Optoacoustic Imaging and Tomography: Reconstruction Approaches and Outstanding Challenges in Image Performance and Quantification

1Nowadays, the terms optoacoustic and photoacoustic are equally used to describe the effect of acoustic wave generation by transient light absorption. Optoacoustic sensing and imaging draws its roots from the discovery of the photophone by Bell and his assistant Tainter in 1880, which the inventors used as the first practical wireless telephony or, in fact, optical communication device. Bell also named the effect the photophonic phenomenon and later suggested the so-called spectrophone. Following its discovery in solid and liquid media, the effect had been further confirmed by Tyndall and Rontgen to occur in gases. However, despite multiple early attempts to put the photophone into practical use, mainly for military communications, it was not until 1938 that Veingerov developed the first widely accepted implementation of the phenomenon, which he called method for gas analysis based on Tyndall-Rontgen optic-acoustic effect. The term photoacoustic was later adopted in 1970s when sensitive spectroscopy of solids was first demonstrated. Rosencwaig was the first to suggest the phenomenon for biological sensing. However, due to the lack of appropriate laser sources and ultrasound detection technologies, theoretical understanding of the underlying phenomena as well as backlogs with the algorithmic and data processing capacities, utilization of optoacoustics for imaging of real biological tissues and organisms has only evolved in the last decade.

2The field of biological optoacoustic imaging is developing tremendously with new technical approaches and applications continuously emerging. By combining now commercially available pulsed laser technology in the nano-second range and sensitive acoustic detectors, it was shown possible to generate opto-acoustic responses from tissue that carry significant spatially-resolved biochemical information. Due to its hybrid nature, i.e., optical excitation and acoustic detection, the optoacoustic imaging technology benefits both from the rich and versatile optical contrast and high (diffraction-limited) spatial resolution associated with low-scattering nature of ultrasonic wave propagation as compared to photon propagation.

3In vascularized tissues, highly absorbing hemoglobin manifests good optoacoustic contrast. It is therefore natural for optoacoustics to attain high fidelity images of vascular anatomy, dynamic microcirculation, tumor neovascularization, as well as blood oxygenation levels deep within highly diffuse tissues without introduction of contrast agents. Besides blood-related contrast, optoacoustics is sensitive to other intrinsic tissue contrast. Nevertheless, imaging of extrinsic photo-absorbing agents may offer important contrast advantage and imaging specificity but would typically require differentiation of these agents on top of spectrally varying background tissue absorption. In response, the so-called multi-spectral optoacoustic tomography (MSOT) technique relies on the spectral identification of chromophoric molecules and particles distributed in tissue over background tissue absorption. Pulses of different wavelengths are used, in a time-shared fashion, whereas the wavelengths are selected to sample a certain spectral characteristic in the absorption spectrum of intrinsic bio-markers and reporter agents of interest. Essentially, when operated at multiple wavelengths, optoacoustics is capable of resolving spectral signatures of chromophores in three dimensions (volumetrically), i.e., it can sense color, or more precisely, listen to color.

4Despite its great promise, like ultrasound, optoacoustic imaging possesses several limitations related to imaging through acoustically mismatched areas, such as lungs and bones. Furthermore, even though spatial resolution performance here does not directly suffer from light scattering in tissue, it is still affected by penetration limits and signal degradation due to light attenuation. Therefore, whole-body imaging with optoacoustics is currently only possible in small animals whereas clinical applications are limited to relatively superficial, low attenuating or otherwise conveniently accessible areas of the body. Optoacoustic imaging is indeed a rapidly evolving area in biomedical imaging sciences. Its in-vivo use puts forth a number of challenging problems demanding intensive investigations, from imaging instrumentation, quantified reconstruction algorithms, spectral processing schemes, detection sensitivity, and other technical issues, to biology-related topics, such as effectiveness of imaging contrast approaches or animal handling. We refer an interested reader to some of the recently published review articles summarizing on progress of optoacoustic imaging applications. On the technical side, while some focused reviews address the mathematical inverse acoustic problem and quantification challenges of multispectral optoacoustic reconstruction methods, the current paper attempts to comprehensively cover the most recent experimentally-driven algorithmic developments in the optoacoustic field and also review some newly introduced image acquisition methodologies and selected practical in vivo imaging approaches.

5Optoacoustic (OA) imaging is based on absorption of light radiation in tissue and conversion of the deposited energy into heat (Figure 1(a)). This results in an expansion of the tissue and mechanical stress which propagates in the form of pressure waves. Typically, light in the visible or near-infrared spectrum (400 nm–1,200 nm) is used for excitation of OA responses due to the relatively weak absorption of biological tissues in this spectral region, also known as the optical window. Alternatively, tissue can be excited with energy in the radiofrequency and microwave spectra, which is however not included in the scope of this paper. A large variety of chromophores absorb at the optical wavelengths of interest, which further leads to a high contrast between different tissues with varying chromophore concentration. For biomedical applications, examples of usable contrast include intrinsic chromophores (e.g., hemoglobin in its oxygenated and deoxygenated form, melanin, fat), extrinsically administered agents (e.g., nanoparticles, fluorophores) or genetically encoded markers, such as fluorescent proteins. Vast majority of the modern OA imaging setups utilize lasers that emit ultra-short pulses, mainly due to a better signal to noise performance and parallelization capacity. In this paper, we therefore limit our discussion to the theory and applications of pulsed OA sources, although the interested reader is indeed encouraged to explore more about OA imaging using modulated continuous-wave sources.

6To simplify the model of the induced pressure variations, two assumptions are usually made on the duration x of the laser pulse: the first is the so-called thermal confinement assuming that the heat conductance during a laser pulse in one image voxel does not affect the neighboring voxels (τ ≪ Δl/4Dt, where Δl is the dimension of the voxels and Dt is the thermal diffusivity). The second is the so-called stress confinement assuming the volume expansion during the laser pulse to be negligible (τ ≪ Δl/c, where c is the speed of sound). With these simplifications, the initially induced pressure p0(r) is generally proportional to the total absorbed optical energy density H(r): where Γ is the dimensionless Grüneisen parameter, β is the isobaric thermal expansion coefficient, and Cp is the isobaric specific heat capacity. In a more general sense, Γ represents the combined OA efficiency of conversion from heat into pressure. In contrast to significant spatial variations in the absorption coefficient, Γ does not normally exhibit dramatic variations among different soft tissues, thus it is most conveniently assumed to be constant. Nevertheless, accounting for spatial variations in Γ may still benefit the efforts on quantification of optoacoustic images. Moreover, some studies attempted utilizing the strong dependence of the Griineisen parameter on temperature in order to sense temperature optoacoustically.

