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

       图像分割是一种重要的图像分析技术。对图像分割的研究一直是图像技术研究中的热点和焦点。医学图像分割是图像分割的一个重要应用领域，也是一个经典难题，至今已有上千种分割方法，既有经典的方法也有结合新兴理论的方法。

       本论文首先介绍了双峰法以及最大类方差自动阈值法，然后重点介绍一种基于小波变换的图像分割方法,该方法先对图像的灰度直方图进行小波多尺度变换,然后从较大的尺度系数到较小的尺度系数逐步定位出灰度阈值。最后，对这几种算法的分割效果进行了比较。实验结果表明, 本设计能够实时稳定的对目标分割提取，分割效果良好。

       医学图像分割是医学图像处理中的一个经典难题。图像分割能够自动或半自动描绘出医学图像中的解剖结构和其它感兴趣的区域，从而有助于医学诊断。

1.1 图像分割技术的现状和发展情况

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

其中最主要的技术是图像分割技术，从图像中将某个特定区域与其它部分进行分离并提取出来的处理。图像分割的方法有许多种，有阈值分割方法，边界分割方法，区域提取方法，结合特定理论工具的分割方法等。早在1965年就有人提出检测边缘算子，边缘检测已产生不少经典算法。越来越多的学者开始将数学形态学、模糊理论、遗传算法理论、分形理论和小波变换理论等研究成果运用到图像分割中，产生了结合特定数学方法和针对特殊图像分割的先进图像分割技术。尤其是近年来迅速发展起来的小波理论为图像处理带来了新的理论和方法。小波变换具有良好局部特性，当小波函数尺度较大时，抗噪声的能力强，当小波函数尺度较小时，提取图像细节的能力强，这样就可以很好地解决抑制噪声和提取图像边缘细节之间的矛盾。

1.2 图像分割主要研究方法

图像分割是图像处理中的一项关键技术，自20世纪70年代起一直受到人们的高度重视，至今已提出了上千种各种类型的分割算法，现提出的分割算法大都是针对具体问题的，并没有一种适合于所有图像的通用分割算法，而且近年来每年都有上百篇相关研究报道发表。然而，还没有制定出选择合适分割算法的标准，这给图像分割技术的应用带来许多实际问题。因此，对图像分割的研究还在不断深入之中，是目前图像处理中研究的热点之一。

图像分割在图像工程中的位置它起着承上启下的作用，可以认为是介于低层次处理和高层次处理的中间层间。最近几年又出现了许多新思路、新方法、或改进算法。下面对一些经典传统方法作简要的概述。

多年来人们对图像分割提出了不同的解释和表述，借助集合概念对图像分割可给出如下定义：令集合R代表整个图像区域，对R的图像分割可以看做是将R分成N个满足以下条件的非空子集R1，R2，R3，…，RN；

（1）在分割结果中，每个区域的像素有着相同的特性；

（2）在分割结果中，不同子区域具有不同的特性，并且它们没有公共特性；

（3）分割的所有子区域的并集就是原来的图像；

（4）各个子集是连通的区域；

图像分割是把图像分割成若干个特定的、具有独特性质的区域并提取出感兴趣目标的技术和过程，这些特性可以是像素的灰度、颜色、纹理等提取的目标可以是对应的单个区域，也可以是对应的多个区域。图像分割方法有许多种分类方式，在这里将分割方法概括为四类：(1)边缘检测方法(2)区域提取方法(3)阈值分割方法(4)结合特定理论工具的分割方法。下面就这些方法展开介绍。

1.2.1 边缘检测法

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

1.2.2 区域提取法

区域提取法有两种基本形式：一种是从单个像素出发，逐渐合并以形成所需的分割区域；另一种是从全图出发，逐渐分裂切割至所需的分割区域。在实际中使用的通常是这两种基本形式的结合。根据以上两种基本形式，区域提取法可以分为区域生长法和分裂合并法。区域生长法的基本思想是将具有相似性质的像素合起来构成区域，具体做法是先给定图像中要分割的目标物体内的一个小块或者说种子区域，再在种子区域的基础上不断将其周围的像素点以一定的规则加入其中，达到最终将代表该物体的所有像素点结合成一个区域的目的。该方法的关键是要选择合适的生长或相似准则。生长准则一般可分为三种：基于区域灰度差准则、基于区域内灰度分布统计性质准则和基于区域形状准则。分裂合并法是先将图像分割成很多的一致性较强的小区域，再按一定的规则将小区域融合成大区域，达到分割图像的目的。区域提取法的缺点是往往会造成过度分割，即将图像分割成过多的区域，因此近年来针对这种方法的研究较少。

1.2.3 阈值分割法

对灰度图像的取阈值分割就是先确定一个处于图像灰度取值范围之中的灰度阈值，然后将图像中各个像素的灰度值都与这个阈值相比较，并根据比较结果将对应的像素分为两类。这两类像素一般分属图像的两类区域，从而达到分割的目的。阈值分割算法主要有两个步骤：

