Bundled-optode implementation for 3D imaging in functional near-infrared spectroscopy

The paper presents a functional near-infrared spectroscopy (fNIRS)-based bundled-optode method for detection of the changes of oxy-hemoglobin (HbO) and deoxyhemoglobin (HbR) concentrations. fNIRS with 32 optodes is utilized to measure five healthy male subjects’ brain-hemodynamic responses to arithmetic tasks. Specifically, the coordinates of 256 voxels in the three-dimensional (3D) volume are computed according to the known probe geometry. The mean path length factor in the Beer-Lambert equation is estimated as a function of the emitter-detector distance, which is utilized for computation of the absorption coefficient. The mean values of HbO and HbR obtained from the absorption coefficient are then applied for construction of a 3D fNIRS image. Our results show that the proposed method, as compared with the conventional approach, can detect brain activity with higher spatial resolution. This method can be extended for 3D fNIRS imaging in real-time applications. © 2016 Optical Society of America OCIS codes: (170.2655) Functional monitoring and imaging; (300.0300) Spectroscopy; (110.0110) Imaging systems; (100.2960) Image analysis. References and links 1. A. Villringer and B. Chance, “Non-invasive optical spectroscopy and imaging of human brain function,” Trends Neurosci. 20(10), 435–442 (1997). 2. L. Gagnon, R. J. Cooper, M. A. Yücel, K. L. Perdue, D. N. Greve, and D. A. Boas, “Short separation channel location impacts the performance of short channel regression in NIRS,” Neuroimage 59(3), 2518–2528 (2012). 3. M. A. Kamran and K.-S. Hong, “Linear parameter-varying model and adaptive filtering technique for detecting neuronal activities: an fNIRS study,” J. Neural Eng. 10(5), 056002 (2013). 4. M. Aqil, K.-S. Hong, M.-Y. Jeong, and S. S. Ge, “Detection of event-related hemodynamic response to neuroactivation by dynamic modeling of brain activity,” Neuroimage 63(1), 553–568 (2012). 5. D. T. Delpy, M. Cope, P. van der Zee, S. Arridge, S. Wray, and J. Wyatt, “Estimation of optical pathlength through tissue from direct time of flight measurement,” Phys. Med. Biol. 33(12), 1433–1442 (1988). 6. H. Santosa, M. J. Hong, S.-P. Kim, and K.-S. Hong, “Noise reduction in functional near-infrared spectroscopy signals by independent component analysis,” Rev. Sci. Instrum. 84(7), 073106 (2013). 7. M. R. Bhutta, K.-S. Hong, B.-M. Kim, M. J. Hong, Y.-H. Kim, and S.-H. Lee, “Note: three wavelengths nearinfrared spectroscopy system for compensating the light absorbance by water,” Rev. Sci. Instrum. 85(2), 026111 (2014). 8. M. J. Khan, M. J. Hong, and K.-S. Hong, “Decoding of four movement directions using hybrid NIRS-EEG brain-computer interface,” Front. Hum. Neurosci. 8, 244 (2014). 9. N. Naseer, M. J. Hong, and K.-S. Hong, “Online binary decision decoding using functional near-infrared spectroscopy for the development of brain-computer interface,” Exp. Brain Res. 232(2), 555–564 (2014). 10. K.-S. Hong, N. Naseer, and Y.-H. Kim, “Classification of prefrontal and motor cortex signals for three-class fNIRS-BCI,” Neurosci. Lett. 587, 87–92 (2015). 11. J. G. Kim and H. Liu, “Variation of haemoglobin extinction coefficients can cause errors in the determination of haemoglobin concentration measured by near-infrared spectroscopy,” Phys. Med. Biol. 52(20), 6295–6322 (2007). 12. X. S. Hu, K.-S. Hong, and S. S. Ge, “fNIRS-based online deception decoding,” J. Neural Eng. 9(2), 026012 (2012). Vol. 7, No. 9 | 1 Sep 2016 | BIOMEDICAL OPTICS EXPRESS 3491


