Research  Open  Published:
A photoacoustic image reconstruction method using total variation and nonconvex optimization
BioMedical Engineering OnLinevolume 13, Article number: 117 (2014)
Abstract
Background
In photoacoustic imaging (PAI), the reduction of scanning time is a major concern for PAI in practice. A popular strategy is to reconstruct the image from the sparseview sampling data. However, the insufficient data leads to reconstruction quality deteriorating. Therefore, it is very important to enhance the quality of the sparseview reconstructed images.
Method
In this paper, we proposed a joint total variation and L _{ p }norm (TVL _{ p }) based image reconstruction algorithm for PAI. In this algorithm, the reconstructed image is updated by calculating its total variation value and L _{ p }norm value. Along with the iteration, an operatorsplitting framework is utilized to reduce the computational cost and the BarzilaiBorwein step size selection method is adopted to obtain the faster convergence.
Results and conclusion
Through the numerical simulation, the proposed algorithm is validated and compared with other widely used PAI reconstruction algorithms. It is revealed in the simulation result that the proposed algorithm may be more accurate than the other algorithms. Moreover, the computational cost, the convergence, the robustness to noises and the tunable parameters of the algorithm are all discussed respectively. We also implement the TVL _{ p } algorithm in the invitro experiments to verify its performance in practice. Through the numerical simulations and invitro experiments, it is demonstrated that the proposed algorithm enhances the quality of the reconstructed images with faster calculation speed and convergence.
Introduction
Photoacoustic imaging (PAI), also known as optoacoustic tomography (OAT) or thermoacoustic tomography (TAT), is a novel hybrid biomedical imaging modality which combines the strengths of both optical and ultrasound imaging [1–6]. Due to its nonionizing nature, it has been considered as a promising imaging technique and developed rapidly during the past decade. PAI reveals physiologically specific optical absorption contrast of the biological tissues, which has great potential in clinic applications such as early tumor detection [7, 8], vessel imaging [9, 10] and brain imaging [11].
PAI is developed based on the photoacoustic effect [1, 2], which is a process describing that the imaging tissues absorb the laser energy and convert it into acoustic waves. In this paper, we focus on the computedtomographic PAI. In this kind of imaging mode, a laser pulse is used to illuminate the imaging tissues from the top. Some of the laser energy is absorbed and converted into heat, leading to thermoelastic expansion and thus wideband ultrasonic wave emission. The generated photoacoustic signals are then detected by a scanning ultrasound transducer or a transducer array to form images. Based on these detected signals, the optical absorption deposition of the imaging tissues can be calculated by using an image reconstruction algorithm.
In PAI, the reconstruction algorithms have become the vital factor of imaging quality. The PAI reconstruction result benefits a lot from a stable, accurate and efficient algorithm. A variety of analytical image reconstruction algorithms have been developed. Reconstruction algorithms based on the inverse spherical radon transform have been proposed in both the timedomain [12, 13] and the frequencydomain [14, 15]. The filtered backprojection (FBP) algorithm proposed by Xu et al. is the most popular algorithm due to its accuracy and convenience [16, 17]. The deconvolution reconstruction algorithm proposed by Zhang et al. has specific advantage under the circumstance of limitedangle sampling and heterogeneous acoustic medium [18, 19]. Several investigations have been made to propose the algorithms in plane geometries for imaging with the linear array of transducer [20, 21]. The analytical reconstruction methods have advantage in the computational cost and implementation convenience. However, the analytical algorithms fail to keep effective when the sampling points are sparse. The ignorance of measurements noises leads to the severely quality decline in the noisy situation. Those drawbacks above limit the applications of the analytical algorithms and impair their performance. Then the iterative image reconstruction methods are proposed to overcome these shortcomings and enhance the image quality of PAI.
The iterative image reconstruction methods usually build up a model to describe the relationship between the detected photoacoustic signals and the optical absorption deposition. So they are also called the modelbased algorithms. Most of them calculate the optical absorption deposition iteratively to get the final reconstructed image. With proper optimization condition setup, the modelbased methods can provide a more accurate and robust image reconstruction compared to the analytical ones [22, 23]. Many methods that proved to be useful in other aspects have been adopted in PAI reconstruction as an optimization condition of the modelbased methods. Some algorithms focus on the compensation of the acoustic inhomogeneous phenomenon [24, 25]. Jose et al. proposes an iterative approach that takes the speedofsound of subject into account. They acquire the 2D speedofsound distributions and use this speedofsound map in their reconstruction algorithm [24]. Huang et al. develop and establish a fullwave iterative reconstruction approach in PAI to deal with the acoustic inhomogeneous and acoustic attenuation problem [25]. The compressed sensing has been involved in PAI reconstruction aiming to reduce the measurements and accelerate the data acquisition [26, 27]. The modelbased algorithm proposed by Rosenthal et al. recovers the image in the wavelet domain with a different strategy [28]. Meng et al. develop a compressed sensing framework by using partially known support [29, 30]. The reported results show some improvement of image qualities. The total variation (TV) coefficient is always used to denoise the image. Some algorithm are proposed by using the total variation minization to PAI image reconstruction [31–34]. Yao et al. propose the total variation minization (TVM) with the TV coefficient involved in the finite elment method to enhance the image quality and overcome the limitangle problem [31, 32]. An adaptive steepestdescentprojection onto convex sets (ASDPOCS) is proposed by Wang et al. with the TV utilized in the iteration [33]. They investigate and employ the TVbased iterative image reconstruction algorithms in threedimensional PAI. Zhang et al. utilizes the TV coefficient along with the gradient descent method in PAI reconstruction to propose the Total variation based gradient descent (TVGD) algorithm [34]. The TVGD method is reported to be stable and efficient under the sparseview circumstance for PAI reconstruction. From the discussion noted above, it can be deduced that the iterative algorithms for PAI reconstruction have advantage in reconstruction qualities and robustness to noise. Now the reduction of scanning time is the main concern of PAI. A popular strategy is to reconstruct the image from sparseview sampling data. Also there exists photoacoustic imaging system which can image the whole area with one laser exposure. These systems usually have large amount of transducers around the imaging area. With the help of sparseview photoacoustic imaging reconstruction method, the transducer amount can be reduced. This reduction benefits the system from two main aspects. First, it helps to manage the system complexity to a lower level. The lower complexity system is more stable and easier to maintain. Second, this reduction also means the reduction of data scale. The data scale reduction can make the acquisition process more simple and flexible. Besides these two aspects, it is also worth to mention that it reduces the cost of the whole system. All those mentioned above is very important for further clinical applications. It is very important to develop a sparseview imaging system. In this situation, the qualities of the iterative reconstructed images have room for improvement. Take the TVGD method for example, it is reported to be an efficient and highquality algorithm in sparseview situation. But the paintinglike artifacts emerge and some detail information is lost in the extremely sparseview reconstruction.
The compressed sensing theories have been adopted in PAI reconstruction, in which the L _{1}norm of the signal is minimized to obtain the reconstructed image. Recently, Chartrand reported that by replacing the L _{1}norm with the L _{ p }norm (0 < p ≤ 1), accurate reconstruction is possible with substantially fewer measurements [35]. This nonconvex optimization setting has been successfully applied to Magnetic Resonance Imaging (MRI) image reconstruction [36, 37]. The results show that the algorithm with L _{ p }norm can provide accurate reconstruction image with fewer measurements comparing to the L _{1}norm based algorithms. To another dimension of this optimization problem, several algorithms have been proposed to get a better performance in image reconstruction through jointly minimizes the TV value and L _{1}norm value [38, 39] in MRI image reconstruction.
In this paper, we present a novel algorithm to the problem of reconstructing the image from sparseview data in PAI. The algorithm is based on the jointly minimization of total variation and nonconvex L _{ p }norm (TVL _{ p }). The reconstructed image is updated by calculating its joint total variation value and L _{ p }norm value. The operatorsplitting framework is used to reduce the computational cost, and the BarzilaiBorwein step size selection method is adopted to obtain the faster convergence. Through the numerical simulation, the image reconstruction in the case of insufficient sampling data was accomplished. The reconstruction result is compared with several other algorithms including the FBP [16], the L _{1}norm [27] and the TVGD method [34]. The computational cost and the convergence of the proposed algorithm is also discussed and compared with other algorithms. The numerical simulations also cover the robustness to the noise and the tunable parameter investigation. Through the numerical simulations and invitro experiments, it is demonstrated that the proposed algorithm enhances the quality of the reconstructed images with the faster calculation speed and convergence. It’s worthwhile to mention that like Ref. [27] and other iteration method, we also used a projection matrix to connect the acoustic pressure measurements with the reconstructed image. But there are some implementation differences between our method and that one. We both use an intermediate variable to simplify our equations. Ref. [27] used the velocity potential as the intermediate variable and we used a linear integration of the initial pressure along an arc whose center is the position of the ultrasound sensor and with a certain radius ct. The Ref. [27] used a sparsifying matrix and minimized the L _{1}norm in sparsifying domain to get the reconstruction. We used the information from sparsifying domain and piecewise continuous behavior to reconstruct the image. Also, we adapted the pnorm minimization into the algorithm, so it can be a more accurate algorithm in sparseview PAI.
The main contribution of this paper is to develop a novel algorithm for solving the problem of reconstructing the image from sparseview data in PAI. Our contributions are threefold. First, we include the nonconvex optimization into the PAI reconstruction. This nonconvex optimization setting can provide more stable and accurate result under the sparseview situation. Second, we combine the nonconvex optimization with TV minimization. The combined method is able to reconstruct more detailed image with sharp edge. Finally, we implement the BarzilaiBorwein method accelerates the reconstruction speed and improves the convergence considerably.
This paper is organized as follows. ‘Theory and method’ describes the theory of the proposed algorithm. The numerical simulation is introduced in ‘Simulation’. The invitro experimental results are shown in ‘Invitro experiments’. The conclusions of this work are drawn in ‘Conclusion’.
Theory and method
Photoacoustic theory
In this paper, the twodimensional PAI is concerned in the simulations and experiments. In 2D PAI, a laser pulse is used to illuminate the imaging tissues from the top. Due to the photoacoustic, the illumination creates an initial acoustical pressure field. The initial acoustical pressure field propagates as ultrasound waves, which can be detected by ultrasound transducers. Based on the physical principle of the photoacoustic effect, assuming that the illumination is spatially uniform, the relationship between the acoustical pressure measurements and the initial pressure rise distribution can be derived as :
where $p\left(\overrightarrow{r},t\right)$ is the acoustic pressure measurements at the position r and the time t, c is the sound speed, C _{ p } is the specific heat, μ is the isobaric expansion coefficient, I(t) is the temporal profile of the laser pulse and $u\left(\overrightarrow{r}\right)$ is the initial pressure rise distribution. In our study and many photoacoustic tomography studies, we employ a laser pulse with a very short duration. Its duration is nano seconds. So here we made an approximation to treat the I(t) as a Diracdelta function.
In order to recover initial data for the wave equation, several inversion formulas have been established to solve this as a filtered backprojection problem [12, 40]. By using the Green’s function [12] to solve equation (1), the acoustic pressure measurements can be deduced as:
where ${\overrightarrow{r}}_{0}$ is the position of the ultrasound transducer.
In PAI experiments, an ultrasound transducer is used to receive the acoustic pressure measurements at different positions, and the image reconstruction is regarded as an inverse problem to obtain the initial pressure rise distribution. In the iteration of the image reconstruction, a projection matrix A is typically established to connect the acoustic pressure measurements with the reconstructed image. The measurements can be calculated based on the reconstructed image, and then the reconstructed image can be repeatedly corrected by minimizing the difference between the calculated measurements and the real ones. In this way, the optimization method can be used for collaboration and then the iteration reconstruction algorithm can be developed.
Compressed sensing for PAI
If the sampling data is insufficient, the projection matrix A is illconditioned. Thus, the matrix A does not have an exact inversion. As a result, it leads to streaking artifacts in the reconstructed image. This problem can be treated by incorporating the compressed sensing theory into PAI.
We define a new variable f as:
Then the equation (2) can be converted as follows:
In practical imaging, the reconstructed image and the measurements are processed discretely, and the image is reshaped into vectors for convenience. If the size of the reconstructed image $u\left(\overrightarrow{r}\right)$ is X pixels × Y pixels, then the total pixel number of the reconstructed image $u\left(\overrightarrow{r}\right)$ is N (N = XY). After vectorization, the reconstructed image $u\left(\overrightarrow{r}\right)$ becomes a vector u with the length of N. If the total number of the detection points is Q, the length of measurement in each detection point is M, the equation (4) can be expressed as:
where f _{ i } is the integration of the $u\left(\overrightarrow{r}\right)$ along the arc that is centered in i th detection point and with a radius of ct, A _{ i } is the projection matrix of the i th detection point, T is the transpose operation of a matrix. The calculation of the projection matrix is as follows:

(a)
Calculate an matrix A _{ i }(j) as:
where $d=\sqrt{{\left(m{m}_{i}\right)}^{2}+{\left(n{n}_{i}\right)}^{2}}$ , (m, n) is the position of the j th point in the reconstructed image,(m _{ i }, n _{ i }) is the position of the i th detection point,dx is the actual length between the two pixels in the reconstructed image,dt is the discretized time step and M is the total sampling points at one detection point.

(b)
Vectorize the matrix A _{ i }(j) as the j th column vector in projection matrix A _{ i }.

(c)
Repeat the calculation M times to get the projection matrix A _{ i }.

(d)
Repeat step(a) to step(c) Q times to get the projection matrix in the different sampling positions (A _{1},A _{2},…,A _{ Q }). Then write the projection matrixes in the forms as follows:
The equation (6) can be expressed as:
where the sizes of f, A and u are MQ pixels × 1 pixel, MQ pixels × N pixels and N pixels × 1 pixel respectively.
To reconstruct the photoacoustic image from incomplete measurements by using the compressed sensing theory, we can solve an optimization problem as follows:
where Ψ is a sparse transform matrix, _{1} and  _{2} are the L _{1}norm and L_{2}norm respectively. By projecting the image onto an appropriate basis set, we can get a sparse representation of the original image. In this domain, most coefficients of the image are small, and a few large coefficients capture most information of the signal. In this way, we can recover a much more accurate image from those undersampled measurements.
In practical applications of PAI, the reconstructed images often show piecewise continuous behavior. The images like this always have small total variation (TV) values, which is defined as follows:
where D _{ i } is a matrix with the size of 2 pixels × N pixels that has two nonzero entries in each row to calculate the finite difference of u at the i th pixel. D is a matrix with the size of 2N pixels × N pixels, and D = (D^{X};D^{Y}),D^{X} and D^{Y} are the horizontal and vertical global finite difference matrixes respectively.
It is reported that the TV based reconstruction algorithm can recover the image accurately from sparse sampling data [34]. Using TV values to reconstruct the image can be expressed mathematically as:
However, the TV minimization still has some limitations that impair its performance. The optimization of the TV value encourages the recovery of images with sparse gradients, thus resulting in the paintinglike staircase artifacts in the reconstructed images.
Recently, some research find out that the nonconvex optimization can reconstruct an accurate image with fewer measurements by replacing the L _{1}norm with the L _{ p }norm (0 < p ≤ 1). Aiming to enhance the reconstruction quality and overcome the problem of TV based algorithm, we joint the L _{ p }norm with TV values to establish a new optimization which can be defined by:
where α and β are parameters corresponding to the weights of the TV value and L _{ p }norm value ,  _{ p } is the L _{ p } – norm in this optimization problem respectively.
Therefore, we can obtain the reconstructed image by solving this new optimization problem in equation (12).
PAI reconstruction algorithm
In this part, we solve the optimization problem in equation (12) to establish a novel photoacoustic image reconstruction algorithm by using the total variation and nonconvex optimization.
We define the finite difference approximations to partial derivatives of u at the i th pixel along the coordinate as variable ω _{ i } = D _{ i }u, the i th pixel’s sparse coefficient as variable z_{ i } = Ψ _{ i } ^{T}u, where Ψ_{ i } is the sparse transform matrix of the i th pixel. The equation (12) can be deduced as:
where ρ is the parameter corresponding to the weight of the constraint condition in this optimization problem.
We form the augemented Lagrangian defined by
where b _{ i } ^{k} is the TV step parameter in k th iteration, c _{ i } ^{k} is the L _{ p }norm step parameter in k th iteration, u^{k} is the vectorized image reconstruction in k th iteration.
This problem can be solved by
where ω^{k+1} is the finite difference approximations to partial derivatives of u in (k + 1)th iteration, z^{k+1} is the sparse coefficient in (k + 1)th iteration, u^{k+1} is the vectorized image reconstruction in (k + 1)th iteration, z _{ i } ^{k+1} is the sparse coefficient of the i th pixel in (k + 1)th iteration, ω_{ i } ^{k+1} is the finite difference approximations to partial derivatives of u at the i th pixel along the coordinate in (k + 1)th iteration; b _{ i } ^{k+1} is the TV step parameter in (k + 1)th iteration, c _{ i } ^{k+1} is the L _{ p }norm step parameter in (k + 1)th iteration.
By using the standard augmented Lagrangian method, the optimization problem in (15) can be deduced as
where ω^{k} is the finite difference approximations to partial derivatives of u in k th iteration, z^{k} is the sparse coefficient in k th iteration, u^{k} is the vectorized image reconstruction in k th iteration; δ_{ k } is the BarzilaiBorwein step parameter in k th iteration.
After using the BarzilaiBorwein method to determine the step size δ, the optimization problem in equation (13) can be transformed into three subproblem as follows:
where ω_{ i } ^{k} is the finite difference approximations to partial derivatives of u at the i th pixel along the coordinate in k th iteration respectively, z _{ i } ^{k} is the sparse coefficient of the i th pixel in k th iteration respectively; δ_{ k+1} is the BarzilaiBorwein step parameter in (k + 1)th iteration.
We use the soft shrinkage operator to obtain the solution to ωsubproblem in equation (17), the operation is as follows:
where a _{1}, a _{2}, t _{1} and t _{2} are the variables used for a succinct expression.
As for the zsubproblem in equation (17), we use the soft pshrinkage operator to solve it. The operator is defined by:
where a _{3}, a _{4}, t _{3} and t _{4} are the variables used for a succinct expression.
The usubproblem in equation (17) is a typical least squares problem. The solution can be easily obtained by:
where F is the Fourier transform matrix.
As a result, the TV L _{ p } algorithm is summarized as follows:

(1)
Initialization: input f, α, β, ϵ ,p and ρ. Set the reconstructed image u^{0} = 0, b = c = 0, δ_{0} = 1, k = 0.

(2)
Apply equation (18) and (19) to update the value of ω and z.

(3)
Apply equation (20) to update the value of u.

(4)
Apply equation (17) to update the value of b, c and δ.

(5)
If the exiting condition is met, end the iterations and output the result. Otherwise repeat the step from (2) to (4). The exiting condition is as follows:
Simulation
To verify the effectiveness of the proposed TV L _{ p } algorithm on PAI reconstructions, the simulations are designed. All the simulations are performed in Matlab v7.14 on a PC with a 3.07 GHz Intel Xeon processor (only 1 core is used in computation) and 32 GB memory. The sparsisfying operator Ψ is set to Haar wavelet transform using Rice wavelet toolbox. The sound speed is set to be consistent in the simulation as 1500 m/s.
Sparseview reconstruction
In the simulation, we choose the SheppLogan phantom to be the initial pressure rise distribution. The forward simulation and inverse reconstruction are all performed in 2D. The phantom is shown in Figure 1. The measurements from the phantom are generated by using equation (2). The size of the phantom is 89.6 mm × 89.6 mm, the radius of the scanning circle is 42 mm and the size of the reconstructed image is 128 pixels × 128 pixels. During the simulation, the scanning circle covers 360°around the imaging phantom. Four different measurements are collected. The scanning step of tomographic angels is set to 2.25°, 4°, 12° and 20° respectively. So the sampling points are 160 views, 90 views, 30 views and 18 views correspondingly.
The parameters α, β, ϵ and ρ are set to be 1 × 10^{2}, 1 × 10^{2}, 1 × 10^{5} and 1 respectively. The influence of these parameters will be discussed later. And the parameter p are set to be two different values as 0.5 and 0.8.
We choose the FBP [16], the L _{1}norm [27] and the TVGD [34] algorithms to be the comparison besides our proposed TV L _{ p } algorithm. The simulation results by using these different algorithms are shown in Figure 2. It’s worthwhile to note that the weight used in the TVGD algorithm is an adaptive parameter, as same as it is reported in [27]. The negative values in FBP reconstructed image are set to be zero.
It is shown in the first column of Figure 2 that, all three iterative algorithms have comparable reconstruction results when the sampling data is sufficient. Moreover, it is shown in Figure 2(a) that the contrast of FBP reconstructed image is not as high as the other three. But its resolution is comparable with the others visually. When the number of sampling points reduces, the qualities of the reconstructed images are strongly affected in the FBP reconstruction. When the sampling point gets sparse in the FBP reconstruction, the arclike artifacts appears due to the backprojection arcs cannot be canceled out with each other. The iterative algorithms can provide better qualities of the reconstructed images than the FBP method in sparseview reconstructions. Among them, the L _{1}norm method struggles to depress the noise. Meanwhile, the TVGD algorithm and the TVL _{ p } algorithm provides highresolution images and have no visually distinguishable decline in qualities of the reconstructed images as the number of sampling points decreases.
As for the extreme sparse sampling points situation (18view and 30view), the image reconstructed by the FBP algorithm, shown in Figure 2(d), has extremely severe artifacts. The L _{1}norm reconstruction and the TVGD algorithm have a decline in image qualities. The noise in the reconstructed images by the L _{1}norm reconstruction, as shown in Figure 2(j), cannot be depressed effectively. As for the TVGD algorithm, the reconstruction produces piecewise artifacts which also make the qualities of the reconstructed images decrease. In the TVL _{ p } image, the noise is depressed more effectively than the above three algorithm. The quality of the reconstructed image is not affected substantially by the insufficient sampling data.
We calculate the peak signaltonoise ratios (PSNR) of the reconstructed images with the original phantom as a gold standard to provide a numeric quantification of the results. The bigger the value of PSNR is, the better quality of the image is. The PSNR is defined as:
where t(i,j) means the grayvalue of the original image, MAX the maximum possible pixel value of the image which in our simulation is 1.
We calculate the value of the PSNR of all images in Figure 2. The quantitative results are shown in Table 1.
From Table 1, it is shown that the PSNR of the FBP algorithm is always in a very low level due to its unsuitability for sparseview sampling condition. As for those three compressed sensing based algorithms, the PSNR value of images reconstructed by the TVL _{ p } algorithm are the highest. The L _{ p }  norm optimization constraint can provide the better performance in the extremely sparse sampling. With this improvement, the TVL _{ p } algorithm is more accurate than other algorithms in the sparseview sampling condition shown in the quantitative results. Between the two different value of p, the parameter p that is set to be 0.5 has a slightly advantage against the other one. Also it is revealed from Table.1 that the 90view shows higher PSNR than 160view. When the sampling points are sufficient, it is possible that the fewerview projection can produce better reconstruction results. But their PSNR is very close with same algorithm. It is fair to say that the results are on the same level of image quality.
From Figure 2, it is shown that the TVGD image and the TVL _{ p } image are very close in the image quality. Here we choose the FORBILD phantom [41], a more complicated and more challenging phantom, to further compare the proposed algorithm with the TVGD algorithms. The phantom is shown in Figure 3. The scanning step of tomographic angels is set to 2.25°, 4°, 6° and 12°. So the sampling points are 160 views, 90 views, 60 views and 30 views correspondingly. The other numerical implementation conditions remain the same with the sheplogan simulation. The simulation results by using the TVGD algorithms and the proposed TV L _{ p } algorithm are shown in Figure 4. It is shown in Figure 4 that, when the sampling data is sufficient, both algorithms can reconstruct the accurate image. When the sampling angles get sparse during the simulation, it is seen that the TVGD reconstruction results in paintinglike staircase artifacts in the smooth regions. Also it fails to give the accurate image in the low contrast regions in the top and left of the phantom. The proposed algorithm provides reasonably good reconstructions in these regions. The PSNR of the reconstructed images are shown in Table 2. From this table, we observe that the TVL _{ p } algorithm provides the better PSNR for all the cases. In the case of more complicated phantom, the TVL _{ p } algorithm shows significant improvement than the TVGD algorithm.
Also, we include a lineplots image of the reconstruction result by the TVGD algorithm and the TVL _{ p }(p = 0.8) algorithm from 30view data. The location of the pixel profile in the image is displayed in Figure 5(a). The comparisons of pixel profiles are displayed in the Figure 5(b).
In Figure 5(b), the solid line and the dotted line represent the pixel profiles of the TVL _{ p } and the TVGD image respectively. It is shown in Figure 5(b) that the TVL _{ p } can reconstruct the image more precisely than the TVGD one. The edges from the TVL _{ p } are sharper than that from the TVGD. The pixel number from 90 to 100 is the high resolution area, the TVL _{ p } image shows the highspeed change of the pixel value while the TVGD fails to do so. In the continuous area, the TVL _{ p } image is smoother.
We continuously decrease the number of the detect points try to find the limit density of the sampling points. During the simulation, we set the criterion of acceptability that the PSNR of the reconstructed image reaches 30 dB. It is found out that the total number of the sampling points is able to be reduced to 15 for TVL _{ p } algorithm in the reconstruction of sheplogan phantom and to 18 for the forbild phantom.
In this part, the TVL _{ p } algorithm is proved to be more accurate and stable than the other algorithms for PAI image reconstruction in the sparse sampling condition.
Convergence and calculation
In this part, we discuss the theoretical calculation complexity and study the convergence of the proposed algorithm. As mentioned above in ‘Theory and method’, in step (2) the update of ω and z is using the soft shrinkages and the computational costs are both O(N). The update of z also includes a wavelet transforms which computational costs are O(N logN). In step (3), the update of u involves two fast Fourier transforms which computational costs are O(N logN) and two operations of A with the computational costs of O(NMQ). The update of the parameters b and c in step (4) are all simple calculations with the computational costs of O(N). As for the paremeter δ, although it involves an operation of A, it can use the result computed in the step (3). So its computational cost is also O(N).
In a nutshell, the calculation complexity of the proposed algorithm in one iteration is 5O(N) + 4O(N logN) + 2O(NMQ). The first two terms is much smaller than the last term in the practical use of photoacoustic imaging and most iterative algorithm is with the operation. In each iteration, we just use the projection matrix twice. So the proposed algorithm has a cheap per iteration computation.
The TVGD algorithm is reported as an efficient and stable iterative algorithm in photoacoustic imaging. In the ‘Sparseview reconstruction’, its reconstruction result is closest to the proposed algorithm. So here we select it to be a comparison with the TVL _{ p } algorithm. We calculate the time cost of those two algorithms in a simulation. The simulation condition is same as in ‘Sparseview reconstruction’. But the iteration ends when the PSNR values reach 30 dB. The result is shown in Table 3. From Table 3, it is shown that the proposed algorithm is faster than TVGD algorithm in the computational time. Based on this result, it could be inferred that the TVL _{ p } algorithm is a more efficient image reconstruction algorithm comparing to the TVGD algorithm. The value of p has also some influence on the time cost. The smaller the p is, the more iteration times are needed to reach the reconstruction result.
Thanks to the use of BarzilaiBorwein step size selection method, the convergence speed can also be significantly improved. For the quantitative analysis, we use a parameter that represents the distance between the reconstructed image and the original phantom image. The parameter d is defined as:
where u is the reconstructed image and t is the original image. The size of the image is X × Y. The smaller the parameter d is, the closer the reconstructed image is with the original phantom. In the TVL _{ p } algorithm, there is a small rate of chance that the optimization will lead to the wrong solution due to its nonconvex nature. So we use the original image to calculate the parameter d to show the image quality and use the parameter d as a reference. We want to show the improvement of the image quality in every iteration step.
The simulation condition is set to be the same as in ‘Sparseview reconstruction’. The sampling view is 60. The parameter p is set to be 0.8. The defined distance d is calculated after each iteration step. If the distance is smaller than 0.05, the iteration will stop. The simulation result is shown in Figure 6. The xaxis is the value of distance and the yaxis is the iteration times. The line ‘·’ refers to the TVGD algorithm and the line ‘*’ represent the TVL _{ p } (p = 0.8)algorithm. The result is shown in Figure 6. The images reconstructed by TVL _{ p } algorithm in each iteration have smaller value of d than the TVGD ones and the TVL _{ p } iteration only takes 9 times as the distance is met the request.
For discussions noted above, it can be surmised that the convergence of TVL _{ p } algorithm is faster and the TVL _{ p } has a cheaper computational cost.
Robustness to the noise
In the practical applications of the photoacoustic tomography, the measurements are usually polluted by those white measuring noises from the ultrasound transducer and the system electronics. Hence, it is very important for an algorithm to maintain stable performance under the noises polluted circumstances. To analyze the robustness of the TVL _{ p } algorithm, we choose 30view simulated photoacoustic signals that we used in ‘Sparseview reconstruction’. The signals are added with white noises of different noise power levels. We use TVL _{ p } algorithm with two different settings of parameter p (p = 0.5 and p = 0.8) and TVGD algorithm to reconstruct images from these white noise polluted measurements.The reconstruction results are shown in Figure 7. In the first row to the last row of the Figure 7, the signal to noise ratio (SNR) of the polluted measurements is 10 dB, 5 dB, 3 dB and 0 dB, respectively. As shown in the image, when the power level of the noises is not very strong (10 dB and 5 dB), the reconstructed images by using the noisy measurements have basically no obviously difference with the ones reconstructed with the noiseless signals. As the noise becomes stronger, the quality of the reconstructed images decreases.
We also plot the profiles of a pixel line in order to show the detail qualities of reconstructed images clearly. In Figure 7, the dotted line and the solid line are the pixel profiles of the reconstructed image and the original image, respectively. From the line plots, it is revealed that the proposed algorithm has better performance in the edge preservation and more accurate in the smooth area. We calculated the PSNR of the reconstructed image. The result is shown in Table 4. Our algorithm outperformed the TVGD algorithm in any noise power level. Giving the credit to the optimal conditions, the reconstructed image is intended to be continuous and sharp. During the iteration, the photoacoustic signals get enhanced and the noise is suppressed.
As we can see from the table that the TVL _{ p } algorithm reconstructed images have a slightly better performance than the TVGD algorithm when it comes to the image qualities. When the noise is extremely strong (0dB), the TVL _{ p } algorithm has a huge advantage than the TVGD one in image quality. As for the different settings of parameter p in the TVL _{ p } algorithm, there is no major difference that can be observed between the two images. The PSNRs of the p = 0.5 setting is about 0.3 dB bigger than the p = 0.8 setting in the first three noise power level. But when the noise getting extremely strong (0 dB), the p = 0.8 setting is 0.1 dB bigger than the p = 0.5 setting in PSNR value.
From this part of simulation we can conclude that our TVL _{ p } algorithm is robust to noise and has a better performance than the TVGD algorithm in the noisy measurement circumstance.
Parameter investigation
As the original optimization problem of the image reconstruction is described in Eq. (14), the TVL _{ p } algorithm contains some parameters that are tunable, which are α, β, ϵ, p and ρ. In those parameters, the choice of ρ does not affect the performance of the TVL _{ p } algorithm theoretically. The result of the simulation also shows that the image quality is not sensitive to the parameter ρ for a large range. Here we set the parameter ρ to a steady value 1. The ϵ is the exiting condition parameter. It can be easily deduced that smaller ϵ will leads to slightly more accurate reconstructed image at the cost of more iteration times. In this part we focus on analyzing the parameter settings of p, α and β.
Parameter setting of p
In the TVL _{ p } algorithm, we replace the L _{1} norm with the L _{ p } norm (0 < p ≤ 1). It is reported in Ref. [35] that theoretically fewer measurements are required for accurate reconstruction in the L _{ p } norm situation. But it also leads to failure in solving the optimization problem. It’s kind of a dilemma for the setting of p. So here we take different values of p to see its influence to the image reconstruction. The parameter α and β are both set to be 1 × 10^{2}.
We choose the 90view and 18view simulated photoacoustic signals that we used in ‘Sparseview reconstruction’. We set the p value as 0.3, 0.5, 0.8 and 1. We calculate the reconstructed images’ PSNR value. It is shown in Table 5. When the p is set as 0.5, it has advantage in the quality of the reconstructed image. However, when the value of p continues to reduce to 0.3, there is no obvious improvement in image quality. But in the same time, the smaller the p is, the higher probability of the solving failure is during the simulation. The reduction of p leads to the increasing of the iteration times in our simulation. So taking these two factors into account, we set the p as 0.8 so that it can provide a great reconstruction performance and stability with a fast convergence.
Parameter settings of α and β
As we describe above in Eq. (14), the parameter α and β are parameters corresponding to the weights of TV value and L _{ p } norm value in this optimization problem respectively. We use these two parameters to balance the terms of the objective function. With different kinds of the objective image, the settings of those two parameters are different. Here we select three different images as the given optical energy deposition to test the universality of our algorithm and present further investigation of the parameter settings. We select a phantom that stand for the vessels and a phantom of dots with different energy degree in the simulation. We also choose a real brain MRI as the original optical energy deposition to demonstrate the performance of the proposed algorithm in reconstructing extremely detailed and complex structured imaging object. Here the TVGD algorithm is used as a comparison. There are four groups of parameter settings, which are (α = 1 × 10^{2}, β = 5 × 10^{3}) , (α = 1 × 10^{2}, β = 1 × 10^{2}), (α = 5 × 10^{2}, β = 5 × 10^{3}) and (α = 5 × 10^{3}, β = 5 × 10^{3}). The reconstructed images are shown in Figure 8. From first row of Figure 6, we can see in the reconstruction of gradient sparse phantom. The TV based algorithm has great performance when the image demonstrate piecewise continuous behavior. All reconstructions are accurate and the background noise is suppressed well. When it comes to the images with the vessel phantom (Figure 8 (g)(l)), those original optical energy depositions are a little bit more complex than the dots. The reconstruction results show that the image reconstructed by TVL _{ p } algorithm is better than the TVGD ones. As in Figure 8(h), TVGD image has some noises in the background and the edge of the vessel is blurred. While the TVL _{ p } images with different parameter settings both have highresolution results. As for the real MRI image, it has very detailed information. As expected, both two groups of parameter setting α = β have the most accurate result among them. The increasing weighting of L _{ p }norm condition can provide more detail information and prevent the reconstructed image emerging plantlike artifacts. The details such as edges and fine structures are well preserved in both reconstructions. The reconstruction results show that the TVGD reconstructed image has severely paintinglike staircase artifacts with some loss in fine details. From our observation, TVL _{ p } algorithm with the parameter setting α = β preserves the fine features better than the TVGD one. α and β are the regularization parameters determining the tradeoff between the data consistency and the sparsity. It is revealed from the above simulation that the parameter setting α = β is a better strategy. In this parameter setting, the TVL _{ p } algorithm provides a 3 dB improvement in the PSNR over the TVGD algorithm based on our calculation.
Limitedview and irregularview simulation
In the real application of PAI, due to the restrains of the shape or the size of the imaging object, a full angular scanning sometimes is hard to achieve. We evaluate the performance of the TVL _{ p } method in limitedview case, lineview case and unequalview case.
The simulation setup and the reconstructed image is shown in Figure 9. In the limitedview simulation (Figure 9(a)), the scanning angular range is set to 150° and the angular step is 3°, so 50view photoacoustic signals are obtained. In the lineview simulation (Figure 9(c)), the transducer array with 60 transducers is placed in the right side of the imaging object and the interval between two transducer elements is 1.49 mm. It is revealed in Figure 9(b) and Figure 9(d) that the quality of the TVL _{ p } reconstruction is not much affected by the limitation of the sampling angle. Because the sampling angle is limited, the information definite is partly missed, yet the TVL _{ p } method can still provide a satisfying reconstruction. In unequal angel step scanning, we randomly choose 30 sampling points from a 60view projection and use these 30view unequal angle step data for image reconstruction. The result is shown in Figure 9(f). As we can see from the image, the reconstruction result can still maintain a very high quality.
Invitro experiments
Experiment setup
We carry out the experiments on invitro signals to demonstrate the proposed TVL _{ p } algorithm’s performance in the practical application.
The framework of the experiment platform is shown in Figure 10. In this platform, an Nd:YAG laser generator (Continum, Surelite I) is used to emit the laser pulse. The wavelength of the laser is 532 nm. A single laser pulse is generated at the frequency of 10 Hz and last 67 ns. The incident laser pulse is emitted towards the top of the phantom through a concave lens with the diameter of 5 cm. The setup of the lens enlarges the illumination area and lead to the pulse energy reduction in the illumination area. The energy is about 6.47 mJcm^{2}, which is lower than the ANSI laser radiation safety standard (20 mJcm^{2}) [1]. Signal acquisition is done by a waterimmersion ultrasound transducer (Panametric, V383SU). The transducer is a linearly unfocused one at 3.5 MHz (6 dB bandwidth at 45%). A digital stepping motor (GCD0301 M) is used to rotate the ultrasound transducer around the phantom placed in water. The scanning radius is 38 mm. The received analog ultrasound signals are amplified by a pulse receiver (Panametric, 5900PR). An oscilloscope (Agilent, 54622D) with the sampling frequency of 16.67 MHz is set to transform the received signals into digital ones. Both the laser generator and the digital motor are controlled by the computer through the serial interface. The transformed digital data is transported to the computer through the general purpose interface bus (GPIB).The imaged phantom we used in the experiment is made by gelatin cylinder. It is shown in Figure 11. There are two different phantoms. The radius of the phantom is 25 mm. The left one is made by two rubber bars with 1 mm diameter that embedded as the optical absorbers. The right one utilizes leaf which pretends as vein and tissue as the optical absorbers.
In the experiment, the transducer tends to measure the photoacoustic signal inplane only, and the reconstruction is also in 2D. The crosssectional image in any plane is mainly determined by the measured data in the same plane, and a set of circular measurement data on the same plane would be sufficient to reconstruct a good image. We use the deconvlution calculation before the reconstruction to eliminate the transducer’s impulse response influence.
Experiment result
In the experiment, 90view and 30view data are collected for reconstruction. The images are constructed by the FBP, the TVGD and the TVL _{ p } algorithms, respectively. The reconstruction results are shown in Figure 12. The left column of the Figure 12 is reconstructed from 90view data. When the sampling data is sufficient, all three algorithms are effective. With respect to the locations and sizes, the optical absorbers are all well reconstructed in the figure. While the FBP reconstructed image is not as clear as the images reconstructed by the iterative algorithms. When we reconstruct the image with a small of sampling angles (right column of Figure 11), the artifacts start to emerge in the FBP reconstructed image and the quality of the image is severely affected. But the TVGD and the TVL _{ p } algorithms can still provide highcontrast images with less noise. In Figure 12(f), it is shown that the image reconstructed by the TVL _{ p } algorithm outperforms other algorithms in image contrast and noise suppression. The structure of optical absorbers is clear and the noise in the background is wellsuppressed. The sparseview of sampling has barely any influence on the quality of the TVL _{ p } reconstructed image.
In vitro imaging of a leaf vein is also performed to further demonstrate the advantages of the TVL _{ p } algorithm. The reconstruction result is shown in Figure 13. As the structure of the phantom is more complex, the FBP are deeply influenced by the artifact and fail to reconstruct the accurate image both under 90view and 30view sampling circumstance. It is shown in the figure, that the TVGD and TVL _{ p } algorithms can still reconstruct the image in a high contrast level. But when the data is insufficient, there is some noise emerging in the background. TVL _{ p } algorithms can suppress the noise better than the TVGD one. The optical absorber in TVL _{ p } one is more distinct than that in the image by TVGD algorithm.
Quantitative comparisons
We use the L _{1}norm algorithm to reconstruct the image of 180view data from the leaf vein phantom. As the sampling view is efficient, the reconstructed image is used as a “standard” one. We calculate the histograms of the difference between the reconstructed one and the “standard” one as shown in Figure 14. Figure 14 (a)(c) are the difference histograms between the standard and the images reconstructed by the FBP, the TVGD and the TVL _{ p }, respectively, with 30view data. In Figure 14, two CSbased algorithms have a large number of pixels with small ranges of difference with the standard one, which suggests that these two algorithms can reconstruct the image more accurately. In the case of the TVL _{ p } algorithm, the major part of the pixel difference is in the range from 0 to 0.1. From this experiment, the results demonstrate that the TVL _{ p } method can outperform the TVGD one in the field of the image quality.
From the experiment result noted above, it is safe to say that the TVL _{ p } algorithm would have better performance in sparseview PAI than other algorithms. It could provide stable and accurate reconstruction in both sufficient data sampling and sparseview sampling situation.
Conclusion
Aiming to reduce the scanning time and enhance the imaging quality of the photoacoustic image reconstruction, we proposed the TVL _{ p } algorithm that applies the total variation method and nonconvex optimization method to the PAI. The main idea of the algorithm is to apply L _{ p }norm nonconvex optimization along with the total variation method. In the proposed algorithm, the BarzilaiBorwein step size selection method is adopted to provide faster convergence and smaller calculation. The effectiveness and universality of the algorithm is demonstrated through the numerical simulations. The numerical simulations show that the TVL _{ p } algorithm provides good imaging quality in sparseview sampling situation. The algorithm convergence, the robustness to noise and the tunable parameters are also discussed. The simulation result reveals that the TVL _{ p } algorithm is a stable image reconstruction method with fast convergence and small computational cost. The TVL _{ p } algorithm is further investigated through some experiments using gelatinmade phantom. Compared with the result of other popular image reconstruction method, the TVL _{ p } imaging algorithm has significant advantage on contrast and noise suppression. From the discussion noted above, it could be concluded that the TVL _{ p } algorithm may be a practical algorithm for sparseview photoacoustic imaging reconstruction.
Abbreviations
 PAI:

Photoacoustic imaging
 TVL _{ p } :

Total variation and Lpnorm
 OAT:

Optoacoustic tomography
 TAT:

Thermoacoustic tomography
 FBP:

Filtered backprojection
 TV:

Total variation
 TVM:

Total variation minization
 ASDPOCS:

Adaptive steepestdescentprojection onto convex sets
 TVGD:

Total variation based gradient descent
 MRI:

Magnetic resonance imaging.
References
 1.
Wang LV: Tutorial on photoacoustic microscopy and computed tomography. IEEE J Sel Top Quantum Electron 2008,14(1):171–179.
 2.
Xu M, Wang LV: Photoacoustic imaging in biomedicine. Rev Sci Instrum 2006,77(4):041101–1041101–22.
 3.
Li C, Wang LV: Photoacoustic tomography and sensing in biomedicine. Phys Med Biol 2009,54(19):R59R97. 10.1088/00319155/54/19/R01
 4.
Wang LV: Prospects of photoacoustic tomography. Med Phys 2008,35(12):5758–5767. 10.1118/1.3013698
 5.
Kruger RA, Reinecke DR, Kruger GA: Thermoacoustic computed tomography—technical considerations. Med Phys 1999,26(9):1832–1837. 10.1118/1.598688
 6.
Kruger RA, Liu P, Fang Y, Appledorn CR: Photoacoustic ultrasound (PAUS)reconstruction tomography. Med Phys 1995,22(10):1605–1609. 10.1118/1.597429
 7.
Guo B, Li J, Zmuda H, Sheplak M: Multyfrequency microwaveinduced thermal acoustic imaging for breast cancer detection. IEEE Trans Ultrason Ferroelectr Freq Control 2007,54(11):2000–2010.
 8.
Pramanik M, Ku G, Li C, Wang LV: Design and evaluation of a novel breast cancer detection system combining both thermoacoustic (TA) and photoacoustic (PA) tomography. Med Phys 2008,35(6):2218–2223. 10.1118/1.2911157
 9.
Zhang EZ, Laufer JG, Pedley RB, Beard PC: In vivo highresolution 3D photoacoustic imaging of superficial vascular anatomy. Phys Med Biol 2009,54(4):1035–1046. 10.1088/00319155/54/4/014
 10.
Niederhauser JJ, Jaeger M, Lemor R, Weber P, Frenz M: Combined ultrasound and optoacoustic system for realtime highcontrast vascular imaging in vivo . IEEE Trans Med Imaging 2005,24(4):436–440.
 11.
Zerda A, Paulus YM, Teed R, Bodapati S, Dollberg Y, KhuriYakub BT, Blumenkranz BS, Moshfeghi DM, Gambhir SS: Photoacoustic ocular imaging. Opt Lett 2010,35(3):270–272. 10.1364/OL.35.000270
 12.
Xu M, Wang LV: Timedomain reconstruction for thermoacoustic tomography in a spherical geometry. IEEE Trans Med Imaging 2002,21(7):814–822. 10.1109/TMI.2002.801176
 13.
Xu M, Xu Y, Wang LV: Timedomain reconstruction algorithms and numerical simulations for thermoacoustic tomography in various geometries. IEEE Trans Biomed Eng 2003,50(9):1086–1099. 10.1109/TBME.2003.816081
 14.
Xu Y, Feng DZ, Wang LV: Exact frequencydomain reconstruction for thermoacoustic tomography—I: Planar geometry. IEEE Trans Med Imaging 2002,21(7):823–828. 10.1109/TMI.2002.801172
 15.
Xu Y, Xu M, Wang LV: Exact frequencydomain reconstruction for thermoacoustic tomography—II: Cylindrical geometry. IEEE Trans Med Imaging 2002,21(7):829–833. 10.1109/TMI.2002.801171
 16.
Xu M, Wang LV: Pulsedmicrowaveinduced thermoacoustic tomography: Filtered backprojection in a circular measurement configuration. Med Phys 2002,29(8):1661–1669. 10.1118/1.1493778
 17.
Xu M, Wang LV: Universal backprojection algorithm for photoacoustic computed tomography. Phys Rev E 2005,71(1):016706–1016706–7.
 18.
