Abstract
Photoacoustic imaging is an innovative imaging technique to image biomedical tissues. The time reversal reconstruction algorithm in which a numerical model of the acoustic forward problem is run backwards in time is widely used. In the paper, a time reversal reconstruction algorithm based on particle swarm optimization (PSO) optimized support vector machine (SVM) interpolation method is proposed for photoacoustics imaging. Numerical results show that the reconstructed images of the proposed algorithm are more accurate than those of the nearest neighbor interpolation, linear interpolation, and cubic convolution interpolation based time reversal algorithm, which can provide higher imaging quality by using significantly fewer measurement positions or scanning times.
1. Introduction
Recently, photoacoustic tomography (PAT) has emerged as a powerful imaging technology for many biomedical applications [1, 2]. PAT is able to provide functional imaging of physiological parameters. Photoacoustic imaging combines electromagnetic and ultrasonic wave synergistically, providing relatively deep speckle-free imaging with high electromagnetic contrast at high ultrasonic resolution. In PAT, the developing reconstruction methods have been an important task [3, 4]. Various reconstruction methods, such as filtered backprojection, Fourier transform, alternative algorithm, time reversal, inversion of the linear Radon transform, and delay and sum beamforming, have been developed under different assumptions and approximation [5–9]. Besides, different shaped measurement surfaces will lead to different shaped final artifacts.
In a word, the image reconstruction progress can be regarded as running a numerical model of forward problem backwards in time domain, at each time step, in which the measured time-varying pressure signals are enforced (in time reversed order) as a Dirichlet boundary condition at their recorded positions that is called time reversal reconstruction method. It has been described as the “least restrictive” imaging algorithm on the basis that it relies on fewer assumptions than many other imaging reconstruction algorithms [10, 11]. The idea of time reversal to PAT was suggested by Fink and Prada [12] and later used for PAT by Xu and Wang [13], and the time reversal method was applied to deriving exact and approximate reconstruction algorithms for PAT in an arbitrary configuration and heterogeneous acoustic media.
In conventional time reversal imaging reconstruction, the recorded pressure time series are enforced in time reversed order as a Dirichlet boundary condition as the position of detectors on the measurement surface. If a sparse array of the detectors points is used to collect the measurement rather than a continuous surface, the enforced time reversed boundary condition will necessarily be discontinuous. This can cause significant blurring in the reconstructed images. To solve the problem, Treeby and Cox [14] improved time reversal image reconstruction by using interpolated sensor data. In the course, the interaction can be avoided by interpolating the recorded data onto a continuous rather than discrete measurement surface within the -space grid used for the reconstruction. The edges of the reconstructed image are considerably sharper and the magnitude has also been improved. After that, Cox and Treeby [15] used the enforced time reversal boundary condition to trap artifacts in the final image, and by truncating the data, or introducing an adaptive threshold boundary condition, this artifact trapping can be mitigated to some extent.
In this paper, what we are concerned with in the method is the artifacts elimination and interaction of image reconstruction when using “time reversal” algorithm. An optimized hybrid interpolation algorithm has been used to create a continuous surface that is spatially equivalent to the original Cartesian measurement surface via interpolation. It is demonstrated that the proposed algorithm can reduce artifact and improve the acoustic magnitude and the signal-to-noise ratio through partial correction for the discontinuous aperture.
2. Theory and Methods
2.1. Acoustic Propagation Theory in PAT
For PAT, in a lossless acoustic medium, the photoacoustic wave equation can be reformulated as an initial value problem. The photoacoustic forward problem may be written aswhere the initial conditions are given bywhere is the acoustic pressure at time and point inside the imaging region , is the acoustic particle velocity, is the sound speed, and , are the acoustic and ambient densities, respectively. Using this framework, it is straightforward to modify the adiabatic equation of state to account for acoustic absorption or nonlinear effects.
2.2. Time Reversal Image Reconstruction
In PAT, the photoacoustic image reconstruction problem is to estimate the initial pressure distribution inside the imaging region given measurement of on arbitrary measurement surface . In time reversal imaging, this estimate is achieved by using the recorded measurements of the acoustic pressure over an arbitrary surface for some time to . In this case, the initial pressure in (2) is set to zero, givingFrom the equation, the pressure in the imaging field from to can be obtained; at the same time, the reverse-time variable from to also can be computed. In a practical sense, the reconstruction is performed by using the acoustic pressure time histories measured on for to in time reversed order as an enforced (time-varying) Dirichlet boundary condition on within a numerical acoustic propagation model. Here, and both and are in silico equivalents to the real world and .
During the time reversal reconstruction, if and vary spatially, the time reversed waves will no longer interfere precisely to reproduce the original wave field, which will lead to the production of additional waves. The extra waves are called vestigial waves, which are the time reversed vestiges of the scattered waves. Vestigial waves are artifact-producing waves, in the sense that any of them remaining in the computational domain would constitute artifacts in the final PAT image. To improve the reconstruction result of the discontinuous aperture, the artifacts and interaction can be largely removed by interpolating the recorded data onto a continuous measurement surface within the -space grid used for the reconstruction. As an optimized hybrid interpolation algorithm, PSO optimized SVM interpolation algorithm can provide higher convergence rate and optimization precision, which has been used to create a continuous surface that is spatially equivalent to the original Cartesian measurement surface via interpolation to solve the problem.
2.3. PSO Optimized SVM Interpolation Algorithm
The essence of the support vector machine (SVM) training parameter selection is a process of optimization search [16]. In the process of interpolation, the index of the training model of the support vector machine, such as the fitting precision and the generalization ability, is directly related to the selection of the nuclear function, nuclear parameter, the penalty coefficient, and other parameters in the parameter optimization algorithm [17].
The particle swarm optimization (PSO) algorithm is a kind of bionic optimization algorithm, which is derived from the approximation behavior simulation of birds and fish population. It is an optimization tool based on iteration, searching for the optimal solution through the collaboration of the individual particles. The algorithm has a strong ability in global optimization. Using the particle swarm optimization algorithm to solve the optimization problem is equal to the search for the space location of a bird. The birds in the space are called “particles.” Each particle adjusts its flight path randomly based on the flight experience of its own so as to close to the optimal point finally. Different particle has different location and speed and individual fitness corresponding to the flight objective function. The flight path is adjusted by tracking two “extreme values.” One of the extreme values is called individual extremum which indicates the optimal solution of the particle itself. The other one is the global extremum, which indicates the optimal solution of the whole swarm. When the two extreme values are found, the particles update their speed and locations as follows:where is the inertia weight; is the current evolutionary iteration time; , ; and are random numbers distributed in the interval ; is the swarm size; and are the accelerating factors; is the -dimension component of the location of the th particle in the iteration; is the -dimension component of the velocity of the th particle in the iteration; is the individual extreme value; is the -dimension component of the best location of the th particle; is the global extreme value; is the -dimension component of the best location of the swarm.
Using the particle swarm optimization algorithm to solve the optimization problem is mainly discussing the optimization problem of the penalty factor , insensitive parameter , and the nuclear parameter . The progress of PSO-SVM method is shown in Figure 1. In general, the parameters of , , and are related with learning samples and practical problems. The penalty factor controls the degree of punishment when samples are misclassified and achieves a compromise between training error and model complexity. The larger is, the higher fitting precision is required. It makes training so difficult and also takes a longer time. But when is smaller, it will lower the fitting precision. Too big or too small value of will deteriorate the generalization capability of the system; the former is called “overlearning,” and the latter is called “less learning.” The insensitive parameter reflects the regression model’s sensitivity of the noise included in input variable. The larger is, the lower model fitting accuracy and complexity is also prone to overfit. The nuclear parameter represents the mean square error of Gaussian function and is the width of function in independent variables direction. When is smaller, the kernel function has better fitting performance, but it will deteriorate the generalization ability. So, using the particle swarm optimization algorithm to solve the optimization problem is mainly discussing the optimization problem of the penalty factor , insensitive parameter , and the nuclear parameter .
In the process of the SVM interpolation, the solution data in solution space cannot be solved by PSO algorithm directly. The right way is to encode it and convert to particle string structure of the searching space to solve the optimization problem. In the process of the optimization, there are three optimization parameters to be solved, which mean that the dimension of the particle location is three, as . Then the real value encoding is used, 3d particles represent the initiate location of the particles swarm, and it is randomly generated according to the optimal interval of parameters. The random number generated in is the initial speed. There exists an updating problem of different range in the parameter optimization of speed and locations of the particles swarm. At this time, we should do the normalized mapping process for the three parameters before the optimization process. And after the updating, the inverse normalized process will be done to the original scope. In the iterating process, fixed evolution algebra is set. In the solving process, the maximum value will be compared with the fixed value in each generation; the operation will be terminated when the two values reach an agreement.
3. Results and Discussion
3.1. Time Reversal Reconstruction Algorithm
The initial pressure distribution of a phantom (simulating tumor tissue) with a 256 × 256 pixel grid is shown in Figure 2(a). The conventional time reversal algorithm is used to reconstruct the PAT phantom, and the result is shown in Figure 2(b). It shows that the reconstructed PAT image has been significantly blurred due to the discontinuous boundary condition, and some artifacts have been caused by vestigial waves. Figure 3 shows a profile of the initial pressure and the conventional reconstruction pressure, which facilitates a large amplitude gap between them.
(a)
(b)
3.2. Time Reversal Reconstruction Based on PSO-SVM Interpolation Method
The algorithm above is slow for real-time applications. Recently fast PSO optimized SVM interpolation algorithm has been developed [18, 19]. These algorithms can be employed for real-time algorithm implementation. In time reversal PAT imaging, an interpolation function is often employed to compute missed Cartesian raster pixel values. The well-known interpolation algorithms are linear interpolation, cubic convolution interpolation, and cubic spline interpolation algorithm, which can be applied to reconstruct the PAT images. To prove the effectiveness and superiority of the proposed method in PAT image reconstruction, a two-dimensional photoacoustic image simulation experiment is given to illustrate the effect of the imaging performance by using various interpolation based time reversal algorithms. Therefore, the initial pressure distribution, the sound heterogeneity, and the measure surface will choose simple geometrical shape; the generation process of the photoacoustic signal will be simulated by K-Wave Matlab toolbox.
The simulation is performed on grid architecture with 5% and 15% Gaussian White Noise (GWN) added to the recorded photoacoustic data, respectively. The PC with Intel P6300 CPU and 4 GB RAM is used in the simulation. The propagation velocity of sound in the medium is 1500 m/s, the medium density is 1000 kg/m3, and the initial pressure distribution is the two locations in the spherical absorber at (). For a central position at (150, 150), the diameter is 1.3 mm and for the other central position at (90, 90), the diameter is 0.8 mm. A single array ultrasonic detector is put around the absorber selection for a round from the grid (0, 0). The collected photoacoustic signal which was measured at multiple points is used for reconstruction; the result images of different interpolations applied to time reversal algorithm using 30 measurement points are displayed in Figures 4(a)–4(f). In Figure 4, Figure 4(a) represents ideal image of light sound absorber, Figure 4(b) represents reconstruction image of conventional time reversal reconstruction algorithm, and Figures 4(c)–4(f) are reconstruction images of using time reversal reconstruction algorithm based on linear interpolation, three cubic interpolations, cubic spline interpolation, and PSO optimized SVM interpolation algorithm. When less measurement points are used, the collected photoacoustic data imperfection between each measurement point is not enough, and pseudowave is generated in the reconstruction process. It demonstrates that PSO optimized SVM interpolation is able to preserve more edges information in detailed regions and remove more artifacts (blurring) than the other approaches.
(a) Ideal image
(b) No interpolation
(c) Linear interpolation
(d) Three cubic interpolations
(e) Cubic spline interpolation
(f) PSO optimized SVM interpolation
The results of 50 measurement points are shown in Figure 5. It demonstrates that the more measurement points are, the better reconstruction image quality we get. It can be seen from the photoacoustic reconstruction images that the PSO optimized SVM interpolation algorithm can preserve more original image details and remove the artifact phenomenon better; thus it will be most effective to apply it to the time inversion algorithm. This is also evident in Figure 6, which shows a profile through the center of the two absorbers. The junction between the absorber (simulated tumor) and the background is noticeably sharper after the interpolation compared to before, especially for the PSO optimized SVM interpolation we utilized. The overall magnitude of the proposed reconstruction has also been improved through partial correction for the attention effect. The quality of the reconstruction images (with 30, 50 measurements added with 5%, 15% GWN, resp.) is measured via mean square error (MSE) and peak signal-to-noise ratio (PSNR), which are shown in Tables 1 and 2. It also demonstrates that the PSO optimized SVM interpolation is the most accurate interpolation for time reversal method to reconstruct PAT images. The experiment shows that the proposed method can significantly improve the image resolution and suppress the artifacts when less measurement points are used in a single-element photoacoustic detector imaging system.
(a) Ideal image
(b) No interpolation
(c) Linear interpolation
(d) Three cubic interpolations
(e) Cubic spline interpolation
(f) PSO optimized SVM interpolation
(a)
(b)
A tissue phantom is built and shown in Figure 7(a). The phantom is irradiated by a laser and observed; the photoacoustic images were reconstructed by filter backprojection (FBP) algorithm, conventional time reversal algorithm, and PSO optimized SVM time reversal algorithm using 50 and 100 measurements. Results depicted in Figures 7(b)–7(g) indicate that the reconstruction of the photoacoustic image from PSO optimized SVM time reversal algorithm represents the best performance, which can obtain more source-image details than conventional time reversal algorithm and has less artifacts than FBP algorithm, and this makes it a suitable algorithm for PAT.
(a) Original image
(b) FBP-50
(c) Conventional-50
(d) PSO-SVM-50
(e) FBP-100
(f) Conventional-100
(g) PSO-SVM-100
In our experiment, the observed sample is 3D polyacrylamide gel (8% acrylamide + 0.5% bisacrylamide) with graphite powder in a standard culture dish (PN: 16235-1SGP, 35 × 12 (mm); diameter: 20 mm). The sample for the optical scanning includes 4 graphite spots (different graphite granule distributions) in different layers. The transducer is one single-element unfocused transducer with 40 MHz center frequency. The image of the sample is shown in Figure 8(a). Figures 8(b)–8(e) show the photoacoustic result reconstructed by conventional time reversal algorithm and PSO optimized SVM time reversal algorithm using 50 and 100 measurements. As expected, the PSO optimized SVM reconstruction result holds the best performance of enhancing spatial resolution for more accurate image details of the absorption of graphite powder distribution.
(a) Image of samples
(b) Conventional-50
(c) PSO-SVM-50
(d) Conventional-100
(e) PSO-SVM-100
4. Discussion
In practical photoacoustic imaging, the measurement surface is often incomplete or irregular. When a discrete measured surface is used, the time reversal algorithm fixes the acoustic pressure at these incomplete points which can make them act as optic point scatters, which may scatter back into the imaging region. This can result in arc-like artifacts across the image. One way to reduce the artifacts is to interpolate the measured data to a complete measurement surface acting as the boundary in time reverse course. The PSO optimized SVM interpolation method is used in the time reversal reconstruction algorithm. Compared with other interpolation methods, the method has higher convergence rate and optimization precision. In the course of interpolation, the optimized method can effectively eliminate the phenomenon of contour jaggies and blurring. With better parameter selection, training course, and more measurements, the method can keep more original image details and remove the artifact phenomenon better. In a word, after time reversal reconstruction, both the image magnitude and resolution are improved, so the proposed method can accurately compensate for the effects of acoustic absorption in incomplete photoacoustic imaging.
5. Conclusions
The primary contribution of this paper is a time reversal algorithm based on PSO optimized SVM interpolation which is proposed to produce more accurate reconstructed PAT images. The proposed method will reduce background interpolation artifacts and blurring in the PAT image at the expense of computational speed. What is more, the algorithm is capable of correcting the attention effect. The effectiveness of the algorithm is verified with the simulation. Also the PSO optimized SVM interpolation based time reversal algorithm can be used for medical PAT imaging system to obtain high resolution.
Conflict of Interests
The authors declared that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This research is supported by the National Natural Science Foundation of China (Grants no. 61201307 and no. 61371045) and the Fundamental Research Funds for the Central Universities (Grant no. HIT. NSRIF. 2013132).