使用OpenCV和C++實現的分水嶺演算法(Watershed)

pamxy發表於2013-05-29
OpenCVC++演算法

轉自：http://blog.csdn.net/twowind/article/details/8988282

分水嶺演算法(watershed)是一種比較基本的數學形態學分割演算法，其基本思想是將灰度影象轉換為梯度影象，將梯度值看作高低起伏的山嶺，將區域性極小值及其鄰域看作一個“集水盆”。設想一個個“集水盆”中存在積水，且水位不斷生長，淹沒低度較低的地方，當水漫過程停止後，影象就可以被分割成幾塊連通區域。

分水嶺演算法有不同的實現方法。本文要實現的是通過人為標註一些種子點，將這些種子點看作集水盆的底部，利用區域增長的方法，完成影象的分割。試圖實現OpenCV中cv::watershed函式的功能，經過測試，與OpenCV相比分割結果相似，但效能差很多。（前者32ms左右，後者8ms左右）。

OpenCV函式的執行結果：（OpenCV函式對分割邊緣也做了處理，我寫的那個程式沒有）


程式執行結果：


參考：

http://wenku.baidu.com/view/d1fde240336c1eb91a375d95.html

http://blog.csdn.net/fdl19881/article/details/6749976

[cpp] view plaincopyprint?
  1. #include <opencv2/imgproc/imgproc.hpp>  
  2. #include <opencv2/objdetect/objdetect.hpp>  
  3. #include <opencv2/highgui/highgui.hpp>  
  4. #include<vector>   
  5. #include<iostream>   
  6. #include<queue>   
  7. #include<fstream>   
  8. cv::Mat marker_mask;  
  9. cv::Mat g_markers;  
  10. cv::Mat img0, img, img_gray,wshed;  
  11. cv::Point_<int> prev_pt(-1,-1);  
  12. using std::vector;  
  13. using std::queue;  
  14. static void my_watershed(cv::Mat img,cv::Mat& markers,int comp_count);  
  15. static void mouse_event(int event,int x, int y,int flags, void*)  
  16. {  
  17.     if(img.rows==0)  
  18.         return;  
  19.      if(event==CV_EVENT_LBUTTONUP||!(flags&CV_EVENT_FLAG_LBUTTON))  
  20.          prev_pt=cv::Point_<int>(-1,-1);  
  21.      else if(event==CV_EVENT_LBUTTONDOWN)  
  22.          prev_pt=cv::Point2i(x,y);  
  23.      else if(event==CV_EVENT_MOUSEMOVE&&(flags&CV_EVENT_FLAG_LBUTTON))  
  24.      {  
  25.          cv::Point2i pt(x,y);  
  26.          if(prev_pt.x<0)  
  27.              prev_pt=pt;  
  28.          cv::line(marker_mask,prev_pt,pt,cv::Scalar(255,255,255),1,8,0);  
  29.          cv::line(img,prev_pt,pt,cv::Scalar(255,255,255),1,8,0);  
  30.          prev_pt=pt;  
  31.          cv::imshow("image",img);  
  32.      }  
  33. }  
  34. int main()  
  35. {  
  36.     img0=cv::imread("Lenna.png",1);  
  37.     img=img0.clone();  
  38.     CvRNG rng = cvRNG(-1);   
  39.     img_gray=img0.clone();  
  40.     wshed=img0.clone();  
  41.     marker_mask=cv::Mat(cv::Size(img0.cols,img0.rows),8,1);  
  42.     g_markers=cv::Mat(cv::Size(img0.cols,img0.rows),CV_32S,1);  
  43.     cv::cvtColor(img,marker_mask,CV_BGR2GRAY);  
  44.     cv::cvtColor(marker_mask,img_gray,CV_GRAY2BGR);  
  45.     for(int i=0;i<marker_mask.rows;i++)  
  46.         for(int j=0;j<marker_mask.cols;j++)  
  47.             marker_mask.at<unsigned char>(i,j)=0;  
  48.     for(int i=0;i<g_markers.rows;i++)  
  49.         for(int j=0;j<g_markers.cols;j++)  
  50.             g_markers.at<int>(i,j)=0;  
  51.     cv::imshow("image",img);  
  52.     cv::imshow("watershed transform",wshed);  
  53.     cv::setMouseCallback("image",mouse_event,0);  
  54.     for(;;)  
  55.     {  
  56.         int c=cv::waitKey(0);  
  57.         if((char)c==27)  
  58.             break;  
  59.         if((char)c=='r')  
  60.         {  
  61.             for(int i=0;i<marker_mask.rows;i++)  
  62.                 for(int j=0;j<marker_mask.cols;j++)  
  63.                     marker_mask.at<unsigned char>(i,j)=0;  
  64.             img0.copyTo(img);  
  65.             cv::imshow("image",img);  
  66.         }  
  67.         if((char)c=='w'||(char)c==' ')  
  68.         {  
  69.             vector<vector<cv::Point>> contours;  
  70.             CvMat* color_tab=0;  
  71.             int comp_count=0;  
  72.             cv::findContours(marker_mask,contours,CV_RETR_CCOMP,CV_CHAIN_APPROX_SIMPLE,cv::Point(0,0));  
  73.             for(int i=0;i<g_markers.rows;i++)  
  74.                 for(int j=0;j<g_markers.cols;j++)  
  75.                         g_markers.at<int>(i,j)=0;  
  76.             vector<vector<cv::Point> >::iterator iter=contours.begin();  
  77.             for(int i=0;i<(int)contours.size();i++)  
  78.             {  
  79.                 cv::drawContours(g_markers,contours,i,cv::Scalar::all(comp_count+1),  
  80.                     1,8,vector<cv::Vec4i>());  
  81.                 comp_count++;  
  82.             }  
  83.           
  84.             if(comp_count==0)  
  85.                 continue;  
  86.             color_tab=cvCreateMat(1,comp_count,CV_8UC3);  
  87.             for(int i=0;i<comp_count;i++)  
  88.             {  
  89.                 uchar* ptr=color_tab->data.ptr+i*3;  
  90.                 ptr[0]=(uchar)(cvRandInt(&rng)%180+50);  
  91.                 ptr[1]=(uchar)(cvRandInt(&rng)%180+50);  
  92.                 ptr[2]=(uchar)(cvRandInt(&rng)%180+50);  
  93.             }  
  94.             cv::Mat temp=g_markers.clone();  
  95.   
  96.             double t=(double)cvGetTickCount();  
  97.             //my_watershed(img0,g_markers,comp_count);  
  98.             cv::watershed(img0,g_markers);  
  99.             t=(double)cvGetTickCount()-t;  
  100.             std::cout<<"exec time= "<<t/(cvGetTickFrequency()*1000.)<<std::endl;  
  101.         for(int i=0;i<g_markers.rows;i++)  
  102.             for(int j=0;j<g_markers.cols;j++)  
  103.             {  
  104.                 int idx=g_markers.at<int>(i,j);  
  105.                 uchar* dst=&wshed.at<uchar>(i,j*3);  
  106.                 if(idx==-1)  
  107.                     dst[0]=dst[1]=dst[2]=(uchar)255;  
  108.                 else if(idx<=0||idx>comp_count)  
  109.                     dst[0]=dst[1]=dst[2]=(uchar)8;  
  110.                 else{  
  111.                     uchar* ptr=color_tab->data.ptr+(idx-1)*3;  
  112.                     dst[0]=ptr[0];dst[1]=ptr[1];dst[2]=ptr[2];  
  113.                 }  
  114.             }  
  115.             cv::addWeighted(wshed,0.5,img_gray,0.5,0,wshed);  
  116.             cv::imshow("watershed transform",wshed);  
  117.             cvReleaseMat(&color_tab);  
  118.         }  
  119.     }  
  120.     return 0;  
  121. }  
  122. static void my_watershed(cv::Mat img0,cv::Mat& markers,int comp_count)  
  123. {  
  124.     cv::Mat gray=cv::Mat(cv::Size(img0.rows,img0.cols),8,1);  
  125.     cv::cvtColor(img0,gray,CV_BGR2GRAY);  
  126.     cv::Mat imge=cv::Mat(cv::Size(img0.rows,img0.cols),8,1);  
  127.     cv::Sobel(gray,imge,CV_8U,1,1);  
  128.     vector<queue<cv::Point2i>*>Labeleddata;//影象中各連通區域的點  
  129.     queue<cv::Point2i>* pque;//某連通區域已包含的點  
  130.     queue<cv::Point2i> quetem; //用於提取某連通區域中輸入種子點中的初始種子點  
  131.     vector<int*> SeedCounts;  
  132.     int* Array;  
  133.     cv:: Point2i temp;  
  134.     int row=imge.rows,col=imge.cols;  
  135.     cv::Mat marker_saved=markers.clone();  
  136.     bool up,down,right,left,uplef,uprig,downlef,downrig;  
  137.     int m,n;  
  138.     for(int i=0;i<comp_count;i++)  
  139.     {  
  140.         Array=new int[256];  
  141.         SeedCounts.push_back(Array);//統計某waterlevel的各個連通區域中種子點的個數  
  142.         pque=new queue<cv::Point2i>[256];   
  143.         Labeleddata.push_back(pque);//儲存該連通區域中種子生長所得的點        
  144.     }  
  145.     for(int i=0;i<row;i++)  
  146.         for(int j=0;j<col;j++)  
  147.         {  
  148.             if(markers.at<int>(i,j)>0)  
  149.             {  
  150.                 temp.x=i;  
  151.                 temp.y=j;  
  152.                 quetem.push(temp);  
  153.                 int num=markers.at<int>(i,j);  
  154.                 markers.at<int>(i,j)=-1;//該點已處理,其他種子點生長時將繞過該點  
  155.                 while(!quetem.empty())  
  156.                 {  
  157.                     up=down=right=left=uplef=uprig=downlef=downrig=false;  
  158.                     temp=quetem.front(); //提取出一個點，在該點的八連通區域內尋找可生長點  
  159.                     m=temp.x;  
  160.                     n=temp.y;  
  161.                     quetem.pop();  
  162.   
  163.                     if(m-1>=0)//若上方可生長則新增為新種子  
  164.                     {  
  165.                         if(markers.at<int>(m-1,n)==num)  
  166.                         {  
  167.                             temp.x=m-1;  
  168.                             temp.y=n;  
  169.                             quetem.push(temp);  
  170.                             markers.at<int>(m-1,n)=-1;  
  171.                         }  
  172.                         else{  
  173.                             up=true;  
  174.                         }  
  175.                     }  
  176.                     if(m-1>=0&&n-1>=0)  
  177.                     {  
  178.                         if(markers.at<int>(m-1,n-1)==num)  
  179.                         {  
  180.                             temp.x=m-1;  
  181.                             temp.y=n-1;  
  182.                             quetem.push(temp);  
  183.                             markers.at<int>(m-1,n-1)=-1;  
  184.                         }  
  185.                         else{  
  186.                             uplef=true;  
  187.                         }  
  188.                     }  
  189.                     if(m+1<=row-1)  
  190.                     {  
  191.                         if(markers.at<int>(m+1,n)==num)  
  192.                         {  
  193.                             temp.x=m+1;  
  194.                             temp.y=n;  
  195.                             quetem.push(temp);  
  196.                             markers.at<int>(m+1,n)=-1;  
  197.                         }  
  198.                         else{  
  199.                             down=true;  
  200.                         }  
  201.                     }  
  202.                     if(m+1<=row-1&&n+1<=col-1)  
  203.                     {  
  204.                         if(markers.at<int>(m+1,n+1)==num)  
  205.                         {  
  206.                             temp.x=m+1;  
  207.                             temp.y=n+1;  
  208.                             quetem.push(temp);  
  209.                             markers.at<int>(m+1,n+1)=-1;  
  210.                         }  
  211.                         else{  
  212.                             downrig=true;  
  213.                         }  
  214.                     }  
  215.                     if(n+1<=col-1)  
  216.                     {  
  217.                         if(markers.at<int>(m,n+1)==num)  
  218.                         {  
  219.                             temp.x=m;  
  220.                             temp.y=n+1;  
  221.                             quetem.push(temp);  
  222.                             markers.at<int>(m,n+1)=-1;  
  223.                         }  
  224.                         else{  
  225.                             right=true;  
  226.                         }  
  227.                     }  
  228.                     if(m-1>=0&&n+1<=col-1)  
  229.                     {  
  230.                         if(markers.at<int>(m-1,n+1)==num)  
  231.                         {  
  232.                             temp.x=m-1;  
  233.                             temp.y=n+1;  
  234.                             quetem.push(temp);  
  235.                             markers.at<int>(m-1,n+1)=-1;  
  236.                         }  
  237.                         else{  
  238.                             uprig=true;  
  239.                         }  
  240.                     }  
  241.                     if(n-1>=0)  
  242.                     {  
  243.                         if(markers.at<int>(m,n-1)==num)  
  244.                         {  
  245.                             temp.x=m;  
  246.                             temp.y=n-1;  
  247.                             quetem.push(temp);  
  248.                             markers.at<int>(m,n-1)=-1;  
  249.                         }  
  250.                         else{  
  251.                             left=true;  
  252.                         }  
  253.                     }  
  254.                     if(m+1<=row-1&&n-1>=0)  
  255.                     {  
  256.                         if(markers.at<int>(m+1,n-1)==num)  
  257.                         {  
  258.                             temp.x=m+1;  
  259.                             temp.y=n-1;  
  260.                             quetem.push(temp);  
  261.                             markers.at<int>(m+1,n-1)=-1;  
  262.                         }  
  263.                         else{  
  264.                             downlef=true;  
  265.                         }  
  266.                     }  
  267.                     //八連通區域中有未標記點，則該點屬於初始種子點  
  268.                     if(up||down||right||left||uplef||downlef||uprig||downrig)  
  269.                     {  
  270.                         temp.x=m;  
  271.                         temp.y=n;  
  272.                         Labeleddata[comp_count-1][imge.at<uchar>(m,n)].push(temp);  
  273.                         SeedCounts[comp_count-1][imge.at<uchar>(m,n)]++;  
  274.                     }  
  275.                 }  
  276.             }  
  277.         }  
  278.         bool active;  
  279.         int waterlevel;  
  280.         for(waterlevel=0;waterlevel<180;waterlevel++)  
  281.         {  
  282.             active=true;  
  283.             while(active) //當1-count_com個連通區域都無可生長點時結束迴圈  
  284.             {  
  285.                 active=false;  
  286.                 for(int i=0;i<comp_count;i++)//將區域i中將waterlevel梯度以下的點用於區域增長  
  287.                 {  
  288.                     if(!Labeleddata[i][waterlevel].empty())//區域增長，經過多次迭代，直至該區域，該waterlevel無可生長點。  
  289.                     {  
  290.                         active=true;      
  291.                         while(SeedCounts[i][waterlevel]>0) //SeedCount中保留了前一輪生長後各區域，各waterlevel的種子點個數，本輪生長結束後，將根據Labeleddata中的元素個數更新  
  292.                         {  
  293.                             SeedCounts[i][waterlevel]--;  
  294.                             temp=Labeleddata[i][waterlevel].front();  
  295.                             Labeleddata[i][waterlevel].pop();  
  296.                             m=temp.x;  
  297.                             n=temp.y;  
  298.                             int num=marker_saved.at<int>(m,n);  
  299.                             if(m-1>=0)  
  300.                             {  
  301.                                 if(!marker_saved.at<int>(m-1,n))//上方點未處理過  
  302.                                 {  
  303.                                     temp.x=m-1;  
  304.                                     temp.y=n;  
  305.                                     marker_saved.at<int>(m-1,n)=num;  
  306.                                     if(imge.at<uchar>(m-1,n)<=waterlevel)  
  307.                                         Labeleddata[i][waterlevel].push(temp);  
  308.                                     else{  
  309.                                         Labeleddata[i][imge.at<uchar>(m-1,n)].push(temp); //本次生長不處理，可能在waterlevel變化到某值時再用於生長  
  310.                                         SeedCounts[i][imge.at<uchar>(m-1,n)]++;  
  311.                                     }  
  312.                                 }  
  313.                             }  
  314.                             if(m+1<=row-1)  
  315.                             {  
  316.                                 if(!marker_saved.at<int>(m+1,n))  
  317.                                 {  
  318.                                     temp.x=m+1;  
  319.                                     temp.y=n;  
  320.                                     marker_saved.at<int>(m+1,n)=num;  
  321.                                     if(imge.at<uchar>(m+1,n)<=waterlevel)  
  322.                                         Labeleddata[i][waterlevel].push(temp);  
  323.                                     else{  
  324.                                         Labeleddata[i][imge.at<uchar>(m+1,n)].push(temp);  
  325.                                         SeedCounts[i][imge.at<uchar>(m+1,n)]++;  
  326.                                     }  
  327.                                 }  
  328.                             }  
  329.                             if(n+1<=col-1)  
  330.                             {  
  331.                                 if(!marker_saved.at<int>(m,n+1))  
  332.                                 {  
  333.                                     temp.x=m;  
  334.                                     temp.y=n+1;  
  335.                                     marker_saved.at<int>(m,n+1)=num;  
  336.                                     if(imge.at<uchar>(m,n+1)<=waterlevel)  
  337.                                         Labeleddata[i][waterlevel].push(temp);  
  338.                                     else{  
  339.                                         Labeleddata[i][imge.at<uchar>(m,n+1)].push(temp);  
  340.                                         SeedCounts[i][imge.at<uchar>(m,n+1)]++;  
  341.                                     }  
  342.                                 }  
  343.                             }  
  344.                             if(n-1>=0)  
  345.                             {  
  346.                                 if(!marker_saved.at<int>(m,n-1))  
  347.                                 {  
  348.                                     temp.x=m;  
  349.                                     temp.y=n-1;  
  350.                                     marker_saved.at<int>(m,n-1)=num;  
  351.                                     if(imge.at<uchar>(m,n-1)<=waterlevel)  
  352.                                         Labeleddata[i][waterlevel].push(temp);  
  353.                                     else  
  354.                                     {  
  355.                                         Labeleddata[i][imge.at<uchar>(m,n-1)].push(temp);  
  356.                                         SeedCounts[i][imge.at<uchar>(m,n-1)]++;  
  357.                                     }  
  358.                                 }  
  359.                             }  
  360.                         }  
  361.                             SeedCounts[i][waterlevel]=Labeleddata[i][waterlevel].size();  
  362.                         }  
  363.                           
  364.                     }  
  365.             }  
  366.         }  
  367.         markers=marker_saved.clone();  
  368. }  

 

相關文章