1. Introduction

The ability to accurately define the absorption and scattering properties of tissue could significant add to the diagnostic sensitivity and specificity of optical measurements. In the case of imaging studies, a common problem is inter–parameter crosstalk [1]. This refers to instances where, for example, localized variations in absorption also appear, in the recovered image, as localized variations in scattering, or vice versa. This issue can arise from fundamental reasons related to the intrinsic information content of a data set, as well as from reasons related to the numerical methods used to reconstruct parameter maps. In the case of DC measurement data, Arridge and Lionheart [2] have claimed to have rigorously proven that there is an underlying non–uniqueness in the inverse problem, and have strongly emphasized that total intensity measurements alone are insufficient to distinguish effects of absorption from scattering.

In this report we describe a new algorithm (i.e., the normalized–constraint algorithm) for DC imaging, and we demonstrate the ability to separate local perturbations in absorption and in scattering, under a wide range of conditions, as evaluated by numerical simulations and experimental laboratory studies. These results clearly demonstrate that, contrary to previous assertions [2,3], DC imaging methods can be applied to characterize spatial variations in the absorption and scattering properties of highly scattering media.

2. Methods

The normalized–constraint algorithm for minimizing inter–coefficient crosstalk employs a three–step process. First, we work with relative detector readings instead of absolute values, based on a previously developed inverse formulation called the normalized difference method [4]. We do this in recognition of the practical limits imposed on obtaining reliable absolute measurement data from arbitrary structures such as tissue. We also do this in recognition of the fact that perturbation methods, in general, tend to yield grossly incorrect solutions if the reference medium used to generate the initial guess differs sufficiently from the actual target background properties. As described in recent reports [4–6], we have shown that selection of insufficiently accurate reference media can severely alter the information content of the data vector. Once corrupted, the recovery of images relatively free of artifact can be very difficult or impossible, even with use of full Newton updates. The second step scales the weight matrix, by normalizing the column vectors to their respective mean values, to obtain a transformed formulation. This makes the weight matrix in the new system equation more uniform and less ill–conditioned, and also serves to suppress numerical errors and accelerate convergence [7,8]. The third step employs a constrained CGD method for solving the resulting system equation, where positivity or negativity constraints are imposed on iteratively computed estimates of the medium’s absorption and diffusion coefficients. Ordinarily, this option is not applicable to the general case wherein the sign of the perturbation is unknown. As is described below, we have found that by adopting a two–step process where solutions are obtained for both signs of the constraints and then summed, satisfactory solutions can be obtained. The details of this methodology are described subsequently.

2.1 Forward model

Light propagation in a scattering medium was modeled as a diffusion process. For a domain having a boundary ∂Λ, this is represented by the expression

where u(r) is the photon intensity at position r, r _s is the position of a DC point source, and D(r) and µ_a (r) are the position–dependent diffusion and absorption coefficients, respectively. Here we define the diffusion coefficient as D(r)=1/{3[µ_a (r)+µ′_s (r)]}, where µ′_s (r)) is the reduced scattering coefficient.

2.2 The inverse formulation

The optical inverse formulation is based on the normalized difference method [3–6] and has the following form:

where δµ_a and δD are the vectors of cross–sectional differences between the optical properties (absorption and diffusion coefficients, respectively) of a target (measured) medium and of a reference (computed or measured) medium used to generate the initial guess; ${\mathbf{W}}_r^{\left({\mu }_a\right)}$ and W^(D)_r are the weight matrices describing the influence that perturbations in the absorption and diffusion coefficients, respectively, of the selected reference medium have on the surface detectors; and δu_r represents a normalized difference between two sets of detector readings, which is defined by the equation

Here, u _r is the computed detector readings corresponding to the selected reference medium, u ₂ and u ₁ represent two sets of measured data (e.g., background vs. target, time-averaged mean vs. a specific time point, etc.), and M is the number of source–detector pairs in the set of measurements.

2.3 Weight matrix scaling

