Maximum entropy based non-negative optoacoustic tomographic image reconstruction
Jaya Prakash Subhamoy Mandal Daniel Razansky Vasilis Ntziachristos
physics.med-ph cs.CV eess.IV physics.optics





Maximum entropy based non-negative optoacoustic 
tomographic image reconstruction
    Jaya Prakash^†, Subhamoy Mandal^†, Member, IEEE, Daniel Razansky,  Member, IEEE, 
and Vasilis Ntziachristos, Senior Member, IEEE 

This article has been accepted for publication in IEEE
Transactions on Biomedical Engineering. DOI: 10.1109/TBME.2019.2892842
^†J.P. and S.M. contributed equally to this work.
J.P. acknowledges support from the Alexander von Humboldt Postdoctoral Fellowship Program. S.M. acknowledges support from DAAD PhD Scholarship Award (A/11/75907) and IEEE Richard E. Merwin Scholarship. D.R. acknowledges funding support from the European Research Council (ERC-2015-CoG-682379), US National Institutes of Health (R21-EY026382-01), Human Frontier Science Program (RGY0070/2016) and Deutsche Forschungsgemeinschaft (RA1848/5-1). V.N. acknowledges funding support from European Research Council (694968, ERC-PREMSOT).
J.P., S.M., D.R., and V.N. are with the Institute of Biological and Medical Imaging, Helmholtz Zentrum Munich, Ingolstaedter Landstr. 1, D-85764 Neuherberg, Germany, and also with the Chair for Biological Imaging, Technical University Munich, Ismaningerstr 22, D-81675, Munich, Germany. (e-mail: v.ntziachristos@tum.de).
    December 30, 2023
=========================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================




Objective: Optoacoustic (photoacoustic) tomography is aimed at reconstructing maps of the initial pressure rise induced by the absorption of light pulses in tissue. In practice, due to inaccurate assumptions in the forward model, noise and other experimental factors, the images are often afflicted by artifacts, occasionally manifested as negative values. The aim of the work is to develop an inversion method which reduces the occurrence of negative values and improves the quantitative performance of optoacoustic imaging. 
Methods: We present a novel method for optoacoustic tomography based on an entropy maximization algorithm, which uses logarithmic regularization for attaining non-negative reconstructions. The reconstruction image quality is further improved using structural prior based fluence correction. 
Results: We report the performance achieved by the entropy maximization scheme on numerical simulation, experimental phantoms and in-vivo samples. 
Conclusion: The proposed algorithm demonstrates superior reconstruction performance by delivering non-negative pixel values with no visible distortion of anatomical structures. Significance: Our method can enable quantitative optoacoustic imaging, and has the potential to improve pre-clinical and translational imaging applications. 



Optical parameters, photoacoustic tomography, inverse problems, image reconstruction, regularization theory.




§ INTRODUCTION


Optoacoustic (OA) imaging detects broadband ultrasound (pressure) waves generated within tissue in response to external illumination with light of transient energy, due to light absorption by tissue elements and thermo-elastic expansion. Using forward models that describe sound propagation in tissue, ultrasound measurements from multiple positions surrounding the object imaged are mathematically reconstructed to resolve the spatial distribution of the initial pressure rise. The reconstructed pressure rise is proportional to the product H=μ_a ϕ, whereby μ_a is the optical absorption coefficient and ϕ is the light fluence <cit.>. The value H has only positive values in biological tissues since both absorption and light fluence are positive. However, the appearance of negative values is common in OA images due to different factors, such as the use of inaccurate forward models, inversion schemes, numerical errors, limited view detection geometry, transducer impulse response, unknown or unpredictable experimental effects or noise in the imaging system. The presence of negative values in the reconstruction does not have physical relevance. Importantly, when spectral techniques are employed, such as Multispectral Optoacoustic Tomography (MSOT) <cit.>, the presence of negative values make spectral quantification problematic. 

It is therefore important to treat the appearance of negative values in the OA tomography problem. Model based reconstruction has been suggested as an alternative to back-projection algorithms to improve the accuracy of OA imaging, further incorporating transducer and laser characteristics into the inversion procedure <cit.>. In principle, accurate inversion can reduce the image artifacts, but errors persist due to different experimental challenges including limited-angle signal collection, limited bandwidth detection, noise and other uncertainties, leading to incomplete data problems and results in the presence of erroneous negative values <cit.>. Consequently, methods to directly treat the problem of negative values have been considered <cit.>. Ding et. al. <cit.> compared the utility of different minimization procedures using non-negative constraints, including steepest descent, conjugate gradient, and quasi-newton based inversion. Typical non-negative constraint schemes truncate the negative values within each step of the gradient iteration, forcing a result containing only positive or zero values. This practice however may bias the solution and generate inaccuracies in the reconstruction.  

An alternative approach to address the problem of negative values is to use image content for image correction. Image features such as the total energy (smoothness), contrast, total variation of an image can be generally employed as prior information to direct the inversion towards pre-determined outcomes, usually based on the assumptions about the nature of the image. For example, ℓ_2- or ℓ_1-norm minimization of the total variation of an image minimizes the edges of the reconstructed image. Using this notion, negative artifacts can then be eliminated by applying an explicit non-negativity constraint along with ℓ_2-norm minimization <cit.>. Another image metric that has been considered for eliminating negative values is the entropy of an image <cit.>. Entropy is the measure of randomness in an image. Randomness of the image implies that information from each subpixel is assumed to be independent of each other and can statistically take any value irrespective of its neighboring subpixel. This becomes very useful in limited data situations; wherein the principle of maximum entropy tries to eliminate all uncertainties within each subpixel (among the different possible solutions) by imposing independent statistical structure on each pixel. Maximization of entropy (i.e. maximizing the term -xlog(x); whereby x is the vectorized image) is equal to minimizing the term xlog(x) and is a method considered in  Positron Emission Tomography (PET) and multi-modal imaging <cit.> or astronomical imaging <cit.>. 
 
In this work, we examine the use of entropy as a prior in OA image inversion, in the context of nonlinear conjugate gradient minimization <cit.>. We hypothesize that the use of an entropy-based prior, which implements an implicit non-negativity constraint, can improve the accuracy of OA inversions over externally imposed non-negativity constraints. To prove this hypothesis, we first theoretically compare a conventional ℓ_2-norm minimization problem using a smoothness constraint to an entropy maximization problem. We show that images reconstructed by entropy maximization cannot take negative values. The reconstructed OA images were further improved by correcting for the fluence, the fluence was estimated using finite volume method after segmenting the imaging domain (phantom or mouse). Thereafter, we compare the performance of inversion (after fluence correction) using entropy maximization and conventional inversion with externally applied non-negativity constraint using numerical simulation, experimental phantoms and small animal imaging. We discuss the performance differences observed and the advantages and limitations of using entropy maximization.     




§ MATERIALS AND METHODS




 §.§ Theoretical  background

