Matrix Compression and Compressed Sensing Reconstruction for Photoacoustic Tomography

Online medical diagnosis is an emerging type of medical services, also it is a new kind of service for electronic commerce. Photoacoustic tomography is a promising imaging modality that can provide high contrast and spatial-resolution images of light-absorption distribution in tissue. This paper presents novel matrix compression methods for reducing required memory and accelerating reconstruction speed for model-based reconstruction. In addition, we integrated the compressed sensing into our proposed reconstruction method for solving the optimization problem. The results of phantom experiment indicate that the proposed method provides a faster and high quality reconstruction with less memory consumption. Therefore, it makes the photoacoustic tomography suitable for online or remote medical diagnosis.


I. INTRODUCTION
With the rapid development of computer, network, and biomedical technology, online or remote medical diagnosis becomes an emerging type of medical service, which is a novel kind of service for electronic commerce.Photoacoustic imaging is a promising image modality which can provide early detection of cancer especially breast cancer.Worldwide, breast cancer comprises 22.9 % of all cancers in women.In 2008, breast cancer caused 458,503 deaths worldwide.Therefore, better online or remote diagnosis will improve the quality of life, and also it is a big market for electronic commerce.In photoacoustic tomography (PAT) [1], by illuminating a target with a short pulse laser, stress waves are produced due to the thermoelastic expansion.Acoustic sensors are placed at surrounding positions to detect PA signals, and then the absorption source is recovered from the detected PA signals through a reconstruction algorithm.Because light energy is converted to ultrasound which has much less absorption and scattering than light, PAT has much better spatial resolution than traditional optical modalities at depths exceeding the optical ballistic regime.Reconstruction algorithms for PAT have been extensively studied in recent decades [1].The simplest reconstruction method is delay and sum (DAS) [2].Filtered back-projection (FBP) method [3] is a common tomography reconstruction method which can also be applied to PAT.Besides the approximation methods, analytic reconstruction methods with fewer assumptions have been proposed to obtain exact reconstruction, for example, universal back-projection (UBP) method [4] and Fourier domain reconstruction method [5].The above reconstruction methods require enclosed detection for circular scanning or an unbounded open surface for planar scanning, which are generally difficult to implement in a clinical situation.To overcome the limited-view problem, model-based reconstruction methods [6]- [12] (also called iterative reconstruction methods) have been investigated.In these methods, the inverse problem is converted into solving an optimization problem by minimizing the error between detected PA signals and calculated signals from a forward model.Although this type of methods needs more computation, it needs fewer detecting sensors and thus less acquisition time.It is also able to model non-ideal physical conditions and measurement environments.Therefore, the degradation of reconstructed images caused by acoustic inhomogeneity and attenuation can be resolved.Generally, model-based reconstruction requires a large amount of memory and a long calculation time.
In order to overcome limitations of model-based reconstruction, we propose novel matrix compression methods to reduce the required memory of system matrix.By utilizing these methods, the matrix can be compressed to 1/250 its size; hence, the proposed reconstruction method can be implemented with a conventional computer at a high speed.Conjugate gradient method is a popular method for solving the optimization problem in model-based reconstruction, however, in some cases this method converges to local minima.In addition, reconstructed images usually is sparse, the compressed sensing method can reconstruct better images when using less scanned data.In order to resolve the problem and improve the quality of reconstruction, we integrated the compressed sensing reconstruction to our proposed method.

