Intraoperative quantitative functional brain mapping using an RGB camera

Abstract. Intraoperative optical imaging is a localization technique for the functional areas of the human brain cortex during neurosurgical procedures. However, it still lacks robustness to be used as a clinical standard. In particular, new biomarkers of brain functionality with improved sensitivity and specificity are needed. We present a method for the computation of hemodynamics-based functional brain maps using an RGB camera and a white light source. We measure the quantitative oxy and deoxyhemoglobin concentration changes in the human brain cortex with the modified Beer–Lambert law and Monte Carlo simulations. A functional model has been implemented to evaluate the functional brain areas following neuronal activation by physiological stimuli. The results show a good correlation between the computed quantitative functional maps and the brain areas localized by electrical brain stimulation (EBS). We demonstrate that an RGB camera combined with a quantitative modeling of brain hemodynamics biomarkers can evaluate in a robust way the functional areas during neurosurgery and serve as a tool of choice to complement EBS.


Introduction
Noninvasive functional brain mapping is an imaging technique used to localize the functional areas of the patient brain. This technique is used during brain tumor resection surgery to indicate to the neurosurgeon the cortical tissues which should not be removed without cognitive impairment. Functional magnetic resonance imaging (fMRI) 1 is widely used to localize patient functional areas. However, after patient craniotomy, a brain shift invalidates the relevance of neuronavigation to intraoperatively localize the functional areas of the patient brain. 2 In order to prevent any localization error, intraoperative MRI has been suggested but it complicates the surgery gesture which makes it rarely used. For these reasons, electrical brain stimulation (EBS) 3 is the gold standard for the identification of brain functional areas. For instance, Roux et al. 4 demonstrated that language fMRI data obtained with naming or verb generation tasks, before and after surgery, were imperfectly correlated with electrical brain mapping. The overall results of their study demonstrated that language fMRI could not be used to make critical surgical decisions in the absence of EBS.
In 1977, Jöbsis 5 demonstrated that the blood and tissue oxygenation changes in the brain can be measured using nearinfrared (NIR) light. Then Chance et al. 6 demonstrated that optical imaging (visible and NIR light) can monitor the brain activity with the determination of concentration changes of oxygenated hemoglobin (Δ½HbO 2 ), deoxygenated hemoglobin (Δ½Hb), blood volume (Δ½Hb þ Δ½HbO 2 ), and cytochrome oxidase. Since then, numerous works have demonstrated the ability of optical imaging to detect functional areas thanks to hemodynamics. [7][8][9][10][11][12][13] The motivation for this paper is derived from the necessity to analyze the hemodynamics in brain tissue following the neuronal activation, which are closely linked to the blood oxygenation level-dependent (BOLD) contrast used in functional MRI studies.
During neurosurgery, the craniotomy gives a direct access to the brain cortex. Intrinsic optical imaging [14][15][16] can be used intraoperatively to localize the patient hemodynamic activity in the cerebral cortex. The intrinsic signal refers to the cortical reflectance changes 14,17 due to hemodynamic response. A hyperspectral camera, 18,19 or a single-wavelength illumination in conjunction with a low-noise CCD camera, 16,20 can be used to acquire the intrinsic signal. The time course of this signal is characterized by the early hemodynamic responses in brain tissue related to neuronal activity (initial dip) followed by a larger response that corresponds to the BOLD signal in fMRI studies. 17 This technique is a powerful tool to understand the cognitive functions at the neural circuit level 15 and to define more precisely the hemodynamic response following a physiological stimulus. 21,22 In some studies, a spectroscopic analysis of the intrinsic signal is computed to assess the cortical hemoglobin concentration changes. New approaches consist in using an RGB 21,23 or hyperspectral 18 camera with a continuous wave white light source 18,21,24 or pulsed narrow bandpass illumination sources. 23,25 These setups have the main advantages of being usable in real time 18,25 and directly in the operative room. Steimers et al. 23 analyzed an exposed rat cortex with an RGB-LED light source and an RGB camera. The results of their work indicate that semiquantitative functional maps (in arbirary units) can be processed with the modified Beer-Lambert law. But the optical mean optical path lengths were not taken into consideration. Bouchard et al. 25 developed an ultrafast device made up of two-pulsed LEDs and a monochromatic camera to assess in real time the hemoglobin concentration changes in blood vessels using the modified Beer-Lambert law. Inhomogeneities of the optical properties were not taken into consideration since an averaged path length provided by Monte-Carlo simulation was used in the model. Pichette et al. 18 used a hyperspectral camera and continuous wave white light illumination to assess in real time the absorbance changes. In their study, the blood vessels and the gray matter pixels were automatically segmented by comparing measured reflectance spectra to blood and gray matter simulated reflectance spectra. This allowed the blood vessels to be masked to reduce their influence on absorbance change measurements. In these works, qualitative or quantitative brain maps are not compared to the patient hemodynamic response. 26 So it remains difficult for the neurosurgeon to efficiently localize the functional brain areas with an RGB camera.
The objective of the present work is to supply the methodological tools for the construction of quantitative functional brain maps based on the analysis of the cortical reflectance spectra acquired by a digital RGB camera. For each camera pixel, the intrinsic signal was converted into oxy-and deoxygenated hemoglobin concentration changes using the modified Beer-Lambert law. This law was computed with estimated mean optical path lengths calculated by Monte-Carlo simulations. These mean optical path lengths are chosen according to the local optical properties of the patient brain. Three Monte-Carlo models have been investigated for the study of the light propagation in cortical tissue such as gray matter, surface blood vessel, and buried blood vessel. We propose in our study to complement the intrinsic optical imaging approaches with statistical analyses inspired by the BOLD fMRI method. The Pearson correlation coefficient is calculated between the expected hemodynamic response 26 and each measured concentration changes time courses to accurately localize the functional areas of the patient brain. The brain areas identified by intraoperative electrical stimulation showed a good correlation with the quantitative brain maps processed with the proposed method. These results could help to get robust intraoperative brain area identification based on RGB imaging.

