# Locating the US continental divide, part 6 – Final computation and visualization 6

Posted by **Steve Eddins**,

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.

- Part 1 - Introduction. Data import and display.
`multibandread`,`imshow`.

- Part 2 - Watershed transform.
`watershed`,`label2rgb`.

- Part 3 - Regional minima.
`imerode`,`imregionalmin`.

- Part 4 - Ocean masks. binary image manipulation,
`bwselect`.

- Part 5 - Minima imposition.
`imimposemin`.

- Part 6 - Computing and visualizing the divide.
`watershed`,`label2rgb`,`bwboundaries`.

- Part 7 - Putting it all together. One script that does everything, from data import through computation and visualization of the divide.

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

**Category:**- Continental divide

## 6 CommentsOldest to Newest

**1**of 6

Dear Steve

could you tell me how can construct with the image overlay, which is widely

used, e.g. in fMRI. Please have a look at these examples,

http://www.sph.sc.edu/comd/rorden/overlay.html

The color areas are the activated areas determined from other images

acquired during activation and these areas are overlaid to gray

images (background).

thank you very much.

**2**of 6

Zhenghui—I’ve written about image display using transparency at least a couple of times. See my 18-Feb-2009 post and my 20-Aug-2008 post. I also noticed that the “Statistical Overlays” page you explains that its overlay method is based Matthew Brett’s “Slice Overlay” MATLAB routines, and it gives a link.

**3**of 6

Hello Steve,

I was inspired by this great post and by your cell segmentation post to try something.

I have a number of images that have features very similar to those cells. They are not cells but the analogy is helpful so i will call them so.

All “cells” have a nucleus, but some cells are isolated, some are partially joined.

I have an example here.

What I would like to do is write code to automatically relabel and separate cells with 2 nuclei from ones with only one nucleus. Then only on the former perform watershed and also use bwboundaries to trace the edge as you do in this post.

I developed a strategy but I am stack with one of the steps.

here’s what i’d do:

1 Perform multiclass Otsu segmentation to label image

2 Relabel differently cells with 2 centers from ones with one centre

3 Convert nucleus of isolated cells to hole and use flood fill operation so they become background

4 Perform watershed (or maybe use distance function) on complement of joined cells (hopefully the only ones left)so that nuclei are minima and cell body is maxima.

5) use bwboundaries to trace the edge

I have below an example result of the Otsu multiclass segmentation from step 1 and it looks promising.

Step 3-5 in theory are not difficult to implement.

But I am really stuck with step 2. Do you have any recommendation/suggestion as to the approach I might use?

Or can you think of a simpler solution altogether.

Thank you very much.

Arni

**4**of 6

Hi Steve,

I left it to stew for a couple of days and think I found a possible solution for step 2.

I could create two separate images using the second otsu class- the green one in my second image, and the highest class, the one with the maxima.

Then I could use the technique you suggested in your post on corresponding labeled objects in two images to find only those cells that have two corresponding highs from the second image.

Makes sense?

I’d still like your comments on my strategy as a whole if you had a chance. Thank you very much. Arni

**5**of 6

Arni—If you can zero-out the segmented center pixels, they will become “holes” in the cells. A cell with two holes has a different Euler number than a cell with only one hole. So try using `regionprops` to get the Euler number of each segmented cell.

**6**of 6

Wow, that’s a great idea. I think mine is implementable but very circuitous. yours is sharp to the point. Thanks Steve!

## Recent Comments