Depth dependence of estimated optical properties of a scattering inclusion by time-resolved contrast functions

In this study we have theoretically and experimentally investigated the behavior of first order approximation contrast function when purely scattering inhomogeneities located at different depths inside a turbid thick slab are considered. Results of model predictions have been compared with Finite element method simulations and tested on phantoms. To this aim, we have developed for the first time to our knowledge a fitting algorithm for estimating both the scattering perturbation parameter and the shift of the inhomogeneity from the middle plane, allowing one to reduce the uncertainties due to depth. This is important for optical mammography because effects of the depth can cause uncertainties in the derived tumor optical properties that are above 20% and the scattering properties of tumors differ from those of the sourrounding healthy tissue by a comparable extent. © 2008 Optical Society of America OCIS codes: (170.3660) light propagation in tissues; (170.5280) photon migration; (170.6510) spectroscopy, tissue diagnostics; (290.1990) diffusion; (290.7050) turbid media. References and links 1. B. Chance and R. R. Alfano, eds., Photon Migration and Imaging in Random Media and Tissues (1993). 2. R. R. Alfano and J. Fujimoto, “Advances in Optical Imaging and Photon Migration,” Optics & Photonics News 7, 37 (1996). 3. G. E. Cohn, W. S. Grundfest, D. A. Benaron, and T. Vo-Dinh, eds., Advanced Biomedical and Clinical Diagnostic Systems II (2004). 4. H. Xu, B. W. Pogue, H. Dehghani, K. D. Paulsen, R. Springett, and J. F. Dunn, “Absorption and scattering imaging of tissue with steady-state second-differential spectral-analysis tomography,” Opt. Lett. 29, 2043–2045 (2004). 5. S. Fantini, M. A. Franceschini, G. Gaida, E. Gratton, H. Jess, W. W. Mantulin, K. T. Moesta, P. M. Schlag, and M. Kaschke, “Frequency-domain optical mammography: Edge effect corrections,” Med. Phys. 23, 149–157 (1996). 6. G. Mitic, J. G. Koelzer, J. Otto, E. Plies, G. Soelkner, and W. Zinth, “Time-resolved transillumination of turbid media,” vol. 2082, pp. 26–32 (SPIE, 1994). 7. R. Cubeddu, C. D’Andrea, A. Pifferi, P. Taroni, A. Torricelli, and G. Valentini, “Effects of the menstrual cycle on the red and near-infrared optical properties of the human breast,” Photochem. Photobiol. 72, 383–391 (2000). 8. A. Cerussi, D. Jakubowski, N. Shah, F. Bevilacqua, R. Lanning, A. Berger, D. Hsiang, J. Butler, R. Holcombe, and B. tromberg, “Spectroscopy enhances the information content of optical mammography,” J. Biomed. Opt. 7, 60–71 (2002). #98536 $15.00 USD Received 10 Jul 2008; revised 6 Sep 2008; accepted 12 Sep 2008; published 17 Oct 2008 (C) 2008 OSA 27 October 2008 / Vol. 16, No. 22 / OPTICS EXPRESS 17667 9. B. W. Pogue, S. Jiang, H. Dehghani, C. Kogel, S. Soho, S. Srinivasan, X. Song, T. D. Tosteson, S. P. Poplack, and K. D. Paulsen, “Characterization of hemoglobin, water, and NIR scattering in breast tissue: analysis of intersubject variability and menstrual cycle changes,” J. Biomed. Opt. 9, 541–552 (2004). 10. L. Spinelli, A. Torricelli, A. Pifferi, P. Taroni, G. M. Danesini, and R. Cubeddu, “Bulk optical properties and tissue components in the female breast from multiwavelength time-resolved optical mammography,” J. Biomed. Opt. 9, 1137–1142 (2004). 11. U. Hampel, E. Schleicher, and R. Freyer, “Volume image reconstruction for diffuse optical tomography,” Appl. Opt. 41, 3816–3826 (2002). 12. Y. Xu, X. Gu, T. Khan, and H. Jiang, “Absorption and Scattering Images of Heterogeneous Scattering Media Can Be Simultaneously Reconstructed by Use of dc Data,” Appl. Opt. 41, 5427–5437 (2002). 13. T. Yates, J. C. Hebden, A. Gibson, N. Everdell, S. R. Arridge, and M. Douek, “Optical tomography of the breast using a multi-channel time-resolved imager,” Phys. Med. Biol. 50, 2503–2517 (2005). 14. J. Ye, K. Webb, R. Millane, and T. Downar, “Modified distorted Born iterative method with an approximate Frchet derivative for optical diffusion tomography,” J. Opt. Soc. Am. A 16, 1814 (1999). 15. F. Gao, Y. Tanikawa, H. Zhao, and Y. Yamada, “Semi-Three-Dimensional Algorithm for Time-Resolved Diffuse Optical Tomography by Use of the Generalized Pulse Spectrum Technique,” Appl. Opt. . 41, 7346 (2002). 16. D. Grosenick, H. Wabnitz, and H. Rinneberg, “Time-resolved imaging of solid phantoms for optical mammography,” Appl. Opt. 36, 221–231 (1997). 17. B. Wassermann, A. Kummrow, K. T. Moesta, D. Grosenick, J. Mucke, H. Wabnitz, M. Moller, R. Macdonald, P. M. Schlag, and H. Rinneberg, “In-vivo tissue optical properties derived by linear perturbation theory for edgecorrected time-domain mammograms,” Opt. Express 13, 8571–8583 (2005). 18. A. Pifferi, P. Taroni, A. Torricelli, F. Messina, R. Cubeddu, and G. Danesini, “Four-wavelength time-resolved optical mammography in the 680–980-nm range,” Opt. Lett. 28, 1138–1140 (2003). 19. J. C. Hebden and S. R. Arridge, “Imaging Through Scattering Media by the Use of an Analytical Model of Perturbation Amplitudes in the Time Domain,” Appl. Opt. 35, 6788 (1996). 20. M. Morin, S. Verrealut, A. Mailloux, J. Fréchette, S. Chatigny, Y. Painchaud, and P. Beaudry, “Inclusion characterization in a scattering slab with time-resolved transmittance measurements: perturbation analysis,” Appl. Opt. 39, 2840–2852 (2000). 21. S. Carraresi, T. S. M. Shatir, F. Martelli, and G. Zaccanti, “Accuracy of a Perturbation Model to Predict the Effect of Scattering and Absorbing Inhomogeneities on Photon Migration,” Appl. Opt. 40, 4622 (2001). 22. L. Spinelli, A. Torricelli, A. Pifferi, P. Taroni, and R. Cubeddu, “Experimental test of a perturbation model for time-resolved imaging in diffusive media,” Appl. Opt. 42, 3145–3153 (2003). 23. H. Rinneberg, D. Grosenick, K. Moesta, J. Mucke, B. Gebauer, C. Stroszczynski, H. Wabnitz, M. Moeller, B. Wassermann, and P. Schlag, “Scanning time-domain optical mammography: detection and characterization of breast tumors in vivo,” Technol Cancer Res Treat 4, 483–496 (2005). 24. R. Esposito, S. De Nicola, M. Lepore, I. Delfino, and P. L. Indovina, “A perturbation approach to characterize absorptive inclusions in diffusing media by time-resolved contrast measurements,” J. Opt. A 6, 736–741 (2004). 25. R. Esposito, S. De Nicola, M. Lepore, and P. L. Indovina, “Perturbation approach to the time-resolved transmittance for a spatially varying scattering inclusion in a diffusive slab,” J. Opt. Soc. Am. A 23, 1937–1945 (2006). 26. R. C. Haskell, L. O. Svaasand, T. T. Tsay, C. Feng, M. S. McAdams, and B. J. Tromberg, “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A 11, 2727–2741 (1994). 27. S. De Nicola, R. Esposito, M. Lepore, and P. L. Indovina, “Time-resolved contrast function and optical characterization of spatially varying absorptive inclusions at different depths in diffusing media,” Phys. Rev. E 69, 031901 (2004).


