## Abstract

The statistical independence between the distributions of different chromophores in tissue has previously been used for linear unmixing with independent component analysis (ICA). In this study, we propose exploiting this statistical property in a nonlinear model-based inversion method. The aim is to reduce the sensitivity of the inversion scheme to errors in the modelling of the fluence, and hence provide more accurate quantification of the concentration of independent chromophores. A gradient-based optimisation algorithm is used to minimise the error functional, which includes a term representing the mutual information between the chromophores in addition to the standard least-squares data error. Both numerical simulations and an experimental phantom study are conducted to demonstrate that, in the presence of experimental errors in the fluence model, the proposed inversion method results in more accurate estimation of the concentrations of independent chromophores compared to the standard model-based inversion.

Published by The Optical Society under the terms of the Creative Commons Attribution 4.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

## 1. Introduction

Photoacoustic tomography is a non-invasive biomedical imaging technique [1] relying on the absorption of optical energy and the generation of ultrasound waves. It has a relatively low cost of implementation and has the advantage of combining the relatively large penetration depth and high resolution of ultrasound imaging with optical absorption based contrast which provides high specificity. In quantitative photoacoustic tomography, the unique optical absorption spectra of the chromophores are exploited in spectroscopic inversions of the multiwavelength photoacoustic images in order to estimate the concentration of each chromophore. The key endogenous chromophores of interest for quantitative photoacoustic tomography are oxy- and deoxyhaemoglobin, because the ratio of their concentrations is related to the blood oxygenation, which is an important physiological parameter. Spectroscopic techniques are also valuable tools for contrast-enhanced photoacoustic molecular imaging applications, where the detection and quantification of the local accumulation of genetically encoded probes and extrinsically administered contrast agents [2] can provide information on biological processes, drug delivery, disease development and treatment response.

Quantifying the chromophore concentrations is a challenging task because the absorbed optical energy density is not linearly related to the chromophore concentrations, but a product of the absorption coefficient and light fluence, which varies both spatially and spectrally and depends on the unknown chromophore concentrations. The model-based inversion scheme [3–8] is a quantification method that accounts for the effect of the fluence by including a numerical model of the fluence distribution in an iterative scheme to solve the inverse problem. Thus, it has the potential to provide accurate estimation of the absolute chromophore concentrations in complex tissue structures. However, in experimental settings, both the approximate nature of the model, and the difficulty in determining all the “known” model parameters, such as the absorption and scattering spectra and the intensity profile of the excitation beam, to high accuracy, lead to model-mismatch. This poses challenges on the practical implementation of model-based inversion schemes for *in vivo* imaging. Previous studies [5,8] have relied on reducing the number of unknown variables by segmenting the images into regions with piece-wise constant optical properties to increase the robustness of the inversion scheme. In this study, we exploit instead the statistical independence between certain chromophores to improve the accuracy and robustness of model-based inversion.

Statistical independence has been utilised to spectrally decompose the chromophores via independent component analysis (ICA) [9]. ICA is a fast and simple unmixing algorithm based on the assumption that the photoacoustic images are linear mixtures of the independent source components representing the chromophores. Glatz *et al* [10] showed that ICA can result in more accurate unmixing than a linear spectroscopic inversion. ICA has subsequently been used to unmix genetic reporters [11] and contrast agents targeted to apoptotic cells [12] or cells expressing selectin [13]. However, ICA is unable to estimate the absolute concentrations of the chromophores and it does not account for the nonlinear fluence distortion.

In this paper, we propose incorporating statistical independence as additional information in the nonlinear model-based inversion method. This involves including a measure of the independence in the error functional in addition to the least-squares data error. The aim is to reduce the quantification errors caused by inaccurate forward modelling of the fluence. The effect of using the statistical independence in the inversion is investigated using both numerical and experimental tissue mimicking phantoms and the quantification results are compared to the model-based inversion using only the data error.

## 2. Mutual information error functional

