About twelve years ago, I was implementing the algorithm that became roifill in the Image Processing Toolbox. In this algorithm, a specified set of pixels is replaced (or filled in) such that each of the replaced pixels equals the average of its north, east, south, and west neighbors. I had implemented an iterative algorithm, which worked but was slow. When I asked MATLAB creator and company cofounder Cleve if he knew a better way, he immediately replied, "It looks like a sparse linear system to me. Why don't you just use backslash?"

Head slap! Of course.

Now here I am again, looking at a very similar set of equations and thinking again about using sparse backslash.

In the Tarboton paper, the upslope area calculation is formulated recursively. That is, the upslope area of a pixel equals 1 plus some fraction of the upslope area of its neighboring pixels. (The specific fraction is a function of the flow direction for the neighbors.) So that's a linear system of equations. Let's see how to set up and solve the system in MATLAB.

We'll start by making a test "DEM" to work with. I like using peaks.

E = peaks;
imshow(E, [], 'InitialMagnification', 'fit')

Next, compute the pixel flow directions, using the function dem_flow from my upslope area contribution on the MATLAB Central File Exchange.

R = dem_flow(peaks);

Each row and column of our sparse system will correspond to an individual pixel. It will be convenient to make a matrix, the same size as E, that "numbers" or labels the individual pixels of E. These labels will be used as the sparse matrix row and column coordinates.

pixel_labels = 1:numel(E);
pixel_labels = reshape(pixel_labels, size(E));
pixel_labels(1:4, 1:6)

ans =

     1    50    99   148   197   246
     2    51   100   149   198   247
     3    52   101   150   199   248
     4    53   102   151   200   249

Let's consider north neighbors. Here's the set of pixels with north neighbors:

pixels_with_north_neighbors = pixel_labels(2:end, :);

And here are the corresponding north neighbors:

north_neighbors = pixel_labels(1:end-1, :);

The angle pointing inward from a pixel's north neighbor to itself is -pi/2:

inward_angle = -pi/2;

A pixel's north neighbor has a flow that points in a particular direction. That direction is contained in R. If the flow direction is the same as the inward angle, then the fraction of flow from that neighbor is 1.0. If the flow direction is pointing more than pi/4 radians away from the inward angle, then the fraction of flow from that neighbor is 0.0. For inbetween angular differences, the fraction varies linearly from 1.0 to 0.0.

I have written simple functions to compute angular differences and the corresponding weights. The angular difference computation looks like this:

 d = mod(theta1 - theta2 + pi, 2*pi) - pi;

and the directional weight computation looks like this:

 w = max(1 - abs(d)*4/pi, 0);

With these functions, we can compute the north neighbor weights for all the pixels at once:

w = directional_weight(angular_difference(R(north_neighbors), ...
    inward_angle));

Ignoring the other neighbors for the moment, the upslope area equation looks like this:

  Ap = 1 + wAn

where Ap is the upslope area of the pixel, An is the upslope area of the north neighbor, and w is the directional weight based on the flow direction of the north neighbor. Rewrite the equation this way to see how to construct the linear equation matrix:

 Ap - wAn = 1

So each row of our matrix is going to have a 1 on the diagonal, and it's going to have the negative of the directional flow weights for each of its neighbors. If we used only the northern weights computed above, we could construct the sparse matrix this way:

i = pixels_with_north_neighbors(:);
j = north_neighbors(:);
w = w(:);

% Add diagonal terms
ii = [i ; pixel_labels(:)];
jj = [j ; pixel_labels(:)];
ww = [w ; ones(numel(pixel_labels), 1)];

% And construct the linear system matrix with a call to sparse
T = sparse(ii, jj, ww);
spy(T)

We have to zoom in to clearly see the off-diagonal terms:

xlim([0 20])
ylim([0 20])

I'm going to call T the flow matrix. To finish computing T, we'd have to repeat some of the steps above for each of the seven other neighbor directions, which is way too tedious to reproduce here. I'll just use the flow_matrix function in my MATLAB Central File Exchange submission:

T = flow_matrix(E, R);

Now we can compute the upslope area of each pixel by solving the matrix equation:

  Ta = b

where T is the sparse matrix we just formed, a is the vector of upslope areas for every pixel, and b is a vector of ones.

a = T \ ones(numel(E), 1);
A = reshape(a, size(E));

imshow(A, [], 'InitialMagnification', 'fit')

