最近需要做一個影象分割的程式，查了opencv的原始碼，發現opencv裡實現的影象分割一共有兩個方法，watershed和mean-shift演算法。這兩個演算法的具體實現都在segmentation.cpp檔案內。

watershed（分水嶺演算法）方法是一種基於邊界點的分割演算法。我想好好的研究一下， 網上找了一些部落格和教程，感覺也就泛泛的解釋了一下實驗的流程，具體演算法的執行過程並不清楚，又把原始論文拿出來看了看，看完了以後也不太清晰，索性把opencv的原始碼挑出來分析一下。

首先，寫一個影象分割的小程式。程式碼如下：

#include "stdafx.h"
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"

#include <iostream>

using namespace cv;
using namespace std;

Vec3b RandomColor(int value);  //生成隨機顏色函式

int main( int argc, char* argv[] )
{
		//src = imread( "E://ylab//DY//ETHZShapeClasses-V1.2//Bottles//green.jpg" );
	//src = imread( "E://ylab//DY//ETHZShapeClasses-V1.2//Mugs//jazzburger.jpg" );
	//src = imread( "E://1.png" );

	//Mat image=imread( "E://ylab//DY//ETHZShapeClasses-V1.2//Bottles//green.jpg" );
	Mat image=imread( "E://ylab//DY//ETHZShapeClasses-V1.2//bottles//green.jpg" );//載入RGB彩色影象
	imshow("Source Image",image);

	//灰度化，濾波，Canny邊緣檢測
	Mat imageGray,imageCanny;
	cvtColor(image,imageGray,CV_RGB2GRAY);//灰度轉換
	GaussianBlur(imageGray,imageGray,Size(5,5),2);   //高斯濾波
	imshow("Gray Image",imageGray); 
	Canny(imageGray,imageCanny,40,100);  
	imshow("Canny Image",imageCanny);

	//查詢輪廓
	vector<vector<Point>> contours;  
	vector<Vec4i> hierarchy;  
	findContours(imageCanny,contours,hierarchy,RETR_LIST,CHAIN_APPROX_SIMPLE,Point());  
	Mat imageContours=Mat::zeros(image.size(),CV_8UC1);  //輪廓	
	Mat marks(image.size(),CV_32S);   //Opencv分水嶺第二個矩陣引數
	marks=Scalar::all(0);
	int index = 0;
	int compCount = 0;
	for( ; index >= 0; index = hierarchy[index][0], compCount++ ) 
	{
		//對marks進行標記，對不同區域的輪廓進行編號，相當於設定注水點，有多少輪廓，就有多少注水點
		drawContours(marks, contours, index, Scalar::all(compCount+1), 1, 8, hierarchy);
		drawContours(imageContours,contours,index,Scalar(255),1,8,hierarchy);  
	}

	//我們來看一下傳入的矩陣marks裡是什麼東西
	Mat marksShows;
	convertScaleAbs(marks,marksShows);
	imshow("marksShow",marksShows);
	imshow("輪廓",imageContours);
	watershed(image,marks);

	//我們再來看一下分水嶺演算法之後的矩陣marks裡是什麼東西
	Mat afterWatershed;
	convertScaleAbs(marks,afterWatershed);
	imshow("After Watershed",afterWatershed);

	//對每一個區域進行顏色填充
	Mat PerspectiveImage=Mat::zeros(image.size(),CV_8UC3);
	for(int i=0;i<marks.rows;i++)
	{
		for(int j=0;j<marks.cols;j++)
		{
			int index=marks.at<int>(i,j);
			if(marks.at<int>(i,j)==-1)
			{
				PerspectiveImage.at<Vec3b>(i,j)=Vec3b(255,255,255);
			}			 
			else
			{
				PerspectiveImage.at<Vec3b>(i,j) =RandomColor(index);
			}
		}
	}
	imshow("After ColorFill",PerspectiveImage);

	//分割並填充顏色的結果跟原始影象融合
	Mat wshed;
	addWeighted(image,0.4,PerspectiveImage,0.6,0,wshed);
	imshow("AddWeighted Image",wshed);

	waitKey();
}

Vec3b RandomColor(int value)    //生成隨機顏色函式
{
	value=value%255;  //生成0~255的隨機數
	RNG rng;
	int aa=rng.uniform(0,value);
	int bb=rng.uniform(0,value);
	int cc=rng.uniform(0,value);
	return Vec3b(aa,bb,cc);
}
 

裡面基本上都是呼叫opencv現成的函式，實現watershed演算法，演算法的執行結果可以看一下：

從程式碼和最終結果裡可以大致的看出演算法的流程：

1.進行灰度化

2.高斯濾波以消除噪聲的干擾

3.用canny運算元檢測邊緣

4.用findcontours查詢輪廓

5.利用輪廓特徵，實現影象分割

然後，把watershed的原始碼單獨拿出來分析一下，我做了一些必要的註釋：

typedef struct CvWSNode
{
    struct CvWSNode* next;
    int mask_ofs;
    int img_ofs;
}
CvWSNode;

typedef struct CvWSQueuev     //這個僅僅是一個頭指標
{
    CvWSNode* first;
    CvWSNode* last;
}
CvWSQueue;
 

watershed的具體實現中，用佇列來實現對畫素點的儲存。（我本來還以為會用什麼很高階的資料結構或者很複雜的演算法，等真正的看完程式碼才發現，用的也就是很基本的資料結構和想法，只是實現起來有些複雜，需要考慮的細節比較多）這是定義的兩個結構體，CvWSNode用來儲存佇列中的點，其中的next元素用來指向佇列中的下一個點，CvWSQueuev是佇列的頭指標。

static CvWSNode*
icvAllocWSNodes( CvMemStorage* storage )     //申請一段連續的記憶體空間，並將記憶體空間連線起來
{
    CvWSNode* n = 0;
    
    int i, count = (storage->block_size - sizeof(CvMemBlock))/sizeof(*n) - 1;

    n = (CvWSNode*)cvMemStorageAlloc( storage, count*sizeof(*n) );
    for( i = 0; i < count-1; i++ )
        n[i].next = n + i + 1;
    n[count-1].next = 0;

    return n;
}
 

這一個函式的作用就是申請一段連續的記憶體空間，用來儲存畫素點

