1 引言

荧光分子断层成像(FMT)是近几年发展起来的一种医学成像方法,因其具有非电离辐射、低成本、高灵敏度等优点,已被成功应用于癌症检测、药物研发、基因表达可视化等领域,成为医学成像技术领域的研究热点^[1-5]。荧光分子断层成像的原理是将荧光生化标记物(如ICG)作为对比剂注射到生物体中,这些标记物将向病变组织聚集。用特定波长的光激发荧光染料,荧光分子吸收能量后跃迁到激发态,经过特定的时间衰减回基态,并发出比激发光波长更长的发射光^[6-7]。通过测量组织边界处的荧光强度,并考虑光在生物组织中的传播模型,重建出组织内部具有荧光光学特性(聚集度、寿命、量子产额、组织对荧光的吸收散射系数等)的图像。

近红外光在生物组织中的传播具有较强的散射性,光子在生物组织中被多次随机散射,这将导致荧光分子断层成像的重建结果与荧光团的深度具有高度相关性。实验表明,随着成像深度的增加,成像的灵敏度将会呈指数下降^[8]。针对荧光成像中这种深度灵敏度的非线性问题,人们提出了一些深度补偿的策略,如深度依赖的色散补偿算法^[9]、基于层的S型调整法^[10]、权值深度补偿法^[11]等,用于补偿因深度不同而造成的成像误差。事实上,通过实验发现,成像灵敏度除了与荧光团深度具有非线性关系,与荧光团的体积大小也有很大关系。对于生物组织体内具有的相同光学参数、相同深度,但不同体积的荧光团,重建出的荧光团光学参数存在较大误差。这就需要有一种基于体积补偿的方法来平衡荧光成像对体积灵敏度的影响。

本文提出了一种基于体积补偿策略的荧光分子断层成像算法,应用于透射式连续波荧光成像模式中。该方法利用迭代自组织数据分析技术(ISODATA),对预迭代的重建图像进行聚类分析,得到各荧光团的体积。再根据荧光团的体积大小设置权值系数,基于对数运算方法,对重建的荧光光学参数进行补偿,使重建结果更加精确。仿真实验结果表明,该补偿算法可以有效修正因荧光团体积大小不一而造成的成像重建误差,从而显著提高重建质量。

2 荧光分子断层成像重建

2.1 正向传播模型

对于激发光为连续波模式的荧光分子断层成像,近红外光的传播模型可以近似表示为两个耦合的扩散方程^[12-13]:

∇·[De(r)∇Φe(r)]-μaerΦe(r)=-Iδ(r-rl)∇·[Dm(r)∇Φm(r)]-μamrΦm(r)=-Φe(r)β(r),(1)

式中r∈Ω,其中Ω为成像区域;Φ_e,m为激发光(下标为e)或发射光(下标为m)的光子密度;μ_ae,am为光吸收系数;D_e,m=1/3(μ_ae,am+μ'_se,sm)为扩散系数(其中μ'_se,sm表示简化的扩散系数);δ(r-r_l)表示在不同L个位置r_l(l=1,2,…,L)处的各向同性激发光源的位置;I为激发点的光的强度;β=ημ_af为待重建的未知荧光团,其中η为量子产额,μ_af为荧光团相对于激发光的吸收系数。

为了求解(1)式,引入(2)式所示的Robin边界条件:

2De,m∂Φe,m/∂n+qΦe,m=0,(2)

式中q=(1+γ)/(1-γ),γ为边界上反映折射率失配的参数,n为边界上外法向量。

扩散方程的求解方法有很多种,主要有解析法、随机统计法和数值方法等。由于处理不规则边界具有独特灵活性,数值方法中的有限元法常被用来解决扩散方程问题,其基本思路就是将满足(1)式的古典解转化为求解相应变分问题的弱解或广义解,并且在有限维空间投影。基于有限元法,(1)式最终可以表示为^[14-15]

Φ=Ax,(3)

式中Ф为相对于所有激发源,在边界上测量到的M维光子密度;A为M×N维灵敏度矩阵;x为N维待重建的未知荧光参数。

2.2 逆向问题

荧光分子断层成像的正向问题是根据给定的光源和假设的生物组织光学参数分布,结合组织体内的结构信息建立的光学传播模型,预测光学测量值。而其逆向问题则是根据已知的光源分布和实际边界测量值,求解生物组织内部光学参数分布。逆向问题可以表示为

x=A-1Φ。(4)

