fNIRS improvements and minor fixes

hypyp/fnirs/dyad.py

@@ -51,6 +51,7 @@ def __init__(self, s1:Recording, s2:Recording, label:str='', is_pseudo:bool=Fals
         self.wtcs = None
         self.df = None
         self.is_pseudo = is_pseudo
+        self.is_intra = s1 == s2
 
         self.label = label
         if self.label == '':
@@ -269,6 +270,7 @@ def compute_wtcs(
         only_time_range: Tuple[float,float] | None = None,
         bin_seconds: float | None = None,
         period_cuts: List[float] | None = None,
+        frequency_cuts: List[float] | None = None,
         verbose: bool = False,
         with_intra: bool = True,
         downsample: int | None = None,
@@ -284,7 +286,8 @@ def compute_wtcs(
                                                 Defaults to None, which means all channels.
             only_time_range (Tuple[float,float] | None, optional): compute only a portion of the signal, defined as time range tuple (start, stop). Defaults to None.
             bin_seconds (float | None, optional): split the resulting WTC in time bins for balancing weights. Defaults to None.
-            period_cuts (List[float] | None, optional): split the resulting WTC in period/frequency bins for balancing weights and finer analysis. Defaults to None.
+            period_cuts (List[float] | None, optional): split the resulting WTC in period/frequency bins for balancing weights and finer analysis. Use this OR frequency_cuts, not both. Defaults to None.
+            frequency_cuts (List[float] | None, optional): split the resulting WTC in period/frequency bins for balancing weights and finer analysis. Use this OR period_cuts, not both. Defaults to None.
             verbose (bool, optional): verbose flag. Defaults to False.
             with_intra (bool, optional): compute intra-subject as well. Defaults to False.
             downsample (int | None, optional): downsample in time the resulting WTC. Useful to save memory space and faster display. Defaults to None.
@@ -296,23 +299,35 @@ def compute_wtcs(
         if wavelet is None:
             wavelet = ComplexMorletWavelet()
 
+        if period_cuts is not None and frequency_cuts is not None:
+            raise RuntimeError('Cannot specify both period_cuts and frequency_cuts')
+
         self.wtcs = []
 
         pairs = self.get_pairs(self.s1, self.s2, ch_match=ch_match, is_pseudo=self.is_pseudo)
 
-        for pair in pairs:
-            if verbose:
-                print(f'Running Wavelet Coherence for dyad "{self.label}" on pair "{pair.label}"')
-            if only_time_range is not None:
-                pair = pair.sub(only_time_range)
-            wtc = wavelet.wtc(pair, bin_seconds=bin_seconds, period_cuts=period_cuts)
-            if downsample is not None:
-                wtc.downsample_in_time(downsample)
+        if self.is_intra:
+            # Force with_intra when subject 1 is the same as subject 2
+            with_intra = True
+        else:
+            for pair in pairs:
+                if verbose:
+                    print(f'Running Wavelet Coherence for dyad "{self.label}" on pair "{pair.label}"')
+                if only_time_range is not None:
+                    pair = pair.sub(only_time_range)
+                wtc = wavelet.wtc(pair, bin_seconds=bin_seconds, period_cuts=period_cuts, frequency_cuts=frequency_cuts)
+                if downsample is not None:
+                    wtc.downsample_in_time(downsample)
 
-            self.wtcs.append(wtc)
+                self.wtcs.append(wtc)
 
         if with_intra:
-            for i, recording in enumerate([self.s1, self.s2]):
+            if self.is_intra:
+                subjects = [self.s1]
+            else:
+                subjects = [self.s1, self.s2]
+
+            for i, recording in enumerate(subjects):
                 recording.intra_wtcs = []
                 # if we have different channels for each recording, the intra wtc should use only the ones for this subject
                 if isinstance(ch_match, tuple):
@@ -327,7 +342,7 @@ def compute_wtcs(
                         print(f'Running Wavelet Coherence intra-subject "{recording.subject_label}" on pair "{pair.label}"')
                     if only_time_range is not None:
                         pair = pair.sub(only_time_range)
-                    wtc = wavelet.wtc(pair, bin_seconds=bin_seconds, period_cuts=period_cuts)
+                    wtc = wavelet.wtc(pair, bin_seconds=bin_seconds, period_cuts=period_cuts, frequency_cuts=frequency_cuts)
                     if downsample is not None:
                         wtc.downsample_in_time(downsample)
                     recording.intra_wtcs.append(wtc)

hypyp/fnirs/preprocessor/implementations/mne_preprocessor_as_is.py

@@ -29,7 +29,7 @@ def read_file(self, path:str, verbose:bool=False):
         if DataBrowser.is_path_snirf(path):
             return mne.io.read_raw_snirf(path, preload=True, verbose=verbose)
 
-        return None
+        raise ValueError(f"No reader for file {path}")
 
     def run(self, raw: mne.io.Raw, verbose: bool = False) -> list[MneStep]:
         if verbose:

hypyp/fnirs/preprocessor/implementations/mne_preprocessor_raw_to_haemo.py

@@ -26,46 +26,38 @@ class MnePreprocessorRawToHaemo(MnePreprocessorAsIs):
     already preprocessed and will be returned as-is.
 
     """
-    def __init__(self):
+    def __init__(self, ppf=6.0, hp_filter=0.02, lp_filter=0.7, quality_sci_threshold=0.1):
         super().__init__()
+        self.ppf = ppf
+        self.hp_filter = hp_filter
+        self.lp_filter = lp_filter
+        self.quality_sci_threshold = quality_sci_threshold
 
     def run(self, raw:mne.io.Raw, verbose:bool=False) -> list[MneStep]:
+        # If have hbo or hbr, it means it is already in concentration
+        if len(mne.pick_types(raw.info, fnirs=['hbo', 'hbr'])) > 0:
+            raise ValueError('Loaded data seems to already be in concentrations. Use MnePreprocessorAsIs instead.')
+
         if verbose:
             print('Using MnePreprocessorRawToHaemo, converting fNIRS raw data to haemoglobin concentrations')
 
         steps = []
         steps.append(MneStep(raw, PREPROCESS_STEP_BASE_KEY, PREPROCESS_STEP_BASE_DESC))
 
-        haemo_picks = mne.pick_types(raw.info, fnirs=['hbo', 'hbr'])
-
-        # If we have haemo_picks, it means it is already preprocessed
-        # TODO: this code flow if confusing
-        if len(haemo_picks) > 0:
-            if verbose:
-                print('Data loaded using MNE seems to already be converted to haemoglobin concentrations')
-
-            steps.append(MneStep(
-                raw.copy().pick(haemo_picks),
-                PREPROCESS_STEP_HAEMO_FILTERED_KEY,
-                PREPROCESS_STEP_HAEMO_FILTERED_DESC
-            ))
-            return steps
-
         raw_od = mne.preprocessing.nirs.optical_density(raw)
         steps.append(MneStep(raw_od, PREPROCESS_STEP_OD_KEY, PREPROCESS_STEP_OD_DESC))
 
         quality_sci = mne.preprocessing.nirs.scalp_coupling_index(raw_od)
-        raw_od.info['bads'] = list(compress(raw_od.ch_names, quality_sci < 0.1))
+        raw_od.info['bads'] = list(compress(raw_od.ch_names, quality_sci < self.quality_sci_threshold))
         picks = mne.pick_types(raw_od.info, fnirs=True, exclude='bads')
         raw_od_clean = raw_od.copy().pick(picks)
         steps.append(MneStep(raw_od_clean, PREPROCESS_STEP_OD_CLEAN_KEY, PREPROCESS_STEP_OD_CLEAN_DESC))
 
-        # TODO: should set Partial Pathlength Factor (PPF) depending on age.
-        # See https://doi.org/10.1117/1.JBO.18.10.105004, Table 1
-        raw_haemo: mne.io.Raw = mne.preprocessing.nirs.beer_lambert_law(raw_od_clean)
+        # For partial pathlength, see https://doi.org/10.1117/1.JBO.18.10.105004, Table 1
+        raw_haemo: mne.io.Raw = mne.preprocessing.nirs.beer_lambert_law(raw_od_clean.copy(), ppf=self.ppf)
         steps.append(MneStep(raw_haemo, PREPROCESS_STEP_HAEMO_KEY, PREPROCESS_STEP_HAEMO_DESC))
-        raw_haemo_filtered = raw_haemo.copy().filter(0.02, 0.7, verbose=verbose)
 
+        raw_haemo_filtered = raw_haemo.copy().filter(self.hp_filter, self.lp_filter, verbose=verbose, method='iir')
         steps.append(MneStep(raw_haemo_filtered, PREPROCESS_STEP_HAEMO_FILTERED_KEY, PREPROCESS_STEP_HAEMO_FILTERED_DESC))
 
         return steps