一種基于dfMC模型的熒光擴散光學斷層圖像重建方法
【技術領域】
[0001] 本發(fā)明屬于生物醫(yī)學工程技術領域，特別是涉及一種基于解耦合熒光蒙特卡羅 dfMC模型的熒光擴散光學斷層圖像重建方法。
【背景技術】
[0002] 近年來，熒光擴散光學斷層成像技術已經(jīng)發(fā)展成為了一種重要的成像工具[1]，并 被廣泛應用于癌癥診斷、藥物研發(fā)以及基因表達可視化的研宄。熒光目標物的精確重建要 依賴于建立一個準確的光子輸運模型（即所謂的前向問題）和光子傳輸模型的反演模型 (逆向問題）。前向問題主要是基于能量守恒的玻爾茲曼輻射傳輸框架下的擴散近似模型
[2] 和Monte Carlo(MC)模型[3]。由于擴散近似模型用于圖像重建時間代價低，被廣泛采用； 但現(xiàn)有的擴散模型僅僅適用于高散射介質(zhì)，在吸收系數(shù)和散射系數(shù)相當，或者源探距離與 平均傳輸自由程相當時，模型就不再適用了。蒙特卡羅模型適應任意光學參數(shù)的活體組織， 和任意復雜的幾何形狀物體，但由于計算代價十分巨大，計算效率低。一般作為前向模型的 "金標準" [3]。
[0003] 現(xiàn)有的適合于逆向問題求解的前向熒光MC模型主要有伴隨法模型（adjoint MC)
[4] 和微擾法模型（perturbation MC) [5]。伴隨方法假定源和探測器是可交換的，即具有對等 性，這種對等性對于光纖式探測系統(tǒng)是成立的，然而對于自由空間探測的情形，光源和CCD 的復雜特性使得這種對等性很難得到保證。同時對于基于自由空間探測情形的重建效率是 非常低的，大量的CCD像素尺寸大小的探測器使得蒙特卡羅模擬次數(shù)劇增。而微擾法主要 是通過對激發(fā)光在組織內(nèi)的MC路徑擾動來計算探測器上接收到得熒光強度。微擾法它不 受成像系統(tǒng)的限制，具有速度和精度的優(yōu)勢。成為研宄的熱點。但是微擾方法忽略熒光團 對激發(fā)光的影響，忽略了激發(fā)光與熒光在組織內(nèi)光學參數(shù)的差異，這諸多的假設使得微擾 法是一種有偏MC模擬，其熒光統(tǒng)計量與實際真實熒光值會有偏差，而且這種偏差會隨著光 學參數(shù)及光子數(shù)的變化而變化，多數(shù)情況下，誤差會很大，而且由于光子路徑信息巨大，在 體素模型情形下將會使得重建效率低下，整個過程耗時很大。
[0004] 參考文獻
[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)有技術的不足，本發(fā)明提供了一種基于解耦合熒光蒙特卡羅 (dfMC)模型的熒光擴散光學斷層圖像重建方法，本發(fā)明的目的是在保持MC模擬"金標準" 的精度，同時改善MC模擬的耗時特性，實現(xiàn)高精度快速熒光斷層圖像重建。
[0011] 本發(fā)明所采用的技術方案是：
[0012] 一種基于dfMC模型的熒光擴散光學斷層圖像重建方法，其特征在于，包括下列步 驟：
[0013] (1)確定探測區(qū)域，選擇光源掃描點，并在不同掃描位置下得到探測器上的熒光漫 射光強分布；
[0014] (2)數(shù)字化待探測對象，建立體素模型，依據(jù)掃描光源位置與方向進行白蒙特卡羅 模擬，記錄被探測器收集到的激發(fā)光子的權重《㈦及在每個體素內(nèi)的路徑信息和相應的物 理量；
[0015] (3)依據(jù)dfMC模型，計算熒光光子權重；
[0016] (4)由激發(fā)光子路徑信息和相應的物理量建立熒光雅可比矩陣熒光雅可比；
[0017] (5)根據(jù)輸出的熒光光子的權重和熒光熒光雅可比矩陣，迭代重建計算出組織內(nèi) 的熒光團位置及其吸收系數(shù)。
[0018] 步驟（2)中，采用GPU集群加速整個白蒙特卡羅過程，記錄到達探測器的激發(fā)光子 路徑上的相關物理量，在體素 V內(nèi)的路徑相關物理量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集群加速整個白蒙特卡羅過程具體如下：
[0023] 2. 1在集群上，依據(jù)光源數(shù)目及集群節(jié)點數(shù)目分配各個節(jié)點上進行MC模擬的次 數(shù)，并把相對應的光源位置和方向信息傳輸?shù)焦?jié)點上；
[0024] 2. 2在每個節(jié)點上，根據(jù)要模擬的光子數(shù)目分配內(nèi)存，可以將組織模型信息分配到 GPU內(nèi)的紋理寄存器，光學參數(shù)信息分配到常數(shù)寄存器來提高讀取速度，并根據(jù)GPU的技術 信息確定最優(yōu)的CUDA運行中并行塊和線程的數(shù)目；
[0025] 2. 3在每個線程上，模擬一個光子包，并跟蹤光子發(fā)射后在組織內(nèi)散射及吸收的過 程，記錄光子在每個體素內(nèi)的相關路徑信息及物理量，光子溢出組織后終止光子模擬。
[0026] 步驟（3)，計算熒光權重，(p)計算公式如下：
【主權項】
1. 一種基于dfMC模型的熒光擴散光學斷層圖像重建方法，其特征在于，包括下列步 驟： (1) 確定探測區(qū)域，選擇光源掃描點，并在不同掃描位置下得到探測器上的熒光漫射光 強分布； (2) 數(shù)字化待探測對象，建立體素模型，依據(jù)掃描光源位置與方向進行白蒙特卡羅模 擬，記錄被探測器收集到的激發(fā)光子的權重<(〃)及在每個體素內(nèi)的路徑信息和相應的物理 量； (3) 依據(jù)dfMC模型，計算熒光光子權重； (4) 由激發(fā)光子路徑信息和相應的物理量建立熒光雅可比矩陣熒光雅可比； (5) 根據(jù)輸出的熒光光子的權重和熒光熒光雅可比矩陣，迭代重建計算出組織內(nèi)的熒 光團位置及其吸收系數(shù)。
2. 根據(jù)權利要求1所述的基于dfMC模型的熒光擴散光學斷層圖像重建方法，其特征在 于步驟（2)中，采用GPU集群加速整個白蒙特卡羅過程，記錄到達探測器的激發(fā)光子路徑上 的相關物理量，在體素V內(nèi)的路徑相關物理量Llv、L2v分別為：
式中m為體素V內(nèi)的散射次數(shù)，毛是體素V內(nèi)第i次散射的單位散射方向矢量，P1是各 向同性方向函數(shù)，？&是各向異性方向函數(shù)，L為體素V內(nèi)光子的第j次散射路徑長度。
3. 根據(jù)權利要求1所述的基于dfMC模型的熒光擴散光學斷層圖像重建方法，其特征在 于：步驟（2)中GPU集群加速整個白蒙特卡羅過程具體如下： 2. 1在集群上，依據(jù)光源數(shù)目及集群節(jié)點數(shù)目分配各個節(jié)點上進行MC模擬的次數(shù)，并 把相對應的光源位置和方向信息傳輸?shù)焦?jié)點上； 2. 2在每個節(jié)點上，根據(jù)要模擬的光子數(shù)目分配內(nèi)存，可以將組織模型信息分配到GPU內(nèi)的紋理寄存器，光學參數(shù)信息分配到常數(shù)寄存器來提高讀取速度，并根據(jù)GPU的技術信 息確定最優(yōu)的CUDA運行中并行塊和線程的數(shù)目； 2. 3在每個線程上，模擬一個光子包，并跟蹤光子發(fā)射后在組織內(nèi)散射及吸收的過程， 記錄光子在每個體素內(nèi)的相關路徑信息及物理量，光子溢出組織后終止光子模擬。
4. 根據(jù)權利要求1所述的基于dfMC模型的熒光擴散光學斷層圖像重建方法，其特征在 于步驟（3)，計算熒光權重W?(p)計算公式如下：
式中：P為位置及單位方向的六維矢量，纟為單位方向矢量，k為光子在組織內(nèi)的路徑所 經(jīng)過的體素數(shù)目，m為光子在組織內(nèi)的散射次數(shù)，Uafv是體素V內(nèi)熒光系數(shù)，/C是體素V內(nèi) 的激發(fā)光散射系數(shù)，n是量子產(chǎn)率，Wrai為激發(fā)光子權重，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ù)權利要求1所述的基于dfMC模型的熒光擴散光學斷層圖像重建方法，其特征在 于步驟（4)中，熒光雅可比矩陣的計算公式如下：
式中：P為位置及單位方向的六維矢量，丨為單位方向矢量，m為光子在組織內(nèi)的散 射次數(shù)，Uafv是體素V內(nèi)熒光系數(shù)，/C是體素V內(nèi)的激發(fā)光散射系數(shù)，n是量子產(chǎn)率， Wrai為激發(fā)光子權重，gv?為體素V內(nèi)熒光的各向異性因子，g/x為體素V內(nèi)激發(fā)光的各 向異性因子，yav?為體素V內(nèi)熒光的吸收系數(shù)，yJx為體素V內(nèi)激發(fā)光的吸收系數(shù)，且 廠(凡-凡，/〇 =哪(AHa)，j為節(jié)點i上的光源編號，d為探測器編號，V為需要重建的模 型體素編號，m為光子在組織內(nèi)的散射次數(shù)，J為熒光雅可比矩陣，心為對應體素V、光源j、 探測器d的熒光雅可比矩陣元。
6. 根據(jù)權利要求1所述的基于dfMC模型的快速熒光擴散光學斷層圖像重建方法，其特 征在于步驟（5)整個熒光重建過程均通過GPU集群加速，各個節(jié)點上計算得到的熒光雅可 比矩陣塊，計算公式如下：
Di= [D y Dlt;…；D yDsJt 式中：s為節(jié)點i上的光源數(shù)目，t為探測器總數(shù)，n為需要重建的模型體素數(shù)目，DiS節(jié)點i上所有光源下的探測器探測到的熒光強度分布與dfMC模擬得到的熒光分布的差值， Dst為節(jié)點i上第s光源下的探測器中第t個探測單元探測到的熒光強度分布與相對應的 dfMC模擬得到的熒光分布的差值。
7. 根據(jù)權利要求6所述的基于dfMC模型的快速熒光擴散光學斷層圖像重建方法，其特 征在于步驟（5)，中間量///X,JiDi被傳輸?shù)街鞴?jié)上再進行共軛梯度迭代求解熒光分布，這 些中間量比熒光雅可比矩陣塊占用的存貯空間少，迭代計算公式如下：
式中：i為節(jié)點編號，N為總節(jié)點數(shù)目，X為熒光團吸收系數(shù)分布矩陣Ji1為矩陣J,的轉(zhuǎn) 置矩陣，f⑴為目標函數(shù)。
【專利摘要】本發(fā)明公開了一種基于解耦合熒光蒙特卡羅(dfMC)模型的熒光擴散光學斷層圖像重建方法，屬于生物醫(yī)學工程技術領域，包括：確定探測區(qū)域，在探測區(qū)域選取若干掃描點，在探測器上獲取相對應的熒光光強分布；建立描述組織光學參數(shù)空間結構的三維數(shù)字模型，根據(jù)光源掃描的位置和方向進行激發(fā)光子的前向白蒙特卡羅(white Monte Carlo)模擬，跟蹤激發(fā)光子并記錄到達探測器的光子沿路徑上相應的物理量；依據(jù)dfMC方法計算熒光光子權重，并計算出熒光雅可比矩陣；最后通過GPU集群迭代重建計算出組織內(nèi)的熒光團位置及其吸收系數(shù)；本發(fā)明的優(yōu)點是，高精度的dfMC模型，基于GPU集群加速的重建迭代過程，為三維熒光斷層成像系統(tǒng)提供了一種精確快速的重建方法。
【IPC分類】G06T17-00, A61B5-00
【公開號】CN104706320
【申請?zhí)枴緾N201510050253
【發(fā)明人】駱清銘, 鄧勇, 羅召洋
【申請人】華中科技大學
【公開日】2015年6月17日
【申請日】2015年1月30日