CN104586363B - 基于图像块稀疏系数的快速光声成像图像重建方法 - Google Patents

Abstract

本发明属于光声成像技术领域，具体为一种基于图像块稀疏系数的快速光声成像的图像重建方法。本发明首先通过离散余弦变换优化采样矩阵，随后计算图像的块稀疏系数，对图像进行迭代修正，随后再根据迭代过的图像重新计算上述参数，进行反复迭代直至达到结束条件，最后输出重建的光声图像。本发明能有效地降低迭代算法的复杂度，同时还具有较好的重建质量，对于光声成像的图像重建具有实际意义。

Description

基于图像块稀疏系数的快速光声成像图像重建方法

技术领域

本发明属于光声成像技术领域，具体涉及一种快速光声成像的图像重建方法。

背景技术

光声成像是一种新型医学成像方法，其理论基础为光声效应。光声成像结合了光学成像和超声成像的优点，具有高对比度和高成像深度的特性。由于其非电离波的特性，不会在检测的过程中对人体产生伤害^[1,2]，在各个应用领域都有着非常大的潜力，现在主要应用领域有乳腺肿瘤检测^[3]、血管成像^[4]和脑损伤探测^[5]等。同时由于所成图像是人体组织的光吸收特性，因而一定程度上反映了成像组织与光吸收特性相关的病理学特性^[2]。利用这一特性可以将光声成像应用于功能成像^[5,6]，对于癌症的早期检测也有着非常好的效果^[1]。

在光声成像的具体应用中，我们使用一个持续时间非常短的激光脉冲照射成像组织。成像组织会吸收一部分光能将它转化为热能，使成像组织发生热弹性膨胀，从而发出超声波。这一过程被称为光声效应。而光声成像就是使用超声换能器检测成像组织发出的超声波，通过超声换能器在不同位置扫描可以采集到光声信号，随后使用图像重建算法就可以计算出组织的光吸收分布。目前针对圆周扫描有逆Radon变换重建方法^[7]、滤波反投影法^[8]、时域重建法^[9]和反卷积重建法^[10]等方法；针对直线扫描有DAS法^[11]和二维重建法^[11]等方法。上述方法都属于解析算法，无需迭代直接得到结果，在采样点较稀疏时重建图像的精度不高，成像质量差，并且受限于特定的扫描方式。现在主流光声成像图像重建算法主要是基于迭代优化的重建方法^[12-15]，这些方法的优势在于成像精度高、不受扫描方式的影响，但缺点是成像速度较慢，算法运算量大。

针对上述图像重建方法存在的问题，本发明中的图像重建方法使用了图像块稀疏系数，同时使用离散余弦变换，降低了迭代算法的运算量和计算时间，提升了重建算法的效率，重建图像质量并没有受到影响。

发明内容

本发明的目的在于提供一种成像精度高、成像速度快、运算量较小的光声成像的图像重建方法。

本发明提出适用于光声成像的图像重建方法，是通过计算图像块稀疏系数，对重建图像进行修正并迭代，结合Barzilai-Borwein梯度下降法，获得最终的重建光声图像。

在光声成像中，用激光短脉冲垂直于待成像平面照射生物组织，超声换能器在待成像平面内进行扫描。通常激光脉冲持续时间远小于组织的热扩散时间，根据光声效应和超声的运动方程和扩散方程，可以得到光声成像的基本方程^[2]：

(1)

其中是位置处的声压，是成像组织的光吸收分布图，t是时间，I(t)是激光脉冲能量函数，c是生物组织中的声速，和分别是生物组织的等压膨胀系数和比热容。光声成像的图像重建技术，就是通过求出。

使用格林函数对方程(1)进行求解^[8]。对于某采样点，有：

(2)

将式(2)进行变形，可得：

(3)

记采样点处实际采样得到的光声信号的积分与采样时间的乘积为：

(4)

实际应用时，可将图像和采样信号积分分别离散化，并记成矢量形式。若重建图像的大小为(X，Y分别为图像的行数和列数)，则重建图像的总像素为N (N=XY)，即矢量化后的图像可记为长度为N的列矢量u。若采样点个数为Q，每个采样点的信号长度为M，可将(3)式写成：

(5)

其中是第i个采样点的光声信号积分与采样时间的乘积矢量；是第i个采样点的采样矩阵，其计算步骤为：

(a) 先计算大小为的矩阵：

(6)

其中，是中的序号，是采样点的坐标，dx是图像相邻像素间的实际距离，dt是离散时间步长；

(b)将矩阵矢量化得到一个N维的列矢量，作为采样矩阵的第j个列矢量。

(c)计算M次(j=1～M)后得到第i点的采样矩阵；

(d)重复步骤(a)~(c)得到Q个采样矩阵 (i=1～Q)。将它们联立起来，得到总的采样矩阵A：

(7)

于是，式(5)可以归纳为：

(8)

其中f、A和u的大小分别为、和。

综上，基于迭代的光声重建方法就是通过总的采样信号积分与采样时间的乘积矢量f和按步骤(a)~(d)得到的采样矩阵A，基于迭代的方法求出重建的光声图像u。

本发明提出了一种实现光声图像重建的有效迭代方法，其具体迭代步骤为：

(1) 输入原始重建图像，设置各参数的初始值；

(2) 设置阈值TH, 对数据和采样矩阵进行筛选；

(3) 根据上一次迭代得到的图像，按照基于块稀疏系数迭代算法计算新的重建图像；

(4) 更新迭代中使用的参数；

(5) 判断是否达到迭代结束条件，若未达到则返回步骤(3)；若达到就结束迭代，得到重建图像。

步骤(1)所述 输入原始重建图像，设置各参数的初始值：设定原始重建图像为，分块数为32块，各参数的初始值设为：残差系数μ=0.8，步长参数r⁰=[r ₁ ⁰, r ₂ ⁰,…, r ₃₂ ⁰]=[0.05，0.05，…,0.05]，迭代终止阈值；

步骤(2) 所述设置阈值TH, 对数据和采样矩阵进行筛选：对于采样数据和采样矩阵进行离散余弦变换：

(9)

其中f’是实际采样得到的光声信号的积分与采样时间的乘积矢量经过离散余弦变换之后得到的值，f是实际采样得到的光声信号的积分与采样时间的乘积矢量，D是离散余弦变换矩阵，A’是光声信号的采样矩阵经过离散余弦变换之后得到的矩阵，A是光声信号的采样矩阵。取出f’中大于TH的值组成新矢量b，从A’中取出与之相对应的行，组成新的采样矩阵W。

步骤(3)所述根据上一次迭代得到的图像，按照基于块稀疏系数迭代算法计算新的重建图像：迭代公式为：

(10)

其中和表示第k次和第k+1次迭代得到以矢量形式表示的重建图像的第a块，r _a ^k 表示第a块图像第k次迭代的步长参数，W[a]表示第a块图像所对应的采样矩阵，u ^k 为第k次迭代得到的以矢量形式表示的重建图像，μ表示迭代残差参数，T表示为矩阵的转置，y表示为了表达简便使用的中间参数。

步骤(4)所述更新迭代中使用的参数，其计算方法为：

(11)

其中，r _a ^k+1表示第a块图像第k+1次迭代的步长参数，u ^k+1为第k+1次迭代得到的以矢量形式表示的重建图像。

步骤(5) 中迭代步数k变为k+1，判断是否达到迭代结束条件，若未达到则返回步骤(3)；若达到就结束迭代，得到重建图像。

具体的判断方式为：

(12)

使用本发明中的方法进行光声图像重建的具体流程图如图1所示。

与现有技术相比，本发明引入了块稀疏系数和离散余弦变换，在图像域和信号域都实现了降维，同时使用了Barzilai-Borwein梯度下降法，算法收敛速度快，降低了运算量，在保证算法质量的基础上提升了迭代算法的运算效率，具有一定的现实意义。

附图说明

图1为本发明基于块稀疏系数的光声成像图像重建方法具体流程图。

图2为待成像组织的光吸收分布图。

图3为对组织进行圆周扫描，采样点数为30个的情况下各种光声图像重建方法的结果比较。其中，(a)为滤波反投影法，(b)为全变分参数梯度下降法，(c)为离散余弦变换法，(d)为本发明方法。其重建结果的峰值信噪比(PSNR)分别为14.34 dB、32.17 dB、23.27dB、28.92 dB。

图4为各扫描方式下的图像重建结果。其中，(a)有限角度圆周扫描，(b)直线扫描。

图5为在离体实验仿体和重建结果。其中，(a)离体组织仿体的照片，(b)显示了用60个探测信号重建的图像。

具体实施方式

对本发明提出的光声图像重建方法在计算机上进行仿真。测试本发明光声图像重建方法的有效性，以及相对于其他方法的优越性。

1、确定组织的光吸收分布图，如图2所示，组织大小为89.6 mm89.6 mm，重建图像大小为128128像素。根据(2)式采集光声信号，圆周扫描半径为42 mm，角度步长为12°，共30个采样点。比较本发明方法和滤波反投影法、全变分参数梯度下降法^[15]和离散余弦变换法^[16]的图像重建结果，结果在图3中给出。选取峰值信噪比PSNR(Peak Signal to NoiseRatio)为量化指标，单位为dB。PSNR值越大，图像的重建效果越好。

(13)

其中为原始图像。

仿真结果表明，在大幅度降低了计算量的情况下，本发明方法的PSNR值高于解析方法（滤波反投影法）和其他高速迭代算法（离散余弦变换法），与全变分参数梯度下降法的结果非常接近。说明本发明在大幅度降低计算量的同时，仍然可以保证非常优秀的图像成像质量。

2、以上述仿真条件，分别验证有限角度扫描和直线扫描方式下算法的有效性。有限角度扫描半径为42 mm，扫描角度为120°，共15个采样点，采样间隔为8°；直线扫描间隔为2.8 mm，共30个采样点。使用本发明中的方法对采集到的光声信号进行重建，得到重建的图像，具体扫描点位置已经在图4中标出。结果在图4中给出。

由仿真实验结果可以看出，本发明的图像重建方法得到的光吸收分布图和组织实际的光吸收分布图非常接近，在有限角度扫描和直线扫描方式下，都能成功地进行光声成像图像重建。

3、以上述的仿真条件，比较在圆周扫描情况下，不同迭代算法（本发明方法，全变分参数梯度下降法和离散余弦变换法）重建相同图像所用时间，迭代终止的条件均设定为重建图像的PSNR值达到20 dB。

结果见表1，本发明方法在运算时间上相比其他迭代算法有比较大的提升，改善了算法效率，降低了运算量。

对本发明提出的光声图像重建方法在实验平台上进行离体实验。验证本实验方法在具体实验环境中的有效性。

离体实验组织是将凝胶加热后冷却制成，直径为50 mm，内部嵌入两段黑色橡胶条作为光吸收体，长度分别为20 mm和12 mm。实验平台采用的激光源为Nd:YAG激光器(Continuum,Surelite I)，波长为532 nm，脉冲重复频率为10 Hz，脉冲宽度为7 ns，单脉冲能量为20 mJ。采用浸入式非聚焦超声换能器(Panametric, V383-SU)接收超声信号，中心频率为3.5 MHz，有效带宽为1.12 MHz，有效直径为9.525 mm。超声信号通过脉冲接收器(Panametric, 5900PR)放大后送入示波器进行采样。示波器采用Agilent的54622D型数字示波器，最高采样率为200 MS/s。步进电机的型号为GCD-0301M型数控转台，精确控制换能器采样的角度。激光器、步进电机和示波器分别通过RS232、USB和PCI-GPIB接口卡与计算机连接。探测器在一个位置采集5个信号后旋转到下一个位置，进行全角度的圆周扫描，扫描半径42mm，一共旋转60次，信号传输到计算机后，用本发明方法进行图像重建。

实验结果表明，本发明的重建方法在实际应用中得到的重建图像和原始图像基本一致，说明本发明方法能在实际离体实验中精确地重建出组织的光吸收分布。

综上所述，本发明中的基于图像块稀疏系数的快速光声成像图像重建方法与现有的其他方法相比，有效地降低迭代算法的复杂度，同时还具有较好的重建质量，对于光声成像的图像重建具有实际意义。

