## Abstract

In this report, we present a method for reducing the inter–coefficient crosstalk problem in optical tomography. The method described is an extension of a previously reported normalized difference method that evaluates relative detector values, and employs a weight matrix scaling technique together with a constrained CGD method for image reconstruction. Results from numerical and experimental studies using DC measurement data demonstrate that the approach can effectively isolate absorption and scattering heterogeneities, even for complex combinations of perturbations in optical properties. The significance of these results in light of recent theoretical findings is discussed.

©2001 Optical Society of America

## 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**
_{2} and **u**
_{1}
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}\xb7\delta {\mu}_{a}$ and $\delta \tilde{\mathbf{D}}={\left[{\mathbf{R}}^{\left(D\right)}\right]}^{-1}\xb7\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)}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}{\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̃**
^{(0)}, compute
**g**
^{(0)}=**A**·δ**x̃**
^{(0)}-**b**;**d**
^{(1)}=-**g**
^{(0)}
; *β*
^{(1)}=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}\xb7{\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̃**
^{(n)} 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 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).

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 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.

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.

*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.

## 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.

## 6. Acknowledgement

This research was support in part by NIH grant no. R01 CA66184. The authors wish to thank C. H. Schmitz for assistance in collecting the experimental data.

## References and links

**1. **T.O. McBride, B. W. Pogue, U.L. Österberg, and K.D. Paulsen, “Separation of absorption and scattering heterogeneities in NIR tomographic imaging of tissue,” in *Biomedical Topical Meetings*, OSA Technical Digest (Optical Society of America, Washington, D.C., 2000), pp. 339–341.

**2. **S. R. Arridge and W. R. B. Lionheart, “Nonuniqueness in diffusion-based optical tomography,” Opt. Lett. **23**, 882–884 (1998). [CrossRef]

**3. **J. C. Hebden, S. R. Arridge, and M. Schweiger, “Investigation of alterative data types for time-resolved optical tomography,” *Trends in Optics and Photonics vol. 21, Advances in Optical Imaging and Photon Migration*, James G. Fujimoto and Michael S. Patterson, eds. (Optical Society of America, Washington, DC1998), pp 162–167.

**4. **Y. Pei, H. L. Graber, and R L. Barbour, “Influence of systematic errors in reference states on image quality and on stability of derived information for DC optical imaging,” Appl. Opt., 2001, in press. [CrossRef]

**5. **Y. Pei, *Optical Tomographic Imaging Using the Finite Element Method*, Ph. D. Thesis (1999), Polytechnic University.

**6. **Y. Pei, F.-B. Lin, and R. L. Barbour, “Model-based imaging of scattering media using relative detector values,” presented at 1999 *OSA Annual Meeting & Exhibit: Optics in High-Tech Industries* (Santa Clara, CA, September 26–30).

**7. **J. Chang, H.L. Graber, R.L. Barbour, and R. Aronson, “Recovery of optical cross-section perturbations in dense-scattering media by transport-theory-based imaging operators and steady-state simulated data,” Appl. Opt. **35**, 3963–3978 (1996) [CrossRef] [PubMed]

**8. **P.E. Gill, W. Murray, and M. H. Wright, *Practical Optimization*, Academic Press, New York (1981).

**9. **R. Aronson and N. Corngold, “Photon diffusion coefficient in an absorbing medium,” J. Opt. Soc. Am. A **14**, 262–266 (1997).

**10. **C. H. Schmitz, M. Löcker, J. Lasker, A. H. Hielscher, and R. L. Barbour, “Performance characteristics of silicon photodiode (SiPD) based instrument for fast function optical tomography,” Proc. SPIE4250, San Jose, CA, *in press*, (2001).

**11. **S. Fantini, M.A. Franceschini, G. Gaida, H. Jess, H. Erdl, W.W. Mantulin, E. Gratton, K.T. Moesta, P.M. Schlag, and M. Kaschke, “Contrast enhancement by edge effect corrections in frequency-domain optical mammography,” in *OSA Trends in Optics and Photonics on Advances in Optical Imaging and Photon Migration* , R.R. Alfano and J.G. Fujimoto, eds., (Optical Society of America, Washington DC, 1996, **Vol. 2**, pp. 160–163.

**12. **S. T. Flock, S. L. Jacques, B. C. Wilson, W. M. Star, and M. J. C. van Gemert, “Optical Properties of Intralipid: A phantom medium for light propagation studies,” Lasers Surg. Med. , **Vol. 12**, 510–519 (1992). [CrossRef]

**13. **H. Dehghani, D. T. Delpy, and S. R. Arridge, “Photon migration in non-scattering tissue and the effects on image reconstruction,” Phys. Med. Biol. **44**, 2897–2906 (1999). [CrossRef]

**14. **R. L. Barbour, H.L. Graber, Y. Pei, S. Zhong, and C.H. Schmitz, “Optical tomographic imaging of dynamic features of dense scattering media,” J. Opt. Soc. Am. A, in press (2001). [CrossRef]

**15. **H. L. Graber, Y. Pei, and R. L. Barbour, “Imaging of spatiotemporal coincident states by dynamic optical tomography,” Proc. SPIE4250, San Jose, CA, in press, (2001).

**16. **G. Landis, S. Blattman, T. Panetta, C. H. Schmitz, H. L. Graber, Y. Pei, and R. L. Barbour, “Clinical applications of dynamic optical tomography in vascular disease,” Proc. SPIE4250, San Jose, CA, in press, (2001).

**17. **N. Iftimia and H. Jiang, “Quantitative optical image reconstruction of turbid media by use of direct-current measurements,” Appl. Opt. **39**, 5256–5261, (2000). [CrossRef]

**18. **Y. Xu, N. Iftimia, H. Jiang, L. Key, and M. Bolster, “Imaging of in vitro and in vivo bones and joints with continuous-wave diffuse optical tomography,” Opt. Express **8**, 447–451 (2001). [CrossRef] [PubMed]

**19. **B.J. Hoenders, “Existence of invisible nonscattering objects and nonradiating sources,” J. Opt. Soc. Am. A **14**, 262–266, (1997). [CrossRef]

**20. **A.H. Hielscher, R.E. Alcouffe, and R.L. Barbour, ‘Comparison of finite-difference transport and diffusion calculations for photon migration in homogeneous and heterogeneous tissues,” Phys. Med. Biol. **43**, 1285–1302 (1998). [CrossRef] [PubMed]