CN103230258B - 基于有限体积法的扩散光学断层成像正向问题处理方法 - Google Patents

Abstract

本发明涉及有限体积法用于处理扩散光学断层成像正向过程的算法，该方法将待求解生物组织区域的光学特征参数离散成多个相互联系的控制体单元，通过每一个控制体单元上物理量的守恒特性及相互关联特性，从而推广到整个组织体区域，然后结合特定的描述光在组织中传播的理论和有限体积法得到边界处光子密度值。本发明在实施过程中将求解域内的计算转化到求解域的边界，这不仅减少了计算量，同时也具有较高的精度，最终不仅使得处理扩散光学成像正向过程的效率得到很大的提升，同时为扩散光学图像重建的逆问题奠定基础。

Description

基于有限体积法的扩散光学断层成像正向问题处理方法

技术领域

本发明属于生物医学成像技术应用和生物医学工程领域，涉及一种基于有限体积法的扩散光学断层成像正向问题处理方法。

背景技术

传统成像手段已经取得了极大的进展，比如CT、核磁共振等。但是伴随着不错的成像效果所带来的风险就是核辐射对人体带来的伤害，也是发展为癌症的诱发因素。光学分子成像手段相对于这些成像方式来说，具有非入侵(non-invasive)、无辐射伤害、灵敏度更高的特点。并且光学分子成像具有更宽的可选荧光探针谱，通过探针结合组织中特定分子靶标可以实现在体动态成像，被认为是未来分子成像技术的突破点。

而以扩散理论为基础的扩散光学断层成像技术(DOT)，它利用穿过组织的漫射光信息可以重建出组织体在某一个断层面甚至三维空间的光学或生理参数分布。这种技术不仅成本低，而且具有较深的成像深度，在乳腺肿瘤的诊断检测等众多领域有重要的地位.

然而现在扩散光断层成像面临的一个重大的瓶颈就是在空间分辨率和时间分辨率上的不足。虽然采用辐射传输方程可以提高精确度，但是，鉴于组织体的复杂性，通常情况下无法得到解析解，正是由于这个原因，数值法，比如有限元、有限差分法等才被广泛地应用于方程的求解，有限体积法也是这众多方法中的一种。

发明内容

本发明的内容在于提供了一种基于有限体积法的扩散光学断层成像正向问题处理方法，该方法可用于二维及三维的扩散光学断层成像正向过程的处理。

本发明提供了基于有限体积法的扩散光学断层成像正向问题处理方法，其步骤如下：

基于有限体积法的扩散光学断层成像正向问题处理方法，其特征在于包括以下步骤：

(1) 创建组织体模型：使用吸收系数、散射系数、折射率和各向异性因子四个参数描述该生物组织不同区域内的光学特性；

(2) 按照精度需求的不同将待求解的生物组织区域按照现有的原始剖分和对偶剖分技术进行剖分，从而将原来连续的生物组织区域转化为离散的三角单元和控制体单元；

(3) 在三角单元和控制体单元的节点处定义试探函数空间Φ_h作为该节点处光子密度的近似值，使用现有的稳态扩散方程作为光子在组织内的传输模型，采用现有的双线性方法，将稳态扩散方程在每一个控制体单元上进行积分；结合边界条件，利用 Green 公式 ，其中Φ表示该节点处的光子密度值，表示控制体单元，表示控制体单元的边界，表示对边界单位外法向的偏导数, 将对控制体单元面积的积分转化成对每个控制体单元上的边界积分的求和；对每一个控制体单元进行相同的稳态扩散方程积分处理，从而得到每一个控制体单元上的稳态扩散方程；

(4)将步骤(3)中得到的所有的控制体单元的稳态扩散方程进行组合，形成总区域上的控制体单元方程，得到与稳态扩散方程相对应的线性方程组；求解该线性方程组即得到生物组织体边界处各个节点处的光子密度值。

步骤(1)中建立组织体模型包括生物组织的背景区域和位于其中的异质体的具体位置、大小及各区域光学特征参数分布情况。