7For biomedical applications, it is the distribution of the absorption coefficient μa(r) that is usually directly related to the concentration ci(r) of the chromophores rather than the totally absorbed energy H(r). Both are related via the expression: where Φ(r,μa,μs,g) is the light fluence (or photon density), μs is the scattering coefficient and g is the anisotropy factor. Although it might seem that H(r) is a simple product, it depends nonlinearly on the absorption coefficient μa because the light fluence generally depends on the underlying optical properties, among them μa itself. Thus, it is important for image quantification purposes to account for the light fluence since it may considerably vary as a function of depth. Figure 1(b) summarizes the main contributions to the OA signal generation, which need to be accounted for when developing accurate image reconstruction algorithms.

8The initial goal of optoacoustic imaging is to retrieve the initial pressure distribution inside the object due to the absorbed laser energy. However, the generated pressure fields can normally only be measured outside the object. Propagation of pressure toward the detection point is described by the optoacoustic wave equation, which, for a non-absorbing homogeneous medium, is written as:

9For tomographic data acquisition, the detectors are placed on a (closed) surface S surrounding the volume of interest. The most common detection patterns are a spherical, a cylindrical and a planar surface or their 2-D counterparts (Figure 2(a)). On the surface the temporal pressure profile pd(rd, t), rd ∈ S, t > 0, is considered known. To reconstruct an OA image from pd(rd, t), the initial pressure p0(r) satisfying: has to be found. To solve this partial differential equation, one usually takes advantage of the free-space Green's function G(r, r′, t, t′). It relates the pressure at two different spatio-temporal points, i.e.,

10Consequently, the pressure profile pd(rd, t) on the detection surface S can be directly related to the initial pressure p0(r) via a Poisson-type integral:

11This type of equation is known as spherical Radon transform or spherical mean transform. This is because contributions of all sources that are located on the surface of a sphere with radius R = c · t = |rd − r| from the detector are integrated (Figure 2(b)).

12A typical optoacoustic setup consists of several key components, an example shown in Figure 1(c). In the pulsed excitation mode, the tissue is illuminated by a laser emitting monochromatic pulses of light with a typical duration of some nanoseconds. For deep tissue imaging applications, optical parametric oscillators are often used to provide a tunable wavelength in the spectrum of interest with pulse repetition rate in the order of a few tens of Hertz and per-pulse energies in the millijoule range. In optoacoustic microscopy and other superficial applications, where such high per-pulse energies are not required, other types of sources in the microjoule and nanojoule range are considered as well, including high repetition passive Q-switched and dye lasers, laser diodes and fiber lasers. For tomographic imaging, the optoacoustically-generated pressure profiles are subsequently captured with detectors surrounding the object whereas the acoustic coupling between object and detector is usually ensured by water or gel medium. In contrast to ultrasonic (US) imaging, OA signal amplitudes are relatively low while their spectral content is broad, spanning the frequencies from several tens of kHz up to a hundred MHz for micron-scale structures. Thus, high sensitivity, ultrawide bandwidth as well as good tomographic coverage are key for ensuring image quality. In addition, for acquisition of spatially resolved data with any type of detector, either a single detector is scanned around the object or, alternatively, multiple detection elements acquire data in parallel. The latter allows fast data acquisition, e.g., for reconstruction of images from a single laser shot with commercially available or tailored detection arrays. Figure 3 illustrates the broad variety of OA setups presented in the literature including OA raster scanning (Figure 3(a,c)), tomography (Figure 3(b,e,f)) and endoscopy systems Figure 3(d). The signals were detected optically with an interferometer for the setups in Figure 3(c,e), whereas in the other systems conventional piezoelectric detection elements were used.

13Two main categories of sensors are used to detect OA signals: The most commonly used type is based on piezoelectric elements, such as piezocomposite- or polymer-based films, which are also widely utilized in US imaging. This well developed technology possesses relatively low cost of implementation and can be readily applied for highly sensitive measurements. However their sensitivity scales with detector's size, making miniaturization challenging. Of the second type are optical resonators and interferometry detectors like microring resonators, Fabry-Perot interferometers, Mach-Zehnder interferometers, or fiber Bragg gratings, which are sensitive to changes in the length of the optical path induced by pressure waves. Interferometric detectors usually possess ultrawideband detection characteristics and further allow all-optical optoacoustic imaging, in some cases also greatly simplifying delivery of the excitation light to the imaged object through transparent optical components. Sensitivity of interferometric detectors is not directly dependent upon the physical size of the detection element, which simplifies miniaturized designs suitable for small scale imaging, such as intravascular or endoscopic applications. On the other hand, interferometric detection approaches are usually difficult to parallelize in order to achieve faster tomography while it is also not always possible to increase detection sensitivity by integration of signals over large areas. Recently, other detection technologies have been introduced and successfully applied for optoacoustic imaging, such as capacitive micromachined ultrasonic transducers (CMUT). The advantages and disadvantages of the various optoacoustic sensor approaches are summarized in Table 1.

14Following their detection, the optoacoustic waveforms are pre-amplified and digitized by a data acquisition system (DAQ). For tomographic data acquisition, a motion stage can be translated or rotated in order to alter the position of the detector/s with respect to the imaged object or vice versa in order to increase the amount of the available data points (projections) around the imaged area. A sufficiently high number of projections and a suitable detection pattern are crucial for good image quality. A PC controls the acquisition process, stores the data and is used to reconstruct and render the images. Prior to image formation, an additional signal processing step normally takes place that may include some of the following: averaging over multiple pulses to enhance signal-to-noise ratio (SNR), denoising the signals, deconvolving the signals from electrical frequency response of the detector, band pass filtering to remove noise and/or enhance visibility of certain spatial frequency components in the images.

15The first step to obtain OA images consists in reconstructing the initial pressure distribution, i.e., solving Equation (4) or inverting Equation (6). Since the early reconstruction approaches were introduced, plenty of research has been done on the problem in the last decade, driven from both theoretical interest and the need to improve images acquired from specific experimental setups.

16The inversion schemes are usually based on idealized assumptions of infinitely small point detectors located on a closed measurement surface and a non-absorbing homogeneous acoustic medium. Modifications for deviations from these idealized assumptions will be further presented below.

17An inversion formula was initially presented by Norton et al. for ultrasonic imaging. The basic idea is to represent the solution p0(r) as a convergent series: where Ψk(r) are eigenfuctions of the negative Laplacian −Δ in the domain of the measurement geometry and αk are the so-called Fourier coefficients. The coefficients αk can be calculated from pd(rd, t). Although it has been shown that series solutions exist for arbitrary closed detection surfaces, only the geometries with known eigenfunctions (sphere, cylinder, plane) are usually of practical experimental interest.

