1.

Introduction

The geometrical and operational characteristics of ultrasonic detectors play a key role in the reconstruction strategy and consequent image quality obtained in optoacoustic (photoacoustic) imaging. Specifically, the size of the active surface (aperture) and the frequency response of the detector(s) employed determine the achievable sensitivity and resolution of the optoacoustic system.¹ Commonly used piezoelectric transducers can be classified into three main groups, namely: flat, spherically focused, and cylindrically focused transducers. Each shape privileges certain detection geometries and thus determines the possible fields of application of the optoacoustic imaging system.

Spherically focused transducers, i.e., detectors with a concave-shaped surface, are used in acoustic resolution optoacoustic microscopy^2–5 to selectively collect time-domain signals providing depth-profiles of the optical absorption at each measuring position. While three-dimensional (3-D) images can be formed by stacking the signals acquired by raster-scanning the transducer, more sophisticated reconstruction procedures accounting for the shape of the detector can render a better image quality.^6–8

Flat detectors on the other hand are used for volumetric tomographic implementations, where optoacoustic signals at multiple distinct positions surrounding the to-be-imaged object are collected.^9–11

Although the frequency-dependent angular acceptance of flat detectors is much higher than in focused transducers, they are still directive due to diffraction, especially for the relatively large sizes required to provide an acceptable signal-to-noise ratio (SNR). Thereby, more advanced reconstruction algorithms have been suggested to take into account the shape of the transducer in the reconstructions, yielding significantly better results than standard reconstruction procedures assuming point detectors.^12,13

High-resolution volumetric optoacoustic imaging is, however, limited by the large acquisition time needed to collect a sufficient number of signals, even if transducer arrays are employed. In order to overcome this limitation, cross sectional optoacoustic tomographic systems have been suggested, which make use of cylindrically focused transducers to selectively collect signals originated in the imaging plane.^14–18 By simultaneously acquiring signals with a cylindrically focused array of transducers, real-time imaging performance has been showcased.^19,20 Furthermore, the high-throughput capacity of such cross sectional imaging systems, combined with the functional and molecular imaging capabilities of the optoacoustic technology, provide new biological insights in small animal research, with promising applications continuously emerging.^21–23 Finally, volumetric imaging can also be achieved by translating the ultrasonic array in the elevational direction.²⁴ However, although cylindrically focused transducers partially reject signals originating outside the imaging plane, the reconstructed images obtained by assuming a two-dimensional (2-D) geometry are generally strongly affected by out-of-plane artifacts. Also, the shape of the transducer influences the in-plane resolution and the quantification performance of the imaging system. Then, the development of a reconstruction procedure to minimize the reconstruction inaccuracies due to the cylindrical shape of the detector becomes an important issue to address in widely used cross sectional optoacoustic imaging systems.

In this work, we investigate the benefits of using a 3-D model-based reconstruction procedure accounting for the shape of the transducer in tomographic optoacoustic systems based on cylindrically focused transducers. First, the effects of the geometrical features of the detectors on the resolution and quantification performance of 2-D model-based reconstructions assuming point detectors are examined. Then, the 3-D model-based methodology is introduced and its reconstruction performance is assessed. Specifically, a 3-D forward model based on an accurate discretization of the Poisson-type integral²⁵ is modified to incorporate the spatial impulse response (SIR) due to the geometrical characteristics of cylindrically focused transducers. Two approaches are suggested. The first one is based on convolving the pressure field at the measuring locations with the spatially dependent SIR of the transducer, as described in Ref. 26 for the 2-D case. The second one consists of discretizing the transducer surface by a set of surface elements and superimposing the corresponding optoacoustic signals. The reconstructed images obtained with 2-D model-based inversion are then contrasted with the respective images obtained by using the full 3-D model. Theoretical predictions are compared with experimental results and the findings and conclusions are finally discussed.

2.

Theory

Under assumption of thermal and stress confinement in a uniform acoustic medium, the optoacoustic wave generation is mathematically described by

An explicit expression for the solution is given by the Poisson type integral

In some cases, the optoacoustic problem is effectively reduced to 2-D by selectively collecting signals generated in the imaging plane. In such case, the surface integral in Eq. (5) is approximately expressed as a line integral multiplied by the height of the imaging slice, i.e., Eq. (5) (in arbitrary units) is reduced to

