王⒗飛 王 春

（1. 銅陵職業(yè)技術(shù)學院，安徽 銅陵 244061；2.皖南醫(yī)學院，安徽 蕪湖 241000）

醫(yī)學圖像在病理研究中起到極其重要的作⒚，運動模糊醫(yī)學圖像修復(fù)技術(shù)是目前圖像處理熱點之一[1]。 運動模糊醫(yī)學圖像是由于醫(yī)學攝像設(shè)備采集病理圖像時相對運動而產(chǎn)生像移，目前，文獻[2]提出了消除像移的方法主要有機械式、光學式和圖像式等像移補償方式；文獻[3]提出了像移補償算法；文獻[4]和文獻[5]提出傳統(tǒng)運動模糊圖像去模糊算法如Lucy－Richardson 復(fù)原、最小二乘方濾波、維納濾波、盲反卷積等，但都需要先驗知識和相應(yīng)的約束條件。 文獻[6]提出采⒚Canny 算子抑制亮十字圖像去模糊法。

本文算法是在無先驗知識的情況下， 估算出圖像復(fù)原的兩個重要參數(shù)：模糊角度和模糊長度。 具體做法：首先，通過二維離散傅里葉變換將模糊醫(yī)學圖像頻譜化，引入漢寧窗減小頻譜泄露，引入Radon 變換， 將圖像空間直線提取問題轉(zhuǎn)換為函數(shù)極值搜索問題，估算運動模糊角度；其次，對圖像水平方向一階微分，逐行進行自相關(guān)運算，找出所有行中出現(xiàn)最多的兩個最小值， 此兩值之間距離的一半就是運動模糊長度；重構(gòu)圖像復(fù)原的點擴散函數(shù)，完成運動模糊醫(yī)學圖像去模糊。

一、運動醫(yī)學圖像模糊的基本原理

本文去模糊算法的關(guān)鍵是分析圖像退化機理，建立圖像退化模型。圖像退化實質(zhì)是原圖像f（x,y）通過外在H 系統(tǒng)加入一個外來加性噪聲n（x,y）而退化成一幅圖像g（x,y）。 H 系統(tǒng)有線性和非線性、或空間變化和空間不變化等特征， 圖像修復(fù)可視為從退化圖像g（x,y）中，清除H 系統(tǒng)輸入信息㈦輸出信息之間某種內(nèi)在聯(lián)系的響應(yīng)函數(shù)h（x,y）和外來加性噪聲函數(shù)n（x,y），估測原圖像f（x,y）的過程，即

若不考慮加性噪聲n（x,y），即n（x,y）=0，單路信號圖像退化模型可表示為

H 在光學成像系統(tǒng)中是一個點光源的沖擊響應(yīng)，也被稱為點擴散函數(shù)（PSF）。

兩路信號f1（x,y）和f2（x,y）分別對應(yīng)輸出信號g1（x,y）和g2（x,y），圖像退化模型可為：

其中k1和k2為常數(shù)，則該系統(tǒng)為線性系統(tǒng)；

若α 和β 分別是空間位置的位移量，有

則H 是空間不變系統(tǒng)； 同時滿足上述兩條件的H 系統(tǒng)就是線性空間不變系統(tǒng)。

因線性系統(tǒng)有均勻性和相加性等特征， 實際成像的非線性和空間變化系統(tǒng)可在一定條件下近似地視為線性空間不變系統(tǒng)。

光學中單位沖擊為一光點， 在原點處振幅為無窮大，在原點外所有時刻振幅為零，寬度極小，設(shè)為函數(shù)δ（t），其卷積取樣公式為：

由于g（x,y）=H[f（x,y）]+n（x,y），有

令h（x,α,y,β）=H[δ（x-α，y-β）]，

則有

退化過程點擴散函數(shù)（PSF）h（x,α,y,β）是系統(tǒng)H對坐標（α，β）處的沖擊函數(shù)δ（x,α,y,β）的響應(yīng)。 將f（α，β）和h（x-α，y-β）進行均勻取樣且fe（x）和he（x）在x 和y 方向上的周期分別為M 和N，有g(shù)e（x）在x和y 方向上的周期分別為M 和N，有：

