 Research
 Open Access
 Published:
Photoacoustic imaging reconstruction using combined nonlocal patch and totalvariation regularization for straightline scanning
BioMedical Engineering OnLinevolume 17, Article number: 105 (2018)
Abstract
Background
For practical straightline scanning in photoacoustic imaging (PAI), serious artifacts caused by missing data will occur. Traditional total variation (TV)based algorithms fail to obtain satisfactory results, with an oversmoothed and blurred geometric structure. Therefore, it is important to develop a new algorithm to improve the quality of practical straightline reconstructed images.
Methods
In this paper, a combined nonlocal patch and TVbased regularization model for PAI reconstruction is proposed to solve these problems. A modified adaptive nonlocal weight function is adopted to provide more reliable estimations for the similarities between patches. Similar patches are searched for throughout the entire image; thus, this model realizes adaptive search for the neighborhood of the patch. The optimization problem is simplified to a common iterative PAI reconstruction problem.
Results and conclusion
The proposed algorithm is validated by a series of numerical simulations and an in vitro experiment for straightline scanning. The results of patchTV are compared to those of two mainstream TVbased algorithms as well as the iterative algorithm only with patchbased regularization. Moreover, the peak signaltonoise ratio, the noise robustness, and the convergence and calculation speeds are compared and discussed. The results show that the proposed patchTV yields significant improvement over the other three algorithms qualitatively and quantitatively. These simulations and experiment indicate that the patchTV algorithm successfully solves the problems of PAI reconstruction and is highly effective in practical PAI applications.
Background
Photoacoustic imaging (PAI), a novel biomedical imaging technique, combines light and ultrasound to detect absorbed photons ultrasonically through the photoacoustic effect [1,2,3]. Compared with traditional imaging techniques, PAI has many advantages. It obtains high image contrast because the photoacoustic images can reflect the laser absorption distribution in the tissue [1]. It is capable of imaging either thicker tissue or deeper organs with better resolution compared to optical imaging because it receives ultrasound signals [3]. What’s more, PAI is also able to provide noninvasive and functional imaging [4, 5]. Because of these advantages, PAI shows great potential in many biomedical applications such as brain imaging [6, 7], tumor detection [8, 9], vessel imaging [10, 11] and molecular imaging [12, 13].
A laser pulse is usually adopted to irradiate the tissue in computedtomographic PAI, which is the main concern of this paper. The light is absorbed by the tissue, and the ultrasound waves are subsequently excited. This process is called the photoacoustic effect [1]. Then, the photoacoustic signals are detected by a scanning transducer or a transducer array. To reconstruct the photoacoustic image from the detected signals, photoacoustic reconstruction algorithms are required, which directly determine the image quality of the reconstruction. Therefore, photoacoustic reconstruction algorithms play an essential role in computedtomographic PAI.
Many efforts have been made to develop photoacoustic reconstruction algorithms. Analytical reconstruction algorithms were first developed, and their techniques are relatively mature [14,15,16,17,18]. The filtered backprojection (FBP) method proposed by Xu et al. was widely used due to its concision and convenience [16]. Zhang et al. proposed the deconvolution reconstruction algorithm, which achieved improved results in the case of both fullview and limitedview scanning [18]. To overcome the strong data dependency of the analytical reconstruction algorithms and improve their performance, the iterative image reconstruction methods were proposed. This kind of reconstruction methods established a forward model from photoacoustic image to photoacoustic signals to calculate the photoacoustic image iteratively [19,20,21,22,23,24,25]. Compressed sensing (CS) theory has been adopted in PAI to reduce the number of samples required and improve the results in sparseview scanning [26,27,28,29,30,31]. Among these algorithms, totalvariation (TV)based reconstruction algorithms have achieved excellent reconstruction quality [32,33,34,35,36,37,38]. The TV minimization can greatly reduce the dependence on data so that images can be accurately recovered from sparse data. Therefore it is potential to improve the performance of the algorithm on limitedview scanning based on TVmethod. An adaptive steepestdescentprojection onto convex sets (ASDPOCS) is proposed by Wang et al. to employ the TVbased iterative image reconstruction algorithms in threedimensional PAI [33]. Zhang et al. proposed a gradient descentbased TV (TVGD) algorithm, which was able to maintain good performance even in sparseview scanning [34]. A joint TV and Lpnorm (TVLp)based algorithm proposed by Zhang et al. was reported to have improved performance especially in the sparseview scanning [39]. Besides, wavelets transform domain [21, 40], total generalised variation [41] as well as deep learning regularization [42, 43] have also been adopted in PAI reconstruction and reported to have successfully addressed some specific problems in PAI. While for wavelets transform domain [21, 40] as well as total generalised variation [41]—based method, there still exists room for improvement in the preservation of structure and detail information particularly under the circumstance of limitedview scanning. As for deep learning based methods [42, 43], the algorithms are too complex and difficult to implement.
The image reconstruction methods at the present stage have worked well with fullview sampled data, but in practical situations, fullview scanning is often unavailable because of the restraint of the body shape or firmware. Under such circumstances, only limitedview projection data can be acquired, which do not conform to the condition of data completeness. In biomedical clinical practice, the linear transducer array is one of the popular ways to collect ultrasound signals. For clinical application, current PAI reconstruction algorithms still have many problems, such as edge blur and serious artifacts [28, 30, 37, 38, 44,45,46,47,48,49]. There is still much room for improvement. It is necessary to develop an image reconstruction method that is effective in clinical applications.
The TV expresses local intensity changes in an image. The classical TVbased reconstruction methods were established based on the assumption that the images are piecewise constant [50]. While the TV model has obtained a good effect in terms of sparseview reconstruction, due to the overinhibition of the highfrequency coefficients, minimizing the TV of an image tends to create oversmoothed geometry construction in the images [50,51,52]. The result is even worse in the case of practical limitedview scanning when some angular projection data are missing, as severe artifacts emerge and detailed information is lost [34, 37, 39]. In recent years, a nonlocal idea involving a priori knowledge that reveals the selfsimilarity of images has been proposed and widely used in image processing and reconstruction [53,54,55,56]. Minimizing TV can be regarded as minimizing the variation between adjacent pixels and can therefore be named local TV. Nonlocal TV extends the spatial neighborhood in the traditional neighborhood filtering to the structured neighborhood with a more generalized geometric meaning [56]. It searches similar patches in a larger area and uses the similarity between patches as the weight. This approach overcomes the limitation of traditional neighborhood weighting and makes better use of the similarities within images. Therefore, the reconstructed images can be improved in terms of texture and structure preservation. By solving the research and clinical problems, the method has obtained better performance in local TV [56,57,58].
In this paper, we propose a novel PAI reconstruction algorithm that incorporates nonlocal patchbased regularization into the TVbased model (patchTV) to improve the reconstruction results for practical straightline scanning. The patch in the image is estimated by weighting the patches in its neighborhood, which are searched throughout the whole image adaptively. The reconstructed image is updated by joint TV and nonlocalpatch regularization. The modified weighting calculation method is adopted with directivity and adaptability to further improve the performance of structure maintenance for the image [59]. Finally, the optimization model is simplified, and efficient variable splitting and the Barzilai–Borweinbased method are adopted to solve the optimization problem [60]. A series of numerical simulations and an in vitro experiment are conducted to validate the proposed patchTV algorithm. The results of the patchTV algorithm are compared to those of TVbased algorithms solved by the gradient descent method (TVGD), the TVLp algorithm as well as the iterative algorithm only with patchbased regularization (PatchRE). The peak signaltonoise ratios (PSNRs), the noise robustness, and the calculation and convergence speeds are also discussed and compared. Both qualitative and quantitative comparisons show that the patchTV algorithm provides superior results to those of TVGD, TVLp and PatchRE. The geometric structures of the images are preserved well, and the quality of the reconstructed images is greatly improved for practical straightline scanning. A series of patch based methods have been applied in imaging, such as [61]. In [61], nonlocal patch was used as a filter to process the image after the updating of each iteration step, which makes the algorithm one kind of image processing rather than image reconstruction. Moreover, the simple and isotropic distance between two blocks is utilized to screen the neighborhood of the block. In the proposed patchTV algorithm, nonlocal patch is used as a constraint item in the optimization problem for reconstruction. The optimization problem is then simplified to a common iterative PAI reconstruction problem so that the complexity of the algorithm is greatly reduced. The modified weighting calculation method which utilizes the modified structure tensor matrix to construct the weight function between two patches with directivity and adaptability is adopts in the proposed algorithm. The screened neighborhood of the patches takes the directivities and geometric structure of the images fully into consideration. It further improves the performance of structure preserving for the image. The nonlocalpatch regularization is combined with TV minimization in the proposed algorithm to obtain better performance in straightline scanning with stability.
There are mainly three points for the contributions of this paper. First, we include the nonlocal patch idea into PAI reconstruction. As far as we know, it is the first time that nonlocal patch ideal is applied to PAI. Second, the combination of the nonlocal patch optimization and TV minimization has been firstly applied into PAI. This combined method is able to solve the problems of PAI reconstruction from straightline scanning. Finally, we simplify the complicated optimization problem to a common iterative PAI reconstruction problem and use efficient variable splitting and the Barzilai–Borweinbased method to solve this problem. The optimization steps are greatly simplified and the convergence is greatly accelerated.
Theory and methods
A. TVbased photoacoustic reconstruction model
The algorithm proposed in this paper mainly targets twodimensional computedtomographic PAI for simple study. The possibility of extending the method to 3D will be discussed in “Discussion and conclusion”. In this imaging mode, laser pulses irradiate perpendicular to the image plane. Assuming that the tissue is irradiated uniformly by the laser, the relationship between the photoacoustic signals and the photoacoustic image can be described by the photoacoustic Equation [1]:
where p(r, t) is the photoacoustic signals at time t and position r, c is the speed of sound, μ is the isobaric expansion coefficient, C_{p} is the specific heat, I(t) is the temporal profile of the laser pulse and A(r) is the light absorption distribution of the tissue.
Assuming I(t) is an impulse signal and the sound velocity and other parameters of tissue are homogeneous, Eq. (1) can be solved by Green’s function [1]:
where r_{0} is the position of the ultrasound transducer.
Now, we establish the forward model from photoacoustic signals to a photoacoustic image. From Eq. (2), it can be derived that:
Define the product of detected photoacoustic signals at sampling points r_{0} and sampling time t, g(r_{0},t), as:
Equation (3) can be rewritten as:
In practical applications, the images and sampling signals tend to be discretized and can be written in the form of a vector [34]:
where A is the matrix of the photoacoustic image of size N_{x}× N_{y}, A′ is a column vector transposing A, l is the number of sampling points and M_{l} is weight matrix for the lth sampling point, g_{l} is the column vector discretized from g(r_{0}, t) for the lth sampling point.
An image’s gray values tend to have no sparsity, while its discrete gradients have more sparsity under some circumstances, such as homogeneous distribution of light in the sample and piecewise constant absorption coefficient.
TV can be expressed as the l_{1} norm of the discrete gradient matrix of the image [62]:
where A_{m,n} is the gray value of the pixel at the position (m, n).
The optimization problem of TVbased photoacoustic reconstruction can be written as:
where α is the parameter corresponding to the weight of TV value in the optimization. Equation (8) can also be written as:
where u_{i}= D_{i}A. D_{i} is a defined matrix that calculates the finite difference of A at the ith pixel.
B. Nonlocal patch regular constraint
There can be many similar patches in an image. In the flat region, most pixels and patches are identical, while the texture and edge regions also show similarities. Buades et al. therefore proposed the nonlocal idea and extended the similarities between pixels to that between patches [53]. For the nonlocal idea, a neighborhood is no longer for pixels in the common sense but is rather the patchset under a certain measure of similarity.
For pixel xi = (xi_{1}, xi_{2}), P_{xi} refers to the patch centered at xi. The selfsimilarity of the image can be represented in terms of the similarities between patches:
where W(x_{i}, x_{j}) is the weight function between P_{xi} and P_{xj}. It measures the similarity degree between the two patches and satisfies \(\sum\nolimits_{{{\mathbf{x}}j \in \delta ({\mathbf{x}}i)}} {W({\mathbf{x}}i,{\mathbf{x}}j)} = 1\). δ(x_{i}) refers to the neighborhood of P_{xi}:
where T is a threshold value to screen the similar patches. If the weight is greater than T, these two patches are considered similar. Otherwise, this patch does not belong to the neighborhood of patch P_{xi}. Equation (11) represents the collection of every pixel whose similarity to patch P_{xi} is greater than T.
There are multiple expressions for the weight function W(x_{i}, x_{j}), and it is usually inversely proportional to the distance between x_{i} and x_{j}. These weight functions failed to maintain the structure and directivity information of the image. So they are not qualified for the adaptive selection of the neighborhood of the patches. Liu et al. proposed the direction adaptive weight function [59], which is adopted in this paper:
where S_{j} is the modified structure tensor matrix. h is the global smoothing parameter and μ_{i} is the local density of samples data. More details can be found in Ref. [59]. The structure tensor matrix S_{j} reflects the information of gray values and gradients for the image. Using this directionadaptive weight function, the neighborhood δ(x_{i}) of patch P_{xi} can be adaptively selected. The selection of the neighborhood takes the directivity and geometric structure of the image fully into consideration, so it can provide more reliable estimations for the weight calculation between patches. Therefore, the structure and directivity information of the image can be well maintained.
The nonlocal patch regular constraint corresponding to the selfsimilarity between patches in Eq. (2) can be written as:
Patch P_{xi} is estimated by using the weights of patches in the neighborhood that have the highest similarities to P_{xi}. It is the first time that nonlocalpatch is applied as the regularized constraint for the reconstruction of image in PAI. By the constraint of the nonlocal patch, the problem concerning the inaccuracy of the similarity estimation through the use of isolated pixel points is surmounted, and the structure information, such as edges and texture, can be well preserved.
C. PatchTV photoacoustic reconstruction algorithm
The TVbased reconstruction model in Eq. (9) has good performance, but it fails to preserve the geometric structure of the image. To solve the problems of TV and make reconstruction algorithms more suitable for practical application, the nonlocal patch regular constraint is incorporated into the TVbased regular term:
where β is the parameter corresponding to the weight of localpatch value in the optimization. Define the nonlocal matrix H consisting of the weight functions W_{s}(x_{i}, x_{j}) [63]:
When x_{j} is in the neighborhood δ(x_{i}) of x_{i}, α_{ij} in H is set to the weight W_{s}(x_{i}, x_{j}). When x_{j} is not in the neighborhood δ(x_{i}) of x_{i}, α_{ij} is set to 0. In this way the summation item in the constraint item of localpatch can be expressed as multiplication between matrix H and A. Define H′ expressing the transversal vector transposing H. The size of H′ is 1 × (N^{2} × M^{2}). The optimization problem in Eq. (14) can be rewritten into the form of a matrix:
where I′ with the same size with that of H′ is the transversal vector transposing the unit matrix I. Combine the first and third terms in Eq. (16) in matrix form:
Using the notation \({\tilde{\mathbf{g}}} = \left[ {\begin{array}{*{20}c} {\mathbf{g}} \\ 0 \\ \end{array} } \right], \, {\mathbf{K}} = \left[ {\begin{array}{*{20}c} {{\mathbf{M}}^{{\mathbf{T}}} } \\ {\beta ({\mathbf{\rm I}}^{'}  {\mathbf{H}}^{'} )} \\ \end{array} } \right],\) Eq. (17) can be simplified as:
The patchTV optimization problem is simplified to a common photoacoustic iterative reconstruction model. The variable splitting and Barzilai–Borweinbased method is employed to solve the optimization problem in Eq. (18) [60]. This method has excellent performance in rapidly solving photoacoustic reconstruction regularized problems. Using the standard augmented Lagrangian method and the Barzilai–Borwein step size to accelerate the convergence speed, Eq. (19) can be deduced as [60, 64]:
where b _{ k} ^{n} is the TV step parameter in the nth iteration and σ_{n} is the defined Barzilai–Borwein step size in the nth iteration. By using the variable splitting method, Eq. (20) can be translated into the following two subproblems:
The two subproblems can be solved using the shrinkage operator method [60]:
where F is the Fourier transform matrix.
The flow of the patchTV photoacoustic reconstruction algorithm can be summarized as follows:

1.
Initialization: Input A, α, β, T. Set the reconstructed image A^{0} = 0, δ_{0} = 1, and b^{0} = 0.

2.
Apply Eq. (21) to update u^{n} for the given A^{n−1′}.

3.
Apply Eq. (22) to update A^{n} for the given u^{n}.

4.
Apply Eq. (22) to update b^{n} and δ_{n}.

5.
If the terminal condition is met, end the iteration. Otherwise, let n = n + 1, and return to steps 2–4. The termination condition is as follows:
$$\frac{{\left\ {u^{n}  u^{n  1} } \right\}}{{\left\ {u^{n} } \right\}} < \varepsilon .$$(23)
Numerical simulation
To verify the reconstruction quality and performance of the proposed patchTV algorithm, a variety of numerical simulations are designed and conducted. To simulate the signal collection in practice, straightline scanning with varying sampling points is executed. Straightline scanning in different directions to the phantom is also tested to validate the universality of the algorithm. The Shepp–Logan phantom, which is widely used in biomedical imaging, and the FORBILD phantom [65], which is more complicated and challenging, are chosen in the simulations. The results for the patchTV algorithm are compared to those of the TVGD and TVLp algorithms. The PSNR, the noise robustness and the convergence of the algorithms are also compared and discussed. The simulations are carried out using Matlab R2013a on a personal computer with a 2.4 GHz Intel(R) Xeon^{®} CPU and 64 GB memory. In the simulations, the sampling frequency is 200 MHz and the recording time of pressure waves is 20 μs for all the cases. The simulations for the signals and reconstructions are all conducted in the same twodimensional plane.
A. Straightline scanning
First, the Shepp–Logan phantom is adopted as the initial pressure rise distribution, which is shown in Fig. 1. The size of the phantom is 76.8 × 76.8 mm, and the reconstructed images size is set to 128 × 128 pixels. The scanning line on the right side of the phantom with the length of 76 mm is also shown in Fig. 1, from which we can see that the scanning line is parallel to the major axis of the ellipse of the phantom. We use the photoacoustic equation (Eq. 3 in paper) for the numerically produced simulated data and the forward projection model we described in the paper to reconstructed the image iteratively under patchTV regulation. Thus the inverse crime is avoided in our method during the generation of simulated signals. The distance from the center of the image to the scanning line is 38 mm. The length of the scanning line remains constant, while the sampling points can be 10, 20, or 50. The iteration number is set to 10 for all algorithms. The parameter settings for patchTV are estimated by testing the values that provides the best performance for the simulations. In this case, α = 0.4, β = 0.35, T = 0.65. The parameters for TVGD and TVLp are set referring [34, 39] to achieve the best performance in the simulations. The parameter settings for these algorithms are also estimated by testing the values that provides the best performance for the simulations.
The reconstruction results for the three algorithms are shown in Fig. 2. The images in this paper are normalized in the same gray level for comparison. The gray values of all pixels are divided by the maximum one in the images to avoid any effect to the quality of the images. In the first row of Fig. 2, the reconstructed images for TVGD have serious artifacts and blurred edges, which severely distort the images, especially in the vertical direction, where the angular information is missing. Regarding TVLp in the second row of Fig. 2, the result is improved over that of TVGD when the sampling points are sufficient. However, the quality of the reconstruction declines rapidly as the number of the sampling points decreases. We can see that for the 10point sparseview reconstruction in Fig. 2f, there is serious vagueness in the perpendicular direction of the image. As for PatchRE, in the third line, the results are even worse than those of TVLp and just slightly better than those of TVGD. It is because without TVoptimization to ensure the quality of the image in each iteration, the effects of the patch regularization will be greatly weaken. The results of patchTV in the third row of Fig. 2 show great improvement over the other two algorithms. The artifacts are effectively suppressed, and the edges of the image are distinct. The geometric structure of the images is preserved well, with almost no blur or distortion. Furthermore, a sharp decrease in the number of sampling points does not have a great effect on the quality of the reconstructed image.
The PSNRs of the reconstruction results for the four algorithms are also calculated and compared as the quantitative criteria for the evaluation of the reconstruction results. The greater the value of PSNR is, the better the reconstruction. The calculation formula of the PSNR is as follows:
where R_{m,n} is the gray value of the original image and MAXI is the maximum possible pixel value of the image. The original images which were not normalized are used for all the PSNR calculations in this paper. The PSNR results are displayed in Table 1.
Table 1 shows that patchTV obtains the highest PSNR values for every case. The PSNR values for TVGD are always low on account of the deficiency of the data for straightline scanning. In fact, the results of TVGD, are poor in all kinds of sampling condition even though when the sampling points are sufficient (50points). We can see that the PSNRs of TVGD are all lower than 20 dB. Under this circumstance, the amount of variation of PSNRs actually does not make much sense. TVLp has a good PSNR for 50point scanning, but the value of the PSNR decreases rapidly as the number of sampling points decreases. The PSNRSs of PatchRE are just slightly higher than that of TVGD. On average, the PSNR of patchTV is approximately 17 dB higher than that of TVGD, 8 dB higher than that of TVLp and 12 dB higher than that of PatchRE.
To test the universality of the algorithm in practical applications, we change the position of the scanning line relative to the phantom. In this case, the scanning line is parallel to the minor axis of the ellipse of the image. Its length and the distance to the center of the image remain unchanged. The numbers of sampling points are again 50, 20 and 10. The diagram of the scanning line is shown in Fig. 3. The parameter settings in this case is α = 0.50, β = 0.42, T = 0.65.
The results of the reconstruction for the three algorithms are shown in Fig. 4. We can see that there are a large number of blurs and distortions in the reconstructed images for TVGD, especially in the horizontal direction. The geometry structure information of the image is destroyed. TVLp and PatchRE fails to obtain ideal results, especially when the sampling points become sparse. Regarding patchTV, the edges and texture structure of the image are better preserved. The artifacts and background noise are effectively suppressed. Even in sparseview scanning, there is almost no blurring in the image.
We also compare the PSNRs of the results for the three algorithms in Table 2. The PSNR of patchTV is approximately 18 dB higher than that of TVGD, 10 dB higher than that of TVLp, on average and 14 dB higher than that of PatchRE.
To further validate the effectiveness of the proposed algorithm, the FORBILD phantom, which is more complex and challenging, is also adopted in the simulation. The phantom and the scanning line are shown in Fig. 5. The size of the phantom and the scanning settings are the same as those in Fig. 1. Fifty, 20, and 10point straightline reconstructions are conducted, and the results of the three algorithms are shown in Fig. 6. The parameter settings in this case is α = 0.65, β = 0.54, T = 0.57. TVGD and PatchRE shows poor performance, yielding bad image quality. The incompleteness of the data has a significant effect on the reconstruction. For TVLp, serious artifacts and blurring occur when the number of sampling points decreases. The contrasts of the images are not high, and the performance is not satisfactory. PatchTV overcomes these problems. The geometric structure of the phantom is distinct, and the artifacts are effectively suppressed.
The PSNR results of the three algorithms are displayed in Table 3. It is obvious that patchTV outperforms the other three algorithms for each sampling status, making the patchTV algorithm superior to the other two algorithms even for a complicated phantom.
B. Noise robustness
In PAI practical applications, it is important that the reconstruction algorithms have excellent noise robustness because the detected photoacoustic signals are usually disturbed by the system noise. The system noise follows a Gaussian distribution. To test the noise robustness of the proposed algorithm, the 20point sampled signals for the FORBILD phantom in “Straightline scanning” are supplemented with white noise and a signaltonoise ratio (SNR) of 10 dB, 5 dB or 0 dB. The parameter settings in this case is α = 0.73, β = 0.60, T = 0.54.
The reconstructed results for the three algorithms for the different SNR signals are shown in Fig. 7. TVGD, TVLp as well as PatchRE fail to maintain high performance, especially at a low SNR. The quality of the images decays seriously, the contrasts of the images decrease, and the artifacts and background noise cannot be suppressed or eliminated. PatchTV shows the highest performance in terms of noise robustness. The geometric structures of the reconstructed images are closer to those of the original image, and the noise is effectively suppressed.
The PSNRs of the reconstruction results are also displayed in Table 4. PatchTV outperforms the other three algorithms, and the advantages are more obvious when the noise energy is stronger.
C. Convergence and calculation
The convergence speed and calculation time are two other important performance indices for a photoacoustic iterative reconstruction algorithm. We define the distance between the reconstructed image and original image d as the quantization parameter:
The smaller d is, the smaller the difference between the reconstructed image and original image. We record d for every iteration step from 10point sampling of the FORBILD phantom in “Straightline scanning” and compare the d values of the four algorithms in each iteration in a line chart in Fig. 8. The results show that in every step, patchTV’s d value is smaller than those of the other three algorithms, and it convergences to the smallest value.
The time costs t of 50, 20, and 10point straightline reconstruction of the Shepp–Logan phantom in “Straightline scanning” for all four algorithms are also compared (Table 5). t calculates the time from input of the simulated data into the reconstruction algorithm to the output of the reconstructed image. The unit of t is second. The Barzilai–Borwe in method used in TVLp greatly accelerates the speed of the algorithm, and TVLp shows a greatly decreased time compared to TVGD. For patchTV, due to the incorporation of the nonlocal patch regularization, the time costs are higher than those of TVGD, TVLp and PatchRE. However, the performance of the algorithm is greatly improved, and the quality of the reconstructed images is enhanced significantly for practical applications.
According to the above simulations and discussion, patchTV is superior to the two popular TVbased algorithms and is a highly efficient photoacoustic image reconstruction algorithm.
Experimental results
To further validate and analyze the performance and practicability of the proposed algorithm, in vitro experiments were conducted. We used a singledetector platform to scan the Gelatin phantom linearly.
The diagram of the singledetector platform is shown in Fig. 9a. It included a Nd:YAG laser device (Surelite I, Continuum, San Jose, California, USA) to emit a laser pulse with a wavelength of 532 nm and a frequency of 10 Hz. The duration of the laser pulse was 4–6 ns. A single transducer (V383SU, Panametrics, Waltham, Massachusetts, USA) with a center frequency of 3.5 MHz and a bandwidth of 1.12 MHz was driven by a stepping motor scanning in the imaging plane. The sampling rate of the system was 16.67 MHz. The sampling frequency of the system is 16.67 MHz and the recording time of pressure waves is 50 μs. The experiment satisfied the American National Standards Institute (ANSI) laser radiation safety standard. The phantom for the straightline scanning is shown in Fig. 9b. The phantom was made of a gelatin cylinder with a black rectangular rubber sheet embedded into it as a light absorber. The radius of the cylinder was 25 mm, and the size of the light absorber was 9 × 14 mm. The scanning line, which was parallel to the longer side of the light absorber, was uniformly distributed with 41 sampling points. The sampling interval was 1 mm. The perpendicular distance from the center of the phantom to the scan line was 45 mm. The radius of the phantom was 25 mm the reconstructed images size was also set to 128 × 128 pixels. The parameter settings in this case is α = 0.55, β = 0.45, T = 0.60.
The reconstructed results for patchTV, TVLp and TVGD are shown in Fig. 10. PatchTV obtained the best imaging quality. There were serious artifacts and blurring in the images for the other two algorithms. Particularly for TVGD, serious distortions occurred in the vertical direction of the light absorber. The edges of the image were hard to recognize. The patchTV result was greatly improved. The edges of the image were distinct, and the distribution of the gray values was relatively uniform. Furthermore, the artifacts and background noise were effectively suppressed. This experiment further validates the effectiveness of the proposed patchTV algorithm. Under the circumstances of limitedview scanning in practice, patchTV outperforms the two mainstream TVbased algorithms and is a practical and efficient reconstruction algorithm for PAI.
Discussion and conclusion
In this paper, nonlocal patch regularization is incorporated into the TVbased photoacoustic imaging reconstruction model to effectively improve the performance in practical limitedview scanning. TVbased optimization minimizes the variation between adjacent pixels. It penalizes the local changes of the image and can therefore be referred to as local total variation. It is based on the assumption that the image is piecewise constant and oversuppresses the highfrequency coefficients. Thus, the geometric structure information of the reconstructed images tends to be oversmoothed. The result is even worse for practical limitedview scanning, in which the data information is insufficient such that serious artifacts and blurring fail to be effectively suppressed in the reconstructed images. However, in the nonlocal idea, the traditional spatial neighborhood is extended to the structured neighborhood in terms of geometric meaning, and the regularization is applied to patches in the whole image instead of only adjacent pixels [43]. Therefore, patchTV shows great improvement in terms of the preservation of the images’ geometric structure and has better results in preclinical applications. The similar patches for weighted calculation for a certain patch Pxi are searched in the entire image according to the value of the weight function W(xi, xj). A threshold value T is set to screen the neighborhood of the patch Pxi. This method overcomes the problems in traditional nonlocal means (NLM) filters, in which the size of the search field is settled and the patch Pxi is estimated by the patches in the determined search field. Thus, for large areas, the calculation costs are increased rapidly, while for small areas, similar patches far apart are missed. Therefore, the size of the neighborhood of the patch Pxi is adaptively controlled. Moreover, the modified weight function is adopted in this paper. It utilizes the anisotropic distance between two patches to adaptively adjust the search of the neighborhood direction. For example, for edge points, their similar patches are searched along the edge direction. In this case, the neighborhood can be an ellipse. The neighborhood of the patches takes the directivities and geometric structure of the images fully into consideration. Therefore, this approach makes more reliable estimations for the weight calculations between patches. The application of this modified weighting calculation method, can better maintain structural and directional information of the images because of its more reliable estimation for the weights between patches. Furthermore, the optimization problem combining nonlocal patch and TV is simplified to a common iterative reconstruction problem. Thus, the solution process is significantly simplified. The variable splitting method and the Barzilai–Borweinbased method are adopted to further accelerate the calculation and convergence speeds.
The proposed patchTV algorithm was validated by a series of simulations and an experiment. The simulations were conducted by means of straightline scanning, which is often used in practical applications. The reconstructed results of patchTV were compared to those of two mainstream TVbased algorithms: TVGD and TVLp. The results show that patchTV is superior to TVGD and TVLp, whether judged visually or in terms of PSNRs. The artifacts caused by the data incompleteness are effectively suppressed, and the geometric structure of the images is well maintained. Furthermore, the noise robustness, the convergence and the calculation speed are also discussed. The experiment carried out on an in vitro phantom adopted traditional straightline scanning with a single transducer. The results show that patchTV outperforms the other two algorithms in each case, with more distinct geometric structure and fewer artifacts.
In this paper, the study is under a systemspecific choice where the circumstance that laser pulses irradiate perpendicular to the image and not a result of having a 2Dreconstruction. While it is considered to be a common case which is easy to study. As for other cases, such as the light irradiate from other angles, we can use the Monte Carlo method in [66] to simulate the optical absorption distribution of the tissue. Actually, these cases mainly lead to the variation of the optical absorption distribution of the tissue yet the way to the algorithm study is the same.
The iteration number is set to 10 in this paper. As reported in [34, 39], the TVGD and TVLp algorithm converged when the number of iterations is 10, which was an appropriate choice for these algorithms. Also as shown in “Convergence and calculation”, the line chart of the distance d in Fig. 8 confirms that the distance versus the iteration curve for these algorithms converges when the number of the iterations is 10, which validates the convergence of these algorithms at 10th iteration.
As for the parameter setting, α is the parameter corresponding to the weight of TV value in the optimization. α with a big value means that the TVterm is dominant and the optimization is expected to have a quicker convergence. But oversized value will break the balance between the two parts of the objective function. The reconstructed images with oversized α will have a great difference from the real images because the data fidelity in the reconstruction is sacrificed to the image regularity. Based on this criterion, α should be set to a value which is neither too large nor too small when compared with the weights of the other part of the objective function to ensure good reconstructions, noise robustness and convergence speed. β is the parameter corresponding to the weight of localpatch value in the optimization. It has similar effects on reconstructions, noise robustness and convergence speed to α. T is a threshold value ranging from 0 to 1 for screening the similar patches. Small value of T means that more patches with smaller similarities will be included into the neighborhood δ(x_{i}) of x_{i}. It will diminish the effect of the constraint of localpatch and increase the time costs. While if T is set to an oversized value, few patches will be qualified for the neighborhood. So it may also degrade the performance of the algorithm. From the simulations and experiments, α can be set between 0.3 and 0.8, β can be set between 0.2 and 0.65, T can be set between 0.55 and 0.80.
It is also worth mentioning that the computation costs of patchTV are higher than those of the other two algorithms due to the incorporation of nonlocal patch regularization. However, the quality of the images is significantly improved, and the convergence speed is greatly accelerated. Additionally, the simplification of the optimization problem and the utilization of variable splitting and the Barzilai–Borweinbased method make the solution efficient and fast.
As for the 3D extension, i.e. 3D PA tomography, the proposed patchTV algorithm can be easily applied to it. The 3D PA tomography have the similar dataset and scanning mode with the 2D one. It’s also worth to mention that the patchTV framework has space independent nature. The implementations can be fulfilled to 3D image reconstructions that use spatial information. But if we want to solving a 3D image volume, further studies need to be carried out. As we mentioned above, the whole converge time and single iteration time of the proposed patchTV algorithm are just slightly more than TVGD and TVLp algorithms, which makes the 3D reconstructions practical.
In conclusion, the proposed patchTV algorithm is an effective and practical PAI reconstruction algorithm.
Abbreviations
 PAI:

photoacoustic imaging
 TV:

total variation
 TVGD:

gradient descentbased TV
 TVLp:

joint TV and Lpnorm
 PatchRE:

the iterative algorithm only with patchbased regularization
 PSNR:

peak signaltonoise ratio
 FBP:

filtered backprojection method
 CS:

compressed sensing
 patchTV:

the combined nonlocal patch the TV regularization
 SNR:

signaltonoise ratio
 NLM:

nonlocal means
References
 1.
Xu M, Wang LV. Photoacoustic imaging in biomedicine. Rev Sci Instrum. 2006;77(4):305–598.
 2.
Wang L. Prospects of photoacoustic tomography. Med Phys. 2008;35(12):5758–67.
 3.
Li C, Wang LV. Photoacoustic tomography and sensing in biomedicine. Phys Med Biol. 2009;54(19):R59–97.
 4.
Kim C, Favazza C, Wang LV. In vivo photoacoustic tomography of chemicals: highresolution functional and molecular optical imaging at new depths. Chem Rev. 2010;110(5):2756–82.
 5.
Stein EW, Maslov K, Wang LV. Noninvasive, in vivo imaging of bloodoxygenation dynamics within the mouse brain using photoacoustic microscopy. J Biomed Opt. 2009;14(2):020502.
 6.