Model-based reconstruction in 2-D reduced systems consists of a numerical discretization followed by an algebraic inversion of Eq. (6), implemented herein with the interpolated model-matrix inversion algorithm (IMMI).²⁸ Discretization based on approximating the line integral in Eq. (6) by equally spaced points has been considered,²⁹ which can also be used for 3-D reconstructions by discretizing Eq. (5).²⁵ In all cases, the pressure at a given point and at a given instant is expressed as a linear combination of the absorption in the pixels (2-D case) or voxels (3-D case) of the reconstructed region of interest (ROI). Thereby, the pressure vector $\tilde{P}$ corresponding to the pressure at a set of positions and instants is given by

In some cases, the inversion procedure needs to be regularized. For example, by employing standard Tikhonov regularization, Eq. (8) is transformed to

Image reconstruction with the model matrix $M$ assumes that the collected signals correspond to the pressure at a given point position (point detector assumption). This hypothesis is not always valid as the acoustic pressure is spatially averaged on the active area of realistic ultrasonic transducers which have a certain extension. In this way, the signal detected by a transducer having a detection surface $S$ can be expressed as

The pressure $P(x_S,t)$ for a set of instants and for a given surface element can be expressed by Eq. (7). Thereby, by combining Eqs. (7) and (11) a model matrix $M_{\mathrm{sum}}$ incorporating the effects of the finite size of the transducer is given by

If the size of the surface elements is the same, the term $\mathrm{\Delta }{x}_S$ in Eqs. (11) and (12) can be dropped for simplicity. The accuracy of this procedure depends on the number of elements $\mathrm{\Delta }{x}_S$ used to approximate the detector shape. Therefore, it can be slow or even impracticable for large detectors. On the other hand, it provides flexibility as any detector shape can be modeled by a set of surface elements.

The second approach considered in this paper is based on calculating the spatial impulse response (SIR) of the transducer.²⁶ The analytical solution of the wave equation for a spatial and temporal Dirac-delta source term, i.e., substituting the right-hand side of Eq. (1) by $\delta (x^{\prime })\delta (t)$, in free space is given by

The integration of Eq. (14) along the active area corresponds to the spatially dependent impulse response of the transducer due to its shape. The SIR of a finite length line transducer, denoted by $h(x,x^{\prime },t)$, can be calculated analytically.²⁶ Here, $t$ is time, $x$ location of the center of the line (detector), and $x^{\prime }$ position of a point source. Then, the cylindrically focused detector can be approximated by $n$ lines (centered at $x_1,\dots ,x_n$) and its impulse response $h_{\det }(x_{\mathrm{c}},x^{\prime },t)$ can be estimated as

Fig. 1

Comparison of pressure signals emitted by four absorbers of different size. (a and d) Middle plane of the region of interest (ROI) with four parabolic absorbers of radius 200 μm and 1 mm placed along the $x$-axis. (b and e) Signals predicted by the model matrices $M_{\det }$ (red continuous) and $M_{\mathrm{sum}}$ (green dotted) together with the analytical signals (black dotted) for a transducer positioned at $(x|y|z)=(2.54\mathrm{cm}|0\mathrm{cm}|0\mathrm{cm})$. (c and f) Fourier transform of the signals in (b) and (e), respectively. (g) Cross section of a mouse. (h and i) The signals predicted by the two matrices and their Fourier transform. (j and k) Approximation of the active surface of the transducer by a set of points and lines. A picture of the actual transducer is shown in (l).

Again combining Eqs. (7) and (16), a numerically efficient convolution of the SIR $h_{\det }$ with the model-matrix $M$ leads to a new matrix $M_{\det }$ (which is equivalent to $M_{\mathrm{sum}}$). This yields a new system for inversion

This latest approach based on convolution may be more suitable for cylindrically focused detector surfaces with a relatively large width as the number of lines needed to approximate the transducer surface is much lower than the number of points required in the first approach.

3.

Simulations

Numerical simulations were employed to examine the accuracy of both models described above. An optical absorption distribution $f(x,y,z)$ was assumed by a 3-D truncated paraboloid with radius $r_0$, i.e.,

