[发明专利]一种基于线性关系的荧光分子断层成像重建方法

摘要

一种基于线性关系的荧光分子断层成像重建方法属于在生物医学中图像重建的技术领域。其特征在于，依次含有以下步骤：一、实验用鼠图谱数据通过图像处理和数值拟合得到外表面轮廓；二、这个外表面轮廓被导入有限元软件包进行网格划分；三、找出每个激发光源以及相应的荧光检测点对应的有限元网格节点；四、针对一个特定的激发光源sk，建立表面检测到的荧光强度同未知荧光探针分布之间的线性关系矩阵W_sk；五、针对所有激发光源，生成总的线性关系矩阵；六、得到这个线性关系后，代数重建方法被用来迭代逼近真实的荧光探针分布。本发明具有可处理非均匀光学参数以及快速建立表面检测到的荧光强度同未知荧光探针分布之间的精确线性关系的特点。

主权项

1.一种基于线性关系的荧光分子断层重建方法，其特征在于，所述方法是用一台计算机依次按以下步骤实现的：步骤(1).把用CT成像设备分割完毕的生物实验用鼠的分层的图谱输入所述计算机中，该计算机通过下述图像处理方法得到外表面轮廓，并把所表述表面轮廓建模成几何体：步骤(1.1).对所述图谱的每一层图像I，用一个相同的阈值thr＝0.5进行二值化：大于所选阈值thr者取1，小于或等于该阈值thr者取0，把该图谱中的组织和背景依次用所述的“1”和“0”分开，得到每一层图像的二值图BI：BI[row,col]=1,I[row,col]>thr,0,I[row,col)]≤thr,]]>其中，row代表该图像I的行，col代表该图像I的列，步骤(1.2).用二值图轮廓线提取的方法得到所述的二值图BI的外轮廓线S：对该二值图BI上的每一个像素点BI[row，col]执行以下操作：若：BI[row，col]为1，而且其相邻的8个像素点至少有一个为0，则该像素点BI[row，col]在所述外轮廓线S上，若：BI[row，col]为0，而且其相邻的8个像素点至少有一个为1，则该像素点BI[row，col]在所述外轮廓线S上，步骤(1.3).把步骤(1.2)得到的所述外轮廓线S从像素坐标对应到笛卡尔坐标，则该外轮廓线S上的任意一个点[row，col]对应到所述笛卡尔坐标[x，y]上为，x＝(row-row₀)*Δ，y＝(col-col₀)*Δ，其中，点[row₀，col₀]为该轮廓线的中心，Δ是每个像素点的几何尺寸，用cm表示，步骤(1.4).把步骤(1.3)得到的该笛卡尔坐标[x，y]对应到极坐标[α，R]上，α是角度，R是极坐标的半径，对所述外轮廓线S上的所有点，把每个点上的所述半径R用角度α按下式做多项式拟合：Rn=Σc=0c=12pcαnc,]]>c为拟合次数，c＝0，1，…12，p_c是拟合得到的多项式系数，用{p_c}表示，n＝1，2，…S_num，S_num是所述轮廓线S上的所有点的数目，在[0，2π]均匀选定α_num＝30个角度{α₁，…，α_{α_num}}，用多项式{pc}拟合这些角度上的半径{R₁，…，R_{α_num}}，使得所述每一层图像的所述每一层轮廓线S表示为：{{α₁，…，α_{α_num}}，{R₁，…，R_{α_num}}，{p₀，…，p₁₂}}，用{S}表示所有层图像上I上得到的轮廓线，步骤(1.5)在所述所有层图像I的轮廓线{S}中，按照高度均匀选定M层，以此表示所述所有层图像I的外轮廓线，则相邻两层H₁<H₂之间，用下式表示任意一层轮廓线的坐标，用半径R表示为：R1=PH1(α),]]>R2=PH2(α),]]>H₁<H₂，R＝R₁*(H₂-H)/(H₂-H₁)+R₂(H-H₁)/(H₂-H₁)，其中，R₁，R₂分别为在所述高度H₁，H₂层角度α所对应的半径，P_H1，分别是在高度H₁，H₂层所对应的12次多项式，步骤(1.6).把步骤(1.4)～(1.5)得到的三维外轮廓线从极坐标[α，R，H]转换到三维笛卡尔坐标系[x，y，z]，得到了所述三维笛卡尔坐标系下的几何模型G；步骤(2).把要进行活体荧光分子断层成像的一个生物实验用鼠放置在一个旋转平台上，激光器发出的激发光从一侧照射所述活体生物实验用鼠，激发光激发出荧光，被激发出的荧光在相对一侧被CCD相机检测，再把检测结果输入计算机中，分角度通过使所述旋转平台旋转，采集不同旋转角度被激发的荧光图像，相当于实现了在不同角度时的激发光点光源激发；步骤(3).把步骤(1)得到的几何模型G导入一个有限元软件包中，同时把所述各个激发光光源所对应的点加入所述几何模型的节点上，再用有限元软件包做网格划分，离散成一个四面体网格Fem_Mesh，该四面体网格的所有节点自动编号为i＝1，…，N，N是总的网格节点数；步骤(4).按以下步骤确定每个激发光点光源s所对应的在所述四面体网格Fem_Mesh上的有限元网格节点：步骤(4.1).把垂直入射到所述活体生物实验用鼠的激光束建模成入射光皮下一个自由传输自由程ltr=1/(3μse′(rs))]]>的位置r_s处的激发光点光源δ(r-r_s)，该δ(r-r_s)代表在位置r_s处的一个冲激函数，是r_s处激发光波长的散射系数，步骤(4.2).对于每个所述四面体网格Fem_Mesh上有限元网格节点i，若r_s＝r_i，则该有限元网格节点i就是激发光点光源s所对应的有限元网格节点；步骤(5).确定所述每个激发光点光源s在所述Fem_Mesh上有限元网格表面上的检测点所对应的所述四面体网格Fem_Mesh上的有限元网格节点i：若：所述表面节点i满足以下关系： r_i∈{angle_low～angle_high}且r_i∈{height_low～height_high}，所述  angle_low，angle_high，height_low，height_high 为设定值，则所述节点i属于所述激发光点光源s对应的网格节点；步骤(6).在得到步骤(5)所述各检测点后，利用所述CCD相机检测得到的荧光图像，得到各检测点发出的荧光图像；步骤(7).针对所述给定位置r_s处的激发光点光源s，按照以下步骤计算对应所述激发光点光源s所述表面检测点{d1，d2，...，dL}上检测到的荧光强度同未知的荧光强度之间的线性关系，其中L是检测点总数：步骤(7.1).按下式计算N×N维刚度矩阵K_e中的元素K_e[i，j]，该刚度矩阵K_e是生物组织在激发光波长在所述四面体网格Fem_Mesh上的有限元节点上的吸收和散射系数的体现，该刚度矩阵K_e的行和列分别用i和j表示，所述的行i对应所述Fem_Mesh上的所有节点N按照编号顺序排列而成的行，所述的列j对应所述Fem_Mesh上的所有节点N按照编号顺序排列而成的列，下同，其中任意一元素K_e[i，j]表示为：Ke[i,j]=∫Ω(De(r)▿ψi(r)·▿ψj(r)+μae(r)ψi(r)ψj(r))dr+12q∫∂Ωψi(r)ψj(r)dr,]]>其中，Ω代表所述的几何模型G内部对应的空间区域，代表所述的几何模型G的表面边界，D_e(r)是位于r位置处激发光波长的漫射系数，所述De(r)=1/(3μse′(r)),]]>为位于r位置处的组织在所述激发光波长的散射系数，为已知值，μ_ae(r)为位于r位置处的组织在所述激发光波长的吸收系数，为已知值，ψ_i(r)和ψ_j(r)是行编号为i，列编号为j的有限元节点对应的Lagrange线性形函数，由所述有限元软件包给出，▿ψi(r)=[∂ψi(r)/∂x,∂ψi(r)/∂y,∂ψi(r)/∂z]]]>和▿ψj(r)=[∂ψj(r)/∂x,∂ψj(r)/∂y,∂ψj(r)/∂z]]]>分别是所述Lagrange线性形函数ψ_i(r)、ψ_j(r)的梯度，q＝0.1511，步骤(7.2).使用Matlab软件包函数PCG按下式计算所述激发光点光源s在所述有限元网格Fem_Mesh节点上的强度分布向量Φ_e，Φ_e是N×1列向量：K_eΦ_e＝B_s，N×1列向量B_s的列编号为j的元素B_s[j]用下式表示：已知了在节点上的强度分布向量Φ_e，得到任意位置r处的强度分步Φe(r)=Σj=1NΦe[j]]]>步骤(7.3).按下式计算N×N矩阵F_s的元素F_s[i，j]，该矩阵F_s是所述激发光强度Φ_e(r)对荧光的激发作用在所述四面体网格Fem_Mesh上的有限元节点上的体现，其中：F_s[i，j]＝∫_ΩΦ_e(r)ψ_i(r)ψ_j(r)dr，步骤(7.4).按下式计算N×N维刚度矩阵K_m中的元素K_m[i，j]，该刚度矩阵K_m表示生物组织在激发光波长在所述四面体网格Fem_Mesh上的有限元节点上的吸收和散射系数的体现，其中：Km[i,j]=∫Ω(Dm(r)▿ψi(r)·▿ψj(r)+μam(r)ψi(r)ψj(r))dr+12q∫∂Ωψi(r)ψj(r)dr,]]>其中，D_m(r)是位于r位置处荧光波长的漫射系数，所述Dm(r)=1/(3μsm′(r)),]]>为位于r位置处的组织在荧光波长的散射系数，为已知值，μ_am(r)为位于r位置处的组织在荧光波长的吸收系数，为已知值，步骤(7.5).假想当把一个荧光波长的点光源放在所述检测点dl处，dl∈{d1，d2，…，dL}，使用Matlab软件包函数PCG计算所产生的荧光场在所有所述四面体网格Fem_Mesh上的有限元节点上的强度分布向量Φ_dl，向量Φ_dl是N×1列向量，其中：K_mΦ_dl＝B_dl，N×1列向量B_dl的列编号为j的元素B_dl[j]用下式表示：步骤(7.6).按照下式生成所述各检测点{d1，d2，…，dL}处检测到的荧光强度同未知的荧光探针分布之间的线性关系，用矩阵W_s表示：Ws=Φd1···ΦdLFs;]]>步骤(8).针对空间顺序分布的激发光点光源sk，k＝1，…，N_s，以及相对应的所述检测点，按照下式构建总的线性关系矩阵W，下述Φ_m，meas是各激发光点光源对应的被检测荧光，为已知值：Φm,meas=Ws1···WsNsU=WU,]]>N×1列向量U是未知的荧光探针在所述有限元网格Fem_Mesh节点上的分布，Φ_m，meas是M×1列向量，M所有激发光点光源对应的总的检测点数目，W是M×N矩阵；步骤(9).按下述代数表达式，根据步骤(8)得到的总的线性关系矩阵W求解位置的荧光探针分布U：U[j]u_iter+1=U[j]u_iter+λΣt1=1MV[t1]-Σt2=1NW[t1,t2]U[t2]u_iterΣt2=1NW[t1,t2]W[t1,t2]W[t1,j],]]>其中：U[j]^u_iter+1、U[j]^u_iter分别是U中第j个元素在第u_iter+1和u_iter次迭代时的值，U[t₂]^u_iter是U中第t₂个元素在第u_iter次迭代时的值，W[t₁，t₂]是矩阵W第t₁行，第t₂列的值，W[t₁，j]是矩阵W第t₁行，第j列的值，Φ_m，meas用V表示，V[t₁]是的V第t₁个元素，λ＝0.1，U迭代初始值设为U＝0。