本發(fā)明屬于生物醫(yī)學(xué)圖像重建技術(shù)領(lǐng)域，具體涉及一種基于TV(Total Variation，全變差)約束的斷層相位顯微鏡重建方法。

背景技術(shù)：

許多活體細(xì)胞和生物樣品呈現(xiàn)一定的透明性，而其折射率本身則通常用于區(qū)分細(xì)胞或生物樣品結(jié)構(gòu)特性的內(nèi)在的對(duì)比源。近年來(lái)，隨著光學(xué)顯微鏡以更高的時(shí)空分辨率不斷發(fā)展，它已成為生物學(xué)和醫(yī)學(xué)研究中不可或缺的工具。其中斷層相位顯微鏡(tomographic phase microscopy，TPM)由于其直接對(duì)生物樣品的折射率分布進(jìn)行成像，而不需要任何添加任何染料，越來(lái)越成為生物學(xué)家的研究熱點(diǎn)。

斷層相位顯微鏡通常采用外差式激光干涉儀，定量獲取樣品在不同角度的光的照射下形成的干涉條紋，并解析出相應(yīng)的相位圖像，隨后基于各角度相位圖像，重建出樣品三維折射率的分布。因而圖像的重建算法在提高折射率圖像的質(zhì)量和時(shí)空分辨率中，有著重要的作用。通常，基于對(duì)激光經(jīng)過(guò)樣品后的衍射場(chǎng)的不同解釋?zhuān)糜陲@微鏡圖像重建算法可以分為兩類(lèi)。第一類(lèi)方法假定激光是直線穿過(guò)樣品，并且所得的相位，是光線穿過(guò)樣品的線積分，這類(lèi)一般利用濾波反投影或者改進(jìn)的濾波反投影算法，對(duì)折射率圖像進(jìn)行重建。該類(lèi)方法直觀，而且在低散射的樣品中表現(xiàn)出很好的效果，但是對(duì)于大樣品或者折射率高的樣品，重建結(jié)果卻不太理想。另一類(lèi)方法考慮到光學(xué)場(chǎng)的衍射效應(yīng)，在頻域空間建立樣品折射率和衍射場(chǎng)的關(guān)系模型，并采用Born近似或Rytov近似，將二者關(guān)系線性化，得到所謂的傅里葉衍射定理，從而簡(jiǎn)化關(guān)系模型；然后，采用三維傅里葉你變換，重建出樣品的折射率分布。

由于以上兩類(lèi)方法都只是在傅里葉空間進(jìn)行線性濾波，忽略了整個(gè)探測(cè)過(guò)程的影響。而且，這種傅里葉理論并不能充分表征圖像數(shù)據(jù)的空間各項(xiàng)異性特征。同時(shí)，許多生物應(yīng)用，都要求對(duì)生物的動(dòng)態(tài)過(guò)程進(jìn)行成像，因而顯微鏡的成像速度，即時(shí)間分辨率，也是需要考慮的重要因素。然而，在通常的斷層相位顯微鏡中需要依次對(duì)樣品進(jìn)行多個(gè)角度的采樣，成像的空間分辨率受制于采集角度的數(shù)目和相鄰采集角度間隔的大小。為了提高成像空間分辨率，就需要對(duì)樣品進(jìn)行密集的、更多角度的采樣，這就需要更多的成像時(shí)間，犧牲了成像的空間分辨率。

雖然Fang-Yen C等在2011年發(fā)表的文章Video-rate tomographic phase microscopy中通過(guò)改進(jìn)儀器設(shè)置減少每一角度的采樣時(shí)間，進(jìn)而是成像的速度達(dá)到視頻級(jí)別。但是這種方法仍然像以往一樣，需要采集很多角度的數(shù)據(jù)。因而，我們認(rèn)為，如果可以在不降低圖像空間分辨率的情況下，進(jìn)一步減小采樣頻率，即減少總得角度采樣次數(shù)，那么我們將進(jìn)一步提高顯微鏡的成像速度，但是隨著采集頻率變小，將很難從采集數(shù)據(jù)從重建出高質(zhì)量的折射率圖像。另一方面，相位顯微鏡中存在的缺失錐的問(wèn)題，即由于顯微鏡有限的數(shù)值孔徑導(dǎo)致光束只能局部掃描部分角度，而另一部分角度數(shù)據(jù)則無(wú)法采集。這種數(shù)據(jù)的缺失，將致使重建后的折射率圖像在軸向拉伸、折射率數(shù)值偏低。由于采集數(shù)據(jù)的不足，利用普通的解析算法難以重建樣品的折射率圖像，這對(duì)圖像重建提出了新的挑戰(zhàn)。

