8000 kdtree (2D) for collider by pca006132 · Pull Request #1160 · elalish/manifold · GitHub
[go: up one dir, main page]
More Web Proxy on the site http://driver.im/
Skip to content

kdtree (2D) for collider #1160

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.

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

Merged
merged 8 commits into from
Feb 23, 2025
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
2 changes: 2 additions & 0 deletions src/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@ set(
smoothing.cpp
sort.cpp
subdivision.cpp
tree2d.cpp
# optional source files
$<$<BOOL:${MANIFOLD_CROSS_SECTION}>:cross_section/cross_section.cpp>
$<$<BOOL:${MANIFOLD_EXPORT}>:meshIO/meshIO.cpp>
Expand All @@ -48,6 +49,7 @@ set(
shared.h
sparse.h
svd.h
tree2d.h
tri_dist.h
utils.h
vec.h
Expand Down
77 changes: 26 additions & 51 deletions src/polygon.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -18,8 +18,8 @@
#include <map>
#include <set>

#include "./collider.h"
#include "./parallel.h"
#include "./tree2d.h"
#include "./utils.h"
#include "manifold/optional_assert.h"

Expand Down Expand Up @@ -143,6 +143,7 @@ void Dump(const PolygonsIdx &polys, double epsilon) {

void PrintFailure(const std::exception &e, const PolygonsIdx &polys,
std::vector<ivec3> &triangles, double epsilon) {
std::cout << std::setprecision(16);
std::cout << "-----------------------------------" << std::endl;
std::cout << "Triangulation failed! Precision = " << epsilon << std::endl;
std::cout << e.what() << std::endl;
Expand Down Expand Up @@ -306,9 +307,8 @@ class EarClip {
double epsilon_;

struct IdxCollider {
Collider collider;
Vec<PolyVert> points;
std::vector<VertItr> itr;
SparseIndices ind;
};

// A circularly-linked list representing the polygon(s) that still need to be
Expand Down Expand Up @@ -502,32 +502,26 @@ class EarClip {
return totalCost;
}

Box earBox = Box{vec3(center.x - radius, center.y - radius, 0),
vec3(center.x + radius, center.y + radius, 0)};
earBox.Union(vec3(pos, 0));
collider.collider.Collisions(VecView<const Box>(&earBox, 1),
collider.ind);
Rect earBox = Rect(vec2(center.x - radius, center.y - radius),
vec2(center.x + radius, center.y + radius));
earBox.Union(pos);
earBox.min -= epsilon;
earBox.max += epsilon;
Copy link
Owner

Choose a reason for hiding this comment

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

Good call - surprised this wasn't a problem before.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

We do need it previously, just that I forgot to add it back. Previously, we enlarge the bbox of each vertex, but now we enlarge the bbox of the query.


const int lid = left->mesh_idx;
const int rid = right->mesh_idx;

totalCost = transform_reduce(
countAt(0), countAt(collider.ind.size()), totalCost,
[](double a, double b) { return std::max(a, b); },
[&](size_t i) {
const VertItr test = collider.itr[collider.ind.Get(i, true)];
if (!Clipped(test) && test->mesh_idx != mesh_idx &&
test->mesh_idx != lid &&
test->mesh_idx != rid) { // Skip duplicated verts
double cost = Cost(test, openSide, epsilon);
if (cost < -epsilon) {
cost = DelaunayCost(test->pos - center, scale, epsilon);
}
return cost;
}
return std::numeric_limits<double>::lowest();
Copy link
Owner

Choose a reason for hiding this comment

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

Is this case handled inside the query function now?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I think for lowest, it just means it will not change the cost (because we take the max).

});
collider.ind.Clear();
QueryTwoDTree(collider.points, earBox, [&](PolyVert point) {
const VertItr test = collider.itr[point.idx];
if (!Clipped(test) && test->mesh_idx != mesh_idx &&
test->mesh_idx != lid &&
test->mesh_idx != rid) { // Skip duplicated verts
double cost = Cost(test, openSide, epsilon);
if (cost < -epsilon) {
cost = DelaunayCost(test->pos - center, scale, epsilon);
}
if (cost > totalCost) totalCost = cost;
}
});
return totalCost;
}

Expand Down Expand Up @@ -836,35 +830,16 @@ class EarClip {
// epsilon_. Each ear uses this BVH to quickly find a subset of vertices to
// check for cost.
IdxCollider VertCollider(VertItr start) const {
Vec<Box> vertBox;
Vec<uint32_t> vertMorton;
ZoneScoped;
std::vector<VertItr> itr;
const Box box(vec3(bBox_.min, 0), vec3(bBox_.max, 0));

Loop(start, [&vertBox, &vertMorton, &itr, &box, this](VertItr v) {
Vec<PolyVert> points;
Loop(start, [&itr, &points, this](VertItr v) {
points.push_back({v->pos, static_cast<int>(itr.size())});
itr.push_back(v);
const vec3 pos(v->pos, 0);
vertBox.push_back({pos - epsilon_, pos + epsilon_});
vertMorton.push_back(Collider::MortonCode(pos, box));
});

if (itr.empty()) {
return {Collider(), itr, {}};
}

const int numVert = itr.size();
Vec<int> vertNew2Old(numVert);
sequence(vertNew2Old.begin(), vertNew2Old.end());

stable_sort(vertNew2Old.begin(), vertNew2Old.end(),
[&vertMorton](const int a, const int b) {
return vertMorton[a] < vertMorton[b];
});
Permute(vertMorton, vertNew2Old);
Permute(vertBox, vertNew2Old);
Permute(itr, vertNew2Old);

return {Collider(vertBox, vertMorton), itr, {}};
BuildTwoDTree(points);
return {std::move(points), std::move(itr)};
}

// The main ear-clipping loop. This is called once for each simple polygon -
Expand Down
46 changes: 46 additions & 0 deletions src/tree2d.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,46 @@
// Copyright 2025 The Manifold Authors.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#include "tree2d.h"

#include "./parallel.h"
#include "./utils.h"

namespace manifold {

// Not really a proper KD-tree, but a kd tree with k = 2 and alternating x/y
// partition.
// Recursive sorting is not the most efficient, but simple and guaranteed to
// result in a balanced tree.
void BuildTwoDTreeImpl(VecView<PolyVert> points, bool sortX) {
using CmpFn = std::function<bool(const PolyVert&, const PolyVert&)>;
CmpFn cmpx = [](const PolyVert& a, const PolyVert& b) {

Check warning on line 28 in src/tree2d.cpp

View check run for this annotation

Codecov / codecov/patch

src/tree2d.cpp#L28

Added line #L28 was not covered by tests
return a.pos.x < b.pos.x;
};
CmpFn cmpy = [](const PolyVert& a, const PolyVert& b) {
< 9E88 div class="Details js-details-container position-relative check-annotation check-annotation-warning py-3" id="annotation_32225202321">

Check warning on line 31 in src/tree2d.cpp

View check run for this annotation

Codecov / codecov/patch

src/tree2d.cpp#L31

Added line #L31 was not covered by tests
return a.pos.y < b.pos.y;
};
manifold::stable_sort(points.begin(), points.end(), sortX ? cmpx : cmpy);
if (points.size() < 2) return;
BuildTwoDTreeImpl(points.view(0, points.size() / 2), !sortX);
BuildTwoDTreeImpl(points.view(points.size() / 2 + 1), !sortX);
}

void BuildTwoDTree(VecView<PolyVert> points) {
ZoneScoped;
// don't even bother...
if (points.size() <= 8) return;
BuildTwoDTreeImpl(points, true);
}
} // namespace manifold
85 changes: 85 additions & 0 deletions src/tree2d.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,85 @@
// Copyright 2025 The Manifold Authors.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#pragma once

#include "manifold/common.h"
#include "manifold/optional_assert.h"
#include "manifold/polygon.h"
#include "manifold/vec_view.h"

namespace manifold {

void BuildTwoDTreeImpl(VecView<PolyVert> points, bool sortX);

void BuildTwoDTree(VecView<PolyVert> points);

template <typename F>
void QueryTwoDTree(VecView<PolyVert> points, Rect r, F f) {
if (points.size() <= 8) {
for (const auto& p : points)
if (r.Contains(p.pos)) f(p);
return;
}
Rect current;
current.min = vec2(-std::numeric_limits<double>::infinity());
current.max = vec2(std::numeric_limits<double>::infinity());

int level = 0;
VecView<PolyVert> currentView = points;
std::array<Rect, 64> rectStack;
std::array<VecView<PolyVert>, 64> viewStack;
std::array<int, 64> levelStack;
int stackPointer = 0;

while (1) {
if (currentView.size() <= 2) {
for (const auto& p : currentView)
if (r.Contains(p.pos)) f(p);
if (--stackPointer < 0) break;
level = levelStack[stackPointer];
currentView = viewStack[stackPointer];
current = rectStack[stackPointer];
continue;
}

// these are conceptual left/right trees
Rect left = current;
Rect right = current;
const PolyVert middle = currentView[currentView.size() / 2];
if (level % 2 == 0)
left.max.x = right.min.x = middle.pos.x;
else
left.max.y = right.min.y = middle.pos.y;

if (r.Contains(middle.pos)) f(middle);
if (left.DoesOverlap(r)) {
if (right.DoesOverlap(r)) {
DEBUG_ASSERT(stackPointer < 64, logicErr, "Stack overflow");
rectStack[stackPointer] = right;
viewStack[stackPointer] = currentView.view(currentView.size() / 2 + 1);
levelStack[stackPointer] = level + 1;
stackPointer++;
}
current = left;
currentView = currentView.view(0, currentView.size() / 2);
level++;
} else {
current = right;
currentView = currentView.view(currentView.size() / 2 + 1);
level++;
}
}
}
} // namespace manifold
25 changes: 24 additions & 1 deletion test/polygons/polygon_corpus.txt
Original file line number Diff line number Diff line change
Expand Up @@ -1670,4 +1670,27 @@ Very_Branchy_With_hole 35 -1.0 2
3
-61.065936 242.624693
-59.65022 241.08199
-61.482003 240.77002
-61.482003 240.77002
CondensedMatter64 19 2.30642e-11 1
21
-0.06842195833766856 9.679939633804798
0.01125349940133374 9.680903111806934
0.005139787287448156 9.690408106108382
0.005139787287445186 9.690408106108382
0.005139787287445186 9.690408106108382
0.007822556637645924 9.684282094785807
-0.06842195833766818 9.679939633804798
-0.01205120016400146 9.694378352997491
-0.01205120016400146 9.694378352997491
-0.01205120016400146 9.694378352997491
-0.02461833537677266 9.699006865448419
-0.02461833537677267 9.699006865448419
-0.03661970614604192 9.704881022991231
-0.03661970614604192 9.704881022991231
-0.04416639052247019 9.709589843848839
-0.06842195833766818 9.679939633804798
-0.04938028753596506 9.713063939287069
-0.05370195201180676 9.716422994315348
-0.06842195833766776 9.6799396338048
-0.06842195833766818 9.679939633804798
-0.06842195833766825 9.679939633804798
Loading
0