Detailed descriptions of statistical independence and its applications in imaging can be found in Refs [14,16,17]. However, as the concept has not yet been widely used in photoacoustic imaging, a brief summary of the essential aspects will be given here. The spatial distributions of certain biomarkers or molecular probes are statistically independent from other tissue chromophores in some cases of practical interest. Examples include fluorescent probes that can be targeted to disease-specific receptors whose spatial distribution is unrelated to that of the blood and the background tissue, and cells that have been genetically modified to express optical absorbers which can be found at locations independent from other absorbers. The statistical independence of two chromophores is related to the histogram of the values of the individual chromophore concentrations, and the joint histogram of the concentrations of both chromophores. The mathematical definition of statistical independence states that two random variables, **y**_{1} and **y**_{2}, are statistically independent if their joint probability density function (PDF), *ρ*_{y1,y2}(*y*_{1}, *y*_{2}), is the product of their marginal PDFs, *ρ*_{y1}(*y*_{1}) and *ρ*_{y2}(*y*_{2}) [14], such that

*y*

_{1}and

*y*

_{2}denote possible values of

**y**

_{1}and

**y**

_{2}. The degree of independence between variables can be measured using the mutual information (MI). MI is an estimate of the amount of information one variable provides on another variable. Variables with higher independence have lower MI, which means that they contain less information about each other. In quantitative photoacoustic tomography, two chromophores are considered independent if the knowledge of the concentration of one chromophore at a location does not affect the estimate of the other chromophore’s concentration at the same location. For example, if a contrast agent is independent of the blood, then the estimation of the contrast agent concentration at a voxel does not in any way predict the blood concentration there. As a counter example, oxy- and deoxyhaemoglobin are not independent chromophores: if a voxel is found to contain a high concentration of deoxyhaemoglobin, then the likelihood that a high concentration of oxyhaemoglobin will be found at the same pixel increases, as the voxel is likely to represent a blood vessel. The MI,

*I*, between

**y**

_{1}and

**y**

_{2}is given by

*ℋ*(

**y**

*) and*

_{k}*ℋ*(

**y**

_{1},

**y**

_{2}) are the entropy and the joint entropy of

**y**

_{1}and

**y**

_{2}respectively. They are defined by

*k*= 1 or 2, and

As indicated in Eqs. (2–4), the MI depends on the probabilities of the variables. We consider the concentrations of *K* independent chromophores **c**_{1}, ..., **c*** _{K}* as continuous random variables with the PDFs

*ρ*

_{c1}, ...,

*ρ*

_{cK}. However, the underlying probabilities of the chromophore concentrations are unknown and therefore

*ρ*

_{ck}needs to be estimated based on the instantiations of the random variable, which are the concentrations of the chromophore at different voxels [

*c*

_{k,1}, ...,

*c*

_{k,M}]. The total number of voxels, which is equal to the total number of instantiations, is denoted with

*M*. Using these data points, the PDF can be approximated with a kernel density estimator. Conceptually, the kernel density estimation method involves placing a smoothly varying spread of values, or kernel, on the value of each data point. The probability is then approximated by the sum of all kernels, such that

*ξ*are the values at which the PDF is evaluated and

_{k}*κ*(

*x*) is a kernel function that satisfies

*κ*(

*x*) ≥ 0 and ${\int}_{-\infty}^{\infty}\kappa (x)\text{d}x=1$. The breve symbol (˘) is used to denote the estimation of a quantity. The motivation for using the kernel density estimator instead of a simple histogram approximation for the probabilities is that the former ensures continuity and differentiability of the estimated PDFs, provided that a continuous and differentiable kernel function is used. To satisfy those criteria, the Gaussian kernel is chosen for this study: where

*h*is the kernel width, which determines the amount of smoothing in the estimations. It is calculated using the expression for the optimal kernel width normally distributed data, given by [15] where

*σ*is the standard deviation of the data.

To estimate the entropies of the chromophore concentrations, the PDF is evaluated at a set of equally spaced points denoted by [*ξ*_{k,1}, ..., *ξ*_{k,Q}], where *Q* is the total number of points. The integral in the definition of the entropy and the joint entropy for continuous variables in Eq. (3) and (4) can be approximated as a sum using the trapezium rule, such that

*q*

_{1}, ...,

*q*and Δ

_{K}*ξ*is the spacing between the sampling points for the PDF. Using the approximated entropies, the MI can be estimated by

_{k}To find the most independent chromophores requires minimising the MI between the chromophores, which can be done efficiently using a gradient-based optimisation approach. The partial derivative of the MI with respect to the chromophore concentration at each voxel is given by [18]

*κ′*denotes the derivative of

*κ*, given by

## 3. Model-based inversion with statistical independence

