Как оценить локальную касательную плоскость для трехмерных точек?

Question

Допустим, у меня есть множество точек,

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

Для каждой точки (р) в R, я определил локальную окрестность с радиусом (г) и высоты (2r) с использованием scipy.cKDTree

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 user guide. Кажется, он обрабатывает трехмерные точки. Любые рекомендации по его внедрению и использованию приветствуются. Я также нашел this similar question.

EDIT: Следует отметить, что данные могут содержать вертикальные или почти вертикальные поверхности, поэтому необходима неявная модель. Я нашел this example in the scipy codebook, но только с данными xy.

Answer 1

Это общий пример того, как установить трехмерную поверхность двумя облаками точек с помощью scipy.odr. Надеюсь, поможет.

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)