基于MATLAB小波变换的医学图像分割的研究

一 图像分割技术的现状和发展

研究图像分割算法已有几十年的历史，一直受到人们的高度重视。国内外发表了许多关于图像分割原理和方法的论文，但没有一种分割方法适用于所有图像分割处理。传统的图像分割方法存在不足，不能满足人们的要求，给进一步的图像分析和理解带来了困难。随着计算机技术的快速发展，以及相关技术的发展和成熟，结合图像增强等技术，可以在计算机上实现图像分割处理。

最重要的技术是图像分割技术，将特定区域与图像中的其他部分分离提取。图像分割的方法有很多，包括阈值分割、边界分割、区域提取、特定理论工具的分割等。早在1965年，就有人提出检测边缘算子，边缘检测就产生了许多经典算法。越来越多的学者开始将数学形态学、模糊理论、遗传算法理论、分形理论和小波变换理论应用于图像分割，产生了结合特定数学方法和特殊图像分割的先进图像分割技术。特别是近年来，小波理论的快速发展为图像处理带来了新的理论和方法。小波变换具有良好的局部特性。当小波函数尺度较大时，抗噪能力较强。当小波函数尺度较小时，提取图像细节的能力较强，可以很好地解决抑制噪声与提取图像边缘细节的矛盾。

二、 图像分割的主要研究方法

2.1 边缘检测法

边缘检测通常是图像分析和理解的第一步。边缘检测是人们研究的一种方法，它通过检测图像中不同区域的边缘来分割图像。边缘检测的本质是用某种算法提取图像中对象与背景问题的交界线。我们将边缘定义为图像中灰度变化迅速的区域边界。图像灰度分布的梯度可以反映图像灰度的变化，因此我们可以使用局部图像微分技术获得边缘检测算子。经典的边缘检测方法是通过对原始图像中像素的小邻域构造边缘检测算子来检测边缘。

2.2 区域提取法

区域提取有两种基本形式：一种是从单个像素开始，逐渐合并形成所需的分割区域；另一种是从整个图片到所需的分割区域。这两种基本形式的结合通常用于实践中。根据上述两种基本形式，区域提取法可分为区域生长法和分裂合并法。区域生长法的基本思想是将具有相似性质的像素结合起来形成一个区域。具体做法是先给定图像中要分割的目标物体中的小块或种子区域，然后在种子区域的基础上，不断以一定的规则添加周围的像素点，最终将代表物体的所有像素点结合成一个区域。该方法的关键是选择适当的生长或类似标准。生长标准一般可分为区域灰度差标准、区域灰度分布统计性质标准和区域形状标准三种。分裂合并法是将图像分割成许多一致性强的小区域，然后按照一定的规则将小区域整合成大区域，达到图像分割的目的。区域提取法的缺点是图像分割过多，因此近年来对该方法的研究较少。

2.3 阈值分割法

灰度图像的阈值分割是首先确定图像灰度值范围内的灰度阈值，然后将图像中每个像素的灰度值与该阈值进行比较，并根据比较结果将相应的像素分为两类。这两种像素通常分为图像的两个区域，以达到分割的目的。阈值分割算法主要有两个步骤：

(1)确定所需阈值；

(2)将分割阈值与像素值进行比较，划分像素。

可见，确定最优阈值是分割的关键。现有的大部分算法都是集中在阈值确定的研究上。根据图像本身的特点，阈值分割方法可分为单阈值分割方法和多阈值分割方法：或基于像素值的阈值分割方法、基于区域性质的阈值分割方法和基于坐标位置的阈值分割方法。也可分为直方图和直方图转换法、最大空间方差法、最小误差法和均匀误差法、共生矩阵法、最大熵法、简单统计法和局部特征法、概率松弛法、模糊集合法等。

2.4 结合特定理论工具的分割方法

