【数字图像处理】3.2:二值图像-形态学处理腐蚀和膨胀

Abstract: 数字图像处理：第10天
Keywords: 形态学，腐蚀，膨胀，边缘

本文最初发表于csdn，于2018年2月17日迁移至此

开篇废话

今天来介绍形态学中最基础也是最重要的两个操作，腐蚀和膨胀，腐蚀和膨胀基本上是所有形态学操作的基础，除此之外还有补集（即二值图全部取反的操作，0变1，1变0），和反射（将所有坐标去反）。

之前使用过腐蚀和膨胀，仅仅是去噪，那时候连结构元（SE）的概念都没有，其实所有形态学操作，核心都是结构元（包括其形状和中心位置，中心位置可以不在结构元区域中），他的变化可以产生千奇百怪的效果，如果你能很好的设计结构元，那么你将能得到你想要的效果。

对于写博客，我觉得坚持下来还是不错的，一是可以反思总结一下学习结果，很多收获都是写博客的时候想到的，虽然写起来浪费很多时间，包括作图之类的工作。二是可以作为资料，以后查看，查漏补缺。三是与别人分享知识，如果有问题，可以有人及时指正，良师益友。

应用

腐蚀和膨胀的应用应该很多，因为其操作简单，属于基础运算，这里简单列举下用途，还是那句话，SE是关键，设计好了可以处理很多二值图像的问题。

膨胀：

桥接缝隙（缝隙点为0，且宽度比SE的宽度小）
消除细小的黑点（二值图像中的0，黑点比SE小）

腐蚀：

消除“桥梁”（细线装的白色条纹，值为1，宽度小于SE的宽度）
消除细小的白点（二值图像中的1，白点比SE小）

此外，经过组合，腐蚀和膨胀将完成基本全部的形态学操作，例如后面介绍的，开操作，闭操作，命中与不命中，提取边缘，提取骨骼，裁剪等等。

数学基础

数学形态学的数学基础是集合论，Milan Sonka etc.的《图像处理、分析与机器视觉》（以下简称IPAMV）中提到的一篇论文【Serra，1982】（此文在结尾附下载地址，版权归IEEE所有，请勿用于商业用途），文中给出了数学形态学的基本操作，和定义，如下：

Center

上图给出了最基本的腐蚀和膨胀操作，以及平移操作，由于上图中，英文较为简单，这里不再过多的解释，基本操作都是集合操作，如，交并补集操作。
首先，我在看这篇文章之前，对于腐蚀膨胀一直停留在SE划过窗口的模式，也就是，SE在图像上扫描，满足某些条件时进行某些操作，但这么做的问题在于，如果SE比窗口大，按照上述将无法操作，而且，我们无法验证膨胀的交换不变性，A膨胀B=B膨胀A，因为除非A，B大小相等，否者大的无法在小的上滑动。
而文中的标准定义是，对图像X进行移动，包括 $b_1,b_2,b_3,b_4\dots b_n$ 即集合B中的所有移动的结构的并集，换句话说，就是集合A使用SE,B集合在膨胀，等于B的子集分别腐蚀的并集（子集的并集等于B）。
上图，做了好久的图：
STEP1：首先观察左侧SE，红色为中心，和明显，SE具有各向同性，这也是一个结构元的重要特点，即各向同性与各向异性有不同的效果，也有不同的应用。

Center 1

STEP2：SE显示，集合X需要向左移1个单位。

Center 2

STEP3：SE显示，集合X需要向右移1个单位。

Center 3

STEP4：SE显示，集合X需要向上移1个单位。

Center 4

STEP5：SE显示，集合X需要向下移1个单位。

Center 5

STEP last：将所有结果取并集，两种颜色的表示两种操作都能产生那个单元。黑色为原始集合。

Center 6

其实这种解释方法比SE滑过窗口的解释更加严谨，也更容易理解，其实SE滑过窗口的原理与着相同，只是用了类似分治的思想。但我实现起来好像两种方法的算法复杂度差不多。

说道复杂度，由于腐蚀和膨胀属于集合操作，所以，不属于线性操作。

腐蚀不是膨胀的逆操作，虽然有时可以经过腐蚀后膨胀来恢复原图，但他们并不是一对互逆的操作。其具有对偶性，后续介绍（因为数学公式不好输入）

