图像处理 腐蚀算法

2015-12-30 16:48:37 minushuang 阅读数 12031

原理:在特殊领域运算形式——结构元素(Sturcture Element),在每个像素位置上与二值图像对应的区域进行特定的逻辑运算。运算结构是输出图像的相应像素。运算效果取决于结构元素大小内容以及逻辑运算性质。

结构元素:膨胀和腐蚀操作的最基本组成部分,用于测试输出图像,通常要比待处理的图像小还很多。二维平面结构元素由一个数值为0或1的矩阵组成。结构元素的原点指定了图像中需要处理的像素范围,结构元素中数值为1的点决定结构元素的邻域像素在进行膨胀或腐蚀操作时是否需要参与计算。

先来定义一些基本符号和关系。

1.         元素

设有一幅图象X,若点aX的区域以内,则称aX的元素,记作aX,如图6.1所示。

2.         B包含于X

设有两幅图象BX。对于B中所有的元素ai,都有aiX,则称B包含于(included in)X,记作B  X,如图6.2所示。

3.         B击中X

设有两幅图象BX。若存在这样一个点,它即是B的元素,又是X的元素,则称B击中(hit)X,记作BX,如图6.3所示。

4.         B不击中X

设有两幅图象BX。若不存在任何一个点,它即是B的元素,又是X的元素,即BX的交集是空,则称B不击中(miss)X,记作BX=Ф;其中∩是集合运算相交的符号,Ф表示空集。如图6.4所示。

6.1     元素

6.2     包含

6.3     击中

6.4     不击中

5.         补集

设有一幅图象X,所有X区域以外的点构成的集合称为X的补集,记作Xc,如图6.5所示。显然,如果BX=Ф,则BX的补集内,即B  Xc

6.5     补集的示意图

6.         结构元素

设有两幅图象BX。若X是被处理的对象,而B是用来处理X的,则称B为结构元素(structure element),又被形象地称做刷子。结构元素通常都是一些比较小的图象。

7.         对称集

设有一幅图象B,将B中所有元素的坐标取反,即令(xy)变成(-x-y),所有这些点构成的新的集合称为B的对称集,记作Bv,如图6.6所示。

8.         平移

设有一幅图象B,有一个点a(x0,y0),将B平移a后的结果是,把B中所有元素的横坐标加x0,纵坐标加y0,即令(xy)变成(x+x0y+y0),所有这些点构成的新的集合称为B的平移,记作Ba,如图6.7所示。

6.6     对称集的示意图

6.7     平移的示意图

好了,介绍了这么多基本符号和关系,现在让我们应用这些符号和关系,看一下形态学的基本运算。

6.1 腐蚀

把结构元素B平移a后得到Ba,若Ba包含于X,我们记下这个a点,所有满足上述条件的a点组成的集合称做XB腐蚀(Erosion)的结果。用公式表示为:E(X)={a| Ba  X}=X  B,如图6.8所示。

6.8     腐蚀的示意图

6.8X是被处理的对象,B是结构元素。不难知道,对于任意一个在阴影部分的点aBa 包含于X,所以XB腐蚀的结果就是那个阴影部分。阴影部分在X的范围之内,且比X小,就象X被剥掉了一层似的,这就是为什么叫腐蚀的原因。

值得注意的是,上面的B是对称的,即B的对称集Bv=B,所以XB腐蚀的结果和X Bv腐蚀的结果是一样的。如果B不是对称的,让我们看看图6.9,就会发现XB腐蚀的结果和X Bv腐蚀的结果不同。

6.9     结构元素非对称时,腐蚀的结果不同

6.8和图6.9都是示意图,让我们来看看实际上是怎样进行腐蚀运算的。

在图6.10中,左边是被处理的图象X(二值图象,我们针对的是黑点),中间是结构元素B,那个标有origin的点是中心点,即当前处理元素的位置,我们在介绍模板操作时也有过类似的概念。腐蚀的方法是,拿B的中心点和X上的点一个一个地对比,如果B上的所有点都在X的范围内,则该点保留,否则将该点去掉;右边是腐蚀后的结果。可以看出,它仍在原来X的范围内,且比X包含的点要少,就象X被腐蚀掉了一层。

6.10   腐蚀运算

6.11为原图,图6.12为腐蚀后的结果图,能够很明显地看出腐蚀的效果。

6.11    原图

6.12   腐蚀后的结果图

下面的这段程序,实现了上述的腐蚀运算,针对的都是黑色点。参数中有一个BOOL变量,为真时,表示在水平方向进行腐蚀运算,即结构元素B  ;否则在垂直方向上进行腐蚀运算,即结构元素B  

腐蚀源码

BOOL Erosion(HWND hWnd,BOOL Hori)

{

       DWORD                             OffBits,BufSize;

LPBITMAPINFOHEADER    lpImgData;

       LPSTR                   lpPtr;

       HLOCAL                  hTempImgData;

       LPBITMAPINFOHEADER    lpTempImgData;

       LPSTR                            lpTempPtr;

       HDC                      hDc;

       HFILE                    hf;

       LONG                    x,y;

       unsigned char              num;

       int                        i;

//为了处理方便,仍采用256级灰度图,不过只用调色板中0和255两项

if( NumColors!=256){  

           MessageBox(hWnd,"Must be a mono bitmap with grayscale palette!",

"Error Message",MB_OK|MB_ICONEXCLAMATION);

return FALSE;

}

OffBits=bf.bfOffBits-sizeof(BITMAPFILEHEADER);

//BufSize为缓冲区大小

       BufSize=OffBits+bi.biHeight*LineBytes;

       //为新的缓冲区分配内存

       if((hTempImgData=LocalAlloc(LHND,BufSize))==NULL)

{

            MessageBox(hWnd,"Error alloc memory!","Error Message",

MB_OK|MB_ICONEXCLAMATION);

return FALSE;

    }

     lpImgData=(LPBITMAPINFOHEADER)GlobalLock(hImgData);    

       lpTempImgData=(LPBITMAPINFOHEADER)LocalLock(hTempImgData);

       //拷贝头信息和位图数据     

       memcpy(lpTempImgData,lpImgData,BufSize);

       if(Hori)

       {   

//在水平方向进行腐蚀运算

              for(y=0;y<bi.biHeight;y++){

                     //lpPtr指向原图数据,lpTempPtr指向新图数据

                     lpPtr=(char *)lpImgData+(BufSize-LineBytes-y*LineBytes)+1;

                     lpTempPtr=(char*)lpTempImgData+

(BufSize-LineBytes-y*LineBytes)+1;

                     for(x=1;x<bi.biWidth-1;x++){ 

//注意为防止越界,x的范围从1到宽度-2

                            num=(unsigned char)*lpPtr;

                            if (num==0){  //因为腐蚀掉的是黑点,所以只对黑点处理

                                   *lpTempPtr=(unsigned char)0;  //先置成黑点

                                   for(i=0;i<3;i++){

                                          num=(unsigned char)*(lpPtr+i-1);

                                          if(num==255){ 

//自身及上下邻居中若有一个不是黑点,则将该点腐

//蚀成白点

                                                 *lpTempPtr=(unsigned char)255;

                                                 break;

                                          }

                                   }

                            }

//原图中就是白点的,新图中仍是白点

                            else *lpTempPtr=(unsigned char)255;  

                            //指向下一个象素

                            lpPtr++; 

                            lpTempPtr++;

                     }

              }

       }

else{ 

//在垂直方向进行腐蚀运算

              for(y=1;y<bi.biHeight-1;y++){ //注意为防止越界,y的范围从1到高度-2

                     //lpPtr指向原图数据,lpTempPtr指向新图数据

                     lpPtr=(char *)lpImgData+(BufSize-LineBytes-y*LineBytes);

                     lpTempPtr=(char *)lpTempImgData+(BufSize-LineBytes-y*LineBytes);

                     for(x=0;x<bi.biWidth;x++){

                            num=(unsigned char)*lpPtr;

                            if (num==0){ //因为腐蚀掉的是黑点,所以只对黑点处理

                                   *lpTempPtr=(unsigned char)0; //先置成黑点

                                   for(i=0;i<3;i++){

                                          num=(unsigned char)*(lpPtr+(i-1)*LineBytes);

                                          if(num==255){

//自身及上下邻居中若有一个不是黑点,则将该点腐

//蚀成白点

                                                 *lpTempPtr=(unsigned char)255;

                                                 break;

                                          }

                                   }

                            }

//原图中就是白点的,新图中仍是白点

                            else *lpTempPtr=(unsigned char)255;

                            //指向下一个象素

                            lpPtr++;

                            lpTempPtr++;

                     }

              }

       }

    if(hBitmap!=NULL)

           DeleteObject(hBitmap);

       hDc=GetDC(hWnd);     

       //产生新的位图

       hBitmap=CreateDIBitmap(hDc,(LPBITMAPINFOHEADER)lpTempImgData,

(LONG)CBM_INIT,

(LPSTR)lpTempImgData+

sizeof(BITMAPINFOHEADER)+

                                         NumColors*sizeof(RGBQUAD),

(LPBITMAPINFO)lpTempImgData, DIB_RGB_COLORS);

       //起不同的结果文件名

       if(Hori)

              hf=_lcreat("c:\\herosion.bmp",0);

       else

              hf=_lcreat("c:\\verosion.bmp",0);

       _lwrite(hf,(LPSTR)&bf,sizeof(BITMAPFILEHEADER)); 

       _lwrite(hf,(LPSTR)lpTempImgData,BufSize);

       _lclose(hf);

       //释放内存及资源

ReleaseDC(hWnd,hDc);

       LocalUnlock(hTempImgData);

       LocalFree(hTempImgData);

       GlobalUnlock(hImgData);

       return TRUE;

}

