{"id":258,"date":"2009-05-01T07:00:42","date_gmt":"2009-05-01T11:00:42","guid":{"rendered":"https:\/\/blogs.mathworks.com\/steve\/2009\/05\/01\/continental-divide-2-watershed\/"},"modified":"2019-10-28T15:30:58","modified_gmt":"2019-10-28T19:30:58","slug":"continental-divide-2-watershed","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/steve\/2009\/05\/01\/continental-divide-2-watershed\/","title":{"rendered":"Locating the US continental divide, part 2 – Watershed transform"},"content":{"rendered":"
\r\n <\/introduction>\r\n

This post is the second in my series<\/a> on computing the location of the US continental divide. Today I'm going to look at the watershed transform<\/i>, including how to compute and visualize it using the Image Processing Toolbox functions watershed<\/tt> and label2rgb<\/tt>.\r\n <\/p>\r\n

The watershed transform of an image treats image pixel values as surface heights. It divides image pixels into catchment basins<\/i>, a term borrowed from geography. A catchment basin is the geographical area draining into a river or reservoir. A watershed\r\n ridge<\/i> divides catchment basins.\r\n <\/p>\r\n

The function watershed<\/tt> identifies catchment basins and ridge pixels in an image. Let me illustrate with a one-dimensional example.\r\n <\/p>

H = [4 3 2 3 6 7 5 4 3 3]<\/pre>
\r\nH =\r\n\r\n     4     3     2     3     6     7     5     4     3     3\r\n\r\n<\/pre>
L = watershed(H)<\/pre>
\r\nL =\r\n\r\n     1     1     1     1     1     0     2     2     2     2\r\n\r\n<\/pre>
The output of watershed<\/tt> is telling us that there are two catchment basins in H<\/tt>, corresponding to the elements of L<\/tt> equal to 1 and 2.  The element of L<\/tt> equal to 0 corresponds to the element of H<\/tt> equal to 7; this element is a local maximum that sits between two catchment basins. It is called a ridge pixel<\/i>. We call L<\/tt> a label matrix<\/i>. The toolbox functions watershed<\/tt>, bwlabel<\/tt>, and bwdist<\/tt> all produce output in this form.\r\n   <\/p>\r\n   
Normally in image processing, the watershed transform is used in situations where the pixel values do not literally represent\r\n      height values.  (See my MATLAB News & Notes article on using the watershed transform for image segmentation.)\r\n   <\/p>\r\n   
For this series, though, we really do have an image whose pixel values represent height, and the watershed transform is what\r\n      we want to use to compute the continental divide. Essentially, we're looking for the Pacific Ocean and Atlantic Ocean catchment\r\n      basins for the continental US.\r\n   <\/p>\r\n   
So let's try computing the watershed transform of our DEM (digital elevation model). I saved the DEM to a MAT-file in my previous\r\n      post in this series.\r\n   <\/p>
s = load('continental_divide'<\/span>);\r\ndem = s.dem_cropped;\r\ndisplay_range = [-500 3000];\r\nimshow(dem, display_range, 'InitialMagnification'<\/span>, 'fit'<\/span>);<\/pre> 
L = watershed(dem);<\/pre>
The toolbox function label2rgb<\/tt> is useful for visualizing label matrices. It converts a label matrix to an image, with different colors assigned to the different\r\n      catchment basins.  I'll use label2rgb<\/tt> to visualize L<\/tt> using the jet<\/tt> colormap, with ridge pixels displayed as white and catchment basin colors assigned randomly.\r\n   <\/p>
rgb = label2rgb(L, 'jet'<\/span>, 'w'<\/span>, 'shuffle'<\/span>);\r\nimshow(rgb, 'InitialMagnification'<\/span>, 'fit'<\/span>)<\/pre> 
Oh boy, what a mess! What does that even mean? Let's zoom in much closer to see what's really going on.<\/p>
axis([5580 5670 2060 2100])<\/pre> 
It looks we have a very large number of catchment basins, most of which contain only a few pixels.  How many catchment basins\r\n      are there? Just look at the maximum value of L<\/tt>.\r\n   <\/p>
max(L(:))<\/pre>
\r\nans =\r\n\r\n      540976\r\n\r\n<\/pre>
So there are more than half a million distinct catchment basins detected by watershed<\/tt>.\r\n   <\/p>\r\n   
The next posts in this series will examine why there are so many, and how to get down to just two.<\/p>\r\n   
About this Series<\/i><\/p>\r\n   
This series<\/a> explores the problem of computing the location of the continental divide for the United States.  The divide separates the\r\n      Atlantic and Pacific Ocean catchment basins for the North American continent.\r\n   <\/p>\r\n   
As an algorithm development problem, computing the divide lets us explore aspects of data import and visualization, manipulating\r\n      binary image masks, label matrices, regional minima, and the watershed transform.\r\n   <\/p>\r\n   
\r\n