表1

算法	运算时间/秒
		全变分参数梯度下降法	120.07
离散余弦变换法	42.32
		本发明方法	29.18

参考文献：

[1] C. Li, and L. V. Wang, “Photoacoustic tomography and sensing inbiomedicine,” Phys. Med. Biol., vol. 5, pp. R59~R97, Sep. 2009.

[2] M. Xu, and L. V. Wang, “Photoacoustic imaging in biomedicine,”Rev. Sci. Instrum., vol. 77, no. 4, pp.041101-1-041101-22, Apr. 2006.

[3] A. Karabutov, V. A. Andreev, B. A. Bell, R. D. Fleming, Z.Gatalica, et al., “Optoacoustic images of early cancer in forward andbackward modes,” In Proc. SPIE, vol. 4434, pp. 13-27, Jun. 2001.

[4] R. G. M. Kolkman, E. Hondebrink, W. Steenbergen, and F. F. M.Mul, “In vivo photoacoustic imaging of blood vessels using an extreme-narrowaperture sensor,” IEEE J. Sel. Top. Quantum electron., vol. 9, no. 2, pp.343-346, Mar. 2003.

[5] X. Wang, Y. Pang, G. Ku, X. Xie, G. Stoica and L. Wang, “Non-invasive laser-induced photoacoustic tomography for structural and functionalimaging of the brain in vivo,” Nat. Biotechno.,vol. 21, no.7,pp. 803–806,Jun. 2003.

[6] H. F. Zhang, K. Maslov, G. Stoica, and L. V. Wang, “Functionalphotoacoustic microscopy for high-resolution and noninvasive in vivoimaging,” Nat. Biotechno., vol. 24, no. 7, pp. 848-851, Jul. 2006.

[7] R. A. Kruger, P. Liu, Y. Fang, and C. R. Appledom, “Photoacousticultrasound (PAUS)-reconstruction tomography,” Med. Phys., vol. 22, no. 10,pp. 1605-1609, Oct. 1995.

[8] M. Xu, and L.V. Wang, “Pulsed-microwave-induced thermoacoustictomography: Filtered back-projection in a circular measurementconfiguration,” Med. Phys., vol. 29, no. 8, pp.1661~1669, Jul. 2002.

[9] M. Xu, and L. V. Wang, “Time-domain reconstruction forthermoacoustic tomography in a spherical geometry,” IEEE Trans. Med. Imaging,vol. 21, no. 7, pp. 814-822, Jul. 2002.

[10] C. Zhang, and Y. Y. Wang, “Deconvolution reconstruction of full-view and limited-view photoacoustic tomography: a simulation study,” J. Opt. Soc. Am. A, vol. 25, no. 10, pp.2436~2443, Sep. 2008.

[11] A. Modgil and P. J. La Rivière, “Implementation and comparisonof reconstruction algorithms for 2D optoacoustic tomography using a lineararray，” In Proc. SPIE, vol. 6856, pp. 13-27, Jan. 2008.

[12] G. Paltauf, J.A. Viator and S.A. Prahl, “Iterativereconstruction algorithm for optoacoustic imaging,” J. Opt. Soc. Am. A vol.112, no. 4, pp.1536–1544, Apr. 2002.

[13] Z. Guo, C. Li, L. Song and L. Wang, “Compressed sensing inphotoacoustic tomography in vivo,” J. of Biomed. Opt. vol.15, no.2, Apr. 2010.

[14] J. Provost and F. Lesage, “The Application of compressed sensingfor photo-acoustic tomography,” IEEE Trans. Med. Imaging, vol.28, no.4,pp.585–594, Apr. 2009.

[15] Y. Zhang, Y. Wang, and C. Zhang, "Total variation based gradientdescent algorithm for sparse-view photoacoustic image reconstruction,"Ultrasonics, vol.52, pp.1046 Aug. 2012.

[16] Y. Zhang, Y. Wang, and C. Zhang, "Efficient discrete cosinetransform model–based algorithm for photoacousticimage reconstruction,"Journal of Biomedical Optics,vol.18, no.6, pp.066008, Jun. 2013。

Claims (5)

