Effective contrast recovery in rapid dynamic near-infrared diffuse optical tomography using ℓ1-norm-based linear image reconstruction method

Calvin B. Shaw; Phaneedra K. Yalavarthy

doi:10.1117/1.JBO.17.8.086009

9 August 2012 Effective contrast recovery in rapid dynamic near-infrared diffuse optical tomography using ℓ₁-norm-based linear image reconstruction method

Calvin B. Shaw, Phaneedra K. Yalavarthy

Author Affiliations +

Journal of Biomedical Optics, Vol. 17, Issue 8, 086009 (August 2012). https://doi.org/10.1117/1.JBO.17.8.086009

Abstract

Traditional image reconstruction methods in rapid dynamic diffuse optical tomography employ ℓ2-norm-based regularization, which is known to remove the high-frequency components in the reconstructed images and make them appear smooth. The contrast recovery in these type of methods is typically dependent on the iterative nature of method employed, where the nonlinear iterative technique is known to perform better in comparison to linear techniques (noniterative) with a caveat that nonlinear techniques are computationally complex. Assuming that there is a linear dependency of solution between successive frames resulted in a linear inverse problem. This new framework with the combination of ℓ1-norm-based regularization can provide better robustness to noise and provide better contrast recovery compared to conventional ℓ2-based techniques. Moreover, it is shown that the proposed ℓ1-based technique is computationally efficient compared to its counterpart (ℓ2-based one). The proposed framework requires a reasonably close estimate of the actual solution for the initial frame, and any suboptimal estimate leads to erroneous reconstruction results for the subsequent frames.

1. Introduction

Near-infrared (NIR) diffuse optical tomography is an emerging imaging modality with applications including breast cancer imaging and brain function assays.¹^–³ The interrogating medium in diffuse optical tomography is NIR light in the spectral range of 600 to 1000 nm. A finite set of boundary measurements are made in NIR tomography that is in turn used to reconstruct the internal distribution of optical properties.² The NIR light is typically delivered through optical fibers and the transmitted light is also collected through the same fibers which are in contact with the external surface of the tissue. The distributions of wavelength-dependent absorption and/or scattering coefficients of the tissue are reconstructed using model-based iterative algorithms that in turn use measured boundary data.² As NIR studies have the advantage of being noninvasive and nonionizing, they are suitable for investigating functional changes in tissue over a prolonged time.³

As the tumor vasculature provides a mechanism to distinguish malignant from benign, NIR imaging has the capability to capture this difference in vasculature due to its high sensitivity to hemoglobin and water.³^,⁴ The rapid dynamic imaging of tumor vasculature using NIR light will also enable characterization of the disease, making diffuse optical imaging highly desirable in the clinic. Dynamic diffuse optical imaging is capable of providing images typically at video-rate ( $> 15 frames / s$ ) to study changes in hemodynamics of the tissue under investigation.⁵^–¹⁰ Even though the data is collected at video-rate, the image reconstruction is typically handled off-line.⁵^–⁸ One of the main bottlenecks to get images at video-rate is the computational complexity of nonlinear image reconstruction techniques,⁵^,⁸ leading to investigations such as linear iterative techniques and singular-value decomposition (SVD) methods.⁹ The best performance for these linear techniques is achieved when the initial guess to the image reconstruction procedure is close to the actual solution.⁹^,¹⁰ Also, the linear techniques that were presented in the literature do not account for the correlation between the frames.⁹^,¹⁰ In this work, a new framework for the linear image reconstruction procedure is presented that takes into account this correlation. The main hypothesis for this framework arises from the fact that the main difference between the two successive dynamic optical images is only in the tumor vasculature, and as tumors are highly localized, the difference in the optical images can be reconstructed effectively. This hypothesis leads to a linearized framework that could be solved using $ℓ_{1}$ -based methods, which are known to promote sparseness and allow sharp changes in the optical property distribution. Using both numerical and experimental studies, it is shown that linearized framework in combination with $ℓ_{1}$ -based technique provides effective quantitation of optical properties for the tissue under investigation compared to traditional $ℓ_{2}$ -based techniques. In the presented work, the discussion has been limited to two dimensions, as the main emphasis is on providing a linearized $ℓ_{1}$ -based framework for the dynamic optical image reconstruction.

2. Dynamic NIR Diffuse Optical Tomography: Forward Problem

Video-rate diffuse optical imaging involves collection of the continuous-wave (intensity alone) data at the boundary of the tissue under investigation.⁵^,⁸ Continuous-wave NIR light propagation in thick biological tissues like breast and brain can be modeled using diffusion equation (DE),²^,¹¹ given as

Eq. (1)

- \nabla \cdot D (r) \nabla Φ (r) + μ_{a} (r) Φ (r) = Q_{o} (r),

where the optical diffusion and absorption coefficients are given by

D (r)

and

μ_{a} (r)

, respectively. The continuous-wave light source, represented by

Q_{o} (r)

, is modeled as isotropic.

Φ (r)

is the photon fluence density at a given position

r

. The diffusion coefficient is defined as

Eq. (2)

D (r) = \frac{1}{3 [μ_{a} (r) + μ_{s}^{'} (r)]},

where

μ_{s}^{'} (r)

is the reduced optical scattering coefficient, which is defined as

μ_{s}^{'} = μ_{s} (1 - g)

with

μ_{s}

as the optical scattering coefficient and

g

as the anisotropy factor. In the present work,

μ_{s}^{'}

is assumed to be known and remains constant throughout the domain. The finite element method (FEM) is used to solve Eq. (1) to generate modeled data for a given distribution of the absorption coefficient

μ_{a} (r)

. A Type III boundary condition is employed to account for the refractive-index mismatch at the boundary.¹² Under the Rytov approximation, the modeled data becomes the natural logarithm of the intensity (A),

G (μ) = \ln (A)

, where the forward model is represented by

G (μ)

and

μ

is the spatially varying

μ_{a}

here. This forward model is used repeatedly in an iterative manner to estimate the optical property of the tissue under investigation.²

3. Dynamic NIR Diffuse Optical Tomography: Inverse Problem

The inverse problem primarily involves the estimation of optical absorption coefficients from the CW boundary measurements [ $\ln (A)$ ] using a model-based approach. This is achieved by matching the experimental measurements with model-based ones iteratively in the least-squares sense over the range of $μ_{a}$ . This minimization problem can be solved using several approaches, the most common one involving computing of repeated solutions of the forward model [including Jacobian ( $J$ )] and solving linear system of equations.¹³

In the present work, the inverse problem needs to be solved for a given set of measurements (also known as frames) that were acquired in a dynamic fashion. Most commonly used techniques estimate distribution of $μ_{a}$ using a set of data acquired at time point $t$ , typically resulting in a nonlinear inverse problem. This type of estimation does not account for correlation among the frames (i.e., between time $t$ and $t + δ t$ , with $δ t$ representing the time step) and independently estimates $μ_{a}$ .⁹^,¹⁰^,¹⁴ In this work, a new framework that assumes that $μ_{a}$ distributions estimated at $t$ and $t + δ t$ are linearly dependent will be presented, which leads to linearization of the estimation problem, resulting in a linear inverse problem.

3.1.

Linearization

As the assumption of two successive frames being linearly dependent is taken into account, the aim of the minimization problem here becomes minimizing the difference in the data–model misfit between successive frames (assuming that the previous frame solution is available). This leads to an objective function defined as

Eq. (3)

Ω = \min_{μ_{2}} {{‖ δ_{2} (μ_{2}) - δ_{1} ‖}_{2}^{2}},

where

δ_{2} (μ_{2})

is the data–model misfit of the second frame,

δ_{2} (μ_{2}) = y_{2} - G (μ_{2})

, and

δ_{1}

is the data–model misfit of the first frame,

δ_{1} = y_{1} - G (μ_{1})

.

y_{2}

and

y_{1}

are the natural logarithms of the amplitude of experimental data of second and first frame, respectively.

μ_{2}

and

μ_{1}

represent the solution of the second and first frame, respectively, and

G (μ_{2})

and

G (μ_{1})

represent the modeled response of the second and first frame, respectively.

As stated earlier, in this framework it is assumed that $μ_{1}$ is available and known, with the unknown being $μ_{2}$ . Equation (3) needs to be minimized over the range of $μ_{2}$ . Note that $μ_{1}$ (corresponding to first frame) could be obtained by any standard estimation technique, here, Levenberg-Marquardt (LM) minimization scheme was utilized.

Using the above information and rewriting Eq. (3) leads to

Eq. (4)

Ω = \min_{μ_{2}} {{‖ y_{2} - G (μ_{2}) - [y_{1} - G (μ_{1})] ‖}_{2}^{2}} .

As

μ_{2}

is the only unknown here, expanding

G (μ_{2})

using Taylor series around

μ_{1}

gives

Eq. (5)

G (μ_{2}) = G (μ_{1}) + G^{'} (μ_{1}) (μ_{2} - μ_{1}) + {(μ_{2} - μ_{1})}^{T} G^{''} (μ_{1}) (μ_{2} - μ_{1}) + \dots,

where

G^{'} (μ_{1}) = {J |}_{μ_{1}}

is the Jacobian [dimension: NMXNN with NM representing number of measurements and NN number of imaging parameters (nodes in the FEM mesh)] evaluated at

μ_{1}

and

G^{''} (μ_{1})

is the Hessian.

Neglecting the higher-order terms, equivalently linearizing the expansion gives

Eq. (6)

G (μ_{2}) ≃ G (μ_{1}) + {J |}_{μ_{1}} (μ_{2} - μ_{1}) .

Using Eq. (6) in Eq. (4) leads to

Eq. (7)

Ω = \min_{Δ μ_{2}} {{‖ Δ y_{2} - {J |}_{μ_{1}} Δ μ_{2} ‖}_{2}^{2}},

where

Δ μ_{2} = μ_{2} - μ_{1}

and

Δ y_{2} = y_{2} - y_{1}

. Now, the solution to the second frame is given by

Eq. (8)

μ_{2} = μ_{1} + Δ μ_{2} .

Note that Eq. (3) represents a nonlinear inverse problem; assuming that

μ_{2}

and

μ_{1}

are linearly dependent [Eq. (8)] makes it a linear inverse problem [Eq. (7)]. Also, the inverse problem is still ill-posed primarily due to the condition number of Jacobian

J

and requires regularization to solve the problem. The linear inverse problem given by Eq. (7) needs to be minimized over the range of

Δ μ_{2}

, whereas the original problem [Eq. (3)] needs to minimized with respect to

μ_{2}

.

Generalizing Eq. (7) for the $n$ th frame, Eq. (9) can be written as

Eq. (9)

Ω = \min_{Δ μ_{n}} {{‖ Δ y_{n} - {J |}_{μ_{n - 1}} Δ μ_{n} ‖}_{2}^{2}},

where

Δ μ_{n} = μ_{n} - μ_{n - 1}

and

Δ y_{n} = y_{n} - y_{n - 1}

. Now, the solution to the

n

th frame is given by

Eq. (10)

μ_{n} = μ_{n - 1} + Δ μ_{n} .

3.2.

$ℓ_{1}$ -Based Linear Reconstruction Method

As this new framework results in a linear inverse problem, it could be solved using $ℓ_{1}$ type of approaches under the assumption that the update $Δ μ_{n}$ need not be a smooth function. This is a valid assumption even in the physiological sense, as tumors tend to be highly localized and the appreciable change in optical properties between successive frames is pulsatile.⁴^,¹⁰^,¹⁵

Even though the regularization term to be added for solving Eq. (9) can take many forms, the $ℓ_{1}$ -based regularization is explored in this work. Under this framework, the regularization term that is added to the objective function is based on $ℓ_{1}$ -norm and in this work it is posed as a basis pursuit denoising (BPDN) problem,¹⁶^,¹⁷ resulting in an additional term to be added to Eq. (9), given by

Eq. (11)

Ω = \min_{Δ μ_{n}} {λ {‖ Δ μ_{n} ‖}_{1} + {‖ Δ y_{n} - {J |}_{μ_{n - 1}} Δ μ_{n} ‖}_{2}^{2}},

where

{‖ Δ μ_{n} ‖}_{1}

denotes the

ℓ_{1}

-norm of

Δ μ_{n}

, i.e.,

{‖ Δ μ_{n} ‖}_{1} = \sum_{i} | Δ μ_{n} |

, and

λ

is the regularization parameter which should be selected (optimized) such that the error between the residual and the image is minimized. The solution obtained with this objective function is substituted in Eq. (10) to get the solution to the

n

th frame. Note that

ℓ_{1}

-framework is shown to provide better noise-tolerance characteristics as well as better convergence properties and promotes sparseness in the solution in other linear estimation problems;¹⁶^,¹⁷ discussion of the same is beyond the scope of this work.

Minimization of Eq. (11) to find solution $Δ μ_{n}$ can be achieved using any of the greedy or basis pursuit algorithms. In this work, we have used YALL1, which is open source and known to provide a lot of versatility in terms of finding $Δ μ_{n}$ . YALL1¹⁸ is derived from the alternating direction method (ADM), which minimizes augmented Lagrangian functions through an alternating minimization scheme and updates the multipliers after each sweep. The YALL1 solver can solve eight equivalent $ℓ_{1}$ -minimization models, although here only one such model was used (details to follow). The construction of these equivalent models consist of two main steps: 1. reformulate an $ℓ_{1}$ problem into one having partially separable objective functions by adding new variables and constraints; and 2. apply an exact or inexact alternating direction method to the resulting problem. This algorithm can be regarded as a first-order primal dual algorithm because both primal and dual variables are updated at each and every iteration. Among the eight $ℓ_{1}$ -minimization models, the model that has the form

Eq. (12)

Ω = \min_{Δ μ_{n}} {{‖ Δ μ_{n} ‖}_{1} + \frac{1}{2 ρ} {‖ Δ y_{n} - {J |}_{μ_{n - 1}} Δ μ_{n} ‖}_{2}^{2}}

was chosen with

ρ

representing the data-fidelity (also known as regularization), chosen as

10^{- 2}

for the studies performed in this work. It could be easily shown that Eqs. (12) and (11) are equivalent with

λ = 2 ρ

. A detailed account of the numerical implementation employed can be found in Ref. 18. Note that the

ℓ_{1}

-based minimization technique employs an iterative approach to find the optimal solution (

Δ μ_{n}

), even though the inverse problem is linear. The choice of number of iterations that were deployed was based on the noise level in the data. For a typical 1% noise case, the number of iterations was chosen to be 60, and as the noise level increased, the number of iterations were also increased (for the case of 5% noise, it was 65).

3.3.

Regularized Minimal Residual Method: $ℓ_{2}$ -Based

The traditional/conventional method of solving an ill-posed linear inverse problem involves the addition of regularization term that is based on $ℓ_{2}$ -norm. As the $ℓ_{1}$ -norm–based regularized problem was solved using an iterative approach, also here an iterative method was employed. Regularized minimal residual (MinRes) method is a computational alternative for conjugate gradient (CG) method which minimizes the $ℓ_{2}$ -norm based optimization problem.¹⁹ In this the objective function has an additional term based on $ℓ_{2}$ -norm [similar to Eq. (11)] and is given by

Eq. (13)

Ω = \min_{Δ μ_{n}} {α {‖ Δ μ_{n} ‖}_{2}^{2} + {‖ Δ y_{n} - {J |}_{μ_{n - 1}} Δ μ_{n} ‖}_{2}^{2}},

where

α

is the regularization parameter. The solution obtained with this objective function is substituted in Eq. (10) to get the solution to the

n

th frame. A computational scheme of the regularized MinRes method¹⁹ for minimizing the objective function [Eq. (13)] is given in Algorithm 1.

Algorithm 1

Regularized MinRes method.

Reconstruction of

μ_{a}

corresponding to frame

n

INPUT:

{J |}_{μ_{n - 1}}

and

Δ y_{n}

; OUTPUT:

Δ μ_{n}

.

Initialize

Δ μ_{n}^{0}

(initial guess),

α

.

for

i = 0, 1, \dots

(representing inner iteration number)

1.

r^{i} = {J |}_{μ_{n - 1}} {Δ μ_{n}}^{i} - Δ y_{n}

2.

l_{α}^{i} = l_{α} (Δ μ_{n}^{i}) = {{J |}_{μ_{n - 1}}}^{T} r^{i} + α Δ μ_{n}^{i}

3.

k_{α}^{i} = \frac{{‖ l_{α} (Δ μ_{n}^{i}) ‖}^{2}}{{‖ {J |}_{μ_{n - 1}} l_{α} (Δ μ_{n}^{i}) ‖}^{2} + α {‖ l_{α} (Δ μ_{n}^{i}) ‖}^{2}}

4. Update equation:

Δ μ_{n}^{i + 1} = Δ μ_{n}^{i} - k_{α}^{i} l_{α}^{i}

5. The iterative process, steps: 1 to 4, is terminated when the misfit reaches the given stopping criterion

(ϵ_{o}) : {‖ r^{i} ‖}^{2} \leq ϵ_{o}

The MinRes method primarily solves the minimization function given by Eq. (13) for a given $α$ and an initial guess for $Δ μ_{n}$ using an iterative approach, without employing any explicit inverse.¹⁹ In this method, for increasing computational efficiency only matrix-vector multiplications are performed. As the name suggests, the aim is to find a solution that gives the minimal residual in an iterative manner. Initially, a guess for $Δ μ_{n}$ (with $n$ representing the frame number) is given as input to the algorithm ( $Δ μ_{n}^{0}$ ), typically chosen to be constant vector (0.001 here). A residual ( $r$ ) is calculated at every step (as given in step 1 of Algorithm 1) with the help of Jacobian (J) and data-model misfit ( $Δ y_{n}$ ) corresponding to that iteration. Subsequently a positively determined vector ( $l_{α}^{i}$ ) that takes into account regularization ( $α$ ) is computed using the residual vector as well as earlier value of $Δ y_{n}$ (either a guess or corresponding to previous iteration). The $l_{α}^{i}$ is also known as the residual of the Euler equation, which is the direct solution to Eq. (13) given by Ref. 13

Eq. (14)

[{{J |}_{μ_{n - 1}}}^{T} * {J |}_{μ_{n - 1}} + α I] Δ μ_{n}^{i} = {{J |}_{μ_{n - 1}}}^{T} Δ y_{n} .

Following this, the step length based on the residual and Jacobian is determined (step 3 of Algorithm 1) before proceeding to the update equation. This process is repeated until the residual (calculated in step 1) becomes small (less than or equal to

ϵ_{o}

). In this work,

ϵ_{o} = 10^{- 4}

and

α = 100

were chosen for all studies performed here, as they provided the best estimates (also shown later). This algorithm is equivalent to the regularized steepest descent method that is applicable to linear inverse problem.

In this work, Jacobian ( $J$ ) corresponding to the initial frame was used for reconstructing rest of the frames in both $ℓ_{1}$ - and $ℓ_{2}$ -based iterative reconstruction schemes, removing the computationally expensive step of computing Jacobian ( $J$ ) for every frame.²⁰ The reconstructions were performed on a Linux workstation that had a 2.4 GHz Intel Quadcore processor along with 8 GB RAM. The details of simulation and experimental studies performed as part of this work are presented in the next section.

4. Simulation and Experimental Evaluation

4.1.

Numerical Experimental Data

As the reconstruction problem in rapid dynamic NIR imaging has been reformulated as a linear problem, initially the performance of the same in recovering high-contrast tumors is taken up. Note that in all simulation cases discussed here, the imaging domain is chosen to be circular in shape having a diameter of 86 mm, mimicking breast imaging. The target is also considered to be circular in shape with a radius of 10 mm. The background optical properties, mimicking the breast, are assumed to be $μ_{a} = 0.01 {mm}^{- 1}$ and $μ_{s}^{'} = 1.0 {mm}^{- 1}$ . The target $μ_{a}$ was varied depending on the numerical experiment, but $μ_{s}^{'}$ was fixed at $1.0 {mm}^{- 1}$ (as it is considered a known parameter). The data collection setup had 16 fibers arranged in equi-spaced fashion on the boundary of the imaging domain, where when one fiber acted as a source, and the rest acted as detectors. This resulted in 240 ( $16 \times 15$ ) measurement points.²¹ The sources were modeled as having a Gaussian profile with full width half maximum of 3 mm to mimic the experimental conditions.²¹ The source was placed at one mean transport length inside the boundary.

For conducting the numerical experiments, two finite element meshes were considered, one for experimental data generation, the other in reconstruction scheme. In the case of experimental data generation, the finite element mesh had 10,249 nodes corresponding to 20,160 triangular elements. For the reconstruction scheme, the nodes were 2728 corresponding to 5362 elements. In all cases (including experimental), the first frame is always reconstructed using nonlinear iterative method [Levenberg-Marquardt (LM) minimization scheme¹³] and the dynamic reconstruction starts from the second frame using both linear methods ( $ℓ_{1}$ and $ℓ_{2}$ ) discussed in the earlier sections.

In the numerical experiments, initially a stationary target located at the center of the domain with increasing contrast initially for four frames and decreasing contrast in subsequent frames is considered. The increase/decrease in $μ_{a}$ in the target at every step is $0.005 {mm}^{- 1}$ , mimicking the hemodynamic response in a tumor. The target distribution of $μ_{a}$ along with the frame numbers are given in top row of Fig. 1(a). The numerical experimental data corresponding to each of these frames were generated using 10,249-node mesh and was calibrated for the 2728-node mesh using the standard calibration routines.²² Note that the calibration procedure uses the analytical solution to DE, where the domain is assumed to be infinite/semi-infinite, making the procedure computationally inexpensive.²² Four cases of noise levels (0%, 1%, 3%, and 5%) in the data were considered to mimic the experimental data. Note that the noise level up to 4% was reported in the literature²³ in diffuse optical tomography systems, the 5% noise forms the worst case.

Fig. 1

Comparison of performance of $ℓ_{1}$ - and $ℓ_{2}$ -based reconstructions for the case of stationary target with data having noise levels of 0% (a), 1% (b), 3%, and 5% (d). The target distribution is given as first row of images in (a). The corresponding frame numbers are given on top of each subfigure. The reconstruction method that was employed has been indicated against each row of the subfigures. The first column gives the reconstruction of first frame using $ℓ_{2}$ -based method that was used as prior image for the second frame.

A similar effort is considered for the case where the target (having a fixed contrast of $2 ∶ 1$ ) placed on the $x$ -axis [with center being at $(20, 0)$ ] moving in 11 steps to reach diagonally opposite location $(- 20, 0)$ . The target distribution for this case is given in Fig. 2(a) top row along with their corresponding frame number given on top.

Fig. 2

Similar to Fig. 1 except the target is moving and the expected contrast is fixed at $2 ∶ 1$ .

In all cases considered here, it was assumed that the reduced scattering coefficient of the imaging domain was uniform and constant, which may not be true in the real cases. To study the influence of scattering on the estimated $μ_{a}$ , a case of ${μ_{s}}^{'}$ having a contrast of $2 ∶ 1$ in the target compared to the background is considered. Note that as it is very unlikely to have a significant variation in the scattering coefficient in the dynamic imaging scenario, it was assumed that there is no change in ${μ_{s}}^{'}$ throughout the frames. For this case, the ${μ_{s}}^{'}$ is an inhomogeneous model with target value being at $2 {mm}^{- 1}$ and background at $1 {mm}^{- 1}$ . The reconstructions were performed for the stationary target located at the center [target distribution is shown in top row of Fig. 1(a)] using 1% numerically generated experimental data assuming a homogenous ${μ_{s}}^{'}$ ( $= 1 {mm}^{- 1}$ ) throughout the domain.

4.2.

Experimental Phantom Data

To assess the capabilities of the proposed method, an experimental phantom dataset that was acquired using a video-rate NIR tomography system⁷ was considered. In this case, the system had an imaging array of 27 mm in diameter (mimicking the small-animal imaging, e.g., rat’s cranium) and 16 equi-spaced channels for collecting the NIR data. The solid phantom with a diameter of 27 mm used here also had a hole of 6.35 mm in diameter for holding the intra-lipid solution. A 785-nm-wavelength laser diode was used as a source and the corresponding optical properties were $μ_{a} = 0.0002 {mm}^{- 1}$ and ${μ_{s}}^{'} = 1.45 {mm}^{- 1}$ . The intra-lipid had matching optical properties of the solid phantom and the undiluted India ink was injected shortly after the data acquisition. The data was acquired at 35 frames per second and the ink was injected using a pipet. Once again, here the data was calibrated to remove biases in the data and to match the response of the numerical model.²²

5. Results

Using both $ℓ_{1}$ - and $ℓ_{2}$ -based linear reconstruction techniques, the contrast recovery study results with varying noise level, 0%, 1%, 3%, and 5%, in data in the case of stationary target located at the center of imaging domain are given in Fig. 1(a), 1(b), 1(c), and 1(d), respectively. Here, the contrast recovery refers to the ratio of the mean value of the reconstructed $μ_{a}$ in the target region (mimicking tumor) to the background. The target distributions are given in the top row of Fig. 1(a). The distributions corresponding to frame 1 (target and reconstructed image using $ℓ_{2}$ -based method) are also given in the respective figures in the first column. Figure 1 results indicate that the contrast recovery was better using $ℓ_{1}$ -based method compared to $ℓ_{2}$ -based method. To assess this quantitatively, the contrast recovered by computing the mean value of $μ_{a}$ in the target region is taken up and the same is plotted in Fig. 3 correspondingly. The results indicate that the contrast recovery is poor when $ℓ_{2}$ -based methods are employed compared to $ℓ_{1}$ -based method.

Fig. 3

Comparative plots of reconstructed mean $μ_{a}$ in the target region (ROI) using $ℓ_{1}$ - and $ℓ_{2}$ -based methods along with expected values (target) for results presented in Fig. 1.

For the case of moving target (fixed contrast of $2 ∶ 1$ ) the reconstruction results using 0%, 1%, 3%, and 5% noisy data are given in Fig. 2(a), 2(b), 2(c), and 2(d), respectively. These results also indicate that the performance of $ℓ_{2}$ -based reconstruction technique is inferior in both qualitative and quantitative nature of reconstructed images compared to $ℓ_{1}$ -based technique.

The results pertaining to the case of having an inhomogeneous ${μ_{s}}^{'}$ and reconstructing the $μ_{a}$ distribution using $ℓ_{1}$ - and $ℓ_{2}$ -based methods assuming a homogeneous ${μ_{s}}^{'}$ are given in Fig. 4(a). Similar to Fig. 3, the contrast recovery plot for this case is given in Fig. 4(b). It is amply clear from this result that having a reasonable ${μ_{s}}^{'}$ distribution close to the actual (expected) one is essential for obtaining quantitatively accurate results.

Fig. 4

(a), Reconstructed distributions of $μ_{a}$ in the case of 1% noise level in the data with different optical scattering coefficient for the background ( $μ_{s}^{'} = 1.0 {mm}^{- 1}$ ) and the target ( $μ_{s}^{'} = 2.0 {mm}^{- 1}$ ) using $ℓ_{1}$ - and $ℓ_{2}$ -based methods. The $μ_{a}$ target distribution is same as the one shown in Fig. 1(a). (b), Recovered mean $μ_{a}$ of the target corresponding to results presented in (a).

The reconstructed $μ_{a}$ distribution using experimental phantom data are presented in Fig. 5(a). Note that even though the NIR data was acquired at $35 frames / s$ , the reconstruction results that were presented here are sampled versions of total available distributions and the corresponding time of acquisition is given on top of Fig. 5(a). The recovered contrast in the target region for these set of results [given in Fig. 5(a)] is also plotted in Fig. 5(b) for providing quantitative assessment. From both Fig. 5(a) and 5(b), it is evident that the contrast recovery is superior in case of $ℓ_{1}$ -based reconstruction technique compared to its counterpart ( $ℓ_{2}$ -based method).

Fig. 5

(a) Reconstructed distributions of $μ_{a}$ in the case of experimental phantom data using $ℓ_{1}$ - and $ℓ_{2}$ -based methods. The time corresponding to each frame is given in top row. (b) Maximum recovered $μ_{a}$ contrast corresponding to results presented in (a).

To understand the computational complexity of the presented reconstruction technique, the total reconstruction time, overhead time, and total number of inner iterations (corresponding time) is recorded (averaged over four runs) for the results of Fig. 2(b) and the same is reported in Table 1.

Table 1

Average computational time recorded for reconstruction of single frame using both ℓ1- and ℓ2-based methods with 1% noise when the target is in motion [results corresponding to Fig. 2(b)].

Task	$ℓ_{1}$	$ℓ_{2}$
Total reconstruction timea	0.5075	0.9001
Overhead timea	0.2070	0.185
Number of iterations	60	670
Timea taken per iteration	0.005	0.00106

^aTime in seconds.

6. Discussion

Rapid dynamic NIR tomography has the potential to reveal the hemodynamic response of the tumor, providing functional data that could potentially characterize/reveal the patho-physiological state of the tissue under investigation. But typically, the image reconstruction is handled off-line (not in sync with data acquisition) due to the computational complexity.⁵^–¹⁰^,¹⁵ Developing image reconstruction methods that can improve the imaging speed to be on par with the data acquisition is highly desirable to take advantage of NIR tomography’s potential. Unfortunately, the linear reconstruction techniques that perform single-step reconstruction can reduce the computational complexity, but are known to have limited capability in terms of contrast recovery,⁵^,⁷^,⁹ which is needed for an accurate estimation of the patho-physiological state of the tissue. $ℓ_{1}$ -Based techniques are known to promote sparse solutions and have been reported to provide better spatial resolution in case of diffuse optical imaging.²⁴^–²⁷ In this work, a new framework that can make use of $ℓ_{1}$ -based regularization in dynamic diffuse optical imaging is developed, and the developed framework also made the reconstruction problem to be linear. Even though it is possible that the solution within this frame work could be sparse, this work did not explicitly take this into account as it was aimed at showing that contrast recovery, in cases resulting in pulsatile absorption,¹⁰^,¹⁵ could be better using $ℓ_{1}$ -based techniques compared to traditional methods ( $ℓ_{2}$ -based). There could be scenarios, especially in small-animal imaging, where the changes in optical properties need not be localized or in the tumor vasculature alone. In this cases, a sparsifying transform (wavelet or Fourier-based) could be employed, which can make the solution sparse to deploy $ℓ_{1}$ -based optimization schemes. This is well-studied for static diffuse optical tomographic case,¹⁴ but requires additional computation for such a transformation.

The hemodynamic changes between the successive frames in rapid dynamic NIR tomography are predominantly in tumors alone, leading to changes in optical properties to be highly localized, resulting in a linear solution between the successive frames. Untill now the reconstruction algorithms that were investigated in the literature for dynamic diffuse optical imaging did not account for the linear dependence.⁵^–¹⁰^,¹⁵ In the $ℓ_{1}$ -based reconstruction method presented here, alternating direction method was deployed, making use of an iterative technique to solve the optimization problem. These iterations are not global in nature and do not require recomputation of either Jacobian ( $J$ ) or forward data. In spirit, they are similar to inner iterations applied in gradient-based methods that were employed earlier in NIR tomography.²^,³ As the $ℓ_{1}$ -based method employed iterative technique to solve the linear inverse problem, an iterative technique based on MinRes method (similar to steepest descent method) was introduced in this work in case of the $ℓ_{2}$ -based method for an effective comparison. The solution obtained from MinRes is the same as the one obtained by direct method [given by Eq. (14)], except that MinRes has the advantage that the order of computation here is $O (n^{2})$ in comparison to $O (n^{3})$ ¹⁹ in case of Eq. (14), with $n$ being the number of columns in Jacobian ( $J$ ). Note that in both $ℓ_{1}$ and $ℓ_{2}$ cases, an unconstrained optimization was performed as the update ( $Δ μ_{n}$ ) could be negative or positive.

The $ℓ_{1}$ -based methods are known to improve the spatial resolution (avoiding over smoothing) as well as depth localization in diffuse optical imaging.²⁴^,¹⁴^,²⁷ Recent works also indicated that they are capable of providing more robustness to noise.¹⁴. The results (Figs. 1 Fig. 3 Fig. 2 Fig. 4–5) obtained using $ℓ_{1}$ -based method in the linear framework for rapid dynamic NIR tomography asserted the same in comparison to $ℓ_{2}$ -based methods. In case of quantitation, the $ℓ_{1}$ -based method is atleast 70% better in terms of error in the contrast recovery compared to $ℓ_{2}$ -based methods (Fig. 3). Also, as the noise level increased, the performance of $ℓ_{2}$ -based method became inferior (Figs. 2 and 3) in comparison to $ℓ_{1}$ -based method.

In case of experimental phantom data, even though the expected contrast was close to infinity (as the ink was undiluted), the $ℓ_{1}$ -based method was at least six times superior in terms of recovered contrast [Fig. 5(b)], making $ℓ_{1}$ -based technique highly desirable in the experimental cases.

It should be noted that the linear inverse problem framework obtained in this work has a limitation that the obtained reconstruction of initial frame (Frame 1) needs to be close to the target/expected distribution, otherwise the quantitation errors will propagate. Having an inaccurate ${μ_{s}}^{'}$ distribution used in the NIR light propagation models induces errors in the forward data calculations. In turn, these errors reflect in the reconstructed $μ_{a}$ values, as it is the only parameter that is assumed to be unknown and allowed to change. The same is observed in Fig. 4. As the NIR light is attenuated by both absorption and scattering, having an inaccurate ${μ_{s}}^{'}$ distribution in the forward model will lead to estimated $μ_{a}$ being higher then the expected value [Fig. 4(b)]. To uniquely estimate the $μ_{a}$ and ${μ_{s}}^{'}$ one requires either frequency-domain or time-domain measurements rather than pure CW measurements (intensity data). Obtaining these frequency or time-domain measurements requires additional resources/instruments and should be very easy to incorporate in the video-rate data acquisition systems. As scattering is unlikely to change with the frames, it is reasonable to assume that the additional data other than pure intensity needs to be obtained only for the initial frame (Frame 1). This could ensure that the errors do not propagate in the time. Interest in most dynamic video-rate optical imaging applications lies with the relative difference in the absorption properties between the frames. Having an inaccurate ${μ_{s}}^{'}$ might change the base value, but the observed relative difference is similar to that of the target in the results shown in Fig. 4. It also could be observed that the inaccurate modeling of ${μ_{s}}^{'}$ only added the bias to the contrast recovered [Fig. 4(b)], with the $ℓ_{1}$ -based method doing better than the $ℓ_{2}$ -based one.

The choice of regularization ( $λ$ in Eq. (11) for $ℓ_{1}$ -based method and $α$ in Eq. (13) for $ℓ_{2}$ -based method) was based on root mean square (RMS) error in the reconstructed $μ_{a}$ distribution. The plot showing this is given in Fig. 6(a) and 6(b) for the choice of $λ$ and $α$ , respectively, showing that lowest RMS errors were obtained for $λ = 0.02$ (in turn $ρ = 10^{- 2}$ ) and $α = 100$ . Note that even though this is shown only for one case [fourth frame of Fig. 1(b)], the same trend corresponding to other frames is observed to confirm that the chosen values of regularization leads to the least amount of errors.

Fig. 6

Plot showing the root mean square (RMS) error in the reconstructed $μ_{a}$ (region of interest) as a function of regularization parameter for $ℓ_{1}$ -based (a) and $ℓ_{2}$ -based (b) methods corresponding to fourth frame as shown in Fig. 1(b).

Note that the total variation regularization will also obey the $ℓ_{1}$ -norm, but is known to be computationally complex in providing a solution to the minimization problem compared to $ℓ_{2}$ -norm based methods.²⁸^,²⁹ Here, the alternating direction method (ADM) that minimizes the $ℓ_{1}$ -norm-based scheme was directly employed through the open-source YALL1.¹⁸ The computational time that was recorded in $ℓ_{1}$ - and $ℓ_{2}$ -based methods (Table 1) also showed that the $ℓ_{1}$ -based method has a distinct advantage not only in providing better quantitation, but also in computational complexity. Also, in the $ℓ_{1}$ -based method the computational time for actual reconstruction (for completing the inner iterations) is on par with the overhead time (Table 1). In the $ℓ_{2}$ -based method, the number of inner iterations is larger compared to its counterpart. Even though the computational time recorded was only reported for the results obtained in Fig. 3, the trends that were observed in general are true for other cases in this work. Moreover, these computations can be accelerated further by using any high-performance computing environments.²⁰

The recent work by Dutta et al.³⁰ in the context of fluorescence molecular tomography (FMT), which has a similar physics to diffuse optical imaging, has shown that $ℓ_{1}$ reconstruction has out performed the $ℓ_{2}$ -based methods in both the reconstructed target (ROI) value and the obtained contrast (signal to background), reconfirming that the trends observed here match with the literature. It is also shown in Ref. 30 that a joint $ℓ_{1}$ and total variation (TV) regularization for the FMT not only provided better contrast in comparison to standard $ℓ_{2}$ regularization, but also reduced the RMS errors in the background, which were higher when used individually. The detailed study in similar lines in the current context will be taken up as a future work.

7. Conclusions

Traditional $ℓ_{2}$ -based techniques are known to provide poor contrast recovery in case of linear reconstruction methods employed in rapid dynamic diffuse optical imaging. In this work, the dynamic diffuse optical image reconstruction problem is reformulated as a linear problem taking into account the linear dependence among the frames. This formulation naturally led to effective usage of $ℓ_{1}$ -based techniques to estimate the corresponding optical distributions for these under-determined linear problems. The $ℓ_{1}$ -based method that was deployed here used alternating direction method (ADM), which is an iterative method. These solutions were compared and contrasted against $ℓ_{2}$ -based methods (iterative in nature) and showed that the $ℓ_{1}$ -based method can provide better quantitation in these dynamic studies and also is more robust to noise. Moreover, these $ℓ_{1}$ -based methods have lesser computational complexity compared to their counterparts ( $ℓ_{2}$ -based). These types of reconstruction techniques that provide better quantification as well as being computationally inexpensive are highly desirable for better characterization of tissue hemodynamics.

Acknowledgments

The authors are thankful to Dr. Daqing Piao for providing the experimental phantom data that was used in this work. The authors also thank Dr. Chandra Murthy and Dr. Namrata Vaswani for their initial discussions on sparse estimation techniques. This work is supported by the Department of Atomic Energy Young Scientist Research Award (No. 2010/20/34/6/BRNS) by Government of India.

References

1.

D. A. Boaset al., “Imaging the body with diffuse optical tomography,” IEEE Signal Process. Mag., 18 (6), 57 –75 (2001). http://dx.doi.org/10.1109/79.962278 ISPRE6 1053-5888 Google Scholar

2.

S. R. ArridgeJ. C. Schotland, “Optical tomography: forward and inverse problems,” Inv. Problems, 25 (12), 123010 (2009). http://dx.doi.org/10.1088/0266-5611/25/12/123010 INPEEY 0266-5611 Google Scholar

3.

T. Durduranet al., “Diffuse optics for tissue monitoring and tomography,” Rep. Prog. Phys., 73 (7), 076701 (2010). http://dx.doi.org/10.1088/0034-4885/73/7/076701 RPPHAG 0034-4885 Google Scholar

4.

B. W. Pogueet al., “Implicit and explicit prior information in near-infrared spectral imaging: accuracy, quantification and diagnostic value,” Phil. Trans. R. Soc. A, 369 (1955), 4531 –4557 (2011). http://dx.doi.org/10.1098/rsta.2011.0228 PTRMAD 1364-503X Google Scholar

5.

C. H. Schmitzet al., “Instrumentation for fast functional optical tomography,” Rev.Sci. Instrum., 73 (2), 429 –439 (2002). http://dx.doi.org/10.1063/1.1427768 RSINAK 0034-6748 Google Scholar

6.

C. H. Schmitzet al., “Dynamic studies of small animals with a four-color diffuse optical tomography imager,” Rev. Sci. Instrum., 76 (9), 094302 (2005). http://dx.doi.org/10.1063/1.2038467 RSINAK 0034-6748 Google Scholar

7.

D. Piaoet al., “Instrumentation for video-rate near-infrared diffuse optical tomography,” Rev.Sci. Instrum., 76 (12), 124301 (2005). http://dx.doi.org/10.1063/1.2149147 RSINAK 0034-6748 Google Scholar

8.

D. Piaoet al., “Video-rate near-infrared optical tomography using spectrally encoded parallel light delivery,” Opt. Lett., 30 (19), 2593 –2595 (2005). http://dx.doi.org/10.1364/OL.30.002593 OPLEDP 0146-9592 Google Scholar

9.

S. Guptaet al., “Singular value decomposition based computationally efficient algorithm for rapid dynamic near-infrared diffuse optical tomography,” Med. Phys., 36 (12), 5559 –5567 (2009). http://dx.doi.org/10.1118/1.3261029 MPHYA6 0094-2405 Google Scholar

10.

Z. Liet al., “Video-rate near infrared tomography to image pulsatile absorption properties in thick tissue,” Opt. Express, 17 (14), 12043 –12056 (2009). http://dx.doi.org/10.1364/OE.17.012043 OPEXFF 1094-4087 Google Scholar

11.

H. Jianget al., “Optical image reconstruction using frequency domain data: simulations and experiments,” J. Opt. Soc. Am. A, 13 (2), 253 –266 (1996). http://dx.doi.org/10.1364/JOSAA.13.000253 JOAOD6 0740-3232 Google Scholar

12.

M. Schweigeret al., “The finite element model for the propagation of light in scattering media: boundary and source conditions,” Med. Phys., 22 (11), 1779 –1792 (1995). http://dx.doi.org/10.1118/1.597634 MPHYA6 0094-2405 Google Scholar

13.

P. K. Yalavarthyet al., “Weight-matrix structured regularization provides optimal generalized least-squares estimate in diffuse optical tomography,” Med. Phys., 34 (6), 2085 –2098 (2007). http://dx.doi.org/10.1118/1.2733803 MPHYA6 0094-2405 Google Scholar

14.

M. SuzenA. GiannoulaT. Durduran, “Compressed sensing in diffuse optical tomography,” Opt. Express, 18 (23), 23676 –23690 (2010). http://dx.doi.org/10.1364/OE.18.023676 OPEXFF 1094-4087 Google Scholar

15.

Z. Liet al., “Rapid magnetic resonance-guided near-infrared mapping to image pulsatile hemoglobin in the breast,” Opt. Lett., 35 (23), 3964 –3966 (2010). http://dx.doi.org/10.1364/OL.35.003964 OPLEDP 0146-9592 Google Scholar

16.

E. CandesM. Wakin, “An introduction to compressive sampling,” IEEE Signal Process. Magaz., 25 (2), 21 –30 (2008). http://dx.doi.org/10.1109/MSP.2007.914731 ISPRE6 1053-5888 Google Scholar

17.

J. Romberg, “Imaging via compressive sampling,” IEEE Signal Process. Magaz., 25 (2), 14 –20 (2008). http://dx.doi.org/10.1109/MSP.2007.914729 ISPRE6 1053-5888 Google Scholar

18.

J. YangY. Zhang, “Alternating direction algorithms for L1-problems in compressive sensing,” SIAM J. Sci. Comput., 33 (1), 250 –278 (2011). Google Scholar

19.

M. S. Zhdanov, Geophysical Inverse Theory and Regularization Problems, 1st ed.Elsevier Science, New York (2002). Google Scholar

20.

J. Prakashet al., “Accelerating frequency-domain diffuse optical tomographic image reconstruction using graphics processing units,” J. Biomed. Opt., 15 (6), 066009 (2010). http://dx.doi.org/10.1117/1.3506216 JBOPFO 1083-3668 Google Scholar

21.

T. O. Mcbrideet al., “A parallel-detection frequency-domain near-infrared tomography system for hemoglobin imaging of the breast in vivo,” Rev. Sci. Instrum., 72 (3), 1817 –1824 (2001). http://dx.doi.org/10.1063/1.1344180 RSINAK 0034-6748 Google Scholar

22.

B. W. Pogueet al., “Calibration of near infrared frequency-domain tissue spectroscopy for absolute absorption coefficient quantitation in neonatal head-simulating phantoms,” J. Biomed. Opt., 5 (2), 185 –193 (2000). http://dx.doi.org/10.1117/1.429985 JBOPFO 1083-3668 Google Scholar

23.

G. Gulsenet al., “Design and implementation of a multifrequency near-infrared diffuse optical tomography system,” J. Biomed. Opt., 11 (1), 014020 (2006). http://dx.doi.org/10.1117/1.2161199 JBOPFO 1083-3668 Google Scholar

24.

N. CaoA. NehoraiM. Jacobs, “Image reconstruction for diffuse optical tomography using sparsity regularization and expectation-maximization algorithm,” Opt. Express, 15 (21), 13695 –13708 (2007). http://dx.doi.org/10.1364/OE.15.013695 OPEXFF 1094-4087 Google Scholar

25.

J. C. Baritauxet al., “Sparsity-driven reconstruction for FDOT with anatomical priors,” IEEE Trans. Med. Imag., 30 (5), 1143 –1153 (2011). http://dx.doi.org/10.1109/TMI.2011.2136438 ITMID4 0278-0062 Google Scholar

26.

O. Leeet al., “Compressive diffuse optical tomography: noniterative exact reconstruction using joint sparsity,” IEEE Trans. Med. Imag., 30 (5), 1129 –1142 (2011). http://dx.doi.org/10.1109/TMI.2011.2125983 ITMID4 0278-0062 Google Scholar

27.

V. C. Kavuriet al., “Sparsity enhanced spatial resolution and depth localization in diffuse optical tomography,” Biomed. Opt. Express, 3 (5), 943 –957 (2012). http://dx.doi.org/10.1364/BOE.3.000943 BOEICL 2156-7085 Google Scholar

28.

K. D. PaulsenH. Jiang, “Enhanced frequency-domain optical image reconstrution in tissues through total-variation minimization,” Appl. Opt., 35 (19), 3447 –3458 (1996). http://dx.doi.org/10.1364/AO.35.003447 APOPAI 0003-6935 Google Scholar

29.

A. Borsicet al., “In vivo impedance imaging with total variation regularization,” IEEE Trans. Med. Imaging, 29 (1), 44 –54 (2010). http://dx.doi.org/10.1109/TMI.2009.2022540 ITMID4 0278-0062 Google Scholar

30.

J. Duttaet al., “Joint

l^{1}

and total variation regularization for fluorescence molecular tomography,” Phys. Med. Biol., 57 (6), 1459 –1476 (2012). http://dx.doi.org/10.1088/0031-9155/57/6/1459 PHMBA7 0031-9155 Google Scholar

Citation Download Citation

Calvin B. Shaw and Phaneedra K. Yalavarthy "Effective contrast recovery in rapid dynamic near-infrared diffuse optical tomography using ℓ₁-norm-based linear image reconstruction method," Journal of Biomedical Optics 17(8), 086009 (9 August 2012). https://doi.org/10.1117/1.JBO.17.8.086009

Published: 9 August 2012

Access the abstract

JOURNAL ARTICLE
11 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 22 scholarly publications.

Explore citations on Lens.org

RIGHTS & PERMISSIONS

Get copyright permission Get copyright permission on Copyright Marketplace

KEYWORDS

Diffuse optical tomography

Near infrared

Data modeling

Inverse problems

Diffuse optical imaging

Image restoration

Tissue optics

1.

Introduction

2.

Dynamic NIR Diffuse Optical Tomography: Forward Problem

Eq. (1)

Eq. (2)

3.

Dynamic NIR Diffuse Optical Tomography: Inverse Problem

3.1.

Linearization

Eq. (3)

Eq. (4)

Eq. (5)

Eq. (6)

Eq. (7)

Eq. (8)

Eq. (9)

Eq. (10)

3.2.

ℓ1-Based Linear Reconstruction Method

Eq. (11)

Eq. (12)

3.3.

Regularized Minimal Residual Method: ℓ2-Based

Eq. (13)

Algorithm 1

Eq. (14)

4.

Simulation and Experimental Evaluation

4.1.

Numerical Experimental Data

Fig. 1

Fig. 2

4.2.

Experimental Phantom Data

5.

Results

Fig. 3

Fig. 4

Fig. 5

Table 1

6.

Discussion

Fig. 6

7.

Conclusions

Acknowledgments

References

Show All Keywords

Keywords/Phrases

Search In:

Publication Years

$ℓ_{1}$ -Based Linear Reconstruction Method

Regularized Minimal Residual Method: $ℓ_{2}$ -Based