分水嶺演算法的詳細介紹（附c代碼）

文章目錄

• 分水嶺演算法（Watershed Algorithm）

分水嶺演算法（Watershed Algorithm）

所謂分水嶺演算法有好多種實現演算法，拓撲學，形態學，浸水類比和降水類比等方式。要搞懂就不容易了。Watershed Algorithm(分水嶺演算法)，顧名思義，就是根據分水嶺的構成來考慮映像的分割。現實中我們可以或者說可以想象有山有湖的景象，那麼那一定是水繞山，山圍水的情形。當然在需要的時候，要人工構築分水嶺，以防集水盆之間的互相穿透。而區分高山（plateaus）與水的界線，以及湖與湖之間的間隔或都是連通的關係，就是我們可愛的分水嶺（watershed）。為了得到一個相對集中的集水盆，那麼讓水漲到都接近周圍的最高的山頂就可以了，再漲就要漏水到鄰居了，而鄰居，嘿嘿，水質不同誒，會混淆自我的。那麼這樣的話，我們就可以用來擷取邊界灰階大，中間灰階小的物體地區了，它就是集水盆。 浸水法，就是先通過一個適當小的閾值得到起點，即集水盆的底；然後是向周圍淹沒也就是浸水的過程，直到得到分水嶺。當然如果我們要一直淹沒到山頂，即是一直處理到映像灰階最高片，那麼，當中就會出現築壩的情況，不同的集水盆在這裡想相遇了，我們要潔身自愛，到這裡為止，因為都碰到邊界了；那麼即使在相遇時沒有碰到最高灰階的地方，也需要人工構築分水嶺，區分不同的地區。不再上山。構築屬於自己的分水嶺在電腦圖形學中，可利用灰階表徵地貌高。映像中我們可以利用灰階高與地貌高的相似性來研究映像的灰階在空間上的變化。這是空域分析，比如還可以通過各種形式的梯度計算以得到演算法的輸入，進行浸水處理。分水嶺具有很強的邊緣檢測能力，對微弱的邊緣也有較好的效果。這與分水嶺擴張的閾值的設定有關係，閾值可以決定集水盆擴張的範圍。但自我構築的能力卻不受影響。為會麼這麼說呢？為什麼有很強的邊緣檢測能力，而又能得到相對集中的連通的集水盆？現實中很好辦，我們在往凹地加水的時候，直到它漲到這一塊緊湊的山嶺邊緣就不加了；但是如果有一條小山溝存在，那沒辦法，在初始閾值分割的時候，也就是山溝與集水盆有同樣的極小值，而且它們之間是以這個高度一直串連的。那沒關係，我們將它連通。在映像上呢？如何?？看看演算法，演算法思想是這樣的：首先準備好山和初始的水。這山就是我們的初始映像了，比如用自然擷取的映像的梯度來表徵山地的每一點的高度吧；而初始的水就是在閾值記為Thre底下，所有的低於這個高度的整個山地都加水，直到這個閾值Thre高度。從而有三個初始量：unsigned char** Ori_image、char** Seed_image和int** Label_image。最後一個是為最終的結果做準備的。當然要做好初始化，比如，Ora_image賦值為原映像（256色灰階圖）的梯度值，Seed_image則是初始狀態下有水的置位，無水的複位，而Label_image則全初始化為0，最終得到的是各點對應的地區號 接下來是考慮將已加的水進行記錄，記錄成連通的地區，也就是看看有多少個互不相關的集水盆，有五個，那麼我們就漲出五個湖，而且儘可能的高，只要大家想到不溢出。在演算法上，有多少個連通的地區就記錄成多少個資料結構，工夫就在於如何將這些連通的地區串連成一塊，並由一個資料結構來表達了。很好，我們準備用一個向量容器來實現初始儲存，儲存所有標記地區種子隊列的數組，裡面放的是種子隊列的指標vector<queue<POINT>*> vque; ，而且這個隊列是由一系列屬於同一個地區的映像點組成，我們來自一個集水盆：）；其儲存方式是這樣的：queue<point> *pque=new queue<point>[256]；vque.push_back(pque)，這樣便將一個成員放進到這個地區來了，即容器－－集水盆的保管都，容器中的每個指標，都指向一個集水盆，也就是我們要的連通地區；所以我們可以方便地由這個容器資料結構直接讀值的方便性進行操作，一個腳標就可以得到一個地區（隊列指標）的指標；而每個隊列還不簡單，並不是一列整形數那麼易搞，所以說啊，這個演算法，真頭痛，這個隊列的一個成員是一個點；而注意到vque裡存放的一256個隊列的的起始指標，真夠殘忍的。也就是說vque [m] [n]就表達了一個隊列，這個隊列裡可以儲存操作一系列的點；顯然容量取256是因為所有的初始或者是最終的地區中可能有0-256之間的不同的灰階的點，那麼我一個地區分用256個隊列來記錄這些成員點啦，很有可能，這裡就只有一個集水盆，那麼，256個灰階的點都存在一個地區就有可能了統計初始連通地區的方法是，八連通鄰域法（還有其他方法：四連通域，六連通域），即從逐一掃描輸入的Seed_image的每個像素點，將所有的標記了的初始集水盆一一納入各自的地區，這是整修映像的掃描，形成外迴圈。先建立一個暫存佇列quetem，用來處理當前初始集水盆的連通串連，將逐一掃描到的屬於一個特定的初始集水盆地區的可生長點暫存，並形成一個內迴圈。對當前掃描點的處理是，首先判斷該點是否為某個初始集水盆的點，如果不是跳過；接下來是，如果是初始集水盆的點，那麼它的八連通域中是否存在不可生長的點（這裡的不可生長是指Seed_image中沒有標記的點），掃描的八連通鄰域中的點是可生長的，即有標記的，則將之加入到暫存佇列中quetem；如果掃描到的連通鄰域中有不可生長的至少一個點存在，那麼加入到種子隊列，記目前範圍號為Num，當前掃描點為(m，n)，從而當前的灰階為Ori_image[m][n]，將當前點添加到種子隊列：vque[Num－1][Ori_image[m][n]].push(POINT(m，n))。這裡有兩個迴圈，一個是quetem，另一個是Seed_image。直到兩個迴圈完整結束，那麼就得到了各個連通初始集水盆的記錄，儲存標記是地區號Num；而我們同時得到了初始的分水嶺，那就放在了儲存地點vque，這裡面標識了它們對應的地區號，和地區裡面對應的點的灰階，即是特定地區特定灰階對應的點的集合；我們可以擷取這些分水嶺的點所在的地區號，可以得到該地區的所有的灰階的點資訊。一句話，統計連通地區的功能有兩個，一是標記初始地區，二是找分水嶺初始的地區標記好了，分嶺也找到了，那麼可以開始“水漫梁山”了。這就是淹沒過程。淹沒過程由也是由一個內嵌迴圈的迴圈來實現的：外迴圈是做水位上升（這裡迴圈次數一定要256以內），waterlevel的上升，原來是已經做過了初始的水位分割，那麼現在可以從Thre開始了，讓水位慢慢上升，讓它原本的湖慢慢擴張，盡量利用其應有的空間，而又不至於淹沒到其它的鄰居湖泊。內迴圈是掃描每個初始地區（當前Num，從而有Num個迴圈）的分水嶺的點（在vque[][]中），按照給定的水位進行擴張。擴張過程是這樣的：掃描到的分水嶺的當前點，對其進行四連通鄰域進行逐一檢查，如果四連通域中有點沒有標記的（那這一定是高度較高的點，較低的前面一定掃描過），那麼先對該點以本地區號做標記Num（注意是當前的Num）；再判斷它在當前水位下是否可生長（灰階是否小於等於waterlevel），如果可生長那麼加入到vque[Num][waterlevel]種子隊列中，將會再次進入內迴圈，否則如果在當前水位下不可生長，則加入到這個鄰域點的分水嶺集合中vque[Num][Ori_image[鄰域點]]隊列中。如此往複迴圈，直到對應地區完成，一個水位，掃描所有的地區的分水嶺，這樣各自同時在一個水位下擴張，保證了不出現跳躍的情況出現（就是一個水位一個地區全域擴張）。最終，所有的地區在每個水位都擴張完畢，得到了分割圖，我們的大湖泊形成了這是分水嶺演算法的一種實現方式。仔細考察不難發現這種實現方式不能產生新的集水盆，也就是說，由於初始集水盆的局限性，很可能會浪費大部分沒有發掘出來的起始點較高的集水盆地。這樣，我們就要對演算法進行修改了。實現方式是：在淹沒的過程中，我們是由閾值Thre的水位開始淹沒的，那麼我們可以對初始地區之外的沒有標記的點（從Seed_image中考察），對之進行標記，條件是先把這一輪的內迴圈做好，然後在剩下的沒標記地區中發掘新的集水盆，並加入到我們的種子隊列中，下一個水位開始，就又多了一個新成員了，讓它們不斷膨脹，成長，擁有自己的小天的成員就會逐一的被分割出來。不過話說回來，我們這裡是採用梯度映像，一般情況下，閾值初始分割能夠滿足我們的要求，把灰階變化平滑的先截取出來，梯度資訊已然足夠強大；而如果採用了新盆地擴張，則比較適用於原始映像。分水嶺演算法主要的分割目的在於找到映像的連通地區。利用梯度資訊作為輸入映像，會有一個矛盾點，如果對原始映像進行梯度計算時不作濾波平滑處理，很容易將物體分割成多個物體，那是因為雜訊的影響；而如果進行濾波處理，又容易造成將某些原本幾個的物體合成一個物體。當然這裡的物體主要還是指映像變化不大或者說是灰階值相近的目的地區域 void WINAPI CDib::Watershed(unsigned char **OriginalImage, char** SeedImage, int **LabelImage, int row, int col){// using namespace std;//標記地區標識號，從1開始int Num=0;int i,j;//儲存每個隊列種子個數的數組vector<int*> SeedCounts;//臨時種子隊列queue<POINT> quetem;//儲存所有標記地區種子隊列的數組，裡面放的是種子隊列的指標vector<queue<POINT>*> vque;int* array;//指向種子隊列的指標queue<POINT> *pque;POINT temp;for(i=0;i<row;i++){ for(j=0;j<col;j++) LabelImage[i][j]=0;}int m,n,k=0;BOOL up,down,right,left,upleft,upright,downleft,downright;//8 directions...//預先處理,提取區分每個標記地區，並初始化每個標記的種子隊列//種子是指標記地區邊緣的點，他們可以在水位上升時向外淹沒（或者說生長）//pan's words：我的理解是梯度值較小的象素點，或者是極小灰階值的點。for(i=0;i<row;i++){ for(j=0;j<col;j++) { //如果找到一個標記地區 if(SeedImage[i][j]==1) { //地區的標識號加一 Num++; //分配數組並初始化為零，表示可有256個灰階 array=new int[256]; ZeroMemory(array,256*sizeof(int)); //種子個數數組進vector，每次掃描則產生一個數組，並用地區標識號來做第一維。灰階級做第二維。 //表示某個盆地地區中某灰階所對應的點的數目。 SeedCounts.push_back(array); //分配本標記號的優先隊列，256個種子隊列， //表示對應一個灰階有一個隊列，並且每個隊列可以儲存一個集合的點資訊 pque=new queue<POINT>[256]; //加入到隊列數組中，對應的是本標記號Num的 vque.push_back(pque); //當前點放入本標記地區的臨時種子隊列中 temp.x=i; temp.y=j; quetem.push(temp); //當前點標記為已處理 LabelImage[i][j]=Num; SeedImage[i][j]=127;//表示已經處理過 //讓臨時種子隊列中的種子進行生長直到所有的種子都生長完畢 //生長完畢後的隊列資訊儲存在vque中，包括地區號和灰階，對應點數儲存在seedcounts中 while(!quetem.empty()) { up=down=right=left=FALSE; upleft=upright=downleft=downright=FALSE; //隊列中取出一個種子 temp=quetem.front(); m=temp.x; n=temp.y; quetem.pop(); //注意到127對掃描過程的影響，影響下面的比較，但是不影響while語句中的掃描 if(m>0) { //上方若為可生長點則加為新種子 if(SeedImage[m-1][n]==1) { temp.x=m-1; temp.y=n; quetem.push(temp);//如果這樣的話，那麼這些標記過的地區將再次在while迴圈中被掃描到，不會，因為值是127 //新種子點標記為已淹沒地區，而且是目前範圍，並記錄地區號到labelImage LabelImage[m-1][n]=Num; SeedImage[m-1][n]=127; } else//否則上方為不可生長 { up=TRUE; } } if(m>0&&n>0) { if(SeedImage[m-1][n-1]==1)//左上方若為可生長點則加為新種子 { temp.x=m-1; temp.y=n-1; quetem.push(temp); //新種子點標記為已淹沒地區，即下一個迴圈中以127來標識不再掃描，而且是目前範圍 LabelImage[m-1][n-1]=Num; SeedImage[m-1][n-1]=127; } else//否則左上方為不可生長 { upleft=TRUE; } } if(m<row-1) { if(SeedImage[m+1][n]==1)//下方若為可生長點則加為新種子 { temp.x=m+1; temp.y=n; quetem.push(temp); //新種子點標記為已淹沒地區，而且是目前範圍 LabelImage[m+1][n]=Num; SeedImage[m+1][n]=127; } else//否則下方為不可生長 { down=TRUE; } } if(m<(row-1)&&n<(col-1)) { if(SeedImage[m+1][n+1]==1)//下方若為可生長點則加為新種子 { temp.x=m+1; temp.y=n+1; quetem.push(temp); //新種子點標記為已淹沒地區，而且是目前範圍 LabelImage[m+1][n+1]=Num; SeedImage[m+1][n+1]=127; } else//否則下方為不可生長 { downright=TRUE; } } if(n<col-1) { if(SeedImage[m][n+1]==1)//右方若為可生長點則加為新種子 { temp.x=m; temp.y=n+1; quetem.push(temp); //新種子點標記為已淹沒地區，而且是目前範圍 LabelImage[m][n+1]=Num; SeedImage[m][n+1]=127; } else//否則右方為不可生長 { right=TRUE; } } if(m>0&&n<(col-1)) { if(SeedImage[m-1][n+1]==1)//右上方若為可生長點則加為新種子 { temp.x=m-1; temp.y=n+1; quetem.push(temp); //新種子點標記為已淹沒地區，而且是目前範圍 LabelImage[m-1][n+1]=Num; SeedImage[m-1][n+1]=127; } else//否則右上方為不可生長 { upright=TRUE; } } if(n>0) { if(SeedImage[m][n-1]==1)//左方若為可生長點則加為新種子 { temp.x=m; temp.y=n-1; quetem.push(temp); //新種子點標記為已淹沒地區 LabelImage[m][n-1]=Num; SeedImage[m][n-1]=127; } else//否則左方為不可生長 { left=TRUE; } } if(m<(row-1)&&n>0) { if(SeedImage[m+1][n-1]==1)//左下方若為可生長點則加為新種子 { temp.x=m+1; temp.y=n-1; quetem.push(temp); //新種子點標記為已淹沒地區 LabelImage[m+1][n-1]=Num; SeedImage[m+1][n-1]=127; } else//否則左方為不可生長 { downleft=TRUE; } } //上下左右只要有一點不可生長，那麼本點為初始種子隊列中的一個 //這裡可否生長是由seedimage中的值來決定的。 if(up||down||right||left|| upleft||downleft||upright||downright) { temp.x=m; temp.y=n; //下面這個向量數組：第一維是標記號；第二維是該映像點的灰階級 //m,n點對應的是while迴圈裡面掃描的像素點。 //Num是當前的地區號 //這樣這個二維資訊就表示了，某個地區中對應某個灰階級對應的成員點的集合與個數 //分別由下面兩個量來表達 vque[Num-1][OriginalImage[m][n]].push(temp);//這兩句中我把Num-1改成了Num...pan's codes... SeedCounts[Num-1][OriginalImage[m][n]]++; } }//while結束，掃描到quetem為空白而止。也就是對應所有的節點都得到不可生長為止（或者是周圍的點要麼不可生長，要麼已生長） }//if結束 // if(Num==5) // return; } }//在上述過程中，如果標記的點為0則表示，沒有掃描到的點，或者表明不是輸入的種子點//這裡相當於是找seedimage傳過來的初始地區的分水嶺界線的所有的點；並且用標號記錄每個地區，同時集水盆的邊緣點進入隊列。//上面是找集水盆的程式。同時也是連通地區。/*************************************///test 這裡測試一下剩下的非水盆地的點數。int seednum;for(i=0;i<row;i++){ for(j=0;j<col;j++) { if(SeedImage[i][j]==0) seednum++; }}CString str;str.Format("pre region num:%d",Num);AfxMessageBox(str);/*************************************/ bool actives;//在某一水位處，所有標記的種子生長完的標誌 int WaterLevel;//淹沒過程開始，水位從零開始上升，水位對應灰階級，採用四連通法for(WaterLevel=0;WaterLevel<180;WaterLevel++)//第二維。。。{ actives=true; while(actives) { actives=false; //依次處理每個標記號所對應的地區，且這個標記號對應的地區的點的個數在SeedCounts裡面 for(i=0;i<Num;i++)//第一維。。。 { if(!vque[i][WaterLevel].empty())//對應的分水嶺不為空白集,i表示地區號，waterlevel表示灰階 { actives=true; while(SeedCounts[i][WaterLevel]>0) { SeedCounts[i][WaterLevel]--;//取出一個點，個數少一 temp=vque[i][WaterLevel].front();//取出該地區的一個點，清空這個邊緣點，表示當前 //灰階級該像素已經處理掉了。 vque[i][WaterLevel].pop(); m = temp.x; n = temp.y;//當前種子的座標 if(m>0) { if(!LabelImage[m-1][n])//上方若未處理，表示沒有標號，應該在輸入前已經作過初始化為0 //本函數中在開頭也作過初始化 { temp.x=m-1; temp.y=n; LabelImage[m-1][n]=i+1;//上方點標記為已淹沒地區 //注意到這個標記是與掃描點的地區號相同，一定在這個標號所屬的地區嗎？是的 //這樣在下一輪至少會掃描到這個點，確保不遺漏，但是下一輪的處理會使它合理 //歸類嗎？問題還有這樣標記並沒有一定將它加入到種子隊列。也就是說它 //只是被淹沒而不能向上淹沒。只有滿足下述可生長條件才行。 if(OriginalImage[m-1][n]<=WaterLevel)//上方若為可生長點則加入當前隊列，當前高度的隊列 { vque[i][WaterLevel].push(temp); } else//否則加入OriginalImage[m-1][n]對應的灰階級的隊列，為什嗎？ { vque[i][OriginalImage[m-1][n]].push(temp); SeedCounts[i][OriginalImage[m-1][n]]++; } } } if(m<row-1) { if(!LabelImage[m+1][n])//下方若未處理 { temp.x=m+1; temp.y=n; LabelImage[m+1][n]=i+1;//下方點標記為已淹沒地區 if(OriginalImage[m+1][n]<=WaterLevel)//下方若為可生長點則加入當前隊列 { vque[i][WaterLevel].push(temp); } else//否則加入OriginalImage[m+1][n]級隊列 { vque[i][OriginalImage[m+1][n]].push(temp); SeedCounts[i][OriginalImage[m+1][n]]++; } } } if(n<col-1) { if(!LabelImage[m][n+1])//右邊若未處理 { temp.x=m; temp.y=n+1; LabelImage[m][n+1]=i+1;//右邊點標記為已淹沒地區 if(OriginalImage[m][n+1]<=WaterLevel)//右邊若為可生長點則加入當前隊列 { vque[i][WaterLevel].push(temp); } else//否則加入OriginalImage[m][n+1]級隊列 { vque[i][OriginalImage[m][n+1]].push(temp); SeedCounts[i][OriginalImage[m][n+1]]++; } } } if(n>0) { if(!LabelImage[m][n-1])//左邊若未處理 { temp.x=m; temp.y=n-1; LabelImage[m][n-1]=i+1;//左邊點標記為已淹沒地區 if(OriginalImage[m][n-1]<=WaterLevel)//左邊若為可生長點則加入當前隊列 { vque[i][WaterLevel].push(temp); } else//否則加入OriginalImage[m][n-1]級隊列 { vque[i][OriginalImage[m][n-1]].push(temp); SeedCounts[i][OriginalImage[m][n-1]]++; } } } }//while迴圈結束 SeedCounts[i][WaterLevel]=vque[i][WaterLevel].size(); }//if結束 }//for迴圈結束 }//while迴圈結束 }//for迴圈結束/**********************************//test whether the origional segmentation num is changed...int nonzeronum=Num;for(m=0;m<Num;m++){ for(i=0;i<row;i++) { for(j=0;j<col;j++) { if(LabelImage[i][j]==m) break; } if(LabelImage[i][j]==m) break; } if(j==col-1 && i==row-1 && LabelImage[i][j]!=m) nonzeronum--;}str.Format("new region num:%d",nonzeronum);//這表示的是可擴張的點AfxMessageBox(str);/*********************************/ while(!vque.empty()){ pque=vque.back(); delete[] pque; vque.pop_back();}while(!SeedCounts.empty()){ array=SeedCounts.back(); delete[] array; SeedCounts.pop_back();}}