Introduction
Imaging with near-infrared radiation (NIR) is currently investigated for medical diagnosis applications [1,2,3]. The complexity of NIR propagation through biological tissues is the main difficulty to realize a noninvasive diagnostic tool for detecting tumors and other abnormalities hidden in thick biological tissues. Particularly, the strong scattering of light by biological tissues dramatically deteriorates image quality and produces images with poor structural details. Although contrast and spatial resolution of optical imaging in turbid media can be improved by time-resolved or frequency-domain techniques, the performance of optical imaging systems are still lower than the x-rays ones [4,5,6]. Nevertheless, the optical images contain exquisite functional information that complements those obtained by other imaging modalities. Particularly, the measurements of diffusely reflected or transmitted NIR light allow the detection of local changes in scattering and absorption properties of tissue, that can be related to composi-tion, vascularization, blood oxygen saturation, and structural properties of normal breast tissue, benign lesions and tumors as well [7,8,9,10].
For mammography, NIR light is used to screen for breast cancer within diffusion approximation theory. Images are generally obtained by means of two approaches: optical tomography and single two-dimensional image projection. The former recovers the spatial distribution of the optical parameters of the breast by using multiple input and output locations on the surface of the sample [11,12,13]. Results are related to the physiological parameters of tissue and are presented as three-dimensional images to clinicians. However, the algorithms used to generate images from measurements are generally expensive due to the large amount of data to handle [14,15]. Also, the measurement set-up is complex with the need to exploit different source-detector combinations.
Alternatively, breast images can be obtained by the simpler approach of two-dimensional image projection. In this case, a source-detector pair synchronously scans the surface of the organ that is slightly compressed between two parallel planes and the infinitely extended slab is usually adopted to model breast [16,17,18]. This approach has the advantage to produce projection images of the compressed breast that can be straightforwardly compared to the standard x-ray mammograms. Since analytical expressions that describe heterogeneous structures cannot be obtained in a closed-form, perturbation techniques are considered. Approximated expressions have been obtained for the time-resolved transmittance and reflectance through a turbid slab and are currently employed for a clinical study on 2-D optical mammography [19,20,21,22]. In this case the optical breast images are constructed by using a priori knowledge of tumor geometrical parameters because the measurements recorded with on-axis geometry do not allow one to separately derive the tumor optical coefficients together with the tumor size and its location along the compression direction. Tumor size is generally inferred from histopathological findings or x-ray mammograms, MR mammograms or ultrasound examinations. On the contrary, the location of a tumor along the compression direction is estimated from transmittance recorded at various projection angles when such measurements are available. Preliminary results of clinical investigations have shown that tumor absorption always exceeds absorption by healthy tissue, on average by about 150%. Tumours may scatter light less or more strongly than the surrounding healthy breast tissue, yet in the majority of cases scattering by tumours is stronger by about 25% [23].
When a priori knowledge of the geometry is not available, the estimates of the optical properties of tumors are affected by uncertainties in tumor shape and location. Inspections on phantoms have evidenced that the influence of shape can be partly compensated by averaging the optical properties of the same tumor derived for mediolateral and craniocaudal projections. On the contrary, the influence of the depth on contrast image is not yet understood at present, although clinical studies have shown that the uncertainty in the depth may causes uncertainties in tumor optical properties derived that are above 20%. Since the scattering differences between anomalous region and sourrounding healthy tissue are comparable with that uncertainties, the dependence of those estimates on the depth must be reduced in order to generate useful scattering images.
For these reasons we have theoretically and experimentally investigated the behavior of the reduced scattering coefficient derived for a defect embedded at different depths in a uniform background by considering a perturbation approach to diffusion equation. On one hand the paper complements the analysis of the perturbation models for 2D imaging reconstruction reported in Refs. [24,25], and on the other it extends it because, the first time to our knowledge we present a fitting algorithm for estimating both the scattering perturbation parameter and the shift of the inhomogeneity from the middle plane in order to reduce uncertainties due to the depth of the defect. Finally, all these results have been experimentally validated by performing  Fig. 1. Geometric scheme assumed for a Gaussian scattering inclusion of cylindrical shape geometry with radius R and height h located at different depths z = z pc inside a turbid slab of thickness d. A pulsed light beam illuminates the front surface of the scattering slab at plane z = 0. The photons are assumed to be initially isotropically scattered at a depth z s = 1/μ s below the front surface. The time-resolved transmittance is measured by a detector at plane z = d coaxial with the source and with the inclusion. measurements on phantoms containing scattering inclusion located at different depths.
In Sec. 2 we will summarize and discuss the basic results of the first order perturbation approach to determine the analytical expression for the time-resolved transmittance through a turbid slab with a Gaussian scattering inclusion in the coaxial detection scheme. In Sec. 3 we analyze the dependence of the perturbation model contrast function on the depth of the inclusion. Particularly, finite-element method (Fem) simulations will be used to validate the perturbation approach for a case of practical interest, that is representative of a slightly compressed breast. In Sec. 4 we describe the novel fitting algorithm that will be used to investigate the accuracy of the perturbation model as different depths and sizes of the inclusion are considered. Results are compared with those obtained with the algorithm that assumes the inclusion fixed at the center of the slab. Furthermore, the accuracy for retrieving the shift of the inclusion from the middle plane of the slab will be discussed. The experimental set up and the phantoms optical characteristics are described in Sec. 5.1, whereas experimental results are discussed in Sec. 5.2. Figure 1 shows a scattering inhomogeneity located at depth z pc inside an infinitely extended homogeneous slab with thickness d. This anomalous region is shaped like a cylinder with a top hat profile that extends from z pc − h/2 to z pc + h/2 along the imaging direction. The scattering changes along the radial direction are assumed to be described by the following Gaussian profile:

Theoretical model
where Δμ s is the maximum deviation, observed on the cylinder axis, of the reduced scattering coefficient of the inclusion from the unperturbed value μ s of the host medium. The radial size of the inhomogeneity is determined by the radius R and the factor f . For f = ln 2 the quantity R is the distance at which the perturbation intensity decreases to Δμ s /2; this is a value that is halved by doubling f . Let us consider a pencil pulsed beam that hits normally the surface z = 0 of the slab at time t s = 0 and at position (x s = 0, y s = 0) (see Fig. 1); this is a configuration that can be always realized through an appropriate choice of the system coordinate. The light that crosses the surface propagates in the medium undergoing many scattering events and, according to diffusion approximation, photon migration takes place as if photons were initially isotropically scattered at a depth z s = 1/μ s below the front surface. Particularly, the density photon fluence rate ϕ(r,t; r s ) at position r and time t is described by the following parabolic differential equation: where r s ≡ (x s = 0, y s = 0, z s = 1/μ s ), v is the speed of light in the medium, μ a is the absorption coefficient of the slab and D = 1/3(μ s + δ μ s ) is the diffusion coefficient that takes into account for the presence of the scattering inclusion at depth z pc . Within the extrapolated boundary conditions, the fluence rate has to vanish at extrapolated flat surfaces outside the turbid slab at a distance z e = 2AD 0 from the physical boundaries, with D 0 = 1/3μ s and A = (1 + R e f f )/(1 − R e f f ) [26]. Consequently, the physical quantity of interest, that is the time-resolved transmittance T (t, z pc ), can be written in the following way: where the position r m is a measurement point on the back z = d surface of the slab. Henceforth, we will consider a coaxial configuration for the source-detector system, that is we will always refer to a measurement point r m ≡ (0, 0, d).
Under the condition |Δμ s |/μ s 1, a first order perturbation approach to diffusion equation (2) can be applied. It results that small changes in scattering properties affect the fluence rate linearly and the time-resolved transmittance can be approximated by a sum of two terms, namely: The quantity T 0 (t) describes the temporal behavior of the time-resolved transmittance through the homogeneous slab and it is given by: where d e = d + 2z e is the extrapolated thickness of the slab. The second term δ T μ s t; z pc on the right-hand side of Eq. (4) is the change, at first-order approximation, in the time-resolved transmittance due to the presence of the scattering inclusion at depth z pc inside the homogeneous turbid slab. When the source, the detector and the Gaussian scattering inclusion, that is modeled by (1), are collinear the scattering time-resolved perturbation δ T μ s t; z pc can be written as follows: The terms R ± k,l (R,t) in Eq. (6) are functions of the width R of the inclusion and the time t, namely: The functions Z ± k,l (z pc , h) in Eq. (6) depend on the depth of the inclusion z pc and on its axial extent h, and can be written in the following form: where the quantities γ ± are given by The analytical expression (6) for the scattering perturbation is slightly different from that reported in Ref. [25], although both formulas are obtained from the same perturbation integral (Eq. (8) of Ref. [25]). Differences arise from a few manipulations that have been performed on this integral in order to make the perturbation formula computationally more efficient. Particularly, the time t 0 = min[t b1 ,t b2 ] in Eqs. (7, 8b) is deduced by considering that no light can be present in the region of the inclusion before the time of flight source-lower side of the defect, Moreover, no light can be received by the detector from the inclusion before the transit time t b2 = (d − z pc − h/2)/v. The quantity ε in Eqs. (8a, 10) has been introduced to avoid the discontinuity of the terms of the series for k = l. Its value has been set to 10 −6 for computational purposes.