18For a planar detection geometry, this approach results in a fairly simple reconstruction algorithm. First, the signals pd(x, y, t) are transformed to Fourier space in x, y and t coordinates, resulting in pd(kx, ky, ω). Second, the temporal component is mapped to its corresponding spatial coordinate with the dispersion relation, resulting in p0(kx, ky, kz). Finally, the reconstruction is obtained by applying the inverse Fourier transform in all three spatial coordinates. In this way, fast reconstructions of 3D volumes can be obtained using an FFT algorithm. A similar FFT-based approach has been presented by Kunyansky for a finite cube used as the detection surface instead of an infinite planar surface which has to be truncated in practical cases. For a spherical measurement geometry, Wang et al. have recently derived a simple FFT-based algorithm in the frequency domain.

19Another approach is based on closed-form analytical filtered back-projection (BP) formulae in the time domain. This type of formula is commonly used in computed x-ray tomography to invert the classical Radon transform. Similar to Equation (7), the initial pressure is obtained by integrating the detected signals over spherical surfaces. Prior to or following the integration, a filter function is applied. The first exact back-projection formulae were obtained by Finch et al. for a sphere in 3-D and later also in 2-D. Xu et al. presented an approximate BP formula for a spherical detection geometry, which was later generalized to the widely used universal back-projection (UBP) formula:

20The latter scheme can be efficiently applied to spherical, cylindrical or planar (with 4π replaced by 2π) geometries. In the far-field approximation with detectors located far from the object, only the derivative term in Equation (8) contributes significantly.

21Time reversal inverse methods are based instead on the Huygens' principle: In three dimensions, the waves from an initial pressure distribution p0(r), confined inside a finite volume V, will leave the volume after a finite duration T while the pressure p(r, t) inside the volume will vanish for t > T. Although Huygens' Principle does not hold in two dimensions, the pressure inside the volume usually decays fast enough to result in good reconstructions for large T. Time reversal methods aim in solving Equation (4) backwards in time: Starting from zero initial condition at time t = T, the pressure distribution is propagated backwards in time step-by-step, reemitting the pressure profile pd(rd, t) on the detection surface for each time point (Figure 2(c)). Arriving at t = 0, the pressure distribution p(r, t) equals the sought-after initial pressure p0(r). Time reversal methods can be implemented by finite differences methods or by k-space methods, which accelerate inversion by using larger time steps, and furthermore offer greater flexibility as they generally work for arbitrary closed measurement surfaces.

22Whereas series solution or the back-projection formulae perform an approximated inversion of Equation (6) analytically, the so-called model-based (MB) algorithms seek instead for an accurate numerical solution. Rosenthal et al. presented an algorithm that is based on a discretized semi-analytical forward model M̃ of the wave propagation that linearly maps the discretized initial pressure p̃0 to the corresponding discretized signals p̃d. Numerical models can also be obtained by finite element (FE) calculations of Equation (4), which can be implemented in the time or frequency domain. The image is then reconstructed by finding a solution that minimizes the difference between the measured data p̃d and the pressure corresponding to the forward-modeled image p̃0:

23In Equation (9), P is an optional regularization term, such as Tikhonov or total variation regularization, which is intended for stabilization of the numerical inversion. The solution can be calculated via the Moore-Penrose pseudo-inverse M̃†. Although calculation of M̃† may become burdensome in terms of CPU time and memory consumption, it needs to be calculated only once for a given experimental system. Then the reconstruction process reduces to a simple and fast matrix-vector multiplication:

24Alternatively, the solutions can be calculated iteratively, for example with the LSQR algorithm that takes advantage of the sparse nature of the model matrix. Benefiting from the advances in computation technology, MB algorithms have drawn considerable interest in the recent years because of their flexibility to account for a broad variety of experimental imperfections in the model. However, their drawbacks are often associated with longer calculation times or a considerable consumption of memory needed to store the model (or its pseudo-inverse), especially when considering the full three dimensional problem.

25Besides image quality and technical feasibility, computational efficiency is yet another important criterion determining applicability of the different methodologies, especially when dealing with large 3D datasets or in applications involving real-time rendering of images. For instance, different approaches can be characterized by their algorithmic complexity in reconstructing 3-D datasets-O(n4) for series solutions, O(n5) for BP algorithms, O(n4) for time reversal methods and O(n9) for calculating the pseudo-inverse in MB algorithms. The series solutions can result in particularly fast algorithms whenever fast methods for summation of the eigenfunctions are applicable (e.g., FFT), resulting in low complexity in the order of O(n3 · logn). In addition to the algorithmic complexity, additional implementation aspects include the different memory consumption requirements between algorithms that only involve standard linear algebra operations for inversion, such as BP, versus the extensive memory size needed to store matrices for some MB schemes. The graphic processing units (GPUs), which are increasingly used in medical imaging in the recent years due to their vast computational power, are expected to mitigate or significantly reduce complications concerned with computational complexities in optoacoustic imaging.

26As mentioned above, MB algorithms offer better reconstruction accuracy and flexibility in accounting various experimental imperfections; however they demand large computational resources. Consequently, methods to overcome this problem have been presented. Wang et al. proposed a method that circumvents storing the large model matrix by approximating voxels in the image by means of spheres, whose corresponding pressure profiles are known analytically. For the GPU implementation of their algorithm they reported a runtime of 313 min per iteration for a 3-D high resolution data set, whereas the corresponding BP reconstruction took only 46 s. Rosenthal et al. developed an efficient framework for model-based inversions based on wavelet packets. In this case, both the image and the signals are decomposed into a two-level wavelet-packet basis. The model matrix is approximately block-diagonal in this basis and thus an approximate inverse matrix can be calculated (and stored) easily. From the obtained inverse, the image is directly calculated in a more efficient and faster manner and can be then further improved iteratively. The authors reported a runtime of 0.1 s for a 2-D data set by direct reconstruction from the approximate inverse, which constituted a 40-fold increase in the reconstruction speed as compared to the regular iterative MB algorithm.

27Due to mathematical and numerical simplicity, most inversion schemes assume infinitely small point detectors. Conversely, in experimental reality, larger detector sizes are often preferable as they lead to a better detection SNR, which in many cases scales with detector size. Furthermore, large area detectors are also required for focused detection geometries. For instance, if the object is placed in the focal zone of a cylindrically focused detector (Figure 4(a)), the signals acquired by the detector originate only from a small region close to the imaging plane. The reconstruction is thus effectively reduced to a 2-D problem. So data acquisition and image rendering can be performed faster or even in real-time. OA microscopy setups often use spherically focused detectors which are in first approximation only collecting signals from a line, in which case images can be formed without complicated reconstruction algorithms. However, simplified assumptions readily lead to out-of-plane (out-of-focus) artifacts in the images, loss of quantification, highly anisotropic resolution and degradation of the overall image quality.

