2015 3.016 Problem Set 6

Watershed Transform and Applications in Image Processing

If I gave you this image and told you to count the number of circular blobs, how would you do it?


Source: http://blogs.mathworks.com/steve/2013/11/19/watershed-transform-question-from-tech-support/

You can either do this by hand, or make a program do it for you automatically. I prefer the latter.

But how will the program know which objects are touching? For example, if some blobs are touching each other, the program will interpret it as one giant blob as opposed to multiple circular blobs.

This is a central issue to image processing called image segmentation. The ability to turn an image into a collection of discrete objects is absolutely vital for downstream data analysis of object size, width, and other properties.

To tackle this problem, the watershed transform is a simple yet powerful method for separating objects in an image.  The watershed transform is commonly used in a variety of applications such as identifying grain boundaries, biomedical imaging, and George’s demos...

Image sources: http://rses.anu.edu.au/highlights/figure.php?id=438                                                     http://www.cellxpress.org/manual/images/Ch2_view-segmentation-results.png                                                               George, Recitation 2 on nucleation & growth

Although WatershedComponents is built into Mathematica, we must understand the segmentation mechanism to take advantage of the function.


Let’s take a look at a simple example. “Consider a spherical cow...”

Before we begin, let’s transform the image into a topographic surface, where pixel intensity is represented by altitude.

blobs=Import["http://blogs.mathworks.com/images/steve/2013/blobs.png","PNG"]

In this binary image, each pixel can be represented by a 0 (black) or 1 (white). The following array is a nested list that represents the above image with each list of numbers being a row of pixels.

arrayBlobs=ImageData[blobs]

The above array is hard to visualize, so we’ll rearrange it into a matrix. If you squint a little, you can see the image emerge again where 0’s are the black background and 1’s are the white blobs.

matrixBlobs = arrayBlobs//MatrixForm

Still hard to see? We’ll plot the array to reconstruct the image. ArrayPlot takes a given number -- 0 or 1 -- and represents it as a colored box, which is conveniently just like a pixel.

ArrayPlot[arrayBlobs,ColorRules->{0->Black, 1->White},Frame->None]

Same as our original image? Yup.

ImageCollage[{blobs,ArrayPlot[arrayBlobs,ColorRules->{0->Black, 1->White},Frame->None]}]

Cool. Now that we’ve confirmed that this is our original image, let’s plot our array. Each element in our array (either 0 or 1) represents a height. We had to use the Reverse function due to ListPlot3D reading the array in a different order.

array3Dplot = ListPlot3D[Reverse[arrayBlobs],Mesh->None,ColorFunction->"LakeColors",ViewPoint->Above,ImageSize->Large]

Note: since manipulating the above graphic kills my laptop, I added an additional input viewpoint which takes in a ListPlot3D option (Front, {Left,Top}, etc.).

Now we’ll write a function that converts an image into a 3D surface. It takes in an input image, converts the image into an array of pixel values, and plots these pixel values on the z-axis.

imageto3Dsurface[image_,viewpoint_]:=ListPlot3D[Reverse[ImageData[image]],Mesh->None,ColorFunction->"LakeColors",ViewPoint->viewpoint,ImageSize->Large]

Let’s move on...what IS a watershed (also known as drainage/catchment basin)?

ImageCollage[{Import["http://www.wm.edu/as/kecklab/images/watershed1.jpg","JPG"],Import["http://gillespiegrade5.weebly.com/uploads/9/9/5/2/9952644/6063990.jpg?342", "JPG", ImageSize->Large],Import["http://www.epa.gov/sites/production/files/styles/large/public/2015-04/cr-watershed-map.jpg","JPG"]}]

Get the idea? The river sits at the local minimum of the watershed. However, this isn’t very useful, as we’ll find out later.

Instead, let’s look at some watersheds containing lakes. The lakes sit at the local minima of each watershed.

lakewatershed=Import["http://lakeeriebasin.org/yahoo_site_admin/assets/images/great_lakes_watershed_Ohio_DNR.147102344_std.gif","GIF"]

Now what if we think of our image as a collection of watersheds? First, we’ll manipulate our original image so that each blob has a local minimum.

