光学学报, 2023, 43 (1): 0111001, 网络出版: 2023-01-06   

校正光通量变化的定量光声内窥成像 下载: 711次

Quantitative Photoacoustic Endoscopic Imaging for Correcting Light Fluence Variation
作者单位
1 华北电力大学电子与通信工程系,河北 保定 071003
2 华北电力大学河北省电力物联网技术重点实验室,河北 保定 071003
摘要
在光声内窥成像中,不均匀或不稳定的照明,以及生物组织的复杂光学特性均会导致成像平面内光通量分布不均匀,从而造成重建图像质量和成像深度的下降。提出了一种校正光通量变化的定量光声内窥成像方法,对光吸收能量分布进行稀疏分解,采用贪婪算法重构得到光吸收系数和光通量的稀疏表示,并通过稀疏矩阵分解实现光吸收系数与光通量分布的联合重建。仿真和体模实验结果表明,与一步法和基于模型的定量重建方法相比,采用所提方法估算光吸收系数的均方根误差(RMSE)可降低约48%,重建图像的归一化平均绝对距离(NMSAD)和结构相似度(SSIM)可分别降低约25%和提高约24%。与其他校正光通量的重建算法相比,所提方法估计光吸收系数的RMSE可降低约22%、NMSAD可降低约20%、SSIM可提高约10%。
Abstract
Objective

Photoacoustic imaging (PAI) is an emerging imaging modality that provides structural and functional information on biological tissue, with the advantages of high contrast of optical imaging and high penetration depth of ultrasound imaging. Photoacoustic endoscopic imaging is an endoscopic application of PAI, which combines non-invasive PAI with endoscopic detection technology. This imaging modality is physically based on the photoacoustic effect of biological tissue illuminated by short laser pulses. The absorbed optical energy density is proportional to the product of the optical absorption coefficient and the local light fluence. Therefore, light fluence is related not only to the optical properties of the tissue but also to the distribution of irradiation intensity. The optical properties of the tissue cannot be accurately reflected in the distribution image of absorption energy density. By solving the optical inverse problem, the optical property parameters can be estimated quantitatively to realize quantitative imaging. Due to the complex optical properties of tissue and non-uniform and non-stationary illumination, photoacoustic image reconstruction is subject to non-uniform light fluence, which leads to a reduction in image quality and imaging depth. The purpose of this work is to solve the problem of reduced accuracy of quantitative imaging due to light fluence variation.

Methods

A quantitative image reconstruction method for photoacoustic endoscopic imaging is presented to correct light fluence variation. It employs a two-step scheme. Firstly, the distribution of absorbed optical energy density on the cross-section of the tubular object is recovered from the photoacoustic signal measured by the ultrasonic detector with conventional image reconstruction methods such as back projection and time reversion. Secondly, the sparse representation of the distribution of absorbed optical energy density is constructed by the weighted sum of a finite number of discrete Curvelet and Haar joint basis functions. The sparse representations of the absorption coefficient and light fluence are then obtained by the greedy algorithm. In addition, the absorption coefficient and light fluence distributions are reconstructed simultaneously by sparse matrix decomposition.

Results and Discussions

The method has been verified by simulation and phantom studies. The comparisons with the one-step method, model-based method, and state-of-the-art fluence compensation method demonstrate the superiority of the proposed method in recovering optical coefficient distribution with high accuracy. The simulations show that the one-step method has a large error when estimating the absorption coefficient, and the reconstructed image cannot distinguish different tissue types clearly. In the images reconstructed by the Bregman method, significant overlap between different tissue regions can be observed. The images reconstructed by the perturbation Monte Carlo method show the low contrast of tissue boundaries in the low-fluence region. The proposed method is superior to other methods in recovering the absorption coefficient closest to the ground truth. In the phantom study, due to similar optical properties of the two materials used to fabricate the phantom, the overall contrast of the image representing the distribution of absorbed optical energy density is low, and it is difficult to distinguish the contours of the embedded targets. In contrast, images showing the distribution of absorption coefficients can clearly distinguish the contours of the targets. In addition, compared with other methods, the proposed method can achieve the highest accuracy in absorption efficient estimation. Simulation and experimental results on the phantom reveal that the root mean square error (RMSE) of the absorption coefficient estimated with this method can be reduced by about 48%, and the normalized mean square absolute distance (NMSAD) and structural similarity (SSIM) of the reconstructed images can be reduced by about 25% and increased by about 24%, respectively, compared with the one-step method and model-based method. The comparisons with the state-of-the-art fluence compensation method show significant improvements in RMSE, NMSAD, and SSIM metrics by about 22%, 20%, and 10%, respectively.

Conclusions

This work presents a joint reconstruction method of the absorption coefficient and light fluence distribution based on sparse decomposition to construct the distribution of absorbed optical energy density. The method shows advantages over the one-step method based on exact solutions, the model-based method, and the state-of-the-art fluence compensation method in recovering absorption coefficient distribution with high accuracy. Furthermore, the quantitative reconstruction accuracy of the proposed method is not sensitive to the recovery of the absorbed optical energy density and the similarity measurement of sparse decomposition. It should be noted that the optimization algorithm for sparse reconstruction influences the reconstruction accuracy, among which the accuracy of multi-candidate orthogonal matching pursuit (MOMP) and Dice orthogonal matching pursuit (DOMP) is higher than that of matching pursuit (MP) and basis pursuit (BP). Future work should involve the following two aspects. One is to verify the feasibility of this method in clinical transplantation through in vivo experiments. The reconstruction accuracy should be further improved with comprehensive consideration of non-ideal factors in practical application scenarios, such as limited-view sparse-sampling measurement, acoustic reflection and scattering, motion artifacts caused by cardiac motion in intracoronary photoacoustic imaging, and characteristics of ultrasonic detectors. The other is to attempt to apply deep learning to eliminate the influence of light fluence variation on PAI quality.

1 引言

光声成像(PAI)是一种基于光声效应的新型复合功能成像方法1。该方法的原理是用短脉冲(ns量级)激光照射组织,组织吸收部分光能量后受热膨胀,激发出宽带(MHz量级)超声波并向组织表面传播,超声换能器接收声压信号后送入计算机,通过求解声学逆问题重建出初始声压分布图或光吸收能量分布图,从而显示出目标的形态结构。其中,常用的方法有反投影2、时间反演(TR)3、傅里叶变换法(FTR)4和反卷积法(DR)5等。由于光吸收能量与光吸收系数和局部光通量的乘积成正比,其中光通量不仅与组织的光学特性有关,还与光强分布有关,因此光吸收能量分布图无法准确反映组织的光学特性。通过求解光学逆问题定量估算组织的光学特性参数(如光吸收系数和散射系数)6、热膨胀系数(Grüneisen 系数)和功能特性参数(如发色团质量分数和血氧饱和度)等,可实现定量成像,获得组织的功能成分信息7

定量光声成像(qPAI)是一个非线性病态问题,一般采用基于模型的最优化方案求解,即建立前向成像物理模型,并迭代求解使声压测量值或光吸收能量重建值与前向模型输出的理论值之间误差最小的非线性最小二乘问题,在每次迭代中更新待估计参数,直到前向模型的输出与测量数据相匹配7-11。实现qPAI的主要困难是光通量未知,为了简化问题,常规方法通常基于理想照明假设和简单光子传输模式,在光源位置固定的情况下,假设入射光在组织表面和内部形成均匀光照。该假设对于较小的成像目标具有一定的合理性,但对于具有复杂结构的大型目标而言,不合理的照明模式会造成光照不均匀,只有靠近和面向入射光的部分组织可被相对均匀地照射。此外,大部分生物组织都是高散射且不透明的混浊介质,组织结构与成分具有高度复杂性和差异性,不同组织成分对光的吸收和散射存在明显差异,也会导致光通量分布不均匀。同时,随着入射光在组织中传播深度的增加,光的穿透能力在减弱,光能量趋于分散,在不同位置处的光通量呈现不同的衰减趋势。上述原因会导致图像质量和成像深度的下降,在低光通量区域中会出现目标模糊和精度降低等问题。

目前解决非均匀光照问题的方法主要包括优化照明设备、校正光通量变化和估计光通量三种方案12。优化照明设备的目的是保证成像目标内各个视角上具有足够的光穿透性,如旋转照明13、定制多模光纤照明14、光捕捉器15和光纤扩散器16等,均需改变成像系统。在图像重建过程中对光通量变化进行校正可提高成像精度,如:Rosenthal等17将图像分解为低空间频率全局分量和高空间频率局域分量,再通过少量基函数的和逼近图像;Bu等18将基于模型的光吸收系数重建与光通量补偿相结合,减少因光通量变化所致的重建误差;Nykänen等19采用基于模型的方法利用初始声压和表面光测量数据同时估计光吸收系数、散射系数和Grüneisen系数,但无法确定最小化问题是否具有唯一最优解。对于光声光谱成像,Deán-Ben等20利用可逆光开关荧光蛋白的波长开关特性校正散射介质内光通量分布的不均匀性,而Zhou等21验证了三种基于模型的光通量校正方法的有效性,即基于一维扩散近似模型、基于扩散偶极子模型和基于蒙特卡罗(MC)模拟的方法。此外,将PAI与其他成像模态相结合,可实现对光通量的估计或测量,如光声-声光(AO)测量22-23、光声-超声成像24、光声-被动超声成像25和光声-扩散光层析成像(DOT)26等。双模态成像方法增加了系统及其操作的复杂度,故其应用范围受限。

光声内窥成像将无创的PAI和内窥检测技术相结合,可用于腔体病变(如消化道病变和冠状动脉粥样硬化等)的早期诊断。提出了一种定量光声内窥成像方法,该方法不依赖于前向光传输模型,而是将Curvelet和Haar联合基与傅里叶基作为基函数,对光吸收能量分布进行稀疏表示,通过矩阵分解实现腔体横截面上光吸收系数和光通量分布图的联合重建,从而校正光通量变化对定量重建精度的影响。

2 原理与方法

图1所示,所提方法基于两步策略:1)采用常规图像重建方法(如反投影或TR)根据超声探测器测量的光声信号重建腔体横截面上的光吸收能量密度(AOED)分布图;2)建立光吸收能量的稀疏表达式,得到光吸收系数和光通量的稀疏表示,并通过稀疏矩阵分解实现光吸收系数与光通量分布的联合重建。

图 1. 联合重建光吸收系数和光通量分布图的流程

Fig. 1. Process for simultaneously reconstructing optical absorption coefficient and optical fluence distribution

下载图片 查看所有图片

在PAI中,光吸收能量、光吸收系数和光通量之间的关系为

Ar=μarΦr

式中:ArμarΦr分别为成像平面中位置r处的光吸收能量、光吸收系数和光通量。成像平面极坐标系中光吸收能量的稀疏表达式为

lgAθ,l=lgμaθ,l+lgΦθ,l=n=1Ncnϕnθ,l+m=1Mdmψmθ,l

式中:θ,l为位置r的极坐标,如图2所示;ϕnθ,lψmθ,l分别为lgμaθ,llgΦθ,l稀疏表示的基函数;cndm为稀疏分解系数,其中n=1,2,,Nm=1,2,,M

图 2. 成像平面极坐标系

Fig. 2. Polar coordinate system of imaging plane

下载图片 查看所有图片

