基于Lanczos双对角化的快速光声成像重建方法

张娜; 冯金超; 李哲; 贾克斌

doi:doi:10.3788/CJL201845.0307018

中国激光, 2018, 45 (3): 0307018, 网络出版: 2018-03-06

基于Lanczos双对角化的快速光声成像重建方法下载： 907次

Fast Photoacoustic Imaging Reconstruction Method Based on Lanczos Double Diagonalization

论文大纲

张娜 ^1,*冯金超 ^1,2李哲 ^1,2贾克斌 ^1,2,3

作者单位

¹ 北京工业大学信息学部，北京 100124

² 先进信息网络北京实验室，北京 100124

³ 北京工业大学未来网络科技高精尖创新中心，北京 100124

AI 词云图 AI一句话精读 AI短摘要

注：本部分内容由 AI 自动生成，请您知悉。

摘要

光声成像结合了光学成像和声成像的优点，是一种具有高空间分辨率、高对比度的无损成像技术，成为当前生物医学成像的研究热点之一。重建光声图像是一个典型的逆问题，具有严重的病态性。针对光声成像的病态性和较大的系统矩阵会导致重建速度慢的问题，提出了一种基于Lanczos双对角化的快速指数滤波重建方法，并通过数值仿真证实了该方法的有效性。仿真结果表明，所提方法在保证重建图像高质量的同时极大地提高了重建速度，其重建时间是指数滤波和后投影方法的1/67~1/47。

Abstract

Photoacoustic imaging, which combines the advantages of optical imaging and acoustic imaging, is a non-invasive imaging technique with high spatial resolution and high contrast ratio. Therefore, photoacoustic imaging is one of the research hot spots in biomedical imaging. However, the reconstruction of photoacoustic image is a typically ill-posed inverse problem. Aiming at the ill-posedness of photoacoustic imaging and the slow velocity of the reconstruction owing to the large size of the system matrix, we present a fast exponential-filtering reconstruction method based on Lanczos double diagonalization and the effectiveness of the proposed method is proved by numerical simulation. The simulation results show that the proposed method can increase the reconstruction velocity greatly, while preserving the high quality of reconstructed images. The reconstruction time is 1/67-1/47 of those using exponential-filtering method and back projection method.

1 引言

光声成像是一种新的生物组织成像方法,具有高分辨率和高对比度^[1-3],有着广阔的应用前景^[4-6]。光声成像是一种基于生物组织内部光学吸收差异、以超声为媒质的无损生物光子成像方法,结合了纯光学成像对比度高和纯超声成像穿透深度大的优点,以超声探测器探测光声信号代替光学成像中的光子检测,从原理上避免了光学散射的影响,在保持生物样本自然状态的前提下能够得到高对比度和高分辨率的组织图像。由于不同生理状态的生物组织对光的吸收不同,光声图像同时反映了组织代谢的差异和病变特征。近年来,光声成像技术被广泛应用于肿瘤检测、血管成像、脑部结构和肾功能成像等领域^[7-10]。Xiang等^[7]利用光声成像技术开展了对乳腺癌和脑胶质瘤的研究,并进一步利用光声技术检测了肿瘤光动力治疗过程中血管的损伤情况。Yang等^[8]将光声成像技术应用于脑出血、脑缺氧和脑血管扩张等脑部疾病和功能活动的研究中。Li等^[9-10]采用基于纳米传感器的光声成像技术监控体内锂离子含量,实现了对针对性药物治疗作用的分析,并且将纳米传感器与治疗性药物结合使用,以实现对新生血管的成像和治疗。

光声成像重建技术就是利用采集的光声信号重建组织的光吸收分布图像。目前各国科研工作者已经发展了多种光声成像方法,如Kruger等^[11]提出了基于逆三维拉东(Radon)变换的反投影重建法。但是,该算法只在探测组织比扫描半径小很多的情况下才能得到较好的重建效果。之后,Kruger等^[12]提出了基于滤波反投影的算法,该算法利用圆周扫描的方式对组织进行信号采集,并对信号进行滤波处理,从而大大提高了成像质量。Xu等^[13-14]从光声信号产生的基础理论出发,利用傅里叶变换得到了不同测量模式下精确的时域重建算法和频域重建算法,但该方法形式繁琐,在实际应用中只能通过近似简化的表达式进行重建,因此存在一定的模型误差。针对该问题,该课题组又提出了一种统一的时域反投影方法^[15]。该方法形式简单,能够方便实现多种扫描形式下的精确成像。此外,科研工作者也提出了一些其他的光声成像重建方法,其中有代表性的方法包括基于傅里叶变换^[16]、有限元法^[17-18]和迭代法^[19]的重建方法。Wang等^[20]在傅里叶变换的基础上提出了一种不需要测量换能器脉冲响应的反卷积重建算法。为提高光声重建的精度,Tick等^[21-22]提出了基于贝叶斯的重建方法。需要注意的是,上述算法对于有限视角的欠采样信号或者非均匀信号不能得到满意的重建质量。为了解决该问题,杨迪武等^[23]在有限角度下采用代数重建算法提高了重建光声图像的分辨率和对比度。此外,Prakash等^[24]提出了一种基于基追踪反卷积的方法,提高了有限数据下光声成像的精度,并采用Lanczos双对角化进行降维,提高了该方法的重建速度。Bhatt等^[25]引入了指数滤波(EF)以重建光声图像,并在数据有限的情况下提高了图像质量,但是该方法需要对系统矩阵进行奇异值分解,系统矩阵较大时奇异值分解的计算时间长,计算过程非常耗时^[25]。

