專利名稱：基于單視圖的多光譜自發(fā)熒光斷層成像重建方法
技術(shù)領(lǐng)域：
本發(fā)明涉及到光學(xué)分子影像成像模態(tài)-自發(fā)熒光斷層成像(BLT)技術(shù)，尤其涉及到一種基于單視圖的多光譜自發(fā)熒光斷層成像重建方法。

背景技術(shù)：
在過去的幾年里，分子影像由于能在體揭示分子和細(xì)胞的信息而受到了越來越多的關(guān)注。因此它已成為一種疾病診斷、藥物療效評價(jià)等的重要工具。作為一種小動(dòng)物分子影像成像模態(tài)，光學(xué)成像技術(shù)特別是自發(fā)熒光斷層成像因?yàn)槠涓咝阅?、低價(jià)格和無創(chuàng)傷等特性已經(jīng)得到了迅速發(fā)展和廣泛的應(yīng)用。自發(fā)熒光斷層成像是最近剛發(fā)展起來的，是一種用來在活體小動(dòng)物體內(nèi)重建自發(fā)熒光光源分布的光學(xué)成像技術(shù)。當(dāng)熒光素與熒光素酶在有氧氣和三磷酸腺苷(ATP)的條件下，就會(huì)發(fā)出熒光。而因?yàn)闊晒馑孛负形灮鹣x素酶(firefly)，磕頭蟲素酶(click beetle)，所以發(fā)出的熒光有不同的光譜，波長一般在400nm-750nm。發(fā)出的熒光穿透生物體而到達(dá)生物體表面，然后用高靈敏度的液氮制冷電荷耦合器件(CCD)獲得。CCD獲得的動(dòng)物表面的數(shù)據(jù)形成了BLT重建的基礎(chǔ)。為采集整個(gè)生物體的表面數(shù)據(jù)，通常將生物體每隔90度旋轉(zhuǎn)一次，用CCD探測器來獲取生物體的前后左右四個(gè)視圖。當(dāng)光源在生物體的位置比較深時(shí)，采集一個(gè)視圖的數(shù)據(jù)需要的時(shí)間大約是5-10分鐘。而當(dāng)采集時(shí)間超過一個(gè)小時(shí)之后熒光強(qiáng)度就會(huì)衰減，所以在多譜的情況下，如果在每個(gè)單譜段都采集四個(gè)視圖的數(shù)據(jù)，那么采集時(shí)間遠(yuǎn)遠(yuǎn)超過1個(gè)小時(shí)，可能就會(huì)使得采集的數(shù)據(jù)不夠準(zhǔn)確。另外，在實(shí)際情況下存在類似像X射線(x-ray)，計(jì)算機(jī)斷層成像(CT)一樣的物理限制，采集數(shù)據(jù)時(shí)就受到角度限制，因此采集到的數(shù)據(jù)僅是動(dòng)物體表面的一部分。多譜采集也會(huì)使得系統(tǒng)矩陣的維數(shù)增加，導(dǎo)致計(jì)算量過大。如何降低測量時(shí)間和減小系統(tǒng)矩陣的維數(shù)是目前的一個(gè)難點(diǎn)。另一方面，在進(jìn)行光源重建時(shí)，因?yàn)锽LT自身具有的不適定特點(diǎn)導(dǎo)致重建結(jié)果不準(zhǔn)確。
BLT的目標(biāo)是能在體實(shí)時(shí)連續(xù)的觀測。因此，發(fā)展快速的重建方法就成為亟待解決的問題。目前的重建方法常用的是有限元方法。理論上有限元網(wǎng)格越細(xì)，得到的結(jié)果越好，但重建的時(shí)間就越長。此外，BLT是病態(tài)問題，如何降低病態(tài)性仍需要進(jìn)一步的探索。


發(fā)明內(nèi)容
本發(fā)明的目的在于克服了自發(fā)熒光斷層成像重建方法重建光源不準(zhǔn)確，重建速度比較慢，不利于實(shí)時(shí)處理以及在多光譜數(shù)據(jù)采集時(shí)可能存在誤差的缺陷，提出了一種基于單視圖的多光譜自發(fā)熒光斷層成像重建方法，該方法基于擴(kuò)散方程模型，并考慮了小動(dòng)物體的非均勻特性，同時(shí)也考慮自發(fā)熒光光源的光譜及實(shí)際應(yīng)用的特點(diǎn)。
為達(dá)到此目的，本發(fā)明提出的基于單視圖的多光譜自發(fā)熒光斷層成像重建方法的實(shí)現(xiàn)步驟主要包括1)數(shù)據(jù)獲??；2)有限元離散化處理；3)優(yōu)化可行光源區(qū)域選擇；4)光源的重建。
1)數(shù)據(jù)獲取是指利用CCD探測器在被測物體表面對逃出的光子進(jìn)行探測，獲得邊界處的光流密度； 有限元離散化處理，是指利用有限元方法將擴(kuò)散方程離散為矩陣方程，并把重建區(qū)域網(wǎng)格剖分為{T1，...，Tk，...}的網(wǎng)格序列； 優(yōu)化可行光源區(qū)域的選擇，主要是在對擴(kuò)散方程離散化處理的基礎(chǔ)上，利用迭代方法快速選擇光源大致存在的區(qū)域，即優(yōu)化可行光源區(qū)域； 光源的重建，主要是在優(yōu)化可行光源區(qū)域選擇的基礎(chǔ)上，確立未知光源密度變量與被測物體表面測量值之間的線性關(guān)系，并利用Tikhonov正則化方法建立目標(biāo)函數(shù)，并最終求解光源密度分布。
本發(fā)明具體包括以下步驟 (1)數(shù)據(jù)獲取 在多光譜下，用濾波片將熒光波長λ離散為m個(gè)波段τ1，...τm，其中τl＝[λl-1，λl)，l＝1，2，...m-1，τl＝[λm-1，λm]，在多光譜情況下m＝5即可提供足夠的先驗(yàn)信息；利用CCD探測器在被測物體的表面