步骤(2)中对待求解生物组织区域进行原始剖分和对偶剖分，形成三角单元和控制体单元，具体步骤为：通过对待求解生物组织区域进行原始的三角剖分和重心对偶剖分，将整个生物组织区域离散为控制体单元的集合；将吸收系数、散射系数、折射率和各向异性因子对应离散到每一个控制体单元的节点上。

步骤(3)中将稳态扩散方程在每一个控制体单元上进行积分，具体步骤为：以试探函数空间Φ_h的取值作为各节点处光子密度的近似值；采用双线性方法，首先定义检验函数空间v_h为分片常数空间，定义如下：

$v_h(x,y)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ψ}}_{p_0}(x,y)v_{p_0}$

其中   ${\mathrm{\ψ}}_{p_0}\left(p\right)=\left\{\begin{array}{cc}1,& p\∈K_{p0}^{*}\\ 0,& p\∉K_{p0}^{*}\end{array}\right.$

其中表示基底，表示该节点处的基函数，p表示节点编号，表示以p₀为中心的控制体单元, N表示节点个数；

在每一个控制体上进行积分，从而得到每一个控制体单元上的稳态扩散方程：

$\underset{k_{p_0}^{*}}{\∫\∫}(-\▿\·D\▿+{\mathrm{\μ}}_a)\mathrm{\Φvdxdy}=\underset{k_{p_0}^{*}}{\∫\∫}\mathrm{qvdxdy}$

其中Φ为光子密度值，μ_a代表吸收系数，D代表扩散系数，表示以p₀为中心的控制体单元， q表示光源，v表示检验函数。

变分形式为：

$a(\mathrm{\Φ},v)=(q,v)$

其中a(Φ,v)表示双线性泛函，Φ为光子密度值，q表示光源，v表示检验函数；

使用Green 公式和Robin边界条件：

$\mathrm{\Φ}\left(r\right)+2D\·\frac{1+R_f}{1-R_f}\·\▿\mathrm{\Φ}\left(r\right)=0,\∀r\∈\∂\mathrm{\Ω}$

其中Φ为光子密度值，D代表扩散系数，R_f为扩散传输内反射系数，r表示节点向量；将对生物组织的区域面积的积分转化为对边界的曲线积分：

$\underset{{\∂k}_{p_0}^{*}}{\∫}-D\frac{\∂\mathrm{\Φ}}{\∂n}\mathrm{vds}+\underset{\∂\mathrm{\Ω}}{\∫}\frac{1}{2\mathrm{\ζ}}\mathrm{\Φvds}+\underset{k_{p_0}^{*}}{\∫\∫}{\mathrm{\μ}}_a\mathrm{\Φvdxdy}=\underset{k_{p_0}^{*}}{\∫\∫}\mathrm{qvdxdy}$ ，

其中其中Φ为光子密度值，D代表扩散系数，μ_a代表吸收系数， 表示以p₀为中心的控制体单元，表示以p₀为中心的控制体单元的边界，表示生物组织边界, 表示对边界单位外法向的偏导数，ζ是一个与边界折射率失配程度有关的常系数，q表示光源，v表示检验函数。

步骤(4)中将所有的控制体单元的稳态扩散方程进行组合形成总的控制体方程，得到与稳态扩散方程相对应的线性方程组；求解该线性方程组即得到组织体边界处各个节点处的光子密度值，具体步骤为：将步骤(3)中得到的每一个控制体单元的稳态扩散方程按照节点顺序进行组合，形成总的矩阵方程，得到与稳态扩散方程相对应的线性方程组，进而将整个生物组织区域的稳态扩散方程写作如下矩阵的形式：

$\mathrm{K\Φ}=(A+B+C)\mathrm{\Φ}=F$

矩阵K代表总的控制体单元刚度矩阵，Φ为光子密度值，A，B，C，F分别对应下式中的矩阵元素，如下所示：

$\begin{array}{c}A=\underset{k_{p_0}^{*}}{\∫\∫}{\mathrm{\μ}}_a\mathrm{\Φvdxdy}\\ B=\underset{{\∂k}_{p_0}^{*}}{\∫}-D\frac{\∂\mathrm{\Φ}}{\∂n}\mathrm{vds}\\ C=\underset{\∂\mathrm{\Ω}}{\∫}\frac{1}{2\mathrm{\ζ}}\mathrm{\Φvds}\\ F=\underset{k_{p_0}^{*}}{\∫\∫}\mathrm{qvdxdy}\end{array}$

其中μ_a代表吸收系数，D代表扩散系数，表示以p₀为中心的控制体单元，表示以p₀为中心的控制体单元的边界，表示对边界单位外法向的偏导数，ζ是一个与边界折射率失配程度有关的常系数，q表示光源表达式，v表示检验函数；

采用求解线性方程组的迭代法求解此方程组，从而得到生物组织区域每一个节点处的光子密度值，并输出边界处的光子密度值。

本发明的技术效果：建立了基于有限体积法的扩散光学断层成像正向问题处理方法，可以在保证精度和准确度的基础上较快地完成对二维生物组织模型的扩散光学断层成像正向过程的计算，即已知生物组织特性参数分布和光源信息前提下得到生物组织模型内及边界处的光子密度值分布。由于扩散光学图像重建的过程需要对一个病态性很高的方程组计算，其中不乏不满足实际情况的解，而采用有限体积法从方程的建立上就是基于守恒的物理思想，从而避免了这种情况，极大减小了重建过程的病态性，提高了扩散光学图像重建的准确性，这是相对于其他算法的一个独特的优势。

附图说明

图1为本发明的基本流程图。

图2 为实施例中进行的原始三角剖分和重心对偶剖分模型图，(a)表示内部控制体单元(阴影部分)；(b)表示包含边界的控制体单元。

图3为实施例下组织体样品模型和异质体位置设置图。

图4为实施例下有限体积法计算扩散光学断层成像正向过程获得的光子密度值分布图(取对数后的结果)。

图5为实施例下样品外边界圆周上光子密度值同蒙特卡罗方法所得结果的比较图。

具体实施方式

本发明提供了一种基于有限体积法的扩散光学断层成像正向问题处理方法，按以下步骤进行：创建组织体模型：使用吸收系数、散射系数、折射率和各向异性因子四个参数描述该生物组织不同区域内的光学特性；通过对待求解生物组织区域进行原始的三角剖分和重心对偶剖分，从而将原来连续的生物组织区域转化为离散的三角单元和控制体单元；在每一个控制体单元上对稳态扩散方程进行积分，形成控制体单元的方程，通过有限体积法守恒特性的思想，将所有不重叠但相互联系的控制体单元方程进行组合形成总的生物组织区域上的方程组来对生物组织区域上的光子传输进行描述，达到了将连续过程离散化的目的。

下面结合附图和实例对本发明作进一步的说明：

如图1所示，本发明的实施步骤如下：

(1)首先建立起组织体模型，比如样品背景和异质体的大小、位置，设置好各生物组织区域的光学特性参数：吸收系数、散射系数、折射率和各向异性因子四个参数；

(2)对待求解生物组织区域进行原始剖分和对偶剖分

对生物组织体进行原始三角剖分和重心对偶剖分，形成三角单元和控制体单元，剖分案例如图2所示，剖分的精细程度决定了求解的精度。将吸收系数和散射系数等参数值对应离散到每一个控制体单元的节点上，此处采用控制体形心算法下的重心对偶剖分算法，剖分的模型如2所示。图2(a)表示的是不包含组织体边界三角形构成的控制体(阴影区), 图2(b)表示的是包含组织体边界的控制体. 其中P_i为三角形的顶点，G_i为三角形的重心，m_i为三角形各边中点，阴影部分是一个控制体。；

(3)各个控制体单元上扩散方程的积分

对于步骤(2)形成的控制体单元分为边界控制体单元和内部控制体单元，边界控制体单元上通过Green 公式 ，其中Φ表示该节点处的光子密度值，表示控制体单元，表示控制体单元的边界，表示对边界单位外法向的偏导数,转化为对控制体边界的积分。此时控制体的边界包括带求解域的边界，所以要结合所使用的Robin边界条件对区域边界进行处理。以试探函数空间Φ_h的取值作为各节点处光子密度的近似值；结合稳态扩散方程，采用双线性的方法，定义检验函数空间v_h为分片常数空间，定义如下：

$v_h(x,y)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ψ}}_{p_0}(x,y)v_{p_0}$