We have previously described a scaling method that serves to improve the conditioning of the absorption weight matrix [7]. Here we extend this method to the case in which absorption and diffusion coefficients are recovered simultaneously. The effect of scaling the weight matrix is to make it more uniform, which can often have the serve to improve its conditioning. Whereas a variety of scaling approaches could be adopted, we have chosen one that scales each column of ${\mathbf{W}}_r^{\left({\mu }_a\right)}$ and ${\mathbf{W}}_r^{\left(D\right)}$ to the average value of the column vector. The form of the resulting new weight matrices is

where, k can be µ_a or D, and R ^(k) is the normalizing matrix whose entries are

in which N is the number of elements used in discretizing the domain Λ. The resulting system equation is

where $\delta {\tilde{\mu }}_a={\left[{\mathbf{R}}^{\left({\mu }_a\right)}\right]}^{-1}·\delta {\mu }_a$ and $\delta \tilde{\mathbf{D}}={\left[{\mathbf{R}}^{\left(D\right)}\right]}^{-1}·\delta \mathbf{D}$.

2.4 Solutions of constrained CGD method

For most measurement geometries the weight matrix W̃ _r is not symmetric positive definite, and the usual practice is to seek a least-squares solution of the system of linear equations shown in Eq. (3). This solution is found by minimizing the mean-squared error E, which is represented as

where δx̃ ^T=[δ${\stackrel{\mathit{~}}{\mathit{\mu }}}_a^{\mathrm{T}}$ δD̃ ^T], ${\tilde{\mathbf{W}}}_r=\left[{\tilde{\mathbf{W}}}_r^{\left({\mu }_a\right)}{\tilde{\mathbf{W}}}_r^{\left(D\right)}\right],$, A=${\stackrel{\mathbf{~}}{\mathbf{W}}}_r^{\mathrm{T}}$ ·W̃ _r and b=${\stackrel{\mathbf{~}}{\mathbf{W}}}_r^{\mathrm{T}}$ ·δu _r . Formally, the least-squares solution is obtained by setting the derivative of E to 0, i.e.,

When the CGD method is applied, instead of explicitly solving Eq. (8), the following iterative formulation is employed for computing δx̃ :

The iterative procedure is:

1) Based on an initial estimate of δx̃ ⁽⁰⁾, compute g ⁽⁰⁾=A·δx̃ ⁽⁰⁾-b;d ⁽¹⁾=-g ⁽⁰⁾ ; β ⁽¹⁾=0;

