Parametric level set reconstruction methods for hyperspectral diffuse optical tomography

Fridrik Larusson; Sergio Fantini; Eric L. Miller

doi:10.1364/BOE.3.001006

1. Introduction

Near-infrared light has proven to be useful for imaging the human body and providing functional information for applications including breast cancer detection and characterization [1–5]. In recent years significant improvements have been made in diffuse optical tomography (DOT) by using data from multiple wavelengths. The use of this kind of hyperspectral information has moved DOT away from the recovery of space and time varying maps of absorption and scattering properties of the breast to the recovery of chromophore concentrations as well as diffusion amplitude. Inverting for multiple chromophores such as hemoglobin, lipids and water has been shown to provide an improved ability to localize tumours or other objects of interest in the breast cancer application [6–8]. While these and other advancements have been critical for moving DOT from the lab into the clinic, there remain significant obstacles to be addressed so that the method can be used to stably recover these quantities.

The imaging of chromophore concentration and diffusion amplitude from DOT data requires the solution of an ill-posed and non-linear inverse problem [3]. The physics associated with the photons travelling through turbid media, such as breast tissue, makes the problem ill-posed and therefore introduces significant challenges for the DOT problem. The ill-posedness make the reconstruction highly sensitive to noise and un-modeled effects which can reduce accuracy in the recovered images [7]. Additional problems are encountered when recovering multiple chromophores such as crosstalk, e.g. when two spatially disjoint chromophores pollute the images of each other.

The ill-posedness of the DOT problem is traditionally solved by putting the recovery of images in a variational context. To achieve this, numerous regularization techniques [9] can be used to reduce artifacts and improve the overall accuracy. Traditionally, methods such as Tikhonov and total variation (TV) regularization have been employed along with L-curve methods to optimally choose regularization parameters [10–12]. Increasing the amount of data used to solve the DOT problem has also proven effective in generating high accuracy images [8, 13]. To this end hyperspectral information is employed, which along with regularization has been shown to reduce the non-uniqueness of the solution to the image formation problem [8, 10, 11].

In this paper, we expand on our previous work in [10] where we employed hyperspectral information for recovery of chromophore concentrations. Implementing the Born approximation along with Tikhonov regularization we demonstrated that hyperspectral information was shown to have utility for a pixel-based reconstruction of chromophores [10]. The large data set introduced significant challenges in constructing and computing with the forward model and great care had to be taken in choosing the regularization parameters. The known limitations of the Born approximation and the high number of unknowns encountered in the pixel-based approach, made reconstructions sensitive to the choice of regularization parameters [14]. This motivated us to move to a geometric reconstruction approach based on a low-order parametric model that has been shown to perform well in the context of ill posed inverse problems [15–17].

The primary contribution of this work is the development of a new approach to the multi-chromophore inverse problem in which the medium is or can be well approximated as piecewise constant. Piecewise constant DOT problems have been considered in the past. In [18] level sets were used in a two-step method for shape estimation assuming that prior information of the absorption parameter was known. Schweiger et al. [19] and Kilmer et al. [20] employed level sets for the DOT problem estimating parameter distributions using a piecewise basis. Arridge et al. investigated shape based methods by estimating level-sets, specifically investigating an explicit method using basis functions and an implicit shape reconstruction to recover absorption and diffusion assuming a known background [21]. For a detailed review of the use of level sets in inverse scattering problems we point the reader to [22].

The approach we consider in this paper is significantly different from those in [18, 19, 21]. In addition to the fact that none of these papers have considered the fully hyperspectral case, some of these methods [18–20] require the recovery of unknown quantities defined on a fine scale pixelated discretization of the region of interest. More specifically in [21] absorption and scattering is estimated using level sets assuming the background known. With the Born approximation we assume the absorption and diffusion is known in the background but here we estimate chromophore concentrations and diffusion of the object of interest as well as chromophore concentration in the background. TV reconstructions recover images while traditional level set methods work with a level set function defined on a pixel-based grid. In both cases, regularization is required to obtain adequate results and one is faced with the corresponding challenge of choosing regularization parameters [11, 19].

In this paper we consider the use of a shape-based approach to the hyperspectral DOT problem based on a newly-developed parametric level set (PaLS) formulation. In [15], a basis function expansion was used to provide a low order representation of the level set function and yielded very strong results for a number of highly ill-posed inverse problems including a restricted form of the DOT problem where a single wavelength was employed to determine only optical absorption. The method in [15] required no explicit regularization and, due to the low-order nature of the model (number of parameters much, much less than number of pixels) was amenable to Newton-type inversion algorithms known to converge more rapidly than gradient-based schemes. Moreover, in [15] our experiments indicated a surprising lack of sensitivity to the initialization of the inversion algorithm.

Although the breast is a highly heterogeneous medium, in the level set method we assume the images to be recovered to be piecewise constant. This approximation is supported in the literature. For example Schweiger et al. assumed anatomical prior information to derive a piecewise constant region basis [4]. The co-located structures of different species of chromophores is discussed further in Section 5. In the future, where we move to more complicated simulations and experiments the goal would be to implement multiple characteristic functions that require the estimation of multiple level sets, discussed further in Section 8. To prove the accuracy of our method, we use simulation data to show parametric level-set reconstructions for multiple chromophores and diffusion amplitude. Experimental data are used to show the improvement of the PaLS method over pixel based methods using previously reported data sets [10].

The remainder of this paper is structured as follows. In Section 2 we discuss our forward model and derive in detail the Born model used. In Section 3 we discuss the PaLS method and in Section 4 we derive the Levenberg Marquardt algorithm used for the inverse problem. In Section 7.1 we present simulation results with multiple chromophore and diffusion amplitude reconstructions, and in Section 7.2 we present experimental results that show the improvement of the PaLS method over pixel based methods.

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 [23, 24] of the form

\nabla \cdot D^{0} (r, λ) \nabla Φ (r, λ) + v μ_{a}^{0} (r, λ) Φ (r, λ) = - v 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 at wavelength λ. For spatially varying scattering we assume that the diffusion coefficient D⁰(r,λ) follows Mie scattering theory where a scattering prefactor Ψ depends on the size and density of scatterers while a scattering exponent b depends on the size of the scatterers. Using this, we write the perturbation in the diffusion coefficient as

Δ D (r, λ) = \frac{v}{3 Δ Ψ} {(\frac{λ}{λ_{0}})}^{b} = v Δ Ψ^{'} d (λ) .

The arbitrarily chosen reference wavelength λ₀ is introduced to achieve a form of the Mie model where Ψ has the units of mm⁻¹ and Ψ′ has units of mm. In this paper, for simplicity we consider the case where the background medium is infinite and homogeneous. Generalization to the more practical case where there are boundaries is straightforward in theory though somewhat more complex in terms of implementation [25]. As the primary objective of this work is the demonstration of a new approach to inversion, we prefer to consider the simpler physical problem holding off on the more realistic implementation for the future. Additionally, we only consider continuous wave experiments since this is what our instrument measures, where Eq. (1) is usually written with a jω/D(λ) term on the right hand side, in our case we consider ω = 0 giving us the form shown in Eq. (1).

We employ the Born approximation by decomposing $μ_{a}^{0} (r, λ)$ and D⁰(r,λ), as the sum of a constant background absorption, μ_a(λ), and a spatially varying perturbation Δμ_a(r,λ) as well as constant background diffusion D(λ) and a diffusion perturbation ΔD(r,λ). To obtain a linear relationship between the measurements and the chromophore concentrations, we subtract Eq. (1) from the perturbed version of the diffusion equation which gives

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

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 solution of Eq. (3) as is the case for an unbounded medium as well as range of bounded problems [26], we can change Eq. (3) into an integral Eq. [27]

\begin{array}{r} Φ_{s} (r, λ) = - \frac{v}{D (λ)} (\int Δ μ (r^{'}, λ) G (r, r^{'}, λ) Φ_{i} (r^{'}, λ) d r^{'} \\ + \int Δ D (r^{'}, λ) \nabla G (r_{d}, r_{j}, λ) \cdot \nabla Φ_{i} (r_{j}, r_{s}, λ) d r^{'}) \end{array}

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. It should also be noted that we obtain this equation under the assumption that the total fluence rate, Φ, can be approximated as the incident field, Φ_i, since Φ_i ≫ Φ_s [2].

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 [28]

Δ μ_{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 4, 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 Eqs. (5) and (6) in Eq. (4) and performing some algebra we obtain

\begin{array}{l} Φ_{s} (λ) & = & a \frac{v}{D (λ)} \sum_{k, j = 1}^{N_{c}, N_{V}} (G (r_{d}, r_{j}, λ) Φ_{i} (r_{j}, r_{s}, λ) ε_{i} (λ) c_{k, j} \\ + & \nabla G (r_{d}, r_{j}, λ) \cdot \nabla Φ_{i} (r_{j}, r_{s}, λ) Δ D_{j} (λ)) . \end{array}

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. (8) we use a as the area of a pixel. This setup is illustrated 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) Definition of domains used for the parametric level-set methods.

Download Full Size | PDF

Setting up the linear algebraic structure associated with Eq. (8) we define c_k ∈ 𝕉^N_v as the vector obtained by lexicographically ordering the unknown concentrations associated with the k^th chromophore and c_k+1 = ΔΨ′ and Φ_s(λ) is 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. (8) is written in matrix-vector notation as

[\begin{matrix} Φ_{s} (λ_{1}) \\ Φ_{s} (λ_{2}) \\ ⋮ \\ Φ_{s} (λ_{N_{λ}}) \end{matrix}] = [\begin{matrix} ε_{1} (λ_{1}) K_{1}^{a} & ε_{2} (λ_{1}) K_{1}^{a} & \dots & ε_{N_{c}} (λ_{1}) K_{1}^{a} & K_{1}^{d} \\ ε_{1} (λ_{2}) K_{2}^{a} & ε_{2} (λ_{2}) K_{2}^{a} & \dots & ε_{N_{c}} (λ_{2}) K_{2}^{a} & K_{2}^{d} \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ ε_{1} (λ_{N_{λ}}) K_{N_{λ}}^{a} & ε_{2} (λ_{N_{λ}}) K_{N_{λ}}^{a} & \dots & ε_{N_{c}} (λ_{N_{λ}}) K_{N_{λ}}^{a} & K_{N_{λ}}^{d} \end{matrix}] [\begin{matrix} c_{1} \\ c_{2} \\ ⋮ \\ \begin{array}{l} c_{N_{c}} \\ Δ Ψ^{'} \end{array} \end{matrix}] \Leftrightarrow Φ_{s} = Kc

It should be noted in Eq. (9) that the matrix has elements which are also the matrices

K_{l}^{a}

and

K_{l}^{d}

. The (m, j)^th element of the

K_{l}^{a}

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

K_{l}^{d}

the matrix elements are given by (v/D(λ_l))∇G(r_m, r_j, λ) ·∇Φ_i(r_j, r_m, λ)ΔD_j(λ).

To evaluate the forward model for realistically sized problems, we compute the N_λ matrices K^a and K^d and store it along with the N_λ × N_c extinction coefficients. This reduces the amount of memory required for the reconstruction.

3. Parametric level-set method

To counter the ill-posedness of the DOT problem we employ a Parametric Level-Set Method (PaLS). For the purpose of our paper we assume that all chromophore concentrations and diffusion coefficient perturbations are co-located. This choice is supported by reports in literature, where decrease in hemoglobin and water concentration along with scattering power are located at the cancer location and the lipid concentrations increase at the same location [29–31]. This means that the geometry of the anomaly in the medium is the same for all chromophores and diffusion amplitude. The domain Ω ⊂ 𝒡 represents the support of the objects of interest, and 𝒡 represents the homogeneous background where the anomaly is located, shown in Fig. 1(b).

Since the same support is used for each object of interest the characteristic function describing the shape is defined as

χ (x, y) = {\begin{array}{l} 1 & if (x, y) \in Ω, \\ 0 & if (x, y) \in ℱ \ Ω . \end{array}

Then each image to be reconstructed can be written as

c_{k} (x, y) = χ (x, y) c_{k}^{a} + [1 - χ (x, y)] c_{k}^{b}

where k = 1, 2,...,N_c + 1. In this formulation the unknown values are the constant concentration values of the anomaly and background,

c_{k}^{a}

and

c_{k}^{b}

respectively.

The characteristic function χ(x,y) is defined to be the zero level set of a Lipschitz continuous object function 𝒪 : 𝒡 → 𝕉 such that 𝒪 > 0 in Ω(x,y), 𝒪 < 0 in Ω\𝒡 and 𝒪(x,y) = 0 in ∂Ω. Using 𝒪(x,y), χ(x,y) is written as

χ (x, y) = H (𝒪 (x, y))

where H is the step function. In practice we use smooth approximations of the step function and the Dirac delta function denoted as H_ε and δ_ε respectively where H_ε is computed as

H_{ε} (x) = {\begin{array}{l} 1 & if x > ε, \\ 0 & if x < - ε . \\ \frac{1}{2} [1 + \frac{x}{ε} + \frac{1}{π} \sin (\frac{π x}{ε})] & if | x | \leq ε . \end{array}

and δ_ε is computed as the derivative of H_ε [32, 33]. As discussed in Section 1 we represent the object function 𝒪(x,y) parametrically, so instead of using a dense collection of pixel or voxel values [34], we represent it by using basis functions

𝒪 (x, y) = \sum_{i = 1}^{L} a_{i} p_{i} (x, y)

where a_i’s are the weight coefficients whereas p_i(x,y) are the functions which belong to the basis set of 𝒫 = {p₁, p₂,..., p_l}. Possible choices for the 𝒫 basis set include polynomial or radial basis functions. For the purpose of this paper we use Gaussian basis function. The width and number of the Gaussians determines how coarse or fine the reconstruction will be, where a choice of few basis functions will, on the one hand, result in a reduced number of unknowns, it will on the other hand, give a coarser estimation of the shape, which can be a problem for imaging finer more complex structures. In the DOT case, where the physics in the forward model will only allow for a coarse reconstruction of the underlying structure, the Gaussian approach is sufficient, especially for the relatively simple geometries and concentrations presented in this paper. When we will move on to a more complicated non-linear model, the choice of the number and type of basis functions will be more important and should be based on a rigid mathematical framework. This is discussed further in Section 8.

In the PaLS formulation, all the parameters of the model are gathered in one vector $θ^{T} = [c_{a}^{1}, c_{a}^{2}, \dots, c_{b}^{1}, a^{T}]$ where a = [a₁,..., a_L]^T. Now our forward model in Eq. (9) can be expressed as

Φ_{s} = K (θ) = Kc (θ)

The forward model has now been parametrized with a vector containing all of the unknowns, which are far fewer then what a pixel based method would attempt to estimate.

4. Inverse problem

The inverse problem, that of using Φ_s to recover the value of c, is solved as a Levenberg-Marquardt optimization problem of the form

\hat{c} = arg min_{c} {‖ W (K (θ) - Φ_{s}) ‖}_{2}^{2}

The W matrix reflects the structure of the noise corrupting the data [2]. While a Poisson model is technically the most appropriate for DOT data, as is frequently done [35] 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 [11]. 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{{S N R}_{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 minimization of the cost function is then achieved by the Levenberg-Marquardt algorithm. For that purpose an error vector is introduced

ε = W (K (θ) - Φ) .

which can be used to write the cost function in term of ε as

M (θ) = ε^{T} ε

In order to employ the Levenberg-Marquardt algorithm, the calculation of the Jacobian matrix J is required. The Jacobian contains derivatives of ε with respect to each element in the parameter vector θ

J = [\frac{\partial ε (θ)}{\partial {c_{a}^{1}, \dots c_{b}^{l}, a}}]

The solution is then obtained by updating θ at each iteration as θⁿ⁺¹ = θⁿ + h where h is the solution to the following linear system

(J^{T} J + ρ I) h = - J^{T} ε with ρ \geq 0

where I is the identity matrix, ρ is the damping parameter affecting the size and direction of h and found via and appropriate line search algorithm [36]. The stopping criteria used when iterating Eq. (21) is the discrepancy principle [37], in that the iterations are stopped when the norm of the residual has reached the noise level within a certain tolerance.

When estimating the parametric vector, we employ a cyclic coordinate decent strategy [38] Essentially this is equivalent to estimating the shape only at even iterations and the concentration values at odd iterations. This is repeated until stopping criteria is reached. This process is expressed in pseudo-code in Algorithm 1, where J_v and J_s denote the Jacobian strictly for the concentration value and shape, respectively, and τ_i represents a tolerance for the stopping critera.

Algorithm 1. Matlab-like code for estimating shape and concentration value

View Table | View all tables in this article

5. Simulation analysis

The arrangement of sources and detectors with respect to the simulated turbid medium 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 [28]. The source-detector separations were set to 5 cm as is shown in Fig. 1(a). In simulations, we reconstruct concentration images of oxygenated and deoxygenated hemoglobin, HbO₂ and HbR respectively, along with lipid and water concentration and diffusion amplitude. These chromophores are chosen since they mainly cause near-infrared absorption in breast tissue [39], and breast cancer tumours have been found to have higher HbO₂ and HbR concentrations than normal tissue [40].

In simulations, values for Ψ and b are obtained from [41] for the female breast. Values for μ_a are calculated from extinction coefficients, which are in the units cm⁻¹/mM and are obtained from data tabulated by Scott Prahl [42]. The concentration in the simulated images are defined in units of millimolars or millimoles per liter, mM for HbO₂ and HbR. For water and lipid the concentrations are in percent by volume and the diffusion amplitude is measured in units of millimeter. The background has HbR concentration of 0.01 mM, HbO₂ concentration of 0.01 mM, lipid concentration of 32%, water concentration of 13% and Ψ′ is set to 1.6 mm. The target concentration of the object of interest is set to 0.015 mM, 0.012 mM, 50 %, 20 % and 0.25 mm for Hb0₂, HbR, lipid, water and ΔΨ′, respectively.

The simulation set is created with all chromophore concentrations and diffusion perturbations in the same location with different target values. The co-located geometry is chosen since it is likely for chromophore concentrations to have similarly located geometries in the real breast [29, 30]. The ground truth images for simulations are shown in Fig. 6. 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 as well as the performance of the shape based approach. The process is initialized with 21 Gaussian basis functions with width of approximately 8 pixels placed uniformly on a grid over the whole medium to be imaged. A representative image of the order of the basis functions is shown in Fig. 2(a). For all experiments presented in this paper, the a_i’s weight coefficients are initialized to 1.

Fig. 6 Reconstruction using the PaLS method. Middle column of images are generated with 8 wavelengths and rightmost images are generated with 176 wavelengths. From top to bottom the rows show HbO₂, HbR, lipid, water and diffusion amplitude, respectively. Concentration units are in mM.

Download Full Size | PDF

Fig. 2 (a) Order of the basis functions used for the PaLS method. (b) 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

To best understand the utility of a hyperspectral data set, we employ the Born model to generate simulated data. Though this may not be realistic, it allows us to avoid the confounding factor of model mismatch in evaluating the inversion method being considered in this paper. Moreover, the shortcomings of this approach are mitigated in Section 7.2, where we consider the processing of experimental data which, obviously, are not the product of the Born model. In a bit more detail, the data we use for our simulation analysis are computed as

Φ = Kc + n

where c are the simulated concentration images for all chromophores and diffusion amplitude, whereas 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. (16), chosen such that a pre-determined signal to noise ratio (SNR) is achieved. This SNR is calculated from Eq. (17).

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 examining the Dice coefficient to judge how well the concentrations are located. For the k^th chromophore, the mean square error is computed by using the following Eq.

{M S E}_{k} = \frac{{‖ c_{k} - {\hat{c}}_{k} ‖}_{2}}{{‖ c_{k} ‖}_{2}}

Along with using MSE to judge the accuracy of the reconstruction, we also employ a Dice Coefficient to quantitatively analyze how well the shape based approach localizes the reconstruction [43, 44]. 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 that belong to the detected segment as well as the ground truth segment, so that if S and G are equal the Dice coefficient is equal to one, indicating an accurate reconstruction. To compute the D(S,G) we use the characteristic function, χ, which essentially works as a binary map of the reconstructed anomaly where the object of interest is represented by 1’s.

6. Experimental analysis

Physical measurements were performed in order to validate the simulation results, where the data used in this paper has previously been used in [10]. We use the same data again since we are using a completely new image formation process which results in a significant improvement in reconstructions. The background medium consists of water and milk in the ratio of 2:1, respectively. Milk, with 2% fat, is used due to the similarities of the optical properties to human breast tissue. Black India ink and blue food dye were added to mimic tissue 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 [45, 46].

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

In order to obtain multi- and hyperspectral reconstruction values for μ_a and D(λ) the background has to be known. In the experimental measurements we assume constant scattering, therefore we are not trying to estimate the perturbation, ΔD(r,λ), as was done in simulations. Therefore we have the unperturbed representaion of the reduced scattering coefficient, μ′_s, is given by

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

The diffusion coefficient relates to the reduced scattering coefficient through, D(λ) = v/3μ′_s. Phase, amplitude and average intensity data are obtained at two wavelengths using a frequency-domain tissue spectrometer to estimate the Ψ and b parameters in Eq. (25) as Ψ = 6.5 cm⁻¹ and b = 0.4. This allows us to extrapolate values for μ′_s at any wavelength [49]. 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.

In our experiments, two phantom inclusions, named set 1 and set 2, are created 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 contrast range is comparable to traditional tumor contrasts reported in literature, which have been close to 3:1 and lower [50]. 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

For both experimental sets 1 and 2, a single 25 cm long transparent cylindrical inclusion containing ink and dye solutions is placed in the background medium. Optical properties are assumed constant along the z-axis. For the light source, an arc lamp is used which delivers light with an average illumination power of 280 mW, that translates into a power density of 3.96 W/cm². For three different locations, at x ± 1 cm, a 5 mm diameter collection optical glass fiber bundle is placed on the opposite side of the inclusions (at a y-axis separation of 5 cm). Experiments are made with the light source placed in succession at 8 positions with 1 cm increments for total of 24 source-detector pairs. To ensure that approximately the same amount of photons are collected for both hyperspectral and multispectral reconstructions, two exposure times are used for the CCD camera, 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. Further details on the experimental setup can be found at [10].

In our experimental setup, we measure the incident field before the perturbation is put into the medium. Then the scattered field is computed as a dataset that has the 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 [51].

A comparison of the absolute concentrations, $\hat{c_{i}}$ and relative concentration, $\hat{c_{i}^{r}}$ to target concentration 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 [52,53]. The relative concentration is calculated from the peak concentration 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.

7. Results

7.1. Simulations

In Fig. 5 reconstruction results using the pixel based method are shown for 8 wavelengths, λ = [660, 734, 760, 808, 826, 850, 930, 980] nm and hyperspectral reconstruction using 176 wavelengths, which are equally spaced over the 650–1000 nm range. In the 8 wavelength case the 6 first wavelengths are optimally chosen according to [54] with two wavelengths added where water and lipids have peak absorption. Reconstructed images created with the PaLS method are shown in Fig. 6. In simulations the SNR is set to 30 dB, as it is defined by Eq. (16) and Eq. (17). When comparing the pixel based reconstruction in Fig. 5 to the PaLS reconstruction in Fig. 6, it is evident that the PaLS method provides superior reconstructions. Examining the PaLS results, the 8 wavelength case shows reasonable accuracy along the x axis but rather diffuse results in y. We also see noticeable artifacts in the reconstructions. Considering the concentration values, the values for HbO₂, HbR and water concentration come close to the actual value. Moving to hyperspectral information, the reconstruction becomes more accurate, estimating the shape close to the ground truth. It should also be noted that the runtime for each reconstruction for the PaLS method is significantly shorter compared to the pixel-based method. A PaLS reconstruction takes around 30 seconds, which is 3–4 times faster than a pixel-based method. Additionally, we do not employ any regularization parameters, freeing us from finding the optimal reconstruction using regularization. This is a major improvement in moving from a pixel-based approach to the PaLS method.

Fig. 5 Reconstruction using a pixel based method. Middle column of images are generated with 8 wavelengths and rightmost images are generated with 176 wavelengths. From top to bottom the rows show HbO₂, HbR, lipid, water and diffusion amplitude, respectively. Concentration units are in mM.

Download Full Size | PDF

The comparison of the Dice coefficient between the PaLS method and pixel-based is tricky, since for the pixel-based method the Dice coefficient is plotted as a function of a threshold. This threshold is required to create a binary map of the location on the anomaly. If the threshold is chosen to only leave extreme peak concentration values in each image, the Dice coefficient would be low due to edge artifacts as in Fig. 7(b). Therefore, in simulations we compare D(S,G) for the pixel based reconstructions using a threshold of 80% to D(S,G) of the PaLS reconstructions. The improvement of the PaLS method is confirmed quantitatively through D(S,G) and MSE displayed in Table 2 and Table 1, respectively. The Dice coefficient, shown in Table 2, gives a clear view of how the shape estimation improves by added wavelengths, where D(S,G) approaches 1 for the hyperspectral case and the PaLS method shows superior performance in the MSE values.

Fig. 7 Pixel 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 1. The MSE is compared for each chromphore for multiple wavelength choices. In each case the reconstructions are done with equally spaced wavelengths over the spectrum except for the 8 wavelength case.

View Table | View all tables in this article

Table 2. D(S,G) is compared for each chromphore for multiple wavelength choices. In each case the reconstructions are done with equally spaced wavelengths over the spectrum except for the 8 wavelength case. D(S,G) is calculated comparing 80% of the target peak to the reconstruction.

View Table | View all tables in this article

7.2. Experimental validation

Pixel based reconstructions of absolute concentrations for both experimental sets are shown in Fig. 7 and PalS approach reconstructions in Fig. 8. As expected, the hyperspectral information provides improved reconstruction for both the pixel-based and PaLS methods. The pixel based results had previously been reported in [10]. Focusing on the PaLS methods, it is evident that forming the reconstruction with shape-based constraints yields improved results. The estimation of relative concentrations and MSE of the absolute values are examined in Table 4 and Table 5 for the pixel-based and PaLS method, respectively. The relative concentration values are better estimated in both cases using the PaLS method, although the hyperspectral method does not show significant improvement for experimental set 2, which was also the case for the pixel based method. Examining the images along with the MSE values for experimental set 1, Fig. 7–8(a) and (c), it is noticeable how the reconstruction does not resolve the structure particurarly well along the x axis. This is somewhat unexpected since in DOT resolving depth information, on the y axis, is usually the more difficult problem. This is noticeable for both the pixel based and PaLS methods, although the PaLS method outperforms the pixel based method, especially in removing edge artifacts. This streaking in the x direction is most likely a combination of how the Gaussian basis are placed within the imaging medium, and measurement error in placing the source and detectors when taking the reference measurement.

Fig. 8 PaLS 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 3. D(S,G) is compared for each chromphore for multiple wavelength choices. In each case the reconstructions are done with equally spaced wavelengths over the spectrum except for the 6 wavelength case where we use optimally chosen wavelengths. D(S,G) is calculated comparing the half maximum of the target peak to the reconstruction.

View Table | View all tables in this article

Table 4. Comparison of $\hat{c_{i}}$ and $\hat{c_{i}^{r}}$ to target concentration values for experimental results, for the pixel-based method. Best performance is highlighted in bold.

View Table | View all tables in this article

Table 5. Comparison of $\hat{c_{i}}$ and $\hat{c_{i}^{r}}$ to target concentration values for experimental results, for the PaLS method. Best performance is highlighted in bold.

View Table | View all tables in this article

For both experimental sets, the PaLS method resolves the location and the shape of the inclusion more accurately, which is verified by the calculation of the Dice coefficient shown in Table 3, the improvement is notable when compared to the pixel-based reconstruction. As discussed in Section 7.1 a choice of a threshold is needed to compare D(S,G) between pixel based reconstructions and the PaLS method. For the experimental reconstructions we use a threshold of 50% to compare D(S,G) of the PaLS reconstructions. This demonstrates the usefulness of the PaLS method correctly and accurately localizing the anomaly. The PaLS method does very well with eliminating edge artifacts that were severe when doing pixel-based reconstructions for the same data set. These effects are very noticeable in Fig. 7(b) and (d), where, especially in the multispectral case, the edge artifacts were significant. Comparing that to the same data in Fig. 8(b) and (d) it is obvious that the improvement is significant.

8. Conclusion

In this paper, using simulations and experimental measurements we have shown that the PaLS method provides more accurate estimation of chromophore concentrations than a regularized pixel-based inversion scheme. Hyperspectral information results in improved performance in terms of both MSE and spatial localization as measured using the Dice coefficient. The parametric approach is shown to give significant improvements to image reconstruction decreasing run time of the iterative process and increasing the quality of reconstructed images. The PaLS method is also easily expandable to more complicated problems where multiple geometries need to be considered.

Physical measurements were also performed to demonstrate these advantages for actual measurement data. Although exact concentration values were not achieved, there is a notable improvement associated with hyperspectral information in conjunction with the PaLS method. Additionally, improved localization of inclusions was observed for both sets when using hyper-spectral information. This emphasizes the advantage of hyperspectral information when doing reconstructions for more than one chromophore.

Based on the results in this paper, we will be extending the work to address more realistic, clinical conditions. We will first implement a fully 3D model employing boundary conditions and consider an inhomogeneous background more similar to what is found in breast imaging. In that setting, estimating multiple level sets for different geometries could prove useful, especially to estimate different regions of the heterogeneous background such as adipose and fibroglandular tissue. As mentioned in Section 3, it is also important to generate a rigid framework to choose different types, sizes and shapes of basis functions to generate the best reconstruction. To be able to estimate all shapes possible, we aim to increase the number of basis functions to include different types. To avoid over complicating the image reconstruction with a high number of basis functions we aim to pose the image reconstruction in a compressed sensing framework where few optimal basis functions estimate a complex shape. Furthermore, relating to the development of a framework for different basis functions, we will develop the Levenberg Marquardt algorithm to have optimally chosen stopping criteria and step size changes. This should increase robustness of our method and ensure accuracy of estimation for different situations and settings in diffuse optical tomography.

Acknowledgments

The authors thank Prof. Angelo Sassaroli and Alireza Aghasi for useful discussion regarding theory and Yang Yu for his help with experimental measurements. The authors would like to acknowledge financial support provided by the National Institutes of Health grant R01 CA154774.

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). [PubMed]

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 techniques for diffuse optical tomographic imaging of absorption coefficient,” Phys. Med. Biol 45, 1051–1069 (2000). [CrossRef] [PubMed]

3. D. A. Boas, D. H. Brooks, E. L. Miller, C. A. DiMarzio, M. Kilmer, R. J. Gaudette, and Q. Zhang, “Imaging the body with diffuse optical tomography,” IEEE Signal Process. Mag. 18(6), 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. S. Kukreti, A. E. Cerussi, W. Tanamai, D. Hsiang, B. J. Tromberg, and E. Gratton, “Characterization of metabolic differences between benign and malignant tumors: high-spectral-resolution diffuse optical spectroscopy,” Radiology 254, 277–284 (2010). [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. 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]

10. F. Larusson, S. Fantini, and E. L. Miller, “Hyperspectral image reconstruction for diffuse optical tomography,” Biomed. Opt. Express 2, 947–965 (2011). [CrossRef]

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

12. K. D. Paulsen and H. Jiang, “Enhanced frequency-domain optical image reconstruction in tissues through total-variation minimization,” Appl. Opt. 35, 3447–3458, (1996). [CrossRef] [PubMed]

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

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

15. A. Aghasi, M. Kilmer, and E. L. Miller “Parametric level set methods for inverse problems,” SIAM J. Imaging Sci. 4, 618–650 (2011). [CrossRef]

16. A. Aghasi, E. L. Miller, and L. M. Abriola “Characterization of source zone architecture: a joint electrical and hydrological inversion approach,” presented at 2011 Fall Meeting, AGU, San Francisco, Calif., 5–9 Dec. 2011.

17. F. Larusson, S. Fantini, and E. L. Miller, “Parametric level-set approach for hyperspectral diffuse optical tomography,” in 2011 IEEE International Symposium on Biomedical Imaging: from Nano to Macro (IEEE, 2011), pp. 949–955.

18. O. Dorn, “A shape reconstruction method for diffuse optical tomography using a transport model and level sets,” in 2002 IEEE International Symposium on Biomedical Imaging, 2002. Proceedings (IEEE, 2002), pp. 1015–1018.

19. M. Schweiger, O. Dorn, and S. R. Arridge, “3-D shape and contrast reconstruction in optical tomography with level sets,” J. Phys.: Conf. Ser. 124, 012043 (2008). [CrossRef]

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

21. S. R. Arridge, O. Dorn, V. Kolehmainen, M. Schweiger, and A. Zacharopoulos, “Parameter and structure reconstruction in optical tomography,” J. Phys.: Conf. Ser. 135, 012001 (2008). [CrossRef]

22. O. Dorn and D. Lesselier, “Level set methods for inverse scattering,” Inverse Probl. 22, R67–R131 (2006). [CrossRef]

23. M. A. O’Leary, 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]

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

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

26. A. Mandelis, Diffusion-Wave Fields: Mathematical Methods and Green Functions (Springer, 2001).

27. B. Brendel, R. Ziegler, and T. Nielsen “Algebraic reconstruction techniques for spectral reconstruction in diffuse optical tomography,” Appl. Opt. 47, 6392–6403 (2008). [CrossRef] [PubMed]

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

29. A. Cerussi, D. Hsiang, N. Shah, R. Mehta, A. Durkin, J. Butler, and B. J. Tromberg “Predicting response to breast cancer neoadjuvant chemotherapy using diffuse optical spectroscopy,” Proc. Nat. Acad. Sci. U.S.A. 104, 4014–4019 (2007). [CrossRef]

30. H. Soliman, A. Gunasekara, M. Rycroft, J. Zubovits, R. Dent, J. Spayne, M. J. Yaffe, and G. Czarnota, “Functional imaging using diffuse optical spectroscopy of neoadjuvant chemotherapy response in women with locally advanced breast cancer,” Clin. Cancer Res. 16, 2605–2614 (2010). [CrossRef] [PubMed]

31. Q. Zhu, P. U. Hegde, A. Ricci, M. Kane, E. B. Cronin, Y. Ardeshirpour, C. Xu, A. Aguirre, S. H. Kurtzman, P. J. Deckers, and S. H. Tannenbaum, “Early-stage invasive breast cancers: potential role of optical tomography with US localization in assisting diagnosis,” Radiology 256, 367–378 (2010). [CrossRef] [PubMed]

32. T. Chan and L. Vese “Active contours without edges,” IEEE Trans. Image Process. 10, 266–277 (2001). [CrossRef]

33. H. K. Zhao, S. Osher, B. Merriman, and M. Kang, “Implicit, nonparametric shape reconstruction from unorganized points using a variational level set method,” Comput. Vision Image Understanding 80, 295–314 (2000). [CrossRef]

34. S. Osher and R. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces (Springer, 2002)

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

36. K. Madsen, H. Bruun, and O. Tingleff “Methods for non-linear least squares problems,” lecture notes (2004).

37. C. R. Vogel, Computational Methods for Inverse Problems (SIAM, 2002) [CrossRef]

38. T. T. Wu and K. Lange, “Coordinate descent algorithms for lasso penalized regression,” Ann. Appl. Stat. 2, 224–244 (2008). [CrossRef]

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

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

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

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

43. 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,” in IEEE International Symposium on Biomedical Imaging: from Nano to Macro, 2009. ISBI ’09 (IEEE, 2009), pp. 366–369 (2009).

44. J. D. Vylder and W. Philips, “A computational efficient external energy for active contour segmentation using edge propagation,” in IEEE 2100 International Conference on Image Processing (ICIP 2010) (IEEE, 2010), pp. 661–664.

45. 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 (2002). [CrossRef] [PubMed]

46. 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 optical mammography,” J. Biomed. Opt. 9, 1137–1142 (2004). [CrossRef] [PubMed]

47. J. P. Culver, R. Choe, M. J. Holboke, L. Zubkov, T. Durduran, A. Slemp, V. Ntziachristos, B. Chance, 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 (2003). [CrossRef] [PubMed]

48. 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, 051602 (2007). [CrossRef] [PubMed]

49. 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. Lond. B. Biol. Sci. 352, 727–735 (1997). [CrossRef] [PubMed]

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

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

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

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

54. 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, 054037 (2008). [CrossRef] [PubMed]

Pixel based method
# λ	MSE HbO₂	MSE HbR	MSE Lipid	MSE H₂O	MSE D
8	0.075	0.030	0.048	0.010	0.052
176	0.062	0.021	0.034	0.015	0.030
PaLS method
# λ	MSE HbO₂	MSE HbR	MSE Lipid	MSE H₂O	MSE D
8	0.07	0.03	0.12	0.06	0.08
176	0.02	0.008	0.01	0.02	0.01

Pixel based method
# λ	D(S,G) HbO₂	D(S,G) HbR	D(S,G) Lipid	D(S,G) H₂O	D(S,G) D
8	0.12	0.088	0.089	0.65	0.8
176	0.554	0.1085	0.043	0.41	0.09

Pixel based method
	D(S,G) Set 1		D(S,G) Set 2
# λ	Ink	Dye	Ink	Dye
6	0.143	0.113	0.139	0.145
126	0.142	0.114	0.145	0.140

	PaLS method
# λ	D(S,G) Set 1	D(S,G) Set 2
6	0.27	0.33
126	0.37	0.80

Experimental set 1, 10% ink and 90% dye
Fig.	# λ	Species	$\hat{c_{i}} [%]$	$\hat{c_{i}^{r}} [%]$	MSE
7(a)	6	Ink	1	4	1.8
7(a)	6	Dye	27	96	1.3
7(c)	126	Ink	17	16	2.8
7(c)	126	Dye	88	84	1.2

Parametric level set reconstruction methods for hyperspectral diffuse optical tomography

Abstract

1. Introduction

2. Forward problem

3. Parametric level-set method

4. Inverse problem

5. Simulation analysis

6. Experimental analysis

7. Results

7.1. Simulations

7.2. Experimental validation

8. Conclusion

Acknowledgments

References and links

Cited By

Figures (8)

Tables (6)

Equations (27)

Biomedical Optics Express