From 6cdc708353ff8a95d2b9019ea359875048459209 Mon Sep 17 00:00:00 2001 From: iyadahmed Date: Mon, 12 Feb 2024 12:47:04 +0200 Subject: [PATCH] Fix TriangleMesh::SamplePointsUniformly not sampling triangles uniformly (#6144) --- cpp/open3d/geometry/TriangleMesh.cpp | 76 +++++++++++----------------- cpp/open3d/geometry/TriangleMesh.h | 3 +- cpp/open3d/utility/Random.h | 46 +++++++++++++++++ 3 files changed, 77 insertions(+), 48 deletions(-) diff --git a/cpp/open3d/geometry/TriangleMesh.cpp b/cpp/open3d/geometry/TriangleMesh.cpp index 48faa3b6a19..88f602770cb 100644 --- a/cpp/open3d/geometry/TriangleMesh.cpp +++ b/cpp/open3d/geometry/TriangleMesh.cpp @@ -433,20 +433,10 @@ std::shared_ptr TriangleMesh::FilterSmoothTaubin( std::shared_ptr TriangleMesh::SamplePointsUniformlyImpl( size_t number_of_points, - std::vector &triangle_areas, - double surface_area, + const std::vector &triangle_areas, bool use_triangle_normal) { - if (surface_area <= 0) { - utility::LogError("Invalid surface area {}, it must be > 0.", - surface_area); - } - - // triangle areas to cdf - triangle_areas[0] /= surface_area; - for (size_t tidx = 1; tidx < triangles_.size(); ++tidx) { - triangle_areas[tidx] = - triangle_areas[tidx] / surface_area + triangle_areas[tidx - 1]; - } + utility::random::DiscreteGenerator triangle_index_generator( + triangle_areas.begin(), triangle_areas.end()); // sample point cloud bool has_vert_normal = HasVertexNormals(); @@ -463,35 +453,30 @@ std::shared_ptr TriangleMesh::SamplePointsUniformlyImpl( if (has_vert_color) { pcd->colors_.resize(number_of_points); } - size_t point_idx = 0; - for (size_t tidx = 0; tidx < triangles_.size(); ++tidx) { - size_t n = size_t(std::round(triangle_areas[tidx] * number_of_points)); - while (point_idx < n) { - double r1 = uniform_generator(); - double r2 = uniform_generator(); - double a = (1 - std::sqrt(r1)); - double b = std::sqrt(r1) * (1 - r2); - double c = std::sqrt(r1) * r2; - - const Eigen::Vector3i &triangle = triangles_[tidx]; - pcd->points_[point_idx] = a * vertices_[triangle(0)] + - b * vertices_[triangle(1)] + - c * vertices_[triangle(2)]; - if (has_vert_normal && !use_triangle_normal) { - pcd->normals_[point_idx] = a * vertex_normals_[triangle(0)] + - b * vertex_normals_[triangle(1)] + - c * vertex_normals_[triangle(2)]; - } - if (use_triangle_normal) { - pcd->normals_[point_idx] = triangle_normals_[tidx]; - } - if (has_vert_color) { - pcd->colors_[point_idx] = a * vertex_colors_[triangle(0)] + - b * vertex_colors_[triangle(1)] + - c * vertex_colors_[triangle(2)]; - } - point_idx++; + for (size_t point_idx = 0; point_idx < number_of_points; ++point_idx) { + double r1 = uniform_generator(); + double r2 = uniform_generator(); + double a = (1 - std::sqrt(r1)); + double b = std::sqrt(r1) * (1 - r2); + double c = std::sqrt(r1) * r2; + size_t tidx = triangle_index_generator(); + const Eigen::Vector3i &triangle = triangles_[tidx]; + pcd->points_[point_idx] = a * vertices_[triangle(0)] + + b * vertices_[triangle(1)] + + c * vertices_[triangle(2)]; + if (has_vert_normal && !use_triangle_normal) { + pcd->normals_[point_idx] = a * vertex_normals_[triangle(0)] + + b * vertex_normals_[triangle(1)] + + c * vertex_normals_[triangle(2)]; + } + if (use_triangle_normal) { + pcd->normals_[point_idx] = triangle_normals_[tidx]; + } + if (has_vert_color) { + pcd->colors_[point_idx] = a * vertex_colors_[triangle(0)] + + b * vertex_colors_[triangle(1)] + + c * vertex_colors_[triangle(2)]; } } @@ -507,12 +492,12 @@ std::shared_ptr TriangleMesh::SamplePointsUniformly( utility::LogError("Input mesh has no triangles."); } - // Compute area of each triangle and sum surface area + // Compute area of each triangle std::vector triangle_areas; - double surface_area = GetSurfaceArea(triangle_areas); + GetSurfaceArea(triangle_areas); return SamplePointsUniformlyImpl(number_of_points, triangle_areas, - surface_area, use_triangle_normal); + use_triangle_normal); } std::shared_ptr TriangleMesh::SamplePointsPoissonDisk( @@ -545,8 +530,7 @@ std::shared_ptr TriangleMesh::SamplePointsPoissonDisk( std::shared_ptr pcl; if (pcl_init == nullptr) { pcl = SamplePointsUniformlyImpl(size_t(init_factor * number_of_points), - triangle_areas, surface_area, - use_triangle_normal); + triangle_areas, use_triangle_normal); } else { pcl = std::make_shared(); pcl->points_ = pcl_init->points_; diff --git a/cpp/open3d/geometry/TriangleMesh.h b/cpp/open3d/geometry/TriangleMesh.h index 8448cccbcc0..ba3105f0d60 100644 --- a/cpp/open3d/geometry/TriangleMesh.h +++ b/cpp/open3d/geometry/TriangleMesh.h @@ -337,8 +337,7 @@ class TriangleMesh : public MeshBase { /// mesh. std::shared_ptr SamplePointsUniformlyImpl( size_t number_of_points, - std::vector &triangle_areas, - double surface_area, + const std::vector &triangle_areas, bool use_triangle_normal); /// Function to sample points uniformly from the mesh. diff --git a/cpp/open3d/utility/Random.h b/cpp/open3d/utility/Random.h index 87a0eca4d5e..1a3411d29a2 100644 --- a/cpp/open3d/utility/Random.h +++ b/cpp/open3d/utility/Random.h @@ -171,6 +171,52 @@ class NormalGenerator { std::normal_distribution distribution_; }; +/// Generate discretely distributed integer values according to a range of +/// weight values. +/// This class is globally seeded by utility::random::Seed(). +/// This class is a wrapper around std::discrete_distribution. +/// +/// Example: +/// ```cpp +/// #include "open3d/utility/Random.h" +/// +/// // Globally seed Open3D. This will affect all random functions. +/// utility::random::Seed(0); +/// +/// // Weighted random choice of size_t +/// std::vector weights{1, 2, 3, 4, 5}; +/// utility::random::DiscreteGenerator gen(weights.cbegin(), +/// weights.cend()); for (size_t i = 0; i < 10; i++) { +/// std::cout << gen() << std::endl; +/// } +/// ``` +template +class DiscreteGenerator { +public: + /// Generate discretely distributed integer values according to a range of + /// weight values. + /// \param first The iterator or pointer pointing to the first element in + /// the range of weights. + /// \param last The iterator or pointer pointing to one past the last + /// element in the range of weights. + template + DiscreteGenerator(InputIt first, InputIt last) + : distribution_(first, last) { + if (first > last) { + utility::LogError("first must be <= last."); + } + } + + /// Call this to generate a discretely distributed integer value. + T operator()() { + std::lock_guard lock(*GetMutex()); + return distribution_(*GetEngine()); + } + +protected: + std::discrete_distribution distribution_; +}; + } // namespace random } // namespace utility } // namespace open3d