diff --git a/doc/examples/plot_pyramid.py b/doc/examples/plot_pyramid.py
new file mode 100644
index 00000000..60067549
--- /dev/null
+++ b/doc/examples/plot_pyramid.py
@@ -0,0 +1,34 @@
+"""
+====================
+Build image pyramids
+====================
+
+The `build_gauassian_pyramid` function takes an image and yields successive
+images shrunk by a constant scale factor. Image pyramids are often used, e.g.,
+to implement algorithms for denoising, texture discrimination, and scale-
+invariant detection.
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+from skimage import data
+from skimage.transform import build_gaussian_pyramid
+
+
+image = data.lena()
+rows, cols, dim = image.shape
+pyramid = tuple(build_gaussian_pyramid(image, downscale=2))
+
+composite_image = np.zeros((rows, cols + cols / 2, 3), dtype=np.double)
+
+composite_image[:rows, :cols, :] = pyramid[0]
+
+i_row = 0
+for p in pyramid[1:]:
+    n_rows, n_cols = p.shape[:2]
+    composite_image[i_row:i_row + n_rows, cols:cols + n_cols] = p
+    i_row += n_rows
+
+plt.imshow(composite_image)
+plt.show()
diff --git a/skimage/transform/__init__.py b/skimage/transform/__init__.py
index a25ce969..fc2b12ef 100644
--- a/skimage/transform/__init__.py
+++ b/skimage/transform/__init__.py
@@ -4,7 +4,8 @@ from .finite_radon_transform import *
 from .integral import *
 from ._geometric import (warp, warp_coords, estimate_transform,
                          SimilarityTransform, AffineTransform,
-                         ProjectiveTransform, PolynomialTransform, 
+                         ProjectiveTransform, PolynomialTransform,
                          PiecewiseAffineTransform)
 from ._warps import swirl, homography, resize, rotate
-
+from .pyramids import (pyramid_reduce, pyramid_expand,
+                       build_gaussian_pyramid, build_laplacian_pyramid)
diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py
index f81e6e22..c3f21db9 100644
--- a/skimage/transform/_geometric.py
+++ b/skimage/transform/_geometric.py
@@ -1016,5 +1016,7 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
     if mode == 'constant' and not (0 <= cval <= 1):
         clipped[out == cval] = cval
 
-    # Remove singleton dim introduced by atleast_3d
-    return clipped.squeeze()
+    if clipped.shape[0] == 1 or clipped.shape[1] == 1:
+        return clipped
+    else: # remove singleton dim introduced by atleast_3d
+        return clipped.squeeze()
diff --git a/skimage/transform/pyramids.py b/skimage/transform/pyramids.py
new file mode 100644
index 00000000..f07fb1c5
--- /dev/null
+++ b/skimage/transform/pyramids.py
@@ -0,0 +1,290 @@
+import math
+import numpy as np
+from scipy import ndimage
+from skimage.transform import resize
+from skimage.util import img_as_float
+
+
+def _smooth(image, sigma, mode, cval):
+    """Return image with each channel smoothed by the gaussian filter."""
+
+    smoothed = np.empty(image.shape, dtype=np.double)
+
+    if image.ndim == 3: # apply gaussian filter to all dimensions independently
+        for dim in range(image.shape[2]):
+            ndimage.gaussian_filter(image[..., dim], sigma,
+                                    output=smoothed[..., dim],
+                                    mode=mode, cval=cval)
+    else:
+        ndimage.gaussian_filter(image, sigma, output=smoothed,
+                                mode=mode, cval=cval)
+
+    return smoothed
+
+
+def _check_factor(factor):
+    if factor <= 1:
+        raise ValueError('scale factor must be greater than 1')
+
+
+def pyramid_reduce(image, downscale=2, sigma=None, order=1,
+                   mode='reflect', cval=0):
+    """Smooth and then downsample image.
+
+    Parameters
+    ----------
+    image : array
+        Input image.
+    downscale : float, optional
+        Downscale factor.
+    sigma : float, optional
+        Sigma for gaussian filter. Default is `2 * downscale / 6.0` which
+        corresponds to a filter mask twice the size of the scale factor that
+        covers more than 99% of the gaussian distribution.
+    order : int, optional
+        Order of splines used in interpolation of downsampling. See
+        `scipy.ndimage.map_coordinates` for detail.
+    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'.
+    cval : float, optional
+        Value to fill past edges of input if mode is 'constant'.
+
+    Returns
+    -------
+    out : array
+        Smoothed and downsampled float image.
+
+    References
+    ----------
+    ..[1] http://web.mit.edu/persci/people/adelson/pub_pdfs/pyramid83.pdf
+
+    """
+
+    _check_factor(downscale)
+
+    image = img_as_float(image)
+
+    rows = image.shape[0]
+    cols = image.shape[1]
+    out_rows = math.ceil(rows / float(downscale))
+    out_cols = math.ceil(cols / float(downscale))
+
+    if sigma is None:
+        # automatically determine sigma which covers > 99% of distribution
+        sigma = 2 * downscale / 6.0
+
+    smoothed = _smooth(image, sigma, mode, cval)
+    out = resize(smoothed, (out_rows, out_cols), order=order,
+                 mode=mode, cval=cval)
+
+    return out
+
+
+def pyramid_expand(image, upscale=2, sigma=None, order=1,
+                   mode='reflect', cval=0):
+    """Upsample and then smooth image.
+
+    Parameters
+    ----------
+    image : array
+        Input image.
+    upscale : float, optional
+        Upscale factor.
+    sigma : float, optional
+        Sigma for gaussian filter. Default is `2 * upscale / 6.0` which
+        corresponds to a filter mask twice the size of the scale factor that
+        covers more than 99% of the gaussian distribution.
+    order : int, optional
+        Order of splines used in interpolation of upsampling. See
+        `scipy.ndimage.map_coordinates` for detail.
+    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'.
+    cval : float, optional
+        Value to fill past edges of input if mode is 'constant'.
+
+    Returns
+    -------
+    out : array
+        Upsampled and smoothed float image.
+
+    References
+    ----------
+    ..[1] http://web.mit.edu/persci/people/adelson/pub_pdfs/pyramid83.pdf
+
+    """
+
+    _check_factor(upscale)
+
+    image = img_as_float(image)
+
+    rows = image.shape[0]
+    cols = image.shape[1]
+    out_rows = math.ceil(upscale * rows)
+    out_cols = math.ceil(upscale * cols)
+
+    if sigma is None:
+        # automatically determine sigma which covers > 99% of distribution
+        sigma = 2 * upscale / 6.0
+
+    resized = resize(image, (out_rows, out_cols), order=order,
+                     mode=mode, cval=cval)
+    out = _smooth(resized, sigma, mode, cval)
+
+    return out
+
+
+def build_gaussian_pyramid(image, max_layer=-1, downscale=2, sigma=None,
+                           order=1, mode='reflect', cval=0):
+    """Yield images of the gaussian pyramid formed by the input image.
+
+    Recursively applies the `pyramid_reduce` function to the image, and yields
+    the downscaled images. Note that the first image of the pyramid will be the
+    original, unscaled image.
+
+    Parameters
+    ----------
+    image : array
+        Input image.
+    max_layer : int
+        Number of layers for the pyramid. 0th layer is the original image.
+        Default is -1 which builds all possible layers.
+    downscale : float, optional
+        Downscale factor.
+    sigma : float, optional
+        Sigma for gaussian filter. Default is `2 * downscale / 6.0` which
+        corresponds to a filter mask twice the size of the scale factor that
+        covers more than 99% of the gaussian distribution.
+    order : int, optional
+        Order of splines used in interpolation of downsampling. See
+        `scipy.ndimage.map_coordinates` for detail.
+    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'.
+    cval : float, optional
+        Value to fill past edges of input if mode is 'constant'.
+
+    Returns
+    -------
+    pyramid : generator
+        Generator yielding pyramid layers as float images.
+
+    References
+    ----------
+    ..[1] http://web.mit.edu/persci/people/adelson/pub_pdfs/pyramid83.pdf
+
+    """
+
+    _check_factor(downscale)
+
+    # cast to float for consistent data type in pyramid
+    image = img_as_float(image)
+
+    layer = 0
+    rows = image.shape[0]
+    cols = image.shape[1]
+
+    prev_layer_image = image
+    yield image
+
+    # build downsampled images until max_layer is reached or downscale process
+    # does not change image size
+    while layer != max_layer:
+        layer += 1
+
+        layer_image = pyramid_reduce(prev_layer_image, downscale, sigma, order,
+                                     mode, cval)
+
+        prev_rows = rows
+        prev_cols = cols
+        prev_layer_image = layer_image
+        rows = layer_image.shape[0]
+        cols = layer_image.shape[1]
+
+        # no change to previous pyramid layer
+        if prev_rows == rows and prev_cols == cols:
+            break
+
+        yield layer_image
+
+
+def build_laplacian_pyramid(image, max_layer=-1, downscale=2, sigma=None,
+                            order=1, mode='reflect', cval=0):
+    """Yield images of the laplacian pyramid formed by the input image.
+
+    Each layer contains the difference between the downsampled and the
+    downsampled plus smoothed image.
+
+    Parameters
+    ----------
+    image : array
+        Input image.
+    max_layer : int
+        Number of layers for the pyramid. 0th layer is the original image.
+        Default is -1 which builds all possible layers.
+    downscale : float, optional
+        Downscale factor.
+    sigma : float, optional
+        Sigma for gaussian filter. Default is `2 * downscale / 6.0` which
+        corresponds to a filter mask twice the size of the scale factor that
+        covers more than 99% of the gaussian distribution.
+    order : int, optional
+        Order of splines used in interpolation of downsampling. See
+        `scipy.ndimage.map_coordinates` for detail.
+    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'.
+    cval : float, optional
+        Value to fill past edges of input if mode is 'constant'.
+
+    Returns
+    -------
+    pyramid : generator
+        Generator yielding pyramid layers as float images.
+
+    References
+    ----------
+    ..[1] http://web.mit.edu/persci/people/adelson/pub_pdfs/pyramid83.pdf
+    ..[2] http://sepwww.stanford.edu/~morgan/texturematch/paper_html/node3.html
+
+    """
+
+    _check_factor(downscale)
+
+    # cast to float for consistent data type in pyramid
+    image = img_as_float(image)
+
+    if sigma is None:
+        # automatically determine sigma which covers > 99% of distribution
+        sigma = 2 * downscale / 6.0
+
+    layer = 0
+    rows = image.shape[0]
+    cols = image.shape[1]
+
+    smoothed_image = _smooth(image, sigma, mode, cval)
+    yield image - smoothed_image
+
+    # build downsampled images until max_layer is reached or downscale process
+    # does not change image size
+    while layer != max_layer:
+        layer += 1
+
+        out_rows = math.ceil(rows / float(downscale))
+        out_cols = math.ceil(cols / float(downscale))
+
+        resized_image = resize(smoothed_image, (out_rows, out_cols),
+                               order=order, mode=mode, cval=cval)
+        smoothed_image = _smooth(resized_image, sigma, mode, cval)
+
+        prev_rows = rows
+        prev_cols = cols
+        rows = resized_image.shape[0]
+        cols = resized_image.shape[1]
+
+        # no change to previous pyramid layer
+        if prev_rows == rows and prev_cols == cols:
+            break
+
+        yield resized_image - smoothed_image
diff --git a/skimage/transform/tests/test_pyramids.py b/skimage/transform/tests/test_pyramids.py
new file mode 100644
index 00000000..9aa236ff
--- /dev/null
+++ b/skimage/transform/tests/test_pyramids.py
@@ -0,0 +1,41 @@
+from numpy.testing import assert_array_equal, run_module_suite
+from skimage import data
+from skimage.transform import (pyramid_reduce, pyramid_expand,
+                               build_gaussian_pyramid, build_laplacian_pyramid)
+
+
+image = data.lena()
+
+
+def test_pyramid_reduce():
+    rows, cols, dim = image.shape
+    out = pyramid_reduce(image, downscale=2)
+    assert_array_equal(out.shape, (rows / 2, cols / 2, dim))
+
+
+def test_pyramid_expand():
+    rows, cols, dim = image.shape
+    out = pyramid_expand(image, upscale=2)
+    assert_array_equal(out.shape, (rows * 2, cols * 2, dim))
+
+
+def test_build_gaussian_pyramid():
+    rows, cols, dim = image.shape
+    pyramid = build_gaussian_pyramid(image, downscale=2)
+
+    for layer, out in enumerate(pyramid):
+        layer_shape = (rows / 2 ** layer, cols / 2 ** layer, dim)
+        assert_array_equal(out.shape, layer_shape)
+
+
+def test_build_laplacian_pyramid():
+    rows, cols, dim = image.shape
+    pyramid = build_laplacian_pyramid(image, downscale=2)
+
+    for layer, out in enumerate(pyramid):
+        layer_shape = (rows / 2 ** layer, cols / 2 ** layer, dim)
+        assert_array_equal(out.shape, layer_shape)
+
+
+if __name__ == "__main__":
+    run_module_suite()
diff --git a/viewer_examples/viewers/collection_viewer.py b/viewer_examples/viewers/collection_viewer.py
index b36ddf82..24c97ab0 100644
--- a/viewer_examples/viewers/collection_viewer.py
+++ b/viewer_examples/viewers/collection_viewer.py
@@ -4,7 +4,7 @@ CollectionViewer demo
 =====================
 
 Demo of CollectionViewer for viewing collections of images. This demo uses
-successively darker versions of the same image to fake an image collection.
+the different layers of the gaussian pyramid as image collection.
 
 You can scroll through images with the slider, or you can interact with the
 viewer using your keyboard:
@@ -21,9 +21,11 @@ home/end keys
 import numpy as np
 from skimage import data
 from skimage.viewer import CollectionViewer
+from skimage.transform import build_gaussian_pyramid
+
 
 img = data.lena()
-img_collection = [np.uint8(img * 0.9**i) for i in range(20)]
+img_collection = tuple(build_gaussian_pyramid(img))
 
 view = CollectionViewer(img_collection)
 view.show()