稀疏角度數據下非視域物體三維重建方法研究

王雪峰,張自豪,陳興穌,王元慶

(1.伊犁師范大學網絡安全與信息技術學院,新疆 伊寧 835000；2.河南工業大學人工智能與大數據學院,河南 鄭州 450001;3.南京大學電子科學與工程學院,江蘇 南京 210093)

1 引 言

非視域成像技術是近年來的研究熱點技術,是剛發展起來的一種間接成像技術,在醫療成像、災難救援、軍事反恐及城市街道安全偵查等都具有重要的研究意義。非視域成像技術通過激光照射中介面,主要是通過光在中介面進行多次反射和散射來獲取隱藏物體的信息,再使用一定的重建算法對采集的數據進行處理,才能完成非視域物體的三維重建。美國麻省理工學院的Velten等人[1-2]利用飛秒激光和條紋相機進行拐角成像,首先使用反投影重建法,實現了非視域物體的三維重建。為了提高重建的精度,對反投影法的優化算法被相繼提出,如使用時間選通門技術及單光子相機來對隱藏物體進行逐層探測[3-5],使用快速反投影算法來提高重建的速度[6],利用橢圓投影特點進行橢圓模型消除重建的偽影[7],文獻[8]提出了兩個權重因子應用在濾波反投影算法中,減少了條紋效應,并能有效應用到噪聲數據中。采用f-k偏移方法[9]為鏡面反射對象產生更高質量的重建；O′Toole等人[10-11]使用光錐變換及共聚焦光路對隱藏物體進行探測；Xin等人[12]使用費馬流算法得到了較高的重建精度。

由于反投影法需要采集足夠的數據(全部角度范圍),才能較好的完成三維重建,所以需要光源進行掃描式數據采集,需要大量時間,但是在實際的應用中,特別是對于較危險的物體,如反恐及災害救援,需要快速、有效的進行數據采集,并能夠快速完成三維重建,這就對非視域成像技術及重建算法提出了新的挑戰。本文針對稀疏角度數據問題,提出了反投影修正MLEM(最大似然期望值最大化)算法,能夠有效提高非視域物體三維重建的精度。

2 基于APD陣列的非視域成像

現有大多數成像方式為反射式成像,也稱為拐角成像,本文采用透射式成像方式,成像設置的示意圖如圖1所示,激光照射中介面一點(如圖中D點),光透射中介面后,中介面選擇擴散膜(不能透過擴散膜直接觀測到后面的物體),繼續散射到擴散膜后面的隱藏物體,光會再次被隱藏物體進行散射和反射,并返回到擴散膜后被探測器APD接收,其中圖1中參考光是檢測光源的出光時間,作為探測器APD接收信號的時間基準。本文提出的基于APD陣列的透射式非視域激光成像系統具有成像結構簡單、操作方便且無需場景掃描、能夠快速完成數據采集的特點。探測器布置在擴散膜附近,直接探測隱藏目標散射回來的光信息。一次激光照射后,APD陣列接收來自隱藏目標多個方向的信息,這相當于進行了多點信息接收,因此不需要進行場景掃描,使數據采集一次完成。

圖1 本文提出的透射式非視域激光成像示意圖

探測器陣列的排布方式采用矩形陣列方式,如圖2所示,以3×3的陣列為例,不同的探測器陣列(如圓環形狀和十字形狀等)采集到的數據會有不同,獲取的隱藏目標的信息也會有所不同,對于重建算法有一定的影響。本文中使用的探測器單元由4個APD單元組成,詳細的APD單元排布方式如圖1所示。4個APD排布在擴散膜的四個對角位置,排布成2行2列的矩陣陣列的方式,四個APD分別標記為APD11、APD12、APD21和APD22。矩形陣列探測的數據表示為:

本文提出的透射式非視域成像的光傳輸過程如圖3所示。

圖2 矩形陣列的探測器排布方式

圖3(a)中,首先激光照射擴散膜D點,光經過擴散膜后,繼續傳輸到隱藏物體上(如點S1、S2和S3),傳輸路徑為D→S1、D→S2和D→S3,傳輸距離為rd1、rd2和rd3；其次光經過隱藏物體后再次傳輸到擴散膜,并被探測器APD接收,如其中一個APD探測單元在A點,則傳輸路徑為S1→A、S2→A和S3→A,傳輸距離為rs1、rs2和rs3；最后在A點的APD探測器單元接收信號的路徑為D→S1→A、D→S2→A和D→S3→A,傳輸距離為rd1+rs1、rd2+rs2和rd3+rs3,則接收信號的三個時間點(ts1、ts2和ts3)如圖3(b)上面所示,對應圖3(a)中隱藏物體的三個點(S1、S2和S3),則ts1=(rd1+rs1)/c、ts2=(rd2+rs2)/c和ts3=(rd3+rs3)/c,其中c為光速。

