Skip to content

Fix #7449: sample_points_poisson_disk is extremely slow #7450

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
nikste wants to merge 4 commits into
base: main
Choose a base branch
from
Open

Fix #7449: sample_points_poisson_disk is extremely slow #7450

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
aa0348e
Fix sample_points_poisson_disk performance (issue #7449)
klaus-kode Mar 1, 2026
b54b266
Update changelog
nikste Mar 1, 2026
70044a4
Move vector declaration outside of loop
nikste Mar 5, 2026
e2433f6
Add unit tests, add optional unsorted parameter to KDtree search
nikste Mar 5, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -73,6 +73,7 @@
- macOS x86_64 not longer supported, only macOS arm64 is supported.
- Python 3.13+3.14 support
- Fix Windows build failure for PyTorch ops due to PyTorch's bundled fmt (v11+) requiring `/utf-8` with MSVC (PR #7447)
- Fix `TriangleMesh::SamplePointsPoissonDisk` performance by incrementally updating neighbor weights instead of recomputing them with additional KD-tree queries (issue #7449)


## 0.13
Expand Down
11 changes: 7 additions & 4 deletions cpp/open3d/geometry/KDTreeFlann.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -115,7 +115,8 @@ template <typename T>
int KDTreeFlann::SearchRadius(const T &query,
double radius,
std::vector<int> &indices,
std::vector<double> &distance2) const {
std::vector<double> &distance2,
bool sorted) const {
// This is optimized code for heavily repeated search.
// Since max_nn is not given, we let flann to do its own memory management.
// Other flann::Index::radiusSearch() implementations lose performance due
Expand All @@ -126,7 +127,7 @@ int KDTreeFlann::SearchRadius(const T &query,
std::vector<nanoflann::ResultItem<Eigen::Index, double>> indices_dists;
int k = nanoflann_index_->index_->radiusSearch(
query.data(), radius * radius, indices_dists,
nanoflann::SearchParameters(0.0));
nanoflann::SearchParameters(0.0, sorted));
indices.resize(k);
distance2.resize(k);
for (int i = 0; i < k; ++i) {
Expand Down Expand Up @@ -187,7 +188,8 @@ template int KDTreeFlann::SearchRadius<Eigen::Vector3d>(
const Eigen::Vector3d &query,
double radius,
std::vector<int> &indices,
std::vector<double> &distance2) const;
std::vector<double> &distance2,
bool sorted) const;
template int KDTreeFlann::SearchHybrid<Eigen::Vector3d>(
const Eigen::Vector3d &query,
double radius,
Expand All @@ -209,7 +211,8 @@ template int KDTreeFlann::SearchRadius<Eigen::VectorXd>(
const Eigen::VectorXd &query,
double radius,
std::vector<int> &indices,
std::vector<double> &distance2) const;
std::vector<double> &distance2,
bool sorted) const;
template int KDTreeFlann::SearchHybrid<Eigen::VectorXd>(
const Eigen::VectorXd &query,
double radius,
Expand Down
3 changes: 2 additions & 1 deletion cpp/open3d/geometry/KDTreeFlann.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -80,7 +80,8 @@ class KDTreeFlann {
int SearchRadius(const T &query,
double radius,
std::vector<int> &indices,
std::vector<double> &distance2) const;
std::vector<double> &distance2,
bool sorted = true) const;

template <typename T>
Copy link

Copilot AI Mar 6, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

KDTreeFlann::SearchRadius’s signature was changed by adding a new parameter. Even with a default value, this is an ABI-breaking change (different mangled symbol) for any downstream C++ code built against the previous Open3D library. To avoid breaking existing binaries, consider keeping the original 4-arg overload and adding a new 5-arg overload that forwards internally (and explicitly instantiate both overloads as needed).

Suggested change
template <typename T>
template <typename T>
int SearchRadius(const T &query,
double radius,
std::vector<int> &indices,
std::vector<double> &distance2) const {
return SearchRadius(query, radius, indices, distance2, true);
}
template <typename T>

Copilot uses AI. Check for mistakes.
int SearchHybrid(const T &query,
Expand Down
24 changes: 15 additions & 9 deletions cpp/open3d/geometry/TriangleMesh.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -588,10 +588,11 @@ std::shared_ptr<PointCloud> TriangleMesh::SamplePointsPoissonDisk(
return std::pow(1 - d / r_max, alpha);
};

std::vector<int> nbs;
std::vector<double> dists2;

auto ComputePointWeight = [&](int pidx0) {
std::vector<int> nbs;
std::vector<double> dists2;
kdtree.SearchRadius(pcl->points_[pidx0], r_max, nbs, dists2);
kdtree.SearchRadius(pcl->points_[pidx0], r_max, nbs, dists2, false);
double weight = 0;
for (size_t nbidx = 0; nbidx < nbs.size(); ++nbidx) {
int pidx1 = nbs[nbidx];
Expand Down Expand Up @@ -635,12 +636,17 @@ std::shared_ptr<PointCloud> TriangleMesh::SamplePointsPoissonDisk(
deleted[pidx] = true;
current_number_of_points--;

// update weights
std::vector<int> nbs;
std::vector<double> dists2;
kdtree.SearchRadius(pcl->points_[pidx], r_max, nbs, dists2);
for (int nb : nbs) {
ComputePointWeight(nb);
// update weights: subtract the contribution of the deleted point from
// its neighbors instead of recomputing each neighbor's weight from
// scratch (which would require an additional KD-tree query per
// neighbor). This matches the reference algorithm in the paper.
kdtree.SearchRadius(pcl->points_[pidx], r_max, nbs, dists2, false);
for (size_t nbidx = 0; nbidx < nbs.size(); ++nbidx) {
int nb = nbs[nbidx];
if (deleted[nb] || nb == pidx) {
continue;
}
weights[nb] -= WeightFcn(dists2[nbidx]);
queue.push(QueueEntry(nb, weights[nb]));
}
}
Expand Down
68 changes: 68 additions & 0 deletions cpp/tests/geometry/TriangleMesh.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1110,6 +1110,74 @@ TEST(TriangleMesh, SamplePointsUniformly) {
}
}

TEST(TriangleMesh, SamplePointsPoissonDisk) {
// Empty mesh should throw
auto mesh_empty = geometry::TriangleMesh();
EXPECT_THROW(mesh_empty.SamplePointsPoissonDisk(100), std::runtime_error);

// Simple triangle mesh in z=0 plane: (0,0,0), (1,0,0), (0,1,0)
std::vector<Eigen::Vector3d> vertices = {{0, 0, 0}, {1, 0, 0}, {0, 1, 0}};
std::vector<Eigen::Vector3i> triangles = {{0, 1, 2}};

auto mesh_simple = geometry::TriangleMesh();
mesh_simple.vertices_ = vertices;
mesh_simple.triangles_ = triangles;

size_t n_points = 100;
auto pcd = mesh_simple.SamplePointsPoissonDisk(n_points);
EXPECT_EQ(pcd->points_.size(), n_points);
EXPECT_TRUE(pcd->colors_.size() == 0);
EXPECT_TRUE(pcd->normals_.size() == 0);

// All points must lie on the triangle surface (z=0, x>=0, y>=0, x+y<=1)
for (size_t i = 0; i < pcd->points_.size(); ++i) {
const auto &p = pcd->points_[i];
EXPECT_NEAR(p(2), 0.0, 1e-10);
EXPECT_GE(p(0), -1e-10);
EXPECT_GE(p(1), -1e-10);
EXPECT_LE(p(0) + p(1), 1.0 + 1e-10);
}

// Poisson disk property: minimum pairwise distance should be > 0
// (points are spread apart, not clustered)
double min_dist = std::numeric_limits<double>::max();
for (size_t i = 0; i < pcd->points_.size(); ++i) {
for (size_t j = i + 1; j < pcd->points_.size(); ++j) {
double d = (pcd->points_[i] - pcd->points_[j]).norm();
min_dist = std::min(min_dist, d);
Comment on lines +1143 to +1147
Copy link

Copilot AI Mar 6, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This test uses std::numeric_limits and std::min, but this file does not include the standard headers that define them. Please add the appropriate includes (e.g., <limits> and <algorithm>) to avoid relying on transitive includes that may change across platforms/compilers.

Copilot uses AI. Check for mistakes.
}
}
// For 100 points on a triangle with area 0.5, the theoretical Poisson disk
// radius is r = 2*sqrt(area/n / (2*sqrt(3))) ≈ 0.057. After sample
// elimination the actual minimum distance should be a meaningful fraction
// of that. We just check it's not degenerate (> 0.01).
EXPECT_GT(min_dist, 0.01);
Comment on lines +1150 to +1154
Copy link

Copilot AI Mar 6, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

SamplePointsPoissonDisk relies on RNG-based sampling; this test doesn’t set a fixed seed, and the hard min_dist > 0.01 threshold may vary with randomness/implementation details, risking flaky CI. Consider seeding (e.g., utility::random::Seed(0)) and/or making the assertion deterministic (e.g., only require min_dist to be above a tiny epsilon, or compare against an expected value with a fixed seed).

Suggested change
// For 100 points on a triangle with area 0.5, the theoretical Poisson disk
// radius is r = 2*sqrt(area/n / (2*sqrt(3))) ≈ 0.057. After sample
// elimination the actual minimum distance should be a meaningful fraction
// of that. We just check it's not degenerate (> 0.01).
EXPECT_GT(min_dist, 0.01);
// The exact minimum distance depends on RNG and implementation details.
// To avoid flaky tests, only check that points are not degenerate
// (no coincident samples).
EXPECT_GT(min_dist, 1e-8);

Copilot uses AI. Check for mistakes.

// With vertex colors and normals
std::vector<Eigen::Vector3d> colors = {{1, 0, 0}, {1, 0, 0}, {1, 0, 0}};
std::vector<Eigen::Vector3d> normals = {{0, 1, 0}, {0, 1, 0}, {0, 1, 0}};
mesh_simple.vertex_colors_ = colors;
mesh_simple.vertex_normals_ = normals;
pcd = mesh_simple.SamplePointsPoissonDisk(n_points);
EXPECT_EQ(pcd->points_.size(), n_points);
EXPECT_EQ(pcd->colors_.size(), n_points);
EXPECT_EQ(pcd->normals_.size(), n_points);

for (size_t pidx = 0; pidx < n_points; ++pidx) {
ExpectEQ(pcd->colors_[pidx], Eigen::Vector3d(1, 0, 0));
ExpectEQ(pcd->normals_[pidx], Eigen::Vector3d(0, 1, 0));
}

// With triangle normals
pcd = mesh_simple.SamplePointsPoissonDisk(n_points, 5, nullptr, true);
EXPECT_EQ(pcd->points_.size(), n_points);
EXPECT_EQ(pcd->normals_.size(), n_points);

for (size_t pidx = 0; pidx < n_points; ++pidx) {
ExpectEQ(pcd->normals_[pidx], Eigen::Vector3d(0, 0, 1));
}
}

TEST(TriangleMesh, FilterSharpen) {
auto mesh = std::make_shared<geometry::TriangleMesh>();
mesh->vertices_ = {{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {-1, 0, 0}, {0, -1, 0}};
Expand Down
Loading