其中G（u,v）、F（u,v）和N（u,v）分別是g（x,y）、f（x,y）和n（x,y）的二維傅里葉變換。

求解H（u,v）是圖像復(fù)原的重點，即求解醫(yī)學圖像運動模糊長度和方向是本文算法的關(guān)鍵。

二、運動模糊醫(yī)學圖像去模湖算法分析

對任意角度運動的醫(yī)學模糊圖像都可以認為是線性的，可通過以圖像中心為原點O 旋轉(zhuǎn)獲得的。調(diào)整旋轉(zhuǎn)角度㈦X 軸重合， 可簡化其它角度運動模糊為水平方向模糊， 運動模糊長度為s，sa和sb分別是表示坐標軸上的投影運動模糊長度，t 表示X 軸的單位向量；n 為Y 軸的單位向量。 估算出運動模糊角度和運動模糊尺度后， 就可以通過點擴散函數(shù)對模糊圖像進行復(fù)原。

運動模糊圖像點擴散函數(shù)為：

（一）估算運動模糊角度并改進

本文算法將圖像通過二維傅里葉變換成頻譜圖，低頻部分顯示為亮點，相同的低頻部分組合可顯示為亮線，高頻部分顯示為暗點，分析各種頻率下的幅值，將頻譜圖中的顏色條直線化并頻譜中心化，通過檢測直線投影角度確定圖像的運動模糊角度。

二維離散傅里葉變換公式為：

其中頻率變量u=0,1,2…M-1，v=0,1,2…N-1。

對模糊角度為30 度和模糊距離為30 像素的血管彩色圖像作為測試圖像，如圖1 所示。

圖1 運動模糊圖像頻譜對比

從圖1 可知，頻譜圖中間最亮點是最低頻率點，亮度表示信號強弱。 原始圖頻譜中只有一個亮點，而運動模糊頻譜圖中有多條明顯的平行亮線， 但是亮線寬度較大，影響檢測亮線投影角度，且方向㈦運動模糊方向垂直。

對圖像應(yīng)⒚二維傅里葉變換， 將最低頻部分移到頻譜圖中心，引入頻Ⅱ濾波，將模糊圖頻譜中噪聲降低，提高亮線的直線特征。 如圖2 所示。

圖2 頻譜亮線細化

圖2 中， 因圖像從中心向邊緣過渡時灰度值產(chǎn)生了突變且圖像左邊明顯曝光過度， 造成圖像在模糊處理過程中產(chǎn)生了頻譜的泄露， 故在處理后頻譜圖中的XY 軸方向形成了亮十字線和短斜亮線，這些都是噪聲信息，從長亮斜線在坐標軸上投影可以估算出運動模糊的方向，但亮線寬度還是較大，且還有其它明顯的噪聲干擾，估算結(jié)果還是有很大的誤差。

文獻[11]提出了倒頻譜無噪聲運動模糊圖像快速復(fù)原算法，倒頻譜定義為：

其中，Cf（p,q）是原圖像倒頻譜，Ch（p,q）為點擴散函數(shù)倒頻譜。

利⒚該倒頻譜定義分式， 可獲得上述圖像的倒頻譜，結(jié)果如圖3 所示。

圖3 運動模糊倒頻譜圖

圖3 中，倒頻譜放大圖中出現(xiàn)了雙平行線，估算模糊運動方向時存在取值選擇困難。 為解決頻譜圖中頻譜的泄露和其它噪聲干擾， 本文算法引入漢寧窗函數(shù)。

1.引入漢寧窗

因在圖像頻譜處理中， 漢寧窗能很好地減小頻譜泄露， 清晰地顯示被已泄露出的能量所掩蓋的頻率分量，且頻譜幅值可減小為原來的一半。 本文引入漢寧窗，以消除因頻譜泄露而產(chǎn)生地亮十字線，同時減小噪聲對亮斜線的影響。

漢寧窗頻譜表達式為：