The propagation of the acoustic pressure wave generated due to the short-pulsed light absorption is governed by the following inhomogenous wave equation <cit.>,

    ∂^2 p(r,t)/∂ t^2 - c^2 ρ∇.(1/ρ∇ p(r,t))= Γ∂ H(r,t)/∂ t,

where the instantaneous light power absorption density in W/m^3 is indicated by H and Γ represents the medium-dependent dimensionless Grueneisen parameter. In Eq. <ref>, the tissue density is represented by ρ while c indicates the speed of sound (SoS). For our experiments, a uniform SoS of 1520 m/sec was heuristically estimated using image autofocusing method <cit.>. The initial pressure rise at position r and time t is given as p(r,t). The solution for the wave equation can then be obtained using a Green's function by assuming H(r,t)=H_r(r) δ(t), which results in <cit.>,

    p(r,t) = Γ/4 π c∂/∂ t∫_R=ctH_r(r')/R dr',

where R=ct represents the radius of the integration circle over a line element given as dr'. The above solution is subsequently discretized into the following matrix equation <cit.>,

    b=Ax,

where b is the boundary pressure measurements, A is the interpolated model matrix and x is the unknown image to be reconstructed, representing the initial pressure rise distribution. The above formulation represents the forward model, i.e. given the initial pressure rise one can estimate the pressure at the boundary locations detected by the transducers. Thus, the acoustic inverse problem involves reconstructing the initial pressure rise given the boundary pressure data. In the ℓ_2-norm formulation, the inverse problem is solved by minimizing a function given as, 

    Ω_ℓ_2 = arg min_x (||Ax-b||_2^2 + λ|| Lx||_2^2),

where λ is the regularization parameter. The term || Ax-b||_2^2 is called the residual term. The term || Lx||_2^2 is a ℓ_2-norm of the second order total-variation of the image x and L indicates the Laplacian operator. The value of the regularization parameter affects the resolution characteristics of the reconstructed image; higher the value of regularization the smoother the reconstructed image. 



 §.§ Entropy Maximization and Non-negative constraint

An alternative method to the minimization problem of Eq. <ref> (optoacoustic reconstruction), is maximization of the entropy of the image. To elaborate on this point, a statistical approach is considered, wherein we assume that the image to be reconstructed follows a Gaussian distribution with estimated mean and standard deviation values. The dimension of the image to be reconstructed is N × N, i.e. a vector of size NN (=N^2). Next, we assume that each pixel j in this image will be formed by a group of subpixels indicated by m_j (>=1) and M=∑_j=1^NN m_j. With these assumptions, let us consider the following experiment: wherein K particles are distributed over all subpixels and let K_i be the number of particles that fall in pixel i. Then the number of combinations to place K particles in NN pixels such that K_j particles are present in pixel j is given as,

    C(K̃) = K!/Π_j=1^NN K_j!,

Further we have m_j^(K_j) ways to put K_j particles into m_j subpixels. Hence, the total number of combinations to create the particle distribution V(K̃) is given as,

    V(K̃) = C(K̃). Π_j=1^NN m_j^(K_j),

The total number of particles in the distribution is given as M^K. Now making the assumption that each particle is equally likely i.e. uniform distribution. We get the probability of distribution of K̃ as,

    p(K̃) = V(K̃)/M^K,

Now using Stirling approximation i.e. K! ≈ K^K e^(-K), we can write,

    log(p(K̃)) = -K ∑_j=1^NN z_j log(z_j/m̂_j),

where z_j=K_j/K and m̂_̂ĵ = m_j/M. The average value inside a pixel x_j will now be proportional to z_j i.e. x_j=S.z_j and m̃_̃j̃ = S.m̂_̂ĵ such that,

    ∑_j=1^NN x_j = ∑_j=1^NNm̂_̂ĵ = S,

with x_j ≥ 0, m̂_̂ĵ>0.


Now let the prior distribution of the image vector be considered as p_A (x), which is given as,

    log(p_A (x)) = -K/S∑_j=1^NN x_j log(x_j/m_j),

which follows the relative entropy definition and is always non-negative (not defined for negative values). Our next assumption is that the error vector or the noise is normally distributed with zero mean and standard deviation σ given as,

    p(r_i) = c.e^-r_i^2/2σ^2,

which can be rewritten as,

    p(y|x) = c.e^|| A x - b||_2^2/2σ^2,

Rewriting the overall expression using Bayes rule we get,

    log(p(x|y)) = -K/S∑_j=1^NN x_j log(x_j/m_j) - 1/2σ^2||A x - b||_2^2,

Neglecting the terms independent of x. We can pose this as an entropy maximization problem which is non-linear convex maximization problem, and this can be solved by minimizing the function,

    Ω_maxent = arg min_x (||Ax-b||_2^2 + λ∑_i=1^NN x_i log(x_i/m_i)),

where -xlog(x/m) indicates the relative entropy function of image x, typically m is assumed to be an arbitrary constant <cit.>. In this work m is assumed to be 1. Detailed mathematical analysis on the use of Eq. <ref> for applying an implicit non-negativity constraint, stability, and convergence of entropy maximization is given in <cit.>. Herein we study how positive values are retained with entropy maximization scheme. 

In ℓ_2-norm minimization (Eq. <ref>), the gradient update equation at iteration i is given as,

    x_i = x_i-1 - ( A^T ( A x_i-1 - b)) - λ L^T  L x_i-1,

The above update equation is obtained by taking the derivative of the objective function in Eq. <ref>. Note that in the above equation all the quantities will always be in real space i.e. (A, x_i-1, x_i, b  ∈  IR), and can take any values due to the absence of any natural non-negativity barrier. Therefore, the ℓ_2-norm based minimization can generate negative values (which can be in IR) during the image reconstruction procedure. In case of entropy maximization (Eq. <ref>), the gradient updated equation at iteration i is given as,

    x_i = x_i-1 - ( A^T ( A x_i-1 - b)) - λ (1 + log(x_i-1/m_j-1)),

The derivation pertaining to applying implicit positivity constraint using entropy maximization is discussed in the Appendix-I.

Choice of regularization plays a key role in reconstructed image quality by defining over-smoothed or under-smoothed approximations in case of ℓ_2-norm based reconstruction. In terms of distance measure, ℓ_2-norm constraint can be considered as Euclidean distance between the prior and the expected image, i.e. || Lx||_2 = < Lx, Lx> ⟺ < Lx, Lx_pr><cit.>, therefore higher regularization will weigh the ℓ_2-norm constraint more and thus resulting in a smoother solution. Similarly entropy maximization can be related to Kullback-Leiber distance, as cross entropy between prior and the expected image, i.e. ∑ x log(x) = ∑ x log(x/x_pr), therefore higher regularization will push the subpixels (i.e. x_pr) in pixel i of image vector x to uniform distribution <cit.>. Thus, low regularization in the entropy maximization scheme will result in minimizing the residual (i.e. noisy reconstruction), whereas choosing higher regularization will result in the initial pressure rise being close to a smooth distribution having intrinsically positive values. The operating range of the regularization parameter in the entropy maximization framework can be found using the L-curve type method, cross-validation based scheme <cit.>. 




 §.§ Choice of regularization parameter - L-curve method

Typically, the regularization parameter (λ) is chosen automatically using the L-curve method <cit.>. 
The L-curve method is a popular method for automatically choosing the regularization parameter for a linear inverse problem and this scheme was earlier used in diffuse optical tomography and OA tomography. In the L-curve method, a graph is plotted between the residual (||Ax-b||_2^2) and the reconstruction (||x_λ||_2^2) as function of regularization parameter (λ). This essentially means that the reconstructed solution (x_λ) is a function of regularization (λ). In an ideal case this curve will be of L-shape. The corner point of this L-shape represents the least distance from the origin, indicating an ideal balance between residual and expected solution.
For the case of entropy maximization the solution norm will be replaced by entropy term i.e. (Σx_λ log(x_λ)).
In this work, we use L-curve type approach to automatically estimate the regularization parameter in both the L2-norm and entropy maximization schemes. 



 §.§ ℓ_2-norm with smoothness and non-negativity constraint

Minimizing the function in Eq. <ref> was performed using a conjugate gradient method (equivalent to iterative least squares QR (LSQR) method), which has a closed form solution as <cit.>,

    x ≈ x_ℓ_2-lsqr = V_k (B_k^T B_k + λS_k^T S_k)^-1β_0 B_k^T e_1,

where B_k, S_k, V_k, β_0, and e_1 can be obtained in the Lanczos diagonalization procedure with [   A; λ L ] and [ b; 0 ]. Here k indicates the number of iterations during the joint bidiagonalization procedure. 

In the ℓ_2-norm formulation with non-negativity constraint, the following minimization is solved, 

    Ω_ℓ_2-NN = arg min_x (||Ax-b||_2^2 + λ|| Lx||_2^2)         s.t.   x>0,

The above minimization is solved using the LSQR solver and then the obtained solution containing negative values are thresholded to 0, as negative values do not have any physical relevance (as optical absorption coefficient in biological tissue is not negative). Eq. <ref> is used to obtain the solution and then the negative values in the solution are thresholded. The regularization parameter was chosen using L-curve method (explained in Sec. II-C)<cit.>.



 §.§ Implementation Steps for Entropy Maximization

Eq. <ref> is minimized using a non-linear conjugate gradient type method and the step-length for the conjugate gradient method is computed using a line search <cit.>. Minimization of the objective function in Eq. <ref> with conjugate gradient requires computing the derivative and then move in independent perpendicular gradient direction. The derivative used in the conjugate gradient scheme for the objective function in Eq. <ref> is computed as,

    ∇Ω_maxent = 2 A^T (A x -b) + λ (1 + log(x_i-1/m_i-1)),

The minimization is presented in more details in the Algorithm-1 section. The regularization parameter was chosen using an L-curve method (as a tradeoff between negative of entropy and residual).






 §.§ Fluence Correction

The image reconstructed in Eq. <ref> (LSQR) and with Algorithm-1 (Entropy Maximization) represents the absorbed energy distribution H_r (r) in tissue, which depends on the fluence distribution and the optical absorption coefficient  μ_a (r) i.e <cit.>,

    x = p_0 (r) = H_r (r) = μ_a (r) Φ (r),

where p_0 (r) is the initial pressure rise distribution and Φ(r) indicates the local light fluence density in mJ/cm^2. To extract the absorption coefficient map, it is therefore critical to estimate the fluence in the medium imaged. Different schemes have been developed for estimating the fluence distribution and quantitatively recover optical absorption coefficient maps, including model-based inversion schemes integrated with fluence compensation <cit.>, wavelet frameworks <cit.>, finite-element implementation of the delta-Eddington approximation to the radiative transfer equation <cit.>, diffusion equation based regularized Newton method <cit.>, or approximations with base spectra <cit.>. Herein we assumed for demonstration purposes a light propagation model based on the diffusion equation, further assuming that scattering dominates over absorption <cit.>, which is a valid approximation for  most biological tissues and NIR measurements, i.e.,

    -∇.[D(r).∇Φ(r)] + μ_a(r) Φ(r) = S_0 (r),

where D(r)=1/(3(μ_a+μ_s^')) is the diffusion coefficient and μ_s^'(r) indicates the reduced scattering coefficient at position r. S_0 (r) indicates the light source at the boundary of the imaging domain. Eq. <ref> is used for fluence estimation, and the diffusion equation is solved using the finite volume method (FVM). Optical properties were based on the known phantom specifications or estimates of absorption and scattering coefficients of tissue from the literature <cit.>. Then, we obtained absorption coefficient maps by normalizing the images with the corresponding calculated fluence distribution <cit.>. Since OA measurements of phantoms were performed in a water bath, we also employed the Beer-Lambert Law (OD = -log(I/I_0) = -μ_a d) to model photon propagation in water. The relative distances in phantom and water were assigned after segmentation of the OA images. The entire workflow of segmentation and fluence correction is integrated with the proposed non-negative entropy maximization algorithm to render improved image quality.




 §.§ Imaging instrumentation and protocol(s):

Experimental data was acquired using the multispectral optoacoustic tomography (MSOT) scanner <cit.> (MSOT256-TF, iThera Medical GmbH, Munich, Germany). The boundary pressure readouts (time-series) were collected at 2,030 discrete time points at 40 Mega samples per second using a 256-element cylindrically focused transducer, resulting in the number of measurements (M) being 2030x256=519,680. The utilized piezocomposite transducer had a central frequency of 5 MHz with a radius of curvature of about 40 mm and an angular coverage of 270^∘. Uniform illumination was achieved with a ring type of light delivery using laser fiber bundles. Numerical simulations were performed with the same configuration as MSOT256-TF system with a realistic breast phantom having spatially varying absorption coefficient (in cm^-1) as shown in Fig. 1(a). Next, we segmented the boundary of the breast region in Fig. 1(a) and estimated the fluence distribution (shown in Fig. 1(b)) by solving the hybrid model (Sec. II.F) with the absorption coefficient and reduced scattering coefficient set to 0.2 cm^-1 and 12 cm^-1 respectively. The initial pressure rise (in kPa) was then estimated by multiplying the fluence distribution (Fig. 1(b)) with the spatially varying optical absorption (Fig. 1(a)), the initial pressure rise distribution (after scaling with acoustic parameters) is shown in Fig. 1(c). Note that we assumed point detector and did not model transducer characteristics in the simulations. The numerical breast phantom was created by using contrast enhanced magnetic resonance imaging <cit.>. Eq. 3 was used to model the acoustic propagation (on a 512×512 grid) and the pressure signals were collected at specific detector locations. The model matrix in Eq. 3 was built using interpolated model matrix method as explained in  Ref. <cit.>. To avoid inverse crime, the simulated data was generated on an imaging grid of size 512x512, while the reconstruction was performed on a grid of size 256x256. The simulated data was added with additive white Gaussian noise, to result in a SNR of 32 dB in the simulated data.

To verify the quantitative reconstruction capabilities of the proposed entropy maximization scheme, a star shaped (irregular) phantom was created. The phantom constituted of a tissue mimicking (7% by volume of Intralipid and pre-computed volume of diluted India ink added) agar core having the optical density of 0.25. Two tubular absorbers made up of India-ink with the absorption coefficient values of 2.5 OD (calibrations done with Ocean Optics USB 4000) were inserted in the phantom. The absorbers were placed at two different depths within the phantom (one at the center and the other at the edge of the imaging domain) to test the sensitivity of the proposed scheme in reconstructing the absorbers at different imaging distances from the sensing arrays. Under normal operating conditions, the fluence at the center of the imaging domain is significantly lower as compared to the boundary of the object imaged, owing to the optical attenuation of the incident irradiation. Hence, performing fluence correction becomes indispensable to assign appropriated intensity to the absorber at the center of the imaging domain. 

The proposed methods were further validated on in-vivo mouse abdomen and brain datasets drawn from a standardized in-vivo murine whole body imaging database (10 mice/30 anatomical datasets) previously developed in Ref. <cit.>. The selected images were obtained at a laser wavelength of 760 nm and 800 nm, and the water (coupling medium) temperature was maintained at 34^∘ C for all experiments. Non-negativity based entropy maximization scheme was further validated using spectral measurements. Spectral measurements were acquired from a tumor bearing nude BALB-C mice with the laser wavelengths running from 680 nm to 900 nm at steps of 20 nm. All animal experiments were conducted under supervision of trained technician in accordance with institutional guidelines, and with approval from the Government of Upper Bavaria.




 §.§ Figure of merit

To develop an objective approach to evaluate imaging performance of different reconstruction methods, we used line plots on the reconstructed image (from phantom and tissue measurements). We also performed quantification using sharpness metric, defined as,

    SM = ∑dI^2/dx^2 + dI^2/dy^2/n,

The sharpness metric indicates the edges in the reconstructed image (I): the higher the value of SM, the sharper the reconstructed image. This figure of metric was used for evaluating the proposed method, as the non-negative constraint tend to introduce zeros in the reconstructed image. The number of non-negative values is also reported for comparing the different reconstruction methods.  Note that the number of negative pixels were calculated from the phantom or mice region (excluding the water region).

Further root mean square error (RMSE) and peak signal to noise ratio (PSNR) was used to evaluate the performance of different reconstruction methods with numerical simulation. RMSE is given as,

    RMSE = √(∑_o (x_o^recon - x_o^true)^2/NN),

is computed for comparing the performance of different algorithm. Here x_o^true is the o^th pixel of ground truth and x_o^recon is the o^th pixel of reconstructed image. PSNR is defined as,

    PSNR = 20 × log(max(x^true)/RMSE),


The calculated sharpness metrices (for phantom and in vivo small animal images), and the RMSE/ PSNR values (for simulations) are given in section <ref>. 



§ RESULTS



Fig. 1(c) shows the initial pressure distribution with the realistic numerical breast phantom used to evaluate the performance of different reconstruction methods. The reconstructed initial pressure rise distribution using the ℓ_2-norm based reconstruction is shown in Fig. 1(d). The solution pertaining to ℓ_2-norm based reconstruction (along with non-negative constraint) is indicated in Fig. 1(e). The reconstructed optoacoustic image using the entropy maximization approach is represented in Fig. 1(f). The reconstructions containing negative values are indicated with a red colormap, hence the negative pixels in Fig 1(d) are shown in red color. From the numerical simulations, it is apparent that the ℓ_2-norm based reconstruction produces negative values by just adding noise to the data and incorporating fluence effects, however these negative values do not appear after thresholding and using entropy maximization scheme as indicated by red arrows in Figs 1(e) and 1(f). Furthermore ℓ_2-norm with thresholding results in a nosier reconstruction with limited structures compared to entropy maximization scheme as shown with red arrows in Figs 1(e) and 1(f). The PSNR values for ℓ_2-norm, ℓ_2-norm with thresholding and entropy maximization reconstruction are 29.9736 dB, 30.2616 dB and 30.3529 dB respectively. The RMSE values for ℓ_2-norm, ℓ_2-norm with thresholding and entropy maximization reconstructions are 0.0453, 0.0451, and 0.0450 respectively. The number of reconstructed negative pixels with ℓ_2-norm reconstruction with numerical breast phantom is 4370. Note that the simulation studies did not model many experimental parameters like impulse response of the transducer, physical dimension of the transducer, pitch of the detector, artifacts arising due to reflections, and these parameters are known to influence the OA measurements in experimental scenarios. Further, we proceeded to study the performance of the proposed entropy maximization scheme with phantom and in-vivo datasets.


Fig. 2 shows reconstructions of the star phantom, which reveal the efficacy of the proposed method vis-a-vis traditional ℓ_2-norm based reconstruction in generating positive values for both the initial pressure rise and absorption coefficient distribution. The reconstructed initial pressure rise and absorption coefficient distribution using the ℓ_2-norm based reconstruction is shown in Figs 2(a) and 2(d) respectively. The reconstructed initial pressure rise and absorption coefficient distribution using the ℓ_2-norm based reconstruction (with non-negative constraint) is indicated in Figs 2(b) and 2(e) respectively. The reconstructed initial pressure rise and absorption coefficient distribution using the entropy maximization based approach is represented in Figs 2(c) and 2(f) respectively. The reconstructions containing negative values are indicated with a red colormap, hence the negative pixels in Figs 2(a) and 2(d) are shown in red color.


The proposed entropy maximization method (Fig. 2(f)) can provide accurate image representation with the ability to reconstruct the absorber (having OD of 2.5) at the center and the edge of the imaging domain along with reconstructing a star shaped background (having OD of 0.25). The negative values obtained using LSQR inversion is shown as red color in Fig. 2(a) and Fig. 2(d). The non-negative based ℓ_2-norm reconstruction is able to generate reconstruction results with positive values, but is not able to correctly reconstruct the internal volume of the star (tissue mimicking agar with 0.25 OD) phantom which is accurately reconstructed using entropy maximization. Fig. 2(g) shows the photograph of the phantom used from front-view (FV) and top-view (TV). Fig. 2(h) indicates the line plot along the vertical red dashed line shown in Fig. 2(b). Fig. 2(i) indicates the line plot along the horizontal blue dashed line shown in Fig. 2(b). The sharpness metric and the number of non-negative values are shown in Table-<ref>. The quantitative metric indicate that the proposed method can provide accurate image representation. Fig. 2(f) and the line plots in Figs 2(h) and 2(i) demonstrate that the maximum entropy based scheme can deliver better contrast while maintaining the background intensity than the standard ℓ_2-norm based reconstructions.The fluence correction was performed by using segmented (boundary) priors obtained automatically using deformable active contour models <cit.>. The results were corroborated with additional phantom (Agar block with 5% intralipid) scans which included India ink insertions of 3 different ODs in tissue relevant concentrations - 0.15, 0.30 and 0.45 OD at 800nm measured using a spectrometer (VIS-NIR; Ocean Optics). The results demonstrate that the signal intensities change proportionately with the changing OD of the insertions, and the values are in agreement with other commonly used inversion algorithm (i.e Tikhonov). The reported signal intensities were obtained by taking the mean of the different ROI's indicated in Table-I of the supplementary. Additionally, the proposed reconstruction scheme recovered higher (absolute) signal intensities while reducing negative values in reconstructed image (see supplementary Table I).

Empirically selecting the regularization biases the reconstruction results.  Therefore an L-curve method was used to automatically choose the regularization parameter for Tikhonov method <cit.> and entropy maximization based scheme. Previous works have used L-curve approach for automatically choosing  the regularization parameter in entropy maximization framework for estimating distance distributions of magnetic spin-pairs <cit.>. Fig. 3 indicates the L-curve criterion used to choose the regularization parameter (details regarding L-curve approach is given in Sec. II-C) as applied to star phantom OA data presented in Fig. 2. Similar approach was used for automatically selecting the regularization parameter with numerical simulations and in-vivo data. Other methods like cross-validation can also be used for automatically choosing the regularization parameter in Tikhonov and entropy based framework <cit.>. Further, we studied the effect of regularization parameter choice on reconstruction image quality. Fig. S1 in supplementary shows maximum entropy reconstruction at different regularization parameter values. It can be seen that at high regularization values, the solution leads to uniform distribution, however maximum entropy scheme seems to have a large operating range from 1 to 10,000.



The maximum entropy based scheme depends on the initial guess used in the non-linear conjugate gradient scheme. The maximum entropy constraint involves a non-linear logarithmic term, and the logarithm of a negative value is not defined, therefore having a large positive value at the initial guess will always generates positive reconstruction distributions and thus plays an important role in intrinsically obtaining non-negative reconstruction. The same is elaborated in the Appendix-I. The reconstruction results corresponding to a backprojection-type initial guess (A^T b containing negative values; A^T indicates transpose of system matrix) is indicated in Fig. 4(a), the image shows the real part of the solution. The reconstruction results corresponding to the initial guess (||b||_2/|| A||_1× ones(NN,1)) is indicated in Fig. 4(b). Fig. 4(a) clearly indicates that the negative values in the entropy maximization reconstructions arises because of initial guess used in the non-linear conjugate gradient scheme i.e. (||b||_2/|| A||_1× ones(NN,1)) gives non-negative results while A^T b results in negative values.  Hence, in all the reconstructions the initial guess was chosen to be (||b||_2/|| A||_1× ones(NN,1)) and the regularization parameter was chosen using the L-curve method. Note that reconstructions in Fig. 4 involve performing additional fluence correction. The colormap in the case of mouse images are normalized to maximum and minimum values and the negative values are indicated in red color. 


Non-negative reconstruction generated with entropy maximization approach was further improved using fluence correction method. Fig. 5(a) shows the performance of segmentation approach in delineating the interface/boundary between the mice body (at the abdominal region) and water. The segmented boundary is used as a source term (after attenuation compensation using Beer-Lambert law in water) for modeling light propagation by solving the diffusion equation. Indeed, this boundary can be a good approximation for source term, as fiber bundle in the MSOT machine are arranged to provide uniform illumination on the sample. The fluence profile obtained after solving diffusion equation is shown in Fig. 5(b), the fluence was estimated with optical properties obtained from the literature <cit.>. Fig. 5(c) represents the initial pressure rise distribution reconstructed with entropy maximization approach. Fig. 5(d) shows the absorption coefficient distribution after normalizing the initial pressure distribution (Fig. 5(c)) with the estimated fluence profile (Fig. 5(b)). It can be clearly seen that signals from deeper regions on the mice gets highlighted more, similar approach was used for other regions of the mice. 




The reconstruction results (corresponding to absorption coefficient distribution) pertaining to the mouse head and mouse abdominal regions using the standard and proposed method are shown in Fig. 6. The reconstruction results corresponding to ℓ_2-norm based scheme (solved using LSQR method) for the mouse head and abdominal region is indicated in Figs 6(a) and 6(e) respectively, and the corresponding results for ℓ_2-norm based non-negative scheme (solved using LSQR method with thresholding) are given by Figs 6(b) and 6(f) respectively. The reconstruction results using the entropy maximization approach (Algorithm-1 with the integrated hybrid fluence correction) for the same anatomical regions is shown in Fig. 6(c) and Fig. 6(g) respectively. The experimental phantom and in-vivo reconstructions were performed on a 200x200 pixel imaging domain which corresponds to a physical field of view of 20mm x 20mm. The optical properties used for fluence estimation was assumed to be homogenous inside the tissue and taken from literature  <cit.>. Figs 6(d) and 6(h) indicate the Fourier domain representation of the reconstructed images (i.e. Fig. 6(f) and 6(g)) using L2-norm with thresholding and entropy maximization schemes respectively. We could clearly see that entropy maximization scheme (Fig. 6(h)) has more low frequency content when compared to L2-norm with thersholding (Fig. 6(d)). Fig. 6(i) indicates the line plot along the red dashed line shown in Fig. 6(b) and Fig. 6(j) shows the line plot along the red dashed line indicated in Fig. 6(f). The sharpness metric and the number of non-negative values for these reconstructions are indicated in Table-<ref>. These metrics show that the proposed method can provide accurate image reconstruction with lesser negative values and increased sharpness. Negative values should not arise during standard OA data acquisition, hence the lesser the number of negative pixels more accurate is the reconstructions. However in some scenarios the presence of negative values might indicate accurate reconstruction like temperature dependent studies <cit.>. However, we are working with standard OA acquisition, and thus more positive values indicate accurate reconstruction. Again, the colormap is normalized to maximum and minimum values, while indicating the negative values in red color.




Finally, we performed a study to check if entropy maximization scheme was able to accurately recover the spectral information. Fig. <ref>(a) shows the reconstruction results pertaining to a tumor bearing mice using L2-norm based scheme with thresholding at 680 nm wavelength. Fig. <ref>(b) shows the recovered mean spectral information using entropy maximization and L2-norm based reconstruction for the red square region shown in  Fig. <ref>(a).  Fig. <ref>(b) indicates that at wavelengths below 700 nm, we have appearance of negative values using L2-norm based reconstruction. Moreover, in some parts of the image, like the one shown using orange arrow in Fig. <ref>(a), the entire recovered spectra turned out to be negative using L2-norm based reconstruction (however maximum entropy scheme was able to recover positive spectral profile). Fig. <ref>(c) shows reconstructed mean spectra information using entropy maximization and L2-norm based reconstruction from the green square region indicated in Fig. <ref>(a). As can be seen from Figs <ref>(b) and <ref>(c), the spectral recovery of maximum entropy scheme is similar to that of L2-norm based reconstruction, however the appearance of negative values in L2-norm based reconstruction will hinder unmixing results in terms of absolute quantification.




§ DISCUSSION AND CONCLUSION


The reconstruction results for the numerical simulations, phantom and in-vivo mouse scans indicate that the proposed entropy maximization scheme renders strictly positive image values that are also close to the a-priori known absorption values in the phantom. Employing a segmented image prior can effectively reduce the aberrations in image contrast by suitably mapping the light propagation pathway in two optically diverse domains (background and tissue), and enhance the performance of (optical) fluence correction methods<cit.>, as demonstrated in Figs 2(f) and 6(g).  Moreover, when a global SoS is attribute to the entire imaging domain, small SoS variation causes aberration at the edge of the surfaces of the imaged object <cit.>, the same two compartment model can be used to remove SoS mismatch. The figure of merits (Table-<ref>), magnitude of Fourier spectrum from the reconstructed images, and the line plots indicate entropy maximization approach provides superior results in comparison with non-negativity constrained reconstructions. Importantly the proposed approach offers an opportunity for exploring a family of differential type non-negative regularization methods (like entropy scheme).


The entropy maximization scheme performed better with experimental data (Figs 2 and 6) compared to numerical simulation (Fig. 1). This is because experimental OA measurements are heavily influenced by experimental factors like laser pulse width, transducer impulse response, pitch and size of the transducer, making the reconstruction problem with experimental OA measurements more challenging. From Figs 2 and 6, it can be observed that the presence of negative pixels is higher in water region and in the center of imaging domain, where the absorption/the fluence is low resulting in lower SNR in time-series OA measurements. Similarly, introduction of noise and fluence effects in simulation studies (Fig. 1) results in large number of negative values in regions where the initial pressure rise is close to 0 and also generating spurious negative values inside the numerical breast phantom.

In recent studies, lot of emphasis has been placed on using ℓ_1-norm based minimizations for performing OA tomographic image reconstruction in different frameworks <cit.>. We have performed ℓ_1-norm based reconstruction as explained in <cit.> and the results pertaining to non-negativity constraint in the ℓ_1-norm minimization is shown in Fig. S2. Fig. S2 also shows the performance comparison of ℓ_1-norm minimization with entropy maximization and Tikhonov reconstruction with printed phantom data. We observe that applying a ℓ_1-norm constraint does not afflict the appearance of negative values and the reconstruction performance is similar to ℓ_2-norm based scheme in terms of reducing negative values. This also demonstrates the superiority of using entropy maximization to generate physically relevant OA reconstructions devoid of negative values. We have not taken up further comparisons with ℓ_1-norm based approach, as our goal was to demonstrate the utility of entropy maximization approach to overcome appearance of pixels with negative values.

Entropy maximization scheme was evaluated with biological datasets acquired from 270^∘ detection angle wherein the acquired dataset consists of highly independent (incoherent) data. While recent developments involve building systems with handheld probes (90^∘ three-dimensional acquisition, or 145^∘ two-dimensional acquisition) with different data-collection geometry. Performing accurate reconstructions with these clinical handheld systems tend to be difficult due to acquisition of limited independent data. Evaluating the performance of the entropy scheme with the limited independent data scenarios can enable utility of OA imaging in different clinical scenarios<cit.>. 

The proposed method preserves the structural integrity (numerical breast phantom and star phantom) and the anatomical structures (mouse data), and was successful in correcting the effects of variations in optical fluence. As part of future work, we aim to integrate the entropy maximization with more accurate light propagation modeling (such as Monte Carlo based schemes) to obtain better representation of the absorption coefficient with the reconstruction process accelerated by means of graphics processing units <cit.>. In this work, we demonstrated a non-negative image reconstruction method with improved image quality using fluence correction step at single acquisition wavelength. Translating the same to multi-wavelength scenario for estimation of quantitative tissue parameters is a fairly complex problem, since the optical properties used for fluence estimation varies nonlinearly with wavelength and is not known beforehand. Combining these problems will lead to generation of infinite possible ways to obtain accurate spatio-spectral representation, and such spectral analysis methods are beyond the scope of the current study. 

In this work, we have shown that entropy maximization is able to accurately recover the spectral information compared to L2-norm based reconstruction (see Fig. 7). However, the ability to resolve intrinsic chromophores like oxyhemoglobin, deoxyhemoglobin, fat, and water by acquiring data at multiple wavelengths is a key benefit of multispectral OA imaging. The unmixing of chromophores is achieved by a solving system of linear equations (direct or non-negatively constrained), or by non-linear unmixing using an integrated fluence correction. All of these approaches use thresholding of negative values, making them suboptimal and error prone. On the other hand, entropy maximization can  purge out the inaccuracies occurring from truncated pixel information, potentially improving the performance of unmixing and image analysis algorithms. Therefore, the future work will involve comparing the different combination of reconstruction (acoustic inverse problem) and unmixing with different solvers like LSQR, non-negative LSQR and entropy maximization to bring out value among these schemes.



§ CONCLUSION


The proposed maximum entropy based OA image reconstruction scheme demonstrates superior reconstruction performance with no visible distortion of anatomical structures associated with delivering of non-negative pixel values. Entropy maximization reconstruction thus tends to be physically relevant and more accurate in resolving the structures (as demonstrated with numerical simulation, experimental phantoms and in-vivo case) in an imaged sample. The developed methodology has the potential to emerge as a suitable data processing tool for OA imaging, and specifically benefiting pre-clinical biomedical <cit.> and translational imaging <cit.>.
 




§ IMPLICIT NON-NEGATIVITY USING ENTROPY MAXIMIZATION


The objective function in the entropy maximization scheme is given as,

    Ω = | |  Ax - b||_2^2 + λx^T log( x/m)

The gradient of the above equation can be written as,

    ∂Ω/∂ x =  A^T(  Ax - b) + λ( 1 + log( x/m)) = 0

Now, we can consider the above minimization problem as minimizing two models in the subspace, one is based on residual i.e. Res = | | Ax - b||_2^2 and the other being relative entropy i.e. Ent = ∑ xlog( x/m). Here the regularization parameter defines the proportion of residual and entropy term in this minimization problem. As in any optimization, the solution is always found using the search directions (these search directions are defined by the gradients). The update equation at i^th gradient iteration will turn out to be,

    x_i = x_i - 1 - α( ∂Ω/∂ x)_x_i - 1

where α is the step length estimated using line search method and is always non-negative. As x_i - 1→ 0, ∇ Ent →  - ∞, the gradient update will be pushed to a very low value using entropy constraint. Also note that as, x_i - 1→ 0, ∇ Ent will reach - ∞ faster, and the ∇ Res →  - A^Tb; importantly ∇ Res cannot reach ∞ as fast as ∇ Ent → - ∞ to nullify the effect of entropy term, therefore the overall gradient will be negative i.e. ( ∂Ω/∂ x)_x_i - 1→  - ve. In any gradient descent method, we traverse in the direction perpendicular to the gradient, therefore the solution will be pushed away from zero to have high positive value, i.e. as x_i - 1→ 0, x_i→  + ve.  Hence, using the entropy constraint will enable the solution to move away from zero and leading to positive real numbers. Since, a natural barrier is created by including the entropy constraint into the optimization framework, this barrier will not allow the solution to take negative values and consequently positive OA reconstructions are generated. In order to converge to positive OA reconstructions, we need to start with a large positive initial guess i.e. when x_0→ IR^ + then ∇ Res →A^TAx_0 and ∇ Ent → IR^ +. Further, using a step-length control i.e. α  = min( - ∇Ω( x_i - 1))^T∇Ω( x_i - 1 - α∇Ω( x_i - 1)) will ensure positive OA reconstructions, because the choice of α (estimated using secant method) would ensure positive solution in next iteration x_i = x_i - 1 - α∇Ω( x_i - 1).

00

b1 V. Ntziachristos, “Going deeper than microscopy: the optical imaging frontier in biology," Nat. Meth., vol. 7, no. 8, pp. 603–14, Aug. 2010.

b2 P. Beard, “Biomedical photoacoustic imaging," Inter. Focus, vol. 1, no. 4, pp. 602–31, Aug. 2011.

b2a S. Mandal, X. L. Dean-Ben, N. C. Burton and D. Razansky, "Extending Biological Imaging to the Fifth Dimension: Evolution of volumetric small animal multispectral optoacoustic tomography.," in IEEE Pulse, vol. 6, no. 3, pp. 47-53, May-June 2015

b3 A. Taruttis and V. Ntziachristos, “Advances in real-time multispectral optoacoustic imaging and its applications," Nat. Photon., vol. 9, no. Apr., pp. 219–227, 2015.

b4 A. Taruttis, A. Rosenthal, M. Kacprowicz, N. C. Burton, and V. Ntziachristos, “Multiscale multispectral optoacoustic tomography by a stationary wavelet transform prior to unmixing," IEEE Trans. Med. Imag., vol. 33, no. 5, pp. 1194–202, May 2014.

b5 A. Rosenthal, V. Ntziachristos, and D. Razansky, “Model-based optoacoustic inversion with arbitrary-shape detectors," Med. Phys., vol. 38, no. 7, pp. 4285–4295, 2011.

b6 K. Sivasubramanian, V. Periyasamy, K. K. Wen, and M. Pramanik, “Optimizing light delivery through fiber bundle in photoacoustic imaging with clinical ultrasound system: Monte Carlo simulation and experimental validation," J. Biomed. Opt., vol. 22, no. 4, p. 041008, 2016.

b6a Y. Zhen and H. Jiang, “Quantitative photoacoustic tomography: Recovery of optical absorption coefficient maps of heterogeneous media," App. Phys. Lett., vol. 88, no. 23, p. 231101, 2006.

b6b P. Shao, B. Cox, and R. J. Zemp, “Estimating optical absorption, scattering, and Grueneisen distributions with multiple-illumination photoacoustic tomography," App. Opt., vol. 50, no. 19, p. 3145-3154, 2011.

b7 A. Buehler, A. Rosenthal, T. Jetzfellner, A. Dima, D. Razansky, and V. Ntziachristos, “Model-based optoacoustic inversions with incomplete projection data," Med. Phys., vol. 38, no. 3, pp. 1694–1704, 2011.

b8 J. Prakash, A. S. Raju, C. B. Shaw, M. Pramanik, and P. K. Yalavarthy, “Basis pursuit deconvolution for improving model-based reconstructed images in photoacoustic tomography," Biomed. Opt. Exp., vol. 5, no. 5, pp. 1363–77, May 2014.

b9 D. Queiros, X. L. Dean-Ben, A. Buehler, D. Razansky, A. Rosenthal, and V. Ntziachristos, “Modeling the shape of cylindrically focused transducers in three-dimensional optoacoustic tomography," J. Biomed. Opt., vol. 18, no. 3, pp. 076014, 2013.

b10 L. Ding, X. Luis Dean-Ben, C. Lutzweiler, D. Razansky, and V. Ntziachristos, “Efficient non-negative constrained model-based inversion in optoacoustic tomography," Phys. Med. Biol., vol. 60, pp. 6733–6750, 2015.

b11 K. Wang, R. Su, A. Oraevsky, and M. Anastasio, “Investigation of iterative image reconstruction in three-dimensional optoacoustic tomography," Phys. Med. Biol., vol. 57, no. 17, pp. 5399–5423, 2012.

b12 Y. Han, L. Ding, X. L. D. Ben, D. Razansky, J. Prakash, and V. Ntziachristos, “Three-dimensional optoacoustic reconstruction using fast sparse representation," Opt. Lett., vol. 42, no. 5, p. 979, 2017.

b13 S. Somayajula, C. Panagiotou, A. Rangarajan, Q. Li, S. R. Arridge, and R. M. Leahy, “PET image reconstruction using information theoretic anatomical priors," IEEE Trans. Med. Imaging, vol. 30, no. 3, pp. 537–549, 2011.

b14 J. Tang and A. Rahmim, “Anatomy assisted PET image reconstruction incorporating multi-resolution joint entropy," Phys. Med. Biol., vol. 60, no. 1, pp. 31–48, 2014.

b15 J. Skilling and R. K. Bryan, “Maximum Entropy Image Reconstruction - General Algorithm," Mon. Not. R. Astron. Soc., vol. 211, p. 111, 1984.

b16 R. Fletcher, Practical Methods of Optimization, vol. 1: Unconst. 1987.

b17 M. Xu and L. V. Wang, “Photoacoustic imaging in biomedicine," Rev. Sci. Instrum., vol. 77, no. 4, p. 41101, 2006.

b18 S. Mandal, E. Nasonova, X. L. Dean-Ben, and D. Razansky, “Optimal self-calibration of tomographic reconstruction parameters in whole-body small animal optoacoustic imaging," PACS, vol. 2, pp. 128–136, Sep. 2014.

b19 A. Rosenthal, D. Razansky, and V. Ntziachristos, “Fast semi-analytical model-based acoustic inversion for quantitative optoacoustic tomography," IEEE Trans. Med. Imag., vol. 29, no. 6, pp. 1275–85, Jun. 2010.

b20 P. C. Hansen and D. P. O’Leary, “The Use of the L-Curve in the Regularization of Discrete Ill-Posed Problems," SIAM J. Sci. Comput., vol. 14, no. 6, pp. 1487–1503, 1993.

b21 C. B. Shaw, J. Prakash, M. Pramanik, and P. K. Yalavarthy, “Least squares QR-based decomposition provides an efficient way of computing optimal regularization parameter in photoacoustic tomography," J. Biomed. Opt., vol. 18, no. 8, p. 80501, 2013.

b22 G. Landl and R. S. Anderssen, “Non-negative differentially constrained entropy-like regularization," Inv. Prob., vol. 12, no. 1, pp. 35–53, 1996.

b23 B. A. Ardekani, M. Braun, B. F. Hutton, I. Kannof, and H. Iida, “Minimum cross-entropy reconstruction of PET images using prior anatomical information," Phys. Med. Biol., vol. 41, no. 11, pp. 2497–2517, 1996.

b24 B. Borden, “Maximum entropy regularization in inverse synthetic aperture radar imagery," IEEE Trans. Sig. Process., vol. 40, no. 4, pp. 969–973, 1992.

b25 E. Levitan and G. T. Herman, “A Maximum a Posteriori Probability Expectation Maximization Algorithm for Image Reconstruction in Emission Tomography," IEEE Trans. Med. Imag., vol. 6, no. 3, pp. 185–92, 1987.

b26 U. Amato and W. Hughes, “Maximum entropy regularization of Fredholm integral equations of the first kind," Inv. Probl., vol. 7, no. 6, pp. 793–808, 1991.

b27 Y. W. Chiang, P. P. Borbat, and J. H. Freed, “Maximum entropy: A complement to Tikhonov regularization for determination of pair distance distributions by pulsed ESR," J. Magn. Reson., vol. 177, no. 2, pp. 184–196, 2005.

b28 C. C. Paige and M. a. Saunders, “LSQR: An Algorithm for Sparse Linear Equations and Sparse Least Squares," ACM Trans. Math. Softw., vol. 8, no. 1, pp. 43–71, Mar. 1982.

b29 P. C. Hansen, “REGULARIZATION TOOLS: A Matlab package for analysis and solution of discrete ill-posed problems," Numer. Algo., vol. 6, no. 1, pp. 1–35, 1994.

b30a S. Bu et al., “Model-Based Reconstruction Integrated With Fluence Compensation for Photoacoustic Tomography," IEEE Trans. Biomed. Engg, vol. 59, no. 5, p. 1354-1363, 2012.

b30b A. Rosenthal, D. Razansky, and V. Ntziachristos, “Quantitative optoacoustic signal extraction using sparse signal representation," IEEE Trans. Med. Imaging, vol. 28, no. 12, pp. 1997–2006, Dec. 2009.

b30c F. M. Brochu, J. Brunker, J. Joseph, M. R. Tomaszewski, S. Morscher and S. E. Bohndiek, “Towards Quantitative Evaluation of Tissue Absorption Coefficients Using Light Fluence Correction in Optoacoustic Tomography," IEEE Trans. Med. Imag., vol. 36, no. 1, p. 322-331, Jan. 2017.

b30d Z. Yuan, Q. Wang, and H. Jiang, “Reconstruction of optical absorption coefficient maps of heterogeneous media by photoacoustic tomography coupled with diffusion equation based regularized Newton method," Opt. Exp.,  vol. 15, no. 26, p. 18076-18081, 2007.

b31 S. Tzoumas, A. Nunes, I. Olefir, S. Stangl, P. Symvoulidis, S. Glasl, C. Bayer, G. Multhoff, and V. Ntziachristos, “Eigenspectra optoacoustic tomography achieves quantitative blood oxygenation imaging deep in tissues," Nat. Commun., vol. 7, no. May, p. 12121, 2016.

b32 S. R. Arridge, “Optical tomography in medical imaging," Inv. Probl., vol. 15, no. 2, pp. R41–R93, 1999.

b33 S. L. Jacques, “Optical Properties of Biological Tissues: A Review," Phys. Med. Biol., vol. 58, no. 11, pp. R37-61, 2013.

b34 T. Jetzfellner, D. Razansky, A. Rosenthal, R. Schulz, K.-H. Englmeier, and V. Ntziachristos, “Performance of iterative optoacoustic tomography with experimental data," Appl. Phys. Lett., vol. 95, no. 1, p. 013703, 2009.

b35D. Razansky, A. Buehler, and V. Ntziachristos, “Volumetric real-time multispectral optoacoustic tomography of biomarkers," Nat. Protoc., vol. 6, no. 8, pp. 1121–9, Aug. 2011.

b36Y. Lou, W. Zhou, T. P. Matthews, C. M. Appleton, and M. A. Anastasio, “Generation of anatomically realistic numerical phantoms for photoacoustic and ultrasonic breast imaging," J. Biomed. Opt., vol. 22, no. 4, p. 041015, 2017.

b37 S. Mandal, X. L. D. Ben, and D. Razansky, “Visual Quality Enhancement in Optoacoustic Tomography using Active Contour Segmentation Priors," IEEE Trans. Med. Imag., vol. PP, no. 99, p. 1, 2016.

b38E. Petrova, A. Liopo, A. A. Oraevsky, and S. A. Ermilov, “Temperature-dependent optoacoustic response and transient through zero Grüneisen parameter in optically contrasted media," PACS, vol. 7, pp. 36–46, 2017.

b39J. Jose, R. G. H. Willemink, W. Steenbergen, C. H. Slump, T. G. van Leeuwen, and S. Manohar, “Speed-of-sound compensated photoacoustic tomography for accurate imaging," Med. Phys., vol. 39, no. 12, pp. 7262–71, Dec. 2012.

b40 H. He, J. Prakash, A. Buehler, and V. Ntziachristos, “Optoacoustic Tomography Using Accelerated Sparse Recovery and Coherence Factor Weighting," Tomography, vol. 2, no. 2, pp. 138–145, Jun. 2016.

b41Y. Han, S. Tzoumas, A. Nunes, V. Ntziachristos, and A. Rosenthal, “Sparsity-based acoustic inversion in cross-sectional multiscale optoacoustic imaging," Med. Phys., vol. 42, no. 9, pp. 5444–5452, 2015.

b42_1 P. K. Upputuri, and M. Pramanik, “Recent advances toward preclinical and clinical translation of photoacoustic tomography: a review," J. Biomed. Opt., vol. 22, no. 4, p. 041006, 2016.

b42 S. L. Jacques, “Coupling 3D Monte Carlo light transport in optically heterogeneous tissues to photoacoustic signal generation," PACS, vol. 2, no. 4, pp. 137–142, 2014.

b43 V. Ermolayev, X. L. Dean-Ben, S. Mandal, V. Ntziachristos, and D. Razansky, “Simultaneous visualization of tumour oxygenation, neovascularization and contrast agent perfusion by real-time three-dimensional optoacoustic tomography,” Eur. Radiol., vol. 26, no. 6, 2016.

b44 M. Heijblom et al., “The state of the art in breast imaging using the Twente Photoacoustic Mammoscope: results from 31 measurements on malignancies,” Eur. Radiol., pp. 1–14, 2016.



