When light is incident on a strongly scattering medium like tissue or clouds, the ballistic light is quickly attenuated through scattering, and imaging capabilities are generally assumed to be lost. Recently, there have been significant research efforts to model the transport of light through strongly scattering media in order to characterize the optical properties of spatially uniform media and to image spatial variations in optical properties [1–5]. Imaging methods have been developed using continuous-wave light (CW domain) [6], pulsed lasers (time domain) [7,8], and intensity modulated light (frequency domain) [9–11]. Many imaging algorithms rely on assumptions which linearize the relationship between the optical contrast of the inhomogeneity and measurements of the scattered light. In this paper, we show that this linearization results in incorrect moments and thus accurate quantitative imaging is not possible. We also discuss the implications that this result has on Diffuse Optical Tomography (DOT); primarily that quantitative DOT requires non-linear reconstruction algorithms which correctly model the moments. Although the results presented here are based on an analysis in the frequency domain, the results also apply to the time domain (as the two are related by a Fourier transform) and the CW domain (which is the zero frequency limit of the frequency domain).

When light traveling through a highly scattering medium is more likely to experience a scattering event than to be absorbed, the transport of light is well-modeled by the photon diffusion equation [12,13]

Here, D=v/(3μ_s ’) is the photon diffusion coefficient [cm²/s], v is the speed of light in the medium [cm/s], μ_s ’ is the reduced scattering coefficient [cm^-1], μ_a is the absorption coefficient [cm^-1], Φ(r,t) is the light fluence rate [W cm^-2], and S(r,t) is the light source term [W cm^-3]. When the source is intensity modulated, the diffuse light generates a Diffuse Photon Density Wave (DPDW) [14,15]. This DPDW behaves like a classical wave in that it refracts [14], diffracts [16], and scatters [17–19] from spatial variations in the optical properties of the medium. From measurements of the scattered DPDW and utilization of an inverse scattering algorithm, it is possible to reconstruct images of the spatially varying optical properties. The angular moments of the scattered DPDW are related to the optical properties of the inhomogeneities scattering the DPDW. It has been shown that full characterization of the inhomogeneity, including size and optical properties, requires detection of at least the monopole, dipole, and quadrupole moments with a signal-to-noise ratio greater than one[20].

An analytic solution of the diffusion equation exists for the scattering of DPDW’s from spherical objects embedded in otherwise homogeneous media. To leading order in the radius (a) and optical properties of the object(μ_s’ and μ_a), the monopole (Φ^l=0(r,t)), dipole (Φ^l=1(r,t)), and quadrupole (Φ^l=2(r,t)) moments of the scattered wave are [20]:

Here, δμ_a is the difference in the absorption coefficient between the object and the background medium (or delta absorption coefficient), δμ_s’ is the difference in the reduced scattering coefficient, k ²=(-vμ_a +iω)/D, ω is the source angular modulation frequency, μ_a is the absorption coefficient of the background medium, D is the photon diffusion coefficient of the background medium, a is the radius of the object, and r_s , r_d , and θ are defined in fig. 1. Note that to leading order the monopole depends on the product of δμ_a and a ³, and the dipole depends on the product of δμ_s ’ and a ³. Thus for an absorbing object, it is not possible to simultaneously determine the delta absorption coefficient and the radius if only the monopole term is detectable. Simultaneous determination requires, in the least, detection of two different moments. Likewise for a scattering object, at least two moments must be detectable in order to simultaneously determine the reduced scattering coefficient and the radius. For a general inhomogeneity with absorption and scattering contrast, at least three moments need to be detectable to permit full characterization of the object. These results are detailed in Boas et al. [20].

 

Fig. 1 The geometry used for calculating the scattered DPDW (see text).

Download Full Size | PDF

Common examples of linearized image reconstruction algorithms for DOT are based on the Born approximation [5,9,11] and Rytov approximation [9,21]. The Born approximation assumes that the DPDW is given by Φ(r _s,r _d) = Φ _inc (r _s,r _d) + Φ _sc(r _s,r _d) while the Rytov approximation assumes that the DPDW is given by Φ(r _s,r _d) = Φ _inc(r _s,r _d) exp(Φ _sc(r _s,r _d)). Φ _inc(r _s,r _d) is the DPDW measured without any inhomogeneities in the optical properties for a light source at position r _s and a detector at position r _d. The scattered DPDW is given by Φ _sc(r _s,r _d). For the Born approximation,

