Novel Approach to Normal Estimation in Surface Reconstruction

Question

Not sure where to post this exactly but I am hoping to find out if my approach is in any way a novel method for estimating normals in a point cloud when performing surface reconstruction (in e.g. ball-pivoting algorithm) or if something already exists in literature that is superior.

Problem

It is my understanding that most normal estimation techniques use normal averaging by creating a triangulation (which form planes) with neighboring points and then simply averaging the normals from the planes formed in the local triangulation. The issue with this approach is that the winding order in the cross products, while maybe consistent, might result in sporadic inversions of the normals.

Solution

In order to avoid this limitation, my thought was the following: if curvature exists on any sample of points in the point cloud then precisely one point will protrude farther than the rest. If the projection of this "protruding" point onto the plane formed by the three nearest points lies within the triangle formed by these points (on the plane they create), then the normal of the "protruding" point is equal to the vector formed from the projection (a point on the plane) toward the point itself. Of course, not all points are guaranteed to have a "successful" projection (one that falls within a triangle). In order to address this, one can loop through neighboring points that might be slightly farther from the protruding point until a projection is successful.

In the event that there is no successful projection, one can identify those points which did not possess an estimation of a normal and then average as one might normally do but using the normals that already exist for neighboring points.

Implementation

In practice, this looks like the following (in Python):

from scipy.spatial import KDTree
from scipy.linalg import solve
from itertools import combinations

points = <ndarray of shape (n, 3)>
kd_tree = KDTree(points)

lookaround = 15
dists, nns = kd_tree.query(points_new, k=lookaround)

normals = []
count = 0

combos = sorted(list(combinations(range(1,lookaround), r=3)), key=lambda E: sum(E))
for group in points[nns]:

    print(count / len(nns) * 100, end='\r')

    # Find three closest points such that the projection of subject point onto the formed plane is within the triangle formed by the same points that create the plane
    # https://stackoverflow.com/questions/25512037/how-to-determine-if-a-point-lies-over-a-triangle-in-3d
    #     Let N = (B-A) X (C-A), N1 = (B-A) X (P-A), N2 = (C-B) X (P-B), N3 = (A-C) X (P-C)
    #     return N1 * N >= 0 and N2 * N >= 0 and N3 * N >= 0; (point lies over triangle)


    try:
        overlap = False
        c = 0

        while not overlap:
            i,j,k = combos[c]  # Preferences proximity somewhat 
            P = group[0, :3]
            A = group[i, :3]
            B = group[j, :3]
            C = group[k, :3]

            N  = np.cross(B-A, C-A)
            N1 = np.cross(B-A, P-A)
            N2 = np.cross(C-B, P-B)
            N3 = np.cross(A-C, P-C)

            if np.dot(N, N1) >= 0 and np.dot(N, N2) >= 0 and np.dot(N, N3) >= 0:
                overlap = True
            else:
                c+=1

        n = np.cross(B - A, C - A)
        c = np.dot(A, n)
        n_x, n_y, n_z = n
        x_0, y_0, z_0 = P
        # This array is created based on Lagrangian multipliers (see https://math.stackexchange.com/questions/723937/find-the-point-on-a-plane-3x-4y-z-1-that-is-closest-to-1-0-1 and https://people.duke.edu/~ccc14/sta-663-2018/notebooks/S09H_Constrained_Optimization.html)
        M = np.array([[  2,   0,   0, -n_x],
                      [  0,   2,   0, -n_y],
                      [  0,   0,   2, -n_z],
                      [n_x, n_y, n_z,    0]])
        b = np.array([2*x_0, 2*y_0, 2*z_0, c])
        x, y, z, _ = solve(M, b)
        # Want to point away from the plane toward the point
        dist = np.array([x_0 - x, y_0 - y, z_0 - z])
        normals.append(dist / np.linalg.norm(dist))

    except IndexError: # No neighbors produced a good projection
        normals.append(np.full(3, np.nan))

    count += 1

normals = np.array(normals)

# Now loop back through and try to correct the normals
nan_idxs = np.unique(np.where(np.isnan(normals))[0])
while len(nan_idxs) > 0:
    # Now conduct the second pass looking for np.nan normals
    nan_idxs = np.unique(np.where(np.isnan(normals))[0])
    print(len(nan_idxs))
    # Get all the these points from points_new
    nan_points = points_new[nan_idxs]
    tmp_tree = KDTree(points_new)
    # Get 4 nearest neaighbors from these points
    _, good_nns = tmp_tree.query(nan_points, k=(pass_+1)*5)

    for i, idx in enumerate(nan_idxs):
        normals[idx, 0] = np.nanmean(normals[good_nns][i, :,0])
        normals[idx, 1] = np.nanmean(normals[good_nns][i, :,1])
        normals[idx, 2] = np.nanmean(normals[good_nns][i, :,2])