膨胀

膨胀(dilation)可以看做是腐蚀的对偶运算,其定义是:把结构元素B平移a后得到Ba,若Ba击中X,我们记下这个a点。所有满足上述条件的a点组成的集合称做XB膨胀的结果。用公式表示为:D(X)={a | BaX}=X  B,如图6.13所示。图6.13X是被处理的对象,B是结构元素,不难知道,对于任意一个在阴影部分的点aBa击中X,所以XB膨胀的结果就是那个阴影部分。阴影部分包括X的所有范围,就象X膨胀了一圈似的,这就是为什么叫膨胀的原因。

同样,如果B不是对称的,XB膨胀的结果和X Bv膨胀的结果不同。

让我们来看看实际上是怎样进行膨胀运算的。在图6.14中,左边是被处理的图象X(二值图象,我们针对的是黑点),中间是结构元素B。膨胀的方法是,拿B的中心点和X上的点及X周围的点一个一个地对,如果B上有一个点落在X的范围内,则该点就为黑;右边是膨胀后的结果。可以看出,它包括X的所有范围,就象X膨胀了一圈似的。

6.13   膨胀的示意图

6.14   膨胀运算

6.15为图6.11膨胀后的结果图,能够很明显的看出膨胀的效果。

6.15   6.11膨胀后的结果图

下面的这段程序,实现了上述的膨胀运算,针对的都是黑色点。参数中有一个BOOL变量,为真时,表示在水平方向进行膨胀运算,即结构元素B  ;否则在垂直方向上进行膨胀运算,即结构元素B  

膨胀源码

BOOL Dilation(HWND hWnd,BOOL Hori)

{

       DWORD                             OffBits,BufSize;

LPBITMAPINFOHEADER    lpImgData;

       LPSTR                   lpPtr;

       HLOCAL                  hTempImgData;

       LPBITMAPINFOHEADER    lpTempImgData;

       LPSTR                     lpTempPtr;

       HDC                     hDc;

       HFILE                    hf;

       LONG                    x,y;

       unsigned char              num;

       int                        i;

//为了处理的方便,仍采用256级灰度图,不过只调色板中0和255两项

if( NumColors!=256){  

            MessageBox(hWnd,"Must be a mono bitmap with grayscale palette!",

"Error Message",MB_OK|MB_ICONEXCLAMATION);

return FALSE;

}

OffBits=bf.bfOffBits-sizeof(BITMAPFILEHEADER);

//BufSize为缓冲区大小

       BufSize=OffBits+bi.biHeight*LineBytes;

//为新的缓冲区分配内存

       if((hTempImgData=LocalAlloc(LHND,BufSize))==NULL)

    {

           MessageBox(hWnd,"Error alloc memory!","Error Message",

MB_OK|MB_ICONEXCLAMATION);

return FALSE;

    }

     lpImgData=(LPBITMAPINFOHEADER)GlobalLock(hImgData);    

       lpTempImgData=(LPBITMAPINFOHEADER)LocalLock(hTempImgData);

       //拷贝头信息和位图数据     

       memcpy(lpTempImgData,lpImgData,BufSize);

       if(Hori)

       {   

//在水平方向进行膨胀运算

              for(y=0;y<bi.biHeight;y++){

                     //lpPtr指向原图数据,lpTempPtr指向新图数据

                     lpPtr=(char *)lpImgData+(BufSize-LineBytes-y*LineBytes)+1;

                     lpTempPtr=(char*)lpTempImgData+

(BufSize-LineBytes-y*LineBytes)+1;

                     for(x=1;x<bi.biWidth-1;x++){ 

//注意为防止越界,x的范围从1到宽度-2

                            num=(unsigned char)*lpPtr;

//原图中是黑点的,新图中肯定也是,所以要考虑的是那些原图

//中的白点,看是否有可能膨胀成黑点

                            if (num==255){

                                   *lpTempPtr=(unsigned char)255; //先置成白点

                                   for(i=0;i<3;i++){ 

                                          num=(unsigned char)*(lpPtr+i-1);

//只要左右邻居中有一个是黑点,就膨胀成黑点

                                          if(num==0){

*lpTempPtr=(unsigned char)0;

                                                 break;

                                          }

                                   }

                            }

//原图中就是黑点的,新图中仍是黑点

                            else *lpTempPtr=(unsigned char)0;

                            //指向下一个象素

                            lpPtr++;

                            lpTempPtr++;

                     }

              }

       }

       else{

//在垂直方向进行腐蚀运算

              for(y=1;y<bi.biHeight-1;y++){ //注意为防止越界,y的范围从1到高度-2

              lpPtr=(char *)lpImgData+(BufSize-LineBytes-y*LineBytes);

                     lpTempPtr=(char *)lpTempImgData+(BufSize-LineBytes-y*LineBytes);

                     for(x=0;x<bi.biWidth;x++){

                            num=(unsigned char)*lpPtr;

                            if (num==255){

                                   *lpTempPtr=(unsigned char)255;

                                   for(i=0;i<3;i++){

                                          num=(unsigned char)*(lpPtr+(i-1)*LineBytes);

//只要上下邻居中有一个是黑点,就膨胀成黑点

                                          if(num==0){

                                                 *lpTempPtr=(unsigned char)0;

                                                 break;

                                          }

                                   }

                            }

                            else *lpTempPtr=(unsigned char)0;

                            lpPtr++;

                            lpTempPtr++;

                     }

              }

       }

    if(hBitmap!=NULL)

           DeleteObject(hBitmap);

       hDc=GetDC(hWnd);     

       //产生新的位图

       hBitmap=CreateDIBitmap(hDc,(LPBITMAPINFOHEADER)lpTempImgData,

(LONG)CBM_INIT,

(LPSTR)lpTempImgData+

sizeof(BITMAPINFOHEADER)+

                                         NumColors*sizeof(RGBQUAD),

(LPBITMAPINFO)lpTempImgData,

DIB_RGB_COLORS);

       //起不同的结果文件名

       if(Hori)

              hf=_lcreat("c:\\hdilation.bmp",0);

       else

              hf=_lcreat("c:\\vdilation.bmp",0);

       _lwrite(hf,(LPSTR)&bf,sizeof(BITMAPFILEHEADER)); 

       _lwrite(hf,(LPSTR)lpTempImgData,BufSize);

       _lclose(hf);

       //释放内存及资源

      ReleaseDC(hWnd,hDc);

       LocalUnlock(hTempImgData);

       LocalFree(hTempImgData);

       GlobalUnlock(hImgData);

       return TRUE;

}

腐蚀运算和膨胀运算互为对偶的,用公式表示为(X  B)c=(Xc  B),即B腐蚀后的补集等于X的补集被B膨胀。这句话可以形象的理解为:河岸的补集为河面,河岸的腐蚀等价于河面的膨胀。你可以自己举个例子来验证一下这个关系。在有些情况下,这个对偶关系是非常有用的。例如:某个图象处理系统用硬件实现了腐蚀运算,那么不必再另搞一套膨胀的硬件,直接利用该对偶就可以实现了。

先腐蚀后膨胀称为开(open),即OPEN(X)=D(E(X))

让我们来看一个开运算的例子(见图6.16)

6.16开运算

在图16上面的两幅图中,左边是被处理的图象X(二值图象,我们针对的是黑点),右边是结构元素B,下面的两幅图中左边是腐蚀后的结果;右边是在此基础上膨胀的结果。可以看到,原图经过开运算后,一些孤立的小点被去掉了。一般来说,开运算能够去除孤立的小点,毛刺和小桥(即连通两块区域的小点),而总的位置和形状不变。这就是开运算的作用。要注意的是,如果B是非对称的,进行开运算时要用B的对称集Bv膨胀,否则,开运算的结果和原图相比要发生平移。图6.17和图6.18能够说明这个问题。

6.17 B膨胀后,结果向左平移了

6.18   Bv膨胀后位置不变

6.17是用B膨胀的,可以看到,OPEN(X)向左平移了。图18是用Bv膨胀的,可以看到,总的位置和形状不变。

6.19为图6.11经过开运算后的结果。

6.19   6.11经过开运算后的结果

开运算的源程序可以很容易的根据上面的腐蚀,膨胀程序得到,这里就不给出了。

先膨胀后腐蚀称为闭(close),即CLOSE(X)=E(D(X))

让我们来看一个闭运算的例子(见图6.20)

6.20   闭运算

在图6.20上面的两幅图中,左边是被处理的图象X(二值图象,我们针对的是黑点),右边是结构元素B,下面的两幅图中左边是膨胀后的结果,右边是在此基础上腐蚀的结果可以看到,原图经过闭运算后,断裂的地方被弥合了。一般来说,闭运算能够填平小湖(即小孔),弥合小裂缝,而总的位置和形状不变。这就是闭运算的作用。同样要注意的是,如果B是非对称的,进行闭运算时要用B的对称集Bv膨胀,否则,闭运算的结果和原图相比要发生平移。

6.21为图6.11经过闭运算后的结果。

6.21   .611经过闭运算后的结果

闭运算的源程序可以很容易的根据上面的膨胀,腐蚀程序得到,这里就不给出了。

你大概已经猜到了,开和闭也是对偶运算,的确如此。用公式表示为(OPEN(X))c=CLOSE((Xc)),或者(CLOSE(X))c=OPEN((Xc))。即开运算的补集等于X的补集的闭运算,或者闭运算的补集等于X的补集的开运算。这句话可以这样来理解:在两个小岛之间有一座小桥,我们把岛和桥看做是处理对象X,则X的补集为大海。如果涨潮时将小桥和岛的外围淹没(相当于用尺寸比桥宽大的结构元素对X进行开运算),那么两个岛的分隔,相当于小桥两边海域的连通(Xc做闭运算)

细化

细化(thinning)算法有很多,我们在这里介绍的是一种简单而且效果很好的算法,用它就能够实现从文本抽取骨架的功能。我们的对象是白纸黑字的文本,但在程序中为了处理的方便,还是采用256级灰度图,不过只用到了调色板中0255两项。

所谓细化,就是从原来的图中去掉一些点,但仍要保持原来的形状。实际上,是保持原图的骨架。所谓骨架,可以理解为图象的中轴,例如一个长方形的骨架是它的长方向上的中轴线;正方形的骨架是它的中心点;圆的骨架是它的圆心,直线的骨架是它自身,孤立点的骨架也是自身。文本的骨架嘛,前言中的例子显示的很明白。那么怎样判断一个点是否能去掉呢?显然,要根据它的八个相邻点的情况来判断,我们给几个例子(如图6.22所示)

6.22   根据某点的八个相邻点的情况来判断该点是否能删除

6.22中,(1)不能删,因为它是个内部点,我们要求的是骨架,如果连内部点也删了,骨架也会被掏空的;(2)不能删,和(1)是同样的道理;(3)可以删,这样的点不是骨架;(4)不能删,因为删掉后,原来相连的部分断开了;(5)可以删,这样的点不是骨架;(6)不能删,因为它是直线的端点,如果这样的点删了,那么最后整个直线也被删了,剩不下什么;(7)不能删,因为孤立点的骨架就是它自身。

总结一下,有如下的判据:(1)内部点不能删除;(2)孤立点不能删除;(3)直线端点不能删除;(4)如果P是边界点,去掉P后,如果连通分量不增加,则P可以删除。

我们可以根据上述的判据,事先做出一张表,从0255共有256个元素,每个元素要么是0,要么是1。我们根据某点(当然是要处理的黑色点了)的八个相邻点的情况查表,若表中的元素是1,则表示该点可删,否则保留。

查表的方法是,设白点为1,黑点为0;左上方点对应一个8位数的第一位(最低位),正上方点对应第二位,右上方点对应的第三位,左邻点对应第四位,右邻点对应第五位,左下方点对应第六位,正下方点对应第七位,右下方点对应的第八位,按这样组成的8位数去查表即可。例如上面的例子中(1)对应表中的第0项,该项应该为0(2)对应37,该项应该为0(3)对应173,该项应该为1(4)对应231,该项应该为0(5)对应237,该项应该为1(6)对应254,该项应该为0(7)对应255,该项应该为0

这张表我已经替大家做好了,可花了我不少时间呢!

static int erasetable[256]={

                                         0,0,1,1,0,0,1,1,          1,1,0,1,1,1,0,1,

                                   1,1,0,0,1,1,1,1,             0,0,0,0,0,0,0,1,

                                          0,0,1,1,0,0,1,1,             1,1,0,1,1,1,0,1,

                                          1,1,0,0,1,1,1,1,             0,0,0,0,0,0,0,1,

                                          1,1,0,0,1,1,0,0,             0,0,0,0,0,0,0,0,

                                          0,0,0,0,0,0,0,0,             0,0,0,0,0,0,0,0,

                                          1,1,0,0,1,1,0,0,             1,1,0,1,1,1,0,1,

                                   0,0,0,0,0,0,0,0,             0,0,0,0,0,0,0,0,

                           0,0,1,1,0,0,1,1,             1,1,0,1,1,1,0,1,

                                          1,1,0,0,1,1,1,1,             0,0,0,0,0,0,0,1,

                                          0,0,1,1,0,0,1,1,             1,1,0,1,1,1,0,1,

                                          1,1,0,0,1,1,1,1,             0,0,0,0,0,0,0,0,

                                          1,1,0,0,1,1,0,0,             0,0,0,0,0,0,0,0,

                                1,1,0,0,1,1,1,1,             0,0,0,0,0,0,0,0,

                                          1,1,0,0,1,1,0,0,             1,1,0,1,1,1,0,0,

                                   1,1,0,0,1,1,1,0,             1,1,0,0,1,0,0,0

                                     };

有了这张表,算法就很简单了,每次对一行一行的将整个图象扫描一遍,对于每个点(不包括边界点),计算它在表中对应的索引,若为0,则保留,否则删除该点。如果这次扫描没有一个点被删除,则循环结束,剩下的点就是骨架点,如果有点被删除,则进行新的一轮扫描,如此反复,直到没有点被删除为止。

实际上,该算法有一些缺陷。举个简单的例子,有一个黑色矩形,如图6.23所示。

6.23经过细化后,我们预期的结果是一条水平直线,且位于该黑色矩形的中心。实际的结果确实是一条水平直线,但不是位于黑色矩形的中心,而是最下面的一条边。

为什么会这样,我们来分析一下:在从上到下,从左到右的扫描过程中,我们遇到的第一个黑点就是黑色矩形的左上角点,经查表,该点可以删。下一个点是它右边的点,经查表,该点也可以删,如此下去,整个一行被删了。每一行都是同样的情况,所以都被删除了。到了最后一行时,黑色矩形已经变成了一条直线,最左边的黑点不能删,因为它是直线的端点,它右边的点也不能删,因为如果删除,直线就断了,如此下去,直到最右边的点,也不能删,因为它是直线的右端点。所以最下面的一条边保住了,但这并不是我们希望的结果。

解决的办法是,在每一行水平扫描的过程中,先判断每一点的左右邻居,如果都是黑点,则该点不做处理。另外,如果某个黑点被删除了,那么跳过它的右邻居,处理下一个点。这样就避免了上述的问题。

