Re-post: Computing line that follows objects orientation in contours

Question

I am using cv2 to compute the line withing the object present in the mask image attached. The orientation/shape of the object might vary (horizontal, vertical, horizontal curved or vertical curved objects) from mask images in my dataset.

I had been suggest by @Tino D, a PCA based approach to compute the points. It worked initially for the image in question but failed to bend the line (to be drawn end to end withing an object) as orientation of an object in mask image is not always straight. it would be helpful if anyone can suggest an much flexible approach so that line is always following the orientation.

Example raw mask image where object is not straight :

here's how line was computed :

Here is what I want :

%matplotlib notebook
import matplotlib.pyplot as plt
import numpy as np
import cv2

im = cv2.imread("mask.png", 0) # read as gray
y, x = np.where(im) # get non-zero elements
centroid = np.mean(x), np.mean(y) # get the centroid of the particle
x_diff = x - centroid[0] # center x
y_diff = y - centroid[1] # center y
cov_matrix = np.cov(x_diff, y_diff) # get the convariance
eigenvalues, eigenvectors = np.linalg.eig(cov_matrix) # apply EVD
indicesForSorting = np.argsort(eigenvalues)[::-1] # sort to get the primary first
eigenvalues = eigenvalues[indicesForSorting]
eigenvectors = eigenvectors[:, indicesForSorting]
plt.figure()
plt.imshow(im, cmap = "gray") # plot image
vecPrimary = eigenvectors[:, 0] * np.sqrt(eigenvalues[0])
plt.plot([centroid[0] - vecPrimary[0], centroid[0] + vecPrimary[0]], 
        [centroid[1] - vecPrimary[1], centroid[1] + vecPrimary[1]])
vecSecondary = eigenvectors[:, 1] * np.sqrt(eigenvalues[1])
plt.plot([centroid[0] - vecSecondary[0], centroid[0] + vecSecondary[0]], 
        [centroid[1] - vecSecondary[1], centroid[1] + vecSecondary[1]])

Answer 1

The method below uses the excellent PCA-based answer by @Tino D as the key starting point, where the principal axes of the image are determined:

The image is then rotated onto its principal axes, and the mid-point of the image is demarcated along the length (green dots). A polynomial is then fitted to those points (magenta):

The curve is interpolated onto a high-resolution axis and then rotated back into the original coordinate space:

The main tunable parameter is polynomial_degree. It defaults to 3 to render a gentle curve. Higher values will track the green points more closely, giving you a more precise (but also more irregular) midline. The example below is with an 8th-degree polynomial.

You could also increase the number of points n_points used for the fit, though it likely won't matter much except if you wanted to experiment with high-degree fits.

---

from matplotlib import pyplot as plt
import numpy as np
from PIL import Image

#Load image as numpy, and binarise
arr = np.array( Image.open('mask4.png') )
if arr.ndim == 3:
    #If it's RGB, flatten the to grayscale
    arr = np.array(arr).mean(axis=2)
arr = np.where(arr > arr.max() / 2, 1, 0)

#Get the non-zero points for calculating centre of mass
y, x = np.where(arr)
xy = np.column_stack([x, y])

com = xy.mean(axis=0)
com_x, com_y = com

#View image, clipped to non-zero bounds
plt.imshow(arr, cmap='binary', interpolation='none', origin='lower')
plt.xlabel('x')
plt.ylabel('y')
plt.axis([x.min(), x.max(), y.min(), y.max()])

#Mark the COM
plt.scatter(com_x, com_y, s=80, edgecolor='white')
plt.axhline(com_y, linewidth=0.5)
plt.axvline(com_x, linewidth=0.5)

#PCA
xy_diff = xy - com
x_diff, y_diff = xy_diff.T
cov_matrix = np.cov(x_diff, y_diff)
eigenvalues, eigenvectors = np.linalg.eig(cov_matrix)
sorting_indices = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[sorting_indices]
eigenvectors = eigenvectors[:, sorting_indices]
pca0 = eigenvectors[:, 0] * eigenvalues[0] ** 0.5 * 1
pca1 = eigenvectors[:, 1] * eigenvalues[1] ** 0.5 * 1

#Overlay PCA axes
plt.plot([com_x, com_x + pca0[0]], [com_y, com_y + pca0[1]], color='tab:red')
plt.plot([com_x, com_x + pca1[0]], [com_y, com_y + pca1[1]], color='tab:red')
plt.show()

#Rotate pixels onto PCA axes, and view
xy_rot = xy_diff @ eigenvectors
x_rot, y_rot = xy_rot.T

plt.scatter(x_rot, y_rot, marker='.', color='black')
plt.gca().set_aspect('equal')

# Get the averaged x/y coordinates along the length. 
n_points = 50
points_x_brackets = np.percentile(x_rot, np.linspace(0, 100, n_points + 1))

points_x = []
points_y = []
for x_lower, x_upper in zip(points_x_brackets[:-1], points_x_brackets[1:]):
    points_x.append( x_rot[(x_rot > x_lower) & (x_rot < x_upper)].mean() )
    points_y.append( y_rot[(x_rot > x_lower) & (x_rot < x_upper)].mean() )

#Add the min and max x values
points_x = [x_rot.min()] + points_x + [x_rot.max()]
points_y = [y_rot[x_rot==x_rot.min()][0]] + points_y + [y_rot[x_rot==x_rot.max()][0]]

#View demarcations
[plt.axvline(x, linewidth=0.5, color='tab:green', linestyle=':') for x in points_x]
plt.scatter(points_x, points_y, marker='.', s=10, color='tab:green')
plt.scatter(0, 0, marker='o', s=80, color='tab:blue', edgecolor='white')

#Fit curve through the midpoints
from numpy.polynomial import Polynomial

polynomial_degree = 8
poly = Polynomial.fit(
    points_x, points_y,
    deg=polynomial_degree
)

#Create a high-resolution axis and get the curve's values over this range
pca0_ax = np.linspace(points_x[0], points_x[-1], num=len(xy_rot))
curve = poly(pca0_ax)

#View curve
plt.plot(pca0_ax, curve, linewidth=1, color='magenta')
plt.show()

#Invert the curve's coordinates onto original axes
curve_xy = np.column_stack([pca0_ax, curve])
curve_xy_inv = curve_xy @ np.linalg.pinv(eigenvectors)

#Overlay curve on the image in native space
plt.imshow(arr, cmap='binary', interpolation='none', origin='lower')
plt.xlabel('x')
plt.ylabel('y')

plt.plot(
    curve_xy_inv[:, 0] + com_x, curve_xy_inv[:, 1] + com_y,
    color='white', linestyle='--', linewidth=1.5
)
plt.axis([x.min(), x.max(), y.min(), y.max()])
plt.show()