For this absorption pattern, the laser-induced pressure wave detected at a point in space can be calculated analytically (see Appendix A). Then, the surface of a finite-size detector was discretized uniformly and very densely with $>2000$ equidistant points, and the analytical signals corresponding to each point were summed by Eq. (11) to determine the response of the entire detector area. Henceforth, we refer to the signals calculated in this manner as the analytical signals.

A cylindrically focused transducer with a circular shape was considered for the simulations (see Fig. 1). The diameter of the circle is 1.3 cm and its focal length is 2.54 cm. The analytical signals were compared to the ones predicted by the two models introduced in the previous section. For this purpose, the pressure vector was calculated from the matrices $M_{\det }$ and $M_{\mathrm{sum}}$ using Eq. (7). $M_{\det }$ was calculated by approximating the surface of the transducer with 21 lines. No significant changes are produced in $M_{\det }$ for a high number of lines. $M_{\mathrm{sum}}$ was calculated by considering 350 equally spaced points on the transducer surface, so that the distance between the points for calculating $M_{\mathrm{sum}}$ is approximately the same as the distance between the lines for calculating $M_{\det }$. The discrete ROI considered consists of $101\times 101\times 21$ voxels, equivalent to $2\times 2\times 0.4{\mathrm{cm}}^3$, resulting in a uniform resolution of 200 μm.

Two different optical absorption distributions containing four absorbers with radii $r_0=200\mu \mathrm{m}$ and $r_0=1\mathrm{mm}$ were considered [Fig. 1(a) and 1(d)]. The signals predicted by the two models show good agreement with the analytical ones ensuring that the models are accurate [Fig. 1(b) and 1(e)]. The signals shown in Fig. 1 correspond to a transducer positioned at $(x|y|z)=(2.54\mathrm{cm}|0\mathrm{cm}|0\mathrm{cm})$. For all the other positions on the detector ring, analytical signals and signals predicted by the models matched equally. As absorbers with different sizes emit pressure waves at different frequencies, it is important that the whole frequency spectrum of broadband signals is accurately calculated. The Fourier transform of the signals [Fig. 1(c) and 1(f)] showcases a good accuracy of the two models throughout the frequency spectrum. There is, however, a discrepancy between the analytical signals and the signals predicted by the models for absorbers in the order of the pixel size, which is mainly due to discretization errors produced by approximating the actual absorption distribution with the discrete ROI. Furthermore, to showcase the equivalent behavior of the two models, a stack of cross sectional images of a mouse [Fig. 1(g)] was considered and the signals predicted by the two matrices were compared. Also here, the signals are almost identical either in time [Fig. 1(h)] or in frequency [Fig. 1(i)] domain. As the two models show the same behavior, the model-matrix $M_{\det }$ was considered in the rest of simulations and experiments due to the lower computational time for wide detectors. Building a model matrix for the entire detection geometry assuming point detectors took about 130 s. Therefore, building the entire matrix $M_{\mathrm{sum}}$ (which is the sum of 350 such point detector models) took over 12 h and required 12.8 GB of memory for storage. Computation of the model–matrix $M_{\det }$ comprised three steps. The SIR $h_{\det }$ and one 3-D model-matrix assuming point transducers had to be calculated and finally convolved. These steps in total took roughly 8 h. Memory requirements for storage of $M_{\det }$ were the same as for $M_{\mathrm{sum}}$.

In a second step, the analytical signals were calculated for a tomographic geometry (Fig. 2). For this, we consider a cylindrically focused transducer scanned along a circumference surrounding the object with an 2.25 deg step (160 projections) and then, additionally scanned linearly along 0.6 cm in the longitudinal direction (31 steps), so that 4960 transducer positions are taken. In this way, we intend to compare the results obtained by reducing the reconstruction problem to 2-D with those obtained by a 3-D reconstruction with modeling the shape of the transducer. In order to show the effects in absorbers of different sizes, we considered the same two absorption distributions depicted in Fig. 1 corresponding to absorbers with radii $r_0=200\mu \mathrm{m}$ and $r_0=1\mathrm{mm}$.

Fig. 2