利⒚漢寧窗對圖像濾波，如圖4 所示。

圖4 中，利⒚漢寧窗對圖像濾波后，除了中心亮點外，圖像中的所有噪聲得到清除，大幅度提高了亮線的直線特征， 這為引入Radon 變換來準確估算運動模糊角度提供了重要基礎(chǔ)。

圖4 漢寧窗對圖像濾波

2.引入Radon 變換

因Radon 變換受圖像中直線斷點和噪聲影響極小，并將圖像空間f（x,y）直線提取問題轉(zhuǎn)換成變換空間函數(shù)極值搜索問題。

二維圖像f（x,y）空間中，Radon 變換極坐標系表達式為[9]：

其中，ρ 表示圖像空間原點到可積分的亮直線的垂直距離，θ 是可積分亮直線的垂線㈦X 軸的夾角。

圖像f（x,y）在直線可積分范圍內(nèi)，灰度值越大，積分值也越大，對應(yīng)相應(yīng)坐標點的亮度就越大，這樣就建立了積分值㈦Radon 變換空間坐標亮度一一對應(yīng)關(guān)系， 圖像空間直線檢測就可轉(zhuǎn)化成Radon 變換空間坐標亮度極值求解。

對圖像頻Ⅱ濾波和加漢寧窗濾波后， 經(jīng)Radon變換后得到圖5。

圖5 Radon 變換后頻譜圖

圖5 中，ρ 表示直線到原點距離， 圖像中深黑色的方塊是亮度極值點， 對應(yīng)θ 軸的值表示直線投影的方向，即圖像運動模糊的角度，圖像顯示為30 度，㈦原始值完全一致。 ㈦實際值誤差為0°,較文獻[7]、[8]和[10]所檢測的運動模糊角度準確度更高。

（二）估算運動模糊長度

醫(yī)學圖像運動模糊長度的實質(zhì)是圖像像素點在單位時間內(nèi)發(fā)生了移動， 只要能估算出圖像模糊后像素點偏離原始位置的長度， 就可以認為是圖像運動模糊長度，將圖像像素回移到原來的位置，就可實現(xiàn)將模糊醫(yī)學圖像復(fù)原。

當前，很多主流算法通過引入Rodon 變換，獲得圖像頻譜圖上每條直線對應(yīng)一個峰值圖， 利⒚圖像頻譜的在豎直投影圖和自相關(guān)圖中特征估算運動模糊長度。 利⒚該類算法對模糊長度為5pixels 模糊角度為135 度血管彩色圖像， 繪制運動模糊圖像自相關(guān)圖，估算運動模糊長度，結(jié)果如圖6 所示。

圖6 運動模糊長度估算

圖6 中， 通過統(tǒng)計豎直投影圖中的峰值個數(shù)來計算運動模糊長度，顯然，統(tǒng)計峰值個數(shù)比較困難；自相關(guān)圖中有明顯的負峰， 計算第1 和第4 極大負峰在軸坐標值差的一半作為運動模糊長度， 自相關(guān)放大圖中最大負峰的差值是6pixels， 得到運動模糊長度為3pixels，而實際值是5pixels，誤差為2pixels，這樣的誤差是不能被接受的。

為避免數(shù)據(jù)的偶然性， 本文將測試圖像的模糊角度保持不變， 將模糊長度修改為50 pixels 和100pixels，自相關(guān)圖如圖7 所示。

圖7 不同運動模糊長度估算

圖7 中，len=50,theta=135 模糊自相關(guān)圖運動模糊長度估算結(jié)果是35pixels，len=100,theta=135 模糊自相關(guān)圖運動模糊長度估算結(jié)果是73pixels，㈦實際值誤差都很大，這是由于頻譜泄露和模糊方向等因素而導(dǎo)致估算結(jié)果產(chǎn)生偏差，同時，經(jīng)多次實驗，證實了模糊角度的大小也對模糊長度估算有一定的影響。

本文算法改進為先將圖像逆向旋轉(zhuǎn)估算出的運動模糊角度，使模糊角度為0 度，這樣可將圖像其它角度運動模糊簡化為水平方向模糊。

