Can this be written in a much better way for fast computations?

2 views (last 30 days)
ropesh goyal on 7 Aug 2019
Commented: Guillaume on 7 Aug 2019
coast=magic(14400)
di_max=100;
mx=14400;
my=14400;
for ix=1:mx
for iy=1:my
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for dx=-di_max:di_max
for dy=-di_max:di_max
jx=ix+dx;
jy=iy+dy;
if((jx>=1) && (jx<=mx) && (jy>=1) && (jy<=my))
dd=sqrt(dx^2+dy^2);
coast(jx,jy) = min([coast(jx,jy),dd]);
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
end
ropesh goyal on 7 Aug 2019
Edited: ropesh goyal on 7 Aug 2019
Well thank you for your help.
This is a sub part (coastline treatment) of a very big code (removing absolute error from a DEM). It is not dependent on any other thing other than a matrix 'water' of same size (14400,14400). The values of water varies from 0-150. So, I don't think the resultant will be all 0. In fact, on running like this, all the values are not zero. The thing is it is very time consuming, and I have to repeat the same on almost 60 tiles.

Joel Handy on 7 Aug 2019
Edited: Joel Handy on 7 Aug 2019
The code below will do what you are asking without needing a loop.
coast=magic(14400);
di_max=100;
mx=14400;
my=14400;
[ix dx] = meshgrid(1:mx, -di_max:di_max);
[iy dy] = meshgrid(1:my, -di_max:di_max);
jx = ix+dx;
jy = iy+dy;
dd= sqrt(dx.^2+dy.^2);
mask = ((jx(:)>=1) & (jx(:)<=mx) & (jy(:)>=1) & (jy(:)<=my));
ropesh goyal on 7 Aug 2019
Thank you very much Joel and Guillaume for sparing your time to provide a solution.

Guillaume on 7 Aug 2019
If Joel hypothesis that "You want to look at each pixel which is mostly water, and assign a distance from the pixel to all surrounding pixels that are not water" is true (we still haven't add a proper explanation!), then as I said in a comment, this can be achieved in just one line. This is called the distance transform and is achieved in matlab with bwdist.
If all zeros pixels in coast are to be replaced by their euclidean distance to the nearest non-zero pixel, it's simply:
dist = bwdist(coast ~= 0);
coast(coast == 0) = dist(coast == 0);
Guillaume on 7 Aug 2019
I'm not clear why you'd want a window and I don't understand how you would want it to work with a window. bwdist will find the nearest non-zero pixel however far it is, and do this very fast
Again, as I've repeatedly asked, explain in words what you're trying to achieve. e.g.: I've got a matrix coast representing _____ and a matrix water representing ____ and I want to replace the (non-zeros?zeros?something?) values by the distance to the nearest _____. If there's no nearest _______ in a window of 100x100, then I want ______.