2) For the n^th iteration, compute ${\alpha }^{\left(n\right)}=\frac{\parallel {\mathbf{g}}^{\left(n-1\right)}{\parallel }^2}{\parallel {\tilde{\mathbf{W}}}_r·{\mathbf{d}}^{\left(n\right)}{\parallel }^2},$, where g ^(n-1)=A·δx̃ ^(n-1)-b (=g ^(n-2)-α ^(n-1) A·d ^(n-1), and ${\mathbf{d}}^{\left(n\right)}=-{\mathbf{g}}^{\left(n-1\right)}+{\beta }^{\left(n\right)}{\mathbf{d}}^{\left(n-1\right)},{\beta }^{\left(n\right)}=\frac{\parallel {\mathbf{g}}^{\left(n-1\right)}{\parallel }^2}{\parallel {\mathbf{g}}^{\left(n-1\right)}{\parallel }^2}.$

In our study, positivity or negativity constraints are separately imposed on the reconstruction results ${\delta \stackrel{\mathit{~}}{\mathit{\mu }}}_a^{\left(n\right)}$ and δD̃ ⁽ⁿ⁾ after each iteration. In some instances we assume prior knowledge of the direction of the perturbation and apply the appropriate constraint. In the more general case, we recognize that more than one approach could be adopted to implement solution constraints. If one assumes no prior knowledge, then arbitrarily imposing a constraint could result in a grossly incorrect solution. We avoid this by adopting a two–step process; first imposing one set of constraints (i.e., positivity on one coefficient and negativity on the other, or positivity or negativity on both coefficients), then reversing these. The sum of the two partial solutions thus obtained results in the vectors δµ̃_a and δD̃.

After the intermediate solutions (δµ̃_a and δD̃) are computed, they are rescaled to get the final results, via the formulas

For all reconstruction results presented below, solutions were limited to a first-order computation involving at most 1000 CGD iterations. The finite-element grid used comprised 1296 elements. The optical properties of the reference medium used for the initial guess were µ_a =0.04 cm^-1 and µ′_s =10 cm^-1 (D≈0.0332 cm) for the numerical studies, and µ_a =0.02 cm^-1 and µ′_s =10 cm^-1 (D≈0.0333 cm) for the experimental studies.

3. Simulation results

Figure 1 shows a sketch of the target medium considered for the simulation studies. As many as four objects, symmetrically positioned and having the dimensions and locations shown, were included. Each tomographic simulation considered 6 sources and 18 detectors, providing 108 source–detector pairs. Computed were the detector responses for the target medium and those corresponding to a selected reference medium having the same external geometry, size, and source–detector configuration. Simulated data from these calculations were subsequently analyzed using the methods described above.

 

Figure 1. Target geometry and source–detector configuration for simulation cases.

Download Full Size | PDF

Figure 2 shows the original and reconstructed profiles for the target medium considered. Rows one and three are the original profiles for δµ_a and δD, respectively; rows two and four are the corresponding reconstructed profiles. The color scale indicates the magnitude and sign of the perturbation. Proceeding from left to right, in the first column we find two objects are present. In the left-hand inclusion, it is only the µ_a value that is increased (δµ_a =0.04 cm^-1) relative to the background; in the right-hand inclusion, only the D value is increased (δD≈0.033 cm). In the second column, we consider four objects. Here only one of the two coefficients is perturbed in any given object. For the left- and right-hand inclusions, the µ_a value is increased (δµ_a =0.04 cm^-1); for the top and bottom inclusions, the D value is decreased (δD≈-0.017 cm). In columns 3–7 we consider the more general case where perturbations in both µ_a and D can occur simultaneously in any object. In column three, two inclusions are present, and the µ_a and D values of both are decreased (δµ_a =-0.02 cm^-1) and increased (δD≈0.033 cm), respectively. In column four, three objects are present. In the left-hand and top inclusions the µ_a value is increased (δµ_a =0.04 cm^-1), while the D value is decreased in the right-hand and top inclusions (δD≈-0.017 cm). Thus, it is only the top inclusion that has perturbations in both coefficients, while the left- and right-hand inclusions have perturbations in only µ_a and only D, respectively. In the case of columns five and six, the locations of the objects and the pattern of their coefficient value perturbations are the same as those in columns three and four, respectively. What differs is the direction of the perturbation for one of the coefficients. For column five, the µ_a of the inclusions was increased by 0.04 cm^-1 relative to the background, while δD is the same as that in column three. For column six, only the D value for the top and right-hand inclusions differs from that in column four. Here, the δD value was ~0.033 cm, while the µ_a value for the right-hand and top inclusions remains unchanged. It is worth noting that these latter two cases are similar to the more difficult examples explored by Arridge and Lionheart [2], in which the influence of an increase or decrease in an object’s absorption coefficient on surface measurements can be offset by a decrease or increase, respectively, in the value of its scattering coefficient (i.e., increase or decrease in D). These are conditions where evidence for solution nonuniqueness using intensity-only data was presented. Finally, in column seven, the µ_a and D values for the inclusions were decreased and increased, respectively, in the left-hand object (δµ_a =-0.02 cm^-1, δD≈0.033 cm), while the opposite trend occurs in the right-hand object (δµ_a =0.04 cm^-1, δD≈-0.017 cm).

 

Figure 2. Original and reconstructed profiles for the target medium considered (7 types). Rows one and two are the original and reconstructed profiles for δµ_a , respectively; rows three and four are the corresponding original and reconstructed D profiles, respectively. Color scale indicates amplitude of perturbation.

Download Full Size | PDF

Inspection of the reconstructed profiles shows that in each case we can effectively isolate perturbations in the absorption and diffusion coefficients, whether or not they are co-located. In the first six cases this was accomplished by assuming prior knowledge of the sign of the perturbation, and applying the appropriate constraint as described in Section 2. In practice, this could correspond to situations where the influence of a particular manipulation of tissue (e.g., inflation of a pressure cuff) would impose an expected response (e.g., venous congestion and hence an increase in absorption by hemoglobin). Note that in columns two and four, a small degree of image contrast is seen in the D image in the locations where only perturbations in µ_a occur. We do not believe this is completely attributable to crosstalk, as a positive perturbation in µ_a will be expected to produce a small negative perturbation in D [9], which is observed. A corresponding result is not seen in columns one and six because of the positivity constraint applied to the D reconstruction.

The more general case, shown in column seven, is that in which the perturbation in either coefficient could be positive or negative. To capture this information we applied a two–step process. First, we imposed a positivity constraint on one coefficient and a negativity constraint on the other. Second, we reversed the directions of these constraints, and summed the two solutions together. The resulting images clearly show that we can achieve parameter isolation with a high degree of fidelity and spatial accuracy. In fact, the spatial localization of the reconstruction is in most cases excellent. This finding is even more notable when it is considered that the results presented are limited to a first-order solution. As commented below, we believe this originates from the considerably different structure of the Jacobian operator that results from matrix scaling.

4. Experimental results

To validate the normalized-constraint method, we performed a series of laboratory studies on vessels containing one or more objects whose optical properties differ from a homogeneous background primarily in either its scattering or scattering and absorption coefficients. In addition to acting as optical perturbations, we also introduced dynamic behavior by translating the objects along circular arcs as indicated (see Figure 3), while all the time acquiring fast tomographic data sets using the DC imager recently described [10].

Sketches of the target media explored are shown in Figure 3. In case one, we introduced three rods (6 mm dia.) composed of white Delrin into a hollow cylindrical vessel, also composed of white Delrin, having an outer diameter of 7.6 cm and filled with 1% (v/v) Intralipid. It has been reported that the optical properties of white Delrin are µ_a ≈0.02 cm^-1, µ′_s ≈12 cm^-1 [11], while those of 1% Intralipid are µ_a ≈0.02 cm^-1, µ′_s ≈10 cm^-1 [12]. Thus, it is only the scattering coefficient of the rods that is increased (decreased D) relative to the background. While we did not independently verify these values, we did observe that the detected light intensity was increased on the same side of the vessel as the source, and decreased on the side opposite the source, in the presence of the rods, a finding consistent with the reported coefficient values.

In case two, a similar experiment was performed, except that one of the Delrin rods was replaced by a black glossy metal rod having a similar diameter. We treated this substitution as introducing perturbations in both the absorption and scattering coefficients. In cases three and four we performed similar studies, except here the complexity of the medium was increased via the introduction of a 3-mm-thick clear plastic layer that served as an optical void. These cases were intended to determine whether our methods could correctly isolate optical perturbations (and dynamic behavior) under conditions where the assumptions of the diffusion equation are strongly violated. Data acquisition involved 16 equally spaced source positions with 16 co-located detectors (256 source–detector pairs). The illumination wavelength was 785 nm (10 mW), and full tomographic scans were acquired at 3 Hz for approximately 25 seconds (75 scans).

 

Figure 3. Experimental phantom cases 1, 2, 3, and 4.

Download Full Size | PDF

Figure 4 shows typical reconstructed profiles of the diffusion (top row) and absorption (bottom row) coefficients for the four cases obtained using the two–step constraint protocol. Column one (left side) shows the reconstructed profiles corresponding to case one. Inspection shows nearly complete separation of the absorption and diffusion coefficients. The absorption image reveals that artifact is restricted to the region near the boundary, while the interior is almost completely featureless. The reconstructed diffusion image shows, with high contrast, the correct locations of the inclusions. Quantitative comparison of the image data reveals that the amplitudes of the computed µ_a perturbations are quite small in absolute terms, and also small compared to the magnitude of the computed perturbation in D, revealing minimal crosstalk.

 

Figure 4. The reconstructed diffusion (top row) and absorption (bottom row) profiles for (left to right) cases one through four. Click on figure with mouse to see movie (<1.2 Mb for each). [Media 1] [Media 2] [Media 3] [Media 4] [Media 5] [Media 6] [Media 7] [Media 8]

Download Full Size | PDF

Another feature we have observed, particularly in solutions obtained using Eq. (2), is that, whereas the recovered δµ_a and δD values scale with the selected reference medium used to compute u _r and W _r , their ratio does not (results not shown). That is, we find that the relative magnitude of crosstalk between the absorption and diffusion images is essentially independent of the optical coefficient values chosen for the reference medium used as the initial guess. By clicking the mouse on the image, a movie of the dynamic behavior can be seen. Note that the periodic pauses seen are real and correspond to the time needed for the operator to adjust his grip to allow for manual rotation of the support structure holding the inclusions.

Column two shows that results of similar quality to those in case one are obtained when both coefficients are perturbed. Note that the inclusion contrast for this case is qualitatively similar to that associated with column four of Figure 2, where one of the three inclusions has contrast in both coefficients. The presence of the plastic rods is revealed only in the reconstructed diffusion images, whereas the glossy black metal rod appears in the both the recovered diffusion and absorption images. As in case one, nearly complete isolation of absorption and diffusion coefficients is obtained for the included plastic rods, with only minimal distortion of object shape and location. We also note with some interest that, in particular for cases 1 and 2, object contrast and resolution actually were improved when the µ′_s value used for the initial guess was lowered from 10 cm^-1 to 3 cm^-1.

Data in columns three and four show results obtained in the presence of a circular optical void (cases three and four). This geometry was considered as a crude representation of the layered structure that occurs in the head. Here, the void is meant to represent the layer of clear cerebrospinal fluid that lies in the subarachnoid space. We have considered this case because it has been reported by Dehghani et al. that reconstruction methods based on the diffusion approximation perform poorly under conditions similar to those examined here [13]. In the case of a single inclusion (white plastic rod, column three), we observe that the presence of the optical void does not prevent nearly complete separation of the coefficient profiles. Data in column four show that qualitatively similar results are obtained using three plastic rods, in terms of minimizing crosstalk, although increased artifacts are present. Note that the images of the rods in this case are more centrally located than those in case one, and correspond closely to their actual spacing. These findings show that, using the methods described here, we do not encounter the difficulties that Dehghani et al. have reported using diffusion-based imaging codes [13]. Finally, it deserves mentioning that we have repeated the above analyses, for both simulated and experimental data, using a wide range of reference media as the initial guess [4–6]. In all cases, the image quality obtained did not differ substantially from that presented here.

Dependence of image quality on use of constraints and matrix scaling methods.

Results shown thus far have corresponded to findings obtained when the CDG method was extended by introducing two modifications. The first is to implement a range constraint; the second is to solve Eq. (6) using a scaled weight matrix. It is instructive to examine the influence that the additional operations separately have on image quality and parameter isolation. Results in Figure 5 show movies of dynamic behavior of the recovered coefficients when CGD only (column one), CGD with range constraints (column two), CGD with matrix scaling (column three) and CGD with range constraints and matrix scaling (column four) methods were used to reconstruct the images. The top row corresponds to the reconstructed µ_a , the bottom row to reconstructed D. The experimental data analyzed are the same as those that produced the results shown in column one of Figure 4 (i.e., three plastic rods). Recall that the input data vectors are quantities normalized to the temporal mean value for each respective source–detector channel. Results in column one (CGD only) show that the three-rod structure is clearly resolved in the µ_a map. While the spatial locations and resolution of the inclusions are notable, the result obtained is almost completely attributable to inter-parameter crosstalk. Inspection of the D map in column one shows that the three-rod structure is also recovered with the correct sign of the perturbation (δD<0), but the image quality is impaired by artifact. Results in column two show that introduction of a range constraint (positivity on µ_a , negativity on D) produces a result dominated by artifact on the boundary, with significant loss of internal structure. Results in column three show that use of a scaled weight matrix alone produces a rather unsatisfactory answer. The D image map is heavily dominated by artifact, and crosstalk is evident in the µ_a map. Results in column four are a reproduction of those shown in column one of Figure 4, and reveal excellent parameter isolation and spatial resolution. Thus we find that it is only in combination that the desired result is obtained.

 

Figure 5. Reconstructed diffusion (top row) and absorption (bottom row) profiles for experimental case one using standard CGD method only (column one), CGD method with range constraints (column two), CGD method with weight matrix scaling (column three) and CGD method with range constraints and matrix scaling (column four). Click on figure with mouse to see movie (<1.2 Mb for each). [Media 9] [Media 10] [Media 11] [Media 12] [Media 13] [Media 14] [Media 15] [Media 16]

Download Full Size | PDF

Influence of Scaling on the Structure of the Weight Matrix

It is apparent that the combined methods provide for a satisfactory answer. In an effort to further understand the influence of matrix scaling, we have compared the structure of scaled and unscaled weight matrices. These are shown as movies in Figure 6, for both optical coefficients, as functions of different source–detector configurations. Panels A and C correspond to unscaled weight matrices for µ_a and D, respectively. Panels B and D are the corresponding scaled weight matrices. Inspection of the unscaled matrices reveals the well-known “banana” structure, with the region of greatest weight occurring in the vicinity of the source and detector. The structure of the scaled matrices, on the other hand, is notably different. Here we see that region of highest normalized weight occurs in the interior. This demonstrates that matrix scaling effectively redistributes regions of importance. The net effect for all source–detector pairs is to produce a more uniform weighting of the contributions of the various regions of the target, rather than the bias towards the surface for unscaled weight matrices. We believe it is this effect that principally explains why: 1) we achieve improved spatial resolution, and 2) together with constraints, we achieve improved parameter isolation.

 

