張尹馨，鄧嘉俊，費建陽

(1. 光電信息教育部重點實驗室(天津大學)，天津300072；2. 天津大學精密儀器與光電子工程學院，天津300072)

在生命科學研究領域，超分辨顯微成像技術蓬勃發(fā)展，并且獲得了超過衍射極限的成像分辨率[1-3].在這些技術中，結構光照明顯微(structured illumination microscopy，SIM)因為其成像速度高、光毒性低和對熒光特性要求少等優(yōu)點而得到廣泛應用[4-7]．在SIM 中，使用周期性的結構光照明樣品后，可以將原本位于光學傳遞函數(shù)(optical transfer function，OTF)支持域外的高頻信息移動到OTF 支持域內，此時被移動的高頻信息與低頻信息混合形成混頻信息．圖像重建算法可分離低頻和高頻信息，并將高頻信息移位，使之與低頻信息合并形成擴展頻譜．最后對拓展頻譜進行傅里葉逆變換得到重建的超分辨圖像[8-10]．

準確測定結構光的相位和空間頻率對于分離和移動高頻信息至關重要，也是重建超分辨率圖像的關鍵[11-13]．有關空間頻率和相位測定算法已有諸多報道，如迭代自相關算法(iterative auto-correlation algorithm，IAC)、迭代正弦相關算法(iterative sinusoidal correlation algorithm，ISC)[14]、峰值相位算法(phase of the peaks algorithm，POP)[15]、自相關算法(autocorrelation algorithm，AC)[16]和頻譜分量互相關算法(spectral-components cross-correlation algorithm，SCC)等[17-18]．在結構光調制度較高時，這些算法均能準確計算結構光的空間頻率和相位，然而并非所有SIM系統(tǒng)都能獲得理想的結構光調制度．投影式SIM 系統(tǒng)[19]因其結構簡單，圖像采集速度快，相比結構復雜的干涉式SIM 系統(tǒng)[20]應用更為廣泛．但是投影式 SIM 系統(tǒng)所生成的結構光調制度會受到物鏡OTF 的影響，且隨著結構光空間頻率和成像深度的增加而減?。送?，當前SIM 系統(tǒng)逐步小型化，從實驗室走向便攜式應用也成為可能，光學結構由于震動產生移位會使得結構光調制度降低．因此，解決低調制度下的結構光空間頻率和相位的準確測算問題不可或缺，算法不僅要保證估算精度，還要具備廣泛適應性．

本文研究了前述幾種結構光空間頻率和相位估算算法在低調制度下的性能，分析發(fā)現(xiàn)SCC 算法能有效保證參數(shù)的估算精度，在不同調制度下具有較為廣泛的適應性．之后進一步研究了SCC 算法，對其進行了簡化．原始的SCC 算法對同一個重疊區(qū)域的相關值進行兩次迭代來分別測定空間頻率和相位，本文簡化后的SCC 算法能夠同時測定空間頻率和相位．通過模擬仿真驗證簡化SCC 算法在低調制度下的性能，相比于其他算法，該算法可以更精確地測定空間頻率和相位．實驗上，利用自建SIM 系統(tǒng)對牛肺動脈內皮細胞(bovine pulmonary artery endothelial cells，BPAE)成像，使用簡化SCC 算法測定了結構光的空間頻率和相位，以此為基礎重建的圖像更為清晰，有效抑制了重建圖像的偽影．

1 SIM的成像原理

在SIM 中，通常使用余弦形式的結構光照明樣本，余弦結構光可以表示為

式中：r 為空間坐標；pd、φ、I0、m 分別為結構光的空間頻率、初始相位、平均強度和調制度；下標d 表示結構光條紋的方向．

SIM 系統(tǒng)觀測到的樣本圖像可以表示為

樣本的熒光分子濃度 C ( r )與結構光Id,φ( r )的乘積，再卷積光學系統(tǒng)的點擴散函數(shù)(point spread function，PSF) H ( r )，即為SIM 原始圖像Ed,φ(r)．原理如圖1 所示，所得到的SIM 原始圖像如圖1(b)、(c)、(d)所示．

對式(2)執(zhí)行傅里葉變換，可以得到SIM 原始圖像的頻譜表達式為

在求解出低頻和高頻分量之后，需要對頻譜分量進行維納濾波[14]，即

式中：上標* 代表復共軛；Ψ0,d、Ψa,d和Ψb,d是的平均噪聲功率，通過和位于OTF 截止頻率外頻譜幅度的平方均值來估算[14]；α和β為常數(shù)，使用迭代非線性回歸的方法計算[15]．濾波后的頻譜分量如圖1(e)、(g)、(i)所示．

圖1 SIM原理示意Fig.1 Schematic diagram of structured illumination microscopy(SIM)

