CN103218813A - 一种基于边界元的dot/xct双模式成像的图像重建算法 - Google Patents

Abstract

本发明涉及一种基于边界元的DOT/XCT双模式成像的图像重建算法，属于分子影像技术领域。所述的方法包括在XCT提供组织边界的基础上，利用边界元法求解DOT的正向问题，即在已知不同组织体光学参数分布的情况下求出光通量在边界处的分布。逆问题中采用Levenberg-Marquardt算法对不同组织体的吸收系数和散射系数进行了重建，即从探测到的光通量数据出发，通过对非线性最小二乘模型的优化过程求出不同组织体的吸收系数和散射系数。与传统的基于有限元的重建相比，本发明在DOT/XCT双模式中利用边界进行重建，实现了对组织体光学参数更加准确，稳定的重建。

Description

一种基于边界元的DOT/XCT双模式成像的图像重建算法

技术领域

本发明属于分子影像技术领域，涉及一种图像重建的算法，特别涉及一种基于边界元的DOT/XCT双模式成像的图像重建算法。

背景技术

扩散光学层析成像技术(Diffuse Optical Tomography，DOT)是一种非入侵(non-invasive)、高灵敏度的光学成像手段，它利用穿过组织的漫射光信息结合组织光传输模型可以重建出组织体的某一个断层面甚至三维空间的光学或生理参数分布（Schweiger, M., et al.,3D level set reconstruction of model and experimentaldata in Diffuse Optical Tomography. Optics express,2010. 18(1): p. 150-164.）。扩散光学层析成像本质上是一个非线性问题，其逆向问题具有严重的病态性。为了提高重建的准确性，近年来多模式成像技术成为一个新的发展趋势（Ale, A., et al., FMT-XCT: in vivo animal studies with hybrid fluorescence molecular tomography-X-ray computed tomography.Nat Meth, 2012. 9(6): p. 615-620.），如微型X射线计算机层析成像（MicroXCT）以及一些其他成像方式可以将结构信息引入到DOT的重建问题中，改善逆问题的病态性，从而提高了重建算法的精度和速度。

现有的DOT/MicroXCT多模式成像中光学参数的重建，一般采用的是有限元法求解DOT的正向问题（Arridge, S.R., Optical tomography in medical imaging. Inverse problems, 1999. 15(2):p. R41.），逆向问题中采用非线性共轭梯度法、牛顿法、非线性Landweber迭代法求解，而MicroXCT的结构信息通过软先验或硬先验方法引入。近几年，由于边界元法对于微分方程中所考虑变量的导数(如光流量)的求解比有限元精确，因而理论上边界元法可以达到更高的正向问题计算精度和更优的重建效果。剖分方式上，边界元法只需分别独立地剖分各包含组织器官的表面，比有限元中所需要的体元剖分更容易实现。并且边界元法的离散处理仅仅涉及到边界，整个区域中不再出现待求参数，使得待求参数的数目大大减小，方程组规模缩小，实现将所求的问题降维，在处理大规模问题时效率也更高，因此目前成为DOT成像中解决正向问题的研究热点（Srinivasan, S., et al., A boundary element approach for image-guided near-infrared absorption and scatter estimation. Medical physics, 2007. 34: p. 4545.）。

发明内容

本发明的目的在于：为了在DOT/XCT双模式成像系统中更为精确和稳定的重建组织体的光学参数，而提供一种基于边界元的DOT/XCT双模式成像的图像重建算法，本发明有别于广泛使用的基于有限元的重建方法，引入边界元法作为计算正向问题的数值解法，特别适用于在XCT提供组织边界的情况下结合Levenberg-Marquardt算法对不同组织体的吸收系数和散射系数进行重建，使得重建的结果更加精确、稳定。

本发明的技术方案是：

