基于分數階螺旋相位片的定量相位顯微成像*

吳迪 蔣子珍 喻歡歡 張晨爽 張嬌 林丹櫻 于斌 屈軍樂

(深圳大學物理與光電工程學院, 光電子器件與系統教育部/廣東省重點實驗室, 深圳 518060)

1 引 言

定量相位成像(quantitative phase imaging,QPI)技術將相位成像和光學顯微相結合, 為觀測透明生物樣品提供了一種無標記、快速、無損、高分辨率的成像手段, 廣泛應用于生物醫學等領域[1].目前, QPI技術主要分為基于干涉法測量和非干涉法測量兩大類, 包括離軸法[2,3]、相移方法[4,5]、共光路法[6,7]、 白光法[8,9]及傅里葉變換光散射法[10]等多種方法.

徑向希爾伯特變換是實現圖像邊緣增強中最常用的方法之一.在4f濾波系統的傅里葉譜平面上使用螺旋相位濾波器(spiral phase filter, SPF)能夠實現徑向希爾伯特變換[11], 從而使振幅和相位樣本產生各向同性的邊緣增強, 可用于觀察透明生物組織和細胞, 因此, SPF倍受研究人員關注[12-14].在SPF的基礎上, Furhapter等[15]發展了螺旋相襯(spiral phase contrast, SPC)顯微成像技術, 通過利用拓撲荷數為1的整數階螺旋相位片進行相位濾波, 再通過三步相移法進行復雜物體相位的定量重構, 但其圖像采集和處理過程相對復雜.當螺旋相位片的拓撲荷數為分數時, SPF系統點擴散函數的徑向對稱被打破, 由此產生了可控的邊緣位錯, 可以通過控制分數階拓撲荷數的大小和初始相位改變其邊緣增強的程度和方向, 實現物體的各向異性邊緣增強[16-18].目前, 基于分數階螺旋相位片的SPF研究主要圍繞圖像的邊緣增強展開,但在QPI方面研究相對較少.

基于迭代的相位恢復方法, 例如Gerchberg-Saxton (GS)[19]、輸入-輸出算法[20]等, 可以從一幅(或多幅)強度測量中定量恢復物體的相位分布, 實現了非干涉法的QPI, 引起研究人員的廣泛關注.然而, 直接利用樣品強度圖像恢復相位是一個病態問題, 存在解的唯一性和迭代耗時的問題, 這往往需要通過對強度信息進行過采樣來解決這一問題.

基于此, 將相位恢復與螺旋相襯顯微相結合,本文提出了一種基于分數階螺旋相位片的定量相位成像方法和系統, 通過一幅經分數階SPF的樣品強度圖像, 利用基于螺旋相襯的GS相位恢復算法(SPC-based GS phase retrieval algorithm,SGSA)實現了樣品相位的定量重構, 簡化了螺旋相襯顯微的實驗過程和相位重構步驟.計算機模擬實驗研究了基于不同拓撲荷數的螺旋相位濾波成像和相位重構過程, 分析了其可行性.最后, 搭建了基于空間光調制器(spatial light modulator,SLM)的分數階螺旋相襯顯微成像系統, 通過對相位型光柵和生物細胞樣品進行了成像和相位重構,驗證了基于分數階螺旋相位片的相襯顯微方法可以有效地提高螺旋相襯顯微成像的對比度, 能夠定量獲其樣品的相位信息, 對于定量相襯顯微術的發展具有重要的研究意義和應用價值.

2 螺旋相襯顯微成像原理

2.1 成像原理

螺旋相襯顯微成像是基于一個4f系統, 其光學系統示意圖如圖1所示.

圖1 基于4f系統的螺旋相襯成像系統示意圖Fig.1.Schematic diagram of 4f system-based spiral phase contrast imaging system.

平行入射光照射樣品后的光場為 g (r,θ) , 經過4f系統的第1個透鏡L1的傅里葉變換, 在其后焦平面得到物光場的頻譜 f (ρ,φ)=FT[g(r,θ)].然后,通過螺旋相位片 H (ρ,φ)=exp(ilφ) 對其進行濾波,l為拓撲荷數, φ為繞螺旋相位片中心的方位角.最后, 螺旋相位片濾波后的光場通過4f系統的第2個透鏡的逆傅里葉變換后, 在其后焦平面, 即探測平面, 得到輸出光場.利用探測器記錄輸出強度圖像 G =|F(r,θ)|2, 其成像過程表示如下:

其中 h (r,θ) 代表螺旋相位片的傅里葉變換, 相當于4f系統的點擴散函數.

當利用整數階螺旋相位片進行相位濾波成像時, 由于整數階螺旋相位片的點擴散函數的中心奇點的振幅為0, 因此, 輸入信息的低頻部分可以被有效地抑制, 并且主瓣各個方向上的振幅分布比較均勻, 所以, 能夠對相位樣品產生各向同性邊緣增強.

由于分數階螺旋相位片可以表示為有限個不同整數階螺旋相位片的線性疊加, 因此, 其透過率函數可以表示為傅里葉級數的形式[21]:

其中ρ為徑向坐標, φ為方位角坐標, n為整數,l取值一般為小于1的分數.

分數階螺旋相位片作為相位濾波器所對應的點擴散函數可以表示為

其中m為最接近l的整數, J|m|為第|m|階貝塞爾函數.

對于拓撲荷數l ＜ 0.5的分數階螺旋相位片的點擴散函數中并沒有渦旋產生, 因此, 并不能對相位物體實現高襯度的邊緣增強.然而, 由于其非對稱性, 仍可以獲得一定的強度梯度圖像[18].利用拓撲荷數l ＜ 0.5分數階螺旋相位片濾波得到的強度圖, 再經過相位恢復算法處理后, 能夠得到比整數階SPF濾波后的強度圖經SGSA處理后對比度更高的相位圖, 更有利于實現定量相位顯微成像.

2.2 定量相位重構

將相位恢復算法與SPC相結合, 利用探測到的螺旋相位濾波圖像, 實現樣品相位的定量重建.在傳統GS算法的基礎上, 在迭代過程中, 在頻譜面引入螺旋相位片的限制條件, 然后, 進行優化.SGSA的算法框圖如圖2所示.

圖2 SGSA框圖Fig.2.Block diagram of the SGSA.

以第n階迭代為例, 具體步驟如下.

1)第n次迭代的目標光場復振幅分布gn(r,θ)經過傅里葉變換獲得其頻譜分布 fn(ρ,φ) :

2)在頻譜面引入SPF對 fn(ρ,φ) 進行限制, 獲得濾波后的頻譜

3)濾波后的頻譜, 經過逆傅里葉變換, 獲得像平面的復振幅分布 Gn(r,θ) :

4)引入像面上的限制條件, 即保持相位不變,但振幅變為預先探測到的像面的振幅分布, 獲得新的復振幅分布

6)在頻譜面引入共軛的SPF對其進行濾波,獲得物體的頻譜

7)將獲得的物體的頻譜, 通過逆傅里葉變換,獲得物體的復振幅分布 g′(r,θ) :

8)再引入物面上的限制條件, 即保持相位不變, 振幅變為已知的物面的振幅分布, 作為下一次迭代的物波函數:

9)反復進行1)至8)的迭代運算, 直到定義的誤差—均方差之和(sum square error, SSE)到達預先設置的精度或達到設置的最大迭代次數為止.SSE的定義如下:

其中 ρ2表示輸出平面探測器接收到的像振幅分布,表示第n次迭代結束后輸出平面像的振幅分布.下面實驗中, 設置迭代終止的 S SE＜10-3.

3 計算機仿真實驗

基于Matlab, 利用快速傅里葉變換算法, 模擬螺旋相襯成像及相位重構過程.以圖像“LENNA”作為相位物體, 像元數為1024 × 1024, 物理尺寸為332.8 μm × 332.8 μm, 波長 λ =640 nm, 相位值在[0 1]弧度之間變化[22], 滿足弱相位成像條件,對于細胞的成像來說, 這是合理的.對于螺旋相位片, 整數階的拓撲荷數取 l =1 , 分數階拓撲荷數取l=0.1.利用方程(1)獲得了經SPF后的模擬成像圖像.考慮到相機的噪聲, 對生成的圖像添加了高斯噪聲(均值為0, 方差為0.01), 以更符合實際情況.然后, 再對所得強度圖利用SGSA算法重構獲得樣品的相位信息.由于SGSA是一個迭代算法,為了進一步說明迭代次數的選擇, 模擬了拓撲荷取0.1時的分數階螺旋相位濾波成像, 計算了SSE和迭代次數的關系, 如圖3所示.另外, 考慮到數值模擬時相位是已知的, 也計算了恢復相位與真實相位之間的均方誤差(mean square error,MSE)與迭代次數的關系, 如圖4所示.

圖3 分數階拓撲荷取0.1時, SSE隨迭代次數的變化曲線Fig.3.SSE error vs.the number of iterations when the fractional topological charge is 0.1.

圖4 恢復相位與真實相位之間的MSE隨迭代次數的變化曲線Fig.4.MSE error between the recovered phase image and ground truth phase image vs.the number of iterations.

MSE定義為

綜合考慮MSE和SSE, 考慮到計算時間, 選擇了迭代次數設置為20次.樣品原圖像, 整數階螺旋相位片濾波與分數階螺旋相位片濾波后的成像強度圖對比及重構出的樣品相位圖對比結果如圖5所示.

圖5 螺旋相位濾波成像及恢復結果對比 (a) 相位型樣品原圖; (b) 傳統整數階螺旋相位片濾波圖像; (c) 對圖(b)用SGSA恢復的樣品相位圖; (d) 分數階螺旋相位片濾波圖像; (e) 對圖(d)用SGSA恢復的樣品相位圖Fig.5.Comparisons of the recorded images and the recovered results: (a) The ground truth phase sample image; (b) the recorded image via traditional integer order spiral phase plate filtering; (c) the recovered phase sample image using SGSA for panel (b);(d) the recorded image via fractional spiral phase plate filtering; (e) the recovered phase sample image using SGSA for panel (d).

為了確定分數階螺旋相位片的拓撲荷數, 模擬分析了拓撲荷數分別取1.0, 0.8, 0.6, 0.5, 0.4, 0.2,0.1和0.08時的分數階螺旋相位濾波成像, 成像結果如圖6所示.從圖6模擬結果顯示, 隨著拓撲荷數取值逐漸減小, SPF圖像的邊緣增強效果逐漸減弱, 但圖像中低頻信息的對比度逐漸提高, 與文獻[16]中的結論一致.在此基礎上, 繼續減小l的取值, 使之小于0.1, 成像對比度與l = 0.1時相比,低頻信息對比度提升不大, 主要原因來自于相位片展開式中的高階項可以忽略不計, 對物的頻譜濾波不起作用.

圖6 分數拓撲荷l取不同值時的直接相位成像結果 ( a)l=1 ; ( b)l=0.8 ; ( c)l=0.6 ; ( d)l=0.5 ; ( e)l=0.4 ; ( f)l=0.2 ;(g)l=0.1 ;(h)l=0.08Fig.6.Direct phase imaging results for different fractional topological charge l :(a)l=1;(b)l=0.8;(c)l=0.6;(d)l=0.5;(e)l=0.4 ; ( f)l=0.2 ; ( g)l=0.1 ; ( h)l=0.08.

根據圖6的模擬成像結果, 選擇拓撲荷數較小的分數階螺旋相位濾波器獲得的相位濾波圖像, 有利于研究人員直接觀察濾波圖像, 獲取樣品的狀態, 選擇樣品感興趣區域進行成像, 然后再利用相位恢復算法做進一步重構處理.

為了說明分數階螺旋相位濾波器的特性, 首先, 比較了不同分數階拓撲荷數下的相位濾波圖像;其次, 在相位濾波強度圖對比度接近的情況下, 再繼續比較分析了拓撲荷數取0.1與拓撲荷小于0.1的恢復相位, 如圖7所示, 并以測量圖像質量的結構相似度(structural similarity, SSIM)指數來定量對比:

圖7 不同拓撲荷值恢復結果對比 (a) 拓撲荷取0.1時經SGSA恢復出的相位圖; (b) 拓撲荷取0.08時經SGSA恢復出的相位圖Fig.7.Comparisons of the recovered results for different topologies: (a) The reovered phase image using SGSA when the topology is 0.1; (b) the revoverd phase image using SGSA when the topology is 0.08.

其中 μX,μY,σX,σY,σXY分別表示圖像X和Y的局部均值、標準差和互補方差.

由圖7可知, 拓撲荷數取值小于0.1時的恢復相位分布圖的對比度比拓撲荷等于0.1時恢復相位分布圖的對比度差, 且表1中的SSIM值也定量地證明了這一結論.因此, 綜合考慮不同拓撲荷下分數階螺旋相位濾波直接成像的強度分布圖對比度與不同拓撲荷下恢復相位分布圖的對比度兩個因素, 分數階螺旋相位片的拓撲荷取0.1較為合適.因此, 在實際實驗中, 選擇拓撲荷數為0.1的螺旋相位片來進行實際樣品成像.

