C++实现分水岭算法（Watershed Algorithm）

2023-09-17 14:00:04 150

分水岭分割方法（WatershedSegmentation），是一种基于拓扑理论的数学形态学的分割方法，其基本思想是把图像看作是测地学上的拓扑地貌，图像中每一点像素的灰度值表示该点的海拔高度，每一个局部极小值及其影响区域称为集水盆，而集水盆的边界则形成分水岭。分水岭的概念和形成可以通过模拟浸入过程来说明。在每一个局部极小值表面，刺穿一个小孔，然后把整个模型慢慢浸入水中，随着浸入的加深，每一个局部极小值的影响域慢慢向外扩展，在两个集水盆汇合处构筑大坝，即形成分水岭。
分水岭的计算过程是一个迭代标注过程。分水岭比较经典的计算方法是L.Vincent提出的。在该算法中，分水岭计算分两个步骤，一个是排序过程，一个是淹没过程。首先对每个像素的灰度级进行从低到高排序，然后在从低到高实现淹没过程中，对每一个局部极小值在h阶高度的影响域采用先进先出（FIFO）结构进行判断及标注。
分水岭变换得到的是输入图像的集水盆图像，集水盆之间的边界点，即为分水岭。显然，分水岭表示的是输入图像极大值点。因此，为得到图像的边缘信息，通常把梯度图像作为输入图像，即：
grad(f(x,y))=((f(x-1,y)-f(x+1,y))^2+(f(x,y-1)-f(x,y+1))^2)^0.5
式中，f(x,y)表示原始图像，grad(.)表示梯度运算。
分水岭算法对微弱边缘具有良好的响应，图像中的噪声、物体表面细微的灰度变化，都会产生过度分割的现象。但同时应当看出，分水岭算法对微弱边缘具有良好的响应，是得到封闭连续边缘的保证的。另外，分水岭算法所得到的封闭的集水盆，为分析图像的区域特征提供了可能。
为消除分水岭算法产生的过度分割，通常可以采用两种处理方法，一是利用先验知识去除无关边缘信息。二是修改梯度函数使得集水盆只响应想要探测的目标。
为降低分水岭算法产生的过度分割，通常要对梯度函数进行修改，一个简单的方法是对梯度图像进行阈值处理，以消除灰度的微小变化产生的过度分割。即：
g(x,y)=max(grad(f(x,y)),gθ)
式中，gθ表示阈值。
程序可采用方法：用阈值限制梯度图像以达到消除灰度值的微小变化产生的过度分割，获得适量的区域，再对这些区域的边缘点的灰度级进行从低到高排序，然后在从低到高实现淹没的过程，梯度图像用Sobel算子计算获得。对梯度图像进行阈值处理时，选取合适的阈值对最终分割的图像有很大影响，因此阈值的选取是图像分割效果好坏的一个关键。缺点：实际图像中可能含有微弱的边缘，灰度变化的数值差别不是特别明显，选取阈值过大可能会消去这些微弱边缘。

下面用C++实现分水岭算法：

#define_USE_MATH_DEFINES

#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include

usingnamespacestd;



typedefvoidGVVoid;
typedefboolGVBoolean;
typedefcharGVChar;
typedefunsignedcharGVByte;
typedefshortGVInt16;
typedefunsignedshortGVUInt16;
typedefintGVInt32;
typedefunsignedintGVUInt32;
typedeflonglongGVInt64;
typedefunsignedlonglongGVUInt64;
typedeffloatGVFloat32;
typedefdoubleGVFloat64;

constGVBooleanGV_TRUE=true;
constGVBooleanGV_FALSE=false;
constGVByteGV_BYTE_MAX=UCHAR_MAX;
constGVInt32GV_INT32_MAX=INT_MAX;
constGVInt32GV_INT32_MIX=INT_MIN;
constGVInt64GV_INT64_MAX=LLONG_MAX;
constGVInt64GV_INT64_MIN=LLONG_MIN;
constGVFloat32GV_FLOAT32_MAX=FLT_MAX;
constGVFloat32GV_FLOAT32_MIN=FLT_MIN;
constGVFloat64GV_FLOAT64_MAX=DBL_MAX;
constGVFloat64GV_FLOAT64_MIN=DBL_MIN;

classGVPoint;



classGVPoint{

public:
GVInt32x;
GVInt32y;

public:
GVPoint():x(0),y(0){}
GVPoint(constGVPoint&obj):x(obj.x),y(obj.y){}
GVPoint(GVInt32x,GVInt32y):x(x),y(y){}

public:
GVBooleanoperator==(constGVPoint&right)const{
return((x==right.x)&&(y==right.y));
}
GVBooleanoperator!=(constGVPoint&right)const{
return(!(x==right.x)||!(y==right.y));
}
};



