Minimal residual method provides optimal regularization parameter for diffuse optical tomography

Ravi Prasad K. Jagannath; Phaneedra K. Yalavarthy

doi:10.1117/1.JBO.17.10.106015

10 October 2012 Minimal residual method provides optimal regularization parameter for diffuse optical tomography

Ravi Prasad K. Jagannath, Phaneedra K. Yalavarthy

Author Affiliations +

Journal of Biomedical Optics, Vol. 17, Issue 10, 106015 (October 2012). https://doi.org/10.1117/1.JBO.17.10.106015

Abstract

The inverse problem in the diffuse optical tomography is known to be nonlinear, ill-posed, and sometimes under-determined, requiring regularization to obtain meaningful results, with Tikhonov-type regularization being the most popular one. The choice of this regularization parameter dictates the reconstructed optical image quality and is typically chosen empirically or based on prior experience. An automated method for optimal selection of regularization parameter that is based on regularized minimal residual method (MRM) is proposed and is compared with the traditional generalized cross-validation method. The results obtained using numerical and gelatin phantom data indicate that the MRM-based method is capable of providing the optimal regularization parameter.

1. Introduction

Diffuse optical tomography (DOT) is one of the promising imaging modalities that provides functional information of the soft biological tissues under investigation using near infrared (NIR) light (600 to 1000 nm) as the probing radiation. This NIR light illuminates the tissue through the optical fibers positioned on the surface, and the same fibers also collect the emerging light at the boundary.¹^–⁶ By using these limited boundary measurements of the exiting light, the optical properties of the tissue are estimated by making use of a model-based image-reconstruction algorithm.⁷ As NIR light is nonionizing, the prolonged monitoring of tissue physiology is feasible and is capable of providing functional parameters of the tissue when multiple wavelengths of NIR light are used.⁸

Reconstructing the internal distribution of optical properties by solving the inverse problem using the limited boundary data is a challenging task and, in general, only approximate solutions are obtained. This inverse problem is nonlinear, ill-posed, and under-determined, which is the result of the multiple scatterings of the light in a dense biological tissue and requires advanced numerical methods to solve it. One of the most powerful methods in existence is the Newton-Raphson-based technique via Tikhonov regularization.⁹^,¹⁰ This involves the inversion of the large Hessian matrix obtained by using the Jacobian (or sensitivity) matrix. This inversion typically requires regularization, and the most often used one is Tikhonov-type regularization. The Tikhonov-type regularization adds a fixed diagonal value to the Hessian, also known as a regularization parameter. The regularization parameter in this case dictates the quality of the reconstructed image. A higher value compared to the optimal one smoothes the reconstructed images resulting in loss of resolution; a lower case results in high-frequency noise in images. This makes the choice of such a regularization parameter critical in the image reconstruction procedure.⁹^–¹¹ The regularization parameter choice theoretically depends on the noise characteristics both in data and image space, and it is highly impractical to obtain the same in real-time, leading to empirical choice of the same.⁹^,¹⁰

Determination of the regularization parameter empirically leads to subjectivity and an unwarranted bias in the solution. Also, such empirical determination varies with the problem at hand and requires prior knowledge about the target images as well as noise in the data. One of the most commonly used automated technique for determining the regularization parameter in Newton-type inversion schemes is generalized cross-validation (GCV),¹²^–¹⁴ which is often used for Tikhonov-type regularization. The GCV method is applicable to linear inverse problems and is widely used in image deblurring.¹²^–¹⁵ As the inverse problem here is nonlinear, which is solved in a series of linear steps, application of GCV to find the optimal regularization parameter will be shown to be suboptimal for the problem under consideration. More importantly, it will be shown that one can obtain the optimal regularization parameter in these cases using the regularized minimal residual method.¹⁶ Such determination of the regularization parameter requires an additional optimization procedure and is not computationally expensive compared to the estimation of optical properties. The cases considered here are both numerical as well as gelatin phantom experiments to prove that the proposed minimal residual method provides an optimal regularization parameter.

As the aim of this work is to introduce a novel way of obtaining the optimal regularization parameter, the discussion is limited to a two-dimensional continuous-wave (CW) case, in which the distribution of coefficient of absorption ( $μ_{a}$ ) is estimated from the limited boundary measurements of the amplitude of exiting light. The implementation of the same is achieved through MATLAB-based open-source NIRFAST,¹⁷ and the concerned code of the proposed algorithm is also available for enthusiastic readers/users as an open-source.¹⁸

2. Methods

2.1.

Diffuse Optical Tomography: Forward Problem

The forward problem in diffuse optical imaging is defined as finding the photon fluence rate through the medium given the optical properties of the medium and source detector locations. The propagation of CW NIR light in a thick nonhomogeneous medium like biological tissues is modeled using an approximation of the radiative transfer equation, called as diffusion equation (DE),⁷^,¹⁷^,¹⁹ given by,

Eq. (1)

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

where

Φ (r)

is the photon fluence rate at a given position

r

and

S_{o} (r)

