Segmenting Python array into unique regions connected by a single cell or less?

Question

I have a numpy array which I wish to segment into discrete regions with unique IDs which looks something like this:

Usually for something like this I would use scipy.ndimage.label to generate unique ids for discrete blobs, but in this case I have several very large continuous regions which I also wish to be segmented into smaller unique regions, ideally when they are only joined to their neighbours by a connection 1 cell wide. To illustrate, here's a sample array, the result I get when running scipy.ndimage.label, and my desired result:

import numpy as np
import scipy.ndimage as ndimage
import matplotlib.pyplot as plt

# Sample study area array
example_array = np.array([[0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0],
                          [0, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 1],
                          [0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1],
                          [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0],
                          [1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1],
                          [1, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1],
                          [1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1],
                          [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0],
                          [1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0],
                          [1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0],
                          [1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0],
                          [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])

# Identify discrete regions and assign unique IDs
current_output, num_ids = ndimage.label(example_array, structure=np.ones((3,3)))

# Plot outputs
plt.imshow(example_array, cmap="spectral", interpolation='nearest')
plt.imshow(current_output, cmap="spectral", interpolation='nearest')

The specific numbering and division of cells isn't consequential; anything resembling the last plot above will do. My best attempt so far has been to use morphological image processing to first erode my array, run scipy.ndimage.label and then dilate, but this has the unfortunate side effect of eliminating all single cell regions or thin linear features (of which there are many).

Would greatly appreciate any thoughts!

Answer 1

Easiest thing to do may be to apply the following prepossessing kernel before the SciPy labeling (with 2,2 as the origin and "?" can be 0 or 1):

IF |? 1 ?| OR |? 0 ?| THEN origin(x,y) == 0
   |0 1 0|    |1 1 1|
   |? 1 ?|    |? 0 ?|

IMPLEMENTATION:

import numpy as np
import scipy.ndimage as ndimage
import matplotlib.pyplot as plt

# Sample study area array
example_array = np.array([[0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0],
                          [0, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 1],
                          [0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1],
                          [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0],
                          [1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1],
                          [1, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1],
                          [1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1],
                          [1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0],
                          [1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0],
                          [1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0],
                          [1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0],
                          [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])

# Get array shape
row_count, column_count = example_array.shape

# Iterate through kernel origins
for (x,y), value in np.ndenumerate(example_array):

    # Skip first and last rows and columns because kernel would be out of bounds
    if not (x == 0 
         or x == row_count-1 
         or y == 0 
         or y == column_count-1):

        # Extract kernel
        kernel = example_array[x-1:x+2, y-1:y+2]
      
        # Apply IF |? 1 ?| OR |? 0 ?| THEN origin(x,y) == 0
        #          |0 1 0|    |1 1 1|
        #          |? 1 ?|    |? 0 ?|
        if ((kernel[1,0:3]).all() == 1 and kernel[0,1] == 0 and kernel[2,1] == 0 
            or (kernel[0:3,1]).all() == 1 and kernel[1,0] == 0 and kernel[1,2] == 0):
            example_array[x,y] = 0
                
# Identify discrete regions and assign unique IDs
current_output, num_ids = ndimage.label(example_array, structure=np.ones((3,3)))

# Plot outputs
plt.imshow(example_array, cmap="spectral", interpolation='nearest')
plt.imshow(current_output, cmap="spectral", interpolation='nearest')

RESULT: