diff --git a/skimage/filter/__init__.py b/skimage/filter/__init__.py
index 1e13a87e..67088b20 100644
--- a/skimage/filter/__init__.py
+++ b/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',
diff --git a/skimage/filter/_gaussian.py b/skimage/filter/_gaussian.py
new file mode 100644
index 00000000..6f0f9a64
--- /dev/null
+++ b/skimage/filter/_gaussian.py
@@ -0,0 +1,91 @@
+from scipy import ndimage
+from skimage.util import img_as_float
+
+__all__ = ['gaussian_filter']
+
+
+def gaussian_filter(image, sigma, output=None, mode='nearest', cval=0,
+                    multichannel=False):
+    """
+    Multi-dimensional Gaussian filter
+
+    Parameters
+    ----------
+
+    image : array-like
+        input image 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: False)
+        Whether the last axis of the image is to be interpreted as multiple
+        channels. 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)
+    """
+    if multichannel:
+        # do not filter across channels
+        ndim = image.ndim
+        sigma = [sigma] * (ndim - 1) + [0]
+    image = img_as_float(image)
+    return ndimage.gaussian_filter(image, sigma, mode=mode, cval=cval)
diff --git a/skimage/filter/tests/test_gaussian.py b/skimage/filter/tests/test_gaussian.py
new file mode 100644
index 00000000..afae5079
--- /dev/null
+++ b/skimage/filter/tests/test_gaussian.py
@@ -0,0 +1,29 @@
+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((3, 3, 3))
+    a[1, 1] = np.arange(1, 4)
+    gaussian_rgb_a = gaussian_filter(a, sigma=1, 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()