Most efficient way to calculate point wise surface normal from a numpy grid-CodePudding

say we have a 2D grid that is projected on a 3D surface, resulting in a 3D numpy array, like the below image. What is the most efficient way to calculate a surface normal for each point of this grid?

CodePudding user response：

I can give you an example with simulated data:

I showed your way, with three points. With three points you can always calculate the cross product to get the perpendicular vector based on the two vectors created from three points. Order does not matter.

I took the liberty to also add the PCA approach using predefined sklearn functions. You can create your own PCA, good exercise to understand what happens under the hood but this works fine. The benefit of the approach is that it is easy to increase the number of neighbors and you are still able to calculate the normal vector. It is also possible to select the neighbors within a range instead of N nearest neighbors.

If you need more explanation about the working of the code please let me know.

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.05)

# Flatten all arrays
X, Y, Z = map(lambda x: x.flatten(), [X, Y, Z])
plt.plot(X, Y, Z, '.')
plt.show()

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))

# map with functools
pca = PCA(n_components=3)
normals3 = list(map(partial(PCA_unit_vector, pca=pca), combinations))

CodePudding user response：

The surface normal for a point cloud is not well defined. One way to define them is from the surface normal of a reconstructed mesh using triangulation (which can introduce artefacts regarding you specific input). A relatively simple and fast solution is to use VTK to do that, and more specifically, vtkSurfaceReconstructionFilter and vtkPolyDataNormals . Regarding your needs, it might be useful to apply other filters.