1. 程式人生 > >分水嶺演算法的詳細介紹(附c…

分水嶺演算法的詳細介紹(附c…

分水嶺演算法(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*> vque; ,而且這個佇列是由一系列屬於同一個區域的影象點組成,我們來自一個集水盆:);其儲存方式是這樣的:queue *pque=new queue[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中考察),對之進行標記,條件是先把這一輪的內迴圈做好,然後在剩下的沒標記區域中發掘新的集水盆,並加入到我們的種子隊列中,下一個水位開始,就又多了一個新成員了,讓它們不斷膨脹,成長,擁有自己的小天的成員就會逐一的被分割出來。不過話說回來,我們這裡是採用梯度影象,一般情況下,閾值初始分割能夠滿足我們的要求,把灰階變化平滑的先截取出來,梯度資訊已然足夠強大;而如果採用了新盆地擴張,則比較適用於原始影象。分水嶺演算法主要的分割目的在於找到影象的連通區域。利用梯度資訊作為輸入影象,會有一個矛盾點,如果對原始影象進行梯度計算時不作濾波平滑處理,很容易將物體分割成多個物體,那是因為噪聲的影響;而如果進行濾波處理,又容易造成將某些原本幾個的物體合成一個物體。當然這裡的物體主要還是指影象變化不大或者說是灰度值相近的目標區域

[cpp] view plaincopyprint?
  1. void WINAPI CDib::Watershed(unsigned char **OriginalImage, char** SeedImage, int **LabelImage, int row, int col)
  2. {
  3. // using namespace std;
  4. //標記區域標識號,從1開始
  5. int Num=0;
  6. int i,j;
  7. //儲存每個佇列種子個數的陣列
  8. vector<<SPAN class=datatypes>int*> SeedCounts;
  9. //臨時種子隊列
  10. queue quetem;
  11. //儲存所有標記區域種子隊列的陣列,裡面放的是種子隊列的指標
  12. vector*> vque;
  13. int* array;
  14. //指向種子隊列的指標
  15. queue *pque;
  16. POINT temp;
  17. for(i=0;i
  18. {
  19. for(j=0;j
  20. LabelImage[i][j]=0;
  21. }
  22. int m,n,k=0;
  23. BOOL up,down,right,left,upleft,upright,downleft,downright;//8 directions...
  24. //預處理,提取區分每個標記區域,並初始化每個標記的種子隊列
  25. //種子是指標記區域邊緣的點,他們可以在水位上升時向外淹沒(或者說生長)
  26. //pan's words:我的理解是梯度值較小的象素點,或者是極小灰度值的點。
  27. for(i=0;i
  28. {
  29. for(j=0;j
  30. {
  31. //如果找到一個標記區域
  32. if(SeedImage[i][j]==1)
  33. {
  34. //區域的標識號加一
  35. Num++;
  36. //分配陣列並初始化為零,表示可有256個灰階
  37. array=new int[256];
  38. ZeroMemory(array,256*sizeof(int));
  39. //種子個數陣列進vector,每次掃描則生成一個數組,並用區域標識號來做第一維。灰度級做第二維。
  40. //表示某個盆地區域中某灰階所對應的點的數目。
  41. SeedCounts.push_back(array);
  42. //分配本標記號的優先佇列,256個種子隊列,
  43. //表示對應一個灰階有一個佇列,並且每個佇列可以儲存一個集合的點資訊
  44. pque=new queue[256];
  45. //加入到佇列陣列中,對應的是本標記號Num的
  46. vque.push_back(pque);
  47. //當前點放入本標記區域的臨時種子隊列中
  48. temp.x=i;
  49. temp.y=j;
  50. quetem.push(temp);
  51. //當前點標記為已處理
  52. LabelImage[i][j]=Num;
  53. SeedImage[i][j]=127;//表示已經處理過
  54. //讓臨時種子隊列中的種子進行生長直到所有的種子都生長完畢
  55. //生長完畢後的佇列資訊儲存在vque中,包括區域號和灰階,對應點數儲存在seedcounts中
  56. while(!quetem.empty())
  57. {
  58. up=down=right=left=FALSE;
  59. upleft=upright=downleft=downright=FALSE;
  60. //佇列中取出一個種子
  61. temp=quetem.front();
  62. m=temp.x;
  63. n=temp.y;
  64. quetem.pop();
  65. //注意到127對掃描過程的影響,影響下面的比較,但是不影響while語句中的掃描
  66. if(m>0)
  67. {
  68. //上方若為可生長點則加為新種子
  69. if(SeedImage[m-1][n]==1)
  70. {
  71. temp.x=m-1;
  72. temp.y=n;
  73. quetem.push(temp);//如果這樣的話,那麼這些標記過的區域將再次在while迴圈中被掃描到,不會,因為值是127
  74. //新種子點標記為已淹沒區域,而且是當前區域,並記錄區域號到labelImage
  75. LabelImage[m-1][n]=Num;
  76. SeedImage[m-1][n]=127;
  77. }
  78. else//否則上方為不可生長
  79. {
  80. up=TRUE;
  81. }
  82. }
  83. if(m>0&&n>0)
  84. {
  85. if(SeedImage[m-1][n-1]==1)//左上方若為可生長點則加為新種子
  86. {
  87. temp.x=m-1;
  88. temp.y=n-1;
  89. quetem.push(temp);
  90. //新種子點標記為已淹沒區域,即下一個迴圈中以127來標識不再掃描,而且是當前區域
  91. LabelImage[m-1][n-1]=Num;
  92. SeedImage[m-1][n-1]=127;
  93. }
  94. else//否則左上方為不可生長
  95. {
  96. upleft=TRUE;
  97. }
  98. }
  99. if(m
  100. {
  101. if(SeedImage[m+1][n]==1)//下方若為可生長點則加為新種子
  102. {
  103. temp.x=m+1;
  104. temp.y=n;
  105. quetem.push(temp);
  106. //新種子點標記為已淹沒區域,而且是當前區域
  107. LabelImage[m+1][n]=Num;
  108. SeedImage[m+1][n]=127;
  109. }
  110. else//否則下方為不可生長
  111. {
  112. down=TRUE;
  113. }
  114. }
  115. if(m<(row-1)&&n<(col-1))
  116. {
  117. if(SeedImage[m+1][n+1]==1)//下方若為可生長點則加為新種子
  118. {
  119. temp.x=m+1;
  120. temp.y=n+1;
  121. quetem.push(temp);
  122. //新種子點標記為已淹沒區域,而且是當前區域
  123. LabelImage[m+1][n+1]=Num;
  124. SeedImage[m+1][n+1]=127;
  125. }
  126. else//否則下方為不可生長
  127. {
  128. downright=TRUE;
  129. }
  130. }
  131. if(n
  132. {
  133. if(SeedImage[m][n+1]==1)//右方若為可生長點則加為新種子
  134. {
  135. temp.x=m;
  136. temp.y=n+1;
  137. quetem.push(temp);
  138. //新種子點標記為已淹沒區域,而且是當前區域
  139. LabelImage[m][n+1]=Num;
  140. SeedImage[m][n+1]=127;
  141. }
  142. else//否則右方為不可生長
  143. {
  144. right=TRUE;
  145. }
  146. }
  147. if(m>0&&n<(col-1))
  148. {
  149. if(SeedImage[m-1][n+1]==1)//右上方若為可生長點則加為新種子
  150. {
  151. temp.x=m-1;
  152. temp.y=n+1;
  153. quetem.push(temp);
  154. //新種子點標記為已淹沒區域,而且是當前區域
  155. LabelImage[m-1][n+1]=Num;
  156. SeedImage[m-1][n+1]=127;
  157. }
  158. else//否則右上方為不可生長
  159. {
  160. upright=TRUE;
  161. }
  162. }
  163. if(n>0)
  164. {
  165. if(SeedImage[m][n-1]==1)//左方若為可生長點則加為新種子
  166. {
  167. temp.x=m;
  168. temp.y=n-1;
  169. quetem.push(temp);
  170. //新種子點標記為已淹沒區域
  171. LabelImage[m][n-1]=Num;
  172. SeedImage[m][n-1]=127;
  173. }
  174. else//否則左方為不可生長
  175. {
  176. left=TRUE;
  177. }
  178. }
  179. if(m<(row-1)&&n>0)
  180. {
  181. if(SeedImage[m+1][n-1]==1)//左下方若為可生長點則加為新種子
  182. {
  183. temp.x=m+1;
  184. temp.y=n-1;
  185. quetem.push(temp);
  186. //新種子點標記為已淹沒區域
  187. LabelImage[m+1][n-1]=Num;
  188. SeedImage[m+1][n-1]=127;
  189. }
  190. else//否則左方為不可生長
  191. {
  192. downleft=TRUE;
  193. }
  194. }
  195. //上下左右只要有一點不可生長,那麼本點為初始種子隊列中的一個
  196. //這裡可否生長是由seedimage中的值來決定的。
  197. if(up||down||right||left||
  198. upleft||downleft||upright||downright)
  199. {
  200. temp.x=m;
  201. temp.y=n;
  202. //下面這個向量陣列:第一維是標記號;第二維是該影象點的灰度級
  203. //m,n點對應的是while迴圈裡面掃描的畫素點。
  204. //Num是當前的區域號
  205. //這樣這個二維資訊就表示了,某個區域中對應某個灰度級對應的成員點的集合與個數
  206. //分別由下面兩個量來表達
  207. vque[Num-1][OriginalImage[m][n]].push(temp);//這兩句中我把Num-1改成了Num...pan's codes...
  208. SeedCounts[Num-1][OriginalImage[m][n]]++;
  209. }
  210. }//while結束,掃描到quetem為空而止。也就是對應所有的節點都得到不可生長為止(或者是周圍的點要麼不可生長,要麼已生長)
  211. }//if結束
  212. // if(Num==5)
  213. // return;
  214. }
  215. }
  216. //在上述過程中,如果標記的點為0則表示,沒有掃描到的點,或者表明不是輸入的種子點
  217. //這裡相當於是找seedimage傳過來的初始區域的分水嶺界線的所有的點;並且用標號記錄每個區域,同時集水盆的邊緣點進入佇列。
  218. //上面是找集水盆的程式。同時也是連通區域。
  219. //test 這裡測試一下剩下的非水盆地的點數。
  220. int seednum;
  221. for(i=0;i
  222. {
  223. for(j=0;j
  224. {
  225. if(SeedImage[i][j]==0)
  226. seednum++;
  227. }
  228. }
  229. CString str;
  230. str.Format("pre region num:%d",Num);
  231. AfxMessageBox(str);
  232. bool actives;//在某一水位處,所有標記的種子生長完的標誌
  233. int WaterLevel;
  234. //淹沒過程開始,水位從零開始上升,水位對應灰度級,採用四連通法
  235. for(WaterLevel=0;WaterLevel<180;WaterLevel++)//第二維。。。
  236. {
  237. actives=true;
  238. while(actives)
  239. {
  240. actives=false;
  241. //依次處理每個標記號所對應的區域,且這個標記號對應的區域的點的個數在SeedCounts裡面
  242. for(i=0;i//第一維。。。
  243. {
  244. if(!vque[i][WaterLevel].empty())//對應的分水嶺不為空集,i表示區域號,waterlevel表示灰階
  245. {
  246. actives=true;
  247. while(SeedCounts[i][WaterLevel]>0)
  248. {
  249. SeedCounts[i][WaterLevel]--;//取出一個點,個數少一
  250. temp=vque[i][WaterLevel].front();//取出該區域的一個點,清空這個邊緣點,表示當前
  251. //灰度級該畫素已經處理掉了。
  252. vque[i][WaterLevel].pop();
  253. m = temp.x;
  254. n = temp.y;//當前種子的座標
  255. if(m>0)
  256. {
  257. if(!LabelImage[m-1][n])//上方若未處理,表示沒有標號,應該在輸入前已經作過初始化為0
  258. //本函式中在開頭也作過初始化
  259. {
  260. temp.x=m-1;
  261. temp.y=n;
  262. LabelImage[m-1][n]=i+1;//上方點標記為已淹沒區域
  263. //注意到這個標記是與掃描點的區域號相同,一定在這個標號所屬的區域嗎?是的
  264. //這樣在下一輪至少會掃描到這個點,確保不遺漏,但是下一輪的處理會使它合理
  265. //歸類嗎?問題還有這樣標記並沒有一定將它加入到種子隊列。也就是說它
  266. //只是被淹沒而不能向上淹沒。只有滿足下述可生長條件才行。
  267. if(OriginalImage[m-1][n]<=WaterLevel)//上方若為可生長點則加入當前佇列,當前高度的佇列
  268. {
  269. vque[i][WaterLevel].push(temp);
  270. }
  271. else//否則加入OriginalImage[m-1][n]對應的灰度級的佇列,為什麼?
  272. {
  273. vque[i][OriginalImage[m-1][n]].push(temp);
  274. SeedCounts[i][OriginalImage[m-1][n]]++;
  275. }
  276. }
  277. }
  278. if(m
  279. {
  280. if(!LabelImage[m+1][n])//下方若未處理
  281. {
  282. temp.x=m+1;
  283. temp.y=n;
  284. LabelImage[m+1][n]=i+1;//下方點標記為已淹沒區域
  285. if(OriginalImage[m+1][n]<=WaterLevel)//下方若為可生長點則加入當前佇列
  286. {
  287. vque[i][WaterLevel].push(temp);
  288. }
  289. else//否則加入OriginalImage[m+1][n]級佇列
  290. {
  291. vque[i][OriginalImage[m+1][n]].push(temp);
  292. SeedCounts[i][OriginalImage[m+1][n]]++;
  293. }
  294. }
  295. }
  296. if(n
  297. {
  298. if(!LabelImage[m][n+1])//右邊若未處理
  299. {
  300. temp.x=m;
  301. temp.y=n+1;
  302. LabelImage[m][n+1]=i+1;//右邊點標記為已淹沒區域
  303. if(OriginalImage[m][n+1]<=WaterLevel)//右邊若為可生長點則加入當前佇列
  304. {
  305. vque[i][WaterLevel].push(temp);
  306. }
  307. else//否則加入OriginalImage[m][n+1]級佇列
  308. {
  309. vque[i][OriginalImage[m][n+1]].push(temp);
  310. SeedCounts[i][OriginalImage[m][n+1]]++;
  311. }
  312. }
  313. }
  314. if(n>0)
  315. {
  316. if(!LabelImage[m][n-1])//左邊若未處理
  317. {
  318. temp.x=m;
  319. temp.y=n-1;
  320. LabelImage[m][n-1]=i+1;//左邊點標記為已淹沒區域
  321. if(OriginalImage[m][n-1]<=WaterLevel)//左邊若為可生長點則加入當前佇列
  322. {
  323. vque[i][WaterLevel].push(temp);
  324. }
  325. else//否則加入OriginalImage[m][n-1]級佇列
  326. {
  327. vque[i][OriginalImage[m][n-1]].push(temp);
  328. SeedCounts[i][OriginalImage[m][n-1]]++;
  329. }
  330. }
  331. }
  332. }//while迴圈結束
  333. SeedCounts[i][WaterLevel]=vque[i][WaterLevel].size();
  334. }//if結束
  335. }//for迴圈結束
  336. }//while迴圈結束
  337. }//for迴圈結束
  338. while(!vque.empty())
  339. {
  340. pque=vque.back();
  341. delete[] pque;
  342. vque.pop_back();
  343. }
  344. while(!SeedCounts.empty())
  345. {
  346. array=SeedCounts.back();
  347. delete[] array;
  348. SeedCounts.pop_back();
  349. }
  350. }