专利名称：基于单视图的多光谱自发荧光断层成像重建方法
技术领域：
本发明涉及到光学分子影像成像模态-自发荧光断层成像(BLT)技术，尤其涉及到一种基于单视图的多光谱自发荧光断层成像重建方法。
背景技术：
在过去的几年里，分子影像由于能在体揭示分子和细胞的信息而受到了越来越多的关注。因此它已成为一种疾病诊断、药物疗效评价等的重要工具。作为一种小动物分子影像成像模态，光学成像技术特别是自发荧光断层成像因为其高性能、低价格和无创伤等特性已经得到了迅速发展和广泛的应用。自发荧光断层成像是最近刚发展起来的，是一种用来在活体小动物体内重建自发荧光光源分布的光学成像技术。当荧光素与荧光素酶在有氧气和三磷酸腺苷(ATP)的条件下，就会发出荧光。而因为荧光素酶含有萤火虫素酶(firefly)，磕头虫素酶(click beetle)，所以发出的荧光有不同的光谱，波长一般在400nm-750nm。发出的荧光穿透生物体而到达生物体表面，然后用高灵敏度的液氮制冷电荷耦合器件(CCD)获得。CCD获得的动物表面的数据形成了BLT重建的基础。为采集整个生物体的表面数据，通常将生物体每隔90度旋转一次，用CCD探测器来获取生物体的前后左右四个视图。当光源在生物体的位置比较深时，采集一个视图的数据需要的时间大约是5-10分钟。而当采集时间超过一个小时之后荧光强度就会衰减，所以在多谱的情况下，如果在每个单谱段都采集四个视图的数据，那么采集时间远远超过1个小时，可能就会使得采集的数据不够准确。另外，在实际情况下存在类似像X射线(x-ray)，计算机断层成像(CT)一样的物理限制，采集数据时就受到角度限制，因此采集到的数据仅是动物体表面的一部分。多谱采集也会使得系统矩阵的维数增加，导致计算量过大。如何降低测量时间和减小系统矩阵的维数是目前的一个难点。另一方面，在进行光源重建时，因为BLT自身具有的不适定特点导致重建结果不准确。
BLT的目标是能在体实时连续的观测。因此，发展快速的重建方法就成为亟待解决的问题。目前的重建方法常用的是有限元方法。理论上有限元网格越细，得到的结果越好，但重建的时间就越长。此外，BLT是病态问题，如何降低病态性仍需要进一步的探索。
发明内容
本发明的目的在于克服了自发荧光断层成像重建方法重建光源不准确，重建速度比较慢，不利于实时处理以及在多光谱数据采集时可能存在误差的缺陷，提出了一种基于单视图的多光谱自发荧光断层成像重建方法，该方法基于扩散方程模型，并考虑了小动物体的非均匀特性，同时也考虑自发荧光光源的光谱及实际应用的特点。
为达到此目的，本发明提出的基于单视图的多光谱自发荧光断层成像重建方法的实现步骤主要包括1)数据获取；2)有限元离散化处理；3)优化可行光源区域选择；4)光源的重建。
1)数据获取是指利用CCD探测器在被测物体表面对逃出的光子进行探测，获得边界处的光流密度； 有限元离散化处理，是指利用有限元方法将扩散方程离散为矩阵方程，并把重建区域网格剖分为{T1，...，Tk，...}的网格序列； 优化可行光源区域的选择，主要是在对扩散方程离散化处理的基础上，利用迭代方法快速选择光源大致存在的区域，即优化可行光源区域； 光源的重建，主要是在优化可行光源区域选择的基础上，确立未知光源密度变量与被测物体表面测量值之间的线性关系，并利用Tikhonov正则化方法建立目标函数，并最终求解光源密度分布。
本发明具体包括以下步骤 (1)数据获取 在多光谱下，用滤波片将荧光波长λ离散为m个波段τ1，...τm，其中τl＝[λl-1，λl)，l＝1，2，...m-1，τl＝[λm-1，λm]，在多光谱情况下m＝5即可提供足够的先验信息；利用CCD探测器在被测物体的表面
对逃出的光子进行探测，获得每个波段τl的光流密度Q(x，τl)，其中CCD探测到的表面为x表示在被测物体中的位置；根据下面公式计算被测物体表面的光子流密度Φmeas(x，τl) D(x，τl)＝1/(3(μa(x，τl)+(1-g)μs(x，τl)))是生物组织的扩散系数，其中μa(x，τl)是生物组织的吸收系数，μs(x，τl)是生物组织散射系数，g是生物组织各向异性参数；v(x)是被测物体边界
的单位法向量；A(x；n，n′)≈(1+R(x))/(1-R(x))，当外界媒质是空气时，R(x)可以近似为R(x)≈-1.4399n-2+0.7099n-1+0.6681+0.0636n，其中n为空气的折射率； (2)有限元离散化处理 光子在生物组织中传输的数学模型用下面的扩散方程表示 其中Ω是被测物体；
是被测物体的表面；Φ(x，τl)是光子流密度；S(x，τl)是光源密度；根据有限元理论，得到下面的弱解形式 H1(Ω)是Sobelev空间，Ψ(x，τl)是任意给定的试探函数。在自适应有限元分析的框架下，基于自适应的网格细分，把重建区域网格剖分表示为{T1，...，Tk，...}的网格序列，其中网格剖分Tk包括
个离散单元。利用有限元方法，对弱解形式进行离散，得到下面矩阵方程 (Kk(τl)+Ck(τl)+Bk(τl))Φk(τl)＝Fk(τl)Sk(τl) 矩阵元素的具体形式为 令Kk(τl)+Ck(τl)+Bk(τl)＝Mk(τl)，考虑到Mk(τl)是稀疏正定矩阵，所以得到 (1)通过删除 [Mk-1(τl)Fk(τl)]中对应非测量点的行得到下面的方程 Φk(τl)＝Ak(τl)Sk(τl) (2) 因此对于m个波段即得到m个单一波段的矩阵方程； (3)优化可行光源区域选择 I)在网格剖分T1上，将上述含有m个单一波段的矩阵方程组合成含有m个波段的一个方程 根据光源谱信息的性质，每个谱段τl上所占整个谱段的能量表示为S(τl)＝ω(τl)·S，其中ω(τl)≥0，ω(τl)可通过实验测定；S表示光源总的能量。在网格剖分T1上，考虑光源谱信息的性质，通过方程(2)得到 AS1＝Φ (3) 其中 Φ1(τl)和A1(τl)都有方程(2)计算得到，S1是网格剖分T1上的未知光源密度变量；在(3)式两边同乘AT，得到 ATAS1＝ATΦ(4) AT是A的转置矩阵，AT·A是一个
的对称矩阵，所以(4)式是一个标准线形方程。
II)在网格剖分T1上，确定优化光源可行区域 具体方法为在网格剖分T1上，假定整个被测物体Ω为光源可行区域，给定未知光源密度的初始值S10，根据公式(4)，网格剖分T1上的未知光源密度变量S1用下面的关系式进行迭代求解 其中步长搜索方向n是迭代次数； 在迭代的过程中，用Φ1meas(τl)代替Φ1(τl)，这是因为噪声在迭代的过程中产生的影响很小；Φ1meas(τl)是网格剖分T1上的边界测量值，由Φmeas(x，τl)通过最近邻域插值得到；当‖βn‖≤δ或n＞Nmax，迭代停止，此时便得到了未知光源密度S1的优化解S*，优化解S*的大致区域为Ω*，这个区域Ω*我们称之为优化的光源可行区域；在每次迭代的过程中，如果光源的密度小于零，则将其置为零；δ取值范围为10-5-10-2，Nmax为最大迭代次数，取值不超过1000。
(4)光源重建 在网格序列{T1，...，Tk，...}上，根据选择好的优化光源可行区域Ω*，重新将含有m个单一波段的矩阵方程组合成含有m个波段的一个方程 考虑未知光源密度Sk(τl)与边界测量值Φkmeas(τl)之间的线性关系，由方程(1)得到 其中Gk(τl)通过移出[Mk(τl)-1Fk(τl)]中对应非测量点的行而得到；考虑到光源谱信息和优化可行光源区域Ω*，得到 其中 Gk表示多谱下的边界测量值与光源位置量之间的关系，Wk为k层网格上的对角矩阵，Φkmeas(τl)是波段τl在网格剖分Tk上的边界测量值，由Φmeas(x，τl)通过最近邻域插值得到；当k＝1时，sk(i)是通过迭代方法得到的优化解S*，skmax是优化解S的最大值；当k≥2时，sk(i)和skmax是网格剖分Tk上的未知光源密度的初始值和初始值的最大值，通过插值网格剖分Tk-1上的重建结果获得，即Sk-1是在网格剖分Tk-1上重建得到的结果，Ik-1k是插值因子，通过子单元继承父单元的重建结果值的方式来实现；γk是比例因子，是随k变化的分段常数；最终通过保留GkWk中对应可行光源区域的列，得到光源可行区域的未知光源密度变量Skp与表面测量值Φkmeas在网格剖分Tk上的关系 II)利用Tikhonov正则化方法建立目标函数并优化从而得到重建结果 上一步虽然得到了未知光源密度变量Skp与表面测量值Φkmeas之间的关系，由于Ak矩阵的病态，不能通过直接求解的方法得到重建结果，因此利用Tikhonov正则的方法，并考虑未知光源密度变量不能为负值的特性，得到在网格剖分Tk上的目标函数 考虑到多光谱以及非匀质模型会极大的增加计算量，因此采用谱梯度大规模优化算法对该目标函数进行优化；Ssupk是光源密度上界，经验取值为不大于1000，‖V‖Λ＝VTΛV，λk正则化系数，经验取值10-9-10-6，Skinit是网格剖分Tk上的未知光源密度初始值，取值范围为10-7-10-4；在判断优化过程是否中止时，采用当前优化梯度
以及优化迭代次数Nki作为评判准则，即当或时，优化停止，得到重建结果，即光源可行区域的光源密度分布；其中εg和Nkmax分别是梯度阈值和最大优化迭代次数，可根据实验经验获得，取值不超过10-6和104； III)当优化完成后，判断重建是否停止 利用优化求得的光源密度求解任意一个谱段τl上的光密度分布Φkc(τl)，并计算‖Φkc(τl)-Φkmeas(τl)‖，当或k≥Kmax，重建停止；Φkc(τl)是通过将优化求得的光源密度代入扩散方程而求得的边界处光子流密度，Φkmeas(τl)是波段τl在网格剖分Tk上边界测量点的测量值，由Φmeas(x，τl)通过最近邻域插值得到；εΦ经验取值范围10-9-10-7；Kmax是网格最大剖分次数；如果不能满足或k≥Kmax，根据重建的结果分别用直接最大值选择方法和等级值亏修正基的误差估计方法选取可行和禁止光源区域；通过以上方法对可行和禁止光源区域内的网格单元的选择，对区域剖分网格进行自适应的细分；并从网格剖分Tk转到Tk+1，然后转到第(4)(I)步，重复前面的步骤，直至重建停止； 本发明基于单视图测量数据，利用多光谱提供的先验信息，通过迭代方法进行可行光源区域的选择，重建光源的不准确问题得到了很好的解决，重建的速度和质量由于光源可行区域的使用而被进一步提高，此方法对自发荧光断层成像的发展具有重要的应用价值。
图1.非均匀仿体 图中(a)非均匀仿体，1代表肌肉、2骨骼、3心脏、4肺、5肝和 6单光源；(b)重建算法用到的初始网格；(c)包含肌肉、骨骼、肺的横截面，光源在右肺中；箭头所示的方向为CCD探测的方向，7前视图方向，8左视图方向，9后视图方向，10右视图方向。
图2.程序主流程图。
图3.确立光源可行区域的子流程图。
图4.确定光源可行区域的算法示意图。
图5.重建结果示意图。
表1.各组织在不同波长下的光学特性参数。
具体实施例方式 下面结合附图详细描述本实施例。
(1)数据 获取 为了验证本方法，我们设计了非匀质模型来近似模仿小鼠腹部内的组织，如图1(a)所示，生物组织分别为肌肉、骨骼、心脏、肺、肝，各个组织的光学特性参数如表1所示。在此实验中，我们将重建自发光源的光源密度分布作为重建目标，在重建实验中，这个光源大小为半径1毫米，光源密度为0.238nano-Watts/mm3，位置为(-3，5，15)。它的波谱范围从500纳米到750纳米，为了进行多光谱实验，我们依据当前的探测水平，把这一波谱范围分成了5个谱段进行分别探测，分别为τ1＝[500，550)nm，τ2＝[550，600)nm，τ3＝[600，650)nm，τ4＝[650，700)nm，τ5＝[700，750)nm。边界测量数据由蒙特卡罗方法产生，并按照图1(c)所示的方法取前视图的边界数据，并求得每个波段的光子流密度Φmeas(x，τl)； (2)有限 元离散化处理 光子在生物组织中传输的数学模型用下面的扩散方程表示 其中Ω是被测物体；
是被测物体的表面；Φ(x，τl)是光子流密度；S(x，τl)是光源密度；根据有限元理论，得到下面的弱解形式 H1(Ω)是Sobelev空间，Ψ(x，τl)是任意的试探函数。在自适应有限元分析的框架下，基于自适应的网格细分，把重建区域网格剖分表示为{T1，...，Tk，...}的网格序列，其中网格剖分Tk包括
个离散单元，网格剖分T1包含6878个离散单元，如图1(b)所示。利用有限元方法，对弱解形式进行离散，在每个波段τl(l＝1，2，3，4，5)上得到下面矩阵方程 (Kk(τl)+Ck(τl)+Bk(τl))Φk(τl)＝Fk(τl)Sk(τl) 矩阵元素的具体形式为 令Kk(τl)+Ck(τl)+Bk(τl)＝Mk(τl)，考虑到Mk(τl)是稀疏正定矩阵，所以得到 (1)通过删除 [Mk-1(τl)Fk(τl)]中对应非测量点的行得到下面的方程 Φk(τl)＝Ak(τl)Sk(τl) (2) 对于5个波段得到5个单一波段的矩阵方程； (3)优化可行光源区域选择 I)在网格剖分T1上，将上述含有5个单一波段的矩阵方程组合成含有5个波段的一个方程 根据光源谱信息的性质，每个谱段τl上所占整个谱段的能量表示为S(τl)＝ω(τl)·S，ω(τl)≥0，其中S表示光源总的光源能量，值为1；ω(τl)通过实验测定分别为ω(τ1)＝0.222，ω(τ2)＝0.167，ω(τ3)＝0.222，ω(τ4)＝0.167，ω(τ5)＝0.222。
在网格剖分T1上，考虑光源谱信息的性质，通过方程(2)得到 AS1＝Φ (3) 其中 Φ1(τl)和A1(τl)都有方程(2)计算得到，S1是网格剖分T1上的未知光源密度变量。在(3)式两边同乘AT，得到 ATAS1＝ATΦ (4) AT是A的转置矩阵，AT·A是一个34390×34390的对称矩阵，所以(4)式是一个标准线形方程。
III)在网格剖分T1上，确定优化光源可行区域 具体方法为在网格剖分T1上，假定整个被测物体Ω为光源可行区域，给定未知光源密度的初始值根据公式(4)，网格剖分T1上的未知光源密度变量S1根据下面的关系式进行迭代求解 其中步长搜索方向n是迭代次数。
在迭代的过程中，用Φ1meas(τl)代替Φ1(τl)，Φ1meas(τl)是τl在网格剖分T1上的边界测量值，由Φmeas(x，τl)通过最近邻域插值得到；参数取值分别为δ＝5×10-3，Nmax＝500。经过5次迭代后，‖βn‖＝4.9e-5，‖βn‖≤δ，所以迭代停止，此时便得到了未知光源密度S1的优化解S*，优化解S*的大致区域为Ω*，这个区域Ω*我们称之为优化的光源可行区域。在每次迭代的过程中，如果光源的密度小于零，则将其置为零。
(4)光源重建 (I)在被测物体Ω的网格剖分{T1，...，Tk，...}的网格序列上，根据选择好的优化光源可行区域Ω*，重新将含有5个波段的矩阵方程组合成含有5个波段的一个方程 考虑未知光源密度Sk(τl)与边界测量值Φkmeas(τl)之间的线性关系，由方程(1)得到 其中Gk(τl)通过移出[Mk(τl)-1Fk(τl)]中对应非测量点的行而得到。考虑到光源谱信息和优化可行光源区域Ω*，得到 其中 Gk表示多谱下的边界测量值与光源位置量之间的关系，Wk为k层网格上的对角矩阵，Φkmeas(τl)是波段τl在网格剖分Tk上的边界测量值。当k＝1时，sk(i)是通过迭代方法得到的优化解S*，skmax是优化解S*的最大值；当k≥2时，sk(i)和skmax是网格剖分Tk上的初始值和初始值的最大值，通过插值网格剖分Tk-1上的重建结果得到，即Sk-1是在网格剖分Tk-1上重建得到的结果，Ik-1k是插值因子，通过子单元继承父单元的重建结果值的方式来实现；γk是比例因子，是随k变化的分段常数；k＝1时，参数γk＝0；k＞1时，γk初始为0.1，并随着k增加，有γk+1＝10γk。最终通过保留GkWk中对应可行光源区域的列，得到光源可行区域的未知光源密度变量Skp与表面测量值Φkmeas在网格剖分Tk上的关系 (II)利用Tikhonov正则化方法建立目标函数并优化从而得到重建结果 上一步虽然得到了未知光源密度变量Skp与表面测量值Φkmeas之间的关系，由于Ak矩阵的病态，不能通过直接求解的方法得到重建结果，利用Tikhonov正则的方法，并考虑未知光源密度变量不能为负值的特性，得到在网格剖分Tk上的目标函数 考虑到多光谱以及非匀质模型会极大的增加计算量，因此采用谱梯度大规模优化算法对该目标函数进行优化。Ssupk是光源密度上界，取值为100，‖V‖Λ＝VTΛV，λk正则化系数，取值1.0×10-8，Skinit是在网格剖分Tk上未知光源密度变量的初始值，取值1.0×10-5。在判断优化过程是否中止时，采用了当前优化梯度
以及优化迭代次数Nki作为评判准则，即当或时，优化停止，得到重建结果，即光源可行区域的光源密度分布，其中εg是梯度阈值，取值1.0×10-8，Nkmax是最大优化迭代次数，取值1.0×104。
(III)当优化完成后，判断重建是否停止 利用优化得到的光源密度分布求解任意一个谱段τl上的光子流密度Φkc(τl)，并计算‖Φkc(τl)-Φkmeas(τl)‖，当(εΦ为正数)或k≥Kmax(Kmax是网格最大剖分次数)，重建停止。Φkc(τl)是通过将优化求得的光源密度代入扩散方程而求得的边界处光子流密度，Φkmeas(τl)是波段τl在第k层网格上边界测量点的测量值；如果不能满足或k≥Kmax，利用重建的结果分别用直接最大值选择方法和等级值亏修正基的误差估计方法选取可行和禁止光源区域。通过以上方法对可行和禁止光源区域内的网格单元的选择，并对区域剖分网格进行自适应的细分，从网格剖分Tk转到Tk+1，然后转到第(4)(I)步，重复前面的步骤，直至重建停止。
停止阈值εΦ和网格最大细分次数Kmax分别是1×10-8和3。当在网格剖分T1上重建完成后，计算在τ1＝[500-550)nm的光密度分布Φkc(τl)，求得此时(k＝1)＜(Kmax＝3)所以重建未停止；利用网格剖分T1上的重建结果分别用直接最大值选择方法和等级值亏修正基的误差估计方法选取可行和禁止光源区域。通过以上方法对可行和禁止光源区域内的网格单元的选择，并对区域剖分网格进行自适应的细分，并从网格剖分T1转到网格剖分T2，然后转到第(4)(I)步，重复前面的步骤；优化完成后计算在τ1＝[500-550)nm的光密度分布Φkc(τl)，求得因此重建停止，重建的结果如图5所示，重建的位置为(-1.53，5.3，14.98)，此时重建的最大光源密度为0.224nano-Watts/mm3。
表1
权利要求
1.一种基于单视图的多光谱自发荧光断层成像重建方法，其特征在于，包括以下步骤
1)数据获取
在多光谱下，用滤波片将荧光波长λ离散为m个波段τ1，...，τm，其中τl＝[λl-1，λl)，l＝1，2，...，m-1，τl＝[λm-1，λm]；利用CCD探测器在被测物体的表面
对逃出的光子进行探测，获得每个波段τl的光流密度Q(x，τl)，其中CCD探测到的表面为γ，x表示在被测物体中的位置；根据下面公式计算被测物体表面的光子流密度Φmeas(x，τl)
D(x，τl)＝1/(3(μa(x，τl)+(1-g)μs(x，τl)))是生物组织的扩散系数，其中μa(x，τl)是生物组织的吸收系数，μs(x，τl)是生物组织散射系数，g是生物组织各向异性参数；v(x)是被测物体边界
的单位法向量；A(x；n，n′)≈(1+R(x))/(1-R(x))，当外界媒质是空气时，R(x)近似为R(x)≈-1.4399n-2+0.7099n-1+0.6681+0.0636n，其中n为空气的折射率；
2)有限元离散化处理
光子在生物组织中传输的数学模型用下面的扩散方程表示
其中Ω是被测物体；
是被测物体的表面；Φ(x，τl)是光子流密度；S(x，τl)是光源密度；根据有限元理论，得到下面的弱解形式
H1(Ω)是Sobelev空间，ψ(x，τl)是任意给定的试探函数；在自适应有限元分析的框架下，基于自适应的网格细分，把重建区域网格剖分表示为{T1，...，Tk，...}的网格序列，其中网格剖分Tk包括
个离散单元；利用有限元方法，对弱解形式进行离散，得到下面矩阵方程
(Kk(τl)+Ck(τl)+Bk(τl))Φk(τl)＝Fk(τl)Sk(τl)
矩阵元素的具体形式为
令Kk(τl)+Ck(τl)+Bk(τl)＝Mk(τl)，考虑到Mk(τl)是稀疏正定矩阵，所以得到
(1)通过删除
[Mk-1(τl)Fk(τl)]中对应非测量点的行得到下面的方程
Φk(τl)＝Ak(τl)Sk(τl)
因此对于m个波段即得到m个单一波段的矩阵方程；
3)优化可行光源区域选择
I)在网格剖分T1上，将上述含有m个单一波段的矩阵方程组合成含有m个波段的一个方程
根据光源谱信息的性质，每个谱段τl上所占整个谱段的能量表示为S(τl)＝ω(τl)·S，其中ω(τl)≥0，ω(τl)通过实验测定；S表示光源总的能量；在网格剖分T1上，考虑光源谱信息的性质，通过方程(2)得到
AS1＝Φ (3)
其中
Φ1(τl)和A1(τl)都有方程(2)计算得到，S1是网格剖分T1上的未知光源密度变量；在(3)式两边同乘AT，得到
ATAS1＝ATΦ (4)
AT是A的转置矩阵，AT·A是一个
的对称矩阵，所以(4)式是一个标准线形方程；
II)在网格剖分T1上，确定优化光源可行区域
具体方法为在网格剖分T1上，假定整个被测物体Ω为光源可行区域，给定未知光源密度的初始值S10，根据公式(4)，网格剖分T1上的未知光源密度变量S1用下面的关系式进行迭代求解
其中步长搜索方向n是迭代次数；
在迭代的过程中，用Φ1meas(τl)代替Φ1(τl)，Φ1meas(τl)是网格剖分T1上的边界测量值，由Φmeas(x，τl)通过最近邻域插值得到；当||βn||≤δ或n＞Nmax，迭代停止，此时便得到了未知光源密度S1的优化解S*，优化解S*的大致区域为Ω*，这个区域Ω*我们称之为优化的光源可行区域；在每次迭代的过程中，如果光源的密度小于零，则将其置为零；δ取值范围为10-5-10-2，Nmax为最大迭代次数，取值不超过1000；
4)光源重建
I)在网格序列{T1，...，Tk，...}上，根据选择好的优化光源可行区域Ω*，重新将含有m个单一波段的矩阵方程组合成含有m个波段的一个方程
考虑未知光源密度Sk(τl)与边界测量值Φkmeas(τl)之间的线性关系，由方程(1)得到
其中Gk(τl)通过移出[Mk(τl)-1Fk(τl)]中对应非测量点的行而得到；考虑到光源谱信息和优化可行光源区域Ω*，得到
其中
Gk表示多谱下的边界测量值与光源位置量之间的关系，Wk为k层网格上的对角矩阵，Φkmeas(τl)是波段τl在网格剖分Tk上的边界测量值，由Φmeas(x，τl)通过最近邻域插值得到；当k＝1时，sk(i)是通过迭代方法得到的优化解S*，skmax是优化解S*的最大值；当k≥2时，sk(i)和skmax是网格剖分Tk上的未知光源密度的初始值和初始值的最大值，通过插值网格剖分Tk-1上的重建结果获得，即Sk-1是在网格剖分Tk-1上重建得到的结果，Ik-1k是插值因子，通过子单元继承父单元的重建结果值的方式来实现；γk是比例因子，是随k变化的分段常数；最终通过保留GkWk中对应可行光源区域的列，得到光源可行区域的未知光源密度变量Skp与表面测量值Φkmeas在网格剖分Tk上的关系
II)利用Tikhonov正则化方法建立目标函数并优化从而得到重建结果
上一步虽然得到了未知光源密度变量Skp与表面测量值Φkmeas之间的关系，由于Ak矩阵的病态，不能通过直接求解的方法得到重建结果，因此利用Tikhonov正则的方法，并考虑未知光源密度变量不能为负值的特性，得到在网格剖分Tk上的目标函数
考虑到多光谱以及非匀质模型会极大的增加计算量，因此采用谱梯度大规模优化算法对该目标函数进行优化；Ssupk是光源密度上界，经验取值为不大于1000，||V||Λ＝VTΛV，λk正则化系数，经验取值10-9-10-6，Skinit是网格剖分Tk上的未知光源密度初始值，取值范围为10-7-10-4；在判断优化过程是否中止时，采用当前优化梯度
以及优化迭代次数Nki作为评判准则，即当或时，优化停止，得到重建结果，即光源可行区域的光源密度分布；其中εg和Nkmax分别是梯度阈值和最大优化迭代次数，根据实验经验获得，取值不超过10-6和104；
III)当优化完成后，判断重建是否停止
利用优化得到的光源密度求解任意一个谱段τl上的光密度分布Φkc(τl)，并计算||Φkc(τl)-Φkmeas(τl)||，当或k≥Kmax，重建停止；Φkc(τl)是通过将优化求得的光源密度代入扩散方程而求得的边界处光子流密度，Φkmeas(τl)是波段τl在网格剖分Tk上边界测量点的测量值，由Φmeas(x，τl)通过最近邻域插值得到；εΦ经验取值范围10-7-10-9；Kmax是网格最大剖分次数；如果不能满足或k≥Kmax，根据重建的结果分别用直接最大值选择方法和等级值亏修正基的误差估计方法选取可行和禁止光源区域；通过以上方法对可行和禁止光源区域内的网格单元的选择，对区域剖分网格进行自适应的细分；并从网格剖分Tk转到Tk+1，然后转到第4)(I)步，直至重建停止。
全文摘要
基于单视图的多光谱自发荧光断层成像重建方法属于光学分子影像成像领域。该方法基于扩散方程模型，并考虑了小动物体的非均匀特性，同时也考虑自发荧光光源的光谱及实际应用的特点。为达到此目的，本发明提出的基于单视图的多光谱自发荧光断层成像重建方法的实现步骤主要包括1)数据获取；2)有限元离散化处理；3)优化可行光源区域选择；4)光源的重建。本发明克服了自发荧光断层成像重建方法重建光源不准确，重建速度比较慢，不利于实时处理以及在多光谱数据采集时可能存在误差的缺陷。
文档编号A61B5/00GK101342075SQ200810116818
公开日2009年1月14日 申请日期2008年7月18日 优先权日2008年7月18日
发明者贾克斌, 冯金超, 捷 田 申请人:北京工业大学