diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 2dff72d90..131b56336 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -29,6 +29,7 @@ set( smoothing.cpp sort.cpp subdivision.cpp + tree2d.cpp # optional source files $<$:cross_section/cross_section.cpp> $<$:meshIO/meshIO.cpp> @@ -48,6 +49,7 @@ set( shared.h sparse.h svd.h + tree2d.h tri_dist.h utils.h vec.h diff --git a/src/polygon.cpp b/src/polygon.cpp index f092785b8..ea7e1077e 100644 --- a/src/polygon.cpp +++ b/src/polygon.cpp @@ -18,8 +18,8 @@ #include #include -#include "./collider.h" #include "./parallel.h" +#include "./tree2d.h" #include "./utils.h" #include "manifold/optional_assert.h" @@ -143,6 +143,7 @@ void Dump(const PolygonsIdx &polys, double epsilon) { void PrintFailure(const std::exception &e, const PolygonsIdx &polys, std::vector &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; @@ -306,9 +307,8 @@ class EarClip { double epsilon_; struct IdxCollider { - Collider collider; + Vec points; std::vector itr; - SparseIndices ind; }; // A circularly-linked list representing the polygon(s) that still need to be @@ -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(&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; 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::lowest(); - }); - 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; } @@ -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 vertBox; - Vec vertMorton; + ZoneScoped; std::vector itr; - const Box box(vec3(bBox_.min, 0), vec3(bBox_.max, 0)); - - Loop(start, [&vertBox, &vertMorton, &itr, &box, this](VertItr v) { + Vec points; + Loop(start, [&itr, &points, this](VertItr v) { + points.push_back({v->pos, static_cast(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 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 - diff --git a/src/tree2d.cpp b/src/tree2d.cpp new file mode 100644 index 000000000..04a500505 --- /dev/null +++ b/src/tree2d.cpp @@ -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 points, bool sortX) { + using CmpFn = std::function; + CmpFn cmpx = [](const PolyVert& a, const PolyVert& b) { + return a.pos.x < b.pos.x; + }; + CmpFn cmpy = [](const PolyVert& a, const PolyVert& b) { + 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 points) { + ZoneScoped; + // don't even bother... + if (points.size() <= 8) return; + BuildTwoDTreeImpl(points, true); +} +} // namespace manifold diff --git a/src/tree2d.h b/src/tree2d.h new file mode 100644 index 000000000..d16b2ab3e --- /dev/null +++ b/src/tree2d.h @@ -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 points, bool sortX); + +void BuildTwoDTree(VecView points); + +template +void QueryTwoDTree(VecView 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::infinity()); + current.max = vec2(std::numeric_limits::infinity()); + + int level = 0; + VecView currentView = points; + std::array rectStack; + std::array, 64> viewStack; + std::array 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 diff --git a/test/polygons/polygon_corpus.txt b/test/polygons/polygon_corpus.txt index a11a637e1..377e00238 100644 --- a/test/polygons/polygon_corpus.txt +++ b/test/polygons/polygon_corpus.txt @@ -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 \ No newline at end of file +-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