分水嶺算法(分水嶺算法分割圖像)

使用C++、opencv進行分水嶺分割圖像

分水嶺概念是以對圖像進行三維可視化處理為基礎的：其中兩個是坐標，另一個是灰度級。基于“地形學”的這種解釋，我們考慮三類點：

a.屬于局部性最小值的點，也可能存在一個最小值面，該平面內的都是最小值點

b.當一滴水放在某點的位置上的時候，水一定會下落到一個單一的最小值點

c.當水處在某個點的位置上時，水會等概率地流向不止一個這樣的最小值點

對一個特定的區域最小值，滿足條件(b)的點的集合稱為這個最小值的“匯水盆地”或“分水嶺”。滿足條件(c)的點的集合組成地形表面的峰線，稱做“分割線”或“分水線”。

分水嶺分割方法，是一種基于拓撲理論的數學形態學的分割方法，目前較著名且使用較多的有2種算法:

(1) 自下而上的模擬泛洪的算法 (2) 自上而下的模擬降水的算法

這里介紹泛洪算法的過程。

我們把圖像看作是測地學上的拓撲地貌，圖像中每一點像素的灰度值表示該點的海拔高度，模擬泛洪算法的基本思想是：假設在每個區域最小值的位置上打一個洞并且讓水以均勻的上升速率從洞中涌出，從低到高淹沒整個地形。當處在不同的匯聚盆地中的水將要聚合在一起時，修建的大壩將阻止聚合。水將達到在水線上只能見到各個水壩的頂部這樣一個程度。這些大壩的邊界對應于分水嶺的分割線。所以，它們是由分水嶺算法提取出來的(連續的)邊界線。

原圖像： 地形俯視圖：

?

原圖像顯示了一個簡單的灰度級圖像，其中“山峰”的高度與輸入圖像的灰度級值成比例。為了阻止上升的水從這些結構的邊緣溢出，我們想像將整幅地形圖的周圍用比最高山峰還高的大壩包圍起來。最高山峰的值是由輸入圖像灰度級具有的最大值決定的。

?

?

?

圖一被水淹沒的第一個階段，這里水用淺灰色表示，覆蓋了對應于圖中深色背景的區域。在圖二和三中，我們看到水分別在第一和第二匯水盆地中上升。由于水持續上升，最終水將從一個匯水盆地中溢出到另一個之中。

?

?

左圖中顯示了溢出的第一個征兆。這里，水確實從左邊的盆地溢出到右邊的盆地，并且兩者之間有一個短“壩”(由單像素構成)阻止這一水位的水聚合在一起。隨著水位不斷上升，如右圖所顯示的那樣。這幅圖中在兩個匯水盆地之間顯示了一條更長的壩，另一條水壩在右上角。這條水壩阻止了盆地中的水和對應于背景的水的聚合。

這個過程不斷延續直到到達水位的最大值(對應于圖像中灰度級的最大值)。水壩最后剩下的部分對應于分水線，這條線就是要得到的分割結果。

?

對于這個例子，分水線在圖中顯示為疊加到原圖上的一個像素寬的深色路徑。注意一條重要的性質就是分水線組成一條連通的路徑，由此給出了區域之間的連續的邊界。

動圖演示了整個分水嶺算法的過程：

?

?

?

分水嶺算法對噪聲等影響非常敏感。所以在真實圖像中，由于噪聲點或者其它干擾因素的存在，使用分水嶺算法常常存在過度分割的現象，這是因為很多很小的局部極值點的存在，比如下面的圖像，這樣的分割效果是毫無用處的。

?

?

為了解決過度分割的問題，可以使用基于標記(mark)圖像的分水嶺算法，就是通過先驗知識，來指導分水嶺算法，以便獲得更好的圖像分段效果。通常的mark圖像，都是在某個區域定義了一些灰度層級，在這個區域的洪水淹沒過程中，水平面都是從定義的高度開始的，這樣可以避免一些很小的噪聲極值區域的分割。下面的動圖很好的演示了基于mark的分水嶺算法過程：

?

上面的過度分割圖像，我們通過指定mark區域，可以得到很好的分段效果：

?

?

以上參考：岡薩雷斯《數字圖象處理（第三版）》和https://www.cnblogs.com/mikewolf2002/p/3304118.html

相關API：

void tMoucallback(const string& winname, MouCallback onMou, void* urdata=0)

winname:窗口的名字onMou:鼠標響應函數，回調函數。指定窗口里每次鼠標時間發生的時候，被調用的函數指針。 這個函數的原型應該為void on_Mou(int event, int x, int y, int flags, void* param);urdate:傳給回調函數的參數