近年来，随着各学科许多新理论和方法的提出，人们还提出了许多结合特定理论工具的分割方法，如基于数学形态学、神经网络、信息理论、模糊集合和逻辑、小波分析和变换、遗传算法等。基于小波分析和变换的分割方法是利用新数学工具小波变换来分割图像的一种方法，现在也是一种非常新的方法。小波变换是一种多尺度、多通道分析工具，更适合多尺度边缘检测图像。例如，高斯函数的一级和二级导数可以用作小波函数Mallat算法分解小波，然后根据马尔算子进行多尺度边缘检测，这里可以控制观察距离的调焦。可以选择检测边缘的细节来改变高斯函数的标准差。小波变换计算复杂度低，抗噪能力强。理论证明，以零点为对称点的对称二进波适用于检测屋顶状边缘，而以零点为对称点的对称二进波适用于检测阶跃状边缘。近年来，多通道小波也开始用于边缘检测。此外，多尺度边缘也可以通过计算和估计图像奇异度来提取，一些边缘类型可以通过正交小波基的小波变换来区分。

三、 预处理图像分割

3.1 图像平滑

平滑图像的目的是减少图像噪声。图像噪声来自系统的外部干扰，如电磁波或电源串入系统引起的外部噪声，以及相机的热噪声、电气机械运动引起的抖动噪声等内部干扰。因此，去除噪声和恢复原始图像是图像处理的重要组成部分。噪声主要来自以下三个方面：

（1）光电子噪声：主要由光的统计本质和图像传感器的光电转换过程引起（如光电管的光量子噪声和电子起伏噪声）；

(2)电子噪声:主要来自电子元件(如电阻引起的热噪声)；

(3)光学噪声:主要由光学现象产生(如胶片颗粒结构产生的颗粒噪声)；

图像在生成和传输过程中受到这些噪声的干扰和影响，使图像处理结果恶化。因此，抑制或消除这些噪声以提高图像质量是图像处理过程中的一个重要预处理，也称为图像的平滑滤波过程。

图2-3显示了图像中值滤波前后的效果比较，图2-3（a）是含噪音的原图，图2-3（b）采用中值滤波处理的图像，滤波窗为3×3.由此可见，中值滤波后的图像不仅滤去了椒盐的噪音，而且保护了边缘。

(a)带噪声图像 (b)消噪后图像

图2-3 噪声图像与中值滤波后图像进行比较

3.2 灰度调整

在成像过程中，扫描系统和光电转换系统中的许多因素，如光强度、光敏度、光学系统不均匀稳定等诸多因素都会导致图像亮度分布不均匀，导致部分亮度和部分暗度。灰度调整是在图像采集系统中纠正图像像素，使整个图像成像均匀。

3.2.1 灰度调整原理

灰度调整是图像增强的重要手段之一，可以增加图像的动态范围，扩大图像对比度，使图像清晰明显。

在曝光不足或过度的情况下，图像灰度可能限制在一个非常小的范围内。此时，显示器上会看到一个模糊的图像，似乎没有灰度水平。线性灰度调整对图像的每个像素灰度进行线性拉伸，将有效提高图像的视觉效果。

3.2.2 灰度调节效果分析

（a）灰度调整前 (b)灰度调整后

(c)原始图像直方图 （d）调整后的直方图

图2-4 灰度调整前后直方图比较

从图2-4可以看出(b)视觉效果较(a)显然，灰度调整前后直方图的比较可以看出，调整后直方图(d)去除原始直方图(c)灰度调整后，图像明显清晰。

四、 基于阈值的图像分割技术

当非灰度图像转换为灰度图像时，图像中目标区域的灰度值会有所不同。如果图像的灰度方形图具有明显的双峰值或多峰值特征，则可以使用阈值法获得最佳阈值，然后合理划分图像。

4.1 阈值分割原理

阈值图像分割是最基本的图像分割方法之一。经过半个多世纪的研究，已经提取了大量的算法。其基本原理是选择灰度图像范围内的一个或多个灰度阈值，然后将图像中每个像素的灰度值与阈值进行比较，并根据比较结果将图像中的对应像素分为两类或两类以上，将图像分为不重叠区域集合，达到图像分割的目的。

在分割阈值图像时，通常需要对图像进行一定的模型假设。尽可能多地使用图像模型来了解图像的几个不同区域。基于图像分割模型的假设：目标或背景中相邻像素之间的灰度值相似，但不同目标或背景的像素在灰度上存在差异。设置原始图像为f(x,y),按照一定的标准f(x,y)在分割过程中找到特征值是阈值T，或者找到合适的区域空间Ω，将图像分割成两个部分，分割后的图像为

