图像分割 区域分割之分水岭算法

学习DIP第60天 转载请标明本文出处： ，出于尊重文章作者的劳动，转载请标明出处！文章代码已托管，欢迎共同开发：https://github.com/Tony-Tan/DIPpro

开篇废话

今天已经是第60篇博客了，这六十篇每一篇平均要两天左右，所以，在过去的四个月学到了这么多知识，想想挺开心，但学的越多就会发现自己不会的越多。从小学到大学，这么多年一直以学习为主要工作但学习又有很多阶段，对于通用知识，比如小学的语文数学此观点不适用，对于一些专业性较强的知识，感觉会有两个很主要的阶段，感觉自己目前处于入门阶段，由于数字图像涉及数学的知识较多，还有信号和信息论的知识，所以，感觉还是比较难学科，不然招聘公司也不会一年几十万的养着图像处理工程师。 下面的图是本人的一点点见解，只是自己总结的，没有实践，也没有科学依据，不喜勿喷：

我感觉图像处理能分成三个阶段，或者更多，第一阶段的人很多，听说这行前景好或者工资高的人，多半会学习点图像的知识，比如彩色空间啊，了解下OpenCV啊等等，还有一些属于纯属上课被逼无奈的，比如我们学校就对电子信息类专业和通信类开数字图像处理这门课，而且我当时考试还挂了。。。。这部分大家会接触一些名词，一些简单算法，有用心的同学可能会实现下代码，这阶段风景不错，而且好多名词可以拿来忽悠HR或者忽悠导师，得到一份不错的工作，或者做点导师的项目。 这个阶段多半使用现成的函数库，了解了基本算法或者听别人说一些算法，自己来跑结果，这个阶段Matlab的用户量较大。 第1阶段后半期也就是平台期，这个阶段是做了一段时间有一定算法使用基础和代码能力的工程师，多半在各企业负责最底层的图像算法编写，在他们面前是一座大山，和一个路口，继续做图像还是转管理。 如果继续选择图像，就会面临一个峭壁，具体是算法底层的数学，说的数学，总是困难的，爬这个峭壁的动力可以是挣更多的钱，或者爱好，因为一旦到达阶段2，就能成为首席图像处理工程师，到任何需要图像处理的公司，都能够独当一面，这些人已经到了不缺钱的地步，而且在行业内一定有一定名气。 第二阶段的一旦到达，，可以说是事业的平稳期，或者巅峰，不缺钱，还能独自决定一些技术层面上的事，指挥手下阶段1的员工工作，这个阶段面临的也是一个选择，就是靠这个吃一辈子饭，绝对没问题，再有就是向更高的境界冲击。 第三阶段没有尽头，因为能促使进入阶段三的动力只有爱好，这部分挣的钱可能没有阶段2挣的多，而且难度更大，看不到尽头，所以这部分属于探索阶段，在这个阶段上看到的一些，都能推动未来学科的发展，所以这个阶段的人多半是在实验室和数学中忙碌一生，然后几十年后出现在各大论文和教材中。 以上属于个人猜想或者愚见，想法幼稚，仅供参考。

算法描述

今天废话太多了，哈哈，其实上面的可以新开一篇博客单独写出来，但觉得，首先自己年轻视野狭窄，第二水平太低，所以当做废话夹在本文中。 今天说分水岭算法，分水岭算法族是一个思想引出的一些列实现方法。 算法内容：首先将二维灰度图像抽象为三维地形，横纵坐标表示经纬度，灰度表示高度，这样就产生了一个地形图，地形中有高山有盆地，我们的任务就是找到盆地并且给盆地划定地盘。 算法过程，将最低点（灰度）作为起始点，开始从下注水，想泉水一样，水位逐渐上升，所有点都是上下透水的，但不会横向流动，每产生一个新的积水盆地时标记这个盆地使之与原有盆地区别开，当两个盆地内的水要汇集的时候，在此点建立一个水坝，将水分开，此水坝无限高，水永远不会溢出，依次建立所有水坝，最后的水坝就是所谓的分水岭或分水线，也就是得到的分割线。 算法步骤：

检测图像中所有的盆地，并依次标记，（此时最终的分割区域数已确定）将盆地周围邻域的像素按照灰度值大小入队（队列，数据结构，先进先出）按照灰度值顺序将像素出队，进行判断 判断1：如果像素周围只存在一个盆地标记，将此像素划分到此盆地 判断2：如果像素周围存在两个以及以上盆地，次像素为分水岭。

此算法描述只是分水岭算法之一，冈萨雷斯的书中使用的是形态学的描述方式，但与上述过程原理一致，感兴趣者可自行查阅。 解释下第一步，第一步是本算法的核心，此步骤首先要确定盆地还是非盆地，可以抽象成下面几种模型：

图像中的红色框内为非初始盆地，或叫做盆底，绿色为盆底，对于第二种盆地的确定需要使用图的深度或者广度优先搜索，判断所有候选点，来判断是第二种模型还是第三种。