一种基于边界元的DOT/XCT双模式成像的图像重建算法，其特征包括以下步骤：(1)通过XCT对小鼠成像，得到不同角度的投影图，利用CT的重建算法重建出不同位置的组织体的CT切片图，将得到的CT切片图进行图像分割，分离出不同组织体的边界，通过边界提取，得到不同组织体的边界的坐标点集；通过DOT系统对同一只小鼠进行成像，采集鼠体边界处不同位置的光通量值，将这些光通量集作为探测数据集；(3)将第（1）提取的各组织体的边界的坐标和第（2）步提取的鼠体边界处不同位置的光通量值导入边界元程序，从一组初始的光学参数出发，利用Levenberg-Marquardt算法对不同组织体的吸收系数和散射系数进行重建，每步迭代中正向问题的计算通过边界元实现，最终得到小鼠不同器官光学参数分布图。

步骤(2)利用DOT成像系统对小鼠进行整体成像，在小鼠体表选取多个光源的位置和多个探测器的位置，得到不同光源和探测器组合下的鼠体边界处不同位置的光通量值，将鼠体边界处不同位置的光通量值作为探测数据集供后续重建使用。

步骤(3) 利用Levenberg-Marquardt算法对不同组织体的吸收系数和散射系数进行重建，具体步骤为：将重建过程建模为一个非线性最小二乘的优化问题：

$\underset{x}{\min }\mathrm{\Ω}={\left|\right|y-F\left(x\right)\left|\right|}^2$

y为探测数据集，F(x)为带入不同组织体的光学参数x（吸收系数μ_a，散射系数μ_s）正向计算得到的预测数据集，Ω为目标函数；正向问题的计算采用边界元方法：首先将第（1）步中提取的边界坐标导入边界元程序，对不同组织体边界进行剖分，带入相应组织体的光学参数值后，利用边界元这种数值解法对一组耦合的扩散方程进行求解：

${\▿}^2{\mathrm{\Φ}}_l\left(r\right)-{\mathrm{\ω}}_l^2{\mathrm{\Φ}}_l\left(r\right)=-q_l,l=\mathrm{1,2}...L$

