In my August 7th post on upslope area, I showed how to construct and solve the flow matrix to determine the upslope area for every pixel in a digital elevation model (DEM). In addition, the flow matrix can be used to compute both the influence map and the dependence map.

The influence map shows where water starting from a particular pixel will drain. The dependence map shows which uphill pixels drain through a particular pixel. Let's look at how to compute these from the flow matrix.

Recall the form of the flow equation for a particular pixel: The upslope area of a pixel equals 1 plus a weighted fraction of the upslope area of each of its neighbors.

The 1 term comes from assuming each pixel contributes an equal volume to the overall flow across the terrain. Think of it as rain falling equally on each pixel. To compute the influence map, then, we just have to modify our set of equations so that the 1 term appears only for the pixel or pixels of interest. Fortunately, that only affects the right-hand side of the equation, not the flow matrix itself. So solving for the influence map a particular pixel looks almost like solving for upslope area for all pixels.

Here's an example from the Milford, Massachusetts DEM that I've used before:

s = load('milford_ma_dem');
E = s.Zc;
imshow(E, [])
limits = [160 220 215 270];
axis(limits);
hold on
plot(169, 227, '*')
hold off
title('Pixel near the top of Peppercorn Hill')

Find the flow directions and flow matrix using dem_flow and flow_matrix from my upslope area toolbox on MATLAB Central:

R = dem_flow(E);
T = flow_matrix(E, R);

Make a vector for the right-hand side that has a single 1 corresponding to the (227,169) pixel:

[M, N] = size(E);
rhs = zeros(numel(E), 1);
idx = (169-1)*M + 227;
rhs(idx) = 1;

Solve for the influence map:

I = T \ rhs;
I = reshape(I, M, N);
imshow(I,[])
axis(limits);
title('Influence map')

The influence_map function does all this for you. I also included a visualization function (vis_map) that superimposes the influence (or dependence) map transparently, in green, on top of the original DEM. The starting pixels are superimposed in blue. Here's what it looks like for our example:

vis_map(I, E, 227, 169);
axis(limits);

To compute the dependence map, you can't use the flow matrix exactly as is; you have to transpose it. Otherwise, it's very similar to the influence map computation. See the code in dependence_map for the details.

Here's an example that computes the dependence map for an entire set of pixels, not just a single pixel. Specifically, it's the dependence map for North Pond:

plateaus = imerode(E, ones(3, 3)) == E;
pond = bwselect(plateaus, 183, 170);
D = dependence_map(E, T, pond);
vis_map(D, E, pond);
xlim([110 240])
ylim([80 260])

Next time I plan to talk about vectorizing the pixel flow direction calculation.

Get the MATLAB code

Published with MATLAB® 7.5

I saw an application recently where someone wanted to handle NaNs in an image differently depending on whether they were "interior." In other words, a set of connected NaN pixels surrounded by non-NaN pixels are handled one way, and a set of connected NaN pixels that touches the image border are handled another way.

Coincidentally, another person recently commented on the need in some measurement situations to remove objects that touch the image border.

So I thought it would be good time to show how to use the Image Processing Toolbox function imclearborder. This function removes all connected components of a binary image that touch any image border. Here's an example. (To run the example, download aug31.png from here.

BW = imread('aug31.png');
imshow(BW)

BW2 = imclearborder(BW);
imshow(BW2)
title('Objects touching image borders removed')

Now back to the NaN problem. Here's a small example:

f = magic(5);
f(2:5,1:2) = NaN;
f(3:4,4) = NaN

f =

    17    24     1     8    15
   NaN   NaN     7    14    16
   NaN   NaN    13   NaN    22
   NaN   NaN    19   NaN     3
   NaN   NaN    25     2     9

The NaNs on the left make up a group touching the border. The NaNs in the fourth column are an interior group. We can distinguish between the two groups by using imclearboder and some logical operations.

nan_mask = isnan(f)

nan_mask =

     0     0     0     0     0
     1     1     0     0     0
     1     1     0     1     0
     1     1     0     1     0
     1     1     0     0     0

border_nans = nan_mask & ~imclearborder(nan_mask)

border_nans =

     0     0     0     0     0
     1     1     0     0     0
     1     1     0     0     0
     1     1     0     0     0
     1     1     0     0     0

interior_nans = nan_mask & ~border_nans

interior_nans =

     0     0     0     0     0
     0     0     0     0     0
     0     0     0     1     0
     0     0     0     1     0
     0     0     0     0     0

You can use also imclearborder on gray scale images. As the documentation describes, the function "suppresses structures that are lighter than their surroundings and that are connected to the image border." You can see the effect using the rice.png image:

I = imread('rice.png');
imshow(I)

J = imclearborder(I);
imshow(J)

The modified image looks darker. As the documentation notes, for gray scale images imclearborder "tends to reduce the overall intensity level in addition to suppressing border structures." Using the display scaling syntax of imshow helps to see the result clearly:

imshow(J, [])

So what happened to those two dark grains (on the left and on the right)? If we zoom into the original image to look closely at the grain on the right, we can see that there are a couple of slightly darker gray pixels connecting the grain to its neighbor on the right.