Introduction
The paper proposes a method for construction of a three-dimensional (3D) brain image using a bundled-optode arrangement in functional near-infrared spectroscopy (fNIRS).The system's multiple channels with different emitter-detector distances are configured for acquisition of the hemodynamic responses (HRs) in a local brain region.The HRs of 256 voxels are utilized to construct a 3D brain image that shows the corresponding brain activity in detail.The proposed method is validated by arithmetic-task experiments.
Brain activity can be measured using non-invasive techniques such as functional magnetic resonance imaging (fMRI), electroencephalography (EEG), and fNIRS.fMRI has the advantage of high spatial resolution, but suffers from low temporal resolution, whereas EEG, contrastingly, has high temporal but low spatial resolution.Meanwhile, fNIRS with acceptable temporal and spatial resolutions can be used for recording of hemoglobin concentration changes [1].Moreover, it is an inexpensive, portable, and noninvasive optical technique [2][3][4].For detection of brain activity, fNIRS normally utilizes multi-wavelength light within the 700-900 nm range that travels through tissues.Portions of the light measured by the detectors from the subject's scalp are assumed to have travelled along the bananashaped path in the cortex.The intensities of the detected light are then converted to the concentration changes of oxy-hemoglobin (HbO) and deoxy-hemoglobin (HbR) using the modified Beer-Lambert law (MBLL) [5][6][7][8][9][10][11][12][13][14][15].Recently, fNIRS has demonstrated the possibility of detecting an initial dip (i.e., a decrease of HbO caused by neuronal firing) [16], which allows an early detection of brain activity.Additionally, several mathematical models were developed for further improvement of the measured HRs [4,17,18] and the spike-firing rate of neurons [19].
Commercial fNIRS equipment now has temporal resolution as high as 250 Hz [20] and spatial resolution of approximately 10 mm [21].For improved spatial resolution, the density of optodes needs to be higher [20,22].The previous work [23] revealed that the multipledistance method can be utilized for an effective removal of motion artifacts.Additionally, overlapping of banana-shaped light paths has demonstrated the improvement of spatial resolution and localization accuracy [20,24].Also, the larger the emitter-detector pair's distance is, the higher the depth sensitivity is [20,25].In our current work, an improved spatial resolution has been achieved with a bundled-optode method (two 16-optode bundles) forming multiple channels with different distances.
In continuous-wave fNIRS, the differential path-length factor (DPF) in the MBLL equation, an unmeasurable optical property, is unknown.In a scattering medium, however, this factor can be approximated as a function of the emitter-detector distance [1,11,26,27]; and at each wavelength, the DPF is considered as a constant.However, it has been shown that a constant DPF can incur errors or otherwise degrades the accuracy of chromophore concentration changes.More seriously, a DPF modulation can lead to an opposite hemoglobin concentration changes [28].To overcome the above problems, in this current work, rather than using the MBLL law with a constant DPF (i.e., the conventional approach), the path lengths (PL) are estimated for specific brain layers.
In fNIRS, optical brain images are constructed based on an inverse model (continuous wave type) or a forward model (frequency domain and time domain types).The inverse model is employed in the case that the optical properties cannot be measured.In this model, the optical image is constructed, instead, from the measured light intensities, y = A × x (where A refers to the partial sensitivity matrix and x denotes the optical properties of the medium) [24,[29][30][31][32][33]: The definition of partial sensitivity is found in [34,35].The optical properties can then be computed by solving the inversion problem using Tikhonov regularization [29][30][31].Use of this method requires that the partial sensitivity be known.However, partial sensitivity measurement for an actual human brain has yet to be achieved.Most previous works have used Monte Carlo simulation to estimate the partial sensitivity of a human-brain model [34,[36][37][38].For a more precise construction of optical brain imaging, the depth compensation algorithm has been employed in [27,[29][30][31].In their work, the sensitivity matrix was modified by adjusting the depth of each layer in the human-brain model.On this basis, the depth compensation algorithm, as combined with adaptive filtering, can improve constructed brain images [27].Meanwhile, in the case that the optical properties can be measured, a forward model is employed.Based on those measured properties, optical brain images are constructed [39].In our work, to simplify brain-image construction, the data set in [34] regarding the partial sensitivity of the layered tissues are used in the estimation of PLs as a function of distance.
For recording the brain-hemodynamic responses to arithmetic tasks, a bundled-optode method is developed.The absorption coefficients are computed from the intensities of the detected light and the changes of the estimated PLs with two-layer brain tissues (i.e., gray matter/GM and white matter/WM).The HbOs obtained from the absorption coefficients are used to construct 3D fNIRS imaging.Our experimental results demonstrated that the proposed method can construct a 3D fNIRS image with an improved spatial resolution.Moreover, our proposed method can detect more precise locations than those obtained by the conventional approach.

Bundled-optode method for a two-layer brain model
According to the continuous-wave fNIRS technique, the incident light emitted into tissues (i.e., scalp, skull, cerebrospinal fluid (CSF), GM, and WM) [40] and the detected light reflected back to the scalp can be measured by a detector placed a few centimeters apart from the given emitter.This phenomenon is described by the Beer-Lambert equation [5,23,[41][42][43] a ( , ) ( ) out in ( , ) ( , ) , where where the superscript i indicates the channel number, in ( , ) i I t λ and out ( , ) i I t λ stand for the incident and detected light intensities at the i-th channel, respectively, and L i denotes the mean path length of the light traveled from the emitter to the detector at the i-th channel.Using Eq. ( 2), the absorbance difference at two time instances, t -1 and t, can be computed as follows.
Actually, the mean path length of the light in a real human brain is yet unavailable.Most researchers have used Monte Carlo simulation to estimate this value using a human brain model [34,36,37].From Eq. (3), the mean path length can be interpreted as the partial sensitivity defined as follows [34,35].a ( , ) .( , ) Provided that the PLs of the five brain layers (i.e., scalp, skull, CSF, GM, and WM) are estimated, Eq. (3) can be rewritten as follows [27,44,45].
The early works [46,47] demonstrated a lower perfusion rate and a longer retention rate in the skin and skull than in the deeper layers.Therefore, assuming that only the absorption in the brain is changed by the brain activation (i.e., Eq. ( 6) can be rewritten as ( ) .
For bundled optodes with N channels, Eq. ( 7) is represented as Solving Eq. ( 8) for two wavelengths ( a 1 ( , ), the following two equations are obtained.
Finally, the changes in hemoglobin concentration at the i-th channel are computed as follows.

Three-dimensional fNIRS imaging
Figure 1(a) illustrates the concept of a bundled-optode configuration for the detection of a brain activity in a local brain region.In this arrangement, the circles and squares indicate emitters and detectors, respectively.Also the (blue) arrows illustrate the direction of light from emitters to detectors.Figure 1(c) depicts the measurement points (in color) in the 3D space for the bundled-optode configuration inside the (blue) dashed rectangle in Fig. 1(a).It is noted that the italic numbers denote the optodes' positions.The main objective of this work is to improve the spatial resolution of the fNIRS technique using densely spaced optodes.To that end, the proposed method is evaluated by plotting the brain activity across all the channels in the 3D space.
Figure 1(b) provides an example showing one emitter and four detectors together with the banana-shaped light paths traveling through the brain tissues at different depths.For the bundled-optode approach, multiple channels formed with different distances can detect a pattern(s) at different depths.Naturally, if more optodes are used, the patterns can be more precisely detected.In this paper, we focus on the analysis to the bundled-optode configuration shown in the blue dashed rectangle in Fig. 1(a).
Let (x, y, z) be the coordinate of an optode in the 3D space.Each emitter in the left-hand side (a circle) can be combined with 16 detectors in the right-hand side (squares).Figure 1(c) illustrates a combination of one emitter and 4 detectors in parallel (the illustration of 16 channels is too messy): For example, in the 2nd row, emitter 6 and detectors 21 -24 form the four channel configuration indicated by the small yellow squares.Optode 13 is chosen as the origin, and the coordinates of the other optodes are defined accordingly.In this work, the coordinate of each channel is computed based on its specific emitter-detector pair.The location of a specific channel is defined as the midpoint of the emitter-detector pair [26].For instance, the positions of optode 6 (emitter) and optode 21 (detector) are (x 1 , y 1 , z 1 ) and (x 2 , y 2 , z 2 ), respectively.Then the position of this channel is computed as (x, y, z) = ((x 1 + x 2 )/2, (y 1 + y 2 )/2, sqrt((x 1 -x 2 ) 2 + (y 1 -y 2 ) 2 )/2) [26].In this case, the z coordinate denotes the depth in which the brain activity has occurred.For our bundled-optode arrangement with 256 channels, it can fully describe the 3D brain image in a local region.
The brain activity of each channel is described by the mean values of HbX (i.e., HbO and HbR) encoded in color.The t-and p-values are calculated for individual trials using the robustfit function available in Matlab (the Mathworks, Inc.) [48].With this function, the averaged HbO in each channel is compared with the expected HR model (or the designed HR) by the least squares method, in which the designed HR has been produced by convolving one trial period (i.e., 10 s task and 20 s rest) with the canonical HR function (i.e., the difference of two gamma functions) [18].If t-and p-values satisfy the threshold criteria (i.e., t > 0 and p < 0.05), the HbXs are included to compute the (trial) average over the 3~11 s window.The averaged HbXs are standardized into the 0~1 range and then encoded in color before plotting them at their 3D coordinates using the scatter3 function (3D scatter plot) available in Matlab.
In order to estimate the partial sensitivity, we use the data set obtained in the work of [34].A fit function available in Matlab is used to fit the data on the GM and WM layers with a chosen function.We used the linear polynomial curve fitting function in Matlab to estimate the partial sensitivities of the GM and WM from the data in [34], which resulted in the following two partial sensitivities of GM and WM.
where d is the emitter-detector distance.
The advantage of a bundled-optode configuration is that, with its multiple channels, it can depict the brain activities at different depths (see Figs. 1(b), 1(c)).Therefore, the proposed method can construct brain images with an improved spatial resolution.Since our method can compute the coordinates of each channel in the 3D space (x, y, z) based on the known probe geometry, the active locations in the measured brain region can be clearly visualized.

Subjects
Five healthy voluntary male students (mean age: 31.6 years) were asked to perform a series of mental computation tasks.None of the subjects had any neurological impairment or mental disorder.Prior to experiments, the procedures were clearly explained to the subjects.In order to eliminate any interference from noises, the experiments were conducted in a dark and quiet room.In addition, the subjects were asked to comfortably sit on a chair and open their eyes during the experiment.Also they were asked not to move their body throughout the entire duration of experiment.
Figure 2 illustrates the experimental paradigm with five trials.Each trial consists of 10 s task and 20 s rest periods.After a 20 s initial rest preparatory to the first trial, each subject is asked to answer a sequence of problems randomly appearing on a 15-inch laptop screen over the span of 10 s.Each problem consists of four operations (i.e., addition and subtraction of two digit numbers, multiplication, and division) presented in a pseudo random order.For instance, (45 -39) × 2 = X, 52 + X = Y, 75 -Y = Z, etc. within the 10 s period, where the bold numbers are the calculated results from the previous questions [10].The total duration of one experiment is 170 s.All the experiments in this study were performed by following the guideline of the Institutional Review Board, Pusan National University, by obtaining the consents from the subjects regarding the experimentation, and were conducted in accordance with the ethical standards encoded in the latest Declaration of Helsinki.For this work, we planned to measure the brain activities in various local brain regions using a bundled-optode arrangement.To strengthen the obtained results, each subject was scanned two times with the same experimental protocol of the arithmetic task shown in Fig. 2. Two different setups of the optode configuration for detecting the active brain regions are described in the following sections.

First scan: configuration covering a wide brain region
In this scan, the optode configuration in Fig. 3 with twelve channels covering a wide brain area in the prefrontal cortex was used to identify a local brain region in the prefrontal cortex.For each subject, the mean values of twelve channels are used to plot a two-dimensional brain map using the griddata function (interpolation of scattered data) available in Matlab.The mean value of the averaged HbO of each channel is computed over the 3~11 s window.It is noted that the active brain regions are composed of the channels with high mean values.

Second scan: configuration for scanning a local brain region
The second scan uses the bundled-optode arrangement with 256 channels illustrated in Fig. 1(c) and Fig. 4 to evaluate the local brain activity in detail.For the same subject participated in the first scan, the optode arrangement in Fig. 4 has been put covering the highly active channels found in the first scan using the optode configuration in Fig. 3.The 3D brain map is then plotted based on the computed mean values of 256 channels (see Subsection 2.2 for more details).According to our bundled-optode configuration with different distances, the brain activities at different depths could be visualized.It means that a pattern can be generated.

Equipment and data conversion
A dual-wavelength (760 and 830 nm) continuous-wave fNIRS instrument with 32 optodes (DYNOT; NIRx Medical Technologies, USA) was utilized to measure the brainhemodynamic responses under constant light intensities.The system has 20mW/wavelength of power.For the fNIRS technique, the previous work [49] revealed that the relative distance for an emitter-detector pair is normally arranged 2.0~5.5 cm to detect the brain activities.Therefore, in our bundled-optode arrangement (see Fig. 1(c) and Fig. 4), the distances have been configured within the 2.0~5.0 cm range.Since an emitter in the left-hand side (circles) can be combined with 16 detectors in the right-hand side (white squares) producing 16 channels, a total of 256 channels in the prefrontal cortex, operating at a sampling rate of 1.81 Hz, can be configured (see Fig. 4).
The measured fNIRS data were converted to hemoglobin concentration changes, HbX, using Eq. ( 10) and our own Matlab codes.The HbX signals were filtered by a low-pass filter with a cut-off frequency of 0.15 Hz to remove the respiration and pulse-related physiological noises.Also, since fNIRS data exhibit some baseline drifts [48,50], a baseline-correction method has been utilized to remove such baseline drifts in the measured fNIRS data: In our work, a 4th order polynomial has been fit to the measured data, then the obtained 4th order polynomial curve has been subtracted from the measured data, which yields the corrected values [51].According to our measured data, in most cases, HbOs increase whereas HbRs decrease during the activation tasks.However, in some channels, HbRs also increase even during the task interval.In relation to this issue, some previous works pointed out that HbO is more sensitive and more reliable than HbR [18,52,53].Therefore, in the present work, only HbOs were used for the analysis.

Regions of interest
The regions of interest (ROIs) are the brain areas in which the t-values are higher than the critical t-value (t crt ) [50].In the present work, t crt is set to 1.674 according to the degree of freedom (M -1 = 53, M is number of data samples of one trial) and the level of statistical significance (i.e., α = 0.05 for one-tailed test).The t-value of each channel is computed by comparing the HbO averaged across five trials with the expected HR using the robustfit function available in Matlab (see Subsection 2.2 for more details).If the HbO response (or the shape) is closer to that of the expected HR, its t-value becomes high.If the response of the signal is different, the t-value becomes low (or opposite).If a channel with the computed tvalue is greater than t crt , the channel is considered active.

Results
In this paper, we proposed a bundled-optode approach for the detection of HbX and construction of the brain activity in the 3D space.The HbX were computed using the Beer-Lambert equation with the estimated PLs.The mean values of HbOs at each channel's computed location in the 3D space were plotted with encoded colors.A total of 256 channels/voxels were used for the construction of a 3D fNIRS image.
In the case of a scattering medium, the PL was assumed as a function of emitter-detector distance [1].In the conventional method, the differential PL is considered as a constant for each wavelength; in our proposed method though, the PL is estimated as a function of the emitter-detector distance.This means that with changes in the distance of an emitter-detector pair, the PL varies accordingly (see Fig. 5).The partial sensitivity (see Eq. ( 5)), L, was estimated as a function of the emitter-detector distance based on the data set obtained in [34].The emitter-detector distances of all of the channels used in this work were arranged from 2.0 to 5.0 cm.
Figure 5 plots the partial sensitivities of GM (blue solid line) and WM (red dashed line).It is noteworthy that both partial sensitivities are affine functions.Moreover, the magnitude of the partial sensitivity of WM is significantly smaller than that of GM.For the 256 channels, their absorption coefficients, µ a (λ, t), were computed from the measured optical densities and the partial sensitivities of GM and WM using Eq. ( 9), and were converted to HbX using Eq.(10).The highly active channels are first identified by using the optode configuration in Fig. 3 covering a wide brain region in the prefrontal cortex.Figure 6 plots the two-dimensional brain activation map of Subject 3 (a representative subject).In this case, all 16 channels were plotted at the same depth (i.e., ~1 cm) because their emitter-detector distances were all 2.2 cm (see Fig. 3).As can be seen, active brain regions were found in Channels 10, 12, 14, and 16.According to Brodmann area (BA), the brain regions of these active channels are BA9 and BA46 [54].highly active brain region of Subjects 2, 3, and 5 was the right dorsolateral prefrontal cortex (DLPFC) whereas that of Subject 1 and 4 was the left DLPFC [6].The bundled-optode arrangement in Fig. 4 was then put at the active brain region (i.e., the black dashed rectangle in Fig. 6 including active channels 10, 12, 14, and 16) for the second scan.The spatial resolution of the optode configuration in Fig. 3 is low (i.e., 16 channels).However, the spatial resolution of the bundled-optode arrangement in Fig. 4 is actually high (i.e., 256 channels).Moreover, the different distances of emitter-detector pairs can make the patterns precisely detected.Therefore, the bundled-optode method gives a precise local brain activity in detail.For this case, the brain activation map can be plotted in the 3D space as in Section 2.2.
To evaluate our proposed method, the optical light intensities measured from the arithmetic-task experiment with 5 subjects were applied to the computation of µ a (λ, t) and HbX using Eq. ( 9) and Eq. ( 10).The PLs of the two layers, GM and WM, in each channel were computed using Eq. ( 11) and Eq. ( 12). Figure 7 illustrates a representative active channel of Subject 3, HbO (Fig. 7(a)) and HbR (Fig. 7(c)).The blue and red thick curves are the unfiltered and filtered HbX, respectively.The gray bars indicate the stimuli of the arithmetic tasks.Figures 7(b) and 7(d) show the averaged HbO and HbR over the five trials in Figs.7(a) and 7(c), respectively.As can be seen, after each stimulus the HbO signal increases whereas the HbR signal decreases.The HbO signal decreases when the stimulus is over.The trend was consistent throughout all five trials (see Fig. 7).It should be noted that the HbOs in this case were computed using the conventional method with constant DPFs of 7.15 and 5.98 at 760 nm and 830 nm, respectively [6,50].Figures 8(b), 8(d), and 8(f) show the 3D brain activity constructed according to the HbOs as computed by our proposed method with Eq. ( 9) and Eq.(10).Figures 8(a As can be seen, the brain activity using our proposed method is distinctive whereas that using the conventional method appears more dispersive.For example, the active locations (i.e., red areas) in Fig. 8(a) (the conventional method) are larger than those in Fig. 8(b) (the proposed method).For the left and right views, the active regions at the bottom in Figs.8(c) and 8(e) are more dispersive than those in Figs.8(d) and 8(f), respectively.
The active locations detected by the proposed method, furthermore, are more diverse than those using the conventional approach (see Figs. 8(c)-8(f)).For instance, the active regions in the proposed method are identified at the top side (see Figs. 8(d) and 8(f), whereas they are not detected by the conventional one (see Figs. 8(c) and 8(e)).Also, this result is consistent across Subjects 1 and 2 (see Fig. 9) and Subject 5 (see Fig. 10).Table 1 shows a comparison of the number of active channels in ROIs between the conventional and the proposed methods.The numbers of active channels detected by the proposed method are higher than those by using the conventional approach.The result demonstrated that our proposed method can improve the detection of active channels.The criterion for identifying active channels has been described in Subsection 3.5.

Discussion
The main objective of this work is to develop a method for presenting a 3D brain image with an improved spatial resolution using a bundled-optode arrangement.We assumed that the estimated PL of two-deeper layers (i.e., GM and WM) in Eq. ( 9) can improve the reconstructed 3D brain image.The proposed method was compared with the conventional MBLL approach to evaluate our hypothesis.For the conventional method, the DPF is considered as constant at each wavelength.However, using a constant DPF can negatively affect the accuracy of chromophore concentration changes [28].In our proposed method, the PL is computed as a function of the emitter-detector distance in two-layer tissues.In this way, the proposed method can describe brain activity more precisely than the conventional approach.Indeed, our results revealed that the proposed method indicates brain activities with improved spatial resolution (see Figs. 8, 9, and 10).Additionally, the proposed method can detect more active brain regions than those using the conventional approach (see Figs.  1).Also, our results demonstrated that, for the same Subject 3, the active locations shown in Fig. 6 are too dispersive (14 optodes were utilized to form 16 channels) whereas the active regions in Fig. 8 are more distinctive in both the conventional method and the proposed method (32 optodes are utilized to form 256 channels).It means that if the greater the density of optodes is, the more improved the spatial resolution is.This issue was argued in some relevant previous works [20][21][22].Especially, Habermehl et al. [21] revealed that an increase of the number of channels can lead a better lateral resolution.In their work, vibro-tactile stimulation was applied to the little and thumb fingers of eight experimental participants.They found that, in five out of eight participants, the active locations of two fingers in the somatosensory cortex could be distinguished.In our work, the bundled-optode arrangement with 256 channels was used to scan the local brain region.Their activity was plotted in detail in the 3D space.In general, our proposed method can show the brain activity more distinctive and detect more diverse brain regions than those using the conventional approach.
The proposed bundled-optode method has certain advantages.Due to the system's multiple-channel configuration (i.e., 256 channels) with different distances arranged in 2.0~5.0 cm, the depth of 3D brain imaging varies in the 1.0~2.5 cm range.Therefore, the proposed method can provide more precise information on the brain activity in a local brain region.Moreover, because the optodes in the proposed method were put side by side, in theory, the channels' banana-shaped light paths overlap, thus enabling highly precise detection of the active locations of a local brain region.This issue has been contested in some of the relevant previous works [20,24].Additionally, the positions of the channels in the 3D plane are computed based on the known probe geometry.It should be noted, however, that with our configuration, only a relatively small brain region (i.e., 5.0 cm × 1.5 cm) could be investigated, due to the equipment's limited number of optodes (32 optodes in total).Another issue is the distance between two optodes in our bundled-optode configuration, 0.5 cm.Therefore, limited number of long (2.0 cm to 5.0 cm) and short (0.5 cm) separation channels were formed.Some previous works revealed that short separation channels could be utilized for removal of any superficial noises [2,[55][56][57].However, we did not pursue this in our current work.Certainly in future study, short separation channels (measuring the extra-cortical layers) will be applied to improve the quality of measured brain signals.

Conclusions
In this paper, the bundled-optode method for construction of a 3D fNIRS brain image was presented.The mean values of HbOs at the 256 voxels were used to project the brain activation in the 3D space.The method was validated by experiments entailing the performance of arithmetic tasks.Our results showed that the proposed method can depict more active locations than those in the conventional approach.Also, the active regions detected by the proposed method are more distinctive than those using the conventional one.Additionally, our proposed method could be used for monitoring 3D brain imaging in realtime.For our future work, a bundled-optode configuration applied to a 5 cm × 5 cm brain region (see Fig. 1(a)) will be investigated.We believe that with this configuration, the spatial resolution of fNIRS brain imaging will be significantly further improved.

Fig. 1 .Fig. 2 .
Fig. 1.(a) Concept of bundled-optode configuration, (b) light traveling through tissues in banana-shaped light path, and (c) locations of channels in the 3D space (x, y, z) computed from a pair of an emitter (circle) and a detector (white square): Italic numbers indicate optodepositions.

Fig. 3 .
Fig. 3.An optode arrangement for identifying a local brain region activated from arithmetic tasks: Fp1 and Fp2 are two reference points in the international 10-20 system, circles and squares denote emitters and detectors, respectively, and numbers in red indicate channels.

Fig. 4 .
Fig. 4. A bundled-optode arrangement for acquisition of HRs in a local brain region.

Fig. 5 .
Fig. 5.The partial sensitivity (L) of GM (blue solid line) and WM (red dashed line) estimated as a function of distance.

Fig. 6 .
Fig. 6.Active brain regions of representative Subject 3 identified using the optode arrangement in Fig. 3: Numbers in white indicate the channels and the color bar on the right-hand side denotes the activation strength according to the mean values of HbOs.

Figure 8
Figure 8 plots the brain activation maps of representative Subject 3 in the 3D space.The small circles mark the optodes' positions (see Fig. 1(c) and Fig. 4 for their arrangements).Figures 8(a), 8(c), and 8(e) show the brain activity reconstructed according to the mean values of HbOs (from the 256 channels).It should be noted that the HbOs in this case were computed using the conventional method with constant DPFs of 7.15 and 5.98 at 760 nm and 830 nm, respectively[6,50].Figures8(b), 8(d), and 8(f) show the 3D brain activity constructed according to the HbOs as computed by our proposed method with Eq. (9) and Eq.(10).Figures8(a)-8(b), 8(c)-8(d), and 8(e)-8(f) provide the top, left and right views, respectively.The color bar at the bottom indicates the activation strength of the mean value of HbOs.The mean values of HbOs were averaged across the five trials prior to the standardization in the 0~1 range.As can be seen, the brain activity using our proposed method is distinctive whereas that using the conventional method appears more dispersive.For example, the active locations (i.e., red areas) in Fig.8(a) (the conventional method) are larger than those in Fig.8(b)(the proposed method).For the left and right views, the active regions at the bottom in Figs.8(c) and 8(e) are more dispersive than those in Figs.8(d) and 8(f), respectively.The active locations detected by the proposed method, furthermore, are more diverse than those using the conventional approach (see Figs.8(c)-8(f)).For instance, the active regions in the proposed method are identified at the top side (see Figs.8(d) and 8(f), whereas they are not detected by the conventional one (see Figs.8(c) and 8(e)).Also, this result is consistent across Subjects 1 and 2 (see Fig.9) and Subject 5 (see Fig.10).
Figure 8 plots the brain activation maps of representative Subject 3 in the 3D space.The small circles mark the optodes' positions (see Fig. 1(c) and Fig. 4 for their arrangements).Figures 8(a), 8(c), and 8(e) show the brain activity reconstructed according to the mean values of HbOs (from the 256 channels).It should be noted that the HbOs in this case were computed using the conventional method with constant DPFs of 7.15 and 5.98 at 760 nm and 830 nm, respectively[6,50].Figures8(b), 8(d), and 8(f) show the 3D brain activity constructed according to the HbOs as computed by our proposed method with Eq. (9) and Eq.(10).Figures8(a)-8(b), 8(c)-8(d), and 8(e)-8(f) provide the top, left and right views, respectively.The color bar at the bottom indicates the activation strength of the mean value of HbOs.The mean values of HbOs were averaged across the five trials prior to the standardization in the 0~1 range.As can be seen, the brain activity using our proposed method is distinctive whereas that using the conventional method appears more dispersive.For example, the active locations (i.e., red areas) in Fig.8(a) (the conventional method) are larger than those in Fig.8(b)(the proposed method).For the left and right views, the active regions at the bottom in Figs.8(c) and 8(e) are more dispersive than those in Figs.8(d) and 8(f), respectively.The active locations detected by the proposed method, furthermore, are more diverse than those using the conventional approach (see Figs.8(c)-8(f)).For instance, the active regions in the proposed method are identified at the top side (see Figs.8(d) and 8(f), whereas they are not detected by the conventional one (see Figs.8(c) and 8(e)).Also, this result is consistent across Subjects 1 and 2 (see Fig.9) and Subject 5 (see Fig.10).

Fig. 7 .
Fig. 7.An active (representative) channel (Ch.234) of Subject 3 (blue and red thick curves represent unfiltered and filtered HbX, respectively) using the proposed method: (a) HbO, (c) HbR, (b) and (d) are HbO and HbR averaged across five trials from (a) and (c), respectively, and gray bars denote stimuli of arithmetic tasks (the boundaries in (b) and (d) denote one standard deviation).

Fig. 8 .
Fig. 8.Comparison of the active locations of two methods (3D brain imaging constructed from the mean values of HbOs of Subject 3): The conventional method ((a), (c), and (e)), the proposed method ((b), (d), and (f)), and the color bar at the bottom denotes the activation strength of the mean values of HbOs standardized in 0~1 range.

Fig. 9 .
Fig. 9. Comparison of the active locations of two methods (left view, Subject 1 (first row) and Subject 2 (second row): The conventional method (a, c) and the proposed method (b, d).

Fig. 10 .
Fig. 10.Comparison of active locations in two methods (left view, Subject 4 (first row) and Subject 5 (second row): The conventional method (a, c) and the proposed method (b, d).