1. INTRODUCTION

The development of new optical techniques for biomedical purposes is an important challenge in recent years because of its potential to lead to non-invasive, low cost clinical tools that probe tissue function. Imaging of hemoglobin concentration, hemoglobin oxygen saturation and exogenous contrast agent uptake are the main target applications of these techniques. The usefulness of optical contrast has already been assessed in breast tumors detection^1–3 and brain function monitoring.⁴

The optical techniques are sensitive to tissue absorption heterogeneities due to spatial variations of blood saturation and blood volume. The accurate detection, localization and characterization of absorption heterogeneities are therefore an important issue for the novel bio-optical methods. Changes in micro-organelle concentration may also yield optical contrast due to local changes in tissue scattering. Diffuse Optical Tomography (DOT) has been developed to locally resolve tissue absorption and scattering. DOT is based on a rigorous forward model that predicts light propagation in diffuse media. By fitting or inverting this model for a set of measurements, one is able to recover absorption and/or scattering tissue contrast. A simple but efficient way to construct the forward model employs a perturbation-based solution of the diffusion equation using the Rytov approximation. This theoretical approach leads to an analytical expression for the forward model.⁵ Such analytical expressions have already given valuable results in the clinical setting.⁶ The key of the success of the Rytov approach is that it has been shown to be more accurate than the Born-type perturbative solution^7,8 and allows performing differential DOT.⁹

An alternative optical technology for detecting heterogeneities is the use of dual interfering sources that is a very sensitive technique to detect heterogeneities in diffuse media.¹⁰ This method has been suggested to be more sensitive to the presence of absorption abnormalities than single-source schemes.^{11 12} The principle of operation is based on a pair of 180°-phase-shifted sources that yield two diffuse photon density waves (DPDW), which propagate as scalar waves. The superposition of these waves in a homogeneous medium creates an ac amplitude null plane located at the mid distance from the two sources, and a 180° sharp transition in the phase at each side of this plane. The presence of heterogeneity modifies this balance pattern, which is expressed by a perturbation of the null-plane and a corresponding phase shift. These characteristics have been shown to lead to a very sensitive detection¹³ of abnormalities with important clinical potential.^14,15

Phased array systems typically scan a set of interfering sources along a surface to produce a two dimensional image of optical density (O.D.) or phase contrast.¹⁵ Despite the sensitive detection, such an approach does not allow for resolving depth. A more elaborate approach that uses single sources and performs the superposition of the detected DPDW by post-processing has been proposed by O’Leary.⁸ This scheme has been shown to enhance the quality of the reconstruction especially when the bulk optical properties of the medium are miss-estimated. On the other hand, this data processing technique is more sensitive to random and biological noise since the “interfering” measurements are not acquired simultaneously. Furthermore this scheme can work only with detectors that are on the null plane; therefore its true tomographic application is limited.

In this work we attempted to combine tomographic principles derived from DOT with dual interfering sources in order to devise a reconstruction scheme where very sensitive detection is combined with accurate localization and quantification. Using the Rytov perturbative method, we derived analytical Greens functions for the dual source scheme. These solutions were used to reconstruct measurements generated using a finite-difference-based simulation of the slab geometry used in our breast clinical examinations.¹⁶ The reconstruction performance was investigated versus object depth, source separation and frequency. Compared to the O’Leary method this approach has the advantage of utilizing detectors at non-symmetric configurations relative to the sources. Furthermore the method suggested here is more robust to amplitude and phase shift noise between the sources.¹⁷ Additionally, from a clinical point of view, the use of two sources illuminating simultaneously the medium allows the reduction of the impact of correlated biological noise such as movement due to respiration or other noise.

2. METHODS

2.1 Theory

Light follows a diffusive propagation in highly scattering-low absorbing media that is well described by the diffusion equation.¹⁹ It is well known that the P₁ approximation^20,21 leads to an accurate model of near infrared (NIR) light propagation in diffuse media for source-detector distances greater than 10 mean free paths.^22,23 This approximation, which is valid for the frequencies of interest (<1GHz) and the optical values encountered in biomedical applications (μ_a << μ′_s), yields a Helmoltz like equation that describes the propagation of the fluence in an infinite homogeneous medium, i.e.:

