用于带异质体组织光学参数重建的快速蒙特卡罗成像方法

摘要

本发明属于生物医学工程技术领域，涉及一种用于带异质体组织光学参数重建的快速蒙特卡罗成像算法，包括：确定探测区域，并对其进行轴向剖分，选取剖分网格线的交点为扫描点，得到不同扫描位置下的漫反射光强分布；做光源位于原点、探测器对称分布在其两侧时的MC模拟，得到光子在每个体元内的平均碰撞次数、平均行走的路程长以及所有探测器各自接受到的光子权值；第三步：还原出源探位于所有扫描点时每个体元内的碰撞次数及行走的路程长，计算输出光子的权重，并求出Jacobi矩阵；第四步：计算出异质体位置及其光学参数。本发明的有益效果为：能够对含有异质体组织的光学参数进行重建，同时优化pMC算法，节约存储空间，并缩短计算时间。

说明书

技术领域

本发明属于生物医学工程技术领域，涉及一种适合于介观、中等空间分辨率的漫射光功能层析成像算法。 

背景技术

癌症不但可导致组织结构和细胞形态的变化，而且可以引起功能代谢活动的改变，发展功能影像技术可为早期宫颈癌诊断提供更为客观、特异性的依据。目前，针对薄层上皮组织如宫颈及皮肤，主要采用介观、中等空间分辨率的漫射光功能层析成像方法。然而在探测距离小于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的计算公式为： 

${\xi }_i={\left(\frac{{\mu }_{\mathrm{si}}^{\prime }/{\mu }_{\mathrm{ti}}^{\prime }}{{\mu }_{\mathrm{si}}/{\mu }_{\mathrm{ti}}}\right)}^{k_i^l}{\left(\frac{{\mu }_{\mathrm{ti}}^{\prime }}{{\mu }_{\mathrm{ti}}}\right)}^{k_i^l}\exp [-({\mu }_{\mathrm{ti}}^{\prime }-{\mu }_{\mathrm{ti}})s_i^l]$

式中：——分别为组成异质体的第i个体元的散射系数和总衰减系数； 

μ_si、μ_ti——分别为组织体均匀部分中第i个体元的散射系数和总衰减系数； 

——分别是探测器位于光源左侧时目标区域中第i个体元内的碰撞次数和行走的路径长； 

——探测器位于光源右侧时目标区域中第i个体元内的碰撞次数和行走的路径长； 

（3）计算每个探测器位置输出光子的权重： 

$w^{\prime }=w·\underset{i=1}{\overset{i=m}{\Pi }}{\xi }_i·\underset{j=1}{\overset{j=n}{\Pi }}{\xi }_j$

式中：ξ_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。然后按照下式计算每个探测器位置输出光子的权重： 

$w^{\prime }=w·\underset{i=1}{\overset{i=m}{\Pi }}{\xi }_i·\underset{j=1}{\overset{j=n}{\Pi }}{\xi }_j---\left(1\right)$

式中：ξ_j与ξ_i求法相同，w为加入异质体前输出光子权值， 

m、n分别为异质体区域中位于光源两侧的体元总数，分别为组成异质体的第i个体元的散射系数和总衰减系数，μ_si、μ_ti分别为组织体均匀部分中第i个体元的散射系数和总衰减系数，当MC模拟的光子数增加到一定数量，关于光源对称的探测器得到同样多的光子权值w。 

当存在多个距离光源不同的探测器时，由于源探距离决定探测深度，因此可以增加对深度方向（即z轴方向）的剖分，从而得到组织体的三维模型。由于剖分的疏密与图像分辨率紧密相关，为了提高横向分辨率，可以对上述xy平面剖分后的每个网格继续进行n次轴向剖分，使体元个数为之前的4n倍，由于此时探测点数量及位置不变，所以源探移动一次体元需要移动2n次。 

对（2）式中衰减率ξ中的分别求偏导得到Jacobi矩阵J(p)： 