Wang D, Wu Y, Xia J. Review on photoacoustic imaging of the brain using nanoprobes. Neurophotonics. 2016;3(1):1.
 7.
Yang X, Wang LV. Monkey brain cortex imaging by photoacoustic tomography. J Biomed Opt. 2008;13(4):044009.
 8.
Zhong J, Wen L, Yang S, Xiang L, Chen Q, Xing D. Imagingguided highefficient photoacoustic tumor therapy with targeting gold nanorods. Nanomed Nanotechnol Biol Med. 2015;11(6):1499–509.
 9.
Kumon RE, Deng CX, Wang X. Frequencydomain analysis of photoacoustic imaging data from prostate adenocarcinoma tumors in a murine mode. Ultrasound Med Biol. 2011;37(5):834–9.
 10.
Lu J, Gao Y, Ma Z, Zhou H, Wang RK, Wang Y. In vivo photoacoustic imaging of blood vessels using a homodyne interferometer with zerocrossing triggering. J Biomed Opt. 2017;22(3):036002.
 11.
Wang B, Su JL, Karpiouk AB, Sokolov KV, Smalling RW, Emelianov SY. Intravascular photoacoustic imaging. IEEE J Sel Top Quantum Electron. 2010;16(3):588–99.
 12.
Pu K, Shuhendler AJ, Jokerst JV. Semiconducting polymer nanoparticles as photoacoustic molecular imaging probes in living mice. Nat Nanotechnol. 2014;9(3):233–9.
 13.