针对光声成像重建时间比较耗时的问题,提出了一种基于Lanczos双对角化的指数滤波重建方法。该方法先对系统矩阵进行降维以降低计算量,并在此基础上利用指数滤波方法进行图像重建。仿真结果表明,所提方法在保证重建图像高质量的前提下,与指数滤波方法相比,可以极大地缩短计算时间,将重建时间缩短了95%~98%。

2 方法

2.1 光声成像的数学模型

当短脉冲激光照射到待测生物组织上时,组织吸收激光能量产生热,若忽略热传导,此时光声成像的近似热方程为^[26]

\begin{matrix} ρ C_{p} \frac{\partial}{\partial t} T (r, t) = H (r, t), (1) \end{matrix}

式中r为三维空间位置矢量,t为时间,H(r,t)表示组织吸收的热量,T(r,t)表示组织升高的温度,ρ为组织的密度,C_p为比热。组织受热膨胀,形成一个初始声场,且声场随着时间的变化而变化。此时,得到光声成像声场的运动方程和扩散方程为

\begin{matrix} \begin{matrix} ρ \frac{\partial^{2}}{\partial t^{2}} u (r, t) = - \nabla p (r, t), (2) \\ \nabla u (r, t) = - \frac{p (r, t)}{ρ c^{2}} + βT (r, t), (3) \end{matrix} \end{matrix}

式中u(r,t)表示声位移,p(r,t)表示声压,c为声速,β为等压膨胀系数。联立(1)~(3)式,并进一步将H(r,t)表示为光吸收分布函数A(r)和入射激光关于时间的分布函数I(t)的乘积,可得

\begin{matrix} \nabla^{2} p (r, t) - \frac{1}{c^{2}} \frac{\partial^{2}}{\partial t^{2}} p (r, t) = - \frac{β}{C_{p}} A (r) \frac{\partial}{\partial t} I (t) 。 (4) \end{matrix}

(4)式是光声成像的基本方程,表示光声信号与组织的光吸收特性之间的关系,是描述光声成像的前向模型。

基于文献[ 25]所述的方法,对(1)~(3)式进行离散后得到的线性方程可表示为

\begin{matrix} AX = b, (5) \end{matrix}

式中A∈R^M×N²为光声场的系统矩阵(每个像素脉冲响应的集合),X∈R^N²^×¹为成像域中每个像素处的初始压力信号矩阵,b∈R^M×¹为光声信号矢量,其中M为探测器的时间步长与探测器个数的乘积,N²为光声图像的大小。为构造系统矩阵A,可以将在生物组织表面采集光声信号的过程看作一个时变的因果系统。系统矩阵A的每一列表示图像中每个像素点的系统响应。将X的第一个元素保留,其余元素置零,便可以得到第一个像素点的系统响应,对应系统矩阵A的第一列元素。为了提高计算效率,可利用光声信号的转换和衰减特性对X中剩余的N²-1个元素的系统响应进行求解。具体构造(5)式的方法可参照文献[ 25]。

2.2 基于Lanczos双对角化的指数滤波方法(Lanczos-EF)

光声成像的图像重建是利用探测器接收到的声压信号作为已知量来重建成像区域的声场分布。由于在实际测量过程中边界数据有限以及测量的光声信号中不可避免地混有噪声,因此不能对(5)式进行直接求解。为重建初始声压分布,最常用的方法是基于Tikhonov的重建方法^[27]。基于Tikhonov重建初始声压分布的方程可表示为

\begin{matrix} Ω = ‖ AX - b_{meas} ‖_{2}^{2} + λ ‖ X ‖_{2}^{2}, (6) \end{matrix}

式中A为(5)式得到的系统矩阵;b_meas为(5)式中与b相对应的边界采集的超声信号矢量;λ为正则化参数;‖·‖₂表示L₂范数。

对A进行奇异值分解,则(6)式的解可表示为

\begin{matrix} \begin{matrix} X_{tik h} = {(V \sum^{T} \sum V^{T} + λI)}^{- 1} V \sum^{T} U^{T} b_{meas} = \\ V \sum_{F} U b_{meas}, (7) \end{matrix} \end{matrix}