Full-view tomographic geometry for the simulations and experiments. The ROI is depicted by the red cuboid. The blue points represent the positions of the centers of the cylindrically focused detectors. All the transducer positions lie on the surface of a cylinder with radius 2.54 cm.

The inversion was performed using two alternative methods. First, a 2-D model matrix representing the 160 detector positions in one plane was calculated. Thereby, each plane was reconstructed separately with this matrix. The inversion was done by means of the LSQR algorithm, resulting in a stack of 2-D reconstructions representing the volumetric ROI. Here, no Tikhonov regularization was employed as the LSQR algorithm converges for full-view acquisition in the 2-D case, i.e., the optimal regularization parameter for Tikhonov regularization [Eq. (9)] is $\lambda =0$. Second, the full 3-D model matrix incorporating the SIR of the transducer was calculated. The inversion in this case was done with the LSQR algorithm and standard Tikhonov regularization.

Four point absorbers with radius $r_0=200\mu \mathrm{m}$ [Fig. 3(a) and 3(d)] and radius $r_0=1\mathrm{mm}$ [Fig. 3(g) and 3(j)] were placed in the middle plane, starting from the center going outward. The reconstruction achieved by inverting a 2-D model for each plane shows the expected smearing in the imaging plane of the absorbers located away from the center of the image (in-plane artifacts) [Fig. 3(b) and 3(h)]. This is due to the assumed point transducers in the model which do not correspond to the actual signals collected by the transducer. Also, the reconstructed absorption values are severely reduced for peripheral absorbers, resulting in significant quantification errors. Smearing over almost the entire ROI is observed in the $z$-direction, corresponding to strong out-of-plane artifacts [Fig. 3(e) and 3(k)].

Fig. 3

Simulation of four absorbers with radius $r_0=200\mu \mathrm{m}$ (a and d) and radius $r_0=1\mathrm{mm}$ (g and j) placed along the $x$-axis in the center of the ROI. (b and h) The maximum intensity projection (MIP) along the $z$-axis of the stack of 2-D reconstructions and (e) and (k) its MIP along the $y$-axis. (c) and (i) depict the MIP along the $z$-axis of the 3-D reconstruction taking the spatial impulse response (SIR) of the transducer into account and its MIP along the $y$-axis (f and l). Absorption values (normalized) along the $x$-axis in the middle plane for four absorbers with radius $r_0=200\mu \mathrm{m}$ (m) and $r_0=1\mathrm{mm}$ (p). Relative improvement of the absorption values in the $z$-direction is shown for the central absorber (n and q) and the outmost absorber (o and r).

On the other hand, the reconstructions retrieved with the full 3-D model including the SIR of the detector significantly reduce the smearing in the out-of-plane [Fig. 3(f) and 3(l)] and in-plane directions [Fig. 3(c) and 3(i)]. Resolution in all spatial dimensions was improved with the 3-D model incorporating the SIR of the transducer.

Also, it is shown [Fig. 3(m) and 3(p)] that the error in the reconstructed absorption value is size dependent. Figure 3(m) and 3(p) shows the horizontal profiles of the four absorbers along the $x$-axis in the middle plane. If 2-D reconstruction is used, the values of the ratio between the retrieved amplitudes for the outer and the inner absorbers are $\sim 10\%$ for the small absorbers and 40% for the big absorbers, respectively. When the full 3-D model is used, these values are 40% and 75%, respectively, which indicates that the quantitative errors are considerably reduced with the latest approach.

It is important to notice that the overall improvement of the reconstruction was also influenced by the length of the scan in the $z$-direction. Longer scans led to a better conditioned model matrix and thereby inversion problem so that better quality reconstructions were generally obtained (see Appendix B).

4.

Experimental Results

The procedure suggested in this work was also tested in experiments with agar phantoms containing microparticles and an ex vivo spleen of a mouse. The cylindrical phantoms with a diameter of 2 cm were prepared using a gel made from distilled water, containing Agar (Sigma-Aldrich, St. Louis, MO, USA) for jellification (1.3% $\mathrm{w}/\mathrm{w}$) and an Intralipid 20% emulsion (Sigma-Aldrich) for light diffusion and more uniform illumination (6% $\mathrm{v}/\mathrm{v}$), resulting in a gel presenting a reduced scattering coefficient of ${\mu }_s^{\prime }\approx 10{\mathrm{cm}}^{-1}$. First, polyethylene microparticles with an approximate diameter of 200 μm (Cospheric BKPMS 180 to 210 μm) have been placed in one plane trying to resemble the setup of the simulations as closely as possible [Figs. 3(a), 3(d) and 4(a), 4(d)]. Additionally, the spleen of a mouse has been embedded in a second agar phantom.