基于模型的迭代逆向图像重建主要分为3步:1)基于一个假定的光学参数空间分布的正向模型来预测探测器的测量值;2)设计一个比较预测值与实际测量值的目标函数;3)调整光学参数的空间分布,直到目标函数达到最小值^[13,16]。重建算法如图1所示。

图1中J(x)为目标函数,满足

minx≥0J(x)=0.5‖Ax-Φ‖22,(5)

故荧光分子断层成像的逆向重建就是要求解一个最小二乘函数问题。

相对于生物组织内的未知光学分布x,在组织表面实际测量的数据Φ_meas的数量是非常有限的,因此属于欠定问题,方程有无限多解。同时,由于光在组织中传输的衰减特性,投影矩阵A一般为一个病态的矩阵,因而,逆向问题是典型的不适定问题。为了降低该问题的病态性,通常采用正则化方法^[17-18]。其基本思路是对解空间加上先验限制条件,来提高目标函数问题的稳定性。Tikhonov正则化方法是在最小化问题中添加了L₂范数的约束条件,限定待重建的x中的元素波动范围,即

minx≥0J(x)=0.5‖Ax-Φ‖22+λ‖x‖22,(6)

式中λ为正则化参数。该公式仍然是标准的二次函数,处理该函数的最优化方法很多,如牛顿法、共轭梯度法、扰动法等。

图 1. FMT重建算法流程

Fig. 1. Flow chart of FMT reconstruction algorithm

下载图片 查看所有图片

2.3 基于改进ISODATA算法的重建数据聚类

基于有限元方法对荧光分子断层成像的两个耦合方程进行求解,可以求得生物体内荧光参数。为了减少荧光团尺寸对重建结果的影响,首先对这些荧光参数进行分类,确定这些荧光参数分别属于哪个荧光团,然后计算各荧光团的体积,再根据体积进行参数补偿和调整,最后计算更加准确的重建图像。重建的荧光参数数据特征相同,且没有已知类别信息,因此需要一种非监督学习方法进行聚类分析。

2.3.1 ISODATA算法及其改进

ISODATA是一种非监督学习的动态聚类方法,可根据样本的特征进行自动聚类并对聚类数进行修正。与K均值算法相比,ISODATA算法以迭代的中间结果为依据,将全部样本按照相似度准则聚类,重新计算样本均值作为新的聚类中心,因此具有自组织性^[19-20]。同时,ISODATA算法在聚类分析过程中可进行动态地合并和分裂,对样本不断调整,因而能得到更加合理的聚类结果。然而,一方面,传统的ISODATA算法对初始聚类中心的位置选择是随机的^[21-23],不合理的初始聚类中心选择会使迭代的收敛速度慢,且易导致运算陷入局部最优解。另一方面,传统ISODATA算法对初始期望分类个数的定义是盲目的,传统的算法不考虑数据集的结构如何,是否具有可分性,已分的聚类数是否合理。

一个优良的初始聚类中心集应满足两个条件:1)初选的聚类中心属于不同的类,即任意两个聚类中心不能属于同一个类;2)选择的初始聚类中心能够尽量接近最终迭代收敛时的聚类中心。对于FMT的重建图像而言,荧光团在生物组织内部是相对稀疏的,重建的荧光团吸收系数或荧光产额比背景光学参数略大。聚类的目的是把具有相似荧光团参数的数据点归为一类,把间隔一定距离的荧光团分成不同类。对ISODATA算法的改进是希望能够根据FMT图像重建的这个特点,选取初始聚类中心,而不是简单地随机选取。

针对传统ISODATA算法的不足,提出的改进算法如下。

1)在图像重建得到的荧光图像数据x_i(i=1,2,…,N)中取最大值,以某个正数d₁为直径,在图像数据空间中画一个球形邻域,将落入该邻域的数据点作为初始聚类。

2)给定一个正数d₂,在离开第一个聚类中心d₂距离之外寻找一个次大的数据点,作为第二个聚类中心。这样可以避免聚类点的过分集中,从而提高聚类效率。

3)以此类推,依次选出C个聚类中心。

考虑到荧光重建数据的实际情况和医学影像的经验,为了提高聚类速度,按照上述方法确定的初始聚类中心,需要在初始期望聚类数量上给出上限。在依次确定距离d₂之外的次大数据点作为聚类中心时,需满足

ξj≥ρ×xb,j=1,2,…,C,(7)

