ENH: adding gabor filters to ndimage

benchmarks/benchmarks/ndimage.py

@@ -0,0 +1,28 @@
+from __future__ import division, absolute_import, print_function
+
+import numpy as np
+import timeit
+
+try:
+    from scipy.ndimage import gabor_filter
+except ImportError:
+    pass
+
+from .common import Benchmark
+
+
+class GaborFilter(Benchmark):
+    param_names = ['truncate']
+    params = [1, 2, 3]
+
+    def setup(self, truncate):
+        np.random.seed(123456)
+        self.image = np.random.randint(0, 256, size=[25, 25])
+        self.volume = np.random.randint(0, 256, size=[25, 25, 25])
+
+    def time_gabor_2d(self, truncate):
+        gabor_filter(self.image, 2.5, 0, 0.1, 0, truncate=truncate)
+
+    def time_gabor_3d(self, truncate):
+        gabor_filter(self.volume, 5, 0, 0.1, 0, truncate=truncate)
+

scipy/ndimage/__init__.py

@@ -19,6 +19,7 @@
    convolve1d - 1-D convolution along the given axis
    correlate - Multidimensional correlation
    correlate1d - 1-D correlation along the given axis
+   gabor_filter
    gaussian_filter
    gaussian_filter1d
    gaussian_gradient_magnitude

scipy/ndimage/_ni_support.py

@@ -79,3 +79,16 @@ def _get_output(output, input, shape=None):
     elif output.shape != shape:
         raise RuntimeError("output shape not correct")
     return output
+
+def _check_axis(axis, rank):
+    if axis < 0:
+        axis += rank
+    if axis < 0 or axis >= rank:
+        raise ValueError('invalid axis')
+    return axis
+
+
+def _check_float(num):
+    return float(num)
+
+

scipy/ndimage/filters.py

@@ -44,7 +44,7 @@
            'uniform_filter1d', 'uniform_filter', 'minimum_filter1d',
            'maximum_filter1d', 'minimum_filter', 'maximum_filter',
            'rank_filter', 'median_filter', 'percentile_filter',
-           'generic_filter1d', 'generic_filter']
+           'generic_filter1d', 'generic_filter', 'gabor_filter']
 
 
 def _invalid_origin(origin, lenw):
