最近需要做一个图像分割的程序，查了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 CvWSNodewatershed的具体实现中，用队列来实现对像素点的存储。（我本来还以为会用什么很高级的数据结构或者很复杂的算法，等真正的看完代码才发现，用的也就是很基本的数据结构和想法，只是实现起来有些复杂，需要考虑的细节比较多）这是定义的两个结构体，CvWSNode用来存储队列中的点，其中的next元素用来指向队列中的下一个点，CvWSQueuev是队列的头指针。
{
struct CvWSNode* next;
int mask_ofs;
int img_ofs;
}
CvWSNode;
typedef struct CvWSQueuev //这个仅仅是一个头指针
{
CvWSNode* first;
CvWSNode* last;
}
CvWSQueue;
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。
后面紧接着一个两层循环的函数，这个两层循环是为了把与边界点相邻的元素都找出来，然后以这些点作为出发点。
后面的算法，紧接着就是一个不断进入队列和离开队列的过程，如果是靠近边界点，就赋值成和边界点一样的像素值，周围的点有非零的点，就赋值成同样的灰度值，如果是边界点。这样，从不同的边界点延伸出的区域，就具有了不同的颜色，也就完成了图像的分割。