（3-1）

对于各种阈值，分割图像可以表示为：一组分割阈值，分割后对应不同区域的图像灰度值，K分割后的区域或目的

标数，。

无论是单阈值分割还是多阈值分割，都选择合理的阈值，以确定图像中的每个像素点应属于目标区域或背景区域，从而产生相应的二值图像。

4.2 图像分割法

阈值分割是设置一个门限（阈值）。如果图像灰度值大于或等于（或小于或等于）门限值，则分为另一类，其中一类为背景，另一类为目标。

4.2.1 图像二值化

基于区域划分的主要方法是二值化。二值化方法对由多个实体和一个对比较强的背景图像组成的场景图像特别有效。二值化方法一般比较快，而且使每个分割出来的物体都具有闭合和连通的边界。图像二值化后信息丢失很严重，由此得到的边界轮廓可能会不精确。因此，可以用速度较快的二值化方法来获得一个关于图像分割结果的较粗略的描述。

4.2.2 双峰法

在一些简单的图像中，对象物的灰度分布较有规律，背景和各个对象物在图像的灰度直方图中各自形成一个波峰，即区域和波峰一一对应。由于每个波峰间形成一个波谷，因为选择双峰间的波谷处所对应的灰度值为阈值，即可将两个区域分离。以此类推，可以在图像背景中分理出各类有意义的区域。

（a）原始图像 （b）原始图像直方图

（c）阈值=25分割图像 （d）阈值=40分割图像

图3-2 二值化双峰分割

图3-2为两个简单阈值分割图，双峰法比较简单，在可能的情况下常常作为首选的阈值确定法，但是图像的灰度直方图形状随着对象、图像输入系统、输入环境等因素的不同而千差万别，当出现双峰间的波谷平坦、各区域直方图的波形重叠等情况时，用双峰法难以确定阈值，必须寻求其他方法，实现自动选择适宜阈值要求。

4.2.3 最大方差自动取阈值（自适应二值化）

图像灰度直方图的形状是多变的，有双峰但是无明显低谷或者是双峰与低谷都不明显，而且两个区域的面积比也难以确定的情况常常出现，采用最大方差自动取阈值往往能得到较为满意的结果。

图像灰度级的集合设为S=（1,2,3,…,i, …L）, 灰度级为i的像素数设为ni，则图像的全部像素数为

(3-3)

将其标准化后，像素数为，其中，i∈S，pi≥0， （3-4）

设有某一图像灰度直方图，t为分离两区域的阈值。由直方图统计可被t分离后的区域1、区域2占整图像的面积比以及整幅图像、区域1、区域2的平均灰度为：

区域1的面积比： ；区域2的面积比 （3-5）

或者

整幅图像平均灰度 ; 区域1的平均灰度；

区域2的平均灰度 (3-6)

式中，G为图像的灰度级数。

整图像平均灰度与区域1、区域2平均灰度值之间的关系为

（3-7）

同一区域常常具有灰度相似特性，而不同区域之间则表现为明显的灰度差异，当被阈值t分离的两个区域之间灰度差较大时，两个区域的平均灰度u1,u2与整图像平均灰度u之差也较大，区域间的方差就是描述这种差异的有效参数，其表达式为：

（3-8）

式中，表示了图像被阈值t分割后的两个阈值之间的方差。显然不同的t值，就会得到不同的区域方差，也就是说，区域方差、区域1均值、区域2均值、区域面积比、区域面积比都是阈值t的函数，因此式（3-8）可写为：

（3-9）

经数学推导，区域间的方差可表示为：

（3-10）

被分割的两区域间的方差达最大时，被认为是两区域的最佳分离状态，由此确定阈定值T：，以最大方差决定阈值不需要认为设定其他参数，是一种自动选择阈值的方法，它不仅适用于两区域的单阈值选择，也可以扩展到多区域的多阈值选择中去。

(a)原始图像 (b)最大方差法分割后图像

图3-11 最大方差自动取阈值法

该方法将图像分成两个类，当类间方差与类内方差的分离度最大时即为最佳阈值．由图3-11表明，该方法能够准确而快速地对图像进行二值化，特别是当对象物和背景的灰度值的差具有一定大小的时候，效果更明显。

五、 基于小波图像阈值分割技术

