Considering sources and detectors distributions for quantitative photoacoustic tomography

Photoacoustic tomography (PAT) is a hybrid imaging modality that takes advantage of high optical contrast brought by optical imaging and high spatial resolution brought by ultrasound imaging. However, the quantification in photoacoustic imaging is challenging. Multiple optical illumination approach has proven to achieve uncoupling of diffusion and absorption effects. In this paper, this protocol is adopted and synthetic photoacoustic data, blurred with some noise, were generated. The influence of the distribution of optical sources and transducers on the reconstruction of the absorption and diffusion coefficients maps is studied. Specific situations with limited view angles were examined. The results show multiple illuminations with a wide field improve the reconstructions. ©2014 Optical Society of America OCIS codes: (170.0170) Medical optics and biotechnology; (170.5120) Photoacoustic imaging. References and links 1. S. A. Ermilov, T. Khamapirad, A. Conjusteau, M. H. Leonard, R. Lacewell, K. Mehta, T. Miller, and A. A. Oraevsky, “Laser optoacoustic imaging system for detection of breast cancer,” J. Biomed. Opt. 14(2), 024007 (2009). 2. M. Heijblom, D. Piras, W. Xia, J. C. G. van Hespen, J. M. Klaase, F. M. van den Engh, T. G. van Leeuwen, W. Steenbergen, and S. Manohar, “Visualizing breast cancer using the Twente photoacoustic mammoscope: What do we learn from twelve new patient measurements?” Opt. Express 20(11), 11582–11597 (2012). 3. D. Razansky, C. Vinegoni, and V. Ntziachristos, “Multispectral photoacoustic imaging of fluorochromes in small animals,” Opt. Lett. 32(19), 2891–2893 (2007). 4. C. Li and L. V. Wang, “Photoacoustic tomography and sensing in biomedicine,” Phys. Med. Biol. 54(19), R59– R97 (2009). 5. H. F. Zhang, K. Maslov, and L. V. Wang, “In vivo imaging of subcutaneous structures using functional photoacoustic microscopy,” Nat. Protoc. 2(4), 797–804 (2007). 6. X. D. Wang, Y. J. Pang, G. Ku, X. Y. Xie, G. Stoica, and L. V. Wang, “Noninvasive laser-induced photoacoustic tomography for structural and functional in vivo imaging of the brain,” Nat. Biotechnol. 21(7), 803–806 (2003). 7. H. F. Zhang, K. Maslov, M. Sivaramakrishnan, G. Stoica, and L. V. Wang, “Imaging of hemoglobin oxygen saturation variations in single vessels in vivo using photoacoustic microscopy,” Appl. Phys. Lett. 90, 053901 (2007). 8. J. Kottmann, J. M. Rey, J. Luginbühl, E. Reichmann, and M. W. Sigrist, “Glucose sensing in human epidermis using mid-infrared photoacoustic detection,” Biomed. Opt. Express 3(4), 667–680 (2012). 9. J. Ripoll and V. Ntziachristos, “Quantitative point source photoacoustic inversion formulas for scattering and absorbing media,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 71(3), 031912 (2005). 10. B. T. Cox, S. R. Arridge, K. P. Köstli, and P. C. Beard, “Two-dimensional quantitative photoacoustic image reconstruction of absorption distributions in scattering media by use of a simple iterative method,” Appl. Opt. 45(8), 1866–1875 (2006). 11. B. Banerjee, S. Bagchi, R. M. Vasu, and D. Roy, “Quantitative photoacoustic tomography from boundary pressure measurements: noniterative recovery of optical absorption coefficient from the reconstructed absorbed energy map,” J. Opt. Soc. Am. A 25(9), 2347–2356 (2008). 12. Y. Zhou, J. Yao, K. I. Maslov, and L. V. Wang, “Calibration-free absolute quantification of particle concentration by statistical analyses of photoacoustic signals in vivo,” J. Biomed. Opt. 19(3), 037001 (2014). 13. J. R. Rajian, P. L. Carson, and X. Wang, “Quantitative photoacoustic measurement of tissue optical absorption spectrum aided by an optical contrast agent,” Opt. Express 17(6), 4879–4889 (2009). 14. K. Daoudi, A. Hussain, E. Hondebrink, and W. Steenbergen, “Correcting photoacoustic signals for fluence variations using acousto-optic modulation,” Opt. Express 20(13), 14117–14129 (2012). (C) 2014 OSA 1 November 2014 | Vol. 5, No. 11 | DOI:10.1364/BOE.5.003960 | BIOMEDICAL OPTICS EXPRESS 3960 #214301 $15.00 USD Received 19 Jun 2014; revised 23 Sep 2014; accepted 13 Oct 2014; published 16 Oct 2014 15. A. Q. Bauer, R. E. Nothdurft, T. N. Erpelding, L. V. Wang, and J. P. Culver, “Quantitative photoacoustic imaging: correcting for heterogeneous light fluence distributions using diffuse optical tomography,” J. Biomed. Opt. 16(9), 096016 (2011). 16. A. Rosenthal, D. Razansky, and V. Ntziachristos, “Quantitative optoacoustic signal extraction using sparse signal representation,” IEEE Trans. Med. Imaging 28(12), 1997–2006 (2009). 17. G. Bal and K. Ren, “On multi-spectral quantitative photoacoustic tomography in diffusive regime,” Inverse Probl. 28(2), 025010 (2012). 18. G. Bal and G. Uhlmann, “Inverse diffusion theory of photoacoustics,” Inverse Probl. 26, 085010 (2010). 19. T. Tarvainen, B. T. Cox, J. P. Kaipio, and S. R. Arridge, “Reconstructing absorption and scattering distributions in quantitative photoacoustic tomography,” Inverse Probl. 28(8), 084009 (2012). 20. G. Bal and K. Ren, “Multiple-source quantitative photoacoustic tomography,” Inverse Probl. 27, 075003 (2011). 21. P. Shao, B. Cox, and R. J. Zemp, “Estimating optical absorption, scattering, and Grueneisen distributions with multiple-illumination photoacoustic tomography,” Appl. Opt. 50(19), 3145–3154 (2011). 22. P. Shao, T. Harrison, and R. J. Zemp, “Iterative algorithm for multiple illumination photoacoustic tomography (MIPAT) using ultrasound channel data,” Biomed. Opt. Express 3(12), 3240–3249 (2012). 23. T. Vo-Dinh, Biomedical Photonics Handbook (CRC Press, 2003). 24. V. Ntziachristos, A. G. Yodh, M. Schnall, and B. Chance, “Concurrent MRI and diffuse optical tomography of breast after indocyanine green enhancement,” Proc. Natl. Acad. Sci. U.S.A. 97(6), 2767–2772 (2000). 25. A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging (IEEE Press, NY, 1987). 26. H. Gao, H. Zhao, and S. Osher, “Bregman methods in quantitative photoacoustic tomography,” Univ. Calif. Los Angel. UCLA Comput. Appl. Math. Rep. Vol. 10 – 42, (2010). 27. T. Saratoon, T. Tarvainen, B. Cox, and S. Arridge, “A gradient-based method for quantitative photoacoustic tomography using the radiative transfer equation,” Inverse Probl. 29(7), 075006 (2013).


