在 3D 数据上拟合圆的方法（Python）

Question

我想做的是在 3D 点云数据的所有横截面上拟合一个类似于管道的圆。

数据类似于圆柱形，中间有不连续的部分，此时可以忽略。这些数据由 X、Y 和 Z 分量组成。

起初我认为最好采用数据的横截面，找到最适合它的圆，然后将其应用于整个数据以获得许多圆。

我首先通过交换X轴和Z轴来获得圆拟合，使数据变得明确（感谢jlandercy回答我之前的问题（用弧形图描述数据点）。这就是这个想法来自这样一个事实：一个圆应该能够适合水平线和垂直线之间 90 度的位置，如下图所示：。

所以我所做的就是沿着横截面数据的曲线找到斜率为 0 的点，并找到另一个斜率为无穷大的点。但由于曲线上没有斜率无限大的点，我假设完美的圆不能沿着整条曲线拟合，但最好的拟合圆可以画在斜率为0的点上。所以我发现斜率为0的点的坐标，并与横截面数据的最低点进行比较。我使用它们的坐标来找到圆的中心，并使用斜率 0 的点与圆心之间的 X 分量之差来找到半径。结果，我能够画出一个看起来最适合横截面曲线的圆，如下图黄色所示。

当时我忽略了Y分量，因为我发现与X和Z分量相比，方差和标准差非常低，这意味着它们并没有真正影响我的分析。

现在我的计划是尝试将其应用于整个数据，以便我可以为每个横截面创建圆圈。有点像下面这样。

最终的目标是连接数据每个横截面上这些圆的中心，获得一条非直线，我相信一旦我能够获得所有横截面上的圆，这不会太困难.

我希望我想要实现的目标是明确的。此时我只是不知道如何在原始数据上应用我的横截面方法。我正在考虑采用 for 循环将其应用于整个数据，但我无法想象如何为其指定循环。我试图获得 Chat-GPT 的帮助，但它并没有真正帮助我。

我还尝试了其他寻找点云中心线的方法，例如另一篇文章Fit Curve-Spline to 3D Point Cloud所建议的Delanay Triangulation和Voronoi Tessellation的组合。但我没能成功。

如果有人能给我一些提示或建议，我将不胜感激。

原始数据太大，所以我在 OneDrive 上链接：原始数据

这是我的代码片段，展示了我如何实现横截面圆拟合的想法。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Load the point cloud data
data = pd.read_csv('data/selected_2.csv')
points = data.to_numpy()
x = points[:, 0]
y = points[:, 1]
z = points[:, 2]

# Calculate polynomial
a = np.polyfit(z, x, 7)
f_poly = np.poly1d(a)

# Find new x and y data along the curve
z_new = np.linspace(z[0], z[-1], 50)
x_new = f_poly(z_new)
data_new = np.vstack((z_new, x_new)).T

# Define the derivative function
def FindingDerivatives(data_new):
    slopes = []
    for i in range(len(data_new) - 1):
        z1, x1 = data_new[i]
        z2, x2 = data_new[i + 1]
        slope = (x2 - x1) / (z2 - z1)
        slopes.append((z1, slope))
    return slopes

# Define the function to find perpendicular line slopes
def FindingPerpendicular(slopes):
    perpendicular = []
    for z1, slope in slopes:
        if slope != 0:
            inverse = -1 / slope # Perpendicular line is the inverse of the original slope
        else:
            inverse = float('inf')  # Perpendicular to a horizontal line is a vertical line
        perpendicular.append((z1, inverse))
    return perpendicular

# Find the maximum slope element
def FindMinMaxSlopeElement(slopes):
    target_slope_elements = []
    max_slope_elements = 0
    min_slope_elements = 0
    max_slope = max(slopes)[1]
    min_slope = min(slopes)[1]
    print(max_slope, min_slope)
    for i in range(len(slopes)):
        if slopes[i][1] == min_slope:
            max_slope_elements = slopes[i]
        if slopes[i][1] == max_slope:
            min_slope_elements = slopes[i]
    target_slope_elements.append(min_slope_elements)
    target_slope_elements.append(max_slope_elements)
    return target_slope_elements

