Here's a quick tip. A user question came in recently that involved a step of finding the pixels adjacent to foreground pixels in a binary image. Suppose you have a binary mask image, like this one:

bw = imread('circles.png');
imshow(bw)

How can you find all the black pixels in bw that are immediately adjacent to a white pixel? You can do this using imdilate and a logical operation.

Use imdilate to "grow" the mask by one pixel:

bw2 = imdilate(bw, ones(3,3));

Now use a logical operation to find which pixels are white in bw2 but black in bw:

bw3 = bw2 & ~bw;
imshow(bw3)
title('Adjacent pixels')

Or do it in one step:

adjacent_pixels = imdilate(bw, ones(3,3)) & ~bw;

Get the MATLAB code

Published with MATLAB® 7.4

Loren wrote about the 2007 MATLAB World Tour in her Art of MATLAB blog recently. Several other MathWorks developers also attended some of the events.

One message that they brought back to us is that many of the good features we put into our products simply have not yet been discovered by many of our users. That has inspired me to think about a new kind of post on this blog - I think of it as the "Didja Know?" post. (Many of the MATLAB Desktop blog posts are Didja Know? posts.) In my Didja Know? posts, I'll try to highlight useful functions or function options that you might not have discovered yet.

First up: the border replication option of imfilter.

When you perform convolution on an image, one question to consider is: How do you compute the convolution sum for pixels along the image border? By default, imfilter assumes that image pixels "outside" the image are zero. Mathematically this is kind of reasonable, but it tends to produce a dark strip along the edge of the output image. Here's an example to show what I mean:

I = imread('cameraman.tif');
imshow(I)
imcredit('Image courtesy of MIT')

J = imfilter(I, ones(5,5)/25);
imshow(J)

But you can tell imfilter to handle image borders by replicating the border pixel values. That often produces a more desirable result:

K = imfilter(I, ones(5,5)/25, 'replicate');
imshow(K)

There are other border options as well. See the imfilter documentation for more information.

Get the MATLAB code

Published with MATLAB® 7.4

I described in my previous connected component labeling post the algorithm used by bwlabeln. This time I'll talk about the variation used by bwlabel. This variation uses run-length encoding as the first step.

Here's what I mean by run-length encoding. Consider this small binary image matrix:

Such a binary image can be represented as a collection of runs, or sequences of 1s. Working columnwise, this image contains 9 runs:

Each run can be represented by its starting pixel and by the number of pixels in the run, which is called the run length. Hence the term run-length encoding.

For a larger binary image containing larger objects, the number of runs can be much smaller than the number of foreground pixels. But the connectivity analysis can be completely determined based on the runs. In our small example:

This set of connectivity pairs can be resolved into equivalence classes using the technique I described earlier based on dmperm. In this example there are only two equivalence classes, corresponding to the two connected components.

The function bwlabel calls computes a run-length encoding first. As it finds each run, it also determines which runs on the previous column are connected, if any. Then it constructs a sparse adjacency matrix and calls dmperm to compute the equivalence classes. Then it labels each run in the output label matrix according to the equivalence class it belongs to.

This is a very efficient method because it scans each pixel in the input image only once, and also because the adjacency matrix is usually relatively small. It is R-by-R, where R is the number of runs.

You can find this method described in Haralick and Shapiro, Computer and Robot Vision Volume 1, Addison-Wesley, 1992, pp. 40-48. (Note that the equivalence class resolution method used by bwlabel is different than the one described in the book.)

Get the MATLAB code

Published with MATLAB® 7.4

Today I'll continue thinking about plateaus in a DEM and what they mean for calculating pixel flow. Here's a very tiny portion of my DEM array, just northwest of North Pond.

a = [120 118 123; 120 118 125; 119 119 126]

a =

   120   118   123
   120   118   125
   119   119   126

The center pixel doesn't have a downhill neighbor, so I'm calling it a plateau pixel. No matter that it's a very small plateau; we can't compute a pixel flow direction for it, and that will cause problems when we get around to computing upslope area.