對逃出的光子進(jìn)行探測，獲得每個(gè)波段τl的光流密度Q(x，τl)，其中CCD探測到的表面為x表示在被測物體中的位置；根據(jù)下面公式計(jì)算被測物體表面的光子流密度Φmeas(x，τl) D(x，τl)＝1/(3(μa(x，τl)+(1-g)μs(x，τl)))是生物組織的擴(kuò)散系數(shù)，其中μa(x，τl)是生物組織的吸收系數(shù)，μs(x，τl)是生物組織散射系數(shù)，g是生物組織各向異性參數(shù)；v(x)是被測物體邊界

的單位法向量；A(x；n，n′)≈(1+R(x))/(1-R(x))，當(dāng)外界媒質(zhì)是空氣時(shí)，R(x)可以近似為R(x)≈-1.4399n-2+0.7099n-1+0.6681+0.0636n，其中n為空氣的折射率； (2)有限元離散化處理 光子在生物組織中傳輸?shù)臄?shù)學(xué)模型用下面的擴(kuò)散方程表示 其中Ω是被測物體；

是被測物體的表面；Φ(x，τl)是光子流密度；S(x，τl)是光源密度；根據(jù)有限元理論，得到下面的弱解形式 H1(Ω)是Sobelev空間，Ψ(x，τl)是任意給定的試探函數(shù)。在自適應(yīng)有限元分析的框架下，基于自適應(yīng)的網(wǎng)格細(xì)分，把重建區(qū)域網(wǎng)格剖分表示為{T1，...，Tk，...}的網(wǎng)格序列，其中網(wǎng)格剖分Tk包括

個(gè)離散單元。利用有限元方法，對弱解形式進(jìn)行離散，得到下面矩陣方程 (Kk(τl)+Ck(τl)+Bk(τl))Φk(τl)＝Fk(τl)Sk(τl) 矩陣元素的具體形式為 令Kk(τl)+Ck(τl)+Bk(τl)＝Mk(τl)，考慮到Mk(τl)是稀疏正定矩陣，所以得到 (1)通過刪除 [Mk-1(τl)Fk(τl)]中對應(yīng)非測量點(diǎn)的行得到下面的方程 Φk(τl)＝Ak(τl)Sk(τl) (2) 因此對于m個(gè)波段即得到m個(gè)單一波段的矩陣方程； (3)優(yōu)化可行光源區(qū)域選擇 I)在網(wǎng)格剖分T1上，將上述含有m個(gè)單一波段的矩陣方程組合成含有m個(gè)波段的一個(gè)方程 根據(jù)光源譜信息的性質(zhì)，每個(gè)譜段τl上所占整個(gè)譜段的能量表示為S(τl)＝ω(τl)·S，其中ω(τl)≥0，ω(τl)可通過實(shí)驗(yàn)測定；S表示光源總的能量。在網(wǎng)格剖分T1上，考慮光源譜信息的性質(zhì)，通過方程(2)得到 AS1＝Φ (3) 其中 Φ1(τl)和A1(τl)都有方程(2)計(jì)算得到，S1是網(wǎng)格剖分T1上的未知光源密度變量；在(3)式兩邊同乘AT，得到 ATAS1＝ATΦ(4) AT是A的轉(zhuǎn)置矩陣，AT·A是一個(gè)