Materials and Methods
A schematic overview of the computation of the functional quantitative maps can be seen in Fig. 1. Once the video was acquired, the following computational steps were applied. The first image was manually segmented into three classes (gray matter, surface blood vessel, and buried blood vessel). This segmentation step aims to associate each camera pixel to the appropriate optical mean path length used in the Beer-Lambert law; see Sec. 2.4. For each frame of the video, the brain repetitive motion was compensated, then data were corrected and filtered. The details of each step are given in Sec. 2.3. A functional model was applied to preprocessed data to compute functional quantitative maps; see Sec. 2.4.

Experimental Setup
The imaging system is composed of an RGB CMOS camera (BASLER acA2000-165uc) in conjunction with an Edmund Optics camera lens (f ¼ 50 mm f∕2 − f∕22), a continuous wave white light source (OSRAM Classic 116-W 230-V light bulb) and a laptop (processor: Intel Core i5-7200U, 2.50 GHz × 4, RAM: 15.3 GiB); see Fig. 2. During data acquisition, the camera also acquired residual light since the operative room lights were on. Data were directly acquired by the laptop via an USB link. A C++ software acquired and processed the images using open source tools such as Qt (v5.9.4), openCV (v3.2.0), 27 FFTW (v3.3.7), 28 and pylon (BASLER library). 8 bits RGB images were acquired every 33 ms (the sampling rate is set to 30 frames per second) with a resolution which at best is 400 ppi (the minimum size of a square pixel is 64 × 64 μm).