void on_Mou(int event, int x, int y, int flags, void* param)

event： CV_EVENT_*變量之一x和y：鼠標指針在圖像坐標系的坐標（不是窗口坐標系） flags：CV_EVENT_FLAG的組合， param是用戶定義的傳遞到tMouCallback函數調用的參數。

附常用的event：CV_EVENT_MOUSEMOVE、CV_EVENT_LBUTTONDOWN 、CV_EVENT_RBUTTONDOWN、 CV_EVENT_LBUTTONUP 、 CV_EVENT_RBUTTONUP

和標志位flags有關的：CV_EVENT_FLAG_LBUTTON

C++: void watershed(InputArray image,InputoutputArray markers)

第一個參數，InputArray類型的src，輸入圖像，即源圖像，填Mat類的對象即可，且需為8位三通道的彩色圖像。

第二個參數，InputOutput Array類型的markers，函數調用后的運算結果存在這里，輸入/輸出32位單通道圖像的標記結果。即這個參數用于存放函數調后的輸出結果，需和源圖片有一樣的尺寸和類型。

代碼實現：

#include "stdafx.h"#include "opencv2/imgproc/imgproc.hpp" #include "opencv2/highgui/highgui.hpp" #include <iostream> #include <fstream> using namespace cv;using namespace std;#define WINDOW_NAME1 "【程序窗口1】" //為窗口標題定義的宏 #define WINDOW_NAME2 "【分水嶺算法效果圖】" //為窗口標題定義的宏 //描述：全局變量的聲明 Mat g_maskImage, g_srcImage;Point prevPt(-1, -1);//描述：全局函數的聲明 static void ShowHelpText();static void on_Mou(int event, int x, int y, int flags, void*);int main(){ //【0】改變console字體顏色 system("color 02"); //【1】載入原圖并顯示，初始化掩膜和灰度圖 g_srcImage = imread("D:\pic-sam\哀.JPG", 1); namedWindow(WINDOW_NAME1, WINDOW_NORMAL); imshow(WINDOW_NAME1, g_srcImage); Mat srcImage, grayImage; g_srcImage.copyTo(srcImage); cvtColor(g_srcImage, g_maskImage, COLOR_BGR2GRAY); cvtColor(g_maskImage, grayImage, COLOR_GRAY2BGR); g_maskImage = Scalar::all(0); //【2】設置鼠標回調函數 tMouCallback(WINDOW_NAME1, on_Mou, 0); //【3】輪詢按鍵，進行處理 while (1) { //獲取鍵值 int c = waitKey(0); //若按鍵鍵值為ESC時，退出 if ((char)c == 27) break; //按鍵鍵值為2時，恢復源圖 if ((char)c == '2') { g_maskImage = Scalar::all(0); srcImage.copyTo(g_srcImage); imshow("image", g_srcImage); } //若檢測到按鍵值為1或者空格，則進行處理 if ((char)c == '1' || (char)c == ' ') { //定義一些參數 int i, j, compCount = 0; vector<vector<Point> > contours; vector<Vec4i> hierarchy; //尋找輪廓 findContours(g_maskImage, contours, hierarchy, RETR_CCOMP, CHAIN_APPROX_SIMPLE); //輪廓為空時的處理 if (contours.empty()) continue; //拷貝掩膜 Mat maskImage(g_maskImage.size(), CV_32S); maskImage = Scalar::all(0); //循環繪制出輪廓 for (int index = 0; index >= 0; index = hierarchy[index][0], compCount++) drawContours(maskImage, contours, index, Scalar::all(compCount + 1), -1, 8, hierarchy, INT_MAX); //compCount為零時的處理 if (compCount == 0) continue; //生成隨機顏色 /*vector<Vec3b> colorTab; for (i = 0; i < compCount; i++) { int b = theRNG().uniform(0, 255); int g = theRNG().uniform(0, 255); int r = theRNG().uniform(0, 255); colorTab.push_back(Vec3b((uchar)b, (uchar)g, (uchar)r)); }*/ //計算處理時間并輸出到窗口中 double dTime = (double)getTickCount(); watershed(srcImage, maskImage); dTime = (double)getTickCount() - dTime; printf(" 處理時間 = %gms ", dTime*1000. / getTickFrequency()); //雙層循環，將分水嶺圖像遍歷存入watershedImage中 Mat watershedImage(maskImage.size(), CV_8UC3); int index1 = 0; for (i = 0; i < maskImage.rows; i++) for (j = 0; j < maskImage.cols; j++) { if(maskImage.at<int>(i, j)>index1) index1 = maskImage.at<int>(i, j); } for (i = 0; i < maskImage.rows; i++) for (j = 0; j < maskImage.cols; j++) { int index = maskImage.at<int>(i, j); //對watershed函數生成的index的規律不是很清楚，經測試，并不是按照標記順序給出index的 //具體每一塊的index是怎么給出的還需要研究源碼 if (index == -1) watershedImage.at<Vec3b>(i, j) = Vec3b(255, 255, 255); el if (index <= 0 || index > compCount) watershedImage.at<Vec3b>(i, j) = Vec3b(0, 0, 0); el if (index ==index1) watershedImage.at<Vec3b>(i, j) = Vec3b(255, 255, 255); el watershedImage.at<Vec3b>(i, j) = Vec3b(index*10, 0, 0);//這里想給不同的物體標記為不同程度的顏色 //方便后面去除背景，顯示目標物體 } //混合灰度圖和分水嶺效果圖并顯示最終的窗口 //watershedImage = watershedImage*0.5 + grayImage*0.5; imshow(WINDOW_NAME2, watershedImage);//直接顯示分水嶺的效果圖 //這里想直接根據index，將背景顯示為黑色，需要分割出來的目標物體直接顯示 //但對index生成的規律還未搞清楚，結果可能不是很穩定 Mat src = imread("D:\pic-sam\哀.JPG", 1); for (int i = 0; i < src.rows; i++) for (int j = 0; j < src.cols; j++) { int a = abs(watershedImage.at<Vec3b>(i, j)[0] - 250) / 150; src.at<Vec3b>(i, j)[0] *= a; src.at<Vec3b>(i, j)[1] *= a; src.at<Vec3b>(i, j)[2] *= a; } namedWindow("dst", WINDOW_NORMAL); imshow("dst", src); } } return 0;}//鼠標消息回調函數 static void on_Mou(int event, int x, int y, int flags, void*){ //處理鼠標不在窗口中的情況 if (x < 0 || x >= g_srcImage.cols || y < 0 || y >= g_srcImage.rows) return; //處理鼠標左鍵相關消息 if (event == CV_EVENT_LBUTTONUP || !(flags & CV_EVENT_FLAG_LBUTTON)) prevPt = Point(-1, -1); el if (event == CV_EVENT_LBUTTONDOWN) prevPt = Point(x, y); //鼠標左鍵按下并移動，繪制出線條 el if (event == CV_EVENT_MOUSEMOVE && (flags & CV_EVENT_FLAG_LBUTTON)) { Point pt(x, y); if (prevPt.x < 0) prevPt = pt; line(g_maskImage, prevPt, pt, Scalar::all(255), 4, 8, 0); line(g_srcImage, prevPt, pt, Scalar::all(255), 4, 8, 0); prevPt = pt; imshow(WINDOW_NAME1, g_srcImage); }}// 描述：輸出一些幫助信息 static void ShowHelpText(){ printf(" 當前使用的OpenCV版本為：" CV_VERSION); printf(" ---------------------------------------------------------------------------- "); //輸出一些幫助信息 printf(" 歡迎來到【分水嶺算法】示例程序~ "); printf(" 請先用鼠標在圖片窗口中標記出大致的區域， 然后再按鍵【1】或者【SPACE】啟動算法。" " 按鍵操作說明: " " 鍵盤按鍵【1】或者【SPACE】- 運行的分水嶺分割算法 " " 鍵盤按鍵【2】- 恢復原始圖片 " " 鍵盤按鍵【ESC】- 退出程序 ");}

源圖像：

?

進行標記的圖像：

?

分水嶺算法得到的圖像：

?

分割后圖像：

?

代碼的第108-122行是對opencv分水嶺算法生成的結果圖進行分析，目前對watershed函數生成的index的規律不是很清楚，經測試，并不是按照標記順序給出index的，具體每一塊的index是怎么給出的還需要研究源碼

代碼第130-138行，目的是想直接根據分水嶺算法生成的圖像中的index，將背景顯示為黑色，需要分割出來的目標物體直接顯示，但對index生成的規律還未搞清楚，結果可能不是很穩定

以上部分參考： 毛星云 《OpenCV3編程入門》

-----------------------------------------------------

2019年4月19日增加：

查閱到opencv分水嶺算法中，在“循環繪制出輪廓”時用到一個參數compCount，這個參數并不是記錄輪廓數目的，它的作用是把每個輪廓設為同一像素值，而maskImage中的像素值就是用1-compcount 的像素值標注的，這樣問題又轉化為不清楚在查找輪廓時，算法是按照什么樣的順序找出輪廓放入vector中的。