while for the Rytov approximation,

L(r) is an operator which depends on the spatially varying absorption coefficient and the spatially varying reduced scattering coefficient[9]. Images of the spatially varying optical properties are reconstructed given measurements of Φ _sc(r _s,r _d), calculations of Φ _inc(r _s,r _d) and G(r _s,r _d) from the known background optical properties, and inverting eq. (3) or eq. (4) to find L(r). G(r _s,r _d) is the Green’s function solution of the photon diffusion equation for the spatial uniform background optical properties. This inversion generally requires discretizing the integral equation to obtain a matrix equation which is then solved by matrix inversion or algebraic reconstruction techniques.

 

Fig. 2 Comparison of the first order Born approximation with the exact solution for a 1 cm diameter absorbing object with different absorption coefficients.

Download Full Size | PDF

Quantitative imaging is possible only if the linearized solutions, eqs. (3) and (4), produce scattered waves which agree well with exact solutions. The fact that sensitivity to the moments of the scattered wave is critical to quantitative characterization of optical inhomogeneities, as described above, suggests that we should analyze the moments of the linear solutions used for image reconstruction algorithms. To do this, we calculated eqs. (3) and (4) for a spherical object embedded in an otherwise uniform, infinite medium, and compared the numerical results with the analytic solution for DPDW’s scattering from spherical objects [17–19]. The geometry used for the calculations is indicated in fig. 1. The optical properties of the infinite medium were μ_s ’=10 cm^-1 and μ_a=0.05 cm^-1. The source was positioned 3 cm along the z axis from the center of the object, which was placed at the origin. The source had a modulation frequency of 200 MHz (or angular modulation frequency of 2π 200 MHz). Detectors were placed every 2° in a circle 2.5 cm from the center of the object. The spherical object had a fixed reduced scattering coefficient, μ_s’=10 cm^-1, and a varying absorption coefficient and radius (0.1 ≤ μ_a ≤ 1.0 cm^-1 and 0.1 ≤ a ≤ 1.0 cm). The integrals in eq. (3) and eq. (4) were evaluated by dividing the bounding cube of the spherical object into 10x10x10 voxels and summing the contribution from each voxel. The moments were calculated from the projection of the detected complex scattered wave onto spherical harmonics Y_l.m(θ) for m=0 (no dependence on azimuthal angle). We do not explore the dependence on the azimuthal index m as our geometry was chosen for the moments with m≠0 in order to have zero amplitude.

 

Fig. 3 Comparison of the first order Born approximation with the exact solution for an absorbing object with an absorption coefficient of 0.3 cm^-1 and different radii.

Download Full Size | PDF

We first consider the Born approximation, eq. (3). In fig. 2 we plot the relative amplitudes and phases of the monopole, dipole, and quadrupole moments of the scattered wave for the Born approximation and the exact solution, to show their differences. Results are graphed versus the absorption coefficient of the object for a 1.0 cm diameter object.

When the absorption coefficient of the object is 0.10 cm^-1, the amplitude (phase) of the exact and Born moments agree within a few percent (1 degree). As the object absorption coefficient increases, the amplitudes disagree by more than 20% and the phases disagree by more than a few degrees.

In fig. 3 we plot the relative amplitudes and phases of the monopole, dipole, and quadrupole moments of the scattered wave versus radius for an object with an absorption coefficient of 0.3 cm^-1. Similar discrepancies are observed. The results in fig. 2 and fig. 3 clearly demonstrate that the first order Born approximation does not accurately predict the amplitude and phase of the scattered DPDW. Thus, image reconstruction algorithms based on the first order Born approximation will not accurately quantify the optical properties of the medium.

 

Fig. 4 Comparison of the Rytov approximation with the exact solution for a 1 cm diameter absorbing object with different absorption coefficients.

Download Full Size | PDF

