CN103356170B - 用于带异质体组织光学参数重建的快速蒙特卡罗成像方法 - Google Patents
用于带异质体组织光学参数重建的快速蒙特卡罗成像方法 Download PDFInfo
- Publication number
- CN103356170B CN103356170B CN201310198484.0A CN201310198484A CN103356170B CN 103356170 B CN103356170 B CN 103356170B CN 201310198484 A CN201310198484 A CN 201310198484A CN 103356170 B CN103356170 B CN 103356170B
- Authority
- CN
- China
- Prior art keywords
- volume elements
- heteroplasmon
- prime
- light source
- photon
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Landscapes
- Investigating Or Analysing Materials By Optical Means (AREA)
Abstract
本发明属于生物医学工程技术领域,涉及一种用于带异质体组织光学参数重建的快速蒙特卡罗成像算法,包括:确定探测区域,并对其进行轴向剖分,选取剖分网格线的交点为扫描点,得到不同扫描位置下的漫反射光强分布;做光源位于原点、探测器对称分布在其两侧时的MC模拟,得到光子在每个体元内的平均碰撞次数、平均行走的路程长以及所有探测器各自接受到的光子权值;第三步:还原出源探位于所有扫描点时每个体元内的碰撞次数及行走的路程长,计算输出光子的权重,并求出Jacobi矩阵;第四步:计算出异质体位置及其光学参数。本发明的有益效果为:能够对含有异质体组织的光学参数进行重建,同时优化pMC算法,节约存储空间,并缩短计算时间。
Description
技术领域
本发明属于生物医学工程技术领域,涉及一种适合于介观、中等空间分辨率的漫射光功能层析成像算法。
背景技术
癌症不但可导致组织结构和细胞形态的变化,而且可以引起功能代谢活动的改变,发展功能影像技术可为早期宫颈癌诊断提供更为客观、特异性的依据。目前,针对薄层上皮组织如宫颈及皮肤,主要采用介观、中等空间分辨率的漫射光功能层析成像方法。然而在探测距离小于1mm时,无论是漫射方程还是更加复杂的混合漫射-P3近似方程在描述光在组织中的传播上都有很大的误差。鉴于蒙特卡洛模拟(MC)[1]在描述光在组织中传播时的广泛适用性,本发明将研究基于蒙特卡罗(MC)模拟[1]的图像重建方法。
MC模拟已经被证明是描述光在任意形状生物组织中传输的有效模型,但其结果的准确性依赖于统计的精度,为了获得准确模型需要追踪大量光子,因此,蒙特卡罗模拟描述光子在组织体中传输的唯一限制是其计算量。不仅如此,光学参数的重构过程多是通过求解非线性问题得以实现,其要求是当选定一个迭代算法之后,正向模型应该给出多个派生信息,由这些信息引导迭代算法的收敛方向和速度,每迭代一次都需要进行一次正向模型运算。因此,正向模型的计算速度问题极大地影响光学参数的重构。
微扰蒙特卡罗(pMC)是一种快速MC方法,利用已知光子在目标区域内碰撞的次数(k)和行走的路程长(s),得到更新光学参数后输出光子的权值,与其他MC模拟方法类似,pMC没有扩散(漫射)方程(diffusion equation,DE)对源探距离和光学参数的严格要求,使其应用领域较为宽广。但由于需要运行一次完整的MC模拟来获得光子在组织中的行走路径,从而求出k和s,使得目前的pMC算法仅限于对均匀介质的光学参数进行重建,或者在异质体位置已知的情况下对异质体光学参数进行重建,而这与现实需求相差甚远。另外,为了保存光子的行走路径,不仅占用大量存储空间,而且严重增加了重建过程的计算时间。
参考文献:
[1]Wang L,Jacques S L,Zheng L.MCML—Monte Carlo modeling of light transport in multi-layered tissues[J].Computer methods and programs in biomedicine,1995,47(2):131-146.
[2]Schweiger M,Arridge S R,I.Gauss–Newton method for image reconstruction in diffuse optical tomography[J].Physics in medicine and biology,2005,50(10):2365.
发明内容
本发明的目的是,克服现有技术的上述不足,提供一种可以对含有异质体组织的光学参数重建的快速层析成像方法。本发明的技术方案如下:
一种用于带异质体组织光学参数重建的快速蒙特卡罗成像方法,包括下列步骤:
第一步:确定探测区域,并对其进行轴向剖分,选取剖分网格线的交点为扫描点,得到不同扫描位置 下的漫反射光强分布;
第二步:做光源位于原点、探测器对称分布在其两侧时的MC模拟,得到光子在每个体元内的平均碰撞次数、平均行走的路程长以及所有探测器各自接受到的光子权值w;
第三步:还原出源探位于所有扫描点时每个体元内的碰撞次数及行走的路程长,计算输出光子的权重,并求出Jacobi矩阵,方法如下:
(1)对MC模拟结果中每一个不全为零的体元按照镜像或者平移关系还原出所有等效的源探-体元关系对;
(2)针对每一个确定的源探-体元关系对,计算相应的衰减率,设某个源探-体元关系对的体元为i,其相应的衰减率ξi的计算公式为:
式中:——分别为组成异质体的第i个体元的散射系数和总衰减系数;
μsi、μti——分别为组织体均匀部分中第i个体元的散射系数和总衰减系数;
——分别是探测器位于光源左侧时目标区域中第i个体元内的碰撞次数和行走的路径长;
——探测器位于光源右侧时目标区域中第i个体元内的碰撞次数和行走的路径长;
(3)计算每个探测器位置输出光子的权重:
式中:ξi——异质体区域中分别位于光源左侧时第i个体元的衰减率;
ξj——异质体区域中分别位于光源右侧时第j个体元的衰减率;
m、n——分别为异质体区域中位于光源两侧的体元总数;
w、w′——分别为加入异质体前后输出光子权值;
(4)通过对散射系数和总衰减系数分别求偏导,得到Jacobi矩阵的系数;
第四步:根据各个输出光子的权重和Jacobi矩阵,计算异质体位置及其光学参数。
作为优选实施方式,第四步中,采用Newton-Raphson法计算出异质体位置及其光学参数。
本发明的有益效果为:能够对对含有异质体组织的光学参数进行重建,特别适合当探测距离小于1mm且存在多种光源分布的情况。同时优化了pMC算法,节约了存储空间,从而缩短计算时间。
附图说明
图1是剖分后的平行平板模型,图中黑点所示为扫描点
图2a是光源点位于原点;图2b是光源点位于目标区域中的某一点;图2c是上述两种情况的合并
图3是体元关于y轴做镜像
图4是体元关于x轴做镜像
其中:1为扫描点,2为光源,3为探测器,4为目标区域,5为体元
具体实施方式
下面结合附图和实施例对本发明进行说明
一、利用互易法提高碰撞次数k及行走的路程长s的计算效率
对于含肿瘤的生物组织体,可以采用光学参数相同的带异质体平板模型(如图(1)所示)近似代替,设光源及探测器位置如图中所示。为了获取模型中光学参数的分布,从而得到异质体位置及其光学参数,需对平板模型进行轴向剖分(见图(1)),划分成一系列匀质体元。选取剖分网格线的交点为扫描点,使光源沿扫描点按图中所示箭头方向逐行扫描,并保持探测器与光源的相对距离,得到不同扫描位置下的漫反射光强分布。将检测到的光强分布与光子在组织体中传播的正向模型反复比较,指导修正体元的光学参数。基于pMC技术建立正向模型时必须得到光子在每个体元内的碰撞次数(k)和行走的路程长(s),虽然通过MC模拟可以计算出特定光源位置下每个体元内的k及s,但是光源位置的不同需要参数配置不同的MC程序,因此如何避免对多个光源位置进行MC模拟成为关键。
考虑图2a中的情况,当光源移至图2b位置时,需要重新计算光源位置然后进行MC模拟。如果将图2b中的光源平移至原点,利用互易法,即体元及源探同时平移或镜像变换后光子传播路径不变,同时考虑xy平面剖分网格与扫描网格重合,因此将探测器及体元通过平移相同的网格数后与图2a合并成图2c的形式,以此类推,对于其他光源不在原点的情况都可以做上述平移,然后通过一次MC模拟就可以得到所有光源位置下每个体元内的k及s,但是此时目标区域为原来的4倍,分布在原点周围的四个象限内。进一步设想,将位于第二、四象限的体元分别做关于y轴、x轴的镜像,对第三象限的体元先做x轴,再坐y轴镜像,同时变换探测器位置(详见具体实施方式),则所有体元回到目标区域内,这样只要对目标区域做光源在原点时的MC模拟即可得到结果。
做一次完整的光源在原点,探测器对称分布在源点左右两边时的MC模拟,分别得到目标区域中第i个体元内的碰撞次数和行走的路程长:及(l、r分别代表探测器位于光源左侧和右侧)。在建立正向模型时,为了得到所有光源位置下每个体元内的k及s,应首先还原所有等效的源探-体元关系对:
(1)当探测器位于原点右侧,对所有s非零的体元按照图4关于x轴做镜像,得到的体元位于第四象限,此时体元在y轴上的位置应该是黑色箭头所标处;当探测器位于原点左侧,需要考虑两种情况:首先关于y轴对所有s非零的体元按照图3做镜像,得到的体元位于第二象限,其x轴坐标应是黑色箭头所标处,为了保持相对关系,探测器也需要同时关于y轴做镜像。其次将位于第二象限的体元再关于x轴做镜像体元,得到位于第三象限的体元。此时,得到了源探相对于体元的四种位置关系(第一象限不需要做镜像)。
(2)在对目标区域逐点探测过程中,源探相对体元的位置关系分为上述四种情况。因此只需将上述四 个象限内的所有源探-体元关系对保持相对距离平移至第一象限,然后在目标区域范围内按扫描步长逐点平移,由于体元内的k、s已知,最终得到每个源探位置下所有体元内的k、s。然后按照下式计算每个探测器位置输出光子的权重:
式中:ξj与ξi求法相同,w为加入异质体前输出光子权值,
m、n分别为异质体区域中位于光源两侧的体元总数,分别为组成异质体的第i个体元的散射系数和总衰减系数,μsi、μti分别为组织体均匀部分中第i个体元的散射系数和总衰减系数,当MC模拟的光子数增加到一定数量,关于光源对称的探测器得到同样多的光子权值w。
当存在多个距离光源不同的探测器时,由于源探距离决定探测深度,因此可以增加对深度方向(即z轴方向)的剖分,从而得到组织体的三维模型。由于剖分的疏密与图像分辨率紧密相关,为了提高横向分辨率,可以对上述xy平面剖分后的每个网格继续进行n次轴向剖分,使体元个数为之前的4n倍,由于此时探测点数量及位置不变,所以源探移动一次体元需要移动2n次。
对(2)式中衰减率ξ中的分别求偏导得到Jacobi矩阵J(p):
式中:代表在剖分网格所有节点上光学参数,ξs(s=1,2,…,S)为表面有限个不同的激励源位置,为与源相对应的有限个表面检测点。
最后,利用得到的Jacobi矩阵及每个探测器位置输出光子的权重,采用Newton-Raphson[2]法计算出异质体位置及其光学参数。
二、优化pMC算法
按照传统的pMC方法,如果探测器接受到n个光子,即生成了n条光子轨迹,根据每条轨迹可以计算出目标区域所有体元内的k、s,然后套用n次(1)式计算光子权重,最后再累加,因此这其中需要存储n条轨迹对应的所有体元内的k、s。
设kin为通过光子的行走路径计算探测器接受到第n个光子时第i个体元的碰撞的次数,按照下式计算该体元内的平均碰撞次数
平均行走的路程长的计算方法与相同。对n个光子的权重进行累加得到:w=w1+w2+...+wn,由于随着发射光子数的增加,左右两侧的探测器接受到的光子权值逐渐趋于相等,因此任选一侧的探测器计算w即可。
这样MC模拟的结果中只需包含目标区域所有体元的k、s及光子权重w,远比存储n次发射的所有光子行走路径节约空间。另外,由于光子集中分布在光源附近,而对于s等于0(当s为0时k一定为0)的大量体元可以舍弃。
Claims (1)
1.一种用于带异质体组织光学参数重建的快速蒙特卡罗成像方法,包括下列步骤:
第一步:确定探测区域,并对其进行轴向剖分,选取剖分网格线的交点为扫描点,得到不同扫描位置下的漫反射光强分布;
第二步:做光源位于原点、探测器对称分布在其两侧时的蒙特卡罗模拟,得到光子在每个体元内的平均碰撞次数平均行走的路程长以及所有探测器各自接受到的光子权值w;
第三步:还原出源探位于所有扫描点时每个体元内的碰撞次数及行走的路程长,计算输出光子的权重,并求出Jacobi矩阵,方法如下:
(1)对蒙特卡罗模拟结果中每一个不全为零的体元按照镜像或者平移关系还原出所有等效的源探-体元关系对;
(2)针对每一个确定的源探-体元关系对,计算相应的衰减率,设某个源探-体元关系对的体元为i,其相应的衰减率ξi的计算公式为:
式中:μ′si、μ′ti——分别为组成异质体的第i个体元的散射系数和总衰减系数;
μsi、μti——分别为组织体均匀部分中第i个体元的散射系数和总衰减系数;
——分别是探测器位于光源左侧时目标区域中第i个体元内的碰撞次数和行走的路径长;
——探测器位于光源右侧时目标区域中第i个体元内的碰撞次数和行走的路径长;
(3)计算每个探测器位置输出光子的权重:
式中:ξi——异质体区域中分别位于光源左侧时第i个体元的衰减率;
ξj——异质体区域中分别位于光源右侧时第j个体元的衰减率;
m、n——分别为异质体区域中位于光源两侧的体元总数;
w、w′——分别为加入异质体前后输出光子权值;
(4)通过对散射系数μ′si、和总衰减系数μ′ti分别求偏导,得到Jacobi矩阵的系数;
第四步:根据各个输出光子的权重和Jacobi矩阵,采用Newton-Raphson法计算出异质体位置及其光学参数。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201310198484.0A CN103356170B (zh) | 2013-05-24 | 2013-05-24 | 用于带异质体组织光学参数重建的快速蒙特卡罗成像方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201310198484.0A CN103356170B (zh) | 2013-05-24 | 2013-05-24 | 用于带异质体组织光学参数重建的快速蒙特卡罗成像方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN103356170A CN103356170A (zh) | 2013-10-23 |
CN103356170B true CN103356170B (zh) | 2015-02-18 |
Family
ID=49359079
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201310198484.0A Expired - Fee Related CN103356170B (zh) | 2013-05-24 | 2013-05-24 | 用于带异质体组织光学参数重建的快速蒙特卡罗成像方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN103356170B (zh) |
Families Citing this family (15)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104706320B (zh) * | 2014-11-14 | 2017-02-22 | 华中科技大学 | 一种基于dfMC模型的荧光扩散光学断层图像重建方法 |
CN104605823B (zh) * | 2014-11-14 | 2017-01-18 | 华中科技大学 | 一种基于pfMC模型的荧光扩散光学断层图像重建方法 |
CN104679946B (zh) * | 2014-11-14 | 2017-12-05 | 华中科技大学 | 一种基于体素的微扰荧光蒙特卡罗模拟方法 |
CN104699993B (zh) * | 2015-03-31 | 2017-10-20 | 西安交通大学 | 辐射能传播蒙特卡洛算法中的局部网格及光子加密方法 |
CN105816151B (zh) * | 2016-03-10 | 2018-08-07 | 天津大学 | 一种基于空间频域测量的均匀组织体光学参数重建方法 |
CN106290199B (zh) * | 2016-09-13 | 2019-02-19 | 天津大学 | 基于稳态辐射率间质测量装置的多波长光学参数反构方法 |
CN110072435B (zh) * | 2016-12-08 | 2022-07-19 | 皇家飞利浦有限公司 | 表面组织跟踪 |
CN107677644B (zh) * | 2017-08-23 | 2019-11-01 | 北京大学 | 一种多层组织体光学参数的检测系统及其检测方法 |
CN107993722B (zh) * | 2017-12-05 | 2021-12-07 | 天津大学 | 一种针对点阵光源下光子在媒质中分布的快速提取方法 |
CN109991194B (zh) * | 2017-12-29 | 2022-04-26 | 天津先阳科技发展有限公司 | 漫反射光谱中抑制温度干扰的方法、光谱分析方法及装置 |
CN108846790A (zh) * | 2018-06-15 | 2018-11-20 | 华中科技大学 | 一种加速图像重建的方法 |
CN108680537A (zh) * | 2018-07-11 | 2018-10-19 | 天津工业大学 | 一种用于检测近似均匀的非透明介质中异质体的快速定位方法 |
CN109342367B (zh) * | 2018-09-30 | 2020-04-10 | 华中科技大学 | 一种基于控制蒙特卡罗方法的扩散光学成像方法及系统 |
CN110200593B (zh) * | 2019-06-19 | 2020-08-07 | 天津工业大学 | 一种头带式异质体检测设备及头部异质体检测方法 |
CN110911007B (zh) * | 2019-12-29 | 2023-07-25 | 杭州科洛码光电科技有限公司 | 基于成像光谱仪的生物组织光学参数重构方法 |
Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US6850656B1 (en) * | 1998-10-07 | 2005-02-01 | Ecole Polytechnique Federale De Lausanne | Method and apparatus for measuring locally and superficially the scattering and absorption properties of turbid media |
CN101286187A (zh) * | 2008-06-10 | 2008-10-15 | 华中科技大学 | 光在生物组织中传输特性的定量蒙特卡罗模拟方法 |
CN101313847A (zh) * | 2008-07-01 | 2008-12-03 | 北京师范大学 | 对人体皮肤病变组织进行无损光学常数成像的装置和方法 |
CN101513343A (zh) * | 2009-01-13 | 2009-08-26 | 华中科技大学 | 获取稳态/瞬态光漫射特性的分析系统及方法 |
CN101856219A (zh) * | 2010-05-13 | 2010-10-13 | 天津大学 | 基于频域近红外光测量的光学参数重构方法 |
CN101966078A (zh) * | 2010-11-09 | 2011-02-09 | 天津析像光电科技有限公司 | 一种近红外漫射光频域信息获取方法 |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2008103486A1 (en) * | 2007-02-23 | 2008-08-28 | Duke University | Scaling method for fast monte carlo simulation of diffuse reflectance spectra |
-
2013
- 2013-05-24 CN CN201310198484.0A patent/CN103356170B/zh not_active Expired - Fee Related
Patent Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US6850656B1 (en) * | 1998-10-07 | 2005-02-01 | Ecole Polytechnique Federale De Lausanne | Method and apparatus for measuring locally and superficially the scattering and absorption properties of turbid media |
CN101286187A (zh) * | 2008-06-10 | 2008-10-15 | 华中科技大学 | 光在生物组织中传输特性的定量蒙特卡罗模拟方法 |
CN101313847A (zh) * | 2008-07-01 | 2008-12-03 | 北京师范大学 | 对人体皮肤病变组织进行无损光学常数成像的装置和方法 |
CN101513343A (zh) * | 2009-01-13 | 2009-08-26 | 华中科技大学 | 获取稳态/瞬态光漫射特性的分析系统及方法 |
CN101856219A (zh) * | 2010-05-13 | 2010-10-13 | 天津大学 | 基于频域近红外光测量的光学参数重构方法 |
CN101966078A (zh) * | 2010-11-09 | 2011-02-09 | 天津析像光电科技有限公司 | 一种近红外漫射光频域信息获取方法 |
Non-Patent Citations (6)
Title |
---|
《Gauss-Newton method for image reconstruction in diffuse optical tomography》;Martin Schweiger et al;《Physics in medicine and Biology》;20051231;第50卷(第20期);第2365-2386页 * |
《MCML-Monte Carlo modeling of light transport in multi-layered tissues》;Wang L etal;《Computer Methods and Programs in Biomedicine》;19951231;第47卷;第131-146页 * |
《Measurement of optical transport properties of normal and malignant human breast tissue》;N.Ghosh etal;《Applied Optics》;20010131;第40卷(第1期);第176-184页 * |
《Perturbation Monte Carlo methods to solve inverse photon migration problems in heterogenous tissues》;Carole K,etal;《Optics Letters》;20010901;第26卷(第17期);第1335-1337页 * |
《基于蒙特卡罗传输模型的组织光学参数重构方法的研究》;关堂兵;《中国优秀硕士学位论文全文数据库 医药卫生科技辑 》;20090415(第4期);E080-61 * |
《生物组织光学断层成像与光学参数提取方法研究》;林林;《中国博士学位论文全文数据库 医药卫生科技辑》;20110915(第9期);E080-10 * |
Also Published As
Publication number | Publication date |
---|---|
CN103356170A (zh) | 2013-10-23 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN103356170B (zh) | 用于带异质体组织光学参数重建的快速蒙特卡罗成像方法 | |
CN103713288B (zh) | 基于迭代最小化稀疏贝叶斯重构线阵sar成像方法 | |
CN102779350B (zh) | 一种锥束ct迭代重建算法投影矩阵构建方法 | |
CN104933737B (zh) | 一种基于共轭梯度法的电离层层析成像混合反演方法 | |
CN101856219B (zh) | 基于频域近红外光测量的光学参数重构方法 | |
Grote et al. | Local nonreflecting boundary condition for time-dependent multiple scattering | |
CN103454695B (zh) | 一种gps电离层tec层析方法 | |
CN102034266B (zh) | 激发荧光断层成像的快速稀疏重建方法和设备 | |
Da Silva et al. | Sub-second pencil beam dose calculation on GPU for adaptive proton therapy | |
CN102988026A (zh) | 基于乘子法的自发荧光断层成像重建方法 | |
CN106097441A (zh) | 基于l1范数与tv范数的复合正则化生物发光断层成像重建方法 | |
Mahariq et al. | On the Attenuation of the Perfectly Matched Layer in Electromagnetic Scattering Problems with the Spectral Element Method. | |
CN104586366B (zh) | 一种激发荧光断层成像重建方法 | |
CN105559750A (zh) | 组织结构引导的复合正则化生物发光断层成像重建方法 | |
Lyubimov et al. | Application of the photon average trajectories method to real-time reconstruction of tissue inhomogeneities in diffuse optical tomography of strongly scattering media | |
Erdoğan et al. | Incorporating topography into 2D resistivity modeling using finite-element and finite-difference approaches | |
Carpio et al. | Noninvasive imaging of three-dimensional micro and nanostructures by topological methods | |
CN104851080A (zh) | 一种基于tv的三维pet图像重建方法 | |
CN101980304A (zh) | 一种三维数字体积图像变形测量方法 | |
Krah et al. | Polynomial modelling of proton trajectories in homogeneous media for fast most likely path estimation and trajectory simulation | |
CN103393410A (zh) | 一种基于交替迭代运算的荧光分子断层成像重建方法 | |
Tadeu et al. | 3D acoustic wave simulation using BEM formulations: Closed form integration of singular and hypersingular integrals | |
CN104706320B (zh) | 一种基于dfMC模型的荧光扩散光学断层图像重建方法 | |
Huang et al. | Evaluation of minimum zone sphericity by combining single-space contraction strategy with multi-directional adaptive search algorithm | |
Afanasiev et al. | Phase fluctuations of radio waves experiencing total reflection from a randomly inhomogeneous plasma layer |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
C14 | Grant of patent or utility model | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20150218 Termination date: 20210524 |
|
CF01 | Termination of patent right due to non-payment of annual fee |