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 :
{\nabla}^{2}p\left(\overrightarrow{r},t\right)\frac{1}{{c}^{2}}\frac{{\partial}^{2}p\left(\overrightarrow{r},t\right)}{\partial {t}^{2}}=\frac{\mu}{{C}_{p}}u\left(\overrightarrow{r}\right)\cdot \frac{\partial I\left(t\right)}{\partial t}
(1)
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:
p\left({\overrightarrow{r}}_{0},t\right)=\frac{\mu}{4\pi \phantom{\rule{0.12em}{0ex}}{C}_{p}}\frac{\partial}{\partial t}{\displaystyle {\u222f}_{\left\overrightarrow{r}{\overrightarrow{r}}_{0}\right=\mathit{ct}}\frac{u\left(\overrightarrow{r}\right)}{t}{\mathrm{d}}^{2}\overrightarrow{r}}
(2)
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:
f\left({\overrightarrow{r}}_{0},t\right)=\frac{4\pi \phantom{\rule{0.12em}{0ex}}{C}_{p}}{\mu}{\displaystyle {\int}_{0}^{t}p\left({\overrightarrow{r}}_{0},t\right)\mathrm{d}t}\cdot t
(3)
Then the equation (2) can be converted as follows:
f\left({\overrightarrow{r}}_{0},t\right)={\displaystyle {\u222f}_{\left\mathbf{r}\text{'}{\mathbf{r}}_{0}\right=\mathit{ct}}u\left(\overrightarrow{r}\right){d}^{2}}\overrightarrow{r}
(4)
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:
{f}_{i}={{A}_{i}}^{\mathrm{T}}\cdot \mathrm{u}\begin{array}{c}\hfill i=1,2,\cdots ,Q\hfill \end{array}
(5)
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}_{i}\left(j\right)=max\left\{1\left\frac{d\cdot \mathrm{d}x}{c\cdot \mathrm{d}t}j\right,\phantom{\rule{0.12em}{0ex}}0\right\}\phantom{\rule{1em}{0ex}}\left(1\le j\le M\right)
(6)

(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.
A=\left(\begin{array}{c}\hfill {{A}_{1}}^{\mathrm{T}}\hfill \\ \hfill {{A}_{2}}^{\mathrm{T}}\hfill \\ \hfill \vdots \hfill \\ \hfill {{A}_{Q}}^{\mathrm{T}}\hfill \end{array}\right)
(7)

(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:
{min}_{u}\left\{{\left{\mathrm{\Psi}}^{\mathrm{T}}\mathrm{u}\right}_{1}+\frac{1}{2}{{\leftA\mathrm{u}f\right}_{2}}^{2}\right\}
(9)
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:
\mathit{TV}\left(\mathrm{u}\right)={\displaystyle \int \leftD\mathrm{u}\right}={\displaystyle \sum _{i=1}^{N}{\left{D}_{i}\mathrm{u}\right}_{2}}\phantom{\rule{1em}{0ex}}i=1,2\dots N
(10)
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:
{min}_{u}\left\{\mathit{TV}\left(\mathrm{u}\right)+\frac{1}{2}{{\leftA\mathrm{u}f\right}_{2}}^{2}\right\}
(11)
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:
{min}_{u}\left\{\alpha \cdot \mathit{TV}\left(\mathrm{u}\right)+\beta {{\left{\mathrm{\Psi}}^{\mathrm{T}}\mathrm{u}\right}_{p}}^{p}+\frac{1}{2}{{\leftA\mathrm{u}f\right}_{2}}^{2}\right\}\phantom{\rule{1em}{0ex}}0<\mathit{p}\le 1
(12)
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:
{min}_{\mathrm{u},\mathrm{z},\omega}\phantom{\rule{0.5em}{0ex}}\left\{\alpha {\displaystyle \sum _{i}{\left{\omega}_{i}\right}_{2}}+\beta {\displaystyle \sum _{i}{\left{\mathrm{z}}_{i}\right}_{p}^{p}}+\frac{\rho}{2}{\leftA\cdot \mathrm{u}\mathrm{f}\right}_{2}^{2}\right\}\phantom{\rule{0.5em}{0ex}}i=1,2\dots N
(13)
where ρ is the parameter corresponding to the weight of the constraint condition in this optimization problem.
We form the augemented Lagrangian defined by
L(\omega ,z,\mathrm{u};{b}^{k},c{}^{k})=\alpha {\displaystyle \sum _{i}\left({\left{\omega}_{i}\right}_{2}+\frac{\rho}{2}{\left{\omega}_{i}{D}_{i}{\mathrm{u}}^{k}{b}_{i}^{k}\right}_{2}^{2}\right)+}\beta {\displaystyle \sum _{i}\left({\left{z}_{i}\right}_{p}^{p}+\frac{\rho}{2}{\left{z}_{i}{\mathrm{\Psi}}_{i}^{\mathrm{T}}{\mathrm{u}}^{k}{c}_{i}^{k}\right}_{1}^{2}\right)+\frac{1}{2}{\leftA\cdot \mathrm{u}\mathrm{f}\right}_{2}^{2}}
(14)
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
\left\{\begin{array}{l}\left({\omega}^{k+1},{z}^{k+1},{\mathrm{u}}^{\mathrm{k}+1}\right)={min}_{\omega ,z,\mathrm{u}}L(\omega ,z,\mathrm{u};{b}^{k},{c}^{k}),\hfill \\ {b}_{i}^{k+1}={b}_{i}^{k}\left({\omega}_{i}^{k+1}{D}_{i}{\mathrm{u}}^{k+1}\right),\hfill \\ {c}_{i}^{k+1}={c}_{i}^{k}\left({z}_{i}^{k+1}{\Psi}_{i}^{\mathrm{T}}{\mathrm{u}}^{k}\right),\hfill \end{array}\right.
(15)
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
\begin{array}{c}\hfill \left({\omega}^{k+1},{z}^{k+1},{\mathrm{u}}^{\mathrm{k}+1}\right)={min}_{\omega ,z,\mathrm{u}}\left\{\alpha {\displaystyle \sum _{i}\left({\left{\omega}_{i}\right}_{2}+\frac{\rho}{2}{\left{\omega}_{i}{D}_{i}{\mathrm{u}}^{k}{b}_{i}^{k}\right}_{2}^{2}\right)+}\beta {\displaystyle \sum _{i}\left({\left{z}_{i}\right}_{p}^{p}+\frac{\rho}{2}{\left{z}_{i}{\Psi}_{i}^{\mathrm{T}}{\mathrm{u}}^{k}{c}_{i}^{k}\right}_{1}^{2}\right)}\right.\hfill \\ \hfill \phantom{\rule{6em}{0ex}}\left(\right)close="\}">+\frac{{\delta}_{k}}{2}\left({\left{\omega}^{\mathrm{k}+1}{\omega}^{\mathrm{k}}\right}_{2}^{2}+{\left{z}^{\mathrm{k}+1}{z}^{\mathrm{k}}\right}_{2}^{2}+{\left\mathrm{u}{\mathrm{u}}^{k}+{\delta}_{k}^{1}{A}^{\mathrm{T}}\left(A{\mathrm{u}}^{k}f\right)\right}_{2}^{2}\right)\hfill \end{array}\n
(16)
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:
\left\{\begin{array}{l}{\omega}_{i}^{k+1}={min}_{{\omega}_{i}}\left\{{\left{\omega}_{i}\right}_{2}+\frac{\rho}{2}{\left{\omega}_{i}{D}_{i}{\mathrm{u}}^{k}{b}_{i}^{k}\right}_{2}^{2}+\frac{{\delta}_{k}}{2\alpha}{\left{\omega}_{i}{\omega}_{i}^{k}\right}_{2}^{2}\right\},\hfill \\ {z}_{i}^{k+1}={min}_{{z}_{i}}\left\{{\left{z}_{i}\right}_{p}^{p}+\frac{\rho}{2}{\left{z}_{i}{\Psi}_{i}^{\mathrm{T}}{\mathrm{u}}^{k}{c}_{i}^{k}\right}_{1}^{2}+\frac{{\delta}_{k}}{2\beta}{\left{z}_{i}{z}_{i}^{k}\right}_{2}^{2}\right\},\hfill \\ {\mathrm{u}}^{\mathrm{k}+1}={min}_{\mathrm{u}}\left\{\phantom{\rule{0.6em}{0ex}}\mathit{\alpha \rho}{\leftD\mathrm{u}{\omega}^{\mathrm{k}+1}\right}_{2}^{2}+\mathit{\beta \rho}{\left{\Psi}^{\mathrm{T}}\mathrm{u}{z}^{\mathrm{k}+1}\right}_{2}^{2}+{\delta}_{k}{\left\mathrm{u}\left({\mathrm{u}}^{k}{\delta}_{k}^{1}{A}^{\mathrm{T}}\left(A{\mathrm{u}}^{k}f\right)\right)\right}_{2}^{2}\right\},\hfill \\ {b}_{i}^{k+1}={b}_{i}^{k}\left({\omega}_{i}^{k+1}{D}_{i}{\mathrm{u}}^{k+1}\right),\hfill \\ {c}_{i}^{k+1}={c}_{i}^{k}\left({z}_{i}^{k+1}{\Psi}_{i}^{\mathrm{T}}{\mathrm{u}}^{k}\right),\hfill \\ {\delta}_{k+1}={\leftA\left({\mathrm{u}}^{\mathrm{k}+1}{\mathrm{u}}^{\mathrm{k}}\right)\right}_{2}^{2}/\left({\left{\omega}^{\mathrm{k}+1}{\omega}^{\mathrm{k}}\right}_{2}^{2}+{\left{z}^{\mathrm{k}+1}{z}^{\mathrm{k}}\right}_{2}^{2}+{\left{\mathrm{u}}^{\mathrm{k}+1}{\mathrm{u}}^{\mathrm{k}}\right}_{2}^{2}\right)\hfill \end{array}\right.
(17)
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:
\begin{array}{l}\left\{\begin{array}{l}{{\omega}_{i}}^{k+1}=max\left\{\u2225\frac{{a}_{1}{t}_{1}+{a}_{2}{t}_{2}}{{a}_{1}+{a}_{2}}\u2225\frac{1}{{a}_{1}+{a}_{2}},0\right\}\frac{1/\left({a}_{1}+{a}_{2}\right)}{\u22251/\left({a}_{1}+{a}_{2}\right)\u2225}\hfill \\ {a}_{1}={D}_{i}{\mathrm{u}}^{k}+{{b}_{i}}^{k}\hfill \\ {a}_{2}={{\omega}_{i}}^{k}\hfill \\ {t}_{1}=\rho \hfill \\ {t}_{2}={\delta}_{k}/\alpha \hfill \end{array}\right.\left(i=1,2\dots N\right)\end{array}
(18)
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:
\left\{\begin{array}{l}{{z}_{i}}^{k+1}=max\left\{\u2225\frac{{a}_{3}{t}_{3}+{a}_{4}{t}_{4}}{{a}_{3}+{a}_{4}}\u2225\frac{1}{{a}_{3}+{a}_{4}}{\u2225\frac{{a}_{3}{t}_{3}+{a}_{4}{t}_{4}}{{a}_{3}+{a}_{4}}\u2225}^{1p},\phantom{\rule{0.12em}{0ex}}0\right\}\frac{1/\left({a}_{3}+{a}_{4}\right)}{\u22251/\left({a}_{3}+{a}_{4}\right)\u2225}\hfill \\ {a}_{3}={{\Psi}_{i}}^{\mathrm{T}}{\mathrm{u}}^{k}+{{c}_{i}}^{k}\hfill \\ {a}_{4}={{z}_{i}}^{k}\hfill \\ {t}_{3}=\rho \hfill \\ {t}_{4}={\delta}_{k}/\beta \hfill \end{array}\right.\phantom{\rule{0.5em}{0ex}}\left(i=1,2\dots N\right)
(19)
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:
{\mathrm{u}}^{k+1}={F}^{\mathrm{T}}\left\{\frac{F\left[\mathit{\alpha \rho}\phantom{\rule{0.12em}{0ex}}{D}^{\mathrm{T}}{\mathrm{\omega}}^{k+1}+\mathit{\beta \rho}\phantom{\rule{0.12em}{0ex}}{\Psi}^{\mathrm{T}}{\mathrm{z}}^{k+1}+{\delta}_{k}{\mathrm{u}}^{k}{A}^{1}\left(A{\mathrm{u}}^{k}\mathrm{f}\right)\right]}{\mathit{\alpha \rho}\phantom{\rule{0.12em}{0ex}}{F}^{\mathrm{T}}{D}^{\mathrm{T}}\mathit{DF}+\mathit{\beta \rho}\phantom{\rule{0.12em}{0ex}}\mathrm{I}+{\delta}_{k}\mathrm{I}}\right\}
(20)
where F is the Fourier transform matrix.
As a result, the TV L
_{
p
} algorithm is summarized as follows:
\frac{\u2225\phantom{\rule{0.12em}{0ex}}{\mathrm{u}}^{k}{\mathrm{u}}^{k1}\phantom{\rule{0.12em}{0ex}}\u2225}{\phantom{\rule{0.12em}{0ex}}\u2225\phantom{\rule{0.12em}{0ex}}{\mathrm{u}}^{k}\u2225}<\u03f5
(21)

(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: