From 56b21969ac95b45dc43666b7fcd6469980aee72b Mon Sep 17 00:00:00 2001 From: eugenio Date: Tue, 23 May 2023 02:44:48 +0200 Subject: [PATCH 01/10] =?UTF-8?q?Modificata=20funzione=20OrientNormals=20i?= =?UTF-8?q?nserendo=20controlli=20sulla=20distanza=20dal=20piano=20-=20due?= =?UTF-8?q?=20metodi:=20cono=20(lambda=3D0,=20cos(a)=20variabile)=20e=20pe?= =?UTF-8?q?nalizzazione=20(lambda!=3D0,=20cos(a)=3D1).=20Possibilit=C3=A0?= =?UTF-8?q?=20di=20mixare=20i=20due=20metodi?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- cpp/open3d/geometry/EstimateNormals.cpp | 65 ++++++++++++++++++++++++- cpp/open3d/geometry/PointCloud.h | 2 +- cpp/pybind/geometry/pointcloud.cpp | 2 +- 3 files changed, 65 insertions(+), 4 deletions(-) diff --git a/cpp/open3d/geometry/EstimateNormals.cpp b/cpp/open3d/geometry/EstimateNormals.cpp index ac084234c25..5bd47a47054 100644 --- a/cpp/open3d/geometry/EstimateNormals.cpp +++ b/cpp/open3d/geometry/EstimateNormals.cpp @@ -8,6 +8,7 @@ #include #include #include +#include #include "open3d/geometry/KDTreeFlann.h" #include "open3d/geometry/PointCloud.h" @@ -359,7 +360,7 @@ void PointCloud::OrientNormalsTowardsCameraLocation( } } -void PointCloud::OrientNormalsConsistentTangentPlane(size_t k) { +void PointCloud::OrientNormalsConsistentTangentPlane(size_t k, const double lambda, const double cos_alpha_tol) { if (!HasNormals()) { utility::LogError( "No normals in the PointCloud. Call EstimateNormals() first."); @@ -380,7 +381,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(); + const auto diff = points_[v0] - points_[v1]; + double dist = diff.squaredNorm() + + lambda * std::abs(diff.dot(normals_[v0])); + + // se lavoriamo con treshold sull'angolo (lambda=0), viene escluso il lato + // che collega due punti che stanno nel cono definito dalla normale e da alpha + // con alpha < alpha_tol + if (lambda == 0) + { + double cos_alpha = std::abs(diff.dot(normals_[v0]))/dist; + if(cos_alpha > cos_alpha_tol) + //return + dist = std::numeric_limits::infinity(); + } delaunay_graph.push_back(WeightedEdge(v0, v1, dist)); graph_edges.insert(edge); } @@ -403,12 +417,33 @@ void PointCloud::OrientNormalsConsistentTangentPlane(size_t k) { edge.weight_ = NormalWeight(edge.v0_, edge.v1_); } + // function to remove outliers from the KNN-> the points that are far from the plane must be removed + auto compute_median = [&](size_t v0, std::vector neighbors) -> double { + std::vector 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 dist_plane_ord = dist_plane; + std::sort(dist_plane_ord.begin(), dist_plane_ord.end()); + double median = dist_plane_ord[static_cast(dist_plane_ord.size()/2)]; + return median; + }; + // Add k nearest neighbors to Riemannian graph KDTreeFlann kdtree(*this); for (size_t v0 = 0; v0 < points_.size(); ++v0) { std::vector neighbors; std::vector dists2; + // Modifications: instead of considering the KNN in terms of euclidean distance + // between points only, we want to include a penalization for the distance from the plane kdtree.SearchKNN(points_[v0], int(k), neighbors, dists2); + double alpha = 2; + double median = lambda == 0 ? std::numeric_limits::quiet_NaN() : compute_median(v0, neighbors); + for (size_t vidx1 = 0; vidx1 < neighbors.size(); ++vidx1) { size_t v1 = size_t(neighbors[vidx1]); if (v0 == v1) { @@ -416,6 +451,28 @@ void PointCloud::OrientNormalsConsistentTangentPlane(size_t k) { } 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() //da ottimizzare + + lambda * normal_dist; + + // se lavoriamo con treshold sull'angolo (lambda=0), viene escluso il lato + // che collega due punti che stanno nel cono definito dalla normale e da alpha + // con alpha < alpha_tol + if (lambda == 0) + { + double cos_alpha = std::abs(diff.dot(normals_[v0]))/dist; + if(cos_alpha > cos_alpha_tol) + continue; + } + else + { + if (normal_dist > alpha * median) + continue; + } double weight = NormalWeight(v0, v1); mst.push_back(WeightedEdge(v0, v1, weight)); graph_edges.insert(edge); @@ -423,6 +480,9 @@ void PointCloud::OrientNormalsConsistentTangentPlane(size_t k) { } } + + + // extract MST from Riemannian graph mst = Kruskal(mst, points_.size()); @@ -457,6 +517,7 @@ void PointCloud::OrientNormalsConsistentTangentPlane(size_t k) { n1 *= -1; } }; + std::cout << "Procedo con tangenti NUOVE open3D" << std::endl; TestAndOrientNormal(Eigen::Vector3d(0, 0, 1), normals_[v0]); while (!traversal_queue.empty()) { v0 = traversal_queue.front(); diff --git a/cpp/open3d/geometry/PointCloud.h b/cpp/open3d/geometry/PointCloud.h index b55840106b6..48df167a5c9 100644 --- a/cpp/open3d/geometry/PointCloud.h +++ b/cpp/open3d/geometry/PointCloud.h @@ -251,7 +251,7 @@ class PointCloud : public Geometry3D { /// /// \param k k nearest neighbour for graph reconstruction for normal /// propagation. - void OrientNormalsConsistentTangentPlane(size_t k); + void OrientNormalsConsistentTangentPlane(size_t k, const double lambda = 0, const double cos_alpha_tol = 1); /// \brief Function to compute the point to point distances between point /// clouds. diff --git a/cpp/pybind/geometry/pointcloud.cpp b/cpp/pybind/geometry/pointcloud.cpp index d9e2ce935c6..695ef1396bd 100644 --- a/cpp/pybind/geometry/pointcloud.cpp +++ b/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,"lamda"_a,"cos_alpha_tol"_a) .def("compute_point_cloud_distance", &PointCloud::ComputePointCloudDistance, "For each point in the source point cloud, compute the " From 7f9b419f3f86cfb9b29e5f9b8e62fdff81cf3415 Mon Sep 17 00:00:00 2001 From: eugenio Date: Wed, 24 May 2023 16:48:06 +0200 Subject: [PATCH 02/10] Cambiata regola di orientazione del primo punto-normale della mesh. Da maxz a minz --- cpp/open3d/geometry/EstimateNormals.cpp | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/cpp/open3d/geometry/EstimateNormals.cpp b/cpp/open3d/geometry/EstimateNormals.cpp index 5bd47a47054..ab5c8453f87 100644 --- a/cpp/open3d/geometry/EstimateNormals.cpp +++ b/cpp/open3d/geometry/EstimateNormals.cpp @@ -496,13 +496,15 @@ void PointCloud::OrientNormalsConsistentTangentPlane(size_t k, const double lamb } // find start node for tree traversal - // init with node that maximizes z - double max_z = std::numeric_limits::lowest(); + // init with node that maximizes z + + //change minimize instead of maximize + double min_z = std::numeric_limits::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; } } @@ -518,7 +520,7 @@ void PointCloud::OrientNormalsConsistentTangentPlane(size_t k, const double lamb } }; std::cout << "Procedo con tangenti NUOVE open3D" << std::endl; - 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(); From 7cc62c35e6580e461d2c5bd00d60cd7e63d6dde4 Mon Sep 17 00:00:00 2001 From: eugenio Date: Sat, 27 May 2023 00:56:48 +0200 Subject: [PATCH 03/10] Miswritngs and deletion of io --- cpp/open3d/geometry/EstimateNormals.cpp | 2 -- cpp/pybind/geometry/pointcloud.cpp | 2 +- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/cpp/open3d/geometry/EstimateNormals.cpp b/cpp/open3d/geometry/EstimateNormals.cpp index ab5c8453f87..d5ae5b2247a 100644 --- a/cpp/open3d/geometry/EstimateNormals.cpp +++ b/cpp/open3d/geometry/EstimateNormals.cpp @@ -8,7 +8,6 @@ #include #include #include -#include #include "open3d/geometry/KDTreeFlann.h" #include "open3d/geometry/PointCloud.h" @@ -519,7 +518,6 @@ void PointCloud::OrientNormalsConsistentTangentPlane(size_t k, const double lamb n1 *= -1; } }; - std::cout << "Procedo con tangenti NUOVE open3D" << std::endl; TestAndOrientNormal(Eigen::Vector3d(0, 0, -1), normals_[v0]); while (!traversal_queue.empty()) { v0 = traversal_queue.front(); diff --git a/cpp/pybind/geometry/pointcloud.cpp b/cpp/pybind/geometry/pointcloud.cpp index 695ef1396bd..9ba3c747209 100644 --- a/cpp/pybind/geometry/pointcloud.cpp +++ b/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,"lamda"_a,"cos_alpha_tol"_a) + "k"_a,"lambda"_a,"cos_alpha_tol"_a) .def("compute_point_cloud_distance", &PointCloud::ComputePointCloudDistance, "For each point in the source point cloud, compute the " From 5220dcb9299dc552b8a1ccb0bfbce9af6d716f29 Mon Sep 17 00:00:00 2001 From: eugenio Date: Wed, 31 May 2023 22:01:09 +0200 Subject: [PATCH 04/10] modified penalization (now use IQR to identify outliers). Added the possibility to mix penalization and treshold approach. --- cpp/open3d/geometry/EstimateNormals.cpp | 75 ++++++++++++++----------- 1 file changed, 42 insertions(+), 33 deletions(-) diff --git a/cpp/open3d/geometry/EstimateNormals.cpp b/cpp/open3d/geometry/EstimateNormals.cpp index d5ae5b2247a..5b01332175b 100644 --- a/cpp/open3d/geometry/EstimateNormals.cpp +++ b/cpp/open3d/geometry/EstimateNormals.cpp @@ -359,7 +359,9 @@ void PointCloud::OrientNormalsTowardsCameraLocation( } } -void PointCloud::OrientNormalsConsistentTangentPlane(size_t k, const double lambda, const double cos_alpha_tol) { +void PointCloud::OrientNormalsConsistentTangentPlane(size_t k, + const double lambda /* = 0*/, + const double cos_alpha_tol /* = 1*/) { if (!HasNormals()) { utility::LogError( "No normals in the PointCloud. Call EstimateNormals() first."); @@ -381,19 +383,16 @@ void PointCloud::OrientNormalsConsistentTangentPlane(size_t k, const double lamb size_t edge = EdgeIndex(v0, v1); if (graph_edges.count(edge) == 0) { const auto diff = points_[v0] - points_[v1]; + // penalization on normal-plane distance double dist = diff.squaredNorm() + lambda * std::abs(diff.dot(normals_[v0])); - // se lavoriamo con treshold sull'angolo (lambda=0), viene escluso il lato - // che collega due punti che stanno nel cono definito dalla normale e da alpha - // con alpha < alpha_tol - if (lambda == 0) - { - double cos_alpha = std::abs(diff.dot(normals_[v0]))/dist; - if(cos_alpha > cos_alpha_tol) - //return - dist = std::numeric_limits::infinity(); - } + // 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]))/dist; + if(cos_alpha > cos_alpha_tol) + dist = std::numeric_limits::infinity(); + delaunay_graph.push_back(WeightedEdge(v0, v1, dist)); graph_edges.insert(edge); } @@ -417,7 +416,7 @@ void PointCloud::OrientNormalsConsistentTangentPlane(size_t k, const double lamb } // function to remove outliers from the KNN-> the points that are far from the plane must be removed - auto compute_median = [&](size_t v0, std::vector neighbors) -> double { + auto compute_q1q3= [&](size_t v0, std::vector neighbors) -> std::array { std::vector dist_plane; for (size_t vidx1 = 0; vidx1 < neighbors.size(); ++vidx1) { @@ -428,8 +427,19 @@ void PointCloud::OrientNormalsConsistentTangentPlane(size_t k, const double lamb } std::vector dist_plane_ord = dist_plane; std::sort(dist_plane_ord.begin(), dist_plane_ord.end()); - double median = dist_plane_ord[static_cast(dist_plane_ord.size()/2)]; - return median; + // calculate quartiles + int q1_idx = static_cast(dist_plane_ord.size() * 0.25); + double q1 = dist_plane_ord[q1_idx]; + + int q3_idx = static_cast(dist_plane_ord.size() * 0.75); + double q3 = dist_plane_ord[q3_idx]; + + std::array q1q3; + q1q3[0] = q1; + q1q3[1] = q3; + + + return q1q3; }; // Add k nearest neighbors to Riemannian graph @@ -437,11 +447,13 @@ void PointCloud::OrientNormalsConsistentTangentPlane(size_t k, const double lamb for (size_t v0 = 0; v0 < points_.size(); ++v0) { std::vector neighbors; std::vector dists2; + // Modifications: instead of considering the KNN in terms of euclidean distance // between points only, we want to include a penalization for the distance from the plane kdtree.SearchKNN(points_[v0], int(k), neighbors, dists2); - double alpha = 2; - double median = lambda == 0 ? std::numeric_limits::quiet_NaN() : compute_median(v0, neighbors); + + const double DEFAULT_VALUE = std::numeric_limits::quiet_NaN(); + std::array q1q3 = lambda == 0 ? std::array{DEFAULT_VALUE, DEFAULT_VALUE} : compute_q1q3(v0, neighbors); for (size_t vidx1 = 0; vidx1 < neighbors.size(); ++vidx1) { size_t v1 = size_t(neighbors[vidx1]); @@ -455,23 +467,23 @@ void PointCloud::OrientNormalsConsistentTangentPlane(size_t k, const double lamb const auto diff = points_[v0] - points_[v1]; double normal_dist = std::abs(diff.dot(normals_[v0])); - double dist = diff.squaredNorm() //da ottimizzare + double dist = diff.squaredNorm() + lambda * normal_dist; - // se lavoriamo con treshold sull'angolo (lambda=0), viene escluso il lato - // che collega due punti che stanno nel cono definito dalla normale e da alpha - // con alpha < alpha_tol - if (lambda == 0) - { - double cos_alpha = std::abs(diff.dot(normals_[v0]))/dist; - if(cos_alpha > cos_alpha_tol) - continue; - } - else + // 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]))/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) { - if (normal_dist > alpha * median) + 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); @@ -480,8 +492,6 @@ void PointCloud::OrientNormalsConsistentTangentPlane(size_t k, const double lamb } - - // extract MST from Riemannian graph mst = Kruskal(mst, points_.size()); @@ -495,9 +505,8 @@ void PointCloud::OrientNormalsConsistentTangentPlane(size_t k, const double lamb } // find start node for tree traversal - // init with node that maximizes z + // init with node that minimizes z - //change minimize instead of maximize double min_z = std::numeric_limits::max(); size_t v0 = 0; for (size_t vidx = 0; vidx < points_.size(); ++vidx) { From 19a28269bca0afdc517fe42d2c31595b8eac91a4 Mon Sep 17 00:00:00 2001 From: eugenio Date: Mon, 5 Jun 2023 09:59:44 +0200 Subject: [PATCH 05/10] given correct default value, updated pybind --- cpp/open3d/geometry/PointCloud.h | 2 +- cpp/pybind/geometry/pointcloud.cpp | 4 +++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/cpp/open3d/geometry/PointCloud.h b/cpp/open3d/geometry/PointCloud.h index 48df167a5c9..62dc8e0a829 100644 --- a/cpp/open3d/geometry/PointCloud.h +++ b/cpp/open3d/geometry/PointCloud.h @@ -251,7 +251,7 @@ class PointCloud : public Geometry3D { /// /// \param k k nearest neighbour for graph reconstruction for normal /// propagation. - void OrientNormalsConsistentTangentPlane(size_t k, const double lambda = 0, const double cos_alpha_tol = 1); + 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. diff --git a/cpp/pybind/geometry/pointcloud.cpp b/cpp/pybind/geometry/pointcloud.cpp index 9ba3c747209..3f406d2bc4b 100644 --- a/cpp/pybind/geometry/pointcloud.cpp +++ b/cpp/pybind/geometry/pointcloud.cpp @@ -139,7 +139,9 @@ void pybind_pointcloud(py::module &m) { &PointCloud::OrientNormalsConsistentTangentPlane, "Function to orient the normals with respect to consistent " "tangent planes", - "k"_a,"lambda"_a,"cos_alpha_tol"_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 " From a125cef892f886ac055691525f2232f8007733c2 Mon Sep 17 00:00:00 2001 From: eugenio Date: Mon, 5 Jun 2023 10:00:39 +0200 Subject: [PATCH 06/10] solved some normalization criticalities. Cleaned the code. --- cpp/open3d/geometry/EstimateNormals.cpp | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/cpp/open3d/geometry/EstimateNormals.cpp b/cpp/open3d/geometry/EstimateNormals.cpp index 5b01332175b..001ab8ae798 100644 --- a/cpp/open3d/geometry/EstimateNormals.cpp +++ b/cpp/open3d/geometry/EstimateNormals.cpp @@ -360,8 +360,8 @@ void PointCloud::OrientNormalsTowardsCameraLocation( } void PointCloud::OrientNormalsConsistentTangentPlane(size_t k, - const double lambda /* = 0*/, - const double cos_alpha_tol /* = 1*/) { + const double lambda /* = 0.0*/, + const double cos_alpha_tol /* = 1.0*/) { if (!HasNormals()) { utility::LogError( "No normals in the PointCloud. Call EstimateNormals() first."); @@ -384,16 +384,16 @@ void PointCloud::OrientNormalsConsistentTangentPlane(size_t k, if (graph_edges.count(edge) == 0) { const auto diff = points_[v0] - points_[v1]; // penalization on normal-plane distance - double dist = diff.squaredNorm() - + lambda * std::abs(diff.dot(normals_[v0])); + 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]))/dist; + double cos_alpha = std::abs(diff.dot(normals_[v0]))/std::sqrt(dist); if(cos_alpha > cos_alpha_tol) dist = std::numeric_limits::infinity(); - delaunay_graph.push_back(WeightedEdge(v0, v1, dist)); + delaunay_graph.push_back(WeightedEdge(v0, v1, dist + penalization)); graph_edges.insert(edge); } }; @@ -415,7 +415,8 @@ void PointCloud::OrientNormalsConsistentTangentPlane(size_t k, edge.weight_ = NormalWeight(edge.v0_, edge.v1_); } - // function to remove outliers from the KNN-> the points that are far from the plane must be removed + // 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 neighbors) -> std::array { std::vector dist_plane; @@ -448,8 +449,6 @@ void PointCloud::OrientNormalsConsistentTangentPlane(size_t k, std::vector neighbors; std::vector dists2; - // Modifications: instead of considering the KNN in terms of euclidean distance - // between points only, we want to include a penalization for the distance from the plane kdtree.SearchKNN(points_[v0], int(k), neighbors, dists2); const double DEFAULT_VALUE = std::numeric_limits::quiet_NaN(); @@ -467,12 +466,11 @@ void PointCloud::OrientNormalsConsistentTangentPlane(size_t k, const auto diff = points_[v0] - points_[v1]; double normal_dist = std::abs(diff.dot(normals_[v0])); - double dist = diff.squaredNorm() - + lambda * normal_dist; + 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]))/dist; + double cos_alpha = std::abs(diff.dot(normals_[v0]))/std::sqrt(dist); if(cos_alpha > cos_alpha_tol) continue; From ff5a7d7009720d380c633edbc9674f367a148277 Mon Sep 17 00:00:00 2001 From: eugenio Date: Tue, 6 Jun 2023 17:43:05 +0200 Subject: [PATCH 07/10] Applied style to the code --- cpp/open3d/geometry/EstimateNormals.cpp | 61 +++++++++++++------------ cpp/open3d/geometry/PointCloud.h | 10 +++- cpp/open3d/t/geometry/PointCloud.cpp | 7 ++- cpp/open3d/t/geometry/PointCloud.h | 10 +++- cpp/pybind/geometry/pointcloud.cpp | 4 +- 5 files changed, 57 insertions(+), 35 deletions(-) diff --git a/cpp/open3d/geometry/EstimateNormals.cpp b/cpp/open3d/geometry/EstimateNormals.cpp index 001ab8ae798..5bc9abc0ba3 100644 --- a/cpp/open3d/geometry/EstimateNormals.cpp +++ b/cpp/open3d/geometry/EstimateNormals.cpp @@ -359,7 +359,8 @@ 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()) { @@ -387,12 +388,14 @@ void PointCloud::OrientNormalsConsistentTangentPlane(size_t k, 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) + // 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::infinity(); - + delaunay_graph.push_back(WeightedEdge(v0, v1, dist + penalization)); graph_edges.insert(edge); } @@ -416,8 +419,11 @@ void PointCloud::OrientNormalsConsistentTangentPlane(size_t k, } // 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 neighbors) -> std::array { + // 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 neighbors) -> std::array { std::vector dist_plane; for (size_t vidx1 = 0; vidx1 < neighbors.size(); ++vidx1) { @@ -439,9 +445,8 @@ void PointCloud::OrientNormalsConsistentTangentPlane(size_t k, q1q3[0] = q1; q1q3[1] = q3; - return q1q3; - }; + }; // Add k nearest neighbors to Riemannian graph KDTreeFlann kdtree(*this); @@ -452,7 +457,10 @@ void PointCloud::OrientNormalsConsistentTangentPlane(size_t k, kdtree.SearchKNN(points_[v0], int(k), neighbors, dists2); const double DEFAULT_VALUE = std::numeric_limits::quiet_NaN(); - std::array q1q3 = lambda == 0 ? std::array{DEFAULT_VALUE, DEFAULT_VALUE} : compute_q1q3(v0, neighbors); + std::array q1q3 = + lambda == 0 + ? std::array{DEFAULT_VALUE, DEFAULT_VALUE} + : compute_q1q3(v0, neighbors); for (size_t vidx1 = 0; vidx1 < neighbors.size(); ++vidx1) { size_t v1 = size_t(neighbors[vidx1]); @@ -461,26 +469,24 @@ void PointCloud::OrientNormalsConsistentTangentPlane(size_t k, } 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; - } + // 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)); @@ -489,7 +495,6 @@ void PointCloud::OrientNormalsConsistentTangentPlane(size_t k, } } - // extract MST from Riemannian graph mst = Kruskal(mst, points_.size()); @@ -503,7 +508,7 @@ void PointCloud::OrientNormalsConsistentTangentPlane(size_t k, } // find start node for tree traversal - // init with node that minimizes z + // init with node that minimizes z double min_z = std::numeric_limits::max(); size_t v0 = 0; diff --git a/cpp/open3d/geometry/PointCloud.h b/cpp/open3d/geometry/PointCloud.h index 62dc8e0a829..beb48fc4361 100644 --- a/cpp/open3d/geometry/PointCloud.h +++ b/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, const double lambda = 0.0, const double cos_alpha_tol = 1.0); + /// \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. diff --git a/cpp/open3d/t/geometry/PointCloud.cpp b/cpp/open3d/t/geometry/PointCloud.cpp index 6ec721ffa4a..21aecff825d 100644 --- a/cpp/open3d/t/geometry/PointCloud.cpp +++ b/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())); diff --git a/cpp/open3d/t/geometry/PointCloud.h b/cpp/open3d/t/geometry/PointCloud.h index 68827fa489c..2bf134b8b43 100644 --- a/cpp/open3d/t/geometry/PointCloud.h +++ b/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. diff --git a/cpp/pybind/geometry/pointcloud.cpp b/cpp/pybind/geometry/pointcloud.cpp index 3f406d2bc4b..b8aa24eda0b 100644 --- a/cpp/pybind/geometry/pointcloud.cpp +++ b/cpp/pybind/geometry/pointcloud.cpp @@ -139,9 +139,7 @@ void pybind_pointcloud(py::module &m) { &PointCloud::OrientNormalsConsistentTangentPlane, "Function to orient the normals with respect to consistent " "tangent planes", - "k"_a, - "lambda"_a = 0.0, - "cos_alpha_tol"_a = 1.0) + "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 " From e8c70efec82406bd09708a2b69281d0c34a2f0e1 Mon Sep 17 00:00:00 2001 From: eugenio Date: Tue, 6 Jun 2023 20:55:00 +0200 Subject: [PATCH 08/10] missing pybind default values t/geometry --- cpp/pybind/t/geometry/pointcloud.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/pybind/t/geometry/pointcloud.cpp b/cpp/pybind/t/geometry/pointcloud.cpp index cc7f98eca59..e0395638e1c 100644 --- a/cpp/pybind/t/geometry/pointcloud.cpp +++ b/cpp/pybind/t/geometry/pointcloud.cpp @@ -304,7 +304,7 @@ infinite value. It also removes the corresponding attributes. &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); pointcloud.def( "estimate_color_gradients", &PointCloud::EstimateColorGradients, py::call_guard(), py::arg("max_nn") = 30, From 0c31f475a5af4adee333d1a9f30931f847ba3fcb Mon Sep 17 00:00:00 2001 From: eugenio Date: Mon, 25 Sep 2023 10:47:24 +0200 Subject: [PATCH 09/10] Added Python doc --- cpp/pybind/t/geometry/pointcloud.cpp | 53 +++++++++++++++++++++++++--- 1 file changed, 48 insertions(+), 5 deletions(-) diff --git a/cpp/pybind/t/geometry/pointcloud.cpp b/cpp/pybind/t/geometry/pointcloud.cpp index e0395638e1c..cc82f872204 100644 --- a/cpp/pybind/t/geometry/pointcloud.cpp +++ b/cpp/pybind/t/geometry/pointcloud.cpp @@ -300,11 +300,54 @@ 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, "lambda"_a = 0.0, "cos_alpha_tol"_a = 1.0); + 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 + as 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. Number of neighbors to use for tangent plane estimation. + lambda. 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. 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. + + # 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::arg("max_nn") = 30, From afd9a5f5343d8709c10e9a0cdb4b88a3733b51e2 Mon Sep 17 00:00:00 2001 From: Benjamin Ummenhofer Date: Tue, 26 Sep 2023 14:13:45 +0200 Subject: [PATCH 10/10] small format changes --- cpp/pybind/t/geometry/pointcloud.cpp | 45 +++++++++++++++++----------- 1 file changed, 28 insertions(+), 17 deletions(-) diff --git a/cpp/pybind/t/geometry/pointcloud.cpp b/cpp/pybind/t/geometry/pointcloud.cpp index cc82f872204..51a61187d0f 100644 --- a/cpp/pybind/t/geometry/pointcloud.cpp +++ b/cpp/pybind/t/geometry/pointcloud.cpp @@ -304,21 +304,37 @@ infinite value. It also removes the corresponding attributes. "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 - as 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). + 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. Number of neighbors to use for tangent plane estimation. - lambda. 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. Cosine threshold angle used to determine the inclusion boundary of neighbors based on the direction of the normal vector. + 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. + 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() @@ -346,8 +362,8 @@ infinite value. It also removes the corresponding attributes. 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])"); + o3d.visualization.draw_geometries([poisson_mesh_robust]) )"); pointcloud.def( "estimate_color_gradients", &PointCloud::EstimateColorGradients, py::call_guard(), py::arg("max_nn") = 30, @@ -633,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."},