腐蚀的具体操作与膨胀类似，但有以下不同：
首先，移动方向，与膨胀不同，如果SE为：
$$
\begin{bmatrix}
0&1&0\\
1&1&1\\
0&1&0
\end{bmatrix}
$$
红色为SE原点，那么绿色代表的位移不是向右移动而是向左移动，即为反方向。
上面数学基础中给出的腐蚀公式并不准确，IPAMV中给出了准确的公式，位移b前应有负号，即-b。
其次，集合操作是交集，这个很关键。

腐蚀和膨胀的性质

这一节先空着，因为公式不好输入，后续会填上。。敬请期待。（广告：本人找工作，211本科，14年毕业，无工作经验，爱好图像处理，有人要的话随时牵走，工作地点限深圳，联系方式：下面留言）。

代码

上代码：

#include <cv.h>
#include <highgui.h>
#include <stdio.h>
#define isSIZEEQU(x,y) (((x)->width)==((y)->width)&&((x)->height)==((y)->height))
typedef int DataType;
struct Position_{
    int x;
    int y;
};
typedef struct Position_ Position;
typedef struct Position_ MoveDirection;
//位移操作，将图像整体移动，如果超出边界舍去
void Translation(IplImage *src,IplImage *dst,MoveDirection *direction){
    int width=src->width;
    int height=src->height;
    //printf("%d,%d\n",direction->x,direction->y);
    IplImage *temp=cvCreateImage(cvSize(width, height), src->depth, src->nChannels);
    cvZero(temp);
    for(int i=0;i<width;i++)
        for(int j=0;j<height;j++){
                if(j+direction->y<height &&
                   i+direction->x<width  &&
                   j+direction->y>=0      &&
                   i+direction->x>=0        )
                cvSetReal2D(temp, j+direction->y, i+direction->x, cvGetReal2D(src, j, i));

        }
    cvCopy(temp, dst, NULL);
    cvReleaseImage(&temp);
}
//将小的图像弄到大的黑色图像中间，或者说是给图像加黑色边框
void Zoom(IplImage *src,IplImage *dst){
    if(dst->width<src->width         ||
       dst->height<src->height       ||
       (dst->height-src->height)%2==1||
       (dst->width-src->width)%2==1){
        if(dst->width<src->width )
            printf("Zoom wrong:dst's width too small!\n");
        if(dst->height<src->height )
            printf("Zoom wrong:dst's height too small!\n");
        if((dst->height-src->height)%2==1||(dst->width-src->width)%2==1)
            printf("Zoom wrong:dst-src not a oushu!\n");
        exit(0);
    }
    MoveDirection m;
    m.x=(dst->width-src->width)/2;
    m.y=(dst->height-src->height)/2;
    cvZero(dst);
    for(int i=m.x,j=0;j<src->width;i++,j++){
        for(int k=m.y,n=0;n<src->height;k++,n++){
            cvSetReal2D(dst, k, i, cvGetReal2D(src, n, j));

        }
    }

}
//逻辑与操作
void And(IplImage *src0,IplImage *src1,IplImage *dst){
    if(!isSIZEEQU(src0,src1)){
        printf("And wrong !\n");
        exit(0);
    }
    if(!isSIZEEQU(src0,dst)){
        printf("And wrong !\n");
        exit(0);
    }
    int width=src0->width;
    int height=src0->height;
    for(int i=0;i<width;i++){
        for(int j=0;j<height;j++){
            if(cvGetReal2D(src0, j, i)>100.0&&
               cvGetReal2D(src1, j, i)>100.0)
                cvSetReal2D(dst, j, i, 255.0);
            else
                cvSetReal2D(dst, j, i, 0.0);
        }
    }
}
//逻辑或操作
void Or(IplImage *src0,IplImage *src1,IplImage *dst){
    if(!isSIZEEQU(src0,src1)){
        printf("And wrong !\n");
        exit(0);
    }
    if(!isSIZEEQU(src0,dst)){
        printf("And wrong !\n");
        exit(0);
    }
    int width=src0->width;
    int height=src0->height;
    for(int i=0;i<width;i++){
        for(int j=0;j<height;j++){
            if(cvGetReal2D(src0, j, i)>100.0||
               cvGetReal2D(src1, j, i)>100.0)
                cvSetReal2D(dst, j, i, 255);

        }
    }
}
//将所有元素设为1
void One(IplImage *src){
    for(int i=0;i<src->width;i++)
        for(int j=0;j<src->height;j++)
            cvSetReal2D(src, j, i, 255.0);


}
//膨胀
void Dilate(IplImage *src,IplImage *dst,IplImage *se,Position *center){
    if(center==NULL){
        Position temp;
        temp.x=se->width/2;
        temp.y=se->height/2;
        center=&temp;
    }
    //printf("%d,%d",center->x,center->y);
    MoveDirection m;
    IplImage *temp=cvCreateImage(cvGetSize(dst), dst->depth,dst->nChannels);
    IplImage *tempdst=cvCreateImage(cvGetSize(dst), dst->depth,dst->nChannels);
    IplImage *realdst=cvCreateImage(cvGetSize(dst), dst->depth,dst->nChannels);
    cvZero(realdst);
    Zoom(src,temp);
    int width=se->width;
    int height=se->height;
    for(int i=0;i<width;i++){
        for(int j=0;j<height;j++){
            if(cvGetReal2D(se, j, i)>100.0){
                m.x=i-center->x;
                m.y=j-center->y;
                Translation(temp,tempdst, &m);
                Or(tempdst, realdst, realdst);
            }
        }
    }
    cvCopy(realdst, dst, NULL);
    cvReleaseImage(&temp);
    cvReleaseImage(&realdst);
    cvReleaseImage(&tempdst);
}
//腐蚀
void Erode(IplImage *src,IplImage *dst,IplImage *se,Position *center){
    if(center==NULL){
        Position temp;
        temp.x=se->width/2;
        temp.y=se->height/2;
        center=&temp;
    }
    MoveDirection m;
    IplImage *temp=cvCreateImage(cvGetSize(dst), dst->depth,dst->nChannels);
    IplImage *tempdst=cvCreateImage(cvGetSize(dst), dst->depth,dst->nChannels);
    IplImage *realdst=cvCreateImage(cvGetSize(dst), dst->depth,dst->nChannels);
    One(realdst);
    Zoom(src,temp);
    int width=se->width;
    int height=se->height;
    for(int i=0;i<width;i++){
        for(int j=0;j<height;j++){
            if(cvGetReal2D(se, j, i)>100.0){
                m.x=center->x-i;
                m.y=center->y-j;
                Translation(temp,tempdst, &m);
                And(tempdst, realdst, realdst);
            }
        }
    }
    cvCopy(realdst, dst, NULL);
    cvReleaseImage(&tempdst);
    cvReleaseImage(&temp);
    cvReleaseImage(&realdst);
}

