腐蝕、膨脹作為一種簡單、基礎的形態學操作，我之前沒有過多的關注，直到最近發現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