Artifact Detection and BCI-Fit Artifact Report (#336)

Author: tab-cmd

bcipy/helpers/demo/demo_report.py

@@ -0,0 +1,131 @@
+from pathlib import Path
+from bcipy.helpers.load import load_json_parameters, load_raw_data, load_experimental_data
+from bcipy.helpers.triggers import trigger_decoder, TriggerType
+from bcipy.config import (
+    BCIPY_ROOT,
+    DEFAULT_PARAMETER_FILENAME,
+    RAW_DATA_FILENAME,
+    TRIGGER_FILENAME,
+    DEFAULT_DEVICE_SPEC_FILENAME)
+
+from bcipy.acquisition import devices
+from bcipy.helpers.acquisition import analysis_channels
+from bcipy.helpers.visualization import visualize_erp
+from bcipy.signal.process import get_default_transform
+from bcipy.signal.evaluate.artifact import ArtifactDetection
+from bcipy.helpers.report import Report, SignalReportSection, SessionReportSection
+
+
+if __name__ == "__main__":
+    import argparse
+
+    parser = argparse.ArgumentParser()
+    parser.add_argument(
+        '-p',
+        '--path',
+        help='Path to the directory with >= 1 sessions to be analyzed for artifacts',
+        required=False)
+
+    args = parser.parse_args()
+    colabel = True
+    # if no path is provided, prompt for one using a GUI
+    path = args.path
+    if not path:
+        path = load_experimental_data()
+
+    trial_window = (0, 1.0)
+
+    positions = None
+    for session in Path(path).iterdir():
+        # loop through the sessions, pausing after each one to allow for manual stopping
+        if session.is_dir():
+            print(f'Processing {session}')
+            prompt = input('Hit enter to continue or type "skip" to skip processing: ')
+            if prompt != 'skip':
+                # load the parameters from the data directory
+                parameters = load_json_parameters(
+                    f'{session}/{DEFAULT_PARAMETER_FILENAME}', value_cast=True)
+
+                # load the raw data from the data directory
+                raw_data = load_raw_data(Path(session, f'{RAW_DATA_FILENAME}.csv'))
+                type_amp = raw_data.daq_type
+                channels = raw_data.channels
+                sample_rate = raw_data.sample_rate
+                downsample_rate = parameters.get("down_sampling_rate")
+                notch_filter = parameters.get("notch_filter_frequency")
+                filter_high = parameters.get("filter_high")
+                filter_low = parameters.get("filter_low")
+                filter_order = parameters.get("filter_order")
+                static_offset = parameters.get("static_trigger_offset")
+
+                default_transform = get_default_transform(
+                    sample_rate_hz=sample_rate,
+                    notch_freq_hz=notch_filter,
+                    bandpass_low=filter_low,
+                    bandpass_high=filter_high,
+                    bandpass_order=filter_order,
+                    downsample_factor=downsample_rate,
+                )
+
+                # load the triggers
+                if colabel:
+                    trigger_type, trigger_timing, trigger_label = trigger_decoder(
+                        offset=parameters.get('static_trigger_offset'),
+                        trigger_path=f"{session}/{TRIGGER_FILENAME}",
+                        exclusion=[TriggerType.PREVIEW, TriggerType.EVENT, TriggerType.FIXATION],
+                    )
+                    triggers = (trigger_type, trigger_timing, trigger_label)
+                else:
+                    triggers = None
+
+                devices.load(Path(BCIPY_ROOT, DEFAULT_DEVICE_SPEC_FILENAME))
+                device_spec = devices.preconfigured_device(raw_data.daq_type)
+                channel_map = analysis_channels(channels, device_spec)
+
+                # check the device spec for any frontal channels to use for EOG detection
+                eye_channels = []
+                for channel in device_spec.channels:
+                    if 'F' in channel:
+                        eye_channels.append(channel)
+                if len(eye_channels) == 0:
+                    eye_channels = None
+
+                artifact_detector = ArtifactDetection(
+                    raw_data,
+                    parameters,
+                    device_spec,
+                    eye_channels=eye_channels,
+                    session_triggers=triggers)
+
+                detected = artifact_detector.detect_artifacts()
+                figure_handles = visualize_erp(
+                    raw_data,
+                    channel_map,
+                    trigger_timing,
+                    trigger_label,
+                    trial_window,
+                    transform=default_transform,
+                    plot_average=True,
+                    plot_topomaps=True,
+                )
+
+                # Try to find a pkl file in the session folder
+                pkl_file = None
+                for file in session.iterdir():
+                    if file.suffix == '.pkl':
+                        pkl_file = file
+                        break
+
+                if pkl_file:
+                    auc = pkl_file.stem.split('_')[-1]
+                else:
+                    auc = 'No Signal Model found in session folder'
+
+                sr = SignalReportSection(figure_handles, artifact_detector)
+                report = Report(session)
+                session = {'Label': 'Demo Session Report', 'AUC': auc}
+                session_text = SessionReportSection(session)
+                report.add(session_text)
+                report.add(sr)
+                report.compile()
+                report.save()

bcipy/helpers/report.py

@@ -1,6 +1,9 @@
+# mypy: disable-error-code="union-attr"
 import io
 from abc import ABC
-from typing import List, Optional
+from typing import List, Optional, Tuple
+
+from matplotlib import pyplot as plt
 
 from matplotlib.figure import Figure
 from reportlab.platypus import SimpleDocTemplate, Paragraph, Image
@@ -10,6 +13,7 @@
 from reportlab.lib.units import inch
 
 from bcipy.config import BCIPY_FULL_LOGO_PATH
+from bcipy.signal.evaluate.artifact import ArtifactDetection
 
 
 class ReportSection(ABC):
@@ -40,9 +44,15 @@ class SignalReportSection(ReportSection):
 
     def __init__(
             self,
-            figures: List[Figure]) -> None:
+            figures: List[Figure],
+            artifact: Optional[ArtifactDetection] = None) -> None:
         self.figures = figures
         self.report_flowables: List[Flowable] = []