28Li et al. introduced the virtual detector concept to OA in order to improve the reconstructions from focused detectors outside the focal point (Figure 4(b)). The detected signals are assumed to originate from a virtual point detector located in the focal spot, delayed by the time corresponding to the focal distance. These signals can then be used to form an image, e.g., by a simple delay-and-sum algorithm for OA microscopy or a back-projection algorithm for tomographic reconstructions.

29A very straightforward reconstruction approach can be applied for flat detectors with dimensions much larger than the object (in contrary to point detectors with dimensions much smaller than the object). The so-called integrating line detectors (Figure 4(c)) simultaneously record signals along a line thus effectively the pressure field of the object integrated along a certain axis. A complete tomographic data set is subsequently obtained by rotating the detectors' orientation. The reconstruction process is then performed in two steps: First, the integrated 2-D images are reconstructed for a given detector orientation with an arbitrary 2-D method. In a second step, the 3-D reconstruction is obtained from the integrated 2-D reconstructions for all detector orientations with the inverse Radon transform. Alternatively to line detectors, detection elements, that are integrating in two dimensions or integrating in one and focused in the second dimension, have been suggested as well.

30An interesting approach, which is also based on detection of integrated pressure fields, was presented by Nuster et al.. Instead of detecting a time-resolved integrated field along a certain line, the approach consists of recording the full field at a given time instant td using a CCD camera. In this way, the integrated acoustic pressure field is detected in a direction perpendicular to the imaging plane. td is chosen in such a way that essential parts of the waves have already left the reconstruction domain. A volumetric image can be then reconstructed from a few simple steps: (1) The 2-D pressure field is integrated (Radon transform) along a fixed direction in the imaging plane. (2) The integrated 1-D signals are back-propagated in time from t = td to t = 0 using D'Alembert's formula, which results an integrated 1-D image. (3) An integrated 2-D image is obtained by the inverse Radon transform from the integrated 1-D images for all in-plane directions. (4) A final 3-D image is reconstructed by the inverse Radon transform from the 2-D images for all orientations of the imaging plane.

31Rosenthal et al. presented a generalized method that allows accounting for arbitrary-shape detectors. It is based on including the effects of the exact detector shape in the forward model of a model-based inversion approach by temporal convolution. As a result, image artifacts related to finite detector size can be readily eliminated using this universal approach, while also improving spatial resolution of the images.

32In real tissues, acoustic properties may significantly deviate from the idealized homogeneous properties initially assumed in the governing equations, which may lead to deterioration of the image quality. Indeed, speed of sound c may vary considerably for different types of tissue within a range of 1,400–1,600 m/s and be also spatially dependent (c = c(r)). Information on the speed of sound within the object can be obtained using various approaches, e.g., from a priori knowledge on the object, ultrasound or hybrid measurements, auto-focusing approaches or can be reconstructed from the signals by a combined FE approach. Although simplified series solutions exist in principle for optoacoustic reconstructions with variable speed of sound, iterative and time reversal methods were further shown efficient for this purpose.

33An additional effect taking place is attenuation of the optoacoustic waves as they propagate through the tissue. Often these effects are also dispersive, i.e., acoustic attenuation increases with frequency, which is usually modeled by exponential attenuation with a power-law of ωn, 1 ≤ n ≤ 2. The penetration depth from which acoustic waves can be detected also limits the imaging depth. Since the maximum OA resolution is determined by the frequency content of the signals, the resolution thus scales with depth. Time-gain compensation applied in US imaging cannot be used for OA signals because they are generally broadband. Instead, the attenuation is accounted for by adding a dissipation term in the time domain for time reversal methods or by a complex frequency-dependent dispersion relation in an iterative approach.

34In media with heterogeneous acoustic properties, scattering and reflection of the acoustic waves may occur at boundaries of highly mismatching media. Examples include bones or air cavities such as the lung. Wang et al. studied the effects of reflection at soft and hard boundaries and were able to considerably improve image quality in the presence of reflections at planar boundaries, which might however not always be applicable to realistic imaging scenarios. Anastasio et al. showed that tomographic OA signals, acquired from multiple view angles (projections), contain complementary information thus images can be equally reconstructed from truncated signals, which may improve image quality in heterogeneous acoustic media due to data redundancy principles. Dean-Ben et al. presented a statistically weighted BP algorithm that efficiently accounted for strong acoustic reflections. The algorithm weights each detected signal by a factor that represents the probability of the signal being undisturbed during propagation. In this way, reconstruction artifacts arising from scattering can be suppressed (Figure 5).

35Ideal tomographic imaging scenarios assume that the object is fully surrounded by detectors on a closed surface, creating the so-called full-view scenario. Full view detection is favorable for the reconstruction process but cannot always be realized in experimental conditions, e.g., because the imaged region is simply not accessible for efficient detection from all directions or cannot be fully immersed in water for efficient acoustic coupling to the detectors. Such scenarios are known as limited view reconstructions and the inversion has to be performed with an incomplete set of data.

36Xu et al. investigated the effects of limited view in an OAT setup with detectors on a truncated circle (Figure 6(a)). It has been found that regions enclosed by the detection surface can be recovered stably. Sharp boundaries of the object are thus visible in the reconstructions if they are facing the detection surface but turn invisible if they are perpendicular to the detection surface.

37Buehler et al. showed that limited view reconstructions suffer from stripe artifacts (Figure 6(b)). The artifacts were reduced using a MB approach with a regularization term that suppresses ripples in the direction of the artifacts. Alternatively, total variation regularization terms can be applied in limited view scenarios to stabilize the inversion and reduce the artifacts. Closed form reconstruction formulae can be also optimized for limited view scenarios. Paltauf et al. presented a BP approach that weights the signals with a smooth function before back-projecting. The weighting function needs to consider whether complementary data from an opposite direction are available or these projections are missing. This approach may reduce artifacts resulting from the limited view while also preserving the simple BP implementation. Kunyansky presented an approach for non-closed detection surfaces based on single-layer potentials which, however, necessitates a non-trivial step of finding and calculating such potential. Haltmeier et al. presented an efficient algorithm for a truncated planar detection geometry based on a nonuniform FFT that reduces the limited view artifacts while keeping computational efficiency of the fast Fourier methods.

