作為光學成像系統中最常見的簡化模型,高斯型點擴展函數(PSF)參數的估計算法是圖像復原領域中人們所研究的熱點問題之一。本文基于文獻中高斯型OTF曲線特征點與OTF參數及圖像頻域表達存在著確定數學關系的思想,提出了一種自動、準確、穩定地估計高斯型PSF參數的改進方法。該方法通過對真實圖像實施高斯卷積模糊退化運算,突出相關計算特征,最終實現顯微圖像PSF參數的自動估計。實驗表明,采用該方法所估計的PSF可以獲得良好的尿沉渣圖像復原效果,對于具有高斯型近似PSF模型的其他類型圖像復原和三維顯微圖像復原研究也有相當的應用和參考價值。
引用本文： 陳兆學, 陳浩. 顯微圖像高斯型點擴展函數參數估計方法. 生物醫學工程學雜志, 2014, 31(1): 53-57. doi: 10.7507/1001-5515.20140011 復制
引言
在顯微圖像成像過程中，一些隨機因素，如光學系統的像差、成像設備與物體之間的相對運動等均可能造成圖像質量的退化。針對不同的情況，可以通過建立合適的點擴展函數(point spread function,PSF)，也即PSF退化模型，從而對退化模型參數進行辨識，最終實現圖像的復原。通常圖像的退化和復原基本流程如圖 1所示[1]。其中f(x,y)代表一幅靜止、二維的顯微圖像，通過退化系統h(x,y)以后，在噪聲n(x,y)的干擾下，其退化圖像為g(x,y)，復原后的圖像為。根據線性系統理論，圖 1所示退化/復原模型可表達為：
| $g\left( x,y \right)=f\left( x,y \right)*h\left( x,y \right)+n\left( x,y \right)$ |
圖1
圖像退化/復原的模型
Figure1.
Image degradation/restoration model
進一步根據卷積定理，空間域上的卷積等同于頻域上的乘積,因此式(1) 所表達的模型可以寫成如下等價的頻域描述，即
| $G\left( u,v \right)=F\left( u,v \right)H\left( u,v \right)+N\left( u,v \right)~,$ |
其中G(u,v)、H(u,v)、F(u,v)和N(u,v)分別為g(x,y)、h(x,y)、f(x,y)和n(x,y)的傅里葉變換。
上式中的H(u,v)也就是通常所謂OTF光學傳遞函數。由于高斯函數是光學成像系統中最常見的PSF簡化模型[2-5]，若假定顯微光學成像系統的H(u,v)是高斯型的，則有
| $H\left( u,v \right)=\frac{1}{\sqrt{2\pi {{\sigma }^{2}}}}{{e}^{\frac{-{{D}^{2}}\left( u,v \right)}{2{{\sigma }^{2}}}}},$ |
其中D(u,v)是距傅里葉變換原點的距離，σ表示曲線擴展的程度。對于具有高斯型OTF的圖像復原研究，主要是能夠找到式(3) 中參數σ，然后進一步根據圖像的退化模型實現圖像的復原。
付煒等[5]根據高斯型OTF曲線的特點,提出了一種簡單、有效的高斯型OTF參數的估計方法。其基本思想是根據高斯型OTF曲線中存在的幾個特征點與降質圖像傅里葉變換的確定數學關系，推導求出OTF參數,進而實現圖像復原。但是，其在仿真實驗過程中直接給出的對Lena圖像附加人為降質操作后所得模糊圖像的實驗結果，未考慮Lena圖像成像時已有系統OTF的存在與影響。而且文中只仿真了高斯型OTF參數，也即式(3) 中σ值較大的情況。雖然其實驗結果和結論基本合理，但顯然不符合實際應用情況。此外，其降質圖像傅里葉變換特征圓半徑的計算需要通過人工觀察，并基于手動操作獲得，難以滿足實際應用過程中客觀、穩定、準確的要求。
本文在付煒等高斯型OTF參數估計基本思想的基礎上，受其實驗過程啟發，研究了人為附加高斯卷積模糊退化操作后的結果圖像和與原圖像所對應系統OTF參數的數學關系，提出了一種適應于實際應用要求、自動、準確、穩定的估計高斯型PSF參數的改進方法，并將其應用于尿沉渣顯微圖像復原，取得了較為良好的效果。本文所提出地改進方法對于具有高斯型近似PSF模型的其它類型圖像的復原和三維顯微圖像復原研究也具有應用和參考價值。
1 PSF參數的自動估計方法
針對高斯型PSF的表達式(3) ，令σ=D0，文獻[5]考察了不同D
1.1 特征圓增強的基本原理
假設為空域高斯型PSF表達模型，n(x,y)為服從N(0,σn2)的白噪聲，即。令，若對式(1) 所表達退化模型進行一次高斯模糊過程，根據卷積性質有[6]：
| $\begin{align} & b\left( x,y \right)=g\left( x,y \right)*{{G}_{b}}\left( x,y \right)=f\left( x,y \right)*h\left( x,y \right)+n\left( x,y \right)* \\ & {{G}_{b}}\left( x,y \right)=f\left( x,y \right)*\left[ h\left( x,y \right)*{{G}_{b}}\left( x,y \right) \right]+n\left( x,y \right)*{{G}_{b}}\left( x,y \right) \\ \end{align}$ |
實際上，式(4) 中h(x,y)*Gb(x,y)和n(x,y)*Gb(x,y)本質上都是兩個高斯函數的卷積，根據高斯函數的性質其結果仍然是高斯函數，即有式(5) ：
| $\begin{align} & h\left( x,y \right)*{{G}_{b}}\left( x,y \right)=\frac{1}{2\pi \left( \sigma _{b}^{2}+\sigma _{h}^{2} \right)}{{e}^{-\frac{{{x}^{2}}+{{y}^{2}}}{2\left( \sigma _{b}^{2}+\sigma _{h}^{2} \right)}}} \\ & n\left( x,y \right)*{{G}_{b}}\left( x,y \right)=\frac{1}{2\pi \left( \sigma _{b}^{2}+\sigma _{h}^{2} \right)}{{e}^{-\frac{{{x}^{2}}+{{y}^{2}}}{2\left( \sigma _{b}^{2}+\sigma _{n}^{2} \right)}}} \\ \end{align}$ |
經過高斯模糊后的圖像可以看成是被具有不同參數的OTF和噪聲模型退化的結果，其中新的OTF和噪聲空域參數分別為。
將式(5) 兩邊取傅里葉變換得OTF頻域表達為，其中。
| ${{D}_{b}}=\frac{1}{\sqrt{\sigma _{b}^{2}+\sigma _{h}^{2}}}$ |
當時，，實際上這正是文獻[5]中在實驗部分所描述的情況。
對于具有高斯型PSF的圖像復原，當其PSF所對應的標準差參數比較小，導致特征圓不明顯時，我們可以通過人為附加一具有已知較大合適參數σb的高斯卷積模糊運算，使得新圖像的傅里葉變換圖像中與OTF相關的特征圓特征變得極為明顯。圖 2所示為附加標準差為σb=2的人工模糊步驟前后，同一顯微圖像傅里葉變換幅度譜的對比，可以看出在附加人工模糊退化步驟后，特征圓特征變得相對明顯。
圖2
顯微圖像人為退化前后的幅度譜變化比較
(a) 原始顯微圖像幅度譜;(b) 人為退化后圖像幅度譜
Figure2. Comparison of the amplitude spectrum of a microscopic image between pre- and post-artificial degradation(a) the amplitude spectrum of an original microscopic image; (b) the amplitude spectrum after degraded
1.2 特征圓半徑的自動檢測
對于高斯模糊退化后進行傅里葉變換所得圖像的特征圓，通過手工方法測量其半徑后，即可計算得到新的OTF對應參數Db[5]，并由式(6) 得到原高斯型OTF參數。需要注意的是，由于傅里葉變換圖像乃頻域采樣圖像，文獻 [5]所述高斯降質和OTF參數實驗都是在頻域比較和度量的，因此在實驗過程中無需考慮頻率分辨率而只取其整數坐標度量即可。但本文涉及頻域和空域參數轉換，需要換算為真實角頻率。假設原始圖像大小為M×N，則其離散傅里葉變換圖像通常也為M×N，但其中橫向和縱向頻率分辨率分別為2π/M和2π/N，則在頻域中度量和空域真實度量之間轉換時需分別乘上相應方向的分辨率數值。
考慮到通過手工方式測量特征圓半徑時很容易引入較大的測量誤差。為了更準確地測量特征圓半徑，盡可能減少隨意性，本文給出了一種自動、穩定的精確測量方法，其基本流程如圖 3所示。
圖3
特征圓半徑自動測量計算流程
Figure3.
Flow chart of automatic radius measurement and computation of the feature circle
由于附加人為模糊后所得圖像的傅里葉變換中，特征圓內外區域對比度相對比較明顯，非常適合基于類間方差最大化原則實現二值化[7]。因此本文采用了大津算法實現圖像閾值化處理，圖 2(b)頻譜圖像的閾值化結果如圖 4(a)所示。可以看出，特征圓區域在二值化圖像中以高亮的類圓形體現，但沿橫縱兩頻率軸可看到一些細長條白色區域干擾。由于這些區域在投影圖中體現為類圓形區域位置兩側的低投影值區域，可以簡單地通過位置投影濾波方式予以排除。由于頻譜圖像二值化結果中的特征圓圓周是由很多隨機白點組成，為了穩定、有效求取特征圓半徑，特提出如下算法：
圖4
顯微圖像幅度譜中特征圓的自動檢測
(a) 幅度譜大津法二值化結果示例；(b) 特征圓檢測結果示例
Figure4. Automatic detection of the feature circle in amplitude spectrum of a microscopic image(a) example of thresholding the amplitude spectrum by Ostu method; (b) example of the result of feature circle detection
假設離散傅里葉變換圖像大小為M×N，則
(1) 求取類圓形區域的質心C(XC,YC),計算
(2) 若Error＞2，轉(7) ，否則轉(3) ；
(3) 以作為特征圓的圓心并設置半徑Rk=1，計算圖像中所有白點數目為NW；
(4) 計算二值圖像中以O為圓心Rk為半徑的圓內白點數目為NWk；
(5) 若，迭代終止；否則，Rk=Rk+1，轉(4) ；
(6) Rk即所求特征圓半徑。
(7) 報錯：特征圓不明顯，不適合用作PSF參數計算。
其最終檢測結果如圖 4(b)所示。
2 實驗與結論
為了驗證本文方法的有效性，基于尿沉渣顯微圖像，進行了相關參數測量和恢復實驗。對一原始512×512尿沉渣圖像，首先分別對其以一系列人為設置的標準差進行高斯模糊退化，然后自動測量特征圓半徑值，并進一步計算出各種人為退化和模糊標準差下對應的原始PSF參數估計值，結果如表 1所示。
表1
不同人為退化和模糊標準差下特征圓半徑和原始PSF參數估計值比較
Table1.
Comparison of feature circle radius and estimated original PSF parameter under different artificial degrading and blurring Standard Deviation
可以看出，人為附加的高斯標準差參數在0.75~2.0之間時，其PSF估計值的均值和標準偏差為1.938 55±0.027 52。可見在此區間內PSF參數估計值相對比較穩定，可以用來計算原始圖像PSF參數的估計值。用其計算所得平均值對圖像進行維納濾波復原[8]所得結果如圖 5(b)所示。可以明顯看出復原圖像對比度更好，細節更加清晰突出，在一定程度上證明了PSF參數估計的正確性。
圖5
復原前后尿沉渣顯微圖像比較
(a) 原始尿沉渣圖像；(b) 基于所估計PSF參數維納濾波復原結果
Figure5. Comparison of a urinary sediments microscopic image before and after restoration(a) original urinary sediments image; (b) the Wiener restoration result based on estimated PSF parameter
本文在文獻[5]基本原理和實驗啟發基礎上，采用附加模糊退化方法對頻譜圖像特征圓予以增強，并通過特征圓半徑自動計算實現了顯微圖像高斯型PSF參數的自動、準確和穩定估計。實驗中對于尿沉渣圖像的復原效果在一定程度上表明了本文所提出改進方法的合理性。由于受實驗條件局限，未能給出本文算法復原結果與更精密光鏡或電鏡等所成圖像的對比，今后將結合相關臨床應用進一步予以驗證。需要特別指出的是，由于高斯型OSF近似模型具有廣泛的適用性，本文OSF參數估計的改進方法不但可以適合顯微圖像的復原，還可以用于其它具有高斯型OSF近似模型的成像系統中(如天文望遠鏡系統和常規視頻攝像系統等)，實現相關圖像的復原。本文的基本思想還可應用于三維顯微圖像復原，實現三維OSF的參數估計。
引言
在顯微圖像成像過程中，一些隨機因素，如光學系統的像差、成像設備與物體之間的相對運動等均可能造成圖像質量的退化。針對不同的情況，可以通過建立合適的點擴展函數(point spread function,PSF)，也即PSF退化模型，從而對退化模型參數進行辨識，最終實現圖像的復原。通常圖像的退化和復原基本流程如圖 1所示[1]。其中f(x,y)代表一幅靜止、二維的顯微圖像，通過退化系統h(x,y)以后，在噪聲n(x,y)的干擾下，其退化圖像為g(x,y)，復原后的圖像為。根據線性系統理論，圖 1所示退化/復原模型可表達為：
| $g\left( x,y \right)=f\left( x,y \right)*h\left( x,y \right)+n\left( x,y \right)$ |
圖1
圖像退化/復原的模型
Figure1.
Image degradation/restoration model
進一步根據卷積定理，空間域上的卷積等同于頻域上的乘積,因此式(1) 所表達的模型可以寫成如下等價的頻域描述，即
| $G\left( u,v \right)=F\left( u,v \right)H\left( u,v \right)+N\left( u,v \right)~,$ |
其中G(u,v)、H(u,v)、F(u,v)和N(u,v)分別為g(x,y)、h(x,y)、f(x,y)和n(x,y)的傅里葉變換。
上式中的H(u,v)也就是通常所謂OTF光學傳遞函數。由于高斯函數是光學成像系統中最常見的PSF簡化模型[2-5]，若假定顯微光學成像系統的H(u,v)是高斯型的，則有
| $H\left( u,v \right)=\frac{1}{\sqrt{2\pi {{\sigma }^{2}}}}{{e}^{\frac{-{{D}^{2}}\left( u,v \right)}{2{{\sigma }^{2}}}}},$ |
其中D(u,v)是距傅里葉變換原點的距離，σ表示曲線擴展的程度。對于具有高斯型OTF的圖像復原研究，主要是能夠找到式(3) 中參數σ，然后進一步根據圖像的退化模型實現圖像的復原。
付煒等[5]根據高斯型OTF曲線的特點,提出了一種簡單、有效的高斯型OTF參數的估計方法。其基本思想是根據高斯型OTF曲線中存在的幾個特征點與降質圖像傅里葉變換的確定數學關系，推導求出OTF參數,進而實現圖像復原。但是，其在仿真實驗過程中直接給出的對Lena圖像附加人為降質操作后所得模糊圖像的實驗結果，未考慮Lena圖像成像時已有系統OTF的存在與影響。而且文中只仿真了高斯型OTF參數，也即式(3) 中σ值較大的情況。雖然其實驗結果和結論基本合理，但顯然不符合實際應用情況。此外，其降質圖像傅里葉變換特征圓半徑的計算需要通過人工觀察，并基于手動操作獲得，難以滿足實際應用過程中客觀、穩定、準確的要求。
本文在付煒等高斯型OTF參數估計基本思想的基礎上，受其實驗過程啟發，研究了人為附加高斯卷積模糊退化操作后的結果圖像和與原圖像所對應系統OTF參數的數學關系，提出了一種適應于實際應用要求、自動、準確、穩定的估計高斯型PSF參數的改進方法，并將其應用于尿沉渣顯微圖像復原，取得了較為良好的效果。本文所提出地改進方法對于具有高斯型近似PSF模型的其它類型圖像的復原和三維顯微圖像復原研究也具有應用和參考價值。
1 PSF參數的自動估計方法
針對高斯型PSF的表達式(3) ，令σ=D0，文獻[5]考察了不同D
1.1 特征圓增強的基本原理
假設為空域高斯型PSF表達模型，n(x,y)為服從N(0,σn2)的白噪聲，即。令，若對式(1) 所表達退化模型進行一次高斯模糊過程，根據卷積性質有[6]：
| $\begin{align} & b\left( x,y \right)=g\left( x,y \right)*{{G}_{b}}\left( x,y \right)=f\left( x,y \right)*h\left( x,y \right)+n\left( x,y \right)* \\ & {{G}_{b}}\left( x,y \right)=f\left( x,y \right)*\left[ h\left( x,y \right)*{{G}_{b}}\left( x,y \right) \right]+n\left( x,y \right)*{{G}_{b}}\left( x,y \right) \\ \end{align}$ |
實際上，式(4) 中h(x,y)*Gb(x,y)和n(x,y)*Gb(x,y)本質上都是兩個高斯函數的卷積，根據高斯函數的性質其結果仍然是高斯函數，即有式(5) ：
| $\begin{align} & h\left( x,y \right)*{{G}_{b}}\left( x,y \right)=\frac{1}{2\pi \left( \sigma _{b}^{2}+\sigma _{h}^{2} \right)}{{e}^{-\frac{{{x}^{2}}+{{y}^{2}}}{2\left( \sigma _{b}^{2}+\sigma _{h}^{2} \right)}}} \\ & n\left( x,y \right)*{{G}_{b}}\left( x,y \right)=\frac{1}{2\pi \left( \sigma _{b}^{2}+\sigma _{h}^{2} \right)}{{e}^{-\frac{{{x}^{2}}+{{y}^{2}}}{2\left( \sigma _{b}^{2}+\sigma _{n}^{2} \right)}}} \\ \end{align}$ |
經過高斯模糊后的圖像可以看成是被具有不同參數的OTF和噪聲模型退化的結果，其中新的OTF和噪聲空域參數分別為。
將式(5) 兩邊取傅里葉變換得OTF頻域表達為，其中。
| ${{D}_{b}}=\frac{1}{\sqrt{\sigma _{b}^{2}+\sigma _{h}^{2}}}$ |
當時，，實際上這正是文獻[5]中在實驗部分所描述的情況。
對于具有高斯型PSF的圖像復原，當其PSF所對應的標準差參數比較小，導致特征圓不明顯時，我們可以通過人為附加一具有已知較大合適參數σb的高斯卷積模糊運算，使得新圖像的傅里葉變換圖像中與OTF相關的特征圓特征變得極為明顯。圖 2所示為附加標準差為σb=2的人工模糊步驟前后，同一顯微圖像傅里葉變換幅度譜的對比，可以看出在附加人工模糊退化步驟后，特征圓特征變得相對明顯。
圖2
顯微圖像人為退化前后的幅度譜變化比較
(a) 原始顯微圖像幅度譜;(b) 人為退化后圖像幅度譜
Figure2. Comparison of the amplitude spectrum of a microscopic image between pre- and post-artificial degradation(a) the amplitude spectrum of an original microscopic image; (b) the amplitude spectrum after degraded
1.2 特征圓半徑的自動檢測
對于高斯模糊退化后進行傅里葉變換所得圖像的特征圓，通過手工方法測量其半徑后，即可計算得到新的OTF對應參數Db[5]，并由式(6) 得到原高斯型OTF參數。需要注意的是，由于傅里葉變換圖像乃頻域采樣圖像，文獻 [5]所述高斯降質和OTF參數實驗都是在頻域比較和度量的，因此在實驗過程中無需考慮頻率分辨率而只取其整數坐標度量即可。但本文涉及頻域和空域參數轉換，需要換算為真實角頻率。假設原始圖像大小為M×N，則其離散傅里葉變換圖像通常也為M×N，但其中橫向和縱向頻率分辨率分別為2π/M和2π/N，則在頻域中度量和空域真實度量之間轉換時需分別乘上相應方向的分辨率數值。
考慮到通過手工方式測量特征圓半徑時很容易引入較大的測量誤差。為了更準確地測量特征圓半徑，盡可能減少隨意性，本文給出了一種自動、穩定的精確測量方法，其基本流程如圖 3所示。
圖3
特征圓半徑自動測量計算流程
Figure3.
Flow chart of automatic radius measurement and computation of the feature circle
由于附加人為模糊后所得圖像的傅里葉變換中，特征圓內外區域對比度相對比較明顯，非常適合基于類間方差最大化原則實現二值化[7]。因此本文采用了大津算法實現圖像閾值化處理，圖 2(b)頻譜圖像的閾值化結果如圖 4(a)所示。可以看出，特征圓區域在二值化圖像中以高亮的類圓形體現，但沿橫縱兩頻率軸可看到一些細長條白色區域干擾。由于這些區域在投影圖中體現為類圓形區域位置兩側的低投影值區域，可以簡單地通過位置投影濾波方式予以排除。由于頻譜圖像二值化結果中的特征圓圓周是由很多隨機白點組成，為了穩定、有效求取特征圓半徑，特提出如下算法：
圖4
顯微圖像幅度譜中特征圓的自動檢測
(a) 幅度譜大津法二值化結果示例；(b) 特征圓檢測結果示例
Figure4. Automatic detection of the feature circle in amplitude spectrum of a microscopic image(a) example of thresholding the amplitude spectrum by Ostu method; (b) example of the result of feature circle detection
假設離散傅里葉變換圖像大小為M×N，則
(1) 求取類圓形區域的質心C(XC,YC),計算
(2) 若Error＞2，轉(7) ，否則轉(3) ；
(3) 以作為特征圓的圓心并設置半徑Rk=1，計算圖像中所有白點數目為NW；
(4) 計算二值圖像中以O為圓心Rk為半徑的圓內白點數目為NWk；
(5) 若，迭代終止；否則，Rk=Rk+1，轉(4) ；
(6) Rk即所求特征圓半徑。
(7) 報錯：特征圓不明顯，不適合用作PSF參數計算。
其最終檢測結果如圖 4(b)所示。
2 實驗與結論
為了驗證本文方法的有效性，基于尿沉渣顯微圖像，進行了相關參數測量和恢復實驗。對一原始512×512尿沉渣圖像，首先分別對其以一系列人為設置的標準差進行高斯模糊退化，然后自動測量特征圓半徑值，并進一步計算出各種人為退化和模糊標準差下對應的原始PSF參數估計值，結果如表 1所示。
表1
不同人為退化和模糊標準差下特征圓半徑和原始PSF參數估計值比較
Table1.
Comparison of feature circle radius and estimated original PSF parameter under different artificial degrading and blurring Standard Deviation
可以看出，人為附加的高斯標準差參數在0.75~2.0之間時，其PSF估計值的均值和標準偏差為1.938 55±0.027 52。可見在此區間內PSF參數估計值相對比較穩定，可以用來計算原始圖像PSF參數的估計值。用其計算所得平均值對圖像進行維納濾波復原[8]所得結果如圖 5(b)所示。可以明顯看出復原圖像對比度更好，細節更加清晰突出，在一定程度上證明了PSF參數估計的正確性。
圖5
復原前后尿沉渣顯微圖像比較
(a) 原始尿沉渣圖像；(b) 基于所估計PSF參數維納濾波復原結果
Figure5. Comparison of a urinary sediments microscopic image before and after restoration(a) original urinary sediments image; (b) the Wiener restoration result based on estimated PSF parameter
本文在文獻[5]基本原理和實驗啟發基礎上，采用附加模糊退化方法對頻譜圖像特征圓予以增強，并通過特征圓半徑自動計算實現了顯微圖像高斯型PSF參數的自動、準確和穩定估計。實驗中對于尿沉渣圖像的復原效果在一定程度上表明了本文所提出改進方法的合理性。由于受實驗條件局限，未能給出本文算法復原結果與更精密光鏡或電鏡等所成圖像的對比，今后將結合相關臨床應用進一步予以驗證。需要特別指出的是，由于高斯型OSF近似模型具有廣泛的適用性，本文OSF參數估計的改進方法不但可以適合顯微圖像的復原，還可以用于其它具有高斯型OSF近似模型的成像系統中(如天文望遠鏡系統和常規視頻攝像系統等)，實現相關圖像的復原。本文的基本思想還可應用于三維顯微圖像復原，實現三維OSF的參數估計。
