Fix TriangleMesh::SamplePointsUniformly not sampling triangles unifor…

cpp/open3d/geometry/TriangleMesh.cpp

@@ -433,20 +433,10 @@ std::shared_ptr<TriangleMesh> TriangleMesh::FilterSmoothTaubin(
 
 std::shared_ptr<PointCloud> TriangleMesh::SamplePointsUniformlyImpl(
         size_t number_of_points,
-        std::vector<double> &triangle_areas,
-        double surface_area,
+        const std::vector<double> &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<size_t> triangle_index_generator(
+            triangle_areas.begin(), triangle_areas.end());
 
     // sample point cloud
     bool has_vert_normal = HasVertexNormals();
@@ -463,35 +453,30 @@ std::shared_ptr<PointCloud> 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<PointCloud> 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<double> 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<PointCloud> TriangleMesh::SamplePointsPoissonDisk(
@@ -545,8 +530,7 @@ std::shared_ptr<PointCloud> TriangleMesh::SamplePointsPoissonDisk(
     std::shared_ptr<PointCloud> 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<PointCloud>();
         pcl->points_ = pcl_init->points_;

cpp/open3d/geometry/TriangleMesh.h

@@ -337,8 +337,7 @@ class TriangleMesh : public MeshBase {
     /// mesh.
     std::shared_ptr<PointCloud> SamplePointsUniformlyImpl(
             size_t number_of_points,
-            std::vector<double> &triangle_areas,
-            double surface_area,
+            const std::vector<double> &triangle_areas,
             bool use_triangle_normal);
 
     /// Function to sample points uniformly from the mesh.

cpp/open3d/utility/Random.h

@@ -171,6 +171,52 @@ class NormalGenerator {
     std::normal_distribution<T> 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<double> weights{1, 2, 3, 4, 5};
+/// utility::random::DiscreteGenerator<size_t> gen(weights.cbegin(),
+/// weights.cend()); for (size_t i = 0; i < 10; i++) {
+///     std::cout << gen() << std::endl;
+/// }
+/// ```
+template <typename T>
+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 <typename InputIt>
+    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<std::mutex> lock(*GetMutex());
+        return distribution_(*GetEngine());
+    }
+
+protected:
+    std::discrete_distribution<T> distribution_;
+};
+
 }  // namespace random
 }  // namespace utility
 }  // namespace open3d