38Recent reconstruction algorithms for OA have been also driven by the rapidly growing field of compressed sensing (CS), which had been so far successfully applied to biomedical imaging in the field of MRI imaging. Provost et al. were the first to suggest a CS scheme in OA imaging, followed by contributions from other groups. The basic idea of CS is that images p̃0 can be represented by only a few large coefficients in a suitable basis instead of many small coefficients. CS exploits sparsity of the images based on a certain model M̃ in that basis. Instead of minimizing the residual for an image in L2-norm (‖·‖2, Euclidian norm) only, CS also enforces the sparsity of the retrieved image by concurrently minimizing the L1-norm (‖·‖1, Taxicab norm) of the image. Thus, instead of Equation (9), one has to minimize: where Ψp̃0 is the representation of the image in the CS basis and λ is a scalar regularization parameter. Equation (11) can be solved by nonlinear minimization methods, which typically require longer calculation times as compared to just solving Equation (9). In comparison to standard image reconstruction methods, CS approaches have further shown reduction of artifacts in images when reconstructing from only a few tomographic projections (Figure 7). Thus, by reducing the number of acquired projections, these approaches hold great promise for dramatic reduction of image acquisition and reconstruction times as well as instrumentation costs.

39The last section has been entirely devoted to the acoustic part of the optoacoustic phenomena, i.e., reconstruction of the initial pressure distribution and the deposited heat H(r) in the imaged object from the detected signals pd(rd, t). Nevertheless, understanding processes behind the optical counterpart governing propagation of the excitation light is equally important, mainly for two reasons. First, since photons are absorbed as they propagate in tissue, the depth at which a sufficient number of photons is left to excite a detectable signal limits the effective penetration depth of the imaging modality. Second, as stated earlier, the eventual goal of optoacoustic tomography is retrieving distribution of the optical absorption coefficient μa(r) in the imaged object, which can be further related to distribution of the various chromophores and bio-markers in tissues. But the inversion process, described in the previous section, only retrieves maps of the deposited energy H(r), not μa(r). Both distributions are related by the light fluence Φ(r, μa, μs, g) (see Equation (2)), which generally depends on the optical properties of the medium. In order to accurately quantify the values of the optical absorption coefficients and the corresponding concentrations of tissue chromophores, one needs to correct for the light fluence, which may vary considerably versus imaging depth. The reader is referred to a review article on that issue by Cox et al., which has been published recently and covers the subject of image quantification in more detail.

40Several approaches exist to describe transport of light in diffusive tissues, which account for both absorption-which is mainly responsible for the contrast-and scattering processes. For visible and near-infrared light propagating in biological tissues, scattering can be typically considered as the dominant process. As depicted in Figure 8(a), one distinguishes between two main regimes of light transport: ballistic photons that have not been scattered and diffuse photons that have undergone multiple scattering events before being absorbed. At depths relevant for deep tissue OA tomography, typically between several millimeters to several centimeters, diffuse photons are dominant.

41One frequently applied approach to model the transport of diffuse light is the Monte Carlo (MC) method, which simulates a random walk of a great number of photons through the scattering tissue. The ensemble of all simulated paths results in the light fluence. Although MC calculations can be easily parallelized and optimized, they generally require an enormous calculation effort, which cannot be performed in real time.

42An alternative description is given by the time-independent radiative transfer equation (RTE) derived from the Boltzmann equation (see Equation (9) in). It models the light transport as integro-differential equation and holds for both ballistic and diffuse photons.

43Because of the dominating scattering process, light transport in OA is most often modeled by a simpler approximation to the RTE, the so-called light diffusion equation (LDE): where q0(r) is a source term and is the diffusion coefficient (with the reduced scattering coefficient expressed via ). For light incident upon a homogeneous scattering half-space, expression for the light fluence can be readily simplified as Φ(z) = Φ0 · exp (−μeff · z) -well known as Beer's law-with an effective absorption coefficient. For instance, if μeff = 1 cm−1, the light fluence will drop to 1/e of its maximal value after 1 cm (Figure 8 (b)). In this case, light penetration depth in tissue will be limited to a few centimeters.

44If collimated light is incident on the object, the LDE is not accurate for depths close to the surface, where the light has not fully diffused. The δ-Eddington approximation is a higher order approximation to the RTE providing better results in the so-called mesoscopic region, typically up to a few millimeters in most tissues.

45Given a certain light transport model and H(r), Equation (2) needs to be solved for retrieving μa. In general, this means solving a non-linear problem because μa also non-linearly and implicitly depends on H(r) via the light fluence Φ(r, μa, μs, g). Various methods aim at obtaining both μa(r) and μs(r) simultaneously or only μa(r) for a given μs(r).

46The solution is trivial for superficial imaging applications where only ballistic photons are of interest, e.g., as in case of optical resolution OA microscopy. Alternatively, simple correction schemes, accounting for the drop of light fluence with depth, can be applied, e.g., analytical solutions to the LDE considering only the average optical properties of the medium.

47A more complex light transport model can be also solved by an iterative fixed point calculation. Here, the unknown μa(r) is first written as an explicit function:

48A better approximation of the unknown μa(r) is then obtained with each additional iteration and is subsequently used as an improved estimate for the right-hand term in the next iteration. Similarly, the problem can be treated in a linearized way by expanding μa and Φ as a Taylor series (μa = μa,0 + δμa, Φ = Φ0 + δΦ0). The problem can be linearly solved then for small perturbations δμa and δΦ from their known solutions μa,0 and Φ0. The fluence deviation δΦ0 is assumed either as unchanged or as being linearly related to δμa in the Born approximation. The solutions can be then obtained by a standard FE method.

49Instead of separately treating the acoustic and the optical inverse problem, the solution can also be obtained in a hybrid manner. Laufer et al. presented a combined algorithm for modeling distributions of both the light (FE method) and the generated pressure field (k-space method). The method directly recovers the optical parameters quantitatively by iterative non-linear minimization of the difference between the computed and the measured pressure profiles.

50Yet, optoacoustics is a high resolution imaging modality, which makes accurate modeling of the light transport and quantification of OA images challenging. If the various models and equations used in the reconstruction process are not too far away from the experimental reality, the estimate for μa(r) is expected to converge to a stable solution and obtain a good approximate after only a few iterations. In some other cases, solution for the absorption coefficient was shown to be divergent or otherwise the solution was found to be non-unique and ill-posed when trying to recover μa and μs simultaneously.

51The importance of knowing the light fluence in real experiments and the challenges with modeling light propagation in optically heterogeneous tissues have led to some alternative approaches attempting to correct for light attenuation.

