From c57722352ca36027b6d58f17c4fead9df76f2d03 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Mon, 21 Jul 2025 16:00:26 +1000
Subject: [PATCH 01/57] Added draft percentile filter in threshold_by_domains.

---
 seas/domains.py | 27 +++++++++++++++++++++------
 1 file changed, 21 insertions(+), 6 deletions(-)

diff --git a/seas/domains.py b/seas/domains.py
index b4351d6..06f7919 100644
--- a/seas/domains.py
+++ b/seas/domains.py
@@ -626,7 +626,8 @@ def threshold_by_domains(components: dict,
                    map_only: bool = True,
                    apply_filter_mean: bool = True,
                    max_loops: int = 2,
-                   ignore_small: bool = True):
+                   ignore_small: bool = True,
+                   thresh_type: str = 'max'):
     '''
     Creates a domain map from extracted independent components.  A pixelwise maximum projection of the blurred signal components is taken through the n_components axis, to create a flattened representation of where a domain was maximally significant across the cortical surface.  Components with multiple noncontiguous significant regions are counted as two distinct domains.
 
@@ -711,12 +712,26 @@ def threshold_by_domains(components: dict,
                 blurred = cv2.GaussianBlur(eigenbrain, (blur, blur), 0)
                 eig_vec.T[index] = blurred.flat
 
-    # This is the money section, return indices across each eig_vec (loading vector for component) where loading is max
-    domain_ROIs_vector = np.argmax(np.abs(eig_vec), axis=1)
-    # Then threshold by clearing eig_vec outside of max indices
     mask = np.zeros_like(eig_vec, dtype=bool)
-    mask[np.arange(eig_vec.shape[0]), domain_ROIs_vector] = True
-    eig_vec[~mask] = 0
+    flipped = components['flipped']
+    match thresh_type:
+        case 'max':
+            # Return indices across each eig_vec (loading vector for component) where loading is max
+            domain_ROIs_vector = np.argmax(np.abs(eig_vec), axis=1)
+            # Then threshold by clearing eig_vec outside of max indices
+            mask[np.arange(eig_vec.shape[0]), domain_ROIs_vector] = True
+            eig_vec[~mask] = 0
+        case 'percentile':
+            # Flip ICs were necessary using flipped from dict
+            flipped_eig_vec = np.multiply(flipped, eig_vec)
+            # Calculate 95 percentile cutoff for each IC
+            cutoff_vector = np.percentile(flipped_eig_vec, 0.95, axis=1)
+            # Mask for all values above cutoff
+            for i in np.arange(eig_vec.shape[0]):
+                mask[i, :] = flipped_eig_vec[i] > cutoff_vector[i]
+            eig_vec[~mask] = 0
+        case _:
+            print("Threshold type is neither max nor percentile.")
 
     output['eig_vec'] = eig_vec
 

From 33868fddd8bb71975798efddae703743482c00a9 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Tue, 22 Jul 2025 10:23:54 +1000
Subject: [PATCH 02/57] Debug flipped assignment.

---
 seas/domains.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/seas/domains.py b/seas/domains.py
index 06f7919..0383af1 100644
--- a/seas/domains.py
+++ b/seas/domains.py
@@ -713,7 +713,7 @@ def threshold_by_domains(components: dict,
                 eig_vec.T[index] = blurred.flat
 
     mask = np.zeros_like(eig_vec, dtype=bool)
-    flipped = components['flipped']
+    
     match thresh_type:
         case 'max':
             # Return indices across each eig_vec (loading vector for component) where loading is max
@@ -722,6 +722,7 @@ def threshold_by_domains(components: dict,
             mask[np.arange(eig_vec.shape[0]), domain_ROIs_vector] = True
             eig_vec[~mask] = 0
         case 'percentile':
+            flipped = components['flipped']
             # Flip ICs were necessary using flipped from dict
             flipped_eig_vec = np.multiply(flipped, eig_vec)
             # Calculate 95 percentile cutoff for each IC

From 08959ac8e4d1a24dd8e40c0a2e0a6151b9da6572 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Tue, 22 Jul 2025 13:54:12 +1000
Subject: [PATCH 03/57] Added mask case to thresh_by_type to yield binary
 activity masks.

---
 seas/domains.py | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/seas/domains.py b/seas/domains.py
index 06f7919..e6b7af2 100644
--- a/seas/domains.py
+++ b/seas/domains.py
@@ -721,6 +721,12 @@ def threshold_by_domains(components: dict,
             # Then threshold by clearing eig_vec outside of max indices
             mask[np.arange(eig_vec.shape[0]), domain_ROIs_vector] = True
             eig_vec[~mask] = 0
+        case 'mask':
+            # Return indices across each eig_vec (loading vector for component) where loading is max
+            domain_ROIs_vector = np.argmax(np.abs(eig_vec), axis=1)
+            # Then threshold by clearing eig_vec outside of max indices
+            mask[np.arange(eig_vec.shape[0]), domain_ROIs_vector] = True
+            eig_vec = mask
         case 'percentile':
             # Flip ICs were necessary using flipped from dict
             flipped_eig_vec = np.multiply(flipped, eig_vec)
@@ -731,7 +737,7 @@ def threshold_by_domains(components: dict,
                 mask[i, :] = flipped_eig_vec[i] > cutoff_vector[i]
             eig_vec[~mask] = 0
         case _:
-            print("Threshold type is neither max nor percentile.")
+            print("Threshold type is neither max, mask, nor percentile.")
 
     output['eig_vec'] = eig_vec
 

From 8d57966804b59714c6e04ab46ae4a1aaca41b024 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Tue, 22 Jul 2025 14:47:34 +1000
Subject: [PATCH 04/57] Added alternative mean addition routine, incl. option
 and masks import from dict.

---
 seas/ica.py | 34 ++++++++++++++++++++++++++++------
 1 file changed, 28 insertions(+), 6 deletions(-)

diff --git a/seas/ica.py b/seas/ica.py
index a00721e..80c08f2 100644
--- a/seas/ica.py
+++ b/seas/ica.py
@@ -304,6 +304,7 @@ def rebuild(components: dict,
             t_start: int = None,
             t_stop: int = None,
             apply_mean_filter: bool = True,
+            apply_masked_mean: bool = False,
             filter_method: str = 'wavelet',
             fps: float = 7.5,
             include_noise: bool = True):
@@ -344,6 +345,7 @@ def rebuild(components: dict,
     roimask = components['roimask']
     shape = components['shape']
     mean = components['mean']
+    masks = components['masks']
     n_components = components['n_components']
     dtype = np.float32
 
@@ -406,14 +408,34 @@ def rebuild(components: dict,
     data_r = np.dot(eig_vec[:, reconstruct_indices],
                     eig_mix[t_start:t_stop, reconstruct_indices].T).T
 
-    if apply_mean_filter:
-        mean_filtered = filter_mean(mean, filter_method, fps=fps)
-        data_r += mean_filtered[t_start:t_stop, None]
+    if apply_masked_mean:
+        assert masks is not None, \
+        "Masks have not been assigned to dictionary"
+        # Apply mean to masks only, zeroing unmasked pixels
+        if apply_mean_filter:
+            combined_mask = np.any(masks[:, reconstruct_indices], axis=1)
+            mean_to_add = np.zeros_like(data_r)
+            mean_filtered = filter_mean(mean, filter_method, fps=fps)
+            mean_to_add[:, combined_mask] = mean_filtered[t_start:t_stop, None]
+            data_r += mean_to_add
 
+        else:
+            print('Not filtering mean')
+            combined_mask = np.any(masks[:, reconstruct_indices], axis=1)
+            mean_to_add = np.zeros_like(data_r)
+            mean_filtered = None
+            mean_to_add[:, combined_mask] = mean[t_start:t_stop]
+            data_r += mean_to_add
     else:
-        print('Not filtering mean')
-        mean_filtered = None
-        data_r += mean[t_start:t_stop, None]
+        # Run original readdition of mean
+        if apply_mean_filter:
+            mean_filtered = filter_mean(mean, filter_method, fps=fps)
+            data_r += mean_filtered[t_start:t_stop, None]
+
+        else:
+            print('Not filtering mean')
+            mean_filtered = None
+            data_r += mean[t_start:t_stop, None]
 
     print('Done!')
 

From 8c1ac66be86cff369872c39f21c8e581c210dcd8 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Tue, 22 Jul 2025 15:10:40 +1000
Subject: [PATCH 05/57] Changed masks output to dictionary rather than encoding
 in eig_vec.

---
 seas/domains.py | 11 +++--------
 1 file changed, 3 insertions(+), 8 deletions(-)

diff --git a/seas/domains.py b/seas/domains.py
index 4f06a18..698f04f 100644
--- a/seas/domains.py
+++ b/seas/domains.py
@@ -721,12 +721,6 @@ def threshold_by_domains(components: dict,
             # Then threshold by clearing eig_vec outside of max indices
             mask[np.arange(eig_vec.shape[0]), domain_ROIs_vector] = True
             eig_vec[~mask] = 0
-        case 'mask':
-            # Return indices across each eig_vec (loading vector for component) where loading is max
-            domain_ROIs_vector = np.argmax(np.abs(eig_vec), axis=1)
-            # Then threshold by clearing eig_vec outside of max indices
-            mask[np.arange(eig_vec.shape[0]), domain_ROIs_vector] = True
-            eig_vec = mask
         case 'percentile':
             flipped = components['flipped']
             # Flip ICs were necessary using flipped from dict
@@ -738,8 +732,9 @@ def threshold_by_domains(components: dict,
                 mask[i, :] = flipped_eig_vec[i] > cutoff_vector[i]
             eig_vec[~mask] = 0
         case _:
-            print("Threshold type is neither max, mask, nor percentile.")
-
+            print("Threshold type is neither max nor percentile.")
+    
+    output['masks'] = mask
     output['eig_vec'] = eig_vec
 
     return output

From 26ad5262a20f15f88bb9639676c234f4e10eafdf Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Tue, 22 Jul 2025 15:50:45 +1000
Subject: [PATCH 06/57] Broadcasting debug on apply_masked_mean.

---
 seas/ica.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/seas/ica.py b/seas/ica.py
index 80c08f2..921c0de 100644
--- a/seas/ica.py
+++ b/seas/ica.py
@@ -424,7 +424,7 @@ def rebuild(components: dict,
             combined_mask = np.any(masks[:, reconstruct_indices], axis=1)
             mean_to_add = np.zeros_like(data_r)
             mean_filtered = None
-            mean_to_add[:, combined_mask] = mean[t_start:t_stop]
+            mean_to_add[:, combined_mask] = mean[t_start:t_stop, None]
             data_r += mean_to_add
     else:
         # Run original readdition of mean

From cd7d91385ba084edbac14bd3a6a4aa2b9f080603 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Wed, 23 Jul 2025 13:18:41 +1000
Subject: [PATCH 07/57] Add temp variable for eig_vec to do blurring for
 threshold (rather than data).

---
 seas/domains.py | 32 ++++++++++++++++----------------
 1 file changed, 16 insertions(+), 16 deletions(-)

diff --git a/seas/domains.py b/seas/domains.py
index 698f04f..05dde92 100644
--- a/seas/domains.py
+++ b/seas/domains.py
@@ -692,6 +692,9 @@ def threshold_by_domains(components: dict,
             signal_indices = np.where(artifact_components == 0)[0]
         # eig_vec = eig_vec[:, signal_indices] # Don't change number of ICs, we're updating back to dict
 
+    threshold_vec = components['eig_vec'].copy()
+    mask = np.zeros_like(eig_vec, dtype=bool)
+    
     if blur:
         print('blurring domains...')
         assert type(blur) is int, 'blur was not valid'
@@ -701,39 +704,36 @@ def threshold_by_domains(components: dict,
         eigenbrain = np.empty(shape)
         eigenbrain[:] = np.NAN
 
-        for index in range(eig_vec.shape[1]):
+        for index in range(threshold_vec.shape[1]):
 
             if roimask is not None:
-                eigenbrain.flat[maskind] = eig_vec.T[index]
+                eigenbrain.flat[maskind] = threshold_vec.T[index]
                 blurred = cv2.GaussianBlur(eigenbrain, (blur, blur), 0)
-                eig_vec.T[index] = blurred.flat[maskind]
+                threshold_vec.T[index] = blurred.flat[maskind]
             else:
-                eigenbrain.flat = eig_vec.T[index]
+                eigenbrain.flat = threshold_vec.T[index]
                 blurred = cv2.GaussianBlur(eigenbrain, (blur, blur), 0)
-                eig_vec.T[index] = blurred.flat
+                threshold_vec.T[index] = blurred.flat
 
-    mask = np.zeros_like(eig_vec, dtype=bool)
-    
     match thresh_type:
         case 'max':
             # Return indices across each eig_vec (loading vector for component) where loading is max
-            domain_ROIs_vector = np.argmax(np.abs(eig_vec), axis=1)
+            domain_ROIs_vector = np.argmax(np.abs(threshold_vec), axis=1)
             # Then threshold by clearing eig_vec outside of max indices
-            mask[np.arange(eig_vec.shape[0]), domain_ROIs_vector] = True
-            eig_vec[~mask] = 0
+            mask[np.arange(threshold_vec.shape[0]), domain_ROIs_vector] = True
         case 'percentile':
             flipped = components['flipped']
-            # Flip ICs were necessary using flipped from dict
-            flipped_eig_vec = np.multiply(flipped, eig_vec)
+            # Flip ICs where necessary using flipped from dict
+            flipped_threshold_vec = np.multiply(flipped, threshold_vec)
             # Calculate 95 percentile cutoff for each IC
-            cutoff_vector = np.percentile(flipped_eig_vec, 0.95, axis=1)
+            cutoff_vector = np.percentile(flipped_threshold_vec, 0.95, axis=1)
             # Mask for all values above cutoff
             for i in np.arange(eig_vec.shape[0]):
-                mask[i, :] = flipped_eig_vec[i] > cutoff_vector[i]
-            eig_vec[~mask] = 0
+                mask[i, :] = flipped_threshold_vec[i] > cutoff_vector[i]
         case _:
             print("Threshold type is neither max nor percentile.")
-    
+
+    eig_vec[~mask] = 0
     output['masks'] = mask
     output['eig_vec'] = eig_vec
 

From c29e041350684660f6cfca8e86f07b67e8f2015b Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Wed, 23 Jul 2025 21:14:45 +1000
Subject: [PATCH 08/57] Added small mask ROI filtering before blurring tp
 threshold_by_domains.

---
 seas/domains.py | 56 ++++++++++++++++++++++++++-----------------------
 1 file changed, 30 insertions(+), 26 deletions(-)

diff --git a/seas/domains.py b/seas/domains.py
index 05dde92..3b24977 100644
--- a/seas/domains.py
+++ b/seas/domains.py
@@ -16,6 +16,7 @@
 from seas.rois import make_mask
 from seas.colormaps import save_colorbar, REGION_COLORMAP, DEFAULT_COLORMAP
 
+from skimage.morphology import remove_small_objects
 
 def get_domain_map(components: dict,
                    blur: int = 21,
@@ -622,6 +623,7 @@ def write_frame(frame):
 
 def threshold_by_domains(components: dict,
                    blur: int = 1,
+                   min_mask_size: int = 64,
                    min_size_ratio: float = 0.1,
                    map_only: bool = True,
                    apply_filter_mean: bool = True,
@@ -691,48 +693,50 @@ def threshold_by_domains(components: dict,
             print('no noise components found')
             signal_indices = np.where(artifact_components == 0)[0]
         # eig_vec = eig_vec[:, signal_indices] # Don't change number of ICs, we're updating back to dict
-
-    threshold_vec = components['eig_vec'].copy()
-    mask = np.zeros_like(eig_vec, dtype=bool)
     
-    if blur:
-        print('blurring domains...')
-        assert type(blur) is int, 'blur was not valid'
-        if blur % 2 != 1:
-            blur += 1
-
-        eigenbrain = np.empty(shape)
-        eigenbrain[:] = np.NAN
-
-        for index in range(threshold_vec.shape[1]):
-
-            if roimask is not None:
-                eigenbrain.flat[maskind] = threshold_vec.T[index]
-                blurred = cv2.GaussianBlur(eigenbrain, (blur, blur), 0)
-                threshold_vec.T[index] = blurred.flat[maskind]
-            else:
-                eigenbrain.flat = threshold_vec.T[index]
-                blurred = cv2.GaussianBlur(eigenbrain, (blur, blur), 0)
-                threshold_vec.T[index] = blurred.flat
+    mask = np.zeros_like(eig_vec, dtype=bool)
 
     match thresh_type:
         case 'max':
             # Return indices across each eig_vec (loading vector for component) where loading is max
-            domain_ROIs_vector = np.argmax(np.abs(threshold_vec), axis=1)
+            threshold_ROIs_vector = np.argmax(np.abs(eig_vec), axis=1)
             # Then threshold by clearing eig_vec outside of max indices
-            mask[np.arange(threshold_vec.shape[0]), domain_ROIs_vector] = True
+            mask[np.arange(eig_vec.shape[0]), threshold_ROIs_vector] = True
         case 'percentile':
             flipped = components['flipped']
             # Flip ICs where necessary using flipped from dict
-            flipped_threshold_vec = np.multiply(flipped, threshold_vec)
+            flipped_threshold_vec = np.multiply(flipped, eig_vec)
             # Calculate 95 percentile cutoff for each IC
-            cutoff_vector = np.percentile(flipped_threshold_vec, 0.95, axis=1)
+            cutoff_vector = np.percentile(flipped, 0.95, axis=1)
             # Mask for all values above cutoff
             for i in np.arange(eig_vec.shape[0]):
                 mask[i, :] = flipped_threshold_vec[i] > cutoff_vector[i]
         case _:
             print("Threshold type is neither max nor percentile.")
 
+    # Filter small mask ROIs and smooth using blur
+    if blur:
+        print('blurring domains...')
+        assert type(blur) is int, 'blur was not valid'
+        if blur % 2 != 1:
+            blur += 1
+
+        eigenbrain = np.empty(shape)
+        eigenbrain[:] = np.NAN
+
+        for index in range(mask.shape[1]):
+
+            if roimask is not None:
+                eigenbrain.flat[maskind] = mask.T[index]
+                filtered = remove_small_objects(eigenbrain, min_size=min_mask_size, connectivity=1)
+                blurred = cv2.GaussianBlur(filtered, (blur, blur), 0)
+                mask.T[index] = blurred.flat[maskind]
+            else:
+                eigenbrain.flat = mask.T[index]
+                filtered = remove_small_objects(eigenbrain, min_size=min_mask_size, connectivity=1)
+                blurred = cv2.GaussianBlur(filtered, (blur, blur), 0)
+                mask.T[index] = blurred.flat
+
     eig_vec[~mask] = 0
     output['masks'] = mask
     output['eig_vec'] = eig_vec

From 1ce6a144db1930d7408db8008e8908d2fca6666c Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Wed, 23 Jul 2025 21:25:02 +1000
Subject: [PATCH 09/57] Debug eigenbrain mask dtype.

---
 seas/domains.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/seas/domains.py b/seas/domains.py
index 3b24977..eece566 100644
--- a/seas/domains.py
+++ b/seas/domains.py
@@ -721,8 +721,8 @@ def threshold_by_domains(components: dict,
         if blur % 2 != 1:
             blur += 1
 
-        eigenbrain = np.empty(shape)
-        eigenbrain[:] = np.NAN
+        eigenbrain = np.zeros(shape, dtype=bool)
+        # eigenbrain[:] = np.NAN
 
         for index in range(mask.shape[1]):
 

From c23ba85d775883e8119cf938ccba50b24f672b55 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Wed, 23 Jul 2025 21:44:14 +1000
Subject: [PATCH 10/57] Debug dtype changes between morph filtering and
 blurring.

---
 seas/domains.py | 12 +++++++-----
 1 file changed, 7 insertions(+), 5 deletions(-)

diff --git a/seas/domains.py b/seas/domains.py
index eece566..d362414 100644
--- a/seas/domains.py
+++ b/seas/domains.py
@@ -721,15 +721,17 @@ def threshold_by_domains(components: dict,
         if blur % 2 != 1:
             blur += 1
 
-        eigenbrain = np.zeros(shape, dtype=bool)
-        # eigenbrain[:] = np.NAN
+        eigenmask = np.zeros(shape, dtype=bool)
+        eigenbrain = np.empty(shape)
+        eigenbrain[:] = np.NAN
 
         for index in range(mask.shape[1]):
 
             if roimask is not None:
-                eigenbrain.flat[maskind] = mask.T[index]
-                filtered = remove_small_objects(eigenbrain, min_size=min_mask_size, connectivity=1)
-                blurred = cv2.GaussianBlur(filtered, (blur, blur), 0)
+                eigenmask.flat[maskind] = mask.T[index]
+                filtered = remove_small_objects(eigenmask, min_size=min_mask_size, connectivity=1)
+                eigenbrain.flat[maskind] = filtered.flat
+                blurred = cv2.GaussianBlur(eigenbrain.flat, (blur, blur), 0)
                 mask.T[index] = blurred.flat[maskind]
             else:
                 eigenbrain.flat = mask.T[index]

From a6c7e8f1cd937ab2552f361ecba566b777dca9d5 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Wed, 23 Jul 2025 22:04:38 +1000
Subject: [PATCH 11/57] Added draft of IC timecourse lpf as part of
 threshold_by_domains.

---
 seas/domains.py | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/seas/domains.py b/seas/domains.py
index d362414..a38bd39 100644
--- a/seas/domains.py
+++ b/seas/domains.py
@@ -15,6 +15,7 @@
 from seas.ica import rebuild_mean_roi_timecourse, filter_mean
 from seas.rois import make_mask
 from seas.colormaps import save_colorbar, REGION_COLORMAP, DEFAULT_COLORMAP
+from seas.signalanalysis import butterworth
 
 from skimage.morphology import remove_small_objects
 
@@ -670,6 +671,7 @@ def threshold_by_domains(components: dict,
     output['domain_blur'] = blur
 
     eig_vec = components['eig_vec'].copy()
+    eig_mix = components['eig_mix'].copy()
     shape = components['shape']
     shape = (shape[1], shape[2])
 
@@ -740,8 +742,16 @@ def threshold_by_domains(components: dict,
                 mask.T[index] = blurred.flat
 
     eig_vec[~mask] = 0
+
+    # Filter component timecourses
+    timecourses = eig_mix.T
+    lpf_timecourses = np.zeros_like(timecourses)
+    for index in range(timecourses.shape[1]):
+        lpf_timecourses[index] = butterworth(timecourses[index], high=1.0)
+    
     output['masks'] = mask
     output['eig_vec'] = eig_vec
+    output['eig_mix'] = lpf_timecourses.T
 
     return output
 

From 3c4ce61596c8acadf5c43ae66fd94b62d909392e Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Wed, 23 Jul 2025 22:07:20 +1000
Subject: [PATCH 12/57] Debug blurring.

---
 seas/domains.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/seas/domains.py b/seas/domains.py
index a38bd39..ec20ad7 100644
--- a/seas/domains.py
+++ b/seas/domains.py
@@ -733,7 +733,7 @@ def threshold_by_domains(components: dict,
                 eigenmask.flat[maskind] = mask.T[index]
                 filtered = remove_small_objects(eigenmask, min_size=min_mask_size, connectivity=1)
                 eigenbrain.flat[maskind] = filtered.flat
-                blurred = cv2.GaussianBlur(eigenbrain.flat, (blur, blur), 0)
+                blurred = cv2.GaussianBlur(eigenbrain, (blur, blur), 0)
                 mask.T[index] = blurred.flat[maskind]
             else:
                 eigenbrain.flat = mask.T[index]

From a22b79622257c13947c477bb7cfe694446ee23dc Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Wed, 23 Jul 2025 22:10:30 +1000
Subject: [PATCH 13/57] Debug timecourse indexing.

---
 seas/domains.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/seas/domains.py b/seas/domains.py
index ec20ad7..85b3f82 100644
--- a/seas/domains.py
+++ b/seas/domains.py
@@ -746,7 +746,7 @@ def threshold_by_domains(components: dict,
     # Filter component timecourses
     timecourses = eig_mix.T
     lpf_timecourses = np.zeros_like(timecourses)
-    for index in range(timecourses.shape[1]):
+    for index in range(timecourses.shape[0]):
         lpf_timecourses[index] = butterworth(timecourses[index], high=1.0)
     
     output['masks'] = mask

From 0c37375292ccc8d0d5ee846f2b8150525ca0ecf6 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Wed, 23 Jul 2025 23:11:30 +1000
Subject: [PATCH 14/57] Commented out lpf code from threshold_by_domains, need
 to rewrite.

---
 seas/domains.py | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/seas/domains.py b/seas/domains.py
index 85b3f82..d0e7b24 100644
--- a/seas/domains.py
+++ b/seas/domains.py
@@ -743,15 +743,18 @@ def threshold_by_domains(components: dict,
 
     eig_vec[~mask] = 0
 
+    # The following filtering isn't working, needs to be in its own function after 
+    # readdition of the mean? Maybe as part of seas.ica.rebuild?
+
     # Filter component timecourses
-    timecourses = eig_mix.T
-    lpf_timecourses = np.zeros_like(timecourses)
-    for index in range(timecourses.shape[0]):
-        lpf_timecourses[index] = butterworth(timecourses[index], high=1.0)
+    # timecourses = eig_mix.T
+    # lpf_timecourses = np.zeros_like(timecourses)
+    # for index in range(timecourses.shape[0]):
+    #     lpf_timecourses[index] = butterworth(timecourses[index], high=1.0)
     
     output['masks'] = mask
     output['eig_vec'] = eig_vec
-    output['eig_mix'] = lpf_timecourses.T
+    # output['eig_mix'] = lpf_timecourses.T
 
     return output
 

From 0d2b19a1d86a32ef601a6338208cc922715b8696 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Thu, 24 Jul 2025 11:04:25 +1000
Subject: [PATCH 15/57] Debug indexing issue.

---
 seas/domains.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/seas/domains.py b/seas/domains.py
index d0e7b24..8bf30c3 100644
--- a/seas/domains.py
+++ b/seas/domains.py
@@ -671,7 +671,7 @@ def threshold_by_domains(components: dict,
     output['domain_blur'] = blur
 
     eig_vec = components['eig_vec'].copy()
-    eig_mix = components['eig_mix'].copy()
+    # eig_mix = components['eig_mix'].copy()
     shape = components['shape']
     shape = (shape[1], shape[2])
 
@@ -732,9 +732,9 @@ def threshold_by_domains(components: dict,
             if roimask is not None:
                 eigenmask.flat[maskind] = mask.T[index]
                 filtered = remove_small_objects(eigenmask, min_size=min_mask_size, connectivity=1)
-                eigenbrain.flat[maskind] = filtered.flat
+                eigenbrain.flat = filtered.flat
                 blurred = cv2.GaussianBlur(eigenbrain, (blur, blur), 0)
-                mask.T[index] = blurred.flat[maskind]
+                mask.T[index] = blurred.flat
             else:
                 eigenbrain.flat = mask.T[index]
                 filtered = remove_small_objects(eigenbrain, min_size=min_mask_size, connectivity=1)

From 2658d3fda4c98efb8c758bd9d925f8bbea42bc25 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Thu, 24 Jul 2025 11:15:16 +1000
Subject: [PATCH 16/57] Revert indexing issue and debug blurring dtype.

---
 seas/domains.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/seas/domains.py b/seas/domains.py
index 8bf30c3..4029edd 100644
--- a/seas/domains.py
+++ b/seas/domains.py
@@ -732,9 +732,9 @@ def threshold_by_domains(components: dict,
             if roimask is not None:
                 eigenmask.flat[maskind] = mask.T[index]
                 filtered = remove_small_objects(eigenmask, min_size=min_mask_size, connectivity=1)
-                eigenbrain.flat = filtered.flat
+                eigenbrain.flat[maskind] = filtered.flat.astype(np.float32)
                 blurred = cv2.GaussianBlur(eigenbrain, (blur, blur), 0)
-                mask.T[index] = blurred.flat
+                mask.T[index] = blurred.flat[maskind]
             else:
                 eigenbrain.flat = mask.T[index]
                 filtered = remove_small_objects(eigenbrain, min_size=min_mask_size, connectivity=1)

From 8c8c9aae5a3f2b41b172eac4b3b60d05004fc070 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Thu, 24 Jul 2025 11:18:44 +1000
Subject: [PATCH 17/57] Debug dtype assignment.

---
 seas/domains.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/seas/domains.py b/seas/domains.py
index 4029edd..22134e7 100644
--- a/seas/domains.py
+++ b/seas/domains.py
@@ -732,7 +732,8 @@ def threshold_by_domains(components: dict,
             if roimask is not None:
                 eigenmask.flat[maskind] = mask.T[index]
                 filtered = remove_small_objects(eigenmask, min_size=min_mask_size, connectivity=1)
-                eigenbrain.flat[maskind] = filtered.flat.astype(np.float32)
+                filtered_float = filtered.astype(np.float32)
+                eigenbrain.flat[maskind] = filtered_float.flat
                 blurred = cv2.GaussianBlur(eigenbrain, (blur, blur), 0)
                 mask.T[index] = blurred.flat[maskind]
             else:

From 1841894d188359d71ced1ad2b505944a13e9b799 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Thu, 24 Jul 2025 11:37:07 +1000
Subject: [PATCH 18/57] Try filtering routine without blurring (which seems to
 be blanking whole video).

---
 seas/domains.py | 20 ++++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/seas/domains.py b/seas/domains.py
index 22134e7..df204d4 100644
--- a/seas/domains.py
+++ b/seas/domains.py
@@ -671,7 +671,7 @@ def threshold_by_domains(components: dict,
     output['domain_blur'] = blur
 
     eig_vec = components['eig_vec'].copy()
-    # eig_mix = components['eig_mix'].copy()
+    eig_mix = components['eig_mix'].copy()
     shape = components['shape']
     shape = (shape[1], shape[2])
 
@@ -732,10 +732,10 @@ def threshold_by_domains(components: dict,
             if roimask is not None:
                 eigenmask.flat[maskind] = mask.T[index]
                 filtered = remove_small_objects(eigenmask, min_size=min_mask_size, connectivity=1)
-                filtered_float = filtered.astype(np.float32)
-                eigenbrain.flat[maskind] = filtered_float.flat
-                blurred = cv2.GaussianBlur(eigenbrain, (blur, blur), 0)
-                mask.T[index] = blurred.flat[maskind]
+                # filtered_float = filtered.astype(np.float32)
+                # eigenbrain.flat[maskind] = filtered_float.flat
+                # blurred = cv2.GaussianBlur(eigenbrain, (blur, blur), 0)
+                mask.T[index] = filtered.flat[maskind]
             else:
                 eigenbrain.flat = mask.T[index]
                 filtered = remove_small_objects(eigenbrain, min_size=min_mask_size, connectivity=1)
@@ -748,14 +748,14 @@ def threshold_by_domains(components: dict,
     # readdition of the mean? Maybe as part of seas.ica.rebuild?
 
     # Filter component timecourses
-    # timecourses = eig_mix.T
-    # lpf_timecourses = np.zeros_like(timecourses)
-    # for index in range(timecourses.shape[0]):
-    #     lpf_timecourses[index] = butterworth(timecourses[index], high=1.0)
+    timecourses = eig_mix.T
+    lpf_timecourses = np.zeros_like(timecourses)
+    for index in range(timecourses.shape[0]):
+        lpf_timecourses[index] = butterworth(timecourses[index], high=1.0)
     
     output['masks'] = mask
     output['eig_vec'] = eig_vec
-    # output['eig_mix'] = lpf_timecourses.T
+    output['eig_mix'] = lpf_timecourses.T
 
     return output
 

From f8609229ee7911f5f9713a9c8c3216c3cdfe12de Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Thu, 24 Jul 2025 16:27:34 +1000
Subject: [PATCH 19/57] Add option apply_component_filter to ica.rebuild and
 domains.threshold_by_domains for lpf on IC timecourses.

---
 seas/domains.py | 20 ++++++++++----------
 seas/ica.py     | 10 ++++++++++
 2 files changed, 20 insertions(+), 10 deletions(-)

diff --git a/seas/domains.py b/seas/domains.py
index df204d4..6d18113 100644
--- a/seas/domains.py
+++ b/seas/domains.py
@@ -628,6 +628,7 @@ def threshold_by_domains(components: dict,
                    min_size_ratio: float = 0.1,
                    map_only: bool = True,
                    apply_filter_mean: bool = True,
+                   apply_component_filter: bool = False,
                    max_loops: int = 2,
                    ignore_small: bool = True,
                    thresh_type: str = 'max'):
@@ -671,7 +672,7 @@ def threshold_by_domains(components: dict,
     output['domain_blur'] = blur
 
     eig_vec = components['eig_vec'].copy()
-    eig_mix = components['eig_mix'].copy()
+
     shape = components['shape']
     shape = (shape[1], shape[2])
 
@@ -744,19 +745,18 @@ def threshold_by_domains(components: dict,
 
     eig_vec[~mask] = 0
 
-    # The following filtering isn't working, needs to be in its own function after 
-    # readdition of the mean? Maybe as part of seas.ica.rebuild?
-
     # Filter component timecourses
-    timecourses = eig_mix.T
-    lpf_timecourses = np.zeros_like(timecourses)
-    for index in range(timecourses.shape[0]):
-        lpf_timecourses[index] = butterworth(timecourses[index], high=1.0)
+    if apply_component_filter:
+        eig_mix = components['eig_mix'].copy()
+        timecourses = eig_mix.T
+        lpf_timecourses = np.zeros_like(timecourses)
+        for index in range(timecourses.shape[0]):
+            lpf_timecourses[index] = butterworth(timecourses[index], high=1.0)
+        output['eig_mix'] = lpf_timecourses.T
     
     output['masks'] = mask
     output['eig_vec'] = eig_vec
-    output['eig_mix'] = lpf_timecourses.T
-
+    
     return output
 
     # if blur:
diff --git a/seas/ica.py b/seas/ica.py
index 921c0de..895bb2e 100644
--- a/seas/ica.py
+++ b/seas/ica.py
@@ -304,6 +304,7 @@ def rebuild(components: dict,
             t_start: int = None,
             t_stop: int = None,
             apply_mean_filter: bool = True,
+            apply_component_filter: bool = False,
             apply_masked_mean: bool = False,
             filter_method: str = 'wavelet',
             fps: float = 7.5,
@@ -388,6 +389,15 @@ def rebuild(components: dict,
 
     eig_mix = components['eig_mix']
 
+    # Filter component timecourses
+    if apply_component_filter:
+        print('Filtering component timecourses using butterworth_lowpass at 0.5Hz...')
+        timecourses = eig_mix.T
+        lpf_timecourses = np.zeros_like(timecourses)
+        for index in range(timecourses.shape[0]):
+            lpf_timecourses[index] = butterworth(timecourses[index], high=0.5)
+        eig_mix = lpf_timecourses.T
+
     if (t_start == None):
         t_start = 0
 

From 592ae07a20ded7f4871becb9dac79f72cb6d8d83 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Thu, 24 Jul 2025 16:29:22 +1000
Subject: [PATCH 20/57] Make default component lpf 0.5Hz.

---
 seas/domains.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/seas/domains.py b/seas/domains.py
index 6d18113..fa736c2 100644
--- a/seas/domains.py
+++ b/seas/domains.py
@@ -751,7 +751,7 @@ def threshold_by_domains(components: dict,
         timecourses = eig_mix.T
         lpf_timecourses = np.zeros_like(timecourses)
         for index in range(timecourses.shape[0]):
-            lpf_timecourses[index] = butterworth(timecourses[index], high=1.0)
+            lpf_timecourses[index] = butterworth(timecourses[index], high=0.5)
         output['eig_mix'] = lpf_timecourses.T
     
     output['masks'] = mask

From a0e246fc38b61094cd7ffa7297e946f51b6aa6bd Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Thu, 24 Jul 2025 16:33:34 +1000
Subject: [PATCH 21/57] Debug masks assignment.

---
 seas/ica.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/seas/ica.py b/seas/ica.py
index 895bb2e..54ffd33 100644
--- a/seas/ica.py
+++ b/seas/ica.py
@@ -346,7 +346,6 @@ def rebuild(components: dict,
     roimask = components['roimask']
     shape = components['shape']
     mean = components['mean']
-    masks = components['masks']
     n_components = components['n_components']
     dtype = np.float32
 
@@ -419,6 +418,7 @@ def rebuild(components: dict,
                     eig_mix[t_start:t_stop, reconstruct_indices].T).T
 
     if apply_masked_mean:
+        masks = components['masks']
         assert masks is not None, \
         "Masks have not been assigned to dictionary"
         # Apply mean to masks only, zeroing unmasked pixels

From fdaa51f96a1ee2897456d0f391f616a4e3f26d5d Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Fri, 25 Jul 2025 14:59:27 +1000
Subject: [PATCH 22/57] Debug and readd combined morph filtering + blur for
 threshold_by_domains.

---
 seas/domains.py | 24 +++++++++++++++---------
 1 file changed, 15 insertions(+), 9 deletions(-)

diff --git a/seas/domains.py b/seas/domains.py
index fa736c2..9e0d121 100644
--- a/seas/domains.py
+++ b/seas/domains.py
@@ -628,7 +628,8 @@ def threshold_by_domains(components: dict,
                    min_size_ratio: float = 0.1,
                    map_only: bool = True,
                    apply_filter_mean: bool = True,
-                   apply_component_filter: bool = False,
+                   apply_component_lpf: bool = False,
+                   chigh: float = 0.5,
                    max_loops: int = 2,
                    ignore_small: bool = True,
                    thresh_type: str = 'max'):
@@ -732,26 +733,31 @@ def threshold_by_domains(components: dict,
 
             if roimask is not None:
                 eigenmask.flat[maskind] = mask.T[index]
+                # Remove small mask objects
                 filtered = remove_small_objects(eigenmask, min_size=min_mask_size, connectivity=1)
-                # filtered_float = filtered.astype(np.float32)
-                # eigenbrain.flat[maskind] = filtered_float.flat
-                # blurred = cv2.GaussianBlur(eigenbrain, (blur, blur), 0)
-                mask.T[index] = filtered.flat[maskind]
+                filtered_float = filtered.astype(np.float64)
+                eigenbrain.flat[maskind] = filtered_float.flat[maskind]
+                # Then blur
+                blurred = cv2.GaussianBlur(eigenbrain, (blur, blur), 0)
+                mask.T[index] = blurred.flat[maskind]
             else:
                 eigenbrain.flat = mask.T[index]
                 filtered = remove_small_objects(eigenbrain, min_size=min_mask_size, connectivity=1)
-                blurred = cv2.GaussianBlur(filtered, (blur, blur), 0)
+                filtered_float = filtered.astype(np.float64)
+                eigenbrain.flat[maskind] = filtered_float.flat
+                blurred = cv2.GaussianBlur(eigenbrain, (blur, blur), 0)
                 mask.T[index] = blurred.flat
 
-    eig_vec[~mask] = 0
+    mask_bool = mask.astype(bool)
+    eig_vec[~mask_bool] = 0
 
     # Filter component timecourses
-    if apply_component_filter:
+    if apply_component_lpf:
         eig_mix = components['eig_mix'].copy()
         timecourses = eig_mix.T
         lpf_timecourses = np.zeros_like(timecourses)
         for index in range(timecourses.shape[0]):
-            lpf_timecourses[index] = butterworth(timecourses[index], high=0.5)
+            lpf_timecourses[index] = butterworth(timecourses[index], high=chigh)
         output['eig_mix'] = lpf_timecourses.T
     
     output['masks'] = mask

From b2b7ad544bc1d4e5e1d102fcf0eafa4ffa67d7ef Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Fri, 25 Jul 2025 16:50:15 +1000
Subject: [PATCH 23/57] Added cfilt and mfilt parameter passthrough for
 seas.ica.rebuild()

---
 seas/ica.py | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/seas/ica.py b/seas/ica.py
index 54ffd33..a5e7b47 100644
--- a/seas/ica.py
+++ b/seas/ica.py
@@ -304,7 +304,10 @@ def rebuild(components: dict,
             t_start: int = None,
             t_stop: int = None,
             apply_mean_filter: bool = True,
+            mlow: float = 0.5,
+            mhigh: float = 1.0,
             apply_component_filter: bool = False,
+            chigh: float = 1.0,
             apply_masked_mean: bool = False,
             filter_method: str = 'wavelet',
             fps: float = 7.5,
@@ -394,7 +397,7 @@ def rebuild(components: dict,
         timecourses = eig_mix.T
         lpf_timecourses = np.zeros_like(timecourses)
         for index in range(timecourses.shape[0]):
-            lpf_timecourses[index] = butterworth(timecourses[index], high=0.5)
+            lpf_timecourses[index] = butterworth(timecourses[index], high=chigh)
         eig_mix = lpf_timecourses.T
 
     if (t_start == None):
@@ -425,7 +428,7 @@ def rebuild(components: dict,
         if apply_mean_filter:
             combined_mask = np.any(masks[:, reconstruct_indices], axis=1)
             mean_to_add = np.zeros_like(data_r)
-            mean_filtered = filter_mean(mean, filter_method, fps=fps)
+            mean_filtered = filter_mean(mean, filter_method, low_cutoff=mlow, high_cutoff=mhigh, fps=fps)
             mean_to_add[:, combined_mask] = mean_filtered[t_start:t_stop, None]
             data_r += mean_to_add
 
@@ -439,7 +442,7 @@ def rebuild(components: dict,
     else:
         # Run original readdition of mean
         if apply_mean_filter:
-            mean_filtered = filter_mean(mean, filter_method, fps=fps)
+            mean_filtered = filter_mean(mean, filter_method, low_cutoff=mlow, high_cutoff=mhigh, fps=fps)
             data_r += mean_filtered[t_start:t_stop, None]
 
         else:

From 6a438a0afd7f4a1da29a9bbf7512378b24c5a0f0 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Tue, 29 Jul 2025 18:47:47 +1000
Subject: [PATCH 24/57] Testing constant as filter_mean method.

---
 seas/ica.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/seas/ica.py b/seas/ica.py
index a5e7b47..418945b 100644
--- a/seas/ica.py
+++ b/seas/ica.py
@@ -543,6 +543,12 @@ def filter_mean(mean: np.ndarray,
         wavelet = waveletAnalysis(mean.astype('float64'), fps=fps)
         mean_filtered = wavelet.noiseFilter(upperPeriod=1 / low_cutoff)
 
+    elif filter_method == 'constant':
+        mean_template = np.zeros_like(mean)
+        meanest_mean = np.mean(mean)
+        mean_filtered = mean_template + meanest_mean
+        print('Mean set as constant: dfof = ' + str(meanest_mean))
+
     else:
         raise Exception("Filter method '" + str(filter_method)\
          + "' not supported!\n\t Supported methods: butterworth, butterworth_bandpass, wavelet")

From dc9f00e6a9c3a1d89580900dea1215026a0caa03 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Tue, 29 Jul 2025 20:47:35 +1000
Subject: [PATCH 25/57] Added z-score threshold to threshold_by_domains.

---
 seas/domains.py | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/seas/domains.py b/seas/domains.py
index 9e0d121..0806898 100644
--- a/seas/domains.py
+++ b/seas/domains.py
@@ -632,7 +632,8 @@ def threshold_by_domains(components: dict,
                    chigh: float = 0.5,
                    max_loops: int = 2,
                    ignore_small: bool = True,
-                   thresh_type: str = 'max'):
+                   thresh_type: str = 'max',
+                   thresh_param: float = None):
     '''
     Creates a domain map from extracted independent components.  A pixelwise maximum projection of the blurred signal components is taken through the n_components axis, to create a flattened representation of where a domain was maximally significant across the cortical surface.  Components with multiple noncontiguous significant regions are counted as two distinct domains.
 