否则停止设置初始聚类中心。式中ξ_j为待设置的初始聚类中心的值,x_b为FMT重建的背景光学参数值,ρ为一正值常数。(7)式的含义在于需要设置一个初始聚类中心的阈值,当待设的初始聚类中心小于该值时,则停止设置聚类中心。

总之,改进的ISODATA算法的步骤如下。

1) 读入N个待分类的重建图像数据样本x_i(i=1,2,…,N)。

2) 设置聚类的一系列控制参数。K为预期的聚类个数;N_{j_min}为一个聚类中允许的最小的样本数,若小于此数,则不能作为独立的聚类;d_S为一个聚类中样本距离分布的标准差,若大于此数,则需要分裂;d_min为两个聚类之间的最小距离,若小于此数,则两个聚类合并;P为每次迭代运算中可以合并的最多类数;Q为允许的迭代的次数。

3) 设置初次聚类半径d₁和聚类中心之间的距离d₂,经过初次筛选,得到C_init个样本ξ_j(j=1,2,…,C_init)作为候选初始聚类中心。

4) 判断初始聚类中心是否满足ξ_j≥ρ×x_b,确定初始聚类中心数量C。

5) 求各个样本到聚类中心的距离,按照欧氏距离最小的原则将各样本分类。若有‖x_i-ξ_j‖<‖x_i-ξ_i‖,i,j=1,2,…,C,j≠i,则x_i∈R_j,R_j为第j个聚类。

6) 修正各聚类中心的值ξ_j=(1/N_j) ∑i=1Njx_i,i=1,2,…,N_j,j=1,2,…,C,其中N_j为第j个聚类中的样本个数。

7) 计算各聚类域中的样本到聚类中心的平均距离 d-j=(1/N_j) ∑i=1Nj‖x_i-ξ_j‖。

8) 计算全部图像重建样本对其相应聚类中心的总平均距离 d-=(1/N_j) ∑j=1CN_jd-j。

9) 分裂操作。计算每个聚类R_j中的样本分布标准差σ_j=(1/Nj)∑xi∈Rj(xi-ξj)2,j=1,2,…,C。若σ_j>d_S且 d-j>d-,且满足C≤K/2,则将聚类R_j分成两个聚类,相应的聚类中心分别为: ξj+=ξ_j+ασ_j, ξj-=ξ_j-ασ_j。其中,α=[0,1]。α的选取应使 ξj+和 ξj-都在R_j的聚类空间里,而其他聚类R_i(i≠j)到 ξj+和 ξj-的距离较远。原来的聚类中心值ξ_j取消,并使聚类数加1。

10) 合并操作。计算聚类中心两两之间的距离d_ij=‖ξ_i-ξ_j‖,i,j=1,2,…C,j≠i。 将d_ij<d_min的值按照从小到大排序,将最小的两个聚类合并,使聚类数减1。新的聚类中心值为ξ_l=(N_iξ_i+N_jξ_j)/(N_i+N_j)。

11) 如果是最后一次迭代或过程收敛,则运行结束,否则转到第5)步,迭代次数加1。

2.3.2 聚类结果评价

FMT重建的荧光团聚类划分之后,需要对聚类结果进行评价。通过定义准则函数可以将聚类结果明确地表达出来,以两个有效的准则函数来评判荧光团聚类划分的质量^[24-25]。

1) 误差平方和准则函数。误差平方和准则函数定义为

JC=∑j=1C∑i=1Nj‖vij-v-j‖2,(8)

式中v_ij为第j类中的第i个样本值, C为聚类的个数, v-j为聚类j中的样本均值,即

v-j=(1/Nj)∑i=1Njvij,j=1,2,…,C。(9)

J_C衡量的是用C个样本均值 v-1, v-2,…, v-C分别代表C类样本时产生的平方和误差。所以,误差平方和聚类准则J_C越小越好。

2) 类间距离平方和准则函数。类间距离平方和准则函数定义为

Jb=∑j=1Cv-j-v-,(10)

式中 v-为全部N个聚类样本的均值,即

v-=∑i=1Nvi。(11)

J_b衡量的是用C个聚类的样本均值相对于整体样本均值的方差,表示不同聚类之间的差别。所以,J_b值越大越好。

2.4 基于体积补偿的FMT图像重建

基于2.3节ISODATA算法的聚类分割,将成像区域Ω上重建的荧光团按照体积大小进行划分。接下来的问题是如何依据体积的大小,对荧光参数进行补偿,以便得到更准确的重建结果。

