基于壓縮采樣匹配追蹤的生物發(fā)光斷層成像

蒲 鑫，侯榆青，朱 力，郭紅波

(西北大學(xué)信息科學(xué)與技術(shù)學(xué)院，陜西西安 710127)

隨著醫(yī)學(xué)影像技術(shù)的發(fā)展，功能成像已為大家所熟知，分子影像可以從細(xì)胞水平對疾病做出提前預(yù)警，并且能夠進(jìn)行在體的醫(yī)學(xué)實驗，實時反映生物體的病理過程，這些優(yōu)點使它成為了醫(yī)學(xué)界研究的熱點。生物發(fā)光斷層成像BLT(bioluminescence tomography)作為分子影像的一種成像方法，具有成本低、無輻射、操作簡單等優(yōu)點［1］。在BLT的研究中，光源重建是個典型的逆向問題，即通過動物體表含噪聲的表面測量光子密度，結(jié)合生物體解剖結(jié)構(gòu)信息和光學(xué)參數(shù)信息等，采用一定重建算法來確定其內(nèi)部的光源分布，定位病灶位置［2］。光在生物組織中傳播時的多重散射，導(dǎo)致了BLT逆問題的病態(tài)性［3－4］，并且需要重建的光源信息量要遠(yuǎn)大于邊界測量的光信息量，使得BLT成為典型的不適定問題。重建算法的核心就是求解這個不適定問題。為了解決這個問題，在文獻(xiàn)［5］中提出了一種光源可行區(qū)域方法，將其作為先驗信息，這種方法在很大程度上降低了其病態(tài)性。然而光源可行區(qū)域是根據(jù)光源表面分布和動物表面的異質(zhì)結(jié)構(gòu)來推斷的，當(dāng)更深層的多光源存在的時候，這種推斷并不一定正確。為此文獻(xiàn)［6］中提出了一種多光譜的方法，將熒光素酶的輻射分為不同波段，增加了外部測量信息，提高了測量效率，在一定程度上降低了其不適定性。然而多光譜的數(shù)據(jù)量大，增加了垃圾信息，給BLT重建的計算效率帶來了一定壓力?；诖?，學(xué)者們?nèi)诤狭丝尚袇^(qū)域、解剖信息以及多光譜等先驗信息以及測量噪聲，將重建問題轉(zhuǎn)換為一個最小二乘問題，利用正則化方法來求解。TIkhonov作為L2正則化的代表［7］首先提供一個相對平滑的解，但其重建效果卻存在偽影，降低了定位的準(zhǔn)確性。后來，L1正則化引起了研究者的注意，由于其可以從較少的邊界測量得到較好的重建圖像，克服了L2范數(shù)求解的過平滑性［8］。L1范數(shù)是一個凸優(yōu)化問題，可以采用現(xiàn)有的優(yōu)化手段求解，比如一些追蹤算法［9］。另一方面研究表明，L1正則化算法即使在測量數(shù)據(jù)非常有限的情況下也能有很好的光源重建［10］。關(guān)于L1正則化的問題在壓縮感知領(lǐng)域已經(jīng)有了廣泛的研究。

本文依據(jù)傳統(tǒng)壓縮感知理論，提出一種新穎的壓縮采樣匹配追蹤(CoSaMP)［11］算法，將光源重建問題看作一個基追蹤問題，并結(jié)合自適應(yīng)有限元思想。應(yīng)用算法本身的回溯特性，在重建過程中，對離散網(wǎng)格上各節(jié)點光信息向量不斷的比對，尋找包含最佳光源信息的點。與一般的L1正則化方法相比，該算法不僅提高了光源重建的質(zhì)量，同時也提高了雙光源的重建速度。

1 基于壓縮采樣匹配追蹤的光源重建

1．1 重建模型

光在組織中的傳播可以由輻射傳輸模型來描述

對其進(jìn)行有限元離散，拉格朗日二次插值，最后就可得到表面測量和動物體內(nèi)光源分布的線性矩陣方程

加上測量噪聲，我們將其轉(zhuǎn)化為一個約束性最小二乘問題

A∈Rm×n是權(quán)值矩陣，b∈Rm是邊界測量的光密度，x∈Rn是生物體內(nèi)光源分布，m?n，e是測量噪聲。

1.2 光源重建

當(dāng)我們進(jìn)行BLT的重建時，受采集設(shè)備采集方式的限制和光源分布的特點，不會對仿體全部表面進(jìn)行探測，只會對感興趣的部位進(jìn)行探測，所采集到的光能量分布相對于需要重建的全部光源能量分布表現(xiàn)出很大的稀疏性，這使得BLT的重建表現(xiàn)出稀疏重建的特性［7］。而壓縮采樣也是以信號的稀疏性和可壓縮性為前提，即信號當(dāng)中只有s個非零系數(shù)，或者小系數(shù)經(jīng)過變換之后都?xì)w零，這是用CoSaMP方法來重建BLT的基礎(chǔ)。設(shè)s為信號的稀疏度，其數(shù)學(xué)定義如下