的對稱矩陣，所以(4)式是一個(gè)標(biāo)準(zhǔn)線形方程。
II)在網(wǎng)格剖分T1上，確定優(yōu)化光源可行區(qū)域 具體方法為在網(wǎng)格剖分T1上，假定整個(gè)被測物體Ω為光源可行區(qū)域，給定未知光源密度的初始值S10，根據(jù)公式(4)，網(wǎng)格剖分T1上的未知光源密度變量S1用下面的關(guān)系式進(jìn)行迭代求解 其中步長搜索方向n是迭代次數(shù)； 在迭代的過程中，用Φ1meas(τl)代替Φ1(τl)，這是因?yàn)樵肼曉诘倪^程中產(chǎn)生的影響很??；Φ1meas(τl)是網(wǎng)格剖分T1上的邊界測量值，由Φmeas(x，τl)通過最近鄰域插值得到；當(dāng)‖βn‖≤δ或n＞Nmax，迭代停止，此時(shí)便得到了未知光源密度S1的優(yōu)化解S*，優(yōu)化解S*的大致區(qū)域?yàn)棣?，這個(gè)區(qū)域Ω*我們稱之為優(yōu)化的光源可行區(qū)域；在每次迭代的過程中，如果光源的密度小于零，則將其置為零；δ取值范圍為10-5-10-2，Nmax為最大迭代次數(shù)，取值不超過1000。
(4)光源重建 在網(wǎng)格序列{T1，...，Tk，...}上，根據(jù)選擇好的優(yōu)化光源可行區(qū)域Ω*，重新將含有m個(gè)單一波段的矩陣方程組合成含有m個(gè)波段的一個(gè)方程 考慮未知光源密度Sk(τl)與邊界測量值Φkmeas(τl)之間的線性關(guān)系，由方程(1)得到 其中Gk(τl)通過移出[Mk(τl)-1Fk(τl)]中對應(yīng)非測量點(diǎn)的行而得到；考慮到光源譜信息和優(yōu)化可行光源區(qū)域Ω*，得到 其中 Gk表示多譜下的邊界測量值與光源位置量之間的關(guān)系，Wk為k層網(wǎng)格上的對角矩陣，Φkmeas(τl)是波段τl在網(wǎng)格剖分Tk上的邊界測量值，由Φmeas(x，τl)通過最近鄰域插值得到；當(dāng)k＝1時(shí)，sk(i)是通過迭代方法得到的優(yōu)化解S*，skmax是優(yōu)化解S的最大值；當(dāng)k≥2時(shí)，sk(i)和skmax是網(wǎng)格剖分Tk上的未知光源密度的初始值和初始值的最大值，通過插值網(wǎng)格剖分Tk-1上的重建結(jié)果獲得，即Sk-1是在網(wǎng)格剖分Tk-1上重建得到的結(jié)果，Ik-1k是插值因子，通過子單元繼承父單元的重建結(jié)果值的方式來實(shí)現(xiàn)；γk是比例因子，是隨k變化的分段常數(shù)；最終通過保留GkWk中對應(yīng)可行光源區(qū)域的列，得到光源可行區(qū)域的未知光源密度變量Skp與表面測量值Φkmeas在網(wǎng)格剖分Tk上的關(guān)系 II)利用Tikhonov正則化方法建立目標(biāo)函數(shù)并優(yōu)化從而得到重建結(jié)果 上一步雖然得到了未知光源密度變量Skp與表面測量值Φkmeas之間的關(guān)系，由于Ak矩陣的病態(tài)，不能通過直接求解的方法得到重建結(jié)果，因此利用Tikhonov正則的方法，并考慮未知光源密度變量不能為負(fù)值的特性，得到在網(wǎng)格剖分Tk上的目標(biāo)函數(shù) 考慮到多光譜以及非勻質(zhì)模型會(huì)極大的增加計(jì)算量，因此采用譜梯度大規(guī)模優(yōu)化算法對該目標(biāo)函數(shù)進(jìn)行優(yōu)化；Ssupk是光源密度上界，經(jīng)驗(yàn)取值為不大于1000，‖V‖Λ＝VTΛV，λk正則化系數(shù)，經(jīng)驗(yàn)取值10-9-10-6，Skinit是網(wǎng)格剖分Tk上的未知光源密度初始值，取值范圍為10-7-10-4；在判斷優(yōu)化過程是否中止時(shí)，采用當(dāng)前優(yōu)化梯度

以及優(yōu)化迭代次數(shù)Nki作為評判準(zhǔn)則，即當(dāng)或時(shí)，優(yōu)化停止，得到重建結(jié)果，即光源可行區(qū)域的光源密度分布；其中εg和Nkmax分別是梯度閾值和最大優(yōu)化迭代次數(shù)，可根據(jù)實(shí)驗(yàn)經(jīng)驗(yàn)獲得，取值不超過10-6和104； III)當(dāng)優(yōu)化完成后，判斷重建是否停止 利用優(yōu)化求得的光源密度求解任意一個(gè)譜段τl上的光密度分布Φkc(τl)，并計(jì)算‖Φkc(τl)-Φkmeas(τl)‖，當(dāng)或k≥Kmax，重建停止；Φkc(τl)是通過將優(yōu)化求得的光源密度代入擴(kuò)散方程而求得的邊界處光子流密度，Φkmeas(τl)是波段τl在網(wǎng)格剖分Tk上邊界測量點(diǎn)的測量值，由Φmeas(x，τl)通過最近鄰域插值得到；εΦ經(jīng)驗(yàn)取值范圍10-9-10-7；Kmax是網(wǎng)格最大剖分次數(shù)；如果不能滿足或k≥Kmax，根據(jù)重建的結(jié)果分別用直接最大值選擇方法和等級值虧修正基的誤差估計(jì)方法選取可行和禁止光源區(qū)域；通過以上方法對可行和禁止光源區(qū)域內(nèi)的網(wǎng)格單元的選擇，對區(qū)域剖分網(wǎng)格進(jìn)行自適應(yīng)的細(xì)分；并從網(wǎng)格剖分Tk轉(zhuǎn)到Tk+1，然后轉(zhuǎn)到第(4)(I)步，重復(fù)前面的步驟，直至重建停止； 本發(fā)明基于單視圖測量數(shù)據(jù)，利用多光譜提供的先驗(yàn)信息，通過迭代方法進(jìn)行可行光源區(qū)域的選擇，重建光源的不準(zhǔn)確問題得到了很好的解決，重建的速度和質(zhì)量由于光源可行區(qū)域的使用而被進(jìn)一步提高，此方法對自發(fā)熒光斷層成像的發(fā)展具有重要的應(yīng)用價(jià)值。



圖1.非均勻仿體 圖中(a)非均勻仿體，1代表肌肉、2骨骼、3心臟、4肺、5肝和 6單光源；(b)重建算法用到的初始網(wǎng)格；(c)包含肌肉、骨骼、肺的橫截面，光源在右肺中；箭頭所示的方向?yàn)镃CD探測的方向，7前視圖方向，8左視圖方向，9后視圖方向，10右視圖方向。
圖2.程序主流程圖。
圖3.確立光源可行區(qū)域的子流程圖。
圖4.確定光源可行區(qū)域的算法示意圖。
圖5.重建結(jié)果示意圖。
表1.各組織在不同波長下的光學(xué)特性參數(shù)。

具體實(shí)施例方式 下面結(jié)合附圖詳細(xì)描述本實(shí)施例。
(1)數(shù)據(jù) 獲取 為了驗(yàn)證本方法，我們設(shè)計(jì)了非勻質(zhì)模型來近似模仿小鼠腹部內(nèi)的組織，如圖1(a)所示，生物組織分別為肌肉、骨骼、心臟、肺、肝，各個(gè)組織的光學(xué)特性參數(shù)如表1所示。在此實(shí)驗(yàn)中，我們將重建自發(fā)光源的光源密度分布作為重建目標(biāo)，在重建實(shí)驗(yàn)中，這個(gè)光源大小為半徑1毫米，光源密度為0.238nano-Watts/mm3，位置為(-3，5，15)。它的波譜范圍從500納米到750納米，為了進(jìn)行多光譜實(shí)驗(yàn)，我們依據(jù)當(dāng)前的探測水平，把這一波譜范圍分成了5個(gè)譜段進(jìn)行分別探測，分別為τ1＝[500，550)nm，τ2＝[550，600)nm，τ3＝[600，650)nm，τ4＝[650，700)nm，τ5＝[700，750)nm。邊界測量數(shù)據(jù)由蒙特卡羅方法產(chǎn)生，并按照圖1(c)所示的方法取前視圖的邊界數(shù)據(jù)，并求得每個(gè)波段的光子流密度Φmeas(x，τl)； (2)有限 元離散化處理 光子在生物組織中傳輸?shù)臄?shù)學(xué)模型用下面的擴(kuò)散方程表示 其中Ω是被測物體；

是被測物體的表面；Φ(x，τl)是光子流密度；S(x，τl)是光源密度；根據(jù)有限元理論，得到下面的弱解形式 H1(Ω)是Sobelev空間，Ψ(x，τl)是任意的試探函數(shù)。在自適應(yīng)有限元分析的框架下，基于自適應(yīng)的網(wǎng)格細(xì)分，把重建區(qū)域網(wǎng)格剖分表示為{T1，...，Tk，...}的網(wǎng)格序列，其中網(wǎng)格剖分Tk包括