Introduction
Photoacoustic tomography (PAT) is a multiwave imaging modality that potentially combines the high optical contrast brought by the optics and high spatial resolution by ultrasound probing.The interests in using this technique for biomedical sensing are growing since 1990s, and its applications are in breast cancer detection [1,2], small animal molecular imaging [3,4], subcutaneous structures imaging [5], functional imaging [6], and other applications.
The principle consists in illuminating a sample with a local absorbing optical contrast (vasculature, presence of a tumor for example).The absorbed energy converts into heat and the temperature of the object increases, causing thermal expansion that generates an acoustic pressure wave.The energy deposited is proportional to the fluence and the local absorption coefficient of the tissues.Hence, the primary aim of PAT is to retrieve the spatial 3D distribution of the absorption coefficient of the probed tissue.This absorption coefficient, proportional to the concentration of the different chromophores composing the tissue, is an intrinsic marker of the physiology and the metabolism.In the red and near infrared, PAT has been applied to mapping hemoglobin concentrations that can reveal the presence of a cancer, through oxy-and deoxy-hemoglobin concentrations which are markers of the oxygen consumption in the tissues [7].Glucose monitoring with PAT has also been reported [8].However, usually the technique is applied to probing tissues at small depths (small animals, subcutaneous examination), for which the fluence can be modelled as a Beer-Lambert's law by considering light scattering is not perturbing too much the quantification.For deep tissues imaging, because of their high level of scattering, diffusion of light has to be taken into account in the evaluation of the fluence that depends non-linearly on both the absorption and diffusion coefficients.Isolating the absorption contribution from the energy deposited is the challenge in quantitative PAT (QPAT), as the fluence does depend itself on the optical parameters (absorption and diffusion coefficients) distributions.
In a first approximation, PAT can be limited to a purely acoustic problem: from timeresolved measurements collected at the periphery, acoustic sources are localized through the reconstruction of the initial pressure distribution, from which one can obtain the absorption coefficient distribution by making the assumption of a simple photon propagation model or a uniform distribution of light.The problem is similar to seismology passive sensing: one seeks at localizing the sources that produced the collected acoustic signal.For all these approaches, a wide field illumination is preferable in order to produce a homogenized initial pressure.To go further in the quantification, different strategies have been adopted taking advantage of specific experimental protocols.
There exist basically four approaches in QPAT which consists in recovering: i. the absorption coefficient directly.Supposing the diffusion coefficient is known, Ripoll et al. [9] recovered the absorption map with a point source accounting for instrumental factors such as the source strength, the shape of the optical pulse, and the impulse response and finite size of the transducers; Cox et al. [10] reconstructed the absorption coefficient assuming the scattering coefficient is known a priori by use of a fixed-point iterative method; Banerjee et al. [11] recovered the absorption coefficient directly from boundary pressure measurements.Very recently, Zhou et al. [12] have introduced a calibration-free method based on the direct measurement of the counts of particles through the statistical analysis of the photoacoustic signal fluctuations.
ii. the absorption coefficient and the fluence.The fluence contribution can be filtered by making use of a contrast agent of known absorption coefficient with photoacoustic signal measurements performed before and after injection [13].Daoudi et al. [14] proposed to combine photoacoustics with acousto-optics allowing also to overcome the contribution of the fluence in the initial pressure measurement.The fluence map can be explicitly measured through Diffuse Optical Tomography [15] (DOT).By introducing prior knowledge on the spatial frequency of the two quantities, Rosenthal et al. [16] extracted both the absorption coefficient (high frequency) and the fluence (low frequency) from the photoacoustic image.As this algorithm is not based on the explicit solution the theoretical light transport equation, it does not require explicit knowledge of the illumination geometry.
iii. the absorption and diffusion coefficients simultaneously by introducing prior knowledge (spectra of known species, diffusion spectral profile) with multiple wavelength illuminations; Bal and Ren [17] have proven that, when multiple wavelength data are available, the absorption, diffusion and Grüneisen coefficients can be reconstructed simultaneously under additional minor prior assumptions.
iv. the absorption and diffusion coefficients simultaneously by introducing active probing with multiple optical illuminations.In the last decade, a lot of efforts have been produced to improve the quantification through the uncoupled reconstructions of the different physical parameters contained in the initial pressure distribution map, that is: absorption and diffusion coefficients and the Grüneisen coefficient accounting self-consistently on the acoustic properties of the tissue [18].Tarvainen et al. [19] reconstructed both the absorption and diffusion coefficients, by considering both the radiative transport and diffusion equations as light transport models; Bal and Ren [20] recovered the absorption, diffusion, and Grüneisen coefficient mathematically when the propagation of radiation is modelled by a second-order elliptic equation; Shao et al. [21] estimated the absorption, scattering and Grüneisen distributions from reconstructed initial pressure map with multiple optical sources; In their updated work [22], they recovered the absorption and diffusion coefficients directly from measured acoustic pressure.With this last approach, the probing becomes active and similar with what occurs in geophysics when searching for buried objects, petroleum layers, cracks...In the present work, this approach is adopted and tested under different sources and detectors geometries.The purpose is to test its robustness under potential experimental situations with single wavelength illumination, point scanning single element transducer and limited view angle examination.The paper is organized as follows: Section 2 presents the theoretical background, the results are presented and discussed in Section 3, Section 4 is the conclusion.