imshow(I)
axis([200 256 81 137])

But that grain on the right touches the border. So when imclearborder removes it, it also has the affect of lowering the gray level of its neighbor.

imshow(J, [])
axis([200 256 81 137])

Have you found imclearborder to be useful? Let me know - I'd like to hear about it.

Get the MATLAB code

Published with MATLAB® 7.5

One the measurements provided by regionprops is 'Centroid'. Here's an example of labeling binary objects, computing the centroid of each object, and plotting the centroid location on top of the image.

I = imread('text.png');
imshow(I)

L = bwlabel(I);
s = regionprops(L, 'Centroid');
hold on
for k = 1:numel(s)
    plot(s(k).Centroid(1), s(k).Centroid(2), 'r*')
end
hold off
xlim([0 60])
ylim([0 60])

Last week, reader Daphne asked how to compute the intensity-weighted centroid. That is, if each labeled region corresponds to a region in a gray scale image, how do you compute the centroid weighted by the gray scale pixel values? The easiest way, I think, is to use the both the PixelIdxList and PixelList properties from regionprops. Let's try it with a simple synthetic image (download the image from here):

I = imread('spheres.png');
imshow(I)

Now let's threshold and label it:

bw = I < 255;
bw = imfill(bw, 'holes');
L = bwlabel(bw);
imshow(label2rgb(L, @jet, [.7 .7 .7]))

Next, get the PixelIdxList and PixelList for each labeled region:

s = regionprops(L, 'PixelIdxList', 'PixelList');

Each row of the PixelList for a region contains the x and y coordinates of a pixel in that region.

s(1).PixelList(1:4, :)

ans =

    48    94
    48    95
    48    96
    48    97

We can get the grayscale value of pixels in a region by indexing into I with PixelIdxList. By weighting those values with the x and y coordinates, we can compute the centroid.

idx = s(1).PixelIdxList;
sum_region1 = sum(I(idx));
x = s(1).PixelList(:, 1);
y = s(1).PixelList(:, 2);

xbar = sum(x .* double(I(idx))) / sum_region1

xbar =

   74.6128

ybar = sum(y .* double(I(idx))) / sum_region1

ybar =

   92.5121

Let's do that in a loop and superimpose the centroid locations on the image.

imshow(I)
hold on
for k = 1:numel(s)
    idx = s(k).PixelIdxList;
    pixel_values = double(I(idx));
    sum_pixel_values = sum(pixel_values);
    x = s(k).PixelList(:, 1);
    y = s(k).PixelList(:, 2);
    xbar = sum(x .* pixel_values) / sum_pixel_values;
    ybar = sum(y .* pixel_values) / sum_pixel_values;

    plot(xbar, ybar, '*')
end
hold off

Get the MATLAB code

Published with MATLAB® 7.5

The functions bwlabel and regionprops are very useful for measuring the properties of shapes in a binary image. There are documentation examples and product demos showing how to do this, and I've shown these functions in action several times in this blog.

But sometimes I get questions about how to process pixel values in the "original" gray scale image. In other words, suppose your process is something like this:

This post is about how to do that last step. Let's start by creating some simple sample images to work with.

I = imread('coins.png');
imshow(I)

bw = im2bw(I, graythresh(I));
bw = imfill(bw, 'holes');
imshow(bw)
title('Thresholded, with holes filled')

L = bwlabel(bw);
imshow(label2rgb(L, @jet, [.7 .7 .7]))
title('Label matrix')

The key for using L to process the pixels in I is to use regionprops to get the 'PixelIdxList' for each labeled region.

s = regionprops(L, 'PixelIdxList')

s = 

10x1 struct array with fields:
    PixelIdxList

s is a struct array. s(k).PixelIdxList is a vector of the linear indices of the k-th region. You can use these values to index directly into I. For example, here is the mean pixel value of the 3rd coin:

idx = s(3).PixelIdxList;
mean(I(idx))

ans =

  172.8769

You can also replace object pixel values in I by indexing on the left side of an assignment statement, like this:

% Turn the 4th coin white and the 5th coin black:
J = I;
J(s(4).PixelIdxList) = 255;
J(s(5).PixelIdxList) = 0;
imshow(J)

Here's a loop that replaces the pixel values for each coin by the mean pixel value for that coin:

F = I;
for k = 1:numel(s)
    idx = s(k).PixelIdxList;
    F(idx) = mean(I(idx));
end
imshow(F)

I'll finish with a question that Mat asked in a blog comment: "How can I find the location of the maximum pixel value in each labeled region?" Here's how to find the locations and then plot them on top of the original image. (In addition to 'PixelIdxList', I'll also ask regionprops for 'PixelList', which gives the x and y coordinates of each pixel in the region.)

s = regionprops(L, 'PixelIdxList', 'PixelList');
imshow(I)
hold on
for k = 1:numel(s)
    idx = s(k).PixelIdxList;
    pixels_in_region_k = I(idx);
    [max_value, max_idx] = max(pixels_in_region_k);
    max_location = s(k).PixelList(max_idx, :);
    plot(max_location(1), max_location(2), '*')
end
hold off

Get the MATLAB code

Published with MATLAB® 7.4