個(gè)離散單元，網(wǎng)格剖分T1包含6878個(gè)離散單元，如圖1(b)所示。利用有限元方法，對弱解形式進(jìn)行離散，在每個(gè)波段τl(l＝1，2，3，4，5)上得到下面矩陣方程 (Kk(τl)+Ck(τl)+Bk(τl))Φk(τl)＝Fk(τl)Sk(τl) 矩陣元素的具體形式為 令Kk(τl)+Ck(τl)+Bk(τl)＝Mk(τl)，考慮到Mk(τl)是稀疏正定矩陣，所以得到 (1)通過刪除 [Mk-1(τl)Fk(τl)]中對應(yīng)非測量點(diǎn)的行得到下面的方程 Φk(τl)＝Ak(τl)Sk(τl) (2) 對于5個(gè)波段得到5個(gè)單一波段的矩陣方程； (3)優(yōu)化可行光源區(qū)域選擇 I)在網(wǎng)格剖分T1上，將上述含有5個(gè)單一波段的矩陣方程組合成含有5個(gè)波段的一個(gè)方程 根據(jù)光源譜信息的性質(zhì)，每個(gè)譜段τl上所占整個(gè)譜段的能量表示為S(τl)＝ω(τl)·S，ω(τl)≥0，其中S表示光源總的光源能量，值為1；ω(τl)通過實(shí)驗(yàn)測定分別為ω(τ1)＝0.222，ω(τ2)＝0.167，ω(τ3)＝0.222，ω(τ4)＝0.167，ω(τ5)＝0.222。
在網(wǎng)格剖分T1上，考慮光源譜信息的性質(zhì)，通過方程(2)得到 AS1＝Φ (3) 其中 Φ1(τl)和A1(τl)都有方程(2)計(jì)算得到，S1是網(wǎng)格剖分T1上的未知光源密度變量。在(3)式兩邊同乘AT，得到 ATAS1＝ATΦ (4) AT是A的轉(zhuǎn)置矩陣，AT·A是一個(gè)34390×34390的對稱矩陣，所以(4)式是一個(gè)標(biāo)準(zhǔn)線形方程。
III)在網(wǎng)格剖分T1上，確定優(yōu)化光源可行區(qū)域 具體方法為在網(wǎng)格剖分T1上，假定整個(gè)被測物體Ω為光源可行區(qū)域，給定未知光源密度的初始值根據(jù)公式(4)，網(wǎng)格剖分T1上的未知光源密度變量S1根據(jù)下面的關(guān)系式進(jìn)行迭代求解 其中步長搜索方向n是迭代次數(shù)。
在迭代的過程中，用Φ1meas(τl)代替Φ1(τl)，Φ1meas(τl)是τl在網(wǎng)格剖分T1上的邊界測量值，由Φmeas(x，τl)通過最近鄰域插值得到；參數(shù)取值分別為δ＝5×10-3，Nmax＝500。經(jīng)過5次迭代后，‖βn‖＝4.9e-5，‖βn‖≤δ，所以迭代停止，此時(shí)便得到了未知光源密度S1的優(yōu)化解S*，優(yōu)化解S*的大致區(qū)域?yàn)棣?，這個(gè)區(qū)域Ω*我們稱之為優(yōu)化的光源可行區(qū)域。在每次迭代的過程中，如果光源的密度小于零，則將其置為零。
(4)光源重建 (I)在被測物體Ω的網(wǎng)格剖分{T1，...，Tk，...}的網(wǎng)格序列上，根據(jù)選擇好的優(yōu)化光源可行區(qū)域Ω*，重新將含有5個(gè)波段的矩陣方程組合成含有5個(gè)波段的一個(gè)方程 考慮未知光源密度Sk(τl)與邊界測量值Φkmeas(τl)之間的線性關(guān)系，由方程(1)得到 其中Gk(τl)通過移出[Mk(τl)-1Fk(τl)]中對應(yīng)非測量點(diǎn)的行而得到?？紤]到光源譜信息和優(yōu)化可行光源區(qū)域Ω*，得到 其中 Gk表示多譜下的邊界測量值與光源位置量之間的關(guān)系，Wk為k層網(wǎng)格上的對角矩陣，Φkmeas(τl)是波段τl在網(wǎng)格剖分Tk上的邊界測量值。當(dāng)k＝1時(shí)，sk(i)是通過迭代方法得到的優(yōu)化解S*，skmax是優(yōu)化解S*的最大值；當(dāng)k≥2時(shí)，sk(i)和skmax是網(wǎng)格剖分Tk上的初始值和初始值的最大值，通過插值網(wǎng)格剖分Tk-1上的重建結(jié)果得到，即Sk-1是在網(wǎng)格剖分Tk-1上重建得到的結(jié)果，Ik-1k是插值因子，通過子單元繼承父單元的重建結(jié)果值的方式來實(shí)現(xiàn)；γk是比例因子，是隨k變化的分段常數(shù)；k＝1時(shí)，參數(shù)γk＝0；k＞1時(shí)，γk初始為0.1，并隨著k增加，有γk+1＝10γk。最終通過保留GkWk中對應(yīng)可行光源區(qū)域的列，得到光源可行區(qū)域的未知光源密度變量Skp與表面測量值Φkmeas在網(wǎng)格剖分Tk上的關(guān)系 (II)利用Tikhonov正則化方法建立目標(biāo)函數(shù)并優(yōu)化從而得到重建結(jié)果 上一步雖然得到了未知光源密度變量Skp與表面測量值Φkmeas之間的關(guān)系，由于Ak矩陣的病態(tài)，不能通過直接求解的方法得到重建結(jié)果，利用Tikhonov正則的方法，并考慮未知光源密度變量不能為負(fù)值的特性，得到在網(wǎng)格剖分Tk上的目標(biāo)函數(shù) 考慮到多光譜以及非勻質(zhì)模型會(huì)極大的增加計(jì)算量，因此采用譜梯度大規(guī)模優(yōu)化算法對該目標(biāo)函數(shù)進(jìn)行優(yōu)化。Ssupk是光源密度上界，取值為100，‖V‖Λ＝VTΛV，λk正則化系數(shù)，取值1.0×10-8，Skinit是在網(wǎng)格剖分Tk上未知光源密度變量的初始值，取值1.0×10-5。在判斷優(yōu)化過程是否中止時(shí)，采用了當(dāng)前優(yōu)化梯度

以及優(yōu)化迭代次數(shù)Nki作為評判準(zhǔn)則，即當(dāng)或時(shí)，優(yōu)化停止，得到重建結(jié)果，即光源可行區(qū)域的光源密度分布，其中εg是梯度閾值，取值1.0×10-8，Nkmax是最大優(yōu)化迭代次數(shù)，取值1.0×104。
(III)當(dāng)優(yōu)化完成后，判斷重建是否停止 利用優(yōu)化得到的光源密度分布求解任意一個(gè)譜段τl上的光子流密度Φkc(τl)，并計(jì)算‖Φkc(τl)-Φkmeas(τl)‖，當(dāng)(εΦ為正數(shù))或k≥Kmax(Kmax是網(wǎng)格最大剖分次數(shù))，重建停止。Φkc(τl)是通過將優(yōu)化求得的光源密度代入擴(kuò)散方程而求得的邊界處光子流密度，Φkmeas(τl)是波段τl在第k層網(wǎng)格上邊界測量點(diǎn)的測量值；如果不能滿足或k≥Kmax，利用重建的結(jié)果分別用直接最大值選擇方法和等級值虧修正基的誤差估計(jì)方法選取可行和禁止光源區(qū)域。通過以上方法對可行和禁止光源區(qū)域內(nèi)的網(wǎng)格單元的選擇，并對區(qū)域剖分網(wǎng)格進(jìn)行自適應(yīng)的細(xì)分，從網(wǎng)格剖分Tk轉(zhuǎn)到Tk+1，然后轉(zhuǎn)到第(4)(I)步，重復(fù)前面的步驟，直至重建停止。
停止閾值εΦ和網(wǎng)格最大細(xì)分次數(shù)Kmax分別是1×10-8和3。當(dāng)在網(wǎng)格剖分T1上重建完成后，計(jì)算在τ1＝[500-550)nm的光密度分布Φkc(τl)，求得此時(shí)(k＝1)＜(Kmax＝3)所以重建未停止；利用網(wǎng)格剖分T1上的重建結(jié)果分別用直接最大值選擇方法和等級值虧修正基的誤差估計(jì)方法選取可行和禁止光源區(qū)域。通過以上方法對可行和禁止光源區(qū)域內(nèi)的網(wǎng)格單元的選擇，并對區(qū)域剖分網(wǎng)格進(jìn)行自適應(yīng)的細(xì)分，并從網(wǎng)格剖分T1轉(zhuǎn)到網(wǎng)格剖分T2，然后轉(zhuǎn)到第(4)(I)步，重復(fù)前面的步驟；優(yōu)化完成后計(jì)算在τ1＝[500-550)nm的光密度分布Φkc(τl)，求得因此重建停止，重建的結(jié)果如圖5所示，重建的位置為(-1.53，5.3，14.98)，此時(shí)重建的最大光源密度為0.224nano-Watts/mm3。
表1