6.23  黑色矩形

6.24  6.23细化后的结果

解决了上面的问题,我们来看看处理后的结果,如图6.24所示。这次变成一小段竖线了,还是不对,是不是很沮丧?别着急,让我们再来分析一下:在上面的算法中,我们遇到的第一个能删除的点就是黑色矩形的左上角点;第二个是第一行的最右边的点,即黑色矩形的右上角点;第三个是第二行的最左边的点;第四个是第二行的最右边的点;……;整个图象处理这样一次后,宽度减少2。每次都是如此,直到剩最中间一列,就不能再删了。为什么会这样呢?原因是这样的处理过程只实现了水平细化,如果在每一次水平细化后,再进行一次垂直方向的细化(只要把上述过程的行列换一下),就可以了。

这样一来,每处理一次,删除点的顺序变成:(先是水平方向扫描)第一行最左边的点;第一行最右边的点;第二行最左边的点;第二行最右边的点;……最后一行最左边的点;最后一行最右边的点;(然后是垂直方向扫描)第二列最上边的点(因为第一列最上边的点已被删除);第二列最下边的点;第三列最上边的点;第三列最下边的点;……倒数第二列最上边的点(因为倒数第一列最上边的点已被删除);倒数第二列最下边的点。我们发现,刚好剥掉了一圈,这也正是细化要做的事。实际的结果也验证了我们的想法。

以下是源程序,黑体字部分是值得注意的地方。

细化源码

BOOL Thinning(HWND hWnd)

{

       DWORD                             OffBits,BufSize;

     LPBITMAPINFOHEADER    lpImgData;

       LPSTR                            lpPtr;

       HLOCAL                  hTempImgData;

       LPBITMAPINFOHEADER    lpTempImgData;

       LPSTR                   lpTempPtr;

       HDC                      hDc;

       HFILE                    hf;

       LONG                    x,y;

       int                                        num;

       BOOL                     Finished;

       int                        nw,n,ne,w,e,sw,s,se;

//为了处理的方便,仍采用256级灰度图,不过只用调色板中0和255两项

       if( NumColors!=256){

MessageBox(hWnd,"Must be a mono bitmap with grayscale palette!",

"Error Message",MB_OK|MB_ICONEXCLAMATION);

return FALSE;

}

OffBits=bf.bfOffBits-sizeof(BITMAPFILEHEADER);

//BufSize为缓冲区大小

       BufSize=OffBits+bi.biHeight*LineBytes;

//为新的缓冲区分配内存

       if((hTempImgData=LocalAlloc(LHND,BufSize))==NULL)

{

            MessageBox(hWnd,"Error alloc memory!","Error Message",

MB_OK|MB_ICONEXCLAMATION);

return FALSE;

}

     lpImgData=(LPBITMAPINFOHEADER)GlobalLock(hImgData);    

       lpTempImgData=(LPBITMAPINFOHEADER)LocalLock(hTempImgData);

       //拷贝头信息和位图数据     

       memcpy(lpTempImgData,lpImgData,BufSize);

       //结束标志置成假

       Finished=FALSE;

while(!Finished){ //还没有结束

              //结束标志置成假

            Finished=TRUE;

       //先进行水平方向的细化

              for (y=1;y<bi.biHeight-1;y++){ //注意为防止越界,y的范围从1到高度-2

                     //lpPtr指向原图数据,lpTempPtr指向新图数据

                     lpPtr=(char *)lpImgData+(BufSize-LineBytes-y*LineBytes);

                     lpTempPtr=(char *)lpTempImgData+(BufSize-LineBytes-y*LineBytes);

                     x=1; //注意为防止越界,x的范围从1到宽度-2

                     while(x<bi.biWidth-1){

                            if(*(lpPtr+x)==0){ //是黑点才做处理

                                   w=(unsigned char)*(lpPtr+x-1);  //左邻点

                                   e=(unsigned char)*(lpPtr+x+1);  //右邻点

                                   if( (w==255)|| (e==255)){ 

//如果左右两个邻居中至少有一个是白点才处理

                                          nw=(unsigned char)*(lpPtr+x+LineBytes-1); //左上邻点

                                          n=(unsigned char)*(lpPtr+x+LineBytes); //上邻点

                                          ne=(unsigned char)*(lpPtr+x+LineBytes+1); //右上邻点

                                          sw=(unsigned char)*(lpPtr+x-LineBytes-1); //左下邻点

                                          s=(unsigned char)*(lpPtr+x-LineBytes); //下邻点

                                          se=(unsigned char)*(lpPtr+x-LineBytes+1); //右下邻点

                                          //计算索引

                            num=nw/255+n/255*2+ne/255*4+w/255*8+e/255*16+

sw/255*32+s/255*64+se/255*128;

                                          if(erasetable[num]==1){ //经查表,可以删除

//在原图缓冲区中将该黑点删除

                                                 *(lpPtr+x)=(BYTE)255; 

//结果图中该黑点也删除

                                                 *(lpTempPtr+x)=(BYTE)255; 

                                                 Finished=FALSE; //有改动,结束标志置成假

                                                 x++; //水平方向跳过一个象素

                                          }

                                   }

                            }

                            x++; //扫描下一个象素

                     }

              }

       //再进行垂直方向的细化

              for (x=1;x<bi.biWidth-1;x++){ //注意为防止越界,x的范围从1到宽度-2

                     y=1; //注意为防止越界,y的范围从1到高度-2

                     while(y<bi.biHeight-1){

                            lpPtr=(char *)lpImgData+(BufSize-LineBytes-y*LineBytes);

                            lpTempPtr=(char*)lpTempImgData+

(BufSize-LineBytes-y*LineBytes);

                            if(*(lpPtr+x)==0){ //是黑点才做处理

                                   n=(unsigned char)*(lpPtr+x+LineBytes);

                                   s=(unsigned char)*(lpPtr+x-LineBytes);

                                   if( (n==255)|| (s==255)){

//如果上下两个邻居中至少有一个是白点才处理

                                          nw=(unsigned char)*(lpPtr+x+LineBytes-1);

                                          ne=(unsigned char)*(lpPtr+x+LineBytes+1);

                                          w=(unsigned char)*(lpPtr+x-1);

                                          e=(unsigned char)*(lpPtr+x+1);

                                          sw=(unsigned char)*(lpPtr+x-LineBytes-1);

                                          se=(unsigned char)*(lpPtr+x-LineBytes+1);

                                          //计算索引

num=nw/255+n/255*2+ne/255*4+w/255*8+e/255*16+

sw/255*32+s/255*64+se/255*128;

                                          if(erasetable[num]==1){ //经查表,可以删除

//在原图缓冲区中将该黑点删除

                                                 *(lpPtr+x)=(BYTE)255; 

//结果图中该黑点也删除

                                                 *(lpTempPtr+x)=(BYTE)255; 

                                                 Finished=FALSE; //有改动,结束标志置成假

                                                 y++;//垂直方向跳过一个象素

                                          }

                                   }

                            }

                            y++; //扫描下一个象素

                     }

              } 

}

     if(hBitmap!=NULL)

           DeleteObject(hBitmap);

       hDc=GetDC(hWnd);     

       //产生新的位图

       hBitmap=CreateDIBitmap(hDc,(LPBITMAPINFOHEADER)lpTempImgData,

(LONG)CBM_INIT,

(LPSTR)lpTempImgData+

sizeof(BITMAPINFOHEADER)+

NumColors*sizeof(RGBQUAD),

(LPBITMAPINFO)lpTempImgData,

DIB_RGB_COLORS);

hf=_lcreat("c:\\thinning.bmp",0);

       _lwrite(hf,(LPSTR)&bf,sizeof(BITMAPFILEHEADER)); 

       _lwrite(hf,(LPSTR)lpTempImgData,BufSize);

       _lclose(hf);

       //释放内存及资源

      ReleaseDC(hWnd,hDc);

       LocalUnlock(hTempImgData);

       LocalFree(hTempImgData);

       GlobalUnlock(hImgData);

       return TRUE;

}

另外补充说明一下,助于理解

腐蚀:删除对象边界的某些像素

膨胀:给图像中的对象边界添加像素

算法:

膨胀算法:用3X3的结构元素,扫描二值图像的每一个像素,用结构元素与其覆盖的二值图像做“与”运算,如果都为0,结构图像的该像素为0,否则为1.结果:使二值图像扩大一圈。