技術(shù)實(shí)現(xiàn)要素：

針對(duì)現(xiàn)有技術(shù)所存在的上述技術(shù)問(wèn)題，本發(fā)明提出了一種基于TV約束的斷層相位顯微鏡重建方法，通過(guò)建立統(tǒng)一的數(shù)據(jù)采集系統(tǒng)模型，然后加入TV約束條件，并利用加速交替方向乘子法快速求解，重建出折射率分布圖像。

一種基于TV約束的斷層相位顯微鏡重建方法，包括如下步驟：

(1)采集放樣品和不放樣品兩種情況下對(duì)應(yīng)多個(gè)角度的激光干涉圖像，并從中解析出純樣品對(duì)應(yīng)各角度的相位圖像；

(2)建立樣品關(guān)于折射率與相位之間的系統(tǒng)模型如下：

f＝Au+e

其中：A為系統(tǒng)矩陣，u為平行于激光旋轉(zhuǎn)方向任一斷層的樣品折射率圖像，f對(duì)應(yīng)為該斷層的多角度樣品相位圖像，e為測(cè)量噪聲向量；

(3)基于上述系統(tǒng)模型加入TV約束，得到顯微鏡圖像重建的目標(biāo)函數(shù)如下，根據(jù)該目標(biāo)函數(shù)遍歷各斷層進(jìn)行優(yōu)化求解，重建得到樣品的三維折射率圖像；

$\underset{u}{\min }\left\{\frac{\mathrm{\μ}}{2}\right||Au-f|{|}_2+\left|\right|u\left|{|}_{tv}\right\}$

其中：||||₂表示2范數(shù)，||||_tv表示TV范數(shù)，μ為權(quán)重系數(shù)。

所述的步驟(1)中利用馬赫-曾德干涉儀采集對(duì)應(yīng)多個(gè)角度的激光干涉圖像。

所述的步驟(1)中解析純樣品對(duì)應(yīng)各角度的相位圖像：對(duì)于任一角度，采用頻域分析法從放樣品情況下對(duì)應(yīng)該角度的激光干涉圖像中提取出包含有樣品和背景的相位信息，同時(shí)從不放樣品情況下對(duì)應(yīng)該角度的激光干涉圖像中提取出只包含純背景的相位信息；進(jìn)而使著兩部分相位信息相減即得到純樣品對(duì)應(yīng)該角度的相位圖像。

所述的多角度樣品相位圖像f由對(duì)應(yīng)該斷層的K組相位信息組成，其中第k組相位信息即為該斷層在純樣品第k個(gè)角度的相位圖像中所對(duì)應(yīng)的一行相位信息，k為自然數(shù)且1≤k≤K，K為角度的總數(shù)量。

所述的系統(tǒng)矩陣A為m×n維矩陣，m為多角度樣品相位圖像f的總像素個(gè)數(shù)，n為樣品折射率圖像u的總像素個(gè)數(shù)，系統(tǒng)矩陣A中的第i行第j列元素值為樣品折射率圖像u中第i個(gè)像素被多角度樣品相位圖像f中第j個(gè)像素所對(duì)應(yīng)光束覆蓋的面積，i和j均為自然數(shù)且1≤i≤n，1≤j≤m。

所述樣品折射率圖像u的TV范數(shù)||u||_tv表達(dá)式如下：

$\left|\right|u|{|}_{tv}=\underset{i=1}{\overset{n}{\Σ}}|\left|D\left(u_i\right)\right||$

其中：D(u_i)為二維向量且該二維向量中的兩個(gè)元素值分別為u_{i_下}-u_i和u_{i_右}-u_i,u_i為樣品折射率圖像u中第i個(gè)像素的折射率，u_{i_下}為樣品折射率圖像u中與第i個(gè)像素相鄰的下邊一個(gè)像素的折射率，u_{i_右}為樣品折射率圖像u中與第i個(gè)像素相鄰的右邊一個(gè)像素的折射率。

所述的步驟(3)中根據(jù)目標(biāo)函數(shù)采用加速交替方向乘子法遍歷各斷層進(jìn)行優(yōu)化求解，重建得到樣品的三維折射率圖像。

本發(fā)明通過(guò)建立顯微鏡成像的系統(tǒng)模型，將數(shù)據(jù)采集過(guò)程中的各種不確定因素納入統(tǒng)一的描述框架，方便耦合成像過(guò)程中的過(guò)程誤差和物理模型誤差；基于圖像的稀疏特性，加上TV的保邊去噪的特性，有效解決了由于缺失錐帶來(lái)的數(shù)據(jù)不完整的問(wèn)題，緩解了重建后圖像軸向拉伸效應(yīng)，提高了顯微鏡的軸向分辨率；同時(shí)也可以通過(guò)稀疏角度采集相位數(shù)據(jù)，減少了總的采集角的數(shù)量，從而減少成像時(shí)間，提高顯微鏡的時(shí)間分辨率，利于觀察生物樣品的動(dòng)態(tài)特性。

附圖說(shuō)明

圖1為本發(fā)明中相位顯微鏡系統(tǒng)的結(jié)構(gòu)示意圖。

圖2為本發(fā)明重建方法的步驟流程示意圖。

圖3(a)為采用約束型FBP方法在1的頻率下采樣數(shù)據(jù)所重建得到微珠折射率圖像(軸向切片)。

圖3(b)為采用約束型FBP方法在1/5的頻率下采樣數(shù)據(jù)所重建得到微珠折射率圖像(軸向切片)。

圖3(c)為采用約束型FBP方法在1/10的頻率下采樣數(shù)據(jù)所重建得到微珠折射率圖像(軸向切片)。

圖3(d)為采用本發(fā)明方法在1的頻率下采樣數(shù)據(jù)所重建得到微珠折射率圖像(軸向切片)。

圖3(e)為采用本發(fā)明方法在1/5的頻率下采樣數(shù)據(jù)所重建得到微珠折射率圖像(軸向切片)。

圖3(f)為采用本發(fā)明方法在1/10的頻率下采樣數(shù)據(jù)所重建得到微珠折射率圖像(軸向切片)。

圖4(a)為采用約束型FBP方法在1的頻率下采樣數(shù)據(jù)所重建得到HeLa細(xì)胞折射率圖像(軸向切片)。

圖4(b)為采用約束型FBP方法在1/5的頻率下采樣數(shù)據(jù)所重建得到HeLa細(xì)胞折射率圖像(軸向切片)。

圖4(c)為采用約束型FBP方法在1/10的頻率下采樣數(shù)據(jù)所重建得到HeLa細(xì)胞折射率圖像(軸向切片)。

圖4(d)為采用本發(fā)明方法在1的頻率下采樣數(shù)據(jù)所重建得到HeLa細(xì)胞折射率圖像(軸向切片)。

圖4(e)為采用本發(fā)明方法在1/5的頻率下采樣數(shù)據(jù)所重建得到HeLa細(xì)胞折射率圖像(軸向切片)。

圖4(f)為采用本發(fā)明方法在1/10的頻率下采樣數(shù)據(jù)所重建得到HeLa細(xì)胞折射率圖像(軸向切片)。

具體實(shí)施方式

為了更為具體地描述本發(fā)明，下面結(jié)合附圖及具體實(shí)施方式對(duì)本發(fā)明的技術(shù)方案進(jìn)行詳細(xì)說(shuō)明。

