多光谱生物发光断层成像的光源形状重建

余景景，相文彬

(陕西师范大学 物理与信息技术学院,陕西 西安 710119)

生物发光断层成像(bioluminescence tomography,BLT)是一种非侵入式的,无电离辐射的光学分子影像技术。结合基因工程,使靶细胞表达出荧光素酶,荧光素酶与荧光素底物相结合,释放出近红外光。通过采集生物体表面的近红外光强度分布,BLT可反演出生物体内靶细胞的能量密度的空间分布信息[1-3],从而为定性、定量地分析细胞的生理活动与新陈代谢服务。因此,BLT在癌症的早期筛查和新型药物的研发方面具有广泛的应用前景[4-5]。

近红外光在生物体内要经历复杂的吸收、散射等过程,辐射传输方程(radiation transfer equation, RTE)能够描述光在生物体内的传输,但由于RTE计算复杂,求解困难,所以，在实际应用中一般采用RTE的近似形式。其中,RTE的一阶近似模型——扩散近似(diffusion approximation, DA)是目前应用最为广泛的光传输模型[6-8]。但是，由于DA模型仅适用于散射为主的介质,因此，基于DA模型的光源重建难免存在较大的模型误差。为平衡模型的精度和效率,研究者们还探索了各种RTE的高阶近似模型,包括简化球谐近似(simplified spherical harmonics equations,SPN)、离散坐标(discrete ordinates,SN)近似和球谐近似(spherical harmonics,PN)等,其中SPN相对效率更高[9]。对于简化球谐近似模型,Yang等人研究表明SP3比DA有着更高的准确度和更为广泛的应用空间[10]。Lu等人在匀质模型上进行了基于SP3的BLT重建,初步展示了SP3模型的优越性[11]。金晨等人探索了各阶次SPN模型的准确性,基于SP3模型提出了变量分离近似稀疏重构的重建方法,提高了重建精度[12]。近年来,光学分子影像在其他领域的研究也表明，结合SP3模型的重建方法可以提高光源或靶目标中心的定位精度[13-14]。

近年来BLT重建的研究致力于进一步提高光源重建的精度或稳定性。这些方法的研究重点在于增加先验信息和优化重建算法，以克服BLT重建问题的高不适定性,例如Liu等人结合多光谱测量数据和水平集策略,提高重建光源对噪声的鲁棒性[15]。Yu等人提出了基于L1/2正则化的非凸稀疏重建算法，进一步提高重建图像质量[16]。但是，目前多数BLT重建方法都侧重于定量评估光源中心位置和能量偏差,对于光源形状的拟合关注较少[17]。实际上,光源形状也是关键信息,例如在BLT引导放疗应用中,重建的光源形状与靶细胞的大小和体积相对应,对于三维适形放疗有重要意义。Gao等人基于RTE前向模型和多级自适应有限元方法,结合TV范数和L1范数正则化进行重建,在二维仿真实验中对于不同形状的光源获得了较好的重建结果[18],但还没见三维成像结果的相关报告。

本文探索了BLT中的光源形状重建问题,为进一步提高光源重建质量及形状拟合程度,采用SP3传输模型减少模型误差,运用L1范数正则化和多光谱测量数据克服BLT逆问题的不适定性，其中，重点是研究光传输模型对光源形状拟合的影响,因此,实验中着重与基于扩散方程的重建方法进行了对比,除了对重建图像的视觉评价外,本文采用重建光源与真实光源交叠部分的体积与重建光源总体积之比,定量评价光源形状重建的优劣。

1 理论方法

RTE的三阶简化球谐近似模型可表示为[9]

(1)

其中，μai=μa+(1-g)iμs,μa为生物组织的吸收系数,μs为生物组织的散射系数,g为各向异性系数，S是光源项,代表内部光源的能量密度分布，φi为辐射度的勒让德矩的线性组合。SP3相应的Robin边界条件为[9]

(2)

其中，n为方向向外的法向量,A1,B1,C1,D1,A2,B2,C2,D2是常数,由文献[9]给出。