基函数需满足稀疏条件,即每个基函数只能单一地表示lgμaθ,llgΦθ,l。由于光吸收系数和光通量分别满足局部高频率和全局低频率特性17,因此本文采用离散Curvelet和Haar联合基27-28表示lgμaθ,l,采用离散傅里叶基29表示lgΦθ,l,构成基函数库。此时,光吸收能量的稀疏分解问题就转化成了寻找满足稀疏表达式的、最少数量系数cndm的问题。定义优化问题为

Y=minαkαk0,s.t. Xk-Ckαk2<εk

式中:k=1,2,,360为成像平面极坐标系中的第k个角度;αk=c1cNd1dMT为由分解系数cndm构成的列向量;εk为重构误差阈值;02分别为零范数和2范数;Ck=Cϕn,kCψm,kL×(N+M)维的稀疏矩阵,其各列分量分别对应基函数库ϕnψmL为一个角度上离散采样点的个数,第L个采样点对应的极径为LΔl(单位为mm),Δl为一个角度上相邻采样点之间的径向距离;Xk为由第k个角度上的光吸收能量对数值组成的列向量。

图 3. 矩阵分解示意图

Fig. 3. Schematic diagram of matrix decomposition

下载图片 查看所有图片

式(3)是非凸优化问题,常用的求解方法有凸松弛算法和贪婪算法30:凸松弛算法是将原目标函数中难以解决的部分用松弛变量代替,使目标函数成为凸函数,如基追踪(BP)算法31;贪婪算法是将求解的问题分成若干子问题,求解子问题的局部最优解并进行合成,得到当前条件下原问题整体最优解的近似解。与凸松弛算法相比,贪婪算法的复杂度低、运算量小、重构效率高且相对稳定,故其在解决稀疏分解问题中的应用最为广泛30。贪婪算法的典型代表是匹配追踪(MP)算法32,它采用迭代的方式根据基函数库重构图像,每次迭代时将基函数库中的一个元素添加到图像的稀疏表示中,当重构精度满足要求或超过预先设定的迭代次数时即停止迭代。由于每个库元素的系数只选择一次,因此无法在当前迭代步骤中修正上一次迭代中的误差。MP算法的改进算法是正交匹配追踪(OMP)算法30,每次迭代都会重新计算所有基函数的最佳系数,可以补偿之前迭代步骤中的误差。同时,OMP算法的搜索过程具有较好的随机性,可有效避免陷入局部最优,但每次迭代只能选择基函数库中的一个元素,故效率较低。OMP算法的改进算法有多候选正交匹配追踪(MOMP)33和基于Dice系数的正交匹配追踪(DOMP)算法34等,这些算法不仅保留了OMP的优点,还具有更高的重构精度和效率。因此,本文采用DOMP算法求解式(3)的优化问题,在第4.1节中将结合仿真实验结果讨论其相对于其他非凸优化求解算法的优势。

DOMP算法的原理是:在每一次迭代中,根据Dice系数匹配准则求解Dice值最大的一个元素,并加入到稀疏矩阵中,更新残差判定重构结果,达到精度要求或允许的迭代次数时停止迭代。求解式(3)的具体步骤如下。

第一步,初始化迭代次数i=1,设置初始稀疏矩阵Ck,0和稀疏分解系数向量αk,0均为空矩阵,则Xk的初始值Xk,0

Xk,0=Ck,0αk,0=

残差Rk,0

Rk,0=Xk-Xk,0=Xk

第二步,计算第i次迭代后的稀疏矩阵。根据Dice系数匹配准则在基函数矩阵库中搜索与残差相似程度最大的库函数作为最匹配的库函数,即

ϕi=argmaxϕDiceϕ,Rk,i-1ψi=argmaxψDiceψ,Rk,i-1

式中:ϕψlgμak,ilgΦk,i的库函数;ϕiψi为第i次迭代后搜索出的最匹配库函数;Rk,i-1为第i-1次迭代的残差;Dice·为Dice系数匹配准则函数,其表达式为

Dicex,y=2i=1Pxi2+yi2i=1Pxiyi

式中: xiyi分别为向量x和向量y中的元素;P为向量中的元素个数。更新稀疏矩阵为

Ck,i=Cϕn,k,i-1ϕiCψm,k,i-1ψi

式中:Cϕn,k,i-1Cψm,k,i-1分别为Ck,i-1中对应基函数库ϕnψm的矩阵。

第三步,计算第i次迭代后的稀疏分解系数矩阵

αk,i=Ck,i+Xk=Ck,iHCk,i-1Ck,iHXk

式中:Ck,i+Ck,iH分别为Ck,i的Moore-Penrose伪逆矩阵和共轭转置矩阵。

第四步,计算第i次迭代后的残差

Rk,i=Xk-Xk,i=Xk-Ck,iαk,i

第五步,判断迭代是否终止。如果第i次迭代后的误差满足Rk,i-Rk,i-120.01Xk2,则终止迭代,输出稀疏矩阵Ck,i和稀疏分解系数矩阵αk,i。否则,更新迭代次数ii+1,并转向第二步。

按照上述步骤求解出稀疏矩阵和稀疏分解系数矩阵,进而得到光吸收系数和光通量对数的向量表示

lg μa,k=Cϕn,kαϕn,klg Φk=Cψm,kαψm,k

式中:Cϕn,kCψm,k分别为稀疏矩阵Ck中与ϕnψm对应的部分;αϕn,kαψm,k分别为系数矩阵αk中与cndm对应的部分;μa,kΦk分别为角度k上的光吸收系数和光通量重建值。

3 结果及分析

由于进行生物组织在体实验测量时,无法获知组织光学特性参数的真实值,因此本文分别采用计算机仿真模型和体模成像数据验证所提方法的可行性并评估其性能。实验中的编程环境是MATLAB R2018a,所用计算机的CPU型号为Intel® CoreTM i7-5500,其主频为2.40 GHz、内存为4 GB,所用计算机的操作系统为Windows 10 64位操作系统。

