20 March 2013 Efficient gradient-free simplex method for estimation of optical properties in image-guided diffuse optical tomography
Author Affiliations +
Abstract
Typical image-guided diffuse optical tomographic image reconstruction procedures involve reduction of the number of optical parameters to be reconstructed equal to the number of distinct regions identified in the structural information provided by the traditional imaging modality. This makes the image reconstruction problem less ill-posed compared to traditional underdetermined cases. Still, the methods that are deployed in this case are same as those used for traditional diffuse optical image reconstruction, which involves a regularization term as well as computation of the Jacobian. A gradient-free Nelder–Mead simplex method is proposed here to perform the image reconstruction procedure and is shown to provide solutions that closely match ones obtained using established methods, even in highly noisy data. The proposed method also has the distinct advantage of being more efficient owing to being regularization free, involving only repeated forward calculations.
Jagannath and Yalavarthy: Efficient gradient-free simplex method for estimation of optical properties in image-guided diffuse optical tomography

Diffuse optical tomography (DOT) is an emerging medical imaging modality that uses nonionizing near-infrared (NIR) light in the range of 600 to 1000 nm as the probing medium with a capability of estimating the pathophysiological state of the tissue under investigation.1 Multiple scattering of the NIR light in the turbid medium inherently limits the resolution of the reconstructed diffuse optical image to around 0.5 cm,1 which is at least five times lower than in traditional medical imaging modalities, such as computed tomography (CT) and magnetic resonance imaging (MRI). Recent advances in developing hybrid optical imaging modalities, which combine the traditional imaging modalities (MRI/CT) with diffuse optical imaging, have proven that the diagnostic/therapeutic information provided by hybrid optical imaging is of greater value than that of the stand-alone individual imaging modalities.1

One of the main aspects of hybrid optical imaging is utilization of structural prior information provided by the traditional imaging modality to either guide or improve the solution. Among methods that guide the solution using the structural information, soft priors and hard priors are prominent.23.4.5 Other approaches, such as use of region boundary information,5 makes the estimation problem into a two-stage approach, first performing the hard priors, and second incorporating the obtained solution as a initial guess into the traditional reconstruction scheme (no-priors). The soft-priors approach uses the structural information in the regularization scheme to effectively guide the solution.23.4 In this, the number of optical parameters to be estimated remain the same as in a traditional approach, except that the regularization matrix consists of the structural information. In the hard-priors approach, the number of optical parameters that needs to be reconstructed is reduced to the number of distinct regions identifiable with traditional imaging, leading to substantial reduction in the problem size. For the case of breast imaging, the typical regions that are identifiable are adipose (fatty), fibro-glandular, and tumor,4,6 resulting in the number of optical parameters to be reconstructed equal to three.

The hard-priors approach makes the diffuse optical image reconstruction problem less ill-posed compared to the traditional approach (including soft priors). Specifically, this converts the problem from underdetermined in nature to overdetermined. Even then, the reconstruction procedures that are adapted in the hard-priors case are the same as those used in the traditional approach, the most popular being Levenberg–Marquardt (LM) minimization scheme.23.4.5 LM scheme requires a regularization parameter, which not only controls the convergence (number of iterations required),7 but also sometimes leads to biased solutions. Moreover, the LM scheme that is typically adapted to diffuse optical imaging requires calculation of the Jacobian (first-order derivative of the model).8

This work aims to show that a gradient-free simplex method,9,10 which does not require regularization or computation of gradient (or its variant Jacobian), is highly efficient compared to the existing methods for solving the image reconstruction problem with the hard-priors approach. The reconstruction results using the proposed method are compared with traditional methods using numerical and gelatin-phantom experimental cases.

Continuous-wave (CW) light propagation through a highly scattering medium such as biological tissues is modeled using an approximation of radiative transport equation called diffusion equation (DE),11 given by

(1)

·[D(r)Φ(r)]+μa(r)Φ(r)=Qo(r),
where Φ(r) represents the photon density (real values) at position r. The isotropic light source is given by Qo(r) at r. The diffusion coefficient D(r) is defined as D(r)=1/{3[μa(r)+μs(r)]}, with μs(r) and μa(r) representing the reduced scattering coefficient (assumed to be known for CW case) and absorption coefficient, respectively. For the above DE, the type III (Robin-type) boundary condition is used, as it models the refractive index mismatch at the boundary of the tissue. This partial differential equation, along with the boundary condition, is solved using finite element method (FEM),12 and Φ(r) is obtained. The modeled data [G(μa)] for the CW case becomes the natural logarithm of the amplitude data.12

The estimation of the optical absorption coefficients (here, equal to the number of regions) involves matching the experimentally measured boundary data (y) iteratively with G(μa) in a least-square sense by varying μa in different regions.2 Explicitly, the objective function that needs to be minimized with respect to μa here is given by

(2)

Ω(μa)=yG(μa)2,
where the dimension of y and μa is given by M (number of measurements) and R (number of regions), respectively. Even though RM, the minimization problem is still ill-posed (also nonlinear), requiring regularization. The popular LM minimization scheme results in an updated equation:2

(3)

Δμa=[JTJ+λI]1JT[yG(μa)],
where JT represents the transposed Jacobian J(=G(μa)/μa, dimension: M×R), I is the identity matrix, Δμa is the update, and the regularization parameter is represented by λ. The λ is systematically decreased with each iteration, here by a factor 100.25. As J is known to be semipositive definite in the case of hard priors,2 the λ affects the convergence properties of the solution.7 After every iteration, the μa is updated by adding Δμa followed by recomputing of J and G(μa) before proceeding to Eq. (3). This procedure is repeated until the change in Ω [Eq. (2)] becomes <2% between successive iterations, ensuring the convergence of the solution.

As the convergence is dependent on λ for the LM minimization, a direct search method that can minimize a scalar-valued nonlinear function, such as Ω [Eq. (2)], is the Nelder and Mead simplex algorithm.9,10 This algorithm is extremely popular for the unconstrained minimization problem, which does not require explicit or implicit calculation of the derivatives (or its variant gradients),9 making it one of the most efficient techniques. The only requirement for this algorithm is that it involves repeated calculation of G(μa). Because the objective function (Ω) has only R number of parameters, the direct search method forms the initial simplex with R+1 points, with each point having a dimension of R×1 using the initial guess of μa. These form the vertices of the simplex. The simplex method of minimization is achieved through a series of steps deployed repeatedly. Let the maximum value of the objective function (Ω) be at the point μaR, which is replaced by a new point μar via reflection, i.e.,

(4)

μar=(1+α)μ¯aαμaR,
where μ¯a represents the centroid of the the simplex and α is the reflection coefficient (positive value here it is 1). If the Ω(μar) lies between Ω(μa1) (lowest value) and Ω(μaR), then μar becomes the μa and we terminate the iteration. Otherwise, if Ω(μar)<Ω(μa1), then we expand using

(5)

μae=(1+χ)μ¯aαμaR,
where χ is the expansion coefficient taking a value of 2. If Ω(μae)<Ω(μar), then we make μa=μae if Ω(μae)Ω(μar), then μa=μar and we terminate the iteration. In case the expand condition is not satisfied, then the contract step is taken, similar to earlier Ω, and using the contraction, μa is compared and updated.10 This follows a shrinkage step.10 The coefficients that are used here are 0.5 (contract) and 0.5 (shrinkage). Note that these coefficients are universally accepted10 as the standard and fixed for the method proposed here. More details of the algorithm, along with the flow chart and tie-breaking rules, are given in Ref. 10. The procedure is repeated until the change in μa between the successive iterations becomes <1012 or the objective function reaches a value <0.1% of the initial objective function value [Ω(μa0)]. The simplex algorithm requires repeated computation of Ω (Eq. 2), which in turn evaluates G(μa) at every step, for converging to solution. This algorithm is known to provide an optimal solution when the parameter space dimensionality is within the range of data space.10 So the hard-priors method becomes a good case for deployment of this method.

To effectively assess the reconstruction performance of the proposed method, a numerical experiment involving a typical MRI-NIR human breast case is considered. The target distribution of the same is shown in the left top corner of Fig. 1. The optical properties of the tissue types are μa=0.01mm1 for the fatty region; μa=0.015mm1 for fibro-glandular; and tumor with μa=0.02mm1. The μs is assumed to be known and kept constant at 1mm1 throughout the imaging domain. The boundary measurements are taken using the 16 light source detector fiber-optic bundles arranged in an equispaced manner along the edge when one fiber delivers the light, then the rest act as detectors, which provides M=240 measurements. The mesh used for the data collection is a fine mesh with 4876 nodes (9567 triangular elements), and 1%, 5%, and 10% Gaussian noise was added to mimic the experimental case. For the reconstruction, a coarse mesh with 1969 nodes (3753 triangular elements) was used. The reconstruction results using fatty region as initial guess is given in Fig. 1, except on the top row where the guess of 0.001mm1, data noise level of 1% used which is far away from the actual solution. The details of the computational cost for the case with 1% noise are noted as follows. The traditional method converges to solution in three iterations (1.9 s) for the initial guess, being close to the actual solution (second row) with λ=0.01. If λ=100 is used in this case, the convergence is achieved in 11 iterations, taking 15.5 s. For the proposed simplex method, total computation time is 2.7 s (47 function evaluations). For the top-row results, the computation time for the traditional LM method is 3.5 s with λ=0.01 and 28.6 s with λ=100. The computational time is 12.8 s for the proposed simplex method. For the cases of 5% and 10% noise, the trend in the observed computational time is the same as for the 1% case. Figure 1 results were obtained using λ=100. The results obtained using the proposed simplex method match within 2% of results that were obtained using the traditional LM method.

Fig. 1

Comparison of reconstruction performance for the case of numerically generated data. The reconstruction techniques are given on top of each figure, with the traditional being LM method and proposed being Nelder–Mead simplex method; target represents actual μa distribution. The noise level in the data is given in the parentheses (middle two rows). The first-row results were obtained using an initial guess of 0.001mm1 (for the rest it was 0.01mm1), which is indicated in the parentheses, and the data noise level was 1%. The one-dimensional cross-section plot along the dotted line on the target image is given in the bottom row.

JBO_18_3_030503_f001.png

Next, experimental data obtained from a cylindrical-shaped three-layered gelatin phantom that mimics the breast structure, of height 25 mm and diameter 86 mm, was used to test the effectiveness of the proposed method. The details of the gelatin phantom are given in Ref. 2 (Sec. 3.3). The reconstructed distributions of μa using both traditional LM and the proposed simplex methods are given in Fig. 2 along with the one-dimensional cross-sectional profiles. The number of iterations required for the traditional method was 6 for λ=1, requiring 3.1 s. For λ=100, the convergence was achieved in 11 iterations, requiring 14.9 s. The proposed method computational time is 9.5 s (180 function evaluations). This case also shows that the proposed method results are in close agreement with the ones obtained using the traditional LM method.

Fig. 2

As in Fig. 1, for the case of experimental gelatin-phantom data. The initial guess here was obtained using the data-calibration procedure.

JBO_18_3_030503_f002.png

In both numerical and gelatin-phantom experiments, the total computational time that was reported for the traditional LM method was dependent on λ. The converged solution is the same for cases when λ lies between 0.01 and 100; outside this range, it is far away from the actual solution. Note that there are automated methods for choosing the optimal λ, but they result in an overhead time that is almost equivalent to performing the reconstruction scheme, making it computationally inefficient.13 The proposed method does not require regularization, in turn removing the unwarranted bias introduced in the solution in the traditional LM method. The proposed method solves the minimization problem in an iterative manner, where the number of iterations acts as indirect regularization of the solution space (similar to least-squares QR14). As the number of iterations is dependent on the stopping criterion, the choice of the same becomes important in this case. More importantly, the proposed method does not require any Jacobian, which might be a time-consuming process even for the image-guided DOT where the detection mechanism is based on CCD cameras resulting in M=1e7.15

In conclusion, we have proposed a new approach for image-guided diffuse optical tomographic image reconstruction that is based on the popular gradient-free Nelder–Mead simplex method. The reconstructed optical property images closely match the ones obtained using traditional image reconstruction methods. The developed algorithm for image-guided DOT is provided as an open source16 for interested readers.

References

1. 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).PTRMAD1364-503X http://dx.doi.org/10.1098/rsta.2011.0228 Google Scholar

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

3. B. Brooksbyet al., “Near-infrared (NIR) tomography breast image reconstruction with a priori structural information from MRI: algorithm development for reconstructing heterogeneities.” IEEE J. Sel. Top. Quant. Electron. 9(2), 199–209 (2003).IJSQEN1077-260X http://dx.doi.org/10.1109/JSTQE.2003.813304 Google Scholar

4. B. Brooksbyet al., “Imaging breast adipose and fibro-glandular tissue molecular signatures using hybrid MRI-guided near-infrared spectral tomography,” Proc. Natl. Acad. Sci. U. S. A. 103(23), 8828–8833 (2006).PNASA60027-8424 http://dx.doi.org/10.1073/pnas.0509636103 Google Scholar

5. M. SchweigerS. R. Arridge, “Optical tomographic reconstruction in a complex head model using a priori region boundary information,” Phys. Med. Biol. 44(11), 2703–2721 (1999).PHMBA70031-9155 http://dx.doi.org/10.1088/0031-9155/44/11/302 Google Scholar

6. R. P. K. JagannathP. K. Yalavarthy, “Approximation of internal refractive index variation improves image guided diffuse optical tomography of breast,” IEEE Trans. Biomed. Eng. 57(10), 2560–2563 (2010).IEBEAX0018-9294 http://dx.doi.org/10.1109/TBME.2010.2053368 Google Scholar

7. R. Asteret al., Parameter Estimation and Inverse Problems, Elsevier, New York (2005). Google Scholar

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

9. J. A. NelderR. Mead, “A simplex method for function minimization,” Comput. J. 7(4), 308–313 (1965).CMPJA60010-4620 http://dx.doi.org/10.1093/comjnl/7.4.308 Google Scholar

10. J. C. Lagariaset al., “Convergence properties of the Nelder-Mead simplex method in low dimensions,” SIAM J. Opt. 9(1), 112–147 (1998).SJOPE81095-7189 http://dx.doi.org/10.1137/S1052623496303470 Google Scholar

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

12. 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).CANMER0748-8025 http://dx.doi.org/10.1002/cnm.v25:6 Google Scholar

13. J. PrakashP. K. Yalavarthy, “A LSQR-type method provides a computationally efficient automated optimal choice of regularization parameter in diffuse optical tomography,” Med. Phys. 40(3), 033101 (2013).MPHYA60094-2405 http://dx.doi.org/10.1118/1.4792459 Google Scholar

14. C. C. PaigeM. A. Saunders, “LSQR: an algorithm for sparse linear equations and sparse least squares,” ACM Trans. Math. Softw. 8(1), 43–71 (1982).ACMSCU0098-3500 http://dx.doi.org/10.1145/355984.355989 Google Scholar

15. S. D. Koneckyet al., “Imaging complex structures with diffuse light,” Opt. Express 16(7), 5048–5060 (2008).OPEXFF1094-4087 http://dx.doi.org/10.1364/OE.16.005048 Google Scholar

© 2013 Society of Photo-Optical Instrumentation Engineers (SPIE)
Ravi Prasad K. Jagannath, Phaneendra K. Yalavarthy, "Efficient gradient-free simplex method for estimation of optical properties in image-guided diffuse optical tomography," Journal of Biomedical Optics 18(3), 030503 (20 March 2013). https://doi.org/10.1117/1.JBO.18.3.030503
JOURNAL ARTICLE
4 PAGES


SHARE
Back to Top