|
// Suppose we have an image of many round or oval objects, such as a
// microscope image of blood cells. After thresholding the image, we
// end up with a binary image like this:
xset('auto clear', 'on');
a = gray_imread(SIPDIR + 'images/disks.bmp');
imshow(a,2);
//
// We want to make the computer count the number of cells in the image,
// but there are some circles that are overlapping, thus forming a single
// connected component. Watershed is classically used for separating
// mingled objects like these.
//
// First, calculate the distance transform:
a = 1-a;
d = bwdist(a);
d = normal(sqrt(d),255); // normalize it to 0-255 range
imshow(d+1, hotcolormap(256));
//
// The latter command shows the distance transform in shades from
// black to red to yellow to white. The brighter the color, the
// greater the distance of a point to the background.
// If you have the ENRICO toolbox, you can nicely plot the distance
// transform in 3D using sadesurf. ENRICO is not necessary for this
// example, but anyway you may download it at:
// http://www.weizmann.ac.il/~fesegre/scistuff.html
//
// Note that the peaks of the distance transform are in the middle of
// each blob. The idea is to run watershed segmentation using these
// peaks as markers. For this, we invert the distance transform so
// that the peaks become the regional minima:
d = 255 - d;
imshow(d+1, hotcolormap(256));
// Now we "and" the distance transform with the original image, so
// that the background remains dark.
d = d .* a;
imshow(d+1,hotcolormap(256));
// Finally, run watershed segmentation. It automatically detects the
// regional minima for us:
w = watershed(d);
imshow(w, rand(256,3));
// 'w' is an image with a unique number for each watershed region.
// The imshow with a random colormap displays each region with a
// unique arbitrary color. Note how the regions were correctly separated by
// watershed, except for the hardest cases. It is extremely easy to
// count the number of regions:
n = maxi(w) - 1 // 26 regions minus the background
// The computer found 25 regions, but there are 20, an error of about 20%
// Let's improve this result. In the cases with many overlapping
// circles, the result would be perfect if it weren't for the small
// spurious regions. These are much smaller than the circles,
// so we can safely eliminate the regions with less than 100 pixels:
w2 = w;
for i=1:n
f = find(w==i); // coordinates of i-th region
if size(f,'*') < 100
w2(f) = 26; // merge small regions with the background
end
end
imshow(w2, rand(256,3)); // note how the small regions are gone
// Now we count again, using a different way:
n = size( unique(w2), '*') - 1 // subtract 1 for the background
// Now it's 100% correct! We have an authomatic method wich is surprisingly
// robust. This is specially useful to deal with bigger images in large
// ammount.
//
// Enjoy!
//
// TIP #1: Another way to improve the results is to do a median
// filtering of the distance transform. This will remove many spurious
// minima. Use mogrify(img, ['-median', '2']). Slight Gaussian smoothing also
// works well.
//
// TIP #2: for grayscale image segmentation, calculate the image
// gradient before watershed. Use edge(img,'sobel',-1) for this.
//
// TIP #3: use xgetpixel in a for loop (or something similar) to select
// the seed pixels to be used with watershed segmentation.
xset('auto clear', 'off');
|