如何对三维点估计局部切平面？

Question

假设我有一组点，

R = [[x1, y1, z1],[x2, y2, z2],...,[xn, yn, zn]]

对于R中的每个点（p），我使用scipy.cKDTree识别了一个局部邻域，其半径为（r），高度为（2r）。

import numpy as np
import scipy.spatial

R = np.array(R)
r = 1 # search radius

xy = R[:,0:2] # make array of ONLY xy
tree = scipy.spatial.cKDTree(xy)

for p in range(len(R)):
    2d_nbr_indices = tree.query_ball_point(xy[p], r) # indices w/in xy neighborhood
    2d_nbr_array = R[2d_nbr_indices] # 3d array of 2d neighbors
    z = R[p][1] # get z value
    zMin = z - r
    zMax = z + r
    # Create boolean array to filter 3d array
    hgt_filter = np.any([2d_nbr_array[:, 2] >= zMin,
                         2d_nbr_array[:, 2] <= zMax], axis=0)
    3d_nbr_array = 2d_nbr_array[hgt_filter] # points in xyz neighborhood

我想为每个邻域计算正交回归平面，确定每个点到平面的距离（正交距离），并计算平面的法向量。有没有人可以建议如何在Python中实现这一点？

编辑：我发现了一个 odr用户指南。它似乎可以处理3D点云。欢迎任何关于它的实现和使用的建议。我还找到了 这个类似的问题。

编辑：我应该提到数据可能包含垂直或接近垂直的表面，因此需要一个隐式模型。我在Scipy代码库中找到了 这个xy数据示例，但只有xy数据。

- Barbarossa

Answer 1

这是一个使用scipy.odr拟合点云的3D表面的通用示例。希望能对你有所帮助。

from scipy.odr import ODR, Model, Data
import numpy as np
import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

def func(beta,data):
    x,y = data
    a,b,c = beta
    return a*x+b*y+c

N = 20
x = np.random.randn(N)
y = np.random.randn(N)
z = func([-3,-1,2],[x,y])+np.random.normal(size=N)


data = Data([x,y],z)
model = Model(func)
odr = ODR(data, model, [0.0,0.0,0.0])
odr.set_job(fit_type = 0)
res = odr.run()

Y,X = np.mgrid[y.min():y.max():20j,x.min():x.max():20j]
Z = func(res.beta, [X,Y])
f = plt.figure()
pl = f.add_subplot(111,projection='3d')
pl.scatter3D(x,y,z)
pl.plot_surface(X,Y,Z,alpha=0.4)

- Christian K.