濾波后再對高頻分量進行移位，使用傅里葉移位定理將高頻分量移回到原來的空間位置[14]．如圖1(f)、(h)、(j)所示．移位之后，根據(jù)式(7)將所有的頻譜分量通過廣義維納濾波進行疊加，得到擴展的超分辨圖像頻譜，如圖1(k)所示．

式中：A ( k) 為切趾函數(shù)；w 為維納濾波常數(shù)．通過對執(zhí)行傅里葉逆變換可獲得重建的超分辨圖像．

2 空間頻率和相位測定算法原理

Lal 等[14]提出了IAC 和ISC 算法分別測定結構光的空間頻率和相位．IAC 算法原理為

ISC 算法的原理為

式中 I1(r)=－cos(2πpdr+φ0)，φ0是一個估計的初始相位．通過迭代變換相位 φ0，并根據(jù)式(10)計算 I1( r)與SIM 原始圖像Ed,φ(r)的相關值，直到迭代求出最大的相關值Cφ，此時的相位 φ0認為是準確值．但是，如果結構光調制度較低，很難將余弦函數(shù) I1( r )與Ed,φ(r)上的結構光條紋進行準確匹配，從而產生誤差.

Shroff 等[15]提出 POP 算法測定結構光的相位．根據(jù)式(3)，當k = pd時，式(3)變?yōu)?/p>

當SIM 原始圖像的信噪比較高，結構光的調制度較高，空間頻率適中且物頻譜的高頻信息衰減得足夠快時，式(11)可以近似為

Wicker[16]也提出AC 算法測定結構光的相位.AC 算法的原理與IAC 算法和POP 算法原理較為相似．首先根據(jù)式(9)計算出自相關值Ck，然后根據(jù)POP 算法的近似條件對 Ck進行近似處理，可以計算出結構光的相位為

對于AC 與POP 算法，當結構光調制度較低時，無法很好地滿足近似條件，此時計算結果存在誤差．

3 簡化頻譜分量互相關算法

式中c 為白噪聲的標準偏差．算法的思路如下：如果已知準確的空間頻率和相位，那么就可以準確地分離出低頻分量和高頻分量，并能將高頻分量移動到正確的空間位置上，由于重疊區(qū)域頻譜信息的相似性，其互相關值較高；如果空間頻率和相位不準確，那么低頻與高頻分量就無法準確地分離，高頻分量也無法準確地移動，此時重疊區(qū)域的互相關值較低．可以利用這種性質，迭代計算出空間頻率和相位．

Gustafsson 等[18]提出的SCC 算法通過迭代優(yōu)化重疊區(qū)域的相關值來計算空間頻率，而Wicker 等[17]使用的SCC 算法也是優(yōu)化同樣的重疊區(qū)域計算相位．本文的研究發(fā)現(xiàn)，兩者是通過優(yōu)化相同的重疊區(qū)域來計算不同的參數(shù)．因此，SCC 算法可以簡化合并，使其能夠同時計算空間頻率和相位．

簡化后SCC 算法流程如圖2 所示．首先，設定空間頻率的迭代初始值pd和相位初始值1φ 、2φ 、3φ ．然后，根據(jù)式(4)分別求出低頻分量和高頻分量，并利用傅里葉移位定理將高頻分量進行移位．接下來根據(jù)式(14)計算重疊區(qū)域的互相關值．之后將計算出的相關值取負，將相關值和迭代的初始點代入Matlab 函數(shù)fminsearch 中，對空間頻率和相位進行迭代優(yōu)化[21]，當－為極小值時，所對應的空間頻率和相位就是準確值．該過程可用式(16)表示．函數(shù)fminsearch 的搜索原理是基于Nelder-Mead 的單純形算法，其具體原理可參考文獻[21]．

圖2 簡化SCC算法流程Fig.2 Flow chart of simplified spectral-component crosscorrelation(SCC) algorithm

圖1(k)中包括兩個重疊的區(qū)域，但是這兩個重疊的區(qū)域所包含的信息在數(shù)學上是相同的，所以只需要計算任意一個重疊區(qū)域的互相關值即可．關于迭代初始點的選擇，可以將實驗中直接測定的空間頻率和相位值作為迭代的初始點．上述是0°方向結構光條紋的空間頻率和相位的計算方法，其他條紋方向的計算方法與之相同．

簡化的SCC 算法是一種迭代的算法，不存在近似的過程，所以理論上迭代出的空間頻率和相位值即為準確值．而且簡化SCC 算法是對分離后的頻譜分量的重疊區(qū)域計算相關值，相比于IAC 算法，其避免了其他頻譜分量對相關值產生影響，所以在低調制度下也能保持準確的結果．