$\left\{\begin{array}{c}{\mathrm{\ω}}_l^2=\frac{{\mathrm{\μ}}_{\mathrm{\α},l}}{{\mathrm{\κ}}_l}\\ q_l=\frac{q_1\left(r\right)}{{\mathrm{\κ}}_1}\·{\mathrm{\δ}}_{i,1}\end{array}\right.$

其中， Φ(r)为r处的光通量大小，μ_α,l为第l个组织体内部的吸收系数，q(r)为光源项，δ_i,1为狄拉克函数。L为子区域的个数，q_l为光源项，κ_l为扩散系数，ω_l为亥姆霍兹方程中的常数项。

求解即可得到边界处的光通量分布；为了重建出不同组织体真实的光学参数，重建算法首先由一组初始猜测的光学参数值，通过边界元进行正向计算从而得到边界处的光通量分布，结合探测到的真实数据，得到光通量关于光学参数的导数信息——雅各比矩阵

其中f_i为第i个测量值，x_i为第i个光学参数

然后利用导数信息更新光学参数：

$x_{i+1}=x_i-{(J^{T}J+\mathrm{\αI})}^{-1}{J}^T(y-F\left(x_i\right))$

其中x_i为第i次迭代中不同组织体的光学参数值，α为一个正实数，y为测量数据集，F(x_i)为以光学参数x_i进行正向计算所得到的光通量数据集，I为单位矩阵。

如此重复，直至从某一步的光学参数值计算所得到的光通量和探测到的光通量之间的误差小于相应范围，则该步的光学参数值就是重建出的最终结果。

与传统的基于有限元的重建相比，本发明在DOT/XCT双模式中利用边界进行重建，实现了对组织体光学参数更加准确、稳定的重建。

附图说明

图1为本发明实施实例提供的一种基于边界元的DOT/XCT双模式成像的图像重建算法流程图。

图2为本发明实施例提供的边界提取示意图。

图3为本发明实施例提供的吸收系数和扩散系数重建结果示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面结合附图对本发明作进一步的详细描述。

本发明提出一种基于边界元的DOT/XCT双模式成像的图像重建算法，首先利用XCT提供的组织边界引入边界元中，利用边界元作为正向问题的数值解法，然后利用Levenberg-Marquardt算法不断优化迭代待重建的光学参数，使每次正向计算的结果逐步逼近真实的探测数据，直至满足一定的终止条件，从而重建出小鼠体内不同器官的光学参数分布。具体的流程如图1所示。其主要步骤如下：

(1)利用XCT对小鼠进行成像：利用XCT对小鼠整体进行成像，得到不同角度的组织体投影图，然后利用相应的CT切片重建算法重建出一组不同位置的小鼠内部组织体的CT切片图，如图2所示。

(2)不同组织体边界坐标的获取：将上一步得到的CT切片图经过图像分割，分离出不同组织体边界的信息，再经过边界读取算法得到各组织体的边界坐标。

(3)获取探测数据：利用DOT成像在小鼠表面指定多个光源和探测器的位置，分别采集不同的光源位置下及探测器位置的光通量大小作为探测数据；或者可指定一组特定的组织体光学参数，然后通过边界元求解正向问题得到一组模拟探测数据。

(4)通过重建算法重建不同组织体光学参数：设定一组初始的光学参数作为猜测值，剖分组织体边界如图2所示，通过边界求解一组耦合的扩散方程，得到正向问题的解——边界处光通量值的大小。耦合的扩散方程描述了不同组织中光子密度波的传播过程：

${\▿}^2{\mathrm{\Φ}}_l\left(r\right)-{\mathrm{\ω}}_l^2{\mathrm{\Φ}}_l\left(r\right)=-q_l,l=\mathrm{1,2}...L$

$\left\{\begin{array}{c}{\mathrm{\ω}}_l^2=\frac{{\mathrm{\μ}}_{\mathrm{\α},l}}{{\mathrm{\κ}}_l}\\ q_l=\frac{q_1\left(r\right)}{{\mathrm{\κ}}_1}\·{\mathrm{\δ}}_{i,1}\end{array}\right.$

其中，Φ(r)为r处的光通量大小，μ_α,l为第l个组织体内部的吸收系数，q(r)为光源项，δ_i,1为狄拉克函数。L为子区域的个数，q_l为光源项，κ_l为扩散系数，ω_l为亥姆霍兹方程中的常数项。

通常采用的边界条件为Robin边界条件，表达式如下：

$\mathrm{\Φ}\left(r\right)+2\mathrm{\α}\·\mathrm{\κ}\left(r\right)\·\▿\mathrm{\Φ}\left(r\right)=0$

其中，

是一个与边界折射率失配程度有关的常系数。R_f为扩散传输内反射系数。κ(r)为扩散系数。

定义如下变量U_l和V_l，其物理意义为光子密度和光流量，在边界处的关系如下所示：

$\left\{\begin{array}{c}{U}_l:={\mathrm{\Φ}}_l{|}_{{\mathrm{\Γ}}_l}={\mathrm{\Φ}}_{l-1}{|}_{{\mathrm{\Γ}}_l}\\ V_l:=-{\mathrm{\κ}}_l\frac{{\∂\mathrm{\Φ}}_l}{\∂{\hat{v}}_l}{|}_{{\mathrm{\Γ}}_l}={\mathrm{\κ}}_{l-1}\frac{{\∂\mathrm{\Φ}}_{l-1}}{\∂{\hat{v}}_{l-1}}\end{array},{|}_{{\mathrm{\Γ}}_l}\≡{\mathrm{\κ}}_{l-1}\frac{{\∂\mathrm{\Φ}}_{l-1}}{\∂\hat{v}}{|}_{{\mathrm{\Γ}}_l}={\mathrm{\κ}}_l\frac{{\∂\mathrm{\Φ}}_l}{\∂\hat{v}}{|}_{{\mathrm{\Γ}}_l}\right.$

其中

代表第l个组织边界的外法向，Γ_l为表示第l条边界。

利用格林第二公式可得积分方程：

式中的c_l是一个与边界处奇异性相关的常数，当边界光滑时

 ，式中G_l为第l个Helmholtz方程的基本解，Ω_l表示第l个子区域。

二维情况下Helmholtz方程基本解和基本解的导数形式如下：

$\left\{\begin{array}{c}G(r,r^{\′})=\frac{1}{4j}{H}_0^{\left(2\right)}(-\mathrm{j\ω}|r-r^{\′}\left|\right)\\ \frac{\∂G(r,r^{\′})}{\∂\hat{v}}=\hat{v}\·\frac{\mathrm{\ω}(r-r^{\′})}{|r-r^{\′}|}{H}_1^{\left(2\right)}(-\mathrm{j\ω}|r-r^{\′}\left|\right)\end{array}\right.$

利用常单元将所有边界积分方程离散化后可得到一个线性方程组：

$\left[\begin{array}{ccccc}\frac{1}{2}I+A_{\mathrm{1,1}}^1+\frac{1}{2\mathrm{\α}}{B}_{\mathrm{1,1}}^1& A_{\mathrm{1,2}}^1& B_{\mathrm{1,2}}^1& ...& B_{1,n}^1\\ A_{\mathrm{2,1}}^1+\frac{1}{2\mathrm{\α}}{B}_{\mathrm{2,1}}^1& \frac{1}{2}I-A_{\mathrm{2,2}}^1& B_{\mathrm{2,2}}^1& ...& B_{2,n}^1\\ ...& ...& ...& ...& ...\\ 0& \frac{1}{2}I+A_{22}^2& B_{22}^2& ...& 0\\ .& & & & .\\ .& ...& ...& ...& .\\ .& & & & .\\ 0& ...& ...& \frac{1}{2}I+A_{\mathrm{nn}}^n& B_{\mathrm{nn}}^n\end{array}\right]\left[\begin{array}{c}{U}_1\\ U_2\\ V_2\\ .\\ .\\ .\\ U_n\\ V_n\end{array}\right]=\left[\begin{array}{c}{Q}_1{|}_{{\mathrm{\Γ}}_1}\\ .\\ .\\ .\\ Q_1{|}_{{\mathrm{\Γ}}_n}\\ 0\\ .\\ .\\ .\\ 0\end{array}\right]$

求解该线性方程组就可以得到边界处的光通量大小。

逆问题建模为一个非线性最小二乘的优化问题：

$\underset{x}{\min }\mathrm{\Ω}={\left|\right|y-F\left(x\right)\left|\right|}^2$

y为测量数据集，F(x)为待重建参数x（吸收系数μ_a，散射系数μ_s）正向计算得到的预测数据集，Ω为目标函数。构建雅各比矩阵如下所示：

m为源个数和探测器个数的乘积，也即测量数据集包含的数据个数，n为带重建参数的个数，边界元法是基于每个异质区域光学参数相同的模型，与基于有限元的重建相比，最大的不同在于这里的重建问题从原来的负定问题(under-dertermined)变成一个过定问题(over-determined)，即只需重建每个异质区域的光学参数，使得带重建参数个数大大的减少。利用雅各比矩阵提供的导数信息，结合Levenberg-Marquardt算法，并采用信赖域策略，可以得到光学参数的迭代更新公式：

$x_{i+1}=x_i-{(J^{T}J+\mathrm{\αI})}^{-1}{J}^T(y-F\left(x_i\right))$

其中I为单位矩阵，该算法的巧妙之处在于通过控制参数α的大小使得搜索的方向在最速下降方向和Gauss-Newton方向之间调整，保证了全局收敛。

通过上述步骤不断迭代更新光学参数，每步迭代中正向问题的求解都采用边界元法，最终直到依据猜测的光学参数所计算的正向数据和探测数据之间的误差低于一定的阈值，或满足其他终止条件后，便可以重建出不同组织体的光学参数分布，绘制成图如图3所示。

总之，以上所述仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。

Claims (3)

1.一种基于边界元的DOT/XCT双模式成像的图像重建算法，其特征包括以下步骤：(1)通过XCT对小鼠成像，得到不同角度的投影图，利用CT的重建算法重建出不同位置的组织体的CT切片图，将得到的CT切片图进行图像分割，分离出不同组织体的边界，通过边界提取，得到不同组织体的边界的坐标点集；(2)通过DOT系统对同一小鼠不同位置组织体的进行成像，采集各组织体边界处的光通量值，将各组织体边界处的光通量值作为探测数据集；(3)将第（1）提取的各组织体的边界的坐标和第（2）步提取的各组织体边界处的光通量值导入边界元程序，从一组初始的光学参数出发，利用Levenberg-Marquardt算法对不同组织体的吸收系数和散射系数进行重建，每步迭代中正向问题的计算通过边界元实现，最终得到小鼠不同器官光学参数分布图。

2.根据权利要求1所述的一种基于边界元的DOT/XCT双模式成像的图像重建算法，其特征在于的：步骤(2)利用DOT成像系统对小鼠进行整体成像，在小鼠体表选取多个光源的位置和多个探测器的位置，得到不同光源和探测器组合下的各组织体边界处光通量值，将各组织体边界处光通量值作为探测数据集供后续重建使用。

3.根据权利要求1所述的一种基于边界元的DOT/XCT双模式成像的图像重建算法，其特征在于：步骤(3) 利用Levenberg-Marquardt算法对不同组织体的吸收系数和散射系数进行重建，具体步骤为：将重建过程建模为一个非线性最小二乘的优化问题：

$\underset{x}{\min }\mathrm{\Ω}={\left|\right|y-F\left(x\right)\left|\right|}^2$

y为探测数据集，F(x)为带入不同组织体的光学参数x（吸收系数μ_a，散射系数μ_s）正向计算得到的预测数据集，Ω为目标函数；正向问题的计算采用边界元方法：首先将第（1）步中提取的边界坐标导入边界元程序，对不同组织体边界进行剖分，带入相应组织体的光学参数值后，利用边界元这种数值解法对一组耦合的扩散方程进行求解：

${\▿}^2{\mathrm{\Φ}}_l\left(r\right)-{\mathrm{\ω}}_l^2{\mathrm{\Φ}}_l\left(r\right)=-q_l,l=\mathrm{1,2}...L$

$\left\{\begin{array}{c}{\mathrm{\ω}}_l^2=\frac{{\mathrm{\μ}}_{\mathrm{\α},l}}{{\mathrm{\κ}}_l}\\ q_l=\frac{q_1\left(r\right)}{{\mathrm{\κ}}_1}\·{\mathrm{\δ}}_{i,1}\end{array}\right.$

其中， Φ(r)为r处的光通量大小，μ_α,l为第l个组织体内部的吸收系数，q(r)为光源项，δ_i,1为狄拉克函数；L为子区域的个数，q_l为光源项，κ_l为扩散系数，ω_l为亥姆霍兹方程中的常数项。

即可得到边界处的光通量分布；为了重建出不同组织体真实的光学参数，重建算法首先由一组初始猜测的光学参数值，通过边界元进行正向计算从而得到边界处的光通量分布，结合探测到的真实数据，得到光通量关于光学参数的导数信息——雅各比矩阵

其中f_i为第i个测量值，x_i为第i个光学参数

然后利用导数信息更新光学参数：

$x_{i+1}=x_i-{(J^{T}J+\mathrm{\αI})}^{-1}{J}^T(y-F\left(x_i\right))$

其中x_i为第i次迭代中不同组织体的光学参数值，α为一个正实数，y为测量数据集，F(x_i)为以光学参数x_i进行正向计算所得到的光通量数据集，I为单位矩阵；

如此重复，直至从某一步的光学参数值计算所得到的光通量和探测到的光通量之间的误差小于相应范围，则该步的光学参数值就是重建出的最终结果。

Images