1. Introduction

Optical diffusion tomography is an exciting new technology that can be used to image an object that is imbedded in a turbid medium. Typically, the turbid medium and object are illuminated with visible or near infrared light that is then scattered by the turbid medium. If the object is imbedded in the medium and has scattering and/or absorption properties different from the background medium, the scattered light in the medium is also perturbed by the object. It is this additional perturbation by the object that is used to reconstruct an image of the object. A schematic of the system geometry for the results in this paper is shown in Fig. 1, where the illuminating light comes from the left, the object is imbedded in an infinite turbid medium, and the detector is a two-dimensional detector located in the z=z_o plane. Other geometries could be considered, such as a semi-infinite medium or an infinite slab. The experimental results in this paper are generated using an essentially infinite slab, but the theory is presented for an infinite geometry to simplify the mathematics. The spatial distribution of the light in the turbid medium, ϕ(x, y, z) , is given by

where, for any point (x,y,z) in the turbid medium, ϕ_hom (x, y,z) is the spatial distribution of the light in the absence of the imbedded object (the homogeneous or background light) and ϕ_pert (x, y, z) represents the perturbation of the light due to the presence of the object. In the case of an absorptive object, under the Born approximation [1], the light perturbed by the object can be written as [2]

where o(x, y,z) is a 3D representation of the absorption properties of the object and h(x, y,z) is the Green’s function solution to the photon diffusion equation for the homogeneous turbid medium and boundary conditions. Equation (2) shows that the light perturbed by the object is given by the 3D convolution of the imbedded object (weighted by the homogeneous light) with the Green’s function. From Eqs.(1) and (2), the signal collected by the detector in the z=z_o plane, ϕ(x,y,z_o), can be considered as the superposition of two components: a blurred component ϕ_pert(x,y,z_o) from the light that is scattered by the object, and a background component ϕ_hom(x,y,z_o) from the incident light scattered by the turbid medium. The inverse problem associated with imaging the imbedded object requires removing both the background component ϕ_hom(x,y,z_o) and the blurring h(x,y,z) from ϕ_pert(x,y,z). Matson and Liu [3] have shown theoretically and experimentally that a backpropagation algorithm can be used to solve the image blur part of the problem when the light scattered by the object is measured in a plane external to the object (experimentally, external to the medium as well), such as is done with CCD cameras. In this case, under the assumption that the detection plane is perpendicular to the z axis and located at z=z_o (see Fig. 1), Eq. (2) can be rewritten as [4]

 

Fig. 1. Schematic showing the geometry of the optical diffusion tomography system for the development in the paper. The imbedded object is represented by the cube.

Download Full Size | PDF

where

and where * is the 2D convolution operator with respect to x and y, z₁ is the largest z value inside the object support, and the depth parameter Δz = z_o-z₁ is the distance between the part of the object closest to the detection plane and the detection plane. The difference between Eq.(2) and Eqs.(3) and (4) is that the integration indicated in Eq.(2), which is over all space, is broken up into two pieces. The first piece is the integration over the object support (Eq.(4)), and the second piece is the propagation of the light perturbed by the object (o_2D(x,y,z₁)) to the detection plane by using the angular spectrum solution to the homogeneous source-free scalar wave propagation [5]. The angular spectrum propagation operation is indicated by the convolution of o_2D(x,y,z₁) with h(x,y,Δz), where the latter term is the two-dimensional inverse Fourier transform of the angular spectrum transfer function for the distance Δz and is referred to as the system point spread function (PSF) throughout this paper. Equation (3) indicates that the image reconstruction problem can be viewed as 2D deconvolution, where the desired reconstructed image is o_2D(x,y,z₁); that is, the desired reconstruction is the projection of the 3D object structure into the z=z₁ plane. To obtain the absorptive properties of the object from the deconvolved data requires solving Eq.(4) and requires multiple views, in general [3]. Matson and Liu [3] also showed that a backpropagation algorithm can be used to determine Δz from the data itself instead of having to know this value a priori. Unfortunately, like many reconstruction methods, the backpropagation method requires prior knowledge of the structure of the mathematical model governing light propagation in turbid media, something that is not always available. A second disadvantage to the backpropagation algorithm is that, in order to tell what value of depth parameter provides the best image, the user needs to have some idea of what the object should look like. Blind deconvolution is an alternate reconstruction method that does not have these restrictions. The method uses generic properties of the measured data, such as its band-limited nature, the convolution relationship between the measured data and the object perturbing the light, and the fact that light intensities are positive, to constrain the reconstruction of the measured data. In this paper we discuss some preliminary results from applying blind deconvolution methods to optical diffusion tomography data.