如圖1所示的斷層相位顯微鏡系統(tǒng)中，He-Ne激光器發(fā)射的激光經(jīng)由分束器分為樣品光束和參考光束兩束光。在樣品光束處，由掃描陣鏡(GM)控制光束的方向，使該光束以設(shè)定的角步旋轉(zhuǎn)，然后光束由透鏡(L1，f＝250mm)聚焦在油浸物鏡(CL，Nikon，1.4NA)的后焦面上。掃描陣鏡到透鏡L1的距離和透鏡L1到油浸物鏡后焦面(BF)的距離相等，以保證掃描陣鏡和樣品(Sample)共軛。樣品放置在載樣容器中，光束由成像物鏡(OL，Nikon，1.4NA)采集，然后經(jīng)由管鏡(L2，f＝200mm)成像到CCD相機(jī)感光片上。在參考光處，光束由聲光調(diào)制器進(jìn)行調(diào)制，產(chǎn)生10Hz的頻移。然后由分光片(BS)將它與樣品光合并，在相機(jī)感光片上產(chǎn)生干涉條紋圖像。掃描陣鏡以0.67°的角步旋轉(zhuǎn)樣品光束，掃描范圍為在樣品處從-67°到67°。由于邊緣角度的成像質(zhì)量差，在實(shí)際重建過(guò)程中，只用到-60°到60°的采集數(shù)據(jù)。這樣總共有180個(gè)采集角度。在每一角步處，我們連續(xù)采集四幀干涉圖像?；谶@四幀干涉圖像，采用相移干涉法計(jì)算出光束穿過(guò)樣品后產(chǎn)生的相位圖像。然后，在不放置樣品的情況下，同樣采集一系列的干涉圖像，并計(jì)算出相應(yīng)的相位圖像，作為背景相位圖像。最后，由放置樣品時(shí)各個(gè)角度的相位圖像減去相應(yīng)的不放置樣品時(shí)背景相位圖像，得到純粹由樣品本身折射率引起的相位圖像，并以此作為重建原始數(shù)據(jù)。

本發(fā)明重建方法的流程如圖2所示，首先建立相位圖像和樣品折射率分布圖像的系統(tǒng)矩陣模型：

f＝Au+e

其中：A為系統(tǒng)矩陣；u為平行于激光旋轉(zhuǎn)方向的任意一層樣品折射率圖像，且為n維向量；f為對(duì)應(yīng)于該層u的多角度的相位圖像，且為m維向量；e為測(cè)量噪聲向量；m和n均為大于1的自然數(shù)。

然后加入TV約束，得到圖像目標(biāo)函數(shù)為：

$min\left\{\frac{\mathrm{\μ}}{2}\right||Au-f|{|}_2+\left|\right|u\left|{|}_{tv}\right\}$

其中：||||₂為二范數(shù)，||||_tv為T(mén)V范數(shù)，μ為權(quán)重系數(shù)，用于平衡上述公式第一項(xiàng)和第二項(xiàng)的權(quán)重。本實(shí)施方式中，μ設(shè)為2⁸。

最后采用加速交替方向乘子法對(duì)目標(biāo)函數(shù)求解，得到以下迭代步驟：

q_k＝(1-α_k-1)q_k-1+α_k-1u_k

p_k＝(1-α_k)q_k+α_ku_k

$u_{k+1}=\arg \min \<\▿G\left(p_k\right),u\>+\frac{{\mathrm{\θ}}_k}{2}\left|\right|w_k-Du|{|}^2+\<{\mathrm{\λ}}_k,D\>+\frac{{\mathrm{\η}}_k}{2}||u-u_k|{|}^2$

$w_{k+1}=\arg \min F\left(w\right)-\<{\mathrm{\λ}}_k,w\>+\frac{{\mathrm{\θ}}_k}{2}\left|\right|w-{\mathrm{Du}}_{k+1}|{|}^2$

v_k+1＝(1-α_k)v_k+α_kw_k+1

λ_k+1＝λ_k-ρ_k(w_k+1-Du_k+1)

其中：arg min表示求解函數(shù)以使該項(xiàng)取得最小值，為梯度算子，＜,＞表示內(nèi)積，D為2×n為T(mén)V差分矩陣，N為圖像像素的總數(shù)量，k為迭代次數(shù)，L_G為L(zhǎng)ipschitz常數(shù)。初始化數(shù)值q₀＝u₀為隨機(jī)初始化生成。

迭代收斂條件如下：

$\frac{\underset{j=1}{\overset{N}{\Σ}}|u_k^j-u_{k-1}^j|}{N}\≤\mathrm{\ρ}$

其中：為第k次迭代后的重建圖像向量的第j元素值，為第k-1次迭代后的重建圖像向量的第j元素值，ρ為給定的收斂閾值，j為自然數(shù)且1≤j≤N。本實(shí)施方式中ρ設(shè)為10^-4。

以下我們采用polystyrene微珠和真實(shí)HeLa細(xì)胞進(jìn)行試驗(yàn)，以驗(yàn)證本實(shí)施方式的有效性。我們?cè)谑褂帽景l(fā)明方法進(jìn)行重建的同時(shí)，也采用約束型濾波反投影重建算法進(jìn)行重建，將兩者的重建結(jié)果進(jìn)行對(duì)比。同時(shí)，也對(duì)比了稀疏角度數(shù)據(jù)的重建結(jié)果，具體操作是：分別以1、1/5、1/10的頻率從180個(gè)角度數(shù)據(jù)抽取數(shù)據(jù)，即分別采用180個(gè)角度數(shù)據(jù)、36個(gè)角度數(shù)據(jù)、18個(gè)角度數(shù)據(jù)進(jìn)行重建。