式中U=(u₁,u₂,…,u_n)、V=(v₁,v₂,…,v_n)为左右奇异正交矩阵,且满足U^TU=V^TV=I_n和∑=diag(σ₁,σ₂,…,σ_n),其中σ₁≥σ₂≥…≥σ_n>0为A的奇异值;∑_F=diag(ϕ_i/σ_i),其中ϕ_i= $\begin{matrix} σ_{i}^{2} \end{matrix}$ /( $\begin{matrix} σ_{i}^{2} \end{matrix}$ +λ)为正则化滤波因子。该方法的实质是在最小二乘解的基础上增加了约束高频分量的滤波因子。但从 (7)式中可以看出,正则化参数对重建的光声图像有着重要的影响。正则化参数值越大,则目标函数中赋予解的范数的权重就越大,从而可以保证所求解的范数较小,但此时以牺牲数据的拟合程度为代价。反之,正则化参数值越小,则(7)式越接近于未正则化的原问题,此时正则化后的问题可能还存在某种程度上的不适应性。

指数滤波因子与Tikhonov正则化奇异值分解的方法类似^[27],不同的是该方法的滤波因子ϕ_i=1-exp(- $\begin{matrix} σ_{i}^{2} \end{matrix}$ /λ)。Tikhonov正则化所用的滤波因子ϕ_i= $\begin{matrix} σ_{i}^{2} \end{matrix}$ /( $\begin{matrix} σ_{i}^{2} \end{matrix}$ +λ),指数滤波正则化所用的滤波因子ϕ_i=1-exp(- $\begin{matrix} σ_{i}^{2} \end{matrix}$ /λ)。

图1为单源仿真中正则化参数变化时两种滤波因子的曲线。图中所用的奇异值是通过对单源实验中系统矩阵A进行奇异值分解后得到的。通过比较这两种滤波因子可以发现,当且仅当σ²≪λ时,Tikhonov正则化才可以看作是指数滤波的近似。

图 1. 单源仿真中正则化参数变化时两种滤波因子的曲线

Fig. 1. Curves of two filter factors versus regularization parameters in single source simulation

下载图片查看所有图片

对系统矩阵A进行奇异值分解时,由于光声成像系统矩阵较大,因此计算比较耗时。为此,可以先对系统矩阵进行Lanczos双对角化处理,以对系统矩阵进行降维。对于M×N²维下的矩阵A,经过Lanczos双对角化处理后可以得到(k+1)×k维的双对角矩阵B和两个标准正交矩阵U_k+₁和V_k,其中B_k可表示为

\begin{matrix} B_{k} = [\begin{matrix} α_{1} \\ β_{2} & α_{2} \\ β_{3} & ⋱ \\ ⋱ & α_{k} \\ β_{k + 1} \end{matrix}], (8) \end{matrix}

式中α_i(i=1,2,…,k)为主对角阵元素,β_j(j=2,3,…,k+1)为辅对角阵元素。U_k+₁=(u₁,u₂,…,u_k+₁), $\begin{matrix} {U^{T}}_{k + 1} \end{matrix}$ U_k+₁=I;V_k=(v₁,v₂,…,v_k), $\begin{matrix} {V^{T}}_{k} \end{matrix}$ V_k=I, 且满足U^TAV=B_k。在重建过程中,输入A_M×N₂,b_M以及迭代次数0<k<min(M,N²);输出为U_M×k+₁,B_k+₁_×k和 $\begin{matrix} V_{k \times N^{2}} \end{matrix}$ 。Lanczos 双对角化处理的算法步骤为:

设置初始值β=norm(b),u=b/β,u=A^Tu,α=norm(v)

v=v/α,B(1,1)=α,U=[u],V=[v];

for j=1,2,…,k

u=A×v-α×u,β=norm(u),u=u/β,

B(j,j-1)=β,U=[U,u],V=A^T×u-β×v,

α=norm(v),v=v/α,B(j∶j)=α,V=[V,v]

end

经过k步迭代后,原始问题的残差可表示为

\begin{matrix} r_{k} = AX - b = A V_{k} y_{k} - β_{1} u_{1} = U_{k + 1} (B_{k} y_{k} - β_{1} e_{1}), (9) \end{matrix}

式中y_k为e,e₁为单位向量。U_k+₁为标准正交矩阵,因此原问题等价于求解以下最小化问题:

\begin{matrix} mi \underset{yk}{n} ‖ B_{k} y_{k} - β_{1} e_{1} ‖_{2}^{2} 。 (10) \end{matrix}

只要求出 $\begin{matrix} y_{k}^{*} \end{matrix}$ ,即可得到原问题的解为

\begin{matrix} X^{*} = V_{k} y^{*}, (11) \end{matrix}

