分水嶺演算法簡單實現

一直沒看分水嶺，vc影像處理的書對該演算法的介紹也很少，cximage也找到相關的演算法。就看了matlab何opencv。matlab對分水嶺的演算法是封裝的也就沒看源碼，而opencv的寫的實在太高升了，水平差在加沒有過opencv，沒辦法看懂。網上找了幾個演算法的源碼，一個老外用template寫的，一個是國人寫的。老外寫的是NB，看不動，國人寫的看了一些也是雲裡霧裡的，最後麼辦法自己想辦法寫了個。

 

首先聲明，這個只是簡單的實現，僅僅針對一些簡單的情況。效果還可以但有不足，一些地方的分割有可能將原本的單一前景給分割掉。和我一樣的菜鳥，可以看看。高人路過可以指導。

 

先介紹下什麼叫分水嶺，

定義：

所謂映像分割是指將映像中具有特殊涵義的不同地區區分開來，這些地區是互相不交叉的，每一個地區都滿足特定地區的一致性

本代碼，沒有完全的實現，只能滿足現行對工作的需要。

分水嶺的原理網上有很多，我僅僅說明下我對這個演算法的簡單理解：該演算法可以解決針對，你可以用分割成前景，但有些幹擾或者其他原因，分割出來的前景可以連成一片，無法區分那個前景是那個。這個時候，就可以用這個分水嶺演算法來幫忙。

（先看幾張分水嶺原理的圖片，這個打字太麻煩，直接看圖片吧）

 

演算法的理解我花了半天多的時間才看明白。我的理解（可能有錯誤或者不足）就是M為灰階值最小（或者最大相對背景）為最低的盆地，其實也就是各個前景最為遠離背景的地方（這個肯定有不足的理解，希望各位看客指導一下）。根據這個盆地我們不斷的加水，使得盆地的水位不斷的提升（也就是不斷的靠近背景），在個提升地盤的時候，可能會出現和其他盆地相衝突的地方。這個時候我們需要本著和諧的態度來防止這種事情的發生，就是在衝突的地方建設堤壩。這個堤壩其實就是我們需要尋找的東西。

提升盆地水位的時候有三種可能

1：盆地上方實力，也就說這個地區的灰階和該盆地相鄰，和其他勢力沒關係。這種情況，肯定是被吞併的。

2：盆地相鄰實力，也就說這個地區的灰階不僅僅和該盆地相鄰，也和其他的盆地有關係。這種情況，就是我們需要尋找的堤壩。

3：當水位提到一定高度的時候出現的新盆地，該地區和原始的盆地均不關。該水位點為該地區最低水位。這個時候就有新的盆地產生，也就是新的惡勢力出現了。

針對上面的理解，我寫了個依據映像距離變化的分水嶺，適合一些小部分地區相交叉的情況。

首先我們需要對映像進行二值化，分割出感興趣的位置地區。二值化，大家隨便，只要能分割進行。

第二進行距離變換，代碼來自一個vc映像演算法原理和實現，大概這個書名吧

//計算二值化映像的距離值<br />bool processing::distance(BYTE *&pdata)<br />{</p><p>// 儲存象素值的數組<br />if (pdata!=NULL)<br />{<br />return false;<br />}<br />int * pnBinary = new int[m_size];<br />int * pnStore = new int[m_size];<br /> pdata = new BYTE[m_size];<br />// 將圖象二值化<br />for (int j = 0; j < m_height; ++j)<br />{<br />for(int i = 0; i < m_biwidth; ++i)<br />{<br />// 白色象素為背景，存成0<br />if(m_pbinary[i+j*m_biwidth])<br />{<br />pnBinary[i+j*m_biwidth] = 0;<br />pnStore [i+j*m_biwidth] = 0;<br />}<br />// 黑象素存成1<br />else<br />{<br />pnBinary[i+j*m_biwidth] = 1;<br />pnStore [i+j*m_biwidth] = 1;<br />}<br />}<br />}</p><p>int s = 1;</p><p>while(s == 1)<br />{<br />s = 0;</p><p>// 進行距離的累加<br />for (int j = 1; j < m_height - 1; ++j)<br />{<br />for(int i = 1; i < m_biwidth - 1; ++i)<br />{</p><p>// 如果是背景，進行下一個迴圈<br />if(pnBinary[i+j*m_biwidth]== 0)<br />continue;</p><p>// 如果當前象素值和四周的值一樣，象素值增加一<br />if(pnBinary[i+j*m_biwidth]==pnBinary[m_biwidth* (j-1) + i] && pnBinary[i+j*m_biwidth]==pnBinary[m_biwidth * (j+1) + i])<br />if(pnBinary[i+j*m_biwidth]==pnBinary[m_biwidth * j + i-1] && pnBinary[i+j*m_biwidth]==pnBinary[m_biwidth* j + i+1])<br />{<br />pnStore[m_biwidth * j + i]++;<br />s=1;<br />}<br />}<br />}</p><p>// 在進行下一輪迴圈前將當前的結果儲存<br />for (int j = 0; j < m_height; j++)<br />for(int i = 1; i < m_biwidth; i++)<br />pnBinary[m_biwidth * j + i] = pnStore[m_biwidth * j + i];<br />}</p><p>for (int j = 0; j < m_height; ++j)<br />{<br />for(int i = 0; i < m_biwidth; ++i)<br />{<br />//pdata[i+j*m_biwidth] =max(0,((25 - pnStore[j * m_biwidth + i]) * 10 + 5)); //<br />pdata[i+j*m_biwidth]=max(0,255 - min(255,pnBinary[i+j*m_biwidth]));<br />}<br />}</p><p>delete pnStore;<br />delete pnBinary;</p><p>return true;<br />}

