Hierarchical Bayesian estimation improves depth accuracy and spatial resolution of diffuse optical tomography

doi:10.1364/OE.20.020427

Hierarchical Bayesian estimation improves depth accuracy and spatial resolution of diffuse optical tomography

Takeaki Shimokawa, Takashi Kosaka, Okito Yamashita, Nobuo Hiroe, Takashi Amita, Yoshihiro Inoue, and Masa-aki Sato »View Author Affiliations

Optics Express, Vol. 20, Issue 18, pp. 20427-20446 (2012)
http://dx.doi.org/10.1364/OE.20.020427

View Full Text Article

Acrobat PDF (1521 KB)

Journal Search
Article Lookup

Article Tools

Citations

Alert me when this article is cited
Export Citation/Save

Article Navigation

▸ Article Outline
– Abstract
1. Introduction
2. Methods
3. Phantom experiment
4. Results
5. Discussion
– References and links

Article Tools
Full Text
Article Info
References
Cited By
Figures
Metrics
Download PDF

Viewing Equations in the
HTML Article

Optimal Viewing Recommendations

Full Text
Article Info
References (42)
Cited By
Figures (5)
Metrics

Abstract

High-density diffuse optical tomography (HD-DOT) is an emerging technique for visualizing the internal state of biological tissues. The large number of overlapping measurement channels due to the use of high-density probe arrays permits the reconstruction of the internal optical properties, even with a reflectance-only measurement. However, accurate three-dimensional reconstruction is still a challenging problem. First, the exponentially decaying sensitivity causes a systematic depth-localization error. Second, the nature of diffusive light makes the image blurred. In this paper, we propose a three-dimensional reconstruction method that overcomes these two problems by introducing sensitivity-normalized regularization and sparsity into the hierarchical Bayesian method. Phantom experiments were performed to validate the proposed method under three conditions of probe interval: 26 mm, 18.4 mm, and 13 mm. We found that two absorbers with distances shorter than the probe interval could be discriminated under the high-density conditions of 18.4-mm and 13-mm intervals. This discrimination ability was possible even if the depths of the two absorbers were different from each other. These results show the high spatial resolution of the proposed method in both depth and horizontal directions.

1. Introduction

Diffuse optical tomography (DOT) is an advanced imaging method to increase spatial resolution of near-infrared (NIR) light measurement. It enables a non-invasive and safe monitoring of biological tissues, such as functional neuroimaging [1

A. Villringer, J. Planck, C. Hock, L. Schleinkofer, and U. Dirnagl, “Near infrared spectroscopy (NIRS): a new tool to study hemodynamic changes during activation of brain function in human adults,” Neurosci. Lett. 154, 101–104 (1993). [CrossRef] [PubMed]

–8

C. Habermehl, S. Holtze, J. Steinbrink, S. P. Koch, H. Obrig, J. Mehnert, and C. H. Schmitz, “Somatosensory activation of two fingers can be discriminated with ultrahigh-density diffuse optical tomography,” NeuroImage 59, 3201–3211 (2011). [CrossRef] [PubMed]

] and breast cancer detection [9

J. P. Culver, R. Choe, M. J. Holboke, L. Zubkov, T. Durduran, A. Slemp, V. Ntziachristos, B. Chance, and A. G. Yodh, “Three-dimensional diffuse optical tomography in the parallel plane transmission geometry: evaluation of a hybrid frequency domain/continuous wave clinical system for breast imaging,” Med. Phys. 30, 235–247 (2003). [CrossRef] [PubMed]

–12

Q. Fang, J. Selb, S. A. Carp, G. Boverman, E. L. Miller, D. H. Brooks, R. H. Moore, D. B. Kopans, and D. A. Boas, “Combined optical and X-ray tomosynthesis breast imaging,” Radiology 258, 89–97 (2011). [CrossRef]

]. In measurements, pairs of light sources and detectors are used. NIR light is emitted into the tissue from the source probes and then the light that passes through the tissue is measured by the detector probes. DOT reconstructs the three-dimensional distribution of the optical properties inside the tissue from the measured light intensities by using knowledge of the NIR photon migration process, which is highly diffusive in biological tissues.

NIR light measurement can be classified from two viewpoints: the arrangement of source/detector probes and the instrumentation used. From the first viewpoint, a distinction is made between transmission measurement and reflectance measurement. If NIR light is strong enough to transmit an object, transmission measurement can be made by setting the source probes on one side of the object and the detector probes on the opposite side. However, this is only applicable to thin biological tissues such as those of hands and arms. Otherwise, reflectance measurement, which detects reflected diffusive light by setting source and detector probes on the same side of the object, can be made. Transmission measurement provides richer depth information than does reflectance measurement. From the second viewpoint, a distinction is made between continuous-wave (CW) systems and time-resolved (TR) systems. While a CW system provides only the light intensity of continuous-wave light, a TR system provides a flight time distribution of the incident pulsed light. The TR system provides richer information than does the CW system [13

F. Gao, H. Zhao, and Y Yamada, “Improvement of image quality in diffuse optical tomography by use of full time-resolved data,” Appl. Opt. 41, 778–791 (2002). [CrossRef] [PubMed]

] because the optical paths can be inferred from the flight time in the TR system.

In this paper, we focus on the development of DOT for reflectance and CW measurement because of its potentially wide-ranging applicability. Reflectance measurement is not restricted to the size of the target object. The CW system is low-cost and has a higher frame rate, which makes it possible to measure broad areas of the object almost simultaneously. However, DOT with reflectance and CW measurement is more difficult because it provides less information on inner optical properties, especially depth information, compared with transmission measurement and TR system measurement. Currently, this approach has been investigated in functional neuroimaging studies, and it might also be applicable to breast cancer detection.