代码/* typedef struct PriQueueNode_ PriQueueHead; typedef struct NLevelPriQueueNode_ * NLevelPriQueue; typedef struct NLevelPriQueueNode_ NLevelPriNode; typedef struct ExPix_ Pix_Label; *___ ____________________ *| P |->| NLevelPriQueue | *|—| |——————–| *| r |->| NLevelPriQueue | *|—| |——————–| *| i |->| NLevelPriQueue | *|—| |——————–| *| Q |->| NLevelPriQueue | *|—| |——————–| *| u |->| NLevelPriQueue | *|—| |——————–| *| e |->| NLevelPriQueue | *|—| |——————–| *| u |->| NLevelPriQueue | *|—| |——————–| *| e |->| NLevelPriQueue | *|___| |____________________| struct NLevelPriNode_{ int x; int y; NLevelPriQueue next; }; struct PriQueueNode_{ int nodeNum; NLevelPriQueue head; NLevelPriQueue tail; }; struct ExPix_{ int grayvalue; int label; }; */inQueue(PriQueueHead* priQueue,int gray_level,int x,int y){priQueue[gray_level].nodeNum++;///malloc new nodeNLevelPriNode* newNode=(NLevelPriNode *)malloc(sizeof(NLevelPriNode));newNode->x=x;newNode->y=y;newNode->next=NULL;if(priQueue[gray_level].head==NULL){priQueue[gray_level].head=newNode;priQueue[gray_level].tail=newNode;}else{priQueue[gray_level].tail->next=newNode;priQueue[gray_level].tail=newNode;}},int x,int y){src[y*width+x]=-value;for(int j=-1;j<2;j++)for(int i=-1;i<2;i++)if(j+y>=0&&i+x>=0&&j+y<height&&i+x<width&&(i!=0||j!=0)){if(src[(j+y)*width+i+x]<value&&src[(j+y)*width+i+x]>0)return 0;}for(int j=-1;j<2;j++)for(int i=-1;i<2;i++)if(j+y>=0&&i+x>=0&&j+y<height&&i+x<width&&(i!=0||j!=0)){if(src[(j+y)*width+i+x]==value){return(isPan(src,width, height, value, i+x, j+y));}}return 1;},int x,int y){src[y*width+x]=-value;for(int j=-1;j<2;j++)for(int i=-1;i<2;i++)if(j+y>=0&&i+x>=0&&j+y<height&&i+x<width&&(i!=0||j!=0)){if(src[(j+y)*width+i+x]==value){repairPan(src, width, height, value, x+i, y+j);}}},int x,int y){dst[y*width+x]=255.0;for(int j=-1;j<2;j++)for(int i=-1;i<2;i++)if(j+y>=0&&i+x>=0&&j+y<height&&i+x<width&&(i!=0||j!=0)){if(src[(j+y)*width+x+i]==value&&dst[(j+y)*width+x+i]==0.0)setMinimal(src,dst,width, height, value, x+i,y+j);}},int x,int y){dst[y*width+x]=127.0;for(int j=-1;j<2;j++)for(int i=-1;i<2;i++)if(j+y>=0&&i+x>=0&&j+y<height&&i+x<width&&(i!=0||j!=0)){if(src[(j+y)*width+x+i]==value&&dst[(j+y)*width+x+i]==0.0)setUnMinimal(src,dst,width, height, value, x+i,y+j);}}//数据类型从double到intvoid Double2Int(double *src,int* dst,int width,int height){for(int i=0;i<width*height;i++)dst[i]=(int)src[i]+1;}//寻找极小值，包括单点极小值和平底锅void findMinimal(double *src,double *dst,int width,int height){int *temp=(int *)malloc(sizeof(int)*width*height);double *dsttemp=(double *)malloc(sizeof(double)*width*height);Zero(dsttemp, width, height);Double2Int(src, temp, width, height);int lessthan=0;int equ=0;double min=findMatrixMin(src, width, height);for(int i=0;i<width*height;i++)if(src[i]==min)dsttemp[i]=255.0;for(int j=0;j<height;j++){for(int i=0;i<width;i++){lessthan=0;equ=0;int pix=temp[j*width+i];if(dsttemp[j*width+i]==0.0){for(int m=-1;m<2;m++)for(int n=-1;n<2;n++)if(j+m>=0&&i+n>=0&&j+m<height&&i+n<width){if(m!=0||n!=0){if(temp[(j+m)*width+i+n]<pix)lessthan=1;if(temp[(j+m)*width+i+n]==pix)equ=1;}}if(equ==1&&lessthan==0){if(isPan(temp, width, height,pix, i, j)){//repairPan(temp,width, height, -pix, i,j);setMinimal(temp,dsttemp, width, height, -pix, i, j);}else {repairPan(temp,width, height, -pix, i,j);setUnMinimal(temp,dsttemp, width, height, pix, i, j);}}if(lessthan==1)dsttemp[j*width+i]=127.0;if(0==lessthan&&0==equ)dsttemp[j*width+i]=255.0;}}}matrixCopy(dsttemp, dst, width, height);free(dsttemp);free(temp);}//标记极小值的label，从-1开始向下增长 -2 -3 -4 -5 -6 -7…..void LabelMinimal(double * src,Pix_Label* dst,int width,int height,int x,int y,int label){dst[y*width+x].label=label;for(int i=-1;i<2;i++){for(int j=-1;j<2;j++)if(x+i>=0&&x+i<width&&y+j>=0&&y+j<height&&(i!=0||j!=0))if(src[(y+j)*width+x+i]==255.0&&dst[(y+j)*width+x+i].label==0){LabelMinimal(src, dst, width, height, x+i, y+j, label);}}}//初始化label数组，此数组与图像数组多加了labelvoid InitLabelMat(double *src,Pix_Label* dst,double *mask,int width,int height){for(int i=0;i<width*height;i++){dst[i].grayvalue=src[i];dst[i].label=0;}int label_minimal=-1;for(int j=0;j<height;j++)for(int i=0;i<width;i++){if(mask[j*width+i]==255.0&&dst[j*width+i].label==0){LabelMinimal(mask, dst, width,height,i,j, label_minimal);label_minimal–;}}}//初始化队列头数组void InitPriQueue(PriQueueHead* priQueue,Pix_Label *srclabel,int width,int height){for(int i=0;i<GRAY_LEVEL;i++){priQueue[i].head=NULL;priQueue[i].nodeNum=0;priQueue[i].tail=NULL;}int inqueue=0;for(int j=0;j<height;j++)for(int i=0;i<width;i++){inqueue=0;if(srclabel[j*width+i].label==0){for(int m=-1;m<2;m++)for(int n=-1;n<2;n++){if(m+j>=0&&m+j<height&&n+i>=0&&n+i<width&&(m!=0||n!=0))if(srclabel[(m+j)*width+n+i].label<0)inqueue=1;}if(inqueue){inQueue(priQueue,srclabel[j*width+i].grayvalue,i,j);srclabel[j*width+i].label=INQUEUE;}}}//for(int i=0;i<GRAY_LEVEL;i++)// printf(“g:%d NumofNode:%d\n”,i,priQueue[i].nodeNum);}int PirQueueisEmpty(PriQueueHead* priqueue){int sum=0;for(int i=0;i<GRAY_LEVEL;i++)sum+=priqueue[i].nodeNum;return !sum;}NLevelPriNode* outQueue(PriQueueHead* priqueue){NLevelPriNode* node=NULL;if(!PirQueueisEmpty(priqueue))for(int i=0;i<GRAY_LEVEL;i++)if(priqueue[i].nodeNum!=0){node=priqueue[i].head;priqueue[i].head=node->next;priqueue[i].nodeNum–;break;}return node;}void findWaterShed(Pix_Label * srclabel,PriQueueHead* priqueue,int width,int height){NLevelPriNode* node=outQueue(priqueue);while(node!=NULL){int y=node->y;int x=node->x;//printf(“x:%d y:%d \n”,x,y);int label=0;int isWatershed=0;for(int j=-1;j<2;j++)for(int i=-1;i<2;i++){if(j+y>=0&&j+y<height&&i+x>=0&&i+x<width&&(j!=0||i!=0)){if(srclabel[(j+y)*width+x+i].label<0){if(label==0)label=srclabel[(j+y)*width+x+i].label;else if(label!=srclabel[(j+y)*width+x+i].label){isWatershed=1;}}}}if(isWatershed)srclabel[y*width+x].label=WATERSHED;else if(label<0){srclabel[y*width+x].label=label;for(int j=-1;j<2;j++)for(int i=-1;i<2;i++){if(j+y>=0&&j+y<height&&i+x>=0&&i+x<width&&(j!=0||i!=0)){if(srclabel[(j+y)*width+x+i].label==0){inQueue(priqueue, srclabel[(j+y)*width+i+x].grayvalue, i+x, j+y);srclabel[(j+y)*width+i+x].label=INQUEUE;}}}}else if(label==0&&isWatershed==0){srclabel[y*width+x].label=0;}free(node);node=outQueue(priqueue);}}//meyer分水岭方法void MeyerWatershed(double *src,double *dst,int width,int height){double *dst_temp=(double *)malloc(sizeof(double)*width*height);Zero(dst_temp, width, height);Pix_Label * srclabel=(Pix_Label*)malloc(sizeof(Pix_Label)*width*height);PriQueueHead priqueue[GRAY_LEVEL];double *minimal=(double *)malloc(sizeof(double)*width*height);Zero(minimal, width, height);findMinimal(src, minimal, width, height);InitLabelMat(src, srclabel, minimal, width, height);InitPriQueue(priqueue, srclabel, width, height);/*for(int i=0;i<GRAY_LEVEL;i++){NLevelPriNode *node=priqueue[i].head;printf(“%d:”,i);while (node!=NULL) {printf(“x:%d,y:%d “,node->x,node->y);node=node->next;}printf(“\n”);}*/findWaterShed(srclabel, priqueue, width, height);for(int i=0;i<width*height;i++)if(srclabel[i].label==WATERSHED)dst_temp[i]=255.0;matrixCopy(dst_temp, dst, width, height);free(dst_temp);free(srclabel);free(minimal);}实现效果

总结积极思考造成积极人生，消极思考造成消极人生。