3.1 仿真实验结果

分别构建了简单点阵仿真模型和包含不同组织成分的冠状动脉血管仿真模型,如图4(a)所示,点阵模型包含25个半径为0.03 mm的点目标,按照内层、中层和外层的顺序排列,光吸收系数分别设定为1.0 cm-1、0.6 cm-1和0.4 cm-1。如图4(b)所示,参照人体冠脉组织切片图建立包含健康管壁组织和不同成分粥样硬化斑块的冠脉血管计算机仿真模型,并参考人体动脉组织学设定各组织成分的光学和声学特性参数35-37。前向仿真中采用MC方法模拟光子在组织中的传输,入射光波长设定为800 nm,光子总数为106,采用MATLAB k-wave工具箱仿真光声信号38。图像平面的Δl设置为9.8×10-2 mm,L设置为180。

图 4. 仿真点阵模型和仿真血管模型的几何结构。(a)仿真点阵模型;(b)仿真血管模型

Fig. 4. Geometry of simulated grid model and simulated vascular model. (a) Simulated grid model; (b) simulated vascular model

下载图片 查看所有图片

分别采用所提方法、文献[6]中的方法、文献[9]中的方法和文献[11]中的方法对仿真模型进行光吸收系数图像重建,结果如图5所示。文献[6]中的方法基于波动方程的精确解从测量的光声信号中直接重建光吸收系数,下文将其简称为一步法。文献[9]和文献[11]是基于模型的方法:文献[9]采用辐射传输方程的扩散近似作为前向光模型,并采用Bregman算法求解优化问题,下文将其简称为Bregman法;文献[11]利用扰动MC(PMC)作为前向光模型,下文将其简称为PMC法。可以发现:一步法重建的图像中无法清晰辨别各组织区域的形状,光吸收系数估算误差较大;Bregman法重建的图像中虽然可以定位斑块,但是不同成分组织区域之间存在重叠,斑块边界处误差较明显;PMC法重建的图像中可以区分斑块形状,但在斑块区域中存在明显误差,且低光通量区域(远离导管的区域)中组织边界的对比度较低;所提方法的重建图像中组织边界最清晰,与真实值最接近。

图 5. 仿真模型的图像结果。(a)点阵模型;(b)血管模型

Fig. 5. Image results of simulation models. (a) Grid model; (b) vascular model

下载图片 查看所有图片

除视觉评估以外,本文采用均方根误差(RMSE)、归一化平均绝对距离(NMSAD)和结构相似度(SSIM)作为量化指标,来客观评价光吸收系数估算精度和重建图像质量。RMSE是衡量估计值与真实值之间偏差的量化指标,定义为

ERMS=1WHx=1Wy=1H[μ^a(x,y)-μa(x,y)]2

式中:WH分别为图像的宽度和高度(单位为pixel);μ^a(x,y)μa(x,y)分别为图像平面中像素(x y)处光吸收系数的估计值和真实值。SSIM和NMSAD是定量评价重建图像与标准图像之间相似度的指标39,分别定义为

SS(G,R)=(2μGμR+c1)(2σGR+c2)(μG2+μR2+c1)(σG2+σR2+c2)DNMSA(G,R)=1x=1Wy=1HR(x,y)x=1Wy=1HR(x,y)-G(x,y)

式中:GR分别为标准图像和重建图像;μGμR分别为图像G和图像R的平均亮度;σGσR分别为图像G和图像R的亮度标准差;σGR为图像G和图像R的亮度协方差;c1c2为极小正数,用于增加计算的稳定性;Gx y)和Rx y)分别为图像G和图像R中像素(x y)处的亮度值。由图6可知,与一步法和基于模型的方法相比,所提方法估算光吸收系数的RMSE可降低约48%,重建图像的NMSAD可降低约25%,SSIM可提高约24%。

图 6. 仿真模型的光吸收系数估算精度和重建图像定量评价指标。(a) RMSE;(b) NMSAD;(c) SSIM

Fig. 6. Quantitative evaluation metrics for estimation accuracy of optical absorption coefficient and reconstructed images of simulation models. (a) RMSE; (b) NMSAD; (c) SSIM

下载图片 查看所有图片

3.2 体模实验结果

采用3D打印技术制作实体模型(简称体模),如图7所示,体模总长为85 mm,横截面直径为18 mm,内含三个嵌入物以模拟不同组织成分。嵌入物1贯穿整个体模且位于横截面中心,其长、宽、高分别为13、6、85 mm。嵌入物2和3分别位于嵌入物1的下方和上方,且形状和尺寸相同(长、宽、高分别为3.5、3.5、50.0 mm)。体模材料是软度为30的白色硅胶,嵌入物的材料是白色树脂,两种材料的光吸收系数分别为0.3 cm-1和0.2 cm-1

图 7. 体模几何结构示意图和实物照片

Fig. 7. Schematic diagram and photo of geometric structure of actual model

下载图片 查看所有图片

利用MSOT inVision 256多光谱光声断层扫描成像系统(iThera Medical GmbH)采集体模的光声数据。该系统由Nd∶YAG近红外激光器泵浦的可调谐光学参量振荡器组成,入射光波长是680~900 nm,重复频率是10 Hz,激光脉冲持续时间是7 ns,最大脉冲能量是120 mJ,激光切换速度是100 ms。采用270°弧形阵列式超声探测装置,共包含256个换能器阵元,探测器中心频率是5 MHz,带宽为60%。图7中体模的截面1和截面2到探测器的距离分别是13 mm和49 mm。图8展示了入射光波长为680 nm时成像系统输出的两个截面的光吸收能量分布图像、重建的光流图,以及采用所提方法、一步法、PMC法和Bregman法重建的光吸收系数分布图。可以发现:由于两种材料的光吸收系数差别不大,故光吸收能量分布图的整体对比度较低,较难区分目标轮廓;在光吸收系数分布图中,可分辨出三个嵌入物的轮廓,且与其他方法相比,所提方法重建的光吸收系数分布图中嵌入物的轮廓最清晰,图像质量与对比度均得到了明显提高。图9中的定量评价结果进一步验证了所提方法的重建精度,与一步法和基于模型的方法相比,所提方法估算光吸收系数的RMSE可降低约33%、重建图像的NMSAD可降低约20%、SSIM可提高约19%。

