任意形狀曲線刃邊的點擴散函數(shù)估計方法

李真?zhèn)ィ迖遥鶑闹蓿瑒?陽

1. 信息工程大學基礎部，河南 鄭州 450001； 2. 上海工程技術大學管理學院，上海 201620

光學系統(tǒng)的點擴散函數(shù)(point spread function,PSF)的精確估計是提高圖像光學質量的必要前提[1]之一。通常，以h(x,y)代表點擴散函數(shù)矩陣。很多情形下，只能利用成像系統(tǒng)的參數(shù)或獲取的圖像對系統(tǒng)退化函數(shù)進行辨識以獲取準確的h(x,y)。常利用遙感圖像中的刀刃目標、線性脈沖目標及方波目標來獲取準確的h(x,y)[2]。但是這些傳統(tǒng)方法均要求嚴格控制靶標的形狀及擺放的角度[3-4]。為了克服對靶標的依賴，很多方法，例如多通道盲反卷積[5-8]、基于聯(lián)合變換的方法[9]、基于全變分正則化盲復原的方法[10]被提出并進行研究測試。這些方法雖然時間開銷小，但是效果很一般。文獻[11]建立了帶有參數(shù)的高斯型點擴散函數(shù)模型，調整參數(shù)獲得合適的h(x,y)；文獻[12]提出了一種利用雙刃邊的方法估算運動模糊，但這些方法適用范圍均較小。近年能適應更多場景的利用刃邊法進行h(x,y)估計的研究逐漸展開[13-15]，主要針對傾斜直線刃邊[16]以及曲線刃邊的情況。

對于直線刃邊，傳統(tǒng)的刃邊法只限于應用在理想刃邊區(qū)域[17]。文獻[18]給出了對齊拐點的方法，但是對于刃邊傾角較大的情況，由于h(x,y)采樣方向與刃邊梯度方向不一致，樣本橫坐標被拉伸，會產生不可避免的誤差[19]。文獻[20]通過模擬離散卷積過程建立樣本h(x,y)與真實h(x,y)的關系方程組，并通過最小二乘法求解真實h(x,y)。雖然該方法從離散的角度在原理上解決了傾斜直線刃邊下獲取真實h(x,y)的方法，但是在無法知道h(x,y)準確尺寸的情況下無法建立樣本h(x,y)與真實的h(x,y)之間關系的方程組。此外，該方法對刃邊邊緣準確度的要求極其苛刻并且對噪聲異常敏感，而采用的離散模型會產生亞像素點位誤差，直接影響h(x,y)的估計精度。傾斜刃邊法(投影法)[21]是一種較為成熟的方法。該方法將像元點投影到刃邊垂直方向上，以獲得較為準確的h(x,y)的采樣數(shù)據。該方法解決了亞像素點位誤差的問題，同時在投影前不需要知道h(x,y)的準確尺寸，是目前較為有效的任意傾角的直線刃邊點擴散函數(shù)估計方法。

盡管最近幾年提出了很多種解決彎曲刃邊獲取h(x,y)的方案，但是大部分方法的思路可以歸結為兩種：基于ISO12233[22]及其一系列改進的方法(傳統(tǒng)刀刃法)；將像元點灰度值及點到彎曲刃邊的距離組成集合后進一步處理的方法(曲線擬合法)。文獻[23]直接通過對邊緣擴散函數(shù)(edge spread function,ESF)曲線族微分得到的線擴散函數(shù)(line spread function,LSF)進行曲線族的對齊及樣本重采樣，最終得到唯一的LSF曲線。但只有在曲線每個點的切線都滿足理想傾角時該方法才較為準確，否則會使橫坐標拉長導致h(x,y)測量不準確。文獻[22]采用二次多項式擬合彎曲刃邊，將像元點灰度值以及點到彎曲刃邊的距離值加入距離—灰度集合中，最后用集合中的元素擬合Fermi函數(shù)作為實際的ESF曲線。一方面，該算法使用了模擬退火算法導致迭代次數(shù)過多，時間開銷大。另一方面刃邊邊緣很難用二次函數(shù)完美擬合。事實上，將像元點的灰度值以及點到彎曲刃邊距離作為ESF樣本的參考本身就值得商榷，這將在后面指出。一般的彎曲刃邊圖像分為左右型和上下型(圖1(a)、(b))，為方便起見，將默認以左右型進行討論，上下型的刃邊以此類推。