在測定出結構光的空間頻率和相位之后，還可以利用重疊區(qū)域內的頻譜信息來測定結構光的調制度m[20]，其計算式為

4 仿 真

通過仿真對算法進行驗證．仿真條件如下：結構光的空間頻率為(165 nm)－1，樣本空間的像素尺寸設定為30 nm，激發(fā)波長488 nm，發(fā)射波長512 nm，物鏡的數(shù)值孔徑為1.49，結構光條紋的方向分別為0°、60°和120°，結構光的理想相位分別為0、2 π /3和4 π /3，并且給相位引入了一個隨機的誤差．

在不同的結構光調制度下，模擬不同的空間頻率和相位測定算法的計算結果．結構光調制度的范圍設定為0.01～1.00．為了分析噪聲的影響，給SIM 的原始圖像添加高斯噪聲和泊松噪聲．高斯噪聲的方差為0.001．對于泊松噪聲，使圖像中最亮的像素總計為1 000 個預期的光子，作為所有單個SIM 原始圖像的總和．在空間頻率的估算仿真上，將簡化的SCC算法與IAC 算法的結果進行比較，結果如圖3 所示．

由于結構光的空間頻率為矢量，通過分別計算空間頻率的長度誤差和方向誤差來表征空間頻率誤差．對于簡化的SCC 算法，在所討論的調制度范圍內空間頻率的長度誤差均低于0.065%，空間頻率的方向誤差低于0.02°．而對于IAC 算法，當調制度為0.18 時，空間頻率的長度誤差分別為0.64%(無噪聲)、1.37%(泊松噪聲)和1.82%(高斯噪聲)．空間頻率的相位誤差分別為0.024 5°(無噪聲)、0.072 7°(泊松噪聲)和0.226 3°(高斯噪聲)．當調制度低于0.18時，IAC 算法的空間頻率誤差還會繼續(xù)增大．當調制度大于0.18 時，這兩種算法的空間頻率誤差都很小．所以，簡化的SCC 算法在低調制度下有更準確的計算結果．

在相位的估算上，將簡化的SCC 算法、AC 算法和ISC 算法的相位誤差進行比較，仿真結果如圖4所示．簡化后的SCC 算法可以同時計算結構光的空間頻率和相位．但是AC 算法和ISC 算法需要提前知道空間頻率才能夠計算相位．在后面的討論中，AC算法和ISC 算法使用IAC 算法所計算出來的空間頻率去計算相位．

在圖4 中，簡化的SCC 算法在無噪聲且調制度低于0.03 時相位誤差低于5°，當調制度大于0.03時，相位誤差低于2.5°．如果存在泊松噪聲，該算法在調制度低于0.04 時相位誤差接近5°，當調制度大于0.04 時，相位誤差稍大于2.5°．假設存在高斯噪聲，在調制度低于0.07 時，該算法的相位誤差大于5°，當調制度大于0.07 時，其相位誤差低于5°，并隨著調制度增大而逐漸減?。畬τ贏C 算法和ISC 算法，當調制度低于0.18 時，由于IAC 算法存在誤差，因此AC 和ISC 算法也受其影響存在較大的誤差．當調制度大于0.18 時，雖然IAC 算法的空間頻率誤差已經(jīng)很小，但是受到噪聲和調制度的影響，AC 算法和ISC 算法的相位誤差仍然大于簡化SCC 算法．所以，在整個調制度范圍內，簡化SCC 算法的相位誤差更小．

為了獲得較高的分辨率，結構光的空間頻率設定約為96.7%的OTF 截止頻率．有學者證明，當結構光的空間頻率很高時，POP 算法會失效[17]．因此，文中沒有討論POP 算法的相位誤差．

圖3 在不同的調制度和噪聲水平下，簡化SCC算法與IAC算法的空間頻率長度與方向誤差Fig.3 Simulation of the magnitude and direction errors of the spatial frequencies of simplified spectral-component crosscorrelation(SCC)and IAC algorithms at different modulation depths and SNRs

在使用不同的算法計算出空間頻率和相位后，利用所計算出的空間頻率和相位值來重建超分辨圖像，結果如圖5 所示．圖5(a)是仿真時使用的樣本圖像，圖5(b)是反卷積寬場圖像，圖5(c)～(k)是使用了不同的空間頻率和相位測定算法并結合式(4)～(8)在不同的調制度下所重建的圖像．可以看到，當調制度為0.12 時，使用IAC、ISC 算法(圖5(f))和IAC、AC算法(圖5(i))重建的圖像存在很大的誤差，而簡化的SCC 算法可以成功地重建超分辨圖像(圖5(c))．當調制度為0.18 時，使用簡化SCC 算法重建的圖像(圖5(d))比使用IAC、ISC 算法和IAC、AC 算法(圖5(g)和圖5(j))重建的圖像質量更好．當調制度為0.80 時，上述算法均能獲得良好的重建效果．