小波变换是近年来得到广泛应用的数学工具,与傅里叶变换、窗口傅里叶变换相比,小波变换是空间(时间) 和频率的局域变换，能有效地从信号中提取信息。

本课题利用小波变换对含噪图像的直方图进行多尺度分解,先在较大的尺度下找出图像分割阈值的粗略值,然后逐渐减小尺度,精确定位分割阈值，算法采用MATLAB 编程仿真。基于小波变换的阈值法图像分割技术则能够有效地避免噪声的影响。该方法的基本思想是首先由二进制小波变换将图像的直方图分解为不同层次的小波系数, 然后依据给定的分割准则和小波系数选择阈值门限, 最后利用阈值标出图像分割的区域。整个分割过程是从粗到细, 由尺度变化来控制, 即起始分割由粗略的L2（R）子空间上投影的直方图来实现, 如果分割不理想, 则利用直方图在精细的子空间上的小波系数逐步细化图像分割。

六、效果示例

（a）小波分解

（b）小波分割 （c）传统分割

图4-10小波阈值分割

传统阈值分割与小波阈值分割方法的比较：

一、全局二值化方法由于采用的是用一个固定门限值来分割，因此门限值的选取十分重要。虽然众多学者提出了许多种选取“最优门限”的方法，但这种分割方法在光照不均匀和需要提取多个复杂特征的物体的时候难以获得较理想的效果。自适应二值化方法由于采用了自动取阈值的方法，避免了采用固定阈值的弊病，但是图像分割后局限性太大，效果不佳。

二、基于灰度直方图小波变换的多阈值分割，在小尺度下受噪声影响较大，但对阈值的定位比较准；大尺度下受噪声影响较小，可以找到确定阈值。峰与峰之间的谷点直接在大尺度下进行，既克服了噪声的影响，又有效的选取到最优阈值。对图像分割的效果明显优于前一种方法。

三、信噪比是信号与噪声的功率谱之比，但通常功率谱难以计算，有一种方法可以近似估计图象信噪比，即信号与噪声的方差之比。首先计算图象所有象素的局部方差，将局部方差的最大值认为是信号方差，最小值是噪声方差，求出它们的比值。信噪比越大，说明混在信号里的噪声越小，信号质量越好。

峰值信噪比一般是用于最大值信号和背景噪音之间的一个工程项目。通常在经过影像压缩之后，输出的影像通常都会有某种程度与原始影像不一样。为了衡量经过处理后的影像品质，我们通常会参考PSNR值来认定某个处理程序够不够令人满意。PSNR值越大，就代表失真越少。

通常一幅图像分割结果的好与坏，以人的主观判断作为评价标准[13]，也就是说是人的视觉决定了分割结果的优良，这样就导致了由于人的视觉差异对图像分割好坏评价的不统一，所以对不同分割方法的结果做一个定量的、定性的评价也是必要的且有意义的。

七、程序代码参考

1、双峰法
tic;
I=imread('医学1.bmp');
A=I;%A=rgb2gray(I);
figure(1);
subplot(1,3,1);
B=imhist(A);
plot(B);
title('直方图');
subplot(1,3,2);
imshow(I);
title('原始图像');
[m,n]=size(I);
for i=1:m
    for j=1:n
        if(I(i,j)<50)
            I(i,j)=255;
          end
    end
end
subplot(1,3,3);
imshow(I);
imwrite(I,'I.bmp');
title('阈值40分割图像');
toc;
%图像的MSE、SNR、PSNR的计算
B=I;
[M,N]=size(A);
a=double(A);b=double(B);
sum1=0;
mseValue=sum1/(M*N);%计算均方误差MSE
sum2=0;
for i=1:M;
    for j=1:N;
        sum2=sum2+(a(i,j))^2;
    end;
end;
P=sum2;
snrValue=10*log10(P/mseValue);%计算信噪比SNR
psnrValue=10*log10(255^2/mseValue);%计算峰值信噪比PSNR
2、最大类间差法
tic;
x=imread('医学1.bmp'); 
subplot(1,2,1);
imshow(x) ;title('原始图像');
count=imhist(x); 
 [m,n]=size(x); 
N=m*n; 
L=256; 
count=count/N;  
for i=1:L 
    if count(i)~=0 
        st=i-1; 
           end 
