openCv學習筆記（九）-自己實現分水嶺演算法並且用openCv測試

   為了驗證上一篇博文中的演算法，自己稍作修改，即利用openCv裡的映像結構。為了簡單，直接採用Mat.at<uchar>(i,j)進行映像操作，效率肯定低，如果感興趣可參考http://blog.csdn.net/caiye917015406/article/details/7791815改進演算法。不過這裡只是測試效果，就將就吧。。。

  在演算法中要手動選取閥值進行映像的二值化，這給分水嶺演算法的效果有很大影響。也曾採用openCv裡的cvAdaptiveThreshold自動選取 尋求,但效果還不如不好，期待改進。。。

 以下是測試代碼

#include<iostream>#include <queue>#include <vector>#include<cv.h>#include<highgui.h>using namespace std;using namespace cv;void Watershed_1(const Mat OriginalImage,Mat SeedImage,Mat LabelImage);/*====================================================================函數名：            Watershed功能：              用標記-分水嶺演算法對輸入映像進行分割演算法實現：          無輸入參數說明：      OriginalImage --輸入映像（灰階圖，0~255）                    SeedImage       --標記映像（二值圖，0-非標記，1-標記）                    LabelImage      --輸出映像（1-第一個分割地區，2-第二個分割地區，...）傳回值說明：        無      ====================================================================*/int main(){//聲明變數IplImage *OriginalImage,*SeedImage,*LabelImage,*LabelImage2;OriginalImage=cvLoadImage("D:\\openCV\\openCVProject\\img\\stuff.jpg",CV_LOAD_IMAGE_GRAYSCALE);//OriginalImage_gray =cvCreateImage(cvGetSize(OriginalImage), 8,1); //cvtColor( OriginalImage, OriginalImage, CV_BGR2GRAY );//轉化為灰階圖cvNamedWindow("sorc");cvNamedWindow("seed");cvNamedWindow("result1");cvNamedWindow("result2");    //建立操作映像LabelImage = cvCreateImage(cvGetSize(OriginalImage), 8,1);SeedImage = cvCreateImage(cvGetSize(OriginalImage), 8,1);LabelImage2= cvCreateImage(cvGetSize(OriginalImage), 8,3);    //映像二值化，要手動取閥值cvThreshold(OriginalImage,SeedImage,195,1,CV_THRESH_BINARY);    //cvAdaptiveThreshold(OriginalImage,SeedImage,100);cvShowImage("sorc",OriginalImage);cvShowImage("seed",SeedImage);    Watershed_1( OriginalImage,SeedImage, LabelImage);    cvShowImage("result1",LabelImage);    //找出處理後的輪廓CvMemStorage* storage = cvCreateMemStorage();CvSeq* first_contour =NULL;int Nc = cvFindContours(LabelImage,storage,&first_contour,sizeof(CvContour),CV_RETR_CCOMP);int n=0;printf("Toatal Contours Detected: %d\n",Nc);cvZero(LabelImage2);//畫出處理後的輪廓for(CvSeq* c=first_contour;c!=NULL;c=c->h_next){//cvCvtColor(img_8uc1,img_8uc3,CV_GRAY2BGR);cvDrawContours(LabelImage2,c,CV_RGB(0,200,0), CV_RGB(255, 0, 0),2,2,8);//printf("Contour #%d\n",n);//cvShowImage("result",LabelImage);//printf("%d elements:\n",c->total);//for(int i=0;i<c->total;++i)//{//CvPoint* p = CV_GET_SEQ_ELEM(CvPoint,c,i);//printf("(%d,%d)\n",p->x,p->y );//}//cvWaitKey(0);n++;}cvShowImage("result2",LabelImage2);//cvShowImage("result2",markers);cvWaitKey(0);cvReleaseImage(&OriginalImage);cvDestroyWindow("sorc");cvReleaseImage(&SeedImage);cvDestroyWindow("seed");cvReleaseImage(&LabelImage);cvDestroyWindow("result1");cvReleaseImage(&LabelImage2);cvDestroyWindow("result2");return 0;}void Watershed_1(const Mat OriginalImage,Mat SeedImage,Mat LabelImage){int row,col;row = OriginalImage.rows;col = OriginalImage.cols;//標記地區標識號，從1開始int Num=0;int i,j;//儲存每個隊列種子個數的數組vector<int*> SeedCounts;//臨時種子隊列queue<POINT> quetem;//儲存所有標記地區種子隊列的數組，裡面放的是種子隊列的指標vector<queue<POINT>*> vque;int* array;//指向種子隊列的指標.at<uchar>(i,j)queue<POINT> *pque;POINT temp;for(i=0;i<row;i++){for(j=0;j<col;j++)LabelImage.at<uchar>(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.at<uchar>(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.at<uchar>(i,j)=Num;SeedImage.at<uchar>(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.at<uchar>(m-1,n)==1){temp.x=m-1;temp.y=n;quetem.push(temp);//如果這樣的話，那麼這些標記過的地區將再次在while迴圈中被掃描到，不會，因為值是127//新種子點標記為已淹沒地區，而且是目前範圍，並記錄地區號到labelImageLabelImage.at<uchar>(m-1,n)=Num;SeedImage.at<uchar>(m-1,n)=127;}else//否則上方為不可生長{up=TRUE;}}if(m>0&&n>0){if(SeedImage.at<uchar>(m-1,n-1)==1)//左上方若為可生長點則加為新種子{temp.x=m-1;temp.y=n-1;quetem.push(temp);//新種子點標記為已淹沒地區，即下一個迴圈中以127來標識不再掃描，而且是目前範圍LabelImage.at<uchar>(m-1,n-1)=Num;SeedImage.at<uchar>(m-1,n-1)=127;}else//否則左上方為不可生長{upleft=TRUE;}}if(m<row-1){if(SeedImage.at<uchar>(m+1,n)==1)//下方若為可生長點則加為新種子{temp.x=m+1;temp.y=n;quetem.push(temp);//新種子點標記為已淹沒地區，而且是目前範圍LabelImage.at<uchar>(m+1,n)=Num;SeedImage.at<uchar>(m+1,n)=127;}else//否則下方為不可生長{down=TRUE;}}if(m<(row-1)&&n<(col-1)){if(SeedImage.at<uchar>(m+1,n+1)==1)//下方若為可生長點則加為新種子{temp.x=m+1;temp.y=n+1;quetem.push(temp);//新種子點標記為已淹沒地區，而且是目前範圍LabelImage.at<uchar>(m+1,n+1)=Num;SeedImage.at<uchar>(m+1,n+1)=127;}else//否則下方為不可生長{downright=TRUE;}}if(n<col-1){if(SeedImage.at<uchar>(m,n+1)==1)//右方若為可生長點則加為新種子{temp.x=m;temp.y=n+1;quetem.push(temp);//新種子點標記為已淹沒地區，而且是目前範圍LabelImage.at<uchar>(m,n+1)=Num;SeedImage.at<uchar>(m,n+1)=127;}else//否則右方為不可生長{right=TRUE;}}if(m>0&&n<(col-1)){if(SeedImage.at<uchar>(m-1,n+1)==1)//右上方若為可生長點則加為新種子{temp.x=m-1;temp.y=n+1;quetem.push(temp);//新種子點標記為已淹沒地區，而且是目前範圍LabelImage.at<uchar>(m-1,n+1)=Num;SeedImage.at<uchar>(m-1,n+1)=127;}else//否則右上方為不可生長{upright=TRUE;}}if(n>0){if(SeedImage.at<uchar>(m,n-1)==1)//左方若為可生長點則加為新種子{temp.x=m;temp.y=n-1;quetem.push(temp);//新種子點標記為已淹沒地區LabelImage.at<uchar>(m,n-1)=Num;SeedImage.at<uchar>(m,n-1)=127;}else//否則左方為不可生長{left=TRUE;}}if(m<(row-1)&&n>0){if(SeedImage.at<uchar>(m+1,n-1)==1)//左下方若為可生長點則加為新種子{temp.x=m+1;temp.y=n-1;quetem.push(temp);//新種子點標記為已淹沒地區LabelImage.at<uchar>(m+1,n-1)=Num;SeedImage.at<uchar>(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.at<uchar>(m,n)].push(temp);//這兩句中我把Num-1改成了Num...pan's codes...SeedCounts[Num-1][OriginalImage.at<uchar>(m,n)]++;}     }//while結束，掃描到quetem為空白而止。也就是對應所有的節點都得到不可生長為止（或者是周圍的點要麼不可生長，要麼已生長）}//if結束}}//在上述過程中，如果標記的點為0則表示，沒有掃描到的點，或者表明不是輸入的種子點//這裡相當於是找seedimage傳過來的初始地區的分水嶺界線的所有的點；並且用標號記錄每個地區，同時集水盆的邊緣點進入隊列。//上面是找集水盆的程式。同時也是連通地區。////////////////////////////////////漫水//////////////////////////////////////////////////////////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.at<uchar>(m-1,n))//上方若未處理，表示沒有標號，應該在輸入前已經作過初始化為0//本函數中在開頭也作過初始化{temp.x=m-1;temp.y=n;LabelImage.at<uchar>(m-1,n)=i+1;//上方點標記為已淹沒地區//注意到這個標記是與掃描點的地區號相同，一定在這個標號所屬的地區嗎？是的//這樣在下一輪至少會掃描到這個點，確保不遺漏，但是下一輪的處理會使它合理//歸類嗎？問題還有這樣標記並沒有一定將它加入到種子隊列。也就是說它//只是被淹沒而不能向上淹沒。只有滿足下述可生長條件才行。if(OriginalImage.at<uchar>(m-1,n)<=WaterLevel)//上方若為可生長點則加入當前隊列，當前高度的隊列{vque[i][WaterLevel].push(temp);}else//否則加入OriginalImage[m-1][n]對應的灰階級的隊列，為什嗎？{vque[i][OriginalImage.at<uchar>(m-1,n)].push(temp);SeedCounts[i][OriginalImage.at<uchar>(m-1,n)]++;}}}if(m<row-1){if(!LabelImage.at<uchar>(m+1,n))//下方若未處理{temp.x=m+1;temp.y=n;LabelImage.at<uchar>(m+1,n)=i+1;//下方點標記為已淹沒地區if(OriginalImage.at<uchar>(m+1,n)<=WaterLevel)//下方若為可生長點則加入當前隊列{vque[i][WaterLevel].push(temp);}else//否則加入OriginalImage[m+1][n]級隊列{vque[i][OriginalImage.at<uchar>(m+1,n)].push(temp);SeedCounts[i][OriginalImage.at<uchar>(m+1,n)]++;}}}if(n<col-1){if(!LabelImage.at<uchar>(m,n+1))//右邊若未處理{temp.x=m;temp.y=n+1;LabelImage.at<uchar>(m,n+1)=i+1;//右邊點標記為已淹沒地區if(OriginalImage.at<uchar>(m,n+1)<=WaterLevel)//右邊若為可生長點則加入當前隊列{vque[i][WaterLevel].push(temp);}else//否則加入OriginalImage[m][n+1]級隊列{vque[i][OriginalImage.at<uchar>(m,n+1)].push(temp);SeedCounts[i][OriginalImage.at<uchar>(m,n+1)]++;}}}if(n>0){if(!LabelImage.at<uchar>(m,n-1))//左邊若未處理{temp.x=m;temp.y=n-1;LabelImage.at<uchar>(m,n-1)=i+1;//左邊點標記為已淹沒地區if(OriginalImage.at<uchar>(m,n-1)<=WaterLevel)//左邊若為可生長點則加入當前隊列{vque[i][WaterLevel].push(temp);}else//否則加入OriginalImage[m][n-1]級隊列{vque[i][OriginalImage.at<uchar>(m,n-1)].push(temp);SeedCounts[i][OriginalImage.at<uchar>(m,n-1)]++;}}}}//while迴圈結束SeedCounts[i][WaterLevel]=vque[i][WaterLevel].size();}//if結束}//for迴圈結束}//while迴圈結束}//for迴圈結束while(!vque.empty()){pque=vque.back();delete[] pque;vque.pop_back();}while(!SeedCounts.empty()){array=SeedCounts.back();delete[] array;SeedCounts.pop_back();}}