Steve on Image Processing

Today is the final post in my continental divide series. In earlier posts I have used the problem of computing the US continental divide as a vehicle for exploring data import, image display scaling, the watershed transform, label matrices, local and regional minima, binary image manipulation, boundary tracing, and visualization techniques.

For this last post, I thought it would useful to gather together all the computational steps in one short script, starting with data import and finishing with the visualization.

% Import the tiles e10g and f10g.
data_size = [6000 10800 1];  % The data has 1 band.
precision = 'int16=>int16';  % Read 16-bit signed ints into a int16 array.
header_bytes = 0;
interleave = 'bsq';          % Band sequential. Not critical for 1 band.
byte_order = 'ieee-le';
E = multibandread('e10g', data_size, precision, header_bytes, ...
    interleave, byte_order);
F = multibandread('f10g', data_size, precision, header_bytes, ...
    interleave, byte_order);
EF = [E, F];

% Crop the data.
dem = EF(1:4000, 6000:14500);

% Form an ocean mask.
ocean_mask = dem == -500;

% Modify the ocean mask so it contains only the two connected components on
% the left and right side.
[M, N] = size(dem);
ocean_mask = bwselect(ocean_mask, [1 N], [1 1]);

% Modify the DEM so that its only regional minima are the ocean_mask
% pixels.
dem_modified = imimposemin(dem, ocean_mask);

% Compute the watershed transform of the modified DEM.
L = watershed(dem_modified);

% Visualize the result. First, convert the output of watershed to a color
% image.
pacific = [0.2 1 0.2];
atlantic = [0.3 0.3 1];
ridge = [1 0 0];
rgb = label2rgb(L, [pacific; atlantic], ridge);

% Superimpose the colored watershed image over the DEM.
imshow(dem, [-500 3000], 'InitialMagnification', 'fit')
hold on
h = imshow(rgb);
set(h, 'AlphaData', 0.2);

% Trace the watershed ridge line and superimpose it as a thick red line.
b = bwboundaries(L == 0, 4);
b1 = b{1};  % There's only 1 boundary found here - the ridge line.
x = b1(:,2);
y = b1(:,1);
plot(x, y, 'r', 'LineWidth', 2);

That's it! I hope you enjoyed this series.

Do you have your own ideas for fun computation and visualization using DEM data? Please comment.

About this Series

This series explores the problem of computing the location of the continental divide for the United States. The divide separates the Atlantic and Pacific Ocean catchment basins for the North American continent.

As an algorithm development problem, computing the divide lets us explore aspects of data import and visualization, manipulating binary image masks, label matrices, regional minima, and the watershed transform.