52One method to increase the amount of available information regarding optical properties of the imaged tissues is to use multi-illumination patterns. The object is illuminated from different angles with light of the same wavelength. Bal et al. presented an approach that can—under some assumptions—recover the optical properties from a multi-illumination pattern in a non-iterative way. Alternatively to multi-illumination patterns, the inversion process can also be done simultaneously for multi-wavelength illumination with knowledge on the spectral dependence of the optical parameters (see next section).

53Rosenthal et al. presented a blind separation approach which does not consider any particular light transport model. Instead, it is based on a fairly realistic assumption that the light fluence Φ(r) is a slowly varying function of space whereas the absorption coefficient μa(r) has mainly high spatial frequency components. The idea is to sparsely represent a logarithmic representation of the reconstructed optoacoustic image log H(r) in a set of two bases, a wavelet basis Λn(r) and a Fourier basis Ψm(r), i.e., where the coefficients αn are attributed to the absorption coefficient μa(r) and βm to the light fluence Φ(r). In this way, both the absorption and the light fluence can be reconstructed from the OA image blindly (Figure 9), without modeling the light propagation in tissues.

54Cox et al. suggested a method that directly ‘measures’ the local light fluence. Its key component is a chromophore whose absorption has a highly non-linear dependence on the light fluence at a certain threshold intensity. With increase of the illumination intensity, the chromophores are ‘switched’ on (or off) at positions where the local light fluence crosses that threshold. The applicability of this approach however requires the availability of a contrast agent with such highly non-linear absorbance.

55Another approach is to experimentally determine the light fluence using a hybrid imaging approach that also accommodates diffuse optical tomography (DOT). Because of their purely optical nature, DOT images are heavily affected by light diffusion thus cannot provide the same spatial resolution as OA. But they can be instead used to improve quantification of OA images. Maslov et al. implanted an absorber with known optical parameter in tissue at depth to investigate the effects of light fluence. Although being an interesting experimental exercise, this invasive method seems impractical for in-vivo imaging applications.

56The contrast in OA stems from differential light absorption and can be of either intrinsic or extrinsic nature. Figure 10 shows the absorption spectra of the most dominant intrinsic absorbers in tissue in the visible and near-infrared range.

57For many applications, the blood-related absorption is of particular interest, providing valuable physiological or functional information. The main challenge arises from the fact that, at any given wavelength, more than one chromophore contributes to the total absorption coefficient μa(r). For instance, at 800 nm, the oxygenated (HbO2) and the deoxygenated (Hb) forms of hemoglobin contribute equally to absorption and thus are indistinguishable using single wavelength measurements. Therefore, it is the concentrations of the different chromophores ci(r), not the total absorption μa(r), which is mainly of interest from the biological point of view. This relation may be expressed via linear superposition of the different chromophores: where ϵi(λ) are their wavelength-dependent molar extinction coefficients and μback(r, λ) is the residual (background) absorption, which might also include noise. To differentiate between contributions of different chromophores, their distinct spectral dependence on wavelength can be assessed. This multi-wavelength approach is known as multispectral optoacoustic tomography (MSOT) or spectroscopic imaging. The process of recovering ci(r) from multispectral measurements is known as unmixing. It can be combined with calculation of the light fluence or treated as a separate image processing step.

58Due to versatility and wide availability of optical molecular agents, sensitive and accurate spectral processing to recover concentration of extrinsically-administered agents may enable longitudinal molecular imaging studies. This can be done by resolving accumulation of agents with specific spectral signatures, such as targeted and activatable fluorescent molecular agents, nanoparticles or genetic markers.

59If spectra of the different components (chromophores) are known, their concentrations can be directly calculated from the OA absorption images. In some special cases this step is rather trivial, e.g., if the certain chromophore of interest is entirely dominant over other absorbers at a given wavelength or if one is only interested in the total blood concentration measured at the isobestic point, for which Hb and HbO2 are equally absorbing.

60The simplest multispectral OA imaging approach was presented by Kruger et al. and consisted of subtracting the images obtained at two different wavelengths. It was used to recover distribution of an optical agent, having sharp variations in its absorption spectrum, over wavelength-independent absorbing background. A similar approach proved efficient for detection of dynamic contrast based on a pump-probe excitation. Here, a phosphorescent chromophore was pumped to an excited state with a laser pulse at the pump wavelength. The transient absorption of the excited chromophore was then probed with a slightly delayed second laser pulse with a different wavelength. The dynamic contrast can be extracted by subtracting the pump-only absorption from the combined pump-probe absorption. Also, it is well known that blood oxygenation level saO2 is independent of the total concentration of Hb and HbO2 but only depends on their ratio, i.e.,

61Consequently, the oxygenation level can in principle be obtained from the ratio of OA images acquired at two different wavelengths.

62Despite its relative simplicity, the image subtraction method cannot provide quantitative measurements of distribution of several chromophores. The approach can be thus generalized to an arbitrary number of chromophores with unknown concentrations and known spectra. In this way, for a given set of wavelengths, the concentrations are linearly fitted to Equation (15) on a per-voxel basis using a linear regression method. This can be efficiently realized with a matrix relation: where ci and μλ are respectively the concentrations of the components and the absorption coefficient in a certain reconstructed image voxel. is the Moore-Penrose pseudo-inverse of the spectra matrix S derived from the components' molar extinction coefficients. In general, the number of measured wavelength should be greater than the number of chromophores.

63In many cases, the absorption spectra of all the absorbing components present in the imaged object may not be exactly known. As a result, the unmixing method needs to simultaneously retrieve both the chromophore concentrations and their spectra. In contrary to the previously presented methods, the so-called blind unmixing methods do not operate on a per-voxel basis but exploit the statistical properties of a set of images without requiring prior knowledge of the spectra.

64Glatz et al. investigated the performance of blind spectral unmixing on OA images using two multivariate methods. Principal component analysis (PCA) is based on the assumption that the different chromophore distributions are statistically uncorrected and yields an orthonormal transformation into a new base, in which the largest data variance is projected onto the first principal component, the largest remaining variance onto the second one, and so on. In this way, the spectrally correlated measurement data are unmixed by being transformed to the uncorrected components. Similarly, independent component analysis (ICA) is based on a more general assumption that the source components are statistically independent. Consequently, the spectrally mixed measurement data are transformed to statistically independent source components. Both methods are able to recover the spectra of the different components via the transformation matrix and their concentrations via the transformed images. Figure 11 summarizes the results obtained with the linear fitting, PCA and ICA algorithms.

65In this section, we shortly highlight several representative examples of the emerging optoacoustic imaging applications while a more comprehensive review on biomedical applications of optoacoustics can be gained from other recently published review articles.