The distance transform is another important image processing algorithm, but we won’t go into that today. Essentially, the brighest areas represent areas that are furthest from the background, which effectively finds the center-ish of each blob.

ImageCollage[{blobs,distancetransformblobs=DistanceTransform[blobs]//ImageAdjust}]

If we let the local minima of the blobs be represented by the Distance Transform values, let’s take the negative of the DistanceTransform image so that the local minima are black instead of white.

colornegateddistance = ColorNegate[distancetransformblobs]

Then we’ll overlay the DistanceTransform image and the original image to create artificial valleys!

ImageCollage[{overlayedimages=ImageCompose[blobs,{colornegateddistance,0.5}],
lakewatershed}]

Now we’ll use the function created earlier (if you have trypophobia now is a good time to look away)

imageto3Dsurface[image_,viewpoint_]:=ListPlot3D[Reverse[ImageData[image]],Mesh->None,ColorFunction->"LakeColors",ViewPoint->viewpoint,ImageSize->Large]

Note how the following graphics have a PlotRange of 0.5 to 1. This is because after overlaying our images, the lowest pixel intensity is now 0.5, which is grey.

imageto3Dsurface[overlayedimages,Above]
imageto3Dsurface[overlayedimages,{Above,Front}]
imageto3Dsurface[overlayedimages,Front]

Imagine that there is a lake at each of these depressions, and each depression is a watershed. As the water level rises in each separate watershed, we want to keep the watersheds separate from each other by constructing a dam at the watershed line.

Let the following plane represent the water level, which can go from a height of 0 to 1.

waterlevel=Animate[ParametricPlot3D[{x,y,z},{x,0,500},{y,0,500},PlotStyle->Blue,PlotRange->{0,1}],{z,0,1},AnimationRunning->False]

Warning: the following takes a while (7+ seconds). Skip to next section for smoother animation.

Animate[Show[imageto3Dsurface[overlayedimages,Front],ParametricPlot3D[{x,y,z},{x,0,500},{y,0,500},PlotStyle->Blue]],{z,0.5,1},AnimationRunning->False]

Animate[Show[imageto3Dsurface[overlayedimages,Above],ParametricPlot3D[{x,y,z},{x,0,500},{y,0,500},PlotStyle->Blue]],{z,0.5,1},AnimationRunning->False]

Since the above function is too slow to obtain any useful information, I created a movie below. There are some missing frames ($Aborted) because when I was exporting the file, that particular frame took too long to export so Mathematica skipped that frame.

Import["http://pruffle.mit.edu/3.016/Data/frontview.avi","Animation"]
Import["http://pruffle.mit.edu/3.016/Data/aboveview.avi","Animation"]

Note how some watersheds merge after water reaches a certain height. We want to prevent this by building a dam at these watershed lines (lines that separate two watersheds).

As we increase the water level, we are “flooding” the image. Imagine a hole punctured at the bottom of each watershed at the local minimum. Next, imagine plunging the surface into a lake (or, you can think of water rising from the ground).

How do we build these dams at the watershed lines? We need to find the local maxima this time of each region in the direction of two adjacent watersheds and record those positions. From the positions of local maxima, we can reconstruct watershed lines that will segment the different watersheds. Let’s take a look at a general 2D case to illustrate this concept.

Manipulate[Show[Plot[(Sin[x]^4+Cos[x]^3)^2-Cos[x]^2+0.6,{x,0,5},ImageSize->Large,PlotLabel->"Construction of Dams for Arbitrary 2D Curve",LabelStyle->Larger,Filling->Bottom,FillingStyle->LightBrown,PlotStyle->Brown,PlotRange->{0,2},Axes->None],Plot[h,{x,0,5},Filling->Bottom],Graphics[{Brown,Line[{{{1.5708,1.6},{1.5708,h}},{{3.14159,0.6},{3.14159,h}},{{6.28319,0.6},{6.28319,h}},{{4.71239,1.6},{4.71239,h}}}]}]],{h,0,2}]

These watershed lines represent building a dam of infinite height between two watersheds.


Source: http://cmm.ensmp.fr/~serra/cours/pdf/en/ch9en.pdf

The actual algorithm for drawing these segmentation lines uses an ordered/priority queue to flood the pixels of lowest “altitude” (or intensity) first. The rest of the pixels are flooded in an order dependent on their relative altitude with surrounding pixels.

This is exactly what we have modeled already with the rising water! The pixels that define the watershed lines are then stored and combined with the original image to produce the final segmented image.

To save time for covering some cool applications, we'll just use the built-in WatershedComponents function. This function should make more sense now, or else I have failed you. :(

I used MorphologicalComponents on the original image to show what a pre-segmented image looks like to a program. The image on the right is the segmented image.

ImageCollage[{Colorize[MorphologicalComponents[blobs],ColorFunction->"RedBlueTones"],segmented=Colorize[WatershedComponents[overlayedimages],ColorFunction->"RedBlueTones"]}]

Uh-oh. What went wrong?

There are two key problems in this image that illustrate the weaknesses of the Watershed Transform:

1) Oversegmentation - too many lines are drawn to segment objects. This occurs if there are multiple local minima within a single object.

WatershedComponents finds basins at each regional minimum, meaning that the enclosed black dots within clusters of blobs are being interpreted as new watersheds, and consequently, a new objects.

2) Identifying markers - the process to create artificial local minima. This wasn’t too bad for us because the distance transform is decent at finding centers of touching circular objects, but the process of identifying markers is often a manual process.

