【摘要】 基于单视图的多光谱自发荧光断层成像重建方法属于光学分子影像成像领域。该方法基于扩散方程模型，并考虑了小动物体的非均匀特性，同时也考虑自发荧光光源的光谱及实际应用的特点。为达到此目的，本发明提出的基于单视图的多光谱自发荧光断层成像重建方法的实现步骤主要包括：1)数据获取；2)有限元离散化处理；3)优化可行光源区域选择；4)光源的重建。本发明克服了自发荧光断层成像重建方法重建光源不准确，重建速度比较慢，不利于实时处理以及在多光谱数据采集时可能存在误差的缺陷。 【专利类型】发明授权 【申请人】北京工业大学 【申请人类型】学校 【申请人地址】100124 北京市朝阳区平乐园100号 【申请人地区】中国 【申请人城市】北京市 【申请人区县】朝阳区 【申请号】CN200810116818.4 【申请日】2008-07-18 【申请年份】2008 【公开公告号】CN101342075B 【公开公告日】2010-06-02 【公开公告年份】2010 【授权公告号】CN101342075B 【授权公告日】2010-06-02 【授权公告年份】2010.0 【IPC分类号】A61B5/00 【发明人】贾克斌; 冯金超; 田捷 【主权项内容】1.一种基于单视图的多光谱自发荧光断层成像重建方法，其特征在于，包括以下步骤： 1)数据获取： 在多光谱下，用滤波片将荧光波长λ离散为m个波段τ1，...，τm，其中τl＝[λl-1，λl]，l＝1，2，...，m-1，τm＝[λ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))，n′为外界媒质的折射系数；当外界媒质是空气时，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)是稀疏正定矩阵，所以得到： 通过删除 [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是网格剖分τ1上的未知光源密度变量；在(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)]中对应非测量点的行而得到；考虑到光源谱信息和优化可行光源区域Ω*，得到： 其中 Wk＝Diag(wk(11)，wk(22)，…wk(ii)，…)， Gk表示多谱下的边界测量值与光源位置量之间的关系，Wk为k层网格上的对角矩阵，wk(ii)是对角矩阵Wk的第i个对角元素；Φ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；在判断优化过程是否中止时，采用当前优化梯度 以及优化迭代次数Nkte作为评判准则，即当 或 时，优化停止，得到重建结果，即光源可行区域的光源密度分布；上标te表示当前的迭代次数；εg是梯度阈值，其经验取值不超过10-6；Nkmax是最大优化迭代次数，经验取值不超过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)步，直至重建停止。 【当前权利人】北京工业大学 【当前专利权人地址】北京市朝阳区平乐园100号 【专利权人类型】公立 【统一社会信用代码】12110000400687411U 【引证次数】6.0 【被引证次数】8 【他引次数】6.0 【被他引次数】8.0 【家族引证次数】6.0 【家族被引证次数】32