The trend observed in fig. 2 for the Born moments is interesting. As the object absorption coefficient increases, the amplitudes of the moments increase proportionally and the phases of the moments remain fixed. This result is expected. Within the Born approximation, a voxel of the discretized integral produces a monopole scattered wave (given a voxel with an absorption coefficient different than the background). The strength of this monopole wave is linearly proportional to the absorption difference between the voxel and background. Doubling the absorption difference will double the monopole moment. An object comprised of a single absorbing voxel will not produce higher order moments (e.g. dipole and quadrupole). It is only from a group of absorbing voxels that these higher moments are produced. Furthermore, as the existence of these moments depends only on the spatial distribution of voxels, increasing the delta absorption coefficient by a multiplicative factor X will increase the amplitudes of the moments by a factor X. Thus, the relative amplitude of the moments is not dependent on the absorption coefficient (contrary to that observed with the exact solution). Furthermore, the phase remains constant because of this linear dependence on the absorption difference.

It is useful to consider what effect this discrepancy between the Born approximation and the exact solution has on image reconstructions. Reconstruction algorithms typically find the image which gives the best fit to the experimental data by means of minimizing the least-squares difference between the experimental measurement (in this case the exact solution) and the model used by the algorithm (in this case the first Born approximation). This is not strictly true as some algorithms use regularization methods that also minimize a norm of the image [7]. Consider a 1.0 cm diameter spherical absorber with an absorption coefficient of 0.3 cm^-1. From fig. 2 and fig. 3 we see that minimizing the least-squares difference of the amplitude can be accomplished by decreasing the size and/or absorption coefficient of the object. Or the size can be overly decreased and compensated by an increase in the absorption coefficient, or vice-versa. To minimize the least-squares difference in the phase, the size of the object must increase, as changing the absorption coefficient has no effect. This suggests that the reconstruction algorithm would reconstruct an object which has a larger diameter and smaller absorption coefficient than the real object parameters. More insight on specific algorithms operating on particular clinical problems could be obtained by correlating the performance of the algorithm with observations like those described in this paragraph.

 

Fig. 5 Comparison of the Rytov approximation with the exact solution for an absorbing object with and absorption coefficient of 0.3 cm^-1 and different radii.

Download Full Size | PDF

We next consider the Rytov approximation, eq. (4). In fig. 4 we plot the relative amplitudes and phases of the monopole, dipole, and quadrupole moments of the scattered wave for the Rytov approximation compared with the exact solution. Results are graphed versus the absorption coefficient of the object for a 1.0 cm diameter object. For an object absorption coefficient of 0.10 cm^-1 there is an amplitude discrepancy of a few percent and a phase difference of 1 to 3 degrees. These differences increase rapidly with an increasing absorption coefficient. In fig. 5 we plot the relative amplitudes and phases of the monopole, dipole, and quadrupole moments of the scattered wave versus radius for an object with an absorption coefficient of 0.3 cm^-1. Similar discrepancies are observed.

The results in fig. 4 and fig. 5 clearly demonstrate that the Rytov approximation does not accurately predict the amplitude and phase of the scattered DPDW and thus image reconstruction algorithms based on the Rytov approximation will not accurately quantify the optical properties of the medium. The major difference between these Rytov results and the Born results presented in fig. 2 and 3 is that the dipole moment calculated by the Rytov approximation has a discrepancy about four times greater than that of the Born approximation. It is interesting to note that the fractional difference in the amplitude is nearly the same for all moments suggesting that the Rytov approximation may be empirically modified.

These results indicate that linear algorithms for diffuse optical tomography cannot accurately quantify spatially varying optical properties. The solution to this problem is to use non-linear reconstruction algorithms which implement a forward problem solver which accurately calculates the moments of the scattered DPDW. Examples of such algorithms include finite-element and finite-difference solutions of the differential equation, perhaps using multi-grid methods [5,22]. Solving the full non-linear problem is computationally expensive such that near real-time image reconstruction would require super-computers. It is desirable to find less costly methods that gain a speed advantage by not calculating the full non-linear forward problem, but still maintain accuracy.

One possible alternative is to use the N^th order Born approximation. This alternative is only desirable if accurate solutions can be found with less computational expense than solving the full non-linear problem. The applet in fig. 6 allows us to determine the accuracy of the N^th order Born approximation by comparing the moments calculated with the N^th order Born approximation and the exact solution. Results are presented for 1 ≤ N ≤ 10, 0.1 ≤ a ≤ 1.0 cm, and 0.1 ≤ μ_a ≤ 1.0 cm^-1. This four dimensional set of data is viewed in two dimensional slices, and the reader can interactively move the slice through the four dimensional space by clicking on the appropriate buttons. The data files used for generating these graphs are available through http links from the applet.