图 8. 体模成像结果

Fig. 8. Image results of actual model

下载图片 查看所有图片

图 9. 体模光吸收系数估算精度和重建图像定量评价指标。(a) RMSE;(b) NMSAD;(c) SSIM

Fig. 9. Quantitative evaluation metrics for estimation accuracy of optical absorption coefficient and reconstructed images of actual model. (a) RMSE; (b) NMSAD; (c) SSIM

下载图片 查看所有图片

4 讨论

本章结合仿真血管模型图像重建结果和定量评价指标,讨论稀疏重构中的相似度测量方法、优化算法和光吸收能量重建方法对光吸收系数重建精度和图像质量的影响,并与其他联合重建方法进行对比。

4.1 影响重建精度的主要因素

4.1.1 相似度测量方法

DOMP是一种根据Dice系数匹配准则计算相似度的贪婪算法。为了讨论所提方法对稀疏重构相似度测量的依赖关系,本节分别采用Dice准则、内积法(IP)29、Jaccard系数法40和余弦相似法(COS法)41求解向量相似性,其余条件不变,重建的光吸收系数分布图和评价指标分别如图10(a)和图11(a)所示。可以发现,当采用不同相似度测量方法时,重建图像的视觉效果和评价指标均无明显差异。出现该现象的主要原因在于虽然采用不同相似度测量的判定结果存在差异,但是当基函数的数量逐渐增加时,这种差异会逐渐缩小。

图 10. 采用不同相似度测量法、优化算法和AOED重建方法时重建的光吸收系数分布图。(a)不同相似度测量法得到的结果;(b)不同优化算法得到的结果;(c)不同AOED重建方法得到的结果

Fig. 10. Optical absorption coefficient diagrams reconstructed by different similarity measurement methods, optimization algorithms and AOED recovery methods. (a) Results obtained by different similarity measurement methods; (b) results obtained by different optimization algorithms; (c) results obtained by different AOED recovery methods

下载图片 查看所有图片

图 11. 采用不同相似度测量法、优化算法和AOED重建方法时重建的光吸收系数分布图的定量评价指标。(a)相似度测量法;(b)优化算法;(c) AOED重建方法

Fig. 11. Quantitative evaluation metrics of optical absorption coefficient diagrams reconstructed by different similarity measurement methods, optimization algorithms and AOED recovery methods. (a) Similarity measurement methods; (b) optimization algorithms; (c) AOED recovery methods

下载图片 查看所有图片

4.1.2 优化算法

分别采用BP31、MP32、MOMP33和DOMP34算法求解式(3)中的优化问题,讨论稀疏重构优化算法对光吸收系数估算精度的影响,其中MP、MOMP和DOMP都属于贪婪算法,BP属于凸松弛算法。光吸收系数分布图重建结果和定量评价指标分别如图10(b)和图11(b)所示。可以发现:贪婪算法的重建精度略优于凸松弛算法;不同优化算法的重建时间存在明显差异,即MOMP耗时最短、BP最长。出现该现象的原因是MOMP每次迭代时选择多个最优元素组成最优集,有效缩短了迭代时间,而其他算法每次迭代时仅选择单个元素。除此之外,MOMP和DOMP每次迭代时会根据相似性选择高匹配度的元素,且迭代中可以实现对前一次迭代误差的有效纠正,故重建图像质量优于其他算法。综合考虑重建质量和时间成本,采用DOMP或MOMP可获得更佳的成像效果。

4.1.3 光吸收能量重建方法

如第2章所述,所提方法是基于两步策略的,故本节讨论光吸收能量分布图的重建质量对光吸收系数重建精度的影响。分别采用反投影2、TR3、FTR4和DR5根据光声信号重建光吸收能量分布图,其余条件不变,最终得到的光吸收系数分布图和定量评价指标如图10(c)和图11(c)所示。可以发现,当采用不同光吸收能量图像重建方法时,光吸收系数图像的视觉效果和评价指标均无明显差异,其微小差异主要来源于光吸收能量稀疏分解过程中引入了与光吸收能量重建结果相关的误差。

4.2 与其他联合重建方法对比

比较所提方法和Rosenthal方法17重建光吸收系数的精度,结果如图12所示。两种方法都是基于稀疏分解的思想,不同之处在于所提方法在原有库函数的基础上引入了Curvelet基函数,将离散Curvelet和Haar小波联合基用于光吸收系数的稀疏表示。这是因为小波函数在复杂边缘问题上具有一定局限性,而Curvelet函数不仅可以有效解决小波函数造成的边界扭曲问题,对曲线和微弱线性结构的重构效果还优于各种小波函数,更适用于具有复杂形态结构的生物组织42。此外,Rosenthal方法采用OMP算法求解式(3)中的优化问题,所提方法则采用了精确度更高的DOMP算法。由图12(a)可知,所提方法重建图像的视觉效果优于Rosenthal方法,图像中的组织轮廓更清晰,更容易区分不同组织成分。由图12(b)可知,与Rosenthal方法相比,所提方法估算光吸收系数的RMSE可降低约22%,重建图像的NMSAD可降低约20%,SSIM可提高约10%。

图 12. 采用不同联合重建方法得到的血管模型的光吸收系数分布图和评价指标。(a)重建图像;(b)评价指标