Patient Inclusion and Experimental Paradigm
The study was conducted at the neurologic center of the Pierre Wertheimer hospital in Bron, France. Three patients presenting a low-grade glioma close to the motor cortex area were included in the study. All experiments were approved by the local ethics committee of Lyon University Hospitals (France). All participating patients signed written consent. The videos were acquired successively after the patient craniotomy and before the brain tumor resection operation. During the acquisition of the videos, the three patients were awake and under anesthesia (awake surgery).
For videos 1, 2, 3, and 5, the stimulation of the motor cortex was achieved through a repetitive and alternative hand opening and closing at ≈1 Hz. For videos 1, 2, and 5, the hand movement was performed by the patient himself. For video 3, the hand movement was induced by an external person. Video 4 was acquired during the stimulation of the sensory cortex through repetitive fingers and palm caresses at ≈1 Hz. These caresses were performed by an external person. For the five videos, the paradigm consisted of three steps: 30 s of rest, followed by 30 s Fig. 1 Overview of the method. Fig. 2 Schematic of the imaging system.
Neurophotonics 045015-2 Oct-Dec 2019 • Vol. 6 (4) of stimulation, and 30 s of rest. All information regarding the patients and the video acquisition is summarized in the Table 1.
The neurosurgeon performed EBS after RGB imaging to localize patient brain motor and sensory areas. This technique stimulates a neural network in the brain through the direct or indirect excitation of its cell membrane by using an electric current. A bipolar electrode (Nimbus Medtronic neurostimulator) was used in this study. The electrodes are 5 mm apart. A biphasic current was used (pulsating frequency: 60 Hz and pulse width: 1 ms). The current intensity varied during the measurements. At the beginning, the current was set to 1 mA, then the current was increased up to 6 mA. EBS introduces an artificial nonphysiologic signal into the brain. For sensorimotor functions, EBS creates a "positive" effect which is mimicking a sensorimotor behavior. When the motor areas were electrically stimulated, the patient twitched his fingers. When the sensory areas were electrically stimulated, the patient expressed that he felt a sensation in his fingers.

Manual image segmentation
The first image of the video sequence was segmented with a semiautomated procedure into three classes: gray matter, surface blood vessel, and buried blood vessel. Pixels were clustered into four clusters using the K-means algorithm from the C++ library OpenCV (v3.2.0). 27 The components of each cluster were manually sorted and attributed to the three classes. The cluster components that have not been selected corresponded to saturated areas (specular reflection) and were discarded. As the video is motion compensated, the segmentation of the first image is valid during the whole video sequence. The segmentation method is able to detect blood vessels with a diameter larger than 500 μm. The blood vessels with a diameter smaller than 500 μm are included in the gray matter class. This segmentation step aims to associate each camera pixel to the appropriate optical mean path length used in the modified Beer-Lambert law; see Sec. 2.4.

Motion compensation
After craniotomy, the brain surface undergoes a repetitive motion due to the patient breath and cardiac pulsation. This motion as well as potential video camera motion prevents accurate video analysis. The motion compensation aims to ensure that each camera pixel corresponds to the same cortical area all along data acquisition. Sdika et al. 29,30 proposed a repetitive motion compensation algorithm. This algorithm is split into two parts. First, the basis of the repetitive motion is learned from few initial frames of the video (50 frames), then, each video frame is registered to the reference image (first image of the video). To ensure that motion compensation worked in normal operation, each registered image is compared to the reference image using the normalized cross covariance (NCC): E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 3 2 6 ; 5 5 7 NCCðtÞ ¼ n P IðtÞ is the image registered at time t and I 0 the reference image (first frame). n denotes the number of pixels of the images. . The NCC dispersion range of the registered videos is defined as the area between the curves y 1 ðtÞ ¼ μ mes ðtÞ − σ mes ðtÞ and y 2 ðtÞ ¼ μ mes ðtÞ þ σ mes ðtÞ. The motion compensation of the five videos is validated if the NCC dispersion range of the registered videos is included in the acceptable NCC variation range.

Data filtering
According to Chance et al., 6 low-frequency modulation of light absorption is linked to cortical activity. These frequencies of interest can be visualized in the expected hemodynamic response spectrum. The hemodynamic response of the brain in relation to the neural activities (cortical stimulus) can be expressed by the hemodynamic impulse response function (HIRF). 26 In our application, the expected hemodynamic response can be obtained by convolving HIRF with the window function P representing the experimental paradigm; see Fig. 3. The expected hemodynamic response is the ideal representation of the hemodynamic response of the patient whose validity is discussed in Sec. 4. A low-pass filtering was implemented by multliplying the signal by a