经过仿真实验研究发现,对于生物组织体内具有相同光学参数、相同深度的荧光团,体积越小的荧光团,重建得到的光学参数比实际值越小;而体积较大的荧光团,重建得到的光学参数相对较准确。为了解决这一问题,设计了一种非线性单调递减函数,使得体积补偿系数准确可控,如

yi=b×logaui+1,i=1,2,…,C,(12)

式中y_i为补偿系数,0<a<1,b为正常数,0<u_i<1为归一化的第i个荧光团的体积权值系数,具体计算表达式为

Vi/ui=V1/1,i=1,2,…,C,(13)

式中V_i和V₁分别为第i个荧光团和最大荧光团的体积,V_i≤V₁。通过设置合适的参数a、b值,调节递减函数的陡度系数,可以较准确地得到体积补偿系数值。

基于体积补偿的FMT图像重建流程如图2所示。

图 2. 基于体积补偿的FMT图像重建流程

Fig. 2. Flow diagram of FMT reconstruction algorithm based on volume compensation

下载图片 查看所有图片

3 仿真实验结果及分析

为了验证基于体积的补偿算法对荧光图像重建的有效性,开展了一系列仿真实验。假设散射系数已知,以荧光团的吸收系数μ_af为重建对象。基于有限元的求解方法,正向数据可通过(1)式获得,逆向求解采用基于雅克比矩阵的扰动法。本文的重建算法实现是基于MATLAB R2014a,计算机的CPU为Intel i7-4790S 3.2 GHz,内存RAM为8 GB。

为了评价图像重建算法的性能,引入了荧光团吸收系数μ_af的真实值和重建值之间的均方误差(MSE,e_MS)和相对误差e_r作为评价指标^[26]。

