Joint sparsity-driven non-iterative simultaneous reconstruction of absorption and scattering in diffuse optical tomography Okkyun Lee and Jong Chul Ye∗ Bio Imaging and Signal Processing Lab., Dept. of Bio and Brain Engineering, Korea Adv. Inst. of Science and Technology (KAIST), 291 Daehak-ro Yuseong-gu, Daejon 305-701, South Korea ∗ [email protected]
Abstract: Some optical properties of a highly scattering medium, such as tissue, can be reconstructed non-invasively by diffuse optical tomography (DOT). Since the inverse problem of DOT is severely ill-posed and nonlinear, iterative methods that update Green’s function have been widely used to recover accurate optical parameters. However, recent research has shown that the joint sparse recovery principle can provide an important clue in achieving reconstructions without an iterative update of Green’s function. One of the main limitations of the previous work is that it can only be applied to absorption parameter reconstruction. In this paper, we extended this theory to estimate the absorption and scattering parameters simultaneously when the background optical properties are known. The main idea for such an extension is that a joint sparse recovery step gives us unknown fluence on the estimated support set, which eliminates the nonlinearity in an integral equation for the simultaneous estimation of the optical parameters. Our numerical results show that the proposed algorithm reduces the cross-talk artifacts between the parameters and provides improved reconstruction results compared to existing methods. © 2013 Optical Society of America OCIS codes: (170.0170) Medical optics and biotechnology; (170.6960) Tomography.
References and links 1. A. Villringer and B. Chance, “Non-invasive optical spectroscopy and imaging of human brain function,” Trends Neurosci. 20, 435–442 (1997). 2. D. A. Boas, D. H. Brooks, E. L. Miller, C. A. DiMarzio, M. Kilmer, R. J. Gaudette, and Q. Zhang, “Imaging the body with diffuse optical tomography,” IEEE Signal Process. Mag. 18, 57–75 (2001). 3. D. A. Benaron, S. R. Hintz, A. Villringer, D. Boas, A. Kleinschmidt, J. Frahm, C. Hirth, H. Obrig, J. C. van Houten, E. L. Kermit, W. F. Cheong, and D. K. Stevenson, “Noninvasive functional imaging of human brain using light,” J. Cereb. Blood Flow Metab. 20, 469–477 (2000). 4. V. Ntziachristos and B. Chance, “Probing physiology and molecular function using optical imaging: Applications to breast cancer,” Breast Cancer Res. 3, 41–46 (2001). 5. R. Weissleder and V. Ntziachristos, “Shedding light onto live molecular targets,” Nat. Med. 9, 123–128 (2003). 6. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15, R41–R93 (1999). 7. A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50, R1–R43 (2005). 8. V. Ntziachristos and R. Weissleder, “Experimental three-dimensional fluorescence reconstruction of diffuse media by use of a normalized Born approximation,” Opt. Lett. 26, 893–895 (2001).
#197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26589
9. V. A. Markel and J. C. Schotland, “Inverse problem in optical diffusion tomography. I. Fourier-Laplace inversion formulas,” J. Opt. Soc. Am. A 18, 1336–1347 (2001). 10. H. Jiang, K. D. Paulsen, U. L. Osterberg, B. W. Pogue, and M. S. Patterson, “Optical image reconstruction using frequency-domain data: Simulation and experiment,” J. Opt. Soc. Am. A 13, 253–266 (1996). 11. Y. Yao, Y. Wang, Y. Pei, W. Zhu, and R. L. Barbour, “Frequency domain optical imaging of absorption and scattering distributions by a Born iterative method,” J. Opt. Soc. Am. A 14, 325–342 (1997). 12. J. C. Ye, K. J. Webb, C. A. Bouman, and R. P. Millane, “Optical diffusion tomography using iterative coordinate descent optimization in a Bayesian framework,” J. Opt. Soc. Am. A 16, 2400–2412 (1999). 13. J. C. Ye, S. Y. Lee, and Y. Bresler, “Exact reconstruction formula for diffuse optical tomography using simultaneous sparse representation,” in International Symposium on Biomedical Imaging: From Nano to Macro, Paris, France (Institute of Electrical and Electronics Engineers, 2008), pp. 1621–1624. 14. O. K. Lee, J. M. Kim, Y. Bresler, and J. C. Ye, “Compressive diffuse optical tomography: Non-iterative exact reconstruction using joint sparsity,” IEEE Trans. Med. Imaging 30, 1129–1142 (2011). 15. J. Chen and X. Huo, “Theoretical results on sparse representations of multiple measurement vectors,” IEEE Trans. Signal Process. 54, 4634–4643 (2006). 16. J. M. Kim, O. K. Lee, and J. C. Ye, “Compressive MUSIC: Revisiting the link between compressive sensing and array signal processing,” IEEE Trans. Inf. Theory 58, 278–301 (2012). 17. A. Profio and G. Navarro, “Scientific basis of breast diaphanography,” Med. Phys. 16, 60–65 (1989). 18. S. Fantini, S. A. Walker, M. A. Franceschini, M. Kaschke, P. M. Schlag, and K. T. Moesta, “Assessment of the size, position, and optical properties of breast tumors in vivo by noninvasive optical methods,” Appl. Opt. 37, 1982–1989 (1998). 19. A. E. Cerussi, D. Jakubowski, N. Shah, F. Bevilacqua, R. Lanning, A. J. Berger, D. Hsiang, J. Butler, R. F. Holcombe, and B. J. Tromberg, “Spectroscopy enhances the information content of optical mammography,” J. Biomed. Opt. 7, 60–71 (2002). 20. N. F. Boyd, J. W. Byng, R. A. Jong, E. K. Fishell, L. E. Little, A. B. Miller, G. A. Lockwood, D. L. Tritchler, and M. J. Yaffe, “Quantitative classification of mammographic densities and breast cancer risk: Results from the Canadian national breast screening study,” J. Natl. Cancer Inst. 87, 670–675 (1995). 21. J. W. Byng, M. J. Yaffe, R. A. Jong, R. S. Shumak, G. A. Lockwood, D. L. Tritchler, and N. F. Boyd, “Analysis of mammographic density and breast cancer risk from digitized mammograms,” Radiographics 18, 1587–1598 (1998). 22. B. W. Pogue, S. Jiang, H. Dehghani, C. Kogel, S. Soho, S. Srinivasan, X. Song, T. D. Tosteson, S. P. Poplack, and K. D. Paulsen, “Characterization of hemoglobin, water, and NIR scattering in breast tissue: Analysis of intersubject variability and menstrual cycle changes,” J. Biomed. Opt 9, 541–552 (2004). 23. M. O’Leary, D. Boas, B. Chance, and A. Yodh, “Experimental images of heterogeneous turbid media by frequency-domain diffusion-photon tomography,” Opt. Lett. 20, 426–428 (1995). 24. J. C. Hebden, F. E. W. Schmidt, M. E. Fry, M. Schweiger, E. M. C. Hillman, D. T. Delpy, and S. R. Arridge, “Simultaneous reconstruction of absorption and scattering images by multichannel measurement of purely temporal data,” Opt. Lett. 24, 534–536 (1999). 25. T. O. McBride, B. W. Pogue, S. Jiang, U. L. Osterberg, and K. D. Paulsen, “Initial studies of in vivo absorbing and scattering heterogeneity in near-infrared tomographic breast imaging,” Opt. Lett. 26, 822–824 (2001). 26. J. Wang, S. D. Jiang, Z. Z. Li, R. M. diFlorio Alexander, R. J. Barth, P. A. Kaufman, B. W. Pogue, and K. D. Paulsen, “In vivo quantitative imaging of normal and cancerous breast tissue using broadband diffuse optical tomography,” Med. Phys. 37, 3715–3724 (2010). 27. B. J. Hoenders, “Existence of invisible nonscattering objects and nonradiating sources,” J. Opt. Soc. Am. A 14, 262–266 (1997). 28. S. Arridge and W. Lionheart, “Nonuniqueness in diffusion-based optical tomography,” Opt. Lett. 23, 882–884 (1998). 29. V. Markel and J. Schotland, “Inverse problem in optical diffusion tomography. II. Role of boundary conditions,” J. Opt. Soc. Am. A 19, 558–566 (2002). 30. L. V. Wang and H. Wu, Biomedical Optics: Principles and Imaging (John Wiley, 2007). 31. J. C. Ye, K. J. Webb, R. P. Millane, and T. J. Downar, “Modified distorted Born iterative method algorithm with an approximate Fr´echet derivative for optical diffusion tomography,” J. Opt. Soc. Am. A 16, 1814–1826 (1999). 32. M. E. Davies and Y. C. Eldar, “Rank awareness in joint sparse recovery,” IEEE Trans. Inf. Theory 58, 1135–1146 (2012). 33. S. F. Cotter, B. D. Rao, K. Engan, and K. Kreutz-Delgado, “Sparse solutions to linear inverse problems with multiple measurement vectors,” IEEE Trans. Signal Process. 53, 2477–2488 (2005). 34. D. Malioutov, M. Cetin, and A. Willsky, “A sparse signal reconstruction perspective for source localization with sensor arrays,” IEEE Trans. Signal Process. 53, 3010–3022 (2005). 35. D. P. Wipf and B. D. Rao, “An empirical Bayesian strategy for solving the simultaneous sparse approximation problem,” IEEE Trans. Signal Process. 55, 3704–3716 (2007). 36. J. M. Kim, O. K. Lee, and J. C. Ye, “Improving noise robustness in subspace-based joint sparse recovery,” IEEE Trans. Signal Process. 60, 5799–5809 (2012).
#197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26590
37. A. Chambolle, “An algorithm for total variation minimization and applications,” J. Math. Imaging Vision 20, 89–97 (2004). 38. M. Afonso, J. Bioucas-Dias, and M. Figueiredo, “An augmented Lagrangian approach to the constrained optimization formulation of imaging inverse problems,” IEEE Trans. Image Process. 20, 681–695 (2011). 39. Y. Pei, H. L. Graber, and R. L. Barbour, “Normalized-constraint algorithm for minimizing inter-parameter crosstalk in DC optical tomography,” Opt. Express 9, 97–109 (2001). 40. Y. Xu, X. Gu, T. Khan, and H. Jiang, “Absorption and scattering images of heterogeneous scattering media can be simultaneously reconstructed by use of DC data,” Appl. Opt. 41, 5427–5437 (2002). 41. Q. Fang and D. A. Boas, “Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units,” Opt. Express 17, 20178–20190 (2009). 42. B. Dogdas, D. Stout, A. F. Chatziioannou, and R. M. Leahy, “Digimouse: A 3D whole body mouse atlas from CT and cryosection data,” Phys. Med. Biol. 52, 577–587 (2007). 43. G. Alexandrakis, F. R. Rannou, and A. F. Chatziioannou, “Tomographic bioluminescence imaging by use of a combined optical-PET (OPET) system: A computer simulation feasibility study,” Phys. Med. Biol. 50, 4225– 4241 (2005). 44. S. A. Prahl, “Optical Properties Spectra,” Oregon Medical Laser Clinic, 2001, http://omlc.ogi.edu/spectra/index.html. 45. D. P. Huttenlocher, G. A. Klanderman, and W. J. Rucklidge, “Comparing images using the Hausdorff distance,” IEEE Trans. Pattern Anal. Mach. Intell. 15, 850–863 (1993).
1.
Introduction
The near-infrared (NIR) optical wavelength range (700 ∼ 1000nm) provides an important opportunity for biological imaging, as tissues are relatively transparent in this regime due to the low absorption rate of the primary absorbers [1, 2]. It is now well-known that NIR light can penetrate tissue up to a depth of several centimeters. Accordingly, the objective of DOT is to exploit this window of opportunity to reconstruct the optical properties of highly scattering biological tissue from scattered and attenuated optical flux measured at the boundary of the medium [3–5]. Extensive studies have been conducted for the theoretical modeling of photon propagation in tissue and to realize accurate reconstructions of the optical parameters [6, 7]. One of the main difficulties of DOT, however, is that the inverse problem is severely ill-posed and non-linear due to the diffusive nature of photon propagation and the non-linear coupling between the unknown optical fluence and parameters [6]. Conventional approaches to address this non-linearity have included linearization methods [8, 9] and iterative methods that update the Green’s function [10–12]. However, the linearization method fails when perturbation in the optical property is beyond the Born approximation limit. Furthermore, the computational overhead of iterative methods is often prohibitive, especially for 3-D imaging because updating the Green’s function requires multiple applications of 3-D partial differential equation (PDE) solvers or Monte Carlo simulation. Recently, a non-iterative reconstruction method was developed by our group based on the theory of joint sparse recovery from compressed sensing [13, 14]. This method exploits the facts that optical parameter variations occur in a relatively small area and their spatial positions are usually invariant during multiple optical illuminations. These properties convert the DOT problem into a joint sparse recovery problem, where we estimate a set of unknown signals with a common sparse support [15, 16]. It has been found that if the background optical parameters are known, the optical parameter perturbations can be reconstructed using a simple two-step procedure in which the unknown support set is initially estimated using a joint sparse recovery algorithm and the optical parameter variations are then calculated using a least-squares fitting method [13, 14]. However, the least-squares step limits the application of the theory to the recovery of only the perturbation in the absorption parameters. While optical absorption parameter perturbation can be observed during brain activity or in the angiogenesis of cancer [1, 4], scattering parameter perturbation is another important property. It has been reported that both absorption and scattering parameters change in breast tumors
#197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26591
[17, 18], and scattering parameter perturbation can provide information about breast composition and physiology changes that may be useful for diagnosing breast malignancy [19]. Moreover, the scattering amplitude and power correlate with the degree of radiographic density, which is correlated to the risk of developing breast cancer [20–22]. Accordingly, simultaneous reconstruction of both absorption and scattering coefficients change has been extensively studied using phantoms [23, 24] and in vivo experiments [25, 26], and theoretical analyses of the existence and uniqueness of solutions have also been reported [9, 27–29]. In this paper, we extended the results in our earlier work [13, 14] and propose a non-iterative method for the simultaneous reconstruction of both absorption and scattering parameter perturbation when the background optical properties are known. One of the key ideas for such an extension is that the least-squares step in the preceding studies [13, 14] can be easily replaced by a linear integral equation of the first kind, thanks to the self-recursive structure in the FoldyLax equation. In particular, as long as we can identify the support set for the absorption or scattering parameter variations accurately, unknown fluence can be estimated accurately, which eliminates the nonlinear terms in the forward integral equation. Moreover, the estimation problem for the support of absorption and scattering parameter variations can be addressed as a joint sparse recovery approach, which is similar to the earlier studies [13, 14]. Based on the new formulation, we demonstrate that the reconstruction quality is improved and, furthermore, that the cross-talk between two optical parameters can be reduced. The paper is organized as follows. After a brief review of conventional methods is given, we describe a theory of simultaneous absorption and scattering parameter estimation using joint sparsity in Section 2. The implementation of the algorithm and the simulation environment are described in detail in Section 3. Section 4 provides numerical results, which is followed by a discussion of the implications of the results in Section 5. Finally, the conclusion is provided in Section 6. 2.
Theory
In this section, a continuous-wave (CW) formulation is described to explain a non-iterative simultaneous reconstruction of perturbation of both absorption and scattering coefficients. An extension of the theory for frequency domain DOT formulation is also straightforward. 2.1.
Diffusion equation
Photon migration within biological tissues is usually modeled by the transport equation. However, when light scattering prevails over absorption, the propagation of light can be modeled by a diffusion equation. Let Ω be a domain filled with some turbid medium with ∂ Ω as its boundary. In a highly scattering medium with low absorption, the optical fluence u(r) at position r can be modeled by the following continuous-wave diffusion equation (DE) [29]: ∇ · D(r)∇u(r) − μa (r)u(r) = −S(r), r ∈ Ω (1) u(r) + nˆ · ∇u(r) = 0, r ∈ ∂Ω Here, μa (r) and D(r) = 1/3(μs (r) + μa (r)) are the absorption and diffusion coefficients, where μs (r) denotes the reduced scattering coefficient, and S(r) is the source intensity profile. Additionally, is a parameter related to the diffusion coefficient, dimension and reflection on the boundary [30], and nˆ denotes a vector normal to the measurement surface. In a DOT problem, the unknown optical parameters are modeled as D(r) = D0 (r) + δ D(r), μa (r) = μa0 (r) + δ μa (r), and μs (r) = μs0 (r) + δ μs (r), where μa0 (r) and μs0 (r) are the known background absorption and reduced scattering parameters, respectively. D0 (r) = 1/3(μs0 (r) + μa0 (r)) and δ D(r), δ μa (r), and δ μs (r) refer to the unknown perturbations, which we need to
#197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26592
estimate. Note that the known background optical parameter does not need to be constant, and our theory holds for any spatially varying known background. The diffusion equation can be then transformed as (∇ · D0 (r)∇ − k(r)2 )u(r) = [δ μa (r) − ∇ · δ D(r)∇]u(r) − S(r), r ∈ Ω (2) u(r) + nˆ · ∇u(r) = 0, r ∈ ∂ Ω, where the diffusion wave number is given by k(r) =
δ μs (r) μa0 (r) and δ D(r) = − δμμa a (r)+ (r)+μ (r) D(r). 0
Then, Eq. (2) can be equivalently represented using the integral equation: u(r) = u0 (r) −
Ω
s0
G0 (r, r )[δ μa (r ) − ∇ · δ D(r )∇]u(r )dr , r ∈ Ω,
(3)
where r denotes a position in volume Ω and the background Green’s function G0 (r, r ) satisfies (∇ · D0 (r)∇ − k(r)2 )G0 (r, r ) = −δ (r − r ), r ∈ Ω (4) r ∈ ∂ Ω, G0 (r, r ) + nˆ · ∇G0 (r, r ) = 0, and the incidence fluence u0 (r) is given by u0 (r) =
Ω
G0 (r, r )S(r )dr .
(5)
Now, assuming that δ D(r) is zero in the neighborhood around the boundary, we have 0
= = =
∂Ω
Ω Ω
nˆ · [δ D(r )G0 (r, r )∇u(r )]dr
∇ · [δ D(r )G0 (r, r )∇u(r )]dr
δ D(r )∇G0 (r, r ) · ∇u(r )dr +
Ω
G0 (r, r )∇ · δ D(r )∇u(r )dr ,
(6)
where we use the divergence theorem to derive the second equality. Therefore, we have u(r) = u0 (r) −
Ω
G0 (r, r )u(r )δ μa (r )dr −
Ω
∇G0 (r, r ) · ∇u(r )δ D(r )dr , r ∈ ∂ Ω.
(7)
Note that the resulting forward integral equation that maps the unknown optical parameter perturbation to the boundary measurement is nonlinear, as u(r ) depends on the unknown optical parameter perturbation (δ μa , δ D). 2.2.
Proposed method
In a conventional iterative approach such as the distorted Born iterative method (DBIM), Eq. (7) is linearized as u(r; x) ≈ u(r; x(k) ) − −
Ω Ω
G0 (r, r ; x(k) )u(r ; x(k) )δ μa (r )dr ∇G0 (r, r ; x(k) ) · ∇u(r ; x(k) )δ D(r )dr ,
(8)
(k)
where x(k) = [μa , D(k) ]T denotes the k-th estimate of the background optical parameter, and [δ μa , δ D]T denotes the unknown perturbation estimated from the k-th estimate of the background optical parameters. Note that the fluence u(r ; x) is linearized with respect to the back#197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26593
ground optical parameter x(k) . After the linearization, [δ μa , δ D]T is estimated and we perform (k+1) (k+1) (k) , D(k+1) ]T , where μa = μa + δ μa , and another linearization of u(r ; x) at x(k+1) = [μa (k+1) (k) = D + δ D. This procedure is repeated until convergence. It has been shown that the D linearized integral equation corresponds to the Fr´echet derivative and the DBIM is equivalent to the Levenberg-Marquardt algorithm [31]. However, the main problem with this type of iterative linearization is that we need to calculate the Green’s function multiple times, which requires multiple applications of computationally extensive forward solvers. Hence, in some conventional algorithms, only one linearization step is performed (i.e., Born approximation), which often degrades the performance. To overcome these drawbacks, we developed a novel non-iterative method for the simultaneous reconstruction of both absorption and scattering coefficients when the background optical properties are known. The main idea comes from the observation that the absorption or scattering properties change locally against the background (e.g., due to localized angiogenesis for cancer or brain activation) and that the reconstruction problem can be converted into a joint sparse recovery. Figure 1 illustrates such a concept. The target area Ωt represents areas where
Fig. 1. Sparsely distributed perturbations in optical properties for simultaneous DOT reconstruction.
the absorption or scattering properties change. We assumed that the area Ωt is sparse compared to the total field of view (FOV) Ω. This is usually true in practice because, for example, tumors are only sparsely distributed. Furthermore, we assumed that there are N distinct illumination patterns and let u(r; l) be the optical fluence from the l-th illumination pattern. Then, Eq. (3) can be rewritten as u(r; l)
= u0 (r; l) − = u0 (r; l) −
Ω
Ωt
G0 (r, r )X(r ; l)dr G0 (r, r )X(r ; l)dr , l = 1, 2, · · · , N,
(9)
where X(r; l) = [δ μa (r) − ∇ · δ D(r)∇]u(r; l) is the induced current from the l-th illumination. In Eq. (9), the second equality comes from the following:
δ μa (r) = 0, δ D(r) = 0, ∇δ D(r) = 0, r ∈ / Ωt .
(10)
Once the support Ωt is estimated, X(r; l) can easily be obtained by means of linear inversion from Eq. (9). Then, the unknown optical fluence u(r; l) can be calculated by restricting r ∈ Ωt . More specifically, we have the following Foldy-Lax equation for fluence estimation, as derived
#197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26594
from Eq. (9) by restricting r ∈ Ωt : u(r; ˆ l) = u0 (r; l) −
Ωt
ˆ ; l)dr , r ∈ Ωt . G0 (r, r )X(r
(11)
ˆ l) denotes the estimated induced current for the l-th illumination. Then, using u(r; Here, X(r; ˆ l), δ μa and δ D can be reconstructed from the following linear integral equation: u(r; l) = u0 (r; l) −
Ωt
G0 (r, r )u(r ˆ ; l)δ μa (r )dr −
Ωt
∇G0 (r, r ) · ∇u(r ˆ ; l)δ D(r )dr ,
(12)
l = 1, 2, · · · , N. ˆ ; l) Note that the nonlinear coupling due to u(r ; l) in Eq. (7) becomes now linear because u(r and ∇u(r ˆ ; l) are calculated from Eq. (11). Note that the major difference from the conventional methods is that as the nonlinearly coupled optical fluence u(r; l) for r ∈ Ωt is now estimated, the iterative procedure for updating Green’s function becomes unnecessary. Therefore, one of the most important steps for the use of Eq. (12) is the joint support estimation of Ωt from multiple illumination patterns to obtain an accurate estimation of the unknown fluence u(r; l). Another important issue is to estimate the optical parameters using Eq. (12). In our previous work [13, 14], we assumed that the perturbation occurs in the absorption parameter alone, which is calculated by least-squares fitting at each position r ∈ Ωt from the estimated ˆ l) = δ μa (r)u(r; ˆ l) and u(r; ˆ l) values for l = 1, 2, · · · , N. However, the presence of the perX(r; turbation in the scattering parameter inhibited the use of this step. In the following section, we describe how to address these issues. 2.3.
Joint support recovery
We assumed that the induced current is described by piecewise constant or spline approximation as: k
X(r; l) = ∑ X(r(i) ; l)b(r, r(i) ), l = 1, 2, · · · , N,
(13)
i=1
where b(r, r(i) ) denotes the basis function. The 3-D locations {r(i) }ki=1 ∈ Ωt of the k-targets are selected from n-possible locations {ri }ni=1 ∈ Ω, where k n. Then, we have the following relationship after substituting Eq. (13) into the induced current term in Eq. (9), u(r; l) = u0 (r; l) −
Ωt
G0 (r, r )X(r ; l)dr , l = 1, 2, · · · , N,
k
= u0 (r; l) − ∑ G˜ 0 (r, r(i) )X(r(i) ; l),
(14)
i=1
where G˜ 0 (r, r(i) ) =
Ωt
G0 (r, r )b(r , r(i) )dr .
(15)
By collecting the measurement of Eq. (14) at detector positions {rdi }m i=1 for all N illuminations, the imaging problem can be represented as the following matrix equation: Y = AX.
(16)
#197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26595
In this equation, Y ∈ Rm×N is scattered fluence, and A ∈ Rm×n and X ∈ Rn×N are the sensing matrix and the induced current matrix, respectively. Let the component of the i-th row and j-th column of matrix A be defined as Ai j . Then, Yi j = u0 (rdi ; j) − u(rdi ; j), Ai j = G˜ 0 (rdi , r j ). Note that the rows of X are only non-zero when ri ∈ {r(i) }ki=1 . Now, we define the row-diversity measure ||X||0 that counts the number of rows in X ∈ Rn×N with nonzero elements. Then, the simultaneous reconstruction of the DOT problem from multiple illumination cases, where the support of the optical parameter variation is sparse, can be stated as the following joint sparse recovery problem [14, 15]: min ||X||0 , subject to ||Y − AX||F ≤ ε .
(17)
X
Here, ||X||F is the Frobenius norm of matrix X. Note that when δ D(r(i) ) = 0 for i = 1, 2, · · · , k, Eq. (17) becomes precisely the same joint sparse recovery problem in our previous work [13, 14]. Even if δ D(r(i) ) = 0, in terms of implementation, the support estimation step is identical to that of the previous work, since δ D is recovered later from the unknown X. One of the main advantages of this joint sparse recovery formulation is that it enables a further reduction in the number of required measurements [15, 16, 32], as the number of measurements required per sensor must account for the minimal features unique to that sensor. There exist a variety of joint sparse recovery algorithms, including ones from our group [16, 33–36]. The central theme of these algorithms is to exploit the diversity of multiple measurements to provide reconstruction from fewer samples or from a low SNR. Among the various algorithms, in this paper we used the M-SBL algorithm [35], given that we did not need to know the sparsity level and because the algorithm is known to be robust to noise. 2.4.
Optical parameter estimation ˆ
ˆ t and the calculated optical fluence u(r ˆ t, ˆ i ) at ri ∈ Ω Based on the estimated supports {r(i) }ki=1 ∈ Ω the unknown optical parameters can be reconstructed by solving the following discretized linear inverse problem derived from Eq. (12): ⎤ ⎡ δ μa,1 ⎤ ⎢ .. ⎥ ⎡ (1) ⎤ ⎡ (1) (1) ⎢ . ⎥ y Aμa AD ⎥ ⎢ ⎥ ⎢ δμ ˆ ⎥ ⎢ .. ⎥ ⎢ . . a, k ⎥ ⎥ = Ax. ⎢ ⎢ .. ⎦ ⎢ y = ⎣ . ⎦ = ⎣ .. (18) δ D1 ⎥ ⎥ ⎢ (N) (N) (N) ⎢ . ⎥ y Aμa AD ⎣ .. ⎦ δ Dkˆ (l)
(l)
(l)
ˆ ( j) ; l)δ v, AD,i j = ∇G0 (rdi , r( j) ) · Here, yi = u0 (rdi ; l) − u(rdi ; l), Aμa,i j = G0 (rdi , r( j) )u(r ∇u(r ˆ ( j) ; l)δ v, and δ v is the volume of the discretized voxel. In addition, δ μa,i = δ μa (r(i) ) and δ Di = δ D(r(i) ) denote the unknown optical parameter perturbations. Due to the restriction of ˆ t , the number of unknowns in the discretized voxel domain is significantly reduced. the ROI to Ω However, the ill-posedness of the problem is still not completely overcome by the restriction since the ill-posedness of the DOT problem originates from the underlying continuous operator. In fact, it is well-known that the linear mapping from DOT is compact mapping, whose singular values are clustered around zero, which makes the inverse mapping very ill-posed. Therefore, to solve Eq. (18), we convert it into the following constrained optimization problem with L1 and the total variation regularization under the assumption of the sparsity and smoothness of
#197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26596
both absorption and scattering perturbations: arg min ||δ μa ||1 + ||δ D||1 + λ ||δ μa ||TV + λ ||δ D||TV x
(19)
||Ax − y||2 ≤ ε , amin ≤ δ μai ≤ amax , Dmin ≤ δ Di ≤ Dmax .
subject to
ˆ
ˆ
In this equation, x = [δ μaT , δ DT ]T ∈ R2k×1 , A = [Aμa , AD ] ∈ RmN×2k , and ||x||1 and ||x||TV are L1 penalty and the total variation of x [37], respectively. Additionally, [amin , amax ] and [Dmin , Dmax ] denote the ranges of the absorption and diffusion parameter variations, respectively. 3. 3.1.
Method Implementation
In this section, we describe the implementation procedure for the proposed method and conventional algorithm in detail. In Eq. (19), λ is set to 1 for the proposed method, whereas it is reduced to 0.03 ∼ 0.05 for conventional methods since we found that the more weight in TV provides inaccurate absorption parameter reconstruction in the conventional methods. Therefore, we wanted to adjust the parameter to make the comparison fair. To solve the aforementioned constrained optimization problem, we exploited a constrained split augmented Lagrangian shrinkage algorithm (C-SALSA) [38] after considering that it is quite effective in dealing with various regularization schemes and constraints. In particular, Chambolle’s algorithm was employed to deal with TV regularization and a soft-threshold function with L1 penalty as in the C-SALSA algorithm [37, 38]. We normalized the sensing matrix in Eqs. (17) and (19) to have the unit norm for each column. For a fair comparison, C-SALSA was also used for conventional methods. More specifically, the same constrained optimization problem in Eq. (19) can also be applied to DBIM or a linearized approach with (l) the sensing matrix A = [Aμa , AD ] ∈ RmN×2n , where Aμa,i j = G0 (rdi , r j ; x(k) )u(r j ; l, x(k) )δ v and (l)
AD,i j = ∇G0 (rdi , r j ; x(k) ) · ∇u(r j ; l, x(k) )δ v for k-th iteration. When k = 1, we use the known background optical parameters as x(1) = [μa0 , D0 ]T . However, note that the size of the sensing matrix for the proposed method is much smaller than that of the conventional methods due to the restriction of the problem within the estimated joint support. The pseudocode implementation of C-SALSA is described in Algorithm 1. In Algorithm 1, Ψτ gi is the Moreau proximal mapping of τ gi [38], where g1 (x) = ||x||1 , g2 (x) = ||x||TV , g3 = ιE(ε ,y) , and g4 = ιB(a,D) . ιE(ε ,y) and ιB(a,D) are the indicator functions of an ε -radius Euclidean ball centered at y and a Box-constraint which depends on the range of x, respectively. As a stopping criterion, we used the relative change of the cost function in Eq. (19) and performed the algorithms until |(Ck −Ck−1 ) /Ck | < 5 × 10−4 is satisfied, where Ck is the cost function at k-th iteration. All simulation parameters are summarized in Table 1. (1) In Table 1, s¯1 (1) is the mean value of the s1 , which is defined in Algorithm 1, and εi and yi are the values of ε and measurements y for the i-th iteration of the DBIM method. The scale factor for ε was forced to increase during iterations in the case of DBIM, since we found that a fixed value often makes the iteration not converge. These parameters were found manually for each method by choosing the optimal parameter among c = {8, 4, 2, 1, 1/2, 1/4, 1/8} for τ = cs¯1 (1) and c = {0.01, 0.03, 0.05, 0.075, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6} for ε1 = c||y1 ||2 , and we did an additional search for ε2 and ε3 from 0.7 ∼ 1.0 to make the DBIM converge. Since the sensing matrices for the proposed method and conventional methods are different, the selected optimal parameters are distinct, as shown in Table 1, for each method.
#197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26597
Algorithm 1 Pseudocode implementation of C-SALSA for DOT problem in Eq. (19). (i)
(i)
Set k = 0, H (1,4) = I, H (2) = λ I, H (3) = A, and choose τ > 0, v0 , d0 , for i = 1, · · · , 4. repeat for i = 1, . . . , 4 do (i) (i) (i) 4: ζk = vk + dk 5: end for −1
T
T 4 4 ( j) ( j) ( j) 6: uk+1 = ∑ H H ∑ H ( j) ζk 1: 2: 3:
j=1
7: 8:
j=1
for i = 1, . . . , 4 do (i) (i) (i) (i) vk+1 = Ψτ gi sk , where sk = H (i) uk+1 − dk (i)
(i)
(i)
dk+1 = dk − H (i) uk+1 + vk+1 end for k = k+1 12: until some stopping criterion is satisfied 9: 10: 11:
Table 1. Parameters used in C-SALSA algorithm.
Proposed Method
Simple Target Simulation
Mouse Phantom Simulation
λ =1
λ =1
τ = 0.25s¯1
(1)
ε = 0.01||y||2
Conventional Methods
3.2.
τ = 2s¯1 (1) ε = 0.01||y||2
λ = 0.05
λ = 0.03
τ = 4s¯1 (1)
τ = 4s¯1 (1)
ε1 = 0.01||y1 ||2
ε1 = 0.03||y1 ||2
ε2 = 0.9||y2 ||2
ε2 = 0.7||y2 ||2
ε3 = ||y3 ||2
ε3 = 0.9||y3 ||2
Simulation geometry I : simple targets
First, we considered two simple targets immersed in a rectangular homogeneous medium, as in Fig. 2. Sixteen line sources with a 2mm space in total were used, and a 16 × 16 array detector with a pitch of 1.33mm measured the optical flux at the bottom of the medium. The medium was discretized with 2 × 2 × 2mm3 voxels and the region of interest (ROI) for reconstruction was 22 × 26 × 22mm3 in size. The optical parameters for the background were set to μa = 0.015mm−1 , and the μs were set to be equal to 1mm−1 . We conducted four simulation studies of this geometry, as described in Table 2, as CASE A to CASE D. In CASE A and CASE B, each target had optical perturbation in both absorption and scattering coefficients, at low and high contrast, respectively. On the other hand, in CASE C and CASE D, target 1 had only absorption changes, whereas target 2 had scattering changes. These cases were designed to evaluate the cross-talk artifacts [39, 40]. For an accurate simulation, we solved the diffusion equation using a GPU-accelerated Monte-Carlo simulation [41]. Noise with a SNR of 40dB was added to the scattered fluence. The results were obtained by averaging 10 runs with independent noise realization for each case.
#197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26598
Fig. 2. Simulation geometry for simple targets.
Table 2. Optical properties for various cases. CASE B
CASE C
CASE D
(low contrast)
(high contrast)
(low contrast)
(high contrast)
δ μa
δ μs
Target 1
0.0015
Target 2
0.0015
Optical parameters in
3.3.
mm−1
CASE A
δ μa
δ μs
δ μs
δ μa
δ μa
0.05
0.03
1.0
0.05
0.03
1.0
δ μs
0.0015
0.0
0.03
0.0
0.0
0.05
0.0
1.0
Simulation geometry II : mouse phantom
Next, we considered more complicated simulation studies using the Digimouse computational mouse phantom [42]. The goal of this study was to localize tumors within a mouse and to recover any optical parameter changes in the tumors. The simulation geometry is described in Fig. 3(a). The mouse phantom was placed in a rectangular shape of matching fluid that has the same optical parameters as mouse skin. Sixteen line sources in total with a space of 2mm were again used, and a 20×20 array detector with a pitch of 1.5mm served to measure the optical flux at the bottom of the medium. The medium was discretized with 2 × 2 × 2mm3 voxels. The ROI for reconstruction is the area between the source and detector with a size of 38 × 30 × 22mm3 , where three tumors were placed, as illustrated in Fig. 3(b). Various optical parameters of the mouse organs with a wavelength of 800nm [43, 44], summarized in Table 3, were used. For simplicity, the optical parameters corresponding to bones, fat, muscle, and other parts were assumed to be the same as skin. The tumor in the liver is denoted as target 1, which has both absorption and scattering coefficients changes, whereas targets 2 and 3, corresponding to tumors in kidneys, have changes only in their absorption and scattering coefficients, respectively, as described in Table 4. Except for the tumors, all spatially varying optical parameters were assumed to be known. We solved the diffusion equation using the GPU-accelerated Monte-Carlo simulation [41], and noise with a SNR of 40dB was added to the scattered fluence. The results were also obtained by averaging 10 runs with independent noise realizations. 4.
Results
We compared the reconstruction results of the proposed method to those of the conventional iterative (DBIM) and linearization methods. The performance was evaluated quantitatively using
#197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26599
Fig. 3. Simulation geometry of the mouse phantom.
Table 3. Optical properties of the mouse [43, 44].
Tissue type
Skin
Heart
μa in mm−1
0.0259
0.0218
μs in mm−1
1.8258
0.748
Lung
Liver
Stomach
Kidney
0.0636
0.126
0.0058
0.0249
1.9789
0.5629
1.2098
1.7237
Table 4. Optical properties of the targets.
optical parameters in mm−1
δ μa
δ μs
Target 1
0.0504
0.8444
Target 2
0.0423
0.0
Target 3
0.0
0.8619
the following relative mean square error: MSE(x) =
||xtrue − xrecon ||2 . ||xtrue ||2
(20)
Although the final equation of Eq. (19) depends on δ D, our practical interest is δ μs . Therefore, δ μs (r) we calculated δ μs using the relationship of δ D(r) = − δμμa a (r)+ (r)+μ (r) D(r) based on the estimated 0
s0
values of δ μa and δ D, as shown in the following: δ D(r) + δ μa (r) . δ μs (r) = − 3D0 (r)[δ D(r) + D0 (r)]
(21)
We also used the Hausdorff distance [45] to measure the morphological difference. The number of iterations for the DBIM method was set to 3.
#197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26600
4.1.
Simple target
The reconstruction results of the simple target simulation are summarized in Table 5. The MSE values in Table 5 are the average values of the MSE values from multiple noise simulation results, and the values in the parenthesis denote the corresponding standard deviation. The relative error for the proposed method is competitive or lower than those of DBIM or the linearized approach. For the case of the linearized approach, the mean MSE values and their standard deviations were the smallest when the perturbation in optical parameters are small (CASE A). However, the algorithm broke down for large optical perturbation (CASE B) or in the presence of the simultaneous perturbation of absorption and scattering (CASE C and CASE D). The results of DBIM were slightly improved compared to those of the linearized approach for absorption coefficient cases. On the other hand, the results of the proposed method for the reconstruction of the absorption coefficients were robust for all cases, and the MSE values for scattering coefficients showed monotonic increasing behaviors from CASE A to CASE D. These phenomena can be easily confirmed from Figs. 4 and 5, which show the cross-sectional image of the reconstructed images for δ μa and δ μs for each case (The horizontal and coronal sections centered on each target are indicated in the inset image in Figs. 4 and 5). Moreover, in the reconstruction images, our method reduced the cross-talk artifacts that were observed in the conventional methods, especially in the reconstructed scattering coefficients. The average reconstruction time of various methods for the simple target simulation are summarized in Table 6 (using a PC with CPU : core i7 sandy bridge and GPU : GTX 560 Geforce series). The support estimation step using the M-SBL algorithm was also added in the calculation of the run time of the proposed method. Even though the proposed method requires this additional step, the resulting computational time for the proposed method is even smaller compared to the linearized approach. In Table 6, we also describe the computation time for each part of the algorithm. Note that the total computation time includes the time required to construct the A matrix, which was about 1.2 seconds. Table 5. The MSE values of the reconstruction using various methods.
proposed method
linearized approach
DBIM
CASE A
MSE(δ μa )
0.149 (± 0.009)
0.360 (± 0.004)
0.359 (± 0.004)
(low contrast)
MSE(δ μs )
0.165 (± 0.009)
0.176 (± 0.006)
0.176 (± 0.006)
CASE B
MSE(δ μa )
0.134 (± 0.007)
0.547 (± 0.069)
0.525 (± 0.076)
(high contrast)
MSE(δ μs )
0.247 (± 0.008)
0.319 (± 0.081)
0.352 (± 0.134)
CASE C
MSE(δ μa )
0.143 (± 0.050)
0.630 (± 0.045)
0.628 (± 0.046)
(low contrast)
MSE(δ μs )
0.328 (± 0.039)
0.297 (± 0.068)
0.301 (± 0.080)
CASE D
MSE(δ μa )
0.123 (± 0.026)
0.715 (± 0.048)
0.696 (± 0.046)
(high contrast)
MSE(δ μs )
0.488 (± 0.028)
0.456 (± 0.052)
0.523 (± 0.045)
4.2.
Mouse phantom
The reconstruction results of the mouse phantom using various methods are summarized in Table 7. The MSE values for the proposed method are lower than those for DBIM or the linearized approach for both the absorption and scattering coefficients. Table 8 summarizes the Hausdorff distance results [45] between the true supports and the estimated supports of each tumor in the liver and kidneys. Both of these quantitative measures show the advantages of the proposed #197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26601
Fig. 4. The cross-sections of the reconstructed images using various methods for simple target simulations (CASE A and CASE B). Table 6. Average run time of the reconstruction using various methods.
proposed method Run time (sec)
linearized approach
M-SBL
C-SALSA
Total
C-SALSA
Total
0.34
0.22
1.92
5.71
7.08
DBIM 1537.5
method. Figure 6 illustrates the location of the tumor volume thresholded with the same number of non-zero supports for each instance of δ μa and δ μs (29 and 30), respectively. The conventional methods exhibited cross-talk between μa and μs , whereas the proposed method did not have such artifacts. Table 7. MSE for the reconstruction results of the mouse phantom using various methods.
proposed method
MSE(δ μa )
MSE(δ μs )
0.326 (± 0.049)
0.624 (± 0.121)
linearized approach
0.491 (± 0.004)
0.770 (± 0.008)
DBIM
0.477 (± 0.005)
0.775 (± 0.008)
#197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26602
Fig. 5. The cross-sections of the reconstructed images using various methods for simple target simulations (CASE C and CASE D).
Fig. 6. The 3-D simultaneous reconstruction results of tumors with optical parameter variations δ μa and δ μs embedded in the mouse phantom.
#197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26603
Table 8. Hausdorff distance in mm for various reconstruction methods. proposed method Hausdorff distance in mm Liver Kidney (left) Kidney (right)
5.
Linearized approach
DBIM
δ μa
δ μs
δ μa
δ μs
δ μa
δ μs
2.0 (± 0.0) 2.0 (± 0.0) -
2.0 (± 0.0) 2.0 (± 0.0)
4.0 (± 0.0) 2.83 (± 0.0) -
2.08 (± 0.263) 2.0 (± 0.0)
3.0 (± 1.054) 2.83 (± 0.0) -
2.08 (± 0.263) 2.0 (± 0.0)
Discussion
In the existing non-iterative method for absorption imaging [14], the second step to calculate the absorption coefficient was performed in a voxel-by-voxel manner using the least-squares fitting approach. However, this step was not valid in the simultaneous reconstruction problem due to the differential operator applied to the diffusion coefficient. In this paper, we demonstrated that the unknown optical fluence can be estimated, regardless of the differential equation form, thanks to the recursive nature of the Foldy-Lax equation. Traditionally, the cross-talk issue was considered unavoidable during the simultaneous reconstruction of DOT problems with CW imaging. Using experimental data, the reconstruction of both absorption and scattering coefficient perturbations was shown to be impossible in CW cases [28]. To overcome this limitation, Pei et al. [39] proposed the normalized-constraint algorithm, and Xu et al. [40] used both normalization and regularized techniques. However, the proposed method is totally different from these works since we utilized the joint sparsity. While the conventional methods require only the optical parameter estimation step, our method needs the additional task of finding the support set on which the optical parameter estimation is conducted. Therefore, the joint sparse recovery step is important for the accurate localization of the target and subsequent calculation of optical parameters. In our algorithm, the main two computational building parts, the M-SBL and C-SALSA step, has complexity that follows the dimension of the sensing matrix. Recall that in the proposed method, the dimension of the matrix for C-SALSA is reduced due to the restriction of the unknown within the estimated support. Moreover, the complexity of M-SBL depends on the number of detectors. Unlike the conventional method, whose complexity depends on the number of source-detector pairs, the significant reduction of the sensing matrix dimension makes the M-SBL and C-SALSA computationally very efficient. Hence, the combination of M-SBL and C-SALSA results in the reduction of the computational burden compared to the linearized approaches. 6.
Conclusion
In this paper, we proposed a non-iterative method for the simultaneous reconstruction of both absorption and scattering coefficients in diffuse optical tomography. By exploiting multiple illumination patterns and the sparseness of the target area, a DOT inverse problem was converted to a joint sparse recovery problem. Moreover, the unknown optical parameters were calculated without an approximation or an iterative update of the Green’s function, thanks to the recursive nature of the Foldy-Lax equation. The proposed method was validated with various simulation studies, showing an improvement over conventional methods. Furthermore, our algorithm reduced the cross-talk artifacts during simultaneous reconstruction from CW data. Acknowledgments This research was supported by a grant from the Korea Science and Engineering Foundation (KOSEF) funded by the Korean government (MEST) (No.2012-0000173).
#197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26604
Abstract: Some optical properties of a highly scattering medium, such as tissue, can be reconstructed non-invasively by diffuse optical tomography (DOT). Since the inverse problem of DOT is severely ill-posed and nonlinear, iterative methods that update Green’s function have been widely used to recover accurate optical parameters. However, recent research has shown that the joint sparse recovery principle can provide an important clue in achieving reconstructions without an iterative update of Green’s function. One of the main limitations of the previous work is that it can only be applied to absorption parameter reconstruction. In this paper, we extended this theory to estimate the absorption and scattering parameters simultaneously when the background optical properties are known. The main idea for such an extension is that a joint sparse recovery step gives us unknown fluence on the estimated support set, which eliminates the nonlinearity in an integral equation for the simultaneous estimation of the optical parameters. Our numerical results show that the proposed algorithm reduces the cross-talk artifacts between the parameters and provides improved reconstruction results compared to existing methods. © 2013 Optical Society of America OCIS codes: (170.0170) Medical optics and biotechnology; (170.6960) Tomography.
References and links 1. A. Villringer and B. Chance, “Non-invasive optical spectroscopy and imaging of human brain function,” Trends Neurosci. 20, 435–442 (1997). 2. D. A. Boas, D. H. Brooks, E. L. Miller, C. A. DiMarzio, M. Kilmer, R. J. Gaudette, and Q. Zhang, “Imaging the body with diffuse optical tomography,” IEEE Signal Process. Mag. 18, 57–75 (2001). 3. D. A. Benaron, S. R. Hintz, A. Villringer, D. Boas, A. Kleinschmidt, J. Frahm, C. Hirth, H. Obrig, J. C. van Houten, E. L. Kermit, W. F. Cheong, and D. K. Stevenson, “Noninvasive functional imaging of human brain using light,” J. Cereb. Blood Flow Metab. 20, 469–477 (2000). 4. V. Ntziachristos and B. Chance, “Probing physiology and molecular function using optical imaging: Applications to breast cancer,” Breast Cancer Res. 3, 41–46 (2001). 5. R. Weissleder and V. Ntziachristos, “Shedding light onto live molecular targets,” Nat. Med. 9, 123–128 (2003). 6. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15, R41–R93 (1999). 7. A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50, R1–R43 (2005). 8. V. Ntziachristos and R. Weissleder, “Experimental three-dimensional fluorescence reconstruction of diffuse media by use of a normalized Born approximation,” Opt. Lett. 26, 893–895 (2001).
#197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26589
9. V. A. Markel and J. C. Schotland, “Inverse problem in optical diffusion tomography. I. Fourier-Laplace inversion formulas,” J. Opt. Soc. Am. A 18, 1336–1347 (2001). 10. H. Jiang, K. D. Paulsen, U. L. Osterberg, B. W. Pogue, and M. S. Patterson, “Optical image reconstruction using frequency-domain data: Simulation and experiment,” J. Opt. Soc. Am. A 13, 253–266 (1996). 11. Y. Yao, Y. Wang, Y. Pei, W. Zhu, and R. L. Barbour, “Frequency domain optical imaging of absorption and scattering distributions by a Born iterative method,” J. Opt. Soc. Am. A 14, 325–342 (1997). 12. J. C. Ye, K. J. Webb, C. A. Bouman, and R. P. Millane, “Optical diffusion tomography using iterative coordinate descent optimization in a Bayesian framework,” J. Opt. Soc. Am. A 16, 2400–2412 (1999). 13. J. C. Ye, S. Y. Lee, and Y. Bresler, “Exact reconstruction formula for diffuse optical tomography using simultaneous sparse representation,” in International Symposium on Biomedical Imaging: From Nano to Macro, Paris, France (Institute of Electrical and Electronics Engineers, 2008), pp. 1621–1624. 14. O. K. Lee, J. M. Kim, Y. Bresler, and J. C. Ye, “Compressive diffuse optical tomography: Non-iterative exact reconstruction using joint sparsity,” IEEE Trans. Med. Imaging 30, 1129–1142 (2011). 15. J. Chen and X. Huo, “Theoretical results on sparse representations of multiple measurement vectors,” IEEE Trans. Signal Process. 54, 4634–4643 (2006). 16. J. M. Kim, O. K. Lee, and J. C. Ye, “Compressive MUSIC: Revisiting the link between compressive sensing and array signal processing,” IEEE Trans. Inf. Theory 58, 278–301 (2012). 17. A. Profio and G. Navarro, “Scientific basis of breast diaphanography,” Med. Phys. 16, 60–65 (1989). 18. S. Fantini, S. A. Walker, M. A. Franceschini, M. Kaschke, P. M. Schlag, and K. T. Moesta, “Assessment of the size, position, and optical properties of breast tumors in vivo by noninvasive optical methods,” Appl. Opt. 37, 1982–1989 (1998). 19. A. E. Cerussi, D. Jakubowski, N. Shah, F. Bevilacqua, R. Lanning, A. J. Berger, D. Hsiang, J. Butler, R. F. Holcombe, and B. J. Tromberg, “Spectroscopy enhances the information content of optical mammography,” J. Biomed. Opt. 7, 60–71 (2002). 20. N. F. Boyd, J. W. Byng, R. A. Jong, E. K. Fishell, L. E. Little, A. B. Miller, G. A. Lockwood, D. L. Tritchler, and M. J. Yaffe, “Quantitative classification of mammographic densities and breast cancer risk: Results from the Canadian national breast screening study,” J. Natl. Cancer Inst. 87, 670–675 (1995). 21. J. W. Byng, M. J. Yaffe, R. A. Jong, R. S. Shumak, G. A. Lockwood, D. L. Tritchler, and N. F. Boyd, “Analysis of mammographic density and breast cancer risk from digitized mammograms,” Radiographics 18, 1587–1598 (1998). 22. B. W. Pogue, S. Jiang, H. Dehghani, C. Kogel, S. Soho, S. Srinivasan, X. Song, T. D. Tosteson, S. P. Poplack, and K. D. Paulsen, “Characterization of hemoglobin, water, and NIR scattering in breast tissue: Analysis of intersubject variability and menstrual cycle changes,” J. Biomed. Opt 9, 541–552 (2004). 23. M. O’Leary, D. Boas, B. Chance, and A. Yodh, “Experimental images of heterogeneous turbid media by frequency-domain diffusion-photon tomography,” Opt. Lett. 20, 426–428 (1995). 24. J. C. Hebden, F. E. W. Schmidt, M. E. Fry, M. Schweiger, E. M. C. Hillman, D. T. Delpy, and S. R. Arridge, “Simultaneous reconstruction of absorption and scattering images by multichannel measurement of purely temporal data,” Opt. Lett. 24, 534–536 (1999). 25. T. O. McBride, B. W. Pogue, S. Jiang, U. L. Osterberg, and K. D. Paulsen, “Initial studies of in vivo absorbing and scattering heterogeneity in near-infrared tomographic breast imaging,” Opt. Lett. 26, 822–824 (2001). 26. J. Wang, S. D. Jiang, Z. Z. Li, R. M. diFlorio Alexander, R. J. Barth, P. A. Kaufman, B. W. Pogue, and K. D. Paulsen, “In vivo quantitative imaging of normal and cancerous breast tissue using broadband diffuse optical tomography,” Med. Phys. 37, 3715–3724 (2010). 27. B. J. Hoenders, “Existence of invisible nonscattering objects and nonradiating sources,” J. Opt. Soc. Am. A 14, 262–266 (1997). 28. S. Arridge and W. Lionheart, “Nonuniqueness in diffusion-based optical tomography,” Opt. Lett. 23, 882–884 (1998). 29. V. Markel and J. Schotland, “Inverse problem in optical diffusion tomography. II. Role of boundary conditions,” J. Opt. Soc. Am. A 19, 558–566 (2002). 30. L. V. Wang and H. Wu, Biomedical Optics: Principles and Imaging (John Wiley, 2007). 31. J. C. Ye, K. J. Webb, R. P. Millane, and T. J. Downar, “Modified distorted Born iterative method algorithm with an approximate Fr´echet derivative for optical diffusion tomography,” J. Opt. Soc. Am. A 16, 1814–1826 (1999). 32. M. E. Davies and Y. C. Eldar, “Rank awareness in joint sparse recovery,” IEEE Trans. Inf. Theory 58, 1135–1146 (2012). 33. S. F. Cotter, B. D. Rao, K. Engan, and K. Kreutz-Delgado, “Sparse solutions to linear inverse problems with multiple measurement vectors,” IEEE Trans. Signal Process. 53, 2477–2488 (2005). 34. D. Malioutov, M. Cetin, and A. Willsky, “A sparse signal reconstruction perspective for source localization with sensor arrays,” IEEE Trans. Signal Process. 53, 3010–3022 (2005). 35. D. P. Wipf and B. D. Rao, “An empirical Bayesian strategy for solving the simultaneous sparse approximation problem,” IEEE Trans. Signal Process. 55, 3704–3716 (2007). 36. J. M. Kim, O. K. Lee, and J. C. Ye, “Improving noise robustness in subspace-based joint sparse recovery,” IEEE Trans. Signal Process. 60, 5799–5809 (2012).
#197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26590
37. A. Chambolle, “An algorithm for total variation minimization and applications,” J. Math. Imaging Vision 20, 89–97 (2004). 38. M. Afonso, J. Bioucas-Dias, and M. Figueiredo, “An augmented Lagrangian approach to the constrained optimization formulation of imaging inverse problems,” IEEE Trans. Image Process. 20, 681–695 (2011). 39. Y. Pei, H. L. Graber, and R. L. Barbour, “Normalized-constraint algorithm for minimizing inter-parameter crosstalk in DC optical tomography,” Opt. Express 9, 97–109 (2001). 40. Y. Xu, X. Gu, T. Khan, and H. Jiang, “Absorption and scattering images of heterogeneous scattering media can be simultaneously reconstructed by use of DC data,” Appl. Opt. 41, 5427–5437 (2002). 41. Q. Fang and D. A. Boas, “Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units,” Opt. Express 17, 20178–20190 (2009). 42. B. Dogdas, D. Stout, A. F. Chatziioannou, and R. M. Leahy, “Digimouse: A 3D whole body mouse atlas from CT and cryosection data,” Phys. Med. Biol. 52, 577–587 (2007). 43. G. Alexandrakis, F. R. Rannou, and A. F. Chatziioannou, “Tomographic bioluminescence imaging by use of a combined optical-PET (OPET) system: A computer simulation feasibility study,” Phys. Med. Biol. 50, 4225– 4241 (2005). 44. S. A. Prahl, “Optical Properties Spectra,” Oregon Medical Laser Clinic, 2001, http://omlc.ogi.edu/spectra/index.html. 45. D. P. Huttenlocher, G. A. Klanderman, and W. J. Rucklidge, “Comparing images using the Hausdorff distance,” IEEE Trans. Pattern Anal. Mach. Intell. 15, 850–863 (1993).
1.
Introduction
The near-infrared (NIR) optical wavelength range (700 ∼ 1000nm) provides an important opportunity for biological imaging, as tissues are relatively transparent in this regime due to the low absorption rate of the primary absorbers [1, 2]. It is now well-known that NIR light can penetrate tissue up to a depth of several centimeters. Accordingly, the objective of DOT is to exploit this window of opportunity to reconstruct the optical properties of highly scattering biological tissue from scattered and attenuated optical flux measured at the boundary of the medium [3–5]. Extensive studies have been conducted for the theoretical modeling of photon propagation in tissue and to realize accurate reconstructions of the optical parameters [6, 7]. One of the main difficulties of DOT, however, is that the inverse problem is severely ill-posed and non-linear due to the diffusive nature of photon propagation and the non-linear coupling between the unknown optical fluence and parameters [6]. Conventional approaches to address this non-linearity have included linearization methods [8, 9] and iterative methods that update the Green’s function [10–12]. However, the linearization method fails when perturbation in the optical property is beyond the Born approximation limit. Furthermore, the computational overhead of iterative methods is often prohibitive, especially for 3-D imaging because updating the Green’s function requires multiple applications of 3-D partial differential equation (PDE) solvers or Monte Carlo simulation. Recently, a non-iterative reconstruction method was developed by our group based on the theory of joint sparse recovery from compressed sensing [13, 14]. This method exploits the facts that optical parameter variations occur in a relatively small area and their spatial positions are usually invariant during multiple optical illuminations. These properties convert the DOT problem into a joint sparse recovery problem, where we estimate a set of unknown signals with a common sparse support [15, 16]. It has been found that if the background optical parameters are known, the optical parameter perturbations can be reconstructed using a simple two-step procedure in which the unknown support set is initially estimated using a joint sparse recovery algorithm and the optical parameter variations are then calculated using a least-squares fitting method [13, 14]. However, the least-squares step limits the application of the theory to the recovery of only the perturbation in the absorption parameters. While optical absorption parameter perturbation can be observed during brain activity or in the angiogenesis of cancer [1, 4], scattering parameter perturbation is another important property. It has been reported that both absorption and scattering parameters change in breast tumors
#197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26591
[17, 18], and scattering parameter perturbation can provide information about breast composition and physiology changes that may be useful for diagnosing breast malignancy [19]. Moreover, the scattering amplitude and power correlate with the degree of radiographic density, which is correlated to the risk of developing breast cancer [20–22]. Accordingly, simultaneous reconstruction of both absorption and scattering coefficients change has been extensively studied using phantoms [23, 24] and in vivo experiments [25, 26], and theoretical analyses of the existence and uniqueness of solutions have also been reported [9, 27–29]. In this paper, we extended the results in our earlier work [13, 14] and propose a non-iterative method for the simultaneous reconstruction of both absorption and scattering parameter perturbation when the background optical properties are known. One of the key ideas for such an extension is that the least-squares step in the preceding studies [13, 14] can be easily replaced by a linear integral equation of the first kind, thanks to the self-recursive structure in the FoldyLax equation. In particular, as long as we can identify the support set for the absorption or scattering parameter variations accurately, unknown fluence can be estimated accurately, which eliminates the nonlinear terms in the forward integral equation. Moreover, the estimation problem for the support of absorption and scattering parameter variations can be addressed as a joint sparse recovery approach, which is similar to the earlier studies [13, 14]. Based on the new formulation, we demonstrate that the reconstruction quality is improved and, furthermore, that the cross-talk between two optical parameters can be reduced. The paper is organized as follows. After a brief review of conventional methods is given, we describe a theory of simultaneous absorption and scattering parameter estimation using joint sparsity in Section 2. The implementation of the algorithm and the simulation environment are described in detail in Section 3. Section 4 provides numerical results, which is followed by a discussion of the implications of the results in Section 5. Finally, the conclusion is provided in Section 6. 2.
Theory
In this section, a continuous-wave (CW) formulation is described to explain a non-iterative simultaneous reconstruction of perturbation of both absorption and scattering coefficients. An extension of the theory for frequency domain DOT formulation is also straightforward. 2.1.
Diffusion equation
Photon migration within biological tissues is usually modeled by the transport equation. However, when light scattering prevails over absorption, the propagation of light can be modeled by a diffusion equation. Let Ω be a domain filled with some turbid medium with ∂ Ω as its boundary. In a highly scattering medium with low absorption, the optical fluence u(r) at position r can be modeled by the following continuous-wave diffusion equation (DE) [29]: ∇ · D(r)∇u(r) − μa (r)u(r) = −S(r), r ∈ Ω (1) u(r) + nˆ · ∇u(r) = 0, r ∈ ∂Ω Here, μa (r) and D(r) = 1/3(μs (r) + μa (r)) are the absorption and diffusion coefficients, where μs (r) denotes the reduced scattering coefficient, and S(r) is the source intensity profile. Additionally, is a parameter related to the diffusion coefficient, dimension and reflection on the boundary [30], and nˆ denotes a vector normal to the measurement surface. In a DOT problem, the unknown optical parameters are modeled as D(r) = D0 (r) + δ D(r), μa (r) = μa0 (r) + δ μa (r), and μs (r) = μs0 (r) + δ μs (r), where μa0 (r) and μs0 (r) are the known background absorption and reduced scattering parameters, respectively. D0 (r) = 1/3(μs0 (r) + μa0 (r)) and δ D(r), δ μa (r), and δ μs (r) refer to the unknown perturbations, which we need to
#197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26592
estimate. Note that the known background optical parameter does not need to be constant, and our theory holds for any spatially varying known background. The diffusion equation can be then transformed as (∇ · D0 (r)∇ − k(r)2 )u(r) = [δ μa (r) − ∇ · δ D(r)∇]u(r) − S(r), r ∈ Ω (2) u(r) + nˆ · ∇u(r) = 0, r ∈ ∂ Ω, where the diffusion wave number is given by k(r) =
δ μs (r) μa0 (r) and δ D(r) = − δμμa a (r)+ (r)+μ (r) D(r). 0
Then, Eq. (2) can be equivalently represented using the integral equation: u(r) = u0 (r) −
Ω
s0
G0 (r, r )[δ μa (r ) − ∇ · δ D(r )∇]u(r )dr , r ∈ Ω,
(3)
where r denotes a position in volume Ω and the background Green’s function G0 (r, r ) satisfies (∇ · D0 (r)∇ − k(r)2 )G0 (r, r ) = −δ (r − r ), r ∈ Ω (4) r ∈ ∂ Ω, G0 (r, r ) + nˆ · ∇G0 (r, r ) = 0, and the incidence fluence u0 (r) is given by u0 (r) =
Ω
G0 (r, r )S(r )dr .
(5)
Now, assuming that δ D(r) is zero in the neighborhood around the boundary, we have 0
= = =
∂Ω
Ω Ω
nˆ · [δ D(r )G0 (r, r )∇u(r )]dr
∇ · [δ D(r )G0 (r, r )∇u(r )]dr
δ D(r )∇G0 (r, r ) · ∇u(r )dr +
Ω
G0 (r, r )∇ · δ D(r )∇u(r )dr ,
(6)
where we use the divergence theorem to derive the second equality. Therefore, we have u(r) = u0 (r) −
Ω
G0 (r, r )u(r )δ μa (r )dr −
Ω
∇G0 (r, r ) · ∇u(r )δ D(r )dr , r ∈ ∂ Ω.
(7)
Note that the resulting forward integral equation that maps the unknown optical parameter perturbation to the boundary measurement is nonlinear, as u(r ) depends on the unknown optical parameter perturbation (δ μa , δ D). 2.2.
Proposed method
In a conventional iterative approach such as the distorted Born iterative method (DBIM), Eq. (7) is linearized as u(r; x) ≈ u(r; x(k) ) − −
Ω Ω
G0 (r, r ; x(k) )u(r ; x(k) )δ μa (r )dr ∇G0 (r, r ; x(k) ) · ∇u(r ; x(k) )δ D(r )dr ,
(8)
(k)
where x(k) = [μa , D(k) ]T denotes the k-th estimate of the background optical parameter, and [δ μa , δ D]T denotes the unknown perturbation estimated from the k-th estimate of the background optical parameters. Note that the fluence u(r ; x) is linearized with respect to the back#197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26593
ground optical parameter x(k) . After the linearization, [δ μa , δ D]T is estimated and we perform (k+1) (k+1) (k) , D(k+1) ]T , where μa = μa + δ μa , and another linearization of u(r ; x) at x(k+1) = [μa (k+1) (k) = D + δ D. This procedure is repeated until convergence. It has been shown that the D linearized integral equation corresponds to the Fr´echet derivative and the DBIM is equivalent to the Levenberg-Marquardt algorithm [31]. However, the main problem with this type of iterative linearization is that we need to calculate the Green’s function multiple times, which requires multiple applications of computationally extensive forward solvers. Hence, in some conventional algorithms, only one linearization step is performed (i.e., Born approximation), which often degrades the performance. To overcome these drawbacks, we developed a novel non-iterative method for the simultaneous reconstruction of both absorption and scattering coefficients when the background optical properties are known. The main idea comes from the observation that the absorption or scattering properties change locally against the background (e.g., due to localized angiogenesis for cancer or brain activation) and that the reconstruction problem can be converted into a joint sparse recovery. Figure 1 illustrates such a concept. The target area Ωt represents areas where
Fig. 1. Sparsely distributed perturbations in optical properties for simultaneous DOT reconstruction.
the absorption or scattering properties change. We assumed that the area Ωt is sparse compared to the total field of view (FOV) Ω. This is usually true in practice because, for example, tumors are only sparsely distributed. Furthermore, we assumed that there are N distinct illumination patterns and let u(r; l) be the optical fluence from the l-th illumination pattern. Then, Eq. (3) can be rewritten as u(r; l)
= u0 (r; l) − = u0 (r; l) −
Ω
Ωt
G0 (r, r )X(r ; l)dr G0 (r, r )X(r ; l)dr , l = 1, 2, · · · , N,
(9)
where X(r; l) = [δ μa (r) − ∇ · δ D(r)∇]u(r; l) is the induced current from the l-th illumination. In Eq. (9), the second equality comes from the following:
δ μa (r) = 0, δ D(r) = 0, ∇δ D(r) = 0, r ∈ / Ωt .
(10)
Once the support Ωt is estimated, X(r; l) can easily be obtained by means of linear inversion from Eq. (9). Then, the unknown optical fluence u(r; l) can be calculated by restricting r ∈ Ωt . More specifically, we have the following Foldy-Lax equation for fluence estimation, as derived
#197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26594
from Eq. (9) by restricting r ∈ Ωt : u(r; ˆ l) = u0 (r; l) −
Ωt
ˆ ; l)dr , r ∈ Ωt . G0 (r, r )X(r
(11)
ˆ l) denotes the estimated induced current for the l-th illumination. Then, using u(r; Here, X(r; ˆ l), δ μa and δ D can be reconstructed from the following linear integral equation: u(r; l) = u0 (r; l) −
Ωt
G0 (r, r )u(r ˆ ; l)δ μa (r )dr −
Ωt
∇G0 (r, r ) · ∇u(r ˆ ; l)δ D(r )dr ,
(12)
l = 1, 2, · · · , N. ˆ ; l) Note that the nonlinear coupling due to u(r ; l) in Eq. (7) becomes now linear because u(r and ∇u(r ˆ ; l) are calculated from Eq. (11). Note that the major difference from the conventional methods is that as the nonlinearly coupled optical fluence u(r; l) for r ∈ Ωt is now estimated, the iterative procedure for updating Green’s function becomes unnecessary. Therefore, one of the most important steps for the use of Eq. (12) is the joint support estimation of Ωt from multiple illumination patterns to obtain an accurate estimation of the unknown fluence u(r; l). Another important issue is to estimate the optical parameters using Eq. (12). In our previous work [13, 14], we assumed that the perturbation occurs in the absorption parameter alone, which is calculated by least-squares fitting at each position r ∈ Ωt from the estimated ˆ l) = δ μa (r)u(r; ˆ l) and u(r; ˆ l) values for l = 1, 2, · · · , N. However, the presence of the perX(r; turbation in the scattering parameter inhibited the use of this step. In the following section, we describe how to address these issues. 2.3.
Joint support recovery
We assumed that the induced current is described by piecewise constant or spline approximation as: k
X(r; l) = ∑ X(r(i) ; l)b(r, r(i) ), l = 1, 2, · · · , N,
(13)
i=1
where b(r, r(i) ) denotes the basis function. The 3-D locations {r(i) }ki=1 ∈ Ωt of the k-targets are selected from n-possible locations {ri }ni=1 ∈ Ω, where k n. Then, we have the following relationship after substituting Eq. (13) into the induced current term in Eq. (9), u(r; l) = u0 (r; l) −
Ωt
G0 (r, r )X(r ; l)dr , l = 1, 2, · · · , N,
k
= u0 (r; l) − ∑ G˜ 0 (r, r(i) )X(r(i) ; l),
(14)
i=1
where G˜ 0 (r, r(i) ) =
Ωt
G0 (r, r )b(r , r(i) )dr .
(15)
By collecting the measurement of Eq. (14) at detector positions {rdi }m i=1 for all N illuminations, the imaging problem can be represented as the following matrix equation: Y = AX.
(16)
#197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26595
In this equation, Y ∈ Rm×N is scattered fluence, and A ∈ Rm×n and X ∈ Rn×N are the sensing matrix and the induced current matrix, respectively. Let the component of the i-th row and j-th column of matrix A be defined as Ai j . Then, Yi j = u0 (rdi ; j) − u(rdi ; j), Ai j = G˜ 0 (rdi , r j ). Note that the rows of X are only non-zero when ri ∈ {r(i) }ki=1 . Now, we define the row-diversity measure ||X||0 that counts the number of rows in X ∈ Rn×N with nonzero elements. Then, the simultaneous reconstruction of the DOT problem from multiple illumination cases, where the support of the optical parameter variation is sparse, can be stated as the following joint sparse recovery problem [14, 15]: min ||X||0 , subject to ||Y − AX||F ≤ ε .
(17)
X
Here, ||X||F is the Frobenius norm of matrix X. Note that when δ D(r(i) ) = 0 for i = 1, 2, · · · , k, Eq. (17) becomes precisely the same joint sparse recovery problem in our previous work [13, 14]. Even if δ D(r(i) ) = 0, in terms of implementation, the support estimation step is identical to that of the previous work, since δ D is recovered later from the unknown X. One of the main advantages of this joint sparse recovery formulation is that it enables a further reduction in the number of required measurements [15, 16, 32], as the number of measurements required per sensor must account for the minimal features unique to that sensor. There exist a variety of joint sparse recovery algorithms, including ones from our group [16, 33–36]. The central theme of these algorithms is to exploit the diversity of multiple measurements to provide reconstruction from fewer samples or from a low SNR. Among the various algorithms, in this paper we used the M-SBL algorithm [35], given that we did not need to know the sparsity level and because the algorithm is known to be robust to noise. 2.4.
Optical parameter estimation ˆ
ˆ t and the calculated optical fluence u(r ˆ t, ˆ i ) at ri ∈ Ω Based on the estimated supports {r(i) }ki=1 ∈ Ω the unknown optical parameters can be reconstructed by solving the following discretized linear inverse problem derived from Eq. (12): ⎤ ⎡ δ μa,1 ⎤ ⎢ .. ⎥ ⎡ (1) ⎤ ⎡ (1) (1) ⎢ . ⎥ y Aμa AD ⎥ ⎢ ⎥ ⎢ δμ ˆ ⎥ ⎢ .. ⎥ ⎢ . . a, k ⎥ ⎥ = Ax. ⎢ ⎢ .. ⎦ ⎢ y = ⎣ . ⎦ = ⎣ .. (18) δ D1 ⎥ ⎥ ⎢ (N) (N) (N) ⎢ . ⎥ y Aμa AD ⎣ .. ⎦ δ Dkˆ (l)
(l)
(l)
ˆ ( j) ; l)δ v, AD,i j = ∇G0 (rdi , r( j) ) · Here, yi = u0 (rdi ; l) − u(rdi ; l), Aμa,i j = G0 (rdi , r( j) )u(r ∇u(r ˆ ( j) ; l)δ v, and δ v is the volume of the discretized voxel. In addition, δ μa,i = δ μa (r(i) ) and δ Di = δ D(r(i) ) denote the unknown optical parameter perturbations. Due to the restriction of ˆ t , the number of unknowns in the discretized voxel domain is significantly reduced. the ROI to Ω However, the ill-posedness of the problem is still not completely overcome by the restriction since the ill-posedness of the DOT problem originates from the underlying continuous operator. In fact, it is well-known that the linear mapping from DOT is compact mapping, whose singular values are clustered around zero, which makes the inverse mapping very ill-posed. Therefore, to solve Eq. (18), we convert it into the following constrained optimization problem with L1 and the total variation regularization under the assumption of the sparsity and smoothness of
#197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26596
both absorption and scattering perturbations: arg min ||δ μa ||1 + ||δ D||1 + λ ||δ μa ||TV + λ ||δ D||TV x
(19)
||Ax − y||2 ≤ ε , amin ≤ δ μai ≤ amax , Dmin ≤ δ Di ≤ Dmax .
subject to
ˆ
ˆ
In this equation, x = [δ μaT , δ DT ]T ∈ R2k×1 , A = [Aμa , AD ] ∈ RmN×2k , and ||x||1 and ||x||TV are L1 penalty and the total variation of x [37], respectively. Additionally, [amin , amax ] and [Dmin , Dmax ] denote the ranges of the absorption and diffusion parameter variations, respectively. 3. 3.1.
Method Implementation
In this section, we describe the implementation procedure for the proposed method and conventional algorithm in detail. In Eq. (19), λ is set to 1 for the proposed method, whereas it is reduced to 0.03 ∼ 0.05 for conventional methods since we found that the more weight in TV provides inaccurate absorption parameter reconstruction in the conventional methods. Therefore, we wanted to adjust the parameter to make the comparison fair. To solve the aforementioned constrained optimization problem, we exploited a constrained split augmented Lagrangian shrinkage algorithm (C-SALSA) [38] after considering that it is quite effective in dealing with various regularization schemes and constraints. In particular, Chambolle’s algorithm was employed to deal with TV regularization and a soft-threshold function with L1 penalty as in the C-SALSA algorithm [37, 38]. We normalized the sensing matrix in Eqs. (17) and (19) to have the unit norm for each column. For a fair comparison, C-SALSA was also used for conventional methods. More specifically, the same constrained optimization problem in Eq. (19) can also be applied to DBIM or a linearized approach with (l) the sensing matrix A = [Aμa , AD ] ∈ RmN×2n , where Aμa,i j = G0 (rdi , r j ; x(k) )u(r j ; l, x(k) )δ v and (l)
AD,i j = ∇G0 (rdi , r j ; x(k) ) · ∇u(r j ; l, x(k) )δ v for k-th iteration. When k = 1, we use the known background optical parameters as x(1) = [μa0 , D0 ]T . However, note that the size of the sensing matrix for the proposed method is much smaller than that of the conventional methods due to the restriction of the problem within the estimated joint support. The pseudocode implementation of C-SALSA is described in Algorithm 1. In Algorithm 1, Ψτ gi is the Moreau proximal mapping of τ gi [38], where g1 (x) = ||x||1 , g2 (x) = ||x||TV , g3 = ιE(ε ,y) , and g4 = ιB(a,D) . ιE(ε ,y) and ιB(a,D) are the indicator functions of an ε -radius Euclidean ball centered at y and a Box-constraint which depends on the range of x, respectively. As a stopping criterion, we used the relative change of the cost function in Eq. (19) and performed the algorithms until |(Ck −Ck−1 ) /Ck | < 5 × 10−4 is satisfied, where Ck is the cost function at k-th iteration. All simulation parameters are summarized in Table 1. (1) In Table 1, s¯1 (1) is the mean value of the s1 , which is defined in Algorithm 1, and εi and yi are the values of ε and measurements y for the i-th iteration of the DBIM method. The scale factor for ε was forced to increase during iterations in the case of DBIM, since we found that a fixed value often makes the iteration not converge. These parameters were found manually for each method by choosing the optimal parameter among c = {8, 4, 2, 1, 1/2, 1/4, 1/8} for τ = cs¯1 (1) and c = {0.01, 0.03, 0.05, 0.075, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6} for ε1 = c||y1 ||2 , and we did an additional search for ε2 and ε3 from 0.7 ∼ 1.0 to make the DBIM converge. Since the sensing matrices for the proposed method and conventional methods are different, the selected optimal parameters are distinct, as shown in Table 1, for each method.
#197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26597
Algorithm 1 Pseudocode implementation of C-SALSA for DOT problem in Eq. (19). (i)
(i)
Set k = 0, H (1,4) = I, H (2) = λ I, H (3) = A, and choose τ > 0, v0 , d0 , for i = 1, · · · , 4. repeat for i = 1, . . . , 4 do (i) (i) (i) 4: ζk = vk + dk 5: end for −1
T
T 4 4 ( j) ( j) ( j) 6: uk+1 = ∑ H H ∑ H ( j) ζk 1: 2: 3:
j=1
7: 8:
j=1
for i = 1, . . . , 4 do (i) (i) (i) (i) vk+1 = Ψτ gi sk , where sk = H (i) uk+1 − dk (i)
(i)
(i)
dk+1 = dk − H (i) uk+1 + vk+1 end for k = k+1 12: until some stopping criterion is satisfied 9: 10: 11:
Table 1. Parameters used in C-SALSA algorithm.
Proposed Method
Simple Target Simulation
Mouse Phantom Simulation
λ =1
λ =1
τ = 0.25s¯1
(1)
ε = 0.01||y||2
Conventional Methods
3.2.
τ = 2s¯1 (1) ε = 0.01||y||2
λ = 0.05
λ = 0.03
τ = 4s¯1 (1)
τ = 4s¯1 (1)
ε1 = 0.01||y1 ||2
ε1 = 0.03||y1 ||2
ε2 = 0.9||y2 ||2
ε2 = 0.7||y2 ||2
ε3 = ||y3 ||2
ε3 = 0.9||y3 ||2
Simulation geometry I : simple targets
First, we considered two simple targets immersed in a rectangular homogeneous medium, as in Fig. 2. Sixteen line sources with a 2mm space in total were used, and a 16 × 16 array detector with a pitch of 1.33mm measured the optical flux at the bottom of the medium. The medium was discretized with 2 × 2 × 2mm3 voxels and the region of interest (ROI) for reconstruction was 22 × 26 × 22mm3 in size. The optical parameters for the background were set to μa = 0.015mm−1 , and the μs were set to be equal to 1mm−1 . We conducted four simulation studies of this geometry, as described in Table 2, as CASE A to CASE D. In CASE A and CASE B, each target had optical perturbation in both absorption and scattering coefficients, at low and high contrast, respectively. On the other hand, in CASE C and CASE D, target 1 had only absorption changes, whereas target 2 had scattering changes. These cases were designed to evaluate the cross-talk artifacts [39, 40]. For an accurate simulation, we solved the diffusion equation using a GPU-accelerated Monte-Carlo simulation [41]. Noise with a SNR of 40dB was added to the scattered fluence. The results were obtained by averaging 10 runs with independent noise realization for each case.
#197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26598
Fig. 2. Simulation geometry for simple targets.
Table 2. Optical properties for various cases. CASE B
CASE C
CASE D
(low contrast)
(high contrast)
(low contrast)
(high contrast)
δ μa
δ μs
Target 1
0.0015
Target 2
0.0015
Optical parameters in
3.3.
mm−1
CASE A
δ μa
δ μs
δ μs
δ μa
δ μa
0.05
0.03
1.0
0.05
0.03
1.0
δ μs
0.0015
0.0
0.03
0.0
0.0
0.05
0.0
1.0
Simulation geometry II : mouse phantom
Next, we considered more complicated simulation studies using the Digimouse computational mouse phantom [42]. The goal of this study was to localize tumors within a mouse and to recover any optical parameter changes in the tumors. The simulation geometry is described in Fig. 3(a). The mouse phantom was placed in a rectangular shape of matching fluid that has the same optical parameters as mouse skin. Sixteen line sources in total with a space of 2mm were again used, and a 20×20 array detector with a pitch of 1.5mm served to measure the optical flux at the bottom of the medium. The medium was discretized with 2 × 2 × 2mm3 voxels. The ROI for reconstruction is the area between the source and detector with a size of 38 × 30 × 22mm3 , where three tumors were placed, as illustrated in Fig. 3(b). Various optical parameters of the mouse organs with a wavelength of 800nm [43, 44], summarized in Table 3, were used. For simplicity, the optical parameters corresponding to bones, fat, muscle, and other parts were assumed to be the same as skin. The tumor in the liver is denoted as target 1, which has both absorption and scattering coefficients changes, whereas targets 2 and 3, corresponding to tumors in kidneys, have changes only in their absorption and scattering coefficients, respectively, as described in Table 4. Except for the tumors, all spatially varying optical parameters were assumed to be known. We solved the diffusion equation using the GPU-accelerated Monte-Carlo simulation [41], and noise with a SNR of 40dB was added to the scattered fluence. The results were also obtained by averaging 10 runs with independent noise realizations. 4.
Results
We compared the reconstruction results of the proposed method to those of the conventional iterative (DBIM) and linearization methods. The performance was evaluated quantitatively using
#197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26599
Fig. 3. Simulation geometry of the mouse phantom.
Table 3. Optical properties of the mouse [43, 44].
Tissue type
Skin
Heart
μa in mm−1
0.0259
0.0218
μs in mm−1
1.8258
0.748
Lung
Liver
Stomach
Kidney
0.0636
0.126
0.0058
0.0249
1.9789
0.5629
1.2098
1.7237
Table 4. Optical properties of the targets.
optical parameters in mm−1
δ μa
δ μs
Target 1
0.0504
0.8444
Target 2
0.0423
0.0
Target 3
0.0
0.8619
the following relative mean square error: MSE(x) =
||xtrue − xrecon ||2 . ||xtrue ||2
(20)
Although the final equation of Eq. (19) depends on δ D, our practical interest is δ μs . Therefore, δ μs (r) we calculated δ μs using the relationship of δ D(r) = − δμμa a (r)+ (r)+μ (r) D(r) based on the estimated 0
s0
values of δ μa and δ D, as shown in the following: δ D(r) + δ μa (r) . δ μs (r) = − 3D0 (r)[δ D(r) + D0 (r)]
(21)
We also used the Hausdorff distance [45] to measure the morphological difference. The number of iterations for the DBIM method was set to 3.
#197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26600
4.1.
Simple target
The reconstruction results of the simple target simulation are summarized in Table 5. The MSE values in Table 5 are the average values of the MSE values from multiple noise simulation results, and the values in the parenthesis denote the corresponding standard deviation. The relative error for the proposed method is competitive or lower than those of DBIM or the linearized approach. For the case of the linearized approach, the mean MSE values and their standard deviations were the smallest when the perturbation in optical parameters are small (CASE A). However, the algorithm broke down for large optical perturbation (CASE B) or in the presence of the simultaneous perturbation of absorption and scattering (CASE C and CASE D). The results of DBIM were slightly improved compared to those of the linearized approach for absorption coefficient cases. On the other hand, the results of the proposed method for the reconstruction of the absorption coefficients were robust for all cases, and the MSE values for scattering coefficients showed monotonic increasing behaviors from CASE A to CASE D. These phenomena can be easily confirmed from Figs. 4 and 5, which show the cross-sectional image of the reconstructed images for δ μa and δ μs for each case (The horizontal and coronal sections centered on each target are indicated in the inset image in Figs. 4 and 5). Moreover, in the reconstruction images, our method reduced the cross-talk artifacts that were observed in the conventional methods, especially in the reconstructed scattering coefficients. The average reconstruction time of various methods for the simple target simulation are summarized in Table 6 (using a PC with CPU : core i7 sandy bridge and GPU : GTX 560 Geforce series). The support estimation step using the M-SBL algorithm was also added in the calculation of the run time of the proposed method. Even though the proposed method requires this additional step, the resulting computational time for the proposed method is even smaller compared to the linearized approach. In Table 6, we also describe the computation time for each part of the algorithm. Note that the total computation time includes the time required to construct the A matrix, which was about 1.2 seconds. Table 5. The MSE values of the reconstruction using various methods.
proposed method
linearized approach
DBIM
CASE A
MSE(δ μa )
0.149 (± 0.009)
0.360 (± 0.004)
0.359 (± 0.004)
(low contrast)
MSE(δ μs )
0.165 (± 0.009)
0.176 (± 0.006)
0.176 (± 0.006)
CASE B
MSE(δ μa )
0.134 (± 0.007)
0.547 (± 0.069)
0.525 (± 0.076)
(high contrast)
MSE(δ μs )
0.247 (± 0.008)
0.319 (± 0.081)
0.352 (± 0.134)
CASE C
MSE(δ μa )
0.143 (± 0.050)
0.630 (± 0.045)
0.628 (± 0.046)
(low contrast)
MSE(δ μs )
0.328 (± 0.039)
0.297 (± 0.068)
0.301 (± 0.080)
CASE D
MSE(δ μa )
0.123 (± 0.026)
0.715 (± 0.048)
0.696 (± 0.046)
(high contrast)
MSE(δ μs )
0.488 (± 0.028)
0.456 (± 0.052)
0.523 (± 0.045)
4.2.
Mouse phantom
The reconstruction results of the mouse phantom using various methods are summarized in Table 7. The MSE values for the proposed method are lower than those for DBIM or the linearized approach for both the absorption and scattering coefficients. Table 8 summarizes the Hausdorff distance results [45] between the true supports and the estimated supports of each tumor in the liver and kidneys. Both of these quantitative measures show the advantages of the proposed #197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26601
Fig. 4. The cross-sections of the reconstructed images using various methods for simple target simulations (CASE A and CASE B). Table 6. Average run time of the reconstruction using various methods.
proposed method Run time (sec)
linearized approach
M-SBL
C-SALSA
Total
C-SALSA
Total
0.34
0.22
1.92
5.71
7.08
DBIM 1537.5
method. Figure 6 illustrates the location of the tumor volume thresholded with the same number of non-zero supports for each instance of δ μa and δ μs (29 and 30), respectively. The conventional methods exhibited cross-talk between μa and μs , whereas the proposed method did not have such artifacts. Table 7. MSE for the reconstruction results of the mouse phantom using various methods.
proposed method
MSE(δ μa )
MSE(δ μs )
0.326 (± 0.049)
0.624 (± 0.121)
linearized approach
0.491 (± 0.004)
0.770 (± 0.008)
DBIM
0.477 (± 0.005)
0.775 (± 0.008)
#197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26602
Fig. 5. The cross-sections of the reconstructed images using various methods for simple target simulations (CASE C and CASE D).
Fig. 6. The 3-D simultaneous reconstruction results of tumors with optical parameter variations δ μa and δ μs embedded in the mouse phantom.
#197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26603
Table 8. Hausdorff distance in mm for various reconstruction methods. proposed method Hausdorff distance in mm Liver Kidney (left) Kidney (right)
5.
Linearized approach
DBIM
δ μa
δ μs
δ μa
δ μs
δ μa
δ μs
2.0 (± 0.0) 2.0 (± 0.0) -
2.0 (± 0.0) 2.0 (± 0.0)
4.0 (± 0.0) 2.83 (± 0.0) -
2.08 (± 0.263) 2.0 (± 0.0)
3.0 (± 1.054) 2.83 (± 0.0) -
2.08 (± 0.263) 2.0 (± 0.0)
Discussion
In the existing non-iterative method for absorption imaging [14], the second step to calculate the absorption coefficient was performed in a voxel-by-voxel manner using the least-squares fitting approach. However, this step was not valid in the simultaneous reconstruction problem due to the differential operator applied to the diffusion coefficient. In this paper, we demonstrated that the unknown optical fluence can be estimated, regardless of the differential equation form, thanks to the recursive nature of the Foldy-Lax equation. Traditionally, the cross-talk issue was considered unavoidable during the simultaneous reconstruction of DOT problems with CW imaging. Using experimental data, the reconstruction of both absorption and scattering coefficient perturbations was shown to be impossible in CW cases [28]. To overcome this limitation, Pei et al. [39] proposed the normalized-constraint algorithm, and Xu et al. [40] used both normalization and regularized techniques. However, the proposed method is totally different from these works since we utilized the joint sparsity. While the conventional methods require only the optical parameter estimation step, our method needs the additional task of finding the support set on which the optical parameter estimation is conducted. Therefore, the joint sparse recovery step is important for the accurate localization of the target and subsequent calculation of optical parameters. In our algorithm, the main two computational building parts, the M-SBL and C-SALSA step, has complexity that follows the dimension of the sensing matrix. Recall that in the proposed method, the dimension of the matrix for C-SALSA is reduced due to the restriction of the unknown within the estimated support. Moreover, the complexity of M-SBL depends on the number of detectors. Unlike the conventional method, whose complexity depends on the number of source-detector pairs, the significant reduction of the sensing matrix dimension makes the M-SBL and C-SALSA computationally very efficient. Hence, the combination of M-SBL and C-SALSA results in the reduction of the computational burden compared to the linearized approaches. 6.
Conclusion
In this paper, we proposed a non-iterative method for the simultaneous reconstruction of both absorption and scattering coefficients in diffuse optical tomography. By exploiting multiple illumination patterns and the sparseness of the target area, a DOT inverse problem was converted to a joint sparse recovery problem. Moreover, the unknown optical parameters were calculated without an approximation or an iterative update of the Green’s function, thanks to the recursive nature of the Foldy-Lax equation. The proposed method was validated with various simulation studies, showing an improvement over conventional methods. Furthermore, our algorithm reduced the cross-talk artifacts during simultaneous reconstruction from CW data. Acknowledgments This research was supported by a grant from the Korea Science and Engineering Foundation (KOSEF) funded by the Korean government (MEST) (No.2012-0000173).
#197011 - $15.00 USD Received 3 Sep 2013; revised 9 Oct 2013; accepted 17 Oct 2013; published 28 Oct 2013 (C) 2013 OSA 4 November 2013 | Vol. 21, No. 22 | DOI:10.1364/OE.21.026589 | OPTICS EXPRESS 26604
