如何估计3D点的局部切平面？

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中进行？

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

编辑：我应该提到，数据可能包含垂直或接近垂直的表面，所以隐式模型是必要的。 我发现 这个例子在scipy的代码手册中是这样的，但只是xy数据。

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)