Open Access
Add to BibSonomy
Abstract
Spectral unmixing designates techniques that allow to decompose measured spectra into linear or non-linear combination of spectra of all targets (endmembers). This technique was initially developed for satellite applications, but it is now also widely used in biomedical applications. However, several drawbacks limit the use of these techniques with standard optical devices like RGB cameras. The devices need to be calibrated and a a priori on the observed scene is often necessary. We propose a new method for estimating endmembers and their proportion automatically and without calibration of the acquisition device based on near separable non-negative matrix factorization. This method estimates the endmembers on spectra of absorbance changes presenting periodic events. This is very common in in vivo biomedical and medical optical imaging where hemodynamics dominate the absorbance fluctuations. We applied the method for identifying functional brain areas during neurosurgery using four different RGB cameras (an industrial camera, a smartphone and two surgical microscopes). Results obtained with the auto-calibration method were consistent with the intraoperative gold standards. Endmembers estimated with the auto-calibration method were similar to the calibrated endmembers used in the modified Beer-Lambert law. The similarity was particularly strong when both cardiac and respiratory periodic events were considered. This work can allow a widespread use of spectral imaging in the industrial or medical field.
© 2023 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
Spectral unmixing is a technique that was initially developed for satellite imaging applications. The spectrum recorded in satellite images is a sum of spectra of all targets collected by the sensor. This implies that the measured spectra are defined as a mix of several constituent spectra (or endmembers), weighted by their proportion in the pixel [1,2]. The mixing process can be considered linear [3] or non-linear [4]. These models are used in remote sensing application to characterize and discriminate materials on the Earth’s surface [4,5], but also to decompose mixtures by spectral features [6,7].
Spectral unmixing is also used in biomedical applications, especially for in vivo optical imaging [8–10]. The objective is to observe non-invasively biological tissues and help in understanding, detect or track metabolic, functional or disease processes in the body. Sowa et al. evaluated the hemodynamic changes in the early post-burn period of skin by near infra-red reflectance spectroscopy [11]. In this work, a calibrated spectrometer and a linear mixing model was used to estimate the proportion of oxy- and deoxyhemoglobin by fitting measured absorbance to the hemoglobin molar extinction spectra (endmembers). Berman et al. implemented an unmixing method with hyperspectral images acquired on cervical tissue [12] to establish a library of biologically interpretable endmembers. Alston et al. have established a linear spectral unmixing model to evaluate the proportion of the two states of 5-ALA induced PpIX fluorescence spectra [13] in low and high grade gliomas [14]. This work has improved fluorescence-guided glioma resection by discriminating healthy tissue from tumor margins. Lu et al. proposed spectral unmxing to decompose hyperspectral images collected on a head and neck cancer animal model into oxy- and deoxygenated hemoglobin concentration maps [15]. The method is based on non-negative matrix factorization and hemoglobin extinction molar spectra were used as initial inputs. This method allows to estimate concentration maps of chromophores, but an initial a priori on the chromophore spectral signature is required.
Non-invasive 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 that should not be removed without cognitive impairment. Functional magnetic resonance imaging (fMRI) [16] is the preoperative gold standard for the identification of the patient brain functional areas. However, after patient craniotomy, a brain shift invalidates the relevance of neuronavigation to intraoperatively localize the functional areas of the patient brain [17]. In order to prevent any localization error, intraoperative MRI has been suggested, but it complicates the surgery gesture, which makes it rarely used. During neurosurgery, electrical brain stimulation [18] (EBS) is the gold standard, but this technique is mainly limited by its low spatial resolution ($\approx 5$ mm [19]) and has the potential risk to trigger epileptic seizures. This technique allows a robust and reliable detection of many functional areas, but could be traumatic for the patient, by inhibiting certain cognitive functions such as speech for example. EBS is also complicated to perform and requires a very strong expertise.
As a complement to EBS, optical imaging provides an ideal solution for intraoperative functional brain mapping. The analysis of the light absorption allows to monitor the brain activity (motor or sensory tasks for example) with quantification of the concentration changes in oxy- ($\Delta C_{HbO_2}$) and deoxygenated hemoglobin ($\Delta C_{Hb}$) in brain cortex [20–28]. In a previous study, we implemented a spectral unmixing method using an RGB camera to quantify hemodynamics in patient brains: the differential form of the modified Beer-Lambert law [21]. This method extracts $\Delta C_{HbO_2}$ and $\Delta C_{Hb}$ by unmixing the measured absorbance changes in a linear combination of calibrated hemoblobin extinction spectra (endmembers). These endmembers were calculated in two steps. First the hemoglobin extinction spectra were multiplied to the optical mean path length of traveled photons, the light source spectra and the spectral sensitivities of the camera. Then, we integrated this product over the spectral range of the camera detectors to obtain the endmembers. This is the standard approach for identifying functional brain areas using near-infrared spectroscopy devices [29–31]. In our previous studies, we also applied this technique during neurosurgery [20–22].
Several drawbacks limit the use of spectral unmixing techniques with standard optical devices like RGB cameras for intraoperative functional brain mapping. The devices must be calibrated, sometimes in complex environments such as operating rooms. For some models such as the modified Beer-Lambert law, a a priori on the observed scene is also necessary: the optical properties of the tissue has to be known for estimating the optical mean path length of traveled photons. Without considering these two aspects, quantification of concentrations or abundances contains large estimation errors [32]. In this paper, we introduce a new method called auto-calibration method for estimating endmembers and their proportion automatically and without calibration of the acquisition device. Endmembers are estimated blindly using a near separable non-negative matrix factorization (NMF) [33]. In traditional approaches, periodic absorbance variations due to heartbeat and respiration decrease the sensitivity to the desired contrast (for example, the low-frequency modulation of light absorption linked to cortical activity [34]). These variations are thus removed from the initial signal with a filtering operation [21,23]. Unlike traditional approaches, this method estimates the endmembers in spectral bands where absorbance changes periodically before applying filtering. Indeed, the a priori on the physiological characteristics of periodic fluctuations due to heartbeat and respiration is very robust compared to that associated with brain activity. This makes it easy to remove their contribution using linear filtering and it simplifies the optical model since it allows to estimate the endmembers automatically based on the periodic temporal variations of their abundance. The algorithm used for the estimation of the endmembers is fast (can run in real time) and robust to small noise levels.
The auto-calibration method was applied to RGB videos collected during neurosurgery for the identification of functional brain areas. Four different setups were used in a clinical study that included $12$ patients and $16$ acquisitions. RGB data were collected with a smartphone, surgical microscopes or the experimental setup we used in our previous studies [20–22]. Functional brain areas were identified with Statistical Parametric Mapping (SPM) technique [35–37] and the biomarkers computed with the modified Beer-Lambert law and the auto-calibration method. These identifications were compared to those obtained with electrical brain stimulation. We show that identifications provided by auto-calibration method are consistent with those obtained with the modified Beer-Lambert law and electrical brain stimulation. The results showed that the endmembers estimated with the auto-calibration method were similar to the calibrated ones used in the modified Beer-Lambert law. In particular, heartbeat-derived signals are suitable for estimating oxygenated hemoglobin changes in tissue and respiration-related events are suitable for estimating oxygenated and deoxygenated hemoglobin changes. The similarity was particularly strong when both cardiac and respiratory periodic events were considered. However, the method was put into default for high noise levels.
2. Material and methods
2.1 Linear spectral unmixing
2.1.1 Modified Beer Lambert law
Let us assume that the absorbance $A$ can be expressed with modified Beer-Lambert law:
Using the Beer law, changes in the absorption coefficient can be linked to the concentration changes $\Delta C$ (in $mol.L^{-1}$) of the chromophore within the medium:
$\epsilon _n$ is the molar extinction coefficient of chromophore $n$ (in $L.mol^{-1}.cm^{-1}$) and $N$ the number of chromophore that compose the medium, that is usually unknown. Using a spectral acquisition device having $K$ spectral channels, Eq. (3) can be expressed by the following matrix system [21,39]:The elements of the matrices $E$, $B$ and $G$ have been calculated as follows:
The spectral sensitivity of the detector $i$ of the spectral acquisition device is represented by $D_i(\lambda )$ and $S(\lambda )$ is the normalized intensity spectrum of the light source. $[\lambda _1;\lambda _2]$ is the wavelength range collected by the spectral acquisition device.
In biomedical applications, especially for functional brain mapping, it is commonly assumed that $L$ does not vary with $\Delta \mu _a$ and $\Delta \mu _s$ (which is not strictly correct [38]). This implies that matrix $B$ in Eq. (5) is null. With this assumption we can write:
In order to estimate $\Delta C$ values, elements of matrix $E$ and $\Delta G$ need to be calculated. As cerebral activity is related to slow events [23], it is common to apply a low-pass filtering on measured $\Delta A$ to remove the physiological noises such as heart rate and breathing rate. The low-pass filtering was operated in the Fourier domain with a Blackman window ($Fc = 0.05$ Hz) [20–22]. Elements of matrix $\Delta G$ depend on the geometry. In functional brain mapping, these factors are mainly due to tissue dessication [40] that are estimated with linear regression. Once matrix $\Delta G$ has been removed from Eq. (7) by band-pass filtering, we obtain Eq. (8):
For the calculation of matrix $E$, we consider that the variations of absorbance due to the changes of absorption are mainly due to two or three chromophores ($N=2$ or $N=3$ in Eq. (5)), whereas these variations are in reality due to a very great unknown number of chromophores.
In functional brain mapping studies, these chromophores are the oxy- and deoxygenated hemoglobin ($HbO_2$ and $Hb$), and sometimes, the oxidative state of cytochrome-c-oxidase ($CCO$) is also considered [41]. The so-called calibration is performed for calculating the matrix $E$, see Eq. (6). In this equation, the optical path length $L \left (\mu _a^{ref},\mu _s^{ref}\right )$ need to be estimated, Monte Carlo simulations [21] or the diffusion approximation [42,43] could be used for this purpose. Finally, chromophores’ concentration changes are obtained by matrix inversion.
2.1.2 Auto-calibration method
The aim of the auto-calibration method is to estimate the matrix $E$ (see Eq. (6)) in a space where Eq. (8) is valid. The term “auto-calibration” refers to the estimation of the $E$ matrix without any prior knowledge (wavelength-dependent parameters $L$, $D$, $S$ or $\epsilon$, see Eq. (6)) and using only measurements.
As we saw in section 1, the modified Beer-Lambert law is an approximate method based on assumptions that are not strictly correct. The principle of the auto-calibration method is to find a model that simplifies the expression of $\Delta A$ into two components (or sources) in a space where Eq. (8) is valid. This space can be found by studying the frequencies of $\Delta A$, for which $\Delta \mu _a$ is very low and $\Delta \mu _s=0cm^{-1}$. These conditions can be met if the concentration of chromophores varies weakly and periodically over time. Under these conditions, matrix $B$ in Eq. (5) is null because $\Delta L(\lambda,t) \approx 0cm$. The factor $\Delta G$ in Eq. (5) may be approximated to $0$ if it is evaluated in a narrow and high frequency bandwidth and not linked to periodic event, which leads us to Eq. (8).
For functional brain mapping, this space could be related to heart rate (narrow frequency band around $1$ Hz) or breathing rate (narrow frequency band around $0.2$ Hz) [44]. It is interesting to note that $\Delta A$ fluctuations related to these physiological events are usually considered as noise and are discarded with a low-pass filtering. Under normal conditions, hemoglobin is conducted from the heart to the brain through arteries with each heartbeat. This hemoglobin is highly oxygenated and is directly related to the oxygen supply to organs. Its normal saturation in oxygen lie between $95{\% }$ and $100{\% }$ [45]. This supply of hemoglobin perfuses the brain tissue and thus oxygenates the functional areas following neuronal activity. After perfusing the brain tissue, blood is returning to the heart through the vein with an oxygen saturation greater than $75{\% }$ [46]. As absorption changes induced by the heartbeat are only due to $Hb$ and $HbO_2$, heartbeat can be roughly modeled as two periodic functions. In this space, it is valid to assume that only two chromophores vary ($N=2$ in Eq. (5)). It is acceptable to assume that there is no change in $\mu _s$ since only concentrations of $Hb$ and $HbO_2$ are periodically varying. These changes in $\mu _s$ are mainly induced by cell conformational changes and swelling [47] and are no larger than $0.4{\% }$ [48] in the low frequency range. $\mu _a$ varies much less than in the case of cerebral activity because the changes in blood flow and volume on the scale of the heartbeat or respiration are only effects of elasticity of the vessels which are thus weakly significant compared to the effects of the physiological cascade related to the cerebral hemodynamic response. Thus, in this space it is valid to assume that the matrix $B$ in Eq. (5) is null. Finally, since we focus on frequencies higher than those related to the matrix $\Delta G$, we can consider their elements as zero.
The auto-calibration method consists of several steps. First spectra of absorbance changes signals were computed (magnitude of the Fourier transform). Negative frequencies of the spectra were discarded and only frequencies centered on the periodic event were kept. This reduced spectra corresponds to the matrix $X$ in Eq. (9). Then, a near separable non-negative matrix factorization (NMF) was computed to estimate the matrix $E$ in Eq. (8). Non-negative matrix factorization is a low rank approximation model in which an input matrix $X$ is factorized into two element-wise non-negative matrices $W$ and $H$:
In Eq. (9), $X$, $W$ and $H$ have dimensions $(m \times n)$, $(m \times k)$ and $(k \times n)$, respectively. Integers $m$ and $n$ are the number of rows and columns of $X$, respectively. The integer $k$ is called the factorization rank and verifies $1 \leq k \leq min(m,n)$. Separable NMF is a variant of NMF, solvable in polynomial time with a robust algorithm [49].
Algorithm 1. Rank -2 near separable NMF
Algorithm 1 can be used to solve separable NMF with factorization rank $k=2$ [33]. Rank-2 separable NMF is a special case of separable NMF where a solution can be computed in closed form. Indeed, the data points (matrix $X$) lie approximately on a segment which the two extreme points are the two sources (matrix $W$). Finding one of these two extreme points is equivalent to finding the data point with the largest norm (first source). In intraoperative optical imaging, data are subject to acquisition noise. Thus, the extreme points may not be associated with sources but with acquisition noise. Indeed, a large amount of noise may introduce outliers in data which will distort the estimation of the sources. To overcome this drawback, we propose to find the sources by considering the mean value of the $P$ extreme points, a method adapted from [50]. The periodic events that may be considered in the auto-calibration method are heartbeat and respiratory. As hemodynamic changes, and especially changes in $\Delta C_{HbO_2}$ are predominant for frequencies related to these physiological events, the first source can be associated with $E_{HbO_2}$ and the second source with $E_{Hb}$ (see Eq. (6)). The algorithm is not iterative and a solution is computed in $O(2m+2n\log _2(n))$, which is compatible with real-time processing.
We applied a few steps of a descent algorithm for non-negative factorization on $X$ with an initial guess $W$ obtained with Algorithm 1. This allows to refine the estimation of the sources by limiting the impact of the acquisition noise. For this purpose, we used the non-negative factorization techniques toolbox in python nn-fac [51].
To compute slow fluctuations of sources $W$, $\Delta A$ signals were low-pass filtered ($\Delta A^{LP}$ of frequencies $F<0.05Hz$ [21]) to isolate the low-frequency modulation of light absorption linked to cortical activity [23]. Finally, we used the $\Delta A^{LP}$ signal to compute slow fluctuations of sources $\Delta C^{auto}$ (in arbitrary units):
with $^+$ the right pseudo-inverse operator.2.2 Clinical application: intraoperative identification of functional brain areas during brain tumor resection operations in human patients
We represented the flowchart of data analysis in Fig. 1. First acquisition and pre-processing steps were executed. Once the RGB video was acquired, the brain repetitive motion was compensated, data were filtered, and low-pass filtered absorbance changes were calculated.
Fig. 1. Flowchart of data analysis, including definitions of the variables used in this manuscript. The flowchart was separated into four parts. First, data were acquired and pre-processed. After video acquisition, motion in images was compensated, RGB reflectance intensities $I_k(p,t)$ measured at time $t$ for the pixel $p$ were filtered and converted into absorbance changes $\Delta A^{LP}_k(p,t)$ ($k$ designate the color channel of the camera: red, green, or blue). The second and third parts are the flowchart are the parallel execution of the modified Beer-Lambert law and auto-calibration method. The modified Beer-Lambert law used the matrix $E$ to convert $\Delta A^{LP}_k(p,t)$ into concentration changes $\Delta C_n(p,t)$, $n$ designates either $HbO_2$ or $Hb$. The auto-calibration method was applied to the unfiltered temporal intensity signal, averaged over the surgical window $I_k(t)$. Spectra of absorbance changes were calculated $|TF(\Delta A_k)|$, specific frequency bands were selected and separable NMF was executed to calculate matrix $W$ (estimate of $E$). The auto-calibration method used the matrix $W$ to convert $\Delta A^{LP}_k(p,t)$ into concentration changes $\Delta C^{auto}_n(p,t)$, $n$ designates the two sources. The last part of the flowchart is the identification of functional brain areas using statistical parametric brain mapping. A linear general modeling leads to the computation of matrices $\beta$ and $Z_{stats}$ for the contrast $c$ ($\Delta C_{HbO_2}$, $\Delta C_{Hb}$, $\Delta C^{auto}_{1}$ or $\Delta C^{auto}_{2}$). $\beta$ represents the parameters of the general linear model that aims to compute the $Z_{stats}$ matrix. Finally, this last matrix was thresholded with the random field theory or the automatic thresholding procedure to identify the functional brain areas. We obtained the binary mask of statistical inferences $SPM_c$.
To estimate the slow hemodynamic fluctuations due to cerebral activity, the modified Beer-Lambert law and the auto-calibration method were applied in parallel, see section 3. The modified Beer-Lambert law converted the low-pass filtered absorbance changes into quantitative concentration changes of $HbO_2$ and $Hb$, see section 1.
The auto-calibration method was applied to the unfiltered temporal intensity signal, averaged over the surgical window, see section 2. This leads to the calculation of matrix $W$ (estimate of matrix $E$). The auto-calibration method converted the low-pass filtered absorbance changes into semi-quantitative contrasts.
A general linear model was separately executed on quantitative an semi-quantitative contrasts by testing the linear association between measurements and theoretical contrasts representing the hemodynamic response following a physiological stimulus, see section 4. Finally, statistical inferences (functional brain areas) were obtained with the random field theory or with the automatic thresholding procedure.
2.2.1 Intraoperative procedure
The study was conducted at the neurological center of the Pierre Wertheimer hospital in Bron, France. Twelve patients presenting a tumor close to the motor cortex area was included in the study. The experiment was approved by the local ethics committee of Lyon University Hospital (France) and the participating patients signed written consent.
After the patient’s craniotomy and before the brain tumor resection surgery, a video was acquired with a wide field optical device, see Fig. 2(A). Four different devices were used depending on the patient, see Table 1. In all cases, the devices were composed of an RGB camera in conjunction with an a continuous wave white light source. The white light source was illuminating the patient exposed cortex and the RGB camera was collecting the reflected light. After data acquisition, videos were processed by a laptop (processor: Intel Core i5-7200U, 2.50GHz $\times$ 4, ram: 15.3GiB).
Fig. 2. A - Schematic of the experimental setup. B - Absorbance changes ($\Delta A$) spectra measured by the optical device 1 (see Table 1) for patient 5 (see Table 2). The blue and green areas highlight the frequencies related to heartbeat and respiratory, respectively. These frequencies can be used as training data by the auto-calibration method. The red areas highlight the frequencies related to cerebral activity ($F\leq 0.05$ Hz).
Table 1. Intraoperative optical setups.
Table 2. Information on patients and acquisitions.
For all setups, videos were acquired with successive periods of rest and motor or sensory stimulation that could be repeated over time. Motor cortex stimulation was performed by the patient or by an external person by repetitive, alternating hand opening and closing at $\approx 1$ Hz. Stimulation of the sensory cortex was performed through repetitive fingers and palm caresses at $\approx 1$ Hz. These caresses were performed by an external person. For each patient, different number of stimulation cycles were performed. Information on patients and acquisitions is summarized in Table 2. For patients 7 and 8, optical setups 2 and 3 acquired videos in the same time. For patient 10, the 3 videos were acquired in a row.
The neurosurgeon performed electrical brain stimulation (EBS) after RGB imaging using a bipolar electrode (Nimbus Medtronic neurostimulator) to identify the functional brain areas during the neurosurgery. A biphasic current was used (pulsating frequency: $60$ Hz, pulse width: $1$ ms). The current was first set to $1$ mA, and increased to $6$ mA. When a functional area was identified by EBS, the neurosurgeon placed a colored pastille on the patient cortex and an RGB image was acquired to store the position of the functional area in the RGB image.
2.2.2 Pre-processing
For each frame of the RGB video, the repetitive brain motion was compensated [54,55]. The slow drift of RGB intensities was corrected due to tissue desiccation [40] and a low-pass filtering was performed to isolate slow hemodynamic fluctuations. The low-pass filtering was operated in the Fourier domain by multiplying the Fourier transform of each intensity time curve by a Blackman window (cut-off frequency: $0.05$ Hz, see the red rectangle in Fig. 2(B). Then, the filtered time curves were obtained by inverse Fourier transform. Finally, absorbance changes curves were computed for each color channel $k$ with low-pass filtered intensity curves:
with $I^{LP}_k(t)$ the filtered reflectance intensity measured at time $t$ by the camera color channel $k$ (red, green or blue). $I^{LP,\text {ref}}_k$ is the reference reflectance intensity measured by the camera color channel $k$ (average of the reflectance intensity over the duration of the first rest step of the experimental paradigm, see section 1).2.2.3 Modified Beer-Lambert law and auto-calibration method
Quantitative and semi-quantitative contrasts were computed with the modified Beer-Lambert law (see section 1) and the auto-calibration method (see section 2), respectively.
Quantitative contrasts $\Delta C_{HbO_2}$ and $\Delta C_{Hb}$ (in $\mu mol.L^{-1}$) were computed for each camera pixel with the modified Beer-Lambert law, see section 1. For the calculation of matrix E, wavelength-dependent parameters $L$, $D$, $S$ and $\epsilon$ need to be known, see Eq. (6). Values of the optical mean path length $L$ were taken from our previous studies [21,52]. The values of $L$ and the details regarding its estimation can be found in Ref. [21]. The molar extinction coefficients $\epsilon _{HbO_2}$ and $\epsilon _{Hb}$ of $HbO_2$ and $Hb$ were taken from Matcher et al. [56]. For the optical device 1 (see table 1), the spectral sensitivities of the detector $D$ were provided by the camera manufacturer and the light source spectra $S$ has been measured with a spectrometer. As we did not know the spectral characteristics of the cameras and white light sources of devices 2, 3 and 4 (see Table 1), the modified Beer-Lambert law could not be properly applied for videos 6, 7, 8.a, 8.b, 11 and 12, see Table 2. For these videos, we used the matrix E (see Eq. (6)) computed for the optical device 1.
Semi-quantitative contrasts $\Delta C_1^{auto}$ and $\Delta C_2^{auto}$ (in arbitrary units) were computed for each camera pixel with the auto-calibration method, see section 2. The first step was to compute unfiltered temporal intensity signal, averaged over the surgical window $I_k(t)$, see Fig. 1. Then, using this signal, absorbance changes were computed:
Spectra of absorbance changes $|TF(\Delta A_k)|$ were calculated, specific frequency bands were selected and separable NMF was executed, see section 2. This leads to the calculation of the matrix $W$, which is an estimate of the matrix $E$. To calculate $W$, we only used measurements, no prior knowledge were included into the model. Finally, low-pass filtered absorbance changes $\Delta A^{LP}$ were converted into semi-quantitative contrasts $\Delta C^{auto}_1$ and $\Delta C^{auto}_2$, see Eq. (10). As we mentioned in section 2, $\Delta A$ fluctuations related to heartbeat and respiratory may be considered in the auto-calibration method. Depending on the signal to noise ratio (SNR) of the acquisitions, we chose frequencies related to heartbeat ($0.95$ Hz $\leq f_H \leq 1.05$ Hz) and/or respiratory ($0.15$ Hz $\leq f_H \leq 0.25$ Hz). Note that these frequency range may change depending on the patient. In the rest of the manuscript, under- or superscripts $_H$, $_R$ or $_{H,R}$ indicate that the auto-calibration method was performed by taking into account heartbeat, respiratory or the concatenation of heartbeat and respiratory frequencies, respectively. For the acquisitions where the two periodic events were distinctly observable on the absorbance change spectra, the auto-calibration method was applied three times. For these cases, we obtained matrices $W_H$, $W_R$ and $W_{H,R}$.
2.2.4 Statistical parametric functional brain mapping
To identify the functional brain areas from the computed contrasts, the statistical parametric mapping (SPM) technique [36] was implemented. This technique was introduced by Karl Friston to process fMRI data. This method was also adapted for functional near infrared applications [37,57]. The basic idea was to test for each camera pixel the linear association between measured contrasts and a theoretical contrast that describes the patient hemodynamic response to a physiological stimulus. This theoretical contrast was obtained by convolving the hemodynamic impulse response function [58] to the window function that represents the patients physiological stimulus. A general linear model was implemented to evaluate the linear association between measured and expected responses. This leads to the computation of a matrix of $Z$ statistics for the contrast $c$, indicated $Z_{stat}(c)$ in the rest of the manuscript.
The random field theory (RFT) [37,59,60] was used to to threshold the $Z_{stat}$ matrices at $5{\% }$ statistical significance level with a correction for multiple comparisons (family-wise error) at the pixel level. In some cases, measured and expected hemodynamic responses did not correspond because the actual hemodynamic impulse response in the patient’s tissue was different from the one used in this paper. This leads to a lack of detection of functional brain areas due to a too high threshold of the $Z_{stat}$ matrix using RFT. In these cases, functional brain areas were identified as portions of patient’s brain where the $Z_{stat}>\mu (Z_{stat}) + 0.75 \sigma (Z_{stat})$ [22], with $\mu$ and $\sigma$ the mean and standard deviation functions, respectively. The binary image obtained after the thresholding operations is indicated $SPM_c$ in the rest of the paper. More details on the computation of the functional brain maps $Z_{stats}$ and $SPM$ can be found in our previous study [35] and in Supplement 1.
Data were processed with our laptop (processor: Intel Core i5-7200U, 2.50GHz $\times$ 4, ram: 15.3GiB) with a C++ software developed with Qt (v5.15.8), OpenCV (v4.6.0) [61] and FFTW (v3.3.10) [62]. For the $100$ s video of patient $5$, the brain motion was compensated in real time (during data loading). Filtering operations were executed in $45$ s after data loading. The auto-calibration method was computed in real time, contrasts were calculated in $8$ s and functional maps were obtained in $8$ s.
2.2.5 Evaluation metrics
For each acquisition, the matrices $W$ ($W_H$, $W_R$ or $W_{H,R}$) computed with the auto-calibration method (see sections 2 and 3) were compared to the matrix $E$ (see section 1) with the normalized cross covariance $NCC$:
We also used the $NCC$ to compare the $Z_{stats}$ matrices computed with the modified Beer-Lambert law with those obtained with the auto-calibration method. For each acquisition, we compared $Z_{stats}(HbO_2)$ to $Z_{stats}(C_1)$ and $Z_{stats}(Hb)$ to $Z_{stats}(C_2)$.
The DICE coefficient [64] was computed between the binary functional maps obtained with MBLL ($X$ in Eq. (14)) and the auto-calibration method ($Y$ in Eq. (14)):
Finally, we tested if the functional brain areas identified by electrical brain stimulation corresponded to the identifications provided by optical functional maps by testing if the electrical brain stimulation results were included in these maps.
3. Results
In Fig. 3, we represented maps of quantitative and semi-quantitative contrasts obtained for patient 5 at time $t=40$ s with the modified Beer-Lambert law and the auto-calibration method, respectively. Maps (a) and (b) corresponded to $\Delta C_{HbO_2}(t=40s)$ and $\Delta C_{Hb}(t=40s)$, respectively. Maps (c) and (d) corresponded to $\Delta C_{1}^{auto}(t=40s)$ and $\Delta C_{2}^{auto}(t=40s)$, respectively. The white dot and the letter M indicate the location of the patient’s motor area that has been identified with the EBS. Temporal quantitative and semi-quantitative contrasts measured at the level of the motor area are plotted at the right side of the figure. Temporal evolution of $\Delta C_{HbO_2}$ and $\Delta C_{Hb}$, $\Delta C_{1}^{auto}$ and $\Delta C_{2}^{auto}$ maps can be observed in Visualization 1.
Fig. 3. Quantitative and semi-quantitative contrasts obtained for patient 5 at time $t=40$ s with the modified Beer-Lambert law and the auto-calibration method, respectively. Maps (a) and (b) corresponded to $\Delta C_{HbO_2}(t=40s)$ and $\Delta C_{Hb}(t=40s)$, respectively. Maps (c) and (d) corresponded to $\Delta C_{1}^{auto}(t=40s)$ and $\Delta C_{2}^{auto}(t=40s)$, respectively. The white dot and the letter M indicate the location of the patient’s motor area that has been identified with the EBS. Temporal quantitative and semi-quantitative contrasts measured at the level of the motor area are plotted at the right side of the figure. On these graphs, solid curves represent quantitative or semi-quantitative contrasts, dashed black curves, the expected hemodyamic response $HR$. The green vertical lines indicate the temporal index used to plot maps (a) to (d), and the blue rectangles indicate periods of patient activity (hand movement).
In Fig. 4, we computed for each acquisition (see Table 2) the NCC (see Eq. (13)) between the matrices $E$ obtained with the modified Beer-Lambert law and matrices $W$ calculated with the auto-calibration method. $NCC(E,W_H)$, $NCC(E,W_R)$ and $NCC(E,W_{H,R})$ values were represented by blue, orange and green bars, respectively. For some acquisitions, no bars are shown. This means that the calculation of $W_H$, $W_R$ or $W_{H,R}$ matrices could not be performed, as frequencies associated with heartbeat or respiratory were not visible in the training signal. For 3 out of 16 acquisitions, the heartbeat frequencies were not visible in the training signal. For 7 out of 16 acquisitions, respiratory event was not visible in te training signal. NCC values computed between matrices $E$ and $W_H$ were $0.66 \pm 0.12$. Those computed between $E$ and $W_R$ were $0.69 \pm 0.13$. When heartbeat and respiratory frequencies were visible in the training signal, $NCC(E,W_R)$ values were $9{\% }$ greater than $NCC(E,W_H)$. $NCC(E,W_R)$ values were in the same order of magnitude than $NCC(E,W_{H,R})$ values ($0.69 \pm 0.15$).
Fig. 4. Normalized cross covariance computed for each acquisition between the matrices $E$ obtained with the modified Beer-Lambert law (see section 1) and matrices $W$ calculated with the auto-calibration method (see section 2).
The functional brain identifications are represented in Figs. 5, 6 and 7. For each acquisition, functional maps ($Z_{stats}$ matrices, see section 4) were plotted in colormap jet on the first image of the video sequence. We plotted functional maps calculated with the modified Beer-Lambert law ($Z_{stats}(HbO_2)$ and $Z_{stats}(Hb)$) and the auto-calibration method ($Z_{stats}(C_1^H)$, $Z_{stats}(C_2^H)$, $Z_{stats}(C_1^R)$, $Z_{stats}(C_2^R)$). Functional brain areas identified by optical imaging with the random field theory and the automatic thresholding procedure were indicated with black and magenta contours, respectively. Those detected with EBS were indicated by white points and letters $M_{i,j}$ and $S_{i,j}$ for motor and sensory areas $j$ of patient $i$.
Fig. 6. Functional brain identifications provided by optics and EBS for acquisitions from 7 to 10.a.
Fig. 7. Functional brain identifications provided by optics and EBS for acquisitions from 10.b to 12.
We indicated the concordance between functional optics and electrical brain stimulation in Table 3. When comparisons could not be performed, cells were filled with the character “-” (heartbeat or respiratory frequencies not visible in the training signal). We checked if identifications provided by EBS (white spots) were included in those detected by functional optics (black and magenta contours), see Figs. 5, 6 and 7. If there was a match, the cell was filled with a Yes, otherwise a No was indicated. Super-scripts $^*$ and $^{**}$ indicate that the functional brain areas were identified with the automatic thresholding procedure and the random field theory, respectively (see section 4).
Table 3. Concordance between functional identifications provided by intraoperative functional brain mapping techniques. A match between functional optics and EBS is indicated by Yes$^{**}$ (identification with the random field theory) or Yes$^*$ (identification with the automatic thresholding procedure).
In Figs. 8 and 9, we compared the functional brain optics maps computed with the modified Beer-Lambert law with those obtained with the auto-calibration method. In both figures, blue, orange and green bars represent values computed when considering frequencies related with heartbeat, respiratory and the concatenation of heartbeat and respiratory frequencies. For some acquisitions, no bars are shown. This means that the calculation of $W_H$, $W_R$ or $W_{H,R}$ matrices could not be performed, as frequencies associated with heartbeat or respiratory were not visible in the training signal.
Fig. 8. Comparison between $Z_{stats}$ matrices obtained with the modified Beer-Lambert law and the auto-calibration method, see colormaps in Figs. 5, 6 and 7.
Fig. 9. Comparison between functional brain areas identified by the modified Beer-Lambert law and the auto-calibration method, see black and magenta contours in Figs. 5, 6 and 7. Bars without hatches represent DICE coefficients computed with functional brain areas obtained with RFT. Hatched bars indicate that the automatic thresholding procedure was used to identify the functional brain areas.
In Fig. 8, we computed for each acquisition the NCC (see Eq. (13)) between $Z_{stats}(HbO_2)$ and $Z_{stats}(C_1)$ and $Z_{stats}(Hb)$ and $Z_{stats}(C_2)$ (colormaps in Figs. 5, 6 and 7). NCC values computed between matrices $Z_{stats}(HbO_2)$ and $Z_{stats}(C_1)$ were $0.69 \pm 0.30$ (for heartbeat frequencies) and $0.69 \pm 0.36$ (for respiratory frequencies). Those computed between matrices $Z_{stats}(Hb)$ and $Z_{stats}(C_2)$ were $0.37 \pm 0.41$ (for heartbeat frequencies) and $0.66 \pm 0.40$ (for respiratory frequencies). When heartbeat and respiratory frequencies were visible in the training signal, $NCC(Z_{stats}(HbO_2),Z_{stats}(C_1^H))$ ($0.79 \pm 0.16$) values were $8{\% }$ greater than $NCC(Z_{stats}(HbO_2),Z_{stats}(C_1^R))$ ($0.73 \pm 0.40$) but were $5{\% }$ lower than $NCC(Z_{stats}(HbO_2),Z_{stats}(C_1^{H,R}))$ ($0.83 \pm 0.17$). On the contrary, $NCC(Z_{stats}(Hb),$ $Z_{stats}(C_2^R))$ ($0.65 \pm 0.49$) values were $49{\% }$ greater than $NCC(Z_{stats}(Hb),Z_{stats}(C_2^H))$ ($0.43 \pm 0.42$) and were $12{\% }$ greater than $NCC(Z_{stats}(Hb),Z_{stats}(C_2^{H,R}))$ ($0.57 \pm 0.47$). $Z_{stats}(C_1)$ maps were more similar to $Z_{stats}(HbO_2)$ when using the concatenation of frequencies related to heartbeat and respiratory. $Z_{stats}(C_2)$ maps were more similar to $Z_{stats}(Hb)$ when using the frequencies associated with respiratory.
In Fig. 9, we computed for each acquisition the DICE coefficient (see Eq. (14)) between $SPM_{HbO_2}$ and $SPM_{C_1}$ and $SPM_{Hb}$ and $SPM_{C_2}$ (see black and magenta contours in Figs. 5, 6 and 7). Bars without hatches represent DICE coefficients computed with functional brain areas obtained with RFT. Hatched bars indicate that the automatic thresholding procedure was used to identify the functional brain areas.
When RFT was used to identify functional brain areas (black contours in Figs. 5, 6 and 7), DICE computed between $SPM_{HbO_2}$ and $SPM_{C_1}$ were $0.22 \pm 0.34$ (for heartbeat frequencies) and $0.25 \pm 0.35$ (for respiratory frequencies). Those computed between $SPM_{Hb}$ and $SPM_{C_2}$ were $0.08 \pm 0.19$ (for heartbeat frequencies) and $0.27 \pm 0.39$ (for respiratory frequencies). When heartbeat and respiratory frequencies were visible in the training signal, $DICE(SPM_{HbO_2},SPM_{C_1^R})$ ($0.30 \pm 0.35$) values were $61{\% }$ greater than $DICE(SPM_{HbO_2},SPM_{C_1^H})$ ($0.18 \pm 0.39$) but were $17{\% }$ lower than $DICE(SPM_{HbO_2},SPM_{C_1^{H,R}})$ ($0.36 \pm 0.39$). $DICE(SPM_{Hb},SPM_{C_2^R})$ values ($0.38 \pm 0.41$) were $121 {\% }$ greater than $DICE(SPM_{Hb},SPM_{C_2^H})$ ($0.17 \pm 0.24$) and $22 {\% }$ greater than $DICE(SPM_{Hb},SPM_{C_2^H})$ ($0.31 \pm 0.33$).
When the automatic thresholding procedure was used to identify functional brain areas (magenta contours in Figs. 5, 6 and 7), DICE computed between $SPM_{HbO_2}$ and $SPM_{C_1}$ were $0.51 \pm 0.35$ (for heartbeat frequencies) and $0.55 \pm 0.35$ (for respiratory frequencies). Those computed between $SPM_{Hb}$ and $SPM_{C_2}$ were $0.35 \pm 0.35$ (for heartbeat frequencies) and $0.60 \pm 0.32$ (for respiratory frequencies). When heartbeat and respiratory frequencies were visible in the training signal, $DICE(SPM_{HbO_2},SPM_{C_1^R})$ ($0.62 \pm 0.35$) values were $8{\% }$ greater than $DICE(SPM_{HbO_2},SPM_{C_1^H})$ ($0.57 \pm 0.32$) and were in the same order of magnitude than $DICE(SPM_{HbO_2},SPM_{C_1^{H,R}})$ ($0.62 \pm 0.36$). $DICE(SPM_{Hb},SPM_{C_2^R})$ ($0.60 \pm 0.33$) values were $26{\% }$ greater than $DICE(SPM_{Hb},SPM_{C_2^H})$ ($0.47 \pm 0.27$) and were $4{\% }$ greater than $DICE(SPM_{Hb},SPM_{C_2^{H,R}})$ ($0.57 \pm 0.31$).
$SPM_{C_1}$ maps were more similar to $SPM_{HbO_2}$ maps when using the concatenation of frequencies related to heartbeat and respiratory rather than just one of these frequency bands. $SPM_{C_2}$ maps were more similar to $SPM_{Hb}$ when using respiratory frequencies rather than heartbeat or the concatenation of heartbeat and respiratory frequencies.
4. Discussion
The auto-calibration method we proposed is based on the training of absorbance changes due to a spectral dependent periodic event. Using a spectral imaging device, our method allows to estimate endmembers in spectral bands where absorbance changes periodically. This technique has several advantages compared to traditional spectroscopic analyses used in diffuse optics, see section 1. First, the auto-calibration method does not require a complex calibration of the imaging setup. For instance, spectral sensitivities of the acquisition device and the light source do not need to be calibrated or even known. Then, the auto-calibration method does not require an estimation of the scattering and absorption properties of the media. We show that this method could be used with data collected from industrial RGB cameras or even smartphones.
The auto-calibration method is a fast and reliable algorithm. Two sources could be retrieved that are directly taken from the data. This allows an easy interpretation of the estimated spectra. In contrast, the estimated sources may not be the right ones, in case there are no columns of $X$ containing $W$, see Eq. (9). We applied the auto-calibration method for the identification of functional brain areas on $16$ acquisitions and saw that this method was capable of identifying functional brain areas, see Figs. 5, 6 and 7 and Table 3.
Although frequencies related to respiratory were not always present in the training signal, these frequencies seem to be particularly suitable for the auto-calibration method. Indeed, we saw in Fig. 4 that matrices $W$ were more similar to $E$ matrices when using the frequencies associated with respiratory ($F_R$) (or the concatenation of heartbeat and respiratory related frequencies $F_{H,R}$) rather than those related to heartbeat ($F_H$). The $NCC$ metric was only computed on $3\times 2$ matrices, which could certainly bias the result. A better comparison between $E$ and $W$ matrices could be operated with data collected by hyperspectral imaging [43,65]. This will allow us to acquire images with a higher spectral resolution and thus to compare the $E$ and $W$ matrices more efficiently. We also saw on Fig. 9 that functional brain areas identified the auto-calibration method were more similar to those obtained with the modified Beer-Lambert law when using $F_R$ and $F_{H,R}$ rather than only $F_H$. It is also interesting to notice that $\Delta C_{HbO_2}$ contrasts seem to be more easily estimated with $F_{H,R}$ or $F_H$ while $F_R$ seem to be preferred for estimating $\Delta C_{Hb}$ contrasts, see Fig. 8.
In a previous study [35], we show that $\Delta C_{Hb}$ contrasts are more suitable for detecting functional brain areas than $\Delta C_{HbO_2}$. Indeed, $\Delta C_{HbO_2}$ measurements are more subject to physiological noises, which does not allow, after filtering, to perfectly recover the hemodynamic response induced by neuronal activity. Since respiratory rate seems to be preferred for estimating $\Delta C_{Hb}$ contrasts, better functional identification is also obtained with these rates than with those related to heartbeat.
We can use the pulse oxymetry litterature to explain the relationship between semi-quantitative contrasts measured with the auto-calibration method ($\Delta C_1^{auto}$ and $\Delta C_2^{auto}$) and those measured with the modified Beer-Lambert law ($\Delta C_{HbO_2}$ and $\Delta C_{Hb}$). We also used this literature to discuss the choice of frequency bands used in the auto-calibration method (heartbeat- or respiration-centered frequencies) for estimating hemodynamic contrasts ($\Delta C_{HbO_2}$ and $\Delta C_{Hb}$). In pulse oxymetry applications, tissue oxygen saturation ($S_{tO_2}$) is measured by assessing arterial and local venous saturation values, denoted $S_{aO_2}$ and $S_{vO_2}$, respectively [66–68]. $S_{aO_2}$ provides an indication about the ventilation and the oxygen exchange in the lungs and $S_{vO_2}$ is a parameter that reflects the local balance between blood flow and oxygen consumption [67]. Near-infrared light has large depth penetration depth inside tissues ($\approx 15mm$ at $793$ nm [69]), which implies that the arterial, venous, and capillary compartments all contribute to the measured optical signal. On the contrary, light collected by RGB cameras has lower penetration depth than near-infrared light ($\approx 1mm$ at $532$ nm [69]) and contribution of arterial, venous, and capillary compartments can be more easily separated in the optical signal.
It has been shown that the contribution of the arterial compartment to absorbance measurements can be isolated [70]. Absorption changes in arteries are associated with the systolic-diastolic blood pressure variation at the heartbeat frequency [71]. This arterial blood is highly oxygenated, its normal saturation in oxygen lie between $95{\% }$ and $100{\% }$ [45]. Thus, after applying the auto-calibration method on heartbeat frequencies, it is possible to estimate a source related to $HbO_2$ since changes in $C_{HbO_2}$ are predominant in this signal. Furthermore, changes in $C_{HbO_2}$ and $C_{Hb}$ due to heartbeat can be considered as co-linear fluctuations since these variations are induced by a change in blood pressure. Thus, using frequencies related to heartbeat with near separable NMF, it is not possible to also have a good estimate of a second source related to $Hb$.
Wolf et al. [72] formulated the hypothesis that hemoglobin oscillation at the respiratory rate are mostly representative in veins. The compliance of the blood vessels, or the capacity of vessels to respond to an increase of blood pressure by distending and increasing its blood volume, is $\approx 20$ larger in veins than in arteries [67]. This means that for a given blood pressure, changes in blood volume will be $\approx 20$ times larger in veins than in arteries. Breathing is a mechanism that facilitates the return of blood from the brain back to the heart. This is also known as respiratory pump, and provoke hemoglobin oscillations at the respiratory frequency, with an increase and decrease of blood volume during inspiration and expiration, respectively. Since blood in veins return back to the heart, a change in $S_{vO_2}$ means a change in the use of oxygenated blood by the tissues [73]. Thus, changes in $C_{HbO_2}$ and $C_{Hb}$ may not be considered as co-linear fluctuations. For these reasons, the auto-calibration method may be more efficient to estimate oxy- and deoxygenated hemoglobin sources using the respiratory frequency. Thus $W$ matrices were more similar to $E$ matrices when using with respiratory rather than heartbeat frequencies. This may also explain the fact that the functional brain areas identified by the auto-calibration method are more similar to those obtained with the modified Beer-Lambert law when using frequencies associated with respiration rather than heart rate.
We plan to explore these hemodynamic oscillations in more detail in a future study. The dynamic model proposed by Sergio Fantini [74] can be used to simulate the temporal evolution of the concentration and oxygen saturation of hemoglobin in tissue with the considering time-dependent hemodynamic and metabolic parameters: blood volume, flow velocity, and oxygen consumption. Using this model, a second step could be to incorporate the changes in concentration and oxygen saturation of hemoglobin in a 3D numerical brain phantom [32] to estimate radiative quantities. Then, the auto-calibration method can be applied on the simulated optical signal to study in more detail the estimation of hemodynamic sources. We also plan to develop a liquid phantom, mimicking both the slow hemodynamic variations associated with brain activity and the periodic hemodynamic variations associated with heartbeat or breathing. The numerical and the liquid phantoms could be used to validate our method.
For 8 acquisitions out of 16 (videos $5$, $6$, $8.a$, $8.b$, $10.a$, $10.b$, $11$ and $12$), identifications obtained by the auto-calibration method were consistent with those provided by EBS (intraoperative gold standard). These videos have been acquired with four different optical setups: two neurosurgical microscopes, a smartphone and the experimental setup we used in our previous studies, see Table 1. For the remaining videos, the method may have been biased due to an important amount of noise in the training signal. The signal to noise ratio (SNR) was in average $2.4$ for videos with correct functional identifications and $1.5$ for videos with inconclusive identifications. The SNR was defined as the ratio between the standard deviation of the training signal computed over the heartbeat and/or respiratory frequency range and the standard deviation calculated over the frequency range which were not related to physiological events. A high level of noise in the training signal has a direct impact on the quality of the source estimation since Algorithm 1 finds directly the sources in data, estimated sources did not correspond to hemodynamic sources but to noise. Moreover, a large amount of noise may introduce outliers in data which will distort the estimation of the sources, which did not allow to identify functional brain areas.
To illustrate, the influence of noise on the estimation of the hemodynamic sources, we plotted the influence of the number of extreme points $P$ used in Algorithm 1 for estimating the matrix $E$ in the video $2$, see Fig. 10. Frequencies related with heartbeat $F_H$ were used to plot the blue curve, those related with respiratory $F_H$ to plot the orange curve and the concatenation of heartbeat and respiratory related frequencies ($F_{H,R}$) to plot the green curve. This explains the difference in length between the three curves. With $F_{H}$ and $F_{H,R}$, thee $NCC$ values increased strongly on the interval $P \in [1;5]$, then increase slightly until $P=\text {Card}(F_H)$ and $P=\text {Card}(F_{H,R})$, respectively. The learning signal seems to be particularly subject to the noise. The use of a number of extreme points $P \geq 5$ aims to increase the quality of the source estimation by reducing the noise impact. Since the variations in heartbeat signal are mainly due to $\Delta C_{HbO_2}$, the use of a greater number of extreme points leads to a better estimation of $HbO_2$ source. For $F_R$, the $NCC$ values oscillate between $0.65$. The number of extreme points has a low impact on the source estimation, since a low noise level is observed for these frequencies.
Fig. 10. Influence of the number of extreme points $P$ used in Algorithm 1 for estimating the matrix $E$ in the video $2$. The vertical red dotted line indicate the value of $P=5$ used in the manuscript.
For those videos with inconclusive identifications, $75{\% }$ (videos $1$ to $4$, $9$ and $10.c$) of the acquisitions were done by our experimental setup and $12.5{\% }$ (video $7$) by a Leica microscope and $12.5{\% }$ with a smartphone. The important noise amount in the collected signal is due to a poor SNR which could be due to multiple reasons. Making acquisitions in the operating room is often complex. The space is limited, and sometimes, it is complex to correctly place the acquisition system at an optimal distance and orientation from the patient’s brain. In many cases, the acquisition device must be moved as far back as possible from the patient, which limits the amount of light captured by the sensor and therefore deteriorates the SNR. The adjustments of the device must also be done quickly in order to minimize the duration of the procedure and the discomfort for the patient.
In order to make the method more robust, several improvements can be implemented. We plan to increase SNR in videos acquired with our experimental setup by connecting the camera to a side-port of the neurosurgical microscope, as proposed in the work of Pichette et al. [43]. Pre-processing steps can be implemented to maximise the heartbeat signal. For example, we can adapt work of Poh et al. who used independent component analysis to measure cardiac pulse beat [75]. In this application, a single training signal was used to retrieve the hemodynamic sources since intensity time curves were averaged over the brain surface. Thus, the heterogeneity of the tissue optical properties were not taken into consideration. Similarly, if a phase shift of heartbeat and respiratory signals are observed at different locations in the cerebral cortex, the average of the time signals could reduce or mask the frequency peaks in the training signal. A solution would be to use a tensor decomposition approach to estimate sources for each pixel independently [76–79]. In this case, the tensor would be an array of dimension $3$ (frequency, wavelength, pixel location).
The modified Beer-Lambert law has been applied on all videos, even for acquisitions with unknown spectral characteristics (videos $6$, $7.a$, $7.b$, $8.a$, $8.b$, $11$ and $12$). 12 out of 16 identifications were consistent with EBS. For the remaining 4 (videos $3$, $4$, $7.a$ and $7.b$), the modified Beer-Lambert coupled with SPM was not able to identify the functional brain areas. Measured hemodynamic responses did not correspond to theoretical responses, since low $Z_{stats}$ values were computed. These theoretical response functions are patient dependent and differ depending on the type of cortical tissue [80]. Moreover, the neurovascular system evolves with age, which implies a change in hemodynamic response [58] and the progression of gliomas over time implies a change in the hemodynamic response [81].
5. Conclusion
In conclusion, we proposed a new linear spectral unmixing method based on near separable NMF for the estimation of constituent spectra (sources) as well as their proportion in measurements: the auto-calibration method. This approach takes advantage of the periodic temporal variations of the sources in order to estimate them automatically and without device calibration. We applied the method for identifying functional brain areas during neurosurgery using RGB cameras (neurosurgical microscopes, industrial camera and smartphone). We have shown that the method can be used with industrial RGB cameras, but that it must be used under certain conditions. Data must be collected with high SNR and frequencies used by the training signal must be centered on the periodic event. Functional identification obtained with the auto-calibration method were consistent with those provided by the intraoperative gold standards. This work can allow a widespread use of spectral imaging in the industrial or medical field. With for example, the use of camera smartphones in telemedicine applications in order to help patients and doctors to carry out a continuous functional follow-up or a diagnostic help. The auto-calibration method could be applied during cardiac surgery to assess myocardial perfusion. It could also be used to monitor tissue reperfusion after heart, kidney or skin transplantation, or to monitor the patency of a vascular anastomosis (e.g. an intra-extra-cranial anastomosis). This work could also have an impact on the development of low-cost imaging devices for widespread use in developing countries.
Funding
LabEx PRIMES (ANR-11-LABX-0063); of Université de Lyon, within the program “Investissements d’Avenir” (ANR-11-IDEX-0007), operated by the French National Research Agency (ANR); Infrastructures d’Avenir en Biologie Santé (ANR-11-INBS-000), within the program “Investissements d’Avenir” operated by the French National Research Agency (ANR); France Life Imaging (ANR-11-INBS-0006); Pulsalys.
Acknowledgments
We want to acknowledge the PILoT facility for the support provided for the image acquisition.
Disclosures
The authors declare no conflicts of interest.
Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.
Supplemental document
See Supplement 1 for supporting content.
References
1. N. Keshava, “A survey of spectral unmixing algorithms,” Lincoln Laboratory Journal 14, 55–78 (2003).
2. C. Quintano and G. Pereira, “Spectral unmixing,” International Journal of Remote Sensing 33(17), 5307–5340 (2012). [CrossRef]
3. B. Hapke, Theory of Reflectance and Emittance Spectroscopy, (Cambridge University Press, 2012).
4. D. A. Roberts, M. O. Smith, and J. B. Adams, “Green vegetation, nonphotosynthetic vegetation, and soils in aviris data,” Remote Sensing of Environment 44(2-3), 255–269 (1993). [CrossRef]
5. M. José Bioucas-Dias, A. Plaza, N. Dobigeon, M. Parente, Q. Du, P. Gader, and J. Chanussot, “Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches,” IEEE J. Sel. Top. Appl. Earth Observations Remote Sensing 5(2), 354–379 (2012). [CrossRef]
6. R. L. Huguenin and J. L. Jones, “Intelligent information extraction from reflectance spectra: Absorption band positions,” J. Geophys. Res.: Solid Earth 91(B9), 9585–9598 (1986). [CrossRef]
7. F. Tsai and W. Philpot, “Derivative analysis of hyperspectral data,” Remote Sensing of Environment 66(1), 41–51 (1998). [CrossRef]
8. Q. Li, X. He, Y. Wang, H. Liu, D. Xu, and F. Guo, “Review of spectral imaging technology in biomedical engineering: achievements and challenges,” J. Biomed. Opt. 18(10), 100901 (2013). [CrossRef]
9. G. Lu and B. Fei, “Medical hyperspectral imaging: a review,” J. Biomed. Opt. 19(1), 010901 (2014). [CrossRef]
10. A. Rehman and S. A. Qureshi, “A review of the medical hyperspectral imaging systems and unmixing algorithms’ in biological tissues,” Photodiagn. Photodyn. Ther. 33, 102165 (2021). [CrossRef]
11. M. G. Sowa, L. Leonardi, J. R. Payette, J. S. Fish, and H. H. Mantsch, “Near infrared spectroscopic assessment of hemodynamic changes in the early post-burn period,” Burns 27(3), 241–249 (2001). [CrossRef]
12. M. Berman, A. Phatak, R. Lagerstrom, and B. R. Wood, “Ice: a new method for the multivariate curve resolution of hyperspectral images,” J. Chemom. 23(2), 101–116 (2009). [CrossRef]
13. B. Montcel, L. Mahieu-Williame, X. Armoiry, D. Meyronet, and J. Guyotat, “Two-peaked 5-ala-induced ppix fluorescence emission spectrum distinguishes glioblastomas from low grade gliomas and infiltrative component of glioblastomas,” Biomed. Opt. Express 4(4), 548–558 (2013). [CrossRef]
14. L. Alston, L. Mahieu-Williame, M. Hebert, P. Kantapareddy, D. Meyronet, D. Rousseau, J. Guyotat, and B. Montcel, “Spectral complexity of 5-ala induced ppix fluorescence in guided surgery: a clinical study towards the discrimination of healthy tissue and margin boundaries in high and low grade gliomas,” Biomed. Opt. Express 10(5), 2478–2492 (2019). [CrossRef]
15. G. Lu, X. Qin, D. Wang, Z. Chen, and B. Fei, “Estimation of tissue optical parameters with hyperspectral imaging and spectral unmixing,” In Medical Imaging 2015: Biomedical Applications in Molecular, Structural, and Functional Imaging, volume 9417, pages 199–205. SPIE, 2015.
16. S. Ogawa, T. M. Lee, A. R. Kay, and D. W. Tank, “Brain magnetic resonance imaging with contrast dependent on blood oxygenation,” Proc. Natl. Acad. Sci. 87(24), 9868–9872 (1990). [CrossRef]
17. I. J. Gerard, M. Kersten-Oertel, K. Petrecca, D. Sirhan, and J. A. Hall, “Brain shift in neuronavigation of brain tumors: A review,” Med. Image Anal. 35, 403–420 (2017). [CrossRef]
18. W. Penfield and E. Boldrey, “Somatic motor and sensory representation in the cerebral cortex of man as studied by electrical stimulation,” Brain 60(4), 389–443 (1937). [CrossRef]
19. J. Pallud, E. Mandonnet, R. Corns, E. Dezamis, E. Parraga, M. Zanello, and G. Spena, “Technical principles of direct bipolar electrostimulation for cortical and subcortical mapping in awake craniotomy,” Neurochirurgie 63(3), 158–163 (2017). [CrossRef]
20. C. Caredda, L. Mahieu-Williame, R. Sablong, M. Sdika, J. Guyotat, and B. Montcel, “Real time intraoperative functional brain mapping based on rgb imaging,” IRBM 42(3), 189–197 (2021). [CrossRef]
21. C. Caredda, L. Mahieu-Williame, S. Raphaíl, S. Michaíl, L. Alston, J. Guyotat, and B. Montcel, “Intraoperative quantitative functional brain mapping using an RGB camera,” Neurophotonics 6(04), 1–14 (2019). [CrossRef]
22. C. Caredda, L. Mahieu-Williame, Sa. Raphaíl, S. Michaíl, F. C. Schneider, J. Guyotat, and B. Montcel, “Intraoperative resting-state functional connectivity based on rgb imaging,” Diagnostics 11(11), 2067 (2021). [CrossRef]
23. B. Chance, Z. Zhuang, C. UnAh, C. Alter, and L. Lipton, “Cognition-activated low-frequency modulation of light absorption in human brain,” Proc. Natl. Acad. Sci. 90(8), 3770–3774 (1993). [CrossRef]
24. T. Meyer, S. B. Sobottka, M. Kirsch, G. Schackert, R. Steinmeier, E. Koch, and U. Morgenstern, “Intraoperative optical imaging of functional brain areas for improved image-guided surgery,” Biomed. Tech. 58(3), 225–236 (2013). [CrossRef]
25. K. A. Morone, J. S. Neimat, A. W. Roe, and R. M. Friedman, “Review of functional and clinical relevance of intrinsic signal optical imaging in human brain mapping,” Neurophotonics 4(3), 031220 (2017). [CrossRef]
26. M. Oelschlägel, M. Tobias, U. Morgenstern, H. Wahl, J. Gerber, G. Reiß, E. Koch, G. Steiner, M. Kirsch, G. Schackert, and S. B. Sobottka, “Mapping of language and motor function during awake neurosurgery with intraoperative optical imaging,” Neurosurgical Focus FOC 48(2), E3 (2020). [CrossRef]
27. M. Oelschlägel, W. H. Polanski, U. Morgenstern, G. Steiner, M. Kirsch, E. Koch, G. Schackert, and S. B. Sobottka, “Characterization of cortical hemodynamic changes following sensory, visual, and speech activation by intraoperative optical imaging utilizing phase-based evaluation methods,” Hum. Brain Mapp. 43(2), 598–615 (2022). [CrossRef]
28. W. H. Polanski, M. Oelschlägel, T. A. Juratli, Hannes Wahl, P. M. Krukowski, U. Morgenstern, E. Koch, G. Steiner, G. Schackert, and S. B. Sobottka, “Topographic mapping of the primary sensory cortex using intraoperative optical imaging and tactile irritation,” Brain Topogr. 36(1), 1–9 (2023). [CrossRef]
29. W.-L. Chen, J. Wagner, N. Heugel, J. Sugar, Y. W. Lee., L. Conant, M. Malloy, J. Heffernan, B. Quirk, A. Zinos, S. A. Beardsley, R. Prost, and H. T. Whelan, “Functional near-infrared spectroscopy and its clinical application in the field of neuroscience: Advances and future directions,” Front. Neurosci. 14, 1 (2020). [CrossRef]
30. M. Ferrari, L. Mottola, and V. Quaresima, “Principles, techniques, and limitations of near infrared spectroscopy,” Can. J. Appl. Physiol. 29(4), 463–487 (2004). [CrossRef]
31. S. Lloyd-Fox, A. Blasi, and C. Elwell, “Illuminating the developing brain: the past, present and future of functional near infrared spectroscopy,” Neurosci. Biobehav. Rev. 34(3), 269–284 (2010). [CrossRef]
32. C. Caredda, L. Mahieu-Williame, S. Raphaíl, S. Michaíl, J. Guyotat, and B. Montcel, “Optimal Spectral Combination of a Hyperspectral Camera for Intraoperative Hemodynamic and Metabolic Brain Mapping,” Appl. Sci. 10(15), 5158 (2020). [CrossRef]
33. N. Gillis, Nonnegative Matrix Factorization (Society for Industrial and Applied Mathematics, Philadelphia, PA, 2020).
34. J. Fandino, S. S. Kollias, H. G. Wieser, A. Valavanis, and Y. Yonekawa, “Intraoperative validation of functional magnetic resonance imaging and cortical reorganization patterns in patients with brain tumors involving the primary motor cortex,” J. Neurosurg. 91(2), 238–250 (1999). [CrossRef]
35. C. Caredda, L. Mahieu-Williame, S. Raphaíl, Sd. Michaíl, J. Guyotat, and B. Montcel, “Intraoperative Statistical Parametric Brain Mapping using RGB Imaging,” In European Conferences on Biomedical Optics 2021 (ECBO) (2021), paper ETu1C.7, page ETu1C.7. Optica Publishing Group, June 2021.
36. K. J. Friston, ed., Statistical parametric mapping: the analysis of funtional brain images, (Elsevier/Academic Press, Amsterdam ; Boston, 1st ed edition, 2007).
37. M. S. Hassanpour, B. R. White, A. T. Eggebrecht, S. L. Ferradal, A. Z. Snyder, and J. P. Culver, “Statistical analysis of high density diffuse optical tomography,” NeuroImage 85, 104–116 (2014). [CrossRef]
38. D. T. Delpy, M. Cope, P. van der Zee, S. Arridge, S. Wray, and J. S. Wyatt, “Estimation of optical pathlength through tissue from direct time of flight measurement,” Phys. Med. Biol. 33(12), 1433–1442 (1988). [CrossRef]
39. M. Kohl-Bareis, B. Ebert, J. P. Dreier, C. Leithner, U. Lindauer, and G. Royl, “Apparatus for measuring blood parameters,” United States of America Patent 20120277559, (WO2011070357A1), 11 2012.
40. M. Oelschlägel, T. Meyer, H. Wahl, S. B. Sobottka, H. Kirsch, G. Schackert, and U. Morgenstern, “Evaluation of intraoperative optical imaging analysis methods by phantom and patient measurements,” Biomed. Tech. 58(3), 1 (2013). [CrossRef]
41. G. Bale, C. E. Elwell, and T. Ilias, “From jöbsis to the present day: a review of clinical near-infrared spectroscopy measurements of cerebral cytochrome-c-oxidase,” J. Biomed. Opt. 21(9), 091307 (2016). [CrossRef]
42. M. Kohl, U. Lindauer, G. Royl, A. Villringer, and U. Dirnagl, “Physical model for the spectroscopic analysis of cortical intrinsic optical signals,” Phys. Med. Biol. 45(12), 3749–3764 (2000). [CrossRef]
43. J. Pichette, A. Laurence, L. Angulo, F. Lesage, A. Bouthillier, Dang Khoa Nguyen, and Frederic Leblond, “Intraoperative video-rate hemodynamic response assessment in human cortex using snapshot hyperspectral optical imaging,” Neurophotonics 3(04), 045003 (2016). [CrossRef]
44. C. Caredda, S. Michaíl, L. Mahieu-Williame, Sa. Raphaíl, and B. Montcel, “Method for processing images of a living biological tissue with auto-calibration,” FR2022050692, Oct. 2022.
45. Z. Kovacsova, G. Bale, and I. Tachtsidis, “Medical utility of nir monitoring,” Encyclopedia of Biomedical Engineering, Elsevier, pages 415–431, 2019.
46. C. Hartog and F. Bloos, “Venous oxygen saturation,” Bailliere’s Best Pract. Res., Clin. Anaesthesiol. 28(4), 419–428 (2014). [CrossRef]
47. H. Radhakrishnan, W. Vanduffel, H. P. Deng, L. Ekstrom, D. A. Boas, and M. A. Franceschini, “Fast optical signal not detected in awake behaving monkeys,” NeuroImage 45(2), 410–419 (2009). [CrossRef]
48. J. Steinbrink, F. C. D. Kempf, A. Villringer, and H. Obrig, “The fast optical signal–robust or elusive when non-invasively measured in the human adult?” NeuroImage 26(4), 996–1008 (2005). [CrossRef]
49. N. Gillis and S. A. Vavasis, “Fast and robust recursive algorithmsfor separable nonnegative matrix factorization,” IEEE Trans. Pattern Anal. Mach. Intell. 36(4), 698–714 (2013). [CrossRef]
50. N. Nadisic, N. Gillis, and C. Kervazo, “Smoothed separable nonnegative matrix factorization,” arXiv, arXiv:2110.05528, (2021). [CrossRef]
51. A. Marmoret and J. Cohen, “NN-FAC: Nonnegative factorization techniques toolbox,” 2020.
52. C. Caredda, L. Mahieu-Williame, R. Sablong, M. Sdika, J. Guyotat, and B. Montcel, “Pixel-wise modified Beer-Lambert model for intraoperative functional brain mapping,” In Q. J. Jones and T. G. van Leeuwen, eds., Clinical and Preclinical Optical Diagnostics II, volume 11073, pages 148–152. International Society for Optics and Photonics, SPIE (2019).
53. C. Caredda, L. Mahieu-Williame, S. Raphaíl, S. Michaíl, J. Guyotat, and B. Montcel, “Intraoperative functional and metabolic brain mapping using hyperspectral imaging,” In S. J. Madsen, V. X. D. Yang, and N. V. Thakor, eds., Clinical and Translational Neurophotonics 2020, volume 11225, pages 24–30. International Society for Optics and Photonics, SPIE, 2020.
54. M. Sdika, L. Alston, L. Mahieu-Williame, J. Guyotat, D. Rousseau, and B. Montcel, “Robust real time motion compensation for intraoperative video processing during neurosurgery,” In 2016 IEEE 13th International Symposium on Biomedical Imaging (ISBI), pages 1046–1049 (2016).
55. M. Sdika, L. Alston, D. Rousseau, J. Guyotat, L. Mahieu-Williame, and B. Montcel, “Repetitive motion compensation for real time intraoperative video processing,” Med. Image Anal. 53, 1–10 (2019). [CrossRef]
56. S. J. Matcher, C. E. Elwell, C. E. Cooper, M. Cope, and D. T. Delpy, “Performance comparison of several published tissue near-infrared spectroscopy algorithms,” Anal. Biochem. 227(1), 54–68 (1995). [CrossRef]
57. I. Tachtsidis, P. H. Koh, C. Stubbs, and C. E. Elwell, “Functional optical topography analysis using statistical parametric mapping (spm) methodology with and without physiological confounds,” In Ta. Eiji and D. F. Bruley, eds. Oxygen Transport to Tissue XXXI, pages 237–243, (Springer US, Boston, MA, 2010.
58. M. Veldsman, T. Cumming, and A. Brodtmann, “Beyond BOLD: Optimizing functional imaging in stroke populations: Optimizing BOLD Imaging in Stroke,” Hum. Brain Mapp. 36(4), 1620–1636 (2015). [CrossRef]
59. J. Cao and K. J. Worsley, Applications of Random Fields in Human Brain Mapping, pages 169–182, (Springer New York, New York, NY, 2001).
60. K. J. Friston, K. J. Worsley, R. S. J. Frackowiak, J. C. Mazziotta, and A. C. Evans, “Assessing the significance of focal activations using their spatial extent,” Hum. Brain Mapp. 1(3), 210–220 (1994). [CrossRef]
61. G. Bradski, “The OpenCV Library,” Dr. Dobb’s Journal of Software Tools (2000).
62. M. Frigo and S. G. Johnson, “The design and implementation of fftw3,” Proc. IEEE 93(2), 216–231 (2005). [CrossRef]
63. J. N. Sarvaiya, P. Suprava, and B. Salman, “Image registration by template matching using normalized cross-correlation,” In 2009 International Conference on Advances in Computing, Control, and Telecommunication Technologies, pages 819–822 (2009).
64. L. R. Dice, “Measures of the amount of ecologic association between species,” Ecology 26(3), 297–302 (1945). [CrossRef]
65. A. Laurence, H. Dè. Toffa, K. Peng, M. Robert, A. Bouthillier, D. K. Nguyen, and F. Leblond, “Multispectral intraoperative imaging for the detection of the hemodynamic response to interictal epileptiform discharges,” Biomed. Opt. Express 13(12), 6245 (2022). [CrossRef]
66. S. Fantini, F. Maria-Angela, J. S. Maier, S. A. Walker, B. B. Barbieri, and E. Gratton, “Frequency-domain multichannel optical detector for noninvasive tissue spectroscopy and oximetry,” Opt. Eng. 34(1), 32–42 (1995). [CrossRef]
67. M. A. Franceschini, D.A. Boas, A. Zourabian, S. G. Diamond, S. Nadgir, D. W. Lin, J. B. Moore, and S. Fantini, “Near-infrared spiroximetry: noninvasive measurements of venous saturation in piglets and human subjects,” J. Appl. Physiol. 92(1), 372–384 (2002). [CrossRef]
68. H. Liu, D. A. Boas, Y. Zhang, A. G. Yodh, and B. Chance, “Determination of optical properties and blood oxygenation in tissue using continuous NIR light,” Phys. Med. Biol. 40(11), 1983–1993 (1995). [CrossRef]
69. C. Caredda, “Imagerie optique spectrale : applications cliniques et précliniques,” PhD thesis, Université Claude Bernard Lyon 1, 2020, ThÊse de doctorat dirigée par Montcel, Bruno et Sablong, Raphaíl Ingénierie pour le vivant Lyon 2020.
70. J. W. Severinghaus, “Takuo Aoyagi: Discovery of Pulse Oximetry,” Anesth. Analg. 105(6), S1–S4 (2007). [CrossRef]
71. M. E. Wagshul, P. K. Eide, and J. R. Madsen, “The pulsating brain: a review of experimental and clinical studies of intracranial pulsatility,” Fluids Barriers CNS 8(1), 5–23 (2011). [CrossRef]
72. M. Wolf, G. Duc, M. Keel, P. Niederer, K. von Siebenthal, and B. Hans-Ulrich, “Continuous noninvasive measurement of cerebral arterial and venous oxygen saturation at the bedside in mechanically ventilated neonates,” Crit. Care Med. 25(9), 1579–1582 (1997). [CrossRef]
73. M. Nitzan, A. Babchenko, B. Khanokh, and H. Taitelbaum, “Measurement of oxygen saturation in venous blood by dynamic near IR spectroscopy,” J. Biomed. Opt. 5(2), 155–162 (2000). [CrossRef]
74. S. Fantini, “Dynamic model for the tissue concentration and oxygen saturation of hemoglobin in relation to blood volume, flow velocity, and oxygen consumption: Implications for functional neuroimaging and coherent hemodynamics spectroscopy (chs),” NeuroImage 85, 202–221 (2014). [CrossRef]
75. M.-Z. Poh, D. J. McDuff, and R. W. Picard, “Non-contact, automated cardiac pulse measurements using video imaging and blind source separation,” Opt. Express 18(10), 10762–10774 (2010). [CrossRef]
76. M. Boizard, R. Boyer, F. Gérard, J. E. Cohen, and P. Comon, “Performance estimation for tensor cp decomposition with structured factors,” In 2015 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 3482–3486. (IEEE, 2015).
77. J. Cohen, R. C. Farias, and P. Comon, “Fast decomposition of large nonnegative tensors,” IEEE Signal Process. Lett. 22(7), 862–866 (2014). [CrossRef]
78. J. E. Cohen and R. Bro, “Nonnegative parafac2: A flexible coupling approach,” In International Conference on Latent Variable Analysis and Signal Separation, pages 89–98. Springer, 2018.
79. M. A. Veganzones, J. E. Cohen, R. C. Farias, J. Chanussot, and P. Comon, “Nonnegative tensor CP decomposition of hyperspectral data,” IEEE Transactions on Geoscience and Remote Sensing 54(5), 2577–2588 (2016). [CrossRef]
80. M. Bruyns-Haylett, Y. Zheng, J. Berwick, and M. Jones, “Temporal coupling between stimulus-evoked neural activity and hemodynamic responses from individual cortical columns,” Phys. Med. Biol. 55(8), 2203–2219 (2010). [CrossRef]
81. M. K. Montgomery, S. H. Kim, A. Dovas, H. T. Zhao, A. R. Goldberg, W. Xu, A. J. Yagielski, M. K. Cambareri, K. B. Patel, A. Mela, N. Humala, D. N. Thibodeaux, M. A. Shaik, Y. Ma, J. Grinband, D. S. Chow, C. Schevon, P. Canoll, and E. M. C. Hillman, “Glioma-induced alterations in neuronal activity and neurovascular coupling during disease progression,” Cell Rep. 31(2), 107500 (2020). [CrossRef]