圖4 非視域成像的橢圓投影原理圖

圖4中兩個定點是光源D點和探測器A點,隱藏物體上的點即是某個橢圓上的點,如橢圓f1上的動點p3；隱藏物體上的兩點p1和p2在同一個橢圓f2上,則p1和p2會在同一個時間點返回信息到探測器上,即光經過D點后傳輸到隱藏物體上兩點p1和p2,并投影到探測器上,這就是非視域成像中的橢圓投影方式。后期的圖像重建也是根據反投影方法的原理,首先將隱藏物體劃分為三維像素塊(稱為體素,如圖4中p1、p2和p3),然后根據橢圓反投影方法對隱藏物體的每個體素進行重建,最終完成整個物體的三維重建。

3 基于BP-MLEM的非視域物體重建算法

3.1 MLEM算法

MLEM算法是最大似然期望值最大化算法(maximum likelihood expectation maximization),它使用期望最大化EM(expectation maximization)算法[13-14]來完成最大似然的估計,是一種統計迭代算法,以“參數估計理論”為基礎,基于觀測數據統計模型的重建算法,廣泛應用于CT及PETCT圖像重建[15-16],CT圖像重建的主要思想也是反投影[17],所以也需要完整的數據才能完成較好的重建,但實際中無法完成180°范圍的數據采集,且需要采集的時間間隔要小、均勻才能更好的重建圖像,針對CT圖像重建中稀疏角度數據的問題,可以使用迭代重建算法。

MLEM算法是根據最大似然估計理論,假設認為每次橢圓投影符合泊松分布,根據似然函數,求投影數據的最大概率,推導得出計算公式,而此公式正好為非線性方程,很難給出解析解,只能通過迭代的方法來求近似解,通常采用期望最大化(EM)來求解。由此,推導得出了MLEM算法,迭代公式[18]如下:

(1)

其中,重建體素設為xj;重建體素個數為M個;投影矩陣中橢圓個數為i個;投影矩陣元素aij表示第j個體素xj對第i個橢圓的貢獻值;pi表示第i個橢圓的投影值。

基于APD陣列的非視域成像MLEM算法的具體步驟如下:

(4)計算第j個重建體素xj的修正值:

(2)

(5)對重建體素xj的值進行修正:

(3)

這里用穿過該體素的所有橢圓對它進行修正,即完成一輪迭代。

(6)k=k+1；將以上迭代的結果作為初值,重復步驟(2)到(5)進行新一輪迭代,直到取得符合收斂條件為止。

3.2 BP-MLEM算法

針對非視域成像的特點,僅能在稀疏角度下獲取數據,提出了BP-MLEM(Back Projection -MLEM)迭代重建算法,更加適用于非視域成像,改進的方法對偽影和噪聲的抑制較好,且具有更快的收斂速度。

MLEM算法中比較重要的迭代過程是修正方法,它決定了每次迭代結果后的修正方向,為了更快更準確的得到收斂結果,需要迭代方向每次都往正確的結果修正,這就需要更好的修正方法,MLEM算法修正方法的修正值為:

采用上次迭代結果和修正值進行相乘的方法進行迭代,即下一次迭代值是上一次迭代值與反投影乘積的結果。這樣迭代后圖像中較大的值增長越大,較小的值增長較慢,圖像的對比度越來越大,即迭代速度快,但有時速度過快,容易造成過收斂,不能得到理想的收斂結果。

本文提出的基于反投影迭代BP-MLEM算法,迭代公式如下:

(4)

其中,重建體素設為xj;重建體素個數為M個;投影矩陣中橢圓個數為i個;投影矩陣元素aij表示第j個體素xj對第i個橢圓的貢獻值;pi表示第i個橢圓的投影值。

基于APD陣列的非視域成像BP-MLEM算法的具體步驟如下:

(4)計算第j個重建體素xj的修正值為:

(5)

(5)對重建體素xj值進行修正,方法為:

