灰度图像--图像分割 区域分割之分水岭算法

学习DIP第60天

转载请标明本文出处：http://blog.csdn.net/tonyshengtan ，出于尊重文章作者的劳动，转载请标明出处！文章代码已托管，欢迎共同开发：https://github.com/Tony-Tan/DIPpro

开篇废话

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

下面的图是本人的一点点见解，只是自己总结的，没有实践，也没有科学依据，不喜勿喷：

我感觉图像处理能分成三个阶段，或者更多，第一阶段的人很多，听说这行前景好或者工资高的人，多半会学习点图像的知识，比如彩色空间啊，了解下OpenCV啊等等，还有一些属于纯属上课被逼无奈的，比如我们学校就对电子信息类专业和通信类开数字图像处理这门课，而且我当时考试还挂了。。。。这部分大家会接触一些名词，一些简单算法，有用心的同学可能会实现下代码，这阶段风景不错，而且好多名词可以拿来忽悠HR或者忽悠导师，得到一份不错的工作，或者做点导师的项目。

这个阶段多半使用现成的函数库，了解了基本算法或者听别人说一些算法，自己来跑结果，这个阶段Matlab的用户量较大。

第1阶段后半期也就是平台期，这个阶段是做了一段时间有一定算法使用基础和代码能力的工程师，多半在各企业负责最底层的图像算法编写，在他们面前是一座大山，和一个路口，继续做图像还是转管理。

如果继续选择图像，就会面临一个峭壁，具体是算法底层的数学，说的数学，总是困难的，爬这个峭壁的动力可以是挣更多的钱，或者爱好，因为一旦到达阶段2，就能成为首席图像处理工程师，到任何需要图像处理的公司，都能够独当一面，这些人已经到了不缺钱的地步，而且在行业内一定有一定名气。

第二阶段的一旦到达，可以说是事业的平稳期，或者巅峰，不缺钱，还能独自决定一些技术层面上的事，指挥手下阶段1的员工工作，这个阶段面临的也是一个选择，就是靠这个吃一辈子饭，绝对没问题，再有就是向更高的境界冲击。

第三阶段没有尽头，因为能促使进入阶段三的动力只有爱好，这部分挣的钱可能没有阶段2挣的多，而且难度更大，看不到尽头，所以这部分属于探索阶段，在这个阶段上看到的一些，都能推动未来学科的发展，所以这个阶段的人多半是在实验室和数学中忙碌一生，然后几十年后出现在各大论文和教材中。

以上属于个人猜想或者愚见，想法幼稚，仅供参考。

算法描述

今天废话太多了，哈哈，其实上面的可以新开一篇博客单独写出来，但觉得，首先自己年轻视野狭窄，第二水平太低，所以当做废话夹在本文中。

今天说分水岭算法，分水岭算法族是一个思想引出的一些列实现方法。

算法内容：首先将二维灰度图像抽象为三维地形，横纵坐标表示经纬度，灰度表示高度，这样就产生了一个地形图，地形中有高山有盆地，我们的任务就是找到盆地并且给盆地划定地盘。

算法过程，将最低点（灰度）作为起始点，开始从下注水，想泉水一样，水位逐渐上升，所有点都是上下透水的，但不会横向流动，每产生一个新的积水盆地时标记这个盆地使之与原有盆地区别开，当两个盆地内的水要汇集的时候，在此点建立一个水坝，将水分开，此水坝无限高，水永远不会溢出，依次建立所有水坝，最后的水坝就是所谓的分水岭或分水线，也就是得到的分割线。

算法步骤：

  1. 检测图像中所有的盆地，并依次标记，（此时最终的分割区域数已确定）
  2. 将盆地周围邻域的像素按照灰度值大小入队（队列，数据结构，先进先出）
  3. 按照灰度值顺序将像素出队，进行判断

    判断1：如果像素周围只存在一个盆地标记，将此像素划分到此盆地

    判断2：如果像素周围存在两个以及以上盆地，次像素为分水岭。

  4. 如果此像素是盆地像素，将该像素邻域像素入队。
  5. 重复步骤3，直到队列为空。
  6. 得出结果