66Naturally, the best intrinsic tissue contrast arises from highly absorbing hemoglobin, thus blood is clearly visible in the images. Functional OA microscopy is thus geared towards investigation of vascular structures, providing an excellent intrinsic contrast and a high spatial resolution in the order of some tens of μm and penetration of several millimeters into scattering tissues. In order to obtain high resolution 3-D images, a single spherically-focused high frequency detection element is mechanically raster-scanned along a plane. Significantly higher spatial resolution is further achieved by using a tightly focused optical illumination spot, resulting in the so-called optical resolution microscopy. Despite providing much better spatial resolution in the sub-micron range, similarly to all the other optical microscopy techniques, the latter method suffers from limited penetration depth of only several hundreds of microns due to intense photon scattering in biological tissues. An example of a typical acoustic resolution in vivo functional optoacoustic microscopy system is shown in Figure 3(a). A tunable laser generated pulses of 6-ns duration which passed through a fiber and formed a ring-shaped illumination pattern to suppress OA signals from the surface. A coaxially aligned focused 50-MHz central frequency detector acquired time-resolved one-dimensional signals. The transducer was scanned in 50 μm steps in the x-y plane to form a 3-D image without an inversion algorithm needed. With this setup, subcutaneous vasculature of a Sprague-Dawley rat was demonstrated. Figure 12(a) shows the structural image obtained at the isobestic wavelength of 584 nm in grayscale and the functional saO2 image obtained by a linear fitting to four different wavelengths in color. The saO2 levels obtained were 0.97 and 0.77 for arterial and venous blood, respectively. Changes in saO2 for hypoxic and hyperoxic conditions of the animal were successfully tracked, offering interesting prospects for non-invasive functional brain imaging.

67Preclinical whole-body imaging of small animals with multispectral optoacoustic tomography (MSOT) systems is yet another key application of OA imaging. These systems can provide tomographic images with a resolution in the order of hundred μm for depths between several millimeters to several centimeters. Possible applications include monitoring of tumor hypoxia, drug response or molecular targets in biological model organisms. The ability of MSOT to visualize deep-seated fluorescent proteins with high resolution has also been demonstrated. The setup utilized a selective-plane light illumination and a confocal ultrasound detection pattern (Figure 3(b)). An OPO laser provided 8-ns laser pulses of tunable wavelength and the signals were captured by rotating a 3.5-MHz central frequency transducer around the object in 3° steps for multiple slices. The images for multiple wavelengths were reconstructed with the UBP algorithm and unmixed with a linear fitting on a per-pixel basis. In this way, a transgenic three-month old zebra fish expressing mCherry fluorescent proteins in the vertebral column was visualized in vivo. Figure 12(b) shows the structural OA image of the hindbrain at 585 nm and its corresponding histological slice (top row). The multispectral OA imaging has accurately attained the mCherry expression, which has also well corresponded with the findings from epi-fluorescent histology (bottom row).

68Buehler et al. presented the capabilities of MSOT for targeted molecular imaging in real time. In the in vivo mouse experiments, a ring-shaped illumination was provided by a tunable OPO (680–950 nm, 10-ns, 10 Hz) and delivered through a fiber bundle. A 64 element transducer array (5-MHz central frequency) covering 172° solid angle detected the signals in parallel without multiplexing. The CD-1 nude mice containing 4T1 tumor allografts were wrapped in a transparent membrane to prevent direct contact of the mouse with water. A near-infrared apoptosis targeting probe (PSS-794, 50-100 nanomoles) was injected in three mice at different time points. The images were reconstructed using a model-based algorithm and unmixed with an ICA spectral unmixing algorithm. Figure 13(a) shows the MSOT image containing both high resolution anatomical and functional information for a tumor-bearing mouse with PSS-794 probe. From the MSOT images (Figure 13(a)), it was clearly determined that the PSS-794 probe mainly accumulated in the blood vessels surrounding the tumors and did not infiltrate into the tumor mass. On the other hand, the use of three-dimensional optical tomography (FMT) has not attained the sufficient spatial resolution that would determine the precise location of the probe (Figure 13(c)), mainly due to its poor spatial resolution resulting from photon scattering.

69At present, the great potential of optoacoustic imaging that was showcased in preclinical research has encouraged translation of this technology into clinical practice. The clinical promise of OA is greatly supported by its non-ionizing radiation, rich contrast, real-time operation, and relatively low cost. It is likely that OA will increasingly enter the clinical imaging segments in areas such as skin, breast, vasculature, and ocular imaging, visualization of tumor and inflammation-related pathology etc.. Indeed, early detection of breast cancer will capture a significant spot in the future developments of OA imaging applications. One example of such clinical system for OA breast angiography without the need for a contrast agent was recently presented. The system employed 128 detection elements (5-MHz central frequency) which were placed on a 100 mm radius hemisphere in a spiral pattern and acquired the signals in parallel (Figure 3(f)). The hemisphere was rotated for higher projection sampling and one complete scan took between 6 to 24 s. The volunteer's breast was immersed into a water-filled bowl from the top and immobilized by an optically and acoustically transparent cup. 3-D images were reconstructed from a filtered BP algorithm with prior frequency filtering and transducers' impulse response deconvolution. The system was used to visualize the vasculature of a 57 year-old volunteer without a contrast agent and vessels from superficial regions up to a considerable depth of several centimeters (Figure 14(a)).

70Naturally, clinical application of OA imaging is limited to regions that can be efficiently illuminated with the excitation light. When the imaged region of interest is not reachable from the outside, certain areas of the human body can also be imaged with minimally invasive catheter- or endoscopy-based techniques. To this end, such systems have been realized for other imaging modalities, such as intra-vascular ultrasound (IVUS) or fluorescent imaging (NIRF). In the recent years, similar approaches have been also attempted with OA for intravascular, esophageal or colonoscopic imaging. One example is shown in Figure 3(d), where a system for simultaneous functional OA and US endoscopy in vivo is presented. The endoscope was able to acquire co-registered pulse-echo US and OA images. A-scans were obtained from an ultrasonic pulse-echo system while two OA pulses, having different wavelengths for spectroscopic imaging, were delivered through a fiber and detected by the focused transducer. 3-D images could be obtained by constant rotation of a scanning mirror and a mechanical pullback of the catheter. Using this system, rabbit esophagus was imaged in vivo, as shown in Figure 14(b).

71Indeed, optoacoustic imaging offers multiple advantages over other imaging modalities. Most favorable is the combination of US-diffraction limited spatial resolution and excellent optical absorption contrast without any ionizing radiation involved. Clearly, much like all other imaging modalities, optoacoustics is limited in certain aspects. As presented in the current review, due to the (non-ideal) highly heterogeneous and wavelength-dependent nature of biological tissues, images reconstructed from experimental data usually contain artifacts, such as negative values in the images which otherwise have no physical meaning. In addition, high noise levels, limited bandwidth of the detection system, inaccurate assumptions on the speed of sound, acoustic attenuation, as well as limited view in detection, can all alter the signals in a way that the reconstructed image accuracy is highly compromised. While some simple algorithms, such as Hilbert transform for envelope detection or thresholding of negative values, can be used to keep image values in the positive region, these non-linear operations in fact only conceal the underlying problem while not assisting with achieving quantitative results. Instead, reliable detection and modeling of signal propagation in tissues is required.

72Imaging at depth brings further challenges and limitations. First, imaging depth is compromised by attenuation of the light fluence in optically opaque tissues. Some of it can be compensated by increasing the amount of the deposited laser energies. However, for in vivo applications in the near-infrared, illumination on the skin surface is limited by the maximum permissible exposure (MPE,) to about 20 mJ/cm2. As a result, imaging depth is usually restricted to regions where the light fluence is sufficiently high to generate detectable pressure variations, typically up to a few centimeters in most soft tissues. In addition, even though spatial resolution of optoacoustic imaging is not affected by light scattering, resolution may still deteriorate with increasing imaging depth due to its dependence on the frequency content of the detected ultrasonic responses. This is because of the dispersive nature of ultrasound waves with high frequency components being rapidly attenuated as they propagate through tissues.

73In multispectral imaging applications, one also faces the problem of so-called ‘spectral coloring’. Due to the non-local and non-linear dependence of light transport on the object's optical properties, spectra of various tissue chromophores and agents, extracted by means of optoacoustics, might be corrupted. For improving quantitative determination of chromophore concentrations using their known spectra, an accurate correction for the light distribution needs to be employed, which, as discussed before, constitutes a challenging problem and an open area of research. Another important aspect influenced by depth is the limits of the imaging system with respect to detection of minimal concentrations of certain intrinsic tissue chromophores and extrinsically-administered contrast agents. To this end, several publications have reported sensitivity in the femtomole to picomole range, depending on the molecular weight of the probe employed, while a recent study has also systematically addressed the question of how the detection limits are affected by imaging depth. In reality, sensitivity limits are affected by multiple additional factors, such as the total volume, spectrum and absorption coefficient of the imaged chromophore, the noise equivalent pressure (NEP) of the detectors, level and spectral dependence of background tissue absorption.

74Instrumentation-related limitations are an additional factor in translating optoacoustic imaging methods into routine application. One of the important advantages, which helps reducing motion-related image artifacts but may also enable powerful applications in in-vivo tracking of dynamic processes, is the possibility for reliable measurements and image rendering in real-time. This requires introducing high repetition pulsed laser technology and fast wavelength tuning capabilities, especially for multi-spectral imaging applications involving real-time visualization of distributions of spectrally distinct contrast agents. On the other hand, the detection technology should become increasingly sensitive but also parallelized to enable real-time image acquisition. Often a trade-off between acquisition speed and quality on the one hand and costs of instrumentation on the other hand has to be found. In this respect, optoacoustic imaging has greatly benefited from the advances in both laser and detection technology over the last years and significant performance enhancements are expected as the technological frontiers advance. Providing information from hybrid US imaging will be of added value for future clinical systems. Besides fast data acquisition, fast image rendering is also of great importance to readily provide findings for efficient real-time guidance and optimization of diagnostic or therapeutic procedures. Recently, extensive research is devoted to the subject of real-time visualization, benefiting both from algorithmic developments and advances in parallel computing technology.

75Quantitative rendering of chromophore concentration is perhaps among the most important but also among the most challenging tasks of the optoacoustic methods. Hemoglobin in its oxygenated and deoxygenated form is omnipresent in tissue and is the most dominant intrinsic absorber in both the visible and most of the near-infrared range of the optical spectrum. Thus, quantitative determination of blood parameters is of great interest as it is related to a variety of physiological parameters. Prominent examples are monitoring angiogenesis or hypoxic states in tumors, functional imaging in the brain or responses to environmental changes. In most cases, the oxygen saturation saO2 (Equation (16)) is of main interest, whose quantitative determination is non-trivial, especially in the presence of strong wavelength dependence of the excitation light fluence in tissues. First attempts to reliably determine saO2 were conducted in vitro whereas accuracy of ±2.5% saO2 and a minimum detectable change of ±1% saO2 were reported. Furthermore, oxygenation saturation of single blood vessels, embedded deep in tissue mimicking phantom, was determined using numerical forward light transport model, reporting an accuracy of ±7% saO2. Most in vivo results were so far obtained from multispectral optical resolution optoacoustic microscopy where light transport issues can be effectively omitted. Zhang et al. reported 4% saO2 systematic error for their ex vivo study with bovine blood, while accuracy for in vivo determination of arterial and venous blood under hypoxic, normal and hyperoxic conditions was not reported. Applications of efficient and practical saO2-determination schemes for experimental data in 3-D in an accurate, robust and computationally feasible manner presents therefore one of the key challenges for the current algorithmic developments. Finally, when imaging other chromophores or extrinsically-administered contrast agents, the blood background needs to be accounted for to enable quantified extraction of biomarker concentrations.

76Owing to its hybrid nature, i.e., optical excitation and ultrasonic detection, optoacoustics benefits from both rich and versatile optical contrast and high (diffraction-limited) spatial resolution, associated with relatively low scattering of ultrasonic waves. Optoacoustic biosensing and imaging provides an excellent platform for multi-scale investigations using the same contrast, from microscopic observations at the single capillary and cell level to whole body imaging of small animals and deep tissue imaging of humans. Much like other optical imaging modalities, optoacoustics is safe for both small animals and clinical use as it utilizes non-ionizing radiation at the visible and near-infrared wavelengths.

77In the last decade, OA imaging is considered to be the fastest growing biomedical imaging modality with multiple already implemented and envisioned applications in biomedical research and clinical practice, from diagnostic applications in cancer research and brain imaging to drug development and treatment monitoring. Clearly, progress of the image reconstruction and quantification methods has been central to those developments and has significantly contributed to creation of new exciting applications of the OA imaging technologies. As discussed in this review, multiple frontiers are still open in the algorithmic areas, where many challenges related to removal of image artifacts, image quantification, reconstruction strategies in the presence of acoustically mismatched areas, real-time operation, and multi-spectral data processing need to be further addressed.