Figure 6 initially shows amplitude versus object absorption coefficient for a 6 mm diameter object. The amplitude difference between the Born and exact calculations is greater than 10% for the monopole and quadrupole moments for an object absorption coefficient greater than 0.3 cm^-1. This is seen by looking at the amplitude ratio (${{\Phi }_{\mathit{\text{Born}}}}^l$ -${{\Phi }_{\mathit{\text{Exact}}}}^l$ )/${{\Phi }_{\mathit{\text{Exact}}}}^l$ . (Note that there is a systematic difference between the exact and Born and Rytov calculations because of the difference in object volume arising from dividing a sphere into a group of cubes. This systematic difference is generally small compared to the intrinsic difference between the Born approximation and the exact solution for small N but does prevent us from obtaining perfect agreement as we increase N.) When N is 5 the difference is less than 10% for absorption coefficients less than 1 cm^-1. The fifth order Born approximation provides a significant improvement. As we increase the diameter of the object to 1 cm the agreement remains good for absorption coefficients less than 0.5 cm^-1, but the agreement quickly deteriorates for larger absorption coefficients. As N increases, the agreement at an absorption coefficient of 1.0 cm^-1 gets exponentially worse. This can be observed by viewing the amplitude ratio versus N. It appears as if an instability occurs in the N^th order Born approximation. Looking at the graph of amplitude ration versus N in fig. 6, we observe that this instability sets in when the product of the radius and absorption coefficient is between 3 and 4. Similar behavior is observed for the phase data. If an imaging algorithm were to be created based on the N^th order Born approximation, then this behavior needs to be better understood. This discrepancy is most likely a result of ignoring the self-interaction term described by Ostermeyer and Jacques (see eq. (13) of [23]). The message from fig. 6 is that an imaging algorithm based on the fifth order Born approximation would be suitable for object diameters less than 8 mm and absorption coefficients less than 0.4 cm^-1.

 

Fig. 6 This is a view of the interactive applet for browsing the four dimensional data set. The data here is for an absorber with a 0.6 cm diameter. Click on the figure to start the applet. [Media 1]

Download Full Size | PDF

An imaging algorithm based on an N^th order Born approximation is likely to follow the same strategy used in other non-linear approaches that use finite-difference or finite-element schemes. That is: (1) an initial guess of the spatial variations in the optical properties is made; (2) the fluence rate is calculated for that initial guess using non-linear method X (where X is finite-difference, finite-element, N^th order Born, etc.); (3) the difference between the calculated and measured fluence rates is used in an integral equation like eq. (3) or eq. (4) to solve for an update in the initial guess of the optical properties; (4) iterate back to step 2) until some convergence criterion is met. Step (2) is performed using a so-called forward problem solver (such as finite-difference, finite-element, or N^th order Born) and step (3) needs a so-called inverse problem solver. Each step is generally computationally expensive. The Nth order Born approximation requires approximately NxM² multiplications where M is the number of voxels in the image. For comparison, multi-grid methods using finite-difference or finite-element schemes generally require M ln (1/ϵ) multiplications where e is the desired accuracy [24]. Solving the forward problem is therefore much more efficient with a multi-grid method than with the N^th order Born approximation.

In summary, we have shown that linearized algorithms for diffuse optical tomography do not accurately model the amplitude and phase of scattered DPDW’s. Quantitative imaging is thus not feasible as it relies on the use of accurate models for predicting the amplitude and phase of the scattered DPDW. The most common linearized algorithms are based on the first order Born approximation and the Rytov approximation. Each approximation poorly models the scattered DPDW. The Rytov approximation calculates the monopole and quadrupole moments scattered from an absorbing object more accurately than the Born approximation, but the Born approximation is much more accurate for the dipole moment. The main consequence of the results presented here is that nonlinear imaging algorithms must be used to achieve quantitative imaging. We have shown that higher order Born approximations can more accurately calculate the scattered wave. When the absorption coefficient is less than 0.4 cm^-1 and the diameter is less than 8 mm, the fifth order Born approximation generally provides sufficient accuracy. However, implementing an N^th order Born approximation in a non-linear imaging algorithm is much more computationally expensive than using multi-grid finite-difference or finite-element methods.