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–8] and breast cancer detection [9–12]. 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] 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], 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,15].

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]. Culver et al. proposed similar spatially variant regularization in which the regularization weight is calculated from the sensitivity matrix [17]. 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, 19]. 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–23]. 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, 25]. 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]. Since it has advantages in expandability, it has recently been applied in multi-modality imaging, such as the combination of MEG and fMRI [27] and that of EEG and fNIRS [28]. 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 ${\Vert \mathit{x}\Vert }_2^2$. The mean value $\frac{1}{T}{\sum }_{t=1}^T{\mathit{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–32]:

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], 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]:

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

Sensitivity-normalized regularization works strongly in superficial layers and weakly in deep layers, since the sensitivity exponentially decays with depth [16]. This weight matrix is similar to that of Culver et al. [17, 37, 38] 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.5mm, where I_≥22.5mm 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:

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, 39]. 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]. 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

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:

The hyperparameter λ is estimated by introducing a hierarchical prior:

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

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

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

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, 39].

In this study, the reliability value of the hierarchical prior was set as γ₀ = (γ₀₁, ⋯ , γ_0N)^T = 0. Thus, the prior of Eq. (11) becomes non-informative, and the value λ_0i does not affect the probabilistic model and the iterative algorithm. In this case, z becomes sparse [35, 39], 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 γ_0i > 0, z does not converge to zeros because of the prior information. However, if γ_0i 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, 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.

Download Full Size | PDF

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\hspace{0.17em}\text{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.

Download Full Size | PDF

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.

Download Full Size | PDF

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.

View Table | View all tables in this article

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

Download Full Size | PDF

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.

View Table | View all tables in this article

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

Download Full Size | PDF

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.

View Table | View all tables in this article

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] 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]. 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,42]. 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.