结合有限元(finite element method,FEM)方法,将SP3方程及其边界条件转化为线性矩阵方程的形式[19]，

(3)

(4)

其中，Miφi+是M矩阵求逆后Miφi对应的分块子矩阵。生物体的边界光流量和内部光源密度的线性关系可以表示为[19]

J+=β1φ1+β2φ2。

(5)

其中，β1,β2是常数[19]。由于在BLT实验中,只能采集到逸出生物体表面的近红外光。除去矩阵中不可测量的部分,可得到表面能量分布和光源能量密度之间的线性关系，

(6)

当采用多光谱测量值时,波长为λ的近红外光与内部光源之间的线性关系可以表示为

(7)

其中，A(λi)代表不同波长近红外光的系数矩阵,η(λi)是不同波长相关频谱的权重[20]。为了避免测量值由单个波长的光主导,采用文献[21]的方法进行归一化，对每个波长的测量值除以该波长的最大测量值，为保证等式两边成立,对等式右边的系数矩阵也作相应处理。可得

(8)

为了克服BLT问题固有的不适定性,需要通过正则化技术在重建中结合尽可能多的先验信息,例如早期常用的各种L2范数正则化方法,以及近些年广泛采用的稀疏正则化方法。本文结合BLT应用中光源分布的稀疏特性,采用L1正则化的方法求解式(8),将式(8)转化为

(9)

其中，τ是正则化参数。

对于式(9)的L1范数最小化问题,已经有众多的求解方法,如迭代收缩阈值(iterative shrinkage and thresholding,IST)算法[22]、变量分离近似稀疏重构(sparse reconstruction by separable approximation,SpaRSA)算法[23]、不完全变量截断共轭梯度(incomplete variables truncated conjugate gradient,IVTCG)算法[24]等。

本文的重点是探讨不同的光传输模型对形状重建的影响,同时，由于IVTCG算法的稳定性和高效性,更重要的是,它已经得到了光学分子影像领域的广泛应用,证实了它的稳定性和有效性[24-25],因此，我们选择用IVTCG算法求解式(9),在下文中，将两种重建方法简记为DA+IVTCG和SP3+IVTCG。

2 实验与结果

为了较为系统地评估两种方法，在非匀质数字鼠模型Digimouse[26]上设计了一系列的单、双光源仿真实验,对比了多种光源设置下两种方法的重建结果。实验中,光源均设置在躯干部分,因此仅截取了该数字鼠模型长为40mm的躯干部分,包括心、胃、肝、肾、肺5种器官还有肌肉组织。为了模拟多光谱测量数据,采用波长分别为670nm，690nm和710nm的单色光。用于计算前向数字鼠表面测量值的光学参数如表1所示[27]。下面的实验中,前向网格随光源设置得略有不同,节点数均在12 000～13 000,同时，为了公平比较,所有重建采用同一个后向网格,后向网格包含6 038个节点和32 880个四面体。

为了评价实验得出的结果,本文采用的量化指标分别为,实际光源中心点到重建光源中心的距离(location error,LE)、重建光源的平均能量密度(density),以及真实光源与重建光源的交叠体积和重建光源总体积之比(ratio between the overlapping and total volume,OVTVR)[28]。其中,OVTVR是用来定量评价重建光源形状优劣的指标。

2.1 单光源仿真实验

分别在低散射高吸收介质(肝)和高散射低吸收介质(肺)中设计了两组单光源重建实验。实验使用了3种尺寸的光源,分别是小光源(r=0.5mm,h=1mm)、中光源(r=1mm,h=1mm)和大光源(r=1mm,h=2mm)，其中r代表真实光源的底面半径,h代表真实光源的高。

第一组实验是重建低散射高吸收介质(肝)中的单光源,为了研究光源尺寸和深度对重建结果的影响,将单个不同尺寸的光源放置在肝脏的不同位置处,以获得不同的光源深度。具体地,光源中心位置分别为(14,8,16.6)mm,(14,10.5,16.6)mm和(14,12,16.6)mm,其中，距离数字鼠表面最深的光源位置为(14,12，16.6)mm,相对数字鼠背部表面深度为9mm,相对数字鼠腹部表面深度为9mm,最浅的光源位置为(14,8，16.6)mm,相对数字鼠背部表面深度为13mm,相对数字鼠腹部表面深度为5mm。真实光源的能量密度均为1nW/mm3。