1.基于图像块稀疏系数的光声成像图像重建方法，其特征在于，通过计算图像块稀疏系数，对重建图像进行修正并迭代，结合Barzilai-Borwein梯度下降法，获得最终的重建光声图像，具体步骤为：

在光声成像中，光声成像的基本方程：

(1)

其中是位置处的声压，是成像组织的光吸收分布图，t是时间，I(t)是激光脉冲能量函数，c是生物组织中的声速，和分别是生物组织的等压膨胀系数和比热容；

使用格林函数对方程(1)进行求解，对于某采样点，有：

(2)

记采样点处实际采样得到的光声信号的积分与采样时间的乘积为：

(3)

将图像和采样信号积分分别离散化，并记成矢量形式；设重建图像的大小为，X、Y分别为图像的行数和列数，则重建图像的总像素为N=XY，即矢量化后的图像可记为长度为N的列矢量u；若采样点个数为Q，每个采样点的信号长度为M，则(3)式可写成：

(4)

其中，是第i个采样点的光声信号积分与采样时间的乘积矢量；是第i个采样点的采样矩阵，其计算步骤为：

(a) 先计算大小为的矩阵：

(5)

其中，是中的序号，是采样点的坐标，dx是图像相邻像素间的实际距离，dt是离散时间步长；

(b)将矩阵矢量化得到一个N维的列矢量，作为采样矩阵的第j个列矢量；j=1～M；

(c)计算M次后得到第i点的采样矩阵；

(d)重复步骤(a)~(c)，得到Q个采样矩阵，i=1～Q；将它们联立起来，得到总的采样矩阵A：

(6)

于是，式(5)归纳为：

(7)

其中f、A和u的大小分别为、和；

基于迭代的方法求出重建的光声图像u 。

2.根据权利要求1所述的基于图像块稀疏系数的光声成像图像重建方法，其特征在于，基于迭代的方法求出重建的光声图像u，具体步骤为：

(1) 输入原始重建图像，设置各参数的初始值；

(2) 设置阈值TH, 对数据和采样矩阵进行筛选；

(3) 根据上一次迭代得到的图像，按照基于块稀疏系数迭代算法计算新的重建图像；

(4) 更新迭代中使用的参数；

(5) 判断是否达到迭代结束条件，若未达到则返回步骤(3)；若达到就结束迭代，得到重建图像。

3.根据权利要求2所述的基于图像块稀疏系数的光声成像图像重建方法，其特征在于，其特征在于步骤(2)中，对于采样数据和采样矩阵进行离散余弦变换：

(8)

其中，f’是实际采样得到的光声信号的积分与采样时间的乘积矢量经过离散余弦变换之后得到的值，f是实际采样得到的光声信号的积分与采样时间的乘积矢量，D是离散余弦变换矩阵，A’是光声信号的采样矩阵经过离散余弦变换之后得到的矩阵，A是光声信号的采样矩阵；

设定一个阈值TH，取出f’中大于TH的值组成新矢量b，从A’中取出与之相对应的行，组成新的采样矩阵W。

4.根据权利要求2所述的基于图像块稀疏系数的快速光声成像图像重建方法，其特征在于步骤(3)中，对于信号矢量b，图像矢量u和采样矩阵W都平均分成K块，对于第a块重建对象的更新公式为：

(9)

其中，u ^k [a]和u ^k+1[a]表示第k次和第k+1次迭代得到以矢量形式表示的重建图像的第a块，r _a ^k 表示第a块图像第k次迭代的步长参数，W[a]表示第a块图像所对应的采样矩阵，u ^k 为第k次迭代得到的以矢量形式表示的重建图像，μ表示迭代残差参数，T表示为矩阵的转置，max{·}表示取集合中最大数的运算符，y表示为了表达简便使用的中间参数；a=1,2,…,K。

5.根据权利要求4所述的基于图像块稀疏系数的快速光声成像图像重建方法，其特征在于步骤(4)中，所述使用的对迭代步长的更新公式为：

(10)

其中，r _a ^k+1表示第a块图像第k+1次迭代的步长参数，u ^k+1为第k+1次迭代得到的以矢量形式表示的重建图像。