Merge pull request #484 from almarklein/3d-selem

Adding predefined structuring elements for 3D morphology.

skimage/morphology/selem.py

@@ -93,7 +93,7 @@ def disk(radius, dtype=np.uint8):
     """
     Generates a flat, disk-shaped structuring element of a given radius.
     A pixel is within the neighborhood if the euclidean distance between
-    it and the origin is no greater than a radius.
+    it and the origin is no greater than radius.
 
     Parameters
     ----------
@@ -114,3 +114,90 @@ def disk(radius, dtype=np.uint8):
     s = X**2
     s += Y**2
     return np.array(s <= radius * radius, dtype=dtype)
+
+
+def cube(width, dtype=np.uint8):
+    """
+    Generates a cube-shaped structuring element (the 3D equivalent of
+    a square). Every pixel along the perimeter has a chessboard distance
+    no greater than radius (radius=floor(width/2)) pixels.
+
+    Parameters
+    ----------
+    width : int
+       The width, height and depth of the cube
+
+    Other Parameters
+    ----------------
+    dtype : data-type
+       The data type of the structuring element.
+
+    Returns
+    -------
+    selem : ndarray
+       A structuring element consisting only of ones, i.e. every
+       pixel belongs to the neighborhood.
+
+    """
+    return np.ones((width, width, width), dtype=dtype)
+
+
+def octahedron(radius, dtype=np.uint8):
+    """
+    Generates a octahedron-shaped structuring element of a given radius
+    (the 3D equivalent of a diamond).  A pixel is part of the
+    neighborhood (i.e. labeled 1) if the city block/manhattan distance
+    between it and the center of the neighborhood is no greater than
+    radius.
+
+    Parameters
+    ----------
+    radius : int
+       The radius of the octahedron-shaped structuring element.
+
+    dtype : data-type
+       The data type of the structuring element.
+
+    Returns
+    -------
+
+    selem : ndarray
+       The structuring element where elements of the neighborhood
+       are 1 and 0 otherwise.
+    """
+    # note that in contrast to diamond(), this method allows non-integer radii
+    n = 2 * radius + 1
+    Z, Y, X = np.mgrid[ -radius:radius:n*1j, 
+                        -radius:radius:n*1j,
+                        -radius:radius:n*1j]
+    s = np.abs(X) + np.abs(Y) + np.abs(Z) 
+    return np.array(s <= radius, dtype=dtype)
+
+
+def ball(radius, dtype=np.uint8):
+    """
+    Generates a ball-shaped structuring element of a given radius (the
+    3D equivalent of a disk). A pixel is within the neighborhood if the
+    euclidean distance between it and the origin is no greater than
+    radius.
+
+    Parameters
+    ----------
+    radius : int
+       The radius of the ball-shaped structuring element.
+
+    dtype : data-type
+       The data type of the structuring element.
+
+    Returns
+    -------
+    selem : ndarray
+       The structuring element where elements of the neighborhood
+       are 1 and 0 otherwise.
+    """
+    n = 2 * radius + 1
+    Z, Y, X = np.mgrid[ -radius:radius:n*1j, 
+                        -radius:radius:n*1j, 
+                        -radius:radius:n*1j]
+    s = X**2 + Y**2 + Z**2
+    return np.array(s <= radius * radius, dtype=dtype)

skimage/morphology/tests/test_selem.py

@@ -36,9 +36,32 @@ class TestSElem():
                 expected_mask = expected_mask[:, np.newaxis]
             assert_equal(expected_mask, actual_mask)
             k = k + 1
-
+    
+    def strel_worker_3d(self, fn, func):
+        matlab_masks = np.load(os.path.join(data_dir, fn))
+        k = 0
+        for arrname in sorted(matlab_masks):
+            expected_mask = matlab_masks[arrname]
+            actual_mask = func(k)
+            if (expected_mask.shape == (1,)):
+                expected_mask = expected_mask[:, np.newaxis]
+            # Test center slice for each dimension. This gives a good
+            # indication of validity without the need for a 3D reference
+            # mask.
+            c = int(expected_mask.shape[0]/2) 
+            assert_equal(expected_mask, actual_mask[c,:,:])
+            assert_equal(expected_mask, actual_mask[:,c,:])
+            assert_equal(expected_mask, actual_mask[:,:,c])
+            k = k + 1
+    
     def test_selem_disk(self):
         self.strel_worker("disk-matlab-output.npz", selem.disk)
 
     def test_selem_diamond(self):
         self.strel_worker("diamond-matlab-output.npz", selem.diamond)
+    
+    def test_selem_ball(self):
+        self.strel_worker_3d("disk-matlab-output.npz", selem.ball)
+
+    def test_selem_octahedron(self):
+        self.strel_worker_3d("diamond-matlab-output.npz", selem.octahedron)