图1展示了第一组实验中,DA+IVTCG和SP3+IVTCG对于中心位置为(14,8，16.6)mm处的不同尺寸光源的3D重建结果。重建光源定义为大于重建最大能量密度的10%的四面体,表2为重建结果的各定量指标的值。由图1以及表2可以看出,SP3+IVTCG不仅在光源中心定位和重建能量密度上优于DA+IVTCG,OVTVR也是DA+IVTCG的1.96～3.12倍。

表1 数字鼠各器官的光学参数[27]Tab.1 Optical parameters of various organs at different wavelengths in digital mouse

图1 DA+IVTCG和SP3+IVTCG对中心位于(14,8，16.6)mm处的不同尺寸单光源的重建结果。(a)～(c)是DA+IVTCG的重建结果在真实光源中心z=16.6mm处的截面视图,其中黑色圆圈为真实光源位置。(d)～(f)是SP3+IVTCG对应的重建结果。Fig.1 The reconstruction results of different single sources centered at (14,8，16.6)mm by DA+IVTCG and SP3+IVTCG.(a)～(c) are the section views at z=16.6mm of the reconstruction results by DA+IVTCG, where the black circles denote the real light sources.(d)～(f) are the corresponding reconstruction results by SP3+IVTCG.

图2和图3分别为第一组实验中光源的中心位置在(14,10.5，16.6)mm和(14,12，16.6)mm处的重建结果,图2,图3的结果也和图1一致,以及由表2也可以看出,SP3+IVTCG的重建光源在形状、定位误差、以及平均能量密度上均优于DA+IVTCG。尤其对(14,12，16.6)mm处的光源优势更明显,OVTVR在小光源重建时达到了DA+IVTCG的4.3倍。由于此时光源深度较深,且处于高吸收低散射的肝脏内部,SP3+IVTCG更加适用于这种情况。

第二组实验对比了两种重建方法在高散射低吸收组织内的光源重建能力。将不同尺寸的单光源置于肺中,中心位置设定为(21.2,12.5,8)mm,相对于数字鼠背部表面深度为8.5mm,相对于数字鼠腹部表面深度为9.5mm,图4展示了重建结果的截面视图,表3对比了两种方法得到的具体定量重建结果。尽管这组光源置于散射占主导地位的生物组织中,且距离外部边界大于几个平均自由程[11],属于DA模型较为适用的情形,但从图4以及表3的重建结果可以看出,SP3+IVTCG仍得到了优于DA+IVTCG的重建结果,只是优势没有在高吸收介质中那么明显。

图2 DA+IVTCG和SP3+IVTCG对中心在(14,10.5，16.6)mm处的不同尺寸单光源的重建结果。(a)～(c)是DA+IVTCG的重建结果在真实光源中心z=16.6mm处的截面视图,其中黑色圆圈为真实光源位置。(d)～(f)是SP3+IVTCG对应的重建结果。Fig.2 The reconstruction results of different single sources centered at (14,10.5，16.6)mm by DA+IVTCG and SP3+IVTCG.(a)～(c) are the section views at z=16.6mm of the reconstruction results by DA+IVTCG, where the black circles denote the real light sources.(d)～(f) are the corresponding reconstruction results by SP3+IVTCG.

图3 DA+IVTCG和SP3+IVTCG对中心位于(14,12，16.6)mm处的不同尺寸单光源的重建结果。(a)～(c)是DA+IVTCG的重建结果在真实光源中心z=16.6mm处的截面视图,其中黑色圆圈为真实光源位置。(d)～(f)是SP3+IVTCG对应的重建结果。Fig.3 The reconstruction results of different single sources centered at (14,12，16.6)mm by DA+IVTCG and SP3+IVTCG.(a)～(c) are the section views at z=16.6mm of the reconstruction results by DA+IVTCG, where the black circles denote the real light sources.(d)～(f) are the corresponding reconstruction results by SP3+IVTCG.