Fig. 12. Optical absorption coefficient diagrams and evaluation metrics of vascular model obtained by different joint reconstruction methods. (a) Reconstructed images; (b) evaluation metrics

下载图片 查看所有图片

4.3 所提方法相对于其他定量成像方案的优势和不足

目前qPAI的主流方案是基于模型的方法,其本质是迭代求解前向问题,即通过使前向成像模型计算的光吸收能量理论值与根据光声信号重建的光吸收能量测量值之间的误差最小化来获得对未知量的最佳估计。所提方案的重建精度取决于前向模型对成像过程描述的准确性,其主要困难是建立用简单方程表示的前向模型来描述光在非均匀散射介质中的传输、热膨胀产生的声压、宽带超声波在声学特性非均匀组织中的传播和超声换能器对光声信号的测量过程。为了使模型更全面、更准确地描述实际成像场景,提高重建精度,应综合考虑各种不理想因素,如非稳定/非均匀的光覆盖、组织随空间变化的声速和密度、超声波在不同类型组织界面处的反射和散射,以及探测器的有限孔径效应、有限带宽效应和有限角度稀疏采样等。由于在迭代过程中需要反复计算前向算子及其伴随算子,故定量成像的计算复杂度取决于前向模型描述成像过程的全面程度,在实际应用中必须在时间成本和重建精度之间进行权衡。为了提高重建效率,通常需要对成像场景进行一些理想化的假设,如组织光学特性分段恒定、背景光学特性均匀分布、光散射系数和Grüneisen系数已知,以及组织中的声速恒定等。然而,在现实应用场景中这些假设多数都不成立,故定量成像精度会下降。

与基于模型的方案不同,所提方法通过对AOED进行稀疏分解估算光吸收系数,可以避免在反复求解前向模型过程中引入的误差,无需依赖前向模型即可实现组织光吸收系数的定量估计,估算结果对成像场景、测量几何和造成光通量不均匀的原因不敏感,且方法的实施不仅限于内窥式的应用中。不足之处在于,所提方法采用两步策略,即先根据探测器采集的光声信号重建光吸收能量分布图(结构成像),再在此基础上估算光吸收系数,从而实现定量成像。两步策略的共同缺点是:光吸收能量的重建误差会累积到光吸收系数的估算结果中,从而间接影响估算精度。特别是,在光声显微43和光声断层成像中,由于激光需要穿过皮肤甚至颅骨,因此组织对于光的传播特性影响更大,这种影响会直接体现在光吸收能量分布图的重建中。在体模实验中,利用光声断层扫描成像系统从外部扫描体模,并根据系统输出的光吸收能量分布图估算光吸收系数。图8中的两个成像截面到探测器的距离分别是13 mm和49 mm,模拟了激光从外部扫描成像目标,对其内部结构进行成像的情况。从图8的结果可以看出,采用所提方法估算光吸收系数的精度高于一步法和基于模型的方法,一定程度上验证了所提方法应用于光声断层成像中的可行性。后续工作中拟采用小鼠的全身光声扫描图像进行实验,进一步验证所提方法的可移植性。

近年来,随着深度学习技术在医学成像领域中应用的日益广泛,基于深度学习的PAI技术也取得了突破性进展44。与解析法和迭代法不同,深度学习以数据驱动的方式实现图像重建,利用大量样本数据对复杂成像模型的参数进行监督或无监督学习,实现对多物理场耦合成像过程的建模,在实现快速、高质量图像重建方面展现出了巨大的潜力45。在过去4年中,深度学习在qPAI领域中也表现出了优越的性能,如:采用深度神经网络根据不同波长入射光照射下获得的初始声压或光吸收能量数据估计发色团的质量浓度或血氧饱和度46;采用双路径卷积神经网络在没有标记实验数据的情况下重建深层组织的光吸收系数图像47。这些方法在理想的计算机模拟环境或特定数据集中均取得了良好的效果,但关于入射激光、探头设计、校准参数或散射特性的理想假设可能会限制它们在常规PAI中的实际应用,其有效性还需要通过在体动物或临床试验来进一步证明。

5 结论

设计了一种基于光吸收能量稀疏分解的光吸收系数与光通量分布图联合重建方法,用有限数量的基函数之和构造光吸收能量分布的稀疏表达式,并通过稀疏矩阵分解得到光吸收系数和光通量的稀疏表示,从而实现二者的联合重建。仿真和体模实验结果表明,采用所提方法重建的光吸收系数分布图的视觉效果和评价指标均优于基于精确解的一步法、基于模型的方法和Rosenthal方法。讨论结果表明,所提方法的定量重建精度对光吸收能量重建值和稀疏分解相似度测量方法不敏感,即其具有较好的稳定性。稀疏重构优化算法对重建精度有一定的影响,其中MOMP和DOMP相对于MP和BP的精度更高。

在未来的工作中,将通过动物在体实验验证所提方法的临床可移植性,并综合考虑实际应用场景中的多种不理想因素,如有限角度稀疏采样测量、组织非均匀声学特性所致的声反射和散射、成像目标的运动状态(如冠状动脉内PAI中由心脏运动所致的伪影)和超声探测器的非理想特性等,进一步提高重建精度。此外,深度学习等人工智能技术在PAI领域中的应用日益广泛,也为解决光照不均匀问题提供了一种新思路。

参考文献

[1] 李娇, 李帅, 陈冀景, 等. 非接触光声成像研究进展及其在生物医学上的应用[J]. 中国激光, 2021, 48(19): 1907001.

    Li J, Li S, Chen J J, et al. Progress and biomedical application of non-contact photoacoustic imaging[J]. Chinese Journal of Lasers, 2021, 48(19): 1907001.

[2] Xu M H, Wang L V. Universal back-projection algorithm for photoacoustic computed tomography[J]. Physical Review E, 2005, 71(1): 016706.