+        self.artifact = artifact
+
+        if self.artifact:
+            assert self.artifact.analysis_done is not False, (
+                "If providing artifact for this report, an analysis must be complete to run this report.")
         self.style = getSampleStyleSheet()
 
     def compile(self) -> Flowable:
@@ -51,9 +61,65 @@ def compile(self) -> Flowable:
         Compiles the Signal Report sections into a flowable that can be used to generate a Report.
         """
         self.report_flowables.append(self._create_header())
+        if self.artifact:
+            self.report_flowables.append(self._create_artifact_section())
         self.report_flowables.extend(self._create_epochs_section())
+
         return KeepTogether(self.report_flowables)
 
+    def _create_artifact_section(self) -> Flowable:
+        """Create Artifact Section.
+
+        Creates a paragraph with the artifact information. This is only included if an artifact detection is provided.
+        """
+        artifact_report = []
+        artifacts_detected = self.artifact.dropped
+        artifact_text = '<b>Artifact:</b>'
+        artifact_section = Paragraph(artifact_text, self.style['BodyText'])
+        artifact_overview = f'<b>Artifacts Detected:</b> {artifacts_detected}'
+        artifact_section = Paragraph(artifact_overview, self.style['BodyText'])
+        artifact_report.append(artifact_section)
+
+        if self.artifact.eog_annotations:
+            eog_artifacts = f'<b>EOG Artifacts:</b> {len(self.artifact.eog_annotations)}'
+            eog_section = Paragraph(eog_artifacts, self.style['BodyText'])
+            artifact_report.append(eog_section)
+            heatmap = self._create_heatmap(
+                self.artifact.eog_annotations.onset,
+                (0, self.artifact.total_time),
+                'EOG')
+            artifact_report.append(heatmap)
+
+        if self.artifact.voltage_annotations:
+            voltage_artifacts = f'<b>Voltage Artifacts:</b> {len(self.artifact.voltage_annotations)}'
+            voltage_section = Paragraph(voltage_artifacts, self.style['BodyText'])
+            artifact_report.append(voltage_section)
+
+            # create a heatmap with the onset values of the voltage artifacts
+            onsets = self.artifact.voltage_annotations.onset
+            heatmap = self._create_heatmap(
+                onsets,
+                (0, self.artifact.total_time),
+                'Voltage')
+            artifact_report.append(heatmap)
+        return KeepTogether(artifact_report)
+
+    def _create_heatmap(self, onsets: List[float], range: Tuple[float, float], type: str) -> Image:
+        """Create Heatmap.
+
+        Creates a heatmap image with the onset values of the voltage artifacts.
+        """
+        # create a heatmap with the onset values
+        fig, ax = plt.subplots()
+        fig.set_size_inches(6, 3)
+        ax.hist(onsets, bins=100, range=range, color='red', alpha=0.7)
+        ax.set_title(f'{type} Artifact Onsets')
+        ax.set_xlabel('Time (s)')
+        # make the label text smaller
+        ax.set_ylabel('Frequency')
+        heatmap = self.convert_figure_to_image(fig)
+        return heatmap
+
     def _create_epochs_section(self) -> List[Image]:
         """Create Epochs Section.
 

