Skimage skeletonize or OpenCV implementation that respects a grid?

Question

Suppose one has lines intersecting each other at right angles.

And you would like to skeletonize it to obtain (you hope) a cross shape. Instead, with sklearn.morphology.skeletonize the following image is obtained:

Let's call it a "holey cross".

On the other hand, you have OpenCV and the OpenCV skeletonize function that is floating around on the internet in several blogs and answers on here:

def skeletonize(bin: numpy.ndarray, erosion_shape=cv2.MORPH_RECT, kernel_sz: Union[int, Tuple[int, int]] = 3):
    kernel_sz = fix_kernel(kernel_sz)
    kernel = cv2.getStructuringElement(erosion_shape, kernel_sz)

    thresh = bin.copy()
    skeleton = numpy.zeros_like(bin)
    eroded = numpy.zeros_like(bin)
    carry = numpy.zeros_like(bin)

    while (True):
        cv2.erode(thresh, kernel, dst=eroded)
        cv2.dilate(eroded, kernel, dst=carry)
        cv2.subtract(thresh, carry, dst=carry)
        cv2.bitwise_or(skeleton, carry, dst=skeleton)
        thresh, eroded = eroded, thresh

        if cv2.countNonZero(thresh) == 0:
            return skeleton

This one produces the following result:

---

So, there is something wrong or off about the basic OpenCV skeletonization function floating around, and the Skimage skeletonization cannot be modified with a structuring shape.

Is there a way to obtain the skeletonized cross/plus sign shape in python?

Answer 1

As I noted in the comments, you can clean up crossover points in a skeletonized image by fitting hough lines:

#!/usr/bin/env python
"""
https://stackoverflow.com/q/66995948/2912349
"""
import numpy as np
import matplotlib.pyplot as plt

from skimage.morphology import skeletonize
from skimage.transform import probabilistic_hough_line
from skimage.draw import line as get_line_pixels

img = np.zeros((20, 20))
img[4:16, 6:14] = 1
img[:, 10] = 1
img[10, :] = 1

skel = skeletonize(img)

lines = probabilistic_hough_line(skel, line_length=10)

# hough_line() returns the start and endpoint of the fitted lines;
# we need all pixels covered by that line;
cleaned = np.zeros_like(img)
for ((r0, c0), (r1, c1)) in lines:
    rr, cc = get_line_pixels(r0, c0, r1, c1)
    cleaned[rr, cc] = 1

fig, axes = plt.subplots(1, 3, sharex=True, sharey=True)
axes[0].imshow(img, cmap='gray')
axes[0].set_title('Raw')
axes[1].imshow(skel, cmap='gray')
axes[1].set_title('Skeleton')
axes[2].imshow(cleaned, cmap='gray')
axes[2].set_title('Hough lines')
plt.show()

If you want to force horizontal or vertical fits, lines can be trivially filtered to exclude non-horizontal and non-vertical lines:

for ((r0, c0), (r1, c1)) in lines:
    if (r0 == r1) or (c0 == c1):
        ...