CN103279966A

CN103279966A - 基于图像稀疏系数p范数和全变分参数的光声成像图像重建方法

Info

Publication number: CN103279966A
Application number: CN2013102130838A
Authority: CN
Inventors: 汪源源; 张晨; 张砚
Original assignee: Fudan University
Current assignee: Fudan University
Priority date: 2013-06-02
Filing date: 2013-06-02
Publication date: 2013-09-04

Abstract

本发明属于光声成像技术领域，具体为一种基于图像稀疏系数p范数和全变分参数的光声成像图像重建方法。本发明通过计算图像的全变分参数和稀疏表示系数的p范数(0<p<1)，对图像进行迭代修正，随后再根据迭代过的图像重新计算上述参数，进行反复迭代直至达到结束条件，最后得到重建的光声图像。本发明方法能提升重建的光声图像的质量，并且收敛速度快，对于光声成像的图像重建具有实际意义。

Description

基于图像稀疏系数p范数和全变分参数的光声成像图像重建方法

技术领域

本发明属于光声成像技术领域，具体涉及一种基于图像稀疏系数p范数和全变分参数的光声成像的图像重建方法。

技术背景

光声成像是一种基于光声效应建立的新型医学成像方法，结合了光学成像的高对比度特性和超声成像的高成像深度的特性，并且由于其非电离波的特性，不会在检测的过程中对人体产生伤害^[1,2]。现在主要应用领域有乳腺肿瘤检测^[3]、血管成像^[4]和脑损伤探测^[5]等。同时由于所成图像是人体组织的光吸收特性，因而一定程度上反映了成像组织与光吸收特性相关的病理学特性^[2]。利用这一特性可以将光声成像应用于功能成像^[5,6]，如血液动力学检测^[5]等领域。

在光声成像中，使用短时脉冲激光照射成像组织。一部分被组织吸收的光能将会被转化为热能，使附近的组织发生热弹性膨胀，从而形成超声波发射。这一超声波可以用超声换能器检测，通过超声换能器在不同位置扫描可以采集到光声信号，随后使用图像重建算法就可以计算出组织的光吸收分布。目前针对圆周扫描有逆Radon变换重建方法^[7]、滤波反投影法^[8]、时域重建法^[9]和反卷积重建法^[10]等方法；针对直线扫描有DAS法^[11]和二维重建法^[11]等方法。上述方法无需迭代就可直接得到结果，但是在采样点较稀疏时重建图像的精度不高，并且受限于特定的扫描方式。针对这一现象，提出了多种迭代重建方法^[12-14]，这些方法的优势在于成像精度高、不受扫描方式的影响，但是缺点是速度较慢。

与上述图像重建方法相比，本发明中的图像重建方法使用了图像稀疏系数p范数和全变分参数，提升了重建图像的质量，同时由于使用了Barzilai-Borwein梯度下降法进行图像重建，成像速度性能相比于其他迭代方法也有优势。

发明内容

本发明的目的在于提出一种图像质量好，成像速度快的光声成像的图像重建方法。

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

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

Figure 2013102130838100002DEST_PATH_IMAGE001

(1)

其中

是位置

处的声压，

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

求出

。

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

，有：

(2)

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

Figure 2013102130838100002DEST_PATH_IMAGE009

(3)

记采样点

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

(4)

实际应用时，可将

和

Figure 2013102130838100002DEST_PATH_IMAGE011

分别离散化，并记成矢量形式。若重建图像的大小为

(X，Y分别为图像的行数和列数)，则重建图像的总像素为N (N=XY)，即矢量化后的图像可记为长度为N的列矢量u。若系统的采样点个数为Q，每个采样点的信号长度为M，则可将(3)式写成：

Figure 2013102130838100002DEST_PATH_IMAGE013

(5)

其中

是第i个采样点的光声信号积分与采样时间的乘积矢量；

Figure 2013102130838100002DEST_PATH_IMAGE015

是第i个采样点的采样矩阵，其计算步骤为：

(a) 先计算大小为

的矩阵

：

(6)

其中

，

Figure 2013102130838100002DEST_PATH_IMAGE019

是

中的序号，

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

(b) 将矩阵

矢量化得到一个N维的列矢量，作为采样矩阵

的第j个列矢量。

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

。

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

(7)

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

(8)

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

Figure 2013102130838100002DEST_PATH_IMAGE023

、

和

Figure 2013102130838100002DEST_PATH_IMAGE025

。

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

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

(1) 设定原始重建图像为

，各参数的初始值设为：全变分参数在优化问题中所占的权重

Figure 2013102130838100002DEST_PATH_IMAGE027

，p范数在优化问题中所占权重

，限制条件在优化问题中所占的权重，初始的联合迭代步长参数、全变分参数迭代步长参数和p范数迭代步长参数，范数

Figure 2013102130838100002DEST_PATH_IMAGE031

，迭代终止阈值

。

(2) 根据上一次迭代得到的图像，分别计算图像的全变分参数和稀疏表示系数的p范数(0<p<1)。全变分参数的迭代公式为：

(9)

其中

是第k次迭代得到以矢量形式表示的重建图像，

Figure 2013102130838100002DEST_PATH_IMAGE035

和

分别是第k次和第k+1次迭代时第i个像素点的全变分参数；

Figure 2013102130838100002DEST_PATH_IMAGE037

是第i个像素点的全变分参数计算矩阵；