$J\left(p\right)=\left(\begin{array}{cccccc}\frac{\partial w^{\prime }({\xi }_1^{\left(1\right)},{\zeta }_1)}{\partial {\mu }_{s1}^{\prime }},& ...,& \frac{\partial w^{\prime }({\xi }_1^{\left(1\right)},{\zeta }_1)}{\partial {\mu }_{\mathrm{sN}}^{\prime }},& \frac{\partial w^{\prime }({\xi }_1^{\left(1\right)},{\zeta }_1)}{\partial {\mu }_{t1}^{\prime }},& ...,& \frac{\partial w^{\prime }({\xi }_1^{\left(1\right)},{\zeta }_1)}{\partial {\mu }_{\mathrm{tN}}^{\prime }}\\ .& & & & & .\\ .& & & & & .\\ .& & & & & .\\ \frac{\partial w^{\prime }({\xi }_d^{\left(s\right)},{\zeta }_s)}{\partial {\mu }_{s1}^{\prime }},& ...,& \frac{\partial w^{\prime }({\xi }_d^{\left(s\right)},{\zeta }_s)}{\partial {\mu }_{\mathrm{sN}}^{\prime }},& \frac{\partial w^{\prime }({\xi }_d^{\left(s\right)},{\zeta }_s)}{\partial {\mu }_{t1}^{\prime }},& ...,& \frac{\partial w^{\prime }({\xi }_d^{\left(s\right)},{\zeta }_s)}{\partial {\mu }_{\mathrm{tN}}^{\prime }}\\ .& & & & & .\\ .& & & & & .\\ .& & & & & .\\ \frac{\partial w^{\prime }({\xi }_D^{\left(S\right)},{\zeta }_S)}{\partial {\mu }_{s1}^{\prime }},& ...,& \frac{\partial w^{\prime }({\xi }_D^{\left(S\right)},{\zeta }_S)}{\partial {\mu }_{\mathrm{sN}}^{\prime }},& \frac{\partial w^{\prime }({\xi }_D^{\left(S\right)},{\zeta }_S)}{\partial {\mu }_{t1}^{\prime }},& ...,& \frac{\partial w^{\prime }({\xi }_D^{\left(S\right)},{\zeta }_S)}{\partial {\mu }_{\mathrm{tN}}^{\prime }}\end{array}\right)$

式中：代表在剖分网格所有节点上光学参数，ξ_s(s=1，2，…，S)为表面有限个不同的激励源位置，为与源相对应的有限个表面检测点。 

最后，利用得到的Jacobi矩阵及每个探测器位置输出光子的权重，采用Newton-Raphson^[2]法计算出异质体位置及其光学参数。 

二、优化pMC算法 

按照传统的pMC方法，如果探测器接受到n个光子，即生成了n条光子轨迹，根据每条轨迹可以计算出目标区域所有体元内的k、s，然后套用n次（1）式计算光子权重，最后再累加，因此这其中需要存储n条轨迹对应的所有体元内的k、s。 

设k_in为通过光子的行走路径计算探测器接受到第n个光子时第i个体元的碰撞的次数，按照下式计算该体元内的平均碰撞次数

${\bar{k}}_i=\frac{k_{i1}{w}_1+k_{i2}{w}_2+···+k_{\mathrm{in}}{w}_n}{w_1+w_2+···+w_n}---\left(1\right)$

平均行走的路程长的计算方法与相同。对n个光子的权重进行累加得到：w＝w₁+w₂+...+w_n，由于随着发射光子数的增加，左右两侧的探测器接受到的光子权值逐渐趋于相等，因此任选一侧的探测器计算w即可。 

这样MC模拟的结果中只需包含目标区域所有体元的k、s及光子权重w，远比存储n次发射的所有光子行走路径节约空间。另外，由于光子集中分布在光源附近，而对于s等于0（当s为0时k一定为0）的大量体元可以舍弃。