2.1.

Multispectral Imaging System

A custom multispectral imaging system connected to a neurosurgical microscope (OPMI Pentero, Zeiss) was used to image the exposed cortex of a patient during epilepsy surgery.¹² The system integrates a snapshot multispectral camera (IMEC, Leuven, Belgium) composed of $4\times 4$ arrays of band-pass filters disposed over a charged-coupled device (CCD) chip of $1024\times 2048\text{pixels}$, leading to raw images of 16 spectral bands (480 to 630 nm with a 15-nm spectral resolution) of $256\times 512\text{pixels}$. The camera was operated at 20 frames per second (fps) to allow sufficient signal-to-noise ratio (SNR) and recordings were performed for 8 min for a total of 9600 images. The white light source (Superlux, 300 W xenon lamp) integrated into the microscope was used to continuously illuminate the brain surface during the imaging session. The microscope working distance was $\sim 25$ cm and the total microscope magnification was set at $1.5\times$. The lenses in the camera adaptor led to a resulting magnification of 4.7, an imaging field of view of $\sim 13{\mathrm{cm}}^2$, and a spatial resolution $\sim 0.1\mathrm{mm}$. Electrocorticography (ECoG) electrodes were placed on the patient’s brain but electrophysiological data were not used as part of the work presented here. Informed consent was obtained from the patient and the Centre Hospitalier de l’Université de Montréal (CHUM) ethics review board approved the research protocol.

The dataset is composed of acquisitions performed during epilepsy surgery in three patients where regions of accumulated blood were noticeable in the field of view. The dataset of patient 1 exhibited a clear drop of blood and was chosen to design the segmentation technique described in Secs. 2.2 and 2.3. Patients 2 and 3 exhibited areas of accumulated blood and were selected to demonstrate the technique’s capabilities to segment blood, blood vessels, and cortex in different patients. The total number of acquisitions performed with the multispectral imaging system reached 15 patients, but only three presented that spontaneous accumulation of blood and were included in the study.

Calibration and biophysical modeling procedures were applied to the spatial–temporal–spectral patient datasets resulting in a time-course of oxygenated hemoglobin (HbO) and deoxygenated hemoglobin (HbR) concentration variations at the surface of the cortex. A detailed description of the algorithms and model parameters can be found in Ref. 12. Briefly, spectral and spatial calibrations were initially performed to account for the system response by normalization using a dataset measured on a reflectance standard (Spectralon, Labsphere). Pixels attaining the detector maximum intensity were associated with specular reflections and were excluded from the analysis. Then, a spatial registration algorithm was implemented using the MATLAB^® Medical Image Registration Toolbox.¹³ That procedure was applied to each image, resulting in a dataset where the mechanical motion of the brain due to breathing and heartbeat was removed. A modified Beer–Lambert law was applied to all images on a pixel-by-pixel basis to compute the relative concentration changes of HbO and HbR.¹⁴ The differential pathlength factor used in the modified Beer–Lambert law was computed using estimates from the literature for the absorption and scattering coefficients $({\mu }_a,{\mu }_s^{\prime })$ in brain.¹⁵ Parameters and equations used to model absorption and scattering coefficients are detailed in Sec. 2.3.1.

2.2.

Hemodynamic Activity Characteristics of Cortex, Blood Vessels, and Bleeding

Following the application of the preprocessing steps, the dataset consisted of videos (9600 images) of brain hemodynamic activity, i.e., time-sequences showing the concentration variations of HbO and HbR for each of the $256\times 512\mathrm{pixels}$. A Fourier transform was applied to the concentration variation of HbO within each pixel to obtain the frequency profile of the hemodynamic activity from 0 to 10 Hz with a resolution of 0.0024 Hz. To avoid confusion, in this paper, the term “temporal frequencies” refers to the frequency spectra resulting from the Fourier transform (Hz). Otherwise, the term “spectra” refers to the optical spectrum of the diffuse reflectance (nm). A standard normal variate (SNV) normalization was applied to temporal frequencies so that the relative distribution of components was independent of absolute power density.