其中   ${\mathrm{\ψ}}_{p_0}\left(p\right)=\left\{\begin{array}{cc}1,& p\∈K_{p0}^{*}\\ 0,& p\∉K_{p0}^{*}\end{array}\right.$

其中表示基底，表示该节点处的基函数，p表示节点编号，表示以p₀为中心的控制体单元, N表示节点个数；在每一个控制体上进行积分，从而得到每一个控制体上的扩散方程：

$\underset{k_{p_0}^{*}}{\∫\∫}(-\▿\·D\▿+{\mathrm{\μ}}_a)\mathrm{\Φvdxdy}=\underset{k_{p_0}^{*}}{\∫\∫}\mathrm{qvdxdy}---\left(8\right)$

其中Φ为光子密度值，μ_a代表吸收系数，D代表扩散系数，表示以p₀为中心的控制体单元， q表示光源，v表示检验函数；

变分形式为：

$a(\mathrm{\Φ},v)=(q,v)$

其中a(Φ,v)表示双线性泛函，Φ为光子密度值，q表示光源，v表示检验函数；使用Green 公式和Robin边界条件：

$\mathrm{\Φ}\left(r\right)+2D\·\frac{1+R_f}{1-R_f}\·\▿\mathrm{\Φ}\left(r\right)=0,\∀r\∈\∂\mathrm{\Ω}$

其中Φ为光子密度值，D代表扩散系数，R_f为扩散传输内反射系数，r表示节点向量；

将对区域面积的积分转化为对边界的曲线积分：

$\underset{{\∂k}_{p_0}^{*}}{\∫}-D\frac{\∂\mathrm{\Φ}}{\∂n}\mathrm{vds}+\underset{\∂\mathrm{\Ω}}{\∫}\frac{1}{2\mathrm{\ζ}}\mathrm{\Φvds}+\underset{k_{p_0}^{*}}{\∫\∫}{\mathrm{\μ}}_a\mathrm{\Φvdxdy}=\underset{k_{p_0}^{*}}{\∫\∫}\mathrm{qvdxdy}$

其中其中Φ为光子密度值，D代表扩散系数，μ_a代表吸收系数， 表示以p₀为中心的控制体单元，表示以p₀为中心的控制体单元的边界，表示生物组织边界, 表示对边界单位外法向的偏导数，ζ是一个与边界折射率失配程度有关的常系数，q表示光源，v表示检验函数；

(4)控制体单元方程的组合形成总生物组织区域上的矩阵方程，求解，输出边界处光子密度值。

将步骤(3)中得到的每一个控制体单元的控制体方程按照节点顺序进行组合，形成总的矩阵方程，得到与稳态扩散方程相对应的线性方程组，进而可以将整个区域的扩散方程写作如下矩阵的形式：

$\mathrm{K\Φ}=(A+B+C)\mathrm{\Φ}=F$

矩阵K代表总的控制体单元刚度矩阵，Φ为光子密度值，A,B,C,F分别对应下式中的矩阵元素，如下所示：

$\begin{array}{c}A=\underset{k_{p_0}^{*}}{\∫\∫}{\mathrm{\μ}}_a\mathrm{\Φvdxdy}\\ B=\underset{{\∂k}_{p_0}^{*}}{\∫}-D\frac{\∂\mathrm{\Φ}}{\∂n}\mathrm{vds}\\ C=\underset{\∂\mathrm{\Ω}}{\∫}\frac{1}{2\mathrm{\ζ}}\mathrm{\Φvds}\\ F=\underset{k_{p_0}^{*}}{\∫\∫}\mathrm{qvdxdy}\end{array}$

其中μ_a代表吸收系数，D代表扩散系数，表示以p₀为中心的控制体单元，表示以p₀为中心的控制体单元的边界，表示对边界单位外法向的偏导数，ζ是一个与边界折射率失配程度有关的常系数，q表示光源表达式，v表示检验函数；

采用求解线性方程组的迭代法求解此方程组，从而得到生物组织区域每一个节点处的光子密度值，并输出边界处的光子密度值，为扩散光学断层成像的图像重建做好准备。

对于使用本发明的技术方案处理扩散光学断层成像正向过程的结果，采用matlab进行模拟、可视化是很好的实例。

下面通过实例进一步阐述本发明。

实施例：

以图2所示高散射圆形介质为生物组织模型，采用有限体积法处理此模型下的扩散光学断层成像正向过程。模型取参数值如下：均匀背景吸收和散射系数分别为0.002mm^-1和1mm^-1，异向性因子g为0.01，介质折射率为1.5，半径为30mm。两个异质体吸收系数和散射系数分别为0.05mm^-1和5 mm^-1，半径分别为3mm和5mm，位置如图3所示。准直光从如图3左方沿径向入射。有限体积法计算结果光子密度值分布的对数结果如图4所示。为了对有限体积法的可行性和准确性进行验证，使用光学分子成像领域公认金标准:蒙特卡罗法模拟来获得参考值。图像重建工作主要关注的是边界上的光子密度值，利用边界上得到的光子密度值进行生物组织区域内部的光学参数重建，为此比较了两种方法在边界上的结果。两种方法获得边界上的结果如图5所示，为边界一周取90个参考点时的对比图，横坐标表示参考点编号，纵坐标是归一化后的光子密度值。

从图4可以看出，在设定的异质体的位置处光子密度值分布很明显改变，与实际情况吻合；从图5可以看出，在组织模型边界处，有限体积法和蒙特卡罗模拟模拟得到的结果吻合的很好。

Claims (4)

1.基于有限体积法的扩散光学断层成像正向问题处理方法，其特征在于包括以下步骤：

(1)创建组织体模型：使用吸收系数、散射系数、折射率和各向异性因子四个参数描述生物组织不同区域内的光学特性；

(2)按照精度需求的不同将待求解的生物组织区域按照现有的原始剖分和对偶剖分技术进行剖分，从而将原来连续的生物组织区域转化为离散的三角单元和控制体单元；

(3)在三角单元和控制体单元的节点处定义试探函数空间Φ_h作为该节点处光子密度的近似值，使用现有的稳态扩散方程作为光子在组织内的传输模型，采用现有的双线性方法，将稳态扩散方程在每一个控制体单元上进行积分；结合边界条件，利用Green公式其中Φ表示该节点处的光子密度值，表示控制体单元，表示控制体单元的边界，表示对边界单位外法向的偏导数,将对控制体单元面积的积分转化成对每个控制体单元上的边界积分的求和；对每一个控制体单元进行相同的稳态扩散方程积分处理，从而得到每一个控制体单元上的稳态扩散方程；

将稳态扩散方程在每一个控制体单元上进行积分，具体步骤为：以试探函数空间Φ_h的取值作为各节点处光子密度的近似值；采用双线性方法，首先定义检验函数空间v_h为分片常数空间，定义如下：

$v_h(x,y)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ψ}}_{p_0}(x,y)v_{p_0}$

