发明名称 组织结构引导的复合正则化生物发光断层成像重建方法
摘要 组织结构引导的复合正则化生物发光断层成像重建方法,本方法属于医学图像处理领域;针对现有生物发光断层成像重建算法中存在的上述问题,本方法可以在已知肿瘤所在的生物组织而不知具体位置时很好的实现光源的重建;另外在不知肿瘤所在生物组织时,通过均分权重系数来实现重建。首先采用有限元方法仿真光在生物组织中的传输,并基于生物发光荧光断层成像中光源的稀疏分布特性,采用结构成像引导的复合正则化方法融入更多先验信息以降低BLT重建的病态性,最后在此基础上实现荧光光源的准确重建。在两个波段设置了双光源进行仿真实验,以验证算法的有效性和准确性。本方法不但可以对荧光光源进行准确重建定位,而且可以极大地提高计算效率。
申请公布号 CN105559750A 申请公布日期 2016.05.11
申请号 CN201510921312.0 申请日期 2015.12.13
申请人 北京工业大学 发明人 冯金超;魏慧军;贾克斌
分类号 A61B5/00(2006.01)I 主分类号 A61B5/00(2006.01)I
代理机构 北京思海天达知识产权代理有限公司 11203 代理人 沈波
主权项 组织结构引导的复合正则化生物发光断层成像重建方法,其特征在于:首先采用有限元方法仿真光在生物组织中的传输,并基于生物发光荧光断层成像中光源的稀疏分布特性,采用结构成像引导的复合正则化方法融入更多先验信息以降低BLT重建的病态性,最后在此基础上实现荧光光源的准确重建;具体包括下述步骤:步骤一,光在生物组织中的传输规律用扩散方程加以表示;根据变分原理,结合Robin边界条件,建立与扩散近似方程和边值问题等价的弱形式:<maths num="0001" id="cmaths0001"><math><![CDATA[<mrow><msub><mo>&Integral;</mo><mi>&Omega;</mi></msub><mrow><mo>(</mo><mrow><mo>(</mo><mi>D</mi><mrow><mo>(</mo><mi>r</mi><mo>)</mo></mrow><mrow><mo>(</mo><mo>&dtri;</mo><mi>&Phi;</mi><mrow><mo>(</mo><mi>r</mi><mo>)</mo></mrow><mo>)</mo></mrow><mo>&CenterDot;</mo><mo>&dtri;</mo><mi>&Psi;</mi><mrow><mo>(</mo><mi>r</mi><mo>)</mo></mrow><mo>)</mo></mrow><mo>+</mo><msub><mi>&mu;</mi><mi>a</mi></msub><mi>&Phi;</mi><mrow><mo>(</mo><mi>r</mi><mo>)</mo></mrow><mi>&Psi;</mi><mrow><mo>(</mo><mi>r</mi><mo>)</mo></mrow><mo>)</mo></mrow><mi>dr</mi><mo>+</mo><msub><mo>&Integral;</mo><mrow><mo>&PartialD;</mo><mi>&Omega;</mi></mrow></msub><mi>&Phi;</mi><mrow><mo>(</mo><mi>r</mi><mo>)</mo></mrow><mi>&Psi;</mi><mrow><mo>(</mo><mi>r</mi><mo>)</mo></mrow><mo>/</mo><mrow><mo>(</mo><mn>2</mn><msub><mi>I</mi><mi>n</mi></msub><mrow><mo>(</mo><mi>r</mi><mo>)</mo></mrow><mo>)</mo></mrow><mi>dr</mi><mo>=</mo><msub><mo>&Integral;</mo><mi>&Omega;</mi></msub><mi>S</mi><mrow><mo>(</mo><mi>r</mi><mo>)</mo></mrow><mi>&Psi;</mi><mrow><mo>(</mo><mi>r</mi><mo>)</mo></mrow><mtext>dr---</mtext><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000876342690000016.GIF" wi="1886" he="118" /></maths>其中D(r)为扩散系数分布,Φ(r)为光子密度分布,μ<sub>a</sub>为光吸收系数,Ψ(r)为任意测试函数,S(r)为光源分布,I<sub>n</sub>(r)为使用Robin边界条件时引入的一个量,<img file="FDA0000876342690000011.GIF" wi="303" he="141" />R<sub>n</sub>≈‑1.4399n<sup>‑2</sup>+0.7099n<sup>‑1</sup>+0.6681+0.636n表示扩算传输内反射系数,n是一个常数,与边界内外的光学折射系数偏差有关;步骤二,有限元离散;对成像体Ω进行网格剖分,实现离散化;本方法采用二维仿体进行仿真,由于三角网格单元能够实现使用较少数量来近似描述不规则区域,所以,本方法使用三角网格剖分成像物体;根据对近似解精度的要求、单元的几何形状、节点个数和节点的自由度信息来选择单元基函数;由于函数均表示为基函数的线性组合,假设成像物体所划分的网格包含N<sub>n</sub>个节点和N<sub>t</sub>个三角单元,光子密度分布Φ(r)表示为:<img file="FDA0000876342690000012.GIF" wi="1305" he="143" />其中,φ<sub>k</sub>(r)表示节点i上的函数值,<img file="FDA0000876342690000013.GIF" wi="115" he="79" />表示节点i处的基函数;同样,光源函数S(r)表示为:<img file="FDA0000876342690000014.GIF" wi="1269" he="140" />其中,s<sub>k</sub>(r)表示节点i上的光源值;步骤三,合成总体刚度矩阵;将(2)、(3)式代入(1)式中,得到:<img file="FDA0000876342690000015.GIF" wi="2015" he="343" />该矩阵方程简单表示为:([K]+[C]+[B])=MΦ=FS   (5)其中,M为刚度矩阵,其他各矩阵元素通过下面公式计算得到:<img file="FDA0000876342690000021.GIF" wi="1212" he="446" />整理(6)式,得到:M<sup>‑1</sup>FS=Φ   (7)进一步地,上式写成如下形式:AS=Φ   (8)其中,A=M<sup>‑1</sup>F;对于生物发光断层成像,探测器所能够采集到的测量数据仅仅是成像物体表面的荧光信息,因此在(7)式中,需要将光子密度向量Φ中的非表面元素去掉,并将系统矩阵A中对应的行去掉,从而得到不含表面荧光信息的光子密度分布Φ<sup>m</sup>和系统矩阵A,此时的光源重建问题表述为:AS=Φ<sup>m</sup>   (9)步骤四,利用结构成像得到组织结构;利用micro‑CT或micro‑MRI等结构成像提供的组织结构信息,通过对生物组织进行分割,得到不同的生物组织,将此作为先验信息融合到后续的重建中;步骤五,采用复合正则化方法进行重建;目标函数设置如下(10)式所示:<maths num="0002" id="cmaths0002"><math><![CDATA[<mrow><mi>m</mi><mi>i</mi><mi>n</mi><mo>{</mo><mi>f</mi><mrow><mo>(</mo><mi>s</mi><mo>)</mo></mrow><mo>}</mo><mo>=</mo><mo>|</mo><mo>|</mo><mi>A</mi><mi>s</mi><mo>-</mo><msup><mi>&Phi;</mi><mi>m</mi></msup><mo>|</mo><msubsup><mo>|</mo><mn>2</mn><mn>2</mn></msubsup><mo>+</mo><mi>&lambda;</mi><mo>|</mo><mo>|</mo><mi>s</mi><mo>|</mo><msub><mo>|</mo><mn>1</mn></msub><mo>+</mo><mi>&beta;</mi><munderover><mo>&Sigma;</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>q</mi></munderover><msub><mi>&omega;</mi><mi>i</mi></msub><mo>|</mo><mo>|</mo><msub><mi>s</mi><mi>G</mi></msub><mo>|</mo><msub><mo>|</mo><mn>2</mn></msub><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>10</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000876342690000022.GIF" wi="1390" he="135" /></maths>上式参数的意义如下:A:利用有限元求解扩散方程得到的系统矩阵;s:是未知量,待求解;Φ<sup>m</sup>:测量值;λ:正则化参数1;β:正则化参数2;ω<sub>i</sub>:不同组织的权重;s<sub>G</sub>:表示其余未知量s分解为的q个不相重合的部分,q可表示将结构图像分割后分成的q个不同组织,ω<sub>i</sub>表示对不同组织的权重;通过求上述最优化问题,可重建出未知量s的分布情况。
地址 100124 北京市朝阳区平乐园100号