【数字图像处理】3.4:二值图像-形态学处理4 其他操作

Abstract: 数字图像处理：第12天
Keywords: 二值图像，形态学，击中，边界提取，孔洞填充，连通分量提取，凸壳，细化，骨架，形态学重建

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

开篇废话

其实写博客是个很痛苦的过程，要准备一些东西，还怕写错会误导别人，但是在总结和准备相关资料的时候会更深入的理解其中的一些算法原理，然后再根据自己的实际操作把过程表述出来，可以发现一些错误，也能理顺思路。我百度找资料的时候（大牛都用google你怎么还用百度？我找不到上Google的方法![）发现，我写的前几篇博文被一个论坛转载了，但图片有些只能显示一半，也没有标明出处，对此我表示很蛋疼，首先，写出来的东西就是给别人看的，我很希望自己的博客被转载，但是，出于礼貌，应该写上原文出处，其次应该尽量把文章保持原状，不要误导别人。
开始正题，形态学基本操作：腐蚀和膨胀，进而演化出了开操作和闭操作，后来有演化出了一些其他操作，能够应用于各种不同的场景，比如最简单的边界提取，稍微复杂的孔洞填充，与空洞填充类似的连通分量提取，凸壳不知道在哪方面应用，但是细化和骨架在很多手势识别，步态识别中应用广泛，形态学重建可以算作一种很简单的图像分割，但高速有效。

数学公式

击中：
边界提取：
孔洞填充：
连通分量提取：
凸壳：
细化：
骨架：
测地腐蚀：
测地膨胀：
（2018年补充：详情可以参考刚萨雷斯数字图像处理，或等作者后续补上）

算法介绍

击中：击中其实是一种模板匹配，如果一个结构元（模板）与图像中的一个连通区域完全相等，那么腐蚀的结果将是一个点，这就算是“击中”了，因为找到了完全一致的模式（模式的含义去查一下吧，这里的模式和模式识别的模式含义相同），但如果模板与联通区域去不完全相同，但腐蚀后的结果还是一个点，我们也称为一种“击中”，但不是严格的击中，比如我们的模板是一个直径为10的实心圆形，连通区域是一个直径为10的实心圆和一个直径为1的实心圆的并列排布，其腐蚀结果也是一个点，但其不是完美的击中，要找到完美的击中，必须在模板外面加个边框（或者是背景），这样的结果就不会出现不完美击中，我们可以选模板外一圈黑色的背景（红色背景）作为边框，与上述连通区域做击中操作，很明显。。击不中。。完美击中和不完美的击中都有很多用途，后面的操作中会继续用到。
上图，看图更明确：

Center

红色表示背景，击中结果是只能有一个像素被点亮

Center 1

一般性的匹配，不关心背景，只寻找一个模式

孔洞填充：最简单的理解，一个白色的圆环，背景为黑，想要的结果是一个白色的实心圆，这就用到了孔洞填充，孔洞填充的一个缺点在于，必须已知一个种子点，从种子点向外膨胀至整个孔洞。
下图中红色为种子点，黑色为背景
种子点和原图

填充后结果
连通分量提取：与孔洞填充类似，孔洞处理的是被连通分量包围的背景，连通分量提取是提取的连通分量，而非中间的孔洞。所以从公式上看也非常相似。连通分量提取也需要种子点，所以也不是很智能。
四个连通分量，红点为示例种子点
凸壳：是为了找到一个凹陷的物体外壳，以不完美命中（腐蚀）为主要操作，通过调整结构元，加上原图，得到相关结果。
原图

凸壳添加结果

添加的凸壳
细化：细化和骨架有点类似，但涉及到一个是否同伦操作，来贴一下同伦的含义，拓扑学的一个概念，是连通性的一个概念定义，如果变换不改变连通结构，则视为同伦，或者同伦树不变，以下来自百度百科对同伦的解释，wiki的数学有点复杂，不好理解：

Center 8

这个解释的很形象，尤其是终结者的实例，很形象。
上图：
Center 9
原图放大好多倍以后，每个方框是一个像素
结果
Center 10

骨架：使用形态学腐蚀减去结果的开运算的骨架与细化相比，缺少的就是同伦性，即骨架操作得到的并不是原图像的同伦变换，而且这种骨架有些地方并不是一个像素，而是多个像素，使得这种骨架算法应用不是很广泛，但骨架却应用相当广泛
形态学重建：提取原始图像中包含某些特定特征的连通区域，需要一个模板和一些种子点，比如重建开操作提取文字中有长竖“I”这种结构的，先用长竖对图像进行几次腐蚀（命中），使之得到一些点，在膨胀这些点，并以原图作为模板，最终得到包含长竖的字母。

代码


#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;
//位移操作，将图像整体移动，如果超出边界舍去
int isEqual(IplImage *src1,IplImage *src2){
    if(!isSIZEEQU(src1, src2))
        return 0;
    int width=src1->width;
    int height=src1->height;
    for(int i=0;i<width;i++)
        for(int j=0;j<height;j++){
            int v0=cvGetReal2D(src1, j, i);
            int v1=cvGetReal2D(src2, j, i);
            if(v0!=v1)
                return 0;
        }
    return 1;
}
//检测图像是否为空
int isEmpty(IplImage *src){
    int width=src->width;
    int height=src->height;
    for(int i=0;i<width;i++)
        for(int j=0;j<height;j++){
            int v=cvGetReal2D(src, j, i);
            if(v!=0.0)
                return 0;
        }
    return 1;
}

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));

        }
    }

}
//逻辑与操作
int And(IplImage *src0,IplImage *src1,IplImage *dst){
    int isChanged=0;
    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);
        }
    }
    return isChanged;
}
//逻辑或操作
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);

        }
    }
}
//取反
void Not(IplImage *src,IplImage *dst){

    if(!isSIZEEQU(src,dst)){
        printf("Not wrong !\n");
        exit(0);
    }
    int width=src->width;
    int height=src->height;
    for(int i=0;i<width;i++){
        for(int j=0;j<height;j++){
            cvSetReal2D(dst, j, i, 255.0-cvGetReal2D(src, j, i));
        }
    }
}
//将所有元素设为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);

}
//击中与击不中
void HitorMiss(IplImage *src,IplImage *se1,IplImage *se2,IplImage *dst,Position *se1center,Position *se2center){
    IplImage *temp1=cvCreateImage(cvGetSize(src), src->depth, src->nChannels);
    IplImage *temp2=cvCreateImage(cvGetSize(src), src->depth, src->nChannels);
    Erode(src, temp1, se1, se1center);
    Not(src, temp2);
    Erode(temp2, temp2, se2, se2center);
    And(temp1, temp2, dst);
    cvReleaseImage(&temp1);
    cvReleaseImage(&temp2);
}
//二值图像，边缘检测
void BinaryEdge(IplImage *src,IplImage* dst){
    IplImage *temp=cvCreateImage(cvGetSize(src), src->depth, src->nChannels);
    Erode(src, temp, NULL, NULL);
    cvSub(src, temp, dst, NULL);
    cvReleaseImage(&temp);
}
//孔洞填充
void FillHole(IplImage *src,IplImage *dst,IplImage *se,Position *seed){
    IplImage * temp=cvCreateImage(cvGetSize(src), src->depth, src->nChannels);
    cvZero(temp);
    IplImage * lasttemp=cvCreateImage(cvGetSize(src), src->depth, src->nChannels);
    IplImage * nsrc=cvCreateImage(cvGetSize(src), src->depth, src->nChannels);
    Not(src, nsrc);
    cvSetReal2D(temp, seed->y, seed->x, 255.0);
    while(!isEqual(lasttemp, temp)){
        cvCopy(temp, lasttemp, NULL);
        Dilate(temp, temp, se, NULL);
        And(temp, nsrc, temp);



    }
    Or(temp, src, dst);
    cvReleaseImage(&temp);
    cvReleaseImage(&lasttemp);
    cvReleaseImage(&nsrc);
}
//连通分量获取
void GetConComponent(IplImage *src,IplImage *dst,IplImage *se,Position *seed){
    IplImage * temp=cvCreateImage(cvGetSize(src), src->depth, src->nChannels);
    cvZero(temp);
    IplImage * lasttemp=cvCreateImage(cvGetSize(src), src->depth, src->nChannels);

    cvSetReal2D(temp, seed->y, seed->x, 255.0);
    while(!isEqual(lasttemp, temp)){
        cvCopy(temp, lasttemp, NULL);
        Dilate(temp, temp, se, NULL);
        And(temp, src, temp);

    }
    cvCopy(temp, dst, NULL);
    cvReleaseImage(&temp);
    cvReleaseImage(&lasttemp);
}
//骨架
void FrameWork(IplImage *src,IplImage *dst,IplImage *se){
    cvZero(dst);
    IplImage *temp=cvCreateImage(cvGetSize(src), src->depth,src->nChannels);
    IplImage *temp_open=cvCreateImage(cvGetSize(src), src->depth,src->nChannels);
    cvCopy(src, temp, NULL);
    while(!isEmpty(temp)){
        Erode(temp, temp, se, NULL);
        cvCopy(temp, temp_open, NULL);
        Open(temp_open, temp_open, se, NULL);
        cvSub(temp, temp_open, temp_open,NULL);
        Or(temp_open, dst, dst);

    }
    cvReleaseImage(&temp);
    cvReleaseImage(&temp_open);
}
//凸壳生成结构元
IplImage* CreateConvexhullSE(int num){
    IplImage *se=cvCreateImage(cvSize(3, 3), 8, 1);
    cvZero(se);
    switch (num) {
        case 0:
        {
            cvSetReal2D(se, 0, 0, 255.0);
            cvSetReal2D(se, 1, 0, 255.0);
            cvSetReal2D(se, 2, 0, 255.0);
        }
            break;
        case 1:
        {
            cvSetReal2D(se, 0, 0, 255.0);
            cvSetReal2D(se, 0, 1, 255.0);
            cvSetReal2D(se, 0, 2, 255.0);

        }
            break;
        case 2:
        {
            cvSetReal2D(se, 0, 2, 255.0);
            cvSetReal2D(se, 1, 2, 255.0);
            cvSetReal2D(se, 2, 2, 255.0);
        }
            break;
        case 3:
        {
            cvSetReal2D(se, 2, 0, 255.0);
            cvSetReal2D(se, 2, 1, 255.0);
            cvSetReal2D(se, 2, 2, 255.0);
        }
            break;
        default:
            break;
    }
    return se;
}
//凸壳
void Convexhull(IplImage *src,IplImage *dst){
    cvCopy(src, dst, NULL);
    IplImage * se[4];
    IplImage *temp=cvCreateImage(cvGetSize(src), src->depth, src->nChannels);
    IplImage *temp_last=cvCreateImage(cvGetSize(src), src->depth, src->nChannels);
    //cvCopy(src, temp, NULL);
    for(int i=0;i<4;i++){
        cvCopy(src, temp, NULL);
        se[i]=CreateConvexhullSE(i);
        while (!isEqual(temp, temp_last)) {
            cvCopy(temp, temp_last, NULL);
            Erode(temp, temp, se[i],NULL);
            Or(temp, dst, temp);


        }
        cvCopy(temp, dst, NULL);
        cvReleaseImage(&se[i]);
    }
    cvReleaseImage(&temp);
    cvReleaseImage(&temp_last);

}
//生成细化结构元
IplImage* CreateThinningSE(int num){
    IplImage *se=cvCreateImage(cvSize(3, 3), 8, 1);
    cvZero(se);
    switch (num) {
        case 0:
        {
            cvSetReal2D(se, 2, 0, 255.0);
            cvSetReal2D(se, 2, 1, 255.0);
            cvSetReal2D(se, 2, 2, 255.0);
            cvSetReal2D(se, 1, 1, 255.0);
        }
            break;
        case 1:
        {
            cvSetReal2D(se, 1, 1, 255.0);
            cvSetReal2D(se, 1, 0, 255.0);
            cvSetReal2D(se, 2, 0, 255.0);
            cvSetReal2D(se, 2, 1, 255.0);

        }
            break;
        case 2:
        {
            cvSetReal2D(se, 1, 1, 255.0);
            cvSetReal2D(se, 1, 0, 255.0);
            cvSetReal2D(se, 2, 0, 255.0);
            cvSetReal2D(se, 0, 0, 255.0);
        }
            break;
        case 3:
        {
            cvSetReal2D(se, 1, 1, 255.0);
            cvSetReal2D(se, 0, 0, 255.0);
            cvSetReal2D(se, 0, 1, 255.0);
            cvSetReal2D(se, 1, 0, 255.0);
        }
            break;
        case 4:
        {
            cvSetReal2D(se, 1, 1, 255.0);
            cvSetReal2D(se, 0, 0, 255.0);
            cvSetReal2D(se, 0, 1, 255.0);
            cvSetReal2D(se, 0, 2, 255.0);
        }
            break;
        case 5:
        {
            cvSetReal2D(se, 0, 1, 255.0);
            cvSetReal2D(se, 0, 2, 255.0);
            cvSetReal2D(se, 1, 1, 255.0);
            cvSetReal2D(se, 1, 2, 255.0);

        }
            break;
        case 6:
        {
            cvSetReal2D(se, 1, 1, 255.0);
            cvSetReal2D(se, 0, 2, 255.0);
            cvSetReal2D(se, 1, 2, 255.0);
            cvSetReal2D(se, 2, 2, 255.0);
        }
            break;
        case 7:
        {
            cvSetReal2D(se, 1, 1, 255.0);
            cvSetReal2D(se, 1, 2, 255.0);
            cvSetReal2D(se, 2, 1, 255.0);
            cvSetReal2D(se, 2, 2, 255.0);
        }
            break;
        default:
            break;
    }
    return se;
}
IplImage* CreateThinningUSE(int num){
    IplImage *se=cvCreateImage(cvSize(3, 3), 8, 1);
    cvZero(se);
    switch (num) {
        case 0:
        {

            cvSetReal2D(se, 0, 1, 255.0);
            cvSetReal2D(se, 0, 2, 255.0);
            cvSetReal2D(se, 0, 0, 255.0);
        }
            break;
        case 1:
        {
            cvSetReal2D(se, 0, 1, 255.0);
            cvSetReal2D(se, 0, 2, 255.0);
            cvSetReal2D(se, 1, 2, 255.0);

        }
            break;
        case 2:
        {

            cvSetReal2D(se, 0, 2, 255.0);
            cvSetReal2D(se, 1, 2, 255.0);
            cvSetReal2D(se, 2, 2, 255.0);
        }
            break;
        case 3:
        {
            cvSetReal2D(se, 1, 2, 255.0);
            cvSetReal2D(se, 2, 1, 255.0);
            cvSetReal2D(se, 2, 2, 255.0);
        }
            break;
        case 4:
        {

            cvSetReal2D(se, 2, 0, 255.0);
            cvSetReal2D(se, 2, 1, 255.0);
            cvSetReal2D(se, 2, 2, 255.0);
        }
            break;
        case 5:
        {

            cvSetReal2D(se, 1, 0, 255.0);
            cvSetReal2D(se, 2, 0, 255.0);
            cvSetReal2D(se, 2, 1, 255.0);

        }
            break;
        case 6:
        {

            cvSetReal2D(se, 0, 0, 255.0);
            cvSetReal2D(se, 1, 0, 255.0);
            cvSetReal2D(se, 2, 0, 255.0);
        }
            break;
        case 7:
        {
            cvSetReal2D(se, 0, 0, 255.0);
            cvSetReal2D(se, 0, 1, 255.0);
            cvSetReal2D(se, 1, 0, 255.0);
        }
            break;
        default:
            break;
    }
    return se;
}

