一種基于dfMC模型的熒光擴(kuò)散光學(xué)斷層圖像重建方法
【技術(shù)領(lǐng)域】
[0001] 本發(fā)明屬于生物醫(yī)學(xué)工程技術(shù)領(lǐng)域，特別是涉及一種基于解耦合熒光蒙特卡羅 dfMC模型的熒光擴(kuò)散光學(xué)斷層圖像重建方法。
【背景技術(shù)】
[0002] 近年來，熒光擴(kuò)散光學(xué)斷層成像技術(shù)已經(jīng)發(fā)展成為了一種重要的成像工具[1]，并 被廣泛應(yīng)用于癌癥診斷、藥物研發(fā)以及基因表達(dá)可視化的研宄。熒光目標(biāo)物的精確重建要 依賴于建立一個(gè)準(zhǔn)確的光子輸運(yùn)模型（即所謂的前向問題）和光子傳輸模型的反演模型 (逆向問題）。前向問題主要是基于能量守恒的玻爾茲曼輻射傳輸框架下的擴(kuò)散近似模型
[2] 和Monte Carlo(MC)模型[3]。由于擴(kuò)散近似模型用于圖像重建時(shí)間代價(jià)低，被廣泛采用； 但現(xiàn)有的擴(kuò)散模型僅僅適用于高散射介質(zhì)，在吸收系數(shù)和散射系數(shù)相當(dāng)，或者源探距離與 平均傳輸自由程相當(dāng)時(shí)，模型就不再適用了。蒙特卡羅模型適應(yīng)任意光學(xué)參數(shù)的活體組織， 和任意復(fù)雜的幾何形狀物體，但由于計(jì)算代價(jià)十分巨大，計(jì)算效率低。一般作為前向模型的 "金標(biāo)準(zhǔn)" [3]。
[0003] 現(xiàn)有的適合于逆向問題求解的前向熒光MC模型主要有伴隨法模型（adjoint MC)
[4] 和微擾法模型（perturbation MC) [5]。伴隨方法假定源和探測(cè)器是可交換的，即具有對(duì)等 性，這種對(duì)等性對(duì)于光纖式探測(cè)系統(tǒng)是成立的，然而對(duì)于自由空間探測(cè)的情形，光源和CCD 的復(fù)雜特性使得這種對(duì)等性很難得到保證。同時(shí)對(duì)于基于自由空間探測(cè)情形的重建效率是 非常低的，大量的CCD像素尺寸大小的探測(cè)器使得蒙特卡羅模擬次數(shù)劇增。而微擾法主要 是通過對(duì)激發(fā)光在組織內(nèi)的MC路徑擾動(dòng)來計(jì)算探測(cè)器上接收到得熒光強(qiáng)度。微擾法它不 受成像系統(tǒng)的限制，具有速度和精度的優(yōu)勢(shì)。成為研宄的熱點(diǎn)。但是微擾方法忽略熒光團(tuán) 對(duì)激發(fā)光的影響，忽略了激發(fā)光與熒光在組織內(nèi)光學(xué)參數(shù)的差異，這諸多的假設(shè)使得微擾 法是一種有偏MC模擬，其熒光統(tǒng)計(jì)量與實(shí)際真實(shí)熒光值會(huì)有偏差，而且這種偏差會(huì)隨著光 學(xué)參數(shù)及光子數(shù)的變化而變化，多數(shù)情況下，誤差會(huì)很大，而且由于光子路徑信息巨大，在 體素模型情形下將會(huì)使得重建效率低下，整個(gè)過程耗時(shí)很大。
[0004] 參考文獻(xiàn)
[0005] [1]V. Ntziachristos? uGoing deeper than microscopy ：the optical imaging frontier in biology，" Nat. Meth.7,603-614(2010).
[0006] [2] X. Song，D. Wang，N. Chen，J. Bai，and H. Wang，"Reconstruction for free-space fluorescence tomography using a novel hybrid adaptive finite element algorithm，" Opt. Express 15,18300-18317 (2007) ·
[0007] [3] J. Swartling? A. Pifferi ? A. Μ. K. Ene jder? and S. Andersson-Engelsi "Accelerated Monte Carlo models to simulate fluorescence spectra from layered tissues，" J. Opt. Soc. Am. A 20, 714-727 (2003) ·
[0008] [4]G. Ma，J. F. Delorme，P. Gallant，and D. A. Boas，"Comparison of simplified Monte Carlo simulation and diffusion approximation for the fluorescence signal from phantoms with typical mouse tissue optical properties, ^Appl.Opt. 46, 1686-1692(2007).
[0009] [5]Sassaroli/iFast perturbation Monte Carlo method for photon migration in heterogeneous turbid media,Opt. Lett. 36, 2095-2097 (2011).

【發(fā)明內(nèi)容】

[0010] 為了克服上述現(xiàn)有技術(shù)的不足，本發(fā)明提供了一種基于解耦合熒光蒙特卡羅 (dfMC)模型的熒光擴(kuò)散光學(xué)斷層圖像重建方法，本發(fā)明的目的是在保持MC模擬"金標(biāo)準(zhǔn)" 的精度，同時(shí)改善MC模擬的耗時(shí)特性，實(shí)現(xiàn)高精度快速熒光斷層圖像重建。
[0011] 本發(fā)明所采用的技術(shù)方案是：
[0012] 一種基于dfMC模型的熒光擴(kuò)散光學(xué)斷層圖像重建方法，其特征在于，包括下列步 驟：
[0013] (1)確定探測(cè)區(qū)域，選擇光源掃描點(diǎn)，并在不同掃描位置下得到探測(cè)器上的熒光漫 射光強(qiáng)分布；
[0014] (2)數(shù)字化待探測(cè)對(duì)象，建立體素模型，依據(jù)掃描光源位置與方向進(jìn)行白蒙特卡羅 模擬，記錄被探測(cè)器收集到的激發(fā)光子的權(quán)重《㈦及在每個(gè)體素內(nèi)的路徑信息和相應(yīng)的物 理量；
[0015] (3)依據(jù)dfMC模型，計(jì)算熒光光子權(quán)重；
[0016] (4)由激發(fā)光子路徑信息和相應(yīng)的物理量建立熒光雅可比矩陣熒光雅可比；
[0017] (5)根據(jù)輸出的熒光光子的權(quán)重和熒光熒光雅可比矩陣，迭代重建計(jì)算出組織內(nèi) 的熒光團(tuán)位置及其吸收系數(shù)。
[0018] 步驟（2)中，采用GPU集群加速整個(gè)白蒙特卡羅過程，記錄到達(dá)探測(cè)器的激發(fā)光子 路徑上的相關(guān)物理量，在體素 V內(nèi)的路徑相關(guān)物理量Llv、L2v分別為：
[0019] Llv =YiP1(SrSm)/Pa (Si-Sm) 1=1
[0020] ?2'=±±1β[α+ι)/ΡΑ[Κ) /=1 j=\
[0021] 式中：η為體素 v內(nèi)的散射次數(shù)，&是體素 v內(nèi)第i次散射的單位散射方向矢量， ?1是各向同性方向函數(shù)，P 4是各向異性方向函數(shù)，L為體素 V內(nèi)的第j次散射路徑長度。
[0022] 步驟（2)中GPU集群加速整個(gè)白蒙特卡羅過程具體如下：
[0023] 2. 1在集群上，依據(jù)光源數(shù)目及集群節(jié)點(diǎn)數(shù)目分配各個(gè)節(jié)點(diǎn)上進(jìn)行MC模擬的次 數(shù)，并把相對(duì)應(yīng)的光源位置和方向信息傳輸?shù)焦?jié)點(diǎn)上；
[0024] 2. 2在每個(gè)節(jié)點(diǎn)上，根據(jù)要模擬的光子數(shù)目分配內(nèi)存，可以將組織模型信息分配到 GPU內(nèi)的紋理寄存器，光學(xué)參數(shù)信息分配到常數(shù)寄存器來提高讀取速度，并根據(jù)GPU的技術(shù) 信息確定最優(yōu)的CUDA運(yùn)行中并行塊和線程的數(shù)目；
[0025] 2. 3在每個(gè)線程上，模擬一個(gè)光子包，并跟蹤光子發(fā)射后在組織內(nèi)散射及吸收的過 程，記錄光子在每個(gè)體素內(nèi)的相關(guān)路徑信息及物理量，光子溢出組織后終止光子模擬。
[0026] 步驟（3)，計(jì)算熒光權(quán)重，(p)計(jì)算公式如下：
【主權(quán)項(xiàng)】
1. 一種基于dfMC模型的熒光擴(kuò)散光學(xué)斷層圖像重建方法，其特征在于，包括下列步 驟： (1) 確定探測(cè)區(qū)域，選擇光源掃描點(diǎn)，并在不同掃描位置下得到探測(cè)器上的熒光漫射光 強(qiáng)分布； (2) 數(shù)字化待探測(cè)對(duì)象，建立體素模型，依據(jù)掃描光源位置與方向進(jìn)行白蒙特卡羅模 擬，記錄被探測(cè)器收集到的激發(fā)光子的權(quán)重<(〃)及在每個(gè)體素內(nèi)的路徑信息和相應(yīng)的物理 量； (3) 依據(jù)dfMC模型，計(jì)算熒光光子權(quán)重； (4) 由激發(fā)光子路徑信息和相應(yīng)的物理量建立熒光雅可比矩陣熒光雅可比； (5) 根據(jù)輸出的熒光光子的權(quán)重和熒光熒光雅可比矩陣，迭代重建計(jì)算出組織內(nèi)的熒 光團(tuán)位置及其吸收系數(shù)。
2. 根據(jù)權(quán)利要求1所述的基于dfMC模型的熒光擴(kuò)散光學(xué)斷層圖像重建方法，其特征在 于步驟（2)中，采用GPU集群加速整個(gè)白蒙特卡羅過程，記錄到達(dá)探測(cè)器的激發(fā)光子路徑上 的相關(guān)物理量，在體素V內(nèi)的路徑相關(guān)物理量Llv、L2v分別為：
式中m為體素V內(nèi)的散射次數(shù)，毛是體素V內(nèi)第i次散射的單位散射方向矢量，P1是各 向同性方向函數(shù)，？&是各向異性方向函數(shù)，L為體素V內(nèi)光子的第j次散射路徑長度。
3. 根據(jù)權(quán)利要求1所述的基于dfMC模型的熒光擴(kuò)散光學(xué)斷層圖像重建方法，其特征在 于：步驟（2)中GPU集群加速整個(gè)白蒙特卡羅過程具體如下： 2. 1在集群上，依據(jù)光源數(shù)目及集群節(jié)點(diǎn)數(shù)目分配各個(gè)節(jié)點(diǎn)上進(jìn)行MC模擬的次數(shù)，并 把相對(duì)應(yīng)的光源位置和方向信息傳輸?shù)焦?jié)點(diǎn)上； 2. 2在每個(gè)節(jié)點(diǎn)上，根據(jù)要模擬的光子數(shù)目分配內(nèi)存，可以將組織模型信息分配到GPU內(nèi)的紋理寄存器，光學(xué)參數(shù)信息分配到常數(shù)寄存器來提高讀取速度，并根據(jù)GPU的技術(shù)信 息確定最優(yōu)的CUDA運(yùn)行中并行塊和線程的數(shù)目； 2. 3在每個(gè)線程上，模擬一個(gè)光子包，并跟蹤光子發(fā)射后在組織內(nèi)散射及吸收的過程， 記錄光子在每個(gè)體素內(nèi)的相關(guān)路徑信息及物理量，光子溢出組織后終止光子模擬。
4. 根據(jù)權(quán)利要求1所述的基于dfMC模型的熒光擴(kuò)散光學(xué)斷層圖像重建方法，其特征在 于步驟（3)，計(jì)算熒光權(quán)重W?(p)計(jì)算公式如下：
式中：P為位置及單位方向的六維矢量，纟為單位方向矢量，k為光子在組織內(nèi)的路徑所 經(jīng)過的體素?cái)?shù)目，m為光子在組織內(nèi)的散射次數(shù)，Uafv是體素V內(nèi)熒光系數(shù)，/C是體素V內(nèi) 的激發(fā)光散射系數(shù)，n是量子產(chǎn)率，Wrai為激發(fā)光子權(quán)重，g/m為體素V內(nèi)熒光的各向異性因 子，g/x為體素V內(nèi)激發(fā)光的各向異性因子，y:為體素V內(nèi)熒光的吸收系數(shù)，y:為體 素V內(nèi)激發(fā)光的吸收系數(shù)，且廠(A = 〇
5. 根據(jù)權(quán)利要求1所述的基于dfMC模型的熒光擴(kuò)散光學(xué)斷層圖像重建方法，其特征在 于步驟（4)中，熒光雅可比矩陣的計(jì)算公式如下：
式中：P為位置及單位方向的六維矢量，丨為單位方向矢量，m為光子在組織內(nèi)的散 射次數(shù)，Uafv是體素V內(nèi)熒光系數(shù)，/C是體素V內(nèi)的激發(fā)光散射系數(shù)，n是量子產(chǎn)率， Wrai為激發(fā)光子權(quán)重，gv?為體素V內(nèi)熒光的各向異性因子，g/x為體素V內(nèi)激發(fā)光的各 向異性因子，yav?為體素V內(nèi)熒光的吸收系數(shù)，yJx為體素V內(nèi)激發(fā)光的吸收系數(shù)，且 廠(凡-凡，/〇 =哪(AHa)，j為節(jié)點(diǎn)i上的光源編號(hào)，d為探測(cè)器編號(hào)，V為需要重建的模 型體素編號(hào)，m為光子在組織內(nèi)的散射次數(shù)，J為熒光雅可比矩陣，心為對(duì)應(yīng)體素V、光源j、 探測(cè)器d的熒光雅可比矩陣元。
6. 根據(jù)權(quán)利要求1所述的基于dfMC模型的快速熒光擴(kuò)散光學(xué)斷層圖像重建方法，其特 征在于步驟（5)整個(gè)熒光重建過程均通過GPU集群加速，各個(gè)節(jié)點(diǎn)上計(jì)算得到的熒光雅可 比矩陣塊，計(jì)算公式如下：
Di= [D y Dlt;…；D yDsJt 式中：s為節(jié)點(diǎn)i上的光源數(shù)目，t為探測(cè)器總數(shù)，n為需要重建的模型體素?cái)?shù)目，DiS節(jié)點(diǎn)i上所有光源下的探測(cè)器探測(cè)到的熒光強(qiáng)度分布與dfMC模擬得到的熒光分布的差值， Dst為節(jié)點(diǎn)i上第s光源下的探測(cè)器中第t個(gè)探測(cè)單元探測(cè)到的熒光強(qiáng)度分布與相對(duì)應(yīng)的 dfMC模擬得到的熒光分布的差值。
7. 根據(jù)權(quán)利要求6所述的基于dfMC模型的快速熒光擴(kuò)散光學(xué)斷層圖像重建方法，其特 征在于步驟（5)，中間量///X,JiDi被傳輸?shù)街鞴?jié)上再進(jìn)行共軛梯度迭代求解熒光分布，這 些中間量比熒光雅可比矩陣塊占用的存貯空間少，迭代計(jì)算公式如下：
式中：i為節(jié)點(diǎn)編號(hào)，N為總節(jié)點(diǎn)數(shù)目，X為熒光團(tuán)吸收系數(shù)分布矩陣Ji1為矩陣J,的轉(zhuǎn) 置矩陣，f⑴為目標(biāo)函數(shù)。
【專利摘要】本發(fā)明公開了一種基于解耦合熒光蒙特卡羅(dfMC)模型的熒光擴(kuò)散光學(xué)斷層圖像重建方法，屬于生物醫(yī)學(xué)工程技術(shù)領(lǐng)域，包括：確定探測(cè)區(qū)域，在探測(cè)區(qū)域選取若干掃描點(diǎn)，在探測(cè)器上獲取相對(duì)應(yīng)的熒光光強(qiáng)分布；建立描述組織光學(xué)參數(shù)空間結(jié)構(gòu)的三維數(shù)字模型，根據(jù)光源掃描的位置和方向進(jìn)行激發(fā)光子的前向白蒙特卡羅(white Monte Carlo)模擬，跟蹤激發(fā)光子并記錄到達(dá)探測(cè)器的光子沿路徑上相應(yīng)的物理量；依據(jù)dfMC方法計(jì)算熒光光子權(quán)重，并計(jì)算出熒光雅可比矩陣；最后通過GPU集群迭代重建計(jì)算出組織內(nèi)的熒光團(tuán)位置及其吸收系數(shù)；本發(fā)明的優(yōu)點(diǎn)是，高精度的dfMC模型，基于GPU集群加速的重建迭代過程，為三維熒光斷層成像系統(tǒng)提供了一種精確快速的重建方法。
【IPC分類】G06T17-00, A61B5-00
【公開號(hào)】CN104706320
【申請(qǐng)?zhí)枴緾N201510050253
【發(fā)明人】駱清銘, 鄧勇, 羅召洋
【申請(qǐng)人】華中科技大學(xué)
【公開日】2015年6月17日
【申請(qǐng)日】2015年1月30日