One technique used in image processing to deal with plateaus is called the lower complete transformation. A lower complete image is one where the only pixels having no downhill neighors are the pixels belonging to regional minima. Roughly speaking, the lower complete transformation operates by elevating plateau pixels. Pixels further from the plateau edge get elevated more. Of course, pixels uphill from the plateau need to be elevated as well so that they remain uphill. I have read about this technique but have never implemented it myself. (For more information, I recommend P. Soille, Morphological Image Analysis: Principles and Applications, 2nd edition, Springer, 2003, pp. 222-226.)

I'm going to use a different technique, one that's based on the algorithm in the Image Processing Toolbox function roifill. This function fills in pixels in a particular region by smoothly interpolating from the pixels surrounding the region. It does so by solving a sparse linear system based on the notion that, in the filled output, each pixel in the filled region equals the average of its north, east, south, and west neighbors. (It's kind of an interesting method. Maybe that'll be a future blog topic.)

Here's the example (slightly modified) from the roifill documentation:

I = imread('eight.tif');
imshow(I)

c = [222 272 300 270 221 194 222];
r = [21 21 75 121 121 75 21];
hold on
plot(c,r,'r','LineWidth',4)

J = roifill(I,c,r);
imshow(J)

And here's what roifill does to our tiny snippet from the DEM:

mask = true(3, 3);
b = roifill(a, mask)

b =

  120.0000  118.0000  123.0000
  120.0000  120.5000  125.0000
  119.0000  119.0000  126.0000

You can see that the center pixel has been replaced by a new value, and that the center pixel now has a downhill neighbor.

Here's a slightly bigger example, where the two center pixels have no downhill neighbors:

a2 = pond(15:17, 19:22)

a2 =

   116   116   116   120
   117   115   115   117
   119   115   115   116

And the output of roifill

mask2 = true(3, 4);
b2 = roifill(a2, mask2)

b2 =

   116   116   116   120
   117   116   116   117
   119   115   115   116

Well, the two center pixels now have downhill neighbors, but WAIT! The (1,2) pixel no longer has a downhill neighbor! OK, time to 'fess up. I didn't really anticipate that this would be a problem. I'll have to give this some thought. My first reaction is: It'll be OK because we can use the original, unmodified DEM values to compute the pixel flow for the nonplateau pixels. We only need the modified DEM values to compute pixel flow for the plateau pixels. But I'll sleep on it.

Assuming I can resolve that problem successfully, we can use roifill to compute pixel flow directions for the midplateau pixels. It won't help us, though, for the plateau pixels that are part of regional maxima or regional minima. I'll figure out what to do about those next week.

Get the MATLAB code

Published with MATLAB® 7.4

In my previous upslope area post, I showed this graphic of pixel flow around North Pond in Milford, Massachusetts:

Notice that pixels in the pond have no arrows. In fact, the pixelFlow M-file I showed previously is unable to compute a pixel flow direction or magnitude for these pixels because they are in a flat region, or plateau. More specifically, pixelFlow returns a flow direction of NaN for any pixel that has no downhill neighbor.

There is a simple way to find all such pixels using morphological erosion. Erosion with a flat structuring element is equivalent to a local minimum operator. If the minimum value of a pixel and its neighbors equals the pixel itself, then we'll call it a plateau pixel.

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

I'd like to classify plateaus into three types:

The Image Processing Toolbox has functions that find regional maxima and regional minima:

max_plateaus = imregionalmax(pond);
imshow(max_plateaus, 'InitialMagnification', 'fit')
title('Max plateaus')

min_plateaus = imregionalmin(pond);
imshow(min_plateaus, 'InitialMagnification', 'fit')
title('Min plateaus')

We can now identify pixels in plateaus with both uphill and downhill neighbors by using logical operators:

mid_plateaus = plateaus & ~(max_plateaus | min_plateaus);
imshow(mid_plateaus, 'InitialMagnification', 'fit')
title('Mid plateaus')

Next time I'll work on a procedure for calculating flow directions for these different kinds of plateau pixels.

Get the MATLAB code

Published with MATLAB® 7.4