式中y^*为y_k的最优解。

所提算法流程图如图2所示。在该算法中,迭代次数k影响着重建结果。经过大量数值仿真,发现当k=25时所提算法在大多数情况下可以获得较好的重建图像质量,同时具有较快的重建速度。因此,除非特别说明,下文中迭代次数k的值设为25。

图 2. 所提算法流程图

Fig. 2. Flow chart of the proposed algorithm

下载图片查看所有图片

2.3 品质因数

2.3.1 皮尔逊相关系数

皮尔逊相关系数(PC)^[24]可用于衡量重建图像与目标图像之间的相关程度,其表达式为

\begin{matrix} C_{PC} (X, X_{recon}) = \frac{cov (X, X_{recon})}{ρ (X) ρ (X_{recon})}, (12) \end{matrix}

式中X为预期的初始压力分布,X_recon为重建的初始压力分布,cov(·)表示协方差,ρ(·)为标准偏差。PC描述的是目标的准确可检测性(即空间保真度),PC值越高则说明重建结果越好。

2.3.2 对比噪声比

对比噪声比(CNR)^[24]是基于对比度的图像质量的评价参数,其定义为

\begin{matrix} R_{CNR} = \frac{μ_{roi} - μ_{back}}{\sqrt[]{ρ_{roi}^{2} a_{roi} + ρ_{back}^{2} a_{back}}}, (13) \end{matrix}

式中μ和ρ分别表示重建初始声压的平均值和标准偏差,下标roi和back分别表示目标区域和背景区域;a_roi和a_back分别表示目标区域和背景区域的面积。CNR值越高,则说明目标区域与背景区域的差异性越大,即图像重建性能越好。

图 3. 光声数据采集示意图

Fig. 3. Diagram of photoacoustic data acquisition

下载图片查看所有图片

3 仿真结果

为了验证方法的有效性,开展了多次数值仿真。在仿真中,假定介质具有均匀的超声特性,并且在传播过程中没有声音的吸收和散射,则将声速设为1500 m·s^-1 (软组织的近似值)。仿真中共设置了40个超声波检测器,其中心频率为2.25 MHz。将这些检测器等距离放置在一个半径为22 mm的圆上,成像域的大小设置为101 pixel×101 pixel(每个像素大小为0.1 mm),光声数据采集示意图如图3所示。在仿真中,光声信号采集的时间间隔是50 ns,采集时间为25 μs。为了证明所提算法(Lanczos-EF)的有效性,将其与目前常用的基于后投影(BP)^[25]、Tikhonov正则化和指数滤波(EF)^[25]的重建方法进行比较。考虑到正则化参数对基于Tikhonov和指数滤波重建方法的结果有着较大的影响,在[0,1]范围内尝试了一组正则化参数,并选取了最优的重建结果。所有仿真都是在处理器为Intel (R)Core(TM)@2.70 GHz、内存4G的个人计算机上进行的。采用的仿真软件为Matlab 2014b,并使用k-Wave工具箱(version 1.2)仿真光声信号的产生、传播和衰减。

这里共开展了3次数值仿真。其中,第一次和第二次仿真分别是对单源和双源的重建,第三次仿真是对血管结构的重建。在仿真中,假设初始压力分布为1 kPa,同时在仿真的超声信号中加入了1%的高斯随机噪声,模拟信噪比为40 dB的信号中的噪声。

图4为单源时不同算法的重建结果。图4(b)~(d)分别为基于BP、Tikhonov和指数滤波方法的重建结果。图4(e)为利用所提方法得到的重建结果。单源时不同方法的量化结果比较如图5和表1所示,表1中λ表示Tikhonov正则化参数,k表示Lanczos双对角化处理中的迭代次数。从图5中可以看出,基于BP和Tikhonov方法得到的PC值(0.69和0.71)和CNR值(22.8和21.2)的差别不大。与BP和Tikhonov方法相比,指数滤波方法可以获得较好的PC值(0.79)和CNR值(28.3)。但从图5中可以看到,这3种方法的计算比较耗时。BP方法相对Tikhonov和指数滤波方法而言,重建速度较快,需要518 s。而利用所提方法则可以明显地提高重建速度,其重建时间与指数滤波方法相比缩短了约98%,降至约23 s。与BP方法相比,所提方法的重建时间也缩短了约95.5%。此外,利用所提方法得到的重建图像质量与基于指数滤波方法得到的结果相比,PC值从0.79提高到0.80,但CNR值没有变化,仍为28.3。这是因为指数滤波方法(可看作一个低通滤波器)能够有效地去除重建光声图像中的高频噪声,同时Lanczos双对角化处理有效地降低了系统矩阵的维数,从而提高了重建速度。因此,所提方法在保证重建图像高质量的同时,可以明显提高重建速度。

