重新发布：计算轮廓中遵循对象方向的线

Question

我正在使用 cv2 来计算所附蒙版图像中存在的对象所在的线。对象的方向/形状可能与我的数据集中的掩模图像不同（水平、垂直、水平弯曲或垂直弯曲对象）。

@Tino D 建议我使用一种基于 PCA 的方法来计算点。它最初适用于有问题的图像，但未能弯曲线（在对象内端到端地绘制），因为蒙版图像中对象的方向并不总是笔直的。如果有人可以建议一种更灵活的方法，以便该线始终遵循方向，那将会很有帮助。

对象不直的原始蒙版图像示例：

以下是直线的计算方式：

这就是我想要的：

%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

下面的方法使用 @TinoD 的基于 PCA 的答案作为起点，其中确定了图像的主轴：

然后将图像旋转到其主轴上，并确定左/右尖端位置。通过尖端位置和旋转空间中的中点拟合二次曲线。

然后将曲线旋转回原始图像的空间：

有一个可调参数，

tip_percentile

。它默认为 1%，但在一系列值下结果非常相似。 “1%”表示左尖端的位置（在旋转空间中）是最左边 1% 像素的平均值（右尖端也类似）。

由于代码只拟合二次方程，因此曲线非常平滑。如果我们拟合更多的点，曲率可以更紧密地遵循图像的轮廓（尽管它会更不规则）。

如果曲线始终太平坦，您可以向 PCA 中点坐标添加一个名为

centre_pt_y

的偏移，以获得更多的曲率。

包含两张原始图像和一张旋转版本的结果：

---

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 x/y coordiantes of the left/right tips.
# Done by finding the extreme x values (e.g. 1%, 99% percentiles) and
#  averaging the resulting coordindates
tip_percentile = 1
xtip_l, xtip_r = np.percentile(x_rot, [tip_percentile, 100 - tip_percentile])
print(f'Averaging {tip_percentile / 100 * x.size:.0f} points for each tip coordinate')

ytip_l = np.median( y_rot[x_rot < xtip_l] )
ytip_r = np.median( y_rot[x_rot > xtip_r] )

#View tip demarcations
[plt.axvline(x, linewidth=0.5, color='black') for x in [xtip_l, xtip_r]]
plt.scatter([xtip_l, xtip_r], [ytip_l, ytip_r], marker='o', s=50, color='tab:green')

centre_pt_x = centre_pt_y = 0 #fit point through the centre of mass
if False:
    #Fit through midpoint in rotated space - not reliable
    centre_pt_x = np.ptp(x_rot) / 2 + x_rot.min() #midpoint in rotated space
    centre_pt_y = y_rot[ np.argmin(abs(x_rot - centre_pt_x) ).ravel()[:1000]].mean()
plt.scatter(0, 0, marker='o', s=80, color='tab:blue', edgecolor='white')
plt.scatter(centre_pt_x, centre_pt_y, marker='o', s=50, color='tab:green')

#fit a quadratic curve to 3 points: the left/right tip heights, and the centre pt
from numpy.polynomial import Polynomial

poly = Polynomial.fit(
    [xtip_l, centre_pt_x, xtip_r], [ytip_l, centre_pt_y, ytip_r],
    deg=2
)

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

#View curve
plt.plot(pca0_ax, curve, linewidth=1, color='tab:green')
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()