where U(r⃗,r⃗_s) is the fluence [W.cm^-2] (isotropic term), D the diffusion coefficient given by D = (3μ′_s)^-1 [cm], μ_a the absorption coefficient [cm^-1], μ′_s the corrected transport scattering coefficient [cm^-1], AS(r⃗_s) the source term and k the wave vector is given by $k^2=-\frac{v{\mu }_a-\mathrm{i\omega }}{\mathrm{vD}}$ (ω angular modulation frequency; v light speed in media). The photon propagation in semi-infinite and slab geometries can be derived using equation (1) and the appropriate boundary conditions.^24,25 In the case of heterogeneities, this equation is implemented with a heterogeneous operator. In this study, we considered only purely absorptive objects. This choice is justified since the main imaging applications pertaining to cancer detection and functional monitoring are realized through mapping of the absorption coefficient with or without enhancement from the contrast agent.^6,26 The heterogeneous diffusion equation can be written as:⁸

where O(r⃗) is the heterogeneous operator, i.e.:

and δμ_a (r⃗) is the spatially varying part of the absorption coefficient μ_a (r⃗) = ${\mathrm{\mu }}_{\mathrm{a}}^0$ + δμ_a (r⃗).

To solve this heterogeneous equation, the perturbative approach has been widely used. Two expansions of the diffuse field are commonly used: the Born and the Rytov.²⁷ Herein we will use the Rytov expansion as it has shown to be more effective with experimental measurements.^7,8 Under the Rytov expansion the diffuse field can be written as:

where U₀(r⃗,r⃗_s) is the solution of the homogeneous equation (1) and Φ_sc(r⃗, r⃗_s) the diffuse Rytov phase. We are going to solve the heterogeneous equation for the specific case of a dual–interfering-source system. In this case the source term is:

when the sources possess the same strength, are out of phase and are located at r⃗_s1 and r⃗_s2 respectively. The method to obtain analytical solution of the heterogeneous equation is derived from Kak.⁷ If we use the Rytov expression (4) in equation (2), we obtain:

as U₀(r⃗,r⃗_s) = e^{Φ₀(r⃗,r⃗_s)} and as the scattered signal is negligible compared to the delta sources functions, the homogeneous equation could be written as:

By subtracting this homogeneous equation from equation (6), one obtains:

or,

Therefore,

This heterogeneous equation is now dependent only on the source term U₀ (r⃗,r⃗_s), i.e. the solution of the homogeneous equation (1)). In equation (10), the source term is a sum of two delta functions. To obtain a solution in this case we use the superposition theorem, which enable us to use the Green function theorem to solve the Helmoltz equation. So, we solve independently two equations for each delta function and the total field is the sum of the two separated fields obtained:

And U₀(r⃗,r⃗_s) = ${\mathrm{U}}_0^1$(r⃗,r⃗_s1) + ${\mathrm{U}}_0^2$ (r⃗, r⃗_s2).

Obviously,

Equation (10) becomes:

By applying the Rytov approximation to equation (13) one finds:

Equation (14) is a Helmoltz type equation which can be converted to an integral equation by convoluting with the Green function solution as the left hand side of the equation (13) is linear:

Equation (15) is the heterogeneous diffusion equation, for infinite media, in the case of dual-interfering sources. This analytical expression for the diffuse Rytov phase function can be used with typical inversion methods for tomographic purposes.

2.2 Simulations

In order to solve the inverse problem, we created data sets using a finite difference solution of the diffusion equation written in the frequency domain. We chose a numerical predictor of the forward measurements to independently assess the performance of the analytical model derived in equation (15). We considered only 2D geometries for computational efficiency. The 2D simulations assume uniformity along the z-axis and suggest a source line and a detector line perpendicular to the computation plane. Nonetheless, the assessment of the analytical solutions obtained above, could be realized with 2D simulations. We further considered transmission geometry modeled after our clinical studies where the breast is slightly compressed between two plates.¹⁶ The dimensions of the medium simulated were 20×5 cm² which correspond to the typical thickness of the compressed breast. This surface was sampled by 401×101 nodes (Δx=0.5mm; Δy =0.5mm). The optical properties of the background were chosen after the typical average values of breast²⁸, i.e. ${\mathrm{\mu }}_{\mathrm{a}}^{\text{background}}$=0.05cm^-1, ${\mu \prime }_{\mathrm{s}}^{\text{background}}$ = 10cm^-1 We further assumed objects, representing tumors after contrast agent enhancement with ${\mathrm{\mu }}_{\mathrm{a}}^{\text{obj}}$=0.5cm^-1 and ${\mu \prime }_{\mathrm{s}}^{\text{obj}}$ =10cm^-1.

2.3 Reconstructions

