From b5731b6c77d24d742ac6f85718b3b751d312fb6e Mon Sep 17 00:00:00 2001
From: Oscar Esteban <code@oscaresteban.es>
Date: Thu, 15 Dec 2022 17:34:56 +0100
Subject: [PATCH] enh: create new filtering submodule

---
 src/eddymotion/data/filtering.py | 64 ++++++++++++++++++++++++++++++++
 src/eddymotion/estimator.py      | 47 ++---------------------
 2 files changed, 68 insertions(+), 43 deletions(-)
 create mode 100644 src/eddymotion/data/filtering.py

diff --git a/src/eddymotion/data/filtering.py b/src/eddymotion/data/filtering.py
new file mode 100644
index 00000000..c84cc02b
--- /dev/null
+++ b/src/eddymotion/data/filtering.py
@@ -0,0 +1,64 @@
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+#
+# Copyright 2022 The NiPreps Developers <nipreps@gmail.com>
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+# We support and encourage derived works from this project, please read
+# about our expectations at
+#
+#     https://www.nipreps.org/community/licensing/
+#
+"""Filtering data."""
+
+
+def advanced_clip(
+    data, p_min=35, p_max=99.98, nonnegative=True, dtype="int16", invert=False
+):
+    """
+    Remove outliers at both ends of the intensity distribution and fit into a given dtype.
+
+    This interface tries to emulate ANTs workflows' massaging that truncate images into
+    the 0-255 range, and applies percentiles for clipping images.
+    For image registration, normalizing the intensity into a compact range (e.g., uint8)
+    is generally advised.
+
+    To more robustly determine the clipping thresholds, spikes are removed from data with
+    a median filter.
+    Once the thresholds are calculated, the denoised data are thrown away and the thresholds
+    are applied on the original image.
+
+    """
+    import numpy as np
+    from scipy import ndimage
+    from skimage.morphology import ball
+
+    # Calculate stats on denoised version, to preempt outliers from biasing
+    denoised = ndimage.median_filter(data, footprint=ball(3))
+
+    a_min = np.percentile(denoised[denoised > 0] if nonnegative else denoised, p_min)
+    a_max = np.percentile(denoised[denoised > 0] if nonnegative else denoised, p_max)
+
+    # Clip and cast
+    data = np.clip(data, a_min=a_min, a_max=a_max)
+    data -= data.min()
+    data /= data.max()
+
+    if invert:
+        data = 1.0 - data
+
+    if dtype in ("uint8", "int16"):
+        data = np.round(255 * data).astype(dtype)
+
+    return data
diff --git a/src/eddymotion/estimator.py b/src/eddymotion/estimator.py
index 6c7a4970..e4fc3d9d 100644
--- a/src/eddymotion/estimator.py
+++ b/src/eddymotion/estimator.py
@@ -91,7 +91,7 @@ def fit(
             ).to_filename(bmask_img)
             kwargs["mask"] = dwdata.brainmask
 
-        kwargs["S0"] = _advanced_clip(dwdata.bzero)
+        kwargs["S0"] = dwdata.bzero
 
         if "num_threads" not in align_kwargs and omp_nthreads is not None:
             align_kwargs["num_threads"] = omp_nthreads
@@ -208,51 +208,12 @@ def fit(
         return dwdata.em_affines
 
 
-def _advanced_clip(
-    data, p_min=35, p_max=99.98, nonnegative=True, dtype="int16", invert=False
-):
-    """
-    Remove outliers at both ends of the intensity distribution and fit into a given dtype.
-
-    This interface tries to emulate ANTs workflows' massaging that truncate images into
-    the 0-255 range, and applies percentiles for clipping images.
-    For image registration, normalizing the intensity into a compact range (e.g., uint8)
-    is generally advised.
-
-    To more robustly determine the clipping thresholds, spikes are removed from data with
-    a median filter.
-    Once the thresholds are calculated, the denoised data are thrown away and the thresholds
-    are applied on the original image.
-
-    """
-    import numpy as np
-    from scipy import ndimage
-    from skimage.morphology import ball
-
-    # Calculate stats on denoised version, to preempt outliers from biasing
-    denoised = ndimage.median_filter(data, footprint=ball(3))
-
-    a_min = np.percentile(denoised[denoised > 0] if nonnegative else denoised, p_min)
-    a_max = np.percentile(denoised[denoised > 0] if nonnegative else denoised, p_max)
-
-    # Clip and cast
-    data = np.clip(data, a_min=a_min, a_max=a_max)
-    data -= data.min()
-    data /= data.max()
-
-    if invert:
-        data = 1.0 - data
-
-    if dtype in ("uint8", "int16"):
-        data = np.round(255 * data).astype(dtype)
-
-    return data
-
-
 def _to_nifti(data, affine, filename, clip=True):
     data = np.squeeze(data)
     if clip:
-        data = _advanced_clip(data)
+        from eddymotion.data.filtering import advanced_clip
+
+        data = advanced_clip(data)
     nii = nb.Nifti1Image(
         data,
         affine,