Data credit: GLOBE Task Team and others (Hastings, David A., Paula K. Dunbar, Gerald M. Elphingstone, Mark Bootz, Hiroshi Murakami, Hiroshi Maruyama, Hiroshi Masaharu, Peter Holland, John Payne, Nevin A. Bryant, Thomas L. Logan, J.-P. Muller, Gunter Schreier, and John S. MacDonald), eds., 1999. The Global Land One-kilometer Base Elevation (GLOBE) Digital Elevation Model, Version 1.0. National Oceanic and Atmospheric Administration, National Geophysical Data Center, 325 Broadway, Boulder, Colorado 80305-3328, U.S.A. Digital data base on the World Wide Web (URL: http://www.ngdc.noaa.gov/mgg/topo/globe.html) and CD-ROMs.

Get the MATLAB code

Published with MATLAB® 7.8

Today's post in my continental divide series demonstrates the final computational step and shows one procedure for visualizing the Pacific and Atlantic catchment basins.

In my previous post, I modified the DEM (digital elevation model) data set so that its only regional minima were the two oceans.

s = load('continental_divide');
dem_modified = s.dem_modified;
imshow(dem_modified, [-500 3000], 'InitialMagnification', 'fit');

Since the modified DEM has only two regional minima, we expect to get only two catchment basins when we compute the watershed transform.

L = watershed(dem_modified);
max(L(:))

ans =

     2

I'll use label2rgb to display the Pacific catchment basin in green, the Atlantic catchment basin in blue, and the watershed ridge pixels (the continental divide) in red.

pacific = [0.2 1 0.2];
atlantic = [0.3 0.3 1];
ridge = [1 0 0];
rgb = label2rgb(L, [pacific; atlantic], ridge);
imshow(rgb, 'InitialMagnification', 'fit')

In a February blog post, I showed how to overlay one image transparently over another. We can use that technique here to achieve an effect similar to the visualization that Brian Hayes used in his blog post on the continental divide, where the ocean catchment basins are shown shaded with two different colors.

imshow(dem, [-500 3000], 'InitialMagnification', 'fit')
hold on
h = imshow(rgb);
set(h, 'AlphaData', 0.2);

It would be nice if we could make the watershed ridge line itself more visible. For example, if we could use plot, then we could use a thick line. But right now we just have the ridge line represented as zero-valued elements of a raster, instead of the x and y vectors we would need to pass to plot. We could use find to get the coordinates of the ridge pixels, but we wouldn't get them in the right order. One solution is to use bwboundaries to "trace" the zero-valued ridge pixels. Note that I specify 4-connectivity in the call to bwboundaries below. I've found that I get better results using 4-connectivity when tracing a single-pixel-wide stroke that's not a closed loop.

b = bwboundaries(L == 0, 4);
b1 = b{1};  % There's only 1 boundary found here - the ridge line.
x = b1(:,2);
y = b1(:,1);
plot(x, y, 'r', 'LineWidth', 2);

I have one more post planned for this series. Next time I'll gather together all the computational steps in a single, short script.

About this Series

This series explores the problem of computing the location of the continental divide for the United States. The divide separates the Atlantic and Pacific Ocean catchment basins for the North American continent.

As an algorithm development problem, computing the divide lets us explore aspects of data import and visualization, manipulating binary image masks, label matrices, regional minima, and the watershed transform.

Data credit: GLOBE Task Team and others (Hastings, David A., Paula K. Dunbar, Gerald M. Elphingstone, Mark Bootz, Hiroshi Murakami, Hiroshi Maruyama, Hiroshi Masaharu, Peter Holland, John Payne, Nevin A. Bryant, Thomas L. Logan, J.-P. Muller, Gunter Schreier, and John S. MacDonald), eds., 1999. The Global Land One-kilometer Base Elevation (GLOBE) Digital Elevation Model, Version 1.0. National Oceanic and Atmospheric Administration, National Geophysical Data Center, 325 Broadway, Boulder, Colorado 80305-3328, U.S.A. Digital data base on the World Wide Web (URL: http://www.ngdc.noaa.gov/mgg/topo/globe.html) and CD-ROMs.

Get the MATLAB code

Published with MATLAB® 7.8

My colleague Greg Wilson, creator of the course Software Carpentry, proposed a couple of programming problems as prerequisites for the course. The idea is that if you can solve the problems using any language at all, then you're ready for the course.

Here's Problem 1:

"Write a function that takes a list or array of numbers as input and return the largest number that is adjacent to a zero. For example, if the input is:

   [1, 5, 3, 0, 2, 7, 0, 8, 9, 1 0]

the output is 8."

When I saw this problem I immediately fired off an e-mail to Greg with a MATLAB one-liner:

   max(a(imdilate(a == 0, [1 1 1])))

a = [1, 5, 3, 0, 2, 7, 0, 8, 9, 1 0];
max(a(imdilate(a == 0, [1 1 1])))

ans =

     8

Of course, one common problem with one-liners tossed off by "experts" is that they are so often wrong. And my answer above is no exception! If all the nonzero elements adjacent to 0s in the input vector are negative, the above solution would incorrectly return 0.

b = [5 4 -1 0 -2 0 -5 8];
max(b(imdilate(b == 0, [1 1 1])))

ans =

     0

Greg raised two questions in response to my e-mail:

Well, I hope by writing this blog to increase the number of MATLAB and Image Processing Toolbox users who could come up with it on their own. But let me tackle Greg's first question - how many users would understand the solution?

I think the answer is "not very many." That's because my code does not clearly express my intent. This is a common problem with one-liners.

I suspect that most MATLAB programmers go through a phase of being overly enamored of the one-liner. (Some never leave that phase!) Now that I've been in the software development business for quite a while, I've (mostly) lost my fondness for the one-liner. I'll still write code that way if I'm experimenting at the MATLAB command prompt, but in code that's going to ship I'm much more likely to break that expression into parts and give each part a meaningful name. I might write it like this (including the bug fix):

   zero_mask = (a == 0);
   adjacent_to_zero_mask = imdilate(zero_mask, [1 0 1]);
   max_value_adjacent_to_zero = max(a(adjacent_to_zero_mask));

zero_mask = (a == 0)

zero_mask =

     0     0     0     1     0     0     1     0     0     0     1

adjacent_to_zero_mask = imdilate(zero_mask, [1 0 1])

adjacent_to_zero_mask =

     0     0     1     0     1     1     0     1     0     1     0

max_value_adjacent_to_zero = max(a(adjacent_to_zero_mask))

max_value_adjacent_to_zero =

     8

Is this solution understandable? Well, it's more understandable than the first version, but it does assume that the code reader is familiar with a few key concepts:

And maybe you need to be comfortable with the idea of applying image processing operators to things that don't appear to be images, such 1-by-11 vectors. :-)

I have used logical indexing and binary image dilation and erosion in countless ways to do useful things in MATLAB, sometimes with images and sometimes not. I strongly recommend that you add these ideas to your bag of tricks.

Get the MATLAB code

Published with MATLAB® 7.8

Today, in part 5 of my series on computing the US continental divide, we'll look at a technique called minima imposition. Previously we looked at data import and display, the watershed transform and label matrices, different kinds of minima, and manipulating binary images to form an ocean mask.

In part 3, we found that our DEM data set for the continental US contained more than half a million regional minima, each of which was turned into a catchment basin by the watershed transform. We need a way to get just two catchment basins.

Minima imposition is technique that modifies surface heights, pushing some height values upward so that the resulting surface has regional minima only at specified locations. Minima imposition performs this modification while preserving the local height ordering relationships between neighboring pixels.

Here's a one-dimensional example to illustrate.

a = [1 2 3 4 3 4 3 2 1]

a =

     1     2     3     4     3     4     3     2     1

The vector a has three different regional minima locations:

imregionalmin(a)

ans =

     1     0     0     0     1     0     0     0     1

Let's make a mask representing just two locations where we want to have regional minima:

mask = [1 0 0 0 0 0 0 0 1]

mask =

     1     0     0     0     0     0     0     0     1

This mask says we only want regional minima at each end of the vector, so we want the regional minima in the middle to go away somehow. We do this using the Image Processing Toolbox function imimposemin.

imimposemin(a, mask)

ans =

      -Inf    2.0030    3.0030    4.0030    4.0030    4.0030    3.0030    2.0030      -Inf

The new result has regional minima only at the desired locations. The local minima in the middle of the vector has been pushed up to the level of the surround elements. (As I write this, I am noticing that the small fractional bump that imimposemin gives to all the middle elements seems possibly unnecessary. This is something I will investigate at another time.)

In my previous post I created a binary image representing an ocean mask:

s = load('continental_divide');
dem = s.dem_cropped;
ocean_mask = s.ocean_mask;
imshow(ocean_mask, 'InitialMagnification', 'fit')

I can use imimposemin, together with the ocean_mask image, to modify the DEM so that the two oceans are the only regional minima.

dem_modified = imimposemin(dem, ocean_mask);

It doesn't look much different:

imshow(dem_modified, [-500 3000], 'InitialMagnification', 'fit')

But we can verify that now there are only two regional minima:

regional_minima_labeled = bwlabel(imregionalmin(dem_modified));
max(regional_minima_labeled(:))

ans =

     2

imshow(label2rgb(regional_minima_labeled), 'InitialMagnification', 'fit')

I'll add the modified DEM to the data I've been collecting in continental_divide.mat:

save continental_divide dem_modified -append

Next time I'll perform the final step in computing the continental divide and show one way to visualize it using transparency.

About this Series

This series explores the problem of computing the location of the continental divide for the United States. The divide separates the Atlantic and Pacific Ocean catchment basins for the North American continent.

As an algorithm development problem, computing the divide lets us explore aspects of data import and visualization, manipulating binary image masks, label matrices, regional minima, and the watershed transform.

Data credit: GLOBE Task Team and others (Hastings, David A., Paula K. Dunbar, Gerald M. Elphingstone, Mark Bootz, Hiroshi Murakami, Hiroshi Maruyama, Hiroshi Masaharu, Peter Holland, John Payne, Nevin A. Bryant, Thomas L. Logan, J.-P. Muller, Gunter Schreier, and John S. MacDonald), eds., 1999. The Global Land One-kilometer Base Elevation (GLOBE) Digital Elevation Model, Version 1.0. National Oceanic and Atmospheric Administration, National Geophysical Data Center, 325 Broadway, Boulder, Colorado 80305-3328, U.S.A. Digital data base on the World Wide Web (URL: http://www.ngdc.noaa.gov/mgg/topo/globe.html) and CD-ROMs.

Get the MATLAB code

Published with MATLAB® 7.8