The first step in the model-based inversion scheme used in this study is to make an initial guess for the unknown chromophore concentrations. This initial guess and the known specific absorption spectra *α* are used to calculate the absorption coefficient, ${\mu}_{a}={\sum}_{k}{\alpha}_{k}{c}_{k}$. The fluence, Φ, is then modelled using the diffusion approximation to the radiative transfer equation for the distribution of the absorption coefficient and the scattering coefficient. The scattering coefficient is assumed to be known in this study. Using the modelled fluence, the absorption coefficient and the Grüneisen parameter Γ, the modelled initial pressure ${p}_{m,{\lambda}_{n}}^{\mathit{model}}$ at voxel *m* and wavelength *λ _{n}* is calculated based on

*S*is the system calibration factor that depends on the acoustic response of the sensors, which can be assumed to be spatially homogeneous and independent of the optical wavelength. The data error,

*ε*, is defined as the sum of squared differences between the modelled and the measured initial pressure, ${p}_{m,{\lambda}_{n}}^{\mathit{meas}}$, and hence the minimisation problem is given by

_{d}*K*,

_{t}*N*, and

*M*respectively, and

**c**

*denotes the*

_{i}*i*

^{th}chromophore. By iteratively updating the values of the chromophore concentrations until

*ε*is minimised, accurate quantification can be achieved in an idealised scenario. However, in practical applications, model-mismatch arises due to the approximations in the model and uncertainty in the model parameters (which are typically determined experimentally). This leads to the minimum of the data error occurring at erroneous chromophore concentrations, and results in inaccurate quantification. Unlike the data error, the statistical independence is a property of the distribution of the chromophores alone, rather than a function of the forward modelling. Therefore, the errors in the fluence model do not affect the MI between the chromophore concentrations, which will always have a minimum at the true solution, provided that the chromophores are statistically independent. Hence, by including a term representing the MI between the independent chromophores in the error functional, the quantification errors can be reduced. The new minimisation problem with both the data error and the MI is given by

_{d}*γ*is the weight parameter for the MI, and total number of chromophores may be larger than or equal to the number of independent chromophores,

*K*≥

_{t}*K*.

## 4. Generating multiwavelength photoacoustic images of tissue mimicking phantoms

The accuracy of the quantification using *ε*_{d+MI} compared to *ε _{d}* was investigated using experimental and numerical tissue mimicking phantoms containing aqueous solutions of copper(II) chloride dihydrate, nickel(II) chloride hexahydrate and black India ink (Winsor & Newton, London, UK), which represent different absorbers in the tissue, and Intralipid, which provides scattering in the medium. The copper(II) chloride dihydrate and nickel(II) chloride hexahydrate will be referred to as CuCl

_{2}and NiCl

_{2}. For both the numerical and the experimental phantom, the distributions of CuCl

_{2}and NiCl

_{2}are arranged such that they are statistically independent of each other.

#### 4.1. Experimentally acquired images

A schematic of the experimental set-up is shown in Fig. 1. The tissue mimicking phantom consists of four polythene tubes with 0.58mm inner diameter and 0.19mm wall thickness (Scientific Laboratory Supplies Ltd, Nottingham, UK) submerged in a background solution of highly diluted India ink and 1% (w/v) Intralipid, which give rise to an absorption and scattering amplitude comparable to that of typical biological tissue [19]. The tubes are arranged in a line at depths of approximately 3.6, 6.1, 8.1 and 9.8mm from the top surface of the phantom, which are all within the diffusive regime. The first and third tube from the top contain 399gL^{−1} NiCl_{2} and the second and the fourth tube contain 36gL^{−1} CuCl_{2}. The absorption spectra of the chromophores and scattering spectrum of Intralipid are shown in Fig. 2. The CuCl_{2} and NiCl_{2} are statistically independent of each other in this phantom, which is clear from the fact that they are contained in distinct regions that are spatially separated. (However, the spatial separation is not a necessary criterion for statistical independence. For example, CuCl_{2} and NiCl_{2} are both also statistically independent from water, even though they can be found at the same locations.)

The phantom is imaged in a V-shaped photoacoustic imaging system [24,25] consisting of two orthogonal Fabry-Perot interferometer sensors [26]. This sensor geometry increases the detection aperture of the photoacoustic signals compared to a single planar detector array and hence reduces the limited-view artefacts [25]. The fibre tip was positioned vertically above the phantom to deliver the pulsed excitation light from a Nd:YAG-pumped optical parametric oscillator (GWU, Spectra-Physics, Santa Clara, USA) with 10Hz repetition rate and a pulse energy of 15–19mJ depending on wavelength. The phantom was imaged at 8 wavelengths with equal spacing between 750nm and 890nm by recording the photoacoustic time series at a 13×13mm^{2} area with 100*μ*m step size for both sensors. A small fraction of the light was directed to an integrating sphere to measure the pulse energy, which was used to normalise the measured signals. The beam position and the spatial distribution of the beam intensity at the surface of the phantom was found by acquiring an additional photoacoustic image at 780nm of a transparent sheet with a printed grid of highly absorbing dots (with 1mm grid spacing) which was resting on the surface. Based on the reconstruction of this image, the illumination source was approximated as a Gaussian beam with a 1/e diameter of 6.6mm.

The iterative time reversal method [27, 28] was used to reconstruct 3D images from the photoacoustic time series. A 2D cross-sectional 12×12mm^{2} region of interest centred at the tubes was used for the optical inversion, and the dimension was reduced to 72×72 pixels to reduce the computational time and memory requirements. The 2D slices are shown for three wavelengths in Fig. 3.

#### 4.2. Numerically simulated images

The numerical phantom has an element spacing of 100*μ*m and represents a 5×5mm^{2} area with six insertions arranged in two columns, as illustrated in Fig. 4. The left column of insertions represents solutions of CuCl_{2} with concentrations 12, 24 and 36gL^{−1} in increasing order, where the top insertion has the lowest concentration. The right column of insertions contain solutions of NiCl_{2} with concentrations 133, 266 and 399gL^{−1}, also in increasing order from the top insertion. The phantom is designed with increasing concentrations at the insertions with depth to improve the signal to noise ratio at the deeper insertions. The absorption of CuCl_{2} and NiCl_{2} are based on the measured spectra shown in Fig. 2(a) and assumed to follow linear dependence on concentration. The concentrations are chosen such that the average absorption of both columns is 0.52mm^{−1} over the wavelength range between 750nm to 890nm, which is similar to the absorption of blood over the same range of wavelengths. Water is present in the whole phantom and the background region outside the insertions represents a solution of India ink and Intralipid, which gives rise to the same absorption and scattering amplitude as shown in Fig. 2(b–c).

The top surface of the domain was illuminated with a radially-symmetric light source with a Gaussian intensity profile with a 1/e width of 3mm. The light fluence distributions for the same 8 wavelengths as the experimental measurement in Sec. 4.1 were simulated using the diffusion approximation with the MATLAB software Toast++ [23]. The system calibration factor and the Grüneisen parameter are assumed to be known and equal to one, and the acoustic reconstruction is assumed to be perfect, such that the simulated photoacoustic images were equal to the product of the fluence and the absorption coefficient. Gaussian noise with an amplitude equal to 10% of the mean of the data, which was comparable to the magnitude of the noise in the experimental images (Sec. 4.1), was added to the simulated images.

## 5. Inverting for the chromophore concentrations

To investigate the effect of incorporating the MI term, the model-based inversion scheme was applied to the simulated and the experimentally acquired multiwavelength 2D photoacoustic images using both *ε _{d}* and

*ε*

_{d+MI}error functionals. The error functionals were minimised with the limited-memory BFGS [29] algorithm which searches for the optimal chromophore concentrations using the functional gradients of the data error term [3] and the MI term [18]. The unknown parameters in the inversion are the concentrations of CuCl

_{2}, NiCl

_{2}, India ink and water, and the known model parameters are the absorption spectra, the scattering distribution, the system calibration factor, the Grüneisen parameter and the light source position and width. The unknown chromophore concentrations were initialised with a spatially homogeneous value equal to the true concentration at the background. The iterative update was run for 300 iterations for the inversion of both the simulated and the experimental images using

*ε*or

_{d}*ε*

_{d+MI}. For all inversions using

*ε*

_{d+MI}, the MI was calculated between CuCl

_{2}and NiCl

_{2}and the weight parameter

*γ*of the MI term was set to zero for the first 200 iterations to avoid the algorithm being trapped in the local minima of the MI term. The difference in computation time for

*ε*