腐蚀算法:用3X3的结构元素,扫描二值图像的每一个像素,用结构元素与其覆盖的二值图像做“与”运算,如果都为1,结构图像的该像素为1,否则为0.结果:使二值图像减小一圈。

 




2020-03-21 20:46:43 Ibelievesunshine 阅读数 247

如果您觉得本文不错!记得点赞哦!


一. 图像形态学简介:

图解图像腐蚀、膨胀 ↑

    经验之谈:形态学操作一般作用于二值图像,来连接相邻的元素(膨胀)或分离成独立的元素(侵蚀)。腐蚀和膨胀是针对图片中的白色(即前景)部分!


二. 图像形态学操作 膨胀和腐蚀的算法:

    膨胀算法:

        对于待操作的像素 f(x,y),不论 f(x,y-1) 、f(x,y+1) 、f(x-1,y) 、f(x+1,y) 哪一个为255,则 f(x,y)=255。

膨胀操作 ↑

 

        换句话说:将待操作的图像像素与以下  4-近邻矩阵 相乘,结果大于255的话,将中心像素设为255。

膨胀:待操作像素 * 上面矩阵 > =255,f(x,y) = 255。 ↑

 

    腐蚀算法:

        对于待操作的像素 f(x,y),只有 f(x,y-1) 、f(x,y+1) 、f(x-1,y) 、f(x+1,y) 都为255,则 f(x,y)=255。

        换句话说:将待操作的图像像素与以下  4-近邻矩阵 相乘,结果小于255*4的话,将中心像素设为0。

腐蚀:待操作像素 * 上面矩阵 < 255*4,f(x,y) = 0 。↑


三. python实现图像膨胀和腐蚀

# Writer : wojianxinygcl@163.com
# Date   : 2020.3.21
import cv2
import numpy as np
import matplotlib.pyplot as plt

# Gray scale
def BGR2GRAY(img):
	b = img[:, :, 0].copy()
	g = img[:, :, 1].copy()
	r = img[:, :, 2].copy()

	# Gray scale
	out = 0.2126 * r + 0.7152 * g + 0.0722 * b
	out = out.astype(np.uint8)

	return out

# Otsu Binalization
def otsu_binarization(img, th=128):
	H, W = img.shape
	out = img.copy()

	max_sigma = 0
	max_t = 0

	# determine threshold
	for _t in range(1, 255):
		v0 = out[np.where(out < _t)]
		m0 = np.mean(v0) if len(v0) > 0 else 0.
		w0 = len(v0) / (H * W)
		v1 = out[np.where(out >= _t)]
		m1 = np.mean(v1) if len(v1) > 0 else 0.
		w1 = len(v1) / (H * W)
		sigma = w0 * w1 * ((m0 - m1) ** 2)
		if sigma > max_sigma:
			max_sigma = sigma
			max_t = _t

	# Binarization
	print("threshold >>", max_t)
	th = max_t
	out[out < th] = 0
	out[out >= th] = 255

	return out


# Morphology Dilate
def Morphology_Dilate(img, Dil_time=1):
	H, W = img.shape

	# kernel
	MF = np.array(((0, 1, 0),
				(1, 0, 1),
				(0, 1, 0)), dtype=np.int)

	# each dilate time
	out = img.copy()
	for i in range(Dil_time):
		tmp = np.pad(out, (1, 1), 'edge')
		for y in range(1, H):
			for x in range(1, W):
				if np.sum(MF * tmp[y-1:y+2, x-1:x+2]) >= 255:
					out[y, x] = 255

	return out


# Morphology Erode
def Morphology_Erode(img, Erode_time=1):
	H, W = img.shape
	out = img.copy()

	# kernel
	MF = np.array(((0, 1, 0),
				(1, 0, 1),
				(0, 1, 0)), dtype=np.int)

	# each erode
	for i in range(Erode_time):
		tmp = np.pad(out, (1, 1), 'edge')
		# erode
		for y in range(1, H):
			for x in range(1, W):
				if np.sum(MF * tmp[y-1:y+2, x-1:x+2]) < 255*4:
					out[y, x] = 0

	return out


# Read image
img = cv2.imread("../paojie.jpg").astype(np.float32)

# Grayscale
gray = BGR2GRAY(img)

# Otsu's binarization
otsu = otsu_binarization(gray)

# Morphology - dilate
erode_result = Morphology_Erode(otsu, Erode_time=2)
dilate_result = Morphology_Dilate(otsu,Dil_time=2)

# Save result
cv2.imwrite("Black_and_white.jpg",otsu)
cv2.imshow("Black_and_white",otsu)
cv2.imwrite("erode_result.jpg", erode_result)
cv2.imshow("erode_result", erode_result)
cv2.imwrite("dilate_result.jpg", dilate_result)
cv2.imshow("dilate_result",dilate_result)
cv2.waitKey(0)
cv2.destroyAllWindows()

四. 实验结果:

二值图像(左),膨胀图像(中),侵蚀图像(右) ↑


五. 参考内容:

    ①  https://www.jianshu.com/p/ba2cec49c981

    ②  https://www.cnblogs.com/wojianxin/p/12542004.html


六. 版权声明:

    未经作者允许,请勿随意转载抄袭,抄袭情节严重者,作者将考虑追究其法律责任,创作不易,感谢您的理解和配合!

2018-08-12 20:01:25 hanshanbuleng 阅读数 2995

用C语言实现图像的膨胀与腐蚀算法

经过几次学习opencv源代码,我决定自己动手写一下膨胀与腐蚀算法,如果具体算法原理不明确的话,可以看看前几篇我总结的膨胀腐蚀算法原理:


腐蚀算法

/*****************************************************
function: achieve the erode algorithm of the binary image,using 3*3 structural elements 
parameter:
  1 IplImage* img : original image
  2 int *elementArray : 3*3 structural elements
return: Run code   0: right  -1: error
*****************************************************/
//在画图中 理解有中心点的位置,并将结果赋给结构元素中心点所对应的位置,但是在实际的编程中是遍历每一个点(i,j),可以理解为对应结构元素的中心点
int ErodeYang01(IplImage* srcImg, IplImage* dstImg, int *elementArray)
{
    int ret = 0;
    int imgHeight, imgWidth;
    int i, j, k, l;
    int rowPosi, colPosi;   //为结构元素在计算图像中的位置 rowPosi:行 colPosi:列

    bool isMatch;

    if (srcImg == NULL || elementArray == NULL)
        return ret = -1;

    imgHeight = srcImg->height;
    imgWidth = srcImg->width;

    memset((void*)dstImg->imageData, 0, dstImg->imageSize);    // memset((void*)dstImg->imageData, 255, dstImg->imageSize); 

    for (i = 1; i < imgHeight; i++)                             //i,j 都从1开始 是防止逐行扫描时,访问越界,四周留出一个像素的宽度
    {
        for (j = 1; j < imgWidth; j++)
        {
            isMatch = true;
            for (k = 0; k < 3; k++)                             //k,l 为结构元素的遍历
            {
                for (l = 0; l < 3; l++)
                {
                    rowPosi = (i - 1 + k)*srcImg->widthStep;
                    colPosi = j - 1 + l;

                    if (elementArray[3 * k + l] == -1)  //此点不关心
                        continue;

                    if (elementArray[3 * k + l] == 1)  //前景
                    {
                        if (srcImg->imageData[rowPosi + colPosi] != -1)   //opencv中 二值化图结果只是存放了 0 与 -1  但正常理解二值化值应该 0与1 或是 0与255
                        {
                            isMatch = false;
                            break;
                        }
                    }   
                    else
                    {
                        printf("structural elements exist illegal values");
                        return ret = -1;
                    }

                }

            }

            if (isMatch)
                dstImg->imageData[i*dstImg->widthStep + j] = 255;   //赋值为图像中位置点

        }
    }

    return ret;
}

膨胀算法