圖1 刃邊圖的分類Fig.1 Classification of edges

為降低對靶標的要求，針對彎曲刃邊處理較為困難的現(xiàn)狀，本文提出了一種彎曲刃邊的點擴散函數(shù)估計方法。對于曲率在0.001～0.01的刃邊，本文算法在較強模糊下，峰值信噪比依舊能保持在35 dB以上，測得點擴散函數(shù)峰值的誤差可以控制在20%以內。試驗結果與現(xiàn)有較新文獻進行了比較，提出的算法峰值信噪比平均能提高10 dB以上，且具有一定抗噪能力。在很大程度降低了對遙感圖靶標的要求，可以適應各種復雜的彎曲刃邊環(huán)境。

1 刃邊法點擴散函數(shù)估計原理

設清晰圖像為e(x,y)，退化圖像為o(x,y)，遙感圖像退化模型可以描述為[24-25]

o(x,y)=e(x,y)*h(x,y)+n(x,y)

(1)

式中，h(x,y)為點擴散函數(shù)矩陣；n(x,y)為加性噪聲；*表示卷積運算。

一般情況下h(x,y)具有可分離性[26]

h(x,y)=hx(x)hy(y)T

(2)

式中，hx(x)是水平采樣方向的h(x,y)，hy(y)為垂直采樣方向上的h(x,y)。hx(x)和hy(y)被稱為LSF。一般情況下，假設h(x,y)是對稱的，即hx(x)和hy(y)默認相等。如果假設U(x)為垂直于刃邊方向上的階躍函數(shù)，有

g(x)=U(x)·hx(x)

(3)

式中，g(x)為該方向上的輸出函數(shù)，即ESF，其與LSF的關系是

(4)

理想傾角下刃邊法獲取PSF使用ISO12233傾斜刃邊法計算PSF[27](見圖2)。基本思想是利用采樣得到的灰度圖計算得到ESF樣本后微分得到LSF樣本，也就是一維點擴散函數(shù)。但是該傳統(tǒng)刃邊方法僅僅適用于采樣間隔為一像元單位的理想刃邊區(qū)域。

圖2 刃邊法Fig.2 Edge method

2 算法原理

2.1 投影法原理

本文基于投影法提出一種點擴散函數(shù)估計的彎曲刃邊法。很多文章中運用投影法在傾斜直線刃邊下估計h(x,y)，但是很少有文章說明投影法的理論依據。本文從傾斜直線出發(fā)分析投影法內在的物理意義。圖像的退化可以表示成離散卷積

(5)

式中，g為退化的圖像；f為原圖像。由式(2)假設PSF分解后長度為5的LSF為hh=[u1,u2,u3,u2,u1]。45°的傾角的刃邊邊緣見圖3(a)，圖中深色的部分表示灰度值較高的部分。其經過上述PSF退化后的像元如圖3(b)所示。

由于顏色區(qū)分度不高，為清楚地表示，圖中的數(shù)字是將灰度值從小到大排序之后的結果，數(shù)字越大，代表灰度值越高。圖中的直線是理論上的刃邊邊緣。其中-3位置的點的灰度值為