_{d+MI}and

*ε*was negligible. The computationally efficient calculation of the MI term and its gradient was achieved using fast Fourier transforms [30], which takes advantage of the fact that Eqs. (6) and (7) have convolution structures.

_{d}#### 5.1. Effect of model-mismatch

Two case studies were conducted using the simulated images to investigate the effect of the uncertainty in different model parameters on the quantification accuracy. In the first case, the beam diameter was set to be up to 75% smaller or larger than the true value in the inversion. In the second case, an error of up to ±75% was included in the scattering amplitude. The average errors of the estimated concentrations of CuCl_{2} and NiCl_{2} at the insertions (ROI) using the erroneous beam diameter are shown in Fig. 5(a), where the circles and asterisks correspond to the inversions using *ε _{d}* or

*ε*

_{d+MI}respectively. The results show that the increase in the percentage error in the beam diameter leads to larger quantification errors, as expected. The quantification errors for the simulated images are relatively small because only one model parameter contains error at the time. In an experimental setting, there is likely to be a combination of modelling errors, resulting in larger quantification errors. Including the MI term results in a reduction in error compared to using only the standard data error for all data points.

The inaccuracies in the scattering amplitude used in the inversion resulted in similar trends for the quantification error, as shown in Fig. 5(b). The errors are generally larger in Fig. 5(b) than 5(a), which suggests that the changes in scattering amplitude have a larger impact on the fluence distribution than changes in the beam diameter for this numerical phantom. Nonetheless, using *ε*_{d+MI} is shown to provide more accurate estimations compared to using *ε _{d}* also for this case. The relative improvement in accuracy varied between 37% and 8% with an average of 22% over all data points for both cases.

#### 5.2. Experimental results

The multiwavelength experimental images were divided by the calibration factor and the spatially varying Grüneisen parameter before the inversion. The calibration factor was determined using a forward simulation with the true concentrations. The Grüneisen parameter of aqueous solutions of CuCl_{2} and NiCl_{2} are both known to increase with concentration and are given by

_{H2O}is the Grüneisen parameter for water, and

*β*are 5.80×10

_{i}^{−3}Lg

^{−1}and 2.25×10

^{−3}Lg

^{−1}for CuCl

_{2}or NiCl

_{2}respectively [31]. In practical applications, where the true concentrations are unknown, the calibration factor can be obtained by measuring the acoustic sensitivity of the sensors [32], and the Grüneisen parameter can be included in the model as a parameter that is linearly dependent on the estimated chromophore concentrations [5,32].

The results from the model-based inversion of the experimental data are presented in Fig. 6. The estimated concentrations of CuCl_{2} (top row) and NiCl_{2} (bottom row) using *ε _{d}* and

*ε*

_{d+MI}are shown in the left and centre columns respectively in Fig. 6(a), while the true concentrations are shown in the right column. The colour scale indicates the concentrations in units of gL

^{−1}. Figure 6(b) compares the estimated with the true concentrations along a line profile across the tubes. The average estimated and expected concentrations for each tube are presented in Table 1. The inversions using

*ε*resulted in high overestimation of CuCl

_{d}_{2}in the second tube and NiCl

_{2}in the first tube, where the estimated concentrations are 94% and 149% larger than the true values respectively. There are also large cross-talk errors in the estimation of both contrast agents. This is most clearly seen for the estimated CuCl

_{2}concentration, where the third tube shows a high false-positive concentration with comparable magnitude to the concentration in the fourth tube. The accuracy of the quantification is significantly improved when the MI term is included in the error functional. The cross-talk errors for both CuCl

_{2}and NiCl

_{2}are almost completely removed when

*ε*

_{d+MI}is used. The absolute concentrations of the CuCl

_{2}is estimated accurately with an error of 3gL

^{−1}on average for the four tubes. The overestimation of the NiCl

_{2}concentration in the top tube remains present with

*ε*

_{d+MI}. However, this overestimation error is reduced when

*ε*

_{d+MI}is used compared to

*ε*.

_{d}## 6. Discussion

