引言
圖像分割中，采用腐蝕、膨脹等形態學操作一般可以對圖像進行一定效果的分割。腐蝕采用的kernel往往決定了分割的效果，核太小則難以分割開來，核太大則會將某些部分腐蝕掉同時對原本圖像形狀產生較大的影響。所以，換了一個思路想，從外邊分割不如從里面擴張。所以分割問題的重點放在了圖像質心的查找。
圖像處理中，針對二值化分割圖可以使用距離變幻得到一張灰度圖，而亮度中心則是質心，查找質心既可以實現簡單的分割。Matlab中查找局部峰值的函數有imregionalmax，而在python中模塊skimage中有peak_local_max函數同樣實現了對灰度圖的峰值查找。但是，在C/C++中，筆者還沒有找到相關的圖像處理庫有這一方面的功能，所以只能對着論文寫了。
1、距離變換
距離變換在圖像處理應用十分廣泛，特別在圖像分割方面。簡單來說，距離變換是針對二值化圖像，其實就是計算非零點到最近零點的距離。所以，圖像經過距離變換后，分割圖像中心則呈現出峰值。距離的定義一般有以下3種類型：
OpenCV中，則有7中類型的距離定義，除此之外用於還可以使用自己定義的距離。（所以，有無限制種吧）
/** Distance types for Distance Transform and M-estimators */ enum { CV_DIST_USER =-1, /**< User defined distance */ CV_DIST_L1 =1, /**< distance = |x1-x2| + |y1-y2| */ CV_DIST_L2 =2, /**< the simple euclidean distance */ CV_DIST_C =3, /**< distance = max(|x1-x2|,|y1-y2|) */ CV_DIST_L12 =4, /**< L1-L2 metric: distance = 2(sqrt(1+x*x/2) - 1)) */ CV_DIST_FAIR =5, /**< distance = c^2(|x|/c-log(1+|x|/c)), c = 1.3998 */ CV_DIST_WELSCH =6, /**< distance = c^2/2(1-exp(-(x/c)^2)), c = 2.9846 */ CV_DIST_HUBER =7 /**< distance = |x|<c ? x^2/2 : c(|x|-c/2), c=1.345 */ };
以下是針對兩個圓形粘連情況下采用CV_DIST_L2（常用歐式距離）作為距離進行距離變換：
2、形態學灰度重建（Morphological Grayscale Reconstruction）
一般對於局部峰值極值點的查找很容易聯想到使用窗口判斷的方法。但是，窗口的大小決定着算法的好壞。更致命的是，窗口這種方法很難應對出現連續峰值的問題（峰值平坦？不太會描述）。所以這種方法很難解決峰值查找的問題。於是，筆者在各種論壇中查找各種峰值查找的方法，最終還是回到Matlab參考文獻。
Matlab中的imregionalmax函數是使用形態學灰度重建的方法實現了局部極大值的查找，具體參考論文：Morphological Grayscale Reconstruction in Image Analysis: Applications and Effcient Algorithms。這里只進行簡單描述。
從上圖可知，對於圖像I，用其偏置I-1進行grayscale reconstruction灰度重建得到重建后的圖像。然后在使用原圖與重建后的圖像相減即可得到局部極大值部分。上圖也可以解釋Matlab中imregionalmax等於原圖與其imreconstruct相減。接下來，則是如何使用算法實現這個過程呢？論文末尾處已經給出了許多算法實現的偽代碼，其中包括並行算法、串行算法等等。本文使用隊列對算法進行實現，該算法在論文中已經給出偽代碼：
特別說明，筆者取J=I-1，R(J)是遍歷每個點取八鄰近最大值來初始化隊列。NG(p)也是取得p點的八鄰近。部分實現代碼如下：
Mat findpeeks(Mat I,Mat J) { queue<Point> que; /*最大值濾波*/ J = regionalmaximaImg(J); /*對滿足條件的點入隊列，初始化隊列*/ for(int i=0;i<I.rows;i++) for(int j=0;j<I.cols;j++) if(J.at<uchar>(i,j) && isexistNBzero(J,Point(j,i))) que.push(Point(j,i)); /*循環隊列，直至隊列為空*/ while(que.empty()) { Point p = que.front(); que.pop(); for(int i=0;i<NEGHTBOR;i++) { /*get the neighbor point*/ Point q; q.x = p.x+x_offset[i]; q.y = p.y+y_offset[i]; if(q.x>I.cols-1) q.x = I.cols-1; if(q.x<0) q.x = 0; if(q.y>I.rows-1) q.y = I.rows-1; if(q.y<0) q.y = 0; if(J.at<uchar>(p.y,p.x)<J.at<uchar>(q.y,q.x) && I.at<uchar>(q.y,q.x)!=I.at<uchar>(q.y,q.x)) { J.at<uchar>(p.y,p.x)=min(J.at<uchar>(p.y,p.x), I.at<uchar>(q.y,q.x)); que.push(p); } } } return (I-J)>0; }
實現效果如下：
由圖像可見峰值點可以找出來，但是不一定為質心，因為連接處也可能存在峰值，這里需要后續的優化。需要注意的是峰值點不是單點，可能出現連續或間斷點；筆者認為后續使用區域生長、主動輪廓模型可以完成實現分割。
結語
形態學灰度重建提供了一種峰值查找的解決方法，論文比較長，我也看得不太懂，由於前階段畢業設計急着需要用就直接按照后面的偽代碼實現了，后續有時間再將畢設后續的內容總結上來。
另外：該算法的並行算法實現詳見我的簡書：OpenCV CUDA 圖像峰值查找，我總覺得搬過來不太好。
本文源碼下載：01_findpeeks
參考文獻：
Morphological Grayscale Reconstruction in Image Analysis: Applications and Effcient Algorithms