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

6 views (last 30 days)
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!

Answers (3)

Wolfgang Schwanghart
Wolfgang Schwanghart on 19 Jul 2012
Hi,
you might be interested in using the function ixneighbors available on the FEX which might make life a little easier for your task.
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.
Regards, W.
  1 Comment
Renee Cammarere
Renee Cammarere on 24 Jul 2012
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.

Sign in to comment.


Image Analyst
Image Analyst on 19 Jul 2012
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);
  1 Comment
Renee Cammarere
Renee Cammarere on 24 Jul 2012
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).

Sign in to comment.


Renee Cammarere
Renee Cammarere on 24 Jul 2012
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.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by