Morphological operations¶

Morphology is the study of shapes. In image processing, some simple operations can get you a long way. The first things to learn are erosion and dilation. In erosion, we look at a pixel’s local neighborhood and replace the value of that pixel with the minimum value of that neighborhood. In dilation, we instead choose the maximum.

import numpy as np
from matplotlib import pyplot as plt, cm
import skdemo

image = np.array([[0, 0, 0, 0, 0, 0, 0],
                  [0, 0, 0, 0, 0, 0, 0],
                  [0, 0, 1, 1, 1, 0, 0],
                  [0, 0, 1, 1, 1, 0, 0],
                  [0, 0, 1, 1, 1, 0, 0],
                  [0, 0, 0, 0, 0, 0, 0],
                  [0, 0, 0, 0, 0, 0, 0]], dtype=np.uint8)
plt.imshow(image);

../_images/3_morphological_operations_3_0.png

The documentation for scikit-image’s morphology module is here.

Importantly, we must use a structuring element, which defines the local neighborhood of each pixel. To get every neighbor (up, down, left, right, and diagonals), use morphology.square; to avoid diagonals, use morphology.diamond:

from skimage import morphology
sq = morphology.square(width=3)
dia = morphology.diamond(radius=1)
print(sq)
print(dia)

[[1 1 1]
 [1 1 1]
 [1 1 1]]
[[0 1 0]
 [1 1 1]
 [0 1 0]]

The central value of the structuring element represents the pixel being considered, and the surrounding values are the neighbors: a 1 value means that pixel counts as a neighbor, while a 0 value does not. So:

skdemo.imshow_all(image, morphology.erosion(image, sq), shape=(1, 2))

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-5-7101991a2e67> in <module>
----> 1 skdemo.imshow_all(image, morphology.erosion(image, sq), shape=(1, 2))

~/work/skimage-tutorials/skimage-tutorials/lectures/skdemo/_skdemo.py in imshow_all(*images, **kwargs)
     83     fig, axes = plt.subplots(nrows=nrows, ncols=ncols, figsize=(width, size))
     84     for ax, img, label in zip(axes.ravel(), images, titles):
---> 85         ax.imshow(img, **kwargs)
     86         ax.set_title(label)
     87 

/opt/hostedtoolcache/Python/3.8.10/x64/lib/python3.8/site-packages/matplotlib/__init__.py in inner(ax, data, *args, **kwargs)
   1359     def inner(ax, *args, data=None, **kwargs):
   1360         if data is None:
-> 1361             return func(ax, *map(sanitize_sequence, args), **kwargs)
   1362 
   1363         bound = new_sig.bind(ax, *args, **kwargs)

/opt/hostedtoolcache/Python/3.8.10/x64/lib/python3.8/site-packages/matplotlib/axes/_axes.py in imshow(self, X, cmap, norm, aspect, interpolation, alpha, vmin, vmax, origin, extent, filternorm, filterrad, resample, url, **kwargs)
   5603             aspect = rcParams['image.aspect']
   5604         self.set_aspect(aspect)
-> 5605         im = mimage.AxesImage(self, cmap, norm, interpolation, origin, extent,
   5606                               filternorm=filternorm, filterrad=filterrad,
   5607                               resample=resample, **kwargs)

/opt/hostedtoolcache/Python/3.8.10/x64/lib/python3.8/site-packages/matplotlib/image.py in __init__(self, ax, cmap, norm, interpolation, origin, extent, filternorm, filterrad, resample, **kwargs)
    898         self._extent = extent
    899 
--> 900         super().__init__(
    901             ax,
    902             cmap=cmap,

/opt/hostedtoolcache/Python/3.8.10/x64/lib/python3.8/site-packages/matplotlib/image.py in __init__(self, ax, cmap, norm, interpolation, origin, filternorm, filterrad, resample, **kwargs)
    255         self._imcache = None
    256 
--> 257         self.update(kwargs)
    258 
    259     def __getstate__(self):

/opt/hostedtoolcache/Python/3.8.10/x64/lib/python3.8/site-packages/matplotlib/artist.py in update(self, props)
   1060                     func = getattr(self, f"set_{k}", None)
   1061                     if not callable(func):
-> 1062                         raise AttributeError(f"{type(self).__name__!r} object "
   1063                                              f"has no property {k!r}")
   1064                     ret.append(func(v))

AttributeError: 'AxesImage' object has no property 'shape'

../_images/3_morphological_operations_7_1.png

and

skdemo.imshow_all(image, morphology.dilation(image, sq))

and

skdemo.imshow_all(image, morphology.dilation(image, dia))

Erosion and dilation can be combined into two slightly more sophisticated operations, opening and closing. Here’s an example:

image = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                  [0, 0, 1, 1, 1, 1, 1, 1, 0, 0],
                  [0, 0, 1, 1, 1, 1, 1, 1, 0, 0],
                  [0, 0, 1, 1, 1, 1, 1, 1, 0, 0],
                  [0, 0, 1, 1, 1, 0, 0, 1, 0, 0],
                  [0, 0, 1, 1, 1, 0, 0, 1, 0, 0],
                  [0, 0, 1, 1, 1, 0, 0, 1, 0, 0],
                  [0, 0, 1, 1, 1, 1, 1, 1, 0, 0],
                  [0, 0, 1, 1, 1, 1, 1, 1, 0, 0],
                  [0, 0, 1, 1, 1, 1, 1, 1, 0, 0],
                  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
                  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], np.uint8)
plt.imshow(image);

What happens when run an erosion followed by a dilation of this image?

What about the reverse?

Try to imagine the operations in your head before trying them out below.

skdemo.imshow_all(image, morphology.opening(image, sq)) # erosion -> dilation

skdemo.imshow_all(image, morphology.closing(image, sq)) # dilation -> erosion

Exercise: use morphological operations to remove noise from a binary image.

from skimage import data, color
hub = color.rgb2gray(data.hubble_deep_field()[350:450, 90:190])
plt.imshow(hub);

Remove the smaller objects to retrieve the large galaxy.

Image analysis in Python

Morphological operations¶