Zerda ADL, Zavaleta C, Keren S, Vaithilingam S, Bodapati S, Liu Z, Levi J, Smith BR, Tejen MA, Oralkan O. Carbon nanotubes as photoacoustic molecular imaging agents in living mice. Nat Nanotechnol. 2008;3(9):557–62.
 14.
Kruger RA, Liu P, Fang YR, Appledorn CR. Photoacoustic ultrasound (PAUS)–reconstruction tomography. Med Phys. 1995;22(10):1605–9.
 15.
Mohajerani P, Kellnberger S, Ntziachristos V. Fast Fourier backprojection for frequencydomain optoacoustic tomography. Opt Lett. 2014;39(18):5455–8.
 16.
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–99.
 17.
Xu Y, Feng D, Wang LV. Exact frequencydomain reconstruction for thermoacoustic tomography–I: planar geometry. IEEE Trans Med Imaging. 2002;21(7):823–8.
 18.
Zhang C, Wang Y. Deconvolution reconstruction of fullview and limitedview photoacoustic tomography: a simulation study. J Opt Soc Am. 2008;25(10):2436–43.
 19.
Huang C, Wang K, Nie L, Wang LV, Anastasio MA. Fullwave iterative image reconstruction in photoacoustic tomography with acoustically inhomogeneous media. IEEE Trans Med Imaging. 2013;32(6):1097–110.
 20.