[3] Sun Z, Han D D, Yuan Y. 2-D image reconstruction of photoacoustic endoscopic imaging based on time-reversal[J]. Computers in Biology and Medicine, 2016, 76: 60-68.

[4] Wang K, Anastasio M A. A simple Fourier transform-based reconstruction formula for photoacoustic computed tomography with a circular or spherical measurement geometry[J]. Physics in Medicine and Biology, 2012, 57(23): N493-N499.

[5] Zhang C, Wang Y Y. Deconvolution reconstruction of full-view and limited-view photoacoustic tomography: a simulation study[J]. Journal of the Optical Society of America A, 2008, 25(10): 2436-2443.

[6] Wang Y, Wang R K. Photoacoustic recovery of an absolute optical absorption coefficient with an exact solution of a wave equation[J]. Physics in Medicine and Biology, 2008, 53(21): 6167-6177.

[7] Cox B T, Laufer J G, Arridge S R, et al. Quantitative spectroscopic photoacoustic imaging: a review[J]. Journal of Biomedical Optics, 2012, 17(6): 061202.

[8] Mohajerani P, Tzoumas S, Rosenthal A, et al. Optical and optoacoustic model-based tomography: theory and current challenges for deep tissue imaging of optical contrast[J]. IEEE Signal Processing Magazine, 2015, 32(1): 88-100.

[9] Zheng S, Lan Z. Reconstruction of optical absorption coefficient distribution in intravascular photoacoustic imaging[J]. Computers in Biology and Medicine, 2018, 97: 37-49.

[10] Javaherian A, Holman S. Direct quantitative photoacoustic tomography for realistic acoustic media[J]. Inverse Problems, 2019, 35(8): 084004.

[11] Leino A A, Lunttila T, Mozumder M, et al. Perturbation Monte Carlo method for quantitative photoacoustic tomography[J]. IEEE Transactions on Medical Imaging, 2020, 39(10): 2985-2995.

[12] 孟琪, 孙正. 生物光声层析成像中不均匀和不稳定照明解决方法[J]. 中国光学, 2021, 14(2): 307-319.

    Meng Q, Sun Z. Solutions to inhomogeneous and unstable illumination in biological photoacoustic tomography[J]. Chinese Optics, 2021, 14(2): 307-319.

[13] Lou Y, Wang K, Oraevsky A A, et al. Impact of nonstationary optical illumination on image reconstruction in optoacoustic tomography[J]. Journal of the Optical Society of America A, 2016, 33(12): 2333-2347.

[14] Mc Larney B, Rebling J, Chen Z Y, et al. Uniform light delivery in volumetric optoacoustic tomography[J]. Journal of Biophotonics, 2019, 12(6): e201800387.

[15] Yu J, Schuman J S, Lee J K, et al. A light illumination enhancement device for photoacoustic imaging: in vivo animal study[J]. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 2017, 64(8): 1205-1211.

[16] Li M C, Lan B X, Liu W, et al. Internal-illumination photoacoustic computed tomography[J]. Journal of Biomedical Optics, 2018, 23(3): 030506.

[17] Rosenthal A, Razansky D, Ntziachristos V. Quantitative optoacoustic signal extraction using sparse signal representation[J]. IEEE Transactions on Medical Imaging, 2009, 28(12): 1997-2006.

[18] Bu S H, Liu Z B, Shiina T, et al. Model-based reconstruction integrated with fluence compensation for photoacoustic tomography[J]. IEEE Transactions on Biomedical Engineering, 2012, 59(5): 1354-1363.

[19] Nykänen O, Pulkkinen A, Tarvainen T. Quantitative photoacoustic tomography augmented with surface light measurements[J]. Biomedical Optics Express, 2017, 8(10): 4380-4395.

[20] Deán-Ben X L, Stiel A C, Jiang Y Y, et al. Light fluence normalization in turbid tissues via temporally unmixed multispectral optoacoustic tomography[J]. Optics Letters, 2015, 40(20): 4691-4694.

[21] Zhou X W, Akhlaghi N, Wear K A, et al. Evaluation of fluence correction algorithms in multispectral photoacoustic imaging[J]. Photoacoustics, 2020, 19: 100181.

[22] Hussain A, Petersen W, Staley J, et al. Quantitative blood oxygen saturation imaging using combined photoacoustics and acousto-optics[J]. Optics Letters, 2016, 41(8): 1720-1723.

[23] Daoudi K, Hussain A, Hondebrink E, et al. Correcting photoacoustic signals for fluence variations using acousto-optic modulation[J]. Optics Express, 2012, 20(13): 14117-14129.

[24] Zhao L Y, Yang M, Jiang Y X, et al. Optical fluence compensation for handheld photoacoustic probe: an in vivo human study case[J]. Journal of Innovative Optical Health Sciences, 2017, 10(4): 1740002.

[25] Jin H R, Zhang R C, Liu S Y, et al. A single sensor dual-modality photoacoustic fusion imaging for compensation of light fluence variation[J]. IEEE Transactions on Biomedical Engineering, 2019, 66(6): 1810-1813.

[26] Bauer A Q, Nothdurft R E, Culver J P, et al. Quantitative photoacoustic imaging: correcting for heterogeneous light fluence distributions using diffuse optical tomography[J]. Journal of Biomedical Optics, 2011, 16(9): 096016.

[27] Candès E J, Demanet L. The curvelet representation of wave propagators is optimally sparse[J]. Communications on Pure and Applied Mathematics, 2005, 58(11): 1472-1528.

[28] Grohs P, Keiper S, Kutyniok G, et al. Cartoon approximation with α‑curvelets[J]. Journal of Fourier Analysis and Applications, 2016, 22(6): 1235-1293.