is the isotropic light source. The optical diffusion and absorption coefficients are given by

D (r)

and

μ_{a} (r)

, respectively and also note that

Eq. (2)

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

where

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

is the reduced scattering coefficient, which is defined as

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

.

μ_{s}

is the scattering coefficient and

g

is the anisotropy factor. A Type-III boundary condition is employed to the DE [Eq. (1)], which accounts for the refractive-index mismatch at the boundary.²⁰ The finite element method (FEM) is used to solve Eq. (1) to generate the modeled data

G (μ_{a})

for a given distribution of absorption coefficient

μ_{a} (r)

,⁷^,¹⁷^,¹⁹ and under the Rytov approximation, the data obtained (

y

) is the natural logarithm of the intensity (

I

) i.e.,

y = [\ln (I)]

.

2.2.

Diffuse Optical Tomography: Inverse Problem

The goal of the inverse problem is to recover the distribution of absorption coefficient using the limited number of boundary measurements of light intensity. The inverse problem is typically posed as a least-square minimization scheme,¹^,⁵^–⁷^,¹⁰^,¹⁷ where the aim is to iteratively match the modeled data obtained using the forward problem with the experimentally measured boundary data. A regularization term is typically added to stabilize the solution. The objective function ( $Ω$ ) to be minimized with respect to $μ_{a}$ is written as follows,

Eq. (3)

Ω = {‖ y - G (μ_{a}) ‖}^{2} + α {‖ μ_{a} - μ_{a 0} ‖}^{2} .

This type of scheme is known to be a Tikhonov regularization scheme, where

y

is the experimental data collected at the boundary of the object, i.e.,

y = \ln {(A)}^{measured}

and

G (μ_{a})

is obtained by solving the diffusion-based forward model, as described in the previous section, for the given distribution of the absorption coefficient (

μ_{a}

).

μ_{a 0}

is the initial guess for the absorption coefficient and

α

is known to be the regularization parameter. The

α

can be shown to be the ratio of variances in the data, and estimated properties¹⁰ and such determination requires prior information about the data noise and image noise characteristics, which are impractical to obtain in real-time.

Writing the Taylor series expansion of $G (μ_{a})$ around the initial distribution of the absorption coefficient $μ_{a 0}$ ,

Eq. (4)

G (μ_{a}) = G (μ_{a 0}) + G^{'} (μ_{a} - μ_{a 0}) + {(μ_{a} - μ_{a 0})}^{t} G^{''} (μ_{a} - μ_{a 0}) + \dots,

where

G^{'} = J = \frac{d G (μ_{a})}{d μ_{a}}

is the Jacobian and similarly

G^{''}

is the second order derivative. Ignoring the higher order terming by assuming the initial guess is very close to the solution gives rise to the linearized inverse problem²¹ [using Eq. (4) in Eq. (3)]

Eq. (5)

Ω = {‖ δ - J Δ μ_{a} ‖}^{2} + α {‖ Δ μ_{a} ‖}^{2},

where

Δ μ_{a} = (μ_{a} - μ_{a 0})

and assuming

δ = y - G (μ_{a 0})

as the data-model misfit. Note that the inverse problem in Eq. (3) is nonlinear, and the above equation is linear, making the nonlinear inverse problem to be solved in a series of linear steps. The direct solution for Eq. (5) is given by¹⁰^,¹⁶

Eq. (6)

(J^{T} J + α I) Δ μ_{a} = J^{T} δ,

which is known to be the Euler equations,¹⁶ which is known to be the direct method for minimizing the Eq. (5). The computational complexity of obtaining a direct solution using the above equation is

O (N N^{3})

with NN being the number of finite element nodes, making it computationally expensive for large problems. After each inversion, the corresponding update in

μ_{a}

is followed by recomputation of Jacobian. The new linear inverse problem is reformulated [Eq. (5)]; this iterative procedure is stopped when the

L^{2}

norm of the data-model misfit does not improve by more than 2%.¹⁰ Typically

α

is chosen empirically at every iteration and is known to influence the reconstructed image quality.

2.2.1.

Regularized minimal residual method

The minimization problem that is given in Eq. (5) could also be solved by using the minimal residual method (MRM).¹⁶^,²² The regularized minimal residual method is equivalent to the regularized steepest descent method and is an iterative technique that minimizes Eq. (5). This method could be seen as solving linear system of equations in Eq. (5) using an iterative method.¹⁶^,²² The formulation for the $n$ update in $μ_{a}$ for a given $α$ is as follows:

Eq. (7)

{(Δ μ_{a})}_{n + 1} = {(Δ μ_{a})}_{n} - k_{n}^{α} l_{n}^{α},

where

l_{n}^{α} = J^{T} r_{n} + α {(Δ μ_{a})}_{n}

with

r_{n} = J {(Δ μ_{a})}_{n} - δ

and

Eq. (8)

k_{n}^{α} = \frac{{∥ l_{n}^{α} ∥}^{2}}{{∥ J l_{n}^{α} ∥}^{2} + α {∥ l_{n}^{α} ∥}^{2}} .

The iteration is stopped when the

{‖ r_{n} ‖}^{2} - {‖ r_{n - 1} ‖}^{2}

is less than

ϵ = 10^{- 6}

(equal within single-precision limits). The computational complexity of the MRM method is

O (m * N N^{2})

, where m is the number of iterations needed to reach the stopping criterion. It will be shown later with a numerical example that minimal residual method is equivalent to solving the Euler equation (direct solution) given in Eq. (6).

2.2.2.

Estimating optimal regularization parameter using minimal residual method

The choice of proper regularization is a key to successful estimation of the optical parameters. The regularization can be obtained optimally using the regularized MRM. For every given $α$ the MRM converges to a suitable distribution of the absorption coefficient ( $μ_{a}$ ) as explained in the previous section. Now this solution of the inverse problem for a given regularization $α$ can be used to find the corresponding data-model misfit, computed using the updated $μ_{a}$ . This ensures that the data-model misfit [ ${‖ y - G (μ_{a}) ‖}^{2}$ ], depends only on the scalar quantity ( $α$ ). So the optimal $α$ is the one that gives the least data-model misfit (or matches the experimental data with the modeled data within a smallest possible neighborhood). This existence of optimal $α$ through the proposed method is discussed in the appendix and shown analytically and graphically that this procedure results in an optimal (suboptimal in noisy data cases) solution for $α$ . Implementation of finding optimal $α$ is algorithmically summarized in the following steps.

1. Compute $J$ , $μ_{a}$ , and $δ$ for the given iteration.
2. Find the $α$ for which ${‖ y - G (μ_{a} + Δ μ_{a}^{α}) ‖}^{2}$ is minimum, wherein $Δ μ_{a}^{α}$ is obtained using regularized MRM described in Sec. 2.2.1 with the inputs computed in step 1.

For finding such an $α$ , a gradient-free simplex method type algorithm²³ is used due to its computational efficiency. Note that the objective function in here (Step 2) requires computing the modeled data (or equivalently solving the forward model) at every iteration of finding an optimal $α$ .

2.2.3.

Estimation optimal regularization using generalized cross-validation

The GCV ¹²^–¹⁴ method is the most popular method for estimating the regularization parameter ( $α$ ). This is based on the principle that an omitted data point could be easily estimated using the regularized solution that is obtained using the reduced dataset. Obtaining an estimate of the regularization is achieved by minimizing a function $G (α)$ defined to be:

Eq. (9)

G (α) = \frac{{‖ {(J J^{T} + α I)}^{- 1} δ ‖}^{2}}{{(trace {(J J^{T} + α I)}^{- 1})}^{2}} .

Efficient evaluation of the

G (α)

is possible by applying the singular value decomposition (SVD) of the Jacobian and doing some algebraic simplification. If the SVD of the Jacobian matrix(

J

) is given to be

J = U Σ V^{t}

, then the above functional is simplified to the following form:

Eq. (10)

G (α) = \frac{\sum_{i = 1}^{rank (J)} {(\frac{u_{i}^{T} δ}{σ_{i}^{2} + α^{2}})}^{2}}{{(\sum_{i = 1}^{rank (J)} \frac{1}{σ_{i}^{2} + α^{2}})}^{2}},

where the

i

th column of the matrix

U

is represented as

u_{i}

, and

σ_{i}

are the singular values of the Jacobian. The function [

G (α)

] thus defined is continuous for every value of

α

, and this is minimized with respect to

α

using the same methods as in the previous section. This results in the required regularization parameter (

α_{GCV}

) for the inverse problem. This is achieved through the usage of MATLAB-based open-source regularization toolbox.¹⁵

3. Experiments

3.1.

Simulation Studies

The effectiveness of the proposed method is assessed on circular meshes with the background optical properties being $μ_{a} = 0.01 {mm}^{- 1}$ , $μ_{s}^{'} = 1 {mm}^{- 1}$ and a uniform refractive index of 1.33. The diameter of all circular meshes used in this work is 86 mm, and the experimental data are generated on a fine FEM mesh with 10,249 nodes was (corresponding to 20,160 linear triangular elements) and the reconstruction was performed on a coarser mesh with 1,785 nodes (corresponding to 3418 linear triangular elements). The source-detector geometry had 16 fibers that were arranged in circular, equi-spaced fashion and each fiber has a dual function of both acting as source and detector. When one fiber acts as a source, the remaining 15 fibers act as detectors and hence we obtain 240 ( $15 \times 16$ ) light amplitude measurements ( $N M$ ). To simulate experimental conditions, the source was modeled as the Gaussian source with a full width at half maximum of 3 mm²⁴ and was placed at the one mean transport length inside the boundary. All meshes were centered around origin.

Two cases of targets were considered here, with the initial case of showing that the solution obtained using MRM is equivalent to standard reconstruction using direct method [Euler equation; Eq. (6)]. For this, a circular target of radius 7.5 mm was placed at (15,0) having $μ_{a} = 0.02 {mm}^{- 1}$ , $μ_{s}^{'} = 1 {mm}^{- 1}$ and the distribution of $μ_{a}$ is shown in first column of Fig. 1. A regularization parameter of 0.01 has been used in this case. A normally distributed Gaussian noise of 1% has been added to the numerically generated data to mimic the experimental case.²⁴ In the second case, two circular targets that were separated by 5 mm centered around origin with a radius of 7.5 mm were considered. The optical properties were similar to the first case, and the target distribution of $μ_{a}$ is given in the first column of Fig. 2. At every iteration, the regularization parameter ( $α$ ) has been obtained using both MRM and GCV methods as described in Secs. 2.2.2 and 2.2.3, respectively. In this case, both 1% and 3% normally distributed Gaussian noise added data was considered to show the effectiveness of the proposed method.

Fig. 1

Reconstructed $μ_{a}$ distributions using standard (direct inversion) and minimal residual methods (iterative inversion) with numerically generated 1% noisy data and $α$ being 0.01. One-dimensional cross-sectional plot of $μ_{a}$ along the line (shown on the target) for the target and reconstructed results presented here are given in the second row.

Fig. 2

Comparison of reconstruction performance using the two methods discussed in this work for the two circular targets of dimension 7.5 mm radius separated by a distance of 5 mm with numerically generated 1% and 3% noisy data. The reconstructed images obtained along with the corresponding method with the percentage of the noise added (given on top of image) are given along with the target image. One-dimensional cross-sectional plot of $μ_{a}$ along the line (shown on the target) for the target and reconstructed results presented here are given in the second row.

3.2.

Gelatin Phantom Studies

For objective verification of the method, the experimental data collected using a gelatin phantom was used. A multilayer cylindrical gelatin phantom of diameter 86 mm, height 25 mm was fabricated following the procedure as explained in Ref. 25. Mixtures of 80% of deionised water and 20% of gelatin (G2625, Sigma Inc.) along with India ink for absorption and Titanium oxide ( ${TiO}_{2}$ ) for scattering was used. These layers of gelatin were successively hardened to form a three-layer phantom to mimic typical breast tissue. The outer layer had optical properties of $μ_{a} = 0.0065 {mm}^{- 1}$ , $μ_{s}^{'} = 0.65 {mm}^{- 1}$ with a thickness of 10 mm, mimicking the fatty layer of the breast. The inner layer, mimicking the fibro-glandular layer, had properties of $μ_{a} = 0.01 {mm}^{- 1}$ , $μ_{s}^{'} = 1.0 {mm}^{- 1}$ . To mimic the tumor characteristics, along the Z direction a cylindrical hole having a radius of 8 mm and height of 24 mm was made and filled with intra-lipid mixed with India ink to act as an absorptive target having the optical properties $μ_{a} = 0.02 {mm}^{- 1}$ , $μ_{s}^{'} = 1.2 {mm}^{- 1}$ . The optical properties were validated using 785 nm wavelength laser diode as the source. The log-amplitude data collection of the exiting light were collected using only one layer of fibers placed at the middle of the cylinder (at $z = 0 mm$ ). The collected data was calibrated using the coarser mesh that was used in the simulation studies using the procedure explained in Ref. 26. The adipose layer optical properties were used as an initial guess for reconstruction procedures.

All computations were carried out on a Linux workstation with dual quad-core Intel Xeon processor of 2.33 GHz with 64 GB RAM. The algorithms were implemented in MATLAB.

4. Results

Initial results to show that regularized MRM is equivalent to standard (direct) reconstruction [Eq. (6)] are given in Fig. 1. The reconstruction methods that were used are given on top of each $μ_{a}$ distribution. A one-dimensional profile plot across the dotted line as shown in Target $μ_{a}$ is given in the second row of Fig. 1. Note that the number of iterations that were needed to converge in both cases were 13 and in case of MRM at each iteration included about 400 inner iterations. From results shown in Fig. 1, it can be clearly seen that the reconstructed $μ_{a}$ distribution using MRM matches with the one obtained using the standard method within 2% error, assuring that the solution obtained using MRM is numerically equivalent to the solution obtained using direct inversion [Euler equation, Eq. (6)]. Using the standard reconstruction method, the computation time needed per iteration was 0.829 s and for MRM, it was 0.688 s.

In the second case of numerical experiments using 1% and 3% noisy data, the regularization parameter was found using both MRM and GCV methods as outlined in Secs. 2.2.2 and 2.2.3. The obtained reconstruction results are given in Fig. 2. The one-dimensional cross-sectional plot for the obtained results similar to Fig. 1 is given in the second row of Fig. 2. The MRM-based estimation of regularization for the 1% noisy data case has resulted in convergence in four iterations with $α$ being 1.3e-2, 1.1e-3, 2.7e-4, and 6.0e-6. For the GCV the convergence was achieved in three iterations with $α$ being 1.4e-2, 1.2e-2, and 1000. Note that for the 3% noisy data case, the obtained $α$ values were in the same numerical range. The results also indicate that the resolution characteristics were better using the MRM-based regularization parameter compared to its counterpart. The number of iterations that were required to find the $α$ in either case were highly dependent on the initial guess.

The results that were obtained using estimated $α$ with MRM and GCV methods for the case of gelatin phantom data are given in Fig. 3. In here as well, the MRM-based regularization parameter selection resulted in convergence in four iterations with $α$ being 1.1e-2, 2.6e-3, 8.1e-4, and 6.0e-6. Correspondingly, the convergence for the GCV-based method was in three iterations, with $α$ being 8.0e-3, 7.8e-3, and 1.69. The one-dimensional cross-sectional plot along the dashed line in the target is given in the second row of Fig. 3. The results indicate that the performance of the MRM to find an optimal $α$ is better compared to the GCV-based method as the reconstructed target $μ_{a}$ value is closer to the expected one. The total computational time for the three cases considered in this work is reported in Table 1 and as mentioned earlier, the MRM-method computational complexity is much higher compared to the GCV-method.

Fig. 3

Similar effort as Fig. 2 for the case of experimental gelatin phantom data.

Table 1

Comparison of total computational time (in seconds) including the overhead, for three reconstruction results (first column) presented in this work. The corresponding method is given in the first row of the table. The total number of iterations taken to converge are given in the parentheses.

Result	MRM method	GCV method
Fig. 2 (1% Noise-case)	55.42 (3)	12.20 (3)
Fig. 2 (3% Noise-case)	104.88 (4)	12.15 (3)
Fig. 3 (Gelatin phantom)	167.75 (4)	12.36 (3)

5. Discussion

The diffuse optical tomographic image reconstruction problem is ill-posed in nature due to multiple scatterings of near infrared light. Many reconstruction methods were proposed for improving the reconstruction accuracy, making it more quantitative in nature, which uses structural and/or functional information.¹⁰^,²⁷^–³⁰ The quantitative accuracy of the results are dependent on the choice of regularization, making the problem of finding optimal $α$ in an automated way highly relevant to the ongoing efforts. Moreover, such an automated method will remove the unwarranted bias in the reconstruction results due to heuristic choice of regularization parameter. Until now, the most successful methods in the literature include the GCV and L-curve methods.³¹^–³³ In this work, the proposed method is compared with the GCV-based approach and is found to be sub-optimal. The L-curve method is not considered in this work as results in overly smooth solutions in DOT³¹ as it can not exhibit a clear corner to find an optimal regularization parameter.³²^,³³

Finding the regularization parameter using the MRM method could also be achieved using Euler equations [Eq. (6)], but this increases the computational complexity associated with it by an order [ $O (N N^{2})$ versus $O (N N^{3})$ ]. In the case of GCV, the regularization parameter increased abruptly before the last global iteration of the reconstruction procedure, resulting in sub-optimal estimates of $μ_{a}$ compared to the MRM-based method. This trend of GCV inability to find an optimal $α$ is also discussed in Ref. 32, where this is caused by either a function being flat or displaying multiple minima, leading to a sub-optimal choice of $α$ , causing over smoothing of the reconstructed images.

To show the inner iterative behavior of the proposed method, the semi-log plot of data-model misfit and the corresponding $α$ was plotted in Fig. 4 for 60 iterations for the set of results presented in Fig. 2 (first global iteration). These plots indicate that there is a clear convergence in finding $α$ and correspondingly the data-model misfit is also a minimum for the optimal $α$ .

Fig. 4

Plot of the estimated regularization parameter ( $α$ ) using MRM and the corresponding data-model misfit $[Ω = {‖ y - G (μ_{a}) ‖}^{2}]$ as a function of the inner iteration number for the first global iteration of the results presented in Fig. 2.

For the case of numerically generated 3% noisy data, the reconstruction results were poor in terms of resolution for the case GCV (Fig. 2), making the targets indistinguishable. For the MRM case, it is clearly evident that the reconstructed image quality is better in comparison to the result obtained using the GCV method (Fig. 2). The results also indicate that the spatial resolution in the reconstructed $μ_{a}$ image is higher by the usage of the MRM method compared to the deployment of the GCV method. The boundary artifacts for the case of the MRM method with 3% noisy data is pronounced compared to the result obtained using the GCV method (Fig. 2) as the last iteration of GCV always leads to higher value of regularization parameter making the reconstructed image appear more smooth. Note that for the MRM-based method, the regularization parameter monotonically decreases with iterations.

Finding the optimal $α$ requires repeated solving of the forward problem, which might be computationally expensive in the case of large problems, especially involving a large number of measurements. The additional step of finding the optimal $α$ through the MRM method for the cases discussed in here, requires about 9.6 s, and for the GCV method it is 2.7 s. This additional computational complexity (Table 1) added by the optimization procedure of finding optimal $α$ is justifiable as it removes the unwarranted bias in the reconstruction results due to either poor or desirable choice of $α$ heuristically and does not require any additional information in finding optimal $α$ . This computation time could be further reduced by the usage of GPU-based computing environments,³⁴ which can speed up the computations in providing the optimal $α$ .

6. Conclusions

The regularization parameter that is typically used in the diffuse optical tomographic image reconstruction procedure provides a fidelity in improving the reconstructed image quality, with a caveat that such fidelity can also bias the results. An automated estimation of regularization parameter that is based on regularized minimal residual method is presented in this work and is compared with the traditional GCV-based method. The reconstruction results using numerical and gelatin phantom data indicate that the proposed MRM-based method can provide an optimal regularization parameter, overcoming the pitfalls of the GCV-based method. Such an automated method requires repeated computations of forward solutions and is not computationally complex compared to the direct image reconstruction step. The computer programs for the developed method that is used in this work are provided as open source.¹⁸

Acknowledgments

This work is supported by the Department of Atomic Energy Young Scientist Research Award (No. 2010/20/34/6/BRNS) by the Government of India. The authors are also thankful to NIR imaging group of Dartmouth for providing necessary experimental data and finite element meshes that were needed to carry out this work. Ravi Prasad K. Jagannath acknowledges the university grants commission (UGC) senior research fellowship (SRF). The authors also acknowledge the IISc Mathematics Initiative (IMI) program.

Appendices

Appendix:

Existence of Optimal Regularization Parameter

In every iterative step of the diffuse optical tomographic image reconstruction, the modeled data is matched with the experimental data in the least-square sense as explained in Sec. 2.2. So the objective function to be minimized is posed as follows,

Eq. (11)

Ω = {‖ y - G (μ_{a}) ‖}^{2} .

Since the iterative procedure of minimizing this starts with a guess of the initial distribution of the absorption coefficient (

μ_{a 0}

), the above equation can be written as

Eq. (12)

Ω = {‖ y - G (μ_{a 0} + Δ μ_{a}) ‖}^{2} .

Applying Taylor’s series expansion, as explained earlier in Sec. 2.2, modifies the above objective function to

Eq. (13)

Ω = {‖ (δ - J Δ μ_{a}) ‖}^{2} = {(δ - J Δ μ_{a})}^{T} (δ - J Δ μ_{a}) .

The update (

Δ μ_{a}

) in the above equation can be obtained using the regularized minimal residual method (MRM), which requires a scalar (

α

) as a regularization parameter (refer to Sec. 2.2.1). This implies that for every given positive real number (

α

) there is a unique update (

Δ μ_{a}^{α}

), in turn making the objective function

Ω

as a continuous function of

α

.

Finding the minima of the objective function, $Ω$ , with respect to $α$ requires the calculation of a first-order derivative, which can be written as

Eq. (14)

\frac{\partial Ω}{\partial α} = 2 {[J^{T} (J Δ μ_{a}^{α} - δ)]}^{T} (\frac{\partial Δ μ_{a}^{α}}{\partial α}) .

Similarly, the second-order derivative is given by

Eq. (15)

\frac{\partial^{2} Ω}{\partial α^{2}} = 2 {[J^{T} (J Δ μ_{a}^{α} - δ)]}^{T} (\frac{\partial^{2} Δ μ_{a}^{α}}{\partial α^{2}}) + {(J^{T} J \frac{\partial Δ μ_{a}^{α}}{\partial α})}^{T} (\frac{\partial Δ μ_{a}^{α}}{\partial α}) .

The optimal regularization parameter (

α^{*}

) is the one that makes the first-order derivative [Eq. (14)] go to zero, simultaneously making the second-order derivative [Eq. (15)] have a positive value, ensuring that the objective function attains the minimum value.

The optimal solution ( $Δ μ_{a}$ ) for the linearized problem makes the objective function, given by Eq. (13), achieve a minimal value. Assuming perfect system characteristics (including noiseless data) ensures that the residue $|| r || = || (J Δ μ_{a}^{α^{*}} - δ) ||$ is zero, where the $α^{*}$ represents the optimal regularization parameter. So in the case of $α^{*}$ , the right-hand side of Eq. (14) becomes zero, as $J Δ μ_{a}^{α^{*}} - δ = 0$ (with the partial derivative being non-zero). Using this in the right-hand side of Eq. (15), results in

Eq. (16)

\frac{\partial^{2} Ω}{\partial α^{2}} {|_{α^{*}} = [{(\frac{\partial Δ μ_{a}^{α}}{\partial α})}^{T} J^{T} J (\frac{\partial Δ μ_{a}^{α}}{\partial α})] |}_{α^{*}},

as

J^{T} J

is symmetric and positive semi-definite matrix,¹⁰ and the partial derivative term is non-zero ensures that Eq. (16) results in a positive quantity. This assures that

α^{*}

exists for the minimization problem and is achievable when the objective function for the image reconstruction problem has a minima.

This argument could also be validated through graphical representation by plotting the objective function ( $Ω$ ) as a function of the regularization parameter $α$ for an example problem. Figure 5 shows such an effort for the results presented in Fig. 2, where the plot shows $α^{*}$ for which the function has a minima, but not zero due to noise in the data and numerical model errors. This makes the choice of $α^{*}$ in the real-time (for a limited number of iterations) being only sub-optimal (optimal in the asymptotic range).

Fig. 5

The semi-logy plot of the objective function $[Ω = {‖ y - G (μ_{a}^{α}) ‖}^{2}]$ versus regularization parameter $α$ for the results presented in Fig. 2 (first iteration of the 1% noisy data case). The optimal $α$ , indicated by $α^{*}$ in the plot, is where $Ω$ takes the minimal value.

References

1.

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

2.

S. L. JacquesB. W. Pogue, “Tutorial on diffuse light transport,” J. Biomed. Opt., 13 (4), 041302 (2008). http://dx.doi.org/10.1117/1.2967535 JBOPFO 1083-3668 Google Scholar

3.

A. GibsonJ. C. HebdenS. R. Arridge, “Recent advances in diffuse optical tomography,” Phys. Med. Biol., 50 (4), R1 –R43 (2005). http://dx.doi.org/10.1088/0031-9155/50/4/R01 PHMBA7 0031-9155 Google Scholar

4.

A. GibsonH. Dehghani, “Diffuse Optical Imaging,” Phil. Trans. R. Soc. A, 367 (1900), 3055 –3072 (2009). http://dx.doi.org/10.1098/rsta.2009.0080 PTRMAD 1364-503X Google Scholar

5.

S. R. Arridge, “Optical tomography in medical imaging,” Inv. Problems, 15 (2), R41 –R93 (1999). http://dx.doi.org/10.1088/0266-5611/15/2/022 INPEEY 0266-5611 Google Scholar

6.

S. R. ArridgeJ. C. Hebden, “Optical imaging in medicine: II. modelling and reconstruction,” Phys. Med. Biol., 42 (5), 841 –853 (1997). http://dx.doi.org/10.1088/0031-9155/42/5/008 PHMBA7 0031-9155 Google Scholar

7.

H. Dehghaniet al., “Numerical modelling and image reconstruction in diffuse optical tomography,” Phil. Trans. R. Soc. A, 367 (1900), 3073 –3093 (2009). http://dx.doi.org/10.1098/rsta.2009.0090 PTRMAD 1364-503X Google Scholar

8.

X. Inteset al., “Diffuse optical tomography with physiological and spatial a priori constraints,” Phys. Med. Biol., 49 (12), N155 –N163 (2004). http://dx.doi.org/10.1088/0031-9155/49/12/N01 PHMBA7 0031-9155 Google Scholar

9.

B. W. Pogueet al., “Spatially variant regularization improves diffuse optical tomography,” Appl. Opt., 38 (13), 2950 –2961 (1999). http://dx.doi.org/10.1364/AO.38.002950 APOPAI 0003-6935 Google Scholar

10.

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

11.

B. W. Pogueet al., “Statistical analysis of nonlinearly reconstructed near-infrared tomographic images: Part I theory and simulations,” IEEE Trans. Med. Imag., 21 (7), 755 –763 (2002). http://dx.doi.org/10.1109/TMI.2002.801155 ITMID4 0278-0062 Google Scholar

12.

G. GolubU. von Matt, “A generalized cross-validation for large-scale problems,” J. Comput. Graph. Statist., 6 (1), 1 –34 (1997). http://dx.doi.org/10.1080/10618600.1997.10474725 1061-8600 Google Scholar

13.

N. NguyenP. MilanfarG. Golub, “A computationally efficient superresolution image reconstruction algorithm,” IEEE Trans. Image Process., 10 (4), 573 –583 (2001). http://dx.doi.org/10.1109/83.913592 IIPRE4 1057-7149 Google Scholar

14.

P. C. HansenJ. G. NagyD. P. O. Leary, Deblurring Images: Matrices, Spectra, and Filtering, SIAM, Philadelphia (2006). Google Scholar

15.

P. C. Hansen, “Regularization Tools Version 4.0 for Matlab 7.3,” Numer. Algorithms, 46 (2), 189 –194 (2007). http://dx.doi.org/10.1007/s11075-007-9136-9 NUALEG 1017-1398 Google Scholar

16.

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

17.

H. Dehghaniet al., “Near infrared optical tomography using NIRFAST: algorithms for numerical model and image reconstruction algorithms,” Commun. Numer. Meth. Eng., 25 (6), 711 –732 (2009). http://dx.doi.org/10.1002/cnm.v25:6 CANMER 0748-8025 Google Scholar

18.

( (2012) https://sites.google.com/site/sercmig/home/regdot September ). 2012). Google Scholar

19.

S. R. ArridgeM. Schweiger, “Photon-measurement density functions. Part 2: finite-element-method calculations,” Appl. Opt., 34 (34), 8026 –8037 (1995). http://dx.doi.org/10.1364/AO.34.008026 APOPAI 0003-6935 Google Scholar

20.

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

21.

M. SchweigerS. R. ArridgeI. Nissila, “Gauss-Newton method for image reconstruction in diffuse optical tomography,” Phys. Med. Biol., 50 (10), 2365 –2386 (2005). http://dx.doi.org/10.1088/0031-9155/50/10/013 PHMBA7 0031-9155 Google Scholar

22.

C. B. ShawP. K. Yalavarthy, “Effective contrast recovery in rapid dynamic near-infrared diffuse optical tomography using

ℓ_{1}

-norm-based linear image reconstruction method,” J. Biomed. Opt., 17 (8), 086009 (2012). http://dx.doi.org/10.1117/1.JBO.17.8.086009 JBOPFO 1083-3668 Google Scholar

23.

J. C. Lagariaset al., “Convergence properties of the nelder-mead simplex method in low dimensions,” SIAM J. Optimization, 9 (1), 112 –147 (1998). http://dx.doi.org/10.1137/S1052623496303470 SJOPE8 1095-7189 Google Scholar

24.

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

25.

B. W. PogueM. Patterson, “Review of tissue simulating phantoms for optical spectroscopy, imaging and dosimetry,” J. Biomed. Opt., 11 (4), 041102 (2006). http://dx.doi.org/10.1117/1.2335429 JBOPFO 1083-3668 Google Scholar

26.

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), 182 –193 (2000). http://dx.doi.org/10.1117/1.429985 JBOPFO 1083-3668 Google Scholar

27.

P. K. Yalavarthyet al., “Structural information within regularization matrices improves near infrared diffuse optical tomography,” Opt. Express, 15 (3), 8043 –8058 (2007). http://dx.doi.org/10.1364/OE.15.008043 OPEXFF 1094-4087 Google Scholar

28.

S. H. KatamreddyP. K. Yalavarthy, “Model-resolution based regularization improves near infrared diffuse optical tomography,” J. Opt. Soc. Am. A, 29 (5), 649 –656 (2012). http://dx.doi.org/10.1364/JOSAA.29.000649 JOAOD6 0740-3232 Google Scholar

29.

M. E. EamesH. Dehghani, “Wavelength dependence of sensitivity in spectral diffuse optical imaging: effect of normalization on image reconstruction,” Opt. Express, 16 (22), 17780 (2008). http://dx.doi.org/10.1364/OE.16.017780 OPEXFF 1094-4087 Google Scholar

30.

F. LarussonS. FantiniE. L. Miller, “Hyperspectral image reconstruction for diffuse optical tomography,” Biomed. Opt. Express, 2 (34), 946 –965 (2011). http://dx.doi.org/10.1364/BOE.2.000946 BOEICL 2156-7085 Google Scholar

31.

J. P. Culveret al., “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 (2), 235 –247 (2003). http://dx.doi.org/10.1118/1.1534109 MPHYA6 0094-2405 Google Scholar

32.

T. Correiaet al., “Selection of regularization parameter for optical topography,” J. Biomed. Opt., 14 (3), 034044 (2009). http://dx.doi.org/10.1117/1.3156839 JBOPFO 1083-3668 Google Scholar

33.

J. Chamorro-Serventet al., “Feasibility of U-curve method to select the regularization parameter for fluorescence diffuse optical tomography in phantom and small animal studies,” Opt. Express, 19 (12), 11490 –11506 (2011). http://dx.doi.org/10.1364/OE.19.011490 OPEXFF 1094-4087 Google Scholar

34.

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

Citation Download Citation

Ravi Prasad K. Jagannath and Phaneedra K. Yalavarthy "Minimal residual method provides optimal regularization parameter for diffuse optical tomography," Journal of Biomedical Optics 17(10), 106015 (10 October 2012). https://doi.org/10.1117/1.JBO.17.10.106015

Published: 10 October 2012

Access the abstract

JOURNAL ARTICLE
7 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 25 scholarly publications.

Explore citations on Lens.org

RIGHTS & PERMISSIONS

Get copyright permission Get copyright permission on Copyright Marketplace

KEYWORDS

Data modeling

Inverse problems

Diffuse optical tomography

Absorption

Optical properties

Near infrared

Tissues

1.

Introduction

2.

Methods

2.1.

Diffuse Optical Tomography: Forward Problem

Eq. (1)

Eq. (2)

2.2.

Diffuse Optical Tomography: Inverse Problem

Eq. (3)

Eq. (4)

Eq. (5)

Eq. (6)

2.2.1.

Regularized minimal residual method

Eq. (7)

Eq. (8)

2.2.2.

Estimating optimal regularization parameter using minimal residual method

2.2.3.

Estimation optimal regularization using generalized cross-validation

Eq. (9)

Eq. (10)

3.

Experiments

3.1.

Simulation Studies

Fig. 1

Fig. 2

3.2.

Gelatin Phantom Studies

4.

Results

Fig. 3

Table 1

5.

Discussion

Fig. 4

6.

Conclusions

Acknowledgments

Appendices

Appendix:

Existence of Optimal Regularization Parameter

Eq. (11)

Eq. (12)

Eq. (13)

Eq. (14)

Eq. (15)

Eq. (16)

Fig. 5

References

Show All Keywords

Keywords/Phrases

Search In:

Publication Years