1. 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 composition, 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 measurements on phantoms containing scattering inclusion located at different depths.

 

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.

Download Full Size | PDF

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.

2. Theoretical model

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:

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=ln2 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 ₀ from the physical boundaries, with D ₀=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 ₀ (t) describes the temporal behavior of the time-resolved transmittance through the homogeneous slab and it is given by:

where d_e=d+2 z_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:

where b=R ²/(fDvt). The terms 𝓡 ^± _k,l(R,t) in Eq. (6) are functions of the width R of the inclusion and the time t, namely:

with

The functions $Z_{k,l}^{\pm }(z_{\mathrm{pc}},h)$ ^± _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 ₀=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, t _b1=(z_pc-h/2-z_s)/v. 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.

3. 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 mmto R=10 mm. The factor f in Eq. 1 has been set to f=2ln2, thus the perturbation intensity δµ′_s(ρ) reduces to a quarter of its maximum value Δµ′_s at the radial distance ρ=R.

 

Fig. 2. (a) The temporal behavior of contrast functions C(t, z_pc,Δµ′_s) and C_num(t, z_pc,Δµ′_s) for an inclusion with radius R=2.5 mm at normalized depths z_pcn=1/2,3/4, 1. The solid curves refer to Fem simulations, whereas the dashed curves are the perturbation model predictions. The relative scattering perturbation parameter Δµ′_s/µ′_s ranges between -20% and +20%. (b) The ratio C_num(t, z_pc)/C_num(t,1/2) of the numerical contrast for different depths of an inclusion with radius R=2.5 mm and perturbation intensity Δµ′_s/µ′_s=0.2.

Download Full Size | PDF

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/ν 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.

4. 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 ₀(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 plane is given by using the relationship Δz _{pc, f it}=|z_{pc, f it} -d/2|. The uniqueness of the derived values of Δµ′_{s, f it} and Δz _{pc, f it} is assured by the uniqueness of the minimum of the chi-squared function (13) in its domain as it is shown, for example, by contour plots of Fig.3 that consider the case of an inclusion with radius R=2.5 mm and perturbation intensity Δµ′_s/µ′_s=0.2. The panel (a) represents the behaviour of the chi-squared function when the inclusion is at the center of the slab (Δz_pc=0 mm), whereas the panel (b) is obtained by shifting the inclusion from the center to Δz_pc=10 mm.

 

Fig. 3. Contour plots of the chi-squared function (13) versus the Δz_pc-Δµ′_s plane for an inclusion with size R=2.5 mm and perturbation intensity Δµ′_s/µ′_s=0.2. (a) the inclusion is at the center of the slab (Δz_pc=0 mm). (b) the inclusion is shifted from the center to Δz_pc=10 mm.

Download Full Size | PDF

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.

 

Fig. 4. The relative error as a function of the normalized depth z_pcn. The solid curves are obtained with the fitting procedure Proc(Δz_pcΔµ′_s), whereas the dashed curves are given by the procedure Proc(Δµ′_s). (a) The inclusion has size R=2.5 mm. (b) The inclusion has size R=10 mm.

Download Full Size | PDF

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.

 

Fig. 5. The error for retrieving the shift Δz_pc of the inclusion normalized with respect to the diameter 2R as a function of the relative perturbation intensity Δµ′_s/µ′_s. Four values of the radius of the Gaussian inclusion have been considered: R=2.5 mm (black curves), R=5 mm (red curves), R=7.5 mm (blue curves), R=10 mm (green curves).

Download Full Size | PDF

5. Experimental measurements

5.1. 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.4mm^-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=-40mm 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 ⁵ 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.

 

Fig. 6. Contrast functions calculated at different depths from experimental time-resolved transmittances concerning two cylindrical scattering inclusions with perturbation parameter Δµ′_s/µ′_s=39%. The background is a turbid slab with thickness d=40 mm, absorption and reduced scattering coefficient equal to µ_a=(8.7±0.2)×10^-3 mm^-1 and µ′_s=1.17±0.01 mm^-1, respectively. (a) The inclusion with radius R=3 mm is located at depths z_pc=20 mm (black curve) and z_pc=28 mm (red curve). (b) The inclusion with radius R=7.5 mm is located at depths z_pc=20 mm (black curve) and z_pc=28 mm (red curve)

Download Full Size | PDF

5.2. 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].

 

Fig. 7. The relative perturbation intensity (Δµ′_s/µ′_s) f it as a function of the scan position for the inclusions. Black points refer to the inclusion at depth z_pc=20 mm whereas red ones to the case zpc=28 mm. Solid curves represent the Gaussian fits. Plots (a) and (c) have been generated by using the fitting procedure Proc(Δµ′_s), wheres plots (b) and (d) have been constructed with the fitting procedure Proc(Δz_pcΔµ′_s).

Download Full Size | PDF

Table 1. The relative perturbation intensity Δµ′_s/µ′_s estimated by means of the Gaussian fits reported in Fig. 7 by using both the procedures Proc(Δµ′_s) and Proc(Δz_pcΔµ′_s).

View Table

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=28mm 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.

6. 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.