分水岭算法实现
1.输入原图像和掩模图像，掩模图像为手动标记的线，每个线对应的在掩模图像中的灰度不一样。（在maskshows图像中可以看出。）
2.统计原图像中除边框外，每个像素与其周围上下左右四个邻居的差别，差别的定义是，两个三通道图像BGR三个通道相减的绝对值的最小值。得出四个邻居差别最小值，根据这个最值放入相应队列中。比如与左邻居在红色通道差别最小为10，则放入10这个队列中。队列最多有0—255 也就是256个
3.队列中每个元素包含，像素点在掩模图像和原图像的位置。遍历整个原图像将除边框外每个像素点情况放入对应队列中。
4.从差别最小的队列依次取出，如果某个像素点在掩模图像值大于0，说明这个像素点是标记点。如果小于等于0，且其四个邻居只有一个大于0，此像素点在掩模图像的值等于邻居在掩模图像中的值。如果此像素点四个邻居中有两个大于0，而且不相同，那么此像素点在掩模图像中的值记为-1，就是分水岭点。
5取出256个队列所有值，分水岭算法完成。
256个队列结构： http://www.xuebuyuan.com/2164320.html
漫水填充过程可以参考：http://blog.csdn.net/tangketan/article/details/39757513
参考《opencv3编程入门》书中的分水岭算法
//--------------------------------------【程序说明】-------------------------------------------
//程序描述：分水岭算法综合示例
//开发测试所用操作系统： Windows 10 64bit
//开发测试所用IDE版本：Visual Studio 2015
//开发测试所用OpenCV版本：3.1.0
//2014年06月 Created by @浅墨_毛星云
//2014年11月 Revised by @浅墨_毛星云
// 2016年3月 changed by yuyangyg
//------------------------------------------------------------------------------------------------
//---------------------------------【头文件、命名空间包含部分】----------------------------
//描述：包含程序所使用的头文件和命名空间
//------------------------------------------------------------------------------------------------
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
#include <iostream>
#include<fstream>
using namespace cv;
using namespace std;
//-----------------------------------【宏定义部分】--------------------------------------------
// 描述：定义一些辅助宏
//------------------------------------------------------------------------------------------------
#define WINDOW_NAME1 "【程序窗口1】" //为窗口标题定义的宏
#define WINDOW_NAME2 "【分水岭算法效果图】" //为窗口标题定义的宏
//-----------------------------------【全局函变量声明部分】--------------------------------------
//描述：全局变量的声明
//-----------------------------------------------------------------------------------------------
Mat g_maskImage, g_srcImage;
Point prevPt(-1, -1);
//-----------------------------------【全局函数声明部分】--------------------------------------
//描述：全局函数的声明
//-----------------------------------------------------------------------------------------------
static void ShowHelpText();
static void on_Mouse(int event, int x, int y, int flags, void*);
void watershed1(InputArray _src, InputOutputArray _markers);
//-----------------------------------【main( )函数】--------------------------------------------
//描述：控制台应用程序的入口函数，我们的程序从这里开始执行
//-----------------------------------------------------------------------------------------------
int main(int argc, char** argv)
{
//【0】改变console字体颜色
system("color 6F");
//【0】显示帮助文字
ShowHelpText();
//【1】载入原图并显示，初始化掩膜和灰度图
g_srcImage = imread("1.jpg", 1);
imshow(WINDOW_NAME1, g_srcImage);
Mat srcImage, grayImage;
g_srcImage.copyTo(srcImage);
cvtColor(g_srcImage, g_maskImage, COLOR_BGR2GRAY);
cvtColor(g_maskImage, grayImage, COLOR_GRAY2BGR);
g_maskImage = Scalar::all(0);//元素初始化为零
//【2】设置鼠标回调函数
setMouseCallback(WINDOW_NAME1, on_Mouse, 0);
//【3】轮询按键，进行处理
while (1)
{
//获取键值
int c = waitKey(0);
//若按键键值为ESC时，退出
if ((char)c == 27)
break;
//按键键值为2时，恢复源图
if ((char)c == '2')
{
g_maskImage = Scalar::all(0);
srcImage.copyTo(g_srcImage);
imshow("image", g_srcImage);
}
//若检测到按键值为1或者空格，则进行处理
if ((char)c == '1' || (char)c == ' ')
{
//定义一些参数
int i, j, compCount = 0;
vector<vector<Point> > contours;
vector<Vec4i> hierarchy;
//寻找轮廓
findContours(g_maskImage, contours, hierarchy, CV_RETR_CCOMP, CV_CHAIN_APPROX_SIMPLE);//找到图像中的边缘轮廓
//轮廓为空时的处理
if (contours.empty())
continue;
//拷贝掩膜
Mat maskImage(g_maskImage.size(), CV_32S);
maskImage = Scalar::all(0);
//循环绘制出轮廓
for (int index = 0; index >= 0; index = hierarchy[index][0], compCount++)
drawContours(maskImage, contours, index, Scalar::all(compCount + 1), -1, 8, hierarchy, INT_MAX);//在和srcImage同样大小的空白图像上画轮廓
//compCount为零时的处理
if (compCount == 0)
continue;
//生成随机颜色
vector<Vec3b> colorTab;
for (i = 0; i < compCount; i++)
{
int b = theRNG().uniform(0, 255);
int g = theRNG().uniform(0, 255);
int r = theRNG().uniform(0, 255);
colorTab.push_back(Vec3b((uchar)b, (uchar)g, (uchar)r));
}
//计算处理时间并输出到窗口中
double dTime = (double)getTickCount();
Mat maskShows;
convertScaleAbs(maskImage*20, maskShows);
imshow("maskShows", maskShows);
watershed1(srcImage, maskImage);//调用分水岭算法
dTime = (double)getTickCount() - dTime;
printf("\t处理时间 = %gms\n", dTime*1000. / getTickFrequency());
//双层循环，将分水岭图像遍历存入watershedImage中
Mat watershedImage(maskImage.size(), CV_8UC3);
for (i = 0; i < maskImage.rows; i++)
for (j = 0; j < maskImage.cols; j++)
{
int index = maskImage.at<int>(i, j);
if (index == -1)
watershedImage.at<Vec3b>(i, j) = Vec3b(255, 255, 255);
else if (index <= 0 || index > compCount)
watershedImage.at<Vec3b>(i, j) = Vec3b(0, 0, 0);
else
watershedImage.at<Vec3b>(i, j) = colorTab[index - 1];
}
imshow("watershedImage", watershedImage);
//混合灰度图和分水岭效果图并显示最终的窗口
watershedImage = watershedImage*0.5 + grayImage*0.5;
imshow(WINDOW_NAME2, watershedImage);
}
}
return 0;
}
//-----------------------------------【onMouse( )函数】---------------------------------------
//描述：鼠标消息回调函数
//-----------------------------------------------------------------------------------------------
static void on_Mouse(int event, int x, int y, int flags, void*)
{
//处理鼠标不在窗口中的情况
if (x < 0 || x >= g_srcImage.cols || y < 0 || y >= g_srcImage.rows)
return;
//处理鼠标左键相关消息
if (event == CV_EVENT_LBUTTONUP || !(flags & CV_EVENT_FLAG_LBUTTON))
prevPt = Point(-1, -1);
else if (event == CV_EVENT_LBUTTONDOWN)
prevPt = Point(x, y);
//鼠标左键按下并移动，绘制出白色线条
else if (event == CV_EVENT_MOUSEMOVE && (flags & CV_EVENT_FLAG_LBUTTON))
{
Point pt(x, y);
if (prevPt.x < 0)
prevPt = pt;
line(g_maskImage, prevPt, pt, Scalar::all(255), 5, 8, 0);
line(g_srcImage, prevPt, pt, Scalar::all(255), 5, 8, 0);
prevPt = pt;
imshow(WINDOW_NAME1, g_srcImage);
}
}
//-----------------------------------【ShowHelpText( )函数】----------------------------------
// 描述：输出一些帮助信息
//----------------------------------------------------------------------------------------------
static void ShowHelpText()
{
//输出欢迎信息和OpenCV版本
//printf("\n\n\t\t\t非常感谢购买《OpenCV3编程入门》一书！\n");
//printf("\n\n\t\t\t此为本书OpenCV2版的第77个配套示例程序\n");
printf("\n\n\t\t\t 当前使用的OpenCV版本为：" CV_VERSION);
printf("\n\n ----------------------------------------------------------------------------\n");
//输出一些帮助信息
printf("\n\n\n\t欢迎来到【分水岭算法】示例程序~\n\n");
printf("\t请先用鼠标在图片窗口中标记出大致的区域，\n\n\t然后再按键【1】或者【SPACE】启动算法。"
"\n\n\t按键操作说明: \n\n"
"\t\t键盘按键【1】或者【SPACE】- 运行的分水岭分割算法\n"
"\t\t键盘按键【2】- 恢复原始图片\n"
"\t\t键盘按键【ESC】- 退出程序\n\n\n");
}
//#include "precomp.hpp"
/****************************************************************************************\
* Watershed *
\****************************************************************************************/
namespace cv
{
// A node represents a pixel to label
struct WSNode //像素节点
{
int next;
int mask_ofs; //相对于mask图像的偏移
int img_ofs; //相对于原图像的偏移
};
// Queue for WSNodes //像素节点队列
struct WSQueue
{
WSQueue() { first = last = 0; }
int first, last;
};
static int
allocWSNodes(std::vector<WSNode>& storage)//分配空间
{
int sz = (int)storage.size();
int newsz = MAX(128, sz * 3 / 2);//空间扩展为128，或者sz的1.5倍？？？
storage.resize(newsz);
if (sz == 0)
{
storage[0].next = 0;
sz = 1;
}
for (int i = sz; i < newsz - 1; i++)
storage[i].next = i + 1;
storage[newsz - 1].next = 0;
return sz;
}
}
void watershed1(InputArray _src, InputOutputArray _markers)
{
// Labels for pixels
const int IN_QUEUE = -2; // Pixel visited 加入到队列记做-2
const int WSHED = -1; // Pixel belongs to watershed 属于mask记做-1
// possible bit values = 2^8
const int NQ = 256; //队列长度，和灰度等级相同
Mat src = _src.getMat(), dst = _markers.getMat();
Size size = src.size();
// Vector of every created node
std::vector<WSNode> storage;
int free_node = 0, node;
// Priority queue of queues of nodes
// from high priority (0) to low priority (255)
WSQueue q[NQ]; //数组每一个元素都是队列
// Non-empty queue with highest priority
int active_queue;
int i, j;
// Color differences
int db, dg, dr;
int subs_tab[513]; //255+256=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])
// Create a new node with offsets mofs and iofs in queue idx
#define ws_push(idx,mofs,iofs) \
{ \
if( !free_node ) \
free_node = allocWSNodes( storage );\
node = free_node; \
free_node = storage[free_node].next;\
storage[node].next = 0; \
storage[node].mask_ofs = mofs; \
storage[node].img_ofs = iofs; \
if( q[idx].last ) \
storage[q[idx].last].next=node; \
else \
q[idx].first = node; \
q[idx].last = node; \
}
// Get next node from queue idx
#define ws_pop(idx,mofs,iofs) \
{ \
node = q[idx].first; \
q[idx].first = storage[node].next; \
if( !storage[node].next ) \
q[idx].last = 0; \
storage[node].next = free_node; \
free_node = node; \
mofs = storage[node].mask_ofs; \
iofs = storage[node].img_ofs; \
}
// Get highest absolute channel difference in diff
#define c_diff(ptr1,ptr2,diff) \
{ \
db = std::abs((ptr1)[0] - (ptr2)[0]);\
dg = std::abs((ptr1)[1] - (ptr2)[1]);\
dr = std::abs((ptr1)[2] - (ptr2)[2]);\
diff = ws_max(db,dg); \
diff = ws_max(diff,dr); \
assert( 0 <= diff && diff <= 255 ); \
} //找出两个像素点在BGR三个通道中相差最大的值记为diff
CV_Assert(src.type() == CV_8UC3 && dst.type() == CV_32SC1);
CV_Assert(src.size() == dst.size()); //输入图像要求是CV_8UC3,格式CV_32SC1
// Current pixel in input image 步长 = 一行字节数 / sizeof(像素数据类型)
const uchar* img = src.ptr(); //函数Mat.ptr(row)其中row表示行索引.
// Step size to next row in input image
int istep = int(src.step / sizeof(img[0])); //原图像到下一行需要步数。
cout << "src.step="<<src.step<<" sizeof(img[0]))="<<sizeof(img[0])<<endl;//img[0]是uchar型故为一个字节
// Current pixel in mask image
int* mask = dst.ptr<int>();//mask[i]表示第i行的头指针
// Step size to next row in mask image
int mstep = int(dst.step / sizeof(mask[0]));//mask图像到下一行需要的步数CV_32SC1为所以占四个字节
cout <<"dst.step="<<dst.step<<" sizeof(mask[0])="<< sizeof(mask[0])<<endl;
for (i = 0; i < 256; i++)
subs_tab[i] = 0;
for (i = 256; i <= 512; i++)
subs_tab[i] = i - 256; //用来max和min
// 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的首行和末行画成分水岭
// initial phase: put all the neighbor pixels of each marker to the ordered queue -
// determine the initial boundaries of the basins
//除边框外，把标记 4 邻域的像素加入队列，q有256个优先级，对应距离0-255，所以具体加入哪个
//队列就看距离最小idx了
for (i = 1; i < size.height - 1; i++)
{
img += istep; mask += mstep;
mask[0] = mask[size.width - 1] = WSHED; // boundary pixels mask的左右边缘记为分水岭
for (j = 1; j < size.width - 1; j++)
{
int* m = mask + j;
if (m[0] < 0) m[0] = 0;//若该点为负数，记做0
if (m[0] == 0 && (m[-1] > 0 || m[1] > 0 || m[-mstep] > 0 || m[mstep] > 0))//分水岭点的四周都不是分水岭点
{
// Find smallest difference to adjacent markers找出与分水岭点在单一像素通道上，相差最小的值
const 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);
}
// Add to according queue
assert(0 <= idx && idx <= 255);
ws_push(idx, i*mstep + j, i*istep + j * 3); //将该点在img和mask中的坐标（一维表示）存储在q[idx]队列中
m[0] = IN_QUEUE; //已经加到队列，下面一行遇到这个点怎么处理？？
}
}
}
// find the first non-empty queue找到第一个非空队列
for (i = 0; i < NQ; i++)
if (q[i].first)
break; //比如 q[0]表示所有标记区附近与标记区距离为0的点
// if there is no markers, exit immediately
if (i == NQ) //与邻居像素差为i,如果为256说明所有队列为空
return;
active_queue = i;
img = src.ptr();
mask = dst.ptr<int>();
// recursively fill the basins递归地填满聚水盆
for (;;)
{
int mofs, iofs; // 将二维图像线性化后图像像素的坐标 mask_offset 和 img_offset 的缩写
int lab = 0, t;
int* m;
const uchar* ptr;
// Get non-empty queue with highest priority
// Exit condition: empty priority queue
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;
// 如果这个灰度上的队列处理完了，就继续找下一个非空队列
}
// Get next node
ws_pop(active_queue, mofs, iofs); // 从q[active_queue]队列中取出一个结点数据
// Calculate pointer to current pixel in input and marker image
m = mask + mofs;
ptr = img + iofs;
// Check surrounding pixels for labels
// to determine label for current pixel
t = m[-1]; // Left 左邻居的标记
if (t > 0) lab = t;
t = m[1]; // Right
if (t > 0)
{
if (lab == 0) lab = t;
else if (t != lab) lab = WSHED; //左右邻居都大于零且左右邻居的标记不一样，这该点为分水岭点若只有一个邻居大于零，该邻居的标记赋给这个点。
} //输入图像时要给mask的每一个标记值不同，可以有maskshow图像看出。
t = m[-mstep]; // Top
if (t > 0)
{
if (lab == 0) lab = t;
else if (t != lab) lab = WSHED;
}
t = m[mstep]; // Bottom
if (t > 0)
{
if (lab == 0) lab = t;
else if (t != lab) lab = WSHED;// 如果该像素点的标记类型和四个邻居中任意一个像素标记类型都 > 0 且不同，则为分水岭
}
// Set label to current pixel in marker image
// 因为标记点要么是初始种子点，要么是初始阶段延伸的种子点的邻接点
// 该点一定存在一个邻接点是标记点，所以lab一定会赋值一次，不为 0
assert(lab != 0);
m[0] = lab;// 若lab > 0 ，则该点被周围的标记点扩充；若lab = -1（WSHED）,则该点定义为分水岭，继续下一个循环
if (lab == WSHED)
continue;
// Add adjacent, unlabeled pixels to corresponding queue
if (m[-1] == 0)
{
c_diff(ptr, ptr - 3, t); //计算最小差值
ws_push(t, mofs - 1, iofs - 3); // 将m[-1]这一未标记的点扩充为标记点，进队
active_queue = ws_min(active_queue, t); // 判断，若t < 当前处理的队列active_queue值，则下一次循环中处理q[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;
}
}
}
CV_IMPL void cvWatershed(const CvArr* _src, CvArr* _markers)
{
cv::Mat src = cv::cvarrToMat(_src), markers = cv::cvarrToMat(_markers);
cv::watershed(src, markers);
}
效果图
参考《opencv3编程入门》
opencv3.10源码
http://blog.csdn.net/tangketan/article/details/39757513
http://blog.csdn.net/u011375993/article/details/46793655
http://www.xuebuyuan.com/2164320.html