Data correction
According to Oelschlägel et al., 24 data preprocessing is mandatory to correct the slow drift of the collected intensity due to the cortical tissue desiccation during the video acquisition.
The assumption is that the beginning and the end of the video corresponded to the same patient physiological state. So we considered that these intensity values had to be identical. Linear regressions were computed on the measured RGB time courses, then the calculated regression lines was subtracted to original data to compensate the slow drift of the collected intensity.

Functional Model
The standard approach for the analysis of the reflected spectra is based on the modified Beer-Lambert law. The assessment of the concentration changes depends on the determination of the mean optical path length of the detected photons; see Eq. (3). For the purpose of measuring more precisely the hemoglobin concentration changes, three different cortical tissues have been modeled; see Fig. 4. Monte-Carlo simulations were processed to estimate their mean optical path length which were assigned to pixels belonging to the segmented classes; see Sec. 2.3.1.

Modified Beer-Lambert law
The modified Beer-Lambert law can be expressed as a matrix system: 31 where E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 3 2 6 ; 3 3 8 ΔA i ðtÞ is the absorbance changes measured at time t: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 3 2 6 ; 2 8 4 where R i ðtÞ is the reflectance intensity measured at time t by the camera color channel i (red, green, or blue), R 0 i is the reference reflectance intensity measured by the camera color channel i (average of the reflectance intensity over the duration of the first rest step of the experimental paradigm; see Sec. 2.2). ε n is the extinction coefficient of the chromophore n (in L mol −1 cm −1 ). Δ½Hb is the deoxygenated hemoglobin molar concentration changes (in mol L −1 ) and Δ½HbO 2 the oxygenated hemoglobin molar concentration changes (in Mol L −1 ). Our model takes into consideration the receiving spectrum of the RGB camera and the emission spectrum of the light source. The spectral sensitivity of the detector i of the RGB camera is represented by D i ðλÞ and SðλÞ is the normalized intensity spectrum of the light source. LðλÞ is the wavelength-dependent mean optical path length of the photons traveled in tissue. Hemoglobin concentration changes are obtained by matrix inversion once the matrix E has been calculated. However, the mean optical path length LðλÞ needs to be determined.

Mean optical path length determination
The mean optical path length LðλÞ was obtained by Monte-Carlo simulations using MCX software. 32 Three cortical tissues were modeled under a homogeneous white light illumination; see Fig. 4. Model 1 represents a 2-mm diameter blood vessel on the surface of the cortical tissue, model 2, a 2-mm diameter blood vessel buried under 1-mm of cortical tissue, and model 3, a cortical tissue without a large blood vessel. The cortical tissue is composed of gray matter perfused by capillaries. Each voxel of the modeled tissues included the information of optical parameters (absorption, reduced scattering, anisotropy coefficients, and refractive index). A homogeneous illumination was achieved by scanning the emission point of the light source along the entire model surface. A white light illumination has been simulated by scanning the optical parameters along the entire illumination spectrum (from 400 to 1000 nm in steps of 10 nm). A total of 50 × 50 × 61 × 3 simulations were computed. 50 × 50 represents the number of simulations used to illuminate the modeled cortical tissues in a homogeneous way. 61 denotes the number of emitted wavelength and 3 the number of modeled cortical tissues. The optical parameters were taken from the literature and correspond to a nominal physiological condition. References 33-35 have been used for the blood vessel parameters, and Refs. 35-38 for the gray matter parameters.
The size of the modeled tissues has been chosen in accordance with the photon sensitivity profile 39 computed for the detector situated at the center of the top face of model 3. A detector collects all the photons reaching the surface in a square area of 1 mm 2 . Model 3 has been chosen because the photons traveling through the cortical tissue without a large blood vessel have a greater probability to have a long trajectory than the photons traveling through the cortical tissue with a large blood vessel. The photon sensitivity profile is roughly represented as a half sphere of radius of 12.5 mm. To avoid any photon loss and inexact results due to the boundary conditions (a simulation of the travel of a packet of photon stops when this packet of photon leaves the volume), the size of the models is set to 50 × 50 × 50 voxels with a resolution of 1 mm 3 . This allows us to precisely compute the mean optical path length of the photons reaching the detectors located at the top face of the volume.
For the three models, the mean optical path length L of all detected photons is calculated; see Fig. 5. For λ < 580 nm, the mean path length of model 1 is on average 39% smaller than the one of model 3. The one of model 2 is on average 4% smaller than the one of model 3. For λ > 580 nm, the mean path length of model 1 is on average 42% smaller than the one of model 3. The one of model 2 is on average 34% smaller than the one of model 3. When the blood vessel is buried under 1 mm of gray matter (model 2), a smaller proportion of photons is absorbed in model 2 than in model 1. However, this blood vessel still has a large impact on the photon propagation.

