CN102306385B - 任意扫描方式下光声成像的图像重建方法 - Google Patents

Abstract

本发明属于光声成像技术领域，具体为一种适用于任意扫描方式光声成像的图像重建方法。本发明通过计算重建图像的投影信号与实际信号的残差，对重建图像进行修正并迭代，在迭代过程中综合考虑了图像全变分稀疏的性质，结合全变分梯度下降法，获得重建图像。本发明能在任意扫描方式情况下精确地重建出光声图像，具有实际使用意义。

Description

任意扫描方式下光声成像的图像重建方法

技术领域

本发明属于光声成像技术领域，具体涉及一种适用于任意扫描方式光声成像的图像重建方法。

背景技术

光声成像是一种新型的生物医学成像技术，近年来发展迅速，并获得了广泛的关注^[1,2]，是当今国际学术界研究的热点前沿课题。光声成像结合了光学成像高对比度和超声成像高穿透深度的优点^[2]，并被尝试用于肿瘤检测^[3]、血管成像^[4]等领域。另外光声成像不仅适用于结构成像，还适用于功能成像，在小鼠的脑部进行了相关实验^[5]。

在光声成像中，利用短脉冲激光照射生物组织，使得组织吸收入射光能量，产生热膨胀进而向外发出超声波，通过超声换能器在不同位置扫描并采集超声信号，使用图像重建方法计算出组织的光吸收分布^[1]。在光声成像中，图像重建是关键技术。目前针对圆周扫描已提出逆Radon变换重建方法^[6]、滤波反投影法^[7]、时域重建法^[8]和反卷积重建法^[9]等；针对直线扫描提出了DAS法^[10]和二维重建法^[10]等。

上述图像重建方法在实际应用中，都局限于特定的扫描方式。本发明提出的图像重建方法，对于任意扫描方式下的光声成像，都能精确地重建出光声图像。

发明内容

本发明的目的在于提出一种可适用于任意扫描方式，并能精确地重建出光声图像的光声成像图像重建方法。

本发明提出的光声成像图像重建方法，具体步骤包括两层迭代：先进行残差迭代，即通过上一次迭代的图像结果获得模拟信号，并计算模拟信号与实际信号的残差，进而修正重建图像，对所有扫描点都进行此迭代。当所有扫描点完成迭代后，通过全变分梯度下降法对图像进行调整；重复上述两部分过程，设置误差大小为迭代结束条件，经过第二层迭代来获得最终的重建图像。下面作具体描述。

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

                             (1)

其中p(r, t)是位置r处的声压，A(r)是电磁波吸收分布，t是时间，I(t)是激光脉冲能量函数，c是生物组织中的声速，β和C _p 分别是生物组织的等压膨胀系数和比热容。光声成像图像重建是一个典型的逆问题，即：如何由

求出

。

通过使用格林函数求解方程(1) ^[7]，得到：

                  (2)

式(2)表示了光声信号与生物组织光吸收系数的关系。将(2)式变形得到：

                  (3)

记：

                                       (4)

将图像和采样信号离散化，重建光吸收分布图像A，大小为

，信号的长度为T，把图像A重排为

维的列向量

，(3)式可以转换为矩阵形式：

                      (5)

其中g、W和A的大小分别为、

和

。

这里W _i 为第i采样点的采样矩阵，对每一个采样点i都计算采样矩阵W _i ，先取W _i (t)为与图像A相同大小的矩阵，然后进行如下计算：

      (6)

其中

，(i, j) 是W _i (t)中坐标，(i ₀, j ₀)是采样点坐标，c是声速，dx是像素之间的实际距离，dt是离散信号的时间步长。将计算结果重排为

维的列向量，作为采样矩阵W _i 的第t个列向量，计算T次得到W _i 。可以看到，各个采样点是相互独立的，且不存在任何几何关系的约束，故本发明中的图像重建方法能适用于任意扫描方式。

迭代初始值取为零矩阵，分别根据(5)式和(6)式计算