此算法描述只是分水岭算法之一，冈萨雷斯的书中使用的是形态学的描述方式，但与上述过程原理一致，感兴趣者可自行查阅。

解释下第一步，第一步是本算法的核心，此步骤首先要确定盆地还是非盆地，可以抽象成下面几种模型：

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

代码

//
//  Watershed
//  [email protected]
//  Created by 谭升 on 15/03/03.
//  Copyright (c) 2015年 谭升. All rights reserved.
//
/*
 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;
 };
 */
#include "watershed.h"
#include <stdio.h>

//将非极小值且与极小值相邻的元素入队
void inQueue(PriQueueHead* priQueue,int gray_level,int x,int y){
    priQueue[gray_level].nodeNum++;
    ///malloc new node
    NLevelPriNode* 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 isPan(int *src,int width,int height,double value,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;

}
//由于判断平底锅时使部分数据损坏，现进行恢复
void repairPan(int *src ,int width,int height,double value,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);
                }
            }
}
//平底锅形状，标记极小值为255
void setMinimal(int *src,double *dst,int width,int height,int value,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);
            }
}
//台阶形状，标记非极小值127
void setUnMinimal(int *src,double *dst,int width,int height,int value,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到int
void 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数组，此数组与图像数组多加了label
void 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);
    //for(int j=0;j<height;j++){
    //    for(int i=0;i<width;i++)
    //        printf("  l:%3d|",srclabel[j*width+i].label);
    //    printf("\n");
    //}
    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);

}

实现效果

总结

以上实现目前还有点缺陷，就是盆底面积很大的时候，递归计算会crash，初步设想的解决办法是使用迭代替换递归，分割结果相对来讲对于前面的算法来说相对较好，但运算速度较慢，对噪声敏感。

灰度图像的内容简单的介绍至此，下篇开始介绍彩色基础知识和算法，欢迎收看。

待续。。。

时间： 2024-08-28 03:05:04

灰度图像--图像分割 区域分割之分水岭算法的相关文章

灰度图像--图像分割 区域分割之区域分离

学习DIP第59天 转载请标明本文出处:http://blog.csdn.net/tonyshengtan ,出于尊重文章作者的劳动,转载请标明出处!文章代码已托管,欢迎共同开发:https://github.com/Tony-Tan/DIPpro 开篇废话 废话开始,今天本来只想写一篇,但晚上觉得还是快把区域分割简单介绍下,后面开始彩色图像类的知识学习和代码实现,下一篇介绍分水岭算法,这才是个头疼的算法,今天的区域分离(合并)相对比较好理解. 算法原理 首先本算法依然是基于区域的,用到的区域的

灰度图像--图像分割 区域分割之区域生长

学习DIP第58天 转载请标明本文出处:http://blog.csdn.net/tonyshengtan ,出于尊重文章作者的劳动,转载请标明出处!文章代码已托管,欢迎共同开发:https://github.com/Tony-Tan/DIPpro 开篇废话 继续说废话,昨天写博客被同事看到了,问我,为什么你每一篇开始都是废话,我说凑字数,在一个可以写点轻松的话,天天在算法的海洋里飘荡,偶尔说几句荒山野岭的废话也算活跃气氛了. 区域分割介绍 今天介绍基于区域的分割方法,前面基于阈值的分割方法暂时

新手学，java使用分水岭算法进行图像分割（一）

最近被图像分割整的天昏地暗的,在此感谢老朋友周洋给我关于分水岭算法的指点!本来打算等彩色图像分割有个完满的结果再写这篇文章,但是考虑到到了这一步也算是一个阶段,所以打算对图像分割做一个系列的博文,于是先写这篇. 啰嗦了这么多!先看效果: 效果一般,存在着很多过分割现象,但比没使用滤波之前的效果好很多,过分割是分水岭算法的通病.这个后续博文会继续解决. 本文用java实现的是基于自动种子区域的分水岭算法,注意本文是基于单色的分割,所以将输入图片首先进行灰度化处理,这个比较简单,不多提了:因此,对于

