diff --git a/skimage/data/__init__.py b/skimage/data/__init__.py
index 2be829ac..1711185c 100644
--- a/skimage/data/__init__.py
+++ b/skimage/data/__init__.py
@@ -10,7 +10,7 @@ import os as _os
 
 from .. import data_dir
 from ..io import imread, use_plugin
-
+from ._binary_blobs import binary_blobs
 
 __all__ = ['load',
            'camera',
diff --git a/skimage/data/_binary_blobs.py b/skimage/data/_binary_blobs.py
new file mode 100644
index 00000000..38977071
--- /dev/null
+++ b/skimage/data/_binary_blobs.py
@@ -0,0 +1,49 @@
+import numpy as np
+from ..filters import gaussian_filter
+
+
+def binary_blobs(length=512, blob_size_fraction=0.1, n_dim=2,
+                 volume_fraction=0.5, seed=None):
+    """
+    Generate synthetic binary image with several blob-like rounded objects.
+
+    Parameters
+    ----------
+    length : int, default 512
+        Linear size of output image.
+    blob_size_fraction : float, default 0.1
+        Typical linear size of blob, as a fraction of ``length``, should be
+        smaller than 1.
+    n_dim : int, default 2
+        Number of dimensions of output image.
+    volume_fraction : float, default 0.5
+        Fraction of image pixels covered by the blobs (where the output is 1).
+        Should be in [0, 1].
+    seed : int, default 0
+        Seed to initialize the random number generator.
+
+    Returns
+    -------
+    blobs : ndarray of bools
+        Output binary image
+
+    Examples
+    --------
+    >>> blobs = binary_blobs(length=256, blob_size_fraction=0.1)
+    >>> # Finer structures
+    >>> blobs = binary_blobs(length=256, blob_size_fraction=0.05)
+    >>> # Blobs cover a smaller volume fraction of the image
+    >>> blobs = binary_blobs(length=256, volume_fraction=0.3)
+    """
+    if seed is None:
+        seed = 0
+    # Fix the seed for reproducible results
+    rs = np.random.RandomState(seed)
+    shape = tuple([length] * n_dim)
+    mask = np.zeros(shape)
+    n_pts = max(int(1. / blob_size_fraction) ** n_dim, 1)
+    points = (length * rs.rand(n_dim, n_pts)).astype(np.int)
+    mask[[indices for indices in points]] = 1
+    mask = gaussian_filter(mask, sigma=0.25 * length * blob_size_fraction)
+    threshold = np.percentile(mask, 100 * (1 - volume_fraction))
+    return np.logical_not(mask < threshold)
diff --git a/skimage/data/tests/test_data.py b/skimage/data/tests/test_data.py
index bfcdc571..139f1ddb 100644
--- a/skimage/data/tests/test_data.py
+++ b/skimage/data/tests/test_data.py
@@ -54,6 +54,14 @@ def test_coffee():
     data.coffee()
 
 
+def test_binary_blobs():
+    blobs = data.binary_blobs(length=128)
+    assert blobs.mean() == 0.5
+    blobs = data.binary_blobs(length=128, volume_fraction=0.25)
+    assert blobs.mean() == 0.25
+    blobs = data.binary_blobs(length=32, volume_fraction=0.25, n_dim=3)
+    assert blobs.mean() == 0.25
+
 if __name__ == "__main__":
     from numpy.testing import run_module_suite
     run_module_suite()