图 4. 真实单源分布和单源时不同算法的重建结果。(a)真实单源分布;基于(b) BP、(c) Tikhonov、(d) EF和(e) Lanczos-EF方法的重建结果

Fig. 4. Distribution of true single source and reconstruction results using different algorithms in case of signal source.(a) Distribution of true single source; reconstruction results obtained by (b) BP, (c) Tikhonov,(d) EF, and (e) Lanczos-EF methods

下载图片查看所有图片

图 5. 单源时不同方法的量化结果比较。(a) PC;(b) CNR;(c)重建时间

Fig. 5. Comparison of quantitative results using different reconstruction methods in case of single source.(a) PC; (b) CNR; (c) reconstruction time

下载图片查看所有图片

图 6. 真实双源分布和双源时不同算法的重建结果。(a)真实双源分布;基于(b) BP、(c) Tikhonov、(d) EF和(e) Lanczos-EF方法的重建结果

Fig. 6. Distribution of true two sources and reconstruction results using different algorithms in case of two sources.(a) Distribution of true two sources; reconstruction results obtained by (b) BP, (c) Tikhonov,(d) EF and (e) Lanczos-EF methods

下载图片查看所有图片

图 7. 双源时不同方法的量化结果比较。(a) PC;(b) CNR;(c)重建时间

Fig. 7. Comparison of quantitative results using different reconstruction methods in case of two sources.(a) PC; (b) CNR; (c) reconstruction time

下载图片查看所有图片

图6为双源时不同算法的重建结果,双源时不同方法的量化结果如图7和表1所示。从图6中可以看到,基于BP方法得到的重建结果[图6(b)]最差。基于Tikhonov方法重建的图像[图6(c)]的PC值最高,约为0.83,CNR值约为13.8。利用指数滤波方法和所提方法也得到了较好的重建结果,两种方法得到的PC值和CNR值都相同,分别为0.82和14.2。与单源情况类似,基于BP、Tikhonov和指数滤波方法都需要较长的重建时间,其中耗时最短的BP方法也需要512 s,其他两种方法的计算时间都超过了1470 s。而所提方法在没有明显降低重建图像质量的前提下,明显提高了重建速度,将重建时间缩短至27 s。

光声成像技术可被广泛应用于可视化内部血管结构,因此为进一步验证算法的有效性,开展了基于血管的重建仿真。图8(a)为血管结构分布,基于不同方法重建的血管图像如图8(b)~(e)所示。与上述仿真类似,基于BP方法得到的重建结果[图8(b)]最差,很难看清血管结构。基于Tikhonov方法得到了较好的重建结果,但是重建的初始压力与真实值有较大的差别。而基于指数滤波的重建方法则获得了最好的重建结果,不仅重建的初始压力与真实值更接近,而且重建的PC值(0.65)和CNR值(4.2)最大。所提方法获得了与指数滤波方法类似的重建图像质量,但明显缩短了重建时间,将重建时间由指数滤波方法需要的1342 s降低到20 s。

从以上结果中可以看出,基于BP和Tikhonov正则化方法重建的光声图像不仅图像质量低(PC值和CNR值较小),而且所需要的计算时间较长。基于指数滤波方法则可以获得较好的光声图像质量(PC值和CNR值较大),但基于指数滤波方法的缺点是重建速度慢(重建时间大于1100 s),限制了该

图 8. 真实血管图像分布和血管仿真结果。(a)真实血管图像分布;基于(b) BP、(c) Tikhonov、(d) EF和(e) Lanczos-EF方法的重建结果

Fig. 8. Distribution of true blood vessel image and simulation results of blood vessel. (a) Distribution of true blood vessel image; reconstruction results obtained by (b) BP, (c) Tikhonov, (d) EF and (e) Lanczos-EF methods

下载图片查看所有图片

图 9. 不同方法重建的血管图像的量化结果比较。(a) PC;(b) CNR;(c)重建时间

Fig. 9. Comparison of quantitative results of blood vessel image obtained by different methods.(a) PC; (b) CNR; (c) reconstruction time

下载图片查看所有图片

表 1. 4种重建方法的量化结果比较

Table 1. Comparison of quantitative results of four reconstruction methods

Method	Single source				Two sources				Blood vessel
Method	λ/k	PC	CNR	Time /s	λ/k	PC	CNR	Time /s	λ/k	PC	CNR	Time /s
BP	-	0.69	22.8	518	-	0.65	4.8	512	-	0.54	3.1	541
Tikhonov	6.2×10^-3	0.71	21.2	1058	5.6×10^-3	0.83	13.8	1480	4.3×10^-3	0.62	4.0	1057
EF	1.2×10^-3	0.79	28.3	1102	1.0×10^-3	0.82	14.2	1470	1.6×10^-3	0.65	4.2	1341
Lanczos-EF	25	0.80	28.3	23	25	0.82	14.2	27	25	0.65	4.2	20