Theory
Hereafter are recalled the definitions of the physical quantities involved in the photoacoustic phenomenon and the main lines leading to the expression of the solution of the forward model.The inverse problem is also presented.It was formulated according to the method proposed by Shao et al. [22] for multiple illuminations and detections situation.

Photoacoustic phenomenon modelling in biological tissues
The amount of heat generated by tissue is proportional to the strength of the radiation.Within the conditions of thermal and stress confinements, light pulse is typically few ns, the heating function can be treated as an instantaneous function proportional to a Dirac () t δ function: Here (,) H t r is the heating function defined as the thermal energy converted at spatial position r and time t by the electromagnetic radiation per unit volume per unit time.() is the absorption coefficient.When large tissues are considered, with illumination in the wavelength range (600-1200 nm) for which light scattering mean free path length is much shorter than the absorption mean free path length, the Diffusion Approximation holds [23] and the spatial distribution of light fluence can be obtained through the resolution of the following diffusion equation: , can be formulated linearly with the following matrix form: and the Jacobian matrix: β ur represents the sensitivity of the energy instantaneously deposited due to the optical perturbation δ u located at positions j r , j r belonging to the reconstruction mesh: rru r r ru r r ru rru r r ru r r (9) In order to avoid infinite values of 0 G and φ at positions j l = rr , the reconstruction mesh pixel surface) for the values of the optical properties was chosen to be slightly different from the modelling mesh ( ).