Quantitative functional brain map
Two quantitative functional brain maps were computed by selecting the Δ½Hb and Δ½HbO 2 time courses which were correlated with the expected hemodynamic response. The Pearson correlation measurement aims to highlight the hemodynamic response associated with the experimental paradigm and ignore any other hemodynamic variations. Only positive correlations were considered, so the Δ½Hb time courses were multiplied by −1. Two quantitative functional brain maps were defined as the Hb and HbO 2 concentration changes averaged over the duration of the patient hand stimulation (see Sec. 2.2). The quantitative maps were only processed for pixels for which the Pearson correlation coefficient was higher than the Pearson correlation coefficient threshold (PCCT) value: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 3 2 6 ; 2 1 7 QMap C ðx; yÞ ¼ Cðx; yÞ; if r C ðx; yÞ ≥ PCCT nonprocessed; otherwise: ; (5) where E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 3 2 6 ; 1 6 1 Cðx; yÞ ¼ QMap C is the chromophore C quantitative functional brain map. C represents either HbO 2 or Hb. ðx; yÞ denotes an image pixel position. n is the frame index, N 1 and N 2 are the patient hand stimulation starting and ending frame indexes. r C ðx; yÞ denotes the Pearson correlation coefficient calculated between   (5) and (6)].
Welch's T tests were computed to score the differences between the mean values of the Cort NF and Cort F areas. These means were denoted hHbi, hHbO 2 i, hr Hb i, and hr HbO 2 i. If all four tests reject the null hypothesis at 1% significance level, the Cort F area is defined as an activated cortical area.

Motion Compensation
In Fig. 6, the blue rectangular area corresponds to the acceptable NCC variation range defined as the area between the horizontal lines y ref 1 5 The red, blue, and gray curves represent the computed wavelength dependent mean optical path length of models 1, 2, and 3, respectively (see Fig. 4).

Quantitative Functional Maps
The Hb and HbO 2 quantitative functional maps of videos 1-5 are represented in Fig. 7. white circles of 10 pixels diameter for video 1 and 30 pixels diameter for videos 2 to 5. For each video, Hb and HbO 2 quantitative functional maps are plotted according to four PCCT values (0, 0.3, 0.5, and 1). Only the Δ½Hb and Δ½HbO 2 time courses whose Pearson correlation coefficient is higher than the PCCT value are considered in Eq. (5). The highest QMap Hb and QMap HbO 2 values are situated at the level of the blood vessels surrounding the motor and sensory areas.

Video 1
The stimulation of the cortex was achieved through a repetitive and alternative hand opening and closing at ≈1 Hz (the movement was performed by the patient). We defined five points of interest within video 1 to explore in more details the hemodynamic time courses in the different brain areas. In Fig. 8(a), the point M 1 is situated at the center of the M 1−1 area. The point S 1 is situated at the center of the S 1−1 area. The point A is situated on the blood vessel above the M 1−1 area. The point B is localized on a blood vessel which is not linked with the functional areas. The point C represents a point gray matter situated far away from the functional areas. The Δ½Hb and Δ½HbO 2 time courses of these points of interest can be compared in Fig. 8(b).

Video 2
The stimulation of the cortex was achieved through a repetitive and alternative hand opening and closing at ≈1 Hz (the movement was performed by the patient). In Fig. 7

Video 3
The stimulation of the cortex was achieved through a repetitive and alternative hand opening and closing at ≈1 Hz (the movement was induced by an external person). In Fig. 7

Video 4
The stimulation of the cortex was achieved through a repetitive fingers and palm caresses at ≈1 Hz (the movement was induced by an external person). In Fig. 7

Video 5
The stimulation of the cortex was achieved through a repetitive and alternative hand opening and closing at ≈1 Hz (the movement was performed by the patient). In Fig. 7

Statistical Analysis
The functional areas of the patient exposed cortex which have been identified by EBS are compared to a reference cortical area C. This area C has not been identified as functional by the EBS. The objective is to compute the activation criterion introduced in Sec. 2.4.4 on each functional area to evaluate our method sensitivity. For each video, an ANOVA and Kruskal-Wallis H-test have been separately computed on the Hb, HbO 2 , r Hb , and r HbO 2 distributions of the areas defined in Fig. 7. These tests reject the null hypothesis that two or more distributions have the same population mean. Each video is studied separately. For each distribution of each video, a Welch's T-test is computed between the sample X i−j and C j (X represents either the motor M or sensory S areas i of the patient j defined in Fig. 7) with the null hypothesis that the two samples have identical mean values. In Figs. 9 and 10, the distribution of the Hb and HbO 2 values of the areas defined in Fig. 7 are represented. In Figs. 11 and 12, the distribution of the r Hb and r HbO 2 values of the areas defined Table 2 Hb, HbO 2 values (in μmol L −1 ) and r Hb , r HbO 2 values of the curves plotted in Fig. 8(b).   Fig. 7. The black diamonds represent the mean values of the distributions (hHbO 2 i) and the half length of the blue lines the standard deviation values. Each video is studied separately. The notation i Ã indicates a T -test's statistical significance for the comparison of the means of the distribution i and C j (j represents the patient id).   Table 3.