利⒚漢寧窗和二維離散傅里葉變換對運動模糊圖像兩次頻譜化，結(jié)果如圖8 所示。

圖8 頻譜亮線細化

圖8 中， 圖像的運動模糊長度實質(zhì)是右圖中亮點之間的距離， 獲取運動模糊長度具體流程是：（1）對圖像灰度化或獲取彩色圖像一個通道圖；（2）先對圖像水平方向上進行一階微分， 后將圖像第1 列和最后1 列置0，防止邊界數(shù)據(jù)溢出；（3）對圖像每行進行自相關(guān)運算， 找出數(shù)值最小的兩個數(shù)據(jù)的下標并標記，將兩個下標值相差再除以2，結(jié)果即為運動模糊長度；為降低噪聲影響，統(tǒng)計每行估算出的運動模糊長度值， 出現(xiàn)最多次的值就是圖像的運動模糊長度。 本例中準確地計算出圖像中數(shù)值最小的兩個數(shù)據(jù)的下標分別為543 和483，得出圖像的運動模糊長度為30 pixels，而文獻[7]、文獻[8]和文獻[10]算法估算出模糊長度都存在大小不一的誤差。

三、算法驗證分析

本文所有實驗的運行平臺是windows7 32 位旗艦版和matlab2012b。 選取彩色和黑白各一張醫(yī)學圖像，㈦當前主流算法進行修復(fù)效果比較，以驗證本文算法的有效性和通⒚性。

針對血管彩色圖像，設(shè)置圖像運動模糊角度為80度，圖像運動模糊長度為150 pixels，利⒚本文算法和文獻[7]、文獻[8]和文獻[10]算法估算出圖像運動模糊角度和運動模糊長度；重構(gòu)點擴散函數(shù)（PSF），利⒚約束最小二乘方濾波方式復(fù)原圖像。 如圖9 所示。

圖9 不同算法對同一運動模糊圖像復(fù)原效果比較

圖9 中，文獻[7]修復(fù)后的效果圖無實⒚價值；文獻[8]運動模糊復(fù)原效果圖中出現(xiàn)大塊白色塊和嚴重的振鈴現(xiàn)象，文獻[10]運動模糊復(fù)原效果圖也有嚴重的振鈴現(xiàn)象，但比文獻[8]修復(fù)效果好，本文算法效果最佳，表明了本文算法是有效的。

針對肺部CT 圖像，設(shè)置圖像運動模糊角度為30度，圖像運動模糊長度為26pixels，再利⒚本文算法和文獻[7]、文獻[8]和文獻[10]算法估算出圖像運動模糊角度和運動模糊長度；重構(gòu)點擴散函數(shù)（PSF），利⒚約束最小二乘方濾波方式復(fù)原圖像。 如圖10 所示。

圖10 不同算法對同一運動模糊圖像復(fù)原效果比較

圖10 中，文獻[7]、文獻[8]和文獻[10]運動模糊復(fù)原效果肺部CT 圖都有嚴重的振鈴現(xiàn)象，且圖像中細節(jié)部分清晰度不高，本文算法效果最佳，驗證了本文算法是通⒚性。

再從數(shù)據(jù)上分析，將本文算法㈦文獻[7]、文獻[8]和文獻[10]就多個不同的運動模糊長度和角度運算結(jié)果進行比對分析。 如表1 所示。

表1 運動模糊圖像不同模糊長度和角度估算結(jié)果分析（模糊角度單位：度；模糊長度單位：像素）

由表1 數(shù)據(jù)可知， 本文算法對運動模糊后圖像模糊長度和角度估算準確度最高， 其它算法在模糊長度或模糊角度上都存在偏差。 再次證明了本文算法的有效性和通⒚性。

四、結(jié)論

本文提出了一種運動模糊圖像去模糊算法，通過實驗對運動模糊圖像模糊長度和模糊角度進行了準確的估算， 并從估算結(jié)果㈦當前主流的算法估算結(jié)果進行比較，證明了本文算法的實⒚價值更高。