Computation of the Jacobian matrix
In the following case studies, the simulations were conducted in 2D.The chosen test object (Fig. 1) was basically the same as in [22]: a 6 cm × 6 cm square, in the center area of which a 2 cm × 2 cm square was selected as the reconstruction area (Fig. 1, left).The background absorption and diffusion coefficients are respectively (dimensions: 0.3cm × 1.1cm and 0.3cm × 0.3cm, 0.6 cm from the center).The perturbations were chosen relatively small (10%) in order to fulfil the Born approximation.The values correspond to optical properties measured in various situations with DOT (for instance in breast cancer imaging [24]).If the perturbations are high compared to the background optical properties, non-linear reconstruction algorithms have to be used [22].The sources are located 0.3 cm away from the reconstruction area, both the number and the length of the sources are variable.Point detectors are located 2 cm from this computation area, right at the periphery of the object.).The FEM mesh typically contains 4300 elements, the modelling mesh is a regular grid of 61 × 61 pixels (pixel size: 0.1 cm), and the reconstruction mesh is 21 × 21 pixels, 1 pixel shifted in horizontal and vertical with the modelling mesh.,0 (() ) β ur is then obtained by assembling the sensitivity matrices, and its dimension is () 2 SL J ×× .α is ( ) The dimension of the Jacobian matrix J is large () 2 TMS J ××× but the Hessian H dimension is reduced 22 J J × .Vector Δ representing the discrepancy between model and measurements has dimension () 1 TMS ××× and b is only 21 J × .The largest matrices J and Δ are computed once and are not involved in the reconstruction process.

