針對現有醫學圖像中存在噪聲干擾與邊緣信號弱等現象,本文通過對二維小波變換進行研究,同時結合圖像的邊緣的方向性與小波系數的相關性,提出一種基于小波特性與邊緣模糊檢測的醫學圖像處理算法。該算法通過改進小波變換與傳統邊緣模糊檢測算法,來提高算法的降噪能力與邊緣優化效果。結果表明,其實驗結果與預測目標基本相符,該算法能夠有效的降低醫學圖像中的噪聲信號同時有效的保留圖像的邊緣信號,具有清晰度高、降噪能力強等優點。
引用本文： 朱柏輝, 萬智萍. 基于小波特性與邊緣模糊檢測的醫學圖像處理. 生物醫學工程學雜志, 2014, 31(3): 493-498. doi: 10.7507/1001-5515.20140091 復制
引言
現代醫學中，圖像的處理被廣泛的運用在診斷、治療等方面，如X光透視、CT、超聲透視等等，為治療提供了寶貴的信息。但由于人體的諸多因素的影響，使得采集回來的圖像信號往往具有信號弱、噪聲強等缺點。因此，如何對醫學圖像降噪與強化邊緣信號成為了研究的關鍵問題。隨圖像處理技術的發展，諸多學者提出了許多有關醫學圖像處理的算法[1-3]，其中常見的算法有插值算法[4]、邊緣模糊檢測[5]、Snake圖像分割算法[6]等，在一定程度上加強了醫學圖像的處理能力。小波變換憑借其變聚焦特性被譽為“數學顯微鏡”，且在圖像的邊緣檢測與消除噪聲方面具有其他算法所無法比擬的優越性。通過近幾年的研究發現，小波變換與邊緣模糊檢測算法仍有需要改進的地方，如：小波變換“過扼殺”小波系數的傾向與不能較好的表示圖像中線與面的奇異性,使得小波變換在圖像降噪與邊緣檢測具有一定的局限性現象。而邊緣模糊算法在邊緣檢測的過程中，將小波系數較小的信號定義為噪聲信號，因而忽略了其中所含有的一些有效信號。
為了解決上述問題，本研究通過引入方向小波變換來提高算法對小波方向特性的利用；在降噪的同時保護圖像中系數較小的有效信號，本文提出小波變換的相關性函數，對邊緣模糊檢測算法進行改進，從而提高算法的檢測效果與降噪能力，實現噪聲信號與邊緣信號的精確分類以及邊緣的快速檢測。
1 基于小波變換的圖像邊緣檢測算法
可以得到f(t)的小波變換[7]為
| ${{W}_{f}}\left( a,b \right)=\frac{1}{\sqrt{{{C}_{\psi }}}}~\frac{1}{\sqrt{\left| a \right|}}\int\limits_{-\infty }^{+\infty }{f\left( t \right)\bar{\Psi }}\left( \frac{t-b}{a} \right)dt,$ |
其中,ψ(ω)為ψ(t)的傅里葉變換，a為伸縮因子，b為平移因子。
設θ(x,y)滿足下列條件
| $\left\{ \begin{align} & \int {{\int }_{{{R}^{2}}}}\theta \left( x,y \right)dxdy=1 \\ & \underset{x,y\to +\infty }{\mathop{lim}}\,\theta \left( x,y \right)=0~ \\ \end{align} \right.,$ |
則稱θ(x,y)為二元平滑函數。且θ(x,y)滿足小波函數為
| $\left\{ \begin{align} & {{\psi }^{1}}\left( x,y \right)=\frac{\partial \theta \left( x,y \right)}{\partial x} \\ & {{\psi }^{2}}\left( x,y \right)=\frac{\partial \theta \left( x,y \right)}{\partial y} \\ \end{align} \right.$ |
將ψ1(x,y)、ψ2(x,y)作為二維小波變換的基函數，即
| $\left\{ \begin{align} & {{\psi }^{1}}\left( x,y \right)=\frac{1}{{{2}^{2}}}{{\psi }^{1}}\left( \frac{x}{{{2}^{j}}},\frac{y}{{{2}^{j}}} \right) \\ & {{\psi }^{2}}\left( x,y \right)=\frac{1}{{{2}^{2}}}{{\psi }^{2}}(\frac{x}{{{2}^{j}}}),\frac{y}{{{2}^{j}}}) \\ \end{align} \right.,$ |
則在尺度為2j時,函數f(x,y)沿水平方向和垂直方向的二進小波變換為
| $\left\{ \begin{align} & {{W}_{{{2}^{j}}}}^{1}f\left( x,y \right)=f*{{\psi }_{{{2}^{j}}}}^{1}f\left( x,y \right) \\ & {{W}_{{{2}^{j}}}}^{2}f\left( x,y \right)=f*{{\psi }_{{{2}^{j}}}}^{2}f\left( x,y \right)\text{ } \\ \end{align} \right.,$ |
記為W2jf(x,y)=W2j1f(x,y)
| $\eqalign{ & {W_{{2^j}}}f\left( {x,y} \right) = \left[ \matrix{ {W_{{2^j}}}^1f\left( {x,y} \right) \hfill \cr {W_{{2^j}}}^2f\left( {x,y} \right) \hfill \cr} \right] = \cr & {2^j}\left[ \matrix{ {\partial \over {\partial x}}f{\theta _{{2^j}}}\left( {x,y} \right) \hfill \cr {\partial \over {\partial y}}f{\theta _{{2^j}}}\left( {x,y} \right) \hfill \cr} \right] = {2^j}\nabla (f{\theta _{{2^j}}})\left( {x,y} \right), \cr} $ |
定義尺度2j時，函數f(x,y)小波變換的小波系數模與幅角表達函數分別為
| $模{{M}_{{{2}^{j}}}}f\left( x,y \right)=\sqrt{|{{W}_{{{2}^{j}}}}^{1}f\left( x,y \right){{|}^{2}}+|{{W}_{{{2}^{j}}}}^{2}f\left( x,y \right){{|}^{2}}}$ |
| $幅角{{A}_{{{2}^{j}}}}f\left( x,y \right)=arctan(\frac{{{W}_{{{2}^{j}}}}^{2}f\left( x,y \right)}{{{W}_{{{2}^{j}}}}^{1}f\left( x,y \right)})$ |
2 傳統的邊緣模糊檢測算法
2.1 傳統算法
通過小波分解，將圖像的多尺度信號描述為模糊矩陣，然后利用模糊集的并-交運算進行信息的合成，最后通過λ截矩陣完成邊緣圖像的獲取。其具體算法為：
(1) 圖像經過多尺度小波分解。
(2) 獲得各尺度下的模值矩陣M2jf(x,y)，濾除兩個尺度間的模值方向相反的噪聲點，得到候選邊緣信息M2jf(x,y)，該算法簡單的將較大閾值定義為邊緣，將較小閾值定義為非邊緣，而介于這兩個閾值之間則定義為其存在一定的模糊性，既可能是邊緣也可能不是。
(3) 建立M2jf(x,y)的隸屬函數[8-9]獲取圖像的有效邊緣信號。
M2jf(x,y)的隸屬函數為
| $U\left[ {{M_{{2^j}}}f\left( {x,y} \right)} \right] = \left\{ \matrix{ 1,{M_{{2^j}}}f\left( {x,y} \right) \ge \beta \hfill \cr {(1 + {(Max - {M_{{2^j}}}f\left( {x,y} \right))^2})^{ - 1}},\alpha {M_{{2^j}}}f\left( {x,y} \right)\beta \hfill \cr 0{M_{{2^j}}}f\left( {x,y} \right) \le \alpha , \hfill \cr} \right.$ |
其中Max為梯度方向上小領域中的最大值，α與β根據圖像的不同而進行改變。
根據M2jf(x,y)的隸屬函數完成候選邊緣信息矩陣M2jf(x,y)的模糊化，并得到相應尺度的模糊矩陣U2jf(x,y)。最后，利用模糊子集的并-交運算合成各尺度信息。其運算函數為：
交運算：
| $Uf\left( x,y \right)=U[({{U}_{{{2}^{i}}}}f\left( x,y \right))\cap ({{U}_{{{2}^{j}}}}f\left( x,y \right))]|i\ne j,i、j\in s$ |
并運算：
| $Uf\left( x,y \right)=U[({{U}_{{{2}^{i}}}}f\left( x,y \right))\cap ({{U}_{{{2}^{j}}}}f\left( x,y \right))]|i\ne j,i、j\in s$ |
最終得到Uf(x,y)為合成后的模糊矩陣，再利用λ截矩陣完成邊緣圖像的獲取，設Uf(x,y)=(Ui,j)m×n，?λ∈[0, 1],記Uλf(x,y)=(Uij(λ))m×n，其中
| ${{U}_{ij}}(\lambda )=\left\{ \begin{align} & 1,{{U}_{ij}}>\lambda ~ \\ & 0,{{U}_{ij}}＜\lambda \\ \end{align} \right.$ |
則Uλf(x,y)為Uf(x,y)的λ截矩陣。得到圖像最終的邊緣信號。
2.2 存在的問題
算法在對圖像的邊緣信號進行邊緣檢測與圖像的降噪處理時，將小波模值較小的系數簡單的定義為圖像噪聲信號，忽略了小波模值較小的系數中含有的有效信息，使一些本來有用的突變點被當成噪聲平滑掉了。而在對邊緣進行模糊的過程中，圖像中的邊緣方向特性在圖像的局部變化與沿著某一方向上的全局中體現，該算法未能利用小波變換的各子帶間的方向特性，如果能夠很好的利用這一特性，將大大提高算法的邊緣檢測效果。
3 本文算法
3.1 方向小波變換
傳統的小波變換能夠將圖像的高頻子帶分解為水平、垂直以及對角線三個方向。傳統的二維小波變換只能描繪出圖像在這三個方向上的屬性而無法反映其方向上的特征。如果能夠很好的利用小波變換的特性，可為圖像的分析與處理提供很大的幫助。 由小波變換定理可知，對于R2中的任一條直線來說，均可以選取適當的(x,y,θ)，記為
| $xcos\theta +ysin\theta -\gamma =0$ |
因此，設f(x,y)是N中的任一函數，引用一個方向角度θ變量反映圖像的方向性質，用于描述圖像在不同方向上的變化。即ψ(x)∈L′(R)且 A G (0)=0 ,令ψs(x,y,θ)=ψs(xcosθ+ysinθ-γ)。由于函數f(x,y)∈L2(R2)，因此，可令其方向小波函數為
| $\begin{align} & W{{\psi }_{s}}f(x,y,\theta )=(f*{{\psi }_{s}})(x,y,\theta )= \\ & {{\int }_{R}}{{\int }_{R}}f\left( x,y \right){{\psi }_{s}}(xcos\theta ,ysin\theta -\gamma )dxdy~ \\ \end{align}$ |
此函數可用于描述小波變換的方向，s是小波變換方向的尺度，此時用θ來描述該函數的變化情況。
設θ(x,y,θ)為平滑函數，則
| $\eqalign{ & \psi (x,y,\theta ) = \left\{ \matrix{ {\partial \over {\partial x}}\Phi (xcos\theta + ysin\theta - \gamma ) = \hfill \cr \Phi \prime (xcos\theta + ysin\theta - \gamma )cos\theta \hfill \cr {\partial \over {\partial y}}\Phi (xcos\theta + ysin\theta - \gamma ) = \hfill \cr \Phi \prime (xcos\theta + ysin\theta - \gamma )sin\theta \hfill \cr} \right.; \cr & \left\{ \matrix{ {W_{{2^j}}}^1f(x,y,\theta ) = f*\left( {{\partial \over {\partial x}}{\Phi _{{2^j}}}} \right)(x,y,\theta ) \hfill \cr = f*\Phi \prime {2^j}(x,y,\theta )cos\theta \hfill \cr {W_{{2^j}}}^2f(x,y,\theta ) = f*\left( {{\partial \over {\partial y}}{\Phi _{{2^j}}}} \right)(x,y,\theta ) \hfill \cr = f*\Phi \prime {2^j}(x,y,\theta )\sin \theta \hfill \cr} \right., \cr} $ |
因此，可以得到小波函數的模函數為：
| ${{M}_{{{2}^{j}}}}f(x,y,\theta )=\sqrt{|{{W}_{{{2}^{j}}}}^{1}f(x,y,\theta ){{|}^{2}}+|{{W}_{{{2}^{j}}}}^{2}f(x,y,\theta ){{|}^{2}}}$ |
圖像中的邊緣方向特性在圖像的局部變化與沿著某一方向上的全局中體現，傳統的小波變換無法體現這一特性，而方向小波變換卻能反映圖像邊緣信號的這一特性，因此，方向小波變換具有傳統小波變換良好的時頻局部化分析能力和方向分析特性。
3.2 改進邊緣模糊檢測
通過小波分解后，圖像的相關性[10]：Corr=Wn*Wn+1，計算可得第n層相關性的總功率(PC=ΣCorr*Corr)和小波系數的總功率(PW=ΣWn*Wn)； 進而得到相關值： Corrn=Corr*PW/PC ，即可得到此相關值的最小絕對值：Corrn min=min{|Corrn min|}和最大絕對值：Corrn max=max{|Corrn max|}。
由于圖像中的噪聲信號是隨機分散的，并不具備圖像的相關性，所以本文通過判斷其相關性的大小來對圖像進行降噪處理，并且將介于最大相關值與最小相關值之間信號定為模糊邊緣信號，得到公式為
| $\eqalign{ & U[{M_{{2^j}}}f(x,y,\theta )] = \cr & \left\{ \matrix{ 1,\left| {Cor{r_n}} \right| > 0.45Cor{r_{nmax}} \hfill \cr {(1 + {(Max - {M_{{2^j}}}f(x,y,\theta ))^2})^{ - 1}},Cor{r_{nmin}}\left| {Cor{r_n}} \right|0.45Cor{r_{nmax}} \hfill \cr 0,\left| {Cor{r_n}} \right|Cor{r_{nmin}} \hfill \cr} \right. \cr} $ |
從而得到候選的邊緣信號M2jf(x,y,θ)的模糊,并得到相應模糊矩陣U2jf(x,y,θ);通過各模糊子集的并-交運算合成最后各尺度信息。本文只運用于不同等級的子帶模糊矩陣間的交運算，根據各子帶小波分解后各子帶間的方向性，得到其運算函數為并運算：
| $Uf(x,y,\theta )=({{U}_{{{2}^{n}}}}f(x,y,\theta ))\cap ({{U}_{{{2}^{n+1}}}}f(x,y,\theta ))~,$ |
交運算：
| $Uf(x,y,\theta )=U({{U}_{{{2}^{n}}}}f(x,y,\theta ))\cap ({{U}_{{{2}^{n+1}}}}f(x,y,\theta ))~,$ |
即Uf(x,y,θ)為新合成的模糊矩陣；最后通過引入閾值來獲取其邊緣圖像，設Uf(x,y,θ)=(Uij)N2,?λ∈[0, 1],記Uλf(x,y,θ)=(Uij(λ))N2 ，
其中
| ${{U}_{ij}}(\lambda )=\left\{ \begin{align} & 1,{{U}_{ij}}>\lambda ~ \\ & 0,{{U}_{ij}}＜\lambda \\ \end{align} \right.$ |
將Uλf(x,y,θ)映射到灰度空間即得到相應的邊緣圖像。
3.3 算法設計
步驟一： 開始。
步驟二： 對大小為N×N的醫學圖像進行3層小波分解，得到圖像在各尺度下的LH、HL、HH高頻子帶和LL低頻子帶。
步驟三： 計算分解后高頻子帶的小波系數的相關性，從而獲取圖像各子帶間的相關性，并通過以上相關值公式計算出各子帶小波系數的相關值。
步驟四： 通過采集獲得的各子帶間的小波相關值，然后對圖像進行方向小波求解，采用公式(16)得到小波系數的模值，最后采用公式(17)對其進行小波邊緣檢測，提取有效的邊緣信號。
步驟五： 采用公式(20)獲取新的邊緣圖像。
步驟六： 圖像重構。
4 仿真實驗
為了驗證算法在醫學圖像處理上的優越性，在MATLAB 7.10.0 的環境下對算法進行仿真，采用傳統的邊緣模糊檢測算法、傳統的小波邊緣算法、Sobel算法、文獻算法[5]以及本文算法分別對同一幅含有噪聲的圖像進行圖像處理，通過對比峰值信噪比(PSNR)與圖像的均方誤差(MSE)作為評價的標準，即。
峰值信噪比函數式
| $PNSR=10ln\frac{\sum\limits_{x=1,y=1}^{{{N}^{2}}}{{{255}^{2}}}}{\sum\limits_{x=1,y=1}^{{{N}^{2}}}{{{\left( f\left( x,y \right)-\text{ }Af-U\text{ }\left( x,y \right) \right)}^{2}}}~},$ |
均方誤差函數式
| $MSE={{\left[ \frac{1}{{{N}^{2}}}\sum\limits_{x=1,y=1}^{{{N}^{2}}}{{{\left( f\left[ x,y \right]-\bar{f}\left( x,y \right) \right)}^{2}}} \right]}^{1/2}}$ |
其中PSNR 是最普遍、最廣泛使用的評鑒畫質的客觀量測法，PSNR的值越大，則說明其邊緣提前的效果越好。而MSE能夠評價數據的變化程度，其中MSE的值越小，則說明預測模型描述實驗數據具有更好的精確度。
本文選取一幅腦部CT圖片作為仿真對象，為了更好的體現算法在降噪方面的優勢，通過改變其圖像的噪聲量的大小來對本文算法進行驗證；經過仿真，得到圖 1為各算法在不同噪聲量下的PNSR曲線圖，圖 2為各算法在不同噪聲量下的MSE曲線圖。
圖1
各算法在不同噪聲量下的PNSR曲線圖
Figure1.
PNSR graph of each algorithm for type of noise amount
圖2
各算法在不同噪聲量下的MSE曲線圖
Figure2.
MSE curves of each algorithm for different amount of noise
由圖 1可以清晰地看到,在噪聲量不斷增加的情況下,傳統的小波邊緣算法與邊緣模糊檢測算法所得到的PNSR曲線很接近。Sobel算法的PSNR曲線始終保持在本文算法曲線的下方，并隨著噪聲量的增大，其差距越來越大。采用文獻算法所得到的PSNR曲線，雖然一開始能夠得到比本文算法更高的PSNR值，但當噪聲量超過1.5時，本文算法的PSNR曲線開始高于文獻算法的PSNR曲線。由圖 2可以看到，本文算法優于小波算法、邊緣模糊檢測算法以及Sobel算法。相比文獻算法，盡管本文得到的最初的MSE值不如文獻算法，但其良好的穩定性使得本文算法所得的MSE曲線的增長幅度較小，并且隨著噪聲量的增大，其文獻算法的MSE曲線的增長幅度大于本文算法。
為了方便分析與對比數據變化的大小，表 1、2列出了各算法在不同噪聲量的PNSR值MSE值。
表1
各算法在不同噪聲量的PNSR值
Table1.
PNSR values of each algorithm in the type of noise amount
表2
各算法在不同噪聲量的MSE值
Table2.
MSE values of each algorithm in the type of noise amount
由表 1可以看出，傳統的邊緣模糊檢測算法比傳統的小波邊緣算法的PNSR值平均高1.5 dB，而本文算法與傳統的邊緣模糊檢測算法從開始高出2.89 dB到最后變為高出5.53 dB。而從表 2可以看出，本文算法所測得的MSE值始終低于傳統的小波邊緣算法、邊緣模糊檢測算法以及Sobel算法。相對于文獻算法，則是先高后低，其中最低低至1.14。為了進一步驗證算法在去噪與邊緣搜索的效果，通過加入噪聲量δ=3.0的腦部CT圖片進行算法仿真，得各算法仿真效果圖，如圖 3所示。
圖3
各算法的仿真效果圖
(a)原始圖像；(b)加噪圖像；(c)小波算法；(d)邊緣模糊檢測算法；(e)Sobel算法；(f)文獻算法；(g)本文算法
Figure3. The simulation effect image of each algorithm(a) original image; (b) noisy images; (c) wavelet algorithm; (d) fuzzy edge detection algorithm; (e) Sobel algorithm; (f) literature algorithm; (g) our algorithm
由圖 3可以看到，小波算法具有很好的去噪能力，但是在去噪的同時，容易將圖像的有效信號一起除去；邊緣模糊檢測算法雖然能較好地提取邊緣同時除去噪聲，但對圖像的弱邊緣檢測不全；Sobel算法雖然能夠得到清晰的圖像，但也存在著邊緣檢測能力弱等問題；文獻算法能夠得到圖像的總體輪廓，但其邊緣清晰度不高；而本文算法所到的仿真圖像優于其他算法。
通過對各算法的數據與仿真圖的分析，可以發現邊緣模糊檢測算法與傳統的小波邊緣算法采用的去噪方法原理相同,都是通過將較小的小波系數定義為噪聲量,而忽略了其中所含的有用的突變點。Sobel算法雖然能夠得到清晰的圖像，但其檢測效果最差,容易出現邊緣的損失。文獻算法雖然能夠得到良好的邊緣圖像，但其對噪聲的處理相對不佳，在不同的噪聲量下無法保持穩定性。本文算法通過提取圖像中各子帶間的相關值,除去相關值較小的小波系數，使得算法能夠在對圖像進行去噪的同時保持良好的邊緣檢測效果，進而有效的證明了算法的穩定性與可行性。
5 結束語
本文在小波變換的基礎上，通過換算提出了方向小波變換以提高算法對邊緣信號的方向性特征的提取。與此同時，在算法中還結合小波分解后各子帶間小波系數的相關性，對邊緣模糊檢測算法進行改進，以提高算法的邊緣提取能力與去噪效果。經實驗，該算法能夠在去噪的同時有效的保護圖像中的邊緣信號。在醫學圖像處理中，具有去噪能力強，邊緣檢測度高的優點。
引言
現代醫學中，圖像的處理被廣泛的運用在診斷、治療等方面，如X光透視、CT、超聲透視等等，為治療提供了寶貴的信息。但由于人體的諸多因素的影響，使得采集回來的圖像信號往往具有信號弱、噪聲強等缺點。因此，如何對醫學圖像降噪與強化邊緣信號成為了研究的關鍵問題。隨圖像處理技術的發展，諸多學者提出了許多有關醫學圖像處理的算法[1-3]，其中常見的算法有插值算法[4]、邊緣模糊檢測[5]、Snake圖像分割算法[6]等，在一定程度上加強了醫學圖像的處理能力。小波變換憑借其變聚焦特性被譽為“數學顯微鏡”，且在圖像的邊緣檢測與消除噪聲方面具有其他算法所無法比擬的優越性。通過近幾年的研究發現，小波變換與邊緣模糊檢測算法仍有需要改進的地方，如：小波變換“過扼殺”小波系數的傾向與不能較好的表示圖像中線與面的奇異性,使得小波變換在圖像降噪與邊緣檢測具有一定的局限性現象。而邊緣模糊算法在邊緣檢測的過程中，將小波系數較小的信號定義為噪聲信號，因而忽略了其中所含有的一些有效信號。
為了解決上述問題，本研究通過引入方向小波變換來提高算法對小波方向特性的利用；在降噪的同時保護圖像中系數較小的有效信號，本文提出小波變換的相關性函數，對邊緣模糊檢測算法進行改進，從而提高算法的檢測效果與降噪能力，實現噪聲信號與邊緣信號的精確分類以及邊緣的快速檢測。
1 基于小波變換的圖像邊緣檢測算法
可以得到f(t)的小波變換[7]為
| ${{W}_{f}}\left( a,b \right)=\frac{1}{\sqrt{{{C}_{\psi }}}}~\frac{1}{\sqrt{\left| a \right|}}\int\limits_{-\infty }^{+\infty }{f\left( t \right)\bar{\Psi }}\left( \frac{t-b}{a} \right)dt,$ |
其中,ψ(ω)為ψ(t)的傅里葉變換，a為伸縮因子，b為平移因子。
設θ(x,y)滿足下列條件
| $\left\{ \begin{align} & \int {{\int }_{{{R}^{2}}}}\theta \left( x,y \right)dxdy=1 \\ & \underset{x,y\to +\infty }{\mathop{lim}}\,\theta \left( x,y \right)=0~ \\ \end{align} \right.,$ |
則稱θ(x,y)為二元平滑函數。且θ(x,y)滿足小波函數為
| $\left\{ \begin{align} & {{\psi }^{1}}\left( x,y \right)=\frac{\partial \theta \left( x,y \right)}{\partial x} \\ & {{\psi }^{2}}\left( x,y \right)=\frac{\partial \theta \left( x,y \right)}{\partial y} \\ \end{align} \right.$ |
將ψ1(x,y)、ψ2(x,y)作為二維小波變換的基函數，即
| $\left\{ \begin{align} & {{\psi }^{1}}\left( x,y \right)=\frac{1}{{{2}^{2}}}{{\psi }^{1}}\left( \frac{x}{{{2}^{j}}},\frac{y}{{{2}^{j}}} \right) \\ & {{\psi }^{2}}\left( x,y \right)=\frac{1}{{{2}^{2}}}{{\psi }^{2}}(\frac{x}{{{2}^{j}}}),\frac{y}{{{2}^{j}}}) \\ \end{align} \right.,$ |
則在尺度為2j時,函數f(x,y)沿水平方向和垂直方向的二進小波變換為
| $\left\{ \begin{align} & {{W}_{{{2}^{j}}}}^{1}f\left( x,y \right)=f*{{\psi }_{{{2}^{j}}}}^{1}f\left( x,y \right) \\ & {{W}_{{{2}^{j}}}}^{2}f\left( x,y \right)=f*{{\psi }_{{{2}^{j}}}}^{2}f\left( x,y \right)\text{ } \\ \end{align} \right.,$ |
記為W2jf(x,y)=W2j1f(x,y)
| $\eqalign{ & {W_{{2^j}}}f\left( {x,y} \right) = \left[ \matrix{ {W_{{2^j}}}^1f\left( {x,y} \right) \hfill \cr {W_{{2^j}}}^2f\left( {x,y} \right) \hfill \cr} \right] = \cr & {2^j}\left[ \matrix{ {\partial \over {\partial x}}f{\theta _{{2^j}}}\left( {x,y} \right) \hfill \cr {\partial \over {\partial y}}f{\theta _{{2^j}}}\left( {x,y} \right) \hfill \cr} \right] = {2^j}\nabla (f{\theta _{{2^j}}})\left( {x,y} \right), \cr} $ |
定義尺度2j時，函數f(x,y)小波變換的小波系數模與幅角表達函數分別為
| $模{{M}_{{{2}^{j}}}}f\left( x,y \right)=\sqrt{|{{W}_{{{2}^{j}}}}^{1}f\left( x,y \right){{|}^{2}}+|{{W}_{{{2}^{j}}}}^{2}f\left( x,y \right){{|}^{2}}}$ |
| $幅角{{A}_{{{2}^{j}}}}f\left( x,y \right)=arctan(\frac{{{W}_{{{2}^{j}}}}^{2}f\left( x,y \right)}{{{W}_{{{2}^{j}}}}^{1}f\left( x,y \right)})$ |
2 傳統的邊緣模糊檢測算法
2.1 傳統算法
通過小波分解，將圖像的多尺度信號描述為模糊矩陣，然后利用模糊集的并-交運算進行信息的合成，最后通過λ截矩陣完成邊緣圖像的獲取。其具體算法為：
(1) 圖像經過多尺度小波分解。
(2) 獲得各尺度下的模值矩陣M2jf(x,y)，濾除兩個尺度間的模值方向相反的噪聲點，得到候選邊緣信息M2jf(x,y)，該算法簡單的將較大閾值定義為邊緣，將較小閾值定義為非邊緣，而介于這兩個閾值之間則定義為其存在一定的模糊性，既可能是邊緣也可能不是。
(3) 建立M2jf(x,y)的隸屬函數[8-9]獲取圖像的有效邊緣信號。
M2jf(x,y)的隸屬函數為
| $U\left[ {{M_{{2^j}}}f\left( {x,y} \right)} \right] = \left\{ \matrix{ 1,{M_{{2^j}}}f\left( {x,y} \right) \ge \beta \hfill \cr {(1 + {(Max - {M_{{2^j}}}f\left( {x,y} \right))^2})^{ - 1}},\alpha {M_{{2^j}}}f\left( {x,y} \right)\beta \hfill \cr 0{M_{{2^j}}}f\left( {x,y} \right) \le \alpha , \hfill \cr} \right.$ |
其中Max為梯度方向上小領域中的最大值，α與β根據圖像的不同而進行改變。
根據M2jf(x,y)的隸屬函數完成候選邊緣信息矩陣M2jf(x,y)的模糊化，并得到相應尺度的模糊矩陣U2jf(x,y)。最后，利用模糊子集的并-交運算合成各尺度信息。其運算函數為：
交運算：
| $Uf\left( x,y \right)=U[({{U}_{{{2}^{i}}}}f\left( x,y \right))\cap ({{U}_{{{2}^{j}}}}f\left( x,y \right))]|i\ne j,i、j\in s$ |
并運算：
| $Uf\left( x,y \right)=U[({{U}_{{{2}^{i}}}}f\left( x,y \right))\cap ({{U}_{{{2}^{j}}}}f\left( x,y \right))]|i\ne j,i、j\in s$ |
最終得到Uf(x,y)為合成后的模糊矩陣，再利用λ截矩陣完成邊緣圖像的獲取，設Uf(x,y)=(Ui,j)m×n，?λ∈[0, 1],記Uλf(x,y)=(Uij(λ))m×n，其中
| ${{U}_{ij}}(\lambda )=\left\{ \begin{align} & 1,{{U}_{ij}}>\lambda ~ \\ & 0,{{U}_{ij}}＜\lambda \\ \end{align} \right.$ |
則Uλf(x,y)為Uf(x,y)的λ截矩陣。得到圖像最終的邊緣信號。
2.2 存在的問題
算法在對圖像的邊緣信號進行邊緣檢測與圖像的降噪處理時，將小波模值較小的系數簡單的定義為圖像噪聲信號，忽略了小波模值較小的系數中含有的有效信息，使一些本來有用的突變點被當成噪聲平滑掉了。而在對邊緣進行模糊的過程中，圖像中的邊緣方向特性在圖像的局部變化與沿著某一方向上的全局中體現，該算法未能利用小波變換的各子帶間的方向特性，如果能夠很好的利用這一特性，將大大提高算法的邊緣檢測效果。
3 本文算法
3.1 方向小波變換
傳統的小波變換能夠將圖像的高頻子帶分解為水平、垂直以及對角線三個方向。傳統的二維小波變換只能描繪出圖像在這三個方向上的屬性而無法反映其方向上的特征。如果能夠很好的利用小波變換的特性，可為圖像的分析與處理提供很大的幫助。 由小波變換定理可知，對于R2中的任一條直線來說，均可以選取適當的(x,y,θ)，記為
| $xcos\theta +ysin\theta -\gamma =0$ |
因此，設f(x,y)是N中的任一函數，引用一個方向角度θ變量反映圖像的方向性質，用于描述圖像在不同方向上的變化。即ψ(x)∈L′(R)且 A G (0)=0 ,令ψs(x,y,θ)=ψs(xcosθ+ysinθ-γ)。由于函數f(x,y)∈L2(R2)，因此，可令其方向小波函數為
| $\begin{align} & W{{\psi }_{s}}f(x,y,\theta )=(f*{{\psi }_{s}})(x,y,\theta )= \\ & {{\int }_{R}}{{\int }_{R}}f\left( x,y \right){{\psi }_{s}}(xcos\theta ,ysin\theta -\gamma )dxdy~ \\ \end{align}$ |
此函數可用于描述小波變換的方向，s是小波變換方向的尺度，此時用θ來描述該函數的變化情況。
設θ(x,y,θ)為平滑函數，則
| $\eqalign{ & \psi (x,y,\theta ) = \left\{ \matrix{ {\partial \over {\partial x}}\Phi (xcos\theta + ysin\theta - \gamma ) = \hfill \cr \Phi \prime (xcos\theta + ysin\theta - \gamma )cos\theta \hfill \cr {\partial \over {\partial y}}\Phi (xcos\theta + ysin\theta - \gamma ) = \hfill \cr \Phi \prime (xcos\theta + ysin\theta - \gamma )sin\theta \hfill \cr} \right.; \cr & \left\{ \matrix{ {W_{{2^j}}}^1f(x,y,\theta ) = f*\left( {{\partial \over {\partial x}}{\Phi _{{2^j}}}} \right)(x,y,\theta ) \hfill \cr = f*\Phi \prime {2^j}(x,y,\theta )cos\theta \hfill \cr {W_{{2^j}}}^2f(x,y,\theta ) = f*\left( {{\partial \over {\partial y}}{\Phi _{{2^j}}}} \right)(x,y,\theta ) \hfill \cr = f*\Phi \prime {2^j}(x,y,\theta )\sin \theta \hfill \cr} \right., \cr} $ |
因此，可以得到小波函數的模函數為：
| ${{M}_{{{2}^{j}}}}f(x,y,\theta )=\sqrt{|{{W}_{{{2}^{j}}}}^{1}f(x,y,\theta ){{|}^{2}}+|{{W}_{{{2}^{j}}}}^{2}f(x,y,\theta ){{|}^{2}}}$ |
圖像中的邊緣方向特性在圖像的局部變化與沿著某一方向上的全局中體現，傳統的小波變換無法體現這一特性，而方向小波變換卻能反映圖像邊緣信號的這一特性，因此，方向小波變換具有傳統小波變換良好的時頻局部化分析能力和方向分析特性。
3.2 改進邊緣模糊檢測
通過小波分解后，圖像的相關性[10]：Corr=Wn*Wn+1，計算可得第n層相關性的總功率(PC=ΣCorr*Corr)和小波系數的總功率(PW=ΣWn*Wn)； 進而得到相關值： Corrn=Corr*PW/PC ，即可得到此相關值的最小絕對值：Corrn min=min{|Corrn min|}和最大絕對值：Corrn max=max{|Corrn max|}。
由于圖像中的噪聲信號是隨機分散的，并不具備圖像的相關性，所以本文通過判斷其相關性的大小來對圖像進行降噪處理，并且將介于最大相關值與最小相關值之間信號定為模糊邊緣信號，得到公式為
| $\eqalign{ & U[{M_{{2^j}}}f(x,y,\theta )] = \cr & \left\{ \matrix{ 1,\left| {Cor{r_n}} \right| > 0.45Cor{r_{nmax}} \hfill \cr {(1 + {(Max - {M_{{2^j}}}f(x,y,\theta ))^2})^{ - 1}},Cor{r_{nmin}}\left| {Cor{r_n}} \right|0.45Cor{r_{nmax}} \hfill \cr 0,\left| {Cor{r_n}} \right|Cor{r_{nmin}} \hfill \cr} \right. \cr} $ |
從而得到候選的邊緣信號M2jf(x,y,θ)的模糊,并得到相應模糊矩陣U2jf(x,y,θ);通過各模糊子集的并-交運算合成最后各尺度信息。本文只運用于不同等級的子帶模糊矩陣間的交運算，根據各子帶小波分解后各子帶間的方向性，得到其運算函數為并運算：
| $Uf(x,y,\theta )=({{U}_{{{2}^{n}}}}f(x,y,\theta ))\cap ({{U}_{{{2}^{n+1}}}}f(x,y,\theta ))~,$ |
交運算：
| $Uf(x,y,\theta )=U({{U}_{{{2}^{n}}}}f(x,y,\theta ))\cap ({{U}_{{{2}^{n+1}}}}f(x,y,\theta ))~,$ |
即Uf(x,y,θ)為新合成的模糊矩陣；最后通過引入閾值來獲取其邊緣圖像，設Uf(x,y,θ)=(Uij)N2,?λ∈[0, 1],記Uλf(x,y,θ)=(Uij(λ))N2 ，
其中
| ${{U}_{ij}}(\lambda )=\left\{ \begin{align} & 1,{{U}_{ij}}>\lambda ~ \\ & 0,{{U}_{ij}}＜\lambda \\ \end{align} \right.$ |
將Uλf(x,y,θ)映射到灰度空間即得到相應的邊緣圖像。
3.3 算法設計
步驟一： 開始。
步驟二： 對大小為N×N的醫學圖像進行3層小波分解，得到圖像在各尺度下的LH、HL、HH高頻子帶和LL低頻子帶。
步驟三： 計算分解后高頻子帶的小波系數的相關性，從而獲取圖像各子帶間的相關性，并通過以上相關值公式計算出各子帶小波系數的相關值。
步驟四： 通過采集獲得的各子帶間的小波相關值，然后對圖像進行方向小波求解，采用公式(16)得到小波系數的模值，最后采用公式(17)對其進行小波邊緣檢測，提取有效的邊緣信號。
步驟五： 采用公式(20)獲取新的邊緣圖像。
步驟六： 圖像重構。
4 仿真實驗
為了驗證算法在醫學圖像處理上的優越性，在MATLAB 7.10.0 的環境下對算法進行仿真，采用傳統的邊緣模糊檢測算法、傳統的小波邊緣算法、Sobel算法、文獻算法[5]以及本文算法分別對同一幅含有噪聲的圖像進行圖像處理，通過對比峰值信噪比(PSNR)與圖像的均方誤差(MSE)作為評價的標準，即。
峰值信噪比函數式
| $PNSR=10ln\frac{\sum\limits_{x=1,y=1}^{{{N}^{2}}}{{{255}^{2}}}}{\sum\limits_{x=1,y=1}^{{{N}^{2}}}{{{\left( f\left( x,y \right)-\text{ }Af-U\text{ }\left( x,y \right) \right)}^{2}}}~},$ |
均方誤差函數式
| $MSE={{\left[ \frac{1}{{{N}^{2}}}\sum\limits_{x=1,y=1}^{{{N}^{2}}}{{{\left( f\left[ x,y \right]-\bar{f}\left( x,y \right) \right)}^{2}}} \right]}^{1/2}}$ |
其中PSNR 是最普遍、最廣泛使用的評鑒畫質的客觀量測法，PSNR的值越大，則說明其邊緣提前的效果越好。而MSE能夠評價數據的變化程度，其中MSE的值越小，則說明預測模型描述實驗數據具有更好的精確度。
本文選取一幅腦部CT圖片作為仿真對象，為了更好的體現算法在降噪方面的優勢，通過改變其圖像的噪聲量的大小來對本文算法進行驗證；經過仿真，得到圖 1為各算法在不同噪聲量下的PNSR曲線圖，圖 2為各算法在不同噪聲量下的MSE曲線圖。
圖1
各算法在不同噪聲量下的PNSR曲線圖
Figure1.
PNSR graph of each algorithm for type of noise amount
圖2
各算法在不同噪聲量下的MSE曲線圖
Figure2.
MSE curves of each algorithm for different amount of noise
由圖 1可以清晰地看到,在噪聲量不斷增加的情況下,傳統的小波邊緣算法與邊緣模糊檢測算法所得到的PNSR曲線很接近。Sobel算法的PSNR曲線始終保持在本文算法曲線的下方，并隨著噪聲量的增大，其差距越來越大。采用文獻算法所得到的PSNR曲線，雖然一開始能夠得到比本文算法更高的PSNR值，但當噪聲量超過1.5時，本文算法的PSNR曲線開始高于文獻算法的PSNR曲線。由圖 2可以看到，本文算法優于小波算法、邊緣模糊檢測算法以及Sobel算法。相比文獻算法，盡管本文得到的最初的MSE值不如文獻算法，但其良好的穩定性使得本文算法所得的MSE曲線的增長幅度較小，并且隨著噪聲量的增大，其文獻算法的MSE曲線的增長幅度大于本文算法。
為了方便分析與對比數據變化的大小，表 1、2列出了各算法在不同噪聲量的PNSR值MSE值。
表1
各算法在不同噪聲量的PNSR值
Table1.
PNSR values of each algorithm in the type of noise amount
表2
各算法在不同噪聲量的MSE值
Table2.
MSE values of each algorithm in the type of noise amount
由表 1可以看出，傳統的邊緣模糊檢測算法比傳統的小波邊緣算法的PNSR值平均高1.5 dB，而本文算法與傳統的邊緣模糊檢測算法從開始高出2.89 dB到最后變為高出5.53 dB。而從表 2可以看出，本文算法所測得的MSE值始終低于傳統的小波邊緣算法、邊緣模糊檢測算法以及Sobel算法。相對于文獻算法，則是先高后低，其中最低低至1.14。為了進一步驗證算法在去噪與邊緣搜索的效果，通過加入噪聲量δ=3.0的腦部CT圖片進行算法仿真，得各算法仿真效果圖，如圖 3所示。
圖3
各算法的仿真效果圖
(a)原始圖像；(b)加噪圖像；(c)小波算法；(d)邊緣模糊檢測算法；(e)Sobel算法；(f)文獻算法；(g)本文算法
Figure3. The simulation effect image of each algorithm(a) original image; (b) noisy images; (c) wavelet algorithm; (d) fuzzy edge detection algorithm; (e) Sobel algorithm; (f) literature algorithm; (g) our algorithm
由圖 3可以看到，小波算法具有很好的去噪能力，但是在去噪的同時，容易將圖像的有效信號一起除去；邊緣模糊檢測算法雖然能較好地提取邊緣同時除去噪聲，但對圖像的弱邊緣檢測不全；Sobel算法雖然能夠得到清晰的圖像，但也存在著邊緣檢測能力弱等問題；文獻算法能夠得到圖像的總體輪廓，但其邊緣清晰度不高；而本文算法所到的仿真圖像優于其他算法。
通過對各算法的數據與仿真圖的分析，可以發現邊緣模糊檢測算法與傳統的小波邊緣算法采用的去噪方法原理相同,都是通過將較小的小波系數定義為噪聲量,而忽略了其中所含的有用的突變點。Sobel算法雖然能夠得到清晰的圖像，但其檢測效果最差,容易出現邊緣的損失。文獻算法雖然能夠得到良好的邊緣圖像，但其對噪聲的處理相對不佳，在不同的噪聲量下無法保持穩定性。本文算法通過提取圖像中各子帶間的相關值,除去相關值較小的小波系數，使得算法能夠在對圖像進行去噪的同時保持良好的邊緣檢測效果，進而有效的證明了算法的穩定性與可行性。
5 結束語
本文在小波變換的基礎上，通過換算提出了方向小波變換以提高算法對邊緣信號的方向性特征的提取。與此同時，在算法中還結合小波分解后各子帶間小波系數的相關性，對邊緣模糊檢測算法進行改進，以提高算法的邊緣提取能力與去噪效果。經實驗，該算法能夠在去噪的同時有效的保護圖像中的邊緣信號。在醫學圖像處理中，具有去噪能力強，邊緣檢測度高的優點。
