Inpart 4of this series I classified DEM plateaus into:

• Max plateaus
• 分钟高原
• Mid plateaus

Inpart 5I outlined the possibility of usingroifillto modify plateau pixels.

Today I'll combine these ideas with a little morphological processing to compute pixel flow directions for all DEM pixels, except for isolated ones in the center of the min plateaus.

First, find all the plateau pixels that have no downhill neighbors. Remember that we can find these using erosion:

s = load('upslope_sample_dem'); pond = s.Zm(140:280, 160:230); plateaus = imerode(pond, ones(3,3)) == pond; imshow(plateaus,'initialmag','fit') title('Plateau pixels')

Next, find the max plateaus usingimregionalmax:

max_plateaus = imregionalmax(pond);

One way to visualize plateau pixels is to overlay them on the original image using another color. Useimoverlay, a utility function I wrote last year. (See myMarch 28, 2006 postfor details.)

vis1 = imoverlay(mat2gray(pond), max_plateaus, [1 0 0]); imshow(vis1,'initialmag','fit') title('Max plateaus in red')

Thebwmorphshrinkoperation reduces connected components to single points (or, for objects with holes, to a ring).

max_plateaus_shrunk = bwmorph(max_plateaus,'shrink', Inf); vis2 = imoverlay(vis1, max_plateaus_shrunk, [1 1 0]); imshow(vis2,'initialmag','fit') title('Shrunken plateau pixels in yellow')

Bump the shrunken plateau pixels to an elevation slightly above their neighbors. We'll bump up by a fixed constant here. To make a completely general version we'll probably have to calculate the bump based on eps and the maximum plateau height.

pond2 = pond; pond2(max_plateaus_shrunk) = pond(max_plateaus_shrunk) + 0.1;

Now let's do the same thing to the min plateaus, except that we'll slightly lower the shrunken plateau pixels.

min_plateaus = imregionalmin(pond); min_plateaus_shrunk = bwmorph(min_plateaus,'shrink', Inf); pond2(min_plateaus_shrunk) = pond2(min_plateaus_shrunk) - 0.1;

As I described inpart 5, I would like to useroifillhere to interpolate the plateau pixels. However, the logic ofroifilldoesn't quite work out right. Specifically, when you pass in a mask toroifill, it only fills in theinteriorpixels of the mask. You can see this in the code:

dbtyperoifill68:74

68 % Find the perimeter pixels of the specified region; these will 69 % be used to form the boundary conditions for the soap film PDE. 70 perimeter = bwperim(mask); 71 72 % Find the interior pixels; these are the pixels that will be 73 % replaced in the output image. 74 interior = mask & ~perimeter;

For this application I want to fill all the pixels specified by the mask. So I'll do what MATLAB users often do: I'll take an existing M-file and modify it to suit my needs. Starting withroifill, I created a new M-file calledregionfill. In the new function,maskspecifies the pixels to be filled, and the boundary pixels are computed from the mask by dilation instead of by usingbwperim:

dbtyperegionfill27:29

27 % Find the N, E, S, and W pixels adjacent to the pixels to be filled. These 28 % will be used to form the boundary conditions for the soap film PDE. 29 perimeter = imdilate(mask, [0 1 0; 1 1 1; 0 1 0]) & ~mask;

Now I can useregionfillto interpolate the plateau pixels. But not all the plateau pixels. We don't want to fill those that got bumped higher or lower. When those bumped pixels get used byregionfillas boundary conditions, they'll help eliminate the other plateau pixels.

mask = plateaus & ~max_plateaus_shrunk & ~min_plateaus_shrunk; pond2 = regionfill(pond2, mask);

Using both the original DEM (pond) and the modified DEM (pond2), we can now compute pixel flow for every pixel. We use the modified DEM only for plateau pixels (see the last two paragraphs ofpart 5for an explanation).

[M,N] = size(pond); r = zeros(M-2, N-2); s = zeros(M-2, N-2);fori = 1:M-2forj = 1:N-2 [r(i,j), s(i,j)] = pixelFlow(pond,i+1,j+1);ifisnan(r(i,j))% Plateau pixel; use the modified DEM.[r(i,j), s(i,j)] = pixelFlow(pond2,i+1,j+1);endendend

As the images below show, this procedure has successfully computed pixel flow directions for all pixels except those belonging to the shrunken min plateaus.

subplot(1,2,1) imshow(isnan(r)) title('No pixel flow direction'次要情节(1、2、2)imshow (min_plateaus_shrunk(2:——结束1,2:end-1)) title('Min plateaus')

Next time I'll start looking at the actual upslope area computation.

Get the MATLAB code

Published with MATLAB® 7.4