（1）确定需要的阈值；

（2）将分割阈值与像素值比较以划分像素。

可以看出，确定一个最优阈值是分割的关键。现有的大部分算法都是集中在阈值确定的研究上。阈值分割方法根据图像本身的特点，可分为单阈值分割方法和多阈值分割方法：也可分为基于像素值的阈值分割方法、基于区域性质的阈值分割方法和基于坐标位置的阈值分割方法．若考虑分割算法所用的特征或准则的特点，还可以分为直方图与直方图变换法、最大类空间方差法、最小误差法与均匀化误差法、共生矩阵法、最大熵法、简单统计法与局部特性法、概率松弛法、模糊集法等。

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

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

       由于图像的直方图可以看作是一维信号，而直方图上的突变点（波峰点和波谷点），往往可以代表图像灰度变化的特征。因此Jean-Christophe Olivo提出了用小波变换对直方图进行处理的方法实现自动阈值提取。Olivo通过检测直方图小波变换的奇异点和区域极值点给出直方图峰值点的特性。而小波变换的波峰和波谷点可以代表图像中灰度代表值和阈值点。利用小波变换多尺度特性实现对图像的阈值分割。又由于小波变换具有多分辨率的特性，因此可以通过对医学图像直方图的小波变换，实现由粗到细的多层次结构的阈值分割。首先在最低分辨率一层进行，然后逐渐向高层推进。小波变换的零交叉点表示了在分辨率2j时低通信号的局部跳变点。当尺度2j减小时，信号的局部微小细节逐渐增多，因此，能够检测出各微小细节的灰度突变点；当尺度2j增大时，信号的局部细节逐渐消失，而结构较大的轮廓却能清晰地反映出来，因而能检测出该结构较大的灰度突变点。因此，可以选择小波为光滑函数的二阶导数，对图像的一维直方图信号进行小波变换，检测出直方图信号的突变点，由此搜索出两峰之问的谷点作为分割阈值点。这就是小波变换用于图像分割的基本原理。对图像的直方图来说，它的各层的小波分解系数表示不同分辨率下的细节信号，它与小波近似信号联合构成直方图的多分辨率小波分解表示。给定直方图，考虑其多分辨率小波分解表示的零交叉点和极值点来确定直方图的峰值点和谷点。

多分辨率阈值选取

      基于直方图和小波变换的图像分割技术由以下几个步骤组成：首先由粗分辨率下的图像直方图细节信息确定分割区域类数；其次，在相邻峰之间自动确定最优阈值；最后用求出的最优阈值分割原图像。

      由于图像的原始直方图一般不够平滑或含有一定的噪声，因此，有必要对原始直方图进行平滑处理，以利于分割目标。其方法为：在空间域中采用保护边缘平滑方法平滑直方图，它既能保留原直方图基本变化特性，又能消除小峰的跳动。或者选取大尺度下的小波变换系数对直方图进行处理，也可以减小噪声的影响。

       在分辨率为2j时，由小波分解后的直方图近似信号的极大值确定初始区域类数，即确定峰的数目。对于灰度级数不多的原始影  像，一个区域类通常对应直方图中的一个峰，然而，对于一幅复杂图像，经小波分解后平滑直方图中的每个峰则不一定都对应一个区域类，它也可能从属于邻近的一个峰，因而有必要通过检查认定哪些峰对应于分割区域类。

峰的独立性判断是为了消除不能成为一类的那些峰，独立峰应满足三个条件：

(1)应具有一定的灰度范围；

(2)应具有一定峰下面积；

(3)应具有一定的峰谷差。

峰与峰之间的谷点选取首先在最低分辨率层进行，然后逐层推进，直到信号的最高分辨率层，并在每一层的阈值选取中，采用前一层的选取结构为引导。这种由粗到精的控制策略，能有效地选取最优分割阈值，同时较好地克服了噪声干扰和搜索空间大的问题。

        设最优阈值分别为为T1，T2，T3，…，Tk(k为正整数)，即可用这些最优阈值分割原始图像f(x，y)，得到分割结果g(x，y)，公式如下：

其中Ck(k=O，1，．．．，K)表示分割后的类别代码。

（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;

        else I(i,j)=0;

        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;

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

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;

        break;

    end

end

for i=L:-1:1

    if count(i)~=0

        nd=i-1;

        break;

    end

end

f=count(st+1:nd+1);  %f是每个灰度出现的概率

p=st;   q=nd-st;

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';

[x,map]=imread(imgname); %x是像素色值，map是色谱

figure(1);

subplot(2,2,1);imshow(x);

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;

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;

snr=10*log10(P/mseValue);%计算信噪比SNR

psnr=10*log10(255^2/mseValue);%计算峰值信噪比PSNR