Electrophysiology chunker MEG improvements

python/loris_ephys_chunker/src/loris_ephys_chunker/chunking.py

@@ -2,7 +2,6 @@
 import math
 import sys
 from collections import OrderedDict
-from collections.abc import Callable
 from pathlib import Path
 from typing import Any, cast
 
@@ -159,8 +158,8 @@ def write_chunks(chunk_dir: Path, channel_chunks_list: list[ChannelArray], chann
 
 def mne_file_to_chunks(
     path: Path,
+    raw: BaseRaw,
     chunk_size: int,
-    loader: Callable[[Path], BaseRaw],
     from_channel_name: str | None,
     channel_count: int | None,
 ) -> tuple[
@@ -171,9 +170,8 @@ def mne_file_to_chunks(
     list[tuple[float, float]],
     list[int],
 ]:
-    parsed = loader(path)
-    time_interval: tuple[np.float64, np.float64] = (parsed.times[0], parsed.times[-1])
-    channel_names = cast(list[str], parsed.info["ch_names"])
+    time_interval: tuple[np.float64, np.float64] = (raw.times[0], raw.times[-1])
+    channel_names = cast(list[str], raw.info["ch_names"])
     channel_ranges: list[tuple[float, float]] = []
     signal_range = (np.inf, -np.inf)
     channel_chunks_list = []
@@ -189,7 +187,7 @@ def mne_file_to_chunks(
 
     for i, channel_name in enumerate(selected_channels, start=1):
         print(f"Processing channel {channel_name} ({i} / {len(selected_channels)})")
-        channel = cast(ChannelArray, parsed.get_data(channel_name))  # type: ignore
+        channel = cast(ChannelArray, raw.get_data(channel_name))  # type: ignore
         channel_min = np.amin(channel)
         channel_max = np.amax(channel)
         channel_ranges.append((channel_min, channel_max))
@@ -215,8 +213,8 @@ def mne_file_to_chunks(
 
 def write_chunk_directory(
     path: Path,
+    raw: BaseRaw,
     chunk_size: int,
-    loader: Callable[[Path], BaseRaw],
     from_channel_index: int = 0,
     from_channel_name: str | None = None,
     channel_count: int | None = None,
@@ -227,7 +225,7 @@ def write_chunk_directory(
     chunk_dir = chunk_dir_path(path, prefix=prefix, destination=destination)
     channel_chunks_list, time_interval, signal_range, \
         channel_names, channel_ranges, valid_samples_in_last_chunk = \
-        mne_file_to_chunks(path, chunk_size, loader, from_channel_name, channel_count)
+        mne_file_to_chunks(path, raw, chunk_size, from_channel_name, channel_count)
 
     if downsamplings is not None:
         channel_chunks_list = channel_chunks_list[:downsamplings]

python/loris_ephys_chunker/src/loris_ephys_chunker/scripts/ctf_to_chunks.py

@@ -11,16 +11,29 @@
 from loris_ephys_chunker.chunking import write_chunk_directory  # type: ignore
 
 
-def load_channels(path: Path) -> RawCTF:
-    return mne.io.read_raw_ctf(  # type: ignore
+def load_ctf_raw(path: Path) -> RawCTF:
+    """
+    Read the CTF acquisition file into an MNE raw object.
+    """
+
+    raw = mne.io.read_raw_ctf(  # type: ignore
         path,
-        preload=False,
         # CTF raw channel names can contain suffixes that causes them to mismatch their
         # corresponding `channels.tsv` entries, the following flag removes these suffixes.
         clean_names=True,
         verbose=False,
     )
 
+    # Apply third-order software gradient compensation to remove environmental noise.
+    # CTF systems use reference sensors to measure ambient magnetic fields (building vibrations,
+    # distant equipment, etc.). This subtraction algorithm cancels this noise from the MEG
+    # channels. Grade 3 is the highest order and standard for analysis/visualization.
+    # Without this, raw channel values reflect environmental noise (millions of fT)
+    # instead of actual brain signals (tens to hundreds of fT).
+    raw.apply_gradient_compensation(3)  # type: ignore
+
+    return raw
+
 
 def main():
     parser = argparse.ArgumentParser(
@@ -43,8 +56,8 @@ def main():
     args = parser.parse_args()
 
     for path in args.files:
-        raw_ctf = load_channels(path)
-        channel_names = cast(list[str], raw_ctf.ch_names)  # type: ignore
+        raw = load_ctf_raw(path)
+        channel_names = cast(list[str], raw.ch_names)  # type: ignore
 
         if args.channel_index < 0:
             print("Channel index must be a positive integer", file=sys.stderr)
@@ -61,10 +74,10 @@ def main():
         print(f'Creating chunks for {path}')
         write_chunk_directory(
             path=path,
+            raw=raw,
             from_channel_index=args.channel_index,
             from_channel_name=channel_names[args.channel_index],  # type: ignore
             channel_count=args.channel_count,
-            loader=load_channels,
             chunk_size=args.chunk_size,
             destination=args.destination,
             prefix=args.prefix

python/loris_ephys_chunker/src/loris_ephys_chunker/scripts/edf_to_chunks.py

@@ -2,7 +2,6 @@
 
 import argparse
 import sys
-from collections.abc import Callable
 from pathlib import Path
 from typing import cast
 
@@ -13,8 +12,12 @@
 from loris_ephys_chunker.chunking import write_chunk_directory
 
 
-def load_channels(exclude: list[str]) -> Callable[[Path], RawEDF]:
-    return lambda path : mne.io.read_raw_edf(path, exclude=exclude, preload=False)  # type: ignore
+def read_edf_raw(path: Path, exclude: list[str]) -> RawEDF:
+    """
+    Read the EDF acquisition file into an MNE raw object.
+    """
+
+    return mne.io.read_raw_edf(path, exclude=exclude)  # type: ignore
 
 
 def main():
@@ -81,9 +84,10 @@ def main():
             # and avoid memory issues
             # we only load the channel at index channel_index+i
             exclude = channel_names[:channel_index] + channel_names[channel_index + 1:]
+            raw = read_edf_raw(path, exclude)
             write_chunk_directory(
                 path=path,
-                loader=load_channels(exclude),
+                raw=raw,
                 from_channel_index=channel_index,
                 from_channel_name=channel_names[channel_index],
                 channel_count=1,

python/loris_ephys_chunker/src/loris_ephys_chunker/scripts/eeglab_to_chunks.py

@@ -12,8 +12,12 @@
 from loris_ephys_chunker.chunking import write_chunk_directory
 
 
-def load_channels(path: Path) -> RawEEGLAB:
-    return mne.io.read_raw_eeglab(path, preload=False)  # type: ignore
+def read_eeglab_raw(path: Path) -> RawEEGLAB:
+    """
+    Read the EEGLAB acquisition file into an MNE raw object.
+    """
+
+    return mne.io.read_raw_eeglab(path)  # type: ignore
 
 
 def main():
@@ -50,12 +54,13 @@ def main():
             sys.exit("Channel count must be a positive integer")
 
         print(f'Creating chunks for {path}')
+        raw = read_eeglab_raw(path)
         write_chunk_directory(
             path=path,
+            raw=raw,
             from_channel_index=args.channel_index,
             from_channel_name=channel_names[args.channel_index],  # type: ignore
             channel_count=args.channel_count,
-            loader=load_channels,
             chunk_size=args.chunk_size,
             destination=args.destination,
             prefix=args.prefix