其中 ${\mathrm{\ψ}}_{p_0}\left(p\right)=\left\{\begin{array}{cc}1,& p\∈K_{p0}^{*}\\ 0,& p\∉K_{p0}^{*}\end{array}\right.$

其中表示基底，表示该节点处的基函数，p表示节点编号，表示以p₀为中心的控制体单元,N表示节点个数；

在每一个控制体上进行积分，从而得到每一个控制体单元上的稳态扩散方程：

$\underset{k_{p_0}^{*}}{\∫\∫}(-\▿\·D\▿+{\mathrm{\μ}}_a)\mathrm{\Φvdxdy}=\underset{k_{p_0}^{*}}{\∫\∫}\mathrm{qvdxdy}$

其中Φ为光子密度值，μ_a代表吸收系数，D代表扩散系数，表示以p₀为中心的控制体单元，q表示光源，v表示检验函数，

变分形式为：

a(Φ,v)＝(q,v)

其中a(Φ,v)表示双线性泛函，Φ为光子密度值，q表示光源，v表示检验函数；使用Green公式和Robin边界条件：

$\mathrm{\Φ}\left(r\right)+2D\·\frac{1+R_f}{1-R_f}\·\▿\mathrm{\Φ}\left(r\right)=0,\∀r\∈\∂\mathrm{\Ω}$

其中Φ为光子密度值，D代表扩散系数，R_f为扩散传输内反射系数，r表示节点向量；将对生物组织的区域面积的积分转化为对边界的曲线积分：

其中其中Φ为光子密度值，D代表扩散系数，μ_a代表吸收系数，表示以p₀为中心的控制体单元，表示以p₀为中心的控制体单元的边界，表示生物组织边界,表示对边界单位外法向的偏导数，是一个与边界折射率失配程度有关的常系数，q表示光源，v表示检验函数；

(4)将步骤(3)中得到的所有的控制体单元的稳态扩散方程进行组合，形成总区域上的控制体单元方程，得到与稳态扩散方程相对应的线性方程组；求解该线性方程组即得到生物组织体边界处各个节点处的光子密度值。

2.根据权利要求1所述的基于有限体积法的扩散光学断层成像正向问题处理方法，其特征在于：步骤(1)中建立组织体模型包括生物组织的背景区域和位于其中的异质体的具体位置、大小及各区域光学特征参数分布情况。

3.根据权利要求1所述的基于有限体积法的扩散光学断层成像正向问题处理方法，其特征在于：步骤(2)中对待求解生物组织区域进行原始剖分和对偶剖分，形成三角单元和控制体单元，具体步骤为：通过对待求解生物组织区域进行原始的三角剖分和重心对偶剖分，将整个生物组织区域离散为控制体单元的集合；将吸收系数、散射系数、折射率和各向异性因子对应离散到每一个控制体单元的节点上。

4.根据权利要求1所述的基于有限体积法的扩散光学断层成像正向问题处理方法，其特征在于：步骤(4)中将所有的控制体单元的稳态扩散方程进行组合形成总的控制体方程，得到与稳态扩散方程相对应的线性方程组；求解该线性方程组即得到组织体边界处各个节点处的光子密度值，具体步骤为：将步骤(3)中得到的每一个控制体单元的稳态扩散方程按照节点顺序进行组合，形成总的矩阵方程，得到与稳态扩散方程相对应的线性方程组，进而将整个生物组织区域的稳态扩散方程写作如下矩阵的形式：

KΦ＝(A+B+C)Φ＝F

矩阵K代表总的控制体单元刚度矩阵，Φ为光子密度值，A,B,C,F分别对应下式中的矩阵元素，如下所示：

$A=\underset{k_{p_0}^{*}}{\∫\∫}{\mathrm{\μ}}_a\mathrm{\Φvdxdy}$

$B=\underset{{\∂k}_{p_0}^{*}}{\∫}-D\frac{\∂\mathrm{\Φ}}{\∂n}\mathrm{vds}$

$F=\underset{k_{p_0}^{*}}{\∫\∫}\mathrm{qvdxdy}$

其中μ_a代表吸收系数，D代表扩散系数，表示以p₀为中心的控制体单元，表示以p₀为中心的控制体单元的边界，表示对边界单位外法向的偏导数，是一个与边界折射率失配程度有关的常系数，q表示光源表达式，v表示检验函数；

采用求解线性方程组的迭代法求解此方程组，从而得到生物组织区域每一个节点处的光子密度值，并输出边界处的光子密度值。