The simulated datasets yielded a N_S × N_D set of measurements that contain the modulated amplitude and the phase for a set of dual-interfering source system at each spatial node sampling the media, where N_s is the number of source pairs and N_d the number of detectors. In each case, 64 detectors and 17 source-pairs positions were considered. The detectors and the source pairs were positioned in the 8cm-wide middle part of the forward mesh to avoid lateral boundary interference. Simulations were performed for both the heterogeneous medium considered and for a homogeneous medium with the same background optical properties as the heterogeneous medium. The diffuse Rytov phase “measured” is then given by:

For inversion, we discretized the media in N_X × N_Y voxels. The inverse problem in its matrix form can be written as:

where Φ_sc (r⃗_si,r⃗_di) is the diffuse Rytov perturbative phase for the i^th source-detector pair, W_ij is the weight function for the j^th voxel and the i^th source-detector pair, and δμ_a (r⃗_j) is the differential absorption coefficient of the j^th voxel.

The analytical weight function, for a infinite medium can be obtained from equation (15), and is expressed as:

Analytical weight functions for slab geometries can be obtained by employing the extrapolated boundary condition and image sources.²⁴ For four couples of positive and negative image source pairs the analytical form of the weight function for the case of dual interfering sources is:

where ⁺ denote positive image and ^- negative image source pairs. Generally, four image source pairs yield satisfactory accuracy.²⁹ The exact expression of the position of the image source pairs can be obtained from Arridge²⁵ or Wang³⁰ for example. Similar to the simulations performed, we considered an air-tissue boundary (n_tissue =1.33), and used the Fresnel reflection coefficient for unpolarized light²⁴ to estimate the position of the extrapolated boundary condition, at z_boundary = 1.69/(μ′_s + μ_a). It has been shown that the position of the extrapolated boundary is a crucial theoretical parameter especially for the transmission-slab geometry.³¹

The system of equations written in equation (17) was inverted using the Algebraic Reconstruction Technique⁷ (ART) with a positive constraint. More precisely, we used the Simultaneous Iterative Reconstruction Technique (SIRT) in this study, i.e.:

where ${\mathrm{x}}_{\mathrm{j}}^{\left(\mathrm{k}\right)}$ is the k^th estimate of j^th element of the object function, b_i the i^th measurement, W_ij the i-j^th element of the weight matrix and λ the relaxation parameter. This last parameter was set to 0.1 for all the reconstructions herein. This technique is known to produce smoother reconstructions but attains slower convergence than ART. To speed-up the SIRT procedure and to improve the reconstructions quality, we used a random projection access order.³²

To select the regularization parameter (i.e. iteration number), we employed numerical estimators, which compare the reconstructions to the model used. The best reconstructions were selected in accordance to these evaluators. Hence, the optimum regularization parameter is selected using prior knowledge.³² The first evaluator that we used was the correlation coefficient between the reconstruction and the model. The mathematical expression of this coefficient is:

where t_i (t¯) and ${\mathrm{x}}_{\mathrm{i}}^{\left(\mathrm{k}\right)}$ (x¯^(k)) each represent the pixel value (average value) in the original and k^th reconstructed images, respectively. The value of this coefficient ranges between 0 (no correlation) and 1 (total correlation).

The second evaluator used is the normalized root mean square error (Euclidian distance), i.e.:

The best regularization parameter was related to the maximum of ${\mathrm{\epsilon }}_1^{\left(\mathrm{k}\right)}$ and the minimum of ${\mathrm{\epsilon }}_2^{\left(\mathrm{k}\right)}$. In clinical applications where no a priori information are available, the convergence rate, which compares the difference between the measurements and the estimated measurements for two consecutive iterations, is used as the criterion of termination of the iterative process.

3. RESULTS

3.1 Synthetic measurements

The measurement vector was generated from the synthetic measurements that were obtained from a finite difference code. An example of generated measurement is presented in Fig. 1. This example shows the perturbation at detectors in transmission geometry when three heterogeneities are embedded in the media (cf. Fig. 5 (a)).

 

Fig. 1: Perturbation reaching the detection plane for 3 positions of sources-systems and for the model of Fig. 5 (a) with sources separated by 2cm and for a 50MHz frequency. The 64 detectors were located at the geometrical markers.

Download Full Size | PDF

The main characteristics of these perturbation profiles are the peaks located on the position of the null-line. For detectors that are close to the plane separating the two out-of-phase sources, the perturbation is quantitatively high. This is striking for the two extreme positions of the null-line. For the central position, the perturbations close to the null-line are smaller. For this case, the central object of Fig. 5 (a) is on the null-line and the perturbations due to the other objects are somehow compensating themselves.