基于BLT的稀疏性，可以把重建過程看作一個信號分解問題，對A的列向量vi(1≤i≤N)進(jìn)行歸一化，其集合就構(gòu)成了一個完備字典D，然后就可以將測量數(shù)據(jù)看作是D中s個向量的線性組合。這s個向量的索引構(gòu)成一個集合Γ，因此，b就可以表示為b=AΓxΓ，AΓ是Γ中A的索引的列向量組成的子陣。xΓ表示Γ中x的索引組成的向量。由于y在字典中的稀疏性，可以尋找最少的一組向量{vi|i∈ˉΓ}，使得殘差r=b－AˉΓxˉΓ最小。

在CoSaMP進(jìn)行迭代計算時，每次迭代過程當(dāng)中同時選取多個向量，而后又根據(jù)一定的條件刪除之前選取的一些向量，保留至少一個用作信號重建，它保證在每次選擇時候選集當(dāng)中的向量數(shù)目都少于3 s個，迭代支撐集有s個，每次迭代的時候刪掉不多于2s個向量。

第k次迭代的結(jié)果，都是一個s－稀疏的近似解ak，它和待重建信號之間滿足下面的關(guān)系

v是由于噪聲原因而產(chǎn)生的能量損耗

1.3 算法描述

下面是用CoSaMP方法進(jìn)行光源重建的迭代過程，對其描述如下:

輸入:剛度矩陣A，含噪聲的測量數(shù)據(jù)向量y，稀疏度s;

輸出:待估信號是x，一個s－稀疏逼近a。

1)初始解a0=0，初始?xì)埐顁0=y，計數(shù)器t=1;

2)計算信號代理值 c=ATrt－1;

3)鑒定成分大小，將支撐集合并 Ω=sup p(c2s)，T= Ω ∪ sup p(at－1);

4)利用最小二乘法估算信號 b|T=A+Ty，b|Tc=0來獲得逼近解at=bm;

5)更新殘差 rt=y－Aat;

6)停止條件滿足，輸出a=at，若否，則令t=t+1，返回第2)步。

2 實驗

2.1 實驗設(shè)置

為了驗證本文提出的方法對復(fù)雜動物模型的重建性能，在此我們使用美國南加州大學(xué)提供的通過冷凍切片、CT和PET數(shù)據(jù)構(gòu)造的三維數(shù)字鼠模型［12］進(jìn)行仿真重建實驗。截取其中高32mm的腹部軀干部分作為研究對象，如圖1(a)所示，此數(shù)字鼠圖譜(mouse atlas)的具體光學(xué)特征信息如表1所示。

表1 各器官光學(xué)參數(shù)Tab．1 Optical properties for the organs

對此小鼠軀干模型進(jìn)行有限元離散、網(wǎng)格剖分，其中四面體單元個數(shù)為112 795，網(wǎng)格點數(shù)為21 277，邊界單元數(shù)為17 942，三角形17 942，邊緣單元數(shù)4 992，邊界元素點2 203，相對于總數(shù)，邊界點很稀疏，如圖1(b)所示。

圖1 三維數(shù)字鼠模型Fig．1 3D digitalmousemodel

2．2 實驗對比

一個好的算法，不僅能在單光源的時候精確的重建出光源，也能在雙光源的時候表現(xiàn)出良好的重建特性，為此我們設(shè)置了單光源和雙光源兩組實驗。在單光源實驗中，一個半徑為0.5 mm，高為1mm的圓柱狀光源被放入肝部，其中心位置為(18.1 mm，6.3 mm，15.4 mm)。初始的光源總能量和密度分別為0.785 nW和1nW/mm3。

我們用L1正則化方法中的截斷牛頓法(L1Ls)和本文所提出的CoSaMP算法對單光源進(jìn)行重建，在不同的高斯噪聲水平下，對重建位置誤差、重建能量密度大小以及算法速度進(jìn)行對比，如圖2所示，(a)(b)為2種算法初次網(wǎng)格上的重建效果，(c)(d)為2種算法自適應(yīng)細(xì)分后的網(wǎng)格上的重建效果，圖中黑圈為真實光源位置。重建數(shù)據(jù)如表格2所示。

圖2 單光源的重建Fig．2 The reconstruction of single source