表1 三種分數拓撲荷下恢復相位圖與原圖的SSIM值Table 1.The SSIM between the recovered phase image and the ground truth phase image for three fractional topological values.

4 實驗結果與討論

4.1 螺旋相位片

為了便于實際應用, 在實驗中, 利用SLM來產生螺旋相位片, 為了避免零級光的影響, 利用相位型的閃耀光柵與螺旋相位片相結合的方式形成新的螺旋相位圖, 如圖8所示.

圖8 (a) 相位型閃耀光柵; (b) 螺旋相位片; (c) 復合螺旋相位圖Fig.8.(a) Phase-type blazed grating; (b) spiral phase plate; (c) composite spiral phase plate.

4.2 實驗系統光路

基于SLM的螺旋相襯定量相位成像系統的實驗光路如圖9所示.

圖9 (a) 基于螺旋相位片濾波的定量相位成像系統光路圖; (b) SLM上加載的整數階叉形光柵; (c) SLM上加載的分數階叉形光柵Fig.9.(a) Optical setup of quantitative phase imaging system based on a spiral phase filter; (b) the integer order fork grating loaded on SLM; (c) the fractional fork grating loaded on SLM.

波長為640 nm的激光器(Coherent, OBIS)發出的光束經衰減片, 針孔濾波器(物鏡, 40/0.65 160/0.17)及準直透鏡(焦距100 mm)產生平行光束, 照射透明相位樣品.再經樣品調制的光束通過成像物鏡(Nikon, Plan Fluor, 20×/0.50), 管鏡(焦距為200 mm)成像在一次像面處, 然后經過4f螺旋相位濾波系統成像在探測器CMOS (Hamamatsu, C11440-22 CU, 像元尺寸 6.5 μm × 6.5 μm)上.其中, 4f透鏡焦距是500 mm, 頻譜面上放置SLM (Holoeye, NIR, 像元尺寸是8 μm × 8 μm),用于螺旋相位片的顯示.

4.3 實驗結果

首先, 利用定做的相位光柵( S iO2玻璃, 折射率1.456, 光柵刻線理論設計深度200 nm, 光柵刻線實際測試深度為150 nm, 光柵周期為6 μm)進行拓撲荷數 l =1 的螺旋相位片濾波的成像實驗,目的是與分數階螺旋相位濾波成像結果進行對比,實驗結果如圖10所示.

圖10 定制相位型光柵的成像 (a) 相位型光柵未濾波明場強度圖; (b) 整數階螺旋相位片濾波成像邊緣增強圖; (c) 恢復相位圖; (d) 恢復深度圖, 橫縱坐標數值為像素值, 每個像素代表0.325 μm, 總長度為332.8 μmFig.10.Imaging of a custom phase gratinig: (a) The unfiltered bright field image of the phase grating; (b) the recorded integer-order spiral phase filtered edge enhancement image; (c) the recovered phase image; (d) the recovered depth image, the abscissa and ordinate values are pixel values, each pixel represents 0.325 μm, and the total length is 332.8 μm.

從圖10可以看出, 整數階螺旋相位濾波的圖像邊緣增強效果顯著, 再對采集到的強度圖利用SGSA進行相位恢復, 可進一步提高圖像的對比度.由SSE隨迭代次數變化曲線圖11可知, 只需兩次迭代即可使得SGSA趨于收斂, 得到樣品相位圖, 耗時1.412 s (Windows10, 內存16 GB, CPU:Intel (R) Core(TM) i5-9400F CPU @2.90 GHz;Matlab R2018b).

圖11 經整數階螺旋相位片濾波后再由SGSA重構的相位型光柵SSE隨迭代次數的變化Fig.11.SSE error vs.the number of iterations for the grating phase reconstruction problem of phase retrieval from a integer-order spiral phase filtering intensity measurement using the SGSA.

再用MATLAB對恢復相位圖進一步分析可知, 得到恢復相位光柵頂部相位值為0.24π, 光柵底部相位值為—0.06π, 平均相位差為0.30π, 結合公式可得, 恢復相位圖的定量深度h平均為216 nm, 與已知光柵刻線深度150 nm的高度值誤差為66 nm 左右.分析其誤差產生的原因, 主要來自于用SLM的像素化所產生的螺旋相位片精度不高; 另外, 加工的光柵刻線深度也有一定的誤差.為了進一步提高測量精度, 可以加工相關高精度相位片.