[29] Zwartjes P M, Gisolf A. Fourier reconstruction with sparse inversion[J]. Geophysical Prospecting, 2007, 55(2): 199-221.

[30] 田金鹏, 刘小娟, 刘燕平, 等. 多候选集广义正交匹配追踪算法[J]. 应用科学学报, 2017, 35(2): 233-243.

    Tian J P, Liu X J, Liu Y P, et al. Multi-candidate set of generalized orthogonal matching pursuit algorithm[J]. Journal of Applied Sciences, 2017, 35(2): 233-243.

[31] 张繁昌, 李传辉. 非平稳地震信号匹配追踪时频分析[J]. 物探与化探, 2011, 35(4): 546-552.

    Zhang F C, Li C H. Matching pursuit time-frequency analysis of non-stationary seismic signals[J]. Geophysical and Geochemical Exploration, 2011, 35(4): 546-552.

[32] Huang S S, Zhu J B. Recovery of sparse signals using OMP and its variants: convergence analysis based on RIP[J]. Inverse Problems, 2011, 27(3): 035003.

[33] Tsaig Y, Donoho D L. Extensions of compressed sensing[J]. Signal Processing, 2006, 86(3): 549-571.

[34] HuangJ J, XuY H, ZhuP, et al. An improved reconstruction algorithm based on multi-candidate orthogonal matching pursuit algorithm[C]//2014 Seventh International Symposium on Computational Intelligence and Design, December 13-14, 2014, Hangzhou, China. New York: IEEE Press, 2014: 564-568.

[35] 陈荣, 谢树森, 陈艳娇, 等. 中国人血液的组织光学参数[J]. 光电子·激光, 2002, 13(1): 92-93, 97.

    Chen R, Xie S S, Chen Y J, et al. Optical parameters of Chinese blood[J]. Journal of Optoelectronics·Laser, 2002, 13(1): 92-93, 97.

[36] Jacques S L. Optical properties of biological tissues: a review[J]. Physics in Medicine and Biology, 2013, 58(11): R37-R61.

[37] 高上凯. 医学成像系统[M]. 2版. 北京: 清华大学出版社, 2010: 88-95.

    GaoS K. Medical imaging system[M]. 2nd ed. Beijing: Tsinghua University Press, 2010: 88-95.

[38] Treeby B E, Cox B T. k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields[J]. Journal of Biomedical Optics, 2010, 15(2): 021314.

[39] 庄天戈. CT原理与算法[M]. 上海: 上海交通大学出版社, 1992: 11-12.

    ZhuangT G. Principle and algorithm of CT[M]. Shanghai: Shanghai Jiao Tong University Press, 1992: 11-12.

[40] JanaC, KaraaslanF. Dice and Jaccard similarity measures based on expected intervals of trapezoidal neutrosophic fuzzy numbers and their applications in multicriteria decision making[M]//Smarandache F, Abdel-Basset M. Optimization theory based on Neutrosophic and Plithogenic sets. Amsterdam: Elsevier, 2020: 261-287.

[41] Ye J. Cosine similarity measures for intuitionistic fuzzy sets and their applications[J]. Mathematical and Computer Modelling, 2011, 53(1/2): 91-97.

[42] 隆刚, 肖磊, 陈学佺. Curvelet变换在图像处理中的应用综述[J]. 计算机研究与发展, 2005, 42(8): 1331-1337.

    Long G, Xiao L, Chen X Q. Overview of the applications of curvelet transform in image processing[J]. Journal of Computer Research and Development, 2005, 42(8): 1331-1337.

[43] 何勇, 廖唐云, 吴俊伟, 等. 基于透明超声换能器的光声显微镜设计[J]. 中国激光, 2022, 49(3): 0307001.

    He Y, Liao T Y, Wu J W, et al. Design of photoacoustic microscope based on transparent ultrasonic transducer[J]. Chinese Journal of Lasers, 2022, 49(3): 0307001.

[44] Yang C C, Lan H R, Gao F, et al. Review of deep learning for photoacoustic imaging[J]. Photoacoustics, 2021, 21: 100215.

[45] 沈康, 刘松德, 施钧辉, 等. 基于双域神经网络的稀疏视角光声图像重建[J]. 中国激光, 2022, 49(5): 0507208.

    Shen K, Liu S D, Shi J H, et al. Dual-domain neural network for sparse-view photoacoustic image reconstruction[J]. Chinese Journal of Lasers, 2022, 49(5): 0507208.

[46] YangC C, LanH R, ZhongH T, et al. Quantitative photoacoustic blood oxygenation imaging using deep residual and recurrent neural network[C]//2019 IEEE 16th International Symposium on Biomedical Imaging, April 8-11, 2019, Venice, Italy. New York: IEEE Press, 2019: 741-744.

[47] Li J, Wang C, Chen T T, et al. Deep learning-based quantitative optoacoustic tomography of deep tissues in the absence of labeled experimental data[J]. Optica, 2022, 9(1): 32-41.

孟琪, 孙正, 侯英飒, 孙美晨. 校正光通量变化的定量光声内窥成像[J]. 光学学报, 2023, 43(1): 0111001. Qi Meng, Zheng Sun, Yingsa Hou, Meichen Sun. Quantitative Photoacoustic Endoscopic Imaging for Correcting Light Fluence Variation[J]. Acta Optica Sinica, 2023, 43(1): 0111001.

本文已被 1 篇论文引用
被引统计数据来源于中国光学期刊网
引用该论文: TXT   |   EndNote

相关论文

加载中...

关于本站 Cookie 的使用提示

中国光学期刊网使用基于 cookie 的技术来更好地为您提供各项服务,点击此处了解我们的隐私策略。 如您需继续使用本网站,请您授权我们使用本地 cookie 来保存部分信息。
全站搜索
您最值得信赖的光电行业旗舰网络服务平台!