Zhang C, Wang Y: Deconvolution reconstruction of fullview and limitedview photoacoustic tomography: a simulation study. J Opt Soc Am A 2008,25(10):2436–2443. 10.1364/JOSAA.25.002436
 19.
Zhang C, Li C, Wang LV: Fast and robust deconvolutionbased image reconstruction for photoacoustic tomography in circular geometry experimental validation. IEEE Photonics J 2010,2(1):57–66.
 20.
Liao CK, Li ML, Li PC: Optoacoustic imaging with synthetic aperture focusing and coherence weighting. Opt Lett 2004,29(21):2506–2508. 10.1364/OL.29.002506
 21.
Modgil D, Rivière PJ: Implementation and comparison of reconstruction algorithms for 2D optoacoustic tomography using a linear array. Proceedings of SPIE 6856:2008; San Jose 2008, 68561D1–68561D12.
 22.
Zhang J, Anastasio MA, Riviere PJ, Wang LV: Effects of different imaging models on leastsquares image reconstruction accuracy in photoacoustic tomography. IEEE Trans Med Imaging 2009,28(11):1781–1790.
 23.
Paltauf G, Viator JA, Prahl SA, Jacques SL: Iterative reconstruction algorithm for optoacoustic imaging. J Acoust Soc Am 2002,112(4):1536–1544. 10.1121/1.1501898
 24.
