diff --git a/skimage/filters/edges.py b/skimage/filters/edges.py
index 501dd8ea..4cf083a3 100644
--- a/skimage/filters/edges.py
+++ b/skimage/filters/edges.py
@@ -14,6 +14,7 @@ from .. import img_as_float
 from .._shared.utils import assert_nD, deprecated
 from scipy.ndimage import convolve, binary_erosion, generate_binary_structure
 
+from ..restoration.uft import laplacian
 
 EROSION_SELEM = generate_binary_structure(2, 2)
 
@@ -37,9 +38,6 @@ ROBERTS_PD_WEIGHTS = np.array([[1, 0],
 ROBERTS_ND_WEIGHTS = np.array([[0, 1],
                                [-1, 0]], dtype=np.double)
 
-LAPLACE_WEIGHTS = np.array([[1,  1,  1],
-                            [1, -8,  1],
-                            [1,  1,  1]]) / 16.0
 
 def _mask_filter_result(result, mask):
     """Return result after masking.
@@ -178,6 +176,7 @@ def hsobel(image, mask=None):
 
     Parameters
     ----------
+
     image : 2-D array
         Image to process.
     mask : 2-D array, optional
@@ -766,33 +765,34 @@ def roberts_negative_diagonal(image, mask=None):
     """
     return np.abs(roberts_neg_diag(image, mask))
 
-def laplace(image, mask=None):
+def laplace(image, ksize=3, mask=None):
     """Find the edges of an image using the Laplace operator.
 
     Parameters
     ----------
-    image : 2-D array
+    image : ndarray
         Image to process.
-    mask : 2-D array, optional
+    ksize : int, optional
+        Define the size of the discrete Laplacian operator such that it
+        will have a size of (ksize,) * image.ndim.
+    mask : ndarray, optional
         An optional mask to limit the application to a certain area.
         Note that pixels surrounding masked regions are also masked to
         prevent masked regions from affecting the result.
 
     Returns
     -------
-    output : 2-D array
+    output : ndarray
         The Laplace edge map.
 
     Notes
     -----
-    We use the following kernel::
-
-      1   1   1
-      1  -8   1
-      1   1   1
+    The Laplacian operator is generated using the function 
+    skimage.restoration.uft.laplacian().
 
     """
-    assert_nD(image, 2)
     image = img_as_float(image)
-    result = convolve(image, LAPLACE_WEIGHTS)
+    # Create the discrete Laplacian operator - We keep only the real part of the filter
+    _, laplace_op = laplacian(image.ndim, (ksize, ) * image.ndim)
+    result = convolve(image, laplace_op)
     return _mask_filter_result(result, mask)
diff --git a/skimage/filters/tests/test_edges.py b/skimage/filters/tests/test_edges.py
index 2754cbe8..1fd5e32f 100644
--- a/skimage/filters/tests/test_edges.py
+++ b/skimage/filters/tests/test_edges.py
@@ -337,15 +337,29 @@ def test_vprewitt_horizontal():
 
 def test_laplace_zeros():
     """Laplace on an array of all zeros."""
-    result = filters.laplace(np.zeros((10, 10)), np.ones((10, 10), bool))
-    assert (np.all(result == 0))
+    # Create a synthetic 2D image
+    image = np.zeros((9,9))
+    image[3:-3] = 1
+    result = filters.laplace(image)
+    res_chk = array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
+                     [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
+                     [ 0.,  0.,  0., -1., -1., -1.,  0.,  0.,  0.],
+                     [ 0.,  0., -1.,  2.,  1.,  2., -1.,  0.,  0.],
+                     [ 0.,  0., -1.,  1.,  0.,  1., -1.,  0.,  0.],
+                     [ 0.,  0., -1.,  2.,  1.,  2., -1.,  0.,  0.],
+                     [ 0.,  0.,  0., -1., -1., -1.,  0.,  0.,  0.],
+                     [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
+                     [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.]])
+    assert_allclose(result, res_chk)
 
 
 def test_laplace_mask():
     """Laplace on a masked array should be zero."""
-    np.random.seed(0)
-    result = filters.laplace(np.random.uniform(size=(10, 10)),
-                             np.zeros((10, 10), bool))
+    # Create a synthetic 2D image
+    image = np.zeros((9, 9))
+    image[3:-3] = 1
+    # Define the mask
+    result = filters.laplace(image, np.zeros((10, 10), bool))
     assert (np.all(result == 0))