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

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

[Guillaume]

Can this be written in a much better way

Yes, some comments explained what it's meant to do would make it a lot better, not in term of speed but certainly in terms of letting the reading understand your code.

I also suspect that the code doesn't do what it's meant to do as all the code does is fill (very slowly) coast with 0s. So, right now the whole lot is equivalent to:

coast = zeros(14400);

[Joel Handy]

@Guillaume where does coast get filled with 0s? I'm missing that part.

[Guillaume]

For each jx and jy, there's always going an iteration where dx and dy are both 0. Hence, for each jx and jy there will be a time when dd is 0. Of course, the minimum of [cost(jx,jy), 0] is 0, since coast starts with positive integers.

[Joel Handy]

There will be AN interation but not ALL interations. Looks like were missing each others calls between these comments and my answer.

[ropesh goyal]

coast = zeros(mx,my)+9999;
coast(water>=50)=0;
for iy=1:my
    for ix=1:mx
        if(coast(ix,iy)==0)
            iscoast = 0;
            for dx=-1:1
                for dy=-1:1
                    jx=ix+dx;
                    jy=iy+dy;
                    if((jx>=1) && (jx<=mx) && (jy>=1) && (jy<=my))
                        if(coast(jx,jy)~=0)
                            iscoast = 1;
                        end
                    end
                end
            end
            
            if (iscoast==1)
                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
    end
end

[Guillaume]

@joel,

For each ix, and iy, there will always be an iteration where dd is 0. In fact, the dx and dy loop don't matter at all. For each index of coast you'll be trying all dd and the smallest one is always when dd is 0.

This is easily tested on smaller arrays:

coast = magic(100);
di_max = 10;
mx = size(coast, 1); my = size(coast,2);
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

>> all(coast(:) == 0)
ans =
  logical
   1

We'll have to see what the new code does.

@ropesh,

Once again posting code with no explanation of what it's meant to do is not useful. Please explain what you want to do in words and with numeric examples, and then we can work out what the code should be.

I suspect that you're trying to compute something similar to a distance transform. If you say it in words, then we can point you to the function that does it in just one line of code.

[Joel Handy]

@Guillaume. I see it now. Good catch.

For anyone missing it like me, the original code sweeps through every "pixel" and assigns the surrounding pixels a distance from that pixel. Since the pixel itself is being included, the minimum distance of a pixel from a neighboring pixel is always zero.

@ropesh, I kind of figured this is where this was going.

There are some bugs in this code. iscoast is only set if coast(ix,iy) == 0. That means if coast(ix,iy) ~=0 0, iscoast is in a leftover state from the last iteration of the two outer loops.

I also suspect this code is going to overwrite individual pixels more than once making results inconsistent.

Is "water"? is it the same size as coast?

Maybe if you can give us the bigger picture we can make some better suggestions. If I were to venture a guess, you have an image of some sort and you know what percentage of each pixel is water. 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.

[ropesh goyal]

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.

Answer 1

[Joel Handy]

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));
idx = sub2ind(size(coast), jx(mask), jy(mask));
coast(idx) = min(coast(idx),dd(mask));

[Guillaume]

No idea if that produces the same result as the original code (which just fills coast with 0), but note that for the given sizes:

coast(jx(mask),jy(mask))

is about 61 TB of memory.

[Joel Handy]

@Guillaume. I did make a bonehead mistake like I thought I might have. Fixed now hopefully.

I dont see coast filled with zeros though and I dont see how that would hapen in the original code. Cost is initialized to non-zero values and dd is rarely zero either.

[ropesh goyal]

Thank you very much Joel and Guillaume for sparing your time to provide a solution.

Answer 2

[Guillaume]

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);

[ropesh goyal]

Thank you Guillaume.

Can we do this in a window of 100x100 around each pixel? Or kindly suggest if the result will be same either way (using a window and not using the window). I think it will not be the same. Please correct if I am wrong.

[Guillaume]

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 ______.