## 2 Theory and Methods  
### Model-Based Photoacoustic  
In computed-tomography PAI, we usually use a short laser pulse to vertically and uniformly incident upon the tissue. Then, the energy of the laser is absorbed by the tissue and the initial ultrasound field is inspired through the photoacoustic effect. The generated photoacoustic signals propagate in the medium and are detected by the ultrasound transducer scanning around the tissue. The relationship between the photoacoustic signals and the photoacoustic images could be described by the following equation [3]:  
\[\
abla^{2}p(\mathbf{r},t)-\frac{1}{c^{2}}\frac{\partial^{2}p(\mathbf{r},t)}{ \partial t^{2}}=-\frac{\beta}{C_{p}}A(\mathbf{r})\cdot\frac{\partial I(t)}{ \partial t} \tag{1}\]  
where \(p(\mathbf{r},\,t)\) is the generated acoustic pressure at the position \(\mathbf{r}\) and the time \(t\), \(A(\mathbf{r})\) is the light absorption deposition of the tissue, \(c\) is the speed of sound, \(\beta\) is the isobaric expansion coefficient and \(C_{p}\) is the specific heat. Assuming that the acoustic property of the medium is homogeneous and the laser pulse approximates the delta function, for the two-dimensional case, which is mainly concerned in this paper, Equation (1) can be solved by [PERSON]'s function [3]:  
\[p(\mathbf{r}_{0},t)=\frac{\beta}{4\pi C_{p}}\frac{\partial}{\partial t}|_{ \mathbf{r}-\mathbf{r}_{0}|=ct}\frac{A(\mathbf{r})}{t}d^{2}\mathbf{r} \tag{2}\]where \(\mathbf{r}_{0}\) is the position of the transducer. In the model-based PAI, a forward model is usually established to express the relationship between photoacoustic image and photoacoustic signals in matrix multiplication. From Equation (2), it can be derived that:  
\[\frac{4\pi C_{p}t}{\beta}{\int_{0}^{t}p(\mathbf{r}_{0},t)dt}=\iint\limits_{| \mathbf{r}-\mathbf{r}_{0}|=ct}A(\mathbf{r})d^{2}\mathbf{r} \tag{3}\]  
Defining the left side of Equation (3) as:  
\[g(\mathbf{r}_{0},t)=\frac{4\pi C_{p}t}{\beta}{\int_{0}^{t}p(\mathbf{r}_{0},t)dt} \tag{4}\]  
where \(g(\mathbf{r}_{0},t)\) is the integral for the sampled signal at \(\mathbf{r}_{0}\) multiplied with \(t\). Then, Equation (3) can be rewritten as:  
\[g(\mathbf{r}_{0},t)=\iint\limits_{|\mathbf{r}_{0}-\mathbf{r}^{\prime}|=ct}A( \mathbf{r}^{\prime})d\mathbf{r} \tag{5}\]  
In a practical reconstruction system, the photoacoustic image and the detected signals usually need to be discretized. The photoacoustic image is discretized to a matrix \(\mathbf{A}\) with the size of \(N_{x}\times N_{y}\) and \(g\) is discretized to the column vector \(\mathbf{g}\). From Equation (5), we can see that \(\mathbf{g}\) is the line integral of \(\mathbf{A}\) so every element in the vector \(\mathbf{g}\) can be expressed as the weighted sum of the elements in matrix \(\mathbf{A}\) and the weight value is related to the position of the integral arc [34]:  
\[\mathbf{g}_{l}=\mathbf{W}_{l}\mathbf{T}\cdot\mathbf{A}^{\prime},\ \ \ l=1,2,3,...,N \tag{6}\]  
where \(l\) refers to the number of the sampling point. (*)\(\mathbf{T}\) means transpose of the matrix. \(N\) is total number of the sampling points. \(\mathbf{A}^{\prime}\) is a column vector which is rearranged from \(\mathbf{A}\). \(\mathbf{W}_{l}\) is the corresponding weight matrix. As a method of the compressed sensing, the TV measures the local variation of an image which can be expressed as the \(L_{1}\) norm of the discrete gradient-matrix of the image [49]:  
\[\text{TV}(\mathbf{A})=\sum_{i,j}{[(A_{i,j}-A_{i-1,j})^{2}+(A_{i,j}-A_{i,j-1}) ^{2}]^{1/2}} \tag{7}\]  
where \(A_{i,j}\) is the gray value of the image at the pixel (\(i\),\(j\)).