## Abstract

The simplicity and computational efficiency of back-projection formulae have made them a popular choice in optoacoustic tomography. Nonetheless, exact back-projection formulae exist for only a small set of tomographic problems. This limitation is overcome by algebraic algorithms, but at the cost of higher numerical complexity. In this paper, we present a generic algebraic framework for calculating back-projection operators in optoacoustic tomography. We demonstrate our approach in a two-dimensional optoacoustic-tomography example and show that once the algebraic back-projection operator has been found, it achieves a comparable run time to that of the conventional back-projection algorithm, but with the superior image quality of algebraic methods.

© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Tomographic inversion problems, which pertain to the determination of 2D or 3D objects from their projections, are common in many imaging fields [1–8]. Most often, the projection data may be described by the Radon transform, or a modification thereof. While the Radon transform is most commonly known in the context of X-ray computerized tomography (X-ray CT) [1–3] it also appears in single-photon emission computed tomography (SPECT) [4], positron emission tomography (PET) [8], optical projection tomography (OPT) [5], and ultrasound tomography (UT) [6]. In the field of optoacoustic tomography (OAT), often referred to as photoacoustic tomography (PAT) [7], as well as in other forms of thermoacoustic tomography (TAT) [9], a variation of the Radon transform is used in which the projections are taken not over lines, but rather over arcs or spherical shells [10].

In OAT, most inversion algorithms fall under one of the following categories [10]: analytical and optimization-based. In analytical algorithms, it is assumed that the projection data are represented by a continuous function and the solution is written as a set of analytical operations that are applied on that function. Known examples of that approach are back-projection (BP) algorithms [11–15] and Fourier-based methods [16,17]. While analytical algorithms are often exact, approximate formulae are also very valuable when they are simple to compute [11]. Approximate BP formulae may also be used to improve the reconstruction quality in non-ideal imaging scenarios, as demonstrated in [18–20]. Since all analytical algorithms are eventually discretized and applied on non-ideal data, small errors in the analytical approximations do not necessarily have a practical importance [1]. In optimization-based algorithms, the operation describing the projection operation is first discretized, and the inversion is performed for the discretized problem using tools from optimization theory. In the literature, such numerical algorithms are known by different names, which emphasize different aspects: algebraic reconstruction technique (ART) [21], model-based algorithms [22–25], and iterative algorithms [26].

The main advantages of analytical algorithms are the relative simplicity in which they may be programmed and the possibility for computationally efficient implementations. Analytical algorithms, and in particular BP algorithms, are especially attractive in the fields of X-ray CT and OAT, where they achieve high-quality reconstructions in various detection geometries. Optimization-based techniques offer several advantages over analytical ones, which include exact compensation for additional effects in the projection formulae [27–30], compatibility with arbitrary detection geometries [10], higher tolerance to discretization errors and additive white noise [26,31], and the ability to use nonlinear tools from compressed-sensing theory such as sparsity in alternative representations [32] and least total variation [26,33]. The two major disadvantages of optimization-based techniques are the complexity of their implementation and their computational burden, especially when non-linear regularizers are used [26,32,33]. In the case of 3D imaging, images with relatively modest resolutions may take several minutes to reconstruct when efficiently implemented on graphical processing units (GPUs) [34]. In the case of volumetric images with high voxel counts, analytical formulae, even when not exact, remain the most practical option for real-time operation [35,36].

In this paper we develop a novel scheme for algebraic calculation of BP operators for OAT, which we term *algebraic back-projection* (ABP). While in standard algebraic inversion, one may invert the projection operator without making any assumptions on the nature of its inverse, in ABP it is assumed that the inverse operation consists of a single BP operator that is applied on all the projections. As we show in this work, this BP operation may be found by an iterative least-square minimization procedure. Thus, similarly to conventional BP formulae, ABPs are most appropriate when the scanning geometry contains symmetries, namely when the scanning surface is either translation or rotation invariant.

To efficiently implement ABP it is required that the discretization of the image and detection surface share the same symmetry, i.e. share the same invariances. Accordingly, in our work, we chose 2D OAT with linear scanning as our geometry for the implementation of ABP, which is translation invariant and thus compatible with image discretization on a Cartesian grid. We demonstrated ABP in numerical simulations, and obtained that it produced inverse operations in which each point in the projection data is back-projected over an arc in the image – the hallmark of optoacoustic BP formulae. In our numerical simulations, the reconstructions of ABP achieved comparable quality to those obtained by the conventional model-based approach, but with the runtime of a BP formula.

In contrast to BP formulae, there is no limitation in ABPs on the mathematical operation that describes the projections. Thus, one may include additional physical effects in the forward model, e.g. the effect of finite detector size in OAT [23], which may not always be accurately included in BP inversion formulae. For example, it has been previously shown that when focused detectors are used [27,28], e.g. in optoacoustic microscopy, approximate BP-based inversion formulae involve significant reconstruction errors. For such geometries, BP operators found via ABP may offer the benefit of accurate reconstruction without the high computational burden of conventional algebraic reconstruction methods.

The paper is organized as follows: In Section II, a short theoretical introduction to the tomographic problem in optoacoustic imaging is given. In Section III, the theory of ABP for the case of 2D OAT with linear scanning is developed. In Section IV, the performance of ABP is demonstrated in numerical simulations, and in Section V the results are discussed.

## 2. Image reconstruction in OAT

In OAT, the image, denoted by $H(r)$, represents the energy-deposition map in the tissue and is indicative of the optical absorption coefficient within the imaged specimen. The projection data relate to a set of 1D acoustic signals obtained at different positions of the ultrasound detector. The relation between the acoustic pressure $p(r,t)$at a given position r and time instant *t* to the optoacoustic image $H(r)$is often described by the following integral relation [10]:

The set of acoustic signals $p(r,t)$obtained over several locations represents the projection data used for reconstructing the image $H(r)$. When $H(r)$ has a finite support, and $p(r,t)$is known over a closed surface that surrounds $H(r)$, Eq. (1) may be uniquely inverted [10]. When the detection surface is either a sphere, a cylinder, or a plane, exact inversion may be achieved by the following BP formula [12]:

While Eq. (2), or variations thereof, is the most commonly used BP formula in optoacoustic tomography, it is important to note that numerous alternative formulae exist. For the case of infinite line detectors, Burgholzer *et al.* developed an exact 2D formula that was successfully tested on experimental data [13]. 2D formulae have also been developed for the case in which the optoacoustic image is restricted to a 2D surface [14,15]. However, these algorithms have not been tested on experimental data.

In algebraic inversion algorithms, both the optoacoustic image and projection data are discretized and represented by vectors, which we denote by$h$and $p$respectively. The discretization of Eq. (1) leads to the following matrix relation [11]:

where $M$is the model matrix that represents the integral in Eq. (1). The calculation of $M$may be performed via the interpolation of $H(r)$at positions that are not on the discretization grid or by assuming that $H(r)$is composed of atoms for which an analytical solution to Eq. (1) is known.In principle, the inversion of Eq. (4) may be performed for any detection surface by solving the following least squares minimization problem [25]:

where $\left|\right|\cdot |{|}_{2},$ is the ${\mathcal{l}}_{2}$ norm. For a given signal $p$, the solution of Eq. (5) may be found by using the pseudo-inverse of $M$ [37–39], or by using iterative methods such as LSQR [40], [41]. The advantage of the pseudo-inverse approach is that it provides a closed-form solution to Eq. (5), given bywhere ${M}^{\u2020}$is the pseudo-inverse matrix, which is given by ${M}^{\u2020}={\left({M}^{T}M\right)}^{-1}{M}^{T}$. Once ${M}^{\u2020}$has been pre-calculated, the reconstruction may be efficiently performed by calculating the matrix-vector product in Eq. (6). However, the computational burden of the pre-calculation of ${M}^{\u2020}$has limited its use. In practice, Eq. (5) is solved by iterative methods that do not require pre-calculations and involve acceptable reconstruction times [41–43], albeit considerably slower than those achieved by BP methods. Both the pseudo-inverse and iterative solutions to Eq. (5) may be extended to include Tikhonov regularization that is required when the projection data is insufficient for an accurate reconstruction of the image on the pre-prescribed grid, leading to an ill-conditioned matrix$M$. In the case of LSQR, regularization may also be performed by stopping the algorithm after a few iterations, before convergence is reached [40,42], and it is the method of regularization chosen in this work. A discussion on the numerical complexity of the different inversion methods is given in Sec. 3.2 and in the Conclusion.## 3. Algebraic calculation of back-projection operators

#### 3.1 Theory

In this section we develop the theory for ABP reconstruction in OAT for a rectangular grid with a linear scan of the acoustic detector. The parameters of the tomographic problem are depicted in Fig. 1(a). Briefly, our optoacoustic image is discretized over a linear 2D grid with ${N}_{x}$ and ${N}_{y}$pixels in the *x* and *y* directions, respectively, where the respective pixel sizes are${\Delta}_{x}$ and ${\Delta}_{y}$. We assume that the projections are obtained over a linear scan in the *x* direction with a step size of ${\Delta}_{x}$ and ${N}_{p}$ scanning positions, where ${N}_{p}>{N}_{x}$. In case the projection data is measured with a step size different than ${\Delta}_{x}$, it may be interpolated to a new grid with a step size of ${\Delta}_{x}$. For each scanning position, the acoustic signal is discretized with ${N}_{t}$ time points, where ${N}_{t}>{N}_{y}$. We further assume that the projections are acquired symmetrically with respect to the center of the image. For each position of the detector, the acoustic signal is discretized by the vector ${p}_{i}$, where the entire projection data are given by the following equation:

Similarly to analytical BP formulae, in ABP the same back-projection operation is applied on each data vector ${p}_{i}$. However, in contrast to analytical formulae, the back-projection operation in ABP is represented by a kernel matrix $K$ that back-projects each ${p}_{i}$over an image with ${N}_{y}$ pixels on the *y* axis and ${N}_{K}={N}_{x}+{N}_{p}-1$pixels on the *x* axis, as illustrated in Fig. 1(b). As shown in Fig. 2, the requirement ${N}_{K}={N}_{x}+{N}_{p}-1$ensures that all the back-projected images overlap with the grid of the image we wish to reconstruct.

The kernel matrix $K$may be represented as follows:

*n*th column in the back-projected image. Accordingly, ${K}_{n}{p}_{i}$represent the

*n*th column of the image back-projected from the

*i*th projection. As illustrated in Fig. 2, the reconstructed image is a sum of the translated back-projected images. Accordingly, the reconstructed image ${h}_{\text{rec}}$may be described by the following matrix relation:where the modified kernel matrix $\widehat{K}$ is given by

#### 3.2 Algorithmic implementation

The implementation of ABP, based on the derivations in the Appendix, may be summarized by the following steps:

- 1) Build the model matrix $M$for the geometry shown in Fig. 1(a).
- 2) Construct the matrix ${M}_{\text{temp}}=\left[{0}_{{N}_{t}{N}_{p}\times {N}_{y}({N}_{K}-{N}_{x})},M\right]$.
- 3) Set ${M}_{s}={M}_{\text{temp}}$.
- 4) Execute ${N}_{p}-1$time the following operations:
- 5) Build the matrix ${I}_{\text{temp}}=\left[{0}_{{N}_{x}^{2}\times {N}_{y}({N}_{K}-{N}_{x})},{I}_{{N}_{x}{N}_{y}}\right]$.
- 6) Set ${I}_{v}={I}_{\text{temp}}$.
- 7) Execute ${N}_{p}-1$ time the following operations:
- 8) Set
*n*= 1. - 9) Repeat ${N}_{t}$ times:
- 10) Construct the kernel matrix $K$by concatenating all the column vectors found in step 9, as shown in Eq. (A7a).
- 11) Build the matrix ${p}_{\text{mat}}$according to Eq. (A10) and calculate ${h}_{\text{mat}}=K{p}_{\text{mat}}$.
- 12) Construct the solution by using Eq. (A11).

Steps 1-10 in the algorithm are considered as pre-calculation and need to be performed only once for a given imaging geometry. In the pre-calculation, steps 1-7 involve merely rearrangement of matrix elements and their execution time is inconsequential when performed with sufficient random access memory (RAM) to store ${M}_{s}$. The computational burden of ABP is mostly in step 9, which involves solving an optimization problem for the matrix ${M}_{s}$. The reason for this high computational burden is that the matrix ${M}_{s}$ has ${N}_{p}$ times more entries than the model matrix $M$. When solving Eq. (A9) in iterations, each iteration of LSQR involves a complexity proportional to calculating ${M}_{s}{K}_{c,n}$ ${N}_{t}$times, i.e. $O\left({N}_{p}^{2}{N}_{t}^{2}{N}_{y}{N}_{K}\right)$. In contrast, step 11, which is individually performed for each dataset, involves a single matrix-matrix multiplication with a complexity of $O\left({N}_{t}{N}_{y}{N}_{K}{N}_{p}\right)$, i.e. ${N}_{t}{N}_{p}$ times less than a single iteration of step 9. Finally, step 12 involves summing ${N}_{p}$ back-projected images and thus has a complexity of $O\left({N}_{p}{N}_{y}{N}_{x}\right)$. Since in most practical application, ${N}_{x},{N}_{y},{N}_{p}$ and **${N}_{t}$** are of comparable magnitude, which we denote by $N$, the complexity of reconstructing an image with ABP is $O\left({N}^{4}\right)$

While the pre-calculation in ABP is high, it is still lower than that of the pseudo-inverse approach (Eq. (6)). The matrix multiplication and inversion steps in Eq. (6) lead to a complexity of $O\left[{\left({N}_{p}{N}_{t}\right)}^{2}{N}_{x}{N}_{y}+{N}_{p}{N}_{t}{\left({N}_{x}{N}_{y}\right)}^{2}+{\left({N}_{x}{N}_{y}\right)}^{3}\right]$. In addition, in the case studied in this work, regularization is required since the projections are obtained with a limited view. As a result, the calculation of the pseudo-inverse matrix requires a Tikhonov regularization term to Eq. (5) [38], and to perform the inversion for several values of the parameter, thus leading to a total computational burden of $O\left[{\left({N}_{p}{N}_{t}\right)}^{2}{N}_{x}{N}_{y}{N}_{\lambda}+{N}_{p}{N}_{t}{\left({N}_{x}{N}_{y}\right)}^{2}{N}_{\lambda}+{\left({N}_{x}{N}_{y}\right)}^{3}{N}_{\lambda}\right]$, where ${N}_{\lambda}$is the number of regularization parameters that are tested.

## 4. Numerical results

Our numerical simulations were performed for the 2D geometry illustrated in Fig. 1 with the parameters ${N}_{x}={N}_{y}=81$ and ${\Delta}_{x}={\Delta}_{y}=0.025$ cm, which correspond to an image size 2 × 2 cm – a typical image size for optoacoustic tomography [43]. The detector was scanned over a line at a distance of 1.5 cm from the center of the image grid and acquired ${N}_{p}=325$ projections. We note that taking ${N}_{p}$ to be considerably larger than ${N}_{x}$is essential for achieving a good tomographic coverage of the object and reducing limited-view artifacts [42]. The speed of sound used in the simulations was *c* = 1500 m/s. The time vector over which the projections were calculated had ${N}_{t}=647$ point with a time step of $\Delta /(3c)$. To avoid using the same model for the forward and inverse problems, the model matrix used for reconstruction assumed acoustic signals with a time step of $\Delta /(2c)$, corresponding to ${N}_{t}=431$ and a matrix size of 140075 × 6561, which occupied 156 MB when saved as a double sparse matrix in Matlab. Before reconstruction, the projections were interpolated from a$\Delta /(3c)$to a $\Delta /(2c)$time step. The simulations were performed without parallelization on a Lenovo ThinkStation P900 workstation with an Intel Xeon E5-2650 CPU operating at 2.3 GHz, 256 GB of RAM, and one terabyte of hard drive space.

To apply ABP, the matrices ${M}_{s}$ and ${I}_{v}$ in Eq. (17) were first calculated. The matrices had a size of 45524375 × 32805 and 45524375 × 431, and occupied 51 GB and 2.2 GB of RAM, respectively. Once the matrices had been calculated, the LSQR algorithm was used to solve the optimization problem in Eq. (A9). The first iteration of LSQR took approximately 17 minutes for each time instant, and approximately 5 days for the entire matrix$K$. Subsequent iterations took only 45 seconds per iteration for a given time instant, or approximately 5 hours for the entire matrix$K$. While in previous studies on conventional algebraic inversion in OAT, it was found that 50 iterations were sufficient for the LSQR algorithm to converge [21,35], in the current study 120 iterations were performed to enable us to properly present the convergence rate of the algorithm.

Figures 3(a)-3(c) show the back-projection operators obtained from $K$ for different times for the 120th iteration. As the figures shows, all the back-projections were on arcs and had a bipolar nature – the hallmarks of an optoacoustic back-projection operation (Eq. (2)). In Fig. 3(d), we show a 1D slice of the back-projection image of Fig. 3(b) taken over the *y* axis at *x =* 0 (solid curve). For comparison, Fig. 3(d) also shows the result obtained for the 1D slice when only a single iteration was used in the LSQR solution of Eq. (A9), in which $K$ was calculated (dashed curve). As the figure clearly shows, the slice from the 1^{st} iteration is non-zero only at two discrete time instants, at which the back-projected signal has opposite signs. This type of back-projected signal corresponds to a derivate operation, which is also the dominant operation in the 3D back-projection formula (Eq. (2)). As Fig. 3(d) shows, when more iterations were performed, a more elaborate back-projection structure emerged. We note that for all the time instants for which $K$ was calculated, the slices of the back-projected arcs had approximately the same temporal profile shown in Fig. 3(d).

Figures 4(a)-4(c) show the three test images used to evaluate the performance of the ABP algorithm. The first image (Fig. 4(a)), composed of identical point-like circles, was chosen to test the reconstruction isotropy, whereas the second image (Fig. 4(b)), which contained circles of various sizes, was chosen to check how different scales in the image affected the reconstructions. In the third image (Fig. 4(c)), the reconstruction was tested for lines, which imitate the structure of blood vessels. Three methods were used to reconstruct the images: the BP formula of Eq. (2), algebraic inversion via the application of LSQR on Eq. (5), and the proposed ABP.

The BP reconstructions are shown in Figs. 4(d)-4(f). Since Eq. (2) does not include the correct scaling factor for 2D images, its reconstructions were normalized. In Fig. 4(d), a good reconstruction is obtained, in which the major error is due to limited-view artifacts and smearing in the lateral direction. As expected, the artifacts become more significant the farther the objects are from the detector. In Fig. 4(e), an additional high-pass effect is observed. Indeed, it has been previously documented that when the BP formula in Eq. (2), which is exact for 3D, is applied on a 2D problem, the resulting reconstruction is high-passed compared to the originating image [25, 28]. The reconstructions in Fig. 4(d)-4(f) suffered from negative values and excess of background noise, similarly to previous studies [22]. The reconstruction time for each image was approximately 0.17 seconds, which is similar to the run times reported in Ref [44]. for a similar image size.

Figures 4(g)-4(i) show the algebraic reconstructions via Eq. (5) when the LSQR algorithm was performed with a single iteration, whereas Figs. 4(j)-4(l) show the results of the 120^{th} iteration. The reconstruction took approximately 0.65 seconds for a single iteration and approximately 43.5 seconds for 120 iterations. When only a single iteration of the LSQR was performed, the reconstructions exhibited the same smearing effect and accentuation of small objects as those observed in the BP reconstructions (Figs. 4(d)-4(f)). In addition, the images also exhibited a depth-dependent attenuation. However, when 120 iterations were performed, the reconstructions were almost identical to the originating images, with only minor limited-view artifacts.

Figures 4(m)-4(o) show the reconstructed images obtained by ABP with a single iteration, and Figs. 4(p)-4(r) show the result of the 120th iterations. Once the matrix $K$had been pre-calculated, the ABP reconstructions required only 0.17 seconds, regardless of how many iterations were used to pre-calculate $K$. In particular steps 11 and 12 in the ABP algorithm (Section IIIB) required 0.12 and 0.05 seconds, respectively.

The first iteration of the ABP reconstructions suffered from the same distortions observed in the first iteration of the algebraic reconstructions (Figs. 4(g)-4(i)). Similarly to the case of the algebraic reconstructions, in the 120th iteration of the ABP Similarly to the case of the algebraic reconstructions, in the 120th iteration of the ABP reconstructions the distortions were significantly reduced, and only weak limited-view artifacts remained (Fig. 4(p)-4(r)), albeit more visible than those obtained in the algebraic reconstructions.

In the last numerical example, we compared the performance of the reconstruction algorithms in the presence of noise for the images in Fig. 4(a)-4(c). To simulate a noisy measurement, we added to each point in the projection data a random variable with a Gaussian distribution, zero mean, and standard-deviation of 20% of the maximum value in the projection data, where all the random variables were uncorrelated.

We used two measures to assess the quality of the reconstructions. The first measure was the root-mean-square error (RMSE) of the images when compared to their respective originating images in Fig. 4(a)-4(c). For the algebraic inversion, the reconstructions with the lowest RMSE were respectively obtained at the 2^{nd}, 8^{th}, and 3^{rd} iterations for the images of Figs. 4(a)-4(c). Figure 5 shows the reconstructions obtained using the BP formula of Eq. (2) (Figs. 5(a)-5(c)), the algebraic reconstruction of Eq. (5) performed with 120 iterations (Figs. 5(d)-5(f)) and with the number of iterations that produced the lowest RMSE (Figs. 5(g)-5(i)), and ABP with 120 iterations (Figs. 5(j)-5(l)). A 1D profile of the reconstructions of Figs. 5(a), 5(d), 5(g) and 5(j) over the horizontal line *y =* 0 is shown in Fig. 6.

The RMSE values for the noiseless and noisy cases are shown in Figs. 7(a)-7(c) as a function of iteration for standard algebraic inversion (Eq. (5)) and ABP. The figure shows that while the algebraic reconstructions achieved a lower RMSE in the noiseless case, they obtained a higher RMSE than ABP in the noisy case. As Figs. 7(a)-7(c) show, the RMSE of ABP in all three examples converges to a value that is larger than zero. The result is expected because of the restriction put on ABP that all the signals are back-projected with the same kernel. In contrast, in conventional algebraic inversion, the noiseless case would yield an RMSE that practically decays to zero, up to rounding errors, assuming the same model matrix is used in the forward and inverse problems. Since our case involved a limited-view reconstruction, the RMSE decay of the algebraic inversion was slow, as expected in an ill-conditioned problem [40,42]. Interestingly, for both ABP and the algebraic reconstruction, the slowest decay was obtained for the image in Fig. 4(b). This result may be explained by the different scales exhibited in the image of Fig. 4(b) in contrast to the images of Figs. 4(a) and 4(c) that exhibit structures of the same scale. In the noisy case, the lowest RMSE of the algebraic reconstruction was obtained by terminating the LSQR before the 120^{th} iteration. We note that in [42], terminating the LSQR algorithm before convergence was used as a way to regularize the algebraic inversion. For ABP, the susceptibility of the later iterations to noise was found to be less significant than in the algebraic reconstructions.

The second measure of quality was the signal-to-noise ratio (SNR) in the noisy reconstructions, calculated by dividing the maximum value of each image by the standard deviation in a 10 × 10 pixel area on the top left of the image, a region in which the originating images are identically equal to zero. Figures 7(d)-7(f) respectively show the SNR of the reconstructions achieved by BP, algebraic inversion, and ABP for the test images of Figs. 4(a)-4(c). For all the test images, ABP achieved the highest SNR, while the BP formula achieved the lowest. This quantitative result can also be observed qualitatively when comparing the images in Fig. 5. We note that in Fig. 7(e), the SNR of ABP improved significantly in the first few iterations, in contrast to the other examples, because of the low signal level obtained in the first iteration (Fig. 5(n)). Thus, while the noise in all the reconstructions increased in the first iterations, only in Fig. 7(e) the increase in the signal was sufficient to improve the SNR.

## 5. Experimental results

To further test the performance of ABP it was tested on experimental optoacoustic data. The optoacoustic excitation was performed with an OPO pulse laser with a repetition rate of 100 Hz and pulse energy of approximately 30 mJ, tuned to λ = 680 nm (SpitLight DPSS 100 OPO, InnoLas Laser GmbH, Krailling, Germany). The imaged object was a dark polyethylene microsphere with a diameter of 0.3 mm (Cospheric LLC, Santa Barbara, California), which was embedded in transparent agar. The detection was performed with a needle hydrophone with a diameter of 1 mm and an acoustic bandwidth of 12 MHz (Precision Acoustics, Dorchester, UK), which was scanned at a distance of approximately 6 mm from the microsphere over a span of 2 cm with a step size of 0.25 mm. The voltage signal was digitized by an oscilloscope with a bandwidth of 70 MHz (Rigol Technologies, Beaverton, OR, USA) and averaged 16 times per position.

Similar to the previous section, the image was reconstructed over a grid of 1 × 1 cm, at a distance of 0.5 cm from the detector. Figures 8(a)-8(c) show the reconstruction results obtained using BP, algebraic inversion, and ABP, respectively. Since the image included only a single microsphere, only a 5 × 5 mm square of the image is presented. The algebraic reconstruction was performed with 5 iteration, where further iterations did not change the magnitude of the microsphere in the image. Figure 8(d) shows a 1D plot of the microsphere reconstruction, obtained on the horizontal line of *y =* -3.72 mm in Figs. 8(a)-8(c). While all three methods obtained a similar reconstruction of the microsphere, the reconstruction of BP suffered from a significantly higher noise level than the algebraic reconstruction and ABP.

## 6. Conclusion

In conclusion, we developed a method for the algebraic determination of back-projection operators. The method, which we term algebraic back-projection (ABP), finds the optimal back-projection operator in the least-square sense under the restriction that the same operator is applied on each of the projections. Once the back-projection operator, represented by the kernel matrix $K$, is found, the reconstruction is achieved by a matrix-multiplication operation followed by a summation operation (Eq. (23)). Since $K$ is pre-calculated for a given geometry and does not depend on the projection data, ABP may be implemented with a numerical efficiency comparable to the one achieved by analytical BP formulae.

While the purpose of this paper is to introduce the general framework of ABP, which represents a new type of solution to tomographic problems, the specific implementation was for the problem of 2D OAT with a linear scanning of the acoustic detector. This specific geometry simplified the implementation of ABP in two ways. First, in this geometry both the projections and the image were discretized in Cartesian coordinates, simplifying the matrix operations required to describe the translation of the acoustic detector or of the image. Second, the use of a 2D example enabled the implementation of ABP without parallel programming. In our numerical simulations, the kernel matrix was calculated iteratively using the LSQR algorithm. When LSQR was applied with a single iteration, the result was a back-projection of the derivative of the acoustic signals onto arcs, which is the well-known far-field approximation to the 3D BP formula (Eq. (2)) [11]. Indeed, when we used the single-iteration ABP in our 2D geometry, we obtained high-passed reconstructions similar to those obtained when the analytical BP formula (Eq. (2)) was applied. When higher iterations were used to calculate **K**, finer details were added to the back-projection operation, though it still maintained its general structure of arcs with derivative operations.

In terms of run time, in our numerical simulations ABP was comparable to BP and considerably faster than iterative algebraic inversion. In iterative algebraic inversion, each iteration of the LSQR algorithm requires multiplying the model matrix $M$by a vector [40], leading to a complexity of $O\left({N}^{4}\right)$per iteration. In contrast, in ABP the total complexity per reconstructed image was $O\left({N}^{4}\right)$. While the complexity of the BP algorithms is only $O\left({N}^{3}\right)$ [40], they involve more complex calculations per pixel than ABP in which the complexity relates to an efficiently implementable matrix-matrix multiplication, and as a result the run times were comparable in our numerical example. We note that future implementations of ABP may achieve a complexity of $O\left({N}^{3}\right)$by utilizing the sparsity of$K$, visible in Fig. 3, and performing the multiplications only for the significant elements of $K$.

The pre-calculation of ABP involves a complexity of $O\left({N}^{5}\right)$per iteration, whereas the pseudo-inverse approach has a complexity of $O\left({N}^{6}\right)$. While efficient methods to calculate the pseudo-inverse have been developed using frequency-domain [45] and combined-space [46] representations, these methods did not reveal the underlying back-projection kernel in time domain shown in Fig. 3. In addition, in contrast to the current work, the methods in Refs [45,46] have not been tested in limited-view scenarios in which the angular coverage was smaller than 180°.

In terms of reconstruction accuracy, the results obtained by ABP in our simulations were visually comparable to those obtained by standard algebraic inversion. While in the noiseless case the standard algebraic inversion obtained a lower reconstruction error for both images tested, it was more sensitive to noise than ABP. Thus, we conclude that the restriction in ABP to a single back-projection operator acts as a regularizer which suppresses noise at the price of moderate loss of image fidelity. In comparison to ABP and standard algebraic inversion, the analytical BP formula obtained the lowest image fidelity and was the most susceptible method to noise. We note that since the geometry studied involved a limited-tomographic view, the reconstruction errors were manifested also by negative image artifacts around the imaged objects and loss of tangential resolution [10].

Since ABP can achieve run times comparable to conventional BP formulae, while potentially offering better noise suppression, it may prove useful even in cases in which exact BP formulae exist. In addition, the generic framework of ABP may enable future implementations of this approach whose benefits go well beyond noise suppression. In principle, ABP may be applied to any linear tomographic problem, regardless of what operation is used to calculate the projection data. The performance of ABP will be determined by whether the discretization of the image and projection data have the same invariants. When they do, symmetry will require that the optimal back-projection operation be the same for all projections, making ABP equivalent to standard algebraic inversion. An example of such a scenario is cross-sectional imaging with ring detectors in OAT [7] with a polar image grid [45]. In 3D OAT, spherical or cylindrical grids may be used to develop ABP for spherical or cylindrical scans, respectively. A cylindrical grid may also be used for developing ABP for the recently introduced spiral scan [36]. When the projection data do not possess the symmetry of the image grid, ABP will find the best compromise, in the least-squares sense, for the back-projection operator.

One possible direction in which ABP may be generalized is the inclusion of complex projection operations, such as those found in practical OAT systems. For example, when the dimensions of the acoustic detector are larger than the acoustic wavelength, a spatial-averaging operation is added to Eq. (1) [10]. This operation leads to projection formulae that exhibit a complex dependency on the detector geometry for which no exact BP formulae are known; a high-fidelity reconstruction in such cases necessitates the use of algebraic methods [27,28]. Two more effects that may be included in the OAT model are non-uniform illumination [29] and signal attenuation [27]. Since all the above-mentioned effects are identical for all the projections, their inclusion is not expected to affect the performance of ABP with respect to the performance of standard algebraic inversion.

One of the major challenges of extending ABP to other imaging scenarios is computational. In our work, for a 2D image with a modest pixel count of 81 × 81, the corresponding matrix ${M}_{s}$ occupied 51 GB of RAM, requiring a workstation to implement the pre-calculations of the algorithm. In addition, 120 iterations performed to calculate $K$ required approximately six weeks when executed on a single processor. Clearly, directly extending the algorithm to higher resolutions or 3D may lead to an impractical computational burden. However, this computational burden may be decreased by minimizing the redundancy in the over-determined inversion problem of Eq. (A9). For example, in the simulations discussed in this work, the number of equations in Eq. (A8) was over 3 orders of magnitude larger than the number of unknowns. Thus, it may be possible to drastically reduce the number of linear equations without sacrificing the accuracy of the reconstruction, facilitating the extension of ABP to images with higher pixel/voxel count.

A different approach for accelerating the pre-calculation in ABP is to use parallel computing. Since the complexity of the pre-calculation stems from performing the ${N}_{t}$ independent calculations in Step 9 in Section IIIB, an improvement on the scale of ${N}_{t}$ may be achieved by using a similar number of processors. Since the pre-calculation is performed only once for a given geometry, a computer cluster may be used to perform the pre-calculation offline. Alternatively, parallelization may be achieved by using graphical processing units (GPUs). Although the matrix ${M}_{s}$ is too large to be saved in the RAM of conventional GPUs, a more efficient representation of its entries may enable one to implement ABP on a GPU. In particular, the matrix ${M}_{s}$ is not composed of unique entries, but rather contains ${N}_{p}$ copies of$M$, which in our example occupied only 156 MB. Thus, the algorithm may be implemented on a GPU with the matrix $M$ stored in the RAM of the GPU, with elements of the matrix ${M}_{s}$ calculated on-the-fly in the iterative solution of Eq. (A9). Moreover, the matrix $M$ itself may be generated from an even smaller set of data, as previously demonstrated in [34,47]. Thus, ABP may be adapted for GPU implementation, which may enable its generalization to 3D.

Finally, we note that the information found in the low-resolution BP operator may be extended to enable image reconstruction at arbitrary scales if an analytical operation is empirically fitted to the numerical BP operator. For example, the result of the first iteration of ABP found in our paper may be approximated by projecting on the arcs the temporal derivatives of the acoustic signals and applying an angle dependent weighting function on the arc. For higher iterations, the signal that is back-projected would need to be a convolution of the acoustic signal with the 1D BP response found in Fig. 3(d). This approach may enable finding empirical filtered back-projection formulae, analogues to those based on the Ram-Lak filter in computerized tomography [48].

## Appendix

In this Appendix, a detailed derivation of the algorithm in Sec. 3.2. is provided. As Eq. (10) shows, the matrix $\widehat{K}$ has a size of ${N}_{x}{N}_{y}\times {N}_{p}{N}_{t}$ and is composed of ${N}_{K}={N}_{x}+{N}_{p}-1$ sub-matrices of size ${N}_{y}\times {N}_{t}$ that appear at several positions in the combined matrix. Thus, although the matrix $\widehat{K}$ contains ${N}_{x}{N}_{y}{N}_{p}{N}_{t}$unknown elements, the number of its unique entries is only${N}_{K}{N}_{y}{N}_{t}$. In comparison, the number of linear equations in Eq. (12) is equal to the number of elements in ${I}_{{N}_{p}{N}_{t}}$, i.e. ${\left({N}_{p}{N}_{t}\right)}^{2}$. Since ${N}_{p}>{N}_{x}\gg 1$ and ${N}_{t}>{N}_{y}$, the number of equations vastly exceeds the number of unknowns.

To determine $\widehat{K}$ via algebraic inversion, Eq. (12) first needs to be rearranged in a way in which the matrix that contains the unknowns has only unique entries. As a first step, we break $\widehat{K}$into ${N}_{p}$ matrices${K}_{v,1}\dots {K}_{v,{N}_{p}}$, each with the size of ${N}_{x}{N}_{y}\times {N}_{t}$,

Similarly to Eq. (12), Eq. (A5) is a matrix representation of ${\left({N}_{p}{N}_{t}\right)}^{2}$linear equations. However, in contrast to the matrix $\widehat{K}$ in Eq. (12), the matrix $K$ has only unique entries. To invert Eq. (A5) and obtain $K$, we first decompose the matrices $K$and ${I}_{\nu}$in Eq. (A5) into a set of ${N}_{t}$column vectors:

Once all the vectors ${K}_{c,n}$ have been found via Eq. (A9), they are used to construct the kernel matrix $K$ via Eq. (A7a). To reconstruct the images, one may construct $K$ via Eq. (10) and then use Eq. (9). Alternatively, one may rearrange the projection data in the following matrix form:

*n*th column vector of ${h}_{\text{mat}}$ by ${h}_{\text{mat,}n}$the translation and summation of the back-projected images may be obtained by the following equation, which is equivalent to the operation shown in Fig. 2:

*i*th to the

*j*th elements of ${h}_{\text{mat,}n}$.

## Funding

The Israel Science Foundation (942/15); the Ministry of Science and Technology and Space (3-12970).

## Acknowledgment

The author acknowledges the assistance of Yoav Hazan and Ahiad Levi with the experimental work.

## Disclosures

The author declares that there are no conflicts of interest related to this article.

## References

**1. **X. Pan, E. Y. Sidky, and M. Vannier, “Why do commercial CT scanners still employ traditional, filtered back-projection for image reconstruction?” Inverse Probl. **25**(12), 123009 (2009). [CrossRef] [PubMed]

**2. **W. A. Kalender, “X-ray computed tomography,” Phys. Med. Biol. **51**(13), R29–R43 (2006). [CrossRef] [PubMed]

**3. **W. Kalender, *Computed Tomography: Fundamentals, System Technology, Image Quality, Applications,* 3rd ed. (Wiley-VCH, 2011).

**4. **R. M. Lewitt and S. Matej, “Overview of methods for image reconstruction from projections in emission computed tomography,” Proc. IEEE **91**(10), 1588–1611 (2003). [CrossRef]

**5. **J. Sharpe, U. Ahlgren, P. Perry, B. Hill, A. Ross, J. Hecksher-Sørensen, R. Baldock, and D. Davidson, “Optical projection tomography as a tool for 3D microscopy and gene expression studies,” Science **296**(5567), 541–545 (2002). [CrossRef] [PubMed]

**6. **C. Li, N. Duric, P. Littrup, and L. Huang, “In vivo breast sound-speed imaging with ultrasound tomography,” Ultrasound Med. Biol. **35**(10), 1615–1628 (2009). [CrossRef] [PubMed]

**7. **A. Taruttis and V. Ntziachristos, “Advances in real-time multispectral optoacoustic imaging and its applications,” Nat. Photonics **9**(4), 219–227 (2015). [CrossRef]

**8. **B. Bendriem and D. W. Townsend, *The Theory and Practice of 3D PET* (Springer, 1998)

**9. **D. Razansky, S. Kellnberger, and V. Ntziachristos, “Near-field radiofrequency thermoacoustic tomography with impulse excitation,” Med. Phys. **37**(9), 4602–4607 (2010). [CrossRef] [PubMed]

**10. **A. Rosenthal, V. Ntziachristos, and D. Razansky, “Acoustic Inversion in optoacoustic tomography: a review,” Curr. Med. Imaging Rev. **9**(4), 318–336 (2014). [CrossRef] [PubMed]

**11. **P. Burgholzer, G. J. Matt, M. Haltmeier, and G. Paltauf, “Exact and approximative imaging methods for photoacoustic tomography using an arbitrary detection surface,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. **75**(4), 046706 (2007). [CrossRef] [PubMed]

**12. **M. Xu and L. V. Wang, “Universal back-projection algorithm for photoacoustic computed tomography,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. **71**(1), 016706 (2005). [CrossRef] [PubMed]

**13. **P. Burgholzer, J. Bauer-Marschallinger, H. Grün, M. Haltmeier, and G. Paltauf, “Temporal back-projection algorithms for photoacoustic tomography with integrating line detectors,” Inverse Probl. **23**(6), S65–S80 (2007). [CrossRef]

**14. **M. Haltmeier, “Universal inversion formulas for recovering a function from spherical means,” SIAM J. Math. Anal. **46**(1), 214–232 (2014). [CrossRef]

**15. **L. A. Kunyansky, “Explicit inversion formulae for the spherical mean radon transform,” Inverse Probl. **23**(1), 373–383 (2007). [CrossRef]

**16. **Y. Xu, D. Feng, and L. V. Wang, “Exact frequency-domain reconstruction for thermoacoustic tomography--I: Planar geometry,” IEEE Trans. Med. Imaging **21**(7), 823–828 (2002). [CrossRef] [PubMed]

**17. **K. Wang and M. A. Anastasio, “A simple Fourier transform-based reconstruction formula for photoacoustic computed tomography with a circular or spherical measurement geometry,” Phys. Med. Biol. **57**(23), N493–N499 (2012). [CrossRef] [PubMed]

**18. **J. Xiao, X. Luo, K. Peng, and B. Wang, “Improved back-projection method for circular-scanning-based photoacoustic tomography with improved tangential resolution,” Appl. Opt. **56**(32), 8983–8990 (2017). [CrossRef] [PubMed]

**19. **H. Huang, G. Bustamante, R. Peterson, and J. Y. Ye, “An adaptive filtered back-projection for photoacoustic image reconstruction,” Med. Phys. **42**(5), 2169–2178 (2015). [CrossRef] [PubMed]

**20. **J. Turner, H. Estrada, M. Kneipp, and D. Razansky, “Improved optoacoustic microscopy through three-dimensional spatial impulse response synthetic aperture focusing technique,” Opt. Lett. **39**(12), 3390–3393 (2014). [CrossRef] [PubMed]

**21. **Y. Zhao, X. Zhao, and P. Zhang, “An extended algebraic reconstruction technique (E-ART) for dual spectral CT,” IEEE Trans. Med. Imaging **34**(3), 761–768 (2015). [CrossRef] [PubMed]

**22. **Z. Yu, J. B. Thibault, C. A. Bouman, K. D. Sauer, and J. Hsieh, “Fast model-based X-ray CT reconstruction using spatially nonhomogeneous ICD optimization,” IEEE Trans. Image Process. **20**(1), 161–175 (2011). [CrossRef] [PubMed]

**23. **A. Rosenthal, V. Ntziachristos, and D. Razansky, “Model-based optoacoustic inversion with arbitrary-shape detectors,” Med. Phys. **38**(7), 4285–4295 (2011). [CrossRef] [PubMed]

**24. **X. L. Deán-Ben, A. Buehler, V. Ntziachristos, and D. Razansky, “Accurate model-based reconstruction algorithm for three-dimensional optoacoustic tomography,” IEEE Trans. Med. Imaging **31**(10), 1922–1928 (2012). [CrossRef] [PubMed]

**25. **A. Rosenthal, D. Razansky, and V. Ntziachristos, “Fast semi-analytical model-based acoustic inversion for quantitative optoacoustic tomography,” IEEE Trans. Med. Imaging **29**(6), 1275–1285 (2010). [CrossRef] [PubMed]

**26. **K. Wang, R. Su, A. A. Oraevsky, and M. A. Anastasio, “Investigation of iterative image reconstruction in three-dimensional optoacoustic tomography,” Phys. Med. Biol. **57**(17), 5399–5423 (2012). [CrossRef] [PubMed]

**27. **M. Á. Araque Caballero, A. Rosenthal, J. Gateau, D. Razansky, and V. Ntziachristos, “Model-based optoacoustic imaging using focused detector scanning,” Opt. Lett. **37**(19), 4080–4082 (2012). [CrossRef] [PubMed]

**28. **G. Drozdov and A. Rosenthal, “Analysis of negatively focused ultrasound detectors in optoacoustic tomography,” IEEE Trans. Med. Imaging **36**(1), 301–309 (2017). [CrossRef] [PubMed]

**29. **T. Jetzfellner, A. Rosenthal, A. Buehler, A. Dima, K.-H. Englmeier, V. Ntziachristos, and D. Razansky, “Optoacoustic tomography with varying illumination and non-uniform detection patterns,” J. Opt. Soc. Am. A **27**(11), 2488–2495 (2010). [CrossRef] [PubMed]

**30. **X. L. Deán-Ben, R. Ma, A. Rosenthal, V. Ntziachristos, and D. Razansky, “Weighted model-based optoacoustic reconstruction in acoustic scattering media,” Phys. Med. Biol. **58**(16), 5555–5566 (2013). [CrossRef] [PubMed]

**31. **X. Jia, Y. Lou, R. Li, W. Y. Song, and S. B. Jiang, “GPU-based fast cone beam CT reconstruction from undersampled and noisy projection data via total variation,” Med. Phys. **37**(4), 1757–1760 (2010). [CrossRef] [PubMed]

**32. **J. Provost and F. Lesage, “The application of compressed sensing for photo-acoustic tomography,” IEEE Trans. Med. Imaging **28**(4), 585–594 (2009). [CrossRef] [PubMed]

**33. **A. Özbek, X. L. Deán-Ben, and D. Razansky, “Optoacoustic imaging at kilohertz volumetric frame rates,” Optica **5**(7), 857 (2018). [CrossRef]

**34. **L. Ding, X. L. Dean-Ben, and D. Razansky, “Efficient 3-D Model-Based Reconstruction Scheme for Arbitrary Optoacoustic Acquisition Geometries,” IEEE Trans. Med. Imaging **36**(9), 1858–1867 (2017). [CrossRef] [PubMed]

**35. **M. Omar, J. Gateau, and V. Ntziachristos, “Raster-scan optoacoustic mesoscopy in the 25-125 MHz range,” Opt. Lett. **38**(14), 2472–2474 (2013). [CrossRef] [PubMed]

**36. **X. L. Deán-Ben, T. F. Fehm, S. J. Ford, S. Gottschalk, and D. Razansky, “Spiral volumetric optoacoustic tomography visualizes multi-scale dynamics in mice,” Light Sci. Appl. **6**(4), e16247 (2017). [CrossRef] [PubMed]

**37. **J. Prakash, H. Dehghani, B. W. Pogue, and P. K. Yalavarthy, “Model-resolution-based basis pursuit deconvolution improves diffuse optical tomographic imaging,” IEEE Trans. Med. Imaging **33**(4), 891–901 (2014). [CrossRef] [PubMed]

**38. **J. Prakash, A. S. Raju, C. B. Shaw, M. Pramanik, and P. K. Yalavarthy, “Basis pursuit deconvolution for improving model-based reconstructed images in photoacoustic tomography,” Biomed. Opt. Express **5**(5), 1363–1377 (2014). [CrossRef] [PubMed]

**39. **M. S. Ždanov, *Geophysical Inverse Theory and Regularization Problems* (Elsevier, 2006).

**40. **C. C. Paige and M. A. Saunders, “LSQR: an algorithm for sparse linear equations and sparse least squares,” ACM Trans. Math. Softw. **8**(1), 43–71 (1982). [CrossRef]

**41. **C. B. Shaw, J. Prakash, M. Pramanik, and P. K. Yalavarthy, “Least squares QR-based decomposition provides an efficient way of computing optimal regularization parameter in photoacoustic tomography,” J. Biomed. Opt. **18**(8), 080501 (2013). [CrossRef] [PubMed]

**42. **A. Buehler, A. Rosenthal, T. Jetzfellner, A. Dima, D. Razansky, and V. Ntziachristos, “Model-based optoacoustic inversions with incomplete projection data,” Med. Phys. **38**(3), 1694–1704 (2011). [CrossRef] [PubMed]

**43. **D. Razansky, A. Buehler, and V. Ntziachristos, “Volumetric real-time multispectral optoacoustic tomography of biomarkers,” Nat. Protoc. **6**(8), 1121–1129 (2011). [CrossRef] [PubMed]

**44. **C. Zhang, C. Li, and L. V. Wang, “Fast and robust deconvolution-based image reconstruction for photoacoustic tomography in circular geometry: experimental validation,” IEEE Photonics J. **2**(1), 57–66 (2010). [CrossRef] [PubMed]

**45. **C. Lutzweiler, X. L. Deán-Ben, and D. Razansky, “Expediting model-based optoacoustic reconstructions with tomographic symmetries,” Med. Phys. **41**(1), 013302 (2013). [CrossRef] [PubMed]

**46. **A. Rosenthal, T. Jetzfellner, D. Razansky, and V. Ntziachristos, “Efficient framework for model-based tomographic image reconstruction using wavelet packets,” IEEE Trans. Med. Imaging **31**(7), 1346–1357 (2012). [CrossRef] [PubMed]

**47. **J. Aguirre, A. Giannoula, T. Minagawa, L. Funk, P. Turon, and T. Durduran, “A low memory cost model based reconstruction algorithm exploiting translational symmetry for photoacustic microscopy,” Biomed. Opt. Express **4**(12), 2813–2827 (2013). [CrossRef] [PubMed]

**48. **T. Taylor and L. R. Lupton, “Resolution, artifacts and the design of computed tomography systems,” Nucl. Instrum. Methods Phys. Res. A **242**(3), 603–609 (1986). [CrossRef]