Paltauf G, Viator JA, Prahl SA, Jacques SL. Iterative reconstruction algorithm for optoacoustic imaging. J Acoust Soc Am. 2002;112(4):1536–44.
 21.
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–57.
 22.
Ding L, Luís DX, Lutzweiler C, Razansky D, Ntziachristos V. Efficient nonnegative constrained modelbased inversion in optoacoustic tomography. Phys Med Biol. 2015;60(17):6733–50.
 23.
Osborne MR, Presnell B, Turlach BA. On the lasso and its dual. J Comput Graph Stat. 2000;9(2):319–37.
 24.
Paltauf G, Nuster R, Haltmeier M, Burgholzer P. Experimental evaluation of reconstruction algorithms for limited view photoacoustic tomography with line detectors. Inverse Probl. 2007;23:S81–94.
 25.
Guo Z, Li C, Song L, Wang LV. Compressed sensing in photoacoustic tomographyin vivo. J Biomed Opt. 2010;15(2):021311(16).
 26.
Zhang J, Anastasio MA, Riviere PJL, Wang LV. Effects of different imaging models on leastsquares image reconstruction accuracy in photoacoustic tomography. IEEE Trans. Med. Imaging. 2009;28(11):1781–90.
 27.
Meng J, Wang LV, Ying L, Liang D, Song L. Compressedsensing photoacoustic computed tomography in vivo with partially known support. Opt Express. 2012;20(20):16510–23.
 28.