MethodActual centerSize/mmRecons. centerLE/mmOVTVRDensity/nW·mm-3r=0.5,h=1(14.28,8.57,16.48)0.651.477.44×10-5(14,8,16.6)r=1,h=1(14.25,8.53,16.26)0.619.154.00×10-4r=1,h=2(14.50,8.22,16.37)0.6017.578.06×10-4r=0.5,h=1(13.22,10.79,16.46)0.851.056.49×10-5DA+IVTCG(14,10.5,16.6)r=1,h=1(13.33,10.73,16.50)0.7111.743.53×10-4r=1,h=2(13.76,10.73,16.19)0.5317.447.36×10-4r=0.5,h=1(13.69,11.72,17.21)0.741.106.20×10-5(14,12,16.6)r=1,h=1(13.65,11.68,17.08)0.678.703.48×10-4r=1,h=2(13.68,11.70,17.20)0.7413.407.34×10-4r=0.5,h=1(14.00,7.85,16.12)0.502.892.35×10-4(14,8,16.6)r=1,h=1(13.82,7.92,16.25)0.4028.601.35×10-3r=1,h=2(13.89,7.93,16.09)0.5238.492.39×10-3r=0.5,h=1(13.50,10.37,16.56)0.521.761.27×10-4SP3+IVTCG(14,10.5,16.6)r=1,h=1(13.89,10.38,16.67)0.1813.796.51×10-4r=1,h=2(13.79,10.45,16.76)0.2922.071.59×10-3r=0.5,h=1(13.87,11.71,16.60)0.324.741.47×10-4(14,12,16.6)r=1,h=1(13.78,11.72,16.52)0.3723.468.14×10-4r=1,h=2(13.61,11.78,16.24)0.5843.621.67×10-3

图4 DA+IVTCG和SP3+IVTCG对中心在(21.2,12.5,8)mm处的不同尺寸单光源的重建结果。(a)～(c)是DA+IVTCG的重建结果在真实光源中心z=8mm处的截面视图,其中黑色圆圈为真实光源位置。(d)～(f)是SP3+IVTCG对应的重建结果。Fig.4 The reconstruction results of different single sources centered at(21.2，12.5，8)mm by DA+IVTCG and SP3+IVTCG.(a)～(c) are the section views at z=8mm of the reconstruction results by DA+IVTCG, where the black circles denote the real light sources.(d)～(f) are the corresponding reconstruction results by SP3+IVTCG.

MethodActual centerSize/mmRecons. centerLE/mmOVTVRDensity/nW·mm-3r=0.5,h=1(21.16,13.36,7.95)0.872.791.03×10-4DA+IVTCG(21.2,12.5,8)r=1,h=1(21.21,13.29,7.95)0.7919.286.22×10-4r=1,h=2(21.16,13.37,7.93)0.8731.691.14×10-3r=0.5,h=1(20.81,12.30,7.49)0.675.601.94×10-4SP3+IVTCG (21.2,12.5,8)r=1,h=1(20.84,12.53,7.53)0.5934.641.09×10-3r=1,h=2(20.87,12.60,7.47)0.7370.482.37×10-3

2.2 双光源仿真实验

为了进一步测试两种重建方法对多目标的分辨能力,本文设计了双光源实验,相同大小的双光源被置于数字鼠的肝脏中,光源边缘间的距离为3mm,光源大小与单光源仿真设置一致。本组实验仍然测试了3种不同尺寸的光源,为了保持双光源边缘距离不变,不同大小光源的中心点位置不同,在光源为r=0.5mm,h=1mm时,双光源的中心点为(14,8,17)mm和(14,12,17)mm,在光源大小为r=1mm,h=1mm以及r=1mm,h=2mm时,双光源的中心点设置为(14,7.5,17)mm和(14,12.5,17)mm。