Depth dependence of contrast functions
Equation (6) gives the change in the time-resolved transmittance due to a Gaussian scattering inclusion of radius R and height h that is hidden inside a turbid slab at depth z pc from the front surface. The analytical expression has been obtained in the framework of the first order perturbation approach to the diffusion equation. It is expected to give accurate results only for local changes of the reduced scattering coefficient that are negligible with respect to the value of the reduced scattering coefficient of the surrounding medium.
According to the procedure adopted in Ref. [27], the accuracy of the perturbation model is now investigated by solving numerically the diffusion equation (2) under the extrapolated boundary conditions for a turbid slab with Fem. Particularly, the numerical time-resolved transmittance T num t, z pc is obtained from the solution ϕ num (r,t; r s ,t s = 0) of the photon fluence rate by using the relation (3). Results refer to an homogeneous slab of thickness d = 40 mm with absorption coefficient μ a = 0.01 mm −1 and reduced scattering coefficient μ s = 1.0 mm −1 . The mismatch of the refractive index at medium-air interface has been set to the value n = 1.4. These values are of practical interest because they are representative of breast tissues. Moreover, the Gaussian scattering inhomogeneity has been shaped with height h equals to the diameter 2R and has been shifted along the probe beam-detector z axis from the front z = 0 to the rear z = d surface of the slab. The relative perturbation parameter Δμ s /μ s has been ranged from −60% to 60% for inclusions that extend in size from R = 2.5 mm to R = 10 mm. The factor f in Eq. 1 has been set to f = 2 ln2, thus the perturbation intensity δ μ s (ρ) reduces to a quarter of its maximum value Δμ s at the radial distance ρ = R.
The first-order-perturbation time-resolved transmittance T (t, z pc ), that is given by Eqs. (4,5,6), is compared with Fem simulations by introducing the following contrast functions The contrast C(t, z pc , Δμ s ) (11a) is the change of transmittance with respect to the homogeneous case when the perturbation model is considered, whereas the Eq. (11b) is the relative change of the transmitted signal as predicted by numerical simulations. For facilitating the comparison between different inclusions a normalized depth z pcn has been defined, as follows: where z min = z s + h/2 is the minimum distance of the center z pc of the inclusion from the front surface and z max = d − z s − h/2 is the maximum depth of the inclusion, that is the minimum distance from the rear surface. In such a way the normalized depth ranges from 0 to 1 whichever the size of the inclusion. The temporal behavior of contrast functions (11a) and (11b) is shown in Fig. 2(a) for a Gaussian inclusion with radius R = 2.5 mm. It results from comparisons with Fem simulations (solid curves) that the perturbation model (dashed curves) underestimates the contrast function. The discrepancies increase as the time-of-flight of transmitted photons tends to the ballistic time d/v and as the inclusion is moved away from the middle plane. Particularly, the differences are less than 16% for a relative perturbation parameter ranging from Δμ s /μ s = −20% to Δμ s /μ s = 20%. It can be also seen that the module of contrast functions increases as the centre of the inclusion is displaced from the central plane, z pcn = 1/2, to a position near to the rear boundary, z pcn = 1. Thus inclusions located close to the medium-air interface affect the migration of detected photons more significantly, which enhances the difference between the perturbed transmitted signal and the unperturbed one.
In order to better investigate the dependence of the time-resolved contrast on the depth z pc , the ratio C num (t, z pc , Δμ s )/C num (t, 1/2, Δμ s ) must be considered because it allows for monitoring the change of the contrast when the inclusion moves from the central plane toward the boundary surface. In Fig. 2(b), it is considered the case of an inclusion with radius R = 2.5 mm and relative perturbation intensity Δμ s /μ s = 0.2 whose depth is changed from the normalized value z pcn = 0.5 to z pcn = 1 with step equal to 0.1. Accordingly to results of Fig. 2(a) the relative change of the contrast increases more and more as the inclusion approaches to the boundary. The curves reported in a semi-logarithmic scale evidence that there is not a merely proportionality between the contrast assessed at depth z pcn and that observed at the middle plane z pcn = 1/2. This behavior is also observed for inclusions with different sizes and perturbation intensities as it follows by further investigations that are not reported here for brevity. It is worth noting that the shift along the z axis of the perturbation significantly alter the shape of the contrast function causing an offset and a deviation at early times ( Fig. 2(a)). Thus, the fact that upon z pcn , the contrast is progressively modified in time (see Fig. 2(b)) is an important pre-requisite for the capability to extract information on z pc out of experimental data.

Simulations
From results of the previous section one can demonstrate that each value of the contrast is observed at two different depths symmetrically placed with respect to the central plane of the scattering slab. Then, the depth z pc of the inclusion cannot be retrieved unambiguously through comparisons of the perturbation model contrast with Fem simulations (or experimental data) at least in the framework of a coaxial source-detector configuration. However, results of Fig. 2(b) suggest that the contrast function analysis can be used for estimating the shift Δz pc = z pc − d/2 of the inclusion from the middle plane as well as its scattering perturbation parameter Δμ s . For achieving this aim we have developed a procedure of comparison between contrast functions that is described in the following.
The approximated time-resolved contrast (11a) is calculated by considering a scattering inclusion hidden inside an homogeneous slab with thickness d = 40 mm and refractive index mismatch n = 1.4. The optical parameters of the slab used for calculating the contrast are obtained by fitting the T 0 (t) function (5) to the Fem simulation of the transmittance through the homogeneous slab. The inclusions have radius R = 2.5, 5, 7.5, 10 mm and height h = 2R. For each value of the radius, a set of time-resolved contrasts is calculated by locating the cylindrical inclusion at different depths along the source-detector axis. The array of the calculated contrast functions for each size R is numerically interpolated to obtain the interpolating function C(t, z pc ) | interp of time t and depth z pc . This function is implemented in a fitting procedure (Proc(Δz pc Δμ s )) that determines the best estimates of the depth, z pc, f it , and the perturbation intensity, Δμ s, f it , by minimizing the chi-squared function defined as: The range of integration [t min ,t max ] in Eq. (13) identifies the points of the homogeneous transmittance with intensity higher than 30% of the peak on the leading edge of the curve and 1% on the trailing edge. The estimation Δz pc, f it for the shift of the inclusion with respect to the middle So far, the fitting procedure for retrieving the perturbation intensity has been performed by assuming that the position of the inclusion was fixed at midway between the source and the detector (Proc(Δμ s )), whichever is the depth [22]. In such a way, the fitting procedure generate an estimate for the scattering perturbation parameter that is affected by that constraint on the depth.
It is evident that the difference between the estimates of the perturbation intensity given by the two procedures Proc(Δz pc Δμ s ) and Proc(Δμ s ) enhances with increasing the shift of the inclusion from the centre of the slab. For investigating how the assessment of the perturbation intensity is affected by the shift of the defect, the relative difference between the fitted value Δμ s, f it and the expected one Δμ s , namely: has been calculated for each fitting procedure. Figure 4 shows the relative error (14) as a function of the normalized depth z pcn for inclusions with a relative perturbation parameter Δμ s /μ s ranging between −60% and 60%. Also in this case the use of the normalised depth z pcn facilitates the comparison between inclusions with different sizes. Particularly, black curves of Fig. 4(a) refer to an inclusion with radius R = 2.5 mm that is moved from z pcn = 0 to z pcn = 1, whereas the green curves of Fig. 4(b) are the results for R = 10 mm. The solid curves describe the relative error ε Δμ s calculated with the fitting procedure Proc(Δz pc Δμ s ) that estimates both the shift and the perturbation intensity of the inclusion. The dashed curves are the relative error ε Δμ s calculated by considering the procedure Proc(Δμ s ). From inspection of the curves it results that the relative error calculated with the procedure Proc(Δz pc Δμ s ) does not appreciably change as the axial position of the inhomogeneity is moved from the front to the rear boundary of the slab for fixed values of the radius R of the inclusion. As expected, the accuracy worsens with increasing the radius R and the worst case is observed for R = 10 mm with | ε Δμ s | 100% and Δμ s /μ s = −0.6. On the contrary, the minimum value of the relative error is | ε Δμ s | 20% and is attained for R = 2.5 mm and Δμ s /μ s = 0.6. As far as it concerns results given by Proc(Δμ s ), it can be noticed that the discrepancies with Proc(Δz pc Δμ s ) curves enhance as the inclusion moves toward the lateral surfaces of the slab and decrease as the size R of the inclusion increases. Particularly, the maximum values of the relative error are observed at normalized depths z pcn = 0, 1 with | ε Δμ s | 600% for R = 2.5 mm and | ε Δμ s | 150% for R = 10 mm with Δμ s /μ s = −0.6. This behaviour can be explained by considering that Proc(Δμ s ) assumes the location of the inclusion is fixed at middle plane of the slab. Hence, the accuracy for retrieving the perturbation intensity worsens when the shift of the inclusion increases. Obviously, this effect reduces as the size of the inclusion increases and vanishes in the limit of a size equal to the thickness of the slab. On the whole, the worst accuracy is obtained when negative inclusions are considered. This behaviour arises from the fact that the discrepancies between Fem simulations and perturbation model contrast functions observed for negative values of the perturbation intensity are more marked than those corresponding to positive values.
As far as it concerns the accuracy in retrieving the shift Δz pc with Proc(Δz pc Δμ s ), it has been investigated by considering the following quantity: Specifically, the deviation from the expected value Δz pc is calculated at different depths for each inclusion radius and each perturbation intensity. The average value of each set is calculated and then it is normalized with respect to the diameter 2R of the inclusion. Figure 5 shows the normalized average error ε Δz pc and its standard deviation as a function of the relative perturbation parameter Δμ s /μ s for four values of the radius R ranging between R = 2.5 mm to R = 10 mm. As it can be seen, the normalized error ε Δz pc is less than ∼ 5% as the relative perturbation intensity Δμ s /μ s ranges from −60% to 60% in the considered range of values for the radius R. Hence, the shift of the inclusions, that have diameter 2R between 5 mm and 20 mm, is retrieved by means of Proc(Δz pc Δμ s ) with a maximum error that goes from 0.2 mm to 1.0 mm, respectively.