I have found that using log(A) works well for visualizing an upslope area matrix:

imshow(log(A), [], 'InitialMagnification', 'fit')

The brighter values have the highest upslope areas.

In my next post in this series, I'll revisit the issue of handling plateaus.

Get the MATLAB code

Published with MATLAB® 7.4

So far in this series I've talked about calculating pixel flow directions and handling plateaus, but I haven't yet discussed the actual upslope area calculation. The Tarboton paper presents a recursive formulation of the upslope area calculation (from page 313):

  Procedure DPAREA(i, j)
    if AREA(i, j) is known
      then
        no action
      else
        AREA(i, j) = 1 (the area of a single pixel)
        for each neighbor (location in, jn)
          p = proportion of neighbor (in, jn) that drains to
              pixel (i, j) based on angle
          if (p > 0) then
            call DPAREA(in, jn) (this is the recursive call to
                 calculate the area for the neighbor)
            AREA(i, j) = AREA(i, j) + p x AREA(in, jn)
  Return

I had to read the surrounding text a couple of times to figure out exactly how the author intends to calculate p. Consider the flow direction for a particular neighbor, (in, jn). If the flow direction points directly at pixel (i, j), then the corresponding weight p is 1. In other words, pixel (i, j) gets all of the flow from (in, jn). If the flow direction is pi/4 or greater away from the direction toward (i, j), then the weight p is 0. Pixel (i, j) gets none of the flow from (in, jn). In between an angular difference of 0 and pi/4, the weight p varies proportionally between 1 and 0.

For example, consider the computation of upslope area for pixel number 5 in this set of 9 pixels:

   1 4 7
   2 5 8
   3 6 9

If the flow direction for pixel 2 is zero radians, its flow is pointing directly at pixel 5, so the corresponding weight is 1.0. If the flow direction for pixel 8 is pi/4 radians, then its flow is pointing directly at pixel 4. The corresponding weight for pixel 5 is 0.0. If the flow direction is pi/8 radians, its flow is pointing between pixels 5 and 4, and the weight for the pixel 5 computation is 0.5.

It also took me a few minutes to figure out how to compute the radial difference between two angles, properly taking into account the 2*pi periodicity of angles. (For example, the radial difference between pi/8 and 2*pi - pi/8 should be pi/4.) Here's what I came up with:

   radial_difference = abs(mod(theta1 - theta2 + pi), 2*pi) - pi)

(Does anyone else have a better way to compute this quantity?) To calculate the weight for a given neighbor of pixel (i,j), first determine the angle that points directly from that neighbor to the pixel (i, j). Call that angle theta_c. Then find the flow direction for that neighbor; call it theta_f. Now compute the radial difference:

   radial_difference = abs(mod(theta_c - theta_f + pi), 2*pi) - pi)

And finally compute the weight:

   p = max(1 - (radial_difference * 4 / pi), 0)

So, are we ready to start writing a recursive MATLAB function based on DPAREA above? Nope! I don't actually want to use a recursive formulation at all.

Here's a different way to think about it:

  AREA(i,j) = 1 + p1*AREA(i1,j1) + p2*AREA(i2,j2) + ... + p8*AREA(i8,j8)

If we write down this equation for each pixel in the image, we'll end up with a system of N linear equations in N unknowns. Although the system is very large (N-by-N, where N is the number of pixels), it is also very sparse, because each equation involves no more than nine of the unknowns.

So we are getting very close to calculating upslope area. We "just" have to set up an extremely large sparse system of equations and then solve it.

Next time I'll tackle that.

Get the MATLAB code

Published with MATLAB® 7.4

Didja know? You can fill holes, or pits, in grayscale images by using the Image Processing Toolbox function imfill.

The other day I was rereading the Tarboton paper on upslope area, trying to decide what to do next in my series on that topic. I noticed that the author described filling in "pits" in digital elevation models (DEMs) as a common preprocessing step. He outlined a method based on pixel flow directions, but didn't give details. I realized that this preprocessing step is very easy to do using the Image Processing Toolbox, but that many users might not know about this capability.

Let me first define a little more precisely what I mean by "pit" or "hole," and then I'll show you how to fill it in. It's a little easier to show in one dimension, so let's start there.

Here's a one-dimensional function:

t = linspace(0, 3*pi, 100);
y = t/2 + sin(t);
plot(t,y)

