發現兩篇數圖像處理的文章。
在此收藏，留待以後細看。
基于Matlab的标記分水嶺分割算法
(2011-05-21 13:57:31)
轉載▼
标簽： matlab标記分水嶺分割代碼it | 分類：技術文章 |
源位址：http://blog.sina.com.cn/s/blog_725866260100rz7x.html
1 綜述
Separating touching objects in an image is one ofthe more difficult image processing operations. The watershedtransform is often applied to this problem. The watershed transformfinds "catchment basins"(集水盆) and "watershed ridge lines"(山脊線) inan image by treating it as a surface where light pixels are highand dark pixels are low.
如果圖像中的目标物體是連接配接在一起的，則分割起來會更困難，分水嶺分割算法經常用于處理這類問題，通常會取得比較好的效果。分水嶺分割算法把圖像看成一幅“地形圖”，其中亮度比較強的區域像素值較大，而比較暗的區域像素值較小，通過尋找“彙水盆地”和“分水嶺界限”，對圖像進行分割。
Segmentation using the watershed transform worksbetter if you can identify, or "mark," foreground objects andbackground locations. Marker-controlled watershed segmentationfollows this basic procedure:
直接應用分水嶺分割算法的效果往往并不好，如果在圖像中對前景對象和背景對象進行标注差別，再應用分水嶺算法會取得較好的分割效果。基于标記控制的分水嶺分割方法有以下基本步驟：
1. Compute a segmentation function. This is animage whose dark regions are the objects you are trying tosegment.
1.計算分割函數。圖像中較暗的區域是要分割的對象。
2. Compute foreground markers. Theseare connected blobs of pixels within each of the objects.
2.計算前景标志。這些是每個對象内部連接配接的斑點像素。
3. Compute background markers. Theseare pixels that are not part of any object.
3.計算背景标志。這些是不屬于任何對象的像素。
4. Modify the segmentation functionso that it only has minima at the foreground and background markerlocations.
4.修改分割函數，使其僅在前景和後景标記位置有極小值。
5. Compute the watershed transformof the modified segmentation function.
5.對修改後的分割函數做分水嶺變換計算。
Use by Matlab Image Processing Toolbox
使用MATLAB圖像處理工具箱
注：期間用到了很多圖像處理工具箱的函數，例如fspecial、imfilter、watershed、label2rgb、imopen、imclose、imreconstruct、imcomplement、imregionalmax、bwareaopen、graythresh和imimposemin函數等。
2 步驟
Step 1: Read in the Color Image and Convert it toGrayscale
第一步：讀入彩色圖像，将其轉化成灰階圖像
clc; clear all; close all;
rgb = imread('pears.png');
if ndims(rgb) == 3
I =rgb2gray(rgb);
else
I = rgb;
end
figure('units', 'normalized', 'position', [0 0 11]);
subplot(1, 2, 1); imshow(rgb); title('原圖');
subplot(1, 2, 2); imshow(I); title('灰階圖');
Step 2: Use the Gradient Magnitude as theSegmentation Function
第2步：将梯度幅值作為分割函數
Use the Sobel edge masks, imfilter, and some simplearithmetic to compute the gradient magnitude. The gradient is highat the borders of the objects and low (mostly) inside theobjects.
使用Sobel邊緣算子對圖像進行水準和垂直方向的濾波，然後求取模值，sobel算子濾波後的圖像在邊界處會顯示比較大的值，在沒有邊界處的值會很小。
hy = fspecial('sobel');
hx = hy';
Iy = imfilter(double(I), hy, 'replicate');
Ix = imfilter(double(I), hx, 'replicate');
gradmag = sqrt(Ix.^2 + Iy.^2);
figure('units', 'normalized', 'position', [0 0 11]);
subplot(1, 2, 1); imshow(I,[]), title('灰階圖像')
subplot(1, 2, 2); imshow(gradmag,[]),title('梯度幅值圖像')
Can you segment the image by using the watershedtransform directly on the gradient magnitude?
可否直接對梯度幅值圖像使用分水嶺算法？
L = watershed(gradmag);
Lrgb = label2rgb(L);
figure('units', 'normalized', 'position', [0 0 11]);
subplot(1, 2, 1); imshow(gradmag,[]),title('梯度幅值圖像')
subplot(1, 2, 2); imshow(Lrgb);title('梯度幅值做分水嶺變換')
No. Without additional preprocessing such as themarker computations below, using the watershed transform directlyoften results in "oversegmentation."
直接使用梯度模值圖像進行分水嶺算法得到的結果往往會存在過度分割的現象。是以通常需要分别對前景對象和背景對象進行标記，以獲得更好的分割效果。
Step 3: Mark the Foreground Objects
第3步：标記前景對象
A variety of procedures could be applied here tofind the foreground markers, which must be connected blobs ofpixels inside each of the foreground objects. In this exampleyou'll use morphological techniques called"opening-by-reconstruction" and "closing-by-reconstruction" to"clean" up the image. These operations will create flat maximainside each object that can be located using imregionalmax.
有多種方法可以應用在這裡來獲得前景标記，這些标記必須是前景對象内部的連接配接斑點像素。這個例子中，将使用形态學技術“基于開的重建”和“基于閉的重建”來清理圖像。這些操作将會在每個對象内部建立機關極大值，使得可以使用imregionalmax來定位。
開運算和閉運算：先腐蝕後膨脹稱為開；先膨脹後腐蝕稱為閉。開和閉這兩種運算可以除去比結構元素小的特定圖像細節，同時保證不産生全局幾何失真。開運算可以把比結構元素小的突刺濾掉，切斷細長搭接而起到分離作用；閉運算可以把比結構元素小的缺口或孔填充上，搭接短的間隔而起到連接配接作用。
Opening is an erosion followed by a dilation, whileopening-by-reconstruction is an erosion followed by a morphologicalreconstruction. Let's compare the two. First, compute the openingusing imopen.
開操作是腐蝕後膨脹，基于開的重建（基于重建的開操作）是腐蝕後進行形态學重建。下面比較這兩種方式。首先，用imopen做開操作。
se = strel('disk', 20);
Io = imopen(I, se);
figure('units', 'normalized', 'position', [0 0 11]);
subplot(1, 2, 1); imshow(I, []); title('灰階圖像');
subplot(1, 2, 2); imshow(Io), title('圖像開操作')
Next compute the opening-by-reconstruction usingimerode and imreconstruct.
接下來，通過腐蝕後重建來做基于開的重建計算。
Ie = imerode(I, se);
Iobr = imreconstruct(Ie, I);
figure('units', 'normalized', 'position', [0 0 11]);
subplot(1, 2, 1); imshow(I, []); title('灰階圖像');
subplot(1, 2, 2); imshow(Iobr, []),title('基于開的重建圖像')
Following the opening with a closing can remove thedark spots and stem marks. Compare a regular morphological closingwith a closing-by-reconstruction. First try imclose:
開操作後，接着進行閉操作，可以移除較暗的斑點和枝幹标記。對比正常的形态學閉操作和基于閉的重建操作。首先，使用imclose：
Ioc = imclose(Io, se);
Ic = imclose(I, se);
figure('units', 'normalized', 'position', [0 0 11]);
subplot(2, 2, 1); imshow(I, []); title('灰階圖像');
subplot(2, 2, 2); imshow(Io, []);title('開操作圖像');
subplot(2, 2, 3); imshow(Ic, []);title('閉操作圖像');
subplot(2, 2, 4); imshow(Ioc, []),title('開閉操作');
Now use imdilate followed by imreconstruct. Noticeyou must complement the image inputs and output of imreconstruct. IM2 = imcomplement(IM) computes thecomplement(補集) of the image IM. IM can be a binary, intensity, orRGB image. IM2 has the same class and size as IM.
現在使用imdilate，然後使用imreconstruct。注意必須對輸入圖像求補，對imreconstruct輸出圖像求補。IM2=imcomplement(IM)計算圖像IM的補集。IM可以是二值圖像，或者RGB圖像。IM2與IM有着相同的資料類型和大小。
Iobrd = imdilate(Iobr, se);
Iobrcbr = imreconstruct(imcomplement(Iobrd),imcomplement(Iobr));
Iobrcbr = imcomplement(Iobrcbr);
figure('units', 'normalized', 'position', [0 0 11]);
subplot(2, 2, 1); imshow(I, []); title('灰階圖像');
subplot(2, 2, 2); imshow(Ioc, []);title('開閉操作');
subplot(2, 2, 3); imshow(Iobr, []);title('基于開的重建圖像');
subplot(2, 2, 4); imshow(Iobrcbr, []),title('基于閉的重建圖像');
As you can see by comparing Iobrcbr with Ioc,reconstruction-based opening and closing are more effective thanstandard opening and closing at removing small blemishes withoutaffecting the overall shapes of the objects. Calculate the regionalmaxima of Iobrcbr to obtain good foreground markers.
通過比較Iobrcbr和loc可以看到，在移除小污點同時不影響對象全局形狀的應用下，基于重建的開閉操作要比标準的開閉重建更加有效。計算Iobrcbr的局部極大來得到更好的前景标記。
fgm = imregionalmax(Iobrcbr);
figure('units', 'normalized', 'position', [0 0 11]);
subplot(1, 3, 1); imshow(I, []); title('灰階圖像');
subplot(1, 3, 2); imshow(Iobrcbr, []);title('基于重建的開閉操作');
subplot(1, 3, 3); imshow(fgm, []);title('局部極大圖像');
To help interpret the result, superimpose(疊加) the foregroundmarker image on the original image.
為了幫助了解這個結果，疊加前景标記到原圖上。
It1 = rgb(:, :, 1);
It2 = rgb(:, :, 2);
It3 = rgb(:, :, 3);
It1(fgm) = 255; It2(fgm) = 0; It3(fgm) = 0;
I2 = cat(3, It1, It2, It3);
figure('units', 'normalized', 'position', [0 0 11]);
subplot(2, 2, 1); imshow(rgb, []);title('原圖像');
subplot(2, 2, 2); imshow(Iobrcbr, []);title('基于重建的開閉操作');
subplot(2, 2, 3); imshow(fgm, []);title('局部極大圖像');
subplot(2, 2, 4); imshow(I2);title('局部極大疊加到原圖像');
Notice that some of the mostly-occluded andshadowed objects are not marked, which means that these objectswill not be segmented properly in the end result. Also, theforeground markers in some objects go right up to the objects'edge. That means you should clean the edges of the marker blobs andthen shrink them a bit. You can do this by a closing followed by anerosion.
注意到大多閉塞處和陰影對象沒有被标記，這就意味着這些對象在結果中将不會得到合理的分割。而且，一些對象的前景标記會一直到對象的邊緣。這就意味着應該清理标記斑點的邊緣，然後收縮它們。可以通過閉操作和腐蝕操作來完成。
se2 = strel(ones(5,5));
fgm2 = imclose(fgm, se2);
fgm3 = imerode(fgm2, se2);
figure('units', 'normalized', 'position', [0 0 11]);
subplot(2, 2, 1); imshow(Iobrcbr, []);title('基于重建的開閉操作');
subplot(2, 2, 2); imshow(fgm, []);title('局部極大圖像');
subplot(2, 2, 3); imshow(fgm2, []);title('閉操作');
subplot(2, 2, 4); imshow(fgm3, []);title('腐蝕操作');
This procedure tends to leave some stray isolated pixels thatmust be removed. You can do this using bwareaopen, which removes all blobs that havefewer than a certain number of pixels. BW2 = bwareaopen(BW,P)removes from a binary image all connected components (objects) thathave fewer than P pixels, producing another binary image, BW2.
這個過程将會留下一些偏離的孤立像素，應該移除它們。可以使用bwareaopen，用來移除少于特定像素個數的斑點。BW2 =bwareaopen(BW,P)從二值圖像中移除是以少于P像素值的連通塊，得到另外的二值圖像BW2。
fgm4 = bwareaopen(fgm3, 20);
It1 = rgb(:, :, 1);
It2 = rgb(:, :, 2);
It3 = rgb(:, :, 3);
It1(fgm4) = 255; It2(fgm4) = 0; It3(fgm4) = 0;
I3 = cat(3, It1, It2, It3);
figure('units', 'normalized', 'position', [0 0 11]);
subplot(2, 2, 1); imshow(I2, []);title('局部極大疊加到原圖像');
subplot(2, 2, 2); imshow(fgm3, []);title('閉腐蝕操作');
subplot(2, 2, 3); imshow(fgm4, []);title('去除小斑點操作');
subplot(2, 2, 4); imshow(I3, []);title('修改局部極大疊加到原圖像');
Step 4: Compute Background Markers
Now you need to mark the background. In thecleaned-up image, Iobrcbr, the dark pixels belong to thebackground, so you could start with a thresholding operation.
第4步：計算背景标記
現在，需要标記背景。在清理後的圖像Iobrcbr中，暗像素屬于背景，是以可以從門檻值操作開始。
bw = im2bw(Iobrcbr, graythresh(Iobrcbr));
figure('units', 'normalized', 'position', [0 0 11]);
subplot(1, 2, 1); imshow(Iobrcbr, []);title('基于重建的開閉操作');
subplot(1, 2, 2); imshow(bw, []);title('門檻值分割');
The background pixels are in black, but ideally we don't wantthe background markers to be too close to the edges of the objectswe are trying to segment. We'll "thin" the background by computingthe "skeleton by influence zones", or SKIZ, of the foreground ofbw. This can be done by computing the watershed transform of thedistance transform of bw, and then looking for the watershed ridgelines (DL == 0) of the result. D = bwdist(BW) computes theEuclidean distance transform of the binary image BW. For each pixelin BW, the distance transform assigns a number that is the distancebetween that pixel and the nearest nonzero pixel of BW. bwdist uses the Euclidean distance metric bydefault. BW can have any dimension. D is the same size as BW.
背景像素在黑色區域，但是理想情形下，不必要求背景标記太接近于要分割的對象邊緣。通過計算“骨架影響範圍”來“細化”背景，或者SKIZ，bw的前景。這個可以通過計算bw的距離變換的分水嶺變換來實作，然後尋找結果的分水嶺脊線（DL==0）。D=bwdist(BW)計算二值圖像BW的歐幾裡得矩陣。對BW的每一個像素，距離變換指定像素和最近的BW非零像素的距離。bwdist預設使用歐幾裡得距離公式。BW可以由任意維數，D與BW有同樣的大小。
D = bwdist(bw);
DL = watershed(D);
bgm = DL == 0;
figure('units', 'normalized', 'position', [0 0 11]);
subplot(2, 2, 1); imshow(Iobrcbr, []);title('基于重建的開閉操作');
subplot(2, 2, 2); imshow(bw, []);title('門檻值分割');
subplot(2, 2, 3); imshow(label2rgb(DL), []);title('分水嶺變換示意圖');
subplot(2, 2, 4); imshow(bgm, []);title('分水嶺變換脊線圖');
Step 5: Compute the Watershed Transform of theSegmentation Function.
The function imimposemin can be used to modify animage so that it has regional minima only in certain desiredlocations. Here you can use imimposemin to modify the gradientmagnitude image so that its only regional minima occur atforeground and background marker pixels.
第5步：計算分割函數的分水嶺變換
函數imimposemin可以用來修改圖像，使其隻是在特定的要求位置有局部極小。這裡可以使用imimposemin來修改梯度幅值圖像，使其隻在前景和後景标記像素有局部極小。
gradmag2 = imimposemin(gradmag, bgm | fgm4);
figure('units', 'normalized', 'position', [0 0 11]);
subplot(2, 2, 1); imshow(bgm, []);title('分水嶺變換脊線圖');
subplot(2, 2, 2); imshow(fgm4, []);title('前景标記');
subplot(2, 2, 3); imshow(gradmag, []);title('梯度幅值圖像');
subplot(2, 2, 4); imshow(gradmag2, []);title('修改梯度幅值圖像');
Finally we are ready to compute the watershed-basedsegmentation.
最後，可以做基于分水嶺的圖像分割計算。
Step 6: Visualize the Result
One visualization technique is to superimpose theforeground markers, background markers, and segmented objectboundaries on the original image. You can use dilation as needed tomake certain aspects, such as the object boundaries, more visible.Object boundaries are located where L == 0.
第6步：檢視結果
一個可視化技術是疊加前景标記、背景标記、分割對象邊界到初始圖像。可以使用膨脹來實作某些要求，比如對象邊界，更加清晰可見。對象邊界定位于L==0的位置。
It1 = rgb(:, :, 1);
It2 = rgb(:, :, 2);
It3 = rgb(:, :, 3);
fgm5 = imdilate(L == 0, ones(3, 3)) | bgm |fgm4;
It1(fgm5) = 255; It2(fgm5) = 0; It3(fgm5) = 0;
I4 = cat(3, It1, It2, It3);
figure('units', 'normalized', 'position', [0 0 11]);
subplot(1, 2, 1); imshow(rgb, []);title('原圖像');
subplot(1, 2, 2); imshow(I4, []);title('标記和對象邊緣疊加到原圖像');
This visualization illustrates how the locations ofthe foreground and background markers affect the result. In acouple of locations, partially occluded darker objects were mergedwith their brighter neighbor objects because the occluded objectsdid not have foreground markers.
可視化說明了前景和後景标記如何影響結果。在幾個位置，部分的較暗對象與它們相鄰的較亮的鄰接對象相融合，這是因為受遮擋的對象沒有前景标記。
Another useful visualization technique is todisplay the label matrix as a color image. Label matrices, such asthose produced by watershed and bwlabel,can be converted to truecolor images for visualization purposes byusing label2rgb.
另外一個有用的可視化技術是将标記矩陣作為彩色圖像進行顯示。标記矩陣，比如通過watershed和bwlabel得到的，可以使用label2rgb轉換到真彩圖像來顯示。
Lrgb = label2rgb(L, 'jet', 'w','shuffle');
figure('units', 'normalized', 'position', [0 0 11]);
subplot(1, 2, 1); imshow(rgb, []);title('原圖像');
subplot(1, 2, 2); imshow(Lrgb);title('彩色分水嶺标記矩陣');
You can use transparency to superimpose thispseudo-color label matrix on top of the original intensityimage.
可以使用透明度來疊加這個僞彩色标記矩陣在原亮度圖像上進行顯示。
figure('units', 'normalized', 'position', [0 0 11]);
subplot(1, 2, 1); imshow(rgb, []);title('原圖像');
subplot(1, 2, 2); imshow(rgb, []); hold on;
himage = imshow(Lrgb);
set(himage, 'AlphaData', 0.3);
title('标記矩陣疊加到原圖像');
參考：http://blueben-chong.blog.sohu.com/141881444.html
總結
代碼：
clc; clear all; close all;
rgb = imread('pears.png');
if ndims(rgb) == 3
I =rgb2gray(rgb);
else
I = rgb;
end
hy = fspecial('sobel');
hx = hy';
Iy = imfilter(double(I), hy, 'replicate');
Ix = imfilter(double(I), hx, 'replicate');
gradmag = sqrt(Ix.^2 + Iy.^2);
L = watershed(gradmag);
Lrgb = label2rgb(L);
se = strel('disk', 20);
Io = imopen(I, se);
Ie = imerode(I, se);
Iobr = imreconstruct(Ie, I);
Ioc = imclose(Io, se);
Iobrd = imdilate(Iobr, se);
Iobrcbr = imreconstruct(imcomplement(Iobrd),imcomplement(Iobr));
Iobrcbr = imcomplement(Iobrcbr);
fgm = imregionalmax(Iobrcbr);
It1 = rgb(:, :, 1);
It2 = rgb(:, :, 2);
It3 = rgb(:, :, 3);
It1(fgm) = 255; It2(fgm) = 0; It3(fgm) = 0;
I2 = cat(3, It1, It2, It3);
se2 = strel(ones(5,5));
fgm2 = imclose(fgm, se2);
fgm3 = imerode(fgm2, se2);
fgm4 = bwareaopen(fgm3, 20);
It1 = rgb(:, :, 1);
It2 = rgb(:, :, 2);
It3 = rgb(:, :, 3);
It1(fgm4) = 255; It2(fgm4) = 0; It3(fgm4) = 0;
I3 = cat(3, It1, It2, It3);
bw = im2bw(Iobrcbr, graythresh(Iobrcbr));
D = bwdist(bw);
DL = watershed(D);
bgm = DL == 0;
gradmag2 = imimposemin(gradmag, bgm | fgm4);
L = watershed(gradmag2);
It1 = rgb(:, :, 1);
It2 = rgb(:, :, 2);
It3 = rgb(:, :, 3);
fgm5 = imdilate(L == 0, ones(3, 3)) | bgm |fgm4;
It1(fgm5) = 255; It2(fgm5) = 0; It3(fgm5) = 0;
I4 = cat(3, It1, It2, It3);
Lrgb = label2rgb(L, 'jet', 'w', 'shuffle');