/*****************************************************
function: achieve the dilate algorithm of the binary image,using 3*3 structural elements
parameter:
1 IplImage* img : original image
2 int *elementArray : 3*3 structural elements
return: Run code   0: right  -1: error
*****************************************************/
int DilateYang(IplImage* img, IplImage* dstImg, int *elementArray /* elementArray[3][3] */)
{
    int ret = 0;
    int imgHeight, imgWidth;
    int i, j;    //图像循环变量
    int i_, j_;  //结构元素的对称集循环变量
    int k, l;    //结构元素循环变量

    uchar temp;

    bool isMatch;

    if (img == NULL || elementArray == NULL)
        return ret = -1;

    imgHeight = img->height;
    imgWidth = img->width;

    memset((void*)dstImg->imageData,0,dstImg->imageSize);    //这是模仿程序写的 没确定是否正确 初始化整幅图都是白色

    //计算结构元素的对称集 --算法中是利用结构元素的对称集进行后续处理的
    for (i = 0; i < 2; i++)
    {
        for (j = 0; j < 3 - i; j++)
        {
            temp = elementArray[i * 3 + j];
            elementArray[i * 3 + j] = elementArray[(2 - i) * 3 + (2 - j)];
            elementArray[(2 - i) * 3 + (2 - j)] = temp;
        }
    }


    for (i = 1; i < imgHeight; i++)
    {
        for (j = 1; j < imgWidth; j++)
        {

            for (k = 0; k < 3; k++)
            {
                for (l = 0; l < 3; l++)
                {
                    int rowImg = (i - 1 + k)*img->widthStep;
                    int colImg = j - l + 1;

                    if (elementArray[k * 3 + l] == -1)
                        continue;

                    if (elementArray[k * 3 + l] == 1)
                    {
                        //if(img->imageData[rowImg + colImg] != 0)
                        //  printf("%d", img->imageData[rowImg + colImg]);

                        if (img->imageData[rowImg + colImg] == -1)   //opencv中 二值化图结果只是存放了 0 与 -1
                        {
                            dstImg->imageData[i*dstImg->widthStep + j] = 255;
                            break;
                        }

                    }
                    else
                    {
                        printf("structural elements exist illegal values");
                        return ret = -1;
                    }
                }
            }
        }
    }

    return ret;
}
2018-07-01 15:11:10 chen134225 阅读数 6581
在图像处理技术中,有一些的操作会对图像的形态发生改变,这些操作一般称之为形态学操作(phology)。数学形态学是基于集合论的图像处理方法,最早出现在生物学的形态与结构中,图像处理中的形态学操作用于图像与处理操作(去噪,形状简化)图像增强(骨架提取,细化,凸包及物体标记)、物体背景分割及物体形态量化等场景中,形态学操作的对象是二值化图像。
有名的形态学操作中包括腐蚀,膨胀,开操作,闭操作等。其中腐蚀,膨胀是许多形态学操作的基础。
  
腐蚀操作:
  顾名思义,是将物体的边缘加以腐蚀。具体的操作方法是拿一个宽m,高n的矩形作为模板,对图像中的每一个像素x做如下处理:像素x至于模板的中心,根据模版的大小,遍历所有被模板覆盖的其他像素,修改像素x的值为所有像素中最小的值。这样操作的结果是会将图像外围的突出点加以腐蚀。如下图的操作过程:

                                                                             图5 腐蚀操作原理

上图演示的过程是背景为黑色,物体为白色的情况。腐蚀将白色物体的表面加以“腐蚀”。在opencv的官方教程中,是以如下的图示说明腐蚀过程的,与我上面图的区别在于:背景是白色,而物体为黑色(这个不太符合一般的情况,所以我没有拿这张图作为通用的例子)。读者只需要了解背景为不同颜色时腐蚀也是不同的效果就可以了。


                                                        图6 腐蚀操作原理2

简单来说在背景为白(1),情景为黑的图像(0)中,核与其覆盖的图像部分做“与”操作,如果全为1,则该像素点为1,否则为0;也就是0容易得到,图像更多的地方变黑了,白色部分被腐蚀了。


膨胀操作:

  膨胀操作与腐蚀操作相反,是将图像的轮廓加以膨胀。操作方法与腐蚀操作类似,也是拿一个矩形模板,对图像的每个像素做遍历处理。不同之处在于修改像素的值不是所有像素中最小的值,而是最大的值。这样操作的结果会将图像外围的突出点连接并向外延伸。如下图的操作过程:


                                                                           图7 膨胀操作原理

下面是在opencv的官方教程中,膨胀过程的图示:


                                                       图8 膨胀操作原理2

简单来说在背景为白(1),前景为黑(0)的图像中,核与其覆盖的图像部分做“与”操作,如果全为0,则该像素点为0,否则为1;也就是1容易得到,图像更多的地方变白了,白色部分膨胀了。


开操作:

作用:放大裂缝和低密度区域,消除小物体,在平滑较大物体的边界时,不改变其面积。消除物体表面的突起。

开操作就是对图像先腐蚀,再膨胀。其中腐蚀与膨胀使用的模板是一样大小的。为了说明开操作的效果,请看下图的操作过程:


                                                                          图9 开操作原理

由于开操作是先腐蚀,再膨胀。因此可以结合图5和图7得出图9,其中图5的输出是图7的输入,所以开操作的结果也就是图7的结果。


闭操作:

作用:排除小型黑洞,突触了比原图轮廓区域更暗的区域,将两个区域连接起来,形成连通域。

  闭操作就是对图像先膨胀,再腐蚀。闭操作的结果一般是可以将许多靠近的图块相连称为一个无突起的连通域。在我们的图像定位中,使用了闭操作去连接所有的字符小图块,然后形成一个车牌的大致轮廓。闭操作的过程我会讲的细致一点。为了说明字符图块连接的过程。在这里选取的原图跟上面三个操作的原图不大一样,是一个由两个分开的图块组成的图。原图首先经过膨胀操作,将两个分开的图块结合起来(注意我用偏白的灰色图块表示由于膨胀操作而产生的新的白色)。接着通过腐蚀操作,将连通域的边缘和突起进行削平(注意我用偏黑的灰色图块表示由于腐蚀被侵蚀成黑色图块)。最后得到的是一个无突起的连通域(纯白的部分)。


                                                                         图10 闭操作原理

4.代码

  在opencv中,调用闭操作的方法是首先建立矩形模板,矩形的大小是可以设置的,由于矩形是用来覆盖以中心像素的所有其他像素,因此矩形的宽和高最好是奇数。

  通过以下代码设置矩形的宽和高。

Mat element = getStructuringElement(MORPH_RECT, Size(m_MorphSizeWidth, m_MorphSizeHeight) );
在这里,我们使用了类成员变量,这两个类成员变量在构造函数中被赋予了初始值。宽是17,高是3.
  设置完矩形的宽和高以后,就可以调用形态学操作了。opencv中所有形态学操作有一个统一的函数,通过参数来区分不同的具体操作。例如MOP_CLOSE代表闭操作,MOP_OPEN代表开操作。
morphologyEx(img_threshold, img_threshold, MORPH_CLOSE, element);


2016-05-17 01:17:08 wozhengtao 阅读数 22623

膨胀与腐蚀算法

  对图像处理有所了解的人都知道图像的形态学处理里最为基础的膨胀和腐蚀算法。二值图像即只有黑白两种颜色组成的图像,一般的白色为内容,黑色为背景。其实简单点理解二值图像的膨胀与腐蚀,腐蚀即是删除对象边界某些像素,也就是让白色的区域瘦一圈;而膨胀则是给图像中的对象边界添加像素,即让白色的区域胖上一圈。而这个“圈”的大小,则是由参数来指定的。下面的表展示了一幅图像经过膨胀和腐蚀算法的结果。可以看出膨胀让白色区域大了一圈,白色区域若有较小的黑色洞,则洞被填上;而腐蚀则让白色区域瘦了一圈,一些面积很小的白色区域则消失。

原图 腐蚀结果 膨胀结果

  腐蚀膨胀的算法原理并不复杂,而且网上有太多的文章都着重于介绍算法的原理思路,对用具体代码实现算法的方式讨论的不多,因而本文专注于几种实现膨胀腐蚀算法的方法。本文介绍了几种不同的腐蚀膨胀算法的实现,每一种实现都各有特点,今后若有更多的方法,也会继续更新加入至本文中。

 