Resolution of the inverse problem
Two different reconstruction methods were tested: -Singular value decomposition (SVD) method for which the Hessian matrix was decomposed into: T =× × HUSV .S is a diagonal matrix with same dimension as H containing nonnegative singular values i σ ordered in decreasing order, U and V are unitary matrices.The inversion of (5) becomes: -The algebraic reconstruction technique (ART) [25] that searches iteratively for the solution δ u such that: hyperparameter that can be tuned in order to regularize or to speed up the reconstruction.In all what follows, n was kept constant and equal to 1. Non negative values constraint is also introduced in the algorithm.
These reconstruction methods are tested for different noise levels.For small noise levels, SVD provides sufficiently good results but is known to be unstable and very sensitive to noise, while the ART behaves better, especially when constrained.

Simulation study on sources and detectors geometries
Sets of synthetic data were generated by solving the forward problem and introducing a fraction of random noise with normal distribution.Hereafter are presented the results of simulations performed on the synthetic phantom (Fig. 1) with the above-described method.The purpose of the study is to highlight the effects of different sources and detectors distributions and geometries on the reconstruction of the perturbation of the optical parameters [] T a δδ δ = u D , in terms of accuracy and robustness to noise.To that end, the following aspects were considered: (i) using different reconstruction methods with different noise level: this preliminary step allowed to select the proper reconstruction algorithm; (ii) number of point sources distribution: the point source illumination has been specifically studied in [22], here the results are shown as reference; (iii) point detectors distributions: the purpose is to understand which transducers distribution would be the best, for a given number of measurements; (iv) using extended sources instead of the point sources: point sources have been studied but, in practice, extended sources are easier to handle experimentally; (v) limited angle examination: this corresponds to actual experimental situations, here the transmission geometry, where sources and detectors belong to different half-spaces, is studied.

Different reconstruction methods
Four point sources and sixty detectors (15 detectors on each side of the object) are considered for this study.A first point that can be stressed is the poor conditioning of the Hessian: after performing the SVD, and examining the spectrum of the singular values i σ of H, the condition number obtained for this configuration is As expected, reconstructions are highly sensitive to noise, a high value of the regularization parameter has to be introduced in the SVD, showing only few elements of H out of noise may be useful for the reconstruction, even at low noise level.With the introduction of non-negative constraint, ART seems to have better potential in the recovery of the optical parameters at higher noise levels.It is interesting to notice that, though ART reconstructs the perturbations with a higher contrast.In all what follows, ART was chosen as the reconstruction algorithm.
Hence, the aim of the following studies is to improve the conditioning of H such that the reconstructions are less sensitive to noise, by optimising the data acquisition geometry.

Sources distributions
In this study, 60 point detectors were distributed around the object and the sources number was varied (2, 4, 8 and 16 sources considered).Examining the condition number of H as a function of the number of sources (Fig. 3), the reconstructions should be indeed more stable as the number of sources increases.Although the inversion is mathematically demonstrated, 2 sources are not enough to obtain uncoupled robust reconstructions.
,, / . The error in the reconstruction values decreases from 1 with 2 sources to 0.4 with 4.An improvement of a factor 10 in the condition number (Fig. 3, Right Top) is obtained by adding only two more sources.However, the improvement in increasing the number of sources might be reduced after 8 sources and the gain in the reconstructions might be small even if the number of measurements (and hence the duration of the examination) is increased.This is all-the-more visible in the cross-plots shown in Fig. 3 Right Bottom: the values reconstructed become closer to the target values, the edges of the objects are sharper and the reconstruction artifacts decrease.

Detectors distributions
Four sources and 60 point detectors are now considered, the number of measurements is now kept constant but the distribution of the detectors is varied.Five situations, schematized in Fig. 4, were examined: a) even and sparse distribution around the object; b) 4-sides examination but short spacing between detectors; c) same situation but with even shorter (compressed) spacing; d) 2-sides examination horizontal detection; e) 2-sides examination vertical detection.
Examining the reconstructions (Fig. 4 Top Left), the quadratic errors (Fig. 4 Bottom Left) and the conditioning of H (Fig. 4 Top Right) obtained for these five situations shows no major improvement in the reconstructed images.Through the values of the condition number and the quadratic errors, one may conclude that distributing the detectors evenly (Fig. 4 situation (a)) may be the most favorable situation.However, reducing the field of view of the detection to the area of interest (Fig. 4 situations (b) and (c)) seems to improve the quantification while not degrading to much the localization.Figure 4 situations (d) and (e) show explicitly the influence of the positioning of the detectors in comparison to the object geometry, especially in the rectangular absorption perturbation: situation (d), when detectors are probing the object through its thinnest width, the emitted acoustic signal is better defined in time and higher in frequency (small objects producing sharper time-resolved acoustic signals), while in situation (e), the resolution is much more degraded in the direction parallel to the detectors (crossplots, Fig. 4 Right Bottom).

