結(jié)構(gòu)光照明顯微的結(jié)構(gòu)光空間頻率和相位測定算法

張尹馨，鄧嘉俊，費(fèi)建陽

(1. 光電信息教育部重點(diǎn)實(shí)驗(yàn)室(天津大學(xué))，天津300072；2. 天津大學(xué)精密儀器與光電子工程學(xué)院，天津300072)

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

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

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

1 SIM的成像原理

在SIM 中，通常使用余弦形式的結(jié)構(gòu)光照明樣本，余弦結(jié)構(gòu)光可以表示為

式中：r 為空間坐標(biāo)；pd、φ、I0、m 分別為結(jié)構(gòu)光的空間頻率、初始相位、平均強(qiáng)度和調(diào)制度；下標(biāo)d 表示結(jié)構(gòu)光條紋的方向．

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

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

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

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

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

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

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

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

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

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

ISC 算法的原理為

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

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

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

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

對于AC 與POP 算法，當(dāng)結(jié)構(gòu)光調(diào)制度較低時，無法很好地滿足近似條件，此時計(jì)算結(jié)果存在誤差．

3 簡化頻譜分量互相關(guān)算法

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

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

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

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

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

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

在測定出結(jié)構(gòu)光的空間頻率和相位之后，還可以利用重疊區(qū)域內(nèi)的頻譜信息來測定結(jié)構(gòu)光的調(diào)制度m[20]，其計(jì)算式為

4 仿 真

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

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

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

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

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

為了獲得較高的分辨率，結(jié)構(gòu)光的空間頻率設(shè)定約為96.7%的OTF 截止頻率．有學(xué)者證明，當(dāng)結(jié)構(gòu)光的空間頻率很高時，POP 算法會失效[17]．因此，文中沒有討論P(yáng)OP 算法的相位誤差．

圖3 在不同的調(diào)制度和噪聲水平下，簡化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

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

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

圖4 在不同的調(diào)制度和噪聲水平下，簡化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 的值越大，說明重建圖像的質(zhì)量越好．首先分析結(jié)構(gòu)光調(diào)制度對重建圖像質(zhì)量的影響．當(dāng)結(jié)構(gòu)光調(diào)制度為0.80 時，無論使用哪種算法，其重建圖像的PSNR 均最大．所以隨著結(jié)構(gòu)光調(diào)制度的增加，重建圖像的質(zhì)量也逐漸提升．對于不同的算法而言，在不同的結(jié)構(gòu)光調(diào)制度下，簡化SCC 算法重建的圖像PSNR 比其他算法重建的圖像PSNR 更大，這說明簡化SCC 算法重建的圖像的質(zhì)量更好．所以由表1 可知，無論結(jié)構(gòu)光調(diào)制度高或低，使用簡化的SCC 算法所計(jì)算的結(jié)構(gòu)光空間頻率和相位所重建的圖像質(zhì)量均優(yōu)于以其他算法為基礎(chǔ)的重建圖像．

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

5 實(shí) 驗(yàn)

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

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

圖6 實(shí)驗(yàn)圖像重建結(jié)果比較Fig.6 Comparison of the experimentally reconstructed images

6 結(jié) 語

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