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]])