update the codebase

src/__init__.py

@@ -0,0 +1,5 @@
+from .data import *
+from .plotting import *
+from .normals import *
+from .orientation import *
+from .metrics import *

src/metrics.py

@@ -1,46 +1,40 @@
 import numpy as np
 
 
-__all__ = ['weight', 'rms_angle_error']
+__all__ = ['rms_angle_error']
 
 
-def weight(p, nbhd, kernel='rbf', sigma=None):
-    """Return scaling weights given distance between a targeted point
-    and its surrounding local neighborhood.
+def _angle_error(n_estim, n_gt, orient=True):
+    """Return the angle error(s) in degrees between the estimated and
+    ground truth normal vector.
     
-    p : numpy_ndarray
-        Targeted point of shape (3, ).
-    nbhd : numpy.ndarray
-        An array of shape (N, 3) representing the local neighborhood.
-    kernel : str, optional
-        The weighting function to use for MLS fitting. Possible options
-        include 'rbf', the Gaussian kernel, 'cosine', the cosine
-        similarity computed as the L2-normalized dot product,
-        'truncated', the truncated quadratic kernel, 'linear', the
-        linear kernel for weighting, 'inverse', the inverse kernel. If
-        not set, all weights will be set to 1.
-    sigma : float, optional
-        A scaling factor for the weighting function. If not given, it
-        is set to 1 / N where N is the total number of points in the
-        local neighborhood.
+    Parameters
+    ----------
+    n_estim : numpy.ndsarray
+        Estimated unit normals of shape (N, 3), where N is the
+        number of points in the point cloud.
+    n_gt : numpy.ndarray
+        Ground truth normals of the corresponding shape.
+    orient : bool, optional
+        If it is set to True, orientation of the normals is taken
+        into account. Otherwise, orientation does not matter.
+    
+    Returns
+    -------
+    float or numpy.ndarray
+        Angle error(s) in degrees.
     """
-    dist = np.linalg.norm(nbhd - p, axis=1)  # squared Euclidian distance
-    if sigma is None:
-        sigma = 1.
-    if kernel == 'rbf':  # Gaussian
-        w = np.exp(-dist ** 2 / (2 * sigma ** 2))
-    elif kernel == 'cosine':
-        w = (nbhd @ p) / np.linalg.norm(nbhd * p, axis=1)
-    elif kernel == 'linear':
-        w = np.maximum(1 - sigma * dist, 0)
-    elif kernel == 'inverse':
-        w = 1 / dist
-    elif kernel == 'truncated':
-        w = np.maximum(1 - sigma * dist ** 2, 0)
+    N = n_gt.shape[0]
+    if orient:
+        dot = np.sum(n_estim * n_gt, axis=1)
     else:
-        w = np.ones_like(dist)
-    return w
-
+        dot = np.sum(np.abs(n_estim) * np.abs(n_gt), axis=1)
+    dot = np.where(dot > 1, 1, dot)
+    dot = np.where(dot < -1, -1, dot)
+    if np.sum(dot) / N == 1:  # all normals matched -> no error
+        return 0
+    return np.arccos(dot) * 180 / np.pi
+
 
 def rms_angle_error(n_estim, n_gt, orient=True):
     """Return the root mean square angle error between estimated and