图5和图6分别展示了3组不同尺寸的双光源重建结果的截面图和3D视图。从图5，图6和表4的重建结果可以看出,SP3+IVTCG在所有的光源设置下,均能得到2个较为均匀的重建目标,中心定位和形状拟合都明显优于DA+IVTCG。尤其是对较大尺寸的光源,SP3+IVTCG比DA+IVTCG更能准确地反应光源的形状。在双大光源时,SP3+IVTCG重建的两个光源的OVTVR分别是DA+IVTCG的1.22和2.46倍。

图5 DA+IVTCG和SP3+IVTCG双光源重建结果在光源中心所在位置的截面视图(z=17mm),其中黑色圆圈为真实光源位置。(a)～(c)是DA+IVTCG的重建结果,(d)～(f)是SP3+IVTCG的重建结果Fig.5 Section views (z=17mm) of the reconstruction results by DA+IVTCG和SP3+IVTCG, where the black circles denote the real light sources.(a)～(c) are the reconstruction results of DA+IVTCG,(d)～(f) are the reconstruction results of SP3+IVTCG.

图6 DA+IVTCG和SP3+IVTCG双光源重建结果的3D视图,其中红色圆柱体为真实光源。(a)～(c)是DA+IVTCG的重建结果,(d)～(f)是SP3+IVTCG的重建结果Fig.6 3D views of the reconstruction result of DA+IVTCG and SP3+IVTCG, where the red cylinder are the real light sources.(a)～(c) are the reconstruction results of DA+IVTCG,(d)～(f) are the corresponding results of SP3+IVTCG

MethodSize/mm#Actual centerRecons. centerLE/mmOVTVR Density/nW·mm-3r=0.5,h=1 S1(14,8,17)(13.95,8.47,16.85)0.492.291.86×10-4S2(14,12,17)(13.76,9.49,16.68)2.541.678.85×10-5DA+IVTCGr=1,h=1 S1(14,7.5,17)(14.59,8.19,16.44)1.0610.901.10×10-3S2(14,12.5,17)(13.74,12.82,16.44)0.697.664.34×10-4r=1,h=2 S1(14,7.5,17)(14.60,8.18,16.52)1.0319.321.82×10-3S2(14,12.5,17)(13.78,12.55,16.33)0.7113.471.12×10-3r=0.5,h=1 S1(14,8,17)(15.53,8.39,16.95)0.612.443.35×10-4S2(14,12,17)(14.03,12.26,17.07)0.273.582.05×10-4SP3+IVTCGr=1,h=1 S1(14,7.5,17)(13.40,8.11,16.83)0.8710.731.90×10-3S2(14,12.5,17)(14.04,12.36,17.00)0.1518.901.57×10-3r=1,h=2 S1(14,7.5,17)(13.98,7.71,16.96)0.2123.663.30×10-3S2(14,12.5,17)(13.99,12.15,17.21)0.4133.203.81×10-3

3 结 论

本文利用多光谱测量数据,结合L1范数正则化IVTCG算法,比较了基于SP3和基于DA的光源重建方法对重建光源形状拟合程度的影响。单光源的仿真实验结果表明,在高散射低吸收的区域,SP3+IVTCG在形状以及重建光源的能量密度上优于DA+IVTCG,在非散射主导的区域内,SP3+IVTCG的优势更加明显,这是由于SP3更加适用于非散射主导的区域,光源深度越深,重建对于光传输模型的准确性的依赖越强,这种优势越明显。双光源实验对比了两种方法的光源分辨能力,SP3+IVTCG重建的两个光源比DA+IVTCG重建的两个光源形状更加近似,平均能量密度更加相当,定位也更加准确。

从仿真实验结果中发现，相对于大尺寸光源,小尺寸光源的形状拟合程度较差。文献[29]的研究表明,对于较小尺寸的光源,基于非凸正则化的重建算法可以得到更为稀疏的解,在未来的工作中,我们拟将SP3模型和非凸稀疏重建算法结合，进一步改善小尺寸光源的形状拟合性能。此外,由于本文的方法结合了多光谱测量数据,所以重建时间较长,未来将结合主成分分析或者GPU硬件加速等，提高方法的效率。