bcipy/helpers/stimuli.py

@@ -412,21 +412,37 @@ def update_inquiry_timing(timing: List[List[float]], downsample: int) -> List[Li
 
 
 def mne_epochs(mne_data: RawArray,
-               trigger_timing: List[float],
                trial_length: float,
-               trigger_labels: List[int],
-               baseline: Optional[Tuple[Any, float]] = None) -> Epochs:
+               trigger_timing: Optional[List[float]] = None,
+               trigger_labels: Optional[List[int]] = None,
+               baseline: Optional[Tuple[Any, float]] = None,
+               reject_by_annotation: bool = False,
+               preload: bool = False) -> Epochs:
     """MNE Epochs.
 
     Using an MNE RawArray, reshape the data given trigger information. If two labels present [0, 1],
     each may be accessed by numbered order. Ex. first_class = epochs['1'], second_class = epochs['2']
     """
-    annotations = Annotations(trigger_timing, [trial_length] * len(trigger_timing), trigger_labels)
-    mne_data.set_annotations(annotations)
-    events_from_annot, _ = mne.events_from_annotations(mne_data)
-    if not baseline:
-        baseline = (None, 0.0)
-    return Epochs(mne_data, events_from_annot, tmax=trial_length, baseline=baseline)
+    old_annotations = mne_data.annotations
+    if trigger_timing and trigger_labels:
+        new_annotations = Annotations(trigger_timing, [trial_length] * len(trigger_timing), trigger_labels)
+        all_annotations = new_annotations + old_annotations
+    else:
+        all_annotations = old_annotations
+
+    tmp_data = mne_data.copy()
+    tmp_data.set_annotations(all_annotations)
+
+    events_from_annot, _ = mne.events_from_annotations(tmp_data)
+    return Epochs(
+        mne_data,
+        events_from_annot,
+        baseline=baseline,
+        tmax=trial_length,
+        tmin=-0.05,
+        proj=False,  # apply SSP projection to data. Defaults to True in Epochs.
+        reject_by_annotation=reject_by_annotation,
+        preload=preload)
 
 
 def alphabetize(stimuli: List[str]) -> List[str]:

bcipy/helpers/triggers.py

@@ -161,6 +161,7 @@ class TriggerType(Enum):
     OFFSET = "offset"
     EVENT = "event"
     PREVIEW = "preview"
+    ARTIFACT = "artifact"
 
     @classmethod
     def list(cls) -> List[str]:

bcipy/helpers/visualization.py

@@ -91,7 +91,7 @@ def visualize_erp(
         baseline = None
 
     mne_data = convert_to_mne(raw_data, channel_map=channel_map, transform=transform)
-    epochs = mne_epochs(mne_data, trigger_timing, trial_length, trigger_labels, baseline=baseline)
+    epochs = mne_epochs(mne_data, trial_length, trigger_timing, trigger_labels, baseline=baseline)
     # *Note* We assume, as described above, two trigger classes are defined for use in trigger_labels
     # (Nontarget=0 and Target=1). This will map into two corresponding MNE epochs whose indexing starts at 1.
     # Therefore, epochs['1'] == Nontarget and epochs['2'] == Target.

bcipy/parameters/devices.json

@@ -33,13 +33,7 @@
         "excluded_from_analysis": [
             "TRG",
             "X1", "X2", "X3",
-            "A2",
-            "T3", "T4",
-            "Fp1", "Fp2",
-            "F7", "F8",
-            "P3", "P4",
-            "F3", "F4",
-            "C3", "C4"
+            "A2"
         ],
         "status": "active",
         "static_offset": 0.1
@@ -67,7 +61,7 @@
         "name": "DSI-Flex",
         "content_type": "EEG",
         "channels": [
-            { "name": "F3", "label": "Fz", "units": "microvolts", "type": "EEG" },
+            { "name": "P4", "label": "Cz", "units": "microvolts", "type": "EEG" },
             { "name": "S2", "label": "Oz", "units": "microvolts", "type": "EEG" },
             { "name": "S3", "label": "P4", "units": "microvolts", "type": "EEG" },
             { "name": "S4", "label": "P3", "units": "microvolts", "type": "EEG" },