Hyperspectral image reconstruction for diffuse optical tomography

Fridrik Larusson; Sergio Fantini; Eric L. Miller

doi:10.1364/BOE.2.000946

1. Introduction

Over the past 15 years diffuse optical tomography (DOT) has received considerable attention as a functional imaging modality for a number of application areas including breast tumour detection and characterization [1–3] and brain imaging [4, 5]. While initial work with DOT focused on the recovery of space and time varying maps of the optical absorption and scattering properties of tissue [6], by moving to instruments in which the medium is probed with multiple wavelengths of light, recent systems have demonstrated the ability to recover more physiologically-relevant parameters; specifically chromophore concentrations of species including oxygenated and deoxygenated hemoglobin, lipids, and water [7–10]. Although these efforts represent important advances for moving DOT from the lab to the clinic, there still remains a variety of challenges in terms of stably and accurately characterizing the distribution of chromophores.

As is well known, the recovery of chromophore concentrations from DOT data requires the solution of an ill-posed, non-linear inverse problem [3]. Roughly speaking, the physics associated with the interaction of light and tissue coupled with the ability to collect limited quantities of data result in an imaging problem that is sensitive to noise in the data and un-modelled physical processes. This sensitivity manifests itself in a number of ways. The reconstructed images can be corrupted by various artifacts such as large amplitude, high frequency oscillations [7] or negative values for quantities such as concentrations that are required to be positive valued [11]. When reconstructing multiple chromophores, cross-talk artifacts can contaminate the various images. Cross-talk arises when there are e.g., two chromophores distributed in spatially disjoint areas, but the reconstruction shows evidence of both species in both regions.

The challenge of ill-posedness is typically addressed in two ways. By posing the image formation problem in a variational context, various regularization techniques [12] can be used to “discourage” the presence of artifacts in the reconstructed images [13–15] and enforce positivity of the recovered concentrations; however regularization does not explicitly address the cross-talk problem. Additionally, much work has been devoted to increasing the quantity of data available for processing by encircling the region to be imaged with sources and detectors [16] and by collecting data at multiple wavelengths and modulation frequencies [8, 9].

Augmenting the data set however is not without its difficulties. In some breast and most brain imaging applications, it is infeasible to place optical sources and detectors around the region of interest and one is constrained to backscatter and/or transmission type geometries [6, 17]. The collection of additional data exacerbates a second fundamental difficulty associated with DOT imaging: computational complexity. The nonlinear relationship between the observed data and the chromophore concentrations requires the use of iterative methods for solving the reconstruction problem. These approaches are based on the underlying physical model of light propagation through tissue. Because these models must be evaluated many times during each iteration, as more types of data are acquired (multi-wavelength or multi-modulation frequency), the computational burden of the algorithms quickly becomes a significant bottleneck. Largely for this reason most multi-spectral systems collect data at fewer than ten wavelengths [8–10].

In this paper, we demonstrate the utility of hyperspectral diffuse optical tomography (Hy-DOT) in which data from a great number of wavelengths (here up to 120) are employed to recover concentration images of multiple chromophores. Of specific interest are problems in which limited data are available for processing such as breast imaging where the breast is placed between compression plates [17] or brain imaging where transmission data cannot even be acquired. For concreteness, the examples in Section 6 focus on the breast imaging problem. In a typical DOT system three measurement schemes are used for measurements: time domain, frequency domain, and continuous wave (cw). The cw method, used in this work, is the simplest, least expensive, and provides the fastest data collection [17].

To address the computational burden associated with the need to model the physics on a wavelength-by-wavelength basis, we employ the Born approximation. Although the limitations of the Born model are well documented [18, 19] for the purposes of establishing the utility of HyDOT, the Born model provides us with a convenient place to start. Indeed extending the effort to consider the full physical model is relatively straightforward in theory though quite time consuming in practice due mostly to the need to consider implementation challenges unrelated to the fundamental issues of interest here; namely establishing the utility of a hyperspectral data set for DOT. Moreover, the results in this paper indicate that such effort may not even be necessary. Using experimental data, we demonstrate that the availability of hyperspectral data in conjunction with a Born model can yield reconstructed chromophore concentrations that are quantitatively quite accurate and with significantly reduced cross-talk relative to images obtained using even well-chosen sets of roughly ten wavelengths.

In a bit more detail, the approach considered here is based on a variational formulation of the image formation problem. The associate optimization functional is comprised of two terms: one penalizing data misfit in which the Born model is embedded and second Tikhonov-type smoothness regularization term. Because we regularize the smoothness in the reconstruction of each chromophore independently, the weight of the penalty associated with each recovered image must be determined. Here we describe a multi-parameter L-curve scheme for addressing this issue [20]. To avoid constructing negative-valued concentrations, a positivity constraint is enforced. The resulting constrained cost function is known to possess a unique minimum which should in theory make the reconstruction results independent of how we initialize the optimization algorithm. In practice however, the results do depend on initialization. Indeed, though there may not be local minima for the cost function, the ill-posedness of the problem appears to yield an objective function that is rather flat in places thereby impeding rapid convergence of the iterates. Hence, we also provide in this paper an approach to initialization that yields the results mentioned above.

The remainder of this paper is structured as follows. In Section 2 we discuss the forward problem for DOT and provide the derivation of the Born model used as the basis for our inversion. The details of the inverse processing are provided in Section 3 where we specifically concentrate on the variational formulation, our strategy for solving the associated optimization problem. In Sections 4 and 5 detail methods used to judge results and for selecting regularization parameters. Section 6 is devoted to presenting simulated and experimental validation demonstrating image reconstruction showing improved accuracy and reduced cross-talk when using hyperspectral information. Finally, conclusions and future efforts are provided in Section 7.

2. Forward problem

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

(\nabla^{2} + \frac{ν μ_{a}^{0} (r, λ)}{D (λ)}) Φ (r, λ) = \frac{- v}{D (λ)} S (r, λ)

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,

μ_{a}^{0} (r, λ)

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 ₀(λ)δ(r − r _s) with S ₀(λ) the source power. Lastly D(λ) is the diffusion coefficient, given by

D (λ) = v / (3 {μ^{'}}_{s} (λ))

where

{μ^{'}}_{s}

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 well-established with frequency-domain or time-domain 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 $μ_{a}^{0} (r, λ)$ , as the sum of a constant background absorption, μ_a(λ), and a spatially varying perturbation Δμ_a(r, λ). Then Eq. (1) becomes

(\nabla^{2} - \frac{v μ_{a} (λ)}{D (λ)} - \frac{v Δ μ_{a} (r, λ)}{D (λ)}) (Φ_{i} (r, λ) + Φ_{s} (r, λ)) = \frac{- v}{D (λ)} S (r, λ) .

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

[\nabla^{2} + k_{0}^{2} (λ)] Φ_{s} (r, λ) = - Δ k^{2} (r, λ) Φ (r, λ)

where

k_{0}^{2} (λ) = - v μ_{a} (λ) / D (λ)

and Δk ²(r, λ) = (v/D(λ)) Δμ_a(r , λ). Assuming the availability of a Green’s function, G(r,r′, λ) for the operator

\nabla^{2} + k_{0}^{2} (λ)

as is the case for an unbounded medium as well as range of bounded problems [25], we rewrite Eq. (3) as

Φ_{s} (r, λ) = \int G (r, r^{'}, λ) Φ (r^{'}, λ) Δ k^{2} (r^{'}, λ) d r^{'} .

Unfortunately, the total field Φ(r, λ) depends implicitly on Δk ² 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

(\nabla^{2} - \frac{v μ_{a} (λ)}{D (λ)}) Φ_{i} (r, λ) = \frac{- v}{D (λ)} S (r, λ)

and thus is not dependent on Δk ². Replacing Φ(r, λ) by Φ_i(r, λ) and writing Δk ² = v/DΔμ_a in Eq. (4) yields the Born model used in this paper which we write as

Φ_{s} (r_{d}, λ) \approx \frac{v}{D (λ)} \int G (r_{d}, r^{'}, λ) Φ_{i} (r^{'}, r_{s}, λ) Δ μ_{a} (r^{'}, λ) d r^{'} .

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 delta-source 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]

Δ μ_{a} (r, λ) = \sum_{k = 1}^{N_{s}} ɛ_{k} (λ) c_{k} (r) .

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)

c_{k} (r) = \sum_{j = 1}^{N_{v}} c_{k, j} φ (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

φ (r) = {\begin{array}{l} 1, & if r \in V_{j} \\ 0, & if r \notin V_{j} . \end{array}

After using Eq. (7) and Eq. (8) in Eq. (4) and performing some of algebra we obtain

Φ_{s} (r_{d}, λ) = a \sum_{k = 1}^{N_{c}} \sum_{j = 1}^{N_{v}} \frac{v}{D (λ)} G (r_{d}, r_{j}, λ) Φ_{i} (r_{j}, r_{s}, λ) ɛ_{i} (λ) c_{k, j} .

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 Fig. 1(a).

Fig. 1 (a) The setup of sources and detectors for simulation reconstructions. Same orientation of axes is used for experimental data. (b) Molar extinction coefficients used in simulations plotted as a function of wavelength.

Download Full Size | PDF

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 ∈ ℝ^N_v 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 source-detector pairs collecting data at wavelength λ. Now, with N_λ the number of wavelengths used in a given experiment, Eq. (10) is written in matrix-vector notation as

[\begin{matrix} Φ_{s} (λ_{1}) \\ Φ_{s} (λ_{2}) \\ ⋮ \\ Φ_{s} (λ_{N_{λ}}) \end{matrix}] = [\begin{matrix} ɛ_{1} (λ_{1}) K_{1} & ɛ_{2} (λ_{1}) K_{1} & \dots & ɛ_{N_{c}} (λ_{1}) K_{1} \\ ɛ_{1} (λ_{2}) K_{2} & ɛ_{2} (λ_{2}) K_{2} & \dots & ɛ_{N_{c}} (λ_{2}) K_{2} \\ ⋮ & ⋮ & ⋮ & ⋮ \\ ɛ_{1} (λ_{N_{λ}}) K_{N_{λ}} & ɛ_{2} (λ_{N_{λ}}) K_{N_{λ}} & \dots & ɛ_{N_{c}} (λ_{N_{λ}}) K_{N_{λ}} \end{matrix}] [\begin{matrix} c_{1} \\ c_{2} \\ ⋮ \\ c_{N_{c}} \end{matrix}] \Leftrightarrow Φ_{s} = Kc

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 source-detector 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_sdN_λ × N_vN_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 Matlab-like pseudo-code in Algorithm 1 with a similar approach possible for evaluating K ^T Φ_s.

Algorithm 1 Matlab-like code for computing Kc product

for l = 1 to N_λ do

for k = 1 to N_c do

Φ_c = Φ_c +ε_k(λ_l)K _k;

end for

Φ = [Φ;Φ_c];

end for

3. Inverse problem

The image reconstruction method is cast as the solution to a non-negative least squares (NNLS) optimization problem of the form

\hat{c} = arg min_{c \geq 0} {‖ W (Φ_{s} - Kc) ‖}_{2}^{2} + {‖ Lc ‖}_{2}^{2}

where for any vector x,

{‖ x ‖}_{2}^{2} \equiv x^{T} x

is the squared two-norm of x. The first term in Eq. (12) requires that the reconstructed concentration images yield simulated data that are consistent with the observations Φ. Following [2], the weight matrix W reflects the structure of the noise corrupting the data. While a Poisson model is technically the most appropriate for DOT data, as is frequently done [26] we employ a Gaussian approximation in which independent, zero mean Gaussian noise is assumed to corrupt each datum. The reason for this is that with a sufficiently large number of detected photons, the Poisson statistics can be approximated by a Gaussian distribution [14]. Letting

σ_{m}^{2}

denote the variance of the noise corrupting the m^th elements of Φ, W is constructed as a diagonal matrix with 1/σ_m the m^th element along the diagonal. For the experimental and simulated data the variance is calculated from

σ_{m}^{2} = Ω (m) 10^{- \frac{{SNR}_{m}}{10}} .

where Ω(m) corresponds to the photon count for each source-detector pair. The SNR for each element of Φ, is then calculated from

{SNR}_{m} = 10 {log}_{10} (Ω (m) / \sqrt{Ω (m)}) .

In experimental data

\sqrt{Ω (m)}

is the standard deviation of the Poisson noise distribution.

The second term on the right hand side in Eq. (12) represents the regularization. As discussed in the Introduction, in this work we use a smoothness-type regularizer in which the amount of regularization is allowed to vary for each chromophore. Due to sensitivity of the reconstruction to the regularization parameters the optimal parameter for one chromophore is not necessarily the optimal value for another. Separating the parameters for each chromophore allows the reconstruction to optimize it for each chromophore and easily include many different species of chromophores. Given the structure of the vector c defined in Eq. (11), the regularization matrix takes the form

L = 𝒟 (α) \otimes [\begin{matrix} \nabla_{x} \\ \nabla_{y} \end{matrix}]

where α ^T = [α ₁ α ₂ ...αN_c ] is a vector of Nc regularization parameters, 𝒟(x) is a diagonal matrix with the elements of the vector x on the main diagonal, ∇_x and ∇_y are matrices representing first difference approximations to the gradient operators in the horizontal and vertical directions respectively, and for matrices A and B, A ⊗ B is the Kronecker product [27] of A and B.

The NNLS problem is solved by using the lsqnonlin algorithm in MATLAB. This algorithm uses a trust-region reflective algorithm that employs matrix-vector products instead of having to compute the value of the sum of squares from Eq. (12) [28]. The K matrix is the Jacobian matrix of the measurements used in our reconstruction scheme. This allows the algorithm to interact with the matrix K only through the matrix vector products Kf and K ^T v, for various vectors f and v. For the case of DOT NNLS becomes highly attractive for its computational efficiency when compared to a direct solution of traditional least squares. This is due in part to the fact that computing K ^T K can require large amounts of computational overhead. The number of voxels in a given solution becomes somewhat limited by the necessity of solving the system defined by K ^T K or some regularized version thereof. Because of the design of K when the number of voxels increases, the size of K ^T K and the computation required for elimination both increase much more rapidly than with NNLS [29].

As discussed in the introduction a good initial guess is important to obtain a good results. The approach we use here is as follows. We start by solving Eq. (1) ignoring the positivity constraint with lsqnonlin and using the the method discussed in Section 4 for determining the optimal regularization parameters. Setting all negative values in the c vector to zero then provides the initial guess for the constrained form of the problem. This initialization process allows us to obtain good results from both simulation and experimental data. Like the unconstrained problem the constrained problem is solved with lsqnonlin and optimal regularization parameters are chosen independently in each case.

4. Simulation analysis

Simulations are done to test the effect of hyperspectral information when doing reconstruction of more than one chromophore in a controlled setting. In this paper we consider 2.5D problems in which 3D delta-type sources are used to illuminate the medium but we assume the chromophore concentrations vary only in two dimensions, x and y, and are constant in the third, z. Referring to Eq. (10) and Eq. (11), this implies that the 3D Green’s function and incident field are used to build the K _i matrices but we need only discretize each chromophore image in two spatial dimensions. The reader is referred to [30, 31] for additional details.

For the infinite boundary problem considered in this paper we use the free-space Green’s function which is [2]

G (r, r^{'}, λ) = \frac{- 1}{4 π | r - r^{'} |} e^{j k_{0} (λ) | r - r^{'} |} .

Putting this into Eq. (10) gives the following relation between measurement and absorption perturbation

Φ_{s} (r_{d}, r_{s}, r_{j}, λ) = \sum_{r_{j}} \frac{- e^{j k_{0} (λ) | r_{d} - r_{j} |}}{4 π | r_{d} - r_{j} |} \frac{- e^{j k_{0} (λ) | r_{j} - r_{s} |}}{4 π | r_{j} - r_{s} |} \frac{v}{D (λ)} Δ μ_{a} (r_{j}, λ)

= \frac{v}{16 π^{2} D (λ)} \sum_{r_{j}} \frac{e^{j k_{0} (λ) (| r_{d} - r_{j} | + | r_{j} - r_{s} |)}}{| r_{d} - r_{j} | | r_{j} - r_{s} |} Δ μ_{a} (r_{j}, λ)

Then for each delta function source we calculate the incident field everywhere in the domain using the Green’s function and then the scattered field present at each detector by Eq. (4). It is assumed that

{μ^{'}}_{s}

follows Mie Scattering theory. A scattering prefactor Ψ depends primarily on the number and size of scatterers, and a scattering exponent b depends on the size of scatterers [32]. This is combined as:

{μ^{'}}_{s} = Ψ {(\frac{λ}{λ_{0}})}^{- b} .

The arbitrarily chosen reference wavelength λ ₀ is introduced to achieve a form of the Mie model where Ψ has the units of cm⁻ ¹. In simulations values for Ψ and b are obtained from [33] for the female breast. Values for μ_a is calculated from extinction coefficients, which are in the units cm⁻¹/mM and are obtained from Scott Prahl of the Oregon Medical Laser Center [34]. The concentration in the simulated images are defined in units of molarity or millimol per liter, mM. The extinction coefficients are shown in Fig. 1(b). The background has HbR concentration of 0.01 mM and HbO₂ concentration of 0.01 mM. In each of the simulated images the target concentrations have concentration of 0.02 mM.

The alignment of sources and detectors with respect to the simulated images is displayed in Fig. 1(a). This arrangement of sources and detectors is chosen to represent a common setup in optical mammography where the breast is compressed between two planes containing sources and detectors [17]. The source-detector separations were set to 10 cm as is shown in Fig. 1(a). This is a bit larger than is traditionally used in experimental setups, but the added distance demonstrates the utility of the approach in generating high quality images. In experimental measurements we use a shorter distance of 5 cm, discussed in Section 5. In simulations we reconstruct concentration images of oxygenated and deoxygenated hemoglobin, HbO₂ and HbR respectively. These chromophores are chosen since they mainly cause absorption in breast tissue [35], and breast cancer tumours have been found to have higher HbO₂ and HbR concentrations than normal tissue [36].

Two different sets of images are created to test the reconstruction. The first set is an image with a concentration for HbO₂ and HbR in separate locations. This is comparable to simulations in [35] where separate locations for HbO₂ and HbR are used to test effects of cross-talk. The second set is more complicated with concentrations for HbO₂ and HbR in the same location with different target values. This image is more realistic where chromophore concentrations are usually co-located. These images are shown in Fig. 5 and Fig. 7. Reconstruction is done for these images to explore effects of adding hyperspectral information to the problem, i.e. the improvement in quantitative accuracy and the reduction of crosstalk where a concentration of one chromophore creates a false concentration in an image for another chromophore. The images are created to test the inverse solution with respect to spectral information, regularization, cross-talk and accuracy. To avoid the “inverse crime” the data are generated using a 40 pixel by 40 pixel discretization of the x – y plane while reconstruction is performed on a 20 × 20 grid.

Fig. 5 Reconstruction for first set. Middle images are generated with 6 wavelengths and rightmost images are done with 126 wavelengths. Upper row is for the HbO₂ chromophore and the lower for HbR. Concentration units are in mM.

Download Full Size | PDF

Fig. 7 Reconstruction for second set. Middle images are generated with 6 wavelengths and rightmost images are done with 126 wavelengths. Upper row is for the HbO₂ chromophore and the lower for HbR. Concentration units are in mM.

Download Full Size | PDF

To best understand the utility of a hyperspectral data set, for the simulations we employ the Born model to generate the data. Though this may not be realistic, it allows us to avoid the confounding factor of model mismatch in evaluating HyDOT. Moreover, the shortcoming of this approach are mitigated in Section 6.2 where we consider the processing of laboratory data which, obviously, is not the product of the Born model. In a bit more detail, the data we use for our simulation analysis is computed as

Φ = K [\begin{matrix} c_{1} \\ c_{2} \end{matrix}] + n

where c ₁ and c ₂ are the simulated concentration images for HbO₂ and HbR respectively on the 40 × 40 grid and n represents additive noise. Specifically, as in [2] n is a vector of zero mean, independent Gaussian random variables with variances

σ_{m}^{2}

, defined in Eq. (13), chosen such that a pre-determined signal to noise ratio (SNR) is achieved. This SNR is calculated from Eq. (14)

The reconstructed images are evaluated in three ways: through visual inspection, using mean square error (MSE) as a measure of overall quantitative accuracy for each chromophore and estimation error. For the k^th chromophore the mean square error is computed by using the following equation

{MSE}_{k} = \frac{{‖ c_{k} - {\hat{c}}_{k} ‖}_{2}}{{‖ c_{k} ‖}_{2}}

The estimation error is computed by

e = {‖ c - \hat{c} ‖}_{2}^{2}

This error is calculated when choosing optimal regularization parameters, discussed further in Section 4.1.

4.1. Using gradient matrix and two parameters

The choice of the optimal regularization parameters is done by inspecting the L-hypersurface, which are plotted in Fig. 6 for the concentrations images shown in Fig. 5 [20]. To construct the L-hypersurface we introduce the following quantity

z (α) = {‖ Φ - K \hat{c} (α) ‖}_{2}^{2}

For a single constraint the L-hypersurface reduces to the conventional L-curve which is simply a plot of the residual norm versus the norm of the regularized solution drawn in an appropriate scale for a set of admissible regularization parameters. This allows us to optimize the regularization to compromise between the minimization of these two quantities. For a hypersurface the optimal regularization parameters then should appear where the curvature is greatest in the surface, in other words in the corner of the surface. This corner in the hypersurface which should correspond to a point where the error estimation is minimal. This curvature is computed as a special case of Gaussian curvature [37] from

H = \frac{rt - s^{2}}{w^{4}}

where we have

w^{2} = 1 + p^{2} + q^{2} .

In Eq. (24) each element is a partial derivative of the surface which we write as

p = \frac{\partial z}{\partial α_{1}}, q = \frac{\partial z}{\partial α_{2}}, r = \frac{\partial^{2} z}{\partial α_{1}^{2}}, t = \frac{\partial^{2} z}{\partial α_{2}^{2}}, s = \frac{\partial^{2} z}{\partial α_{1} \partial α_{2}}

Fig. 6 (a) L-hypersurface, defined by (23) plotted against regularization parameters. (b) H curvature, defined by (24), computed for the L-hypersurface. (c) Error estimation surface, defined by (22), plotted against regularization parameters. The optimal regularization parameters are marked in each case with a red arrow.

Download Full Size | PDF

Because we know the ground truth for these simulations, it is possible to determine optimal values (i.e., the one that minimized the mean square error) for α ₁ and α ₂. For a simple chromophore concentrations like in the first set, choosing the regularization parameters is an easy problem. The reason for separating the regularization parameters in this case is that the MSE for HbR reaches a lower value for a slightly different parameter than HbO₂.

The importance of separating the regularization parameters becomes even more evident when regularizing more complex concentration sets. When doing reconstruction for more complicated sets the lowest MSE values for HbO₂ and HbR occur at two very different values. For this set the separation of the regularization parameters is very important. Using only one regularization parameter in this set and more complicated ones, would result in a trade off between reconstructions of chromophores. To reduce that trade off the separation of the chromophores becomes very important. This separation becomes even more important when dealing with data sets with low SNR values such as true measurement data.

5. Experimental analysis

In order to validate the simulation results, physical measurements were performed. The background medium is constructed using milk and water. Milk, with 2% fat, is used due to the similarities of the optical properties to human skin. Similarly, black India ink and blue food dye were used to represent two different chromophores. The ink and dye are mixed into the background of milk and water to achieve μ_a = 0.029 cm⁻¹, at 600 nm, which is in the range of optical absorption of the female breast [38, 39].

The absorption spectra for the ink and dye inclusions have the most significant effect in the 450–700 nm range, shown in Fig. 2. We chose these chromophores because the spectral shapes of their absorption are similar to those of HbO2 and HbR and have been widely used in literature [6, 50]. Therefore the 6 specifically chosen wavelengths were selected from this region where these chromophores had the most significant absorption.

Fig. 2 Absorption spectra of the ink and dye solutions chromophores used in experimental measurements. Specifically chosen wavelengths are marked with an asterisk.

Download Full Size | PDF

In order to obtain multi- and hyperspectral reconstruction values for μ_a and ${μ^{'}}_{s}$ the background has to be known. Phase, amplitude and average intensity data are obtained at two wavelengths using a frequency-domain tissue spectrometer (OxiplexTS, ISS inc., Champaign, IL) to estimate the Ψ and b parameters in Eq. (19) as Ψ = 6.5 cm⁻¹ and b = 0.4. This allows us to have values for ${μ^{'}}_{s}$ at any wavelength [40]. Determining spectrally extrapolated values for the absorption coefficient is harder. Since μ_a does not follow a law like ${μ^{'}}_{s}$ , values are estimated using extinction coefficient data for ink, dye, milk and water. These extinction coefficient are measured in a standard spectrophotometer (Lambda 35, Perkin Elmer Instruments, Shelton, CT).

Two phantom inclusions, named set 1 and set 2 are mixed for different absorption contrasts relative to the background in the range of 3:1 to 1:1. The inclusion in set 1 contains 10% ink and 90% dye and the inclusion for set 2 contains 70% dye and 30% ink. This range of contrast is in the range of traditional tumour contrasts reported in literature, which have been close to 3:1 and lower [41]. Reconstructions are done for 126 wavelengths equally spaced over the whole spectrum and 6 specifically chosen wavelengths as λ = [480, 550, 610, 630, 650, 690] nm. The wavelengths are chosen around the isosbestic point, where the contrast between the chromophores is the highest and where each chromophore has highest absorption. The absorption spectra and the contrast over the spectrum for set 1 and set 2 are shown in Fig. 3 and Fig. 4, respectively.

Fig. 3 (a) Absorption spectra for the background, μ_a, and the inclusion, μ_a + Δμ_a, in experimental set 1, containing 10% ink and 90% dye. (b) Contrast between the background and the inclusion for experimental set 1.

Download Full Size | PDF

Fig. 4 (a) Absorption spectra for the background, μ_a, and the inclusion, μ_a + Δμ_a, in experimental set 2, containing 70% ink and 30% dye. (b) Contrast between the background and the inclusion for experimental set 2.

Download Full Size | PDF

In experimental sets 1 and 2 one cylindrical inclusion containing ink and dye solutions are placed in the background medium. These inclusions are 25 cm long transparent tubes so that optical properties are assumed constant along the z-axis. The light source is an arc lamp (Model No.6258, Oriel Instrument, Stratford, CT) whose emission is first spectrally filtered (400–1000 nm) to reject ultraviolet and infrared light, and then focused onto a 3-mm-diameter illumination optical glass fiber bundle, which delivers light with an average illumination power of 280 mW, which translates into a power density of 3.96 W/cm². A 5 mm diameter collection optical glass fiber bundle is located at three positions, at x± 1 cm, on the opposite side of the inclusions (at a y-axis separation of 5 cm) and linearly scanned. Experiments are made with the light source placed in succession at 8 positions with 1 cm increments for total of 24 source-detector pairs. The collection optical fiber delivers light to a spectrograph (Model No. SP-150, Acton Research Corp., Acton, MA), which disperses the light onto the detector array of a charge coupled device (CCD) camera (Model No. DU420A-BR-DD, Andor Technology, South Windsor, CT). Two exposure times are used for the CCD camera to ensure that approximately the same number of photons are collected for reconstructions using 6 wavelengths and for 126 wavelengths. To this end, two exposure times were used, a longer one of 10 s for the 6 wavelength case and 500 ms for the 126 wavelength case. Since the goal of this paper is to demonstrate the improvement of including hyperspectral information we present an ideal case where the signal to noise ratio is large, thereby providing a best-case scenario for the few-wavelength reconstruction against which we compare our approach as well as using realistic absorption contrasts for the inclusions. The spectrograph features a grating blazed at 700 nm with 350 g/mm, resulting in a dispersion of 20 nm/mm at the exit port. The size of the CCD camera pixels of 26 μm×26 μm results in a spectral sampling rate of two data points per nanometer, even though the spectral resolution is not as high because of the size of the entrance slit (2 mm) used to accommodate the large collection optical fiber bundle. From the data we only retain the wavelength band 650–900 nm where the signal-to-noise ratio is adequate.

In our experiment the incident field is a data set taken before the perturbation is put into the medium. The scattered field is then computed as a dataset that has the original unperturbed dataset subtracted from it. For in vivo measurements it is possible to use a priori structural information from other modalities, e.g. MRI, to estimate the incident field by determining the optical properties of the assumed piecewise constant chromophore distribution over these segments. This method has been employed in literature, as in [45]. In that paper the authors used a very reduced order model to estimate the incident field and then imaged the perturbations about that model. Some methods avoid estimating the incident field by processing measurement data to generate relative changes measured by the detector. Barbour et. al. employ this method by computing the temporal mean value of the readings from detectors to be used as reference values [46]. If obtaining or estimating the incident field proves to challenging in a practical setting the need to move to a nonlinear inverse model could prove itself useful even though the computational intensity is higher than for a linear mode.

A comparison of the absolute concentrations, $\hat{c_{i}}$ and relative concentration, $\hat{c_{i}^{r}}$ to target values is done to test the accuracy of the reconstructions. The relative concentrations for ink are calculated as

\hat{c_{ink}^{r}} = \hat{c_{ink}} / (\hat{c_{ink}} + \hat{c_{dye}})

and similarly for dye [42, 43]. The relative concentration is calculated from the peak value in each reconstruction. This allows us to inspect how well our approach manages to separate and estimate each species of chromophores in the process.

To quantitatively analyse the localization of the reconstruction, we employ the Dice coefficient to judge how well the reconstruction locate the inclusion for the experimental measurements [48,49]. If S is the reconstructed image and G is the ground truth created for each set the Dice coefficient between S and G can be defined as

D (S, G) = \frac{2 | S \cap G |}{| S | + | G |}

|S ∩ G| contains all pixels which both belong to the detected segment as well as the ground truth segment, so if S and G are equal the Dice coefficient is equal to one, indicating an accurate reconstruction.. To compute the D(S,G) the reconstructed images need to be converted to binary maps, where reconstructed concentrations are denoted by 1’s. To do this we threshold our reconstructed images at 10% to 90% of the maximum value of each reconstructed image and inspect D(S,G) for each threshold.

6. Results

6.1. Simulations

In Fig. 5 reconstruction results for the first set are shown for 6 wavelengths, λ = [660, 734, 760, 808, 826, 850] nm, which are optimally chosen according to [44] and also for hyperspectral reconstruction using 126 wavelengths, which are equally spaced over the 650–900 nm range. In simulations the SNR is set to 40 dB, as it is defined by Eq. (13) and (14). For the 6 wavelength case the reconstruction of the HbR chromophore almost totally fails. The concentration is diffused over the whole left half of the region not achieving a good localization. The reconstruction for the HbO₂ is somewhat better but the localization is again not successful. The size of the reconstructed anomaly approaches the true shape but its actual location is off. The values of the concentrations in both cases is underestimated by a significant amount. By employing 126 wavelength, equally spaced over the spectrum, the reconstruction successfully localizes the concentrations and the estimate of the values is close to the ground truth, resulting in the MSE values of 0.17 and 0.16 for HbO₂ and HbR, respectively. It should be noted that crosstalk in the 126 wavelength case is not noticeable generating a good estimation of the chromophore concentration.

The optimal choice of regularization parameters is shown on the L-hypersurface in Fig. 6, using 126 wavelengths, where it is evident that the best parameters used for the reconstruction shown in Fig. 5 occur at a corner in the hypersurface where the estimation error is minimized. The optimal values are α ₁ = 0.1045 and α ₂ = 0.101 resulting in a estimation error of 0.067.

For the second more complicated set reconstructions using 6 and 126 wavelengths are shown in Fig. 7. When using 6 wavelengths the reconstruction is diffused and exact localization of the concentrations is hard to achieve. Noticeably the concentration for the HbO₂ is better defined than for the HbR, where it is almost overshadowed by the concentration in the background. Increasing the number of wavelengths to 126 the reconstruction is much more accurate, localizing the chromophore concentration to the correct areas and reach the target values resulting in the MSE values of 0.17 and 0.07 for HbO₂ and HbR, respectively.

The choice of optimal parameters, for the second set, is done in the same way as for the first set. In this more complex and realistic reconstruction the benefit to separating the regularization parameters is clear. The optimal choice of parameters is in the corner of the hypersurface which corresponds to a low estimation error. The optimal values are α ₁ = 0.14 and α ₂ = 0.37 resulting in a estimation error of 0.05.

To emphasize the advantage of the hyperspectral information, the MSE error can be inspected when different number of wavelengths are used in the reconstruction. In Table 1 the MSE is shown as a function of wavelengths for different sets of equally spaced wavelengths, except in the 6 wavelength case where they are optimally chosen according to [44]. In each case optimal regularization parameters are used to obtain the best reconstruction. This shows that the benefit of hyperspectral information is significant, even though the HbO₂ species shows a small decrease from using 63 wavelengths to 126 wavelengths the benefit for HbR is larger. This is promising considering more complicated reconstruction schemes which have to take into account a higher number of chromophores. The Dice coefficients are shown for both simulation sets in Figs. 8(a) and 8(b). The improvement is evident in both simulation sets, especially for the first set, where the multispectral reconstructed concentration were diffuse over the medium,. Combined with the MSE results and the Dice coefficient, hyperspectral information shows significant improvement over multispectral reconstructions.

Fig. 8 (a) Dice coefficients for the first simulation set as a function of threshold. (b) Dice coefficients for the second simulation set as a function of threshold.

Download Full Size | PDF

Table 1. MSE Compared to the Number of Wavelengths used in the Reconstruction^{^a}

View Table | View all tables in this article

6.2. Experimental validation

Reconstructions of absolute concentrations are shown in Fig. 10. In both cases the hyperspectral reconstructions show better results than reconstructions using 6 wavelengths. The improvement from using hyperspectral information is especially notable in the localization of the inclusions where the reconstruction is diffuse for the 6 wavelength images. When the relative concentration values are compared in Table 2 in both cases the hyperspectral information estimates the value more accurately, where the best estimation is highlighted in bold in the table. For experimental set 2 where the phantom contains 70% and 30% the MSE is 0.18 and 0.29 for ink and dye respectively. For experimental set 1 the hyperspectral reconstruction shows significant improvement both in terms of the quantitative accuracy of the recovered relative concentrations and the localization of the objects as measured by the Dice coefficient. Examining values for the Dice coefficient, shown in Fig. 9, it is evident that the inclusions are localized better when using hyperspectral information. In experimental set 1 the localization of the Ink chromophore completely fails for the multispectral reconstruction, shown in Fig. 9(a), where the hyperspectral information gives good results. It should be noted that the ink chromophore has an absorption peak around 580 nm where experimental set 1 has a very low contrast, as is shown in Fig. 3(b). This causes difficulty in recovering the concentration values. The reconstructions shown in Fig. 10(a) and 10(c), show that the hyperspectral information results in a better reconstruction although it is still significantly diffused.

Fig. 9 (a) Dice coefficients for the first experimental set as a function of threshold. (b) Dice coefficients for the second experimental set as a function of threshold.

Download Full Size | PDF

Fig. 10 Reconstruction from both experimental sets, set 1 containing 10% ink and 90% dye and set 2 70% ink and 30% dye.

Download Full Size | PDF

Table 2. Comparison of $\hat{c_{i}}$ and $\hat{c_{i}^{r}}$ to Target Values for Experimental Results^{^a}

View Table | View all tables in this article

7. Conclusion

We have shown through simulations and measurements that using hyperspectral information for the DOT problem can greatly help in creating concentration images of multiple chromophores. Not only can images be created simultaneously for many chromophores, including information for tens or even hundreds of wavelength information results in increased accuracy in the image and reduces cross-talk.

Optimizing the regularization parameter associated with the individual chromophores is an important element to obtaining a good reconstruction. The effects of changing the regularization constants were examined and chosen to get the best reconstruction. Together, these two factors eliminate artifacts in the image and evoke interest in using this kind of technique for physical measurements. Simulated reconstructions showed significant improvements of including hyperspectral data along with the importance of using multiple regularization parameters for each chromophore.

Physical measurements were also performed to demonstrate these advantages for actual measurement data. Although exact values of concentrations were not achieved there is a notable improvement shown when using hyperspectral information. Additionally, improved localization of inclusions was observed for both sets when using hyperspectral information. This is especially for experimental set 1 in Fig. 10(a) and (c). This emphasizes the advantage of hyperspectral information when doing reconstructions for more than one chromophore.

In this paper we have demonstrated that the hyperspectral information improves reconstructions. One area of future work is to explore these gains in a more quantitative manner. Nonlinear forward models can be employed to improve quantitative accuracy, although this would introduce computational issues. This could be addressed by recycling Krylov subspace information for the most efficient solution [51]. Under the Born model at least, one way to accomplish quantitative accuracy is through analysis of the singular value decomposition of the overall K matrix as is done in [2]. There Gaudette used the singular value spectrum to quantify how the addition of a second wavelength reduced the ill-posedness of the inverse problem. This method could be considered for HyDOT although significant numerical challenges would need to be overcome due to the size of the full hyperspectral K matrix as is discussed in Section 2.

The results are encouraging and demonstrate the potential of HyDOT for multiple chromophore reconstruction. For future work we will also consider spatially varying scattering. Imaging scattering has some challenges and it has been stated that the cw method lacks the capability of separating absorption from scattering in the DOT image reconstruction [17, 44], but it has been shown that preconditioning of data, regularization and use of multispectral information can separate absorption from scattering coefficient [16] and result in good image reconstructions. We anticipate considering the use of a level set approach [47] to estimate the chromophore concentrations as well as scattering properties.

Performing hyperspectral reconstructions results in longer computational time, therefore it would be interesting to optimize the code by moving from MATLAB code and implementing it in a lower level language. This approach could take advantage of hardware acceleration(multi-core or graphical processing units) to lower computational time. The use of hardware acceleration has shown promise in literature [52, 53].

Additionally, implementing semi-infinite boundaries would be interesting to represent a more realistic setting. The switch from infinite to semi-infinite boundaries should benefit our method by generating more accurate reconstructions of different chromophores.

Acknowledgments

The authors thank Prof. Angelo Sassaroli, Yang Yu, Ning Liu and Debbie Chen for their help in the experimental measurements and for all their input. The authors would like to acknowledge financial support provided by the National Institute of Health grant R01 NS059933.

References and links

1. S. Fantini, E. L. Heffer, V. E. Pera, A. Sassaroli, and N. Liu, “Spatial and spectral information in optical mammography,” Technol. Cancer Res. Treat. 4, 471–482 (2005).

2. R. J. Gaudette, D. H. Brooks, C. A. DiMarzio, M. E. Kilmer, E. L. Miller, T. Gaudette, and D. A. Boas, “A comparison study of linear reconstruction techinques for diffuse optical tomographic imaging of absorption coefficient,” Phys. Med. Biol. 45, 1051–1069 (2000). [CrossRef] [PubMed]

3. D. Boas, D. Brooks, E. Miller, C. DiMarzio, M. Kilmer, R. Gaudette, and Q. Zhang, “Imaging the body with diffuse optical tomography,” IEEE Trans. Signal Process. 18, 57–75 (2001). [CrossRef]

4. M. Schweiger and S. R. Arridge, “Optical tomographic reconstruction in a complex head model using a priori region boundary information,” Phys. Med. Biol. 44, 2703–2721 (1999). [CrossRef] [PubMed]

5. S. Fantini, D. Hueber, M. A. Franceschini, E. Gratton, W. Rosenfeld, P. G. Stubblefield, D. Maulik, and M. R. Stankovic, “Non-invasive optical monitoring of the newborn piglet brain using continuous-wave and frequency-domain spectroscopy,” Phys. Med. Biol. 44, 1543–1563 (1999). [CrossRef] [PubMed]

6. J. P. Culver, R. Choe, M. J. Holboke, L. Zubkov, T. Durduran, A. Slemp, and A. G. Yodh, “Three-dimensional diffuse optical tomography in the parallel plane transmission geometry: Evaluation of a hybrid frequency domain/continuous wave clinical system for breast imaging,” Med. Phys. 30, 235–247 (2002). [CrossRef]

7. A. Li, G. Boverman, Y. Zhang, D. Brooks, E. L. Miller, M. E. Kilmer, Q. Zhang, E. M. C. Hillman, and D. A. Boas, “Optimal linear inverse solution with multiple priors in diffuse optical tomography,” Appl. Opt. 44, 1948–1956 (2005). [CrossRef] [PubMed]

8. H. Dehghani, B. W. Pogue, S. P. Poplack, and K. D. Paulsen, “Multiwavelength three-dimensional near-infrared tomography of the breast: initial simulation, phantom, and clinical results,” Appl. Opt. 42, 135–145 (2004). [CrossRef]

9. G. Boverman, Q. Fang, S. A. Carp, E. L. Miller, D. H. Brooks, J. Selb, R. H. Moore, D. B. Kopans, and D. A. Boas, “Spatio-temporal imaging of the hemoglobin in the compressed breast with diffuse optical tomography,” Phys. Med. Biol. 52, 3619–3641 (2007). [CrossRef] [PubMed]

10. A. Corlu, T. Durduran, R. Choe, M. Schweiger, E. M. C. Hillman, S. R. Arridge, and A. G. Yodh, “Uniqueness and wavelength optimization in continuous-wave multispectral diffuse optical tomography,” Opt. Lett. 28, 2339–2341 (2003). [CrossRef] [PubMed]

11. S. Lam, F. Lesage, and X. Intes, “Time domain fluorescent diffuse optical tomography: analytical expressions,” Opt. Express 13, 2263–2275 (2005). [CrossRef] [PubMed]

12. J. P. Culver, V. Ntziachristos, M. J. Holboke, and A. G. Yodh, “Optimization of optode arrangements for diffuse optical tomography: a singular-value analysis,” Opt. Lett. 26, 701–703 (2001). [CrossRef]

13. A. Li, Q. Zhang, J. P. Culver, E. L. Miller, and D. A. Boas, “Reconstructing chromosphere concentration images directly by continuous diffuse optical tomography,” Opt. Lett. 29, 256–258 (2004). [CrossRef] [PubMed]

14. B. W. Pogue, T. O. McBride, J. Prewitt, U. L. Osterberg, and K. D. Paulsen, “Spatially variant regularization improves diffuse optical tomography,” Appl. Opt. 38, 2950–2961 (1999). [CrossRef]

15. P. K. Yalavarthy, B. W. Pogue, H. Dehghani, C. M. Carpenter, S. Jiang, and K. D. Paulsen, “Structural information within regularization matrices improves near infrared diffuse optical tomography,” Opt. Express 15, 8043–8058 (2007). [CrossRef] [PubMed]

16. C. Li, S. R. Grobmyer, L. Chen, Q. Zhang, L. L. Fajardo, and H. Jiang, “Multispectral diffuse optical tomography with absorption and scattering spectral constraints,” Appl. Opt. 46, 8229–8236 (2007). [CrossRef] [PubMed]

17. A. Corlu, R. Choe, T. Durduran, K. Lee, M. Schweiger, S. R. Arridge, E. M. Hillman, and A. G. Yodh, “Diffuse optical tomography with spectral constraints and wavelength optimization,” Appl. Opt. 44, 2082–2093 (2005). [CrossRef] [PubMed]

18. D. A. Boas, “A fundamental limitation of linearized algorithms for diffuse optical tomography,” Opt. Express 1, 404–413 (1997). [CrossRef] [PubMed]

19. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15, R41–R93 (1999). [CrossRef]

20. M. Belge, M. E. Kilmer, and E. L. Miller, “Efficient determination of multiple regularization parameters in a generalized l-curve framework,” Inverse Probl. 18, 1161–1183 (2002). [CrossRef]

21. M. A. OLeary, D. A. Boas, B. Chance, and A. G. Yodh, “Experimental images of heterogeneous turbid media by frequency-domain diffusing-photon tomography,” Opt. Lett. 20, 426–428 (1995). [CrossRef]

22. B. Brooksby, B. W. Pogue, S. Jiang, H. Dehghani, S. Srinivasan, C. Kogel, T. D. Tosteson, J. Weaver, S. P. Poplack, and K. D. Paulsen, “Imaging breast adipose and fibroglandular tissue molecular signatures by using hybrid MRI-guided near-infrared spectral tomography,” Proc. Natl. Acad. Sci. U.S.A. 103, 8828–8833 (2006). [CrossRef] [PubMed]

23. J. Wang, S. C. Davis, S. Srinivasan, S. Jiang, B. W. Pogue, and K. D. Paulsen, “Spectral tomography with diffuse near-infrared light: inclusion of broadband frequency domain spectral data,” J. Biomed. Opt. 13, 1–10 (2008). [CrossRef]

24. D. Boas, M. A. O’Leary, B. Chance, and A. Yodh, “Scattering of diffuse photon density waves by spherical inhomogeneities within turbid media: analytic solution and applications,” Proc. Natl. Acad. Sci. U.S.A. 91, 4887–4891 (1994). [CrossRef] [PubMed]

25. A. Mandelis, Diffusion-Wave Fields: Mathematical Methods and Green Functions, 1st ed.(Springer, 2001).

26. M. Guven, B. Yazici, X. Intes, and B. Chance, “Diffuse optical tomography with a priori anatomical information,” Phys. Med. Biol. 50, 2837–2858 (2005). [CrossRef] [PubMed]

27. A. J. Laub, Matrix Analysis for Scientists and Engineers, 1st ed. (SIAM: Society for industrial and applied mathematics, 2004).

28. T. Coleman and Y. Li, “On the convergence of reflective newton methods for large-scale nonlinear minimization subject to bounds,” Math. Prog. 61, 189–224 (1994). [CrossRef]

29. C. C. Paige and M. A. Saunders, “Lsqr: An algorithm for sparse linear equations and sparse least squares,” ACM Trans. Math. Softw. 8, 43–71 (1982). [CrossRef]

30. F. Larusson, “Hyperspectral imaging for diffuse optical tomography,” M.S. Thesis (Department of Electrical Engineering, Tufts University, Medford, 2009).

31. B. Chen, J. J. Stamnes, and K. Stamnes, “Reconstruction algorithm for diffraction tomography of diffuse photon density waves in a random medium,” J. Eur. Opt. Soc. A 7, 1161–1180 (1998). [CrossRef]

32. S. D. Konecky, R. Choe, A. Corlu, K. Lee, R. Wiener, S. M. Srinivas, J. R. Saffer, R. Freifelder, J. S. Karp, N. Hajjioui, F. Azer, and A. G. Yodh, “Comparison of diffuse optical tomography of human breast with whole-body and breast-only positron emission tomography,” Med. Phys. 35, 446–455 (2008). [CrossRef] [PubMed]

33. P. Taroni, A. Pifferi, A. Torricelli, D. Comelli, and R. Cubeddu, “In vivo absorption and scattering spectroscopy of biological tissues,” Photochem. Photobiol. Sci. 2, 124–129 (2003). [CrossRef] [PubMed]

34. S. Prahl, “Tabulated molar extinction coefficient for hemoglobin in water,” http://omlc.ogi.edu/spectra/hemoglobin/summary.html(Oregon Medical Laser Center, 2007).

35. B. Brendel and T. Nielsen, “Selection of optimal wavelengths for spectral reconstruction in diffuse optical tomography,” J. Biomed. Opt. 14, 1–10 (2009). [CrossRef]

36. S. Srinivasan, B. W. Pogue, S. Jiang, H. Dehghani, C. Kogel, S. Soho, J. J. Gibson, T. D. Tosteson, S. P. Poplack, and K. D. Paulsen, “Interpreting hemoglobin and water concentration, oxygen saturation, and scattering measured in vivo by near-infrared breast tomography,” Proc. Natl. Acad. Sci. U.S.A. 100, 12349–12354 (2003). [CrossRef] [PubMed]

37. H. FlandersDifferential Forms with Applications to the Physical Sciences, 1st ed. (Dover, 1989).

38. T. Durduran, R. Choe, J. P. culver, L. Zubkov, M. J. Holboke, J. Giammarco, B. Chance, and A. G. Yodh, “Bulk optical properties of healthy female breast tissue,” Phys. Med. Biol. 47, 2847–2861 (2004). [CrossRef]

39. L. Spinelli, A. Torricelli, A. Pifferi, P. Taroni, G. M. Danesini, and R. Cubeddu, “Bulk optical properties and tissue components in the female breast from multiwavelength time-resolved mammography,” J. Biomed. Opt. 9, 1137–1142 (2004). [CrossRef] [PubMed]

40. E. Gratton, S. Fantini, M. A. Franceschini, C. Gratton, and M. Fabiani, “Measurements of scattering and absorption changes in muscle and brain,” Philos. Trans. R. Soc. London, Ser. B 352, 727–735 (1997). [CrossRef] [PubMed]

41. B. W. Pogue, S. Jiang, H. Dehghani, C. Kogel, S. Soho, S. Srinivasan, X. Song, T. D. Tosteson, S. P. Poplack, and K. D. Paulsen, “Characterization of hemoglobin, water, and NIR scattering in breast tissue: analysis of intersubject variability and menstrual cycles changes,” J. Biomed. Opt. 9, 541–552 (2004). [CrossRef] [PubMed]

42. S. Fantini, M. A. Franceschini, J. S. Maier, S. A. Walker, B. Barbieri, and E. Gratton, “Frequency-domain multi-channel optical detector for noninvasive tissue spectroscopy and oximetry,” Opt. Eng. 34, 32–42 (1995). [CrossRef]

43. M. A. Franceschini, V. Toronov, M. E. Filiaci, E. Gratton, and S. Fantini, “On-line optical imaging of the human brain with 160-ms temporal resolution,” Opt. Express 6, 49–57 (2000). [CrossRef] [PubMed]

44. M. E. Eames, B. W. Pogue, and H. Dehghani, “Wavelength band optimization in spectral near-infrared optical tomography improves accuracy while reducing data acquisition and computational burden,” J. Biomed. Opt. 13, 1–9 (2008). [CrossRef]

45. G. Boverman, E. L. Miller, A. Li, Q. Zhang, T. Chaves, D. H. Brooks, and D. A. Boas, “Quantitative spectroscopic diffuse optical tomography of the breast guided by imperfect a priori structural information,” Phys. Med. Biol. 50, 3941–3956 (2005). [CrossRef] [PubMed]

46. 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 18, 3018–3036 (2001). [CrossRef]

47. M. E. Kilmer, E. L. Miller, A. Barbaro, and D. Boas, “Three-dimensional shape-based imaging of absorption perturbation for diffuse optical tomography,” Appl. Opt. 42, 3129–3144 (2003). [CrossRef] [PubMed]

48. A. A. Joshi, A. J. Chaudhari, D. W. Shattuck, J. Dutta, R. M. Leahy, and A. W. Toga, “Posture matching and elastic registration of a mouse atlas to surface topography range data,” Proc IEEE Intl. Symp. Biomedical Imaging, 366–369 (2009).

49. J. D. Vylder and W. Philips, “A computational efficient external energy for active contour segmentation using edge propagation,” Intl. con. on Im. Processing, (2010).

50. N. Liu, A. Sassaroli, and S. Fantini, “Paired-wavelength spectral approach to measuring the relative concentrations of two localized chromophores in turbid media: an experimental study,” J. Biomed. Opt. 12, 05160 (2007). [CrossRef]

51. M. E. Kilmer and E. de Sturler, “Recycling subspace information for diffuse optical tomography,” SIAM J. Sci. Comput. 27, 2140–2166 (2004). [CrossRef]

52. S. Bélanger, M. Abran, X. Intes, C. Casanova, and F. Lesage, “Real-time diffuse optical tomography based on structured illumination,” J. Biomed. Opt. 15, 016006 (2010). [CrossRef] [PubMed]

53. Q. Fang and D. a. Boas, “Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units,” Opt. Express 17, 20178–20190 (2009). [CrossRef] [PubMed]

(a) First simulation set			(b) Second simulation set
# λ	MSE HbO₂	MSE HbR	# λ	MSE HbO₂	MSE HbR
6 [44]	0.217	0.203	6 [44]	0.20	0.11
9	0.20	0.20	9	0.18	0.11
26	0.18	0.16	26	0.19	0.09
63	0.17	0.18	63	0.18	0.09
126	0.17	0.16	126	0.17	0.07

(a) Experimental set 1, 10% ink and 90% dye						(b) Experimental set 2, 70% ink and 30% dye
Fig.	# λ	Species	$\hat{c_{i}} [%]$	$\hat{c_{i}^{r}} [%]$	MSE	Fig.	# λ	Species	$\hat{c_{i}} [%]$	$\hat{c_{i}^{r}} [%]$	MSE
10a	6	Ink	1	4	0.64	10b	6	Ink	56	82	0.18
10a	6	Dye	27	96	0.07	10b	6	Dye	12	18	0.41
10c	126	Ink	17	16	0.62	10d	126	Ink	65	61	0.18
10c	126	Dye	88	84	0.07	10d	126	Dye	41	39	0.29

(a) First simulation set			(b) Second simulation set
# λ	MSE HbO₂	MSE HbR	# λ	MSE HbO₂	MSE HbR
6 [44]	0.217	0.203	6 [44]	0.20	0.11
9	0.20	0.20	9	0.18	0.11
26	0.18	0.16	26	0.19	0.09
63	0.17	0.18	63	0.18	0.09
126	0.17	0.16	126	0.17	0.07

(a) Experimental set 1, 10% ink and 90% dye						(b) Experimental set 2, 70% ink and 30% dye
Fig.	# λ	Species	$\hat{c_{i}} [%]$	$\hat{c_{i}^{r}} [%]$	MSE	Fig.	# λ	Species	$\hat{c_{i}} [%]$	$\hat{c_{i}^{r}} [%]$	MSE
10a	6	Ink	1	4	0.64	10b	6	Ink	56	82	0.18
10a	6	Dye	27	96	0.07	10b	6	Dye	12	18	0.41
10c	126	Ink	17	16	0.62	10d	126	Ink	65	61	0.18
10c	126	Dye	88	84	0.07	10d	126	Dye	41	39	0.29

Hyperspectral image reconstruction for diffuse optical tomography

Abstract

1. Introduction

2. Forward problem

3. Inverse problem

4. Simulation analysis

4.1. Using gradient matrix and two parameters

5. Experimental analysis

6. Results

6.1. Simulations

6.2. Experimental validation

7. Conclusion

Acknowledgments

References and links

Cited By

Figures (10)

Tables (2)

Equations (28)

Biomedical Optics Express