/*
*
*imagedata;
*imagewidth;
*imageheight;
*imageoflabeledwatersheds.
*/
GVVoidgvWatershed(
constGVByte*image,
GVInt32width,
GVInt32height,
GVInt32*label)
{

//Localconstants
constGVInt32WSHD=0;
constGVInt32INIT=-1;
constGVInt32MASK=-2;
constGVPointFICT_PIXEL=GVPoint(~0,~0);

//Imagestatisticsandsorting
GVInt32size=width*height;
GVInt32*image_stat=newGVInt32[GV_BYTE_MAX+1];
GVInt32*image_space=newGVInt32[GV_BYTE_MAX+1];
GVPoint*image_sort=newGVPoint[size];
::memset(image_stat,0,sizeof(GVInt32)*(GV_BYTE_MAX+1));
::memset(image_space,0,sizeof(GVInt32)*(GV_BYTE_MAX+1));
::memset(image_sort,0,sizeof(GVPoint)*size);
for(GVInt32i=0;!(i==size);++i){
image_stat[image[i]]++;
}
for(GVInt32i=0;!(i==GV_BYTE_MAX);++i){
image_stat[i+1]+=image_stat[i];
}
for(GVInt32i=0;!(i==height);++i){
for(GVInt32j=0;!(j==width);++j){
GVBytespace=image[i*width+j];
GVInt32index=image_stat[space]-(++image_space[space]);
image_sort[index].x=j;
image_sort[index].y=i;
}
}
for(GVInt32i=GV_BYTE_MAX;!(i==0);--i){
image_stat[i]-=image_stat[i-1];
}

//Watershedalgorithm
GVPoint*head=image_sort;
GVInt32space=0;
GVInt32*dist=newGVInt32[size];
GVInt32dist_cnt=0;
GVInt32label_cnt=0;
std::queuequeue;
::memset(dist,0,sizeof(GVInt32)*size);
::memset(label,~0,sizeof(GVInt32)*size);
for(GVInt32h=0;!(h==(GV_BYTE_MAX+1));++h){
head+=space;
space=image_stat[h];
for(GVInt32i=0;!(i==space);++i){
GVInt32index=head[i].y*width+head[i].x;
GVInt32index_l=((head[i].x-1)<0)?-1:((head[i].x-1)+head[i].y*width);
GVInt32index_r=!((head[i].x+1)>width)?-1:((head[i].x+1)+head[i].y*width);
GVInt32index_t=((head[i].y-1)<0)?-1:(head[i].x+(head[i].y-1)*width);
GVInt32index_b=!((head[i].y+1)>height)?-1:(head[i].x+(head[i].y+1)*width);
label[index]=MASK;
if((!(index_l<0)&&!(label[index_l]width)?-1:((top.x+1)+top.y*width);
GVInt32index_t=((top.y-1)<0)?-1:(top.x+(top.y-1)*width);
GVInt32index_b=!((top.y+1)>height)?-1:(top.x+(top.y+1)*width);
queue.pop();
if(top==FICT_PIXEL){
if(queue.empty())break;
else{
++dist_cnt;
queue.push(FICT_PIXEL);
top=queue.front();
queue.pop();
}
}
if(!(index_l<0)){
if((dist[index_l]WSHD){
if((label[index]==MASK)||(label[index]=WSHD))label[index]=label[index_l];
elseif(!(label[index]==label[index_l]))label[index]=WSHD;
}elseif(label[index]==MASK){
label[index]=WSHD;
}
}elseif((label[index_l]==MASK)&&(dist[index_l]==0)){
dist[index_l]=dist_cnt+1;
queue.push(GVPoint(top.x-1,top.y));
}
}
if(!(index_r<0)){
if((dist[index_r]WSHD){
if((label[index]==MASK)||(label[index]=WSHD))label[index]=label[index_r];
elseif(!(label[index]==label[index_r]))label[index]=WSHD;
}elseif(label[index]==MASK){
label[index]=WSHD;
}
}elseif((label[index_r]==MASK)&&(dist[index_r]==0)){
dist[index_r]=dist_cnt+1;
queue.push(GVPoint(top.x+1,top.y));
}
}
if(!(index_t<0)){
if((dist[index_t]WSHD){
if((label[index]==MASK)||(label[index]=WSHD))label[index]=label[index_t];
elseif(!(label[index]==label[index_t]))label[index]=WSHD;
}elseif(label[index]==MASK){
label[index]=WSHD;
}
}elseif((label[index_t]==MASK)&&(dist[index_t]==0)){
dist[index_t]=dist_cnt+1;
queue.push(GVPoint(top.x,top.y-1));
}
}
if(!(index_b<0)){
if((dist[index_b]WSHD){
if((label[index]==MASK)||(label[index]=WSHD))label[index]=label[index_b];
elseif(!(label[index]==label[index_b]))label[index]=WSHD;
}elseif(label[index]==MASK){
label[index]=WSHD;
}
}elseif((label[index_b]==MASK)&&(dist[index_b]==0)){
dist[index_b]=dist_cnt+1;
queue.push(GVPoint(top.x,top.y+1));
}
}
}
for(GVInt32i=0;!(i==space);++i){
GVInt32index=head[i].y*width+head[i].x;
dist[index]=0;
if(label[index]==MASK){
label_cnt++;
label[index]=label_cnt;
queue.push(head[i]);
while(!queue.empty()){
GVPointtop=queue.front();
GVInt32index_l=((top.x-1)<0)?-1:((top.x-1)+top.y*width);
GVInt32index_r=!((top.x+1)>width)?-1:((top.x+1)+top.y*width);
GVInt32index_t=((top.y-1)<0)?-1:(top.x+(top.y-1)*width);
GVInt32index_b=!((top.y+1)>height)?-1:(top.x+(top.y+1)*width);
queue.pop();
if(!(index_l<0)&&(label[index_l]==MASK)){
queue.push(GVPoint(top.x-1,top.y));
label[index_l]=label_cnt;
}
if(!(index_r<0)&&(label[index_r]==MASK)){
queue.push(GVPoint(top.x+1,top.y));
label[index_r]=label_cnt;
}
if(!(index_t<0)&&(label[index_t]==MASK)){
queue.push(GVPoint(top.x,top.y-1));
label[index_t]=label_cnt;
}
if(!(index_b<0)&&(label[index_b]==MASK)){
queue.push(GVPoint(top.x,top.y+1));
label[index_b]=label_cnt;
}
}
}
}
}

//Releaseresources
delete[]image_stat;
delete[]image_space;
delete[]image_sort;
delete[]dist;
}

以上就是本文的全部内容，希望对大家的学习有所帮助，也希望大家多多支持毛票票。

C++实现分水岭算法（Watershed Algorithm）

热门推荐

随机推荐