# Find the corresponding coordinates to the maximum slope element
def FindCorrespondingCoordinates(target_slope_elements, data):
    min_coordinates = []
    max_coordinates = []
    target_coordinates = []
    for i in range(len(data)):
        if data[:, 0][i] == target_slope_elements[0][0]:
            min_coordinates.append(data[:, 0][i])
            min_coordinates.append(data[:, 1][i])
        if data[:, 0][i] == target_slope_elements[1][0]:
            max_coordinates.append(data[:, 0][i])
            max_coordinates.append(data[:, 1][i])
    print(min_coordinates)
    print(max_coordinates)
    target_coordinates.append(min_coordinates)
    target_coordinates.append(max_coordinates)
    return target_coordinates

# Find the circle based on the coordinate values corresponding to the maximum and minimum slopes
def FindCircle(target_coordinates):
    circle_info = []
    center_z = target_coordinates[0][1] 
    center_x = target_coordinates[1][0]
    center = [center_z, center_x]
    radius = target_coordinates[0][0] - target_coordinates[1][0]
    print(center, radius)
    circle_info.append(radius)
    circle_info.append(center)
    return circle_info

# Draw the circle
def CircleDraw(circle, test_data, z_new, x_new):
    radius = circle[0]
    fig, ax = plt.subplots()
    ax.scatter(test_data[:, 1], test_data[:, 0])
    ax.plot(z_new, x_new, 'y')
    cir = plt.Circle((circle[1][0], circle[1][1]), radius, color='r',fill=False)
    center_label = f"center: \n {round(circle[1][0], 2), round(circle[1][1], 2)}"
    plt.annotate(center_label, xy=(circle[1][0], circle[1][1]), ha='center') 
    ax.set_aspect('equal', adjustable='datalim')
    ax.add_patch(cir)
    ax.set_xlabel('Z Data')  # X axis data label
    ax.set_ylabel('X Data')  # Y axis data label
    plt.show()


# Organize the test data (only taking 'X' and 'Z' coordinates of the original data)
test_data = points[:, [0, 2]]
# Find the slope of the original data
slopes = FindingDerivatives(test_data)
# Find perpendicular lines of the original data
perpendicular = FindingPerpendicular(slopes)
# Find minimum and maximum slope elements 
target_slope_elements = FindMinMaxSlopeElement(slopes) 
print("Target Slopes", target_slope_elements)
# Find the target coordinates
target_coordinates = FindCorrespondingCoordinates(target_slope_elements, test_data) 
print("Target Coordinates: ", target_coordinates)
# Find the circle information
circle = FindCircle(target_coordinates)
print("radius: %s center: %s" % (circle[0], circle[1]))
# Call the circle function
CircleDraw(circle, test_data, z_new, x_new)

# Here is the cross-section of the data
# ('Y' component is not really used as I decided that they don't really affect my analysis)
X,Y,Z
1702113.3399999738,2434223.5100000203,9080.16
1702132.699999988,2434227.2700000107,9079.86
1702172.7199999988,2434233.629999995,9106.22
1702153.0400000215,2434230.4699999997,9087.98
1702143.6899999976,2434227.5699999933,9083.15
1702167.079999983,2434224.6899999976,9100.11
1702159.8799999952,2434236.3600000143,9093.52
1702177.9900000098,2434225.169999987,9113.97
1702192.0399999917,2434220.0899999742,9155.37
1702186.7199999988,2434228.7599999905,9128.54
1702191.330000013,2434221.2299999893,9174.77
1702182.9099999962,2434234.8899999857,9120.3
1702182.580000013,2434217.7899999917,9121.94
1702192.1200000048,2434219.9599999785,9164.910000000002
1702189.5399999917,2434224.2199999997,9137.12
1702191.229999989,2434221.3899999857,9146.08
1702188.9300000072,2434225.100000024,9185.21
1702183.5,2434216.300000012,9198.22
1702175.5,2434229.150000006,9220.28
1702178.099999994,2434224.9799999893,9209.75
1702188.080000013,2434226.4600000083,9223.87
1702184.669999987,2434231.9600000083,9196.33
1702214.9900000098,2434218.3599999845,9220.34

如果我的问题不够清楚，请评论。我将根据需要更新和修改问题。

非常感谢！

Answer 1

您没有指定编程语言，因此如果 C++ 适合您，那么 PCL 可能就是您正在寻找的语言。 PCL 具有 RANSAC 方法的实现，用于在噪声数据中查找几何模型。支持许多不同的模型，您最感兴趣的是 CIRCLE2D、CIRLCE3D 和 CYLINDER（请参阅 https://pointclouds.org/documentation/group__sample__consensus.html）

这里是有关如何安装圆柱体的教程：https://pcl.readthedocs.io/projects/tutorials/en/master/Columns_segmentation.html