@@ -706,6 +707,12 @@ def threshold_by_domains(components: dict,
             threshold_ROIs_vector = np.argmax(np.abs(eig_vec), axis=1)
             # Then threshold by clearing eig_vec outside of max indices
             mask[np.arange(eig_vec.shape[0]), threshold_ROIs_vector] = True
+        case 'z-score':
+            mean_ROIs_vector = np.nanmean(eig_vec, axis=1)
+            std_ROIs_vector = np.nanstd(eig_vec, axis=1)
+            z_ROIs_vector = (eig_vec - mean_ROIs_vector)/std_ROIs_vector
+            for i in np.arange(eig_vec.shape[0]):
+                mask[i, :] = np.abs(z_ROIs_vector[i]) > thresh_param
         case 'percentile':
             flipped = components['flipped']
             # Flip ICs where necessary using flipped from dict

From d52373bbe19245b5d22bf8bd38dc462990b91348 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Tue, 29 Jul 2025 21:21:01 +1000
Subject: [PATCH 26/57] Debug z-score for threshold_by_domains

---
 seas/domains.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/seas/domains.py b/seas/domains.py
index 0806898..47978b8 100644
--- a/seas/domains.py
+++ b/seas/domains.py
@@ -708,8 +708,8 @@ def threshold_by_domains(components: dict,
             # Then threshold by clearing eig_vec outside of max indices
             mask[np.arange(eig_vec.shape[0]), threshold_ROIs_vector] = True
         case 'z-score':
-            mean_ROIs_vector = np.nanmean(eig_vec, axis=1)
-            std_ROIs_vector = np.nanstd(eig_vec, axis=1)
+            mean_ROIs_vector = np.nanmean(eig_vec, axis=0)
+            std_ROIs_vector = np.nanstd(eig_vec, axis=0)
             z_ROIs_vector = (eig_vec - mean_ROIs_vector)/std_ROIs_vector
             for i in np.arange(eig_vec.shape[0]):
                 mask[i, :] = np.abs(z_ROIs_vector[i]) > thresh_param

From 3ab58bfb9943304bb016a272cbd65e6d571d54cb Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Wed, 30 Jul 2025 16:59:44 +1000
Subject: [PATCH 27/57] Cleaned and consolidated apply_component_filter option
 for rebuilding. Consolidated new options in docstrings.

---
 seas/ica.py | 21 ++++++++++++++++-----
 1 file changed, 16 insertions(+), 5 deletions(-)

diff --git a/seas/ica.py b/seas/ica.py
index 418945b..10bdee3 100644
--- a/seas/ica.py
+++ b/seas/ica.py
@@ -309,7 +309,7 @@ def rebuild(components: dict,
             apply_component_filter: bool = False,
             chigh: float = 1.0,
             apply_masked_mean: bool = False,
-            filter_method: str = 'wavelet',
+            filter_method: str = 'butterworth_highpass',
             fps: float = 7.5,
             include_noise: bool = True):
     '''
@@ -329,8 +329,20 @@ def rebuild(components: dict,
             The frame to stop rebuilding the movie at.  If none is provided, the rebuilt movie ends at the last frame
         apply_mean_filter:
             Whether to apply a filter to the mean signal.
-        filter_method:;
-            The filter method to apply (see filter_mean function).
+        mlow:
+            A float determining the highpass cutoff for the mean filter, if used.
+        mhigh:
+            A float determining the lowpass cutoff for the mean filter, if used.
+        apply_component_filter:
+            Whether to apply a butterworth_lowpass filter to IC timecourses before rebuild.
+        chigh:
+            A float determining the lowpass cutoff for the component filter, if used.
+        apply_masked_mean:
+            If True, only re-adds the mean signal to pixels where at least one IC is defined. To be used for thresholded ICs.
+        filter_method:
+            The filter method to apply to the mean. Choose from 'butterworth_bandpass', 'butterworth_lowpass', 'butterworth_highpass', or 'constant'. Behaviour for 'wavelet' as yet undefined.
+        fps:
+            A float determining the fps for the source video.
         include_noise:
             Whether to include noise components when rebuilding.  If noise_components should not be included in the rebuilt movie, set this to False
 
@@ -389,11 +401,10 @@ def rebuild(components: dict,
         assert eig_vec[:,0].size == maskind[0].size, \
         "Eigenvector size is not compatible with the masked region's size"
 
-    eig_mix = components['eig_mix']
-
     # Filter component timecourses
     if apply_component_filter:
         print('Filtering component timecourses using butterworth_lowpass at 0.5Hz...')
+        eig_mix = components['eig_mix']
         timecourses = eig_mix.T
         lpf_timecourses = np.zeros_like(timecourses)
         for index in range(timecourses.shape[0]):

From 0ade1919fd3560e764588b47b84661fdd3749b8e Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Wed, 30 Jul 2025 17:00:21 +1000
Subject: [PATCH 28/57] Consolidated docstrings for threshold_by_domains, and
 removed component filtering (now in seas.ica.rebuild).

---
 seas/domains.py | 58 ++++++++++++++-----------------------------------
 1 file changed, 16 insertions(+), 42 deletions(-)

diff --git a/seas/domains.py b/seas/domains.py
index 47978b8..82ba502 100644
--- a/seas/domains.py
+++ b/seas/domains.py
@@ -625,49 +625,31 @@ def write_frame(frame):
 def threshold_by_domains(components: dict,
                    blur: int = 1,
                    min_mask_size: int = 64,
-                   min_size_ratio: float = 0.1,
-                   map_only: bool = True,
-                   apply_filter_mean: bool = True,
-                   apply_component_lpf: bool = False,
-                   chigh: float = 0.5,
-                   max_loops: int = 2,
-                   ignore_small: bool = True,
                    thresh_type: str = 'max',
                    thresh_param: float = None):
     '''
-    Creates a domain map from extracted independent components.  A pixelwise maximum projection of the blurred signal components is taken through the n_components axis, to create a flattened representation of where a domain was maximally significant across the cortical surface.  Components with multiple noncontiguous significant regions are counted as two distinct domains.
+    Function based on modified get_domain_map(). Thresholds ICs using a variety of methods for selective rebuild.
 
     Arguments:
         components: 
-            The dictionary of components returned from seas.ica.project.  Domains are most interesting if artifacts has already been assigned through seas.gui.run_gui.
+            The dictionary of components returned from seas.ica.project.  ROIs are most interesting if artifacts has already been assigned through seas.gui.run_gui.
         blur: 
-            An odd integer kernel Gaussian blur to run before segmenting.  Domains look smoother with larger blurs, but you can lose some smaller domains.
-        map_only:
-            If true, compute the map only, do not rebuild time courses under each domain.
-        apply_filter_mean:
-            Whether to compute the filtered mean when calculating ROI rebuild timecourses.
-        min_size_ratio:
-            The minimum size ratio of the mean component size to allow for a component.  If a the size of a component is under (min_size_ratio x mean_domain_size), and the next most significant domain over the pixel would result in a larger size domain, this next domain is chosen.
-        max_loops:
-            The number of times to check if the next most significant domain would result in a larger domain size.  To entirely disable this, set max_loops to 0.
-        ignore_small:
-            If True, assign undersize domains that were not reassigned during max_loops to np.nan.
+            An odd integer kernel Gaussian blur to run before segmenting.  ROIs look smoother with larger blurs, but you can lose some smaller domains.
+        min_mask_size:
+            An integer determining the minimum ROIs passed from each thresholded IC.
+        thresh_type:
+            A string used to determine IC threshold method. Choose from either 'max', 'z-score' or 'percentile'.
+        thresh_param:
+            A float used to determine the parameter for the given thresh_type. For 'z-score', this is the z-score threshold (eg; 2.0 for 2std). For 'percentile' this is the percentile used to threshold (eg; 95th percentile = 0.95).
 
     Returns:
         output: a dictionary containing the results of the operation, containing the following keys
             domain_blur:
                 The Gaussian blur value used when generating the map
-            component_assignment: 
-                A map showing the index of which *component* was maximally significant over a given pixel.  Here, 
-                This is in contrast to the domain map, where each domain is a unique integer.  
-            domain_ROIs: 
-                The computed np.array domain map (x,y).  Each domain is represented by a unique integer, and represents a discrete continuous unit.  Values that are masked, or where large enough domains were not detected are set to np.nan.
-
-        if not map_only, the following are also included in the output dictionary:
-            ROI_timecourses: 
-                The time courses rebuilt from the video under each ROI.  The frame mean is not included in this calculation, and must be re-added from mean_filtered.
-            mean_filtered: 
-                The frame mean, filtered by the default method.
+            eig_vec: 
+                The thresholded eigenvectors (ICs).  
+            thresh_masks: 
+                The boolean masks used to threshold eig_vec.
     '''
     print('\nExtracting Domain ROIs\n-----------------------')
     output = {}
@@ -718,7 +700,7 @@ def threshold_by_domains(components: dict,
             # Flip ICs where necessary using flipped from dict
             flipped_threshold_vec = np.multiply(flipped, eig_vec)
             # Calculate 95 percentile cutoff for each IC
-            cutoff_vector = np.percentile(flipped, 0.95, axis=1)
+            cutoff_vector = np.percentile(flipped, thresh_param, axis=0)
             # Mask for all values above cutoff
             for i in np.arange(eig_vec.shape[0]):
                 mask[i, :] = flipped_threshold_vec[i] > cutoff_vector[i]
@@ -757,17 +739,9 @@ def threshold_by_domains(components: dict,
 
     mask_bool = mask.astype(bool)
     eig_vec[~mask_bool] = 0
-
-    # Filter component timecourses
-    if apply_component_lpf:
-        eig_mix = components['eig_mix'].copy()
-        timecourses = eig_mix.T
-        lpf_timecourses = np.zeros_like(timecourses)
-        for index in range(timecourses.shape[0]):
-            lpf_timecourses[index] = butterworth(timecourses[index], high=chigh)
-        output['eig_mix'] = lpf_timecourses.T
     
-    output['masks'] = mask
+    output['thresh_masks'] = mask
+    # output['thresh_vec'] = eig_vec
     output['eig_vec'] = eig_vec
     
     return output

From 8b5113e6310be7f9c82e0133da093b391ec83995 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Thu, 7 Aug 2025 11:22:12 +1000
Subject: [PATCH 29/57] Refactor filtering ICs into its own function
 seas.ica.filter_components().

---
 seas/ica.py | 36 ++++++++++++++++++++++++++++++------
 1 file changed, 30 insertions(+), 6 deletions(-)

diff --git a/seas/ica.py b/seas/ica.py
index 10bdee3..23266b7 100644
--- a/seas/ica.py
+++ b/seas/ica.py
@@ -403,13 +403,9 @@ def rebuild(components: dict,
 
     # Filter component timecourses
     if apply_component_filter:
-        print('Filtering component timecourses using butterworth_lowpass at 0.5Hz...')
         eig_mix = components['eig_mix']
-        timecourses = eig_mix.T
-        lpf_timecourses = np.zeros_like(timecourses)
-        for index in range(timecourses.shape[0]):
-            lpf_timecourses[index] = butterworth(timecourses[index], high=chigh)
-        eig_mix = lpf_timecourses.T
+        lpf_eig_mix = filter_components(eig_mix, fps=fps, high_cutoff=chigh)
+        eig_mix = lpf_eig_mix
 
     if (t_start == None):
         t_start = 0
@@ -567,6 +563,34 @@ def filter_mean(mean: np.ndarray,
     return mean_filtered
 
 
+def filter_components(eig_mix: np.ndarray,
+                      fps: float = 7.5,
+                      high_cutoff: float = 0.5):
+    '''
+    Applies a butterworth low pass filter to the ica component timecourses.
+
+    Arguments:
+        eig_mix: 
+            The mixing matrix containing IC timecourses.
+        fps:
+            Sampling rate of the video.
+        high_cutoff:
+            The frequency cutoff to apply the low pass filter at.
+
+    Returns:
+        lpf_eig_mix: The filtered IC timecourses reconstructed as the eig_mix matrix.
+    '''
+    print('Filtering component timecourses using butterworth_lowpass at '+ high_cutoff +'Hz...')
+
+    timecourses = eig_mix.T
+    lpf_timecourses = np.zeros_like(timecourses)
+    for index in range(timecourses.shape[0]):
+        lpf_timecourses[index] = butterworth(timecourses[index], fps=fps, high=high_cutoff)
+    lpf_eig_mix = lpf_timecourses.T
+
+    return lpf_eig_mix
+
+
 def rebuild_mean_roi_timecourse(components: np.ndarray,
                                 mask: np.ndarray,
                                 include_zero: bool = True,

From 43af207204f1b1d0821dd87bfe89eb1571d5ce3e Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Thu, 7 Aug 2025 11:42:49 +1000
Subject: [PATCH 30/57] dtype debug

---
 seas/ica.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/seas/ica.py b/seas/ica.py
index 23266b7..e588484 100644
--- a/seas/ica.py
+++ b/seas/ica.py
@@ -580,7 +580,7 @@ def filter_components(eig_mix: np.ndarray,
     Returns:
         lpf_eig_mix: The filtered IC timecourses reconstructed as the eig_mix matrix.
     '''
-    print('Filtering component timecourses using butterworth_lowpass at '+ high_cutoff +'Hz...')
+    print('Filtering component timecourses using butterworth_lowpass at '+ str(high_cutoff) +'Hz...')
 
     timecourses = eig_mix.T
     lpf_timecourses = np.zeros_like(timecourses)

From 86daec88fca74061e04d29653436346e827a735b Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Wed, 13 Aug 2025 17:48:55 +1000
Subject: [PATCH 31/57] Debug eig_mix assignment order.

---
 seas/ica.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/seas/ica.py b/seas/ica.py
index 10bdee3..7a57806 100644
--- a/seas/ica.py
+++ b/seas/ica.py
@@ -358,6 +358,7 @@ def rebuild(components: dict,
     assert type(components) is dict, 'Components were not in format expected'
 
     eig_vec = components['eig_vec']
+    eig_mix = components['eig_mix']
     roimask = components['roimask']
     shape = components['shape']
     mean = components['mean']
@@ -404,7 +405,6 @@ def rebuild(components: dict,
     # Filter component timecourses
     if apply_component_filter:
         print('Filtering component timecourses using butterworth_lowpass at 0.5Hz...')
-        eig_mix = components['eig_mix']
         timecourses = eig_mix.T
         lpf_timecourses = np.zeros_like(timecourses)
         for index in range(timecourses.shape[0]):

From 3a65829ef2ee5427f6004c6d3e9b6be80cd1c185 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Mon, 1 Sep 2025 12:21:55 +1000
Subject: [PATCH 32/57] Update all instances of np.NAN to np.nan for numpy 2.x
 compatability.

---
 seas/domains.py | 8 ++++----
 seas/gui.py     | 2 +-
 seas/ica.py     | 4 ++--
 3 files changed, 7 insertions(+), 7 deletions(-)

diff --git a/seas/domains.py b/seas/domains.py
index 82ba502..27e2f52 100644
--- a/seas/domains.py
+++ b/seas/domains.py
@@ -97,7 +97,7 @@ def get_domain_map(components: dict,
             blur += 1
 
         eigenbrain = np.empty(shape)
-        eigenbrain[:] = np.NAN
+        eigenbrain[:] = np.nan
 
         for index in range(eig_vec.shape[1]):
 
@@ -117,7 +117,7 @@ def get_domain_map(components: dict,
 
     if roimask is not None:
         domain_ROIs = np.empty(shape)
-        domain_ROIs[:] = np.NAN
+        domain_ROIs[:] = np.nan
         domain_ROIs.flat[maskind] = domain_ROIs_vector
 
     else:
@@ -716,7 +716,7 @@ def threshold_by_domains(components: dict,
 
         eigenmask = np.zeros(shape, dtype=bool)
         eigenbrain = np.empty(shape)
-        eigenbrain[:] = np.NAN
+        eigenbrain[:] = np.nan
 
         for index in range(mask.shape[1]):
 
@@ -751,7 +751,7 @@ def threshold_by_domains(components: dict,
 
     # if roimask is not None:
     #     domain_ROIs = np.empty(shape)
-    #     domain_ROIs[:] = np.NAN
+    #     domain_ROIs[:] = np.nan
     #     domain_ROIs.flat[maskind] = domain_ROIs_vector
 
     # else:
diff --git a/seas/gui.py b/seas/gui.py
index 76b214a..ec70389 100644
--- a/seas/gui.py
+++ b/seas/gui.py
@@ -469,7 +469,7 @@ def update(self, component_id):
 
                     if component_id is None:  # Clear image.
                         im = np.empty((x, y))
-                        im[:] = np.NAN
+                        im[:] = np.nan
                         self.imgplot = self.ax.imshow(im)
                         self.canvas.draw()
                         return ()
diff --git a/seas/ica.py b/seas/ica.py
index dd9df93..2fbde55 100644
--- a/seas/ica.py
+++ b/seas/ica.py
@@ -713,7 +713,7 @@ def rebuild_eigenbrain(eig_vec: np.ndarray,
         else:
             eigenbrains = np.empty(
                 (roimask.shape[0], roimask.shape[1], eig_vec.shape[1]))
-            eigenbrains[:] = np.NAN
+            eigenbrains[:] = np.nan
             eigenbrains[x, y, :] = eig_vec
         eigenbrains = np.swapaxes(eigenbrains, 0, 2)
         eigenbrains = np.swapaxes(eigenbrains, 1, 2)
@@ -730,7 +730,7 @@ def rebuild_eigenbrain(eig_vec: np.ndarray,
             eigenbrain = eigenbrain.reshape(eigb_shape)
         else:
             eigenbrain = np.empty(roimask.shape)
-            eigenbrain[:] = np.NAN
+            eigenbrain[:] = np.nan
             eigenbrain.flat[maskind] = eig_vec.T[index]
 
         return eigenbrain

From 2ae6341aac783ef1e4bb2a01411b96f7b5b2a851 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Wed, 5 Nov 2025 14:04:14 +1000
Subject: [PATCH 33/57] Debug include_noise logic in ica.rebuild()

---
 seas/ica.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/seas/ica.py b/seas/ica.py
index 2fbde55..b1a9bf3 100644
--- a/seas/ica.py
+++ b/seas/ica.py
@@ -376,7 +376,8 @@ def rebuild(components: dict,
     elif artifact_components == 'none':
         print('including all components')
         artifact_components = np.zeros(n_components)
-    elif ((not include_noise) and ('noise_components' in components.keys())):
+    
+    if ((not include_noise) and ('noise_components' in components.keys())):
         print('Not rebuilding noise components')
         artifact_components += components['noise_components']
         artifact_components[np.where(artifact_components > 1)] = 1

From 974398231dd40b38b9052427b9ea680de096ebdc Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Tue, 18 Nov 2025 13:22:48 +1000
Subject: [PATCH 34/57] Debug use of thresh_masks in ica rebuild.

---
 seas/ica.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/seas/ica.py b/seas/ica.py
index b1a9bf3..b3f4dd7 100644
--- a/seas/ica.py
+++ b/seas/ica.py
@@ -429,7 +429,7 @@ def rebuild(components: dict,
                     eig_mix[t_start:t_stop, reconstruct_indices].T).T
 
     if apply_masked_mean:
-        masks = components['masks']
+        masks = components['thresh_masks']
         assert masks is not None, \
         "Masks have not been assigned to dictionary"
         # Apply mean to masks only, zeroing unmasked pixels

From 2faadf8bbe1ce48686fcd13734b19b3d20546528 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Tue, 18 Nov 2025 14:05:40 +1000
Subject: [PATCH 35/57] Added method to threshold component timeseries' and
 associated option in rebuild.

---
 seas/ica.py | 41 +++++++++++++++++++++++++++++++++++++++--
 1 file changed, 39 insertions(+), 2 deletions(-)

diff --git a/seas/ica.py b/seas/ica.py
index b3f4dd7..47fb29a 100644
--- a/seas/ica.py
+++ b/seas/ica.py
@@ -308,6 +308,8 @@ def rebuild(components: dict,
             mhigh: float = 1.0,
             apply_component_filter: bool = False,
             chigh: float = 1.0,
+            apply_component_threshold: bool = False,
+            cthresh: float = 2.0,
             apply_masked_mean: bool = False,
             filter_method: str = 'butterworth_highpass',
             fps: float = 7.5,
@@ -337,6 +339,10 @@ def rebuild(components: dict,
             Whether to apply a butterworth_lowpass filter to IC timecourses before rebuild.
         chigh:
             A float determining the lowpass cutoff for the component filter, if used.
+        apply_component_threshold:
+            Whether to apply a z-score threshold on the component timeseries.
+        cthresh:
+            A float determining the z-score threshold for the component threshold, if used.
         apply_masked_mean:
             If True, only re-adds the mean signal to pixels where at least one IC is defined. To be used for thresholded ICs.
         filter_method:
@@ -408,6 +414,11 @@ def rebuild(components: dict,
         lpf_eig_mix = filter_components(eig_mix, fps=fps, high_cutoff=chigh)
         eig_mix = lpf_eig_mix
 
+    # Threshold component timecourses
+    if apply_component_threshold:
+        thresh_eig_mix = threshold_components(eig_mix, thresh_param=cthresh)
+        eig_mix = thresh_eig_mix
+
     if (t_start == None):
         t_start = 0
 
@@ -568,7 +579,7 @@ def filter_components(eig_mix: np.ndarray,
                       fps: float = 7.5,
                       high_cutoff: float = 0.5):
     '''
-    Applies a butterworth low pass filter to the ica component timecourses.
+    Applies a butterworth low pass filter to the IC timecourses.
 
     Arguments:
         eig_mix: 
@@ -581,8 +592,8 @@ def filter_components(eig_mix: np.ndarray,
     Returns:
         lpf_eig_mix: The filtered IC timecourses reconstructed as the eig_mix matrix.
     '''
+    
     print('Filtering component timecourses using butterworth_lowpass at '+ str(high_cutoff) +'Hz...')
-
     timecourses = eig_mix.T
     lpf_timecourses = np.zeros_like(timecourses)
     for index in range(timecourses.shape[0]):
@@ -591,6 +602,32 @@ def filter_components(eig_mix: np.ndarray,
 
     return lpf_eig_mix
 
+def threshold_components(eig_mix: np.ndarray,
+                         thresh_param: float):
+    '''
+    Applies a z-score threshold to the IC timecourses.
+
+    Arguments:
+        eig_mix: 
+            The mixing matrix containing IC timecourses.
+        thresh_param:
+            Z-score thresholding parameter (standard deviations).
+
+    Returns:
+        thresh_eig_mix: The thresholded IC timecourses reconstructed as the eig_mix matrix.
+    '''
+
+    print('Thresholding component timecourses using z-score: >' + str(thresh_param) +'s.d.')
+    timecourses = eig_mix.T
+    thresh_timecourses = np.zeros_like(timecourses)
+    for index in range(timecourses.shape[0]):
+        mean = np.mean(timecourses[index])
+        std = np.std(timecourses[index])
+        thresh_timecourses = timecourses[index][timecourses > mean + thresh_param*std]
+    thresh_eig_mix = thresh_timecourses.T
+
+    return thresh_eig_mix
+
 
 def rebuild_mean_roi_timecourse(components: np.ndarray,
                                 mask: np.ndarray,

From c3e956b1adb33fa2105e13a3c05eb8967e416b06 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Tue, 18 Nov 2025 22:11:06 +1000
Subject: [PATCH 36/57] Added argument ica.project to optionally retain excess
 noise.

---
 seas/ica.py | 36 ++++++++++++++++++++----------------
 1 file changed, 20 insertions(+), 16 deletions(-)

diff --git a/seas/ica.py b/seas/ica.py
index 47fb29a..a8c3d08 100644
--- a/seas/ica.py
+++ b/seas/ica.py
@@ -17,6 +17,7 @@ def project(vector: np.ndarray,
             shape: Tuple[int, int, int],
             roimask: np.ndarray = None,
             n_components: int = None,
+            crop_excess_noise: bool = True,
             svd_multiplier: float = 5,
             calc_residuals: bool = True,
             max_iter: int = 1000):
@@ -183,23 +184,26 @@ def project(vector: np.ndarray,
         reduced_n_components = int((noise.size - noise.sum()) * 1.25)
 
         print('reduced_n_components:', reduced_n_components)
-
-        if reduced_n_components < n_components:
-            print('Cropping', n_components, 'to', reduced_n_components)
-
-            ev_sort = np.argsort(eig_mix.std(axis=0))
-            eig_vec = eig_vec[:, ev_sort][:, ::-1]
-            eig_mix = eig_mix[:, ev_sort][:, ::-1]
-            noise = noise[ev_sort][::-1]
-
-            eig_vec = eig_vec[:, :reduced_n_components]
-            eig_mix = eig_mix[:, :reduced_n_components]
-            n_components = reduced_n_components
-            noise = noise[:reduced_n_components]
-
-            components['lag1_full'] = components['lag1_full'][ev_sort][::-1]
+        
+        if crop_excess_noise:
+            if reduced_n_components < n_components:
+                print('Cropping', n_components, 'to', reduced_n_components)
+
+                ev_sort = np.argsort(eig_mix.std(axis=0))
+                eig_vec = eig_vec[:, ev_sort][:, ::-1]
+                eig_mix = eig_mix[:, ev_sort][:, ::-1]
+                noise = noise[ev_sort][::-1]
+
+                eig_vec = eig_vec[:, :reduced_n_components]
+                eig_mix = eig_mix[:, :reduced_n_components]
+                n_components = reduced_n_components
+                noise = noise[:reduced_n_components]
+
+                components['lag1_full'] = components['lag1_full'][ev_sort][::-1]
+            else:
+                print('Less than 75% signal.  Not cropping excess noise.')
         else:
-            print('Less than 75% signal.  Not cropping excess noise.')
+            print('Noise retention enabled. Not cropping excess noise.')
 
         components['noise_components'] = noise
         components['cutoff'] = cutoff

From 93a38bdcc794628194afeb1f6f1b1b2c374a15bc Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Tue, 18 Nov 2025 22:25:04 +1000
Subject: [PATCH 37/57] Debug calc_residuals fork in ica.project.

---
 seas/ica.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/seas/ica.py b/seas/ica.py
index a8c3d08..098c434 100644
--- a/seas/ica.py
+++ b/seas/ica.py
@@ -264,8 +264,7 @@ def project(vector: np.ndarray,
         try:
             vector = vector.astype('float64')
             rebuilt = rebuild(components,
-                              artifact_components='none',
-                              vector=True).T
+                              apply_mean_filter=False).T
 
             rebuilt -= rebuilt.mean(axis=0)
             vector -= vector.mean(axis=0)

From 190114cc3ea1d8a54a6acdc5088933b43b1ecf2c Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Tue, 18 Nov 2025 22:27:20 +1000
Subject: [PATCH 38/57] Replace missing artifact_components argument in
 calc_residuals fork.

---
 seas/ica.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/seas/ica.py b/seas/ica.py
index 098c434..b43ddf6 100644
--- a/seas/ica.py
+++ b/seas/ica.py
@@ -264,6 +264,7 @@ def project(vector: np.ndarray,
         try:
             vector = vector.astype('float64')
             rebuilt = rebuild(components,
+                              artifact_components='none',
                               apply_mean_filter=False).T
 
             rebuilt -= rebuilt.mean(axis=0)

From 0cc2a4058108caa92cbfde326b38953c419f8b23 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Sat, 22 Nov 2025 13:39:35 +1000
Subject: [PATCH 39/57] Add new functions to export event mask compilation tifs
 and thresholded video.

---
 seas/domains.py    |   1 -
 seas/experiment.py |  56 +++++++++++++++++-
 seas/ica.py        | 137 +++++++++++++++++++++++++++++++++++++++++++--
 3 files changed, 186 insertions(+), 8 deletions(-)

diff --git a/seas/domains.py b/seas/domains.py
index 27e2f52..fd7a913 100644
--- a/seas/domains.py
+++ b/seas/domains.py
@@ -621,7 +621,6 @@ def write_frame(frame):
     print('Saving Colorbar to:' + cbarpath)
     save_colorbar(scale, cbarpath, colormap=colormap)
 
-
 def threshold_by_domains(components: dict,
                    blur: int = 1,
                    min_mask_size: int = 64,
diff --git a/seas/experiment.py b/seas/experiment.py
index 85124dd..6cea737 100644
--- a/seas/experiment.py
+++ b/seas/experiment.py
@@ -8,11 +8,12 @@
 from seas.filemanager import sort_experiments, get_exp_span_string, read_yaml
 from seas.rois import roi_loader, make_mask, get_masked_region, insert_masked_region, draw_bounding_box
 from seas.hdf5manager import hdf5manager
-from seas.ica import project, filter_mean
+from seas.ica import project, filter_mean, rebuild_eigenbrain, threshold_by_domains, filter_components, threshold_components, rebuild
 from seas.signalanalysis import sort_noise, lag_n_autocorr
 from seas.waveletAnalysis import waveletAnalysis
 
 from typing import List
+import tifffile as tif
 
 
 class Experiment:
@@ -112,7 +113,7 @@ def __init__(self,
 
         if np.any(np.isnan(movie)):
             # If the video was already masked
-            roimask = np.zeros(movie[0].shape, dtype='uisnt8')
+            roimask = np.zeros(movie[0].shape, dtype='uint8')
             roimask[np.where(~np.isnan(movie[0]))] = 1
             self.roimask = roimask
 
@@ -461,3 +462,54 @@ def ica_project(self,
             f.print()
 
         return components
+
+def export_event_masks(components: dict,
+                       outpath: str,
+                       blur: int = 3,
+                       thresh_type: str = 'z-score',
+                       thresh_param: float = 7) -> None:
+    components_copy = components.copy()
+    threshold = threshold_by_domains(components_copy, 
+                                     blur = blur, 
+                                     thresh_type = thresh_type, 
+                                     thresh_param = thresh_param)
+    components_copy.update(threshold)
+    artifacts_bool = components_copy['artifact_components'].astype(np.bool)
+    event_components = components_copy['eig_vec'][:, ~artifacts_bool]
+    event_masks = rebuild_eigenbrain(event_components,
+                                     roimask = components_copy['roimask'], 
+                                     bulk = True)
+    event_masks = np.where(np.abs(event_masks) > 0, 255, 0)
+    event_masks = np.where(np.mean(event_masks, axis = 0) == 255, 0, event_masks)
+    tif.imwrite(outpath, event_masks.astype(np.float32),imagej=True)
+
+def export_event_video(components: dict,
+                       outpath: str,
+                       artifact_components: np.ndarray = None,
+                       t_start: int = None,
+                       t_stop: int = None,
+                       apply_mean_filter: bool = True,
+                       cthresh: float = 2.0,
+                       apply_masked_mean: bool = False,
+                       filter_method: str = 'constant',
+                       include_noise: bool = True) -> None:
+    components_copy = components.copy()
+    threshold = threshold_by_domains(components_copy, 
+                                     blur = 3, 
+                                     thresh_type = 'z-score', 
+                                     thresh_param = 7)
+    eig_mix = filter_components(components_copy['eig_mix'])
+    eig_mix = threshold_components(eig_mix, 
+                                   thresh_param = cthresh)
+    components_copy.update(threshold)
+    components_copy['eig_mix'] = eig_mix
+    rebuilt = rebuild(components_copy,
+                      artifact_components = artifact_components,
+                      t_start = t_start, 
+                      t_stop = t_stop,
+                      apply_mean_filter = apply_mean_filter,
+                      cthresh = cthresh,
+                      apply_masked_mean = apply_masked_mean,
+                      filter_method = filter_method,
+                      include_noise = include_noise)
+    tif.imwrite(outpath, rebuilt.astype(np.float32), imagej=True)
\ No newline at end of file
diff --git a/seas/ica.py b/seas/ica.py
index b43ddf6..4e2b404 100644
--- a/seas/ica.py
+++ b/seas/ica.py
@@ -12,6 +12,9 @@
 from seas.hdf5manager import hdf5manager
 from seas.video import rotate, save, rescale, play, scale_video
 
+import cv2
+from skimage.morphology import remove_small_objects
+
 
 def project(vector: np.ndarray,
             shape: Tuple[int, int, int],
@@ -302,7 +305,6 @@ def project(vector: np.ndarray,
     print('\n')
     return components
 
-
 def rebuild(components: dict,
             artifact_components: np.ndarray = None,
             t_start: int = None,
@@ -625,13 +627,139 @@ def threshold_components(eig_mix: np.ndarray,
     timecourses = eig_mix.T
     thresh_timecourses = np.zeros_like(timecourses)
     for index in range(timecourses.shape[0]):
-        mean = np.mean(timecourses[index])
-        std = np.std(timecourses[index])
-        thresh_timecourses = timecourses[index][timecourses > mean + thresh_param*std]
+        timecourse = timecourses[index]
+        mean = np.mean(timecourse)
+        std = np.std(timecourse)
+        threshold = mean + thresh_param*std
+        timecourse[np.abs(timecourse) < np.abs(threshold)] = 0
+        thresh_timecourses[index] = timecourse
     thresh_eig_mix = thresh_timecourses.T
 
     return thresh_eig_mix
 
+def threshold_by_domains(components: dict,
+                   blur: int = 1,
+                   min_mask_size: int = 64,
+                   thresh_type: str = 'max',
+                   thresh_param: float = None):
+    '''
+    Function based on modified get_domain_map(). Thresholds ICs using a variety of methods for selective rebuild.
+
+    Arguments:
+        components: 
+            The dictionary of components returned from seas.ica.project.  ROIs are most interesting if artifacts has already been assigned through seas.gui.run_gui.
+        blur: 
+            An odd integer kernel Gaussian blur to run before segmenting.  ROIs look smoother with larger blurs, but you can lose some smaller domains.
+        min_mask_size:
+            An integer determining the minimum ROIs passed from each thresholded IC.
+        thresh_type:
+            A string used to determine IC threshold method. Choose from either 'max', 'z-score' or 'percentile'.
+        thresh_param:
+            A float used to determine the parameter for the given thresh_type. For 'z-score', this is the z-score threshold (eg; 2.0 for 2std). For 'percentile' this is the percentile used to threshold (eg; 95th percentile = 0.95).
+
+    Returns:
+        output: a dictionary containing the results of the operation, containing the following keys
+            domain_blur:
+                The Gaussian blur value used when generating the map
+            eig_vec: 
+                The thresholded eigenvectors (ICs).  
+            thresh_masks: 
+                The boolean masks used to threshold eig_vec.
+    '''
+    print('\nExtracting Domain ROIs\n-----------------------')
+    output = {}
+    output['domain_blur'] = blur
+
+    eig_vec = components['eig_vec'].copy()
+
+    shape = components['shape']
+    shape = (shape[1], shape[2])
+
+    if 'roimask' in components.keys() and components['roimask'] is not None:
+        roimask = components['roimask']
+        maskind = np.where(roimask.flat == 1)[0]
+    else:
+        roimask = None
+
+    if 'artifact_components' in components.keys():
+        artifact_components = components['artifact_components']
+
+        print('Switching to signal indices only for domain detection')
+
+        if 'noise_components' in components.keys():
+            noise_components = components['noise_components']
+
+            signal_indices = np.where((artifact_components +
+                                       noise_components) == 0)[0]
+        else:
+            print('no noise components found')
+            signal_indices = np.where(artifact_components == 0)[0]
+        # eig_vec = eig_vec[:, signal_indices] # Don't change number of ICs, we're updating back to dict
+    
+    mask = np.zeros_like(eig_vec, dtype=bool)
+
+    match thresh_type:
+        case 'max':
+            # Return indices across each eig_vec (loading vector for component) where loading is max
+            threshold_ROIs_vector = np.argmax(np.abs(eig_vec), axis=1)
+            # Then threshold by clearing eig_vec outside of max indices
+            mask[np.arange(eig_vec.shape[0]), threshold_ROIs_vector] = True
+        case 'z-score':
+            mean_ROIs_vector = np.nanmean(eig_vec, axis=0)
+            std_ROIs_vector = np.nanstd(eig_vec, axis=0)
+            z_ROIs_vector = (eig_vec - mean_ROIs_vector)/std_ROIs_vector
+            for i in np.arange(eig_vec.shape[0]):
+                mask[i, :] = np.abs(z_ROIs_vector[i]) > thresh_param
+        case 'percentile':
+            flipped = components['flipped']
+            # Flip ICs where necessary using flipped from dict
+            flipped_threshold_vec = np.multiply(flipped, eig_vec)
+            # Calculate 95 percentile cutoff for each IC
+            cutoff_vector = np.percentile(flipped, thresh_param, axis=0)
+            # Mask for all values above cutoff
+            for i in np.arange(eig_vec.shape[0]):
+                mask[i, :] = flipped_threshold_vec[i] > cutoff_vector[i]
+        case _:
+            print("Threshold type is neither max nor percentile.")
+
+    # Filter small mask ROIs and smooth using blur
+    if blur:
+        print('blurring domains...')
+        assert type(blur) is int, 'blur was not valid'
+        if blur % 2 != 1:
+            blur += 1
+
+        eigenmask = np.zeros(shape, dtype=bool)
+        eigenbrain = np.empty(shape)
+        eigenbrain[:] = np.nan
+
+        for index in range(mask.shape[1]):
+
+            if roimask is not None:
+                eigenmask.flat[maskind] = mask.T[index]
+                # Remove small mask objects
+                filtered = remove_small_objects(eigenmask, min_size=min_mask_size, connectivity=1)
+                filtered_float = filtered.astype(np.float64)
+                eigenbrain.flat[maskind] = filtered_float.flat[maskind]
+                # Then blur
+                blurred = cv2.GaussianBlur(eigenbrain, (blur, blur), 0)
+                mask.T[index] = blurred.flat[maskind]
+            else:
+                eigenbrain.flat = mask.T[index]
+                filtered = remove_small_objects(eigenbrain, min_size=min_mask_size, connectivity=1)
+                filtered_float = filtered.astype(np.float64)
+                eigenbrain.flat[maskind] = filtered_float.flat
+                blurred = cv2.GaussianBlur(eigenbrain, (blur, blur), 0)
+                mask.T[index] = blurred.flat
+
+    mask_bool = mask.astype(bool)
+    eig_vec[~mask_bool] = 0
+    
+    output['thresh_masks'] = mask
+    # output['thresh_vec'] = eig_vec
+    output['eig_vec'] = eig_vec
+    
+    return output
 
 def rebuild_mean_roi_timecourse(components: np.ndarray,
                                 mask: np.ndarray,
@@ -777,7 +905,6 @@ def rebuild_eigenbrain(eig_vec: np.ndarray,
 
         return eigenbrain
 
-
 def filter_comparison(components: dict,
                       downsample: int = 4,
                       savepath: str = None,

From cbd2936b6167cbd0febcc3bb6b4f2d883226b72c Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Sun, 23 Nov 2025 12:18:42 +1000
Subject: [PATCH 40/57] Added experiment.sort_components() to handle unsorted
 outputs in Balbi data. Added experiment.flip_negative_components() to convert
 older dictionaries to all-positive. Corrected deprecated flip adjustment in
 seas.gui.

---
 seas/experiment.py | 57 +++++++++++++++++++++++++++++++++++++++++++++-
 seas/gui.py        |  8 +++----
 2 files changed, 60 insertions(+), 5 deletions(-)

diff --git a/seas/experiment.py b/seas/experiment.py
index 6cea737..e337e59 100644
--- a/seas/experiment.py
+++ b/seas/experiment.py
@@ -11,6 +11,7 @@
 from seas.ica import project, filter_mean, rebuild_eigenbrain, threshold_by_domains, filter_components, threshold_components, rebuild
 from seas.signalanalysis import sort_noise, lag_n_autocorr
 from seas.waveletAnalysis import waveletAnalysis
+from seas.domains import get_domain_map
 
 from typing import List
 import tifffile as tif
@@ -512,4 +513,58 @@ def export_event_video(components: dict,
                       apply_masked_mean = apply_masked_mean,
                       filter_method = filter_method,
                       include_noise = include_noise)
-    tif.imwrite(outpath, rebuilt.astype(np.float32), imagej=True)
\ No newline at end of file
+    tif.imwrite(outpath, rebuilt.astype(np.float32), imagej=True)
+
+def sort_components(components: dict):
+    eig_vec = components['eig_vec']
+    eig_mix = components['eig_mix']
+    #lag1 = components['lag1']
+    lag1_full = components['lag1_full']
+    noise = components['noise_components']
+
+    ev_sort = np.argsort(eig_mix.std(axis=0)) # Sorting by timecourse standard deviation.
+    #ev_sort = np.argsort(lag1) # Sorting by lag1 auto-correlation
+    eig_vec = eig_vec[:, ev_sort][:, ::-1]
+    eig_mix = eig_mix[:, ev_sort][:, ::-1]
+    #lag1 = lag1[ev_sort][::-1]
+    #noise = noise[ev_sort][::-1]
+    lag1_full = lag1_full[ev_sort][::-1]
+            
+    noise, cutoff = sort_noise(eig_mix.T)
+    components['noise_components'] = noise
+    components['cutoff'] = cutoff
+
+    components['eig_mix'] = eig_mix
+    components['timecourses'] = eig_mix.T
+
+    components['eig_vec'] = eig_vec
+    components['lag1'] = lag_n_autocorr(components['timecourses'], 1)
+    components['lag1_full'] = lag1_full
+
+    # Recalculate domain map
+    domain_map = get_domain_map(components, map_only = False)
+    components.update(domain_map)
+
+    return components
+
+def flip_negative_components(components: dict):
+    n_components = components['n_components']
+    eig_vec = components['eig_vec']
+    eig_mix = components['eig_mix']
+
+    # Track component orientation and ensure positive spatial patterns
+    flipped = np.ones(n_components)
+    for i in range(n_components):
+        # Find the index of maximum absolute value
+        max_idx = np.argmax(np.abs(eig_vec[:, i]))
+        # If that maximum value is negative, flip the component
+        if eig_vec[max_idx, i] < 0:
+            eig_vec[:, i] *= -1
+            eig_mix[:, i] *= -1
+            flipped[i] = -1
+
+    components['flipped'] = flipped
+    components['eig_vec'] = eig_vec
+    components['eig_mix'] = eig_mix
+
+    return components
\ No newline at end of file
diff --git a/seas/gui.py b/seas/gui.py
index ec70389..fb7cad4 100644
--- a/seas/gui.py
+++ b/seas/gui.py
@@ -115,11 +115,11 @@ def run_gui(components: dict,
         print('initializing artifact_components toggle')
         toggle = np.zeros((n_components,), dtype='uint8')
 
-    if 'flipped' in components:
-        flipped = components['flipped']
+    # if 'flipped' in components:
+    #     flipped = components['flipped']
 
-        timecourses = timecourses * flipped[:, None]
-        eig_vec = eig_vec * flipped
+    #     timecourses = timecourses * flipped[:, None]
+    #     eig_vec = eig_vec * flipped
 
     if 'domain_ROIs' in components:
         domain_ROIs = components['domain_ROIs']

From 97388e7356362c66cb20e721fd1fe96d15277173 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Mon, 24 Nov 2025 12:32:09 +1000
Subject: [PATCH 41/57] Replace np.bool with bool for new numpy compatability.

---
 seas/experiment.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/seas/experiment.py b/seas/experiment.py
index e337e59..6dbf098 100644
--- a/seas/experiment.py
+++ b/seas/experiment.py
@@ -475,7 +475,7 @@ def export_event_masks(components: dict,
                                      thresh_type = thresh_type, 
                                      thresh_param = thresh_param)
     components_copy.update(threshold)
-    artifacts_bool = components_copy['artifact_components'].astype(np.bool)
+    artifacts_bool = components_copy['artifact_components'].astype(bool)
     event_components = components_copy['eig_vec'][:, ~artifacts_bool]
     event_masks = rebuild_eigenbrain(event_components,
                                      roimask = components_copy['roimask'], 

From 8fb47034e1eebe847156d0b24f8c8c48858cd492 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Mon, 24 Nov 2025 12:49:03 +1000
Subject: [PATCH 42/57] Added extra thresholding step to
 ica.threshold_by_domains as 'schematic' mode and implementation in
 experiment.export_event_masks. This is to give poster activity maps a more
 contracted/schematic-like appearance.

---
 seas/experiment.py | 6 ++++--
 seas/ica.py        | 6 +++++-
 2 files changed, 9 insertions(+), 3 deletions(-)

diff --git a/seas/experiment.py b/seas/experiment.py
index 6dbf098..c839c51 100644
--- a/seas/experiment.py
+++ b/seas/experiment.py
@@ -468,12 +468,14 @@ def export_event_masks(components: dict,
                        outpath: str,
                        blur: int = 3,
                        thresh_type: str = 'z-score',
-                       thresh_param: float = 7) -> None:
+                       thresh_param: float = 7,
+                       schematic: bool = False) -> None:
     components_copy = components.copy()
     threshold = threshold_by_domains(components_copy, 
                                      blur = blur, 
                                      thresh_type = thresh_type, 
-                                     thresh_param = thresh_param)
+                                     thresh_param = thresh_param,
+                                     schematic = schematic)
     components_copy.update(threshold)
     artifacts_bool = components_copy['artifact_components'].astype(bool)
     event_components = components_copy['eig_vec'][:, ~artifacts_bool]
diff --git a/seas/ica.py b/seas/ica.py
index 4e2b404..72a319b 100644
--- a/seas/ica.py
+++ b/seas/ica.py
@@ -641,7 +641,8 @@ def threshold_by_domains(components: dict,
                    blur: int = 1,
                    min_mask_size: int = 64,
                    thresh_type: str = 'max',
-                   thresh_param: float = None):
+                   thresh_param: float = None,
+                   schematic: bool = False):
     '''
     Function based on modified get_domain_map(). Thresholds ICs using a variety of methods for selective rebuild.
 
@@ -710,6 +711,9 @@ def threshold_by_domains(components: dict,
             z_ROIs_vector = (eig_vec - mean_ROIs_vector)/std_ROIs_vector
             for i in np.arange(eig_vec.shape[0]):
                 mask[i, :] = np.abs(z_ROIs_vector[i]) > thresh_param
+                if schematic:
+                    schem_thresh = np.percentile(np.abs(z_ROIs_vector[i])[mask[i,:]],75) 
+                    mask[i, :] = np.abs(z_ROIs_vector[i]) > schem_thresh
         case 'percentile':
             flipped = components['flipped']
             # Flip ICs where necessary using flipped from dict

From 0debc34e8c8e1c3a70ac33f1da801ba1b0c82c98 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Mon, 24 Nov 2025 12:56:40 +1000
Subject: [PATCH 43/57] Debug schematic thresholding.

---
 seas/ica.py | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/seas/ica.py b/seas/ica.py
index 72a319b..4b694f3 100644
--- a/seas/ica.py
+++ b/seas/ica.py
@@ -710,10 +710,12 @@ def threshold_by_domains(components: dict,
             std_ROIs_vector = np.nanstd(eig_vec, axis=0)
             z_ROIs_vector = (eig_vec - mean_ROIs_vector)/std_ROIs_vector
             for i in np.arange(eig_vec.shape[0]):
-                mask[i, :] = np.abs(z_ROIs_vector[i]) > thresh_param
+                abs_z = np.abs(z_ROIs_vector[i])
+                mask[i, :] = abs_z > thresh_param
                 if schematic:
-                    schem_thresh = np.percentile(np.abs(z_ROIs_vector[i])[mask[i,:]],75) 
-                    mask[i, :] = np.abs(z_ROIs_vector[i]) > schem_thresh
+                    abs_z = abs_z[mask[i, :]]
+                    schem_thresh = np.percentile(abs_z, 75) 
+                    mask[i, :] = abs_z > schem_thresh
         case 'percentile':
             flipped = components['flipped']
             # Flip ICs where necessary using flipped from dict

From 06f66087284e72d6f3849a5ccbf662c763a6a1ca Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Mon, 24 Nov 2025 13:31:59 +1000
Subject: [PATCH 44/57] Further debug of schematic option.

---
 seas/ica.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/seas/ica.py b/seas/ica.py
index 4b694f3..351e5e4 100644
--- a/seas/ica.py
+++ b/seas/ica.py
@@ -697,7 +697,7 @@ def threshold_by_domains(components: dict,
             signal_indices = np.where(artifact_components == 0)[0]
         # eig_vec = eig_vec[:, signal_indices] # Don't change number of ICs, we're updating back to dict
     
-    mask = np.zeros_like(eig_vec, dtype=bool)
+    mask = np.zeros_like(eig_vec, dtype = bool)
 
     match thresh_type:
         case 'max':
@@ -713,8 +713,8 @@ def threshold_by_domains(components: dict,
                 abs_z = np.abs(z_ROIs_vector[i])
                 mask[i, :] = abs_z > thresh_param
                 if schematic:
-                    abs_z = abs_z[mask[i, :]]
-                    schem_thresh = np.percentile(abs_z, 75) 
+                    event = abs_z[mask[i, :]]
+                    schem_thresh = np.percentile(event, 75) 
                     mask[i, :] = abs_z > schem_thresh
         case 'percentile':
             flipped = components['flipped']

From 009e11f90e3fc2b156e729ea493b1a93beeea09c Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Mon, 24 Nov 2025 13:41:31 +1000
Subject: [PATCH 45/57] Further debug on schematic option.

---
 seas/ica.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/seas/ica.py b/seas/ica.py
index 351e5e4..762cf5a 100644
--- a/seas/ica.py
+++ b/seas/ica.py
@@ -712,7 +712,7 @@ def threshold_by_domains(components: dict,
             for i in np.arange(eig_vec.shape[0]):
                 abs_z = np.abs(z_ROIs_vector[i])
                 mask[i, :] = abs_z > thresh_param
-                if schematic:
+                if schematic and abs_z.size != 0:
                     event = abs_z[mask[i, :]]
                     schem_thresh = np.percentile(event, 75) 
                     mask[i, :] = abs_z > schem_thresh

From 3abf1df18531aa3da7fa7725a1d0d986d988db6e Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Mon, 24 Nov 2025 13:44:33 +1000
Subject: [PATCH 46/57] Further debug of schematic option.

---
 seas/ica.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/seas/ica.py b/seas/ica.py
index 762cf5a..30c6d94 100644
--- a/seas/ica.py
+++ b/seas/ica.py
@@ -712,8 +712,8 @@ def threshold_by_domains(components: dict,
             for i in np.arange(eig_vec.shape[0]):
                 abs_z = np.abs(z_ROIs_vector[i])
                 mask[i, :] = abs_z > thresh_param
-                if schematic and abs_z.size != 0:
-                    event = abs_z[mask[i, :]]
+                event = abs_z[mask[i, :]]
+                if schematic and event.size != 0:
                     schem_thresh = np.percentile(event, 75) 
                     mask[i, :] = abs_z > schem_thresh
         case 'percentile':

From 63ec4c1e445754bbddebdba4d76be7d75a4f5e1c Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Mon, 24 Nov 2025 15:34:57 +1000
Subject: [PATCH 47/57] Added schematic fork for ica.threshold_by_domains.

---
 seas/ica.py | 29 +++++++++++++++++++++++++----
 1 file changed, 25 insertions(+), 4 deletions(-)

diff --git a/seas/ica.py b/seas/ica.py
index 30c6d94..1e67ec0 100644
--- a/seas/ica.py
+++ b/seas/ica.py
@@ -14,6 +14,8 @@
 
 import cv2
 from skimage.morphology import remove_small_objects
+from skimage import draw, measure
+from scipy import ndimage
 
 
 def project(vector: np.ndarray,
@@ -712,10 +714,11 @@ def threshold_by_domains(components: dict,
             for i in np.arange(eig_vec.shape[0]):
                 abs_z = np.abs(z_ROIs_vector[i])
                 mask[i, :] = abs_z > thresh_param
-                event = abs_z[mask[i, :]]
-                if schematic and event.size != 0:
-                    schem_thresh = np.percentile(event, 75) 
-                    mask[i, :] = abs_z > schem_thresh
+                # event = abs_z[mask[i, :]]
+                # Deprecated but produced an interesting result
+                # if schematic and event.size != 0:
+                #     schem_thresh = np.percentile(event, 75) 
+                #     mask[i, :] = abs_z > schem_thresh
         case 'percentile':
             flipped = components['flipped']
             # Flip ICs where necessary using flipped from dict
@@ -758,6 +761,24 @@ def threshold_by_domains(components: dict,
                 blurred = cv2.GaussianBlur(eigenbrain, (blur, blur), 0)
                 mask.T[index] = blurred.flat
 
+    if schematic:
+        event_schematics = np.zeros(shape, dtype=np.uint8)
+
+        for i in range(mask.shape[1]):
+            event_schematics.flat[maskind] = mask.T[i]
+            labelled, num_features = ndimage.label(eigenbrain, 
+                                                    structure = [[0,1,0],
+                                                                 [1,1,1],
+                                                                 [0,1,0]])
+            for j in range(num_features):
+                centroids = ndimage.center_of_mass(labelled, 
+                                                   labels = range(num_features))
+                event_size = np.sum(labelled, where = labelled == j)/j
+                schem_radius = np.sqrt(event_size/np.pi)
+                rr, cc = draw.disk(centroids[i], schem_radius, shape = shape)
+                event_schematics[rr, cc] = 255
+                mask.T[i] = event_schematics.flat
+
     mask_bool = mask.astype(bool)
     eig_vec[~mask_bool] = 0
     

From 74b08400f9c64399dad8affc89d8c78fcccc3ad6 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Mon, 24 Nov 2025 15:51:51 +1000
Subject: [PATCH 48/57] Debug schematic fork.

---
 seas/ica.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/seas/ica.py b/seas/ica.py
index 1e67ec0..9e1af59 100644
--- a/seas/ica.py
+++ b/seas/ica.py
@@ -773,9 +773,10 @@ def threshold_by_domains(components: dict,
             for j in range(num_features):
                 centroids = ndimage.center_of_mass(labelled, 
                                                    labels = range(num_features))
+                int_centroids = [tuple(int(x) for x in y) for y in centroids]
                 event_size = np.sum(labelled, where = labelled == j)/j
                 schem_radius = np.sqrt(event_size/np.pi)
-                rr, cc = draw.disk(centroids[i], schem_radius, shape = shape)
+                rr, cc = draw.disk(int_centroids[i], schem_radius, shape = shape)
                 event_schematics[rr, cc] = 255
                 mask.T[i] = event_schematics.flat
 

From 5383a36eb2945655f77fe10b37017cf6666952bb Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Mon, 24 Nov 2025 15:57:55 +1000
Subject: [PATCH 49/57] Further debug on schematic fork.

---
 seas/ica.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/seas/ica.py b/seas/ica.py
index 9e1af59..9795cf8 100644
--- a/seas/ica.py
+++ b/seas/ica.py
@@ -763,9 +763,11 @@ def threshold_by_domains(components: dict,
 
     if schematic:
         event_schematics = np.zeros(shape, dtype=np.uint8)
+        eigenbrain = np.empty(shape)
+        eigenbrain[:] = np.nan
 
         for i in range(mask.shape[1]):
-            event_schematics.flat[maskind] = mask.T[i]
+            eigenbrain.flat[maskind] = mask.T[i]
             labelled, num_features = ndimage.label(eigenbrain, 
                                                     structure = [[0,1,0],
                                                                  [1,1,1],

From ec15638df4eb0666c3a8f32e042e81d3fe6a807f Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Mon, 24 Nov 2025 16:06:06 +1000
Subject: [PATCH 50/57] Further debug schematic fork.

---
 seas/ica.py | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/seas/ica.py b/seas/ica.py
index 9795cf8..c3944ab 100644
--- a/seas/ica.py
+++ b/seas/ica.py
@@ -763,18 +763,21 @@ def threshold_by_domains(components: dict,
 
     if schematic:
         event_schematics = np.zeros(shape, dtype=np.uint8)
+        eigenmask = np.zeros(shape, dtype=bool)
         eigenbrain = np.empty(shape)
         eigenbrain[:] = np.nan
 
         for i in range(mask.shape[1]):
-            eigenbrain.flat[maskind] = mask.T[i]
-            labelled, num_features = ndimage.label(eigenbrain, 
+            eigenmask.flat[maskind] = mask.T[i]
+            eigenbrain.flat[maskind] = eig_vec.T[i]
+            labelled, num_features = ndimage.label(eigenmask, 
                                                     structure = [[0,1,0],
                                                                  [1,1,1],
                                                                  [0,1,0]])
             for j in range(num_features):
-                centroids = ndimage.center_of_mass(labelled, 
-                                                   labels = range(num_features))
+                centroids = ndimage.center_of_mass(eigenbrain, 
+                                                   labels = labelled,
+                                                   index = j)
                 int_centroids = [tuple(int(x) for x in y) for y in centroids]
                 event_size = np.sum(labelled, where = labelled == j)/j
                 schem_radius = np.sqrt(event_size/np.pi)

From eaec6f4997e619f6a5380fd7732a5ff3c26222c3 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Mon, 24 Nov 2025 16:09:22 +1000
Subject: [PATCH 51/57] Further debug schematic fork.

---
 seas/ica.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/seas/ica.py b/seas/ica.py
index c3944ab..941bed9 100644
--- a/seas/ica.py
+++ b/seas/ica.py
@@ -780,7 +780,7 @@ def threshold_by_domains(components: dict,
                                                    index = j)
                 int_centroids = [tuple(int(x) for x in y) for y in centroids]
                 event_size = np.sum(labelled, where = labelled == j)/j
-                schem_radius = np.sqrt(event_size/np.pi)
+                schem_radius = int(np.sqrt(event_size/np.pi))
                 rr, cc = draw.disk(int_centroids[i], schem_radius, shape = shape)
                 event_schematics[rr, cc] = 255
                 mask.T[i] = event_schematics.flat

From a737a68ac72420e27726bf3a7677587661d5d3c1 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Mon, 24 Nov 2025 16:36:49 +1000
Subject: [PATCH 52/57] Further debug schematic fork.

---
 seas/ica.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/seas/ica.py b/seas/ica.py
index 941bed9..c809630 100644
--- a/seas/ica.py
+++ b/seas/ica.py
@@ -775,13 +775,13 @@ def threshold_by_domains(components: dict,
                                                                  [1,1,1],
                                                                  [0,1,0]])
             for j in range(num_features):
-                centroids = ndimage.center_of_mass(eigenbrain, 
+                centroid = ndimage.center_of_mass(eigenbrain, 
                                                    labels = labelled,
                                                    index = j)
-                int_centroids = [tuple(int(x) for x in y) for y in centroids]
+                int_centroid = tuple(int(x) for x in centroid)
                 event_size = np.sum(labelled, where = labelled == j)/j
                 schem_radius = int(np.sqrt(event_size/np.pi))
-                rr, cc = draw.disk(int_centroids[i], schem_radius, shape = shape)
+                rr, cc = draw.disk(int_centroid[i], schem_radius, shape = shape)
                 event_schematics[rr, cc] = 255
                 mask.T[i] = event_schematics.flat
 

From 35c9c69894323509ef0c100603f7e1348fbb582e Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Tue, 25 Nov 2025 21:49:59 +1000
Subject: [PATCH 53/57] Brittle but working schematic fork. Use only blur=1 to
 avoid labelling padded boundary from blur operation.

---
 seas/ica.py | 45 +++++++++++++++++++++++++++------------------
 1 file changed, 27 insertions(+), 18 deletions(-)

diff --git a/seas/ica.py b/seas/ica.py
index c809630..453b380 100644
--- a/seas/ica.py
+++ b/seas/ica.py
@@ -16,6 +16,7 @@
 from skimage.morphology import remove_small_objects
 from skimage import draw, measure
 from scipy import ndimage
+import tifffile as tif
 
 
 def project(vector: np.ndarray,
@@ -762,32 +763,40 @@ def threshold_by_domains(components: dict,
                 mask.T[index] = blurred.flat
 
     if schematic:
-        event_schematics = np.zeros(shape, dtype=np.uint8)
-        eigenmask = np.zeros(shape, dtype=bool)
+        eigenmask = np.zeros(shape, dtype=np.uint8)
         eigenbrain = np.empty(shape)
         eigenbrain[:] = np.nan
 
         for i in range(mask.shape[1]):
+            event_schematic = np.zeros(shape, dtype=np.uint8)
             eigenmask.flat[maskind] = mask.T[i]
             eigenbrain.flat[maskind] = eig_vec.T[i]
-            labelled, num_features = ndimage.label(eigenmask, 
-                                                    structure = [[0,1,0],
-                                                                 [1,1,1],
-                                                                 [0,1,0]])
-            for j in range(num_features):
-                centroid = ndimage.center_of_mass(eigenbrain, 
-                                                   labels = labelled,
-                                                   index = j)
-                int_centroid = tuple(int(x) for x in centroid)
-                event_size = np.sum(labelled, where = labelled == j)/j
-                schem_radius = int(np.sqrt(event_size/np.pi))
-                rr, cc = draw.disk(int_centroid[i], schem_radius, shape = shape)
-                event_schematics[rr, cc] = 255
-                mask.T[i] = event_schematics.flat
-
+            # print("i is:", i)
+            # print(eigenmask)
+            if eigenmask.any():
+                # tif.imwrite("/home/apluff/dev/test_data/eigenmasks/sub-070_eigenmask"+str(i)+".tif", eigenmask, imagej=True)
+                labelled, num_features = ndimage.label(eigenmask)
+                # print(labelled)
+                # print("labelled contains values:", np.unique(labelled))
+                # print("num_features is:", num_features)
+                for j in range(1, num_features + 1):
+                    centroid = ndimage.center_of_mass(eigenmask, 
+                                                      labels = labelled,
+                                                      index = j)
+                    # print("j is:", j)
+                    # print("centroid is:", centroid)
+                    int_centroid = tuple(int(x) for x in centroid)
+                    event_size = np.sum(labelled, where = labelled == j)/j
+                    schem_radius = int(np.sqrt(event_size/np.pi))
+                    rr, cc = draw.disk(int_centroid, schem_radius, shape = shape)
+                    event_schematic[rr, cc] = 255
+                    # print("mask shape is:", mask.shape)
+                    # print("event_schematics shape is:", event_schematic.shape)
+                    mask.T[i] = event_schematic.flat[maskind]
+    
     mask_bool = mask.astype(bool)
     eig_vec[~mask_bool] = 0
-    
+
     output['thresh_masks'] = mask
     # output['thresh_vec'] = eig_vec
     output['eig_vec'] = eig_vec

From f3ce117ba84031ef7920f23d0947cfffe3d32616 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Wed, 10 Dec 2025 09:58:22 +1100
Subject: [PATCH 54/57] Refactor and update experiment.sort_components to use
 lag1 by default.

---
 seas/experiment.py | 38 ++++++++++++++++++++++----------------
 1 file changed, 22 insertions(+), 16 deletions(-)

diff --git a/seas/experiment.py b/seas/experiment.py
index c839c51..72ea08c 100644
--- a/seas/experiment.py
+++ b/seas/experiment.py
@@ -517,32 +517,38 @@ def export_event_video(components: dict,
                       include_noise = include_noise)
     tif.imwrite(outpath, rebuilt.astype(np.float32), imagej=True)
 
-def sort_components(components: dict):
+def sort_components(components: dict, sort_by_noise: bool = True):
     eig_vec = components['eig_vec']
     eig_mix = components['eig_mix']
-    #lag1 = components['lag1']
+    lag1 = components['lag1']
     lag1_full = components['lag1_full']
     noise = components['noise_components']
 
-    ev_sort = np.argsort(eig_mix.std(axis=0)) # Sorting by timecourse standard deviation.
-    #ev_sort = np.argsort(lag1) # Sorting by lag1 auto-correlation
+    if sort_by_noise:
+        ev_sort = np.argsort(lag1) # Sorting by lag1 auto-correlation
+    else:
+        ev_sort = np.argsort(eig_mix.std(axis=0)) # Sorting by timecourse standard deviation.
+    
     eig_vec = eig_vec[:, ev_sort][:, ::-1]
     eig_mix = eig_mix[:, ev_sort][:, ::-1]
-    #lag1 = lag1[ev_sort][::-1]
-    #noise = noise[ev_sort][::-1]
+    lag1 = lag1[ev_sort][::-1]
     lag1_full = lag1_full[ev_sort][::-1]
-            
-    noise, cutoff = sort_noise(eig_mix.T)
-    components['noise_components'] = noise
-    components['cutoff'] = cutoff
-
-    components['eig_mix'] = eig_mix
-    components['timecourses'] = eig_mix.T
-
+    noise = noise[ev_sort][::-1]
+    
+    # Recalculation calls (how PySEAS does it originally)
+    #noise, cutoff = sort_noise(eig_mix.T)
+    #components['cutoff'] = cutoff
+    #components['lag1'] = lag_n_autocorr(components['timecourses'], 1)
+    
+    # Save sorted values
     components['eig_vec'] = eig_vec
-    components['lag1'] = lag_n_autocorr(components['timecourses'], 1)
+    components['eig_mix'] = eig_mix
+    components['lag1'] = lag1
     components['lag1_full'] = lag1_full
-
+    components['noise_components'] = noise
+    # Derived
+    components['timecourses'] = eig_mix.T
+    
     # Recalculate domain map
     domain_map = get_domain_map(components, map_only = False)
     components.update(domain_map)

From 624651e93acdfee1fde1b7a80918f977cf1c9ae7 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Thu, 11 Dec 2025 11:07:07 +1000
Subject: [PATCH 55/57] Refine sort_components for Balbi data. Sort
 artifact_components and don't run domain_map.

---
 seas/experiment.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/seas/experiment.py b/seas/experiment.py
index 72ea08c..a70760e 100644
--- a/seas/experiment.py
+++ b/seas/experiment.py
@@ -523,6 +523,7 @@ def sort_components(components: dict, sort_by_noise: bool = True):
     lag1 = components['lag1']
     lag1_full = components['lag1_full']
     noise = components['noise_components']
+    artifacts = components['artifact_components']
 
     if sort_by_noise:
         ev_sort = np.argsort(lag1) # Sorting by lag1 auto-correlation
@@ -534,6 +535,7 @@ def sort_components(components: dict, sort_by_noise: bool = True):
     lag1 = lag1[ev_sort][::-1]
     lag1_full = lag1_full[ev_sort][::-1]
     noise = noise[ev_sort][::-1]
+    artifacts = artifacts[ev_sort][::-1]
     
     # Recalculation calls (how PySEAS does it originally)
     #noise, cutoff = sort_noise(eig_mix.T)
@@ -546,12 +548,13 @@ def sort_components(components: dict, sort_by_noise: bool = True):
     components['lag1'] = lag1
     components['lag1_full'] = lag1_full
     components['noise_components'] = noise
+    components['artifact_components'] = artifacts
     # Derived
     components['timecourses'] = eig_mix.T
     
     # Recalculate domain map
-    domain_map = get_domain_map(components, map_only = False)
-    components.update(domain_map)
+    # domain_map = get_domain_map(components, map_only = False)
+    # components.update(domain_map)
 
     return components
 

From c738e41c8932fee887dfa179bce4fc557c2d5f15 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Thu, 11 Dec 2025 11:16:32 +1000
Subject: [PATCH 56/57] Make sorting of artifact_components conditional to
 generalise to newly decomposed pyseas dicts.

---
 seas/experiment.py | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/seas/experiment.py b/seas/experiment.py
index a70760e..34dd5ad 100644
--- a/seas/experiment.py
+++ b/seas/experiment.py
@@ -523,7 +523,6 @@ def sort_components(components: dict, sort_by_noise: bool = True):
     lag1 = components['lag1']
     lag1_full = components['lag1_full']
     noise = components['noise_components']
-    artifacts = components['artifact_components']
 
     if sort_by_noise:
         ev_sort = np.argsort(lag1) # Sorting by lag1 auto-correlation
@@ -535,7 +534,11 @@ def sort_components(components: dict, sort_by_noise: bool = True):
     lag1 = lag1[ev_sort][::-1]
     lag1_full = lag1_full[ev_sort][::-1]
     noise = noise[ev_sort][::-1]
-    artifacts = artifacts[ev_sort][::-1]
+    
+    if 'artifact_components' in components:
+        artifacts = components['artifact_components']
+        artifacts = artifacts[ev_sort][::-1]
+        components['artifact_components'] = artifacts
     
     # Recalculation calls (how PySEAS does it originally)
     #noise, cutoff = sort_noise(eig_mix.T)
@@ -548,7 +551,7 @@ def sort_components(components: dict, sort_by_noise: bool = True):
     components['lag1'] = lag1
     components['lag1_full'] = lag1_full
     components['noise_components'] = noise
-    components['artifact_components'] = artifacts
+    
     # Derived
     components['timecourses'] = eig_mix.T
     

From 027b49fb97f7ba60e1a02d5711d51920590c6d43 Mon Sep 17 00:00:00 2001
From: Adam Luff <aluff2@gmail.com>
Date: Thu, 11 Dec 2025 13:39:58 +1000
Subject: [PATCH 57/57] Tidying sort_components.

---
 seas/experiment.py | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/seas/experiment.py b/seas/experiment.py
index 34dd5ad..20bfcc4 100644
--- a/seas/experiment.py
+++ b/seas/experiment.py
@@ -540,11 +540,6 @@ def sort_components(components: dict, sort_by_noise: bool = True):
         artifacts = artifacts[ev_sort][::-1]
         components['artifact_components'] = artifacts
     
-    # Recalculation calls (how PySEAS does it originally)
-    #noise, cutoff = sort_noise(eig_mix.T)
-    #components['cutoff'] = cutoff
-    #components['lag1'] = lag_n_autocorr(components['timecourses'], 1)
-    
     # Save sorted values
     components['eig_vec'] = eig_vec
     components['eig_mix'] = eig_mix
@@ -552,10 +547,15 @@ def sort_components(components: dict, sort_by_noise: bool = True):
     components['lag1_full'] = lag1_full
     components['noise_components'] = noise
     
-    # Derived
+    # Derive from sorted values
     components['timecourses'] = eig_mix.T
     
-    # Recalculate domain map
+    # Recalculation calls (how PySEAS does it originally)
+    #noise, cutoff = sort_noise(eig_mix.T)
+    #components['cutoff'] = cutoff
+    #components['lag1'] = lag_n_autocorr(components['timecourses'], 1)
+
+    # Recalculate domain map (doesn't work for some reason)
     # domain_map = get_domain_map(components, map_only = False)
     # components.update(domain_map)