Point sources versus wide field illumination
QPAT has been demonstrated by using point sources illumination [22].However, in practice, a wide field illumination is easier and safer to handle experimentally because the irradiation dose can be decreased by distributing the energy over a larger area.For this study, four line sources illuminating the four sides of the reconstruction area were considered (Fig. 5).Their lengths were varied from zero, corresponding to the point source situation, up to the maximum length allowing a full angle probing of the reconstruction area.Increasing the length of the source brings major improvements in the conditioning of H (Fig. 5 Right Top) and, hence, in the reconstructions (Fig. 5 Left Top).These results show that the quantification (see QE Fig. 5 Left Bottom) and localization (see crossplots Fig. 5 Right Bottom) are both definitively improved by using wider field illumination.Hence, the number of measurements required for obtaining accurate reconstructions in terms of localization and quantification can be reduced.When comparing reconstruction results shown in Fig. 3 for situation (c) and those in Fig. 5 for situations (b) to (e), with two times less measurements, the reconstructions obtained are better with 4 line sources illuminations than with 8 point sources (QEs are sensibly the same).However, the medium has still to be probed under multiple points of views: as shown for comparison on Fig. 5 situation (f), using one single wide source does not bring sufficient information to discriminate absorption and diffusion.In this last test case, the absorption coefficient is well reconstructed while the diffusion coefficient is far underestimated, showing that this illumination protocol senses essentially absorption abnormalties (see cross-plots on Fig. 5 Right Bottom).

Example of simulations for experimental situations with limited angle examination
Examining the object under all angles may be experimentally complicated or simply not feasible.Here are reported reconstructions obtained under constrained examination geometry: the medium is probed exclusively under transmission geometry, sources and detectors can rotate and probe the object on its four sides.The purpose of these calculations is to test if satisfying reconstructions can be obtained under less favorable situations corresponding nevertheless to simpler experimental protocols.Figure 6 Top Left shows the five different synthetic experiments considered here.For all of them, the four sides of the object are probed in transmission, sequentially (this corresponds to 4 measurements with an object rotating around its center of gravity with rotation angles 0°, 90°, 180° and 270°), with (a) 1 point source illumination on one side and 15 point detectors evenly distributed on the other side; from this reference situation, the illumination source was enlarged, (b) 1 line source (length: 26 mm) and 15 detectors; then the number of detectors was reduced up to a very limited number (c) same line source and 1 detector only (in this situation, sources and detectors are interchanged compared to situation (a)); (d) same as (c) but with 2 point detectors; (e) same as (c) with 3 point detectors.The condition numbers and quadratic errors obtained for the corresponding reconstructions are also reported Fig. 6.
A first remark is the reconstructions obtained with point or line sources illuminations (Fig. 6 situations (a) and (b)) are improved compared to those for which the detectors are collecting signal from all the sides of the object: when the object is probed in transmission, although the condition numbers are several orders of magnitude higher, still one can get reconstructions of satisfying quality (QE≤0.4) with highly reduced number of measurements (60 measurements only).Figure 6(c) is the oversimplified situation, a single wide source and a single point detector: the absorption coefficient is satisfyingly reconstructed while the diffusion coefficient is not.However, slightly increasing the number of detectors (2 or 3 point detectors per view Fig. 6(d) and 6(e)) allows recovering reconstructions of improved quality.
This result shows the possibility of using simplified experimental PAT setups.Reducing the number of detectors points, up to 3 only here, corresponding to 12 measurements only, of course worsens the quality of the reconstructions but still allows getting better results than with complete angle detection: in comparison with Fig. 5(e), obtained with 240 measurements, the number of measurements is reduced by a factor 20. Additionally, one may think in combining the measurements in order to improve the quality of the reconstructions: the absorption coefficient seems to be fairly well reconstructed with the simplest situation (1 source, 1 detector, Fig. 6 (c)), these results could be used as prior knowledge in a second step reconstruction using 2 or 3 additional measurements (Fig. 6(d) or 6(e)).

Conclusion
In the present work, a multiple illumination photoacoustic tomography algorithm has been implemented according to the method suggested in [22].The method is quantitative in the sense that it allows reconstructing both the absorption and diffusion coefficient while the conventional approaches in PAT, with a single wide field illumination give the access to the initial pressure distribution map.A description with explicit practical implementation details was provided.Different linear reconstruction algorithms were tested: the ART with nonnegative constraints was shown to be less sensitive to noise than SVD or LSQR.Even if not involved in these reconstruction loops, the computation of the Jacobian matrix is still a burden, especially if one seeks to extend the approach to the non-linear case.Gradient-based methods may be used instead because more memory-efficient as they require only gradient information to approximate the Hessian matrix [17,26,27].
With this tool, a 2D simulation study on a variety of sources and detectors geometries situations was conducted.The sets of simulations were conducted through the examination of the conditioning of the Hessian matrix.The reconstructions highlighted several important results: -Optical parameters are reconstructed with improved quality by using wide field illuminations: this study demonstrates that it is preferable to use enlarged sources instead of multiplying the number of point sources illuminations.
-The position of the detectors with respect to the perturbations to be reconstructed is crucial and deserves additional prospections.When comparing the influence of sources and detectors, the role of the illumination and detection schemes are nested and complex.The illumination scheme seems to play a major role in the quantification, while the detection in the localization with high resolved time-offlight measurement.An important result is that probing the sample under transmittance geometry only provides better reconstructions than with full angular coverage detection.
-Different simple experimental situations were tested and good quality reconstructions were obtained with a number of measurements extremely reduced.Experimental validations may support these findings that may influence future experimental setups designs.The present study has highlighted the importance of the data acquisition protocol.A smart combination between different types of measurements will undoubtedly improve the quality of the reconstructions.An example of such a process would be to combine one measurement taken under a full angle illumination, with a field of illumination as wide as possible, and a set of data acquired under different illumination angles.The first measurement would serve in a first reconstruction loop that would provide initialization and constraint values for a second reconstruction loop, performed with the second series of data.
the Hessian.The multiplication by T J allows handling an inversion process with smaller dimensions of the matrices involved (reduced to the dimension of the chosen reconstruction mesh).The main problem remains in expressing J .The medium is first meshed into L voxels of volume rr , the model for the considered measurable quantity at time evolution transfer matrix propagating the pressure created at initial time at position l r to point detector located at d r ; ,0 (() ) ls j

δG
is the Kronecker symbol.Within the perturbation approach, knowing the Green's function 0 of the diffusion Eq. (

2. 3 . 2 α
Time evolution transfer matrix propagating the pressure , dl τ This matrix represents the acoustic propagation between each pixel of the modelling mesh and the transducer.It carries essentially the information on the distance between detectors and the different pixels of the modelling mesh dl d l R =− rr .It is calculated by considering the points located at the center of each pixel of the modelling mesh { } l r .l V Δ is the surface of each pixel, The Grüneisen coefficient is chosen to be a classical one in biological tissues 0steps are considered with a time interval of 0.2 µs.The dimension of , dl τ

1 −
b .A Tikhonov regularization was adopted to avoid instabilities caused by poor conditioning of H : the values of the diagonal of are expected to be unstable and highly sensitive to noise.The two different reconstruction methods, SVD (with regularization when indicated) and ART (with non-negative constraint, 1000 iterations, and 1, n n λ =∀) were tested under different noise levels (0 up to 10 −2 ). Figure 2 (Left) shows the maps of the reconstructed perturbations on the absorption ( a δ ) and diffusion ( δ D ) coefficients through the different reconstruction methods with increasing noise levels.On Fig. 2 (Right) are represented crossplots extracted from these maps along a horizontal line through the middle of the reconstruction area.