@@ -303,6 +303,140 @@ def gaussian_filter(input, sigma, order=0, output=None,
     return output
 
 
+@_ni_docstrings.docfiller
+def gabor_filter(input, sigma, phi, frequency, offset=0.0, output=None,
+                 mode="reflect", cval=0.0, truncate=3.0):
+    """Multidimensional Gabor filter. A gabor filter
+    is an elementwise product between a Gaussian
+    and a complex exponential.
+
+    Parameters
+    ----------
+    %(input)s
+    sigma : scalar or sequence of scalars
+        Standard deviation for Gaussian kernel. The standard
+        deviations of the Gaussian filter are given for each axis as a
+        sequence, or as a single number, in which case it is equal for
+        all axes.
+    phi : scalar or sequence of scalars
+        Angles specifying orientation of the periodic complex
+        exponential. If the input is n-dimensional, then phi
+        is a sequence of length n-1. Convention follows
+        https://en.wikipedia.org/wiki/N-sphere#Spherical_coordinates.
+    frequency : scalar
+        Frequency of the complex exponential. Units are revolutions/voxels.
+    offset : scalar
+        Phase shift of the complex exponential. Units are radians.
+    output : array or dtype, optional
+        The array in which to place the output, or the dtype of the
+        returned array. By default an array of the same dtype as input
+        will be created. Only the real component will be saved
+        if output is an array.
+    %(mode)s
+    %(cval)s
+    truncate : float
+        Truncate the filter at this many standard deviations.
+        Default is 3.0.
+
+    Returns
+    -------
+    real, imaginary : arrays
+        Returns real and imaginary responses, arrays of same
+        shape as `input`.
+
+    Notes
+    -----
+    The multidimensional filter is implemented by creating
+    a gabor filter array, then using the convolve method.
+    Also, sigma specifies the standard deviations of the
+    Gaussian along the coordinate axes, and the Gaussian
+    is not rotated. This is unlike
+    skimage.filters.gabor, whose Gaussian is
+    rotated with the complex exponential.
+    The reasoning behind this design choice is that
+    sigma can be more easily designed to deal with
+    anisotropic voxels.
+
+    Examples
+    --------
+    >>> from scipy.ndimage import gabor_filter
+    >>> a = np.arange(50, step=2).reshape((5,5))
+    >>> a
+    array([[ 0,  2,  4,  6,  8],
+           [10, 12, 14, 16, 18],
+           [20, 22, 24, 26, 28],
+           [30, 32, 34, 36, 38],
+           [40, 42, 44, 46, 48]])
+    >>> gabor_filter(a, sigma=1, phi=[0.0], frequency=0.1)
+    (array([[ 3,  5,  6,  8,  9],
+            [ 9, 10, 12, 13, 14],
+            [16, 18, 19, 21, 22],
+            [24, 25, 27, 28, 30],
+            [29, 31, 32, 34, 35]]),
+    array([[ 0,  0, -1,  0,  0],
+            [ 0,  0, -1,  0,  0],
+            [ 0,  0, -1,  0,  0],
+            [ 0,  0, -1,  0,  0],
+            [ 0,  0, -1,  0,  0]]))
+    >>> from scipy import misc
+    >>> import matplotlib.pyplot as plt
+    >>> fig = plt.figure()
+    >>> plt.gray()  # show the filtered result in grayscale
+    >>> ax1 = fig.add_subplot(121)  # left side
+    >>> ax2 = fig.add_subplot(122)  # right side
+    >>> ascent = misc.ascent()
+    >>> result = gabor_filter(ascent, sigma=5, phi=[0.0], frequency=0.1)
+    >>> ax1.imshow(ascent)
+    >>> ax2.imshow(result[0])
+    >>> plt.show()
+    """
+    input = numpy.asarray(input)
+    sigmas = _ni_support._normalize_sequence(sigma, input.ndim)
+    phi = _ni_support._normalize_sequence(phi, input.ndim - 1)
+    frequency = _ni_support._check_float(frequency)
+    offset = _ni_support._check_float(offset)
+
+    limits = [numpy.ceil(truncate * sigma).astype(int) for sigma in sigmas]
+    ranges = [range(-limit, limit + 1) for limit in limits]
+    coords = numpy.meshgrid(*ranges, indexing='ij')
+    filter_size = coords[0].shape
+    coords = numpy.stack(coords, axis=-1)
+
+    new_shape = numpy.ones(input.ndim)
+    new_shape = numpy.append(new_shape, -1).astype(int)
+    sigmas = numpy.reshape(sigmas, new_shape)
+
+    g = numpy.zeros(filter_size, dtype=numpy.complex)
+    g[:] = numpy.exp(-0.5 * numpy.sum(numpy.divide(coords, sigmas)**2, axis=-1))
+
+    g /= (2 * numpy.pi)**(input.ndim / 2) * numpy.prod(sigmas)
+
+    orientation = numpy.ones(input.ndim)
+    for i, p in enumerate(phi):
+        orientation[i + 1] = orientation[i] * numpy.sin(p)
+        orientation[i] = orientation[i] * numpy.cos(p)
+    orientation = numpy.flip(orientation, axis=0)
+
+    rotx = coords @ orientation
+
+    g *= numpy.exp(1j * (2 * numpy.pi * frequency * rotx + offset))
+
+    if isinstance(output, (type, numpy.dtype)):
+        otype = output
+    elif isinstance(output, str):
+        otype = numpy.typeDict[output]
+    else:
+        otype = None
+
+    output = convolve(input, weights=numpy.real(g),
+                      output=output, mode=mode, cval=cval)
+    imag = convolve(input, weights=numpy.imag(g),
+                    output=otype, mode=mode, cval=cval)
+
+    result = (output, imag)
+    return result
+
+
 @_ni_docstrings.docfiller
 def prewitt(input, axis=-1, output=None, mode="reflect", cval=0.0):
     """Calculate a Prewitt filter.

scipy/ndimage/tests/test_filters.py

@@ -315,6 +315,21 @@ def test_gaussian_truncate():
     n = nonzero_indices.ptp() + 1
     assert_equal(n, 15)
 
+    # Test gabor_filter
+    num_nonzeros_2 = (np.abs(sndi.gabor_filter(arr, 5, 0, 5, truncate=2))[0] > 0).sum()
+    assert_equal(num_nonzeros_2, 21**2)
+    num_nonzeros_5 = (np.abs(sndi.gabor_filter(arr, 5, 0, 5, truncate=5))[0] > 0).sum()
+    assert_equal(num_nonzeros_5, 51**2)
+
+    f = np.abs(sndi.gabor_filter(arr, [0.5, 2.5], 0, 5, truncate=3.5)[0])
+    fpos = f > 0
+    n0 = fpos.any(axis=0).sum()
+    # n0 should be 2*int(2.5*3.5 + 0.5) + 1
+    assert_equal(n0, 19)
+    n1 = fpos.any(axis=1).sum()
+    # n1 should be 2*int(0.5*3.5 + 0.5) + 1
+    assert_equal(n1, 5)
+
 
 class TestThreading(object):
     def check_func_thread(self, n, fun, args, out):

scipy/ndimage/tests/test_ndimage.py

@@ -692,6 +692,54 @@ def derivative(input, axis, output, mode, cval, a, b):
             extra_keywords={'b': 2.0})
         assert_array_almost_equal(tmp1, tmp2)
 
+    def test_gabor_filter01(self):
+        input = numpy.array([[1, 2, 3],
+                        [2, 4, 6]], numpy.float32)
+        real, imag = ndimage.gabor_filter(input, 1.0, 0, 5)
+        assert_equal(input.dtype, real.dtype)
+        assert_equal(input.shape, real.shape)
+        assert_equal(input.dtype, imag.dtype)
+        assert_equal(input.shape, imag.shape)
+
+    def test_gabor_filter02(self):
+        input = numpy.arange(100 * 100).astype(numpy.float32)
+        input.shape = (100, 100)
+        otype = numpy.float64
+        real, imag = ndimage.gabor_filter(input, [1.0, 1.0], 0, 5, output=otype)
+        assert_equal(real.dtype.type, numpy.float64)
+        assert_equal(input.shape, real.shape)
+        assert_equal(imag.dtype.type, numpy.float64)
+        assert_equal(input.shape, imag.shape)
+
+    def test_gabor_filter03(self):
+        # Tests that inputting an output array works
+        input = numpy.arange(100 * 100).astype(numpy.float32)
+        input.shape = (100, 100)
+        output = numpy.zeros([100, 100]).astype(numpy.float32)
+        ndimage.gabor_filter(input, [1.0, 1.0], 0, 5, output=output)
+        assert_equal(input.dtype, output.dtype)
+        assert_equal(input.shape, output.shape)
+
+    def test_gabor_filter04(self):
+        # Tests that list of sigmas is same as single
+        input = numpy.arange(100 * 100).astype(numpy.float32)
+        input.shape = (100, 100)
+        otype = numpy.float64
+        real1, imag1 = ndimage.gabor_filter(input, [1.0, 1.0], 0, 5, output=otype)
+        real2, imag2 = ndimage.gabor_filter(input, 1.0, 0, 5, output=otype)
+        assert_array_almost_equal(real1, real2)
+        assert_array_almost_equal(imag1, imag2)
+
+    def test_gabor_filter05(self):
+        # Tests that list of phi works
+        input = numpy.arange(100 * 100 * 100).astype(numpy.float32)
+        input.shape = (100, 100, 100)
+        otype = numpy.float64
+        real1, imag1 = ndimage.gabor_filter(input, 1.0, [0.5, 0.5], 5, output=otype)
+        real2, imag2 = ndimage.gabor_filter(input, 1.0, 0.5, 5, output=otype)
+        assert_array_almost_equal(real1, real2)
+        assert_array_almost_equal(imag1, imag2)
+
     def test_uniform01(self):
         array = numpy.array([2, 4, 6])
         size = 2