Moved band-pass filters to dedicated module in misc

Signed-off-by: Nicola VIGANÒ <nicola.vigano@esrf.fr>

doc/apidoc/nabu.misc.fourier_filters.rst

+nabu.misc.fourier_filters module
+========================
+
+.. automodule:: nabu.misc.fourier_filters
+   :members:
+   :undoc-members:
+   :show-inheritance:

nabu/misc/fourier_filters.py

+# -*- coding: utf-8 -*-
+"""
+Fourier filters.
+"""
+
+import numpy as np
+
+try:
+    import scipy.special as spspe
+    __have_scipy__ = True
+except ImportError:
+    import math
+    __have_scipy__ = False
+
+
+def get_lowpass_filter(img_shape, cutoff_par):
+    """Computes a low pass filter with the erfc.
+
+    Parameters
+    ----------
+    img_shape: tuple
+        Shape of the image
+    cutoff_par: float or sequence of two floats
+        Position of the cut off in pixels, if a sequence is given the second float expresses the
+        width of the transition region which is given as a fraction of the cutoff frequency.
+        When only one float is given for this argument a gaussian is applied whose sigma is the
+        parameter.
+        When a sequence of two numbers is given then the filter is 1 ( no filtering) till the cutoff
+        frequency while a smooth erfc transition to zero is done
+
+    Raises
+    ------
+    ValueError
+        In case of malformed cutoff_par
+
+    Returns
+    -------
+    numpy.array_like
+        The computed filter
+    """
+    if isinstance(cutoff_par, (int, float)):
+        cutoff_pix = cutoff_par
+        cutoff_trans_fact = None
+    else:
+        try:
+            cutoff_pix, cutoff_trans_fact = cutoff_par
+        except ValueError:
+            raise ValueError(
+                "Argument cutoff_par  (which specifies the pass filter shape) must be either a scalar or a"
+                " sequence of two scalars"
+            )
+        if (not isinstance(cutoff_pix, (int, float))) or (not isinstance(cutoff_trans_fact, (int, float))):
+            raise ValueError(
+                "Argument cutoff_par  (which specifies the pass filter shape) must be  one number or a sequence"
+                "of two numbers"
+            )
+
+    coords = [np.fft.fftfreq(s, 1) for s in img_shape]
+    coords = np.meshgrid(*coords, indexing="ij")
+
+    r = np.sqrt(np.sum(np.array(coords) ** 2, axis=0))
+
+    if cutoff_trans_fact is not None:
+        k_cut = 0.5 / cutoff_pix
+        k_cut_width = k_cut * cutoff_trans_fact
+        k_pos_rescaled = (r - k_cut) / k_cut_width
+        if __have_scipy__:
+            res = spspe.erfc(k_pos_rescaled) / 2
+        else:
+            res = np.array(list(map(math.erfc, k_pos_rescaled))) / 2
+    else:
+        res = np.exp(- (np.pi ** 2) * (r ** 2) * (cutoff_pix ** 2) * 2)
+
+    return res
+
+
+def get_highpass_filter(img_shape, cutoff_par):
+    """Computes a high pass filter with the erfc.
+
+    Parameters
+    ----------
+    img_shape: tuple
+        Shape of the image
+    cutoff_par: float or sequence of two floats
+        Position of the cut off in pixels, if a sequence is given the second float expresses the
+        width of the transition region which is given as a fraction of the cutoff frequency.
+        When only one float is given for this argument a gaussian is applied whose sigma is the
+        parameter, and the result is subtracted from 1 to obtain the high pass filter
+        When a sequence of two numbers is given then the filter is 1 ( no filtering) above the cutoff
+        frequency and then a smooth  transition to zero is done for smaller frequency
+
+    Raises
+    ------
+    ValueError
+        In case of malformed cutoff_par
+
+    Returns
+    -------
+    numpy.array_like
+        The computed filter
+    """
+    return 1 - get_lowpass_filter(img_shape, cutoff_par)

nabu/preproc/alignment.py

 import numpy as np
-import math
+
 import logging
 from numpy.polynomial.polynomial import Polynomial
 
 from nabu.utils import previouspow2
+from nabu.misc import fourier_filters
 
 try:
     from scipy.ndimage.filters import median_filter
-    import scipy.special as spspe
 
     __have_scipy__ = True
 except ImportError:
@@ -16,87 +16,6 @@ except ImportError:
     __have_scipy__ = False
 
 
-def _get_lowpass_filter(img_shape, cutoff_par):
-    """Computes a low pass filter with the erfc.
-
-    Parameters
-    ----------
-    img_shape: tuple
-        Shape of the image
-    cutoff_par: float or sequence of two floats
-        Position of the cut off in pixels, if a sequence is given the second float expresses the
-        width of the transition region which is given as a fraction of the cutoff frequency.
-        When only one float is given for this argument a gaussian is applied whose sigma is the
-        parameter.
-        When a sequence of two numbers is given then the filter is 1 ( no filtering) till the cutoff
-        frequency while a smooth erfc transition to zero is done
-
-    Raises
-    ------
-    ValueError
-            In case cutoff_par is not properly given
-
-    Returns
-    -------
-    numpy.array_like
-        The computed filter
-    """
-    if isinstance(cutoff_par, (int, float)):
-        cutoff_pix = cutoff_par
-        cutoff_trans_fact = None
-    else:
-        try:
-            cutoff_pix, cutoff_trans_fact = cutoff_par
-        except ValueError:
-            raise ValueError("Argument cutoff_par  (which specifies the pass filter shape) must be either a scalar or a"
-                             " sequence of two scalars")
-        if (not isinstance(cutoff_pix, (int, float))) or (not isinstance(cutoff_trans_fact, (int, float))):
-            raise ValueError("Argument cutoff_par  (which specifies the pass filter shape) must be  one number or a sequence"
-                             "of two numbers")
-
-    coords = [np.fft.fftfreq(s, 1) for s in img_shape]
-    coords = np.meshgrid(*coords, indexing="ij")
-
-    r = np.sqrt(np.sum(np.array(coords) ** 2, axis=0))
-
-    if cutoff_trans_fact is not None:
-        k_cut = 0.5 / cutoff_pix
-        k_cut_width = k_cut * cutoff_trans_fact
-        k_pos_rescaled = (r - k_cut) / k_cut_width
-        if __have_scipy__:
-            res = spspe.erfc(k_pos_rescaled) / 2
-        else:
-            res = np.array(list(map(math.erfc, k_pos_rescaled))) / 2
-    else:
-        res = np.exp(- np.pi*np.pi*r*r*cutoff_pix*cutoff_pix*2)
-
-    return res
-
-
-def _get_highpass_filter(img_shape, cutoff_pix, cutoff_par):
-    """Computes a high pass filter with the erfc.
-
-    Parameters
-    ----------
-    img_shape: tuple
-        Shape of the image
-    cutoff_par: float or sequence of two floats
-        Position of the cut off in pixels, if a sequence is given the second float expresses the
-        width of the transition region which is given as a fraction of the cutoff frequency.
-        When only one float is given for this argument a gaussian is applied whose sigma is the
-        parameter, and the result is subtracted from 1 to obtain the high pass filter
-        When a sequence of two numbers is given then the filter is 1 ( no filtering) above the cutoff
-        frequency and then a smooth  transition to zero is done for smaller frequency
-
-    Returns
-    -------
-    numpy.array_like
-        The computed filter
-    """
-    res = 1 - _get_lowpass_filter(img_shape, cutoff_pix, cutoff_par)
-    return res
-
-
 class AlignmentBase(object):
     @staticmethod
     def refine_max_position_2d(f_vals: np.ndarray, fy=None, fx=None):
@@ -332,20 +251,13 @@ class AlignmentBase(object):
         img[np.isinf(img)] = invalid_val
 
         if high_pass is not None or low_pass is not None:
-            myfilter = np.ones_like(img)
+            img_filter = np.ones(img.shape[-2:], dtype=img.dtype)
             if low_pass is not None:
-                myfilter[:] *= _get_lowpass_filter(img.shape[-2:], low_pass)
+                img_filter[:] *= fourier_filters.get_lowpass_filter(img.shape[-2:], low_pass)
             if high_pass is not None:
-                myfilter[:] *= _get_highpass_filter(img.shape[-2:], high_pass)
-            img_shape = img.shape
-            if len(img_shape) == 2:
-                img = np.fft.ifft2(np.fft.fft2(img) * myfilter).real
-            elif len(img_shape) > 2:
-                # if dealing with a stack of images, we have to do them one by one
-                img = np.reshape(img, (-1,) + tuple(img_shape[-2:]))
-                for ii in range(img.shape[0]):
-                    img[ii, ...] = np.fft.ifft2(np.fft.fft2(img[ii, ...]) * myfilter).real
-                img = np.reshape(img, img_shape)
+                img_filter[:] *= fourier_filters.get_highpass_filter(img.shape[-2:], high_pass)
+            # fft2 and iff2 use axes=(-2, -1) by default
+            img = np.fft.ifft2(np.fft.fft2(img) * img_filter).real
 
         if median_filt_shape is not None:
             img_shape = img.shape