Normal orientation issues in complex geometries

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

@@ -695,12 +695,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

@@ -300,11 +300,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,
@@ -590,11 +649,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."},