DistanceTransform was used to create our markers, but it isn’t perfect. Since the objects aren’t perfectly circular, some of the objects have multiple local minima, resulting in multiple watersheds per object and too many watershed lines being drawn.

This leads us to one of the central problems of image segmentation -

Image segmentation is difficult because the processing techniques are often case-specific in order to account for oversegmentation and automatic marker identification.

The watershed transform can’t segment an image on its own. Combined with various other algorithms, the process of image segmentation that uses the watershed transform is called watershed segmentation. It can be used with edge detection and other concepts in mathematical morphology and computer vision to tackle real-life applications.

Let’s take a look at some applications introduced earlier:

Figure 1 							Figure 2
“The relationship between rock fabric and acoustic anisotropy in reservoir sandstones
Figure 1. Perspective view of a cubic portion of the tomogram, segmented into grain space and pore space.
Figure 2. Segmented dataset after application of Euclidean Distance Map and Watershed Transform, which locate grain boundaries and distinguish individual grains (colours are arbitrary).”

Although the authors of this paper used the distance and watershed transforms on a 3D visualization, in many cases, the 3D representation must first be reconstructed with segmented 2D slices stacked on top of each other.


2.4 Cell segmentation - cellXpress: Multimarker watershed segmentation algorithm
Identifies regions of fluorescent DNA to use as markers for the watershed transform. The Gaussian blur is used to create gradient image of the markers within the constraints of the cell areas. Thresholding is used to identify foreground objects vs. background pixels.


Recitation #2 - Growth and Nucleation
“Finally here’s a little visualization of how spherulite growth looks like for homogeneous nucleation.”

George uses the watershed transform using the initial dots as markers for the segmentation process. Each marker grows until it contacts another watershed, and a line is drawn at the watershed line. The different color blocks represent separate watersheds.

“This toy visualization is surprisingly insightful. For example, the straight lines obtained after impingement suggest that the particles were formed at the same instance.” 

Since each marker is growing at the same rate, the watershed lines will be equidistance from the nearest two points - which describes the Voronoi diagram. Interesting connection - the watershed transform on a set of points (“markers”, “seeds”) creates a general Voronoi diagram.

For future studies:

How would you go about segmenting the grains in this SEM image of molybdenum foil?

Source: http://www.lv-em.com/images/web-gallery/SEM/Material-Science/Surfac-Morphology/Grain-Boundries/molybdenum-foil.jpg

Hint:

Morphological Image Analysis: Principles and Applications

The key takeaway is that the watershed transform cannot handle perfect, automated image segmentation on its own. However, understanding of the watershed transform is essential to developing more complex algorithms to tackle segmentation problems in real applications.

Additional Readings

The Morphological Approach to Segmentation: The Watershed Transformation
S. Beucher and F. Meyer

2D Image Segmentation of Metallic Slices
Jarrell W. Waggoner, Jeff Simmons, and Song Wang