Binary erosion without using SciPy ndimage

One of the best parts about Python is all of the great add-on modules that allow you to do very complex operations in just a few lines of code.  That being said, I also greatly enjoy when I can’t use a specific module or function and I am forced to really learn an algorithm in order to recreate it.  I recently needed to do a binary erosion of a NumPy array without using the binary_erosion() function in the SciPy ndimage package (the ArcGIS Python install only includes NumPy and matplotlib).

The following is an example of binary erosion using the SciPy ndimage function

>>> import numpy as np
>>> from scipy import ndimage
>>> input_array = np.ones((7,7), dtype=np.int)
>>> input_array[3,0] = 0
>>> input_array[1,3:5] = 0
>>> input_array
array([[1, 1, 1, 1, 1, 1, 1],
       [1, 1, 1, 0, 0, 1, 1],
       [1, 1, 1, 1, 1, 1, 1],
       [0, 1, 1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1, 1, 1]])
>>> structure = np.ones((3,3), dtype=np.int)
>>> ndimage.binary_erosion(input_array, structure)
array([[0, 0, 0, 0, 0, 0, 0],
       [0, 1, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0],
       [0, 0, 1, 1, 1, 1, 0],
       [0, 0, 1, 1, 1, 1, 0],
       [0, 1, 1, 1, 1, 1, 0],
       [0, 0, 0, 0, 0, 0, 0]])

I searched stackoverflow for awhile, but most of the responses were to use ndimage.  Since the algorithm itself is quite simple (set the output pixel to the minimum value of all the pixels in the input pixel’s structure/neighborhood) I decided to start from scratch, instead of borrowing from stackoverflow like I usually do.  Initially, I was going to just iterate over every cell in the input array, select the subset using the structure element as a mask, and then calculate the minimum.  The drawback to this approach is that I would need quite a bit of logic to handle the edge cases because the structure element would also need to be sliced.  Instead, I decided to just pad the input array with zeros so the structure element would fit without modification.

>>> pad_array = np.zeros((9,9), dtype=np.int)
array([[0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0]])
>>> pad_array[1:-1,1:-1] = input_array
>>> pad_array
array([[0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 1, 1, 1, 1, 1, 1, 1, 0],
       [0, 1, 1, 1, 0, 0, 1, 1, 0],
       [0, 1, 1, 1, 1, 1, 1, 1, 0],
       [0, 0, 1, 1, 1, 1, 1, 1, 0],
       [0, 1, 1, 1, 1, 1, 1, 1, 0],
       [0, 1, 1, 1, 1, 1, 1, 1, 0],
       [0, 1, 1, 1, 1, 1, 1, 1, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0]])

Padding with zeros also made it trivial to propagate the edge cells into the array like occurs with ndimage.binary_erosion().  One potential downside of this approach is that I am creating an additional copy of the input array in order to pad it.

Below is the numpy binary erosion function. I am sure that this is a very inefficient approach, but it worked!

def np_binary_erosion(input_array, structure=np.ones((3,3)).astype(np.bool)):
    '''NumPy binary erosion function
    
    No error checking on input array (type)
    No error checking on structure element (# of dimensions, shape, type, etc.)
    
    Args:
    input_array: Binary NumPy array to be eroded. Non-zero (True) elements
        form the subset to be eroded
    structure: Structuring element used for the erosion. Non-zero elements
        are considered True. If no structuring element is provided, an
        element is generated with a square connectivity equal to two 
        (square, not cross).
    Returns:
        binary_erosion: Erosion of the input by the stucturing element
    '''
    rows, cols = input_array.shape
    
    ## Pad output array (binary_erosion) with extra cells around the edge
    ## so that structuring element will fit without wrapping.
    ## A 3x3 structure, will need 1 additional cell around the edge
    ## A 5x5 structure, will need 2 additional cells around the edge
    pad_shape = (
        input_shape[0] + structure_shape[0] - 1, 
        input_shape[1] + structure.shape[1] - 1)
    input_pad_array = np.zeros(pad_shape).astype(np.bool)
    input_pad_array[1:rows+1,1:cols+1] = input_array
    binary_erosion = np.zeros(pad_shape).astype(np.bool)
    
    ## Cast structure element to boolean
    struc_mask = structure.astype(np.bool)
    
    ## Iterate over each cell
    for row in xrange(rows):
        for col in xrange(cols):
            ## The value of the output pixel is the minimum value of all the
            ## pixels in the input pixel's neighborhood.
            binary_erosion[row+1,col+1] = np.min(
                input_pad_array[row:row+3, col:col+3][struc_mask])
    return binary_erosion[1:rows+1,1:cols+1]

And the result using the function

>>> np_binary_erosion(input_array, structure).astype(np.int)
array([[0, 0, 0, 0, 0, 0, 0],
       [0, 1, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0],
       [0, 0, 1, 1, 1, 1, 0],
       [0, 0, 1, 1, 1, 1, 0],
       [0, 1, 1, 1, 1, 1, 0],
       [0, 0, 0, 0, 0, 0, 0]])

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s