圖像分割是一種重要的圖像分析技術。對圖像分割的研究一直是圖像技術研究中的熱點和焦點。醫學圖像分割是圖像分割的一個重要應用領域，也是一個經典難題，至今已有上千種分割方法，既有經典的方法也有結合新興理論的方法。
本論文首先介紹了雙峰法以及最大類方差自動門檻值法，然後重點介紹一種基于小波變換的圖像分割方法,該方法先對圖像的灰階直方圖進行小波多尺度變換,然後從較大的尺度系數到較小的尺度系數逐漸定位出灰階門檻值。最後，對這幾種算法的分割效果進行了比較。實驗結果表明, 本設計能夠實時穩定的對目标分割提取，分割效果良好。
醫學圖像分割是醫學圖像進行中的一個經典難題。圖像分割能夠自動或半自動描繪出醫學圖像中的解剖結構和其它感興趣的區域，進而有助于醫學診斷。
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值越大，就代表失真越少。
| 客觀評價 分割方法 分割方法 | TIME | SNR | PSNR | Threshold |
| 雙峰法 | 0.859649s | 35.48 | 1.56 | 40 |
| 最大類方差法 | 0.242231s | 35.18 | 1.26 | 62 |
| 小波變換 | 0.841680s | 48.22 | 14.30 | 253 |
通常一幅圖像分割結果的好與壞，以人的主觀判斷作為評價标準[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