Those profiles demonstrate previous observations that the phased array enhances the perturbations close to the null-plane and outline the possibility to perform highly sensitive 2D localization.^10–13

3.2 Sensitivity profiles

In order to elucidate the spatial sensitivity of the method in detecting small perturbations we plotted the weights for a single detector and source pair employing equation 19. Examples of sensitivity maps derived from our analytical models are presented in Fig. 2 (ν: modulation frequency of the sources, d: separation between the two out of phases sources constituting a pair.).

The cases presented here use the same configuration than the one presented in section 2.2. The sensitivity is displayed in terms of natural logarithm of the weight function for a specific source detector configuration.

 

Fig. 2: Sensitivity profile for a sources-detector pair: (a) ν=50MHZ – d=2cm; (b) ν=50MHZ – d=2cm; (c) ν=200MHZ – d=1cm; (d) ν=200MHZ – d=1cm

Download Full Size | PDF

In all the cases presented in Fig. 2, the sensitivity maps exhibit some common features such as the null-sensitivity line corresponding to the plane separating the two out of phase sources. This plane is characteristic to the phased array configuration. It corresponds to the plane of 180° transition of the phase.^11,33 No perturbation can be detected along this plane since the excitation pattern is cancelled.

The results further demonstrate that when the detector is placed outside the null-line the system loses sensitivity. In this case we transit from a specific “heart” shape to a more classical banana shape for the sensitivity of the technique suggesting that when the detector is shifted from the null-line, the phased array systems sensitivities are close to one of the single source systems.³⁴ This is also seen in the profiles of Fig. 1.

The sensitivity also depends on the frequency (ν) and source separation (d), as shown in Fig. 2 (a) and (d). Generally, the higher the frequency and the lower the separation, the more spatially confined is the sensitivity. However, the sensitivity retains a dissymmetric spatial distribution in contrast to the banana shape of the single source technique.

3.3 Single objects

A 4×4mm² object was embedded in the center of the medium. Its optical values were set to have a 10:1 contrast versus the background. As described above, 17 positions of the phased sources were used with 64 detectors. The sources were separated by 2cm and modulated at 50MHz. Reconstruction for the object embedded in the center of the media is presented in Fig. 3.

The reconstruction presented corresponds to the iteration giving the best reconstructions in terms of the mathematical evaluators of section 2.3 (equation (21) and (22)). After this iteration, the reconstructions are getting worse as long as SIRT is a semi-convergent technique.⁸

 

Fig. 3: One object reconstruction for a 50MHz modulation and a 2cm-separation between the sources. 17 couples of sources and 64 detectors were considered with 80×50 voxels: (a) model – (b) reconstruction.

Download Full Size | PDF

From Fig. 3, we see that the inhomogeneity was well recovered in terms of spatial and quantitative information. The reconstruction is artifact free demonstrating that the analytical model is able to recover the information.

 

Fig. 4: Correlation coefficient ε₁ and Euclidian distance ε₂ for (a)-(b) d=1cm versus frequency and (c)-(d) υ=50MHz versus sources separation.

Download Full Size | PDF

Reconstruction for several modulations frequencies (ν=50-100-200 MHz), several sources separation (d=1-2-3 cm) and several object depths (1-1.5-2-2.5-3-3.5-4 cm) were performed. Those values represent the ones used experimentally.10-13 In all cases, the heterogeneity was well recovered. However, differences between the different configurations were recorded in the quality of the reconstructions. Two examples of mathematical evaluators for different sources configuration and object depth are presented in Fig. 4.

3.4 Multiple objects

In this section we present reconstructions for imaging three objects embedded at several depths.

Three identical objects of 0.4×0.4 cm² were embedded at 1–2.5–4 cm respectively as shown in Fig. 5 (a). They were separated by 2cm in the longitudinal dimension and 1.5cm in the transversal dimension. Fig. 5 display reconstructions obtained with 17 sources and 64 detectors for different modulation frequencies and separation between the sources pair.

 

Fig. 5: Differential absorption maps: (a) Simulated – (b) d=1cm; υ=50MHz, (c ) d=1cm; υ=200MHz – (d) d=2cm; υ=200MHz.

Download Full Size | PDF

The three objects are resolved in all the reconstructions presented. That was also the case for reconstructions using the same frequencies and same source separations as described in the previous section (results not shown here). However, the three objects are recovered with different quality, which is discussed in the next section..

