从NumPy网格计算逐点表面法线的最有效方法

Question

从NumPy网格计算逐点表面法线的最有效方法

pythonnumpycomputational-geometrysurface

3

假设我们有一个二维网格，投影在三维表面上，得到一个三维numpy数组，如下图所示。最有效的方法是什么，以计算此网格每个点的表面法线？

- azerila

1

如何定义表面法线？我可以基于三个点或N个最近邻计算法线。对于边缘情况，你想怎么处理？也许更重要的是，你尝试了什么，哪些方法对你不起作用？ - 3dSpatialUser

@3DspatialUser 很好的问题，对我来说3个点已经足够了。但是如果不会增加太多计算量的话，有选择n个点的灵活性会很有用。 - azerila

仍然是哪三个点？最近的三个点？还是包括该点本身的三个点？我建议您结合最近邻和主成分分析来查看 https://stats.stackexchange.com/questions/2691/making-sense-of-principal-component-analysis-eigenvectors-eigenvalues。 - 3dSpatialUser

@3DspatialUser 包括该点本身在内共有3个点。PCA的复杂度可能很高，不是吗？我想象一下，如果我们可以为每个点找到2个最近的点，并且计算量非常低，我们就可以更快地完成它。 - azerila

1

对于所有点，使用kD树找到所有点的第一个最近邻（比如8个），并拟合最小二乘平面。如果它们的距离看起来异常，请考虑丢弃这些点。请注意，方向将是不确定的。 - user1196549

2个回答

1

点云的表面法线并不是很明确。一种定义方法是使用三角化重建网格的表面法线（这可能会引入与您特定输入相关的伪影）。一个相对简单和快速的解决方案是使用VTK，更具体地说，使用vtkSurfaceReconstructionFilter和vtkPolyDataNormals 。根据您的需求，应用其他过滤器可能会有所帮助。

- Jérôme Richard

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- 3dSpatialUser · Accepted Answer

我可以用模拟数据给你举个例子：

我展示了一种方法，使用三个点。使用三个点，您始终可以计算叉积，以基于从三个点创建的两个向量获得垂直向量。顺序无关紧要。

我还自作主张添加了使用预定义的sklearn函数的PCA方法。您可以创建自己的PCA，这是理解底层发生的事情的好练习，但这也很好用。该方法的好处是很容易增加邻居的数量，而且仍然能够计算法向量。还可以选择在范围内选择邻居，而不是N个最近的邻居。

如果您需要更多有关代码工作原理的解释，请告诉我。

from functools import partial
import numpy as np
from sklearn.neighbors import KDTree

from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt

from sklearn.decomposition import PCA


fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# Grab some test data.
X, Y, Z = axes3d.get_test_data(0.25)

X, Y, Z = map(lambda x: x.flatten(), [X, Y, Z])

plt.plot(X, Y, Z, '.')
plt.show(block=False)

data = np.array([X, Y, Z]).T

tree = KDTree(data, metric='minkowski') # minkowki is p2 (euclidean)

# Get indices and distances:
dist, ind = tree.query(data, k=3) #k=3 points including itself

def calc_cross(p1, p2, p3):
    v1 = p2 - p1
    v2 = p3 - p1
    v3 =  np.cross(v1, v2)
    return v3 / np.linalg.norm(v3)

def PCA_unit_vector(array, pca=PCA(n_components=3)):
    pca.fit(array)
    eigenvalues = pca.explained_variance_
    return pca.components_[ np.argmin(eigenvalues) ]

combinations = data[ind]

normals = list(map(lambda x: calc_cross(*x), combinations))

# lazy with map
normals2 = list(map(PCA_unit_vector, combinations))


## NEW ##

def calc_angle_with_xy(vectors):
    '''
    Assuming unit vectors!
    '''
    l = np.sum(vectors[:,:2]**2, axis=1) ** 0.5
    return np.arctan2(vectors[:, 2], l)

    

dist, ind = tree.query(data, k=5) #k=3 points including itself
combinations = data[ind]
# map with functools
pca = PCA(n_components=3)
normals3 = list(map(partial(PCA_unit_vector, pca=pca), combinations))

print( combinations[10] )
print(normals3[10])


n = np.array(normals3)
n[calc_angle_with_xy(n) < 0] *= -1

def set_axes_equal(ax):
    '''Make axes of 3D plot have equal scale so that spheres appear as spheres,
    cubes as cubes, etc..  This is one possible solution to Matplotlib's
    ax.set_aspect('equal') and ax.axis('equal') not working for 3D.

    Input
      ax: a matplotlib axis, e.g., as output from plt.gca().
    
    FROM: https://dev59.com/oGYr5IYBdhLWcg3wb5rc
    '''

    x_limits = ax.get_xlim3d()
    y_limits = ax.get_ylim3d()
    z_limits = ax.get_zlim3d()

    x_range = abs(x_limits[1] - x_limits[0])
    x_middle = np.mean(x_limits)
    y_range = abs(y_limits[1] - y_limits[0])
    y_middle = np.mean(y_limits)
    z_range = abs(z_limits[1] - z_limits[0])
    z_middle = np.mean(z_limits)

    # The plot bounding box is a sphere in the sense of the infinity
    # norm, hence I call half the max range the plot radius.
    plot_radius = 0.5*max([x_range, y_range, z_range])

    ax.set_xlim3d([x_middle - plot_radius, x_middle + plot_radius])
    ax.set_ylim3d([y_middle - plot_radius, y_middle + plot_radius])
    ax.set_zlim3d([z_middle - plot_radius, z_middle + plot_radius])


u, v, w = n.T

fig = plt.figure()
ax = fig.add_subplot(projection='3d')
# ax.set_aspect('equal')

# Make the grid
ax.quiver(X, Y, Z, u, v, w, length=10, normalize=True)
set_axes_equal(ax)
plt.show()