g(-3)=SUM(PSF·

(6)

-4位置的點的灰度值為

(7)

離散條件下有以下結論：

在離散條件下，退化后像元的灰度值等于退化前以該點為中心的h(x,y)大小的區(qū)域的灰度值與h(x,y)對應位置值點乘后的和(見圖4(a))。

圖3 退化前后的45°傾角刃邊Fig.3 Edge of 45° before and after degeneration

在離散條件下投影法的原理可以理解為：保持像元點到刃邊邊緣的距離不變，退化前以該點為中心的h(x,y)大小的區(qū)域的灰度值與h(x,y)對應位置值點乘后求和的結果不變。進一步，由于h(x,y)的值以及退化前刃邊邊緣兩端灰度值均為常數(shù)，可以理解為：保持像元點到刃邊邊緣的距離不變，退化前以該點為中心的h(x,y)大小的區(qū)域內灰度值的和保持不變(見圖4(b))。但是在離散的情況下投影法的理論依據并不滿足，直觀上可以從圖4(b)中觀察到，保持像元點到刃邊邊緣的距離h不變，刃邊邊緣的兩端面積卻發(fā)生了改變。

圖4 離散條件下投影法的理論依據Fig.4 The theory evidence of projection method under discrete conditions

造成這種矛盾的原因是真實的h(x,y)以方陣來描述是建立在方陣足夠大的前提下的，更準確的描述是以某點為圓心向外擴散的圓，并且h(x,y)上的點的值僅與該點距離h(x,y)中心的距離有關(見圖5(a))。在真實亞像素條件下，退化后亞像素像元點的灰度值g為

(8)

式中，Ω1、Ω2分別為退化前以像元點為中心的h(x,y)大小的圓所覆蓋位置被刃邊邊緣所分割的兩塊區(qū)域(見圖5(b))；g1、g2分別為退化前Ω1、Ω2區(qū)域所對應的灰度值；F代表歸一化后h(x,y)各位置灰度值函數(shù)。

在直線刃邊的情況下，由于刃邊同到刃邊的最短距離的直線必定垂直，所以只需要通過h(x,y)的值以及亞像素像元點到刃邊的距離，就可以由式(8)計算出退化后該像元點的灰度值g。同時，無論將刃邊旋轉到任何角度，退化前以該點為中心的h(x,y)大小的圓內Ω1、Ω2部分的積分值不會發(fā)生改變(見圖5(c))，也就是說該像元點到刃邊直線的距離可以認為是理想傾角下距離刃邊邊緣的橫坐標距離。這就是投影法進行h(x,y)樣本采樣的理論依據。

圖5 真實的亞像素條件下投影法的理論依據Fig.5 The theory evidence of projection method under real sub-pixel conditions

2.2 曲線刃邊投影法理論依據

在上節(jié)基礎上考慮彎曲刃邊時就會發(fā)現(xiàn)：單純尋找每一點到曲線刃邊的垂直距離作為h(x,y)樣本的橫坐標的參考會產生巨大的誤差，因為這個距離反映的僅僅是該距離下直線刃邊的情況(見圖6(a))，這也驗證了本文算法的問題。并且由于曲線方程難以找到合適的表達式，很難通過計算來確定這些覆蓋的部分。事實上，只需要得到精確的部分h(x,y)樣本，并不需要每個退化后的點都為h(x,y)樣本的采樣點。如果能夠選出曲線刃邊中較為平緩(曲率半徑較大)的部分周圍的點(見圖6(b))，用直線刃邊投影代替曲線刃邊投影，在這些點處，點到刃邊直線的距離就可以近似反映出該點為中心的h(x,y)尺寸的圓所覆蓋的部分，進一步就可以通過插值和重采樣較為精確地得到較好的ESF樣本。

基于這種思想，本文提出了移動窗口投影的曲線刃邊求解方法，每次選取一定的行擬合邊緣曲線為直線并采用任意傾角直線刃邊投影法進行h(x,y)樣本的采樣，然后向下移動步長個采樣行，繼續(xù)h(x,y)樣本的采樣，直到完成對所有樣本行的采樣(圖6(c))，最后將那些明顯不適合作為h(x,y)樣本的采樣點去除后得到最終h(x,y)樣本并進行進一步的處理。

圖6 曲線刃邊投影法理論依據以及算法Fig.6 The theory evidence of projection method of curve edge and the schematic diagram of algorithm

3 算法步驟

利用移動窗口投影法，對左右型刃邊的PSF計算的總體流程見圖7。本節(jié)將詳細說明流程中的關鍵步驟。

圖7 算法流程Fig.7 The flow chart of algorithm

3.1 ESF樣本的采樣

ESF樣本的采樣是彎曲刃邊法的核心部分。目標圖像的3個超參數(shù)：移動窗口即每次采樣選取的行數(shù)大小L，每次窗口移動的步長S，每次的投影中心C的選取會在試驗部分給出。

單次需要處理的部分就是某一窗口中選取的行，見圖6(c)。先采用最小二乘法對選定行數(shù)的邊緣點的坐標進行直線擬合。其次按照任意傾角刃邊投影法對選取部分進行ESF樣本的采樣：設擬合直線與垂直方向傾角為θ(見圖6(c))，θ為小于等于45°的任意角度。采樣的部分某點的坐標為(i,j)，其中i為采樣點所在列數(shù)，j為所在行數(shù)，則該點投影位置的坐標P=j*cosθ±i*sinθ。(對于左下到右上的直線，符號取負，對于左上到右下的直線符號取正)。

采樣完單次的ESF樣本之后將窗口向下移動S個采樣行繼續(xù)采樣，直到完成對所有樣本行的采樣。在每個窗口內計算投影位置坐標應默認邊緣點擬合得到的直線上的點位于投影的中心位置，所以在最終得到的ESF樣本中需要將投影中心對齊。為此將所有點的投影位置坐標P加上一個偏差值，設每次窗口中第一行刃邊邊緣對應的坐標為(io,jo)，偏差值為C-(jo*cosθ±io*sinθ)，此時對應窗口中某點(i,j)的投影位置坐標值就應當修正為P*=j*cosθ±i*sinθ+C-(jo*cosθ±io*sinθ)。將P*和對應點灰度值作為一組ESF的采樣值存儲。值得注意的是同一行的數(shù)據可能會按照不同的傾角采樣多次，但是在刀刃邊緣曲線連續(xù)且曲率不是太大的情況下，同一行得到的ESF樣本投影位置坐標幾乎相同。

3.2 異常與噪聲ESF樣本的剔除

圖8 容忍域Fig.8 Tolerance domain

3.3 ESF樣本的插值與重采樣

剔除粗差之后，一方面刃邊邊緣檢測存在誤差，另一方面ESF樣本中會出現(xiàn)缺失值，這將影響進一步的LSF樣本采樣。因此對剔除粗差之后的ESF樣本進行插值后重采樣的操作。具體步驟如下：第1步對剔除粗差之后的樣本進行插值(插值的目的是為了將剔除異常點和噪聲點之后的ESF樣本缺失值補上)。先按照投影位置坐標P*的大小將ESF樣本進行排序，然后采用線性插值對投影位置坐標P*的取值區(qū)間中的ESF樣本灰度值進行插值，移動步長為0.05像素。第2步將插值得到的點同原ESF樣本一起重新按照投影位置的值P*進行排序，得到新的ESF樣本。第3步對新的ESF樣本進行分段擬合，移動窗口獲取ESF重采樣值，移動窗口的尺寸為1像素，窗口內對ESF樣本進行三次多項式擬合，將窗口中心處擬合值作為重采樣值，移動步長的設定比較隨意，這里選擇0.05像素。

3.4 LSF采樣與擬合

由式(4)可知，LSF可以對ESF直接微分獲得。對于離散樣本，LSF可以通過ESF樣本差分得到

(9)

式中，L為LSF樣本；E為ESF樣本。在對LSF樣本進行擬合處理之前，算法先采用文獻[28]的方法估算了LSF尺寸的大小。文獻[28]中只考慮了4個方向，對于不是這幾個方向的直線會出現(xiàn)偏差。提出的算法中處理的數(shù)據是投影差分后的LSF樣本，因此并不會出現(xiàn)文獻[28]中的問題。噪聲會干擾輪廓值(LSF尺寸)的準確性，因此本文設定了一個閾值加以限定。為確定輪廓大小，從左端(右端)開始遍歷獲得的LSF樣本(采樣間隔為0.05個單位長度)。當某一點連同其相鄰的兩點都大于某個閾值的時候將其坐標記錄為端點1(端點2)。閾值視噪聲大小而定。由端點2坐標減去端點1坐標并向下取整得到輪廓值大小。

下面對LSF進行尾部截取以減弱LSF尾部抖動影響，為了保證LSF擬合的準確性，最終截取部分的長度應小于之前得到的輪廓值的大小。LSF的形狀近似于高斯型，為了減少估計的誤差，對得到的LSF樣本進行一維高斯函數(shù)擬合。擬合模型可以表示為

h(x)=Aexp[-(x-u)2/2σ2]

(10)

式中，A為擬合后的峰值；u為擬合中心點的位置；x為LSF的位置；σ為擬合高斯函數(shù)的標準差。可依據最小二乘原理對上述模型進行求解。最后以u作為截取的中心截取輪廓值的大小得到最終的LSF向量。由式(2)，h(x,y)可以通過兩個方向的LSF向量乘積得到，歸一化之后就是最終的h(x,y)。

4 試驗測試

4.1 參數(shù)的選取

試驗測試需要確定3個超參數(shù)，分別是投影的中心位置常數(shù)項C，每次采樣的行數(shù)L以及框架移動的步長S。投影中心的常數(shù)項C選擇相對隨意，僅需要保證每個樣本投影后的坐標為正即可。在比較刃邊算法性能的時候，常用峰值信噪比(peak signal to noise ratio,PSNR)作為算法性能的客觀評價指標[29]。每次選取的行數(shù)大小L與圖片或者h(x,y)的尺寸有著很大的關系。本文先以128×128的模擬圖為例(見圖9(b)和圖9(c))，以選取的行數(shù)為自變量，以所估h(x,y)與真實值的峰值信噪比作為因變量，最終得到的結果如圖9(a)所示。圖中的3條曲線分別是以2作為步長h(x,y)尺寸為5×5，9×9，以及15×15的每次取樣行數(shù)與峰值信噪比的關系圖。圖中可以看出，最優(yōu)的選取行數(shù)與h(x,y)的大小無關，每次采樣行數(shù)在10左右的時候峰值信噪比會比較穩(wěn)定，采樣行數(shù)過小或者過大都會對結果產生影響。確定每次采樣的行數(shù)與h(x,y)的尺寸無關以后本文又以64×64的模擬圖進行同樣的試驗，最優(yōu)選取行數(shù)在5左右。由此最優(yōu)選取行數(shù)一般為采樣圖片尺寸的十分之一左右。

圖9 參數(shù)選取及模擬圖Fig.9 Parameters selection and simulation diagram

4.2 不同曲率刃邊試驗結果

為驗證本文方法在理論上的可行性及測試算法的極限，本文針對不同曲率的曲線進行試驗，利用半徑為R的圓曲率為1/R的性質，可以在理論上構建任何曲率的曲線，但是一方面受到h(x,y)尺寸的限制曲率不能過大，一方面曲率過小曲線將接近直線，這里僅僅考慮曲率在0.001～0.01的情況：對不同曲率的曲線圖分別進行尺寸為15×15,方差為0.5和1的高斯濾波，PSNR與曲率的關系圖分別見圖10的模擬結果1和2。在結果2中選取曲率為0.07的刃邊，獲得的LSF曲線與真實LSF曲線的比較圖見圖10模擬結果3，獲得的LSF曲線與真實LSF曲線基本吻合。對于曲率在0.001～0.01的刃邊，提出的算法在方差為0.5的較弱模糊下，PSNR保持在40 dB以上，而在方差為1的較強模糊下，PSNR依舊能保持在35 dB以上，測得PSF峰值的誤差可以控制在20%以內。

圖10 模擬試驗結果Fig.10 Results of simulation experiment

4.3 不同曲線刃邊算法的比較

選取不同的曲線刃邊以及方法進行對比試驗，分別以兩個角度的曲線刃邊圖(如圖9(b)(c))對傳統(tǒng)刃邊法、二次曲線擬合法以及提出的方法測試，最終的試驗結果記錄在表1。對于兩種曲線刃邊，傳統(tǒng)刃邊法都有比較穩(wěn)定的效果，但是曲線擬合法卻在測量結果上出現(xiàn)了很大幅度的波動，一些情況下甚至會出現(xiàn)測量嚴重不準確的情況。為對比效果，對曲線刃邊2進行尺寸為15×15,方差為1的高斯濾波后的結果，分別用3種方法測得的PSF進行維納濾波。用真實h(x,y)，本文測得的h(x,y),傳統(tǒng)刃邊法測得的h(x,y)復原結果,曲線擬合法測得的h(x,y)進行維納濾波的結果分別見圖11(c)、(d)、(e)、(f)。由于測量偏差大，傳統(tǒng)刀刃法測得h(x,y)復原結果出現(xiàn)了嚴重的振鈴效應，而曲線擬合法測得的h(x,y)復原結果則和復原前模糊圖幾乎相同。提出的算法在各種情況下都要遠遠優(yōu)于其他兩種算法，PSNR平均可以在最優(yōu)效果的基礎上提高10 dB左右。

圖11 不同方法測得的PSF維納濾波復原后效果Fig.11 Images after restoration by Wiener filtering using PSF measured by different methods

4.4 算法的抗噪性研究

為測試噪聲對于算法的影響，分別對曲線刃邊1用尺寸15×15，方差為1的PSF高斯濾波后的模擬圖中加入高斯噪聲以及隨機噪聲，來測試算法的抗噪能力，結果分別見圖12(a)、(b)。經測試，對于隨機噪聲，雖然測量精度有所下降，但是仍能較準確獲得h(x,y)估計值。算法對于高斯噪聲較為敏感，對于標準差在10以上的高斯噪聲無法獲得預期的效果。

4.5 遙感圖復原效果

對于真實場景，應當選取邊界為階躍變化且未受到強噪聲污染的彎曲刃邊作為靶標圖，同時應當保證刃邊曲率在0.1以下。本節(jié)中選取某衛(wèi)星遙感圖一部分進行測試，對場景進行大小為9×9，方差為1的高斯濾波后，截取中間64×64大小的作為靶標，采用本文方法進行彎曲刃邊法測試如圖13(a)，最終測得的LSF曲線與真實LSF曲線對比圖見圖13(b)。圖13(c)和圖13(d)分別是退化圖采用算法估計的h(x,y)進行維納濾波前后截取的同一部分的結果。可以看到，圖像所能獲得的刃邊區(qū)域邊緣經常會有不連續(xù)的點，在這種情況下，算法仍可以較好地估計出真實的h(x,y)，并且獲得不錯的復原效果。

圖13 遙感圖像測試Fig.13 Test of remote sensing images

5 結 論

針對大量靶標受到限制條件下的遙感圖像估計PSF的問題，本文提出的方法可以利用任意形狀的刃邊進行PSF的估計，從理論上分析了現(xiàn)有算法的不足以及提出方法的可行性。對于曲率在0.001～0.01的刃邊，本文算法在較強模糊下，PSNR依舊能保持在35 dB以上，測得PSF峰值的誤差可以控制在20%以內。相比傳統(tǒng)刃邊法和曲線擬合法，本文方法測量結果的PSNR平均可以在最優(yōu)效果的基礎上提高10 dB左右，對于實際遙感圖有不錯的效果。算法對隨機噪聲有一定抗性，對于高斯噪聲較為敏感，只適用于標準差在10以下的場景，有待在未來的工作中改進。