II. RECONSTRUCTION METHOD AND MATRIX COMPRESSION
A PA wave generated in an acoustically homogenous and non-viscous medium can be described as [13] Matrix Compression and Compressed Sensing Reconstruction for Photoacoustic Tomography Shuhui Bu 1 , Zhenbao Liu 1 , Tsuyoshi Shiina2 , Kazuhiko Fukutani where p(r,t) is pressure, H(r,t) is a heating function defined as the thermal energy converted, C p is the isobaric specific heat, β is the coefficient of isobaric volume expansion, and v s is acoustic speed.Under the conditions of both thermal and stress confinement, the heating time can be treated as a delta function, such as , where A(r,t) denotes absorbed energy.Wave equation ( 1) can be solved by using a Green function approach [1], and the solution is The above equation describes a relationship in which the detected PA pressure at position r and time t comes from sources over a spherical surface centered at r with a radius of |r-r'|.Based on the equation ( 2), the forward model for calculating PA signals can be expressed in a matrix form as where H denotes the system matrix representing the geometric relationship between the initial pressure and the detected PA signals, f represents absorbed energy, and p denotes the detected PA signals.
Although model-based reconstruction for PAT has many advantages, it has a critical limitation that the system matrix consumes huge amounts of memory without matrix compression.In order to apply the model-based reconstruction method to clinical situations, we researched matrix-compression methods to reduce the memory requirement.First, matrix H is a sparse matrix in which more than 99% of the entries are zero.The required memory can be greatly reduced by using the compressed sparse row (CSR) format, where non-zero entries are stored in continuous memory locations and corresponding column indices are stored in an integer array.We use another integer array to store the index of first non-zero entry in each row.
In PAT, the ultrasound transducer is generally assumed to detect pressure waves in directions from 0 to 90 degrees with the same sensitivity, however, in real situations the sensitivity is not uniform due to the finite-sized unfocused transducer.The spatial impulse response of the transducer introduces limited detection angles, waveform distortion, and time-delay errors for reconstruction.For example, the angular sensitivity of a 2x2 mm 2 unfocused rectangular transducer is plotted in Fig. 1(a).When the incident angle is 45 degrees, the sensitivity is -50dB, which is less than the detectable signal level for a conventional sensor.Hence, the incident wave can be ignored.In order to compress the matrix H, if the incident angle exceeds a given threshold, its value is set to zero so that the entry is not stored.With this approach, matrix H can be further decreased to about half its original size when the threshold value is 45 degrees.
The matrix can be further compressed to one-fourth its original size based on the symmetry of wave propagation [10], as illustrated in Fig. 1(b).This figure depicts only four sensors and eight voxels.The PA signals detected by sensor 1 are calculated from the inner product H(1,:)f.Because sensors 2 and 1 are symmetrical about the Y axis, the second row of H can be obtained by changing the X-axis order of the first row.Similarly, the third row of H can be obtained by changing the Y-axis order of the first row, and the fourth row of H can be obtained by changing the X-and Y-axis orders of the first row.Therefore, only one-fourth of the rows of H are necessary for reconstruction, and other rows' data can be obtained during the multiplication of sparse matrix and vector.In total, the required memory can be reduced to approximately 1/250 its original size by using the above methods.The DREAM toolbox [14] was used to calculate the result.(b) Matrix compression method using the symmetry of wave propagation.For simplicity, in this figure only four sensors and eight voxels are illustrated.Only the first 1/4 row of H is stored; coefficients in other rows can be obtained by rearranging the order of the first 1/4 row's coefficients.
In the system matrix m x n  HR , usually the number of scanned data m is not large enough compared to the number of voexl n, simple least-squares method leads to under-determinate.In order to overcome the problem, Tikhonov regularization can be applied.However, this method does not consider the sparseness of the imaging objects.The compressed sensing theory provides a better solution which can eliminate under-determinate by incorporating sparsity constraints.In the CC method, the image is projected onto an appropriate basis set W, such as discrete cosine transformation (DCT) or wavelet basis.And therefore, the reconstruction of image f is obtained by solving the following constrained optimization problem where f*=Wf denote the coefficients of image data f.In this research we adopt l 1 -regularized least squares [15] to solve the convex optimization problem.The optimized regularization parameter  is initialized set to 0.05, and optimally updated by the method Guo et al., [16] proposed.All programs were implemented in parallel and run on a cluster system containing six workstations.Each workstation had two quad-core Xeon 2.4GHz CPUs.The floating point operation per second (FLOPS) of the cluster system was 251 G.After the simulated PA signals were obtained, the reconstruction program was used to reconstruct the absorption distribution.The required memory for the matrix H before compression was 1.043 TB; after compression the size was reduced to 3.250 GB.The calculation time for one iteration was 0.9 seconds, and there were 50 iterations in total.Fig. 4 presents the maximum intensity projection (MIP) images reconstructed by the proposed method, compared with the conventional UBP method.Here, MIP X-Y is the projection along the Z axis, MIP X-Z is the projection along the Y axis, and MIP Z-Y is the projection along the X axis.Fig. 4(a) present MIP images reconstructed by the UBP method.The MIP images reconstructed by the proposed method are presented in Figs.4(b).The resulting images indicate that the absorption distribution reconstructed using the proposed method is better than that using the UBP method.Noise and artifacts were produced in the UBP reconstruction.In contrast, these were greatly reduced by the proposed method.The amount of required memory was greatly reduced due to the proposed matrix compression methods.In addition, because imaging object sparse constraint was considered and l 1 -regularized least squares was adopted for solving the optimization problem, better image quality is achieved.

IV. CONCLUSIONS
In this paper, we presented novel system matrix compression methods for 3-D planar PAT model-based reconstruction.Matrix compression methods made the proposed method applicable to 3-D reconstruction under planar measurement conditions, and the calculation speed was also accelerated.In addition, the quality of reconstructed image had significant improvement due to the integrating the compressed sensing method.The phantom experiment results indicated that the proposed method reconstructs better quality images with a high speed and less memory requirement.We expect that the proposed method will benefit the development of PAT for clinical diagnosis which could provide a new type of service for electronic commerce.

Fig. 1 .
Fig. 1.Matrix compression methods.(a) Angular sensitivity of a 2x2 mm 2 unfocused rectangular transducer.The distance between transducer and target is L=10 mm.The horizontal axis represents the wave incident angle θ (in degree), and the vertical axis represents the sensitivity of the transducer.The DREAM toolbox [14] was used to calculate the result.(b) Matrix compression method using the symmetry of wave propagation.For simplicity, in this figure only four sensors and eight voxels are illustrated.Only the first 1/4 row of H is stored; coefficients in other rows can be obtained by rearranging the order of the first 1/4 row's coefficients.

Fig. 4 .
Fig. 4. Results of image reconstructed from experiment data.(a) images reconstructed by the UBP method.(b) images reconstructed by the proposed method.
Manuscript received March 19, 2012; accepted May 15, 2012.This work is partly supported by the Innovative Techno-Hub for Integrated Medical Bio-imaging Project of the Special Coordination Funds, from the MEXT, Japan.This work is also partly funded by NSFC (61003137, 61202185), Natural Science Basic Research Plan in Shaanxi Province of China (Program No. 2012JQ8037), SRF for ROCS, and NWPU fundamental funds (JC201202, JC201220).