Fig. 4

Positioning of four microspheres in the experiment; MIP along the $z$-axis (a) and along the $y$-axis (d). (b) and (e) The corresponding reconstruction via inversion of the 2-D model-matrix system. (c) The reconstruction obtained via inversion of the 3-D model matrix with the SIR convolved and its MIP along the $y$-axis (f). In (g), the absorption values of the four microspheres projected along the $y$-axis for both reconstructions are shown. The relative improvement of the absorption values in the $z$-direction is shown for the central absorber (h) and the outmost absorber in (i). (j) The position of the four absorbers in the middle plane with a threshold from 0 to 0.1. (k) The middle plane of the 2-D reconstruction with the same threshold. (l) The central plane of the reconstruction obtained by the full 3-D model with the same threshold from 0 to 0.1.

The layout of the tomographic optoacoustic system employed to image the phantom containing microparticles is depicted in Fig. 5. The optoacoustic signals were measured by a standard cylindrically focused piezoelectric immersion transducer (Panametrics V320-SU, Olympus NDT Inc., Waltham, MA, USA) with the same dimensions as in the simulations, i.e., a diameter of 1.3 cm and a focal length of 2.54 cm. The central frequency of the transducer is 7.5 MHz with a bandwidth of $\sim 70\%$. Optoacoustic pressure waves were excited with a tunable (680 to 950 nm) optical parametric oscillator laser (Phocus, Opotek Inc., Carlsbad, CA, USA), delivering $<10\mathrm{ns}$ duration pulses with repetition frequency of 10 Hz. The laser was set at a wavelength of 760 nm (corresponding to the maximum laser power) and the acquired signals at each projection were averaged 10 times and bandpass filtered from 0.5 to 12 MHz.

Fig. 5

Layout of the experimental setup.

The laser beam was guided with a fiber bundle into the water tank and (due to spatial limitations) deflected with a mirror in order to uniformly illuminate the phantom from the bottom. The transducer was kept at a fixed position while the sample was rotated and moved vertically by means of rotation and translation stages. The Q-switch output of the laser was used to trigger the data acquisition card (Spectrum, M3i.4121) embedded in the personal computer controlling the stages.

The transducer was positioned at a distance from the center of rotation equal to its focal length, namely, 2.54 cm. In order to cover the entire volumetric ROI, the phantom was scanned in 31 positions along the $z$-axis with a step of 200 μm covering a total distance of 6 mm (Fig. 2). For each vertical position of the sample, 160 projections over 360 deg were recorded (angular step of 2.25 deg).

Much as in the simulations, the ROI was a volume with dimensions of $2\times 2\times 0.4{\mathrm{cm}}^3$ discretized by $101\times 101\times 21$ voxels resulting in a uniform resolution of 200 μm. It was centered inside the detector ring and the scanning range. The reconstructions with the 2-D model matrix with no impulse response and with the 3-D model matrix including the SIR of the transducer were computed on a workstation computer $2\times$ Intel Xeon DP X5650 ($6\times 2.67\mathrm{GHz}$) with 144 GB of random-access memory (RAM). The LSQR algorithm was executed with MATLAB (Mathworks, Natick, MA, USA). For inversion in the 2-D case, the LSQR algorithm shows convergence, so no Tikhonov regularization was employed (or equivalently, Tikhonov regularization was used with $\lambda =0$). The 3-D reconstruction was always performed by means of the LSQR algorithm with standard Tikhonov regularization by optimizing the value of $\lambda$ to give the best possible image quality.

Figure 4 displays the reconstruction results. The stack of 2-D reconstructions shows the effects showcased in the simulations. The reconstructed microparticles are elongated in the azimuthal direction with an increasing distance from the center of rotation. The absorbers also become smeared throughout almost the entire ROI in the elevational direction [Fig. 4(b) and 4(e)], resulting in strong out-of-plane artifacts.