Recently, high-density diffuse optical tomography (HD-DOT) has received much attention as a unique solution of DOT for reflectance and CW measurement. This permits three-dimensional reconstruction by using a large number of spatially-overlapping and multi-distance measurement channels by means of high-density probe arrays. Cortical maps of angle and eccentricity in the visual field were measured with a 13-mm probe-interval arrangement [7

B. W. Zeff, B. R. White, H. Dehghani, B. L. Schlaggar, and J. P. Culver, “Retinotopic mapping of adult human visual cortex with high-density diffuse optical tomography,” Proc. Natl. Acad. Sci. U.S.A. 104, 12169–12174 (2007). [CrossRef] [PubMed]

], and a somatosensory activation of two fingers was discriminated with a 7.5-mm probe-interval arrangement [8

Although the application of HD-DOT is promising, there are still problems to solve. First, the difficulty in depth localization causes large quantitative errors in the reconstructed optical properties. Second, the reconstructed image becomes blurred due to the nature of diffusive light. Although the above problems can be mitigated by using a higher-density probe arrangement [8

], it is also necessary to improve the reconstruction method. If these problems were effectively dealt with by improving this method, we could reduce, to a certain extent, the requirements on probe density. This would reduce the cost and preparation time for experiments. The use of spatial variant regularization is considered an effective approach for the first problem, and the introduction of sparsity is considered an effective approach for the second problem. In this study, we incorporate these approaches in the hierarchical Bayesian estimation method [14

M. Guven, B. Yazici, X. Intes, and B. Chance, “Diffuse optical tomography with a priori anatomical information,” Phys. Med. Biol. 50, 2837–2858 (2005). [CrossRef] [PubMed]

,15

F. Abdelnour, C. Genovese, and T. Huppert, “Hierarchical Bayesian regularization of reconstructions for diffuse optical tomography using multiple priors,” Biomed. Opt. Express 1, 1084–1103 (2010). [CrossRef]

Let us explain the first problem in more detail. NIR light intensity exponentially decreases with increased distance from the source at the surface. Accordingly, the amplitude of sensitivity, which relates the inner optical properties to the measured light intensity, is highly non-uniform in the reflectance measurements. It is strong in the neighborhood of the surface where the light sources and detectors are placed but weak in the deep layers. This leads to the reconstructed images becoming biased toward the surface. This systematic depth-localization error causes large quantitative errors because the sensitivity varies exponentially with depth direction. Several approaches have been proposed to overcome this problem. Pogue et al. proposed spatially variant regularization, in which the regularization is weighted by exponential function according to the distance from the surface [16

B. W. Pogue, T. O. McBride, J. Prewitt, U. L. Österberg, and K. D. Paulsen, “Spatially variant regularization improves diffuse optical tomography,” Appl. Opt. 38, 2950–2961 (1999). [CrossRef]

]. Culver et al. proposed similar spatially variant regularization in which the regularization weight is calculated from the sensitivity matrix [17

J. P. Culver, T. Durduran, D. Furuya, C. Cheung, J. H. Greenberg, and A. G. Yodh, “Diffuse optical tomography of cerebral blood flow, oxygenation, and metabolism in rat during focal ischemia,” J. Cereb. Blood Flow Metab. 23, 911–924 (2003). [CrossRef] [PubMed]

]. Niu et al. proposed a depth-compensation algorithm in which the sensitivity matrix in each layer is normalized with the maximum singular value of the sensitivity matrix in that layer [18

H. Niu, F. Tian, Z. J. Lin, and H. Liu, “Development of a compensation algorithm for accurate depth localization in diffuse optical tomography,” Opt. Lett. 35, 429–431 (2010). [CrossRef] [PubMed]

, 19

H. Niu, Z. J. Lin, F. Tian, S. Dhamne, and H. Liu, “Comprehensive investigation of three-dimensional diffuse optical tomography with depth compensation algorithm,” J. Biomed. Opt. 15, 046005 (2010). [CrossRef] [PubMed]

]. These methods are mathematically similar to each other, and they improve the depth localization. However, they are still insufficient if the Tikhonov regularization is used for reconstruction because it cannot localize the solution, which leads to an unwanted compensation to irrelevant regions. The introduction of sparsity, which is explained in the next paragraph, is also required for accurate depth localization.

In regard to the second problem, sparse regularization has recently been applied in optical imaging [20

N. Cao, A. Nehorai, and M. Jacobs, “Image reconstruction for diffuse optical tomography using sparsity regularization and expectation-maximization algorithm,” Opt. Express 15, 13695–13708 (2007). [CrossRef] [PubMed]

–23

S. Okawa, Y. Hoshi, and Y. Yamada, “Improvement of image quality of time-domain diffuse optical tomography with lp sparsity regularization,” Biomed. Opt. Express 2, 3334–3348 (2011). [CrossRef] [PubMed]

]. Previously, the Tikhonov regularization, which is also called the L₂ norm regularization, has been conventionally used in optical imaging because of its ease of implementation and the robustness of this solution. However, the L₂ norm penalty makes the solution overly smoothed because it excessively suppresses strong values. In contrast, the localized solution is obtained by the L₁ norm regularization or another sparse estimation method, which favors the solution with a small number of non-zero elements. Functional neuroimaging and breast cancer detection are considered potential applications of DOT. Task-related brain activities are localized. Tumors in the breast, especially in the early stage, are also localized. For these purposes, accurate localization of activations or tumors is the most crucial issue. However, the solution of the conventional L₂ norm approaches results in localization bias and blurred reconstruction. The sparse estimation method can detect accurate positions with high resolution by using prior information that favors spatially localized patterns. This sparse estimation technique has been investigated for the source-localization problem of electroencephalography (EEG) and magnetoencephalography (MEG) [24

D. Wipf and S. Nagarajan, “A unified Bayesian framework for MEG/EEG source imaging,” Neuroimage 44, 947–966 (2009). [CrossRef]

, 25

F. Lucka, S. Pursiainen, M. Burger, and C. H. Wolters, “Hierarchical Bayesian inference for the EEG inverse problem using realistic FE head models: depth localization and source separation for focal primary currents,” Neuroimage 61, 1364–1382 (2012). [CrossRef] [PubMed]

]. A hierarchical Bayesian model that incorporates an automatic relevance determination (ARD) prior is another way to introduce sparsity. It has been argued that the ARD approach is superior to separable regularization constraints such as L₁ norm regularization [26

D. Wipf and S. Nagarajan, “A new view of automatic relevance determination,” Adv. Neural Inf. Process. Syst. 20, 1625–1632 (2008).

]. Since it has advantages in expandability, it has recently been applied in multi-modality imaging, such as the combination of MEG and fMRI [27

M. Sato, T. Yoshioka, S. Kajihara, K. Toyama, N. Goda, K. Doya, and M. Kawato, “Hierarchical Bayesian estimation for MEG inverse problem,” NeuroImage 23, 806–826 (2004). [CrossRef] [PubMed]

] and that of EEG and fNIRS [28

T. Aihara, Y. Takeda, K. Takeda, W. Yasuda, T. Sato, Y. Otaka, T. Hanakawa, M. Honda, M. Liu, M. Kawato, M. Sato, and R. Osu, “Cortical current source estimation from electroencephalography in combination with near-infrared spectroscopy as a hierarchical prior,” NeuroImage 59, 4006–4021 (2012). [CrossRef]

]. In this study, we introduce sparsity by using the ARD prior in the hierarchical Bayesian model.

We developed a three-dimensional reconstruction method by introducing the above two approaches, that is, spatial variant regularization and sparsity. Then, we examined the proposed method by applying it to the data obtained by a phantom experiment measured with a commercially available near-infrared imaging system. We evaluated the depth accuracy and the spatial resolution of the proposed method by placing one or two absorbers in the phantom. The evaluation was done in the three interval conditions of the probe arrangement (i.e. 26 mm, 18.4 mm, and 13 mm) to investigate the relationship between probe density and estimation accuracy. As a result, we found that the depth of the absorber can be estimated accurately using the proposed method in the high-density arrangements of 18.4-mm and 13-mm probe intervals. It was also possible to discriminate two absorbers whose distance is shorter than the probe interval. These results demonstrate that the proposed method allows three-dimensional reconstruction with high resolution in both depth and horizontal directions.

The following notations are used throughout this paper. Scalars are represented in normal-type letters as x. Vectors are represented in lower-case boldface letters as x. Matrices are represented in upper-case boldface letters as X. The element in the i-th row and the j-th column of matrix X is denoted as (X)_ij = X_ij. The transpose, trace, and determinant of matrix X are denoted as X^T, tr(X), and |X|, respectively. A diagonal matrix with diagonal elements x is denoted as diag(x). A vector consisting of the diagonal elements of matrix X is denoted as diag(X). The L₂ norm of vector x is denoted as

{‖ x ‖}_{2}^{2}

. The mean value

\frac{1}{T} \sum_{t = 1}^{T} x_{t}

is denoted as mean(x). The N-dimensional vector consisting of 1s and the identity matrix of size N are denoted as 1_N and I_N, respectively. N(x; μ, Σ) denotes that a vector random variable x follows a multivariate Gaussian distribution with mean vector μ and covariance matrix Σ. Γ(x;μ, γ) denotes that a scalar random variable x follows a gamma distribution with mean μ and shape parameter γ (note that the scale parameter corresponds to μγ⁻¹).

2. Methods

In this DOT study, we reconstruct a three-dimensional distribution of the absorption change in a scattering medium from the observed change in light intensity. First, we describe the relationship between the absorption change and light intensity change in the forward problem. Then, we estimate the distribution of the absorption change given the observed light intensity change in the inverse problem using the relationship calculated in the forward problem.

2.1. Forward problem

In this work, we used continuous-wave intensity data to reconstruct absorption changes. The photon fluence in the baseline condition Φ₀ changes to Φ in the presence of absorption changes in the scattering medium. Under the condition that the absorption changes are small, the measured changes in photon fluence ϕ_pert (r_s, r_d) = − ln(Φ(r_s, r_d)/Φ₀(r_s, r_d)) are linearly related to the spatial variations in the absorption coefficient δμ_a(r) by the Rytov approximation [29

A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging (IEEE Press, New York, 1988).

–32

S. R. Arridge and J. C. Schotland, “Optical tomography: forward and inverse problems,” Inverse Probl. 25, 123010 (2009). [CrossRef]

ϕ_{pert} (r_{s}, r_{d}) = \int d r \frac{Φ_{0} (r_{s}, r) Φ_{0} (r, r_{d})}{Φ_{0} (r_{s}, r_{d})} δ μ_{a} (r),

(1)

where r_s and r_d denote the positions of source and detector, respectively. Throughout this study, we discretized the scattering medium into 2.5mm × 2.5mm × 2.5mm voxels. Thus, Eq. (1) can be written in a matrix form as y = Ax, where x = (x₁, ⋯ ,x_N)^T represents the absorption changes x_i = δμ_a(r_i), y = (y₁, ⋯ ,y_M)^T represents the measured light intensity changes y_j = ϕ_pert (r_s₍_j₎, r_d₍_j₎), and A represents the sensitivity matrix

A_{j i} = \frac{Φ_{0} (r_{s (j)}, r_{i}) Φ_{0} (r_{i}, r_{d (j)})}{Φ_{0} (r_{s (j)}, r_{d (j)}))}

. Here, the index i denotes the ith voxel (i = 1, ⋯ ,N) and the index j denotes the jth measurement logical channel (j = 1, ⋯ ,M). The value of Φ₀ can be calculated by using the Monte Carlo simulation [33

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). [CrossRef] [PubMed]

] or solving the diffusion equation. In this study, we obtained the value by solving the diffusion equation for a one-layer model [34

R. C. Haskell, L. O. Svaasand, T. T. Tsay, T. C. Feng, M. S. McAdams, and B. J. Tromberg, “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A 11, 2727–2741 (1994). [CrossRef]

In a measurement duration, a number of sample data were obtained. We summarized the observed light intensity changes in a matrix form as Y = (y₁, ⋯ , y_t, ⋯ , y_T), where T denotes the number of samples. Likewise, the absorption changes are summarized as X = (x₁, ⋯ , x_t, ⋯ , x_T).

2.2. Inverse problem

The inverse problem in which the absorption change x is estimated from the observed light intensity change y is ill-posed, since the number of unknowns N is typically much larger than the number of data M. Furthermore, three-dimensional reconstruction is very difficult due to the characteristics of the sensitivity matrix. Sensitivity is exponentially decaying with depth [16

B. W. Pogue, T. O. McBride, J. Prewitt, U. L. Österberg, and K. D. Paulsen, “Spatially variant regularization improves diffuse optical tomography,” Appl. Opt. 38, 2950–2961 (1999). [CrossRef]

], and thus absorption changes tend to be estimated in superficial layers. In order to overcome this difficulty, we introduce sensitivity-normalized regularization and sparsity for better localization. The proposed estimation method consists of two steps. First, we obtain the solution with the sensitivity-normalized Tikhonov regularization. Second, we refine the solution with the hierarchical Bayesian estimation method.

2.2.1. Sensitivity-normalized Tikhonov regularization

The Tikhonov regularization is conventionally used to solve this inverse problem, in which the following cost function is minimized [35

C. M. Bishop, Pattern Recognition and Machine Learning (Springer, New York, 2006).

C (x_{t}; λ) = {‖ Σ_{y}^{- \frac{1}{2}} (y_{t} - A x_{t}) ‖}_{2}^{2} + λ {‖ x_{t} ‖}_{2}^{2},

(2)

where Σ_y is the measurement noise covariance matrix. The first term denotes the sum of squared error, and the second term

{‖ x_{t} ‖}_{2}^{2} = \sum_{i = 1}^{N} x_{i t}^{2}

denotes L₂ norm regularization to avoid overfitting. λ is the regularization parameter controlling the relative importance of the regularization term compared with the sum of squared error term. In this study, we determined the appropriate λ value by minimizing Akaike’s Bayesian information criterion (or maximizing the marginal likelihood) [36

H. Akaike, “Likelihood and the Bayes procedure,” in Bayesian Statistics , J. M. Bernardo, M. H. De Groot, D. V. Lindley, and A. F. M. Smith, eds. (Univ. Press, Valencia, 1980), 143–166.

Here, we propose sensitivity-normalized Tikhonov regularization that introduces the normalization weight matrix D in the regularization term to balance the error term and the regularization term in each voxel as

C_{D} (x_{t}; λ) = {‖ Σ_{y}^{- \frac{1}{2}} (y_{t} - A x_{t}) ‖}_{2}^{2} + λ {‖ D^{\frac{1}{2}} x_{t} ‖}_{2}^{2},

(3)

where D is a diagonal matrix defined as the sum of two terms: the sensitivity normalization term, which is the diagonal matrix of the product of the sensitivity matrix and the measurement noise covariance matrix, and the constant term

D = diag (ρ + β 1_{N}),

(4)

ρ_{i} = {(A^{T} Σ_{y}^{- 1} A)}_{i i},

(5)

where ρ = (ρ₁, ⋯ ,ρ_N)^T is the sensitivity normalization term. β is a positive control constant. D has positive diagonal elements because the noise covariance matrix Σ_y is positive definite. The constant term is introduced in order to make the regularization penalty remain constant in deep layers.

Sensitivity-normalized regularization works strongly in superficial layers and weakly in deep layers, since the sensitivity exponentially decays with depth [16

B. W. Pogue, T. O. McBride, J. Prewitt, U. L. Österberg, and K. D. Paulsen, “Spatially variant regularization improves diffuse optical tomography,” Appl. Opt. 38, 2950–2961 (1999). [CrossRef]

]. This weight matrix is similar to that of Culver et al. [17

, 37

J. Selb, A. M. Dale, and D. A. Boas, “Linear 3D reconstruction of time-domain diffuse optical imaging differential data: improved depth localization and lateral resolution,” Opt. Express 15, 16400–16412 (2007). [CrossRef] [PubMed]

, 38

H. Dehghani, B. R. White, B. W. Zeff, A. Tizzard, and J. P. Culver, “Depth sensitivity and image reconstruction analysis of dense imaging arrays for mapping brain function with diffuse optical tomography,” Appl. Opt. 48, D137–D143 (2009). [CrossRef] [PubMed]

] except for the measurement noise covariance matrix term. This is also similar to the weight matrix of Niu et al. [18

], but they introduced normalization for each layer. The proposed sensitivity normalization compensates for the non-uniformity of the sensitivity matrix, not only for the depth direction but also for the horizontal direction.

In this study, the control constant β was determined with reference to the maximum value of sensitivity for the region deeper than 22.5 mm from the surface: β = max(ρ_i), i ∈ I_≥22.5_mm, where I_≥22.5_mm denotes the set of voxel indices deeper than 22.5 mm, i.e. maximum is taken over the voxels deeper than 22.5 mm.

In the case of L₂ norm regularization, the solution x̂_D can be obtained in an explicit form:

{\hat{x}}_{t, D} = {(A^{T} Σ_{y}^{- 1} A + λ D)}^{- 1} A^{T} Σ_{y}^{- 1} y_{t} .

(6)

The second term in brackets represents the regularization term. By introducing the constant term β, the weight D is almost the same as β in layers deeper than 22.5 mm: D_ii ∼ β ≫ ρ_i. Consequently, a solution for layers deeper than 22.5 mm is suppressed by the regularization term. In contrast, the weight D is almost the same as the normalization term in layers shallower than 22.5 mm: D_ii ∼ ρ_i ≫ β. Therefore, the amplitudes of the first and second terms in the brackets of Eq. (6) are balanced. Accordingly, the solution is appropriately regularized regardless of the exponentially non-uniform sensitivity.

2.2.2. Hierarchical Bayesian estimation

In this section, we introduce the hierarchical Bayesian method in order to refine the solution x̂_D, using it as an initial value. As already noted in the previous section, the Tikhonov regularization penalizes the solution uniformly at every position. Sensitivity-normalized regularization improves the solution by using a weighted penalty according to the sensitivity value at each position. However, it is still insufficient to localize the solution, which is biased as examined in section 4.1. The hierarchical Bayesian estimation method introduces the regularization parameters λ_i at each voxel position by an automatic relevance determination (ARD) prior [35

C. M. Bishop, Pattern Recognition and Machine Learning (Springer, New York, 2006).

, 39

A. C. Faul and M. E. Tipping, “Analysis of sparse Bayesian learning,” Adv. Neural Inf. Process. Syst. 14, 383–389 (2002).

]. These parameters λ_i, which control the degree of penalty, are automatically adjusted from the data. As a result, the solution becomes sparse and high-resolution. However, the estimation of λ_i values becomes a nonlinear optimization, which requires an appropriate initial value. Therefore, we set the solution of the sensitivity-normalized regularization x̂_D as an initial value of λ_i. Note that the result of the hierarchical Bayesian method becomes biased toward the surface if the solution of the conventional Tikhonov regularization is set as an initial value.

In the hierarchical Bayesian estimation method, observation and regularization are described as hierarchical probabilistic models. The following models and algorithm were originally developed in the context of the electro-magnetic source localization problem [27

M. Sato, T. Yoshioka, S. Kajihara, K. Toyama, N. Goda, K. Doya, and M. Kawato, “Hierarchical Bayesian estimation for MEG inverse problem,” NeuroImage 23, 806–826 (2004). [CrossRef] [PubMed]

]. In this paper, we, for the first time, introduce the ARD prior for CW-DOT and propose combining its sparse-promoting nature with the initial value obtained by the sensitivity-normalized regularization method. This combination makes it possible to obtain tomography with high spatial resolution and high depth accuracy.

First, the observation model is constructed assuming a Gaussian noise in measuring light intensity change. This assumption is validated by the phantom experimental data. We also assume that the noise covariance matrix in the absorption change measurement is proportional to the noise covariance matrix measured under the baseline condition Σ_y. Then, the conditional probability for the observed light-intensity change Y given the absorption change X is given by

P (Y | X, σ) = \prod_{t = 1}^{T} N (y_{t}; A x_{t}, σ^{- 1} Σ_{y}) .

(7)

where σ is the scale parameter inversely proportional to the measurement noise variance. The scale parameter σ is estimated by introducing a hierarchical prior, as

P (σ) \propto σ^{- 1} .

(8)

Second, the regularization model is constructed by adopting the ARD prior. In order to prevent the solution from becoming too sparse, a spatially localized smoothing filter W is incorporated. We assume that the solution x is represented as x = Wz, where z denotes the sparse expression of the solution. The ARD prior is adopted for z. The probability distribution of z is given by a Gaussian distribution with a zero mean vector and a diagonal covariance matrix whose diagonal elements are adjustable hyperparameters:

P (Z | λ) = \prod_{t = 1}^{T} N (z_{t}; 0, Λ^{- 1}),

(9)

where Λ = diag(λ), and λ = (λ₁, λ₂, ⋯ ,λ_N)^T represent the adjustable hyperparameters. In this study, a 5-mm FWHM Gaussian filter was used as W. Equation (9) can be equivalently rewritten as

P (X | λ) = \prod_{t = 1}^{T} N (x_{t}; 0, W Λ^{- 1} W^{T}) .

(10)

The hyperparameter λ is estimated by introducing a hierarchical prior:

P (λ) = \prod_{i = 1}^{N} Γ (λ_{i}; {\bar{λ}}_{0 i}, γ_{0 i}),

(11)

where λ̄₀_i represents the mean value of λ_i in the hierarchical prior, and γ₀_i represents the reliability (or confidence) of the hierarchical prior [27

M. Sato, T. Yoshioka, S. Kajihara, K. Toyama, N. Goda, K. Doya, and M. Kawato, “Hierarchical Bayesian estimation for MEG inverse problem,” NeuroImage 23, 806–826 (2004). [CrossRef] [PubMed]

In the Bayesian estimation, the posterior probability of the unknown X, λ, σ given the observation data Y is written as

P (X, λ, σ | Y) = \frac{P (Y, X, λ, σ)}{P (Y)} = \frac{P (Y | X, σ) P (X | λ) P (λ) P (σ)}{\int d X d λ d σ P (Y | X, σ) P (X | λ) P (λ) P (σ)} .

(12)

The calculation of the posterior cannot be done analytically. Accordingly, the posterior calculation problem is converted to the maximization problem of the free energy function. The free energy for a trial distribution Q(X, λ, σ), which is the approximation of the posterior P(X, λ, σ|Y), is defined by

\begin{array}{l} F (Q) & = & \int d X d λ d σ Q (X, λ, σ) log \frac{P (Y, X, λ, σ)}{Q (X, λ, σ)}, \\ = & log P (Y) - KL [Q (X, λ, σ) | | P (X, λ, σ | Y)], \end{array}

(13)

where

KL [P (θ) | | Q (θ)] = \int d θ P (θ) log \frac{P (θ)}{Q (θ)}

denotes the Kullback-Leibler distance. Equation (13) implies that the maximization of the free energy F(Q) is equivalent to the minimization of the Kullback-Leibler distance between the posterior and the trial distribution.

The maximization of the free energy can be solved using a factorization approximation restricting the solution space [40

H. Attias, “Inferring parameters and structure of latent variable models by variational Bayes,” Proc. 15th Conf. on Uncertainty in Artificial Intelligence , Morgan Kaufmann, 21–30 (1999).

, 41

M. Sato, “Online model selection based on the variational Bayes,” Neural Comput. 13, 1649–1681 (2001). [CrossRef]

Q (X, λ, σ) = Q_{X} (X) Q_{λ} (λ, σ) .

(14)

Under the factorization assumption, the maximum free energy is obtained by alternately maximizing the free energy with respect to Q_X and Q_λ. The convergence of the free energy was guaranteed by the previous theoretical studies [26

D. Wipf and S. Nagarajan, “A new view of automatic relevance determination,” Adv. Neural Inf. Process. Syst. 20, 1625–1632 (2008).

, 40

H. Attias, “Inferring parameters and structure of latent variable models by variational Bayes,” Proc. 15th Conf. on Uncertainty in Artificial Intelligence , Morgan Kaufmann, 21–30 (1999).

, 41

M. Sato, “Online model selection based on the variational Bayes,” Neural Comput. 13, 1649–1681 (2001). [CrossRef]

]. In the first step (X-step), the free energy F(Q) is maximized with respect to Q_X while Q_λ is fixed. The solution Q_X (X) becomes a normal distribution [27

M. Sato, T. Yoshioka, S. Kajihara, K. Toyama, N. Goda, K. Doya, and M. Kawato, “Hierarchical Bayesian estimation for MEG inverse problem,” NeuroImage 23, 806–826 (2004). [CrossRef] [PubMed]

log Q_{X} (X) = \int d λ d σ Q_{λ} (λ, σ) log P (Y, X, λ, σ) = \sum_{t = 1}^{T} log N (x_{t}; {\bar{x}}_{t}, Σ_{x}),

(15)

where x̄_t and Σ_x denote the expectation value and covariance matrix of x_t, respectively. In the second step (λ-step), the free energy F(Q) is maximized with respect to Q_λ while Q_X is fixed. The solution Q_λ (λ, σ) becomes a gamma distribution:

log Q_{λ} (λ, σ) = \int d X Q_{X} (X) log P (Y, X, λ, σ) = \sum_{i = 1}^{N} log Γ (λ_{i}; {\bar{λ}}_{i}, γ_{i}) + log Γ (σ; \bar{σ}, γ_{σ}),

(16)

where λ̄_i and γ_i denote the expectation and shape parameter values of λ, respectively, and σ̄ and γ_σ denote the expectation and shape parameter values of σ, respectively. The above X- and λ-steps are repeated until the free energy converges (see Appendix).

In this procedure, the free energy function can be expressed as

F (Q) = \sum_{t = 1}^{T} F_{t} ({\bar{x}}_{t}, \bar{λ}, \bar{σ}) + f (\bar{λ}, \bar{σ}, Σ_{x}),

(17)

F_{t} ({\bar{x}}_{t}, \bar{λ}, \bar{σ}) = - \frac{1}{2} [\bar{σ} {‖ Σ_{y}^{- \frac{1}{2}} (y_{t} - A {\bar{x}}_{t}) ‖}_{2}^{2} + {‖ {\bar{Λ}}^{\frac{1}{2}} W^{- 1} {\bar{x}}_{t} ‖}_{2}^{2}] .

(18)

In Eq. (17), f (λ̄, σ̄, Σ_x) represents a sparse-promoting prior function. Its detailed form is explained in the Appendix.

In the Tikhonov regularization of Eq. (2), the solution is uniformly suppressed by the single parameter λ shared by all voxel positions. In contrast, the regularization parameters λ = (λ₁, λ₂, ⋯ ,λ_N)^T are appropriately adjusted at every voxel position in this hierarchical Bayesian estimation method. Therefore, a sparse solution is obtained [35

C. M. Bishop, Pattern Recognition and Machine Learning (Springer, New York, 2006).

, 39

A. C. Faul and M. E. Tipping, “Analysis of sparse Bayesian learning,” Adv. Neural Inf. Process. Syst. 14, 383–389 (2002).

In this study, the reliability value of the hierarchical prior was set as γ₀ = (γ₀₁, ⋯ , γ₀_N)^T = 0. Thus, the prior of Eq. (11) becomes non-informative, and the value λ₀_i does not affect the probabilistic model and the iterative algorithm. In this case, z becomes sparse [35

C. M. Bishop, Pattern Recognition and Machine Learning (Springer, New York, 2006).

, 39

A. C. Faul and M. E. Tipping, “Analysis of sparse Bayesian learning,” Adv. Neural Inf. Process. Syst. 14, 383–389 (2002).

], which means that the elements of z converge to zeros in most of the voxels and non-zeros in the remaining small number of voxels. Therefore, the solution x becomes a mixture of a small number of local activities smoothed by W. In numerical calculation, z becomes very small values in most of the voxels if the algorithm is terminated by some convergence criteria. In the case where informative prior γ₀_i > 0, z does not converge to zeros because of the prior information. However, if γ₀_i is small, the values of z become numerically very small in most of the voxels. We incorporated the result of the sensitivity-normalized Tikhonov regularization as an initial value of the iterative algorithm as λ̄_init = (10 · mean(x̂_D))⁻² (see Appendix). The ten-times value of mean(x̂_D) was used because the solution of the Tikhonov regularization is overly smoothed.

3. Phantom experiment

3.1. Experimental conditions

We evaluated the estimation accuracy of the proposed method using a liquid phantom experiment as shown in Fig. 1. The tank was filled with a solution of 1.4% intralipid and ink to make the optical properties of the liquid similar to that of gray and white matter [6

D. A. Boas and A. M. Dale, “Simulation study of magnetic resonance imaging-guided cortically constrained diffuse optical tomography of human brain function,” Appl. Opt. 44, 1957–1968 (2005). [CrossRef] [PubMed]

, 7

]. The estimated absorption and reduced scattering coefficients are μ_a = 0.019mm⁻¹ and μ_s′ = 1.1mm⁻¹, respectively, based on measurement with a spectral photometer and a near-infrared imaging system.

Fig. 1 Setup of liquid phantom experiment. (a) Schematic of the phantom. (b) Experimental device. (c) Source and detector probes are on the bottom of the tank, z = 0. (d) The x–y positions of the absorber in measurement: (i) center (0, 0), (ii) midpoint (

- \frac{l}{2}

, 0), (iii) source (

- \frac{l}{2}

- \frac{l}{2}

), (iv) detector (

- \frac{l}{2}

+ \frac{l}{2}

), (v) second absorber position (

- \frac{l}{2} + d, 0

). Here, l denotes the probe interval and d denotes the horizontal distance between the centers of the two absorbers. The red circle represents the source and the blue square represents the detector.

We used a continuous-wave near-infrared imaging system (FOIRE-3000, Shimadzu Corporation) to acquire the data for the reconstruction of the absorption change. We used only 805-nm wavelength data. The source and detector probes are placed on the bottom alternately in a rectangular arrangement as shown in Fig. 1(c). The sampling rate is 130 ms. We measured for 20 seconds each time, with and without absorbers. We tested the estimation results under three conditions of probe interval: l = 13 mm, 18.4 mm (

= 13 mm \times \sqrt{2}

), 26 mm (=13 mm×2).

In experiments, one or two spherical absorbers of 5 mm in diameter were placed in the liquid phantom. These are black absorbers made of polyacetal (Fig. 1(b)). First, we evaluated the positional accuracy of the proposed method using one absorber placed in four x–y positions (see section 4.2): (i) center of the four probes, (ii) midpoint of the source and the detector, (iii) beneath the source, (iv) beneath the detector, as shown in Fig. 1(d). Second, we evaluated the spatial resolution of the proposed method using two absorbers. The first absorber is placed in (ii) and the second absorber is placed in (v), separated by d from (ii) (see section 4.3). Finally, we evaluated the depth accuracy in a more general setting where the two absorbers are placed in (ii) and (v) at different depths (see section 4.4). One example of this experiment is shown first in section 4.1 for understanding the main idea of the proposed method.

The number of voxels N was 7500 = 25 × 25 × 12 (we assumed that absorption changes occur inside the region of 62.5mm × 62.5mm × 30mm). The number of measurement logical channels M was 56 (l = 13 mm, up to 3rd nearest-neighbour(NN) channels), 48 (l = 18.4 mm, up to 2nd NN channels), or 24 (l = 26 mm, up to 1st NN channels).

4. Results

4.1. Example of 3D reconstruction

In this section, we first explain the main idea of the proposed method by applying it to one example of a 3D reconstruction in the phantom experiment. The example data was obtained by the liquid phantom experiment with an arrangement of 18.4-mm probe intervals. In this case, two absorbers were used to generate the absorption change. The positions of the centers of the spherical absorbers were (ii) (−9.2, 0, 10) [mm] and (v) (3.3, 0, 15) [mm] as shown in Fig. 2(a), where the horizontal distance d is 12.5 mm and the vertical distance is 5 mm. Therefore, the distance between the two absorbers is shorter than the probe interval of 18.4 mm. Note that the coordinate system of Fig. 2 is flipped upside down to that of Fig. 1 in order to make it similar to that of functional neuroimaging, in which the probes are usually placed on top of the surface.

Fig. 2 Example of 3D reconstruction image. The three-dimensional image (tomography) is represented in gray scale and the x–y maximum intensity projection (topography) is represented in green scale on the z = 0 surface. (a) True distribution of the absorption change. The positions of the centers of the two absorbers are (−9.2, 0, 10) [mm] and (3.3, 0, 15) [mm]. The red cylinder represents the source and the blue square prism represents the detector. (b) Solution of the Tikhonov regularization with spatially uniform regularization. (c) Solution of the Tikhonov regularization with sensitivity-normalized regularization. (d) Solution of the hierarchical Bayesian estimation method using the result of (c) as an initial value.

Figure 2(b) shows the estimation result of the Tikhonov regularization with spatially uniform regularization, where D is the identity matrix. The solution is biased toward superficial layers because the sensitivity is large in superficial layers while the effect of the regularization suppression is uniform and independent of the depth. This depth localization error makes the reconstructed absorption change too small. Moreover, the solution is overly smoothed and the two absorbers cannot be distinguished.

Figure 2(c) shows the estimation result of the sensitivity-normalized Tikhonov regularization. Although the solution looks biased toward deep layers in contrast to the previous result of Fig. 2(b), this seems to be a reasonable result. Typically, a weak absorption change in superficial layers and a strong absorption change in deep layers generate similar light intensity changes. Since the Tikhonov regularization method cannot localize the solution, the solution might become a mixture state of a superficial weak change and a deep strong change. In the figure, the deep strong values are more visible than the superficial weak values. Values deeper than 22.5 mm are suppressed by the introduction of the control constant β as explained in section 2.2.1.

Figure 2(d) shows the estimation result of the hierarchical Bayesian estimation method starting from the sensitivity-normalized Tikhonov regularization solution. By using the Bayesian algorithm incorporating the automatic relevance determination prior, the appropriate solution is selected from the mixture state solution. The two absorbers can be accurately distinguished and the positional error is smaller than the voxel size of 2.5 mm. This result is summarized again in section 4.4 (l = 18.4 mm, d = 12.5 mm, z = 10 & 15 mm).

4.2. One-absorber experiment

In this section, we evaluate the horizontal and depth estimation accuracies of the proposed method by using a one-absorber setting. We first evaluated the estimation results at various depths of a spherical absorber’s center, z = 7.5 mm, 10 mm, 12.5 mm, 15 mm, 17.5 mm, 20 mm, 22.5 mm, while the horizontal position was fixed. The analysis was done under three conditions of probe interval: l = 13 mm, 18.4 mm, 26 mm.

First, the absorber was placed in center position (i) (Fig. 1(d)). Figure 3(a) shows the estimation results for various depth and probe-interval conditions. Here, the 3D estimation results are represented as the 2D x–z projection images to clearly visualize the estimation accuracy. In the figure, the maximum intensity of the three central voxel regions in the y direction (7.5-mm-thick layer in Fig. 3(b)), including the true position of the 5-mm-diameter spherical absorber, was projected onto the x–z plane in order to check the positional accuracy in the y direction. If a positional error occurs in the y direction, the estimation value in this region becomes nearly zero. The white circle in the projection images in Fig. 3(a) represents the true position of the absorption change. We judged the estimation results as successful when the positional error between the center position of the true absorption change and the voxel center of the maximum estimation value was within one voxel (2.5 mm) in each x, y, z direction and when the maximum estimation value was greater than 0.025mm⁻¹. The estimation image is surrounded by a bold frame if the estimation result was judged as successful. Figure 3(c) represents the depth of the maximum estimation value exceeding 0.025mm⁻¹ as a function of the true depth of the absorption change.

Fig. 3 Estimation results of one-absorber experiment. (a) Estimation results under various depth and probe-interval conditions. From left, results under the depth conditions of the absorber: z = 7.5 mm, 10 mm, 12.5 mm, 15 mm, 17.5 mm, 20 mm, 22.5 mm. From top, results under three conditions of the probe interval: l = 13 mm, 18.4 mm, 26 mm. The white circle represents the true position of the absorption change. The estimation result is surrounded by a bold frame if the positional error was within one voxel (2.5 mm) in each x, y, z direction and the maximum estimation value was greater than 0.025mm⁻¹. (b) The projection region is the central 7.5-mm-thick layer including the true distribution of the absorption change. The estimation images of (a) are obtained by projecting the maximum estimation value within the projection region in the y direction onto the x–z space. (c) Estimated depth as a function of the true depth.

The first line from the top of Fig. 3(a) shows the estimation results for the 26-mm probe-interval condition. The depths of the maximum estimation value are the same regardless of the variations in the true depth conditions from z = 7.5 mm to z = 15 mm (see also Fig. 3(c)). Furthermore, the amplitudes of the estimated absorption change strongly vary depending on the true depth conditions. These results indicate that the depth of the absorption change cannot be estimated at all in the 26-mm probe-interval condition. In the 26-mm interval condition, we could use only 1st nearest-neighbour (NN) channels because 2nd NN separation, which is 58.1 mm, is too distant to measure the light intensity changes [38

]. Therefore, the depth information of the absorption change could not be obtained without overlapping measurements [38

]. Moreover, no absorption change deeper than z = 17.5 mm could be detected.

The second line from the top of Fig. 3(a) shows the estimation results for the 18.4-mm probe-interval condition. The absorption changes are accurately estimated from z = 7.5 mm to z = 20 mm, where the positional errors are within one voxel (see also Fig. 3(c)). The amplitudes of the estimated absorption change are maintained at a roughly constant level. In the 18.4-mm interval condition, we could use the overlapping measurement channels of the 1st and 2nd NN channels, whose separations are 18.4 mm and 41.1 mm, respectively. No absorption change deeper than z = 22.5 mm could be detected.

The third line from the top of Fig. 3(a) shows the estimation results for the 13-mm probe-interval condition. The absorption changes are accurately estimated from z = 7.5 mm to z = 22.5 mm, where the positional errors are within one voxel (see also Fig. 3(c)). The amplitudes of the estimated absorption change are maintained roughly constant. In the 13-mm interval condition, we could use the 1st, 2nd, and 3rd NN channels, whose separations are 13 mm, 29.1 mm and 39 mm, respectively. No absorption change deeper than z = 25 mm could be detected. These results show that, with overlapping measurement channels, the depth of the absorption change is accurately estimated.

Next, the estimation accuracy with respect to horizontal position was evaluated by applying the same analyses in four conditions of the x–y absorber position: (i) center, (ii) midpoint, (iii) source, and (iv) detector (Fig. 1(d)). By deepening the spherical absorber by 2.5 mm from z = 7.5 mm, we investigated the depth limit where successful estimation is possible, which means that the positional error is within one voxel (2.5 mm) in each x, y, z direction and the maximum detected value is greater than 0.025mm⁻¹. The summarized data are shown in Table 1 under the various conditions of the x–y absorber position and the probe interval. Note that in the 26-mm probe-interval condition, the estimation failed in every case because the depth could not be estimated at all. The value in the parentheses represents the depth limit where positional accuracy is evaluated with respect to only the x, y direction for reference.

Table 1 Depth limit of spherical absorber’s position where accurate estimation is possible.

probe interval	26mm	18.4mm	13mm

absorber position	depth limit [mm]
(i) center	− (^*15.0)	20.0	22.5
(ii) midpoint	− (^*17.5)	20.0	22.5
(iii) source	− (^*17.5)	17.5	20.0
(iv) detector	− (^*17.5)	17.5	20.0

^* Value in parentheses represents the depth limit where positional accuracy is evaluated with respect to only the x, y direction for reference.

In the 26-mm interval condition, sensitivity is weakest in the center (i). In contrast, in the 18.4-mm and 13-mm interval conditions, in which overlapping measurement is possible, estimation accuracy is relatively weak at just below the source (iii) and the detector (iv).

As the probe interval becomes shorter, a deeper absorption change can be detected. This is because there are measurement channels in suitable positions owing to the high-density arrangement and the increasing number of simultaneous independent observations by the high-density overlapping measurements.

4.3. Two-absorber experiment 1 (at the same depth)

In this section, we evaluate the spatial resolution of the proposed method by a two-absorber experiment. In this phantom experiment, two spherical absorbers were placed at the same depth in (ii) (

- \frac{l}{2}

, 0, z) and (v) (

- \frac{l}{2} + d

, 0, z) (Fig. 1(d)), where d denotes the horizontal distance, z denotes the depth of the two absorbers, and l denotes the probe interval. We evaluated the estimation results under two conditions of probe interval, 18.4 mm and 13 mm, because the 26-mm probe-interval condition cannot obtain good spatial resolution due to the lack of overlapping measurements.

Figure 4 shows the estimation results for the 18.4-mm probe-interval condition. The estimations are tested under four conditions of the horizontal distance, d = 17.5 mm, 15 mm, 12.5 mm, 10 mm, and under four conditions of the depth, z = 12.5 mm, 15 mm, 17.5 mm, 20 mm. The 3D estimation results are represented as the 2D x–z intensity projection images in the same way as Fig. 3(a). The white circles represent the true positions of the two absorbers. We judged the estimation results as successful if the positional errors of both absorbers were within one voxel (2.5 mm) in each x, y, z direction and both maximum estimation values were greater than 0.025mm⁻¹. The estimation image is surrounded by a bold frame if the estimation result was judged as successful.

Fig. 4 Estimation results for the 18.4-mm probe-interval condition. From top, the results under four conditions of horizontal distance: d = 17.5 mm, 15 mm, 12.5 mm, 10 mm. From left, the results under four conditions of depth: z = 12.5 mm, 15 mm, 17.5 mm, 20 mm. The white circles represent the true positions of the absorption change. The estimation image is surrounded by a bold frame if the positional errors of both absorbers were within one voxel (2.5 mm) in each x, y, z direction and both maximum estimation values were greater than 0.025mm⁻¹.

In the case of the horizontal distance d = 17.5 mm, which is almost the same as the probe interval, the estimation was accurate until the depth condition of z = 17.5 mm as shown in Fig. 4 (d = 17.5 mm, z = 17.5 mm). In the depth condition of z = 20 mm, which corresponds to the detectable depth limit as summarized in Table 1, the two absorbers could not be distinguished.

Even in the case of the horizontal distance d = 10 mm, which is nearly half the probe interval, and where the clearance distance of the two spherical absorbers is 5 mm, the estimation was accurate until the depth condition of z = 15 mm as shown in Fig. 4 (d = 10 mm, z = 15 mm).

The same analysis was also done under the 13-mm probe-interval condition. Table 2 summarizes the depth limits where successful estimation is possible.

Table 2 Depth limit where the two absorbers are accurately estimated.

probe interval	18.4mm	13mm

horizontal distance [mm]	depth limit [mm]
15.0	17.5	17.5
12.5	15.0	17.5
10.0	15.0	17.5

The estimation results of the 13-mm probe-interval condition were better than those of the 18.4-mm condition. However, the depth limit z = 17.5 mm in the horizontal distance d = 15 mm seems unsatisfactory compared to the detectable depth limit z = 22.5 mm as summarized in Table 1 (13-mm probe interval). In the 4 × 4 rectangular arrangement we used, the numbers of 1st NN, 2nd NN, and 3rd NN channels are 24, 24, and 8, respectively. The number of the 3rd NN channels, which are sensitive to deep absorption changes, would not be high enough to achieve a high spatial resolution in deep layers. Additionally, the estimation results may suffer from the nonlinear effect, especially under the 13-mm probe-interval condition due to the use of the strong absorber, which will be discussed in section 5.

4.4. Two-absorber experiment 2 (at different depths)

Next, we conducted a two-absorber experiment at different depths to confirm the depth accuracy of the proposed method in more general settings. The two spherical absorbers were placed at depths differing by 5 mm in (ii) (

- \frac{l}{2}

, 0, z −5.0) and (v) (

- \frac{l}{2} + d

, 0, z) (Fig. 1(d)), where d denotes the horizontal distance and l denotes the probe interval. We evaluated the estimation results under two conditions of probe interval, 18.4 mm and 13 mm.

Figure 5 shows the estimation results under the 18.4-mm probe-interval condition. The estimations are tested under four conditions of horizontal distance, d = 17.5 mm, 15 mm, 12.5 mm, 10 mm, and under four conditions of depth, z = 7.5 & 12.5 mm, 10 & 15 mm, 12.5 & 17.5 mm, 15 & 20 mm. The 3D estimation results are represented as the 2D x–z intensity projection images in the same way as Fig. 3(a) and Fig. 4. The white circles represent the true positions of the two absorbers. We judged the estimation results as successful if the positional errors of both absorbers were within one voxel (2.5 mm) in each x, y, z direction and both maximum estimation values were greater than 0.025mm⁻¹. The estimation image is surrounded by a bold frame if the estimation result was judged as successful.

Fig. 5 Estimation results for the 18.4-mm probe-interval condition. From top, the results under four conditions of horizontal distance: d = 17.5 mm, 15 mm, 12.5 mm, 10 mm. From left, the results under four conditions of depth: z = 7.5&12.5 mm, 10&15 mm, 12.5&17.5 mm, 15&20 mm. The white circles represent the true positions of the absorption change. The estimation image is surrounded by a bold frame if the positional errors of both absorbers were within one voxel (2.5 mm) in each x, y, z direction and both maximum estimation values were greater than 0.025mm⁻¹.

In the case of the horizontal distance d = 17.5 mm, which is almost the same as the probe interval, the estimation was successful until the depth condition of z = 15 mm & 20 mm as shown in Fig. 5 (d = 17.5 mm, z = 15 & 20 mm). The absorber at z = 20 mm, which corresponds to the detectable depth limit as summarized in Table 1, could be estimated accurately undisturbed by the other absorber.

In the case of the horizontal distance d = 12.5 mm, which is much shorter than the probe interval, the estimation was successful until the depth condition of z = 12.5 & 17.5 mm as shown in Fig. 5 (d = 12.5 mm, z = 12.5 & 17.5 mm). However, in the case of the horizontal distance d = 10 mm, the estimation was difficult but successful until the depth condition of z = 7.5 & 12.5 mm as shown in Fig. 5 (d = 10 mm, z = 7.5 & 12.5 mm).

The same analysis was also done for the 13-mm probe-interval condition. Table 3 summarizes the depth limits where successful estimation is possible.

Table 3 Depth limit where two absorbers at different depths are accurately estimated.

probe interval	18.4mm	13mm

horizontal distance [mm]	depth limit [mm]
15.0	15.0 & 20.0	15.0 & 20.0
12.5	12.5 & 17.5	12.5 & 17.5
10.0	7.5 & 12.5	7.5 & 12.5

The estimation results under the 13-mm probe-interval condition were the same as those under the 18.4-mm condition in spite of the higher density. The first reason would be due to the probe arrangement we use. The number of 3rd NN channels is small as mentioned in section 4.3, and the 1st NN channels of the 13-mm interval condition are sensitive in only superficial layers compared to those of the 18.4-mm interval condition. The second reason would be the nonlinear effect as mentioned in section 4.3.

5. Discussion

We proposed a three-dimensional reconstruction method for diffuse optical tomography. The main feature of the method is its high resolution in both depth and horizontal directions. The key concept that achieves this feature is the introduction of spatially variant regularization and sparsity. These are introduced in two steps. First, a provisional solution is obtained by the sensitivity-normalized Tikhonov regularization. Second, the solution is refined by the hierarchical Bayesian estimation method.

The proposed method was evaluated by applying it to the data obtained from the phantom experiment under three conditions of probe interval: 26 mm, 18.4 mm, 13 mm. This phantom experiment is consistent with our prior assumption for sparsity. Under the 26-mm interval condition, three-dimensional reconstruction failed because there was no overlapping channel due to the large distance of 2nd NN channels. In contrast, three-dimensional reconstruction was successful in the 18.4-mm and 13-mm probe-interval conditions. Under these conditions, the two absorbers whose distance is shorter than the probe interval were accurately estimated. Even when the horizontal distance was 10 mm, which is almost half the probe interval 18.4 mm, the two absorbers were accurately distinguished. Moreover, as the probe interval became shorter, deeper absorption change could be detected. The depth limit where accurate estimation was possible was around 20 mm under the 18.4-mm or 13-mm probe-interval condition as summarized in Table 1. This 20-mm depth corresponds to a depth of 24 mm with respect to the light intensity beneath the source when the medium is a head-like multilayer [4

E. Okada, M. Firbank, M. Schweiger, S. R. Arridge, M. Cope, and D. T. Delpy, “Theoretical and experimental investigation of near-infrared light propagation in a model of the adult head,” Appl. Opt. 36, 21–31 (1997). [CrossRef] [PubMed]

] consisting of a 3-mm scalp (μ_a = 0.018mm⁻¹ and μ_s′ = 0.7mm⁻¹), a 10-mm skull (μ_a = 0.013mm⁻¹ and μ_s′ = 0.9mm⁻¹), 2-mm CSF (μ_a = 0.004mm⁻¹ and μ_s′ = 0.3mm⁻¹), and 20-mm gray & white matter (μ_a = 0.019mm⁻¹ and μ_s′ = 1.1mm⁻¹).

For the sensitivity-normalized regularization, an appropriate control constant value β is required in Eq. (4). We determined the β value with reference to the sensitivity of a 22.5-mm layer. The small change of the β value does not affect the results of the hierarchical Bayesian estimation method. We confirmed that the depth limits of accurate estimation summarized in Tables 1, 2, 3 do not change even if the β value is increased by half or decreased by half. However, if β is too large, deep activities are excessively suppressed in a regularized solution. Then, the deep activities cannot be estimated by the hierarchical Bayesian estimation. If β is too small, false solutions appear in deep layers in a regularized solution because the regularization penalty becomes too small. This also leads to false solutions in the hierarchical Bayesian estimation. The appropriate value of β should be determined by considering the target region of absorption changes and the SN ratio of the measurement devices.

Finally, it should be noted that the linear Rytov approximation Eq. (1) may be violated in this phantom experiment due to the use of a strong absorber. Hence, we checked the nonlinear effect by using the Monte Carlo simulations [33

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). [CrossRef] [PubMed]

]. We compared the light intensity change calculated by linear Rytov approximation and that calculated by Monte Carlo simulation, which includes the nonlinear effect. We found the following two phenomena. First, the estimated values of the absorption change were smaller than the true value as investigated in [30

M. A. O’Leary, “Imaging with diffuse photon density waves,” Ph.D. Thesis, Unversity of Pennsylvania (1996).

,42

M. A. O’Leary, D. A. Boas, B. Chance, and A. G. Yodh, “Experimental images of heterogeneous turbid media by frequency-domain diffusing-photon tomography,” Opt. Lett. 20, 426–428 (1995). [CrossRef]

]. Second, the nonlinearity effects of the observed light intensity change depend on the positions of the source, the detector, and the absorber. These effects tended to be stronger in measurement channels where the light intensity changes were large, and stronger under the 13-mm probe-interval condition than under the 18.4-mm and 26-mm conditions. In order to take account of these nonlinear effects in the linear model, we considered nonlinear effects as noise terms that are proportional to the squared observation values. This procedure can be considered a simplified form of the second-order perturbation expansion of the nonlinear model. This was done only in the 13-mm probe-interval condition (see Appendix). The estimation results became inaccurate without these noise terms. Of course, when the linear approximation is valid, accurate estimation results are obtained without these terms, which we confirmed through a simulation study. Although we treated the nonlinear effect as noise terms, the estimation results may improve if we can apply the nonlinear effect appropriately in the forward model.

Appendices

Appendix

Noise model for nonlinear effects

Since the nonlinear effect is relatively large in the 13-mm probe-interval condition, the value proportional to the squared observation value, mean(y_j)², is added in the diagonal element of the noise covariance matrix Σ_y. This means that we considered the nonlinear effect as a noise in this study. These squared values are added only in the 13-mm probe-interval condition and not in the 18.4-mm or 26-mm interval condition.

Free energy function

The free energy function can be calculated under the factorization assumption,

F (Q) = - \frac{1}{2} \sum_{t = 1}^{T} [\bar{σ} {‖ Σ_{y}^{- \frac{1}{2}} (y_{t} - A {\bar{x}}_{t}) ‖}_{2}^{2} + {‖ {\bar{Λ}}^{\frac{1}{2}} W^{- 1} {\bar{x}}_{t} ‖}_{2}^{2}] + f (\bar{λ}, \bar{σ}, Σ_{x}),

(19)

\begin{array}{l} f (\bar{λ}, \bar{σ}, Σ_{x}) & = & \sum_{i = 1}^{N} γ_{0 i} [log ({\bar{λ}}_{i} / λ_{0 i}) - ({\bar{λ}}_{i} / λ_{0 i}) + 1] + \frac{T}{2} [log | Σ_{x} | + log | \bar{Λ} | + log | \bar{σ} Σ_{y}^{- 1} |] \\ - \frac{T}{2} [tr (Σ_{x} (\bar{σ} A^{T} Σ_{y}^{- 1} A + W^{- 1 T} \bar{Λ} W^{- 1}))] + const . \end{array}

(20)

where, x̄_t, λ̄, and σ̄ represent the expectation values of x_t, λ, and σ with respect to a trial distribution Q, respectively. Σ_x represents the covariance matrix of x with respect to a trial distribution Q.

Estimation algorithm

Hierarchical prior values and initial values are set.
In this study, we set hierarchical prior values as γ₀ = 0 and initial values as λ̄ = λ̄_init = (10 · mean(x̂_D))⁻², σ̄ = 1.
The free energy is maximized by repeating the X-step and λ-step alternately. The leading order of the computational cost is o(M²N) for each iteration step when N ≫ M. In the following equations, := denotes the assignment operator, namely, the value of the left-hand side variable is updated by the current value of the right-hand side equations.
We repeated these steps until the following two criteria were met. The change ratio of the free energy is smaller than 10⁻⁵. The iteration count reaches 1000 times.
[X-step]
In the X-step, the following values are calculated:

$Σ : = {\bar{σ}}^{- 1} + Σ_{y} + AW {\bar{Λ}}^{- 1} W^{T} A^{T},$
(21)

$h : = diag [{\bar{Λ}}^{- 1} W^{T} A^{T} Σ^{- 1} Y Y^{T} Σ^{- 1} AW + T (I_{N} - {\bar{Λ}}^{- 1} W^{T} A^{T} Σ^{- 1} AW)],$
(22)

$S : = tr [{\bar{σ}}^{- 1} Σ_{y} Σ^{- 1} Y Y^{T} Σ^{- 1} + T {\bar{Λ}}^{- 1} W^{T} A^{T} Σ^{- 1} AW] .$
(23)
The expectation and covariance of X can be calculated using these values as X̄ = WΛ̄⁻¹W^TA^TΣ⁻¹Y and $Σ_{x}^{- 1} = \bar{σ} A^{T} Σ_{y}^{- 1} A + W^{- 1 T} \bar{Λ} W^{- 1}$ , respectively.
[λ-step] In the λ-step, the expectations and shape parameters of Λ and σ are calculated as

$γ_{i} : = γ_{0 i} + \frac{T}{2},$
(24)

$γ_{i} {\bar{λ}}_{i}^{- 1} : = γ_{0 i} {\bar{λ}}_{0 i}^{- 1} + \frac{1}{2} h_{i} {\bar{λ}}_{i}^{- 1},$
(25)

$γ_{σ} : = \frac{1}{2} M T,$
(26)

$γ_{σ} {\bar{σ}}^{- 1} : = \frac{1}{2} S {\bar{σ}}^{- 1} .$
(27)
After the convergence criteria are met, X-step is performed one more step. Then, the estimation result is obtained as

${\hat{X}}_{B} = W {\bar{Λ}}^{- 1} W^{T} A^{T} Σ^{- 1} Y .$
(28)

Acknowledgments

The authors thank Makoto Fukushima for helpful comments. This research was supported by a contract with the National Institute of Information and Communications Technology entitled, ’Multimodal integration for brain imaging measurements’. A part of this research was supported by SRPBS, MEXT.

References and links

1.	A. Villringer, J. Planck, C. Hock, L. Schleinkofer, and U. Dirnagl, “Near infrared spectroscopy (NIRS): a new tool to study hemodynamic changes during activation of brain function in human adults,” Neurosci. Lett. 154, 101–104 (1993). [CrossRef] [PubMed]
2.	Y. Hoshi and M. Tamura, “Detection of dynamic changes in cerebral oxygenation coupled to neuronal function during mental work in man,” Neurosci. Lett. 150, 5–8 (1993). [CrossRef] [PubMed]
3.	A. Maki, Y. Yamashita, Y. Ito, E. Watanabe, Y. Mayanagi, and H. Koizumi, “Spatial and temporal analysis of human motor activity using noninvasive NIR topography,” Med. Phys. 22, 1997–2005 (1995). [CrossRef] [PubMed]
4.	E. Okada, M. Firbank, M. Schweiger, S. R. Arridge, M. Cope, and D. T. Delpy, “Theoretical and experimental investigation of near-infrared light propagation in a model of the adult head,” Appl. Opt. 36, 21–31 (1997). [CrossRef] [PubMed]
5.	A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50, R1–R43 (2005). [CrossRef] [PubMed]
6.	D. A. Boas and A. M. Dale, “Simulation study of magnetic resonance imaging-guided cortically constrained diffuse optical tomography of human brain function,” Appl. Opt. 44, 1957–1968 (2005). [CrossRef] [PubMed]
7.	B. W. Zeff, B. R. White, H. Dehghani, B. L. Schlaggar, and J. P. Culver, “Retinotopic mapping of adult human visual cortex with high-density diffuse optical tomography,” Proc. Natl. Acad. Sci. U.S.A. 104, 12169–12174 (2007). [CrossRef] [PubMed]
8.	C. Habermehl, S. Holtze, J. Steinbrink, S. P. Koch, H. Obrig, J. Mehnert, and C. H. Schmitz, “Somatosensory activation of two fingers can be discriminated with ultrahigh-density diffuse optical tomography,” NeuroImage 59, 3201–3211 (2011). [CrossRef] [PubMed]
9.	J. P. Culver, R. Choe, M. J. Holboke, L. Zubkov, T. Durduran, A. Slemp, V. Ntziachristos, B. Chance, and A. G. Yodh, “Three-dimensional diffuse optical tomography in the parallel plane transmission geometry: evaluation of a hybrid frequency domain/continuous wave clinical system for breast imaging,” Med. Phys. 30, 235–247 (2003). [CrossRef] [PubMed]
10.	A. Li, E. L. Miller, M. E. Kilmer, T. J. Brukilacchio, T. Chaves, J. Stott, Q. Zhang, T. Wu, M. Chorlton, R. H. Moore, D. B. Kopans, and D. A. Boas, “Tomographic optical breast imaging guided by three-dimensional mammography,” Appl. Opt. 42, 5181–5190 (2003). [CrossRef] [PubMed]
11.	H. Dehghani, M. E. Eames, P. K. Yalavarthy, S. C. Davis, S. Srinivasan, C. M. Carpenter, B. W. Pogue, and K. D. Paulsen, “Near infrared optical tomography using NIRFAST: Algorithm for numerical model and image reconstruction,” Commun. Numer. Methods Eng. 25, 711–732 (2008). [CrossRef] [PubMed]
12.	Q. Fang, J. Selb, S. A. Carp, G. Boverman, E. L. Miller, D. H. Brooks, R. H. Moore, D. B. Kopans, and D. A. Boas, “Combined optical and X-ray tomosynthesis breast imaging,” Radiology 258, 89–97 (2011). [CrossRef]
13.	F. Gao, H. Zhao, and Y Yamada, “Improvement of image quality in diffuse optical tomography by use of full time-resolved data,” Appl. Opt. 41, 778–791 (2002). [CrossRef] [PubMed]
14.	M. Guven, B. Yazici, X. Intes, and B. Chance, “Diffuse optical tomography with a priori anatomical information,” Phys. Med. Biol. 50, 2837–2858 (2005). [CrossRef] [PubMed]
15.	F. Abdelnour, C. Genovese, and T. Huppert, “Hierarchical Bayesian regularization of reconstructions for diffuse optical tomography using multiple priors,” Biomed. Opt. Express 1, 1084–1103 (2010). [CrossRef]
16.	B. W. Pogue, T. O. McBride, J. Prewitt, U. L. Österberg, and K. D. Paulsen, “Spatially variant regularization improves diffuse optical tomography,” Appl. Opt. 38, 2950–2961 (1999). [CrossRef]
17.	J. P. Culver, T. Durduran, D. Furuya, C. Cheung, J. H. Greenberg, and A. G. Yodh, “Diffuse optical tomography of cerebral blood flow, oxygenation, and metabolism in rat during focal ischemia,” J. Cereb. Blood Flow Metab. 23, 911–924 (2003). [CrossRef] [PubMed]
18.	H. Niu, F. Tian, Z. J. Lin, and H. Liu, “Development of a compensation algorithm for accurate depth localization in diffuse optical tomography,” Opt. Lett. 35, 429–431 (2010). [CrossRef] [PubMed]
19.	H. Niu, Z. J. Lin, F. Tian, S. Dhamne, and H. Liu, “Comprehensive investigation of three-dimensional diffuse optical tomography with depth compensation algorithm,” J. Biomed. Opt. 15, 046005 (2010). [CrossRef] [PubMed]
20.	N. Cao, A. Nehorai, and M. Jacobs, “Image reconstruction for diffuse optical tomography using sparsity regularization and expectation-maximization algorithm,” Opt. Express 15, 13695–13708 (2007). [CrossRef] [PubMed]
21.	P. Mohajerani, A. A. Eftekhar, J. Huang, and A. Adibi, “Optimal sparse solution for fluorescent diffuse optical tomography: theory and phantom experimental results,” Appl. Opt. 46, 1679–1685 (2007). [CrossRef] [PubMed]
22.	Y. Lu, X. Zhang, A. Douraghy, D. Stout, J. Tian, T. F. Chan, and A. F. Chatziioannou, “Source reconstruction for spectrally-resolved bioluminescence tomography with sparse a priori information,” Opt. Express 17, 8062–8080 (2009). [CrossRef] [PubMed]
23.	S. Okawa, Y. Hoshi, and Y. Yamada, “Improvement of image quality of time-domain diffuse optical tomography with lp sparsity regularization,” Biomed. Opt. Express 2, 3334–3348 (2011). [CrossRef] [PubMed]
24.	D. Wipf and S. Nagarajan, “A unified Bayesian framework for MEG/EEG source imaging,” Neuroimage 44, 947–966 (2009). [CrossRef]
25.	F. Lucka, S. Pursiainen, M. Burger, and C. H. Wolters, “Hierarchical Bayesian inference for the EEG inverse problem using realistic FE head models: depth localization and source separation for focal primary currents,” Neuroimage 61, 1364–1382 (2012). [CrossRef] [PubMed]
26.	D. Wipf and S. Nagarajan, “A new view of automatic relevance determination,” Adv. Neural Inf. Process. Syst. 20, 1625–1632 (2008).
27.	M. Sato, T. Yoshioka, S. Kajihara, K. Toyama, N. Goda, K. Doya, and M. Kawato, “Hierarchical Bayesian estimation for MEG inverse problem,” NeuroImage 23, 806–826 (2004). [CrossRef] [PubMed]
28.	T. Aihara, Y. Takeda, K. Takeda, W. Yasuda, T. Sato, Y. Otaka, T. Hanakawa, M. Honda, M. Liu, M. Kawato, M. Sato, and R. Osu, “Cortical current source estimation from electroencephalography in combination with near-infrared spectroscopy as a hierarchical prior,” NeuroImage 59, 4006–4021 (2012). [CrossRef]
29.	A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging (IEEE Press, New York, 1988).
30.	M. A. O’Leary, “Imaging with diffuse photon density waves,” Ph.D. Thesis, Unversity of Pennsylvania (1996).
31.	S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15, R41–R93 (1999). [CrossRef]
32.	S. R. Arridge and J. C. Schotland, “Optical tomography: forward and inverse problems,” Inverse Probl. 25, 123010 (2009). [CrossRef]
33.	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). [CrossRef] [PubMed]
34.	R. C. Haskell, L. O. Svaasand, T. T. Tsay, T. C. Feng, M. S. McAdams, and B. J. Tromberg, “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A 11, 2727–2741 (1994). [CrossRef]
35.	C. M. Bishop, Pattern Recognition and Machine Learning (Springer, New York, 2006).
36.	H. Akaike, “Likelihood and the Bayes procedure,” in Bayesian Statistics , J. M. Bernardo, M. H. De Groot, D. V. Lindley, and A. F. M. Smith, eds. (Univ. Press, Valencia, 1980), 143–166.
37.	J. Selb, A. M. Dale, and D. A. Boas, “Linear 3D reconstruction of time-domain diffuse optical imaging differential data: improved depth localization and lateral resolution,” Opt. Express 15, 16400–16412 (2007). [CrossRef] [PubMed]
38.	H. Dehghani, B. R. White, B. W. Zeff, A. Tizzard, and J. P. Culver, “Depth sensitivity and image reconstruction analysis of dense imaging arrays for mapping brain function with diffuse optical tomography,” Appl. Opt. 48, D137–D143 (2009). [CrossRef] [PubMed]
39.	A. C. Faul and M. E. Tipping, “Analysis of sparse Bayesian learning,” Adv. Neural Inf. Process. Syst. 14, 383–389 (2002).
40.	H. Attias, “Inferring parameters and structure of latent variable models by variational Bayes,” Proc. 15th Conf. on Uncertainty in Artificial Intelligence , Morgan Kaufmann, 21–30 (1999).
41.	M. Sato, “Online model selection based on the variational Bayes,” Neural Comput. 13, 1649–1681 (2001). [CrossRef]
42.	M. A. O’Leary, D. A. Boas, B. Chance, and A. G. Yodh, “Experimental images of heterogeneous turbid media by frequency-domain diffusing-photon tomography,” Opt. Lett. 20, 426–428 (1995). [CrossRef]

OCIS Codes
(100.3010) Image processing : Image reconstruction techniques
(100.3190) Image processing : Inverse problems
(170.3880) Medical optics and biotechnology : Medical and biological imaging

ToC Category:
Image Processing

History
Original Manuscript: April 9, 2012
Revised Manuscript: August 12, 2012
Manuscript Accepted: August 15, 2012
Published: August 21, 2012

Virtual Issues
Vol. 7, Iss. 10 Virtual Journal for Biomedical Optics

Citation
Takeaki Shimokawa, Takashi Kosaka, Okito Yamashita, Nobuo Hiroe, Takashi Amita, Yoshihiro Inoue, and Masa-aki Sato, "Hierarchical Bayesian estimation improves depth accuracy and spatial resolution of diffuse optical tomography," Opt. Express 20, 20427-20446 (2012)
http://www.opticsinfobase.org/vjbo/abstract.cfm?URI=oe-20-18-20427

Sort: Author | Year | Journal | Reset

References

A. Villringer, J. Planck, C. Hock, L. Schleinkofer, and U. Dirnagl, “Near infrared spectroscopy (NIRS): a new tool to study hemodynamic changes during activation of brain function in human adults,” Neurosci. Lett.154, 101–104 (1993). [CrossRef] [PubMed]
Y. Hoshi and M. Tamura, “Detection of dynamic changes in cerebral oxygenation coupled to neuronal function during mental work in man,” Neurosci. Lett.150, 5–8 (1993). [CrossRef] [PubMed]
A. Maki, Y. Yamashita, Y. Ito, E. Watanabe, Y. Mayanagi, and H. Koizumi, “Spatial and temporal analysis of human motor activity using noninvasive NIR topography,” Med. Phys.22, 1997–2005 (1995). [CrossRef] [PubMed]
E. Okada, M. Firbank, M. Schweiger, S. R. Arridge, M. Cope, and D. T. Delpy, “Theoretical and experimental investigation of near-infrared light propagation in a model of the adult head,” Appl. Opt.36, 21–31 (1997). [CrossRef] [PubMed]
A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol.50, R1–R43 (2005). [CrossRef] [PubMed]
D. A. Boas and A. M. Dale, “Simulation study of magnetic resonance imaging-guided cortically constrained diffuse optical tomography of human brain function,” Appl. Opt.44, 1957–1968 (2005). [CrossRef] [PubMed]
B. W. Zeff, B. R. White, H. Dehghani, B. L. Schlaggar, and J. P. Culver, “Retinotopic mapping of adult human visual cortex with high-density diffuse optical tomography,” Proc. Natl. Acad. Sci. U.S.A.104, 12169–12174 (2007). [CrossRef] [PubMed]
C. Habermehl, S. Holtze, J. Steinbrink, S. P. Koch, H. Obrig, J. Mehnert, and C. H. Schmitz, “Somatosensory activation of two fingers can be discriminated with ultrahigh-density diffuse optical tomography,” NeuroImage59, 3201–3211 (2011). [CrossRef] [PubMed]
J. P. Culver, R. Choe, M. J. Holboke, L. Zubkov, T. Durduran, A. Slemp, V. Ntziachristos, B. Chance, and A. G. Yodh, “Three-dimensional diffuse optical tomography in the parallel plane transmission geometry: evaluation of a hybrid frequency domain/continuous wave clinical system for breast imaging,” Med. Phys.30, 235–247 (2003). [CrossRef] [PubMed]
A. Li, E. L. Miller, M. E. Kilmer, T. J. Brukilacchio, T. Chaves, J. Stott, Q. Zhang, T. Wu, M. Chorlton, R. H. Moore, D. B. Kopans, and D. A. Boas, “Tomographic optical breast imaging guided by three-dimensional mammography,” Appl. Opt.42, 5181–5190 (2003). [CrossRef] [PubMed]
H. Dehghani, M. E. Eames, P. K. Yalavarthy, S. C. Davis, S. Srinivasan, C. M. Carpenter, B. W. Pogue, and K. D. Paulsen, “Near infrared optical tomography using NIRFAST: Algorithm for numerical model and image reconstruction,” Commun. Numer. Methods Eng.25, 711–732 (2008). [CrossRef] [PubMed]
Q. Fang, J. Selb, S. A. Carp, G. Boverman, E. L. Miller, D. H. Brooks, R. H. Moore, D. B. Kopans, and D. A. Boas, “Combined optical and X-ray tomosynthesis breast imaging,” Radiology258, 89–97 (2011). [CrossRef]
F. Gao, H. Zhao, and Y Yamada, “Improvement of image quality in diffuse optical tomography by use of full time-resolved data,” Appl. Opt.41, 778–791 (2002). [CrossRef] [PubMed]
M. Guven, B. Yazici, X. Intes, and B. Chance, “Diffuse optical tomography with a priori anatomical information,” Phys. Med. Biol.50, 2837–2858 (2005). [CrossRef] [PubMed]
F. Abdelnour, C. Genovese, and T. Huppert, “Hierarchical Bayesian regularization of reconstructions for diffuse optical tomography using multiple priors,” Biomed. Opt. Express1, 1084–1103 (2010). [CrossRef]
B. W. Pogue, T. O. McBride, J. Prewitt, U. L. Österberg, and K. D. Paulsen, “Spatially variant regularization improves diffuse optical tomography,” Appl. Opt.38, 2950–2961 (1999). [CrossRef]
J. P. Culver, T. Durduran, D. Furuya, C. Cheung, J. H. Greenberg, and A. G. Yodh, “Diffuse optical tomography of cerebral blood flow, oxygenation, and metabolism in rat during focal ischemia,” J. Cereb. Blood Flow Metab.23, 911–924 (2003). [CrossRef] [PubMed]
H. Niu, F. Tian, Z. J. Lin, and H. Liu, “Development of a compensation algorithm for accurate depth localization in diffuse optical tomography,” Opt. Lett.35, 429–431 (2010). [CrossRef] [PubMed]
H. Niu, Z. J. Lin, F. Tian, S. Dhamne, and H. Liu, “Comprehensive investigation of three-dimensional diffuse optical tomography with depth compensation algorithm,” J. Biomed. Opt.15, 046005 (2010). [CrossRef] [PubMed]
N. Cao, A. Nehorai, and M. Jacobs, “Image reconstruction for diffuse optical tomography using sparsity regularization and expectation-maximization algorithm,” Opt. Express15, 13695–13708 (2007). [CrossRef] [PubMed]
P. Mohajerani, A. A. Eftekhar, J. Huang, and A. Adibi, “Optimal sparse solution for fluorescent diffuse optical tomography: theory and phantom experimental results,” Appl. Opt.46, 1679–1685 (2007). [CrossRef] [PubMed]
Y. Lu, X. Zhang, A. Douraghy, D. Stout, J. Tian, T. F. Chan, and A. F. Chatziioannou, “Source reconstruction for spectrally-resolved bioluminescence tomography with sparse a priori information,” Opt. Express17, 8062–8080 (2009). [CrossRef] [PubMed]
S. Okawa, Y. Hoshi, and Y. Yamada, “Improvement of image quality of time-domain diffuse optical tomography with lp sparsity regularization,” Biomed. Opt. Express2, 3334–3348 (2011). [CrossRef] [PubMed]
D. Wipf and S. Nagarajan, “A unified Bayesian framework for MEG/EEG source imaging,” Neuroimage44, 947–966 (2009). [CrossRef]
F. Lucka, S. Pursiainen, M. Burger, and C. H. Wolters, “Hierarchical Bayesian inference for the EEG inverse problem using realistic FE head models: depth localization and source separation for focal primary currents,” Neuroimage61, 1364–1382 (2012). [CrossRef] [PubMed]
D. Wipf and S. Nagarajan, “A new view of automatic relevance determination,” Adv. Neural Inf. Process. Syst.20, 1625–1632 (2008).
M. Sato, T. Yoshioka, S. Kajihara, K. Toyama, N. Goda, K. Doya, and M. Kawato, “Hierarchical Bayesian estimation for MEG inverse problem,” NeuroImage23, 806–826 (2004). [CrossRef] [PubMed]
T. Aihara, Y. Takeda, K. Takeda, W. Yasuda, T. Sato, Y. Otaka, T. Hanakawa, M. Honda, M. Liu, M. Kawato, M. Sato, and R. Osu, “Cortical current source estimation from electroencephalography in combination with near-infrared spectroscopy as a hierarchical prior,” NeuroImage59, 4006–4021 (2012). [CrossRef]
A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging (IEEE Press, New York, 1988).
M. A. O’Leary, “Imaging with diffuse photon density waves,” Ph.D. Thesis, Unversity of Pennsylvania (1996).
S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl.15, R41–R93 (1999). [CrossRef]
S. R. Arridge and J. C. Schotland, “Optical tomography: forward and inverse problems,” Inverse Probl.25, 123010 (2009). [CrossRef]
Q. Fang and D. A. Boas, “Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units,” Opt. Express17, 20178–20190 (2009). [CrossRef] [PubMed]
R. C. Haskell, L. O. Svaasand, T. T. Tsay, T. C. Feng, M. S. McAdams, and B. J. Tromberg, “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A11, 2727–2741 (1994). [CrossRef]
C. M. Bishop, Pattern Recognition and Machine Learning (Springer, New York, 2006).
H. Akaike, “Likelihood and the Bayes procedure,” in Bayesian Statistics, J. M. Bernardo, M. H. De Groot, D. V. Lindley, and A. F. M. Smith, eds. (Univ. Press, Valencia, 1980), 143–166.
J. Selb, A. M. Dale, and D. A. Boas, “Linear 3D reconstruction of time-domain diffuse optical imaging differential data: improved depth localization and lateral resolution,” Opt. Express15, 16400–16412 (2007). [CrossRef] [PubMed]
H. Dehghani, B. R. White, B. W. Zeff, A. Tizzard, and J. P. Culver, “Depth sensitivity and image reconstruction analysis of dense imaging arrays for mapping brain function with diffuse optical tomography,” Appl. Opt.48, D137–D143 (2009). [CrossRef] [PubMed]
A. C. Faul and M. E. Tipping, “Analysis of sparse Bayesian learning,” Adv. Neural Inf. Process. Syst.14, 383–389 (2002).
H. Attias, “Inferring parameters and structure of latent variable models by variational Bayes,” Proc. 15th Conf. on Uncertainty in Artificial Intelligence, Morgan Kaufmann, 21–30 (1999).
M. Sato, “Online model selection based on the variational Bayes,” Neural Comput.13, 1649–1681 (2001). [CrossRef]
M. A. O’Leary, D. A. Boas, B. Chance, and A. G. Yodh, “Experimental images of heterogeneous turbid media by frequency-domain diffusing-photon tomography,” Opt. Lett.20, 426–428 (1995). [CrossRef]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Click here to see a list of articles that cite this paper

Figures


Fig. 1	Fig. 2	Fig. 3


Fig. 4	Fig. 5

« Previous Article | Next Article »

OSA is a member of CrossRef.

Select a Conference
Session or Paper #
Year

OpticsInfoBase

Virtual Journal for Biomedical Optics

| EXPLORING THE INTERFACE OF LIGHT AND BIOMEDICINE

Hierarchical Bayesian estimation improves depth accuracy and spatial resolution of diffuse optical tomography

Article Tools

Article Navigation

Abstract

1. Introduction

2. Methods

2.1. Forward problem

2.2. Inverse problem

2.2.1. Sensitivity-normalized Tikhonov regularization

2.2.2. Hierarchical Bayesian estimation

3. Phantom experiment

3.1. Experimental conditions

4. Results

4.1. Example of 3D reconstruction

4.2. One-absorber experiment

4.3. Two-absorber experiment 1 (at the same depth)

4.4. Two-absorber experiment 2 (at different depths)

5. Discussion

Appendices

Appendix

Noise model for nonlinear effects

Free energy function

Estimation algorithm

Acknowledgments

References and links

References

Adv. Neural Inf. Process. Syst.

Appl. Opt.

Biomed. Opt. Express

Commun. Numer. Methods Eng.

Inverse Probl.

J. Biomed. Opt.

J. Cereb. Blood Flow Metab.

J. Opt. Soc. Am. A

Med. Phys.

Neural Comput.

NeuroImage

Neurosci. Lett.

Opt. Express

Opt. Lett.

Phys. Med. Biol.

Proc. Natl. Acad. Sci. U.S.A.

Radiology

Other

Cited By

Figures