Haltmeier M, Berer T, Moon S, Burgholzer P. Compressed sensing and sparsity in photoacoustic tomography. J Opt. 2016;18:114004.
 29.
Betcke MM, Cox BT, Huynh N, Zhang EZ, Beard PC, Arridge SR. Acoustic wave field reconstruction from compressed measurements with application in photoacoustic tomography. IEEE Trans Comput Imaging. 2017;3(4):710–21.
 30.
AuthorSandbichler M, AuthorKrahmer F, AuthorBerer T, AuthorBurgholzer P, AuthorHaltmeier M. A novel compressed sensing scheme for photoacoustic tomography. SIAM J Appl Math. 2015;75(6):2475–94.
 31.
Hammernik K, Pock T, Nuster R. Variational photoacoustic image reconstruction with spatially resolved projection data. In: Proceedings of the SPIE 10064, photons plus ultrasound: imaging and sensing 10064, 100643I. International Society for Optics and Photonics; 2017.
 32.
Arridge S, Beard P, Betcke M, Cox B, Huynh N, Lucka F, Ogunlade O, Zhang E. Accelerated highresolution photoacoustic tomography via compressed sensing. Phys Med Biol. 2016;61(24):8908–40.
 33.
Wang K, Su R, Oraevsky AA, Anastasio MA. Investigation of iterative image reconstruction in threedimensional optoacoustic tomography. Phys Med Biol. 2012;57(17):5399–423.
 34.
Zhang Y, Wang Y, Zhang C. Total variation based gradient descent algorithm for sparseview photoacoustic image reconstruction. Ultrasonics. 2012;52(8):1046–55.
 35.
Wang K, Anastasio MA, Sidky EY, Pan X, Oraevsky AA. Limited data image reconstruction in optoacoustic tomography by constrained total variation minimization. In: Proceedings of the SPIE. 7899; 78993U1; 2011.
 36.
Lei Y, Jiang H. Photoacoustic image reconstruction from fewdetector and limitedangle data. Biomed Opt Express. 2011;2(9):2649–54.
 37.
Jiang H, Yao L. Enhancing finite elementbased photoacoustic tomography using total variation minimization. Appl Opt. 2011;50(25):5031–41.
 38.
Arridge SR, Betcke MM, Cox BT, Lucka F, Treeby BE. On the adjoint operator in photoacoustic tomography. Inverse Probl. 2016;32(11):115012(1–19).
 39.
