Algebraic determination of back-projection operators for optoacoustic tomography

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


Introduction
Tomographic inversion problems, which pertain to the determination of 2D or 3D objects from their projections, are common in many imaging fields [1][2][3][4][5][6][7][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][2][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 backprojection (BP) algorithms [11][12][13][14][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][19][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 optimizationbased 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][23][24][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][28][29][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 nonlinear 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 realtime 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 BPbased 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.

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 t r at a given position r and time instant t to the optoacoustic image ( ) H r is often described by the following integral relation [10]: where, Γ is the Grüneisen coefficient, c is the speed of sound in the medium, and t is time. The projection formula in Eq. (1) represents 3D integration (or projection) over a sphere with a radius of ct whose origin is at r . In 2D problems in which it is assumed ( ) H r is approximately confined to a plane, the integral in Eq. (1) is calculated over a circle.
The set of acoustic signals ( , ) p t r obtained over several locations represents the projection data used for reconstructing the image ( ) H r . When ( ) H r has a finite support, and ( , ) p t r 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]: where ( ') r dΩ r is the solid angle element corresponding to the detector surface element ' dS when viewed from r , and is given by the following equation [11]: where ( ') n r denotes the outward-pointing normal of the integration surface at ' r . In Eq. (2), 0 Ω is the solid angle covered by the detection surface and is equal to 4π for spherical and cylindrical surfaces, and to 2π for planar surfaces. When the detection surface is in the far field of the imaged object, Eq. (2) represents an approximate solution to the inverse problem for any closed-surface [11]. In that case, the first term in Eq. (2) may be omitted since ( ', ) ( ', ). p t t t p t ∂ ∂ r r  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 and h 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 2 || || , ⋅ is the 2  norm. For a given signal p , the solution of Eq. (5) may be found by using the pseudo-inverse of M [37][38][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 by † rec h = M p (6) where † M is the pseudo-inverse matrix, which is given by ( )

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 . 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 i p , where the entire projection data are given by the following equation: ay be represent K the back-proje cordingly, n K ction. As illust ted images. A matrix relation matrix K is giv projection operato shown by the arc on is the same fo ted so it is alway overlap with the re square), where the ted as follows: Multiplying Eq. (9) by M and substituting Eq. (4), we obtain recˆ.

p = MKp
An exact reconstruction would fulfill rec p = p for any vector p . Using Eq. (11), a kernel matrix that achieves an exact reconstruction should fulfill ˆ,

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 temp 3) Set  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 s M . The computational burden of ABP is mostly in step 9, which involves solving an optimization problem for the matrix ( ) In contrast, step 11, which is individually performed for each dataset, involves a single matrix-matrix multiplication with a complexity of ( ) 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 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 where N λ is the number of regularization parameters that are tested.

Numerical results
Our numerical simulations were performed for the 2D geometry illustrated in Fig. 1 with the parameters 81 x y N N = = and 0.025 x y Δ = Δ = 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 325 p N = projections. We note that taking p N to be considerably larger than x N 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 647 . 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 / (2 ) c Δ , corresponding to 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 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 s M and v I 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 nonzero 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 reconstruction run times repo hen the LSQR a the results of a single itera ration of the L d accentuation In addition, th rations were pe with only minor [22]. The ilar to the s e = e g e a n a algorithm the 120 th ation and SQR was n of small he images erformed, r limited-

Figures 4( and Figs. 4(p calculated, th iterations wer (Section IIIB)
The first i in the first ite the algebraic algebraic reco were significa albeit more vi compared the the images in in the project andard-deviatio variables were quality of the re f the images the algebraic in d at the 2 nd , 8 th structions obta tion of Eq. (5) ns that produce l)). A 1D profi y = 0 is shown tained by ABP ns. Once the m 17 seconds, re teps 11 and 12 ely. ed from the sam gs. 4(g)-4(i)). S the ABP Simi ABP reconstr w artifacts rem econstructions from noisy data. Th a (BP) in Eq. (2 ber of iterations th 20 iterations in th e performance n Fig. 4(a)-4(c tion data a ra on of 20% of uncorrelated. econstructions when compar nversion, the r h , and 3 rd itera ained using the ) performed wi ed the lowest R ile of the recon n in Fig. 6.   Fig. 4(b). Fig. 4  Similar to distance of 0 obtained using a single micr reconstruction magnitude of reconstruction three methods suffered from

Conclusio
In conclusion operators. Th back-projectio is applied on kernel matrix operation foll given geomet numerical effi While the represents a n for the proble geometry sim projections an operations re Second, the programming the LSQR alg projection of field approxim iteration ABP obtained whe used to calcul maintained its In terms considerably iteration of th 8. (a-c) Reconstruc using (a) BP (b) alg terations. For the a e microsphere rec struction through sphere reconstruct the previous s 0.5 cm from g BP, algebraic rosphere, only n was perform f the microsphe n, obtained on s obtained a si m a significantly on n, we develope he method, wh on operator in each of the pr x K , is foun lowed by a su try and does no ficiency compa e purpose of th new type of so em of 2D OA mplified the imp nd the image w quired to des use of a 2D . In our numer gorithm. When the derivative mation to the 3 P in our 2D ge n the analytica late K, finer d s general struct of run time, i faster than ite he LSQR algor ction of a microsp gebraic inversion w algebraic inversion construction was o the horizontal li tions. section, the im the detector. c inversion, an y a 5 × 5 mm med with 5 iter ere in the imag the horizontal imilar reconstru y higher noise l ed a method f hich we term a the least-squar rojections. On nd, the recon ummation oper ot depend on th arable to the on his paper is to olution to tomo AT with a linea plementation o were discretize cribe the tran example enab rical simulation n LSQR was ap e of the acoust 3D BP formul eometry, we ob al BP formula details were ad ture of arcs wit in our numeri erative algebra rithm requires phere with a diam with 5 iterations an n, no further impr obtained with mo ine y = 3.72 mm mage was recon   [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 ( ) 3 O N 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 ( ) 5 O N per iteration, whereas the pseudo-inverse approach has a complexity of ( ) 6 O N . 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 leastsquares 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 spatialaveraging 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 s M 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 t N independent calculations in Step 9 in Section IIIB, an improvement on the scale of t N 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 precalculation offline. Alternatively, parallelization may be achieved by using graphical processing units (GPUs). Although the matrix s M 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 s M is not composed of unique entries, but rather contains p N 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 s M 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 K has a size of To determine 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 K The matrix v K in Eq. (A3) is connected to the matrix K in Eq. (8) by the following transformation:  N N N N N N N N N N N N N N N N N N N N N N N N N N N N N N N N N N N N N N N N N N N N N N N N N N N We denote the first matrix on the right hand side of Eq. (A4) by s I . As Eq. (A4) shows, the size of the matrix s I is ( 1 )  Kp . Each column in mat h represents a single backprojected image with the size of K N and y N in the x and y directions, respectively, as illustrated in Fig. 1(b). Denoting the nth column vector of mat h by mat,n h 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: ( )

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