The 3-D reconstruction on the other hand shows the positive effect of the incorporation of the impulse response. In-plane, the absorbers are clearly visible showing four distinct peaks and its shape is almost perfectly recovered. Out-of-plane, their extension is also reduced to make their localization much more precise [Fig. 4(c) and 4(f)]. For the central and the outmost absorber, the resolution improvement in the $z$-direction is depicted in Fig. 4(h) and 4(i), respectively. Moreover, the overall signal-to-noise level is higher in the 3-D reconstruction [Fig. 4(g)]. Figure 4(k) depicts the middle plane of the stack of 2-D reconstructions with a threshold between 0 and 0.1 in order to make the noise floor more visible. Clearly, the same plane obtained via the full 3-D reconstruction and the same threshold shows visibly less noise [Fig. 4(l)]. As a measure for the noise level, the standard deviation of the reconstructed pixel values excluding the absorber region has been calculated ($y$-values between $[-1,-0.5]$ and [0.5, 1]). For the 2-D model reconstruction the value was 1.1655 whereas in the case of the full 3-D model reconstruction the calculation yielded 0.5373 corroborating the assertion of a lower noise level in the 3-D reconstruction.

The agar phantom containing the spleen of a mouse was imaged with a high-throughput optoacoustic tomographic system as described in detail in Ref. 19. In this case, the signals were collected with a transducer array consisting of 64 cylindrically focused elements covering 172 deg. The size of the elements was approximately 2 and 15 mm in the azimuthal and elevational directions, respectively. Due to the small width of the elements as compared to their height, the main effect in the reconstructions was the out-of-plane spreading of the absorbers. Thereby, the modeling of the transducer in this case was simplified by discretizing it with 150 surface elements in the elevational direction and neglecting azimuthal extensions.

A comparison of the results obtained with the 2-D model and the 3-D simplified model is displayed in Fig. 6. An improvement in the elevational resolution is obtained with the 3-D model as shown in the maximum intensity projection (MIP) along the $y$-direction [Fig. 6(c) and 6(d)]. The out-of-plane artifacts are especially significant for the background absorption, which generates mainly low-frequency acoustic waves, and the focusing capacity of the transducer is lower in this case. The reduction of the out-of-plane artifacts corresponding to low-spatial frequencies also improves the visual quality of the MIP along the $z$-direction. The weak SNR observed in the 2-D reconstruction [Fig. 6(a)] was improved in the 3-D reconstruction [Fig. 6(b)]. This could be observed even when Tikhonov regularization was also employed for the 2-D reconstructions as the signal-to-noise level did not change significantly for different values of the regularization parameter $\lambda$.

Fig. 6

Reconstruction of a mouse spleen. (a and c) Maximum intensity projections obtained via a 2-D model. (b and d) The corresponding reconstructions from a 3-D model with the detector properties incorporated.

5.

Discussion

In this work, we adapted a model-based reconstruction algorithm to account for the effects of cylindrically focused transducers in commonly used cross sectional tomographic optoacoustic systems. Two different approaches based, respectively, on approximating the active surface of the transducer by a set of surface elements and by considering the SIR of a discretization of the transducer by lines were presented and compared. It was shown that the signals predicted by both models are equivalent, so that application of any of them yields a 3-D model matrix incorporating the effects of the transducer shape. Then, optoacoustic reconstruction was performed by inverting such a model matrix with the LSQR algorithm.

The reconstruction performance of the model-matrix $M_{\mathrm{sum}}$ depends strongly on the number of surface elements used to discretize the transducer surface. As the time required for calculating $M_{\mathrm{sum}}$ is proportional to the number of elements used, it can lead to very large computational times for wide transducers being modeled with uniformly spaced surface elements. In contrast to that, a relatively small number of lines is needed to calculate the model-matrix $M_{\det }$ for flat detectors. Typically, $<30$ lines are sufficient to accurately model the optoacoustic wave detection, whereas the number of elements needed to obtain the equivalent $M_{\mathrm{sum}}$ is in the order of several hundreds to thousands. However, the model approximating the transducer area by surface elements may be more convenient for narrow transducers where the effect of the detector breadth is negligible. Beyond that, the discretization of the active area of the transducer with a set of surface elements provides more flexibility to model an arbitrary curved surface.