2. Blind deconvolution algorithm

Blind deconvolution is a method for determining the characteristics of both the object being imaged and the system PSF from a single blurred image of the object. In its most basic form, this is an ill-posed, undetermined problem that has a multitude of solutions, not all of which are physically meaningful. However, by applying prior knowledge about the imaging process to the problem, it is possible to constrain the set of solutions to a subset that are feasible and thus useful.

There are many different mathematical techniques for blind deconvolution, each with their own strengths and weaknesses. Here we use an iterative approach that is based on a forward model of the imaging process and a maximum-residual-likelihood [6] determination of the optimal model parameters; that is, we find the object [o_2D(x,y)] and PSF [h(x, y)] intensity distributions that minimize the cost function

where r(x, y) are weighted residuals as given by

Here d(x, y) is the recorded image, s_d(x, y) is a binary function that masks out bad pixels or regions of low signal-to-noise ratio in the recorded image, n(x, y) is the detector read noise, and ⊗ is the correlation operator. The weighting term (the denominator) allows for both Poisson and Gaussian noise [7]. We have dropped the term z₁ in the argument list of o_2D and the term Δz in the argument list of h for clarity. The iterative approach uses a conjugate gradient routine to find the intensity distributions that minimize ε.

When there is little or no information available about the object or PSF intensity distributions, it is necessary to model them in the most basic way: i.e., sample the distributions at the same resolution as the measurement image, and treat each sample (or pixel) as a free parameter. The standard constraints used by most blind deconvolution algorithms, including ours, are that: 1) the object and PSF are positive functions that extend, or at best are observable, only over a finite extent of image space, and 2) the PSF is band limited. We enforce these properties through the reparameterizations

and

Here s_o(x,y)and s_h(x,y)are binary functions that designate the spatial extent of o_2D(x,y)and h(x,y),and θ(x,y) and ψ(x, y) are the parameters to be optimized. The PSF reparameterization includes a correlation (band-limit) function, τ(x,y), [8] that is a low-pass filter whose cutoff frequency is chosen to exclude all the noise-dominated Fourier components of the measured data.

In some situations, however, we may have sufficient additional prior information to model the object and/or PSF with substantially fewer free parameters than are required by the pixel-by-pixel method. This is the case for imagery obtained through turbid media where the diffusion equation can be used for modeling the scattered light propagation [1,2]. We thus use a turbid medium optical transfer function (OTF), which is based on the diffusion equation in an infinite medium, as a means for modeling the PSF. For continuous-wave (CW) illumination, the turbid medium OTF is given by [9]

where α depends on the scattering and absorption characteristics of the turbid medium. The incoherent, band-limited PSF is given by the inverse Fourier transform

where T(u, v) is the Fourier transform of τ(x, y) and N² is the number of pixels in the image. In practice, we model the PSF via the turbid medium OTF during the early part of the iteration, and then change to the pixel-by-pixel PSF model when the object reconstruction appears to be fairly well established. The reason for this is that Eq. (9) is quantitatively valid only for an infinite and homogeneous medium, and we are inherently dealing with a finite, inhomogeneous medium (due to the presence of an object). By switching to the pixel-by-pixel model during the later iterations, we help avoid artifacts in the reconstructed object due to shortcomings in the PSF model.

Equation (9) shows that the theoretical OTF has zero Fourier phase, implying that the Fourier phase of the measured data is a noisy realization of the object’s Fourier phase. We have thus found it beneficial to: 1) use the measured data as the initial estimate of the object, and 2) replace the Fourier phase of the reconstructed object, after a few iterations, with the Fourier phase of the data (over the region where the measured Fourier phase signal is above the noise level) and then continue the minimization.

3. Experimental setup

The experimental data were collected using the transmission geometry¹ setup described in Matson and Liu [11] and shown in Fig. 2. Basically, a CW laser was used to illuminate a target that was either set against a turbid medium phantom or was imbedded in one. The turbid medium phantoms were created using a resin mixed with TiO₂ to give scattering properties, and a dye to give absorption properties [11]. The absorption (μ_a) and reduced scattering (μ’_s) coefficients for the phantoms were 0.01cm^-1 and 18 cm^-1, respectively. The scattered light was recorded using a SpectraSource 16 bit 512×512 CCD camera and stored in a computer. The targets are shown in Fig. 3. The blurred images of the targets were produced by subtracting from the blurred target measurement a second measurement of the same turbid media without a target present. We realize that it is generally not practical to measure and remove the background component of the measured signal, either in the field (surveillance), or in vivo (biomedical). As such, this approach has limited applicability. However, it does allow us to make a preliminary assessment of blind deconvolution as a potential tool for post processing imagery obtained through a turbid medium.

 