CV_IMPL void
cvWatershed( const CvArr* srcarr, CvArr* dstarr )
{
    const int IN_QUEUE = -2;
    const int WSHED = -1;
    const int NQ = 256;
    cv::Ptr<CvMemStorage> storage;
    
    CvMat sstub, *src;
    CvMat dstub, *dst;
    CvSize size;
    CvWSNode* free_node = 0, *node;
    CvWSQueue q[NQ];
    int active_queue;
    int i, j;
    int db, dg, dr;
    int* mask;
    uchar* img;
    int mstep, istep;
    int subs_tab[513];

    // MAX(a,b) = b + MAX(a-b,0)
    #define ws_max(a,b) ((b) + subs_tab[(a)-(b)+NQ])
    // MIN(a,b) = a - MAX(a-b,0)
    #define ws_min(a,b) ((a) - subs_tab[(a)-(b)+NQ])

    #define ws_push(idx,mofs,iofs)  \
    {                               \
        if( !free_node )            \
            free_node = icvAllocWSNodes( storage );\
        node = free_node;           \
        free_node = free_node->next;\
        node->next = 0;             \
        node->mask_ofs = mofs;      \
        node->img_ofs = iofs;       \
        if( q[idx].last )           \
            q[idx].last->next=node; \
        else                        \
            q[idx].first = node;    \
        q[idx].last = node;         \
    }

    #define ws_pop(idx,mofs,iofs)   \
    {                               \
        node = q[idx].first;        \
        q[idx].first = node->next;  \
        if( !node->next )           \
            q[idx].last = 0;        \
        node->next = free_node;     \
        free_node = node;           \
        mofs = node->mask_ofs;      \
        iofs = node->img_ofs;       \
    }

    #define c_diff(ptr1,ptr2,diff)      \
    {                                   \
        db = abs((ptr1)[0] - (ptr2)[0]);\
        dg = abs((ptr1)[1] - (ptr2)[1]);\
        dr = abs((ptr1)[2] - (ptr2)[2]);\
        diff = ws_max(db,dg);           \
        diff = ws_max(diff,dr);         \
        assert( 0 <= diff && diff <= 255 ); \
    }

    src = cvGetMat( srcarr, &sstub );
    dst = cvGetMat( dstarr, &dstub );

    if( CV_MAT_TYPE(src->type) != CV_8UC3 )
        CV_Error( CV_StsUnsupportedFormat, "Only 8-bit, 3-channel input images are supported" );

    if( CV_MAT_TYPE(dst->type) != CV_32SC1 )
        CV_Error( CV_StsUnsupportedFormat,
            "Only 32-bit, 1-channel output images are supported" );
    
    if( !CV_ARE_SIZES_EQ( src, dst ))
        CV_Error( CV_StsUnmatchedSizes, "The input and output images must have the same size" );

    size = cvGetMatSize(src);                //圖片的大小size   有height和width分量
    storage = cvCreateMemStorage();

    istep = src->step;
    img = src->data.ptr;
    mstep = dst->step / sizeof(mask[0]);      //step為每行元素的位元組數，
    mask = dst->data.i;

    memset( q, 0, NQ*sizeof(q[0]) );    //void *memset(void *s,int c,size_t n)  總的作用：將已開闢記憶體空間 s 的首 n 個位元組的值設為值 c。

    for( i = 0; i < 256; i++ )
        subs_tab[i] = 0;
    for( i = 256; i <= 512; i++ )
        subs_tab[i] = i - 256;

    // draw a pixel-wide border of dummy "watershed" (i.e. boundary) pixels
    for( j = 0; j < size.width; j++ )
        mask[j] = mask[j + mstep*(size.height-1)] = WSHED;    //把mask的上邊界和下邊界畫素賦值為-1

    // initial phase: put all the neighbor pixels of each marker to the ordered queue -
    // determine the initial boundaries of the basins
    for( i = 1; i < size.height-1; i++ )
    {
        img += istep; mask += mstep;
        mask[0] = mask[size.width-1] = WSHED;     //把mask的左邊界和右邊界畫素賦值為-1

        for( j = 1; j < size.width-1; j++ )     //把和邊界點相鄰的點連成連結串列，進行儲存
        {
            int* m = mask + j;
            if( m[0] < 0 ) m[0] = 0;
            if( m[0] == 0 && (m[-1] > 0 || m[1] > 0 || m[-mstep] > 0 || m[mstep] > 0) )
            {
                uchar* ptr = img + j*3;
                int idx = 256, t;
                if( m[-1] > 0 )
                    c_diff( ptr, ptr - 3, idx );
                if( m[1] > 0 )
                {
                    c_diff( ptr, ptr + 3, t );
                    idx = ws_min( idx, t );
                }
                if( m[-mstep] > 0 )
                {
                    c_diff( ptr, ptr - istep, t );
                    idx = ws_min( idx, t );
                }
                if( m[mstep] > 0 )
                {
                    c_diff( ptr, ptr + istep, t );
                    idx = ws_min( idx, t );
                }
                assert( 0 <= idx && idx <= 255 );
                ws_push( idx, i*mstep + j, i*istep + j*3 );
                m[0] = IN_QUEUE;    //  IN_QUEUE=-2  在序列裡，就把Mark矩陣邊緣點相鄰點的位置賦值為-2，
            }
        }
    }

    // find the first non-empty queue
    for( i = 0; i < NQ; i++ )
        if( q[i].first )
            break;

    // if there is no markers, exit immediately
    if( i == NQ )
        return;

    active_queue = i;
    img = src->data.ptr;
    mask = dst->data.i;

    // recursively fill the basins     遞迴的填充盆地
    for(;;)
    {
        int mofs, iofs;
        int lab = 0, t;
        int* m;
        uchar* ptr;
        
        if( q[active_queue].first == 0 )
        {
            for( i = active_queue+1; i < NQ; i++ )
                if( q[i].first )
                    break;
            if( i == NQ )
                break;
            active_queue = i;
        }

        ws_pop( active_queue, mofs, iofs );

        m = mask + mofs;
        ptr = img + iofs;
        t = m[-1];
        if( t > 0 ) lab = t;
        t = m[1];
        if( t > 0 )
        {
            if( lab == 0 ) lab = t;
            else if( t != lab ) lab = WSHED;
        }
        t = m[-mstep];
        if( t > 0 )
        {
            if( lab == 0 ) lab = t;
            else if( t != lab ) lab = WSHED;
        }
        t = m[mstep];
        if( t > 0 )
        {
            if( lab == 0 ) lab = t;
            else if( t != lab ) lab = WSHED;
        }
        assert( lab != 0 );
        m[0] = lab;
        if( lab == WSHED )
            continue;

        if( m[-1] == 0 )
        {
            c_diff( ptr, ptr - 3, t );
            ws_push( t, mofs - 1, iofs - 3 );
            active_queue = ws_min( active_queue, t );
            m[-1] = IN_QUEUE;
        }
        if( m[1] == 0 )
        {
            c_diff( ptr, ptr + 3, t );
            ws_push( t, mofs + 1, iofs + 3 );
            active_queue = ws_min( active_queue, t );
            m[1] = IN_QUEUE;
        }
        if( m[-mstep] == 0 )
        {
            c_diff( ptr, ptr - istep, t );
            ws_push( t, mofs - mstep, iofs - istep );
            active_queue = ws_min( active_queue, t );
            m[-mstep] = IN_QUEUE;
        }
        if( m[mstep] == 0 )
        {
            c_diff( ptr, ptr + istep, t );
            ws_push( t, mofs + mstep, iofs + istep );
            active_queue = ws_min( active_queue, t );
            m[mstep] = IN_QUEUE;
        }
    }
}

