3.5分水岭算法
3.5.1分水岭算法概述
分水岭分割方法。它是一种基于拓扑理论的数学形态学的分割方法，其基本思想是把图像看作是测地学上的拓扑地貌，图像中每一点像素的灰度值表示该点的海拔高度，每一个局部极小值及其影响区域称为集水盆，而集水盆的边界则形成分水岭。分水岭的概念和形成可以通过模拟浸入过程来说明。在每一个局部极小值表面，刺穿一个小孔，然后把整个模型慢慢浸入水中，随着浸入的加深，每一个局部极小值的影响域慢慢向外扩展，在两个集水盆汇合处构筑大坝，即形成分水岭。
分水岭比较经典的计算方法是L．Vincent于1991年在PAMI上提出的《An efficientalgorithms based on immersion simulation》。传统的分水岭分割方法，是一种基于拓扑理论的数学形态学的分割方法，其基本思想是把图像看作是测地学上的拓扑地貌，图像中每一像素的灰度值表示该点的海拔高度，每一个局部极小值及其影响区域称为集水盆地，而集水盆地的边界则形成分水岭。分水岭的概念和形成可以通过模拟浸入过程来说明。在每一个局部极小值表面，刺穿一个小孔，然后把整个模型慢慢浸人水中，随着浸入的加深，每一个局部极小值的影响域慢慢向外扩展，在两个集水盆汇合处构筑大坝如下图所示，即形成分水岭。
一般的分水岭算法会对微弱边缘，图像中的噪声，物体表面细微的灰度变化造成过度的分割。OpenCV中的分水岭算法对此进行了改进，它使用预定义的一组标记来引导对图像的分割。理解这一节，关键是理解这个标记。
然而基于梯度图像的直接分水岭算法容易导致图像的过分割，产生这一现象的原因主要是由于输入的图像存在过多的极小区域而产生许多小的集水盆地，从而导致分割后的图像不能将图像中有意义的区域表示出来。所以必须对分割结果的相似区域进行合并。
传统的分水岭算法的步骤：
第一步：找到图像的局部最低点，这个方法很多了，可以用一个内核去找，也可以一个一个比较，实现起来不难。
第二步：从最低点开始注水，水开始网上满（图像的说法就是梯度法），其中那些最低点已经被标记，不会被淹没，那些中间点是被淹没的。
第三步：找到局部最高点，就是图中3位置对应的两个点。
第四步：这样基于局部最小值，和找到的局部最大值，就可以分割图像了。
因为传统分水岭算法存在过分割的不足，OpenCV提供了一种改进的分水岭算法，使用一系列预定义标记来引导图像分割的定义方式。使用OpenCV的分水岭算法cv::wathershed，需要输入一个标记图像，图像的像素值为32位有符号正数（CV_32S类型），每个非零像素代表一个标签。它的原理是对图像中部分像素做标记，表明它的所属区域是已知的。分水岭算法可以根据这个初始标签确定其他像素所属的区域。传统的基于梯度的分水岭算法和改进后基于标记的分水岭算法示意图如下图所示。
从上图可以看出，传统基于梯度的分水岭算法由于局部最小值过多造成分割后的分水岭较多。而基于标记的分水岭算法，水淹过程从预先定义好的标记图像（像素）开始，较好的克服了过度分割的不足。本质上讲，基于标记点的改进算法是利用先验知识来帮助分割的一种方法。因此，改进算法的关键在于如何获得准确的标记图像，即如何将前景物体与背景准确的标记出来。
基于标记点的分水岭算法的步骤：
第一步：封装分水岭算法类
第二步：获取标记图像
获取前景像素，并用255标记前景
获取背景像素，并用128标记背景，未知像素，使用0标记
合成标记图像
第三步：将原图和标记图像输入分水岭算法
第四步：显示结果
3.5.2实现分水岭算法：watershed()函数
watershed()函数讲解
C++: void watershed(InputArray image,
InputOutputArray markers)【参数】
第一个参数，image – Input 8-bit 3-channel image.
第二个参数，markers – Input/output 32-bit single-channel image (map) of markers. It should have the same size as image .
watershed()函数源代码
/*【watershed ( )源代码】*************************************************************
* @Version：OpenCV 3.0.0（Opnencv2和Opnencv3差别不大，Linux和PC的对应版本源码完全一样，均在对应的安装目录下）
* @源码路径：…\opencv\sources\modules\imgproc\src\ segmentation.cpp
* @起始行数：88行
********************************************************************************/
void cv::watershed( InputArray _src, InputOutputArray _markers )
{
// Labels for pixels
const int IN_QUEUE = -2; // Pixel visited
const int WSHED = -1; // Pixel belongs to watershed
// 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];
// 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 ); \
}
CV_Assert( src.type() == CV_8UC3 && dst.type() == CV_32SC1 );
CV_Assert( src.size() == dst.size() );
// Current pixel in input image
const uchar* img = src.ptr();
// Step size to next row in input image
int istep = int(src.step/sizeof(img[0]));
// Current pixel in mask image
int* mask = dst.ptr<int>();
// Step size to next row in mask image
int mstep = int(dst.step / sizeof(mask[0]));
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;
// 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; // boundary pixels
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) )
{
// 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 );
m[0] = IN_QUEUE;
}
}
}
// 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.ptr();
mask = dst.ptr<int>();
// recursively fill the basins
for(;;)
{
int mofs, iofs;
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 );
// 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;
}
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;
}
// Set label to current pixel in marker image
assert( lab != 0 );
m[0] = lab;
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 );
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;
}
}
}3.5.3分水岭算法实例
代码参看附件【demo1】。
【注】原图中的白色曲线是用鼠标勾勒的大致轮廓。
参考：
英文
扫一扫
2805
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
