腐蚀膨胀的快速实现

腐蚀、膨胀作为一种简单、基础的形态学操作，我之前没有过多的关注，直到最近发现OpenCV的实现要比自己的实现快几十倍，才进行了深入研究，发现这个操作也并没有想象中的那么简单。

0.准备工作

一般来说，腐蚀和膨胀都是基于二值图像做的，因此我把经典的lena.jpg转换成了二值图像，用于测试效果和性能。代码如下：

    //convert a RGB image to binary
    Mat image=imread("/Users/lena.jpg");
    Mat gray(image.cols,image.rows,CV_8UC1);
    cvtColor(image,gray,CV_BGR2GRAY);
    Mat binary(image.cols,image.rows,CV_8UC1);
    for(int i=0;i<gray.cols*gray.rows*gray.channels();++i)
    {
        if(gray.data[i]>128)
            binary.data[i]=255;
        else
            binary.data[i]=0;
    }
    imshow("lena",binary);
    waitKey(0);

转换后的图片如下：

1.OpenCV函数接口及其实现方式

OpenCV中有dilate和erode两个函数用于膨胀和腐蚀，调用方法如下：

    //opencv dilate&erode functions
    Mat dilateImg(image.cols,image.rows,CV_8UC1);
    Mat erodeImg(image.cols,image.rows,CV_8UC1);
    double dur;
    clock_t start,end;
    
    start = clock();
    dilate(binary,dilateImg,Mat(),cv::Point(-1,-1),5);
    end = clock();
    dur = (double)(end - start);
    printf("OpenCV dilate Use Time:%f ms
",(dur/CLOCKS_PER_SEC*1000));
    
    start = clock();
    erode(binary,erodeImg,Mat(),cv::Point(-1,-1),5);
    end = clock();
    dur = (double)(end - start);
    printf("OpenCV erode Use Time:%f ms
",(dur/CLOCKS_PER_SEC*1000));
    
    imshow("dilate",dilateImg);
    waitKey(0);
    imshow("erode",erodeImg);
    waitKey(0);

在腐蚀膨胀时，均使用了半径／迭代次数＝5。在Release模式下，在我的电脑上运行时间大概在0.35-0.5ms之间。

实际上，OpenCV提供的erode和dilate的实现均在/opencv-2.4.11/modules/imgproc/src/morph.cpp这个文件中，腐蚀和膨胀操作只是针对的像素不同而已，原理是一样的。最终调用的是morphOp这个函数，这个函数又使用

MorphologyRunner来完成真正的操作，而MorphologyRunner又使用了Ptr<FilterEngine>。总之，代码写的错综复杂，但基本可以确定的是，OpenCV是利用多线程和SSE优化来完成的腐蚀膨胀操作。

 

2.自己实现v1

假设做的是膨胀（腐蚀其实同理），考虑最简单的思路，有两种：

1.遍历每个像素，如果该像素值不是0，则其k＊k临域内的所有像素，将被设置成和该像素相同的值（边缘像素需要单独考虑）。

2.遍历每个像素，如果该像素值不是0，或者如果该像素k＊k临域内，有像素值不为0的，则将该像素值设置成非0值。

设图像宽高分别为w，h，膨胀半径／膨胀迭代次数为k。则以上思路时间复杂度为：o(w＊h＊k＊k)。空间复杂度o(w*h)

下面是实现代码，因为膨胀和腐蚀的原理相同，因此后面的代码都是实现膨胀算法。

void myDilateV1(Mat&input,Mat&output,int radius)
{
    for(int i=0;i<input.rows;++i)
    {
        for(int j=0;j<input.cols;++j)
        {
            if(input.data[i*input.cols+j]==255)
            {
                int startX=0>j-radius?0:j-radius;
                int startY=(i-radius)>0?(i-radius):0;
                int endX=input.cols-1<j+radius?input.cols-1:j+radius;
                int endY=(i+radius)<input.rows?(i+radius):(input.rows-1);
                for(int m=startY;m<=endY;++m)
                {
                    for(int n=startX;n<=endX;++n)
                    {
                        output.data[m*input.cols+n]=255;
                    }
                }
            }
        }
    }
}

在我的电脑上运行大概13-15ms之间，比openCV的实现大概慢了30倍。

3.自己实现v2

把每个像素膨胀的区域先设置成另一个值，例如128，遍历完整个图像后，再把这些128值得像素设置成255，这样就可以在原始图像上做膨胀了，代码如下：

void myDilateV2(Mat&input,int radius)
{
    for(int i=0;i<input.rows;++i)
    {
        for(int j=0;j<input.cols;++j)
        {
            if(input.data[i*input.cols+j]==255)
            {
                int startX=0>j-radius?0:j-radius;
                int startY=(i-radius)>0?(i-radius):0;
                int endX=input.cols-1<j+radius?input.cols-1:j+radius;
                int endY=(i+radius)<input.rows?(i+radius):(input.rows-1);
                for(int m=startY;m<=endY;++m)
                {
                    for(int n=startX;n<=endX;++n)
                    {
                        if(input.data[m*input.cols+n]==0)
                            input.data[m*input.cols+n]=128;
                    }
                }
            }
        }
    }
    
    for(int i=0;i<input.rows;++i)
    {
        for(int j=0;j<input.cols;++j)
        {
            if(input.data[i*input.cols+j]==128)
                input.data[i*input.cols+j]=255;
        }
    }
}

时间复杂度没变，运行耗时减慢到了18-20ms，空间复杂度变成了o(1)

4.自己实现v3

换个思路，以k为半径做膨胀，可以转换为连续k次半径为1的膨胀，基于v2版本，修改代码如下：

void myDilateBy1(Mat&input)
{
    for(int i=0;i<input.rows;++i)
    {
        for(int j=0;j<input.cols;++j)
        {
            if(input.data[i*input.cols+j]==255)
            {
                int startX=0>j-1?0:j-1;
                int startY=(i-1)>0?(i-1):0;
                int endX=input.cols-1<j+1?input.cols-1:j+1;
                int endY=(i+1)<input.rows?(i+1):(input.rows-1);
                for(int m=startY;m<=endY;++m)
                {
                    for(int n=startX;n<=endX;++n)
                    {
                        if(input.data[m*input.cols+n]==0)
                            input.data[m*input.cols+n]=128;
                    }
                }
            }
        }
    }
    
    for(int i=0;i<input.rows;++i)
    {
        for(int j=0;j<input.cols;++j)
        {
            if(input.data[i*input.cols+j]==128)
                input.data[i*input.cols+j]=255;
        }
    }
}

void myDilateV3(Mat&input,int radius)
{
    for(int i=0;i<radius;++i)
    {
        myDilateBy1(input);
    }
}

时间复杂度还是没变，运行时间14-16ms

5.自己实现v4

看到这里，大家应该发现，一直按照之前的思路做下去是不行的，必须的换个角度思考问题了。假设这么一种情况：如果我们预先知道每个像素和离它最近的值为255的像素间的距离，那么我们只需要遍历每个像素，判断距离是否大于给定的radius，就可以完成膨胀算法，时间复杂度只有o(w*h)。那么这个预先知道的信息能不能获得呢？

先考虑一维情况，假设有一个长度为n的一维数组，我先从左到右遍历一遍，寻找某个元素和它所有左边值为255像素的最近距离，并保存下来，然后我们再从右到左遍历一遍，寻找某个元素和它所有右边值为255像素的最近距离，如果比第一遍的值小，那么我们就更新。这样两遍下来，上面说的预先知道的信息就得到了。推广到二维情况，代码如下：

void myDilateV4(Mat&src,Mat&dst,int radius)
{
    int i, j;
    unsigned char maxDist=254;
    //step1 forward
    //the top left pixel
    dst.data[0]=(src.data[0]==255)?0:maxDist;
    //the first row
    for (i=1;i<src.cols;++i)
    {dst.data[i]=(src.data[i]==255)?0:MIN(maxDist,dst.data[i-1]+1);}
    //rest pixles
    for (i=1;i<src.rows;++i) // traverse from top left to bottom right
    {
        //left-most pixel
        dst.data[i*src.cols]=(src.data[i*src.cols]==255)?0:MIN(MIN(dst.data[(i-1)*src.cols],dst.data[(i-1)*src.cols+1])+1,maxDist);
        for (j=1;j<src.cols-1;++j)
        {
            if (src.data[i*src.cols+j]==255)
            {
                dst.data[i*src.cols+j]=0;//first pass and pixel was on, it gets a zero
            }
            else
            {
                // pixel was off
                dst.data[i*src.cols+j] = MIN(maxDist, dst.data[(i-1)*src.cols+j]+1);
                dst.data[i*src.cols+j] = MIN(MIN(dst.data[i*src.cols+j-1]+1,dst.data[i*src.cols+j]),MIN(dst.data[(i-1)*src.cols+j-1]+1, dst.data[(i-1)*src.cols+j+1]+1));
            }
        }
        //right-most pixle
        dst.data[(i+1)*src.cols-1]=(src.data[(i+1)*src.cols-1]==255)?0:MIN
        (MIN(dst.data[i*src.cols-1]+1,dst.data[i*src.cols-2]+1),
         MIN(dst.data[(i+1)*src.cols-2]+1,maxDist));
    }
    //step2 backward
    //    the bottom right pixel
    dst.data[src.cols*src.rows-1]=(src.data[src.cols*src.rows-1]==255)?0:dst.data[src.cols*src.rows-1];
    //the last row
    for (i=src.cols-2;i>=0;--i)
    {dst.data[(src.rows-1)*src.cols+i]=(src.data[(src.rows-1)*src.cols+i]==255)?0:MIN(dst.data[(src.rows-1)*src.cols+i],dst.data[(src.rows-1)*src.cols+i+1]+1);}
    //rest pixles
    for (i=src.rows-2; i>=0; --i){
        //right-most pixel
        dst.data[(i+1)*src.cols-1]=(src.data[(i+1)*src.cols-1]==255)?0:MIN(MIN(dst.data[(i+2)*src.cols-1],dst.data[(i+2)*src.cols-2])+1,dst.data[(i+1)*src.cols-1]);
        for (j=src.cols-2; j>0; --j)
        {
            if(src.data[i*src.cols+j]!=255)
            {
                dst.data[i*src.cols+j] = MIN(dst.data[i*src.cols+j], dst.data[i*src.cols+j+1]+1);
                dst.data[i*src.cols+j] = MIN(
                                     MIN(dst.data[i*src.cols+j], dst.data[(i+1)*src.cols+j]+1),
                                     MIN(dst.data[(i+1)*src.cols+j+1]+1, dst.data[(i+1)*src.cols+j-1]+1)
                                     );
            }
        }
        //left-most pixel
        dst.data[i*src.cols]=(src.data[i*src.cols]==255)?0:MIN
        (MIN(dst.data[i*src.cols+1],dst.data[(i+1)*src.cols])+1,
         MIN(dst.data[(i+1)*src.cols+1]+1,dst.data[i*src.cols]));
    }
    for (i=0;i<src.rows;++i)
    {
        for (j=0;j<src.cols;++j)
        {
            dst.data[i*src.cols+j] = ((dst.data[i*src.cols+j]<=radius)?255:src.data[i*src.cols+j]);
        }
    }
}

在我的电脑上耗时大概2-3ms，空间复杂度o(w*h)。虽然还是比opencv慢了5-6倍，但比之前几个版本的实现快了5倍以上。

参考资料：

http://blog.ostermiller.org/dilate-and-erode