ws_max和我說ws_min這兩個函式說白了就是求最大值和最小值，但是這裡用來求最大最小值的方法還真是奇特，我還是第一次見到，其中使用了subs_tab，然後利用這個陣列的值就可以求極大值和極小值。

ws_push和ws_pop這兩個函式，從名字上就可以看出來，一個是往佇列的末端加入元素，另一個是從佇列的頂端取出元素。

c_diff是為了計算圖片RGB三個分量的梯度，並選出一個梯度最大的。

傳入watershed函式共有兩個引數，第一個引數是原始的3通道彩色圖片，第二個引數是Mark矩陣，這個矩陣裡，每一條輪廓都被表示了出來，而且，每一條輪廓的灰度值是不同的，這些輪廓，就被該演算法成為種子，演算法就是從這些輪廓上的點開始計算並分割成不同的區域的。除了邊界點外的其他區域，畫素值都賦值為0。

後面緊接著一個兩層迴圈的函式，這個兩層迴圈是為了把與邊界點相鄰的元素都找出來，然後以這些點作為出發點。

後面的演算法，緊接著就是一個不斷進入佇列和離開佇列的過程，如果是靠近邊界點，就賦值成和邊界點一樣的畫素值，周圍的點有非零的點，就賦值成同樣的灰度值，如果是邊界點。這樣，從不同的邊界點延伸出的區域，就具有了不同的顏色，也就完成了影象的分割。