How can I perform a comparison function (of surrounding pixels) on a DEM using a starting pixel?

I would like to write a function that would take as input, a starting pixel in a Digital Elevation Model raster, that would identify all pixels connected with this pixel that have the same (or lower) elevation. I'm not sure if I should use a simple region-growing algorithm for this. Thanks!

Answer 1

[Wolfgang Schwanghart]

Hi,

you might be interested in using the function ixneighbors available on the FEX which might make life a little easier for your task.

http://www.mathworks.de/matlabcentral/fileexchange/16991-neighbor-indexing

Say, you are interested in the pixel with the linear index ix you can proceed as follows

[ic,icd] = ixneighbors(dem,ix);
I = dem(icd)<=dem(ic);
ixn = icd(I);

In addition, you might be interested in TopoToolbox that has various functions for terrain analysis in Matlab.

http://physiogeo.unibas.ch/topotoolbox

Regards, W.

[Renee Cammarere]

I tried this and it worked OK. Thank you so much for your input. Since I need to eventually convert this code to C#, it was brought to my attention that I should probably seek for a way that relies less on Matlab special methods and functions.

Answer 2

[Image Analyst]

If you have the Image Processing Toolbox you could do it in 4 lines without writing your own region growing code by using morphological reconstruction. Here's some untested code just off the top of my head. You could just threshold your image

binaryImage = yourImage <= pixelValue;

Then use imreconstruct to grab only that one region connected to your specified pixel.

% Create a marker image.
markerImage = false(size(yourImage));  % Initialize to all false.
markerImage(row, column) = true; % Set the pixel that you're interested in.
% Extract the connected pixels.
connectedBlob = imreconstruct(binaryImage, markerImage, 8);

[Renee Cammarere]

I have taken note of this method because I may be able to use it eventually. Thank you so much for your input. Since I need to (at some point) convert this code to C#, it was brought to my attention that I should probably seek for a way that relies less on Matlab special methods and functions (although that's the reason we love Matlab).

Answer 3

[Renee Cammarere]

I ended up using a neighbor footprint, and then calculating the neighbor coordinate, and comparing from that . . .

            endneigb=[-1 0; 1 0; 0 -1; 0 1; -1 -1; 1 -1; 1 1; -1 1];
            for j=1:8,
                  % Calculate the neighbor coordinate
                  xn = xcol + neigb(j,1);%column 
                  yn = yrow + neigb(j,2);%row
                  % Check if neighbor is inside image
                    ins=(xn>=1)&&(yn>=1)&&(yn<=size(sink0data,1))&&(xn<=size(sink0data,2));        
                  % Add if inside, not already part of region and if neighbor is < or = maxVal
                  if(ins&&(PondRegion(yn,xn)==0)&&nk_e066_n30data(yn,xn)<=maxElev)
                       % IF <=, then keep xn and yn
                       PondRegion(yn,xn)=current_sink_value;
                  end
             end

I know that it's best to stay away from for loops in Matlab, but again, I'm anticipating having to convert the code later on.