查看所有表

方法在实际中的应用。因此,在利用指数滤波方法的基础上基于Lanczos双对角化对系统矩阵进行降维,在保证高重建图像质量的同时降低了重建速度。仿真结果表明,所提方法不仅可以获得较好的重建图像质量,而且能够极大地提高计算效率。此外,结果表明基于指数滤波方法可以获得较好的重建图像质量。这是由于:随着奇异值数量的增加,指数滤波方法采用递减权重对奇异值进行过滤,因此可以将其看作一个低通滤波器,可有效去除重建光声图像中的高频噪声。但需要注意的是,Lanczos双对角化中的迭代次数k和指数滤波因子中的正则化参数λ都会对重建结果产生较大的影响。虽然实验中将迭代次数k设为25,并取得了较好的结果,但如何自适应地选择k仍需要进一步研究。考虑到正则化参数λ对重建结果的影响,在[0,1]范围内尝试了大量正则化参数,并选取了重建质量最好的光声重建图像作为最终的重建结果,这不仅会耗费大量的计算时间,而且不够客观。因此,如何准确地对指数滤波方法中的正则化参数进行自适应选择也需要进一步考虑。

4 结论

光声成像的应用潜力依赖于重建算法研究的新进展。针对光声成像重建算法展开研究,提出了一种重建图像的高质量、快速的光声成像重建方法,并对该方法的有效性进行了验证。结果表明,基于Lanczos双对角化的指数滤波方法在保证重建图像质量较高的同时,可以大大缩短重建时间;所提方法的重建时间是指数滤波和BP方法的1/67~1/47。但基于Lanczos双对角化的指数滤波方法对于成像区域较简单情况的重建效果较好,但是对于目标区域复杂情况的重建效果较差。这是由于指数滤波方法相当于一个低通滤波器,它在有效去除高频噪声的同时也去除了重建图像中的部分细节(即图像中的高频分量),从而导致重建的图像边缘模糊。今后将研究并改进基于Lanczos双对角化的指数滤波,以提高复杂目标区域的成像效果。

参考文献

[1] Xiang LZ, Gu HM. Photoacoustic imaging with focused ultrasonic probe beam[C]. SPIE, 2006, 6047: 60470H.

[2] Xiang LZ, XingD, Gu HM, et al. Pulse laser integrated photodynamic therapy and photoacoustic imaging[C]. SPIE, 2007, 6437: 64372B.

[3] 曾光, 石岩, 宋亮, 等. 新型声学分辨率光声显微镜系统照明设计[J]. 中国激光, 2016, 43(2): 0204002.

Zeng G, Shi Y, Song L, et al. Novel illumination design of acoustic resolution photoacoustic microscopy system[J]. Chinese Journal of Lasers, 2016, 43(2): 0204002.

[4] Wang X D, Pang Y J, Ku G, et al. Noninvasive laser-induced photoacoustic tomography for structural and functional in vivo imaging of the brain[J]. Nature Biotechnology, 2003, 21(7): 803-806.

[5] Kruger R A, Kiser W L, Reinecke D R, et al. Thermoacoustic computed tomography using a conventional linear transducer array[J]. Medical Physics, 2003, 30(5): 856-860.

[6] 任重, 刘国栋, 黄振. 基于可调谐脉冲激发的血糖浓度光声无损检测研究[J]. 中国激光, 2016, 43(2): 0204001.

Ren Z, Liu G D, Huang Z. Study on photoacoustic noninvasive detection for blood glucose concentration based on tunable pulsed[J]. Chinese Journal of Lasers, 2016, 43(2): 0204001.

[7] Xiang L Z, Xing D, Gu H M, et al. In vivo monitoring of neovascularization in tumour angiogenesis by photoacoustic tomography[J]. Chinese Physics Letters, 2007, 24(3): 751-754.

[8] Yang S H, Xing D, Lao Y, et al. Noninvasive monitoring of traumatic brain injury and post-traumatic rehabilitation with laser-induced photoacoustic imaging[J]. Applied Physics Letters, 2007, 90(24): 243902.

[9] Li L, Jun X. Terence T W W, et al. In vivo deep brain imaging of rats using oral-cavity illuminated photoacoustic computed tomography[J]. Journal of Biomedical Optics, 2015, 20(1): 016019.

[10] Yao J J, Jun X, Wang L V. Multiscale functional and molecular photoacoustic tomography[J]. Ultrasonic Imaging, 2016, 38(1): 44-62.

[11] Kruger R A, Liu P, Fang Y R, et al. Photoacoustic ultrasound (PAUS) - Reconstruction tomography[J]. Medical Physics, 1995, 22(10): 1605-1609.

[12] Kruger R A, Reinecke D R, Kruger G A. Thermoacoustic computed tomography-Technical considerations[J]. Medical Physics, 1999, 26(9): 1832-1837.