//细化操作
void Thinning(IplImage *src,IplImage *dst){
    IplImage *temp=cvCreateImage(cvGetSize(src), src->depth, src->nChannels);
    IplImage *temp_last=cvCreateImage(cvGetSize(src), src->depth, src->nChannels);
    IplImage *temp_com=cvCreateImage(cvGetSize(src), src->depth, src->nChannels);
    cvZero(temp_last);
    cvCopy(src, temp, NULL);

    while(!isEqual(temp, temp_com)){
        cvCopy(temp, temp_com, NULL);
        for(int i=0;i<8;i++){
            cvCopy(temp, temp_last, NULL);
            IplImage *se1=CreateThinningSE(i);
            IplImage *se2=CreateThinningUSE(i);
            HitorMiss(temp, se1, se2, temp, NULL, NULL);
            cvSub(temp_last, temp, temp, NULL);
            cvReleaseImage(&se1);
            cvReleaseImage(&se2);
        }

    }
    cvCopy(temp, dst, NULL);
    cvReleaseImage(&temp);
    cvReleaseImage(&temp_com);
    cvReleaseImage(&temp_last);
}
//重建开操作
void reBuildOpen(IplImage *src,IplImage *dst,IplImage *ground,IplImage *dilateSE,IplImage *erodeSE,int eroden){
    IplImage *temp=cvCreateImage(cvGetSize(src), src->depth, src->nChannels);
    IplImage *temp_last=cvCreateImage(cvGetSize(src), src->depth, src->nChannels);
    cvCopy(src, temp, NULL);
    for(int i=0;i<eroden;i++){
        Erode(temp, temp, erodeSE, NULL);
    }

    while(!isEqual(temp, temp_last)){
        cvCopy(temp, temp_last, NULL);
        Dilate(temp, temp, dilateSE, NULL);
        And(temp, ground, temp);

    }
    cvCopy(temp, dst, NULL);
    cvReleaseImage(&temp);
    cvReleaseImage(&temp_last);

}

int main(){

    return 0;
}