對(duì)于微珠實(shí)驗(yàn)，我們采用直徑為4.5μm、折射率為1.588的微珠，放置在折射率為1.518的油中進(jìn)行成像。我們采集到180個(gè)角度的相位圖像，每個(gè)圖像的大小為196×196。隨后對(duì)這些數(shù)據(jù)重組成196層的2D切片sinogram圖像。然后依次對(duì)各層進(jìn)行重建，得到196層196×196的2D圖像，并合并成196×196×196的3D圖像，其重建結(jié)果如圖3所示，其中各圖像均為中間一層的軸向截面。圖3(a)～圖3(c)分別為約束型FBP方法由1、1/5、1/10的頻率采樣數(shù)據(jù)重建得到的結(jié)果。圖3(d)～圖3(f)分別為本發(fā)明重建方法由1、1/5、1/10的頻率采樣數(shù)據(jù)重建得到的結(jié)果。結(jié)果顯示，約束型FBP方法重建的結(jié)果，折射率圖像在軸向有明顯的拉伸，而本發(fā)明方法的重建結(jié)果很好地保持了微珠的球形形狀，減小了缺失錐帶來(lái)的影響。同時(shí)，隨著采樣頻率的減少，約束型FBP方法的圖像質(zhì)量明顯下降，而本發(fā)明方法的重建即使在采樣頻率減少到1/10時(shí)(相當(dāng)于，只需要原來(lái)1/10采集時(shí)間)，圖像質(zhì)量依舊很好。

對(duì)于HeLa細(xì)胞實(shí)驗(yàn)，我們將HeLa細(xì)胞浸在折射率為1.335的細(xì)胞培養(yǎng)液磷酸鹽緩沖鹽水中。同樣，我們采集到180個(gè)角度的相位圖像，每個(gè)圖像的大小為256×256。隨后對(duì)這些數(shù)據(jù)重組成256層的2D切片sinogram圖像。然后依次對(duì)各層進(jìn)行重建，得到256層256×256的2D圖像，并合并成256×256×256的3D圖像，其重建結(jié)果如圖4所示，其中各圖像均為中間一層的橫向截面。圖4(a)～圖4(c)分別為約束型FBP方法由1、1/5、1/10的頻率采樣數(shù)據(jù)重建得到的結(jié)果。圖4(d)～圖4(f)分別為本發(fā)明重建方法由1、1/5、1/10的頻率采樣數(shù)據(jù)重建得到的結(jié)果。結(jié)果同樣顯示，隨著采樣頻率的減少，約束型FBP方法重建的圖像質(zhì)量逐漸下降，并慢慢淹沒(méi)在噪聲之中，而本發(fā)明方法的卻很好的保持了圖像的質(zhì)量。

本實(shí)施方式對(duì)斷層相位顯微鏡圖像進(jìn)行重建，與傳統(tǒng)的約束型FBP方法相比，本實(shí)施方式既有效緩解了缺失錐問(wèn)題帶來(lái)的圖像軸向拉伸的問(wèn)題，又可以在稀疏采樣數(shù)據(jù)的情況下重建出高質(zhì)量的圖像，這也意味著，只需要采集少量的角度數(shù)據(jù)就可以恢復(fù)折射率分布圖像，大大提高了顯微鏡的成像速度。

上述對(duì)實(shí)施例的描述是為便于本技術(shù)領(lǐng)域的普通技術(shù)人員能理解和應(yīng)用本發(fā)明。熟悉本領(lǐng)域技術(shù)的人員顯然可以容易地對(duì)上述實(shí)施例做出各種修改，并把在此說(shuō)明的一般原理應(yīng)用到其他實(shí)施例中而不必經(jīng)過(guò)創(chuàng)造性的勞動(dòng)。因此，本發(fā)明不限于上述實(shí)施例，本領(lǐng)域技術(shù)人員根據(jù)本發(fā)明的揭示，對(duì)于本發(fā)明做出的改進(jìn)和修改都應(yīng)該在本發(fā)明的保護(hù)范圍之內(nèi)。