eMS=(1/n)∑i=1n(xi-x'i)2,i=1,2,…,n,(14)er=(1/n)∑i=1nxi-x'ixi,i=1,2,…,n,(15)

式中,n为重建区域内的待求目标点个数,x_i为设定的真实荧光参数值,x'_i为重建得到的荧光参数值。

3.1 双目标体体积补偿

二维圆形数字仿体的直径设置为25 mm,两个深度相同、光学参数相同的目标体分布在中心的两侧,如图3(a)所示。基于有限元的求解方法,数字仿体区域Ω可以被剖分为三角形网格,共生成1885个三角元和993个节点,其中边界上的节点数为99个(探测器可以获取的数据在边界节点上)。如图3(b)所示,红色区域为两个目标体的位置及大小。两个目标体的横截面面积的比值为22∶11。

10个各向同性的光源和30个探测器等间距地分布在圆形仿体的边缘。单个光源每激发一次,30个探测器就可以获得30个发射光数据。所以,总共可以获得300组边界上的测量数据。仿真模型中的背景光学参量和两个目标体的荧光参量如表1 所示。由于激发光和发射光的波长相差较小,可忽略它们两者之间的光学参数的差别。

在改进的ISODATA聚类算法中,取d₁=0.5 mm,d₂=0.8 mm,ρ=4,FMT图像预迭代次数为10次。图4和图5为预迭代图像重建之后,基于ISODATA算法的聚类分割结果评价。x轴为迭代次数,y轴为相应的J_C和J_b值。由图4和

图 3. 双目标体的图像重建。(a)双目标数字仿体;(b)有限元三角剖分

Fig. 3. Image reconstruction of double target phantom. (a) Double target phantom; (b) triangular mesh for FEM

下载图片 查看所有图片

图 4. 算法改进前后类内误差平方和JC对比图

Fig. 4. Convergence performance of different clustering algorithm in terms of the cluster uniformity JC

下载图片 查看所有图片

图 5. 算法改进前后类间距离平方和Jb对比图

Fig. 5. Convergence performance of different clustering algorithm in terms of the cluster uniformity Jb

下载图片 查看所有图片

图 6. 双目标体的图像重建结果。(a)未经补偿的重建图像;(b)体积补偿后的重建图像

Fig. 6. Reconstructed results for double target phantom. Reconstructed results (a) before and (b) after volume compensation

下载图片 查看所有图片

图5可以看出,虽然改进前后的ISODATA算法都得到了较好的收敛结果,但是改进后的算法在J_C和J_b的初始值和收敛速度上有较大优势,经过5次迭代,改进的算法就可基本收敛到最优解。

图6(a)为基于Tikhonov正则化,利用扰动法进行重建的结果(未经补偿)。其中,迭代次数为60次,正则化参数λ=0.00001。左侧体积较小目标体的荧光参数明显小于右侧,重建误差明显。图6(b)为基于对数函数,利用非线性体积补偿重建算法得到的图像,(12)式中系数取a=1/2,b=1。由图6可以看出,本文的补偿策略显著提高了较小荧光团的重建参数准确度。

为了测试算法的稳健性,在正向数据中加入不同信噪比(SNR,R_SN)的高斯噪声,分析比较本文算法在不同噪声数据下的聚类能力。如图7所示,图7(a)、(b)、(c)分别为在模拟正向数据中加入信噪比为40 dB,30 dB,20 dB的高斯噪声,体积补偿前后的图像重建结果。从图7中可以看出,即使正向数据的信噪比为20 dB时,本文聚类算法也能准确地确定目标,并进行参数补偿。当正向数据的噪声进一步增大时,图像的重建质量越来越差,甚至算法不能收敛,无法完成图像重建。

图 7. 加入高斯噪声后体积补偿前后的重建图像对比图。(a)加入信噪比40 dB的噪声;(b)加入信噪比30 dB的噪声;(c)加入信噪比20 dB的噪声

Fig. 7. Reconstruction images before and after the volume compensation when adding Gaussian noises of (a) RSN=40 dB, (b) RSN =30 dB, and (c) RSN =20 dB

下载图片 查看所有图片

为了进一步定量分析基于体积补偿的重建算法性能,对重建图像的均方误差、相对误差和计算时间等指标作进一步分析比较,如表2所示。由表2可以看出,经过体积补偿的重建图像的均方误差和相对误差都得到了明显降低。由于体积补偿算法中增加了基于ISODATA的聚类划分过程和函数补偿运算,整体的重建迭代过程将占用更多的CPU运算时间。

3.2 三目标体体积补偿

为了进一步评价本文算法的性能,设计了3个荧光团的数字仿体实验。如图8(a)所示,在3.1节中原有的双目标体的基础上,再增加一个荧光目标体。3个目标体设置为相同荧光参数和相同的深度,其体积比为22∶11∶6。图8(b)为基于有限元法的三角剖分图,图8(c)、(d)分别为体积补偿前后的重建图像。补偿函数及ISODATA聚类算法中的参数选取与上述两个荧光团图像重建中的设置一致。从主观视觉上比较,本文的非线性体积补偿算法较好地克服了小目标体重建荧光参数误差较大的问题,使重建图像的质量显著提高。

表3列出了上述3个目标体图像重建的算法性能比较,从均方误差和相对误差的量化指标上看,体积补偿后的重建图像质量明显优于补偿前的图像。

4 结论

针对在荧光分子断层成像过程中,具有相同荧光参数和相同深度,但不同体积的荧光团在图像重建后存在较大误差的问题,提出了基于体积补偿的图像重建方法。该方法首先利用改进的ISODATA算法对初步的重建图像聚类划分,设计了初始聚类中心选取和聚类个数确定的方法,以提高ISODATA算法的聚类效率。根据聚类算法得到不同荧光团之间的体积比率,基于一种体积权值系数的对数运算,对重建荧光参数进行非线性补偿。仿真实验结果表明,本文提出的基于对数函数的非线性体积补偿算法,能够有效减少图像重建误差,显著提高重建质量。

图 8. 3个目标体的图像重建。(a) 3个目标体;(b)有限元三角剖分图;(c)未经补偿的重建图像; (d)体积补偿后的重建图像

Fig. 8. Reconstruction results of triple target phantom. (a) Triple targets phantom; (b) mesh for FEM; reconstructed results (c) before and (d) after volume compensation

下载图片 查看所有图片

本文的体积补偿算法是在以扰动法为逆向重建算法的基础上设计的,当目标体之间的距离小于2 mm或目标体深度大于6 mm时,重建的目标体图像融合到一起,难以聚类分开,本文补偿策略就不再适用。今后将重点研究扰动法以外的更高效的逆向重建算法,以提高图像的分辨率和成像深度,统筹考虑荧光目标体的深度和体积等因素,给出综合补偿策略。另外,本文算法中的参数是基于经验值选取的,如何选取最佳的算法参数以达到最准确的数据补偿,如何设定一个优良的补偿标准,这也将在今后的研究中加以解决。