为第k次迭代时第i个像素点全变分参数迭代步长参数，为第k次迭代时联合迭代步长参数；α为表示全变分参数在优化问题中所占的权重，ρ为限制条件在优化问题中所占的权重；

、

、

和是为了表达简洁所引入的中间变量。

图像的稀疏系数的p范数的更新方法为：

(10)

其中

和

为第k次和第k+1次迭代时第i个像素点的稀疏表示系数；

是第i个像素点的稀疏变换矩阵，T表示矩阵的转置；

为第k次迭代时p范数的迭代步长参数，p为范数，β为p范数在优化问题中所占权重；

、

、

和

是为了表达简洁所引入的中间变量。

(3) 根据得到的全变分参数和稀疏表示系数的p范数，计算新的重建图像。其计算方法为：

(11)

其中

是第k+1次迭代得到以矢量形式表示的重建图像；

为第k+1次迭代得到以矢量形式表示的整幅图像的全变分参数，由每个像素点的全变分参数组合而成；

表示第k+1次迭代得到以矢量形式表示的整幅图像的稀疏表示系数，由每个像素点的稀疏表示系数组合而成；F为傅里叶变换矩阵，D为全变分参数计算矩阵，A为光声信号的采样矩阵，f是实际采样得到的光声信号的积分矢量与采样时间的乘积组成的矢量；I为单位阵。

(4) 更新迭代中使用的参数。其计算方法为：

(12)

其中

和分别为第k+1次迭代时第i个像素点全变分参数迭代步长参数和p范数迭代步长参数，

为第k+1次迭代时的联合步长参数。

(5) 迭代步数k变为k+1，判断是否达到迭代结束条件，若未达到则返回步骤(2)；若达到就结束迭代，得到重建的光声图像。具体的判断方式为：

(13)

图1给出了使用本发明中的方法进行光声图像重建的具体流程图。

与现有技术相比，本发明引入了全变分参数和稀疏系数的p范数，提升了图像的重建质量，同时使用了Barzilai-Borwein梯度下降法，算法收敛速度快，具有很高实际应用价值。

附图说明

图1、本发明基于图像稀疏系数p范数和全变分参数的光声成像图像重建方法具体流程图。

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

图3、各扫描方式下的图像重建结果，(a)圆周扫描，(b)直线扫描。

图4、对组织进行圆周扫描，采样点数为20个的情况下各种光声图像重建方法的结果比较。(a) 滤波反投影法，(b) L1范数法，(c) 全变分参数梯度下降法，(d) 本发明方法。其重建结果的峰值信噪比(PSNR)分别为14.34 dB、30.23 dB、33.17 dB、35.86 dB。

图5、圆周扫描采样角情况下，L1范数法、全变分参数梯度下降法和本发明方法在相同迭代次数时的MSE值。

具体实施方式

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

1、确定组织的光吸收分布图，如图2所示，组织大小为89.6 mm

89.6 mm，重建图像大小为128

128像素。根据(2)式采集光声信号，分别采用了圆周扫描和直线扫描两种扫描方式：圆周扫描半径为42 mm，角度步长为

，共20个采样点；直线扫描间隔为4.2 mm，共20个采样点。使用本发明中的方法对采集到的光声信号进行重建，得到重建的图像。结果在图3中给出。

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

2、以上述的仿真条件，比较本发明方法和滤波反投影法、L1范数法和全变分参数梯度下降法的图像重建结果，结果在图4中给出。选取峰值信噪比PSNR(Peak Signal to Noise Ratio)为量化指标，单位为dB。PSNR值越大，图像的重建效果越好。

(14)

其中

为给定的原始图像，

为重建的图像。

仿真结果表明，本发明方法的PSNR值高于其他方法，说明本发明中的重建方法的成像质量要优于其他发法。

3、以上述的仿真条件，比较相同迭代次数时，本发明方法与L1范数法和全变分参数梯度下降法的图像重建结果，以此来比较不同方法间的收敛速度，结果在图5中给出。选取平均均方误差MSE(mean square error)为量化参数，MSE值越小，光声图像的重建效果越好。

(15)

仿真结果表明，在相同的迭代次数下，本算法的MSE值低于其他算法，说明本发明中的重建方法的收敛速度要优于其他算法。

综上所述，本发明中的基于图像稀疏系数p范数和全变分参数的光声成像图像重建方法与现有的其他方法相比，可以有效地提高重建光声图像的成像质量，并且提升图像重建方法的收敛速度。

参考文献：

[1] C. Li, and L. V. Wang, “Photoacoustic tomography and sensing in biomedicine,” 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 and backward 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-narrow aperture 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 functional imaging 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, “Functional photoacoustic microscopy for high-resolution and noninvasive in vivo imaging,” Nat. Biotechno., vol. 24, no. 7, pp. 848-851, Jul. 2006.

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

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

[9] M. Xu, and L. V. Wang, “Time-domain reconstruction for thermoacoustic 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 comparison of reconstruction algorithms for 2D optoacoustic tomography using a linear array，” In Proc. SPIE, vol. 6856, pp. 13-27, Jan. 2008.

[12] G. Paltauf, J.A. Viator and S.A. Prahl, “Iterative reconstruction 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 in photoacoustic tomography in vivo,” J. of Biomed. Opt. vol.15, no.2, Apr. 2010.

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

X. Ye, Y. Chen and F. Huang, “Computational Acceleration for MR Image Reconstruction in Partially Parallel Imaging,” IEEE Trans. Med. Imaging, vol.30, no.5, pp.1055–1063, May. 2011。