We present also in Fig. 6, reconstructions for the same configuration but with noise added to the synthetic measurements (c.f. Fig. 6). This uniformly distributed noise was added independently on each source of a phased pair. The reconstructions presented are for a modulation frequency of 200MHz and a separation of 2cm between the sources pair.

4. DISCUSSION and CONCLUSIONS

In this paper, analytical solutions of the heterogeneous equation for dual-interfering sources illuminating simultaneously a diffuse medium were derived. The solutions obtained with the first order Rytov perturbation approach were used to address DOT with phased-array system. Reconstructions of data simulating abnormalities embedded in human breast tissue were performed to check the validity of the analytical expressions.

 

Fig. 6: Differential absorption maps: (a) ±0.5% – ± 0.25 ° – (b) ±2.5% – ± 1 °.

Download Full Size | PDF

The accurate retrieval of the modeled maps using a separate forward model measurement generator was performed. The results obtained from these simulations in the case of single or multiple heterogeneities showed that the theoretical approach presented here to perform DOT with phased array system was able to retrieve the objects but produced irregular shapes.

In all cases, the central object was more accurately recovered in terms of size and shape. The two adjacent objects displayed irregular shapes with an evident preference in direction. We found that the higher the frequency, the more irregular the shape, especially for the object close to the source plane. On the other hand, it was this object that was more accurately quantified. In Fig. 5 (a), the central object and the one closest to the detector plane were underestimated. For higher frequencies (cf. Fig. 5 (d)), this non-uniform estimation is less apparent. For this case, the three objects were well quantified.

This effect could be due to mismatches of the two forward estimators used but also due to the non-uniformity of the weight function. Especially, the weight functions are highly dissymmetric with high sensitivity close to the boundaries. Thus as seen in Fig. 5 and Fig. 6, the object close to the boundaries suffer from this heart-shape sensitivity profiles.

The quality of the reconstructions is dependent on several parameters. Phased array techniques are based on the destructive summation of diffuse photon density wave. As interfering patterns, the destructive summation is dependent on the frequency and on the separation of the sources. The higher the frequency and the closest the source, the higher the destructive summation was. This has been seen in the sensitivity profiles. Thus, modulation frequency and source separation within each pair have an impact on the reconstruction quality. This overall statement is important in regard to the optimization of the existing apparatus. Generally, the correlation coefficient and the Euclidian distance give similar trends but slightly contradictory results. The correlation coefficient indicates that the best reconstructions were obtained for the central positions of the object. For positions around the middle of the medium, the reconstructions are worse. On the other hand, the Euclidian distance indicates that the best reconstructions are obtained for the object close to the source surface. The reconstructions are getting worse as the object is buried deeper in the media until it improves after 3cm.

One should note that these two estimators are sensitive to two different parameters. The correlation coefficient is sensitive to the spatial quality of the reconstruction and the Euclidian distance is sensitive to the quantification quality of the reconstruction. Therefore these results imply that the best location of the object is obtained for the centrals position of the object; conversely, the best estimates of the differential absorption are obtained for the position close to the boundaries. These findings were in agreement with human operator checking.

Overall, our results indicated that the higher the frequency, the better the reconstruction. However, the separation of the sources did not yield a definitive outcome. Reconstruction quality was dependent on both the relative position of the object to the source pairs and the actual source pair separation. The latter finding is related to the spatial dependence of the sensitivity of the technique and the regularization parameter chosen. In our case, the regularization parameter was the iteration number, which is a global estimate and does not discriminate the different depths. Since the sensitivity of the technique varies with depth, differences in the reconstruction quality as a function of depth are expected if the regularization parameter is fixed. This is important for clinical applications where there is not always a-priori information on the location of the suspicious mass or on the activation region. Thus, variant spatial regularization techniques³⁵ would be desirable for phased array DOT to produce more accurate distributions of the object function at various depths simultaneously.

An alternative solution could capitalize on more advanced features of the phased-array technique, i.e., the modification of the relative strength³⁶ or relative phase difference¹⁶ of the sources in order to bend the null-line. Experimental validation of this method to produce a scanning-null line has already been performed.³⁴ With this approach, one could produce an apparatus that will always optimally work with a detector on the null line, enabling more uniform sensitivity across the medium of investigation

Extension of this work to the reconstruction of scattering and fluorescent properties (lifetime and quantum yield) could be made. Also comparison to classical DOT configuration is necessary to estimate the benefit of dual-interfering sources.