[13] Xu Y, Xu M, Wang L V. Exact frequency-domain reconstruction for thermoacoustic tomography. II. Cylindrical geometry[J]. IEEE Transactions on Medical Imaging, 2002, 21(7): 829-833.

[14] Xu M H, Wang L V. Time-domain reconstruction for thermoacoustic tomography in a spherical geometry[J]. IEEE Transactions on Medical Imaging, 2002, 21(7): 814-822.

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

[16] Köstli K P, Beard P C. Two-dimensional photoacoustic imaging by use of Fourier-transform image reconstruction and a detector with an anisotropic response[J]. Applied Optics, 2003, 42(10): 1899-1908.

[17] Jiang H B, Yuan Z, Gu X J. Spatially varying optical and acoustic property reconstruction using finite-element-based photoacoustic tomography[J]. Journal of the Optical Society of America A, 2006, 23(4): 878-888.

[18] Yuan Z, Jiang H B. Three-dimensional finite element-based photoacoustic tomography: Initial results[J]. SPIE, 2007, 6437: 64371A.

[19] Paltauf G, Viator J A, Prahl S A, et al. Iterative reconstruction algorithm for optoacoustic imaging[J]. Journal of the Acoustical Society of America, 2002, 112(4): 1536-1544.

[20] Wang Y, Xing D, Zeng Y G, et al. Photoacoustic imaging with deconvolution algorithm[J]. Physics in Medicine & Biology, 2004, 49(14): 3117-3124.

[21] Tick J, Pulkkinen A, Tarvainen T. Image reconstruction with uncertainty quantification in photoacoustic tomography[J]. Journal of the Acoustical Society of America, 2016, 139(4): 1951-1961.

[22] Pulkkinen A, Cox B T, Arridge S R, et al. Direct estimation of optical parameters from photoacoustic time series in quantitative photoacoustic tomography[J]. IEEE Transactions on Medical Imaging, 2016, 35(11): 2497-2508.

[23] 杨迪武, 邢达, 王毅, 等. 基于代数重建算法的有限角度扫描的光声成像[J]. 光学学报, 2005, 25(6): 772-776.

Yang D W, Xing D, Wang Y, et al. Limited-view scanning photoacoustic imaging based on algebraic reconstruction techniques[J]. Acta Optica Sinica, 2005, 25(6): 772-776.

[24] Prakash J, Raju A S, Shaw C B, et al. Basis pursuit deconvolution for improving model-based reconstructed images in photoacoustic tomography[J]. Biomedical Optics Express, 2014, 5(5): 1363-1377.

[25] Bhatt M, Gutta S, Yalavarthy P K. Exponential filtering of singular values improves photoacoustic image reconstruction[J]. Journal of the Optical Society of America A, 2016, 33(9): 1785-1792.

[26] Wang L V. Tutorial on photoacoustic microscopy and computed tomography[J]. IEEE Journal of Selected Topics in Quantum Electronics, 2008, 14(1): 171-179.

[27] Martin D R, Reichel L. Projected tikhonov regularization of Large-Scale discrete Ill-Posed problems[J]. Journal of Scientific Computing, 2013, 56(3): 471-493.

1 引言

2 方法

2.1 光声成像的数学模型

2.2 基于Lanczos双对角化的指数滤波方法(Lanczos-EF)

2.3 品质因数

3 仿真结果

4 结论

张娜, 冯金超, 李哲, 贾克斌. 基于Lanczos双对角化的快速光声成像重建方法[J]. 中国激光, 2018, 45(3): 0307018. Zhang Na, Feng Jinchao, Li Zhe, Jia Kebin. Fast Photoacoustic Imaging Reconstruction Method Based on Lanczos Double Diagonalization[J]. Chinese Journal of Lasers, 2018, 45(3): 0307018.

基于Lanczos双对角化的快速光声成像重建方法下载： 907次

1 引言

2 方法

2.1 光声成像的数学模型

2.2 基于Lanczos双对角化的指数滤波方法(Lanczos-EF)

图 1. 单源仿真中正则化参数变化时两种滤波因子的曲线

Fig. 1. Curves of two filter factors versus regularization parameters in single source simulation

图 2. 所提算法流程图

Fig. 2. Flow chart of the proposed algorithm

2.3 品质因数

图 3. 光声数据采集示意图

Fig. 3. Diagram of photoacoustic data acquisition

3 仿真结果

图 4. 真实单源分布和单源时不同算法的重建结果。(a)真实单源分布;基于(b) BP、(c) Tikhonov、(d) EF和(e) Lanczos-EF方法的重建结果

Fig. 4. Distribution of true single source and reconstruction results using different algorithms in case of signal source.(a) Distribution of true single source; reconstruction results obtained by (b) BP, (c) Tikhonov,(d) EF, and (e) Lanczos-EF methods