其次, 對拓撲荷數 l =0.1 分數階螺旋相位濾波進行實驗, 結果如圖12所示.

圖12 定制相位型光柵的成像 (a) 相位型光柵未濾波明場強度圖; (b) 拓撲荷l取0.1時分數階螺旋相位片濾波成像強度圖;(c) 恢復相位圖; (d) 恢復深度圖(坐標同圖10(d))Fig.12.Imaging of a custom phase grating: (a) The unfiltered phase grating bright field image; (b) the fractional spiral phase plate filtered image when the topological charge l is 0.1; (c) the recovered phase image; (d) the recovered depth image (the coordinates are the same as Fig.10.(d)).

從圖13可以看出, 當使用l = 0.1的分數階螺旋相位片直接對相位型樣品進行螺旋相位濾波成像時, 成像圖的對比度相對于未經濾波得到的成像圖對比度明顯得到提升, 效果與模擬結果得出的結論一致.因此, 為了方便觀察樣品, 可以直接利用分數階螺旋相位片進行螺旋相位濾波成像.然后對單次拍攝采集到的單幅螺旋相位濾波強度圖, 利用SGSA進行相位恢復進一步重構樣品的相位信息.實驗結果證明, 分數階螺旋相位濾波可以得到保留更多低頻信息的高對比度強度分布圖, 有利于對相位物體進行實時顯微觀測.從圖13的誤差曲線可知, 需迭代10次, 耗時2.839 s, 誤差幾乎減為0.進一步分析可知, 恢復相位光柵頂部平均相位值為0.25π, 光柵底部平均相位值為—0.03π, 均值相位差為0.28π, 則恢復相位圖的定量深度h為198 nm, 與已知光柵實際刻蝕深度150 nm的高度值誤差為48 nm.與整數階SPC相比, 測量誤差減小27%.

圖13 經分數階螺旋相位片濾波后再由SGSA重構的相位光柵SSE隨迭代次數的變化Fig.13.SSE error vs.the number of iterations for the grating phase reconstruction problem of phase retrieval from a fractional spiral phase filtering intensity measurement using the SGSA.

為了進一步驗證分數階SPC成像效果, 對SHSY5Y人神經母細胞瘤細胞進行成像, 如圖14所示.

圖14 SH-SY5Y細胞成像 (a) SH-SY5Y細胞未濾波明場強度圖; (b) 拓撲荷l取0.1時分數階螺旋相位片濾波神經元細胞成像強度圖; (c) 恢復相位圖; (d) 定量相移圖(坐標同圖10(d))Fig.14.SH-SY5Y cell imaging: (a) The unfiltered SH-SY5Y cell bright field image; (b) the intensity image of the neuron cell using the fractional spiral phase filter when the topological charge l is 0.1; (c) the recoverd phase image; (d) the quantitative phase image(the coordinates are the same as Fig.10.(d)).

從圖15誤差隨迭代次數變化曲線可知, 只需迭代5次即可, 耗時1.782 s.由于并不清楚細胞樣品內部的具體折射率分布, 因此, 只給出細胞的定量相移圖.恢復相位圖的最大相位值為 0.6141π , 最小相位值為 - 0.3234π.從圖15可以看出, 當使用分數階螺旋相位片對生物細胞進行成像時, 成像效果明顯, 對比度高, 再經過SGSA恢復后得到的樣品的相位信息, 與模擬結果相一致.

圖15 SH-SY5Y細胞經分數階螺旋相位濾波后再由SGSA重建的SSE隨迭代次數的變化Fig.15.SSE error vs.the number of iterations for the SHSY5Y cells reconstruction problem of phase retrieval from a fractional spiral phase filtering intensity measurement using the SGSA.

5 結 論

本文實現了一種基于空間光調制器的分數階螺旋相襯顯微成像系統, 基于單幅螺旋相位濾波圖像, 通過SGSA實現了定量相位顯微成像, 簡化了螺旋相襯顯微的圖像采集和相位重構過程.計算機仿真和實驗結果表明, 基于分數階拓撲荷數的螺旋相位成像, 圖像低頻信息得到保留, 直接成像對比度高, 更加接近樣品真實相位信息, 有利于在實驗過程中及時觀察樣品的形態和選取感興趣區域進行相位濾波成像, 為進一步開展其在活體細胞快速定量相位成像中的應用奠定了技術基礎.