和W- _i  (i=1,2,…m)，迭代公式为：

        (7)

其中m是采样点个数，

是重建图像，

是第i个采样点的采样矩阵，

是第i点的实际信号；

在全部采样点完成迭代之后，先将

中所有的负值置0，然后用全变分梯度下降法对图像进行处理，方法如下：

计算

，为上一步迭代调整前和调整后图像之间的欧氏距离，将

恢复成图像矩阵A，使用梯度下降法对图像进行调整，计算公式为：

                       (8)

其中a为常数，在实际应用中可取为0.2，用作调节残差迭代和梯度下降法的相互影响。全变分偏导的计算公式^[11]如下：

                (9)

其中

为一很小的正数，避免出现分母为0的情况，一般取为10^-8即可。

全变分梯度下降法的结果作为式(7)的初始值继续进行迭代，设置误差大小为迭代结束条件，通过这两层迭代计算可以得到重建图像。

本发明计算的流程图如图1所示。

附图说明

图1、本发明计算方法流程图。

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

图3、各扫描情况下的图像重建结果，(a)圆周扫描，(b)直线扫描，(c)不规则扫描。

具体实施方式

在计算机上进行了仿真实验。首先建立待成像组织的模型，设定已知的光吸收分布图，并根据(2)式分别计算圆周扫描、直线扫描和不规则扫描情况下采集得到的光声信号，圆周扫描半径为42 mm，角度步长为18°，共20个采样点；直线扫描间隔为4.42 mm，共20个采样点；不规则扫描一共有20个采样点，分布在组织的周围，不符合任何几何形状。用仿真的超声信号重建出待成像组织的光吸收系数图。

图2显示了待成像组织的光吸收系数图，组织大小为89.6 mm×89.6 mm，重建图像大小为128×128像素。

图3显示了用本发明重建的光吸收系数图，图3(a)为在圆周扫描情况下的重建图像，图3 (b)为再直线扫描情况下的重建图像，图3(c)为在不规则扫描情况下的重建图像。

由实验结果可见，本发明的重建图像和原始图像非常接近，说明本发明能精确地重建出光吸收分布图，并且在任意扫描方式下，都能成功地进行光声成像图像重建，使得图像重建不再受到扫描方式的限制，具有实际意义。

本发明受国家自然科学基金(10974035)和上海市优秀学科带头人计划(10XD1400600)资助。

参考文献

(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. 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)    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.

(6)    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.

(7)    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.

(8)    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.

(9)    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.

(10) D. 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.

 E. Y. Sidky, C .M. Kao, and X. Pan, “Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT”, J. X-ray Sci. Technol. Vol. 14, no. 2, Jun. 2006.。

Claims (1)

1.一种适用于任意扫描方式光声成像的图像重建方法，其特征在于：通过上一次迭代的图像结果获得模拟光声信号，并计算模拟光声信号与实际信号的残差，利用此残差修正重建图像，对所有扫描点都进行此迭代；当所有扫描点完成迭代后，计算重建图像的全变分，并通过全变分梯度下降法对图像进行调整；设置误差大小为结束条件，重复上述迭代和梯度下降法两部分过程，重建光声图像；其中：

所述迭代的算式为：

         (1)

其中，m是采样点个数，

是重建图像，

是第i个采样点的采样矩阵，是第i点的实际信号；其中

的算式为：

                         (5)

W _i 的算式为： 

      (6)

式（6）中，，(i, j) 是W _i (t)中坐标，(i ₀, j ₀)是采样点坐标，c是声速，dx是像素之间的实际距离，dt是离散信号的时间步长；将计算结果重排为

维的列向量，作为采样矩阵W _i 的第t个列向量，计算T次得到W _i ；

将

恢复成图像矩阵A，所述的全变分梯度下降法算式为：

                     (2)

式（2）中a为常数，

，为所述迭代的重建初始值与重建结果的欧氏距离；

公式(2)中全变分偏导的计算公式为：

             (3)；

其中，

取10^-8。

Images