Discussion
In videos 1, 2, and 5 (see Table 1), the statistical results of Table 3 indicate that the motor and sensory areas identified by EBS correspond to the highlighted cortical areas of Fig. 7. This demonstrates the ability of our functional model to identify in a robust manner these functional areas. In videos 3 and 4, the motor and sensory cortex was explored in a more subtle way. In video 3, the hand movement was induced by an external person. All motor and sensory areas are defined as activated cortical areas except the S 2−2 area. In video 4, the stimulation of the cortex was achieved through a repetitive fingers and palm caresses. Only the S 1−2 area is defined as an activated cortical area. These results seem to indicate that the S 1−2 area could be linked to somatosensory functions, whereas the S 2−2 and S 3−2 areas to sensorimotor functions. This indicates that our model seems to be rather sensitive to subtle differences in the physiological stimuli of the hand.
These results should be taken with caution because our functional model needs to be improved. In particular, a physiological a priori on the hemoglobin concentration change values could be incorporated in the model. However, there is still no consensus in the literature on such threshold values linked with a physiological stimulus. In video 3, when the Hb and HbO 2 distributions of the S 2−2 area are statistically significant (see Table 3), the hHbi and hHbO 2 i values of this area differ from the area C 2 in the sense of a deactivation of the cortical area, which is difficult to interpret (the hHbi value of the S 2−2 area is higher than the one of the C 2 area and the hHbO 2 i value of the S 2−2 area is smaller than the one of the C 2 area). The S 2−2 area is defined as nonactive because its hr Hb i value is not statistically different from the one of the C 2 area. This result shows that our functional model is able to eliminate this nonphysiologic hemodynamic event by comparing the similarity of the hemoglobin changes time courses to the expected hemodynamic response. This confirms that our model is more robust than models with a nonphysiological a priori. The same phenomenon can be observed for the M 1−2 area of video 4.
Our functional model has to be used with the EBS. Indeed, it requires a reference cortical area which is defined as a nonfunctional area by the EBS. The main limits are that the results of the statistical analysis depend of the choice of the reference and the definition of the cortical activation is not obtained for each camera pixel but for a group of pixels. A more robust definition of cortical areas will be explored in future works, notably with the implementation of a SPM 40 like analysis.
The Hb, HbO 2 , r Hb , and r HbO 2 noise levels in our measurements were calculated by taking the standard deviation of the mean values of 10 × 5 cortical areas which have not been identified as functional areas by the EBS (10 measurements were realized in each of the five videos). The Hb and HbO 2 noise levels are equal to 0.26 μmol L −1 , and the r Hb and r HbO 2 noise levels are equal to 0.07. In video 4, the S 2−2 area is not defined as an activated area only because the hHbi value is not statistically different from the one of the C 2 area. The mean values of the other distributions are statistically different from the ones of the C 2 areas and the sign of the differences corresponds to the sense of an activation of the cortical area. When comparing the hHbi, hHbO 2 i, hr Hb i, and hr HbO 2 i values of the S 2−2 and C 2 areas, the differences between these values are approximately equal to the noise levels. Thus by increasing the signal-to-noise ratio, the activated cortical areas could be determined in a more robust way.
The five videos have been studied separately, so no general PCCT value was proposed to localize the activated cortical areas in the five videos. The Hb, HbO 2 , r Hb , and r HbO 2 distributions of the functional areas identified by EBS have different mean and standard deviation values. For instance, the r HbO 2 values of the M 1−2 area in video 2 range between 0.42 AE 0.25 and the ones of the same motor area range between 0.53 AE 0.22 in video 3. The r HbO 2 values of video 3 are higher than the ones of video 2. This may be due to a more constant hand movement in video 3 than in video 2. Indeed, in video 2, the patient has moved himself his Oct-Dec 2019 • Vol. 6 (4) hand, whereas the movement was induced by an external person in video 3. In our study, the stimulation step of the experimental paradigm (see Sec. 2.2) is not exactly constant during 30 s (patient sleepness and irregular movement). This implies that the experimental paradigm function plotted in Fig. 3 does not exactly fit a window function. Since an ideal expected hemodynamic response (convolution of the HIRF 26 with a window function) is compared to the Δ½Hb and Δ½HbO 2 time courses, the r Hb and r HbO 2 distributions may differ depending on the quality of the stimulus. Furthermore, habituation is known to interfere with the linearity of the hemodynamic response during a 30-s plateau. This would further modify the actual stimulus from the ideal window function.
The model of the HIRF may be improved. The positive BOLD signal used in MRI studies occurs in superficial cortical layers (0 to 1 mm). This positive BOLD signal is often accompanied with a negative BOLD signal occurring in deeper cortical layers (1 to 2 mm). 41 The HIRF may be a conjunction of these two BOLD signals. In our study, we consider that the theoretical Hb impulse response function correspond to the opposite of the HbO 2 impulse response function. This assumption may introduce some errors since it has been shown that these HIRFs do not correspond to an opposite operation. 42 It also has been shown that the HIRF is different in the gray matter and in the arterioles. 42 To improve our functional model, several HIRFs may be determined for each camera pixel depending on its biological attribution (gray matter, surface blood vessel, and buried blood vessel). For this purpose, the segmentation step (see Sec. 2.3.1) could be used.
The highest Hb and HbO 2 values are localized at the level of the blood vessels (see Table 2). Indeed, the partial volume effect has less impact on the blood vessels than on the gray matter. So the underestimation of the concentration changes is lower at the level of the blood vessel than at the level of the gray matter. It is also interesting to notice that the concentration changes of the blood vessels which are not directly linked to the functional areas are higher than those of gray matter pixels associated with the sensory and motor areas; see Fig. 8. Similar observations can be found in the literature. 18,25 Pichette et al. 18 proposed to mask the blood vessels in order to render more clearly the smaller hemodynamic changes in gray matter. The computation of the quantitative maps in conjunction with the measurement of the Pearson correlation coefficient between the concentration changes time courses with the expected hemodynamic response allows us to visualize the activated cortical areas without masking the blood vessels. With the results of this study, it is complicated to understand the role of the blood vessels. Indeed, it is impossible to get the direction of the blood flow and thus conclude that a blood vessel is associated with a functional area. A complementary imaging modality such as speckle imaging 43 could help to interpret the role of these blood vessels.
In Sec. 2.4.2, three simple models for the light transport of photons in cortical tissue are defined for the analysis of reflectance spectra. These models are applied to each camera pixel depending on its biological attribution (gray matter, surface blood vessel, and buried blood vessel). This segmentation allows us to approach the optical properties inhomogeneities of a cortical tissue. If the mean path length of the gray matter model (see the curve of model 3 in Fig. 5) was used for the computation of the Δ½Hb and Δ½HbO 2 time courses of the surface blood vessel point A, this would imply a 40% underestimation of the Hb value and 35% underestimation of the HbO 2 value.
In our model, the distance of a gray matter pixel to the nearest blood vessel is not taken into consideration. We consider that the mean path length associated with this pixel is model 3 curve in Fig. 5. This is an approximation since the blood vessel still has an impact on the photon propagation. To solve this problem, a pixel-wise determination of the optical mean path length can be implemented. This could be realized by fitting theoretical reflectance spectra (computed for a great number of Monte-Carlo models) with measured reflectance spectra. This comparison appears not to be an easy task since only three spectral values can be obtained which are integrated over the broad wavelength ranges covered by each detector. A hyperspectral camera seems to be a more appropriate solution to answer this problem with the acquisition of high-resolution reflectance spectra. A similar method has been implemented by Pichette et al. 18 using a hyperspectral imaging device and a two class image segmentation but does not take account of the distance between cortical tissues and blood vessel.
The quantitative maps are displayed after the end the data acquisition. Unlike Bouchard et al.'s 25 work, our method is not processed in real time. But it would be possible to develop a real-time processing using a multithreading software architecture and a code optimization using GPU (graphics processing units).