Phantoms and Experimental Setup
The measurements were carried out by using a d = 40 mm thick cell made of black polyvinyl chloride with two small transparent apertures for reducing the inner boundary reflections. The media were prepared using an aqueous solution of Intralipid and black ink with optical properties μ a = (8.7 ± 0.2) × 10 −3 mm −1 and μ s = 1.17 ± 0.01 mm −1 . As for the inclusions, they were made of 2% agar, Intralipid, and black ink in distilled water. Two cylindrical samples were prepared with radius R = 3 mm and R = 7.5 mm, respectively, whereas the thickness h were equal to the diameter 2R. The reduced scattering coefficient of the inclusion were the only optical parameter changed with respect to the background and corresponded to a perturbation intensity Δμ s = 0.46 ± 0.4 mm −1 .
The experimental data were recorded with a setup for time-resolved transmittance based on a picosecond solid state laser at 785 nm (mod. PDL 800, PicoQuant GmbH, Germany) and an electronic chain for time-correlated single photon counting (mod. SPC-300, Becker & Hickl GmbH, Germany). The instrument response function of the setup is about 140 ps. Point measurements were performed by moving the inclusion inside the host medium along the x direction from x = −40 mm to x = 40 mm every 1 mm, transversally to the source-detector axis that is fixed at position x = 0.0 ± 0.5 mm (see Fig. 1). The acquisition time for each point has been chosen so as to obtain an homogeneous transmittance with peak value of ∼ 3× 10 5 counts. In this way the 1% limit considered to calculate the chi-squared function of Eq. 13 identifies points at the trailing edge whose intensity has a 2% relative error about and consequently are not affected by fluctuations of the baseline level.

