Joint sparsity-driven non-iterative simultaneous reconstruction of absorption and scattering in diffuse optical tomography

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échet 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/in


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][4][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][11][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 [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][21][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][28][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 Foldy-Lax 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.

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.

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]: 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) = D 0 (r) + δ D(r), μ a (r) = μ a 0 (r) + δ μ a (r), and μ s (r) = μ s 0 (r) + δ μ s (r), where μ a 0 (r) and μ s 0 (r) are the known background absorption and reduced scattering parameters, respectively.D 0 (r) = 1/3(μ s 0 (r) + μ a 0 (r)) and δ D(r), δ μ a (r), and δ μ s (r) refer to the unknown perturbations, which we need to 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 where the diffusion wave number is given by k(r) = μ a 0 (r) and δ D(r) = − δ μ a (r)+δ μ s (r) μ a 0 (r)+μ s 0 (r) D(r).Then, Eq. ( 2) can be equivalently represented using the integral equation: where r denotes a position in volume Ω and the background Green's function G 0 (r, r ) satisfies and the incidence fluence u 0 (r) is given by Now, assuming that δ D(r) is zero in the neighborhood around the boundary, we have where we use the divergence theorem to derive the second equality.Therefore, we have 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).

Proposed method
In a conventional iterative approach such as the distorted Born iterative method (DBIM), Eq. ( 7) is linearized as where 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- ground optical parameter x (k) .After the linearization, [δ μ a , δ D] T is estimated and we perform another linearization of u(r ; x) at a + δ μ a , and D (k+1) = D (k) + δ D. This procedure is repeated until convergence.It has been shown that the linearized integral equation corresponds to the Fréchet 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 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 where In Eq. ( 9), the second equality comes from the following: 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 from Eq. ( 9) by restricting r ∈ Ω t : Here, X(r; l) denotes the estimated induced current for the l-th illumination.Then, using û(r; l), δ μ a and δ D can be reconstructed from the following linear integral equation: Note that the nonlinear coupling due to u(r ; l) in Eq. ( 7) becomes now linear because û(r ; l) and ∇ û(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 X(r; l) = δ μ a (r) û(r; l) and û(r; l) values for l = 1, 2, • • • , N. However, the presence of the perturbation in the scattering parameter inhibited the use of this step.In the following section, we describe how to address these issues.

Joint support recovery
We assumed that the induced current is described by piecewise constant or spline approximation as: where b(r, r (i) ) denotes the basis function.The 3-D locations {r (i) } k i=1 ∈ Ω t of the k-targets are selected from n-possible locations {r i } n i=1 ∈ Ω, where k n.Then, we have the following relationship after substituting Eq. ( 13) into the induced current term in Eq. ( 9), where By collecting the measurement of Eq. ( 14) at detector positions {r d i } m i=1 for all N illuminations, the imaging problem can be represented as the following matrix equation: In this equation, Y ∈ R m×N is scattered fluence, and A ∈ R m×n and X ∈ R n×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 A i j .Then, Y i j = u 0 (r d i ; j) − u(r d i ; j), A i j = G0 (r d i , r j ).Note that the rows of X are only non-zero when r i ∈ {r (i) } k i=1 .Now, we define the row-diversity measure ||X|| 0 that counts the number of rows in X ∈ R n×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]: Here, ||X|| F is the Frobenius norm of matrix X.Note that when δ D(r 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][34][35][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.

Optical parameter estimation
Based on the estimated supports {r (i) } k i=1 ∈ Ωt and the calculated optical fluence û(r i ) at r i ∈ Ωt , the unknown optical parameters can be reconstructed by solving the following discretized linear inverse problem derived from Eq. ( 12): y (1)  . . .
Here, y and δ v is the volume of the discretized voxel.In addition, δ μ a,i = δ μ a (r (i) ) and δ D i = δ D(r (i) ) denote the unknown optical parameter perturbations.Due to the restriction of the ROI to Ωt , the number of unknowns in the discretized voxel domain is significantly reduced.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 L 1 and the total variation regularization under the assumption of the sparsity and smoothness of both absorption and scattering perturbations: (19) subject to ||Ax − y|| 2 ≤ ε, and ||x|| 1 and ||x|| TV are L 1 penalty and the total variation of x [37], respectively.Additionally, [a min , a max ] and [D min , D max ] denote the ranges of the absorption and diffusion parameter variations, respectively.

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 L 1 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 the sensing matrix A = [A μ a , A D ] ∈ R mN×2n , where A (l) μ a,i j = G 0 (r d i , r j ; x (k) )u(r j ; l, x (k) )δ v and A (l) D,i j = ∇G 0 (r d i , 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) = [μ a 0 , D 0 ] 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, Ψ τg i is the Moreau proximal mapping of τg i [38], where g 1 (x) = ||x|| 1 , g 2 (x) = ||x|| TV , g 3 = ι E(ε,y) , and g 4 = ι 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 is satisfied, where C k is the cost function at k-th iteration.All simulation parameters are summarized in Table 1.In Table 1, s1 (1) is the mean value of the s 1 , which is defined in Algorithm 1, and ε i and y i 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 τ = c s1 (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||y 1 || 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.19).

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 × 2mm 3 voxels and the region of interest (ROI) for reconstruction was 22 × 26 × 22mm 3 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.

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 × 2mm 3 voxels.The ROI for reconstruction is the area between the source and detector with a size of 38 × 30 × 22mm 3 , 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.

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 Although the final equation of Eq. ( 19) depends on δ D, our practical interest is δ μ s .Therefore, we calculated δ μ s using the relationship of δ D(r) = − δ μ a (r)+δ μ s (r) μ a 0 (r)+μ s 0 (r) D(r) based on the estimated values of δ μ a and δ D, as shown in the following: We also used the Hausdorff distance [45] to measure the morphological difference.The number of iterations for the DBIM method was set to 3.

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.

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

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.

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.

Fig. 4 .
Fig. 4. The cross-sections of the reconstructed images using various methods for simple target simulations (CASE A and CASE B).

Fig. 5 .
Fig. 5.The cross-sections of the reconstructed images using various methods for simple target simulations (CASE C and CASE D).

Fig. 6 .
Fig. 6.The 3-D simultaneous reconstruction results of tumors with optical parameter variations δ μ a and δ μ s embedded in the mouse phantom.

Table 2 .
Optical properties for various cases.

Table 4 .
Optical properties of the targets.

Table 5 .
The MSE values of the reconstruction using various methods.

Table 6 .
Average run time of the reconstruction using various methods.

Table 7 .
MSE for the reconstruction results of the mouse phantom using various methods.

Table 8 .
Hausdorff distance in mm for various reconstruction methods.