In this paper we restrict our attention to problems in which the transport physics [
2] associated with the propagation of light at wavelength
λ through tissue can be adequately approximated using a diffusion model [
21,
22] of the form
where Φ(
r,
λ) is the photon fluence rate at position
r due to light of wavelength
λ injected into the medium,
v is the electromagnetic propagation velocity in the medium,
is the spatially varying absorption coefficient, and
S(
r,
λ) is the photon source with units of optical energy per unit time per unit volume. For the work in this paper the sources are considered to be delta sources in space and can be written as
S(
r,
λ) =
S
_{0}(
λ)
δ(
r −
r
_{s}) with
S
_{0}(
λ) the source power. Lastly
D(
λ) is the diffusion coefficient, given by
where
is the reduced scattering coefficient. We assume that the reduced scattering coefficient is spatially constant and known, and we focus solely on reconstructing the chromophore concentrations. Though recent work has demonstrated the possibility of using cw data for the recovery of chromophore concentrations and scattering parameters [
10], for simplicity, here we concentrate exclusively on the problem of recovering concentration information from hyperspectral data. It should be noted that the first assumption of a known reduced scattering coefficient, and effects associated with a wrong choice, can be addressed in practice by a preliminary measurement of scattering properties of tissue. Such measurements of average tissue properties are wellestablished with frequencydomain or timedomain techniques and do not constitute a basic limitation to the proposed imaging approach. The second assumption of a uniform scattering coefficient is also justified by a much larger contrast (both in terms of value and spectral shape) typically provided by absorption versus scattering in a large number of cases (cancer, functional activation, localized hemorrhage, etc.). However, this assumption can also be avoided by estimating space varying scattering properties in addition to chromophore structure. In Section 7 we discuss our ongoing efforts in this area. Finally, we note that
Eq. (1) is often written in the frequency domain with a term
jω/
D(
λ) included in the parentheses on the left hand side where
ω is the modulation frequency of the light intensity [
23]. Here we consider exclusively problems for which
ω = 0 so that the diffusion equation takes the form shown in
Eq. (1).
We consider an infinite medium problem. In medical imaging, measurements are typically made by placing the source and detector on the tissue surface. Treating such system as infinite will obviously lead to some discrepancies between theory and experiment. Considering that wavefronts have the same general shape except near the boundaries and object sensitivity is unaffected by the boundaries we only do reconstruction with no boundaries [
24]. In Section 7 we briefly discuss alterations to our approach needed to deal with finite geometries.
The Born approximation is derived by decomposing
, as the sum of a constant background absorption,
μ_{a}(
λ), and a spatially varying perturbation Δ
μ_{a}(
r,
λ). Then
Eq. (1) becomes
In
Eq. (2) the total fluence rate, Φ(
r,
λ) is written as the sum of the fluence rate Φ
_{i}(
r,
λ) due to the source acting on the background medium and a scattered fluence rate, Φ
_{s}(
r,
λ), due to the inhomogeneities. As explained more fully in Section 4, we assume the data we have for inversion are in fact samples of Φ
_{s}. To obtain a linear relationship between the scattered fluence rate and the chromophore concentrations, we first need a linear mapping between the perturbations to the material properties and Φ
_{s} which is derived by subtracting
Eq. (1) from
Eq. (2) giving
where
and Δ
k
^{2}(
r,
λ) = (
v/D(
λ)) Δ
μ_{a}(
r
, λ). Assuming the availability of a Green’s function,
G(
r,
r′,
λ) for the operator
as is the case for an unbounded medium as well as range of bounded problems [
25], we rewrite
Eq. (3) as
Unfortunately, the total field Φ(
r,
λ) depends implicitly on Δ
k
^{2} via
Eq. (2) thereby resulting in a nonlinear relationship between the scattered field and the absorption perturbation. The Born linearization is achieved under the assumption that the total fluence rate, Φ, appearing in the integrand of
Eq. (4) can be approximated as the incident field, Φ
_{i}, which satisfies
and thus is
not dependent on Δ
k
^{2}. Replacing Φ(
r,
λ) by Φ
_{i}(
r,
λ) and writing Δ
k
^{2} =
v/DΔ
μ_{a} in
Eq. (4) yields the Born model used in this paper which we write as
where
r
_{d} is the location of the detector and (with a small abuse of notation) Φ
_{i}(
r,
r
_{s},
λ) is used here to denote the incident field at position
r and wavelength
λ due to a deltasource located at
r
_{s}.
Equation (4) provides a linear relationship between the scattered fluence rate and the absorption perturbation. To relate the scattered fluence rate to concentrations of chromophores, Δ
μ_{a} is decomposed as follows [
17]
where
N_{s} is the number of absorbing species for the problem under investigation,
ε_{k}(
λ) is the extinction coefficient for the
k^{th} species at wavelength
λ, and
c_{k}(
r) is the concentration of species
k at location
r. To obtain the fully discrete form of the Born model used in Section 3, we expand each
c_{k}(
r)
where
c_{k, j} is the value of the concentration for species
k in
V_{j}, the
j^{th} “voxel”. The
(
r) function is an indicator function where
After using
Eq. (7) and
Eq. (8) in
Eq. (4) and performing some of algebra we obtain
We approximate
Eq. (4) as the value at the center of each pixel multiplied by the area or volume of each pixel or voxel, so in
Eq. (10) we use
a as the area of a pixel. This setup is demonstrated in .
The computational tractability of the inversion scheme we describe in Section 3 arises from the linear algebraic structure associated with
Eq. (10). We start by defining
c
_{k} ^{Nv} as the vector obtained by lexicographically ordering the unknown concentrations associated with the
k^{th} chromophore and Φ
_{s}(
λ) to be the vector of observed scattered fluence rate associated with all sourcedetector pairs collecting data at wavelength
λ. Now, with
N_{λ} the number of wavelengths used in a given experiment,
Eq. (10) is written in matrixvector notation as
It should be noted in
Eq. (11) that the matrix has elements which are also the matrices
K
_{l}. The (
m,
j)
^{th} element of the
K
_{l} matrix is given by the product (
v/D(
λ_{l}))
G(
r
_{m},
r
_{j},
λ_{l})Φ
_{i}(
r
_{j},
r
_{m},
λ_{l}), where
m represents the
m^{th} sourcedetector pair and as before
j represents the
j^{th} voxel. Assuming that for a given experiment
N_{sd} source detector pairs operate at all
N_{λ} wavelengths, then each
K
_{l} has
N_{sd} rows and
N_{v} columns so that the whole matrix
K is of size
N_{sd}N_{λ} ×
N_{v}N_{c}. If, for example, in an experimental setup where
N_{sd} = 57,
N_{c} = 2, and image reconstruction is done for 1800 pixels,
N_{v} = 1800, and
N_{λ} = 126 results in a
K matrix of size 7182
× 3600.
While for realistically sized problems, it is difficult to store the full K matrix in memory, the processing methods developed in Section 3 require only the result of multiplying K or K
^{T} (the transpose of K) by appropriately sized vectors. Hence, we need only compute and store the N_{λ} matrices K
_{l} as well as the N_{λ} × N_{c} array of extinction coefficients. Then computation of the product Kc can be carried out using the Matlablike pseudocode in Algorithm 1 with a similar approach possible for evaluating K
^{T} Φ_{s}.
Algorithm 1 Matlablike code for computing Kc product

forl = 1 to N_{λ}do 
fork = 1 to N_{c}do 
Φ_{c} = Φ_{c} +ε_{k}(λ_{l})K_{k}; 
end for 
Φ = [Φ;Φ_{c}]; 
end for 