Figure 6. Weight functions corresponding to source–detector pairs with source 1 and detectors 1 to 16 for absorption (A and B) and diffusion (C and D) coefficients, before (A and C) and after (B and D) applying matrix scaling. Click on figure with mouse to see movie (<0.5 Mb for each). [Media 17] [Media 18]

Download Full Size | PDF

5. Discussion and Conclusions

In this report, we have explored the problem of inter-parameter crosstalk involving simultaneous reconstruction of images of spatial distributions of the absorption and diffusion coefficients of highly scattering media. By imposing range constraints and scaling the weight matrix, we can effectively isolate perturbations in one coefficient from those in the other.

This capability, verified by experiment, is in striking contrast to the findings that might be expected on the basis of theoretical studies by Arridge and Lionheart [2], who demonstrated an underlying nonuniqueness between the optical properties of a scattering medium and the surface intensities. Specifically, it was shown that a “null space” state can exist, in which the effects of pairs of absorption and scattering coefficient distributions on the surface intensity effectively cancel each other. The inference drawn in Refs. [2,3] and by reviewers of our recent work who point to it is that, given a set of detector responses, essentially any spatial map can be derived because of the “null space” phenomenon. From this it would seem to logically follow that DC imaging methods should not be capable of providing useful information. As shown by results presented here and in several recent reports from our group [10,14–16] and others [17,18], this is not the collective experience.

Before commenting on these seemingly contradictory findings, we do wish to point out one observation pertaining to the use of constraints. Whereas parameter isolation is best achieved when the direction of the constraint is known, as noted, use of an inappropriate constraint can lead to a grossly incorrect solution. We mention this because in instances where we apply the two–step constraint, it is the case that only half of the applied constraints can be correct. In these circumstances, we observe that the two–step process can lead to increased artifact levels in the extreme boundary (cf. column seven, Fig. 2). In fact, it has been our experience that these artifacts are seen almost exclusively in the region of the extrapolated boundary. As this extends beyond the physical boundary, this region has been excluded from the image results shown for the experimental studies.

Influence of Nonuniqueness on DC imaging studies.

Here we attempt to reconcile our results with the reported theoretical findings [2]. We do so, mainly to bring clarity to those who have objected to our numerical and experimental findings. As a matter of established methodological principle, however, we stress that empirical facts have the right-of-way; if a theoretical derivation yields a conclusion that is at odds with experimental results, the reconciliatory burden falls on the theorist, not on the experimentalist. A general explanation we offer is based on a consideration of the dimensionality of infinite sets. For instance, whereas the number of points in 3D space and a line are both infinite, the likelihood that a randomly chosen point in space would fall on a particular pre-selected line is vanishingly small. Likewise, it is our view that whereas there may be an infinite number of absorption-scattering pairs that satisfy the null space criteria for any given set of detector values obtained at the surface, the subset that yield these surface measurements and also are physically (i.e., both µ_a and D must be everywhere positive) and biologically (i.e., values of both µ_a and D must everywhere lie within the ranges that one can expect to find in tissue) possible may have a dimensionality that is small relative to that of the full null space. We would further suggest that many of the absorption-scattering distributions that do lie within the plausible subset are qualitatively similar to each other, and to the physical medium from which the measurements were, in fact, taken. If these hypotheses are correct for a given target medium, then it is possible to recover accurate, or at least reasonable, estimates of that medium’s optical properties from intensity-only measurements.

This consideration brings into view the comments made by Hoenders [19] who had earlier examined this issue of uniqueness and noted that the real-world solution of inverse problems depends on many more things than just a nonuniqueness property. One factor we have considered, and one emphasized in Ref. [19] is the use of prior knowledge. As we have demonstrated here, this can be applied in ways that seemingly retain general applicability. No doubt other sensible strategies can likewise be adopted to provide useful measures using DC methods.

Of greater importance for real-world problems is a clear statement of precisely what information one seeks to infer from measurement of the light-field distribution. Thus it may turn out, paraphrasing Hoenders [19], that the property of nonuniquness has hardly any influence on one’s ability to extract the particular information sought. One class of problem for which we believe this almost certainly holds concerns the characterization of dynamic states [14]. In fact, in a recent report [4] we demonstrated that the accuracy with which temporal variations in the optical properties of inclusions can be localized and correctly characterized is much greater than are those associated with the absolute coefficient values, given the expected constraints of any real experiment.

As a corollary to this, we call attention to two findings in the present report wherein the results obtained, although not completely correct, nevertheless provide highly useful information. Specifically, we point to the absorption coefficient images shown in column one of Figure 5. As indicated, while the presence of inclusion contrast in these images is notably incorrect (evidence of crosstalk), the image quality, in terms of the identified location of the inclusions, artifact levels and the observed dynamic behavior, is nonetheless striking. A similar assignment applies to the results presented in Figure 4, for experimental measures obtained in the presence of an optical void. Notably absent from the recovered images is evidence of the void itself. Surely the absence of this information in no way invalidates our ability to distinguish variations in absorption from those in scattering, to determine object location or to observe dynamic behavior. The flip side to this is that the methodology described herein clearly is not well suited to characterize these properties and those of the void. In this instance it may be that fundamental limitations stemming from properties of uniqueness, modeling error or other factors will render inverse solutions obtained using DC illumination techniques unsuitable. That such limitations may exist, however, in no way supports the suggestion made in Refs. [2,3], that DC illumination methods are somehow inherently unsuitable for characterizing the optical properties of highly scattering media.