@@ -62,14 +56,5 @@ def rms_angle_error(n_estim, n_gt, orient=True):
     float
         Root mean square angle error in degrees.
     """
-    N = n_gt.shape[0]
-    if orient:
-        rel = np.sum(n_estim * n_gt, axis=1)
-    else:
-        rel = np.sum(np.abs(n_estim) * np.abs(n_gt), axis=1)
-    rel = np.where(rel > 1, 1, rel)
-    rel = np.where(rel < -1, -1, rel)
-    if np.sum(rel) / N == 1:  # all normals matched -> no error
-        return 0
-    theta = np.arccos(rel) * 180 / np.pi
-    return np.sqrt(np.mean(theta ** 2))
+    err = _angle_error(n_estim, n_gt, orient)
+    return np.sqrt(np.mean(err ** 2))

src/normals.py

@@ -0,0 +1,257 @@
+import numpy as np
+from scipy import interpolate
+from scipy import spatial
+from sklearn.decomposition import KernelPCA
+
+from .utils import apply_weight, polyfit2d
+
+
+__all__ = ['estimate_normals_pca',
+           'estimate_normals_kpca',
+           'estimate_normals_spline',
+           'estimate_normals_poly']
+
+
+def estimate_normals_pca(xyz, k, kernel=None, **kwargs):
+    """Return the unit normals by fitting local tangent plane at each
+    point in the point cloud.
+    
+    Ref: Hoppe et al., in proceedings of SIGGRAPH 1992, pp. 71-78,
+         doi: 10.1145/133994.134011
+    
+    Parameters
+    ----------
+    xyz : numpy.ndarray
+        The point cloud of shape (N, 3), N is the number of points.
+    k : float
+        The number of nearest neighbors of a local neighborhood around
+        a current query point.
+    kernel : string, optional
+        Kernel for computing distance-based weights.
+    kwargs : dict, optional
+        Additional keyword arguments for computing weights. For details
+        see `apply_weight` function.
+    
+    Returns
+    -------
+    numpy.ndarray
+        The unit normals of shape (N, 3), where N is the number of
+        points in the point cloud.
+    """
+    # create a kd-tree for quick nearest-neighbor lookup
+    tree = spatial.KDTree(xyz)
+    n = np.empty_like(xyz)
+    for i, p in enumerate(xyz):
+        # extract the local neighborhood
+        _, idx = tree.query([p], k=k, eps=0.1, workers=-1)
+        nbhd = xyz[idx.flatten()]
+
+        # compute the kernel function and create the weights matrix
+        if kernel:
+            w = apply_weight(p, nbhd, kernel, **kwargs)
+        else:
+            w = np.ones((nbhd.shape[0], ))
+        W = np.diag(w)
+
+        # extract an eigenvector with smallest associeted eigenvalue
+        X = nbhd.copy()
+        X = X - X.mean(axis=0)
+        C = (X.T @ (W @ X)) / (nbhd.shape[0] - 1)
+        U, S, VT = np.linalg.svd(C)
+        n[i, :] =  U[:, np.argmin(S)]
+    return n
+
+
+def estimate_normals_kpca(xyz, k, kernel=None, **kwargs):
+    """Return the unit normals by fitting local tangent plane at each
+    point in the point cloud by using kernel PCA approach to better
+    handle non-linear patterns.
+    
+    Parameters
+    ----------
+    xyz : numpy.ndarray
+        The point cloud of shape (N, 3), N is the number of points.
+    k : float
+        The number of nearest neighbors of a local neighborhood around
+        a current query point.
+    kernel : string, optional
+        Kernel for computing distance-based weights.
+    kwargs : dict, optional
+        Additional keyword arguments for the initialization of the
+        `sklearn.decomposition.KernelPCA` class.
+    
+    Returns
+    -------
+    numpy.ndarray
+        The unit normals of shape (N, 3), where N is the number of
+        points in the point cloud.
+    """
+    # create a kd-tree for quick nearest-neighbor lookup
+    tree = spatial.KDTree(xyz)
+    n = np.empty_like(xyz)
+    for i, p in enumerate(xyz):
+        # extract the local neighborhood
+        _, idx = tree.query([p], k=k, eps=0.1, workers=-1)
+        nbhd = xyz[idx.flatten()]
+
+        # normalize the data
+        X = nbhd.copy()
+        X = X - X.mean(axis=0)
+
+        # create an instance of the kernel PCA class
+        kpca = KernelPCA(n_components=3, kernel=kernel, **kwargs)
+        kpca.fit(X)
+        U = X.T @ kpca.eigenvectors_
+        ni = U[:, np.argmin(kpca.eigenvalues_)]
+        n[i, :] = ni / np.linalg.norm(ni)
+    return n
+
+
+def estimate_normals_spline(xyz,
+                            k,
+                            deg=3,
+                            s=None,
+                            unit=True,
+                            kernel=None, 
+                            **kwargs):
+    """Return the (unit) normals by constructing smooth bivariate
+    B-spline at each point in the point cloud considering its local
+    neighborhood.
+    
+    Parameters
+    ----------
+    xyz : numpy.ndarray
+        The point cloud of shape (N, 3), N is the number of points.
+    k : float
+        The number of nearest neighbors of a local neighborhood around
+        a current query point.
+    deg : float, optional
+        Degrees of the bivariate spline.
+    s : float, optional
+        Positive smoothing factor defined for smooth bivariate spline
+        approximation.
+    unit : float, optional
+        If true, normals are normalized. Otherwise, surface normals are
+        returned.
+    kernel : string, optional
+        Kernel for computing distance-based weights.
+    kwargs : dict, optional
+        Additional keyword arguments for computing weights. For details
+        see `apply_weight` function.
+    
+    Returns
+    -------
+    numpy.ndarray
+        The (unit) normals of shape (N, 3), where N is the number of
+        points in the point cloud.
+    """
+    # create a kd-tree for quick nearest-neighbor lookup
+    tree = spatial.KDTree(xyz)
+    n = np.empty_like(xyz)
+    for i, p in enumerate(xyz):
+        _, idx = tree.query([p], k=k, eps=0.1, workers=-1)
+        nbhd = xyz[idx.flatten()]
+
+        # change the basis of the local neighborhood
+        X = nbhd.copy()
+        X = X - X.mean(axis=0)
+        C = (X.T @ X) / (nbhd.shape[0] - 1)
+        U, _, _ = np.linalg.svd(C)
+        Xt = X @ U
+
+        # compute weights given specific distance function
+        if kernel:
+            w = apply_weight(p, nbhd, kernel, **kwargs)
+        else:
+            w = np.ones((nbhd.shape[0], ))
+
+        # create a smooth B-Spline representation of the "height" function
+        h = interpolate.SmoothBivariateSpline(*Xt.T, w=w, kx=deg, ky=deg, s=s)
+
+        # compute normals as partial derivatives of the "height" function
+        ni = np.array([-h(*Xt[0, :2], dx=1).item(),
+                       -h(*Xt[0, :2], dy=1).item(),
+                       1])
+
+        # convert normal coordinates into the original coordinate frame
+        ni = U @ ni
+
+        # normalize normals by considering the magnitude of each
+        if unit:
+            ni = ni / np.linalg.norm(ni, 2)
+        n[i, :] = ni
+    return n
+
+
+def estimate_normals_poly(xyz,
+                          k,
+                          deg=1,
+                          unit=True,
+                          kernel=None,
+                          **kwargs):
+    """Return the (unit) normals by fitting 2-D polynomial at each
+    point in the point cloud considering its local neighborhood.
+    
+    Parameters
+    ----------
+    xyz : numpy.ndarray
+        The point cloud of shape (N, 3), N is the number of points.
+    k : float
+        The number of nearest neighbors of a local neighborhood around
+        a current query point.
+    deg : float, optional
+        Degrees of the polynomial.
+    unit : float, optional
+        If true, normals are normalized. Otherwise, surface normals are
+        returned.
+    kernel : string, optional
+        Kernel for computing distance-based weights.
+    kwargs : dict, optional
+        Additional keyword arguments for computing weights. For details
+        see `apply_weight` function.
+    
+    Returns
+    -------
+    numpy.ndarray
+        The (unit) normals of shape (N, 3), where N is the number of
+        points in the point cloud.
+    """
+    # create a kd-tree for quick nearest-neighbor lookup
+    n = np.empty_like(xyz)
+    tree = spatial.KDTree(xyz)
+    for i, p in enumerate(xyz):
+        _, idx = tree.query([p], k=k, eps=0.1, workers=-1)
+        nbhd = xyz[idx.flatten()]
+
+        # change the basis of the local neighborhood
+        X = nbhd.copy()
+        X = X - X.mean(axis=0)
+        C = (X.T @ X) / (nbhd.shape[0] - 1)
+        U, _, _ = np.linalg.svd(C)
+        X_t = X @ U
+
+        # compute weights given specific distance function
+        if kernel:
+            w = apply_weight(p, nbhd, kernel, **kwargs)
+        else:
+            w = np.ones((nbhd.shape[0], ))
+
+        # fit parametric surface by usign a (weighted) 2-D polynomial
+        X_t_w = X_t * w[:, np.newaxis]
+        c = polyfit2d(*X_t_w.T, deg=deg)
+
+        # compute normals as partial derivatives of the "height" function
+        cu = np.polynomial.polynomial.polyder(c, axis=0)
+        cv = np.polynomial.polynomial.polyder(c, axis=1)
+        ni = np.array([-np.polynomial.polynomial.polyval2d(*X_t_w[0, :2], cu),
+                       -np.polynomial.polynomial.polyval2d(*X_t_w[0, :2], cv),
+                       1])
+
+        # convert normal coordinates into the original coordinate frame
+        ni = U @ ni
+
+        # normalize normals by considering the magnitude of each
+        if unit:
+            ni = ni / np.linalg.norm(ni, 2)
+        n[i, :] = ni
+    return n