權(quán)利要求
1.一種基于單視圖的多光譜自發(fā)熒光斷層成像重建方法，其特征在于，包括以下步驟
1)數(shù)據(jù)獲取
在多光譜下，用濾波片將熒光波長λ離散為m個(gè)波段τ1，...，τm，其中τl＝[λl-1，λl)，l＝1，2，...，m-1，τl＝[λm-1，λm]；利用CCD探測器在被測物體的表面
對逃出的光子進(jìn)行探測，獲得每個(gè)波段τl的光流密度Q(x，τl)，其中CCD探測到的表面為γ，x表示在被測物體中的位置；根據(jù)下面公式計(jì)算被測物體表面的光子流密度Φmeas(x，τl)
D(x，τl)＝1/(3(μa(x，τl)+(1-g)μs(x，τl)))是生物組織的擴(kuò)散系數(shù)，其中μa(x，τl)是生物組織的吸收系數(shù)，μs(x，τl)是生物組織散射系數(shù)，g是生物組織各向異性參數(shù)；v(x)是被測物體邊界
的單位法向量；A(x；n，n′)≈(1+R(x))/(1-R(x))，當(dāng)外界媒質(zhì)是空氣時(shí)，R(x)近似為R(x)≈-1.4399n-2+0.7099n-1+0.6681+0.0636n，其中n為空氣的折射率；
2)有限元離散化處理
光子在生物組織中傳輸?shù)臄?shù)學(xué)模型用下面的擴(kuò)散方程表示
其中Ω是被測物體；
是被測物體的表面；Φ(x，τl)是光子流密度；S(x，τl)是光源密度；根據(jù)有限元理論，得到下面的弱解形式
H1(Ω)是Sobelev空間，ψ(x，τl)是任意給定的試探函數(shù)；在自適應(yīng)有限元分析的框架下，基于自適應(yīng)的網(wǎng)格細(xì)分，把重建區(qū)域網(wǎng)格剖分表示為{T1，...，Tk，...}的網(wǎng)格序列，其中網(wǎng)格剖分Tk包括
個(gè)離散單元；利用有限元方法，對弱解形式進(jìn)行離散，得到下面矩陣方程
(Kk(τl)+Ck(τl)+Bk(τl))Φk(τl)＝Fk(τl)Sk(τl)
矩陣元素的具體形式為
令Kk(τl)+Ck(τl)+Bk(τl)＝Mk(τl)，考慮到Mk(τl)是稀疏正定矩陣，所以得到
(1)通過刪除
[Mk-1(τl)Fk(τl)]中對應(yīng)非測量點(diǎn)的行得到下面的方程
Φk(τl)＝Ak(τl)Sk(τl)
因此對于m個(gè)波段即得到m個(gè)單一波段的矩陣方程；
3)優(yōu)化可行光源區(qū)域選擇
I)在網(wǎng)格剖分T1上，將上述含有m個(gè)單一波段的矩陣方程組合成含有m個(gè)波段的一個(gè)方程
根據(jù)光源譜信息的性質(zhì)，每個(gè)譜段τl上所占整個(gè)譜段的能量表示為S(τl)＝ω(τl)·S，其中ω(τl)≥0，ω(τl)通過實(shí)驗(yàn)測定；S表示光源總的能量；在網(wǎng)格剖分T1上，考慮光源譜信息的性質(zhì)，通過方程(2)得到
AS1＝Φ (3)
其中
Φ1(τl)和A1(τl)都有方程(2)計(jì)算得到，S1是網(wǎng)格剖分T1上的未知光源密度變量；在(3)式兩邊同乘AT，得到
ATAS1＝ATΦ (4)
AT是A的轉(zhuǎn)置矩陣，AT·A是一個(gè)
的對稱矩陣，所以(4)式是一個(gè)標(biāo)準(zhǔn)線形方程；
II)在網(wǎng)格剖分T1上，確定優(yōu)化光源可行區(qū)域
具體方法為在網(wǎng)格剖分T1上，假定整個(gè)被測物體Ω為光源可行區(qū)域，給定未知光源密度的初始值S10，根據(jù)公式(4)，網(wǎng)格剖分T1上的未知光源密度變量S1用下面的關(guān)系式進(jìn)行迭代求解
其中步長搜索方向n是迭代次數(shù)；
在迭代的過程中，用Φ1meas(τl)代替Φ1(τl)，Φ1meas(τl)是網(wǎng)格剖分T1上的邊界測量值，由Φmeas(x，τl)通過最近鄰域插值得到；當(dāng)||βn||≤δ或n＞Nmax，迭代停止，此時(shí)便得到了未知光源密度S1的優(yōu)化解S*，優(yōu)化解S*的大致區(qū)域?yàn)棣?，這個(gè)區(qū)域Ω*我們稱之為優(yōu)化的光源可行區(qū)域；在每次迭代的過程中，如果光源的密度小于零，則將其置為零；δ取值范圍為10-5-10-2，Nmax為最大迭代次數(shù)，取值不超過1000；
4)光源重建
I)在網(wǎng)格序列{T1，...，Tk，...}上，根據(jù)選擇好的優(yōu)化光源可行區(qū)域Ω*，重新將含有m個(gè)單一波段的矩陣方程組合成含有m個(gè)波段的一個(gè)方程
考慮未知光源密度Sk(τl)與邊界測量值Φkmeas(τl)之間的線性關(guān)系，由方程(1)得到
其中Gk(τl)通過移出[Mk(τl)-1Fk(τl)]中對應(yīng)非測量點(diǎn)的行而得到；考慮到光源譜信息和優(yōu)化可行光源區(qū)域Ω*，得到
其中
Gk表示多譜下的邊界測量值與光源位置量之間的關(guān)系，Wk為k層網(wǎng)格上的對角矩陣，Φkmeas(τl)是波段τl在網(wǎng)格剖分Tk上的邊界測量值，由Φmeas(x，τl)通過最近鄰域插值得到；當(dāng)k＝1時(shí)，sk(i)是通過迭代方法得到的優(yōu)化解S*，skmax是優(yōu)化解S*的最大值；當(dāng)k≥2時(shí)，sk(i)和skmax是網(wǎng)格剖分Tk上的未知光源密度的初始值和初始值的最大值，通過插值網(wǎng)格剖分Tk-1上的重建結(jié)果獲得，即Sk-1是在網(wǎng)格剖分Tk-1上重建得到的結(jié)果，Ik-1k是插值因子，通過子單元繼承父單元的重建結(jié)果值的方式來實(shí)現(xiàn)；γk是比例因子，是隨k變化的分段常數(shù)；最終通過保留GkWk中對應(yīng)可行光源區(qū)域的列，得到光源可行區(qū)域的未知光源密度變量Skp與表面測量值Φkmeas在網(wǎng)格剖分Tk上的關(guān)系
II)利用Tikhonov正則化方法建立目標(biāo)函數(shù)并優(yōu)化從而得到重建結(jié)果
上一步雖然得到了未知光源密度變量Skp與表面測量值Φkmeas之間的關(guān)系，由于Ak矩陣的病態(tài)，不能通過直接求解的方法得到重建結(jié)果，因此利用Tikhonov正則的方法，并考慮未知光源密度變量不能為負(fù)值的特性，得到在網(wǎng)格剖分Tk上的目標(biāo)函數(shù)
考慮到多光譜以及非勻質(zhì)模型會(huì)極大的增加計(jì)算量，因此采用譜梯度大規(guī)模優(yōu)化算法對該目標(biāo)函數(shù)進(jìn)行優(yōu)化；Ssupk是光源密度上界，經(jīng)驗(yàn)取值為不大于1000，||V||Λ＝VTΛV，λk正則化系數(shù)，經(jīng)驗(yàn)取值10-9-10-6，Skinit是網(wǎng)格剖分Tk上的未知光源密度初始值，取值范圍為10-7-10-4；在判斷優(yōu)化過程是否中止時(shí)，采用當(dāng)前優(yōu)化梯度
以及優(yōu)化迭代次數(shù)Nki作為評判準(zhǔn)則，即當(dāng)或時(shí)，優(yōu)化停止，得到重建結(jié)果，即光源可行區(qū)域的光源密度分布；其中εg和Nkmax分別是梯度閾值和最大優(yōu)化迭代次數(shù)，根據(jù)實(shí)驗(yàn)經(jīng)驗(yàn)獲得，取值不超過10-6和104；
III)當(dāng)優(yōu)化完成后，判斷重建是否停止
利用優(yōu)化得到的光源密度求解任意一個(gè)譜段τl上的光密度分布Φkc(τl)，并計(jì)算||Φkc(τl)-Φkmeas(τl)||，當(dāng)或k≥Kmax，重建停止；Φkc(τl)是通過將優(yōu)化求得的光源密度代入擴(kuò)散方程而求得的邊界處光子流密度，Φkmeas(τl)是波段τl在網(wǎng)格剖分Tk上邊界測量點(diǎn)的測量值，由Φmeas(x，τl)通過最近鄰域插值得到；εΦ經(jīng)驗(yàn)取值范圍10-7-10-9；Kmax是網(wǎng)格最大剖分次數(shù)；如果不能滿足或k≥Kmax，根據(jù)重建的結(jié)果分別用直接最大值選擇方法和等級值虧修正基的誤差估計(jì)方法選取可行和禁止光源區(qū)域；通過以上方法對可行和禁止光源區(qū)域內(nèi)的網(wǎng)格單元的選擇，對區(qū)域剖分網(wǎng)格進(jìn)行自適應(yīng)的細(xì)分；并從網(wǎng)格剖分Tk轉(zhuǎn)到Tk+1，然后轉(zhuǎn)到第4)(I)步，直至重建停止。
全文摘要
基于單視圖的多光譜自發(fā)熒光斷層成像重建方法屬于光學(xué)分子影像成像領(lǐng)域。該方法基于擴(kuò)散方程模型，并考慮了小動(dòng)物體的非均勻特性，同時(shí)也考慮自發(fā)熒光光源的光譜及實(shí)際應(yīng)用的特點(diǎn)。為達(dá)到此目的，本發(fā)明提出的基于單視圖的多光譜自發(fā)熒光斷層成像重建方法的實(shí)現(xiàn)步驟主要包括1)數(shù)據(jù)獲?。?)有限元離散化處理；3)優(yōu)化可行光源區(qū)域選擇；4)光源的重建。本發(fā)明克服了自發(fā)熒光斷層成像重建方法重建光源不準(zhǔn)確，重建速度比較慢，不利于實(shí)時(shí)處理以及在多光譜數(shù)據(jù)采集時(shí)可能存在誤差的缺陷。
文檔編號A61B5/00GK101342075SQ200810116818
公開日2009年1月14日 申請日期2008年7月18日 優(yōu)先權(quán)日2008年7月18日
發(fā)明者賈克斌, 馮金超, 捷 田 申請人:北京工業(yè)大學(xué)