為了定量地比較使用不同的算法在不同的調制度下所重建的圖像的質量，對重建圖像的峰值信噪比(peak-signal-to-noise ratio，PSNR)進行計算．計算結果如表1 所示．

圖4 在不同的調制度和噪聲水平下，簡化SCC 算法、AC算法和ISC算法的相位誤差Fig. 4 Phase errors of the simplified spectral-component cross-correlation(SCC)，AC，and ISC algorithms at different modulation depths and SNRs

表1 重建圖像的峰值信噪比Tab.1 PSNRs of the reconstructed images

PSNR 的值越大，說明重建圖像的質量越好．首先分析結構光調制度對重建圖像質量的影響．當結構光調制度為0.80 時，無論使用哪種算法，其重建圖像的PSNR 均最大．所以隨著結構光調制度的增加，重建圖像的質量也逐漸提升．對于不同的算法而言，在不同的結構光調制度下，簡化SCC 算法重建的圖像PSNR 比其他算法重建的圖像PSNR 更大，這說明簡化SCC 算法重建的圖像的質量更好．所以由表1 可知，無論結構光調制度高或低，使用簡化的SCC 算法所計算的結構光空間頻率和相位所重建的圖像質量均優(yōu)于以其他算法為基礎的重建圖像．

圖5 圖像重建仿真Fig.5 Simulations of reconstructed images

5 實 驗

在自建的 SIM 系統(tǒng)上對牛肺動脈內皮細胞(bovine pulmonary artery endothelial cells，BPAE)成像，檢驗簡化SCC 算法在超分辨圖像重建上的效果.SIM 系統(tǒng)所使用的物鏡是 Nikon CFI APOTRIF 100XH，數(shù)值孔徑1.49．顯微鏡主體的中間放大倍率為1.5 倍，總的放大倍率為150 倍．顯微鏡的探測器是靈敏的sCOMS(Andor Zyla 4.2 Plus)，像素尺寸為6.5 μm. 樣品是三色染色的 BPAE 細胞，用MitoTracker?Red CMXRos 標記線粒體，藍色熒光DNA 染料DAPI 標記細胞核，Alexa Fluor 488 鬼筆環(huán)肽標記F-肌動蛋白．激發(fā)波長488 nm，發(fā)射波長512 nm．結構光的調制度約為0.17．在重建圖像之前，需要對原始圖像進行預處理．首先，通過Matlab函數(shù) adapthisteq 來提高原始圖像的局部對比度．然后，選擇其中一張原始圖像作為參考基準，調整其他原始圖像的強度直方圖與參考圖像的直方圖匹配，可使用Matlab 函數(shù)imhistmatch 來進行[12]．

用不同算法計算出的空間頻率和相位所重建的圖像如圖6 所示．在結構光調制度為0.17 時，使用簡化SCC 算法所計算出的空間頻率和相位重建的圖像(圖6(b))明顯優(yōu)于以其他算法為前提的重建圖像(圖6(c)和(d))．與反卷積后的寬場圖像(圖6(a))相比，圖6(b)中重建圖像呈現(xiàn)的結構更為精細，如白色箭頭標注的位置能清晰展現(xiàn)3 個分支，而反卷積寬場圖像則無法分辨．此處的歸一化強度曲線如圖6(e)所示，重建圖像的分辨率明顯提升．實驗結果表明，簡化SCC 算法在低調制度下有助于更清晰地重建超分辨圖像．

圖6 實驗圖像重建結果比較Fig.6 Comparison of the experimentally reconstructed images

6 結 語

本文研究了結構光在低調制度下空間頻率和相位的測定算法，通過對多種算法的分析和對比發(fā)現(xiàn)SCC 算法能在不同調制度下具有更廣泛的適應性，精確地測算相關參數(shù)．文中對SCC 算法做了進一步研究，通過合并空間頻率和相位的迭代過程實現(xiàn)了SCC 算法的簡化，簡化后的SCC 算法可以同時精確地測定空間頻率和相位．仿真結果表明，該算法在結構光調制度較低時仍能準確地計算參數(shù)，相比于其他算法，其空間頻率和相位誤差更?。谧越ǖ腟IM 系統(tǒng)中對BPAE 細胞樣本成像，在低調制度下，用簡化的SCC 算法測定空間頻率和相位，更為清晰地重建了BPAE 細胞圖像，分辨率明顯提升，并抑制了重建圖像的偽影．簡化的SCC 算法能適應不同SIM 系統(tǒng)對結構光空間頻率和相位精確測定的需求，也能提高SIM 系統(tǒng)的魯棒性．