Minimising the MI as well as the data error led to improved quantification accuracy for both simulated and experimental multiwavelength photoacoustic images, compared to using the data error alone. In the numerical study, despite the significant fractional decrease of the quantification errors using *ε*_{d+MI}, in absolute terms, the error was decreased only by a few per cent. This was mainly due to the fact that only one model parameter was erroneous in each inversion, while all other assumptions in the model were accurate, which resulted in relatively small quantification errors, even when only *ε _{d}* was used. In the experimental study, on the other hand, a combination of different types of modelling errors was likely to have been present simultaneously, leading to poorer quantification results in the absence of the MI term. The main causes of model-mismatch may be due to the limited size of the modelled domain, which does not account for the backscattered light from outside of the domain, and the 2D modelling of the light fluence, which assumes that the light source is constant in the direction along the tubes, while in the experiment, the beam was of circular cross-section. Other possible errors may include uncertainty in the scattering spectra and amplitude, as different values have been reported in the literature [21,33]. These errors in the model affect the calculation of ${p}_{m,{\lambda}_{n}}^{\mathit{model}}$, which consequently also affect

*ε*. The MI is not affected by the fluence modelling errors, because MI is calculated based on only the distribution of the estimated chromophore concentrations in each iteration, and does not require the forward modelling of ${p}_{m,{\lambda}_{n}}^{\mathit{model}}$. Therefore, the quantification errors were greatly reduced when

_{d}*ε*

_{d+MI}was used in the experimental study. These results suggest that incorporating the statistical independence can improve the robustness of model-based inversion schemes for independent chromophores and thus potentially enhance their applicability to pre-clinical or clinical imaging studies.

In order to obtain accurate results with the inversion using *ε*_{d+MI}, it is necessary to use an appropriate weight parameter *γ* for the MI term. The weight was determined through manual trial and error. The solution was continuously dependent on the weight parameter, in the sense that small changes in the weight parameter resulted in small changes in the solution. The same weight was used for all inversions of the simulated and the experimental data, despite the differences in the data and/or the model parameters. This suggests that the concentration estimates were not highly sensitive to small variations of the weighting of the MI term around this value and, although non-trivial [34, 35], it may be possible to develop a general method for finding the optimal weight parameter for different types of applications.

The 2D fluence model based on the diffuse approximation assumes that the features are constant in the third dimension and located at depths within the diffusive regime. These assumptions are appropriate for the phantom geometry used in this study. However, full 3D fluence modelling will be required for applications of the model-based inversion in biological tissue with complex structures. More accurate modelling of the fluence for the superficial layer can be achieved by incorporating the delta-Eddington approximation [5] or using the radiative transfer equation [36]. The calculation of the MI can be straightforwardly extended to 3D without causing significant increase in computation time. Furthermore, using the MI term in the error functional is also compatible with other regularisation methods such as the total-variation regulariser [37–39].

## 7. Conclusion

We proposed exploiting the statistical independence between certain chromophores in the model-based inversion method by minimising the MI between the independent chromophores in addition to the data error. The improvement in the accuracy of the estimated chromophore concentrations was demonstrated using both numerical simulations and an experimental phantom. The results suggest that the sensitivity of the model-based inversion to model-mismatch can be reduced by incorporating the additional information of statistical independence. Thus, the robustness and hence usefulness of the inversion scheme can potentially be improved for *in vivo* imaging experiments.

## Funding

This work was funded by the UCL EPSRC Centre for Doctoral Training in Medical Imaging.

## Disclosures

The authors declare that there are no conflicts of interest related to this article.

## References and links

**1. **P. Beard, “Biomedical photoacoustic imaging,” Interface Focus **1**, 602–631 (2011). [CrossRef]

**2. **J. Weber, P. C. Beard, and S. E. Bohndiek, “Contrast agents for molecular photoacoustic imaging,” Nat. Methods **13**, 639–650 (2016). [CrossRef] [PubMed]

**3. **B. T. Cox, S. R. Arridge, and P. C. Beard, “Estimating chromophore distributions from multiwavelength photoacoustic images,” J. Opt. Soc. Am. A **26**(2), 443–455 (2009). [CrossRef]

**4. **Y Sun, E Sobel, and H Jiang, “Quantitative three-dimensional photoacoustic tomography of the finger joints: an in vivo study,” J. Biomed. Opt. **14**(6), 064002 (2009). [CrossRef]

**5. **J. Laufer, B. Cox, E. Zhang, and P. Beard, “Quantitative determination of chromophore concentrations from 2D photoacoustic images using a nonlinear model-based inversion scheme,” Appl. Opt. **49**(8), 1219–1233 (2010). [CrossRef] [PubMed]

