Merge pull request #712 from emmanuelle/gaussian_filter

Add a wrapper around `scipy.ndimage.gaussian_filter` with useful default behaviors.

skimage/filter/__init__.py

@@ -1,5 +1,6 @@
 from .lpi_filter import inverse, wiener, LPIFilter2D
 from .ctmf import median_filter
+from ._gaussian import gaussian_filter
 from ._canny import canny
 from .edges import (sobel, hsobel, vsobel, scharr, hscharr, vscharr, prewitt,
                     hprewitt, vprewitt, roberts , roberts_positive_diagonal,
@@ -16,6 +17,7 @@ __all__ = ['inverse',
            'wiener',
            'LPIFilter2D',
            'median_filter',
+           'gaussian_filter',
            'canny',
            'sobel',
            'hsobel',

skimage/filter/_gaussian.py

@@ -0,0 +1,105 @@
+import collections as coll
+import numpy as np
+from scipy import ndimage
+import warnings
+
+from ..util import img_as_float
+from ..color import guess_spatial_dimensions
+
+__all__ = ['gaussian_filter']
+
+
+def gaussian_filter(image, sigma, output=None, mode='nearest', cval=0,
+                    multichannel=None):
+    """
+    Multi-dimensional Gaussian filter
+
+    Parameters
+    ----------
+
+    image : array-like
+        input image (grayscale or color) to filter.
+    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.
+    output : array, optional
+        The ``output`` parameter passes an array in which to store the
+        filter output.
+    mode : {'reflect', 'constant', 'nearest', 'mirror', 'wrap'}, optional
+        The `mode` parameter determines how the array borders are
+        handled, where `cval` is the value when mode is equal to
+        'constant'. Default is 'nearest'.
+    cval : scalar, optional
+        Value to fill past edges of input if `mode` is 'constant'. Default
+        is 0.0
+    multichannel : bool, optional (default: None)
+        Whether the last axis of the image is to be interpreted as multiple
+        channels. If True, each channel is filtered separately (channels are
+        not mixed together). Only 3 channels are supported. If `None`,
+        the function will attempt to guess this, and raise a warning if
+        ambiguous, when the array has shape (M, N, 3).
+
+    Returns
+    -------
+
+    filtered_image : ndarray
+        the filtered array
+
+    Notes
+    -----
+
+    This function is a wrapper around :func:`scipy.ndimage.gaussian_filter`.
+
+    Integer arrays are converted to float.
+
+    The multi-dimensional filter is implemented as a sequence of
+    one-dimensional convolution filters. The intermediate arrays are
+    stored in the same data type as the output. Therefore, for output
+    types with a limited precision, the results may be imprecise
+    because intermediate results may be stored with insufficient
+    precision.
+
+    Examples
+    --------
+
+    >>> a = np.zeros((3, 3))
+    >>> a[1, 1] = 1
+    >>> a
+    array([[ 0.,  0.,  0.],
+           [ 0.,  1.,  0.],
+           [ 0.,  0.,  0.]])
+    >>> gaussian_filter(a, sigma=0.4)  # mild smoothing
+    array([[ 0.00163116,  0.03712502,  0.00163116],
+           [ 0.03712502,  0.84496158,  0.03712502],
+           [ 0.00163116,  0.03712502,  0.00163116]])
+    >>> gaussian_filter(a, sigma=1)  # more smooting
+    array([[ 0.05855018,  0.09653293,  0.05855018],
+           [ 0.09653293,  0.15915589,  0.09653293],
+           [ 0.05855018,  0.09653293,  0.05855018]])
+    >>> # Several modes are possible for handling boundaries
+    >>> gaussian_filter(a, sigma=1, mode='reflect')
+    array([[ 0.08767308,  0.12075024,  0.08767308],
+           [ 0.12075024,  0.16630671,  0.12075024],
+           [ 0.08767308,  0.12075024,  0.08767308]])
+    >>> # For RGB images, each is filtered separately
+    >>> from skimage.data import lena
+    >>> image = lena()
+    >>> filtered_lena = gaussian_filter(image, sigma=1, multichannel=True)
+    """
+    spatial_dims = guess_spatial_dimensions(image)
+    if spatial_dims is None and multichannel is None:
+        msg = ("Images with dimensions (M, N, 3) are interpreted as 2D+RGB" +
+                   " by default. Use `multichannel=False` to interpret as " +
+                   " 3D image with last dimension of length 3.")
+        warnings.warn(RuntimeWarning(msg))
+        multichannel = True
+    if multichannel:
+        # do not filter across channels
+        if not isinstance(sigma, coll.Iterable):
+            sigma = [sigma] * (image.ndim - 1)
+        if len(sigma) != image.ndim:
+            sigma = np.concatenate((np.asarray(sigma), [0]))
+    image = img_as_float(image)
+    return ndimage.gaussian_filter(image, sigma, mode=mode, cval=cval)

skimage/filter/tests/test_gaussian.py

@@ -0,0 +1,42 @@
+import numpy as np
+from skimage.filter._gaussian import gaussian_filter
+
+
+def test_null_sigma():
+    a = np.zeros((3, 3))
+    a[1, 1] = 1.
+    assert np.all(gaussian_filter(a, 0) == a)
+
+
+def test_energy_decrease():
+    a = np.zeros((3, 3))
+    a[1, 1] = 1.
+    gaussian_a = gaussian_filter(a, sigma=1, mode='reflect')
+    assert gaussian_a.std() < a.std()
+
+
+def test_multichannel():
+    a = np.zeros((5, 5, 3))
+    a[1, 1] = np.arange(1, 4)
+    gaussian_rgb_a = gaussian_filter(a, sigma=1, mode='reflect',
+                                     multichannel=True)
+    # Check that the mean value is conserved in each channel
+    # (color channels are not mixed together)
+    assert np.allclose([a[..., i].mean() for i in range(3)],
+                        [gaussian_rgb_a[..., i].mean() for i in range(3)])
+    # Test multichannel = None
+    gaussian_rgb_a = gaussian_filter(a, sigma=1, mode='reflect')
+    # Check that the mean value is conserved in each channel
+    # (color channels are not mixed together)
+    assert np.allclose([a[..., i].mean() for i in range(3)],
+                        [gaussian_rgb_a[..., i].mean() for i in range(3)])
+    # Iterable sigma
+    gaussian_rgb_a = gaussian_filter(a, sigma=[1, 2], mode='reflect',
+                                     multichannel=True)
+    assert np.allclose([a[..., i].mean() for i in range(3)],
+                        [gaussian_rgb_a[..., i].mean() for i in range(3)])
+
+
+if __name__ == "__main__":
+    from numpy import testing
+    testing.run_module_suite()