The patient 1 dataset showed a blood drop within the imaging field of view that slowly accumulated in the top-left corner of the image and started dripping toward the bottom of the image. The time frames clearly exhibiting the blood drop (#4800 to 8000) were selected to evaluate the temporal and spectral characteristics of the three different components of the exposed brain: blood vessels, cortex, and bleeding. For patients 2 and 3, no apparent movement of blood was noticed and 4000 frames were selected for data processing. Manual segmentation of the three physiological components was performed within those time frames to identify some regions for which the nature of the imaged tissue could be ascertained. Bleeding was identified from areas of dense accumulated blood; blood vessels with clear boundaries including arteries and veins were selected; clear areas of cortex tissue were circled. The temporal frequency profiles of all pixels within each segmented area were averaged to visualize the temporal frequency profile of each component.

The temporal frequency interval of 0.03 to 0.3 Hz was identified as exhibiting large differences between blood vessels, bleeding, and cortex mainly due to the respiration rate that has a noticeable effect in HbO concentration trends and was set by the medical team to 0.23, 0.2, and 0.22 Hz for the three patients, respectively. Temporal frequency spectra were averaged over the 0.03 to 0.3 Hz range for each pixel of the segmented components and a Kruskal–Wallis one-way analysis of variance between regions was performed for the frequency range. A threshold value corresponding to the 75th percentile of the bleeding component was selected for each patient to determine the lower limit value of hemodynamic activity intensity, used in a semiautomatic segmentation algorithm.

2.3.

Segmentation Algorithm Using Reflectance Spectra and Hemodynamic Activity

Using information obtained from Sec. 2.2, a semiautomatic segmentation technique was developed to identify blood vessels, bleeding, and cortex in the image. Compared to existing methods detecting the presence of blood vessels based on the presence of hemoglobin absorption signature in acquired reflectance spectra, this new method combines the spectral detection of hemoglobin with the temporal frequency characteristics of HbO time-response. This allows more specific discrimination between dynamic blood (blood vessels supplied in oxygen) and static blood that has accumulated.

The segmentation technique is performed in three steps, detailed in the following sections. For step one, pixels with a reflectance spectrum indicating the presence of hemoglobin are identified as blood, including here bleeding and blood vessels. For step two, pixels with HbO temporal frequency profiles indicating strong low-frequency components are identified as dynamic (cortex, blood vessels). Other pixels are identified as static (bleeding, skull, and other structures). Step 3 consists of combining the information of steps 1 and 2 to categorize pixels into bleeding, blood vessels, or cortex. Each of the three segmentation steps is detailed in the following sections of the paper.

3.1.

Hemodynamic Activity of Cortex, Blood Vessels, and Bleeding

Figure 1 shows the manual segmentation [Fig. 1(a)] and average temporal frequency spectrum [Fig. 1(b)] for patient 1. Averaged HbO concentration variations for each region of the brain exhibited an inverse trend ($1/f$) with respect to the temporal frequency $f$, which has been observed to be associated with spontaneous hemodynamic activity,¹⁷ combined with sharper and distinctive peaks at specific frequencies [Fig 1(b)]. The peak at $\sim 0.23\mathrm{Hz}$ corresponds to the respiratory rate for this patient, which was set by the medical team at 14 breathing cycles per minute during the surgical procedure. This peak was present in the cortex and vessel regions, but it had noticeably smaller amplitude within the blood drop. The second harmonic of the respiratory rate is also visible in the temporal frequency spectrum at $\sim 0.46\mathrm{Hz}$, as well as the heartbeat at $\sim 1.3\mathrm{Hz}$ ($79\text{beats}/\min$), which are consistent with the values recorded by the physiological monitoring device during the procedure.

Fig. 1

Hemodynamic activity profiles of cortex, blood vessels, and bleeding: (a) manual segmentation of three visually distinct regions: cortex (green), blood vessels (blue), and bleeding (red). ECoG electrodes, not used in this study, can be observed in the top right corner of the image. (b) Average SNV temporal frequency spectrum of HbO concentration for each of the three regions. Identified peaks correspond to respiration rate (R1), its second harmonic (R2), and heartbeat (HB).

Whisker plots were produced to display the distribution of frequency intensity within all three segmented regions and for two frequency ranges exhibiting the most differences between regions: 0.03 to 0.3 Hz [Fig. 2(a)] associated with respiration rate and 1 to 1.7 Hz [Fig. 2(b)] associated with the heartbeat. Outside of these ranges ($f<0.03$ and $f>1.7\mathrm{Hz}$), data showed little differences in amplitude and less SNR. Amplitudes were noticeably lower for bleeding in all frequency ranges, except for the lowest frequencies (0 to 0.03 Hz) where the differences between regions were less apparent. A Kruskal–Wallis one-way analysis of variance between regions showed a $p$-value $<0.01$ for all frequency ranges. The temporal frequency amplitudes of pixels for cortex and blood vessels also appeared to be differentiable, with the largest observed differences in the 0.03 to 0.3 Hz range. The outliers in the whisker plots corresponded to pixels detecting high reflectivity (specular reflection) and were mostly associated with CCD saturation. Those pixels represented $<0.7\%$ of all data points in the 0.03 to 0.3 Hz range and $<4\%$ in the 1.0 to 1.7 Hz range. The performance of the brain motion correction algorithm was diminished for those high intensity pixels, which can explain the larger number of outliers in the frequency range containing the signal form heartbeat. For each patient, the bleeding region was the one with the lowest signal amplitude, and the cortex region appeared to have the highest signal amplitude. Average signal intensity varies from one patient to another, with patient 2 having the lowest. A threshold specific to each patient was then chosen to separate the bleeding component from cortex and blood vessels based on these intensity differences.

Fig. 2

Distribution of the intensity of the SNV-normalized temporal frequency HbO spectra of each pixel comprised within the segmented components of bleeding (B), cortex (C), and blood vessels (V) for each patient (labeled 1, 2, and 3): (a) 0.03 to 0.3 Hz and (b) 1.0 to 1.7 Hz.

3.2.

Segmentation of Blood Vessels, Bleeding, and Cortex

Figure 3(a) shows the resulting segmented regions with hemodynamic activity (in black) and without activity (in light red) for patient 1. The blood drop at the image top in Fig. 3(c) and the blood accumulation toward the image bottom are identified as inactive by the algorithm and thus correspond to bleeding. Figure 3(b) displays the blood (light blue) and cortex (black) segmentation using the reflectance spectra-based segmentation algorithm. Both bleeding and blood vessels are identified as blood components by the algorithm. Figure 3(d) shows the final segmentation step, combining the results of Figs. 3(a) and 3(b). Bleeding is presented in red, blood vessels are in blue, and cortex in green. Gray pixels indicate saturated pixels due to specular reflections and black pixels were not classified in any category.

Fig. 3

Semiautomatic segmentation of bleeding, blood vessels, and brain cortex in patient 1: (a) segmentation of regions with hemodynamic activity (black) and regions without hemodynamic activity (light red) using a method based on the temporal frequency of reflectance signals; (b) segmentation of cortex (black) and blood (light blue) using a method based on the reflectance spectrum; (c) color visualization of the region of interest; and (d) bleeding, blood vessels, and cortex segmentation using logical intersections of images in (a) and (b). Gray pixels indicate specular reflections.

Color visualization of the region of interest is presented for patients 2 [Fig. 4(a)] and 3 [Fig 4(b)], with the final segmentation results of blood vessels, cortex, and bleeding [Figs. 4(c)–4(d)]. Areas of blood pooling are identified as bleeding by the algorithm. Blood vessels are identified for patient 2, but only one deep blood vessel is present for patient 3’s region of interest and was not detected as a blood vessel by the algorithm. The areas corresponding to the plastic of the electrodes are associated with the largest errors in terms of performance for the registration algorithm used in the preprocessing phase. This implies that residues of the movement of the brain due to heartbeat and respiration rate appear in HbO time sequences and thus the areas are classified as “dynamic” by the algorithm.

Fig. 4

Bleeding, blood vessels, and cortex segmentation on two patients: (a), (b) color visualization of the region of interest for patient 2 and patient 3; (c), (d) bleeding, blood vessels, and cortex segmentation. Gray pixels indicate specular reflections.

3.3.

Minimal Acquisition Time Analysis

The analysis presented in Fig. 3 for patient 1 was performed using 3200 frames acquired at 20 fps, corresponding to a total acquisition time of 160 s. However, such large acquisition times are not practical for surgical applications, which typically require much faster identification of vessels or blood pooling in a matter of seconds if not milliseconds. Considering that the lowest frequency within the interval displaying the most discriminatory potential (0.03 Hz) has a period of 33.3 s, the technique theoretically requires an acquisition time of around 30 s. An evaluation of the minimal acquisition time necessary to segment the image using the 0.03 to 0.3 Hz frequency range was performed by extracting sections of the dataset with a different number of frames and performing the segmentation algorithm on every section. Sections were formed by selecting $n$ consecutive frames, with $n$ ranging from 100 to 3600 by steps of 100 frames. The results were compared to reference values obtained with the full acquisition, i.e., with 3200 frames. The average temporal frequency spectrum intensities within the 0.03 to 0.3 Hz range were similar ($<10\%$ variation) to the reference values for acquisition times associated with $>500\text{frames}$ (25 s), and variations increased to $>20\%$ for an acquisition of 200 frames (10 s). The two-step automatic segmentation algorithm for these limited time-sequences further confirmed significant degradation of segmentation results at 500 frames.