**6. **G. Bal and K. Ren, “On multi-spectral quantitative photoacoustic tomography in diffusive regime,” Inverse Probl. **28**(2), 025010 (2012). [CrossRef]

**7. **T. Tarvainen, B. T. Cox, J. P. Kaipio, and S. R. Arridge, “Reconstructing absorption and scattering distributions in quantitative photoacoustic tomography,” Invserse Prob. **28**(8), 084009 (2012). [CrossRef]

**8. **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. **36**(1), 322–331 (2017). [CrossRef]

**9. **A. Hyvarinen and E. Oja, “Independent component analysis: algorithms and applications,” Neural Netw. **13**, 411–430 (2000). [CrossRef] [PubMed]

**10. **J. Glatz, N. C. Deliolanis, A. Buehler, D. Razansky, and V. Ntziachristos, “Blind source unmixing in multi-spectral optoacoustic tomography,” Opt. Express **19**(4), 3175–3184 (2011). [CrossRef] [PubMed]

**11. **N. C. Deliolanis, A. Ale, S. Morscher, N. C. Burton, K. Schaefer, K. Radrich, D. Razansky, and V. Ntziachristos, “Deep-tissue reporter-gene imaging with fluorescence and optoacoustic tomography: A performance overview,” Mol. Imaging Biol. **16**(5), 652–660 (2014). [CrossRef] [PubMed]

**12. **A. Buehler, E. Herzog, A. Ale, B. D. Smith, V. Ntziachristos, and D. Razansky, “High resolution tumor targeting in living mice by means of multispectral optoacoustic tomography,” EJNMMI Research **2**(1), 1–6 (2012). [CrossRef]

**13. **A. Taruttis, M. Wildgruber, K. Kosanke, N. Beziere, K. Licha, R. Haag, M. Aichler, A. Walch, E. Rummeny, and V. Ntziachristos, “Multispectral optoacoustic tomography of myocardial infarction,” Photoacoustics **1**(1), 3–8 (2013). [CrossRef] [PubMed]

**14. **A. Hyvarinen, J. Karhunen, and E. Oja, *Independent Component Analysis* (John Wiley & Sons, 2004), Chap. 2.

**15. **B. W. Silverman, *Density Estimation for Statistics and Data Analysis* (Chapman and Hall1986). [CrossRef]

**16. **J. Wang and C. I. Chang, “Applications of independent component analysis in endmember extraction and abundance quantification for hyperspectral imagery,” IEEE Trans. Geosci. Remote Sens. **44**(9), 2601–2616 (2006). [CrossRef]

**17. **V. D. Calhoun, J. Liu, and T. Adal, “A review of group ICA for fMRI data and ICA for joint inference of imaging, genetic, and ERP data,” NeuroImage **45**(1, Supplement 1), S163–S172 (2009). [CrossRef]

**18. **C. Panagiotou, S. Somayajula, A. P. Gibson, M. Schweiger, R. M. Leahy, and S. R. Arridge, “Information theoretic regularization in diffuse optical tomography,” J. Opt. Soc. Am. **26**(5), 1277–1290 (2009). [CrossRef]

**19. **S. L. Jacques, “Optical properties of biological tissues: a review,” Phys. Med. Biol. **58**, R37–R61 (2013). [CrossRef] [PubMed]

**20. **L. Kou, D. Labrie, and P. Chylek, “Refractive indices of water and ice in the 0.65-to 2.5-μm spectral range,” Appl. Opt. **32**(19), 3531–3540 (1993). [CrossRef] [PubMed]

**21. **H. J. Van Staveren, C. J. Moes, J. van Marie, S. A. Prahl, and M. J. C. van Gemert, “Light scattering in lntralipid-10% in the wavelength range of 400–1100nm,” Appl. Opt. **30**(31), 4507–4514 (1991). [CrossRef] [PubMed]

**22. **L. An, T. Saratoon, M. Fonseca, R. Ellwood, and B. Cox, “Exploiting statistical independence for quantitative photoacoustic tomography,” Proc. SPIE **10064**, 1006419 (2017). [CrossRef]

**23. **M. Schweiger and S. Arridge, “The Toast++ software suite for forward and inverse modeling in optical tomography,” J. Biomed. Opt. **19**(4), 040801 (2014). [CrossRef] [PubMed]

