Research on the migration of light through highly scattering random media has led to the realization that diffuse light can be used to probe and image spatial variations in macroscopic optical properties [1–5]. This possibility is now driving major research efforts to develop a new medical imaging technique that offers new contrast possibilities with an inexpensive, non-invasive technology. The vast literature on solving inverse problems has accelerated the development of diffuse optical tomography (DOT) algorithms for medical imaging. Most DOT algorithms are essentially matrix inversion techniques that utilize singular value decomposition with regularization [6], or algebraic techniques for minimizing a least squares difference [3, 6, 7]. More recently, non-linear methods have been developed to overcome systematic errors due to the linear approximation [2, 4]. Sufficient comparison has yet to be performed to quantify the systematic errors that non-linear methods are purported to improve upon [8]. Furthermore, it is not clear that quantitative reconstruction algorithms are absolutely necessary for all applications. For this reason, the development of new and improved linear algorithms is still an important objective. In addition, most algorithm development has looked at reconstructing images with data collected in transmission geometry, although many important applications (such as imaging brain hemorrhage and function within the adult head) only have access to reflectance data. Here we demonstrate that diffuse optical diffraction tomography, developed in transmission mode by Li et al. [9, 10] and Matson et al. [11, 12] can be extended to reconstruct psuedo-3D images from reflectance data. We present images of absorbing objects from continuous-wave (CW) experimental data collected with a single source and CCD camera. In addition to being the first demonstration of images reconstructed from CW data using diffuse optical diffraction tomography, this is an important advance for DOT as it is the first experimental demonstration that psuedo-3D images can be reconstructed from reflectance data.

The diffusion of light through a highly scattering medium on a macroscopic level can be treated as a scalar wave, called a diffuse photon density wave (DPDW), whether or not the source of light is intensity modulated [1, 13]. In a medium with spatially uniform optical properties, this wave propagates unobstructed. However, the presence of spatial variations in the optical properties results in a perturbation of the DPDW, which is accurately modeled as a scattered DPDW [14–16]. Therefore, light that diffuses through a heterogeneous medium can be described as a superposition of the incident wave and those waves scattered from local variations in the optical properties. The first Born approximation of the photon diffusion equation indicates that the wave scattered from spatial variations in the absorption coefficient is given by

ϕ_sc (r _d) is the scattered photon fluence measured at the detector at position r _d. ϕ_inc (r _s ,r) is the incident fluence at position r generated in the medium by a source at position r _s. G(r,r _d) is the Green’s function solution of the diffusion equation for the given medium. Both ϕ_inc (r _s ,r) and G(r,r _d) are calculated given the spatially uniform, average background optical properties of the medium: μ_s′, the reduced scattering coefficient, and μ_a, the absorption coefficient. D = ν/(3μ_s′) is the photon diffusion coefficient, ν is the speed of light in the medium, and δμ_a(r) is the spatial variation of the absorption coefficient at position r from the average background value μ_a. Diffuse optical tomography is essentially based on reconstructing an image of δμ_a(r) from measurements of ϕ_sc (r _d) at multiple source, r _s, and detector, r _d, positions.

Recognizing that eq. (1) resembles the Somerfield - Kirchoff scalar diffraction integral, we can use diffraction tomography to reconstruct an image of δμ_a(r) from measurements of ϕ_sc (r _d) in a plane . This diffraction model is valid in transmission mode. Neglecting this requirement, we proceed to derive an equation that is valid for reflection mode. The resulting expression should be described as an expression for reflection tomography.

We start by noting that for imaging applications where measurements are made in reflection mode, we can usually approximate the highly scattering medium as semi-infinite. In this case, placing the free-space boundary in the xy-plane at z = 0, using a collimated pencil beam incident on the medium at x = y = z = 0, and making measurements in the same plane with the source fixed, then eq. (1) can be re-written as

where we have taken the Fourier transform with respect to x_d and y_d, and the object function A(ω_x, ω _y, z), which we wish to reconstruct an image of, is given by