end 
for i=L:-1:1 
    if count(i)~=0 
        nd=i-1; 
            end 
end 
f=count(st+1:nd+1);  %f是每个灰度出现的概率 
u=0; 
for i=1:q 
    u=u+f(i)*(p+i-1);  %u是像素的平均值  
    ua(i)=u;           %ua（i）是前i个像素的平均灰度值 
end; 
 
for i=1:q 
    w(i)=sum(f(1:i));  %w（i）是前i个像素的累加概率 
end; 
 
d=(u*w-ua).^2./(w.*(1-w)); 
[y,tp]=max(d);  %可以取出数组的最大值及取最大值的点 
th=tp+p; 
 
for i=1:m 
    for j=1:n 
        if x(i,j)>th 
            x(i,j)=0; 
        else 
            x(i,j)=255; 
        end 
    end 
end  
subplot(1,2,2);
imshow(x); title('最大类方差法');
toc;
%图像的MSE、SNR、PSNR的计算
A=imread('医学1.bmp');
B=x;
[M,N]=size(A);
a=double(A);b=double(B);
sum1=0;
for i=1:M;
    for j=1:N;
        sum1=sum1+(a(i,j)-b(i,j))^2;
    end;
end;
mseValue=sum1/(M*N);%计算均方误差MSE
sum2=0;
for i=1:M;
    for j=1:N;
        sum2=sum2+(a(i,j))^2;
    end;
end;
P=sum2;
snrValue=10*log10(P/mseValue);%计算信噪比SNR
psnrValue=10*log10(255^2/mseValue);%计算峰值信噪比PSNR
3、小波分割
clear; 
tic;
imgname='医学1.bmp'; 
wavename='haar';   %用haar小波来分解
mode='per'; 
title('原图');
r=x;  %r=rgb2gray(x) 三维图像变成二维图像               
xx=histeq(r); %直方图均衡化
subplot(2,2,2);imshow(xx); 
title('增强后的图像'); 
deccof=struct('ca',[],'ch',[],'cv',[],'cd',[]); %三层分解
reccof=struct('rx',[]); 
sx=size(xx); 
nbcol=size(map,1); 
dx=xx; 
deccof(1).ca=xx; 
[deccof(2).ca,deccof(2).ch,deccof(2).cv,deccof(2).cd]=dwt2(dx,wavename,'mode',mode); %使用小波函数对图像dx进行二维小波变换
figure(2);imshow([deccof(2).ca/255,deccof(2).ch/255;deccof(2).cv/255,deccof(2).cd/255;]); 
title('小波分解'); 
am=deccof(2).ca/255; 
figure(1); 
subplot(2,2,3);imshow(am); 
title('近似分量'); 
%找阈值 
L=size(am,1); 
W=size(am,2); 
N=L*W; 
[Count Ret]=imhist(am); 
Pi=Count'./N; 
g=[]; 
for t=0:255 
    w0=sum(Pi(1:t+1)); 
    w1=1-w0; 
    mu0=sum(Pi(1:t+1).*(0:t))/w0; 
    mu1=sum(Pi(t+2:256).*(t+1:255))/w1; 
    mu=sum(Pi(1:256).*(0:255)); 
    g=[g w0*w1*(mu0-mu1)^2/(w1*(mu0-mu)^2+w0*(mu1-mu)^2)]; 
end 
[Ret1 T ]=max(g); 
T1=T; 
T=num2str(T-1); 
disp(['最佳灰度threhold:   ' T]); 
%分割 
I1=im2bw(am,T1/255); 
figure(1),subplot(2,2,4);imshow(I1); 
title('阈值分割后的图像'); 
u=idwt2(I1,deccof(2).ch/127,deccof(2).cv/127,deccof(2).cd/127,wavename,'mode',mode); 
figure(3),imshow(u);title('小波分割图像') ;
toc;
 
%图像的MSE、SNR、PSNR的计算

A=imread('医学1.bmp');
B=u;
[M,N]=size(A);
a=double(A);b=double(B);
sum1=0;
P=sum2;
snr=10*log10(P/mseValue);%计算信噪比SNR
psnr=10*log10(255^2/mseValue);%计算峰值信噪比PSNR