图像分割之分水岭算法

理论 任何灰度图像都可以看作是一个地形表面,其中高强度表示山峰和丘陵,而低强度表示山谷.用不同颜色的水(标签)填充每个孤立的山谷(局部极小值).当水上升时,根据附近的峰(梯度),不同山谷不同的颜色的水,显然会开始融合.为了避免这种情况,你在水就要融合的地方及时增加屏障(增高水坝).你继续填满水,建造屏障,直到所有的山峰都被淹没.然后,您创建的屏障会给出分割结果.这就是分水岭背后的“哲学”.你可以访问分水岭的CMM网页,里面有动画帮助理解. 但是这种方法会由于图像中的噪声或其他不规则性因素而导致过

新手学，java使用分水岭算法进行图像分割（二）

最近要考试了,所以现在不写,怕这段时间都没空写了.(算法的效率暂时不考虑,有时间了再来解决彩色信息和效率问题) 继上一篇的算法:http://blog.csdn.net/abcd_d_/article/details/41218549,本文对分水岭算法进行了区域合并,合并准则采用hsv颜色空间的区域特征的直方图相似度进行合并.且看效果:图一是原图,图二是采用之前的文章算法的效果,图三为进行了区域合并后的效果.(大小被我调整过) (图一) (图二) (图三) 在第一篇的基础之上,增加了区域合并算法

二十一 分水岭算法

一.原理 分水岭算法主要用于图像分段,通常是把一副彩色图像灰度化,然后再求梯度图,最后在梯度图的基础上进行分水岭算法,求得分段图像的边缘线. 任意的灰度图像可以被看做是地质学表面,高亮度的地方是山峰,低亮度的地方是山谷.给每个孤立的山谷(局部最小值)不同颜色的水(标签),当水涨起来,根据周围的山峰(梯度),不同的山谷也就是不同的颜色会开始合并,要避免这个,你可以在水要合并的地方建立障碍,直到所有山峰都被淹没.你所创建的障碍就是分割结果,这个就是分水岭的原理,但是这个方法会分割过度,因为有噪点,或

python数字图像处理（19）：骨架提取与分水岭算法

骨架提取与分水岭算法也属于形态学处理范畴,都放在morphology子模块内. 1.骨架提取 骨架提取,也叫二值图像细化.这种算法能将一个连通区域细化成一个像素的宽度,用于特征提取和目标拓扑表示. morphology子模块提供了两个函数用于骨架提取,分别是Skeletonize()函数和medial_axis()函数.我们先来看Skeletonize()函数. 格式为:skimage.morphology.skeletonize(image) 输入和输出都是一幅二值图像. 例1: from s

opencv分水岭算法对图像进行切割

先看效果 说明 使用分水岭算法对图像进行切割,设置一个标记图像能达到比較好的效果,还能防止过度切割. 1.这里首先对阈值化的二值图像进行腐蚀,去掉小的白色区域,得到图像的前景区域.并对前景区域用255白色标记 2.相同对阈值化后的图像进行膨胀,然后再阈值化并取反.得到背景区域. 并用128灰度表示 3.将前景和背景叠加在一起在同一幅图像中显示. 4.用标记图和原图,利用opencv的watershed对图像进行切割. 源代码 class WatershedSegment{ private: cv

OpenCV 源码中分水岭算法 watershed 函数源码注解

为了研究分水岭算法,阅读了OpenCV 2.4.9 中watershed函数的源码实现部分,代码位于 opencv\sources\modules\imgproc\src\segmentation.cpp 文件中.先贴出加了注解的代码,以后补充对分水岭算法的解释. #include "precomp.hpp" /******************************************************* Watershed **********************