//开操作
void Open(IplImage *src,IplImage *dst,IplImage *se,Position *center){
    Erode(src, dst, se, center);
    Dilate(dst, dst, se, center);

}
//关操作
void Close(IplImage *src,IplImage *dst,IplImage *se,Position *center){
    Dilate(src, dst, se, center);
    Erode(dst, dst, se, center);

}

int main(){
    IplImage *se=cvLoadImage("/Users/Tony/Binary_Image/mask6.jpg",0);

    IplImage *src=cvLoadImage("/Users/Tony/lena/lena_BW.jpg", 0);
    IplImage *dst=cvCreateImage(cvGetSize(src), 8, 1);
    IplImage *subdst=cvCreateImage(cvGetSize(src), 8, 1);



    Close(src, dst, se, NULL);
    cvSub(dst,src,subdst, NULL);
    cvNamedWindow("SRC", 1);
    cvShowImage("SRC", src);
    cvNamedWindow("DST", 1);
    cvShowImage("DST", dst);
    cvNamedWindow("SUB", 1);
    cvShowImage("SUB", subdst);
    cvSaveImage("/Users/Tony/Binary_Image/lena_close_sub.jpg", subdst, 0);
    cvSaveImage("/Users/Tony/Binary_Image/lena_close.jpg", dst, 0);
    cvWaitKey(0);
    return 0;
}

结果

膨胀：第一行SE各向同性，后两行SE各向异性
膨胀结果,结构元,与原图的差
Center 7

腐蚀：第一行SE各向同性，后两行SE各向异性
腐蚀结果,结构元,原图的差
Center 8

同一结构元不同中心位置的不同结果：
Center 9
Center 10
Center 11
结构元为简单，各向同性：
$$
\begin{bmatrix}
0&1&0\\
1&1&1\\
0&1&0
\end{bmatrix}
$$
实际中的应用，lena图，灰度图100为阈值后的二值图：
原图
Center 12

腐蚀与原图的差,膨胀与原图的差
Center 13