(6)

這里用穿過該體素的所有橢圓對重建體素xj進行修正,并將反投影BPΔi加入到修正方法中,這樣就完成一輪迭代。

(6)k=k+1；將以上迭代的結果作為初值,重復步驟(2)到(5)進行新一輪迭代,直到取得符合收斂條件為止。

4 實驗結果

實驗設置如圖1所示,隱藏目標選擇“十”字形狀,距離探測器220 cm,4個APD探測器的位置如圖1所示。分別使用反投影算法、MLEM算法和BP-MLEM算法進行實驗,為了更好的進行對比實驗,本文進行了多次迭代重建,分別給出不同迭代次數的對比結果。

圖5顯示了使用傳統的反投影算法得到重建結果,反投影算法中使用了R-S-L濾波算子,從圖中可以看出,反投影的重建結果存在較大的偽影,重建結果精度不高。

圖5 濾波反投影算法重建結果圖

圖6顯示的是使用MLEM算法對非視域物體“十”字的重建結果,分別顯示了4次不同迭代次數(分別是第3、4、9、20次迭代),從圖中可以看出,當第3次迭代(圖6(a))時,重建結果圖的偽影在減少,但是重建結果圖并不理想,圖像中心點強度很高,而邊緣強度很低(接近背景色),不能較好的重建“十”字形狀的結果；再次進行6次迭代后,即第9次迭代(圖6(c)),重建結果的偽影比第3次迭代少,但是重建的結果仍然是中心點強度過大,而邊緣強度太小；當達到第20次迭代時,重建結果反而沒有最初的迭代結果好,幾乎無法重建“十”形狀結果,這就是MLEM算法迭代速度過快,從而產生過收斂的結果。

本文提出的BP-MLEM算法的重建結果如圖7所示,同時也給出了4個不同迭代次數(分別是第3、4、9、20次迭代),總體來看,4次迭代的重建結果都能夠重建出“十”字形狀圖,并且隨著迭代次數的增加,重建偽影在不斷減少,重建的精度也在不斷增加；當迭代次數達到4次時(圖7(b)),物體“十”字形狀就得到了較好的重建結果圖；當迭代次數不斷增加時,重建結果沒有出現MLEM算法的過收斂結果,而是達到了更好的重建效果,重建結果圖像趨于穩定。實驗表明,BP-MLEM算法能夠在較少的迭代次數下,得到較好的重建結果；并且當迭代次數增加時,重建算法能夠得到較穩定的重建結果圖,不會出現過收斂等問題。

為了更加準確的對比實驗結果,本文使用了SSIM(Structural SIMilarity)結構相似度來評價重建結果圖的圖像質量,對比實驗結果顯示在圖8中,分別對使用BP-MLEM算法的重建結果圖(圖7)和MLEM算法的重建結果圖(圖6)計算SSIM值,圖8中分別計算了第3、4、9、20次迭代的SSIM值,其中第一個值(橫坐標標識為BP)是使用反投影算法重建結果圖(圖5)。從圖8中可以看出,利用BP-MLEM算法重建結果圖的SSIM值都高于MLEM算法重結果圖的SSIM值,而且隨著迭代次數的增加,BP-MLEM算法的重建結果圖的SSIM也在不斷增加,表明重建結果圖的質量在不斷增加；而MLEM算法的重建結果圖的SSIM值雖然剛開始也在不斷增加,但增加的趨勢較慢,且當迭代次數為第9次和20次時,SSIM值基本沒有增加,表明重建結果圖的質量并沒有隨著迭代次數的增加而不斷增加。

圖8 BP-MLEM算法與MLEM算法重建結果圖在 不同迭代次數下的SSIM值對比圖

5 結 語

圖像重建算法在非視域成像系統中具有重要的作用,它關乎著整個系統的實現結果,對實驗系統的驗證起到了關鍵作用。本文構建了基于APD陣列的透射式成像系統,快速進行隱藏物體探測,一次完成數據采集,不需要使用光源進行掃描成像；并針對非視域成像系統僅能在稀疏角度下獲取不完整數據的特點,提出了BP-MLEM迭代算法,實驗結果表明,相比反投影算法,能夠更好的去除偽影；并且能在較少的迭代次數下,得到較好的重建結果圖；當重建次數不斷增加時,重建結果趨于穩定,表明BP-MLEM算法較穩定,不會出現過收斂等問題。