matlab 从 3D 数据点拟合椭圆

Question

我在给定平面 (3D) 中有一组 3D (X,Y,Z) 数据点。我希望在这些点上拟合一个椭圆。

我找到了很多关于如何在 2D 点中拟合椭圆的答案。所以更准确地说，我的问题是如何转换 3D 数据(x,y,z) 点 -> 2D 数据(x,y)？

Answer 1

这是我解决这个问题的 Python 代码。 这个链接帮助我完成了我的实现：https://meshlogic.github.io/posts/jupyter/curve-fitting/fitting-a-circle-to-cluster-of-3d-points/

import numpy as np
from skimage.measure import EllipseModel

#-------------------------------------------------------------------------------
# RODRIGUES ROTATION
# - Rotate given points based on a starting and ending vector
# - Axis k and angle of rotation theta given by vectors n0,n1
#   P_rot = P*cos(theta) + (k x P)*sin(theta) + k*<k,P>*(1-cos(theta))
#-------------------------------------------------------------------------------
def rodrigues_rot(P, n0, n1):
    
    # If P is only 1d array (coords of single point), fix it to be matrix
    if P.ndim == 1:
        P = P[np.newaxis,:]
    
    # Get vector of rotation k and angle theta
    n0 = n0/np.linalg.norm(n0)
    n1 = n1/np.linalg.norm(n1)
    k = np.cross(n0,n1)
    k = k/np.linalg.norm(k)
    theta = np.arccos(np.dot(n0,n1))
    
    # Compute rotated points
    P_rot = np.zeros((len(P),3))
    for i in range(len(P)):
        P_rot[i] = P[i]*np.cos(theta) + np.cross(k,P[i])*np.sin(theta) + k*np.dot(k,P[i])*(1-np.cos(theta))

    return P_rot

def fit_an_ellipse(P):
    P_mean = P.mean(axis=0)
    P_centered = P - P_mean
    
    # Fitting plane by SVD for the mean-centered data
    U,s,V = np.linalg.svd(P_centered, full_matrices=False)
    
    # Normal vector of fitting plane is given by 3rd column in V
    # Note linalg.svd returns V^T, so we need to select 3rd row from V^T
    # normal on 3d plane
    normal = V[2,:]
    
    # Project points to coords X-Y in 2D plane
    P_xy = rodrigues_rot(P_centered, normal, [0,0,1])
    
    # Use skimage EllipseModel to fit an ellipse to set of 2d points
    ell = EllipseModel()
    ell.estimate(P_xy[:, :2])
    
    # Generate n 2D points on the fitted elippse
    n = 100
    xy = ell.predict_xy(np.linspace(0, 2 * np.pi, n))
    
    # Convert the 2D generated points to the 3D space
    points = []
    for i in range(len(xy)):
        points.append([xy[i, 0], xy[i, 1], 0])
    points = np.array(points)
    ellipse_points_3d = rodrigues_rot(points, [0,0,1], normal) + P_mean
    
    return ellipse_points_3d

要测试代码，您可以运行它并检查输出结果：

import numpy as np
import plotly
import plotly.graph_objs as go

P = np.array([[52.21818786,  7.86337722, 57.83456389],
       [30.55316226, 32.36591494, 14.35753359],
       [59.77387002, 14.29531811, 53.6462596 ],
       [42.85677086, 32.67223954, -5.95323959],
       [44.46449002,  1.43144171, 54.0253186 ],
       [27.6464027 , 19.80836045, -1.5754063 ],
       [63.48591069,  6.88329618, 57.55556516],
       [44.19484831, 28.32302575,  6.01730042],
       [46.09443886,  2.71782362, 57.98617489],
       [22.55050927, 30.28315605, 42.5642505 ],
       [20.16244533, 18.55944689, 34.06871328],
       [69.4591254 , 33.62256919, 40.91996533],
       [24.42183439,  5.95578526, 35.80224431],
       [70.09161495, 24.03152634, 45.77915268],
       [28.68122335, -6.64788396, 37.53577535],
       [59.84340586, 23.47833222, 60.01530894],
       [23.98376474, 14.23114661, 32.43676647],
       [73.28044481, 29.29426891, 39.28801852],
       [28.48679585, -5.33220296, 36.04206575],
       [54.66351746, 15.7561502 , 51.20981383],
       [38.33444206, -0.08003422, 41.2639318 ],
       [57.27722964, 39.91662965, 20.63778872],
       [43.24856256,  7.79042068, 50.95451935],
       [64.68788661, 31.78841088, 27.19632274],
       [41.67377653, -0.18313508, 49.56081237],
       [60.577958  , 35.8138609 , 28.9005053 ]])
ellipse_points = fit_an_ellipse(P)
lines = []
lines.append(go.Scatter3d(x=P[:, 0], y=P[:, 1], \
                        z=P[:, 2],name='P'\
                          ,opacity=1, ))
lines.append(go.Scatter3d(x=ellipse_points[:, 0], y=ellipse_points[:, 1], \
                        z=ellipse_points[:, 2],name='ellipse_points'\
                          ,opacity=1, ))
plotly.offline.iplot(lines)

输出结果： the 3d fitting result

您可以在 colab 中自己尝试代码：https://colab.research.google.com/drive/1snI2-_S1CY8iUtszRP1bzsFULEQYdJym?usp=sharing

Answer 2

在标准投影中，椭圆（圆是椭圆，a = b = r）将被投影为椭圆或直线。所以我将使用这种行为，这样你想要的 3D 椭圆将由你可以计算的不同椭圆来定义。

我不会显示代码，但方法可以是：

1. 将数据定义为 M x 3 矩阵，其中 M 为点数，格式为[x,y,z]
2. 以z=f(x,y) 的形式定义平面
3. 在数据矩阵中搜索等于或类似于[x,y,f(x,y)] 向量的行
4. 假设结果点是椭圆形间隔的。然后使用答案how to fit ellipse to [x,y] 对在搜索结果的矩阵中（忽略z 部分可以理解为对x-y 平面的投影）。
5. 将拟合数据转换为[x_fit,y_fit]形式的N x 2矩阵
6. 将最后一个矩阵展开为[x_fit,y_fit,f(x_fit,y_fit)]的形式
7. 瞧 - 我们在这里安装了省略号。