图 5. 单源时不同方法的量化结果比较。(a) PC;(b) CNR;(c)重建时间

Fig. 5. Comparison of quantitative results using different reconstruction methods in case of single source.(a) PC; (b) CNR; (c) reconstruction time

图 6. 真实双源分布和双源时不同算法的重建结果。(a)真实双源分布;基于(b) BP、(c) Tikhonov、(d) EF和(e) Lanczos-EF方法的重建结果

Fig. 6. Distribution of true two sources and reconstruction results using different algorithms in case of two sources.(a) Distribution of true two sources; reconstruction results obtained by (b) BP, (c) Tikhonov,(d) EF and (e) Lanczos-EF methods

图 7. 双源时不同方法的量化结果比较。(a) PC;(b) CNR;(c)重建时间

Fig. 7. Comparison of quantitative results using different reconstruction methods in case of two sources.(a) PC; (b) CNR; (c) reconstruction time

图 8. 真实血管图像分布和血管仿真结果。(a)真实血管图像分布;基于(b) BP、(c) Tikhonov、(d) EF和(e) Lanczos-EF方法的重建结果

Fig. 8. Distribution of true blood vessel image and simulation results of blood vessel. (a) Distribution of true blood vessel image; reconstruction results obtained by (b) BP, (c) Tikhonov, (d) EF and (e) Lanczos-EF methods

图 9. 不同方法重建的血管图像的量化结果比较。(a) PC;(b) CNR;(c)重建时间

Fig. 9. Comparison of quantitative results of blood vessel image obtained by different methods.(a) PC; (b) CNR; (c) reconstruction time

表 1. 4种重建方法的量化结果比较

Table 1. Comparison of quantitative results of four reconstruction methods

4 结论

Article Outline

关于本站 Cookie 的使用提示

全站搜索

基于Lanczos双对角化的快速光声成像重建方法 下载： 907次

1 引言

2 方法

2.1 光声成像的数学模型

2.2 基于Lanczos双对角化的指数滤波方法(Lanczos-EF)

图 1. 单源仿真中正则化参数变化时两种滤波因子的曲线

Fig. 1. Curves of two filter factors versus regularization parameters in single source simulation

图 2. 所提算法流程图

Fig. 2. Flow chart of the proposed algorithm

2.3 品质因数

图 3. 光声数据采集示意图

Fig. 3. Diagram of photoacoustic data acquisition

3 仿真结果

图 4. 真实单源分布和单源时不同算法的重建结果。(a)真实单源分布;基于(b) BP、(c) Tikhonov、(d) EF和(e) Lanczos-EF方法的重建结果

Fig. 4. Distribution of true single source and reconstruction results using different algorithms in case of signal source.(a) Distribution of true single source; reconstruction results obtained by (b) BP, (c) Tikhonov,(d) EF, and (e) Lanczos-EF methods

图 5. 单源时不同方法的量化结果比较。(a) PC;(b) CNR;(c)重建时间

Fig. 5. Comparison of quantitative results using different reconstruction methods in case of single source.(a) PC; (b) CNR; (c) reconstruction time

图 6. 真实双源分布和双源时不同算法的重建结果。(a)真实双源分布;基于(b) BP、(c) Tikhonov、(d) EF和(e) Lanczos-EF方法的重建结果

Fig. 6. Distribution of true two sources and reconstruction results using different algorithms in case of two sources.(a) Distribution of true two sources; reconstruction results obtained by (b) BP, (c) Tikhonov,(d) EF and (e) Lanczos-EF methods

图 7. 双源时不同方法的量化结果比较。(a) PC;(b) CNR;(c)重建时间

Fig. 7. Comparison of quantitative results using different reconstruction methods in case of two sources.(a) PC; (b) CNR; (c) reconstruction time

图 8. 真实血管图像分布和血管仿真结果。(a)真实血管图像分布;基于(b) BP、(c) Tikhonov、(d) EF和(e) Lanczos-EF方法的重建结果

Fig. 8. Distribution of true blood vessel image and simulation results of blood vessel. (a) Distribution of true blood vessel image; reconstruction results obtained by (b) BP, (c) Tikhonov, (d) EF and (e) Lanczos-EF methods

图 9. 不同方法重建的血管图像的量化结果比较。(a) PC;(b) CNR;(c)重建时间

Fig. 9. Comparison of quantitative results of blood vessel image obtained by different methods.(a) PC; (b) CNR; (c) reconstruction time

表 1. 4种重建方法的量化结果比较

Table 1. Comparison of quantitative results of four reconstruction methods

4 结论

Article Outline

相关论文

相关资讯

关于本站 Cookie 的使用提示

全站搜索

基于Lanczos双对角化的快速光声成像重建方法下载： 907次