1. Introduction

Diffuse optical tomography (DOT) is a non-invasive imaging modality known for the past two decades, which uses non ionizing radiation (650-900 nm) to provide functional information about tissue hemodynamics and vascular oxygenation levels [1–3]. Specifically, it allows quantification of hemoglobin concentrations and thus estimation of brain functions or tumor angiogenesis. For functional brain imaging by DOT, reflectance geometry is a typical setup [3], where the light sources and detectors are placed on the same side of a human head to be measured. The measurement sites are carefully selected such that the diffuse reflected photons between the sources and detectors interrogate specific cortical regions that are activated by a chosen task. The measured perturbation in diffuse photon patterns, induced by the cortical activations, is used to locate or image the stimulated areas in the brain.

To achieve DOT, a variety of mathematical methods have been developed in the last 20-25 years; two recent, topical reviews provide very comprehensive descriptions in this regard [2,3]. Overall, the image reconstruction involves three major steps: First, a forward model is selected to solve light propagation distribution in tissues based on the diffusion equation or radiative transport equation. Second, the Jacobian (also called sensitivity matrix) [1,2] is generated for a chosen measurement setup with an assumption that changes in optical properties of the imaged objects are regarded as perturbations from a homogeneous “background” medium. This assumption simplifies DOT and makes it a linear inversion problem, since the light distribution under perturbation is set to be its first-order Taylor series expansion about its background distribution [3]. The last step is to inversely solve the optical properties of the imaged objects to form DOT based on all measurement readings.

Given the principle of forming DOT, there exist two technical challenges in order to obtain high quality DOT. First, since the number of measurements is much fewer than the number of pixels or voxels to be reconstructed, the inverse problems in DOT are highly ill-posed. To solve this under-determined problem in image reconstruction, a commonly used method is to introduce a regularization parameter in the inverse problem. This regularization is a theoretical key “that can be adjusted to smooth image artifacts from experimental noise and other errors at the cost of decreasing the spatial resolution” [3]. In the current field of DOT, a frequently utilized approach for regularization is based on Euclidian or L2 norm by minimizing least squares errors [4]. L2 regularization reduces high frequency noises but appears to over-smooth the reconstructed images. Thus, it is difficult to reconstruct sharp-edge profiles for DOT images with L2 regularization.

The second challenge faced by the DOT research community results from the fact that it is problematic to reconstruct accurate depth locations of an object located deeper than 2 cm due to exponential decay of light intensities with depth [1,5,6]. Several techniques have been developed to tackle this problem including (a) introduction of spatially variant regularization [7,8] besides spatially uniform regularization, and (b) usage of hard spatial priors for functional DOT in human brain imaging to enforce the reconstructed images at the correct depths or locations [1]. Recently, our group has developed a depth compensation algorithm (DCA) [6,9] to significantly improve the accuracy of depth localization in DOT. In DCA, the compensation for depth is achieved by introducing a weight matrix within L2 regularization to counter-balance the reduced sensitivity of measurement at deeper depth. After adding this weight matrix to modify the sensitivity matrix, we can still utilize the conventional reconstruction procedure. However, DCA has a limited ability to significantly improve the spatial resolution of DOT images.

Close inspection on DOT-based functional brain images reveals that most of reconstructed images appear to be in a sparse format, where most of image elements are zeros with respect to the unchanged background; only a few of them are nonzero due to changes induced by cortical activations. Given this special characteristic of DOT images, L1 regularization may be an alternative in minimizing least squares errors so as to improve the spatial resolution of DOT since L1 regularization promotes sparse solutions [10,11]. For example, two recent studies have demonstrated that L1 regularization can successfully improve the spatial resolution of reconstructed DOT images [11,12]. Then, one critical question was whether the spatial resolution and depth localization of DOT could be further improved if we combine DCA with L1 regularization.

In this study, we explored a combined approach of DCA with L1 regularization (DCA-L1) to improve the spatial resolution and depth localization of DOT; namely, reconstructed images are obtained by using L1 regularization technique after modification of the Jacobian or sensitivity matrix by DCA. Specifically, to validate the proposed DCA-L1 approach, we conducted laboratory phantom experiments using (1) a fiber-based, multichannel DOT system and (2) a camera-based DOT imaging system. Then, we compared the reconstructed DOT images using DCA-L1 and DCA with L2 regularization (DCA-L2), and estimated the performance of the two methods using the volume ratio and contrast to noise ratio. After the validation, we further demonstrated the usefulness of the DCA-L1 method, using actual experimental results taken from a human brain measurement under a finger-tapping protocol. By the end of this paper, it is concluded that L1 regularization outperforms L2 regularization for DOT image reconstruction.

2. Methods

2.1. Theory of DOT

A comprehensive review for DOT including its theoretical basis and derivation was recently published in Ref. [3]. In particular, for functional brain imaging using DOT, the theoretical derivation is clearly laid out in Refs. [1,13]. Here, we provide a brief summary that covers the theoretical foundation needed for this study.

For the forward model, the propagation of light in a highly scattering medium or biological tissue can be modeled by the diffusion approximation, which is given by [1–3]

where _{$\phi \left(r,t\right)$} is the light fluence rate given at point r and time point t. _{$S\left(r,t\right)$} is the isotropic light source, and c is the speed of light in the medium or tissue. D is the diffusion coefficient defined by the absorption coefficient,_{${\mu }_a$}, and reduced scattering coefficient, _{${\mu }_s^{\text{'}}$}, written as_{$D={\left[3\left({\mu }_a+{\mu }_s^{\text{'}}\right)\right]}^{-1}$}. Then, light propagation distribution _{$\phi \left(r,t\right)$} in Eq. (1) can be solved for a semi-infinite medium with an extrapolated zero boundary condition when either a pulsed laser [14] or continuous wave light [15,16] is utilized.

To inversely reconstruct an embedded object within tissue by DOT, namely, to identify localized heterogeneities within a homogenous medium, theoretical approaches are needed so that the measured optical quantities can be associated with the heterogeneities that perturb the homogeneous background. The commonly used forms to solve the inverse problem are either the Born or the Rytov approach, with the assumption that changes in _{${\mu }_a$} and _{${\mu }_s^{\text{'}}$}of the heterogeneities are regarded as perturbations from those of the background medium [3].

Specifically, Rytov approximation [3,13,17] assumes that the diffuse light fluence φ can be written as_{$\phi =\phi {}_{0}exp({\phi }_{scat})$}, where _{${\phi }_0$}is the background field with optical properties of the medium, and _{${\phi }_{scat}$} is the perturbed field induced by the heterogeneities. Quantity of φ can be linearly related to the spatial variations (including the heterogeneities) of the optical properties (i.e., _{$\delta {\mu }_a$}and_{$\delta {\mu }_s^{\text{'}}$}), if we assume that changes in optical properties of the imaged heterogeneities are relatively small with respect to those of the homogeneous “background”. This assumption simplifies DOT and makes it a linear inversion problem since the diffuse light fluence φ will be equal to its first-order Taylor series expansion about its background distribution, _{${\phi }_0$} [3].

In particular, for functional brain imaging, a further assumption is practically used: brain activations only perturb light absorption without changing the uniform scattering background across the medium or the brain (i.e., _{$\delta {\mu }_s^{\text{'}}$} = 0). With all given assumptions and conditions, the Rytov formulation leads to a matrix notation for functional brain imaging by DOT [1,3,13,17]:

where matrix elements of y are given by _{$y_{i\text{}}=-\log \frac{\phi (r_{s,i},r_{d,i})}{{\phi }_0(r_{s,i},r_{d,i})}$} ,with r_s,i and r_d,_i representing positions of source, s, and detector, d, respectively, for the ith measurement; matrix elements of x are _{$x_j=\delta {\mu }_{a,j}$}, which signify the perturbation in absorption, _{${\mu }_a$}, at the jth voxel within the imaging medium. Matrix A is the Jacobian or sensitivity matrix and can be derived from the photon diffusion equation, Eq. (1), using the Rytov approximation [13,17], which gives rise to

where _{${\phi }_0(r_{s,i},\text{\hspace{0.17em}}{r}_j)$}and _{${\phi }_0(r_j,\text{\hspace{0.17em}}{r}_{d,i})$} represent the light fluence rate at point r_j resulting from source, s, and contributing to detector, d, for the ith measurement, under the background medium conditions.

Overall, in Eq. (2), y represents the vector of measured optical densities, as defined by -log(φ/φ₀), between all possible pairs of sources and detectors at the measurement boundary and has the size of NM × 1, where NM is the number of measurements (or all possible combined pairs between sources and detectors). Vector x in Eq. (2) represents changes in absorption in the three dimensional (3D) image space and has the size of NV × 1, where NV is the number of total voxels in the 3D space. Matrix A has dimensions of the number of measurements by number of voxels, namely, NM × NV. Equivalently, Matrix A is often called sensitivity matrix since it reflects the measurement sensitivity to the perturbation heterogeneity within each voxel in the medium. Matrix A is also the linear transformation between the measurement space and the voxel-based image space. Based on Eq. (3), matrix A can be analytically calculated by solving the forward model of light diffusion equation, Eq. (1), with the extrapolated boundary conditions [18] for a semi-infinite medium [15], provided that optical properties of the background medium and source-detector (S-D) geometry are given.

2.2. Regularization methods for DOT

2.2.1. Depth compensation method (DCA)

It is well-known that the number of photons decreases dramatically with the increase in propagation depth, leading to the measurement sensitivity in deep tissue significantly lower than that in superficial tissue. The lower measurement sensitivity for deeper layers results in poor depth resolution and biases reconstructed images towards the superficial layers. In order to overcome this problem, weighted matrix M [6,9] was introduced and calculated, providing a pseudo-exponential increase in magnitude with depth, to compensate the sensitivity of A in deeper layers. Unlike other SVR methods [7,19] which modify the penalty term of regularization, the weighted matrix M is introduced to directly compensate the sensitivity matrix A. While the comprehensive details on DCA can be found in [6,9], we briefly review it here. The weight matrix M is formed as

where _{$M\left(A_i\right)$} represents the maximum singular values for measurement sensitivities within the particular layer i = 1,2,L, which is decomposed from the forward matrix A; γ is an adjustable power and varies between 0 and 3. Notice that the maximum singular values are arranged inversely with respect to the matrix A, namely, by the order from the bottom to surface, providing the maximum counterbalance for the deepest layer and vice versa. According to the previous studies, γ = 1.2–1.6 is considered to be appropriate for high-quality DOT images to recover embedded objects in deep tissue. In this study, a medium γ value of 1.3 was used. The adjusted sensitivity matrix A^# is defined as A^# = AM; the modified inverse problem is given by

2.2.2. Combination of DCA with L2 regularization

Similar to conventional matrix A, A^# is also under-determined and ill-posed, because the number of measurements are usually much fewer than the number of voxels to be reconstructed, as given in Eq. (5) [17]. Regularization techniques are often needed to stabilize the inversion of Eq. (5) and to overcome the ill-posed inverse problem.

In general, regularization techniques involve an addition of a second term that can be adjusted to minimize image artifacts from experimental noise by controlling a regularization parameter at the cost of reducing image spatial resolution. In DOT, the conventional form of regularization used is the L2 norm. Specifically, to solve Eq. (5), L2 least squares formulation for DOT can be given as

where λ>0 is the regularization parameter and _{${\Vert \Vert }_2^2$} denotes L2 norm. Equation (6) has an analytical solution, which can be solved directly or iteratively [2,3,13] as given by

where I is the identity matrix, S_max is the maximum eigenvalue of _{$A^{\#}{A}^{{\#}^T}$}, and α is usually set in the range of 10⁻³ to 10⁻¹ to suppress the measurement noise and stabilize the solution. While L2 norm regularization is an effective means of achieving stable solutions for the inverse problem and increasing predictive performance, it doesn’t promote sparse, sharp-edge solutions.

2.2.3. Combination of DCA with L1 regularization

On the other hand, L1 regularization promotes sparse solutions and has been reported for its uses [10,11,20]. As mentioned earlier, L1 regularization has also been studied for DOT, showing improvements in spatial resolution for sharper-edge images [11,12]. The objective function of L1-regularized least squares is given by

where _{${\Vert \Vert }_1^{}$} denotes L1 norm. In general, Eq. (8) does not have any analytical solution; the quality of the regularized solution depends on the choice of the regularization parameter, which was often selected manually. Also, the quality of reconstructed images depends on the user’s judgment. Several automatic methods, such as L-curve method, generalized cross validation method, and Morozov disperency principle, all were reported in [20,21] for this particular task.

In our work, we do not intend to develop any new methodology for L1 regularization. Instead, we applied the already developed knowledge and methodology of L1 regularization to DOT image reconstruction. Specifically, we utilized the same approach as that reported by Ref. [22] to solve the objective function with L1 regularization. While the details can be found in [22], we explain the basis of L1 regularization and how we executed it briefly as follows.

In Eq. (8), A^# is the modified sensitivity matrix or Jacobian after incorporating DCA into the objective function, y is the measurement matrix with dimension of NM × 1, and λ is the regularization parameter. Equation (8) doesn’t have an analytical solution but can be transformed into a convex quadratic form, which can be solved by standard convex optimization methods, such as interior point methods [22], as given below:

where the new variable_{$u\in R^n$} provides constraints on x. Next, adding logarithmic barrier penalties results in

As t varies from 0 to ∞, Eq. (10) converges to an optimal point. Equation (10) reaches the optimal point by utilizing Newton’s steps, NS_Newton, and searching directions by pre-conjugate gradient (PCG) method [23]. Both the number of iterations in PCG, NI_PCG, and NT_Newton have a significant impact on the reconstructed images.

2.2.4. Implementation of combined DCA with L1 regularization

Given all the needed mathematical information, our steps to implement the combined DCA-L1 algorithm are described below:

For step 4, with some modification, the L1-LS function can be expressed as x = l1_ls (A^#, y, λ, NI_PCG, NS_Newton), where x is a vector of NV × 1 to cover the 3D image volume, y is again the measurement vector containing the observed data, and λ is the regularization parameter. The reason we chose utilizing L1-LS toolbox was that it has been developed, tested, and supported by its publication [22], as well as it has the capability to handle a large set of 3D data and to have a fast computational speed.

In practice, we had to address how to determine critical empirical parameters during the regularization process: they were (1) the regularization parameter, λ (2) the value of gamma, γ, for the weight matrix M, (3) the number of iterations in PCG, NI_PCG, and (4) Newton’s steps, NS_Newton. The optimal selection of these four parameters determined the final quality of reconstructed DOT images. Based on literature, λ values of 0.1-0.01 were usually chosen, depending on experimental noise levels. Based on our own studies [6,9,26], γ values between 1.2 and 1.6 were appropriate in order to accurately recover embedded objects in deep tissue. The key issue in DCA-L1 algorithm was how to choose NI_PCG and NS_Newton. In this study, we finally selected NI_PCG and NS_Newton to be 60 and 15, based on trial and error. The ranges used to choose appropriate values of NI_PCG and NS_Newton in the trials were set 20-100 and 10-20, with an increment of 10 and 1, respectively. Specifically, the volume ratio (VR) between the reconstructed and actual objects was calculated for all the trials. A VR of unity served as a good performance criterion since VR was ideally expected to be close to “1”. In this way, the selected values of NI_PCG = 60 and NS_Newton = 15 provided us with an optimal VR in our current phantom study. Note that while running the trials to determine optimal values for NI_PCG and NS_Newton, we fixed values of λ and γ to be 0.01 and 1.3, respectively.

2.3. Tissue phantom experiments

We utilized laboratory tissue phantoms in order to assess the performance of both DCA-L1 and DCA-L2 regularizations. Two absorbers embedded inside the tissue phantom were imaged by an optical fiber-based and a CCD-camera-based imager. Volumetric image reconstruction was performed using both regularizations. Finally, the reconstructed images were compared and quantified on the basis of VR [27,28] and Contrast to Noise Ratio (CNR) [27,29,30].

2.3.1. Fiber-based DOT imager

The laboratory experiment was performed by utilizing a 32-channel, continuous-wave DOT imager [31] (DYNOT, NIRx Medical Technologies). The system delivered and collected two wavelengths of Laser at 760 nm and 830 nm, sequentially from each optical fiber. For the study, 25 bifurcated optodes were utilized and arranged as a square array of 5x5 (with a separation of 1.4 cm in both x and y direction), which was placed on the surface of the phantom (see Fig. 1 ). The data was selected from the first to sixth nearest S-D pairs (188 measurements) and used for image reconstruction. Our DOT measurement results were wavelength independent since the embedded objects were made with a low concentration of diluted black ink with a flat absorption spectrum. Thus, we utilized the data taken only from 830 nm for image reconstruction.

 

Fig. 1 (a) Experimental setup; (b) the bifurcated source-detector configuration used.

Download Full Size | PDF

A liquid tissue-mimicking phantom was prepared by filling a container of dimensions of 15x10x10 cm³ with 1% Intralipid solution. This solution served as the homogeneous background medium with an absorption coefficient (µ_a) of 0.08 cm⁻¹ and reduced scattering coefficient (_{${\mu }_s^{\text{'}}$}) of 8.8 cm⁻¹. Two spherical absorbers (µ_a = 0.3 cm⁻¹) of 1-cm diameter were placed at 3-cm depth around the center of optode array from the surface of container and separated by 3 cm, as shown in Fig. 1(a).

Volumetric image reconstruction was performed with the dimensions of x = −4 cm to 4 cm, y = −4 cm to 4 cm, and z = 0 to −5 cm. The voxel size of the reconstructed images was set to be 0.1x0.1x0.1 cm³. After reconstruction, the resultant images were sliced along both lateral cross section and depth cross section separately to show the locations of the absorbers. The dotted lines in Fig. 1(a) outline the slices of both lateral (XY plane at Z = −3) and vertical cross section (XZ plane). All reconstructed images were normalized between 0 and 1 for comparison.

2.3.2. Camera-based DOT imager

Using a fiber-based DOT system has several advantages, such as being compatible with different geometry and shape of a measured organ as well as having low noise because of direct contact of fibers on the tissue surface. However, a limited surface area on small animals is often a constraint to place many fiber optodes, and thus restricts the spatial resolution of reconstructed images. In recent years, CCD cameras have been commonly used as multichannel detectors [32,33] since they can serve as a detector array with possible thousands of virtual detectors and cover a wide field of view (FOV). Moreover, a CCD-camera-based DOT system is simpler and more portable with lower cost, as compared to a fiber-based, multichannel, DOT system.

Our CCD-camera-based DOT system consisted of an optical multiplexor (Avantes Inc. Multiplexor Channels 1x16) to deliver light at multiple locations and a 12-bit CCD camera (SamBa Q34 with Navitar Zoom 7000 lens) to serve as the detector. The FOV of the camera was set to be 13 × 11 mm². A broadband white light source (Illumination Technologies, Inc. Model 3900, quartz-tungsten halogen lamp) was connected to the multiplexor. Eight source fibers coming out from the multiplexor were arranged on the phantom surface (as shown in Fig. 2 ) to deliver the optical signals; the diffuse reflectance signals were captured by the CCD camera placed above the phantom surface. The total FOV was divided into 143 (13 columns × 11 rows) virtual detectors, while each virtual detector had 38 × 38 pixels. Total 8 (sources) x143 (virtual detectors) measurements were grouped and used to perform DOT image reconstructions.

 

Fig. 2 (a) Experimental setup; (b) light sources and CCD-camera configuration used.

Download Full Size | PDF

Similar to the fiber-based DOT experiment, an intralipid solution was used to create the liquid tissue phantom with background optical properties of μ_a = 0.1 cm⁻¹ and _{${\mu }_s^{\text{'}}$} = 10 cm⁻¹. Two spherical absorbers of 8-mm in diameter were embedded at a depth of 2 cm below the liquid surface and separated by 2.5 cm; the two absorbers had a 3:1 contrast ratio in absorption between the absorbers and background.

The volumetric images were reconstructed with dimension of x = −2 to 2 cm, y = −2 to 2 cm, and z = 0 to −3 cm. The voxel size of the reconstructed images was 0.1x0.1x0.1 cm³, being the same as that in the fiber-based imaging case. Then, we sliced reconstructed images along lateral cross section (XY-plane at Z = −2mm) and depth cross section (XZ plane) to show the locations of the reconstructed absorbers. Reconstructed images were normalized between 0 and 1 for comparison.

2.3.3. Measurement metrics

The reconstruction performances using both L1 and L2 regularizations were quantified by two measurement metrics: (1) VR [27,28] and (2) CNR [27,29,30]. Specifically, VR is the ratio of the reconstructed volume of absorber to the true volume of absorber. The volume of the reconstructed absorber was defined as the total volume of the voxels whose reconstructed µ_a values are above 50% of the maximum µ_a in the reconstructed image.

CNR indicates whether the reconstructed object can be clearly distinguished from the background. To calculate CNR, two regions, which are the volume of interest (VOI) and volume of background (VOB), were derived from the reconstructed image. VOI was defined by the location and size of the actual reconstructed object. VOB was defined by the remaining volume of the image. The CNR can be calculated by

where w_VOI and w_VOB are the weight factor of the VOI and VOB relative to the entire volume (i.e., VOI or VOB divided by the entire volume), μ_VOI and μ_VOB are the mean values of µ_a in the object and background volumes in a 3D reconstructed image, and σ_VOI and σ_VOB are the standard deviations of the two regions. In general, a high-quality reconstructed image possesses a VR value close to 1 and a high CNR value.

2.4. Applications of DCA-L1 method for human brain in vivo measurements

The phantom study provided us with an accurate and quantitative means to compare the performance between DCA-L1 and DCA-L2 and to validate that DCA-L1 outperforms DCA-L2. Then, we wished to demonstrate improved quality of DOT images by DCA-L1 using actual functional brain imaging data taken from human in vivo measurements, as an example. Specifically, we chose to image the motor cortex with our DOT while having the human subject perform a motor task (i.e., finger-tapping task) as a brain stimulation protocol. The reason to choose the motor task for assessing DOT was that this protocol has been studied intensively with either single-channel or multichannel near infrared spectroscopy (NIRS) by many research groups over the last decade [34]. According to Ref. [34], there were more than 180 papers published in this area. Thus, temporal and spatial patterns of DOT images in response to finger tapping tasks are well known and adequately published. Such knowledge could help us determine and verify whether our reconstructed DOT images were accurate with improved spatial resolution.

In the actual human measurements, we followed the protocol which was previously reported in [35] for the human brain measurements. Briefly, the subject was instructed to simultaneously tap four fingers (except thumb) up and down without moving the wrist and arm. For reliability, the subject watched a video clip of finger tapping at a frequency of ~1.5 Hz while being asked to follow the same rhythm. Specifically, an epoch of 15 seconds of tapping and 25 seconds of rest was repeated 10 times in each session. A 30-second, pre-session baseline and a 20-second, post-session baseline were also recorded. Before a formal session, the subject was instructed to perform a practice session. The experimental protocol has been approved by the Institution Review board of the University of Texas at Arlington and the University of Texas Southwestern Medical Center at Dallas.

For data acquisition from the sensorimotor cortex, a multichannel, continuous-wave NIRS system (CW-5, Techen Inc., Milford, MA) [36] was used. As shown in Fig. 3 , eight sources and sixteen detectors were used for a bilateral imaging scan. The sources were designed to emit light at 690 nm and at 830 nm, since two wavelengths were required to calculate changes in concentrations of oxy- and deoxy- hemoglobin (ΔHbO and ΔHbR) [1–3]. The source and detector optodes were arranged in such a way that they covered an area of ~8 × 5.2 cm² on each lateral side of the subject’s head and provided a total of 28 nearest S-D channels at a nearest S-D distance of 3.0 cm, as marked in Fig. 3(b).

 

Fig. 3 (a) Placement of optodes on the subject’s head; (b) S-D (Source-detector) configuration showing optode separations.

Download Full Size | PDF

It has been reported that a high-density probe array used in DOT would greatly improve DOT spatial resolution [37]. In the meantime, our recent study using laboratory phantom measurements has revealed that the quality of reconstructed DOT images depends on the measurement density asymptotically, having an optimal point for measurement density beyond which more overlapping measurements would not significantly improve the quality of reconstructed images [27]. Based on the conclusion from Ref. [27], we decided to utilize the probe geometry shown in Fig. 3 since it would provide us with a moderate spatial resolution while the setup time to place the probe array on the subject’s head with good optical contacts was reasonable (about 10-15 minutes).

The actual data acquisition rate was at 100 Hz, which was later down-sampled to 10 Hz for data processing and image reconstruction. The data was high-pass and low-pass filtered at 0.01 Hz and 0.3 Hz, respectively, to remove the baseline drift and interference due to arterial pulsations. Changes in optical density were calculated as a function of time at each wavelength. Then, the data from each S-D channel was block-averaged across time. Consequently, the block-averaged time profiles were further temporally averaged within the stimulation (i.e., finger tapping) period for each S-D channel, all of which served as inputs of y_i for Eq. (5). Reconstructed hemodynamic images of motor activation were obtained after following the steps given in Section 2.2 (more specifically in Section 2.2.4). The sensitivity matrix, A, was generated assuming μ_a (background) = 0.1 cm⁻¹ and _{${\mu }_s^{\text{'}}$}(background) = 10 cm⁻¹ for both wavelengths. In the process of DOT image reconstruction, volumetric imaging space of 20.32 × 5.84 × 4 cm³ was created with a voxel size of 0.2x0.2x0.1 mm³. The constructed images were sliced at 2.5-cm depth along the lateral cross section (i.e., in XY plane at Z = −2.5 cm) in order to compare the performance of L1 and L2 regularizations.

3. Results

3.1. Results from tissue phantom experiments

The reconstructed DOT images from the fiber-based measurement are shown in Fig. 4 . Figures 4(a) and 4(b) show the results using DCA-L1, along the vertical cross section (XZ plane at Y = 0) and lateral cross section (XY plane at Z = −3 cm), respectively. Figures 4(c) and 4(d) show the results using DCA-L2, along the same vertical (XZ plane at Y = 0) and lateral cross section (XY plane at Z = −3 cm), respectively. Figures 4(a) and 4(b) clearly reveal that L1-regularization makes the reconstructed absorbers well matched with the real absorbers in size and location. On the other hand, Figs. 4(c) and 4(d) clearly depict that the reconstructed images with L2-regularization are more blurry and diffused, as compared with the images obtained using DCA-L1. These two figures also exhibit that the size of the reconstructed absorbers are larger than their expected sizes (1cm diameter).

 

Fig. 4 Reconstructed images of two objects placed symmetrically around the center of X-Y plane and at Z = −3 cm (See Fig. 1 for geometry details). (a) and (b) were obtained with DCA-L1 and plotted in X-Z and X-Y plane, respectively; (c) and (d) were obtained with DCA-L2 and also plotted in X-Z and X-Y plane, respectively; (e) and (f) were obtained with L1 only and also plotted in X-Z and X-Y plane, respectively. The dashed circles in each panel indicate the true size and location of the absorbers to be reconstructed. The reconstructed images are normalized between 0 and 1.

Download Full Size | PDF

For comparison, we also reconstructed volumetric DOT images using L1 regularization without DCA, as shown in Figs. 4(e) and 4(f). These two figures noticeably illustrate that L1 regularization alone does not compensate the severe attenuation of measurement sensitivity with increased depth while it does greatly reduce blurry effects on the lateral plane (i.e., XY plane) of the image. This set of results are expected since L1-regularization promotes sparsity enhanced image reconstruction, but does not specifically counterbalance the sensitivity along depth.

The data collected by the CCD-camera-based imager was also analyzed to form DOT images, as shown in Fig. 5 . Figures 5(a) and 5(b) were obtained using DCA-L1, while Figs. 5(c) and 5(d) were reconstructed with DCA-L2. It is clearly seen that the two reconstructed absorbers in Figs. 5(a) and 5(b) are completely separate, while the reconstructed absorbers by L2 regularization [Figs. 5(c) and 5(d)] are more blurry and diffused. Moreover, the vertical cross-section plot in Fig. 5(c) illustrates that the two reconstructed objects generated by L2 regularization are more distorted in their shape and also pulled toward the center. This can be explained by the fact that the number of overlapping sets of measurements is relatively higher at the center than at periphery regions. Namely, the measurement sensitivities at the center are relatively higher. However, Fig. 5(a) depicts better shapes and locations of the reconstructed absorbers with respect to the true objects (dashed circles in the figure). Note that the artifacts seen in Fig. 5(c) near the superficial layers (from −0.5 cm to −1 cm) are significantly reduced in Fig. 5(a).

 

Fig. 5 Reconstructed images of two objects placed symmetrically around the center of X-Y plane and at Z = −2 cm (see Fig. 2 for the measurement geometry setup). (a) and (b) were obtained with DCA-L1 regularization, in X-Z and X-Y plane, respectively; (c) and (d) were obtained with DCA-L2 regularization and also plotted in X-Z and X-Y plane, respectively. The dashed circles in each panel indicate the true size and location of the absorbers to be reconstructed. The reconstructed images are normalized between 0 and 1.

Download Full Size | PDF

In addition, Table 1 lists VRs and CNRs between the two reconstructed absorbers with respect to the actual ones. Those ratios were calculated for both DCA-L1 and DCA-L2 methods and for each of the fiber-based and camera-based DOT imaging systems. In case of the fiber-based imager, it is clearly seen that the VRs of the two absorbers resulting from DCA-L2 are very high, 8-12 times bigger than “1”, indicating that the reconstructed absorbers are 8-10 times larger in their volume than the actual absorbers. On the contrary, DCA-L1 gives rise to the VRs of both absorbers to be close to 1 (i.e., 0.86 and 1.25). Moreover, DCA-L1 leads to a CNR to be 2-3 times better than DCA-L2. In case of camera-based DOT measurement, DCA-L1 still provides 2 times better in VRs and CNRs, as compared to DCA-L2 regularization.

Table 1. Comparison of DCA-L1 versus DCA-L2 algorithm

View Table

Overall, the results taken from tissue phantom studies with both fiber-based and camera-based DOT systems confirm that DCA-L1 surely outperforms DCA-L2, improving greatly the spatial resolution and depth localization in volumetric DOT. Next, we wish to demonstrate the usefulness of DCA-L1 using actual human brain measurement data, as an example.

3.2. Results from functional human brain imaging

The original study was reported in Ref. [35] having a total number of 8 human control subjects measured. The corresponding data across all controls were analyzed to conclude the study. For this paper, we randomly selected only one subject’s data as a representative from the 8 human subjects and performed 2D DOT images using both DCA-L1 and DCA-L2, as follows.

As expected, the subject had brain activation on the right side of the brain due to the left hand finger tapping (contra-lateral activation). Four panels in Fig. 6 show 2D slices from the reconstructed volumetric images at 2.5 cm depth from the scalp surface (i.e., XY plane at Z = −2.5 cm). The images were normalized between −1 and 1 for comparison. Figures 6(a) and 6(b) present the reconstructed ΔHbO and ΔHbR images obtained using L1 regularization; Figs. 6(c) and 6(d) show the reconstructed ΔHbO and ΔHbR images obtained using L2 regularization. In comparison, Figs. 6(a) and 6(b) show much sharper and more localized reconstructed images with L1 regularization, as compared to the blurred L2-dervied images in Figs. 6(c) and 6(d). These figures clearly demonstrate that L1 regularization can be valuably applied to functional human brain studies, and can greatly improve the spatial resolution of in vivo human brain images.

 

Fig. 6 2D slices (2.5 cm below the scalp surface) of reconstructed human brain images induced by finger tapping tasks. It shows a localized area with (a) an increase in oxy-hemoglobin concentration and (b) a decrease in deoxy-hemoglobin concentration when DCA-L1 is applied. In contrast, a larger or more diffused region is observed with (c) an increase in oxy-hemoglobin and (d) a decrease in deoxy-hemoglobin concentration when DCA-L2 is utilized for image reconstruction.

Download Full Size | PDF

4. Discussion and conclusion

In this study, we initially proposed to combine two previously published techniques, (1) DCA and (2) L1-regularization, for enhancing or improving the quality of reconstructed DOT images. Investigations on these two methods have been individually reported [6,12]: the DCA method compensates the loss of measurement sensitivity in depth and has been used with L2 regularization for volumetric DOT imaging. Such a combination, however, still cannot greatly improve the diffuse nature of DOT with over-smoothed image edges [see Figs. 4(c) and 4(d)]. On the other hand, L1 regularization has the ability to provide sharper reconstructed images and to reduce the blurry effect on the images. However, L1 regularization alone doesn’t compensate for the severe attenuation of measurement sensitivity with increased depth, as demonstrated in Figs. 4(e) and 4(f).

The key idea of this study is to make good use of both DCA and L1-regularization jointly so as to yield high-quality reconstructed images with improved depth localization and spatial resolution. We tested and validated the combined approach using laboratory tissue phantoms, with two typical scenarios of DOT, namely, a fiber-based and camera-based imaging system. The phantom studies yielded conclusive results based on two measurement metrics of VR and CNR. The VRs from both imaging systems indicate clearly that DCA-L1 algorithm outperforms DCA-L2 algorithm, at least 2-3 times better for volumetric DOT imaging. Furthermore, CNR from L1 regularization is 2 times larger than that from L2 regularization.

Throughout the study, we have also observed that the reconstructed images with the fiber-based DOT system appear to be less noisy and more accurate than the images obtained with the camera-based system. This can be mainly attributed to the bifurcated optode geometry used in the fiber-based system. If non-bifurcated fibers were used, both imaging systems would generate similar quality of DOT images regardless of the locations of absorbers [28,38].

The usefulness of the DCA-L1 approach has been also examined using actual functional brain imaging data taken from human in vivo measurements. With DCA-L1, the reconstructed human brain images from a randomly selected human subject show significant improvement in depth localization and spatial resolution of the imaged activation region/volume in the brain. Specifically, reconstructed ΔHbO and ΔHbR changes derived from DCA-L1 are more localized and concentrated in the specific or expected region [see Figs. 6(a) and 6(b)], as compared to those resulting from DCA-L2 method. In contrast, the blurry effect of L2 regularization is clearly seen in reconstructed images [see Figs. 6(c) and 6(d)]. The reason that we believe DCA-L1 produces better DOT images in response to finger tapping is given as follows.

It is well accepted by the biomedical optics community that functional brain images derived from multi-channel NIRS or DOT suffer from low-spatial resolution, as compared to fMRI, due to the scattering nature of light when it interacts with tissue. Without setting up a high cut-off threshold in ΔHbO or ΔHbR amplitude, the cortical activation regions imaged by DOT are usually bigger and less localized than those given by fMRI [9,39,40] unless a high-density DOT is utilized [37,41]. In the latter case, the reconstructed DOT images may have comparable spatial resolutions with respect to fMRI. According to Ref. [34], the motor cortex often exhibits prompt hemodynamic responses (i.e., ΔHbO) in the hemisphere contralateral to the performing limb. Consistently, we did observe contralateral activation in this study, which was approximately located near the motor cortical strip. While we did not have a quantitative measure to determine the width of motor strip in the current study, we estimated that the size of motor activation should be smaller than 2–3 cm across the motor strip (i.e., across the y-axis). This estimation was based on a recent, similar study that was performed over several human subjects using a similar DOT system, but with a high-density probe array [41]. The data in Ref. [41] allowed us to deduce that the activation size near the motor cortex evoked by a simple finger-tapping task was about ~2 × 3 cm. The reconstructed images in Fig. 6 clearly illustrate that DCA-L1 leads to the reconstructed brain activation size more localized and accurate than DCA-L2, with respect to the expected activation region. Note that in this human brain study, the optimized empirical parameters (i.e., NI_PCG and NS_Newton) were derived from our phantom results. For quantitative or rigorous validation of DCA-L1, we plan to conduct a joint fMRI-DOT study in order to make volumetric DOT possible for human brain imaging.

In conclusion, based on tissue-phantom studies, we have validated that the combination of DCA with L1-regularization can offer significant improvement in depth localization and spatial resolution for DOT images. We have further demonstrated the applicability and usefulness of this method using in vivo measurements of functional human brain imaging. In general, this DCA-L1 approach can be extended to other applications of DOT, such as breast and prostate cancer detection, and can be further explored to improve quantification of tumor optical properties because of more localized targets identified.