在雙光源實驗中，我們在數(shù)字鼠體內(nèi)設(shè)置兩個光源，光源的各項參數(shù)和單光源相同。光源位置分別為(11.6，10.8，16.4)和 (11.6，6.3，16.4)，2種算法重建效果如圖3所示，(a)(b)為初次網(wǎng)格上的重建效果，(c)(d)為自適應(yīng)細(xì)分后的網(wǎng)格上的重建效果，圖中最內(nèi)側(cè)黑斑為真實光源位置，重建數(shù)據(jù)如表2表3所示。

表2 單光源重建數(shù)據(jù)Tab．2 The data of single source reconstruction

圖3 雙光源重建Fig．3 The reconstruction of single source

我們發(fā)現(xiàn)，在不同的噪聲水平下，在初次和細(xì)分后的網(wǎng)格上，CoSaMP方法的重建能量密度均要比L1Ls方法好，重建位置偏差更小，且在計算速度方面我們發(fā)現(xiàn)CoSaMP方法明顯要優(yōu)于L1Ls方法。

從圖和數(shù)據(jù)中可以看出，在初次網(wǎng)格和自適應(yīng)細(xì)分一次后，CoSaMP方法對2個光源位置的定位比較準(zhǔn)確，重建的雙光源能量密度和重建位置都要優(yōu)于L1Ls方法，這是因為CoSaMP算法受約束等距特性的啟發(fā)，在每次迭代中都要計算信號代理中前2s個元所對應(yīng)的支撐集以及前次計算的逼近解，將它們合并，作為本次候選支撐集，最后用最小二乘法求出其系數(shù)，這種不斷“刪除”和“添加”的回溯思想提高了計算的精確性。

3 結(jié)語

在BLT問題中，由于生物組織的影響，使得光子的傳播過程變得異常復(fù)雜，不可避免地添加了很多噪聲，因此，選取一種既能穩(wěn)定重建光源位置，又能很好抑制噪聲的方法是重建光源算法的核心。本文為了解決BLT重建中存在的困難，基于傳統(tǒng)的壓縮感知理論，提出了一種新穎的壓縮采樣匹配追蹤算法。并通過數(shù)字鼠數(shù)值仿真實驗得以驗證，結(jié)果顯示:在不同噪聲水平下，不論是單光源還是雙光源，這種計算方法都能夠快速而準(zhǔn)確地實現(xiàn)光源重建，重建精度和速度都較L1正則化方法有所提高，并且實驗證明CoSaMP方法具有很好的魯棒性。

表3 雙光源重建數(shù)據(jù)Tab．3 The data of double source reconstruction

［1］ WILLMANN JK，van BRUGGEN N，DINKELBORG L M，et al．Molecular imaging in drug development［J］．Nature Reviews Drug Discovery，2008，7(7):591-607．

［2］ GU X，ZHANG Q，LARCOM L，et al．Three-dimensional bioluminescence tomography with model-based reconstruction［J］．Optics Express，2004，12(17):3996-4000．

［3］ CONGW，WANG G，KUMAR D，et al．Practical reconstruction method for bioluminescence tomography［J］．Opt Express，2005，13(18):6756-6771．

［4］ WANG G，Li Y，JIANG M．Uniqueness theorems in bioluminescence tomography［J］．Medical Physics，2004，31:2289．

［5］ SHEN H，GOLDSTEIN A S，WANG G．Biomedical imaging and image processing in tissue engineering［M］//Tissue Engineering，Berlin Heidelberg:Springer，2011:155-178．

［6］ DEHGHANIH，DAVISSC，JIANG S，et al．Spectrally resolved bioluminescence optical tomography［J］．Optics Letters，2006，31(3):365-367．

［7］ LU Y，ZHANG X，DOURAGHY A，et al．Source reconstruction for spectrally-resolved bioluminescence tomography with sparse a priori information［J］．Optics Express，2009，17(10):8062-8080．

［8］ GAO H，ZHAO H．Multilevel bioluminescence tomography based on radiative transfer equation Part 1:l1 regularization［J］．Opt Express，2010，18(3):1854-1871．

［9］ DONOHO D L，HUO X．Uncertainty principles and ideal atomic decomposition［J］．Information Theory，IEEE Transactions on，2001，47(7):2845-2862．

［10］ CANDESE J，ROMBERG JK，TAO T．Stable signal recovery from incomplete and inaccuratemeasurements［J］．Communications on pure and applied mathematics，2006，59(8):1207-1223．

［11］ NEEDELL D，TROPP JA．CoSaMP:Iterative signal recovery from incomplete and inaccurate samples［J］．Applied and Computational Harmonic Analysis，2009，26(3):301-321．

［12］ DOGDAS B，STOUT D，CHATZIIOANNOU A F，et al．Digimouse:a 3D whole body mouse atlas from CT and cryosection data［J］．Physics in Medicine and Biology，2007，52:577-587．