**24. **R. Ellwood, O. Ogunlade, E. Z. Zhang, P. C. Beard, and B. T. Cox, “Orthogonal Fabry-Pérot sensors for photoacoustic tomography,” Proc. SPIE **9708**, 97082N (2016).

**25. **R. Ellwood, O. Ogunlade, E. Zhang, P. Beard, and B. Cox, “Photoacoustic tomography using orthogonal Fabry-Pérot sensors,” J. Biomed. Opt. **22**(4), 041009 (2017). [CrossRef]

**26. **E. Zhang, J. Laufer, and P. Beard, “Backward-mode multiwavelength photoacoustic scanner using a planar Fabry-Perot polymer film ultrasound sensor for high-resolution threedimensional imaging of biological tissues,” Appl. Opt. **47**, 561–577 (2008). [CrossRef] [PubMed]

**27. **G. Paltauf, J. A. Viator, S. A. Prahl, and S. L. Jacques, “Iterative reconstruction algorithm for optoacoustic imaging,” J. Acoust. Soc. Am. **112**(4), 1536–1544 (2002). [CrossRef] [PubMed]

**28. **P. Stefanov and G. Uhlmann, “Thermoacoustic tomography with variable sound speed,” Inverse Probl. **25**(7), 075011 (2009). [CrossRef]

**29. **J. Nocedal and S. Wright, *Numerical Optimization* (Springer Science & Business Media, 2006).

**30. **S. Shwartz, M. Zibulevsky, and Y. Y. Schechner, “ICA using kernel entropy estimation with NlogN complexity,” in *Independent Component Analysis and Blind Signal Separation: Fifth International Conference, ICA 2004, Granada, Spain, September 22–24, 2004. Proceedings*, C. G. Puntonet and A. Prieto, eds. (Springer, 2004), pp. 422–429. [CrossRef]

**31. **J. Laufer, E. Zhang, and P. Beard, “Evaluation of absorbing chromophores used in tissue phantoms for quantitative photoacoustic spectroscopy and imaging,” IEEE J. Sel. Top. Quantum Electron. **16**(3), 600–607 (2010). [CrossRef]

**32. **M. Fonseca, E. Malone, F. Lucka, R. Ellwood, L. An, R. Arridge, P. Beard, and B. Cox, “Three-dimensional photoacoustic imaging and inversion for accurate quantification of chromophore distributions,” Proc. SPIE **10064**, Photons Plus Ultrasound: Imaging and Sensing 2017, 1006415 (2017). [CrossRef]

**33. **S. T. Flock, S. L. Jacques, B. C. Wilson, W. M. Star, and M. J. van Gement, “Optical properties of intralipid: A phantom medium for light propagation studies,” Lasers. Surg. Med. **12**(5), 510–519 (1992). [CrossRef] [PubMed]

**34. **C. R. Vogel, *Computational Methods for Inverse Problems* (Society for Industrial and Applied Mathematics (SIAM): Philadelphia, PA: 2002). Chap. 7. [CrossRef]

**35. **A. M. Thompson, J. C. Brown, J. W. Kay, and D. M. Titterington, “A study of methods of choosing the smoothing parameter in image restoration by regularization,” IEEE Trans. Pattern Anal. Mach. Intell. **13**(13),326–339 (1991). [CrossRef]

**36. **T. Saratoon, T. Tarvainen, B. T. Cox, and S. R. Arridge, “A gradient-based method for quantitative photoacoustic tomography using the radiative transfer equation,” Invserse Prob. **29**(7), 075006 (2013). [CrossRef]

**37. **L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Physica D. **60**(1–4), 259–268 (1992). [CrossRef]

**38. **C. Huang, K. Wang, L. Nie, L. Wang, and M. Anastasio, “Full-wave iterative image reconstruction in photoacoustic tomography with acoustically inhomogeneous media,” IEEE Trans. Med. Imaging. **32**(6), 1097–1110 (2013). [CrossRef] [PubMed]

**39. **S. Arridge, P. Beard, M. Betcke, B. Cox, N. Huynh, F. Lucka, O. Ogunlade, and E. Zhang, “Accelerated high-resolution photoacoustic tomography via compressed sensing,” Phys. Med. Biol. **61**(24), 8908 (2016). [CrossRef] [PubMed]