The plot above has a local minimum in its interior. This is a hole, or pit. If we are careful to define our connectivity appropriately for one dimension, we can use imfill with the 'holes' option to fill in the hole.

% Define a one-dimensional connectivity in the horizontal direction.
conn = [0 0 0; 1 1 1; 0 0 0];
y2 = imfill(y, conn, 'holes');

plot(t,y2)

The filled curve has this key property: From any interior point, you can now travel to a boundary minimum without ever going uphill. Here's what it looks for an image:

I = imread('tire.tif');
I2 = imfill(I,'holes');  % the default connectivity is fine
subplot(1,2,1)
imshow(I)
title('Original image')
subplot(1,2,2)
imshow(I2)
title('Filled image')

The hole-filling algorithm in imfill is based on morphological reconstruction (imreconstruct). If you are interested in the theoretical details, see section 6.3.7 of Morphological Image Analysis: Principles and Applications, 2nd ed., Springer, by P. Soille.

Get the MATLAB code

Published with MATLAB® 7.4

L = membrane(1, 100);
imagesc(L), axis equal, colormap(hot)
axis equal

Note that L has negative values.

min(L(:))

ans =

   -0.3338

max(L(:))

ans =

     1

imagesc automatically scales the values of L "into" the colormap, so that the minimum value is displaying using the first colormap color and the maximum value is displayed using the last colormap color. (See my post about scaled indexed images.)

Although Matt's specific suggestion was about getting imoverlay to work with this data, there's a more general concept here about how to save images like the one above to image file formats, such as JPEG or PNG. Image file formats have generally much less flexibility than MATLAB about displaying pixel values as colors. In particular, the most common formats don't have the notion of scaling an arbitrary range of data values into a colormap. In addition, the JPEG format doesn't even have the notion of using a colormap at all.

We could try writing the L values directly into a JPEG file, like this:

imwrite(L,'membrane1.jpg')
imshow membrane1.jpg

But we lost all the color information, and there's no automatic data scaling. All the values of L less than 0 were simply clipped to black in the JPEG file.

So how can we turn L and the hot colormap into something that can be stored in a JPEG file?

We can do it in three steps:

gray = mat2gray(L);
X = gray2ind(gray, 256);
rgb = ind2rgb(X, hot(256));

imwrite(rgb, 'membrane2.jpg')
imshow membrane2.jpg

Now that the membrane image is in RGB form, we can use imoverlay. For example, let's overlay all the negative pixels with a light blue color.

mask = L < 0;
light_blue = [.6 .6 1];
overlay = imoverlay(rgb, mask, light_blue);

imwrite(overlay, 'membrane3.jpg')
imshow membrane3.jpg

Of course, every time you launch MATLAB you see this membrane shape in a more familiar form that looks something like this:

membrane(1,50)
colormap(hot)
camlight
shading interp

The L-shaped membrane has played a fascinating role in the history of MATLAB and The MathWorks. If you are interested in the origins and history of MATLAB, you should check out the Origins of MATLAB video by MATLAB creator and company co-founder Cleve Moler.

Get the MATLAB code

Published with MATLAB® 7.4

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 using imregionalmax:

max_plateaus = imregionalmax(pond);

One way to visualize plateau pixels is to overlay them on the original image using another color. Use imoverlay, a utility function I wrote last year. (See my March 28, 2006 post for details.)

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

The bwmorph shrink operation 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 in part 5, I would like to use roifill here to interpolate the plateau pixels. However, the logic of roifill doesn't quite work out right. Specifically, when you pass in a mask to roifill, it only fills in the interior pixels of the mask. You can see this in the code:

dbtype roifill 68: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 with roifill, I created a new M-file called regionfill. In the new function, mask specifies the pixels to be filled, and the boundary pixels are computed from the mask by dilation instead of by using bwperim:

dbtype regionfill 27: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 use regionfill to 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 by regionfill as 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 of part 5 for an explanation).

[M,N] = size(pond);
r = zeros(M-2, N-2);
s = zeros(M-2, N-2);
for i = 1:M-2
    for j = 1:N-2
        [r(i,j), s(i,j)] = pixelFlow(pond,i+1,j+1);
        if isnan(r(i,j))
            % Plateau pixel; use the modified DEM.
            [r(i,j), s(i,j)] = pixelFlow(pond2,i+1,j+1);
        end
    end
end

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

subplot(1,2,2)
imshow(min_plateaus_shrunk(2:end-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