Conclusion
A quantitative and noninvasive method for imaging the oxygenated hemoglobin concentration changes (Δ½HbO 2 ) and deoxygenated hemoglobin concentration changes (Δ½Hb) of a human exposed cortex using a digital RGB camera was demonstrated in the present report. The results showed that the method can be used by the neurosurgeon to localize the motor and sensory areas before a tumor removal surgery. The brain areas defined as activated cortical areas by our functional model are well correlated with the functional areas localized by EBS. Our model is based on quantitative measurements and physiological and statistical comparisons. Each step of the model can be improved, such as the noise level reduction, the consideration of a physiological a priori in the statistical comparisons, or the precision of the quantitative measurements. Indeed, the estimation of three different path lengths associated with three different cortical tissues was used to process the Beer-Lambert law. The Monte-Carlo models used in our study remain an approximation of the optical properties inhomogeneities of the patient brain. To increase the accuracy of our method, a detailed tissue structure that best matches the cortical pattern of the acquired images can be implemented in the Monte-Carlo simulations. However, this would imply the determination of a three-dimensional cortical tissue from an RGB image and would require considerably more computational time and a large amount of simulated data. These directions will be explored in future works. The current work demonstrates that an RGB camera combined with a quantitative modeling of brain hemodynamics biomarkers could evaluate in a robust way the functional areas during neurosurgery. This strengthens the relevance of using a classical RGB camera for functional intraoperative brain imaging.

Disclosures
No conflicts of interest, financial or otherwise, are declared by the authors.