Fig. 2. Schematic of the experimental setup used for the data acquisition.

Download Full Size | PDF

 

Fig. 3. The different targets used to generate the experimental data. The first two targets, from left to right, were positioned behind a 2.5 cm thick phantom. The third target (far right) was imbedded in the center of a 5.5 cm thick phantom. From left to right: F-111 model aircraft, F-18 model aircraft, and Boeing 747 model aircraft. The F-111 and F-18 aircraft models, when illuminated and viewed through the turbid phantom, produce similar images (see Fig. 4). This provides a stringent test of the resolving capability of the two restoration algorithms considered in this work.

Download Full Size | PDF

4. Reconstruction results

Using the algorithm described in section 2, along with different input guesses for Δz and α, we have found that there are several {Δz, α} pairs that produce almost identical PSF profiles, and thus essentially the same values for the cost function. This degeneracy in the PSF is reminiscent of what Carasso [12] found when using Levy distributions to describe the system OTF in his blind deconvolution algorithm. However, for high signal-to-noise ratio imagery (such as that presented here), we find that the degeneracy can be reduced by penalizing solutions with spatial frequencies beyond the cutoff frequency observed for the measured data. We do this by adding the penalty term [13]

to the cost function ε in Eq.(5). Figure 4 shows the reconstructions for the targets shown in Fig. 3 using our blind deconvolution (BD) algorithm, and the backpropagation (BP) algorithm of Matson and Liu [11]. We note that the OTFs for the measured data fall off extremely quickly with increasing spatial frequency. Thus, the best reconstructions that can be expected are low-frequency versions of the truth. In addition, because the targets can be considered as purely absorptive, we will only measure their shadow. Figure 4 shows that the morphology of the blind deconvolution restorations is better (visually) than that of the BP restorations. In particular, although both the BP and BD reconstructions allow the user to distinguish between the F-111 and F-18 targets, the BD result more clearly differentiates the tail structures of the two aircraft. The strong artifacts and distortion of the morphology in the BP reconstructions are a manifestation of both the type of spatial filter and the cutoff frequency used to regularize the inversion process [3]. By comparison, the BD algorithm is able to reduce the artifacts and distortion through the inclusion of prior knowledge about the imaging process into the reconstruction procedure.

From their work using BP techniques, Matson and Liu [3] concluded that the quality of the reconstructed image is directly related to the accuracy of the depth parameter used in the BP inverse filter. The best image quality is obtained when the measured data are inverted using the OTF calculated for the correct depth. This suggests that the OTF parameters estimated by the BD algorithm have the potential to provide information on the object location (depth) and the properties of the scattering medium. Our mean estimated value for Δz for the case of the targets located behind the turbid medium phantom, <Δz> = 2.3±0.2 cm, is close to its expected value of 2.5cm (the thickness of the phantom). In addition, the estimated value of Δz for the (imbedded) Boeing 747 data, Δz=2.6 cm, is in good agreement with its true location [11]. These results are encouraging as the depth values that are obtained using BP are strongly dependent on the choice of regularization filter [3]. However, our estimated values of α are very different from the expected value (3μ’_sμ_a). This is most likely due to the fact that Eq. (9) is based upon the assumption of an infinite background medium, while the data were generated using a finite slab for the background medium. The OTF model for use with a finite background medium requires additional terms to account for the boundaries of the medium.

 

Fig. 4. Columns, left to right: raw data, reconstructed object, o_2D(x,y;z₁), using BP, and BD respectively. Top row: F-111 model airplane, second row: F-18 model airplane, bottom row: 747 model.

Download Full Size | PDF

5. Discussion

Our first results from applying BD techniques to imagery obtained through a turbid medium are encouraging. Blind deconvolution appears to be capable of providing a high-quality reconstruction of the object, and a good estimate of its location, without prior knowledge of either the characteristics of the turbid medium, or of the nature of the object. However, more work needs to be done to see if the reconstructed PSF can be a source of information about the properties of the turbid medium. In particular, the model for the turbid medium OTF needs to be improved.

As noted in the section on the experimental setup, before attempting any image reconstruction, we preprocessed the measured data to remove the component of the detected signal that was due to the light scattered by the turbid medium so as to leave just the light scattered from the target. We realize that it is generally not practical to do this, either in the field (surveillance), or in vivo (biomedical). In principle, simultaneous estimation of the homogeneous wave should be possible as it only requires the estimation of one extra parameter. Thus, we are currently investigating how to adapt our BD algorithm to estimate the background signal as well.