Steve on Image Processing

August 23rd, 2007

Upslope area - handling plateaus, revisited

Earlier in this series I have discussed the problem of plateaus. Specifically, how do you assign a flow direction to a DEM pixel that has no downhill neighbors? I proposed in part 5 that the problem could be addressed by using the interpolation technique in roifill to modify plateau values so that they are no longer equal, and then recomputing the flow direction based on the modified values.

OK, true confessions time—I could never get this idea to really work! Here's a small set of DEM values to illustrate the typical problem:

The pixels labeled A and B in the right column have no downhill neighbors. Modifying those values using roifill so they were no longer flat had the effect of raising pixel A above the value of pixel C. As a result, my recomputed flow directions had the amusing property that pixel C flowed to pixel A, and pixel A flowed to pixel C!

I discovered this problem while trying to understand why, in some cases, I was getting flow matrices that were singular.

So it was back to the drawing board. I eventually found a good solution in F. Meyer, "Skeletons and Perceptual Graphs," Signal Processing 16 (1989) 335-363. This paper has a couple of algorithms for determining ordering relationships among pixels on a plateau. One of the algorithms, called arrowing, turned to be suitable for computing flow matrix weights for plateau pixels. (Meyer credits Maisonneuve for the original idea.)

First scan. During the first scan every pair of neighbors (x,y) is considered. If the pixel x is higher than the pixel y, it is given an arrow towards y. During this scan all descending borders of all plateaus receive arrows.

Next scans. Until convergence, the following operation is repeated: if x and y are two neighbours sharing the same altitude, and if x possesses an arrow but not y, then y is given an arrow pointing towards x.

–pages 360-361

Here's an example showing how to compute flow directions using arrowing. Below, on the left, is a 4-by-4 matrix of DEM pixels. Some of these DEM pixels have no downhill neighbors. These pixels are indicated via NaNs in the matrix on the right.

With our earlier method for calculating flow directions, the descending borders of plateaus already have flow directions computed for them, which means we can skip the first scan step of the arrowing algorithm. In each subsequent scan, we identify the pixels that do not have a flow direction yet, and which have equal-altitude neighbors that do. Here is the first set of such pixels:

We'll say that each shaded pixel flows equally to all of its equal-altitude neighbors that already have flow directions established, like this:

On the second scan, we find only one pixel without a flow direction, and which has equal-altitude neighbors that do have flow directions:

And so we assign flow directions like this:

The third scan finds no pixels satisfying the criteria, so we know we're done.

Revisiting the example I started with, here's how the arrowing procedure assigns flow directions to pixels A and B:

The MATLAB code is a little too involved to show in detail here. If you care to see exactly how it works, look at the plateau_flow_weights subfunction of flow_matrix.m in my upslope area code on the MATLAB Central File Exchange.

2 Responses to “Upslope area - handling plateaus, revisited”

  1. Urs (us) Schwarz replied on :

    warning: OT

    an all-new mug-shot - ready for image processing…
    :-)

    urs

  2. Steve replied on :

    Urs—Yep. I hated that old picture. :-)

Leave a Reply

Wrap code fragments inside <pre> tags, like this:

<pre class="code">
a = magic(3);
sum(a)
</pre>

If you have a "<" character in your code, either follow it with a space or replace it with "&lt;" (including the semicolon).


Steve Eddins manages the Image & Geospatial development team at The MathWorks and coauthored Digital Image Processing Using MATLAB. He writes here about image processing concepts, algorithm implementations, and MATLAB.

  • Sana: hi steve, could you explain to me how i would be able to use the dir function, to do a loop through a directory...
  • Nishtha: Sir, I have preprocessed the image in following steps: [1] adaptive histogram equalization [2] thresholding...
  • Kristof: I also strongly support the idea. I have just recently bumped into the problem that im2single was not...
  • Steve: David—I’ m glad you found it useful!
  • David Lalejini: I found your example very useful for finding connected nodes in a large set of input pairs. I start...
  • tommy: Dear Steve, I have a question,please if you are kind to help me regarding the accumulator array dimensions of...
  • Steve: Abc—I don’t know how to distinguish the faces. You might try posting your question in the MATLAB...
  • Manju: well if we have a few ovals within each other like in a cell how do we measure the distance from the center...
  • Steve: Manju—What do you mean? How is each region defined?
  • Manju: if we have 2-3 regions within each other how do we measure the regions of each one?

These postings are the author's and don't necessarily represent the opinions of The MathWorks.