Jose J, Willemink R, Steenbergen W, Slump C, van Leeuwen T, Manohar S: Speedofsound compensated photoacoustic tomography for accurate imaging. Med Phys 2012, 39: 7262–7281. 10.1118/1.4764911
 25.
Huang C, Wang K, Nie L, Wang LV, Anastasio M: FullWave Iterative Image Reconstruction in Photoacoustic Tomography With Acoustically Inhomogeneous Media. IEEE Trans Med Imaging 2013,32(6):1097–1110.
 26.
Provost J, Lesage F: The application of compressed sensing for photoacoustic tomography. IEEE Trans Med Imaging 2009,28(4):585–594.
 27.
Guo Z, Li C, Song L, Wang LV: Compressed sensing in photoacoustic tomography in vivo . J Biomed Opt 2010,15(2):021311–1021311–6.
 28.
Rosenthal A, Jetzfellner T, Razansky D, Ntziachristos V: Efficient framework for modelbased tomographic image reconstruction using wavelet packets. IEEE Trans Med Imaging 2012,31(7):1346–1357.
 29.
Meng J, Wang LV, Ying L, Liang D, Song L: Compressedsensing photoacoustic computed tomography in vivo with partially known support. Opt Express 2012,20(15):16510–16523. 10.1364/OE.20.016510
 30.
Meng J, Wang LV, Ying L, Liang D, Song L: In vivo opticalresolution photoacoustic computed tomography with compressed sensing. Opt Lett 2012,37(22):4573–4575. 10.1364/OL.37.004573
 31.
Yao L, Jiang H: Photoacoustic image reconstruction from fewdetector and limitedangle data. Biomed Opt Express 2011,2(9):2649–2654. 10.1364/BOE.2.002649
 32.
Yao L, Jiang H: Enhancing finite elementbased photoacoustic tomography using total variation minimization. Appl Optics 2011,50(25):5031–5041. 10.1364/AO.50.005031
 33.
Wang K, Su R, Oraevsky AA, Anastasio M: Investigation of iterative image reconstruction in threedimensional optoacoustic tomography. Phys Med Biol 2012,57(17):5399–5423. 10.1088/00319155/57/17/5399
 34.
Zhang Y, Wang Y, Zhang C: Total variation based gradient descent algorithm for sparseview photoacoustic image reconstruction. Ultrasonics 2012,52(8):1046–1055. 10.1016/j.ultras.2012.08.012
 35.
Chartrand R: Exact reconstruction of sparse signals via nonconvex minimization. IEEE Trans Signal Proc Let 2007,14(10):707–710.
 36.
Majumdar A, Ward RK: An algorithm for sparse MRI reconstruction by Schatten p norm minization. Magn Reson Imaging 2011,29(3):408–417. 10.1016/j.mri.2010.09.001
 37.
Majumdar A: Improved dynamic MRI reconstruction by exploiting sparsity and rankdeficiency. Magn Reson Imaging 2013,31(5):789–795. 10.1016/j.mri.2012.10.026
 38.
Ma S, Yin W, Zhang Y, Chakraborty A: An efficient algorithm for compressed MR imaging using total variation and wavelets. Proceedings IEEE Conference on Computer Vision Pattern Recongnition: 2008; Anchorage 2008, 1–8.
 39.
Ye X, Chen Y, Huang F: Computational acceleration for MR image reconstruction in partially parallel imaging. IEEE Trans Med Imaging 2011,30(5):1055–1063.
 40.
Finch D, Haltmeier M, Rakesh : Inversion of spherical means and the wave equation in even dimensions. SIAM J Appl Math 2007,68(2):392–412. 10.1137/070682137
 41.
Yu Z, Noo F, Dennerlein F, Wunderlich A, Lauritsch G, Hornegger J: Simulation tools for twodimensional experiments in xray computed tomography using the FORBILD head phantom. Phys Med Biol 2012,57(13):237–252. 10.1088/00319155/57/13/N237
Acknowledgment
This work was supported by the National Natural Science Foundation of China (No. 61271071 and No. 11228411), the National Key Technology R&D Program of China (No. 2012BAI13B02) and Specialized Research Fund for the Doctoral Program of Higher Education of China (No. 20110071110017).
Author information
Additional information
Competing interests
The authors declared that they have no competing interests.
Authors’ contributions
Study concept and design (CZ); drafting of the manuscript (CZ); critical revision of the manuscript for important intellectual content (CZ, YZ and YW); obtained funding (YW); administrative, technical, and material support (CZ and YZ); study supervision (YW). All authors read and approved the final manuscript.
Authors’ original submitted files for images
Rights and permissions
About this article
Received
Accepted
Published
DOI
Keywords
 Photoacoustic imaging
 Image reconstruction
 Total variation
 L _{ p }norm
 Nonconvex optimization