结构元素与窗口形态

  在介绍算法逻辑之前,有必要先介绍跟腐蚀膨胀算法关系密切的结构元素。结构元素是形态学的基本算子,合理选取结构元素直接影响图像处理的效果和质量。结构元素的选择在于结构元素的形状和尺寸。结构元素可以有不同的形状,圆形、正方形、菱形、六边形、线段形都是可以选择的形状。圆形结构元素,由于各向同性,因此可以得到与方向无关的运算结果,正方形、菱形可以看作是圆盘形的变异。不同形状的结构元素运算结果会有差异,应针对待处理图像的几何形状进行选择。下表对比了正方形,圆形和正菱形三种结构元素形态。

预览
ElementSize 121 98 61
WindowSize 11×11 (r=5) 11×11 (r=5) 11×11 (r=5)
非空点个数计算方式 |x|<=r&&|y|<=r x2+y2<=r2 |x|+|y|<=r

  在算法实现过程中,往往会将卷积窗口中所有像素相对中心像素的偏移存在一个数组之中,这样在对不规则形状的卷积窗口进行处理时,可以不重复判断结构元素中哪些位置为有效位置,能减少计算次数。在实现之前首先贴上基本数据结构的实现,其中visit_count用来记录像素的访问次数:

复制代码
#define byte unsigned char
struct IntDouble
{
    int X;
    int Y;
    IntDouble(int x,int y)
    {
        this->X=x;
        this->Y=y;
    }
    IntDouble()
    {
        this->X=0;
        this->Y=0;
    }
};
class Bitmap2d
{
private:
    byte* data;
    int width;
    int height;
public:
    int visit_count;
    Bitmap2d(int width,int height,byte v)
    {
        this->data=new byte[width*height];
        memset(data,v,width*height*sizeof(byte));
        this->width=width;
        this->height=height;
        this->visit_count=0;
    }
    Bitmap2d(Bitmap2d& bitmap)
    {
        this->width=bitmap.Width();
        this->height=bitmap.Height();
        this->data=new byte[width*height];
        memcpy(this->data,bitmap.data,sizeof(byte)*Length());
        this->visit_count=0;
    }
    ~Bitmap2d()
    {
        delete[] data;
    }
    inline byte GetValue(int x,int y)
    {
        visit_count++;
        return data[x+y*width];
    }
    inline void SetValue(int x,int y,byte v)
    {
        visit_count++;
        data[x+y*width]=v;
    }
    inline int Width()
    {
        return width;
    }
    inline int Height()
    {
        return height;
    }
    inline int Length()
    {
        return width*height;
    }
    inline bool InRange(int x,int y)
    {
        return x>=0&&x<width&&y>=0&&y<height;
    }
    void ReadRaw(const char* filename)
    {
        FILE* file=fopen(filename,"rb");
        fread(data,sizeof(byte),Length(),file);
        fclose(file);
    }
    void SaveRaw(const char* filename)
    {
        FILE *const file = fopen(filename,"wb");
        fwrite(data,sizeof(byte),Length(),file);
        fclose(file);
    }
};
复制代码

 

实现思路1—根据定义直接算

  首先最为简单的思路是按算法基本原理直接正向求取输出图片的像素值:

  • 膨胀:对于输出图像的所有像素点P,调查原图像中对应窗口中的像素集合S,若S中至少有一个255,则P为255。
  • 腐蚀:对于输出图像的所有像素点P,调查原图像中对应窗口中的像素集合S,若S中至少有一个0,则P为0。

  该算法的腐蚀实现函数(膨胀的类似,就不重复贴,代码工程里有)如下:

复制代码
Bitmap2d* Execute()
{
    Bitmap2d* newBmp=new Bitmap2d(bmp.Width(),bmp.Height(),0);
    for(int j=0;j<bmp.Height();j++)
    {
        for(int i=0;i<bmp.Width();i++)
        {
            if(HasBlackInWindow(this->bmp,i,j))
                newBmp->SetValue(i,j,0);
            else
                newBmp->SetValue(i,j,255);
        }
    }
    return newBmp;
}
复制代码
复制代码
bool HasBlackInWindow(Bitmap2d& bmp,int i,int j)
{
    for(size_t k=0;k<winOffsets.size();k++)
    {
        int tx=i+winOffsets[k].X;
        int ty=j+winOffsets[k].Y;
        if(!bmp.InRange(tx,ty))
            continue;
        if(bmp.GetValue(tx,ty)==0)
        {
            return true;
        }
    }
    return false;
}
复制代码

  膨胀腐蚀算法的主要时间开销来至于对像素的访问,从上述实现不难得该算法对于n*n的位图和k*k大小正方形结构元素访问像素的次数为k2n2。事实上这是实现腐蚀膨胀算法最直接但也是最慢的方式。下图是Engine数据的一个切片二值化之后分别用正方形、圆形和菱形为结构元素膨胀和腐蚀的效果图:

腐蚀(正方形) 腐蚀(圆形) 腐蚀(菱形)
原图预览 膨胀(正方形) 膨胀(圆形) 膨胀(菱形)

 

实现思路2—跳过一些内部区域

  考虑到思路1的算法逻辑耗费在访问黑色区域和白色区域内部的时间较多,我们可以考虑只对黑白交界的边界考察窗口像素。这样的过程我们就可以想象成一把具有尺寸的刷子,膨胀算法刷子为白色,腐蚀算法刷子为黑色,然后让刷子在黑白交界的地方刷过,这样的过程生成的结果等价于思路1的结果。

  其优化的部分是针对远离边界的内部区域的涂刷,这样就能很大程度上减少像素的访问次数。不难想象出,对远离边界的内部区域的涂刷是不起效果的,这就是思路2对思路1改进的主要原因。设算法对n*n的图像按k*k的结构元素腐蚀,则访问像素的次数为a+4b1+(4+k2)b2,其中a为白色像素个数,b1为黑色内部像素个数,b2为黑色边界像素个数,且有a+b1+b2=n2按思路2实现的算法代码如下:

复制代码
Bitmap2d* Execute2()
{
    Bitmap2d* newBmp=new Bitmap2d(bmp);
    for(int j=0;j<bmp.Height();j++)
    {
        for(int i=0;i<bmp.Width();i++)
        {
            if(bmp.GetValue(i,j)==0&&HasWhiteAdjacencyPixel(i,j))
            {
                SetWindowValue(*newBmp,i,j,0);
            }
        }
    }
    return newBmp;
}
复制代码
复制代码
bool HasWhiteAdjacencyPixel(int i,int j)
{
    if(i>0&&bmp.GetValue(i-1,j)==255)
        return true;
    if(i<bmp.Width()-1&&bmp.GetValue(i+1,j)==255)
        return true;
    if(j>0&&bmp.GetValue(i,j-1)==255)
        return true;
    if(j<bmp.Height()-1&&bmp.GetValue(i,j+1)==255)
        return true;
    return false;
}
复制代码
复制代码
void SetWindowValue(Bitmap2d& bmp,int i,int j,byte v)
{
    for(size_t k=0;k<winOffsets.size();k++)
    {
        int tx=i+winOffsets[k].X;
        int ty=j+winOffsets[k].Y;
        if(!bmp.InRange(tx,ty))
            continue;
        bmp.SetValue(tx,ty,v);
    }
}
复制代码

 

基于结构元素分解的算法

  对于一些具有规则形状的结构元素,可以利用矩阵分解的原理降低计算次数,例如3*3的正方形结构元素,等价于一个3*3的矩阵,这个矩阵可以为解为{1,1,1}与{1,1,1}-1的乘积。这样使用3*3的矩阵对图像进行卷积等价于先使用{1,1,1}进行卷积,再将结果使用{1,1,1}-1进行卷积。

  由于膨胀腐蚀算法本质上属于卷积的一种特殊形式,这样,正方形结构元素的膨胀腐蚀可以使用如下的方式实现:

复制代码
Bitmap2d* Execute4()
{
    Bitmap2d* newBmp=new Bitmap2d(bmp);
    Bitmap2d* newBmp2=new Bitmap2d(bmp);
    if(this->mode==SQUARE)
    {
        winOffsets.clear();
        for (int i = 0; i < 2 * radius + 1; i++)
        {
            IntDouble t(i-radius,0);
            this->winOffsets.push_back(t);
        }
        for(int j=0;j<bmp.Height();j++)
        {
            for(int i=0;i<bmp.Width();i++)
            {
                if(HasBlackInWindow(this->bmp,i,j))
                    newBmp->SetValue(i,j,0);
                else
                    newBmp->SetValue(i,j,255);
            }
        }
        winOffsets.clear();
        for (int j = 0; j < 2 * radius + 1; j++)
        {
            IntDouble t(0,j-radius);
            this->winOffsets.push_back(t);
        }
        for(int j=0;j<newBmp->Height();j++)
        {
            for(int i=0;i<newBmp->Width();i++)
            {
                if(HasBlackInWindow(*newBmp,i,j))
                    newBmp2->SetValue(i,j,0);
                else
                    newBmp2->SetValue(i,j,255);
            }
        }
    }
    newBmp2->visit_count+=newBmp->visit_count;
    delete newBmp;
    return newBmp2;
}
复制代码

  经过测试可以知道这种方式可以大大减少像素访问次数,以k*k的结构元素腐蚀n*n的图像为例,用思路1的方法需要至少访问k2n2次像素,经过分解再处理两次只需要2kn2次访问。这个思路的详细数学原理可以参考链接

  下图是分解的方法与思路1的方法的结果对比,可以看出这两个算法的结果确实是完全等价的。

思路1 分解的方法

 

基于曼哈顿距离的算法

  上述思路1思路2可以适用于任意形状的处理窗口。还有一种基于曼哈顿距离的实现方式,来源于链接,这种方式主要是实现了基于菱形窗口的膨胀腐蚀。这里简单介绍一下曼哈顿距离,曼哈顿距离(Manhattan Distance)是种使用在几何度量空间的几何学用语,用以标明两个点在标准坐标系上的绝对轴距总和。其计算公式为:

  这个距离简单点理解就是“格子距离”,如下图所示:A到B的走格子的最少步数是4,那么AB的曼哈顿距离就是4。

  设我们需要膨胀的图像是下图左这样一个背景为0,内容为1的二值图像。假如我们能够求得所有0像素到离自己最近的1像素的距离的话,我们便做成了一张曼哈顿距离图(下图右)。曼哈顿距离图中像素标的数字代表该像素在左图中寻找最近的1的曼哈顿距离。假如这个像素在左图中本来就是1,则该像素处的曼哈顿距离为0。可以看出,01边界处的0像素的曼哈顿距离较小,而原理边界的0像素曼哈顿距离很大。

原图 原图得到的曼哈顿距离图

  对于二值图像I,若能够一定处理计算得到他的曼哈顿距离图D,则想求取他的菱形结构元素膨胀结果会非常容易。不难想到,对D进行一个阈值化既可以达到结果。若将曼哈顿图D中曼哈顿距离大于等于1与小于1的像素区分开,则等于原二值图像;若将曼哈顿距离大于等于2与小于2的像素区分开,则等价于对原二值图像进行一个尺寸为1的菱形元素膨胀;若将曼哈顿距离大于等于k(k>1)与小于k的像素区分开,则等价于对原二值图像进行一个尺寸为k的菱形元素膨胀。

  而腐蚀同样可以使用这个思路来完成,前面介绍的曼哈顿距离图是适用与膨胀的,求取的是每个0像素与距离最近的1的距离。在腐蚀的场合下,我们可以求取所有1像素与距离最近的0像素距离的曼哈顿图,这样再进行阈值化,也就完成了腐蚀操作。利用曼哈顿图的好处还体现在需要使用对很多组不同大小的结构元素对相同图像进行膨胀或腐蚀的场合。一旦计算出曼哈顿距离图,就可“一次预处理,多次复用”,预处理的开销只在初次处理产生,之后的所有操作都是阈值化的过程,而阈值化我们知道只需要width*height的访问开销。

  所以问题的关键在与如何实现对二值图像I求取其曼哈顿距离图D。这里以求取膨胀的曼哈顿距离图为例进行说明。其实我们可以利用一种类似于动态规划的思想来解决这个问题。不难发现这个问题是能够分解为规模更小并且可以复用的小型子问题的和。这基于如下的事实:

  1. 对于所有I中为1的像素,D中他们为0。因为他们自己就是1像素显然到自己最近,所以不需要走格子。
  2. 对于I中的0像素p,若其四邻域像素在D中为d0、d1、d2、d3,则D(p)=min(d0,d1,d2,d3)+1。不难看出p到离其最近的1像素的通路必然经过了其四邻域像素。所以0像素p到最近的1的像素的曼哈顿距离可以基于其四邻域的曼哈顿距离求得。

  要实现这个思路,可以使用递归,但也可以使用更加直接的方式,下面的代码使用两次双循环来求得D。首先每个像素d值默认值为最大值width+height,第一次双循环,对每一个像素实际上是考察了上方和左方的像素,经过这一次循环,其d值不一定正确,仅是能够保证每个像素处的d值是相对与上方和左方的最小值加1;但第二次双循环是逆向,从下方和右方访问像素,依次再改变之前的d值,这样就实现了d值确实为min(d0,d1,d2,d3)+1。

行序正向赋值,每个像素参考了两个父方向的d值 第二次迭代行序逆向复制,每个像素参考4个方向的d值

  采用这个思路实现的一个演示程序如下(不能跑刷新几次试试..):

  其实现的代码如下:

复制代码
class DistenceMap
{
private:
    int* data;
    int width;
    int height;
public:
    int visit_count;
    DistenceMap(int width,int height,int v)
    {
        this->data=new int[width*height];
        for(int i=0;i<width*height;i++)
            data[i]=v;

        this->width=width;
        this->height=height;
        this->visit_count=0;
    }
    ~DistenceMap()
    {
        delete[] data;
    }
    inline int GetValue(int x,int y)
    {
        visit_count++;
        return data[x+y*width];
    }
    inline void SetValue(int x,int y,int v)
    {
        visit_count++;
        data[x+y*width]=v;
    }
    inline int Width()
    {
        return width;
    }
    inline int Height()
    {
        return height;
    }
    inline int Length()
    {
        return width*height;
    }
};
复制代码
复制代码
Bitmap2d* Execute3()
{
    Bitmap2d* newBmp=new Bitmap2d(bmp);
    DistenceMap* dmap=GetDistenceMap();
    for (int i=0; i<bmp.Width(); i++)
    {
        for (int j=0; j<bmp.Height(); j++)
        {
            byte v=dmap->GetValue(i,j)<=radius?0:255;
            newBmp->SetValue(i,j,v);
        }
    }
    newBmp->visit_count+=dmap->visit_count;
    delete dmap;
    return newBmp;
}
复制代码
复制代码
DistenceMap* GetDistenceMap()
{
    DistenceMap* distenceMap=new DistenceMap(this->bmp.Width(),this->bmp.Height(),0);
    for (int i=0; i<bmp.Width(); i++)
    {
        for (int j=0; j<bmp.Height(); j++)
        {
            if (bmp.GetValue(i, j) == 0)
            {
                distenceMap->SetValue(i,j,0);
            } 
            else
            {
                distenceMap->SetValue(i,j, bmp.Width()+bmp.Height());
                if (i>0) 
                    distenceMap->SetValue(i,j,Min(distenceMap->GetValue(i,j),distenceMap->GetValue(i-1,j)+1));
                if (j>0) 
                    distenceMap->SetValue(i,j,Min(distenceMap->GetValue(i,j), distenceMap->GetValue(i,j-1)+1));
            }
        }
    }

    for (int i=bmp.Width()-1; i>=0; i--)
    {
        for (int j=bmp.Height()-1; j>=0; j--)
        {
            if (i+1<bmp.Width())
                distenceMap->SetValue(i,j,Min(distenceMap->GetValue(i,j), distenceMap->GetValue(i+1,j)+1));
            if (j+1<bmp.Height()) 
                distenceMap->SetValue(i,j,Min(distenceMap->GetValue(i,j), distenceMap->GetValue(i,j+1)+1));
        }
    }
    return distenceMap;
}
复制代码

 

总结

  本文介绍的实现方式,思路1和思路2是基本方法,其中思路2是对思路1的极大改进;矩阵分解方法适用于一些特殊形状的结构元素,其核心是把结构元素所代表的矩阵分解成两个更简单的矩阵的乘积,然后再使用这两个更简单的矩阵作为结构元素。这个思路同样能与思路1和2相配合使用;曼哈顿距离法使用一步预处理先计算出曼哈顿距离图,之后再对这个图进行阈值化,等价于使用菱形结构元素进行的膨胀腐蚀的结果,对于需要多次膨胀腐蚀的场合,这个方法非常适用。

膨胀与腐蚀算法

阅读数 8999