Results and discussions
Figure 6(a) shows the experimental time-resolved contrast meas k curve refers to the inclusion located at center of the slab, that is z pc = 20 mm, whereas the red curve is obtained when the depth z pc = 28 mm is considered. Similarly, Fig. 6(b) reports plots obtained by considering the inclusion with radius R = 7.5 mm. In agreement with results discussed in section 3, the experimental contrast enhances by decreasing the time of flight and by shifting the defect away from the center of the slab. Particularly, the displacement of the inclusion causes a rise in contrast of roughly 30% for R = 3 mm and of 50% when R is equal to 7.5 mm.
For all the phantoms, the time-resolved contrast function has been calculated for each point of the linear scan in order to reconstruct the profile of the inclusion. To this purpose, the perturbed contrast function (11a) has been fitted to the experimental data by using both procedures Proc(Δμ s ) and Proc(Δz pc Δμ s ). The former admits only the perturbation intensity as fitting parameter, whereas the depth is fixed to the value of the center of the cell. The latter, instead, considers both the perturbation intensity and the depth of the inclusion as fitting parameters. To eliminate the system response, the fitting has been performed by adopting the model (11a) with expressions (5) and (6) convolved with system response. According to these procedures, we expect the fitted value Δμ s, f it (x) obtained at the scan position x corresponds to the scattering perturbation Δμ s of the inclusion only when it is actually located on the source-detector axis. Otherwise, the fitting detects an effective inclusion located on the axis with a perturbation intensity value between zero and Δμ s . For these reasons, a maximum should be observed in the map of the optical perturbation at position of the center of the inclusion. The ratio between this maximum and the background intensity, Δμ s, f it (0)/μ s , is called contrast image and it is usually estimated by fitting the image of the perturbation with a Gaussian curve [22].   In Fig. 7 we have plotted the estimated values of the relative perturbation intensity Δμ s, f it (x)/μ s versus the position x for the inclusions with radius R = 3 mm and R = 7.5 mm.
The procedure Proc(Δμ s ) produces the plots shown in 7(a) and (c), whereas plots 7(b) and (d) show results of the procedure Proc(Δz pc Δμ s ). Furthermore, in each panel black points refer to the inclusion at depth z pc = 20 mm, red points to the case z pc = 28 mm, whereas solid curves represent the Gaussian fits. Moreover, in table 1 they are reported the amplitude values of the fitted Gaussians in the different cases.
As it can be seen from the first row of the table, the procedure Proc(Δμ s ) images the inclusions with a contrast that depends on the depth. Particularly, the shift of the smallest inclusion from the center of the cell to the depth z pc = 28 mm causes an enhancement in contrast image of 41%, which gives rise to two different images as it is shown in panels (a) and (c) of Fig. 7. Data reported in row 1, column 3 and 4 of table 1 evidence that the effect of contrast enhancement is also observed for the inclusion with R = 7.5 mm, although the contrast difference reduces to 13% in this case, in agreement with results of Fig. 4 discussed in section 4. Hence experimental results confirm investigations on Fem simulations that have evidenced a marked sensitivity of perturbation intensity to the depth when the fitting procedure Proc(Δμ s ) is used. However, theoretical investigations have also shown that this behavior can be reduced by adopting the procedure Proc(Δz pc Δμ s ). To this regard, panels (b) and (d) of Fig. 7 show the images of the two inclusions produced by this procedure. As it can be seen from an inspection of the curves and from data reported in row 2 of table 1, the contrast image calculated at different depths has a value that is very close to that observed in the central plane. Particularly, the difference in contrast reduces to 10% for the smallest inclusion and to 4% for the largest one. These values are in good agreement with predictions than are deduced from theoretical analysis on the accuracy of the first order perturbation model. In fact, it follows from results reported in Fig. 4 that 8 mm shift of the inclusion with R = 3 mm causes a 13% enhancement in contrast that decreases to 3% when the size of the defect increases to R = 7.5 mm.

Conclusions
We have investigated the behavior of time-resolved contrast functions for a scattering inclusion when its depth inside a turbid slab is shifted from the central plane to a boundary surface. The theoretical model has been derived by considering a Gaussian first order perturbation model to the diffusion equation. The accuracy of the predictions of this model has been studied through comparisons with Fem simulations by developing for the first time to our knowledge a fitting procedure that retrieves both the perturbation intensity Δμ s and the shift Δz pc of the inclusion from the center of the slab (Proc(Δz pc Δμ s )). The effectiveness of the proposed procedure has been tested against the usual fitting algorithm that allows only the perturbation intensity to vary, whereas the depth of the inclusion is fixed at the center of the slab (Proc(Δμ s )).
The Proc(Δz pc Δμ s ) provides an accuracy in the estimation of the scattering perturbation intensity that does not change appreciably as the inhomogeneity moves from the front to the rear boundary. Particularly, the relative error in retrieving the scattering perturbation parameter is less than 20% for R ≤ 2.5 mm and |Δμ s /μ s | ≤ 60%. Moreover, the shift of the inclusion is retrieved with an error less than 8% of the value of the radius R of the inclusion. These figures must be compared with a maximum relative error of | ε Δμ s | 600% for R = 2.5 mm obtained from an inclusion near the boundary surfaces with the procedure Proc(Δμ s ).
We have also test experimentally our perturbation analysis by considering two scattering inclusions with different radius (R = 3 and 7.5 mm), the same perturbation intensity (Δμ s = 0.46 mm −1 ), located at different positions inside a tank filled with a scattering solution: one is at the center and the other is at about a quarter of the tank thickness. To this regard, a Fem based simulation investigation was preliminarily used to assess the amount of discrepancies in retrieving the perturbation intensity, when using a theoretical model that refers to a spatially varying scattering coefficient inside the defect compared to the case of homogeneous objects employed in the experiments. Discrepancies have been found less than 2%, thus validating the data analysis.
To conclude, we have observed that the time-resolved contrast functions for an inclusion out of the center of the slab are not merely proportional to those for the same inclusion located at the center. This fact allows one to recover both the perturbation intensity and the inclusion displacement from the center by exploiting transmittance on-axis measurements, that are typically done in time-resolved 2-D projective optical mammography. In this way we increase the accuracy of the retrieved scattering properties of inhomogeneities hidden inside a diffusing medium at an unknown depth. Then, improvements in physiological information on the tumor and the surrounding tissue will be gained, thus increasing the diagnostic value of optical mammography. However, a closer characterization of the tumor could be achieved by extending our procedure to the case of objects that differ not only in scattering but also in absorption from the surrounding tissue. For this reason, we aim in the future to develop such an algorithm, by paying particular attention to the technique of data analysis. Indeed, it is not clear whether it is more efficient to retrieve the optical parameters simultaneously or through a separation of the scattering and the absorptive contribution.