发明名称 一种用于自发荧光成像的空间加权的有限元重建方法
摘要 本发明提出一种用于自发荧光成像的空间加权的有限元重建方法。本发明使用扩散光学层析成像diffuse optical tomography(DOT)技术来精确重建主要生物组织的光学参数。本发明对可能的真实光源区域进行了限制,将整个重建区域分为:光源可能区域和光源不可能区域来提高重建问题的数值稳定性和有效性,并且减轻了BLT重建问题的病态性。为了避免著名的”Inverse Crime”问题,本发明采用Monte Carlo(MC)随机的方法来模拟光在生物组织中传输过程。然后,本发明提出了一种基于空间加权单元的有限元方法并采用了一种线性约束优化问题的容忍算法来重建光源信息。
申请公布号 CN101539518B 申请公布日期 2011.06.15
申请号 CN200810102314.7 申请日期 2008.03.20
申请人 中国科学院自动化研究所 发明人 田捷;石金;杨鑫;徐敏
分类号 G01N21/64(2006.01)I;A61B5/00(2006.01)I;G06T11/00(2006.01)I 主分类号 G01N21/64(2006.01)I
代理机构 中科专利商标代理有限责任公司 11021 代理人 周国城
主权项 1.一种用于自发荧光成像的空间加权的有限元重建方法,其特征在于,包括:重建感兴趣区域的光学参数;使用MC方法模拟前向问题;将MC仿真结果与重建网格进行空间距离最小化的匹配;将整个网络分为光源可行区和不可行区;使用空间加权单元的有限元方法将扩散方程离散化;建立内部未知光源S<sup>pr</sup>和边界处通亮强度Φ<sup>bound</sup>之间的线性关系;将BLT重建问题转化为一个有简单边界约束的优化问题;使用线性约束的容忍算法迭代求解所述优化问题,从而求得重建的光源信息;所述有限元重建方法的具体步骤如下:步骤1:扩散光学层析成像技术使用有限元方法来模拟光在组织中的传输,并得到其时域相应信号的积分,然后通过Levenberg-Marquardt方法来迭代重建感兴趣区域的组织光学特性参数的准确的空间分布;步骤2:根据CCD上捕获的光学信号分布对真实光源的分布区域进行了限制,将整个区域分成了光源可行区和光源不可行区;步骤3:使用Monte Carlo(MC)随机的方法来模拟光在生物组织中传输过程,并将MC网格上的仿真结果与重建网格通过空间距离最近原则进行匹配;步骤4:通过使用有限元分析的方法,稳态扩散方程被描述为相应的弱解形式:<maths num="0001"><![CDATA[<math><mrow><mo>-</mo><mo>&dtri;</mo><mo>&CenterDot;</mo><mrow><mo>(</mo><mi>D</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>&dtri;</mo><mi>&Phi;</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>)</mo></mrow><mo>+</mo><msub><mi>&mu;</mi><mi>a</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mi>&Phi;</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>=</mo><mi>S</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mrow><mo>(</mo><mi>x</mi><mo>&Element;</mo><mi>&Omega;</mi><mo>)</mo></mrow></mrow></math>]]></maths><maths num="0002"><![CDATA[<math><mrow><mi>&Phi;</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>+</mo><mn>2</mn><mi>A</mi><mrow><mo>(</mo><mi>x</mi><mo>;</mo><mi>n</mi><mo>,</mo><mi>n</mi><mo>&prime;</mo><mo>)</mo></mrow><mi>D</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mrow><mo>(</mo><mi>v</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>&CenterDot;</mo><mo>&dtri;</mo><mi>&Phi;</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>)</mo></mrow><mo>=</mo><mn>0</mn><mrow><mo>(</mo><mi>x</mi><mo>&Element;</mo><mo>&PartialD;</mo><mi>&Omega;</mi><mo>)</mo></mrow></mrow></math>]]></maths><maths num="0003"><![CDATA[<math><mrow><msub><mo>&Integral;</mo><mi>&Omega;</mi></msub><mrow><mo>(</mo><mi>D</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mrow><mo>(</mo><mo>&dtri;</mo><mi>&Phi;</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>)</mo></mrow><mo>&CenterDot;</mo><mo>(</mo><mo>&dtri;</mo><mi>&Psi;</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>)</mo><mo>+</mo><msub><mi>&mu;</mi><mi>a</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mi>&Phi;</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mi>&Psi;</mi><mo>(</mo><mi>x</mi><mo>)</mo><mo>)</mo></mrow><mi>dV</mi></mrow></math>]]></maths><maths num="0004"><![CDATA[<math><mrow><mo>+</mo><msub><mo>&Integral;</mo><mrow><mo>&PartialD;</mo><mi>&Omega;</mi></mrow></msub><mfrac><mn>1</mn><mrow><mn>2</mn><mi>A</mi><mrow><mo>(</mo><mi>x</mi><mo>;</mo><mi>n</mi><mo>,</mo><mi>n</mi><mo>&prime;</mo><mo>)</mo></mrow></mrow></mfrac><mi>&Phi;</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mi>&Psi;</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mi>dA</mi><mo>=</mo><msub><mo>&Integral;</mo><mi>&Omega;</mi></msub><mi>S</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mi>&Psi;</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mi>dV</mi><mo>,</mo><mo>&ForAll;</mo><mi>&Psi;</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>&Element;</mo><msup><mi>H</mi><mn>1</mn></msup><mrow><mo>(</mo><mi>&Omega;</mi><mo>)</mo></mrow></mrow></math>]]></maths>其中,Ω和<img file="FSB00000467965400015.GIF" wi="73" he="45" />分别是感兴趣区域和其边界;Φ(x)代表光子流量密度的分布;S(x)代表了荧光光源的密度分布;D(x)=1/[3(μ<sub>a</sub>(x)+(1-g)μ<sub>s</sub>(x))]是扩散系数;μ<sub>a</sub>(x)是吸收系数,而μ<sub>s</sub>(x)是散射系数,g是各向异性参数;v是边界点x处的外法线方向;A=(1+R)/(1-R),R依赖于媒介的反射性质并且能够被近似为:R≈-1.4399n<sup>-2</sup>+0.7099n<sup>-1</sup>+0.6681+0.0636n;步骤5:根据标准的有限元方法,对于任意一个分段连续的测试函数Ψ(x)∈C<sup>0</sup>(Ω),则Φ(x)按如下离散化:<maths num="0005"><![CDATA[<math><mrow><mi>&Phi;</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>&ap;</mo><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>Np</mi></munderover><msub><mi>&phi;</mi><mi>i</mi></msub><msub><mi>&Psi;</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow></mrow></math>]]></maths>其中{ψ1,...,ψnp}是离散化后的空间A的基函数,φi是相应的第i个节点值;令{ζ1,...,ζNs}为插值基函数,则S(x)近似为:<maths num="0006"><![CDATA[<math><mrow><mi>S</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>&ap;</mo><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>Ns</mi></munderover><msub><mi>S</mi><mi>i</mi></msub><msub><mi>&zeta;</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow></mrow></math>]]></maths>其中Si为相应的插值节点值,Ns为插值基函数的个数;步骤6:将步骤4中积分式等号左侧的三个积分项转化成矩阵乘积的形式MΦ,而步骤4中积分式等号右面的积分项如下处理,对于任意一个四面体有:∫S(x)Ψ(x)dV=[ω<sub>l1</sub> ω<sub>l2</sub> ω<sub>l3</sub> ω<sub>l4</sub>]*[ψ<sub>i1</sub> ψ<sub>i2</sub> ψ<sub>i3</sub> ψ<sub>i4</sub>]<sup>T</sup>*S<sub>i</sub>其中S<sub>i</sub>代表的第i个重建单元的光源重建值;ωi1、ωi2、ωi3、ωi4是该单元上的节点值;ψi1、ψi2、ψi3、ψi4是相应的基函数,用以描述该单元的空间位置和形状;步骤7:整个有限元网格上建立荧光光源密度分布和光子通量密度分布的线性关系:MΦ=FS通过删除内部的节点并且选择光源可行区域的方法,建立内部未知光源S<sup>pr</sup>和边界处通量强度Φ<sup>bound</sup>之间的线性关系:M<sup>mod</sup>Φ<sup>bound</sup>=F<sup>mod</sup>S<sup>pr</sup>其中,M<sup>mod</sup>为M矩阵按Φ<sup>bound</sup>所在行删除,按S<sup>pr</sup>所在列删除,剩余的部分,F<sup>mod</sup>为F矩阵按S<sup>pr</sup>所在行删除,按Φ<sup>bound</sup>所在列删除,剩余的部分;步骤8:将BLT重建问题转化为一个有简单界约束的优化问题:<maths num="0007"><![CDATA[<math><mrow><msup><mi>S</mi><mi>low</mi></msup><mo>&le;</mo><mover><msup><mi>S</mi><mi>pr</mi></msup><mi>min</mi></mover><mo>&lt;</mo><mo>&lt;</mo><msup><msup><mi>S</mi><mi>up</mi></msup><mrow><mo>{</mo><mo>|</mo><mo>|</mo><msup><mi>M</mi><mi>mod</mi></msup><msup><mi>&Phi;</mi><mi>bound</mi></msup><mo>-</mo><msup><mi>F</mi><mi>mod</mi></msup><msup><mi>S</mi><mi>pr</mi></msup><msub><mrow><mo>|</mo><mo>|</mo></mrow><mo>^</mo></msub><mo>+</mo><mi>&lambda;&xi;</mi><mrow><mo>(</mo><msup><mi>S</mi><mi>pr</mi></msup><mo>)</mo></mrow><mo>}</mo></mrow></msup><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>10</mn><mo>)</mo></mrow></mrow></math>]]></maths>其中,S<sup>low</sup>和S<sup>up</sup>分别是荧光光源强度的上下限;∧为权重矩阵,||V||∧=V<sup>T</sup>∧V;λ为正则参数;ξ为惩罚函数。
地址 100080 北京市海淀区中关村东路95号