腐蚀膨胀的快速实现
腐蚀、膨胀作为一种简单、基础的形态学操作，我之前没有过多的关注，直到最近发现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\n",(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\n",(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