Fig. 2 .
Fig. 2. Comparison of different reconstruction methods at different noise levels: Left, results of the reconstructions of the perturbations on the absorption a δ and diffusion δ D coefficient maps; Right: cross-plots of the values extracted from a horizontal line in the middle of the vertical axis in the reconstruction areas.

Figure 3
Left Top shows the reconstructions obtained: with 2 sources the diffusion coefficient perturbation is almost not detected, but interestingly, the absorption perturbations are well recovered.The gain obtained in the reconstructions can be quantified through the quadratic error (QE, Fig.3Left Bottom) between the reconstructed values  

Fig. 3 .
Fig. 3. Influence of the number of sources.Left Top: reconstructed perturbation absorption a δ and diffusion δ D coefficient maps with (a) 2, (b) 4, (c) 8 and (d) 16 point sources illuminations.Left Bottom: Quadratic errors on the reconstructions of the absorption (black circles) and diffusion (blue squares) coefficients perturbations as a function of the number of point sources.Right Top: Normalized singular values (Left) and condition number of (H) as a function of the number of sources (Right).Right Bottom: cross-plots of the values extracted from a horizontal line in the middle of the vertical axis in the reconstruction areas.

Fig. 4 .
Fig. 4. Influence of the detectors distributions.Left Top: reconstructed perturbation absorption a δ and diffusion δ D coefficient maps.(a) the detectors are evenly distributed around the object (15 detectors on each side); (b) on each side, the detectors cover a length that is half the length of the object; (c) the length covered by the detectors is half of (b); (d) 30 detectors evenly distributed on each of the two vertical sides of the object; (e) 30 detectors evenly distributed on each of the two horizontal sides of the object.Left Bottom: Quadratic errors on the reconstructions of the absorption (black circles) and diffusion (blue squares) coefficients perturbations for the five different configurations.Right Top: Normalized singular values (Left) and condition number of (H) as a function of the number of sources (Right).Right Bottom: cross-plots of the values extracted from a horizontal line in the middle of the vertical axis in the reconstruction areas.

Fig. 5 .
Fig. 5. Influence of the sources shape.Left Top: reconstructed perturbation absorption a δ and diffusion δ D coefficient maps with (a) 4 points sources; 4 line sources with a length of (b) 0.2 mm; (c) 4 mm; (d) 8 mm; (e) 26 mm; (f) full field of view probing with one single source composed of lines.Left Bottom: Quadratic errors on the reconstructions of the absorption (black circles) and diffusion (blue squares) coefficients perturbations as a function of the length of the four sources ("Integral" represents situation (f), reconstructions obtained with a single measurement but with full field of view source illumination).Right Top: Normalized singular values (Left) and condition number of (H) as a function of the number of sources (Right).Right Bottom: cross-plots of the values extracted from a horizontal line in the middle of the vertical axis in the reconstruction areas.

Fig. 6 .
Fig. 6.Top Left: Schemas of the experimental situations.For each situation, only 4 measurements were performed at 0°, 90°, 180° and 270°, the rotation angles corresponding to rotations around the axis located at point O (cross marker), and perpendicular to the object plan.Schemas correspond to the measurement taken at 0°: (a) 4 points sources, 15 detectors; (b) 4 line sources with length 26 mm, 15 detectors; (c) same line sources, 1 detector; (d) same line sources, 2 detectors; (e) same line sources, 3 detectors; for each situations, four synthetic measurements were considered: these configurations (0°), and object rotated by 90°, 180° and 270°.The condition numbers of the corresponding (H) are reported for each situation.Top Right: reconstructed perturbation absorption a δ and diffusion δ D coefficient maps.Bottom: Quadratic errors (QE) on the reconstructions of the absorption (black circles) and diffusion (blue squares) coefficients perturbations for the five different configurations.