Normal orientation issues in complex geometries (#6198)

Given an input point cloud, the algorithm processes the data primarily based on Euclidean distance, using the K-nearest neighbors (KNN) algorithm as a basis. In the case of specific geometries, this can be problematic when orienting two distinct surfaces of the same volume that should have opposite orientations (e.g., cusps, stenosis, etc.).

Therefore, we provide the user with the ability to manipulate two parameters in the OrientNormalsConsistentTangentPlane function (or orient_normals_consisten_tangent_plane in Python), namely lambda and cos_alpha_tol, which respectively allow for a penalty and a threshold on the distance of a point x from the plane defined by a point z and its estimated normal
(calculated in the EstimateNormal method).

Report: https://eugeniovaretti.github.io/meshreco/Piazza_Valentini_Varetti_MeshReconstructionFromPointCloud_2023.pdf

---------

Co-authored-by: Benjamin Ummenhofer <[email protected]>

cpp/open3d/geometry/EstimateNormals.cpp

@@ -359,7 +359,10 @@ void PointCloud::OrientNormalsTowardsCameraLocation(
     }
 }
 
-void PointCloud::OrientNormalsConsistentTangentPlane(size_t k) {
+void PointCloud::OrientNormalsConsistentTangentPlane(
+        size_t k,
+        const double lambda /* = 0.0*/,
+        const double cos_alpha_tol /* = 1.0*/) {
     if (!HasNormals()) {
         utility::LogError(
                 "No normals in the PointCloud. Call EstimateNormals() first.");
@@ -380,8 +383,20 @@ void PointCloud::OrientNormalsConsistentTangentPlane(size_t k) {
         v1 = pt_map[v1];
         size_t edge = EdgeIndex(v0, v1);
         if (graph_edges.count(edge) == 0) {
-            double dist = (points_[v0] - points_[v1]).squaredNorm();
-            delaunay_graph.push_back(WeightedEdge(v0, v1, dist));
+            const auto diff = points_[v0] - points_[v1];
+            // penalization on normal-plane distance
+            double dist = diff.squaredNorm();
+            double penalization = lambda * std::abs(diff.dot(normals_[v0]));
+
+            // if cos_alpha_tol < 1 some edges will be excluded. In particular
+            // the ones connecting points that form an angle below a certain
+            // threshold (defined by the cosine)
+            double cos_alpha =
+                    std::abs(diff.dot(normals_[v0])) / std::sqrt(dist);
+            if (cos_alpha > cos_alpha_tol)
+                dist = std::numeric_limits<double>::infinity();
+
+            delaunay_graph.push_back(WeightedEdge(v0, v1, dist + penalization));
             graph_edges.insert(edge);
         }
     };
@@ -403,19 +418,76 @@ void PointCloud::OrientNormalsConsistentTangentPlane(size_t k) {
         edge.weight_ = NormalWeight(edge.v0_, edge.v1_);
     }
 
+    // The function below takes v0 and its neighbors as inputs.
+    // The function returns the quartiles of the distances between the neighbors
+    // and a plane defined by the normal vector of v0 and the point v0.
+    auto compute_q1q3 =
+            [&](size_t v0,
+                std::vector<int> neighbors) -> std::array<double, 2> {
+        std::vector<double> dist_plane;
+
+        for (size_t vidx1 = 0; vidx1 < neighbors.size(); ++vidx1) {
+            size_t v1 = size_t(neighbors[vidx1]);
+            const auto diff = points_[v0] - points_[v1];
+            double dist = std::abs(diff.dot(normals_[v0]));
+            dist_plane.push_back(dist);
+        }
+        std::vector<double> dist_plane_ord = dist_plane;
+        std::sort(dist_plane_ord.begin(), dist_plane_ord.end());
+        // calculate quartiles
+        int q1_idx = static_cast<int>(dist_plane_ord.size() * 0.25);
+        double q1 = dist_plane_ord[q1_idx];
+
+        int q3_idx = static_cast<int>(dist_plane_ord.size() * 0.75);
+        double q3 = dist_plane_ord[q3_idx];
+
+        std::array<double, 2> q1q3;
+        q1q3[0] = q1;
+        q1q3[1] = q3;
+
+        return q1q3;
+    };
+
     // Add k nearest neighbors to Riemannian graph
     KDTreeFlann kdtree(*this);
     for (size_t v0 = 0; v0 < points_.size(); ++v0) {
         std::vector<int> neighbors;
         std::vector<double> dists2;
+
         kdtree.SearchKNN(points_[v0], int(k), neighbors, dists2);
+
+        const double DEFAULT_VALUE = std::numeric_limits<double>::quiet_NaN();
+        std::array<double, 2> q1q3 =
+                lambda == 0
+                        ? std::array<double, 2>{DEFAULT_VALUE, DEFAULT_VALUE}
+                        : compute_q1q3(v0, neighbors);
+
         for (size_t vidx1 = 0; vidx1 < neighbors.size(); ++vidx1) {
             size_t v1 = size_t(neighbors[vidx1]);
             if (v0 == v1) {
                 continue;
             }
             size_t edge = EdgeIndex(v0, v1);
             if (graph_edges.count(edge) == 0) {
+                const auto diff = points_[v0] - points_[v1];
+                double normal_dist = std::abs(diff.dot(normals_[v0]));
+
+                double dist = diff.squaredNorm();
+
+                // if cos_alpha_tol < 1 some edges will be excluded. In
+                // particular the ones connecting points that form an angle
+                // below a certain threshold (defined by the cosine)
+                double cos_alpha =
+                        std::abs(diff.dot(normals_[v0])) / std::sqrt(dist);
+                if (cos_alpha > cos_alpha_tol) continue;
+
+                // if we are in a penalizing framework do not consider outliers
+                // in terms of distance from the plane (if any)
+                if (lambda != 0) {
+                    double iqr = q1q3[1] - q1q3[0];
+                    if (normal_dist > q1q3[1] + 1.5 * iqr) continue;
+                }
+
                 double weight = NormalWeight(v0, v1);
                 mst.push_back(WeightedEdge(v0, v1, weight));
                 graph_edges.insert(edge);
@@ -436,13 +508,14 @@ void PointCloud::OrientNormalsConsistentTangentPlane(size_t k) {
     }
 
     // find start node for tree traversal
-    // init with node that maximizes z
-    double max_z = std::numeric_limits<double>::lowest();
+    // init with node that minimizes z
+
+    double min_z = std::numeric_limits<double>::max();
     size_t v0 = 0;
     for (size_t vidx = 0; vidx < points_.size(); ++vidx) {
         const Eigen::Vector3d &v = points_[vidx];
-        if (v(2) > max_z) {
-            max_z = v(2);
+        if (v(2) < min_z) {
+            min_z = v(2);
             v0 = vidx;
         }
     }
@@ -457,7 +530,7 @@ void PointCloud::OrientNormalsConsistentTangentPlane(size_t k) {
             n1 *= -1;
         }
     };
-    TestAndOrientNormal(Eigen::Vector3d(0, 0, 1), normals_[v0]);
+    TestAndOrientNormal(Eigen::Vector3d(0, 0, -1), normals_[v0]);
     while (!traversal_queue.empty()) {
         v0 = traversal_queue.front();
         traversal_queue.pop();

cpp/open3d/geometry/PointCloud.h

@@ -248,10 +248,18 @@ class PointCloud : public Geometry3D {
     /// \brief Function to consistently orient estimated normals based on
     /// consistent tangent planes as described in Hoppe et al., "Surface
     /// Reconstruction from Unorganized Points", 1992.
+    /// Further details on parameters are described in
+    /// Piazza, Valentini, Varetti, "Mesh Reconstruction from Point Cloud",
+    /// 2023.
     ///
     /// \param k k nearest neighbour for graph reconstruction for normal
     /// propagation.
-    void OrientNormalsConsistentTangentPlane(size_t k);
+    /// \param lambda penalty constant on the distance of a point from the
+    /// tangent plane \param cos_alpha_tol treshold that defines the amplitude
+    /// of the cone spanned by the reference normal
+    void OrientNormalsConsistentTangentPlane(size_t k,
+                                             const double lambda = 0.0,
+                                             const double cos_alpha_tol = 1.0);
 
     /// \brief Function to compute the point to point distances between point
     /// clouds.

cpp/open3d/t/geometry/PointCloud.cpp

@@ -731,12 +731,15 @@ void PointCloud::OrientNormalsTowardsCameraLocation(
     }
 }
 
-void PointCloud::OrientNormalsConsistentTangentPlane(size_t k) {
+void PointCloud::OrientNormalsConsistentTangentPlane(
+        size_t k,
+        const double lambda /* = 0.0*/,
+        const double cos_alpha_tol /* = 1.0*/) {
     PointCloud tpcd(GetPointPositions());
     tpcd.SetPointNormals(GetPointNormals());
 
     open3d::geometry::PointCloud lpcd = tpcd.ToLegacy();
-    lpcd.OrientNormalsConsistentTangentPlane(k);
+    lpcd.OrientNormalsConsistentTangentPlane(k, lambda, cos_alpha_tol);
 
     SetPointNormals(core::eigen_converter::EigenVector3dVectorToTensor(
             lpcd.normals_, GetPointPositions().GetDtype(), GetDevice()));

cpp/open3d/t/geometry/PointCloud.h

@@ -517,10 +517,18 @@ class PointCloud : public Geometry, public DrawableGeometry {
     /// \brief Function to consistently orient estimated normals based on
     /// consistent tangent planes as described in Hoppe et al., "Surface
     /// Reconstruction from Unorganized Points", 1992.
+    /// Further details on parameters are described in
+    /// Piazza, Valentini, Varetti, "Mesh Reconstruction from Point Cloud",
+    /// 2023.
     ///
     /// \param k k nearest neighbour for graph reconstruction for normal
     /// propagation.
-    void OrientNormalsConsistentTangentPlane(size_t k);
+    /// \param lambda penalty constant on the distance of a point from the
+    /// tangent plane \param cos_alpha_tol treshold that defines the amplitude
+    /// of the cone spanned by the reference normal
+    void OrientNormalsConsistentTangentPlane(size_t k,
+                                             const double lambda = 0.0,
+                                             const double cos_alpha_tol = 1.0);
 
     /// \brief Function to compute point color gradients. If radius is provided,
     /// then HybridSearch is used, otherwise KNN-Search is used.

cpp/pybind/geometry/pointcloud.cpp

@@ -139,7 +139,7 @@ void pybind_pointcloud(py::module &m) {
                  &PointCloud::OrientNormalsConsistentTangentPlane,
                  "Function to orient the normals with respect to consistent "
                  "tangent planes",
-                 "k"_a)
+                 "k"_a, "lambda"_a = 0.0, "cos_alpha_tol"_a = 1.0)
             .def("compute_point_cloud_distance",
                  &PointCloud::ComputePointCloudDistance,
                  "For each point in the source point cloud, compute the "

cpp/pybind/t/geometry/pointcloud.cpp

@@ -320,11 +320,70 @@ infinite value. It also removes the corresponding attributes.
                    "Function to orient the normals of a point cloud.",
                    "camera_location"_a = core::Tensor::Zeros(
                            {3}, core::Float32, core::Device("CPU:0")));
-    pointcloud.def("orient_normals_consistent_tangent_plane",
-                   &PointCloud::OrientNormalsConsistentTangentPlane,
-                   "Function to orient the normals with respect to consistent "
-                   "tangent planes.",
-                   "k"_a);
+    pointcloud.def(
+            "orient_normals_consistent_tangent_plane",
+            &PointCloud::OrientNormalsConsistentTangentPlane, "k"_a,
+            "lambda"_a = 0.0, "cos_alpha_tol"_a = 1.0,
+            R"(Function to consistently orient the normals of a point cloud based on tangent planes.
+
+The algorithm is described in Hoppe et al., "Surface Reconstruction from Unorganized Points", 1992. 
+Additional information about the choice of lambda and cos_alpha_tol for complex
+point clouds can be found in Piazza, Valentini, Varetti, "Mesh Reconstruction from Point Cloud", 2023 
+(https://eugeniovaretti.github.io/meshreco/Piazza_Valentini_Varetti_MeshReconstructionFromPointCloud_2023.pdf).
+
+Args:
+    k (int): Number of neighbors to use for tangent plane estimation.
+    lambda (float): A non-negative real parameter that influences the distance 
+        metric used to identify the true neighbors of a point in complex 
+        geometries. It penalizes the distance between a point and the tangent 
+        plane defined by the reference point and its normal vector, helping to 
+        mitigate misclassification issues encountered with traditional 
+        Euclidean distance metrics.
+    cos_alpha_tol (float): Cosine threshold angle used to determine the 
+        inclusion boundary of neighbors based on the direction of the normal 
+        vector.
+
+Example:
+    We use Bunny point cloud to compute its normals and orient them consistently.
+    The initial reconstruction adheres to Hoppe's algorithm (raw), whereas the 
+    second reconstruction utilises the lambda and cos_alpha_tol parameters. 
+    Due to the high density of the Bunny point cloud available in Open3D a larger
+    value of the parameter k is employed to test the algorithm.  Usually you do 
+    not have at disposal such a refined point clouds, thus you cannot find a 
+    proper choice of k: refer to 
+    https://eugeniovaretti.github.io/meshreco for these cases.::
+
+        import open3d as o3d
+        import numpy as np
+        # Load point cloud
+        data = o3d.data.BunnyMesh()
+
+        # Case 1, Hoppe (raw):
+        pcd = o3d.io.read_point_cloud(data.path)
+
+        # Compute normals and orient them consistently, using k=100 neighbours
+        pcd.estimate_normals()
+        pcd.orient_normals_consistent_tangent_plane(100)
+
+        # Create mesh from point cloud using Poisson Algorithm
+        poisson_mesh = o3d.geometry.TriangleMesh.create_from_point_cloud_poisson(pcd, depth=8, width=0, scale=1.1, linear_fit=False)[0]
+        poisson_mesh.paint_uniform_color(np.array([[0.5],[0.5],[0.5]]))
+        poisson_mesh.compute_vertex_normals()
+        o3d.visualization.draw_geometries([poisson_mesh])
+
+        # Case 2, reconstruction using lambda and cos_alpha_tol parameters:
+        pcd_robust = o3d.io.read_point_cloud(data.path)
+
+        # Compute normals and orient them consistently, using k=100 neighbours
+        pcd_robust.estimate_normals()
+        pcd_robust.orient_normals_consistent_tangent_plane(100, 10, 0.5)
+
+        # Create mesh from point cloud using Poisson Algorithm
+        poisson_mesh_robust = o3d.geometry.TriangleMesh.create_from_point_cloud_poisson(pcd_robust, depth=8, width=0, scale=1.1, linear_fit=False)[0]
+        poisson_mesh_robust.paint_uniform_color(np.array([[0.5],[0.5],[0.5]]))
+        poisson_mesh_robust.compute_vertex_normals()
+
+        o3d.visualization.draw_geometries([poisson_mesh_robust]) )");
     pointcloud.def(
             "estimate_color_gradients", &PointCloud::EstimateColorGradients,
             py::call_guard<py::gil_scoped_release>(), py::arg("max_nn") = 30,
@@ -612,11 +671,6 @@ The implementation is inspired by the PCL implementation. Reference:
             m, "PointCloud", "orient_normals_towards_camera_location",
             {{"camera_location",
               "Normals are oriented with towards the camera_location."}});
-    docstring::ClassMethodDocInject(
-            m, "PointCloud", "orient_normals_consistent_tangent_plane",
-            {{"k",
-              "Number of k nearest neighbors used in constructing the "
-              "Riemannian graph used to propagate normal orientation."}});
     docstring::ClassMethodDocInject(
             m, "PointCloud", "crop",
             {{"aabb", "AxisAlignedBoundingBox to crop points."},