一种基于线性关系的荧光分子断层成像重建方法

专利名称：：一种基于线性关系的荧光分子断层成像重建方法
技术领域：
：本发明属于在生物医学中图像重建的
技术领域：
，涉及一种三维图像重建方法。技术背景分子成像是近年來生物医学领域的一个热门方向。分子成像，广义上可以定义为在分子和细胞水平对活体状态下的生物过程进行的定性和定量的测量[R.Weissleder，Radiology,2001]。荧光分子断层成像(FMT)是分子成像中的一个重要方面。它被广泛的应用在生物，医学研究领域，使得荧光分子成像已经从一种定性的用作切片检査的工具变成了一种能够用于准确的，3维成像的工具。荧光分子断层成像的关键在于获得多个角度的投影数据用于重建。同CT等基于高能射线的设备有所不同，光(400nm900nm)在生物组织中的传播不是遵循直线传播。这个光传播路径可以近似用漫射方程来描述。利用基于漫射方程模型的重建方法，从多个角度投影数据就可以断层地重建出荧光探针的浓度分布。R.B.Schulz等人[R.B.Schulz，正EETransactionsonMedicalImaging,2004]提出了一个基于解析解的重建方法，能够得到表面检测到荧光强度同内部荧光体浓度分布之间的线性关系。这种方法假设成像物体的光学参数是均匀的。然而，实际成像的实验用鼠的光学参数是非均匀的。因而，光学参数估计的不准确性会导致重建结果的很大失真。S.Bjoem等人[S.Bjoem，AdvancesinMedicalEngineering,2007]对光学参数估计不准确造成失真做了研究。白净教授领导的医学成像实验室在2007年发表了基于有限元方法的重建方法[X.Song，OpticsExpress,2007]，能够处理非均匀的光学参数，同时在方法中计算得到表面检测到的荧光强度同未知的荧光探针分布之间的线性关系，进而进行重建。这些方法具有准确，在迭代过程中不需要反复求解梯度以及二届梯度的特点，因而计算速度以5-10分钟计。然而，在计算线性关系的过程中，有针对大规模稀疏矩阵求逆的操作，限制了方法的计算速度地进一步提高
发明内容本发明的目的是为荧光分子断层成像系统提供一种重建方法，其特点是能够快速准确地建立实验用鼠表面检测到的荧光强度同未知的荧光探针分布之间的线性关系，实现重建。在这个过程中，方法能够处理不规则的成像物体以及非均匀的光学参数。方法中在其它成像模式，如CT，得到的分割好的实验用鼠图谱的基础上，利用二值图像轮廓线提取以及数值拟合的方法得到实验用鼠的不规则外表面轮廓；在实验用鼠图谱上，通过査表的方法确定实验用鼠体内的非均匀光学参数，包括吸收和散射参数；针对每个激发光源-荧光检测对，利用一个向量-稀疏矩阵乘积的方法得到这个检测点处的荧光强度同未知的荧光探针分布之间的线性关系。本发明的特征在于，所述方法是用一台计算机依次按以下步骤实现的-.步骤(l).把用CT成像设备分割完毕的生物实验用鼠的分层的图谱输入所述计算机中，该计算机通过下述图像处理方法得到外表面轮廓，并把所表述表面轮廓建模成几何体步骤(l.l).对所述图谱的每--层图像/，用一个相同的阈值*=0.5进行二值化大于所选阈值Ar者取l，小于或等于该阈值Ar者取O，把该图谱中的组织和背景依次用所述的"1"和"0"分开，得到每一层图像的二值图BI:<formula>formulaseeoriginaldocumentpage8</formula>其中，row代表该图像/的行，co/代表该图像/的列，步骤(1.2).用二值图轮廓线提取的方法得到所述的二值图S/的外轮廓线S:对该二值图W上的每一个像素点5/[raw，co/]执行以下操作若5/[row,co/]为1,而且其相邻的8个像素点至少有一个为0，则该像素点万/[raw，co/]在所述外轮廓线S上，若S/[row,co/]为0，而且其相邻的8个像素点至少有一个为1，则该像素点B/[row,co/]在所述外轮廓线S上，步骤(1.3).把步骤(1.2)得到的所述外轮廓线S从像素坐标对应到笛卡尔坐标，则该外轮廓线S上的任意一个点[raw,co/]对应到所述笛卡尔坐标[je，W上为<formula>formulaseeoriginaldocumentpage8</formula>其中，点[row。，co/。]为该轮廓线的中心，A是每个像素点的几何尺寸，用cw表示，步骤(1.4).把步骤(1.3)得到的该笛卡尔坐标[x，y]对应到极坐标[cr,i]上，《是角度，i是极坐标的半径，对所述外轮廓线S上的所有点，把每个点上的所述半径i用角度a按下式做多项式拟合<formula>formulaseeoriginaldocumentpage9</formula>c为拟合次数，c=0，l,'.'12，A是拟合得到的多项式系数，用{八.}表示，"二l,2,…&m^，S—m^是所述轮廓线S上的所有点的数目，在均匀选定a—"臓=30个角度{,，...,}，用多项式OJ拟合这些角度上的半径{《,，&,}，使得所述每一层图像的所述每一层轮廓线S表示为""p…，"c"V,、U"，…'A2"，用表示所有层图像上/上得到的轮廓线，步骤(1.5).在所述所有层图像/的轮廓线{5}中，按照高度均匀选定M层，以此表示所述所有层图像/的外轮廓线，则相邻两层//,<//2之间，用下式表示任意一层轮廓线的坐标，用半径i表示为<formula>formulaseeoriginaldocumentpage9</formula>，其中，A,A分别为在所述高度/Zp/^层角度a所对应的半径，i^,/^分别是在高度巧,/^层所对应的12次多项式，步骤(1.6).把步骤(1.4)(1.5)得到的三维外轮廓线从极坐标[,凡//]转换到三维笛卡尔坐标系[x,乂z]，得到了所述三维笛卡尔坐标系下的几何模型G;步骤(2).把要进行活体荧光分子断层成像的一个生物实验用鼠放置在一个旋转平台上，激光器发出的激发光从一侧照射所述活体生物实验用鼠，激发光激发出荧光，被激发出的荧光在相对一侧被CCD相机检测，再把检测结果输入计算机中，分角度通过使所述旋转平台旋转，采集不同旋转角度被激发的荧光图像，相当于实现了在不同角度时的激发光点光源激发；步骤(3).把歩骤(1)得到的几何模型G导入一个有限元软件包中，同时把所述各个激发光光源所对应的点加入所述几何模型的节点上，再用有限元软件包做网格划分，离散成一个四面体网格Fem—ikfe^，该四面体网格的所有节点自动编号为!、1,…,W，iV是总的网格节点数；步骤(4).按以下步骤确定每个激发光点光源s所对应的在所述四面体网格i^m—M^/7上的有限元网格节点X说明书第4/14页步骤(4.1).把垂直入射到所述活体生物实验用鼠的激光束建模成入射光皮下一个自由传输自由程=1/(3/^(0)的位置r、.处的激发光点光源5(r-O，该5(r-O代表在位置rs,处的一个冲激函数，//,Ur、)是^处激发光波长的散射系数，歩骤(4.2).对于每个所述四面体网格Fem—7kfei上有限元网格节点/，若^,=。则该有限元网格节点/就是激发光点光源^所对应的有限元网格节点；步骤(5).确定所述每个激发光点光源s在所述Few—iWe^上有限元网格表面上的检测点所对应的所述四面体网格Fem—上的有限元网格节点/:若所述表面节点/满足以下关系^e{""g/e/<w""g7e—且a;e{/ze/g/z/1—/譜—，所述a"g7e—/，，"wg7e—/n.沐—/，,/e妙/1—/z妙为设定值，则所述节点/属于所述激发光点光源^对应的网格节点；步骤(6).在得到步骤(5)所述各检测点后，利用所述CCD相机检测得到的荧光图像，得到各检测点发出的荧光图像；步骤(7).针对所述给定位置^;处的激发光点光源5，按照以下歩骤计算对应所述激发光点光源s所述表面检测点W1^2，…，t/W上检测到的荧光强度同未知的荧光强度之间的线性关系，其中丄是检测点总数步骤(7.1).按下式计算7VxiV维刚度矩阵《中的元素《[/,y]，该刚度矩阵&是生物组织在激发光波长在所述四面体网格i^w一7kfe^上的有限元节点上的吸收和散射系数的体现，该刚度矩阵《的行和列分别用,'和J'表示，所述的行/对应所述Fem—7kfe^上的所有节点按照编号顺序排列而成的行，所述的列y'对应所述Fem一iWeA上的所有节点TV按照编号顺序排列而成的列，下同，其中任意一元素《[/，/]表示为&["刀=J"。(A(>)▽W,0)Vy乂(/)+//"e0》A+由(r)^(r)A，其中，Q代表所述的几何模型G内部对应的空间区域，3Q代表所述的几何模型G的表面边界，A(/)是位于H立置处激发光波长的漫射系数，所述马(0=1/(3&)，；^(r)为位于r位置处的组织在所述激发光波长的散射系数，为已知值，/^(r)为位于r位置处的组织在所述激发光波长的吸收系数，为已知值，V,(。和^(r)是行编号为"列编号为y的有限元节点对应的Lagrange线性形函数，由所述有限元软件包给出，VV,(r)=[3(r)/3y/,0)/,d(r)/3z]和Vy,(/)=[30)/5y,(/)/*，3y/,0)/&]分别是所述Lagrange线性形函数y,O)、^.(r)的梯度，《=0.1511，步骤(7.2).使用Matlab软件包函数PCG按下式计算所述激发光点光源s在所述有限元网格Few—M;W节点上的强度分布向量^，(^是7Vxl列向《Ot'=S、，TVx1列向量5、的列编号为y的元素5,[y]用下式表示,fl,有限元节点_/位于所述;;.位置上，已知了在节点上的强度分布向量O,，得到任意位置^处的强度分步0￡(0=1]^[_/]步骤(7.3).按下式计算7VxN矩阵《的元素《[/,刀，该矩阵《是所述激发光强度0)￡(0对荧光的激发作用在所述四面体网格Few—MeW上的有限元节点上的体现，其中步骤(7.4).按下式计算WxiV维刚度矩阵i^中的元素i^[/,/1，该刚度矩阵i^表示生物组织在激发光波长在所述四面体网格Fe附—Me^上的有限元节点上的吸收和散射系数的体现，其中l['.，/]=J"。(化(。▽K0)'▽^0)+//。m(r》,(,(r》A+*J"3ny,(。&0)&'其中，DJ。是位于r位置处荧光波长的漫射系数，所述=1/(3/^0))，&(。为位于r位置处的组织在荧光波长的散射系数，为己知值，/^(r)为位于H立置处的组织在荧光波长的吸收系数，为已知值，步骤(7.5).假想当把一个荧光波长的点光源放在所述检测点W处，t//eWl,"2,…,d"，使用Matlab软件包函数PCG计算所产生的荧光场在所有所述四面体网格尸em—il^A上的有限元节点上的强度分布向量CE^，向量(D』是iVxl列向量，其中"'u」否则步骤(7.6).按照下式生成所述各检测点^1，J2,…，^^处检测到的荧光强度同未知的荧光探针分布之间的线性关系，用矩阵R表示步骤(8).针对空间顺序分布的激发光点光源W,A^1,…,A^，以及相对应的所述检测点，h=^，ivxi列向量a,的列编号为y的元素a,l/]用下式表示:按照下式构建总的线性关系矩阵『，下述。已知值.是各激发光点光源对应的被检测荧光，为<formula>formulaseeoriginaldocumentpage12</formula>成Wx1列向量是未知的荧光探针在所述有限元网格Few—M^/z节点上的分布，Om,_是Mxl列向量，M所有激发光点光源对应的总的检测点数目，『是Mx7V矩阵；步骤(9).按下述代数表达式，根据步骤(8)得到的总的线性关系矩阵『求解位置的荧光探针分布f/:<formula>formulaseeoriginaldocumentpage12</formula>分别是f/中第)个元素在第"—/f"+l和"—!'^次迭代时的值，其中<formula>formulaseeoriginaldocumentpage12</formula>t/[a'，是中第^个元素在第"次迭代时的值，是矩阵^第^行，第^列的值，『[/,，/1是矩阵W第^行，第y列的值，。,画,用7表示，^,]是的r第A个元素，/1=0.1，C7迭代初始值设为^/=0。所述的朋g/e—/wv=-48°，""g/e—/z/g/z=48°，所述的/z".g/z/1—=-1.5c附,/7"'g似—/妙=1.5cw。为了验证本重建方法的可靠性，实验选择了一套在实验用鼠图谱上得到的仿真数据。从重建结果可以看出，重建的荧光体具有良好的定位和定量效果。图l荧光分子成像的机理l代表的是荧光滤光片。图2荧光分子断层成像的机理。图3实验用鼠图谱的外表面轮廓。图4非接触荧光分子断层成像系统结构图l代表的是荧光滤光片。图5实验用鼠外表面轮廓划分得到的有限元网格的部分显示。图6实验用鼠图谱模型上的仿真数据重建结果。图7计算机程序的流程框图。具体实施方式荧光分子成像的机理可以用图1说明激发光光源发出的激发光激发组织内的荧光探针(标记的蛋白，肿瘤，荧光染料等)发射荧光，传播到表面的荧光被CCD照相机检测到。在图l中，位于位置r、,的激发光点光源发出的激发光在实验用鼠体内传播，用实线箭头表示。激发光遇到荧光探针后，会被荧光探针吸收，然后荧光探针被激发，发射出更高波长的荧光。荧光在实验用鼠体内传播，用虚线箭头表示。到达激发光源对侧实验用鼠表面的荧光可以用CCD照相机检测到。CCD照相机前有荧光滤光片，保证只有被激发的荧光才被CCD检测到。由于生物组织具有强散射特性，因而激发光和被激发的荧光在实验用鼠体内都不是直线传播。这个成像的机理以及光在生物组织中的传播可以用两个耦合的称为漫射方程的偏微分方程组来表示其中下标e代表激发光波长，下标m代表荧光波长。①e(。和OJr)分别代表激发光和被激发荧光的强度分布。r代表整个几何模型内的空间位置坐标。在本文中，(...)表示函数操作，而[…]表示向量、图像、或者矩阵操作。例如，(D。(0表示激发光光强场在r处的值，而o丄/]则代表激发光强场向量c^的第y个元素。光在组织中的传播受到两个因素的影响，称为吸收和散射。吸收特性会吸收光的一部分能量，而散射特性会改变光的传播方向，使得光在组织中的不能直线传播。组织的吸收和散射特性在不同波长会不相同。A(。和Z)Jr)分别是位于r处的激发光和荧光波长的漫射系数，对应为D。(r)=l/(3Mre(0)和A>)=l/(3^>))，其中&(r)和&(0分别是组织在激发光和荧光波长的是散射系数。圪6("和/^一)分别是组织在激发光和荧光波长的吸收系数。在荧光分子断层成像中，把垂直入射到实验用鼠表面的激光束建模成入射点皮下一个传输自由程^=1/(3//〕处的位置r、.的激发光点光源50-O。-O代表在位置。处的一个冲激函数。t/(r)是位于r处的荧光探针的分布。荧光分子断层成像是利用已知的在表面位置处检测到的荧光强度，来重建未知的荧光探针的浓度分布C/(r)。其原理用图2简要说明分布在实验用鼠体内的所有荧光探针对检测点"l和d2处的检测到的荧光强度都有贡献，但是贡献的权重会有不一样。在图中，用荧光探针1，2，3简要代表实验用鼠体内的所有荧光探针分布。以检测点^1为例，荧光探针1，2，3对"l贡献权重分别为f/lxwl，t/2xw2，t/3xw3。其中/71，f/2，f/3分别是荧光探针l，2，3的浓度，而wl，w2，w3是位于荧光探针l，2，3位置处的单位浓度的荧光探针对检测点dl的贡献权重，称为单位贡献权重。单位贡献权重M，w2，w3可以通过荧光分子断层重建方法求解得到。对不同的激发光源-荧光检测点对，可以计算得到相应的单位贡献权重。利用所有计算得到的单位贡献权重以及检测到的荧光强度，就可以重建出未知的荧光浓度f/1，f/2，t/3。这个单位贡献权重在本文中被称为检测到的荧光强度同未知的荧光探针分布之间的线性关系。实际的重建方法计算有限元网格节点上的单位荧光浓度的荧光探针对检测到的荧光强度的单位贡献权重，然后可以用优化方法重建得到有限元网格节点上的荧光探针浓度分布。在荧光分子断层成像中，多个角度的激发光点光源顺序激发，同时检测到相应的荧光投影数据，从多个角度投影数据就可以断层重建出未知的荧光探针浓度分布。在荧光分子断层成像系统中，用CT成像设备先得到实验用鼠的解剖结构信息，然后利用医学图像分割的方法分割CT图像，得到实验用鼠的器官图谱。这里提出的方法可以在得到的实验用鼠图谱的基础上，利用测量到的多个角度的荧光投影数据，重建出未知的荧光探针浓度分布。方法的主要步骤包括步骤(l)，利用分割完毕的实验用鼠图谱，通过图像处理得到外表面轮廓，并将外表面轮廓建模成几何体，这个方法依次按以下步骤实现步骤(l.l).在已经得到的实验用鼠图谱中，图谱中组织的灰度值都是大于等于1，而背景的灰度值是0。因而，选定一个阈值^"=0.5，就能够分开图谱中的组织和背景。对图谱的每一层图像/，利用阈值Ar进行二值化，大于阈值的取l，小于阈值的取0，得到二值图5/:其中ww代表图像的行，而"/代表图像的列。步骤(1.2).用二值图轮廓线提取的方法得到二值图S/的外轮廓线S。方法如下对二值图3/的每一个像素点A/[tow,co/]:(1)如果5/[raw,co/]为1，而且其相邻的8个像素点至少有一个为0，那么这个像素点在外轮廓线S上；(2)如果6/[raw，co/]为0，而且其相邻的8个像素点至少有一个为1，那么这个像素点在外轮廓线S上。这样就得到了图谱中一层图像的外轮廓线。对图谱中所有的层都执行二值化，然后进行轮廓线提取，就能得到所有层的外轮廓线，用{5}表示。步骤(1.3).依次把外轮廓线5*从像素坐标对应到笛卡尔坐标，给这些轮廓线选定一个中心[tow。，co/。]，每个像素点的几何尺寸是Acm.那么以[row。,co/。]为原点，轮廓线上的任意一个点[tow,co/]都对应到几何坐标[xj]:_y=(co/—co/0)*A然后把笛卡尔坐标[x，y]对应到极坐标[a,i]，其中《是极坐标系中的角度，A是极坐标系中的半径。对每一条外轮廓线S上的所有点，每个点的半径R用角度"按下式做12次多项式拟合凡二Za《，其中c为拟合次数，c=0,l，".12。a是拟合得到的多项式系数，用{&}表示。"=l,2,...Smw7，其中S—m/m是轮廓线S上的所有点的数目。在均匀选定a—mw^30个角度&,…,aa—用多项式{/^}拟合在这些角度上的半径{^,'..,&■}。这样，每一层图像/所对应的轮廓线S表示为U",，…,"。—…A—步骤(1.4).每一层图像都能够得到一条轮廓线S，用{5}表示所有层图像所得到的轮廓线。在{5}中，按照高度均匀选定M层，用这M层轮廓线来代表整个外轮廓。在相邻两层之间，可以通过插值的方法得到轮廓坐标。也就是说任意指定一个角度a和高度i/，都能够相应计算得到其半径i。设高度//处于高度为//,和//2的两层(//1<//2)之间，那么其中/,,&分别是在Hp/7,层角度a所对应的半径，/^和i^分别是在巧，A所对应的12次多项式。步骤(1.5).得到的外轮廓从三维极坐标[a,i，/Z]转换到三维笛卡尔坐标[x，y,z]，这样就得到了在笛卡尔坐标系下的几何模型G，如图3。步骤(2).在图4所示的成像系统中，实验用鼠被放置在一个旋转平台上。激光从一侧照射实验用鼠，激发实验用鼠体内的荧光，被激发出的荧光在对侧被CCD照相机检测到。通过旋转平台，然后采集图像，就能在不同角度实现激发，也就对应多个激发光点光源。把步骤(1)得到的几何模型G导入有限元软件包中，同时把各个激发光点光源所对应的点加到这个几何模型G的节点上，然后用有限元软件做网格划分，离散成四面体网格i^附—MeA。图5部分显示了该四面体网格。该四面体网格的所有节点自动编号为/=1,一,^，iV是总的网格节点数。步骤(3).确定每个激发光点光源s以及其所有检测点所对应的有限元节点。在离散的四面体网格&w一ikfei中，每个激发光点光源都对应一个有限元节点，通过以下方式确定激发光点光源对应的节点对每个有限元节点/，如果^、.=。那么/就是激发光点光源s所对应的有限元节点，其中^是激发光点光源s所对应的位置坐标，而r是有限元节点/所对应的位置坐标。同时，对每个激发光点光源5，在其对侧对应的角度""g/e—/owtwg/e—Wg/以及高度/^'g似一/ow~/ze/g&—Wg/7范围内检测到荧光。在方法中，角度范围选为-48'~48°，而高度范围选为相对激发光点光源高度-1.51.5cm内。在方法中，要求检测点放置在有限元网格的表面节点，因而，可以用以下方式确定检测点对应的所有有限元节点对每一个表面节点/，如果ce{""g/e—/cw""g/e—/n.g/z}而且r,.e{/ze/g/^—/ow~—/z/g/}，另卩么节点/属于这个光源对应的检测点，否则不属于。对特定激发光点光源，找到相应检测点后，利用拍摄到的荧光图象，得到在这些检测点处的测量到的荧光强度。步骤(4).描述荧光分子断层成像中光传播的漫射方程组在有限元网格Fem—MeA上可以离散成一对线性方程组<formula>formulaseeoriginaldocumentpage16</formula>其中Wx1列向量cDt,和iVx1列向量。分别是是激发光和被激发荧光在网格节点上的强度分布；Wxl列向量"代表了荧光探针在网格节点上的分布；A^l列向量5、,代表了激发光点光源5(r-r)在网格节点的分布。在有限元方法中，己知了在离散的网格节点上的分布，可以通过有限元节点所对应的Lagmnge线性形函数来得到在任意位置r上的分布。以激发光强度为例，在任意位置r上的激发光强度(^(。，可以表示为0)=^0力']^(0，其中①』刀是列向量A的第7'个元素。①jy]也是激发光强度在编号为y的网格节点上的值。^,(r)是编号为的网格节点所对应的Lagrange线性形函数。&和《m在有限元方法中称为刚度矩阵。矩阵^是组织在激发光波长的吸收和散射特性在网格节点上的体现。矩阵《是组织在荧光波长的吸收和散射特性在网格节点上的体现。矩阵《是激发光强度分布OJ。对荧光探针的激发作用在网格节点上的体现。这三个矩阵都是WxiV的稀疏矩阵。对于矩阵《，/^和《，它们的行和列分别用/和y表示，所述的行/对应Fem—T^fey/z上的所有节点7V按照编号顺序排列而成的行，所述的列_/对应Few_ikfe^上的所有节点W按照编号顺序排列而成的列。针对位于位置^;的激发光源s,其表面检测节点WU2,…,^^上检测到的荧光强度同未知的荧光探针分布之间的线性关系可以按照以下步骤得到步骤(4.1).计算刚度矩阵&。^的第/行，第y列元素&[/j]可以用下式求得，其中Q代表了整个几何模型G内部所对应的空间区域；3Q则代表几何模型G的表面边界；下标e代表激发光波长；r是位置坐标，A(r)是位于r处的漫射系数，对应为化0)=1/(3&(r》；《和组织和空气的折射系数有关，一般取g=0.1511;y,0)和^0)是行编号为/，列编号为j'的有限元节点对应的Lagrange线性形函数，由所述有限元软件包给出。Vw,(r)=[3y,(r)/ax,3y,0)/^y,3y,(/)/&]和V(^(r)=[3<//,0)/&，3(^,0)//&]分别是所述Lagrange线性形函数^,(/)、^.(/)的梯度。稀疏矩阵《通过在有限元软件包中设置好组织在激发光波长所对应的光学参数z^,00和mUO的接口函数，利用有限元软件包求解得到。光学参数/^(r)和a。O0是这样获取的首先，获取位置r的图谱的灰度值，确定位置/"处所属的组织器官；然后，查表得到该组织器官在激发光波长对应的光学参数。下表就是不同组织器官在激发光和荧光波长的光学参数值<table>tableseeoriginaldocumentpage18</column></row><table>其中，光学参数的单位是cnT1。歩骤(4.2).计算激发光点光源s所对应的激发光在有限元网格Few—ikfe^节点上的强度分布Ot,，《A=A。^的第_/个元素0)丄/]代表了激发光光强在编号为./的有限元节点上的值。列向量A的第./个元素5、[/l按下式表示，p,有限元节点y位于所述/;位置上力]^0，否则°利用雅克比预处理的共轭梯度法来迭代求解这个线性方程组，使用Matlab软件包函数PCG来实现，输入矩阵^和向量i,，输出向量d^。步骤(4.3).计算矩阵F、。《的第z'行，第J'列元素《[/，)]表示为:其中ojo是激发光点光源^的光强场在空间位置r处的值，表示为o=^>。[刀1//，这个累加计算可以由有限元软件包来完成。F、通过在有限元软件包中设置好接口函数，利用有限元软件包求解得到。这个接口函数是这样的，输入位置r，输出r处的激发光强强度，。步骤(4.4).计算刚度矩阵《。矩阵尺第/行，第y列元素^[/，刀表示为，i["y]=J"q(ao)vw,(。.vw乂o)+(r)w,o》&+*j^y,o，乂。其中，下标m代表被激发光波长。这个矩阵的求解类似这个计算矩阵《，通过在有限元软件包中设置好组织在被激发光波长的光学参数^m(r)和/^(r)接口函数，利用有限元软件求解得到，光学参数Am(0和&(0是这样获取的首先，判断位置r处所属的组织器官；然后，査表得到该器官在荧光波长对应的光学参数。步骤(4.5).假设把一个荧光波长的点光源放在检测节点^^/€{^1^2,...,^}处，计算产生的荧光场在有限元节点上的分布向量a^，其中iVxl列向量Sw的第J个元素5J/^'":n|。这个线性方程组通过雅克比预处理的共wJL0，否贝U轭梯度法来迭代求解，见步骤(4.2).步骤(4.6).生成检测点{^/1^2,...,"}处检测到的荧光同未知的荧光探针分布之间的线性关系矩阵『、.，可以看到，对于每个激发光源-荧光检测对，例如卜J1，在检测点dl检测到的荧光同未知的荧光探针分布之间的线性关系通过一个向量-稀疏矩阵乘积。a尸,得到。步骤(5).针对多个顺序激发的光源AJ=l...iV,以及其相应的检测点，叠加生成总的线性关系矩阵『。按照歩骤(4)，分别建立每个激发光源A和其相应检测点的线性关系矩阵『、.p然后把这些矩阵叠加可以得到总的线性关系矩阵，如下所示①■si其中，iVxl列向量t/是未知的荧光探针在有限元网格Few—7kfeW节点上的分布。Mxl列向量0),,。,是所有检测点所对应的荧光强度值，其中M所有激发光点光源对应的总的检测点数目。『是MxiV矩阵。步骤(6).根据线性关系矩阵重建未知的荧光探针分布，利用代数重建法(ART)来实现重建，把0^,，用向量r表示，求解线性方程组『17=「时，按照下式迭代<formula>formulaseeoriginaldocumentpage20</formula>其中"w、t/[y]",是;y中第乂个元素在第"—/^+i和"」^次迭代时的值。是t/中第^个元素在第"—/^次迭代时的值。吗^,g是矩阵r第f,行，第^列的值。是矩阵ff第^行，第乂列的值。r[U是F第6个元素。义取值为义=0.1.迭代的初值取[/=0。示例的重建结果如图6所示。图6是重建的荧光探针分布在激发光点光源那个高度层上，坐标x,y各点画出的等浓度线分布。权利要求1.一种基于线性关系的荧光分子断层重建方法，其特征在于，所述方法是用一台计算机依次按以下步骤实现的步骤(1).把用CT成像设备分割完毕的生物实验用鼠的分层的图谱输入所述计算机中，该计算机通过下述图像处理方法得到外表面轮廓，并把所表述表面轮廓建模成几何体步骤(1.1).对所述图谱的每一层图像I，用一个相同的阈值thr＝0.5进行二值化大于所选阈值thr者取1，小于或等于该阈值thr者取0，把该图谱中的组织和背景依次用所述的“1”和“0”分开，得到每一层图像的二值图BI其中，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-row0)*Δ，y＝(col-col0)*Δ，其中，点[row0，col0]为该轮廓线的中心，Δ是每个像素点的几何尺寸，用cm表示，步骤(1.4).把步骤(1.3)得到的该笛卡尔坐标[x，y]对应到极坐标[α，R]上，α是角度，R是极坐标的半径，对所述外轮廓线S上的所有点，把每个点上的所述半径R用角度α按下式做多项式拟合c为拟合次数，c＝0，1，…12，pc是拟合得到的多项式系数，用{pc}表示，n＝1，2，…S_num，S_num是所述轮廓线S上的所有点的数目，在均匀选定α_num＝30个角度{α1，…，αα_num}，用多项式{pc}拟合这些角度上的半径{R1，…，Rα_num}，使得所述每一层图像的所述每一层轮廓线S表示为{{α1，…，αα_num}，{R1，…，Rα_num}，{p0，…，p12}}，用{S}表示所有层图像上I上得到的轮廓线，步骤(1.5)在所述所有层图像I的轮廓线{S}中，按照高度均匀选定M层，以此表示所述所有层图像I的外轮廓线，则相邻两层H1<H2之间，用下式表示任意一层轮廓线的坐标，用半径R表示为H1<H2，R＝R1*(H2-H)/(H2-H1)+R2(H-H1)/(H2-H1)，其中，R1，R2分别为在所述高度H1，H2层角度α所对应的半径，PH1，RH2分别是在高度H1，H2层所对应的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).把垂直入射到所述活体生物实验用鼠的激光束建模成入射光皮下一个自由传输自由程的位置rs处的激发光点光源δ(r-rs)，该δ(r-rs)代表在位置rs处的一个冲激函数，是rs处激发光波长的散射系数，步骤(4.2).对于每个所述四面体网格Fem_Mesh上有限元网格节点i，若rs＝ri，则该有限元网格节点i就是激发光点光源s所对应的有限元网格节点；步骤(5).确定所述每个激发光点光源s在所述Fem_Mesh上有限元网格表面上的检测点所对应的所述四面体网格Fem_Mesh上的有限元网格节点i若所述表面节点i满足以下关系ri∈{angle_low～angle_high}且ri∈{height_low～height_high}，所述angle_low，angle_high，height_low，height_high为设定值，则所述节点i属于所述激发光点光源s对应的网格节点；步骤(6).在得到步骤(5)所述各检测点后，利用所述CCD相机检测得到的荧光图像，得到各检测点发出的荧光图像；步骤(7).针对所述给定位置rs处的激发光点光源s，按照以下步骤计算对应所述激发光点光源s所述表面检测点{d1，d2，...，dL}上检测到的荧光强度同未知的荧光强度之间的线性关系，其中L是检测点总数步骤(7.1).按下式计算N×N维刚度矩阵Ke中的元素Ke[i，j]，该刚度矩阵Ke是生物组织在激发光波长在所述四面体网格Fem_Mesh上的有限元节点上的吸收和散射系数的体现，该刚度矩阵Ke的行和列分别用i和j表示，所述的行i对应所述Fem_Mesh上的所有节点N按照编号顺序排列而成的行，所述的列j对应所述Fem_Mesh上的所有节点N按照编号顺序排列而成的列，下同，其中任意一元素Ke[i，j]表示为其中，Ω代表所述的几何模型G内部对应的空间区域，代表所述的几何模型G的表面边界，De(r)是位于r位置处激发光波长的漫射系数，所述为位于r位置处的组织在所述激发光波长的散射系数，为已知值，μae(r)为位于r位置处的组织在所述激发光波长的吸收系数，为已知值，ψi(r)和ψj(r)是行编号为i，列编号为j的有限元节点对应的Lagrange线性形函数，由所述有限元软件包给出，和分别是所述Lagrange线性形函数ψi(r)、ψj(r)的梯度，q＝0.1511，步骤(7.2).使用Matlab软件包函数PCG按下式计算所述激发光点光源s在所述有限元网格Fem_Mesh节点上的强度分布向量Φe，Φe是N×1列向量KeΦe＝Bs，N×1列向量Bs的列编号为j的元素Bs[j]用下式表示已知了在节点上的强度分布向量Φe，得到任意位置r处的强度分步步骤(7.3).按下式计算N×N矩阵Fs的元素Fs[i，j]，该矩阵Fs是所述激发光强度Φe(r)对荧光的激发作用在所述四面体网格Fem_Mesh上的有限元节点上的体现，其中Fs[i，j]＝∫ΩΦe(r)ψi(r)ψj(r)dr，步骤(7.4).按下式计算N×N维刚度矩阵Km中的元素Km[i，j]，该刚度矩阵Km表示生物组织在激发光波长在所述四面体网格Fem_Mesh上的有限元节点上的吸收和散射系数的体现，其中其中，Dm(r)是位于r位置处荧光波长的漫射系数，所述为位于r位置处的组织在荧光波长的散射系数，为已知值，μam(r)为位于r位置处的组织在荧光波长的吸收系数，为已知值，步骤(7.5).假想当把一个荧光波长的点光源放在所述检测点dl处，dl∈{d1，d2，…，dL}，使用Matlab软件包函数PCG计算所产生的荧光场在所有所述四面体网格Fem_Mesh上的有限元节点上的强度分布向量Φdl，向量Φdl是N×1列向量，其中KmΦdl＝Bdl，N×1列向量Bdl的列编号为j的元素Bdl[j]用下式表示步骤(7.6).按照下式生成所述各检测点{d1，d2，…，dL}处检测到的荧光强度同未知的荧光探针分布之间的线性关系，用矩阵Ws表示步骤(8).针对空间顺序分布的激发光点光源sk，k＝1，…，Ns，以及相对应的所述检测点，按照下式构建总的线性关系矩阵W，下述Φm，meas是各激发光点光源对应的被检测荧光，为已知值N×1列向量U是未知的荧光探针在所述有限元网格Fem_Mesh节点上的分布，Φm，meas是M×1列向量，M所有激发光点光源对应的总的检测点数目，W是M×N矩阵；步骤(9).按下述代数表达式，根据步骤(8)得到的总的线性关系矩阵W求解位置的荧光探针分布U其中U[j]tt_iter+1、U[j]tt_iter分别是U中第j个元素在第u_iter+1和u_iter次迭代时的值，U[t2]tt_iter是U中第t2个元素在第u_iter次迭代时的值，W[t1，t2]是矩阵W第t1行，第t2列的值，W[t1，j]是矩阵W第t1行，第j列的值，Φm，meas用V表示，V[t1]是的V第t1个元素，λ＝0.1，U迭代初始值设为U＝0。2."A'…'—…'W、"細〉，(/V…，Pi2"'用{S}表示所有层图像上/上得到的轮廓线，步骤(1.5).在所述所有层图像/的轮廓线{^中，按照高度均匀选定M层，以此表示所述所有层图像/的外轮廓线，则相邻两层//,<//2之间，用下式表示任意一层轮廓线的坐标，用半径i表示为W=《*(//2-77)/(//2—//,)+/2(//—/(//2—H》，其中，《,^2分别为在所述高度//1,//2层角度"所对应的半径，/^,i^分别是在高度//,，//2层所对应的12次多项式，步骤(1.6).把步骤(1.4)(1.5)得到的三维外轮廓线从极坐标[a,i,/f]转换到三维笛卡尔坐标系[;c,少,zl，得到了所述三维笛卡尔坐标系下的几何模型G;步骤(2).把要进行活体荧光分子断层成像的一个生物实验用鼠放置在一个旋转平台上，激光器发出的激发光从一侧照射所述活体生物实验用鼠，激发光激发出荧光，被激发出的荧光在相对一侧被CCD相机检测，再把检测结果输入计算机中，分角度通过使所述旋转平台旋转，采集不同旋转角度被激发的荧光图像，相当于实现了在不同角度时的激发光点光源激发；步骤(3).把步骤(1)得到的几何模型G导入一个有限元软件包中，同时把所述各个激发光光源所对应的点加入所述几何模型的节点上，再用有限元软件包做网格划分，离散成一个四面体网格Few—ikfe^，该四面体网格的所有节点自动编号为z、l,…,iV，iV是总的网格节点数；步骤(4).按以下步骤确定每个激发光点光源s所对应的在所述四面体网格i^^^M^/z上的有限元网格节点步骤(4.1).把垂直入射到所述活体生物实验用鼠的激光束建模成入射光皮下一个自由传输自由程=1/(3&(O)的位置r、.处的激发光点光源50-O，该-O代表在位置。处的一个冲激函数，/^(。是r、处激发光波长的散射系数，步骤(4.2).对于每个所述四面体网格&m—Ma/7上有限元网格节点/，若^=;，则该有限元网格节点/就是激发光点光源^所对应的有限元网格节点；步骤(5).确定所述每个激发光点光源s在所述Fem一Me^上有限元网格表面上的检测点所对应的所述四面体网格Few—M^/z上的有限元网格节点/:3.则所述节点/属于所述激发光点光源s对应的网格节点；歩骤(6).在得到歩骤(5)所述各检测点后，利用所述CCD相机检测得到的荧光图像，得到各检测点发出的荧光图像；步骤(7).针对所述给定位置r、处的激发光点光源"按照以下步骤计算对应所述激发光点光源^所述表面检测点{"1^2,...,^:}上检测到的荧光强度同未知的荧光强度之间的线性关系，其中Z是检测点总数步骤(7.1).按下式计算A^7V维刚度矩阵&中的元素A[/,/1，该刚度矩阵《是生物组织在激发光波长在所述四面体网格Few—MeA上的有限元节点上的吸收和散射系数的体现，该刚度矩阵《的行和列分别用/和_/表示，所述的行/对应所述&m—M^/z上的所有节点iV按照编号顺序排列而成的行，所述的列/对应所述Fem—Afe^上的所有节点iV按照编号顺序排列而成的列，下同，其中任意一元素《[/,y]表示为d["刀=J"q(AK(,)'▽0)+(0)^+由J"a^,(,V,，其中，Q代表所述的几何模型G内部对应的空间区域，3Q代表所述的几何模型G的表面边界，A(。是位于r位置处激发光波长的漫射系数，所述AO)"/(3&(")，A。(。为位于H立置处的组织在所述激发光波长的散射系数，为已知值，/^(。为位于r位置处的组织在所述激发光波长的吸收系数，为已知值，y,(o和^.oo是行编号为/，列编号为y的有限元节点对应的Lagmnge线性形函数，由所述有限元软件包给出，V^,.(r)=[3^,(。/&，3^/,(。/^,3^(。/&]禾口V=[5y/,(r)/分别是所述Lagrange线性形函数y/,0)、W,O)的梯度，^=0.1511,步骤(7.2).使用Matlab软件包函数PCG按下式计算所述激发光点光源s在所述有限元网格Few—^fe^节点上的强度分布向量0。，d^是vVxl列向量已知了在节点上的强度分布向量^,,得到任意位置r处的强度分步Oe(。-Z气[j']〖A=5、，Wxl列向量A的列编号为y的元素S、,[/l用下式表示nf.，fi,有限元节点y位于所述r,位置上，4.步骤(7.3).按下式计算WxN矩阵F、,的元素F、.[/,/1，该矩阵《是所述激发光强度^(r)对荧光的激发作用在所述四面体网格Few_Mm/2上的有限元节点上的体现，其中步骤(7.4).按下式计算A^W维刚度矩阵i^中的元素^J/,/1，该刚度矩阵/^表示生物组织在激发光波长在所述四面体网格/^w—M^/7上的有限元节点上的吸收和散射系数的体现，其中&['■，./]=J"。("m(。▽k0)▽+(Ok0》^+iJ"3。k0V,0)&'其中，DJO是位于r位置处荧光波长的漫射系数，所述^=1/(3^)，z^O)为位于r位置处的组织在荧光波长的散射系数，为已知值，为位于r位置处的组织在荧光波长的吸收系数，为已知值，步骤(7.5).假想当把一个荧光波长的点光源放在所述检测点W处，^eW1^2,…^Z)，使用Matlab软件包函数PCG计算所产生的荧光场在所有所述四面体网格i^w—7kfe^上的有限元节点上的强度分布向量0^，向量0^,是iVxl列向量，其中kA=A,，A^xl列向量^,的列编号为y'的元素A,[/l用下式表示，=1，7.=t//Lo，否则步骤(7.6).按照下式生成所述各检测点Wl,d2,…,^^处检测到的荧光强度同未知的荧光探针分布之间的线性关系，用矩阵R表示—①①F、,；步骤(8).针对空间顺序分布的激发光点光源A,^:=1,...，乂，以及相对应的所述检测点，按照下式构建总的线性关系矩阵『，下述o^,、是各激发光点光源对应的被检测荧光，为已知值「l①.w、.iVxl列向量t/是未知的荧光探针在所述有限元网格Few7kfey/2节点上的分布，0>5.是Mxl列向量，M所有激发光点光源对应的总的检测点数目，PF是MxiV矩阵；步骤(9).按下述代数表达式，根据步骤(8)得到的总的线性关系矩阵『求解位置的荧光探针分布"M厂w-iy^"2]"K]"-齡t/[乂]"--十1二"[乂]"-'旨+义￡^-『[qJ]射f/[刀"」"、"[刀"」""分别是C/中第/个元素在第"—+1和"—次迭代时的值，t/[g",是[/中第^个元素在第"—ter次迭代时的值，吗^g是矩阵『第U亍，第^列的值，W[f,J]是矩阵W第f,行，第/列的值，气，，用F表示，化]是的K第/,个元素，义=0.1，"迭代初始值设为^/=0。2.根据权利要求l所述的一种基于线性关系的荧光分子断层重建方法，其特征在于，所述的""g/e—/ow=-48°，""g/e—/z/g//=48°，所述的/ze/g//"—/ow=-1.5cm,/ze/g/^—/z/g/=1.5c附。全文摘要一种基于线性关系的荧光分子断层成像重建方法属于在生物医学中图像重建的
技术领域：
。其特征在于，依次含有以下步骤一、实验用鼠图谱数据通过图像处理和数值拟合得到外表面轮廓；二、这个外表面轮廓被导入有限元软件包进行网格划分；三、找出每个激发光源以及相应的荧光检测点对应的有限元网格节点；四、针对一个特定的激发光源sk，建立表面检测到的荧光强度同未知荧光探针分布之间的线性关系矩阵W_sk；五、针对所有激发光源，生成总的线性关系矩阵；六、得到这个线性关系后，代数重建方法被用来迭代逼近真实的荧光探针分布。本发明具有可处理非均匀光学参数以及快速建立表面检测到的荧光强度同未知荧光探针分布之间的精确线性关系的特点。文档编号G01N21/64GK101396262SQ20081022544公开日2009年4月1日申请日期2008年10月31日优先权日2008年10月31日发明者欣刘,李明泽,汪待发,王洪凯,净白,燕金,陈延平,陈欣潇申请人:清华大学