這個函數是寫在一個類裡面的所以，一些映像的基本參數為類的成員。對這個函數我就不再說明了。

//分水嶺處理二值化距離映像<br />bool processing::Watershed(BYTE* Data,int * hist)<br />{<br />//修改成try機制<br />if (Data==NULL||hist==NULL)<br />{<br />return false;<br />}<br />queue<POINT > *statistic = new queue<POINT>[255];//記錄同級盆地隊列，背景點不處理<br />int * wdata = new int[m_size];<br />memset(wdata,255,m_size*sizeof(int));<br />bool cue = true;//變化標記<br />int start_queue=0;//標記起始隊列<br />int lab_cout=256;//盆地標號</p><p>POINT pt;</p><p>const int WSHED = 255;//水壩<br />//所有資料壓入隊列//邊間點不進行處理<br />for (int i=1;i<m_biwidth-1;++i)<br />{<br />for (int j=1;j<m_height-1;++j)<br />{<br />wdata[i+j*m_biwidth]=Data[i+j*m_biwidth];<br />if (Data[i+j*m_biwidth]<255)<br />{<br />pt.x=i;<br />pt.y=j;<br />statistic[Data[i+j*m_biwidth]].push(pt);<br />}<br />}<br />}<br />//完成尋找第一低位盆地起始水位 代碼說明1<br />for(;start_queue<256&&hist[start_queue]==0;++start_queue);</p><p>//完成尋找所有起始水位盆地 代碼說明2<br />while (hist[start_queue]>0)<br />{<br />pt=statistic[start_queue].front();<br />statistic[start_queue].pop();<br />wdata[pt.x+pt.y*m_biwidth] = lab_cout;<br />--hist[start_queue];<br />cue = true;<br />while(cue)<br />{<br />cue = false;<br />int size = (int)statistic[start_queue].size();<br />for (int i=0;i<size;++i)<br />{<br />int t;<br />pt=statistic[start_queue].front();<br />statistic[start_queue].pop();<br />//8個方向尋找lab<br />t = wdata[pt.x+(pt.y-1)*m_biwidth];<br />if (t>255)<br />{<br />wdata[pt.x+pt.y*m_biwidth]=lab_cout;<br />--hist[start_queue];<br />cue = true;<br />continue;<br />}<br />t = wdata[pt.x-1+(pt.y-1)*m_biwidth];<br />if (t>255)<br />{<br />wdata[pt.x+pt.y*m_biwidth]=lab_cout;<br />--hist[start_queue];<br />cue = true;<br />continue;<br />}<br />t = wdata[pt.x-1+pt.y*m_biwidth];<br />if (t>255)<br />{<br />wdata[pt.x+pt.y*m_biwidth]=lab_cout;<br />--hist[start_queue];<br />cue = true;<br />continue;<br />}<br />t = wdata[pt.x-1+(pt.y+1)*m_biwidth];<br />if (t>255)<br />{<br />wdata[pt.x+pt.y*m_biwidth]=lab_cout;<br />--hist[start_queue];<br />cue = true;<br />continue;<br />}<br />t = wdata[pt.x+(pt.y+1)*m_biwidth];<br />if (t>255)<br />{<br />wdata[pt.x+pt.y*m_biwidth]=lab_cout;<br />--hist[start_queue];<br />cue = true;<br />continue;<br />}<br />t = wdata[pt.x+1+(pt.y+1)*m_biwidth];<br />if (t>255)<br />{<br />wdata[pt.x+pt.y*m_biwidth]=lab_cout;<br />--hist[start_queue];<br />cue = true;<br />continue;<br />}<br />t = wdata[pt.x+1+(pt.y)*m_biwidth];<br />if (t>255)<br />{<br />wdata[pt.x+pt.y*m_biwidth]=lab_cout;<br />--hist[start_queue];<br />cue = true;<br />continue;<br />}<br />t = wdata[pt.x+1+(pt.y-1)*m_biwidth];<br />if (t>255)<br />{<br />wdata[pt.x+pt.y*m_biwidth]=lab_cout;<br />--hist[start_queue];<br />cue = true;<br />continue;<br />}<br />statistic[start_queue].push(pt);//沒有處理到繼續壓入隊列<br />}</p><p>}<br />++lab_cout;<br />}<br />//完成尋找盆地和漲水過程//背景不處理 </p><p>for (int k=start_queue+1;k<255;++k)<br />{</p><p>while (hist[k]>0)<br />{<br />cue = true;<br />while(cue)<br />{<br />cue = false;<br />int size = (int)statistic[k].size();<br />for (int i=0;i<size;++i)<br />{<br />int t;<br />int lab=0;<br />pt=statistic[k].front();<br />statistic[k].pop();<br />//8個方向尋找lab<br />t = wdata[pt.x+(pt.y-1)*m_biwidth];<br />if (t>255)<br />{<br />if (lab == 0) lab = t;<br />else if( t != lab ) lab = WSHED;</p><p>}<br />t = wdata[pt.x-1+(pt.y-1)*m_biwidth];<br />if (t>255)<br />{<br />if (lab == 0) lab = t;<br />else if( t != lab ) lab = WSHED;</p><p>}<br />t = wdata[pt.x-1+pt.y*m_biwidth];<br />if (t>255)<br />{<br />if (lab == 0) lab = t;<br />else if( t != lab ) lab = WSHED;</p><p>}<br />t = wdata[pt.x-1+(pt.y+1)*m_biwidth];<br />if (t>255)<br />{<br />if (lab == 0) lab = t;<br />else if( t != lab ) lab = WSHED;</p><p>}<br />t = wdata[pt.x+(pt.y+1)*m_biwidth];<br />if (t>255)<br />{<br />if (lab == 0) lab = t;<br />else if( t != lab ) lab = WSHED;</p><p>}<br />t = wdata[pt.x+1+(pt.y+1)*m_biwidth];<br />if (t>255)<br />{<br />if (lab == 0) lab = t;<br />else if( t != lab ) lab = WSHED;</p><p>}<br />t = wdata[pt.x+1+(pt.y)*m_biwidth];<br />if (t>255)<br />{<br />if (lab == 0) lab = t;<br />else if( t != lab ) lab = WSHED;</p><p>}<br />t = wdata[pt.x+1+(pt.y-1)*m_biwidth];<br />if (t>255)<br />{<br />if (lab == 0) lab = t;<br />else if( t != lab ) lab = WSHED;</p><p>}<br />if (lab == WSHED)<br />{<br />--hist[k];<br />wdata[pt.x+pt.y*m_biwidth]=WSHED;<br />continue;<br />}<br />else<br />{<br />if (lab == 0)<br />{<br />statistic[k].push(pt);//沒有處理到繼續壓入隊列<br />}<br />if (lab>255)<br />{<br />wdata[pt.x+pt.y*m_biwidth]=lab;<br />--hist[k];<br />cue = true;<br />}<br />}</p><p>}<br />}<br />if (!statistic[k].size())<br />{<br />pt=statistic[k].front();<br />statistic[k].pop();<br />wdata[pt.x+pt.y*m_biwidth] = ++lab_cout;;<br />--hist[k];<br />}</p><p>}</p><p>}</p><p>for (int i=0;i<m_biwidth;++i)<br />{<br />for (int j=0;j<m_height;++j)<br />{<br />if (wdata[i+j*m_biwidth]>255)<br />{<br />Data[i+j*m_biwidth] = 0;<br />}<br />else<br />Data[i+j*m_biwidth] = wdata[i+j*m_biwidth];<br />}<br />}<br />delete []statistic;<br />return true;<br />}

這個分水嶺有二個參數，第一是有二值得到的距離資料，第二個是這個距離資料的長條圖統計。

對長條圖的統計來決定迴圈的次數，當一個灰階被處理，如被水位吞沒或者成為堤壩，均認為被處理。

說明1：首先找到長條圖統計中最低的灰階做為起點水位

說明2：根據這個水位來尋找所以的盆地。此時的盆地就2中情況。一獨立，二和相關的盆地相串連。由於是根據長條圖統計來做的所以不知道是否連通。所以會出現如下情況，第一在8個方向能找到起始的盆地，那麼我們將這個點淹沒做為盆地，並將記錄標記設定為真。二8方向沒有盆地標記，那麼其有可能是此刻盆地的勢力範圍，但暫時沒有發展過來，或者和此刻的盆地無關。那麼當標記為真的情況均處理完畢的時候，說明低水位的第一個盆地已經找完了。剩下的是這個水位下的其他盆地。此時我們增加一個盆地標記，開始尋找第二個盆地。如何迴圈直到所有的初始水位盆地均被找出來。

...（太晚了明天再接著寫）