One other point deserves emphasis regarding the practical utility of methods for characterizing the optical properties of highly scattering media. For the current report the diffusion approximation to the radiation transport equation was used when computing the weight matrices and reference detector readings that enter the system equation (Eq. (6)). It is widely appreciated that this formulation is not well suited for predicting light distributions in media containing voids [20]. As noted above, however, using this approach we are nevertheless fully capable of separating absorption from scattering, defining object location and observing dynamic behavior, even in the presence of a 3-mm-thick clear circular layer. These findings stand in marked contrast to a recent report by Dehghani et al. [13] who performed a largely similar experiment involving numerical simulations on a static medium. They found that in cases where the thickness of the circular void approaches 1 mm, resolution of internal structure is completely lost. Common to both studies is a pronounced mismatch between the conditions that propagating photons experience in the medium and those that are assumed to exist in the model used for image recovery. In the computations presented in Ref. 13, this was modeled by considering a two–dimensional medium containing a non-scattering ring wherein the inverse problem was solved using the diffusion approximation, while the forward problem was solved using the transport equation. We believe the apparent discrepancy can be attributed mainly to the different types of detector data considered. In our case we examine a differential measure, while Dehghani et al. do not. This consideration is consistent with a recent report wherein we demonstrated that the quality of information derived from examination of light-field distributions can vary greatly depending on the specifics of the how the data vectors are treated [4]. It would seem that the practical limits on inferring useful information from highly scattering media are not as restrictive as some have perceived.