It was shown in the simulations that the in-plane smearing observed in 2-D reconstructions assuming point transducers can be corrected when employing the complete 3-D model, even for absorbers located at a relatively large distance from the center of the ROI. Likewise, out-of-plane smearing was also reduced, although the out-of-plane resolution is limited by the acquisition geometry, which is not optimal for 3-D imaging. Generally, the smearing effect when using 2-D reconstruction is frequency dependent. In-plane, the effect is lower for large absorbers, since signals emitted by large objects have longer wavelengths (low frequencies) and are, therefore, less distorted by the SIR of the detector, which performs as a low-pass filter. Out-of-plane, spreading of the absorbers is, however, stronger for large absorbers as the wavelength-dependent focusing capacity of the transducer is higher for short wavelengths (high frequencies). Both effects could be significantly reduced, improving the resolution in all spatial dimensions, by using the 3-D model including the SIR of the transducer.

Much like the smearing effect, degradation of the reconstructed absorption values was also shown to be dependent on the size of the absorbers and therefore, on the frequencies of the emitted ultrasonic waves. Smaller absorbers emit ultrasonic waves with higher frequencies and get corrupted more strongly in the reconstruction. Such abatement leads to significant quantification errors in the 2-D reconstructions assuming point transducers, which are substantially reduced when considering the full 3-D model. The maximum achievable improvement, however, is also dependent on the scanning length in the $z$-direction, as this influences the conditioning of the inversion (see Appendix B). Also, as recently shown in Ref. 13, different regularization techniques affect the 3-D reconstruction performance.

The improvement achieved with the full 3-D model shown in the simulations was also confirmed experimentally with microparticles. Spreading of the absorbers in azimuthal and elevational directions could clearly be seen in the 2-D reconstructions. The 3-D reconstruction algorithm including the SIR of the detector showed improvements in all spatial dimensions. Moreover, the noise level in the images was generally reduced with the 3-D model.

The positive results obtained in a biological specimen (mouse spleen) demonstrate the applicability and convenience of the method to reconstruct actual biological tissues. The blurry images obtained by 2-D reconstruction could be significantly enhanced by using the full 3-D model. As out-of-plane artifacts and blurring are mainly due to background absorption, they correspond to low-frequency optoacoustic waves. The focusing capacity of the transducer, however, is weak in that frequency range. Therefore, the full 3-D model enhanced significantly the reconstruction quality. The simplification of the model neglecting the azimuthal extension of the transducer also led to acceleration of the inversion due to the corresponding sparser matrix obtained.

The main disadvantage of the method proposed is its memory requirements, as all algebraic reconstruction methods, also the IMMI algorithm employed here, become very time- and memory-consuming with an increasing resolution. Extending it to 3-D only aggravates the problem. Both models presented herein, $M_{\det }$ and $M_{\mathrm{sum}}$, needed each 12.8 GB of memory for storage. Calculation time for building the matrices was $\sim 8\mathrm{h}$ for $M_{\det }$ and $>12\mathrm{h}$ for $M_{\mathrm{sum}}$. Depending on the size of the considered detector surface the number of nonzero elements of the model matrix increases significantly. For transducers typically used in optoacoustics, a 10-fold increase in nonzero elements is produced with respect to the model matrix for point detectors. Then, storage of the model matrix is an important issue to address in order for the method to be a viable tool in algebraic image reconstruction. This can be partially alleviated by simplifications of the model (as done herein for imaging the spleen). However, larger matrices inevitably require the development of other strategies such as computation of the action of the model matrix in each step of the iterative reconstruction procedure or matrix-compression schemes. These issues will be addressed in future work.

Overall, the showcased good performance of the methods anticipates its convenience in those cases where high-resolution and quantitative optoacoustic tomographic imaging is wanted. Typical artifacts in optoacoustic tomography resulting from data collected by commonly used cylindrically focused ultrasonic transducers could be efficiently reduced, leading to a consequent improvement in the spatial resolution of the imaging system.