三、数学形态学在图像处理中的主要应用
近年来,数学形态学在图像处理方面得到了日益广泛的应用。下面主要就数学形态学在边缘检测、图像分割、图像细化以及噪声滤除等方面的应用做简要介绍。
(1)边缘检测
边缘检测是大多数图像处理必不可少的一步,提供了物体形状的重要信息。对于二值图像,边缘检测是求一个集合A的边界,记为B(A):
对于灰度图像,边缘检测是求一幅图像的形态学梯度,记为g:
数学形态学运算用于边缘检测,存在着结构元素单一的问题。它对与结构元素同方向的边缘敏感,而与其不同方向的边缘(或噪声)会被平滑掉,即边缘的方向可以由结构元素的形状确定。但如果采用对称的结构元素,又会减弱对图像边缘的方向敏感性。所以在边缘检测中,可以考虑用多方位的形态结构元素,运用不同的结构元素的逻辑组合检测出不同方向的边缘。
梁勇等人构造了8个方向的多方位形态学结构元素,应用基本形态运算,得到8个方向的边缘检测结果,再把这些结果进行归一化运算、加权求和,得到最终的图像边缘。该算法在保持图像细节特征和平滑边缘等方面,取得了较好的效果。
原图: 边界提取:
边缘检测源代码:
/********************************************************************
形态学基本操作采用二值图像
***********************************************************************/
#include<cv.h>
#include <highgui.h>
int main(){
IplImage * image,*image2,*image3;
image = cvLoadImage("E:\\image\\mapleleaf.tif", 0);
cvNamedWindow("image",1);
cvShowImage("image",image);
/*边界提取*/
image2 = cvCreateImage(cvSize(image->width, image->height),image->depth ,1);
image3 = cvCreateImage(cvSize(image->width, image->height),image->depth ,1);
int i , j ;
int left,right, up , down;
int n = 2;//窗口大小为5*5
int r,s,flag;
unsigned char * ptr, *dst;
for (i = 0 ; i< image->height; i++)
{
for (j = 0 ; j< image->width; j++)
{
//窗口设定
left = j - n;
right = j +n;
up = i -n;
down = i+n;
//窗口出界处理
if(left< 0){left = 0;}
if(right >= image->width){right = image->width-1;}
if(up< 0){up =0;}
if(down >= image->height){down = image->height -1;}
//腐蚀处理
dst = (unsigned char *)image2->imageData + image2->widthStep*i + j;
flag = 1;
for (r = up;r <= down;r++)
{
for (s = left ; s<= right; s++)
{
ptr = (unsigned char *)image->imageData + r*image->widthStep + s;
if(*ptr != 255){
flag = 0;
}
}
}
if (flag == 1)
{
*dst = 255;
}
else{
*dst = 0;
}
}
}
cvSub(image,image2,image3,0);
cvNamedWindow("image3",1);
cvShowImage("image3",image3);
cvSaveImage("mapleleafboard.bmp", image3);
cvWaitKey(0);
cvReleaseImage(&image);
cvReleaseImage(&image2);
cvReleaseImage(&image3);
cvDestroyAllWindows();
return 0 ;
}
(2)图像分割
基于数学形态学的图像分割算法是利用数学形态学变换,把复杂目标X分割成一系列互不相交的简单子集X1,X2,…,XN,即:
对目标X的分割过程可按下面的方法完成:首先求出X的最大内接“圆”X1,然后将X1从X中减去,再求X-X1的最大内接“圆”X2,…,依此类推,直到最后得到的集合为空集为止。下面以二值图像为例,介绍用数学形态学方法求解子集X1,X2,…,XN的过程。
设B为结构元素,B可以是圆、三角形、正方形等简单的几何基元,那么“简单”形状集合Xi可以用下面的公式来定义:
式中ni为一整数,用上式定义Xi分割目标,有时会产生分割过程不唯一的现象。为此可采用下面公式来定义简单集合Xi:
其中Li为一个点或一条线,当Li为点时,则与(12)式定义等价。(13)式定义的简单形状Xi可由niB沿线Li移动而产生。即将“产生器”niB的中心沿“脊骨”Li移动产生。如果niB为圆,则得到的Xi称Blum带。它具有一些特殊的性质,如Xi的边界是光滑的,Xi的最大圆与其边界相切,Xi的脊骨与产生器都是唯一的等等。
有了简单形状集合Xi的定义,则目标X可按下面方法分割。首先按式(14)求出X的最大内切结构元素Xi:
数学形态学用于图像分割的缺点是对边界噪声敏感。为了改善这一问题,刘志敏等人提出了基于图像最大内切圆的数学形态学形状描述图像分割算法和基于目标最小闭包结构元素的数学形态学形状描述图像分割算法,并使用该算法对二值图像进行了分割,取得了较好的效果。邓世伟等人提出一种基于数学形态学的深度图像分割算法。作者首先利用形态学算子获得分别含有阶跃边缘与屋脊边缘的凸脊和凹谷图像,然后利用控制区域生长过程得到最终的分割结果。与传统方法相比,该方法速度快,抗噪性能好。
(3) 形态骨架提取
形态骨架描述了物体的形状和方向信息。它具有平移不变性、逆扩张性和等幂性等性质,是一种有效的形状描述方法。二值图像A的形态骨架可以通过选定合适的结构元素B,对A进行连续腐蚀和开启运算来求取,设S(A)代表A的骨架,定义为:
蒋刚毅等人运用数学形态学方法,对交通标志的内核形状提取形态骨架函数,将其作为用于模式匹配的形状特征。A的形态骨架函数SKF(A)表示为:
SKF(X)中值较大的点对应大的n,并代表了形态骨架的主要成分,即表达了形状的主体结构;而SKF(X)中值较小的点对应小的n,是形态骨架的细节成分,与形状的边缘信息相联系。
形态骨架函数完整简洁地表达了形态骨架的所有信息,因此,根据形态骨架函数的模式匹配能够实现对不同形状物体的识别。算法具有位移不变性,因而使识别更具稳健性。
原图: 提取骨架: DNA原图: 骨架:
骨架提取原代码:
/************************************************************************/
/* 骨架提取*/
/************************************************************************/
#include<cv.h>
#include <highgui.h>
int main(){
IplImage* image = cvLoadImage("E:\\image\\bone.tif",0);
cvNamedWindow("image",1);
cvNamedWindow("image2",1);
cvNamedWindow("image3",1);
cvNamedWindow("image4",1);
cvNamedWindow("image5",1);
cvNamedWindow("image6",1);
cvNamedWindow("result",1);
cvShowImage("image", image);
//cvWaitKey(0);
//当前图片image2 当前被腐蚀后的图片image3 被腐蚀开操作之后的图片 image5
//image4 作为开操作的中间值 作差后的图片image6 并之后的图片result
IplImage *image2, *image3 , *image4, *image5,*image6,*result;
image2 = cvCreateImage(cvSize(image->width,image->height),IPL_DEPTH_8U,1);
image3 = cvCreateImage(cvSize(image->width,image->height),IPL_DEPTH_8U,1);
image5 = cvCreateImage(cvSize(image->width,image->height),IPL_DEPTH_8U,1);
image4 = cvCreateImage(cvSize(image->width,image->height),IPL_DEPTH_8U,1);
image6 = cvCreateImage(cvSize(image->width,image->height),IPL_DEPTH_8U,1);
result = cvCreateImage(cvSize(image->width,image->height),IPL_DEPTH_8U,1);
//循环标志 flag
bool flag = true;
//腐蚀判断标志 flag2
bool flag2 = true;
int i,j,r,s;
unsigned char * ptr,*dst;
unsigned char B[9] = {255 ,255,255,255,255,255,255,255,255};
//对result进行赋值 全部为0
for (i = 0 ; i< result->height; i++)
{
for (j = 0 ; j< result->width ; j++)
{
dst = (unsigned char *)(result->imageData + result->widthStep *i +j);
*dst = 0;
}
}
image2 = cvCloneImage(image);
cvShowImage("image2", image2);
//cvWaitKey(0);
while (flag)
{
//flag = false;
cvShowImage("image2",image2);
//进行腐蚀操作,开环操作 作差 并
for (i = 0 ; i< image3->height; i++)
{
for (j = 0 ; j< image3->width ; j++)
{
dst = (unsigned char *)(image3->imageData + i*image3->widthStep + j);
if ((i == 0 )|| (j == 0) ||( i == image->height -1 ) || (j == image->width -1 ))
{
*dst = 0;
//break;
}
else{
flag2 = true;
for (r = i-1 ; r<= i+1 ; r++)
{
for (s = j -1 ; s<= j+1 ; s++)
{
ptr = (unsigned char *)(image2->imageData + r*image2->widthStep + j);
if(*ptr != 255){
flag2 =false;
}
}
}
if (flag2)
{
*dst = 255;
}
else {*dst = 0;}
}
}
}
cvShowImage("image3",image3);
//开操作 先腐蚀 后膨胀
for (i = 0 ; i< image4->height; i++)
{
for (j = 0 ; j< image4->width ; j++)
{
dst = (unsigned char *)(image4->imageData + i*image4->widthStep + j);
if ((i == 0 )|| (j == 0) ||( i == image->height -1 ) || (j == image->width -1 ))
{
*dst = 0;
//break;
}
else{
flag2 = true;
for (r = i-1 ; r<= i+1 ; r++)
{
for (s = j -1 ; s<= j+1 ; s++)
{
ptr = (unsigned char *)(image3->imageData + r*image3->widthStep + s);
if(*ptr != 255){
flag2 =false;
}
}
}
if (flag2)
{
*dst = 255;
}
else {*dst = 0;}
}
}
}
cvShowImage("image4",image4);
//膨胀操作
for (i = 0 ; i< image5->height; i++)
{
for (j = 0 ; j< image5->width ; j++)
{
dst = (unsigned char *)(image5->imageData + i*image5->widthStep + j);
if ((i == 0 )|| (j == 0) ||( i == image5->height -1 ) || (j == image5->width -1 ))
{
*dst = 0;
//break;
}
else{
flag2 = false;
for (r = i-1 ; r<= i+1 ; r++)
{
for (s = j -1 ; s<= j+1 ; s++)
{
ptr = (unsigned char *)(image4->imageData + r*image4->widthStep + s);
if(*ptr == 255){
flag2 = true;
}
}
}
if (flag2)
{
*dst = 255;
}
else {*dst = 0;}
}
}
}
cvShowImage("image5",image5);
//作差
cvSub(image3,image5,image6,0);
//并运算
for (i = 0 ; i< result->height; i++)
{
for (j = 0 ; j< result->width ; j++)
{
dst = (unsigned char *)(result->imageData + result->widthStep *i +j);
ptr = (unsigned char *)(image6->imageData + image6->widthStep *i +j);
if (*ptr == 255)
{
*dst = 255;
}
}
}
cvShowImage("image6",image6);
//将腐蚀后的图像复制给当前图像image2
image2 = cvCloneImage(image3);
//循环标志判定
flag = false;
for (i = 0 ; i< image2->height; i++)
{
for (j = 0 ; j< image2->width ; j++)
{
ptr = (unsigned char *)(image2->imageData + image2->widthStep *i +j);
if (*ptr == 255)
{
flag = true;
}
}
}
cvShowImage("result", result);
cvWaitKey(40);
}
cvShowImage("image2", image2);
cvShowImage("image3", image3);
//cvShowImage("image4", image4);
//cvShowImage("image5", image5);
//cvShowImage("image6", image6);
cvShowImage("result",result);
cvSaveImage("E:\\image\\bonegujia.bmp",result);
cvWaitKey(0);
cvReleaseImage(&image);
cvDestroyAllWindows();
return 0;
}
形态学算法对于提取交叉的物体,会产生断裂,一般会在提取之后紧跟连接操作。
原图: 提取的角点:
角点提取源代码:采用击中击不中
/********************************************************************
形态学基本操作采用二值图像 形态学方法在边界获取和形状检测中有很多应用 寻找角点
***********************************************************************/
#include<cv.h>
#include <highgui.h>
int main(){
IplImage * image,*image2,*image3,*image4,*result;
image = cvLoadImage("E:\\image\\jiaodian.bmp", 0);
cvNamedWindow("image",1);
cvShowImage("image",image);
/*击中击不中变换 寻找角点*/
/************************************************************************/
/* 最终结果为(79, 85 )和 (129 ,134)两个符合要求的角点 在结果图像中可以看到两个白点*/
/************************************************************************/
image2 = cvCreateImage(cvSize(image->width, image->height),image->depth ,1);
image3 = cvCreateImage(cvSize(image->width, image->height),image->depth ,1);
image4 = cvCreateImage(cvSize(image->width, image->height),image->depth ,1);
result = cvCreateImage(cvSize(image->width, image->height),image->depth ,1);
int i,j,r,s ;
unsigned char B1[9] ={ 255,255,255,255,0,0,255,0,0};
unsigned char B2[9] = {0,0,0,0,255,255,0,255,255 };
int flag;
//iamge2 是对image进行求补的结果
unsigned char *ptr, *dst,*ptr2;
for (i = 0 ; i< image2->height; i++)
{
for (j = 0 ; j< image2->width; j++)
{
ptr = (unsigned char*)image->imageData + i * image->widthStep+ j;
dst = (unsigned char*)image2->imageData + i*image2->widthStep + j;
*dst = 255- (*ptr);
}
}
//对源图像进行腐蚀
for (i = 0 ; i< image3->height; i++)
{
for (j = 0 ; j< image3->width; j++)
{
flag = 1;
dst = (unsigned char*)image3->imageData + i*image3->widthStep +j;
//边界判断
if (( i == 0) || (j == 0) || (i == image3->height-1) || (j == image3->width -1 ))
{
*dst = 0;
}
else{
for (r = i -1 ; r<= i+1; r++)
{
for (s = j-1 ; s <= j +1 ; s++)
{
ptr = (unsigned char*)image->imageData + r * image->widthStep+ s;
if (*ptr != B1[3*(r-i+1) + s-j+1])
{
flag = 0;
break;
}
}
}
if (flag == 1)
{
ptr = (unsigned char*)image->imageData + i * image->widthStep+ j;
*dst = 255;
}
else{
*dst = 0;
}
}
}
}
//显示腐蚀结果
for (i = 0 ; i< image2->height; i++)
{
for (j = 0 ; j< image2->width; j++)
{
ptr = (unsigned char*)image->imageData + i * image->widthStep+ j;
ptr2 = (unsigned char*)image3->imageData + i*image3->widthStep + j;
if (*ptr2 == 255)
{
printf("x, %d y: %d %d\n",j, i, *ptr2 - *ptr);
}
}
}
//对补图像进行腐蚀
for (i = 0 ; i< image4->height; i++)
{
for (j = 0 ; j< image4->width; j++)
{
flag = 1;
dst = (unsigned char*)image4->imageData + i*image4->widthStep +j;
//边界判断
if (( i == 0) || (j == 0) || (i == image4->height-1) || (j == image4->width -1 ))
{
*dst = 0;
}
else{
for (r = i -1 ; r<= i+1; r++)
{
for (s = j-1 ; s <= j +1 ; s++)
{
ptr = (unsigned char*)image2->imageData + r * image2->widthStep+ s;
if ((*ptr) != B2[3*(r- i+1) + s-j +1])
{
flag = 0;
}
}
}
if (flag == 1)
{
ptr = (unsigned char*)image2->imageData + i * image2->widthStep+ j;
*dst = 255;
}
else{
*dst = 0;
}
}
}
}
//显示腐蚀结果
for (i = 0 ; i< image4->height; i++)
{
for (j = 0 ; j< image4->width; j++)
{
ptr = (unsigned char*)image2->imageData + i * image->widthStep+ j;
ptr2 = (unsigned char*)image4->imageData + i*image3->widthStep + j;
if (*ptr2 == 255)
{
printf("x, %d y: %d %d\n",j,i, *ptr2 - *ptr);
}
}
}
//二者求交集
for (i = 0 ; i< result->height; i++)
{
for (j = 0 ; j< result->width; j++)
{
ptr = (unsigned char *)image3->imageData + image3->widthStep * i + j;
ptr2 = (unsigned char *)image4->imageData + image4->widthStep * i + j;
dst = (unsigned char *)result->imageData + result->widthStep*i + j;
if (((*ptr) == 255) && ((*ptr2) == 255))
{
*dst = 255;
printf("x : %d\n" , j);
printf("y : %d\n" , i);
}
else{
*dst = 0;
}
}
}
cvNamedWindow("image3",1);
cvNamedWindow("image4",1);
cvNamedWindow("result",1);
cvShowImage("image3",image3);
cvShowImage("image4",image4);
cvShowImage("result",result);
//cvSaveImage("mapleleafboard.bmp", image3);
cvWaitKey(0);
cvReleaseImage(&image);
cvReleaseImage(&image2);
cvReleaseImage(&image3);
cvDestroyAllWindows();
return 0 ;
}
(4)噪声滤除
对图像中的噪声进行滤除是图像预处理中不可缺少的操作。将开启和闭合运算结合起来可构成形态学噪声滤除器。
对于二值图像,噪声表现为目标周围的噪声块和目标内部的噪声孔。用结构元素B对集合A进行开启操作,就可以将目标周围的噪声块消除掉;用B对A进行闭合操作,则可以将目标内部的噪声孔消除掉。该方法中,对结构元素的选取相当重要,它应当比所有的噪声孔和噪声块都要大。
对于灰度图像,滤除噪声就是进行形态学平滑。实际中常用开启运算消除与结构元素相比尺寸较小的亮细节,而保持图像整体灰度值和大的亮区域基本不变;用闭合运算消除与结构元素相比尺寸较小的暗细节,而保持图像整体灰度值和大的暗区域基本不变。将这两种操作综合起来可达到滤除亮区和暗区中各类噪声的效果。同样的,结构元素的选取也是个重要问题。