For the described geometry, G(r,r _d) is given by

where k=(3μ_s ' μ_a)^½ and

where R_eff accounts for Fresnel reflections at the boundary [17]. The 2D Fourier transform of G(r,r _d) with respect to coordinates x_d and y_d is thus

ϕ_inc (r _s ,r) is equal to G(r,r _d) (eq. (4)) multiplied by the source amplitude.

An analytic inversion of eq. (2) can now be obtained by assuming that the object function (eq. (3)) is zero everywhere except in a slice of width h at position z. Thus, the analytic solution for the spatially varying absorption coefficient becomes

where the inverse Fourier transform has been taken to obtain δμ_a in real coordinates and to separate the source fluence from the absorption coefficient. The assumption that the object is contained within a slice of width h is necessary, but, as we show below, this does not prevent us from reconstructing psuedo-3D images.

 

Fig. 1. Schematic illustration of the set-up for diffuse optical reflectance tomography. The phantom is 30 × 30 × 20 cm with μ_s′ = 10.0 cm^-1 and μ_a = 0.05 cm^-1. The properties of the objects are given in the text.

Download Full Size | PDF

We used the experimental geometry depicted in fig. 1 to demonstrate the utility of diffuse optical reflection tomography, as given by eq. (6). An approximately semi-infinite phantom of titanium dioxide and India ink suspended in water was created in a 30 by 30 by 20 cm container. The optical properties were determined to be μ_s′ = 10 cm^-1 and μ_a = 0.05 cm^-1 using the method described by Farrell et al. [18]. Two absorbing blocks made from polyester resin mixed with titanium dioxide and absorbing dye were placed in the phantom as absorbing objects. The blocks were 6 × 6 mm with a thickness of 3 mm aligned with the z-axis (the xy-plane is parallel with the surface of the sample). The reduced scattering coefficient of each was the same as the background medium. The absorption coefficient of the block at a depth of z = 5 to 8 mm and centered at (x,y) = (-8, 8 mm) was 0.2 cm^-1. For the block at a depth of z = 10 to 13 mm and (x,y) = (-4, -12 mm) the absorption coefficient was 0.5 cm^-1. A HeNe laser provided a collimate beam of 632 nm light at the center of the phantom at an angle of 10° from the phantom surface normal. This small angle permitted a 12 bit CCD camera with the image plane parallel to the phantom surface to image a 60 × 60 mm area with the source at the center (x,y) = (0,0).

 

Fig. 2 Difference of the images obtained by the CCD camera without and with the absorbing objects. The absorbing objects cause the measured diffuse reflectance to decrease. The color scale is linear where red corresponds to maximum change and blue corresponds to zero change. The peak attenuation for the first (second) object corresponds to a 24% (8%) change in the signal.

Download Full Size | PDF

 

Fig.3. Images reconstruction at different depths from the raw data presented in fig. 2. The objects can be localized in depth by minimizing the size of the object in the X and Y directions. The color scale is linear. Reconstruct absorption coefficients are given in the text.

Download Full Size | PDF

In fig. 2 we plot the difference of raw images obtained by the camera with and without the absorbing objects. From the raw data without the absorbing objects, we observe the expected exponential decay of the diffuse reflectance versus radial distance from the source. Fitting this experimental data with the solution of the photon diffusion equation for a semi-infinite medium, we obtain the optical properties of the medium reported above. The presence of the absorbing objects decreases the measured intensity in localized regions. The decreased intensity, indicated by the data in fig. 2, is the scattered photon fluence, ϕ_sc (r _d), used for reconstructing an image. Although the presence of two absorbing objects is clear from fig. 2, and their x and y position are accurate, we do not know the 3D position nor do we know their optical properties or size. For instance, it appears that the 2^nd (lower center) object is less absorbing than the 1^st object (upper left) when in fact the 2^nd object is more strongly absorbing but is at a greater depth.

To reconstruct an image from fig. 2, the 2D fast Fourier transform was calculated and images reconstructed using eq. (6) with h = 3 mm. The images reconstructed at depths of 2, 5, 9, and 15 mm are shown in fig. 3 (actually this is the object function A rather than the absorption coefficient δμ_a). Notice that the depth of the two objects can be determined by tuning them in and out of focus by adjusting the depth at which the image is reconstructed. The best focus is determined by minimizing the size of the object. The 1^st object is focussed best at a depth of 5 mm while the 2^nd object is best focussed at a depth of 0.9 cm. At the appropriate depth, the full width-half maximum of the reconstructed objects is 4.1 and 6.6 mm for the 1^st and 2^nd objects respectively, agreeing well with the expected value of 6 mm. The reconstructed absorption coefficients for the 1^st and 2^nd objects are 0.21 cm^-1 and 0.41 cm^-1 respectively, agreeing well with the actual values of 0.2 cm^-1 and 0.5 cm^-1. This demonstrates that diffuse optical refraction tomography enables us to fully characterize the 3D location, lateral size, and absorption coefficient of absorbing objects. We note that each image in fig. 3 requires approximately 1 second to reconstruct on a 200 MHz Pentium computer.

The images displayed in fig. 3 were reconstructed from the difference in the data collected with and without the object. In most imaging applications, it is generally not possible to obtain data without the object present. In this case, it is necessary to approximate the background data that would have been collected if the object were not present. In fig. 4 we demonstrate a method for approximating the background data which allows us to reconstruct an image from a single measurement of the sample with the object present. We start by plotting the measured reflected intensity versus radial distance from the source. In fig. 4 we plot the intensity along a line that passes through the object in the upper left, through the source at the center at the origin, and through the lower right quadrant. Given that the sample is well approximated by a semi-infinite medium, we fit the data from the lower right quadrant with eq. (4) to determine μ_s′ and μ_a. This fit gives μ_s′ =10 cm^-1 and μ_a = 0.05 cm^-1. We then calculate the diffuse reflectance for a homogeneous medium with those optical properties, subtract it from the measured data, and reconstruct images using the determined μ_s′ and μ_a for calculating ϕ_inc (r _s, r) and G(r, r _d). The images reconstructed at the depth of best focus are shown in fig 5. Note that in fig. 5a the contribution of the lower center object was neglected by replacing the raw data in that quadrant with the theoretical expectation for a homogeneous medium. The same was done in fig. 5b, but for the upper left object. Using this method, we find the absorption coefficient of the first (second) object to be 0.16 cm^-1 (0.47 cm^-1), agreeing well with the known values. Once again, we observed that the correct depth for each object could be determined by optimizing the focus of each object.

 

Fig. 4. Measured reflectance versus radial position from the source. Experimental data with the absorbing objects present given by symbols. Theoretical fit for a semi-infinite homogeneous medium is given by the solid line. The object is to the left resulting in the difference between theory and experiment.

Download Full Size | PDF

 

Fig. 5. A) Reconstruction of the upper left object by subtracting a theoretical background from the experimental data. B) Reconstruction of the lower center object by subtracting the theoretical background from the experimental data.

Download Full Size | PDF

In summary, we have shown that diffuse optical diffraction tomography can be used in reflectance mode to reconstruct psuedo-3D images of isolated absorbing objects from CW data. This method is simple and fast, as it does not require regularization parameters or multiple iterations, and images can be reconstructed from measurements made with a single source. The major disadvantage of the method is that absolute absorption images require knowledge of the thickness of the object. Although the 3D position of objects can be determined by pseudo-focussing the object, it is not yet clear how to extract the 3D size of the object. In addition, the present method cannot resolve an object that is shadowed by another object nor is it likely to image structure more complicated than isolated objects. Future improvements include incorporating measurements with multiple sources into the algorithm to enable better 3D imaging. Extension of this method to imaging spatial variations in the reduced scattering coefficient and index of refraction should be straight-forward.