Zhang C, Zhang Y, Wang Y. A photoacoustic image reconstruction method using total variation and nonconvex optimization. Biomed Eng Online. 2014;13(1):117.
 40.
Provost J, Lesage F. The application of compressed sensing for photoacoustic tomography. IEEE Trans Med Imaging. 2009;28(4):585–94.
 41.
Boink YE, Lagerwerf MJ, Steenbergen W, Van SG, Manohar S, Brune C. A framework for directional and higherorder reconstruction in photoacoustic tomography. Phys Med Biol. 2018;63(4):045018(118).
 42.
Hauptmann A, Lucka F, Betcke M, Huynh N, Adler J, Cox B, Beard P, Ourselin S, Arridge S. Model based learning for accelerated, limitedview 3D photoacoustic tomography. Imaging: IEEE Trans Med; 2018.
 43.
Antholzer S, Haltmeier M, Nuster R, Schwab J. Photoacoustic image reconstruction via deep learning. In: Proceedings of the SPIE, photons plus ultrasound: imaging and sensing, 10494, 104944U; 2018.
 44.
Gamelin JK, Aguirre A, Zhu Q. Fast, “limiteddata photoacoustic imaging for multiplexed systems using a frequencydomain estimation technique”. Med Phys. 2011;38(3):1503–18.
 45.
Wu D, Tao C, Liu X. Photoacoustic tomography in scattering biological tissue by using virtual time reversal mirror. J Appl Phys. 2011;109(8):084702.
 46.
Xu Y, Wang LV, Ambartsoumian G, Kuchment P. Reconstructions in limitedview thermoacoustic tomography. Med Phys. 2004;31(4):724–33.
 47.
Tao C, Liu X. Reconstruction of high quality photoacoustic tomography with a limitedview scanning. Opt Express. 2010;18(3):2760–6.
 48.
Gao H, Feng J, Song L. Limitedview multisource quantitative photoacoustic tomography. Inverse Probl. 2015;31(6):065004(1–24).
 49.
Feng J, Zhou W, Gao H. Multisource quantitative photoacoustic tomography with detector response function and limitedview scanning. J Comput Math. 2016;34(6):588–670.
 50.
Fei X, Wei Z, Xiao L. Iterative directional total variation refinement for compressive sensing image reconstruction. IEEE Signal Process Lett. 2013;20(11):1070–3.
 51.
Candes EJ, Wakin MB, Boyd SP. Enhancing sparsity by reweighted L1 minimization. J Fourier Anal Appl. 2008;14(5–6):877–905.
 52.
Tian Z, Jia X, Yuan K, Pan T, Jiang SB. Lowdose CT reconstruction via edgepreserving total variation regularization. Phys Med Biol. 2011;56(18):5949–67.
 53.
Buades A, Coll B, Morel JM. A review of image denoising algorithms, with a new one. Multiscale Model Sim. 2005;4(2):490–530.
 54.
Lou Y, Zhang X, Osher S, Bertozzi A. Image recovery via nonlocal operators. J Sci Compt. 2010;42(2):185–97.
 55.
Huang J, Ma J, Liu N, Zhang H, Bian Z, Feng Y, Feng Q, Chen W. Sparse angular CT reconstruction using nonlocal means based iterativecorrection POCS. Comp Bio Med. 2011;41(4):195–205.
 56.
LikformanSulem L, Darbon J, Smith EHB. Enhancement of historical printed document images by combining total variation regularization and nonlocal means filtering. Image Vision Comput. 2011;29(5):351–63.
 57.
Ertas M, Yildirim I, Kamasak M, Akan A. An iterative tomosynthesis reconstruction using total variation combined with nonlocal means filterin. Biomed Eng Online. 2014;13(1):65.
 58.
Liu HY, Wei ZH, Zhang ZR. Adaptive nonlocal patch regularization for image restoration. Chin J Electron. 2012;40(3):512–7.
 59.
Liu H, Wei Z. An edgeadaptive structure tensor kernel regression for image interpolation. In: Proceedings of the 2nd international confrence on future computer and communication. Wuhan: IEEE Press. 2010; 2:681–4.
 60.
Ye X, Chen Y, Huang F. Computational acceleration for MR image reconstruction in partially parallel imaging. IEEE Trans Med Imaging. 2011;30(5):1055–63.
 61.
Li X, Chen Z, Xing Y. Multi‐segment limitedangle CT reconstruction via a BM3D filter. In: IEEE nuclear science symposium and medical imaging conference (NSS/MIC); 2012. p. 2390–4.
 62.
Rudin LI, Osher S, Fatemi E. Nonlinear total variation based noise removal algorithms. Phys D. 1992;60:259–68.
 63.
Dong W, Zhang L, Shi G, Wu X. Image deblurring and superresolution by adaptive sparse domain selection and adaptive regularization. IEEE Trans Image Process. 2011;20(7):1838–57.
 64.
Afonso MV, BioucasDias JM, Figueiredo MAT. An augmented lagrangian approach to the constrained optimization formulation of imaging inverse problems. IEEE Trans Image Process. 2011;20(3):681–95.
 65.
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–52.
 66.
Wang L, Jacques SL, Zheng L. MCMLMonte Carlo modeling of light transport in multilayered tissues. Comput Method Program Biomed. 1995;47(2):131–46.
Authors’ contributions
Study concept and design (JW); drafting of the manuscript (JW); critical revision of the manuscript for important intellectual content (JW and YW); obtained funding (YW); administrative, technical, and material support (JW); study supervision (YW). Both authors read and approved the final manuscript.
Acknowledgements
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Avaliability of data and materials
The datasets used and/or analyzed during in current study are available from the corresponding author on reasonable requests.
Consent for publication
Not applicable.
Ethics approval and consent to participate
Not applicable.
Fundings
This work was supported by the National Natural Science Foundation of China (No. 11474071).
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Author information
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 Photoacoustic imaging
 Straightline reconstruction
 Nonlocal patch
 Total variation