From ddb3cbeb3a7c7f5bd91c6828c6c50ca9cefd4098 Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Sun, 7 Jul 2024 08:37:43 -0700 Subject: [PATCH] Remove fancy iterators (#852) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * removed zip from collider * removed zip from boolean result * removed zip from boolean3 * refactored constructors * refactored mesh_fixes * refactored impl * ≈removed reduce_by_key * more zip removal * refactored properties * refactored shared * refactored smoothing * refactored sort * refactored subdivision * refactored polygon * removed zip completely * cleanup * transform iterator * strided range * format * fix pointer/reference to strided range * avoid thrust copy * fix iter * correct location --------- Co-authored-by: pca006132 --- src/collider/src/collider.cpp | 35 ++-- src/manifold/src/boolean3.cpp | 95 +++++------ src/manifold/src/boolean_result.cpp | 100 ++++++----- src/manifold/src/constructors.cpp | 73 +++----- src/manifold/src/csg_tree.cpp | 37 ++--- src/manifold/src/impl.cpp | 148 ++++++----------- src/manifold/src/manifold.cpp | 24 +-- src/manifold/src/mesh_fixes.h | 25 ++- src/manifold/src/properties.cpp | 45 ++--- src/manifold/src/shared.h | 27 +-- src/manifold/src/smoothing.cpp | 95 +++++------ src/manifold/src/sort.cpp | 205 +++++++++-------------- src/manifold/src/subdivision.cpp | 166 ++++++++----------- src/polygon/src/polygon.cpp | 17 +- src/utilities/include/iters.h | 248 ++++++++++++++++++++++++++++ src/utilities/include/par.h | 29 +++- src/utilities/include/sparse.h | 71 ++++---- src/utilities/include/utils.h | 74 ++------- 18 files changed, 774 insertions(+), 740 deletions(-) create mode 100644 src/utilities/include/iters.h diff --git a/src/collider/src/collider.cpp b/src/collider/src/collider.cpp index 321818b1f..02c1c35bb 100644 --- a/src/collider/src/collider.cpp +++ b/src/collider/src/collider.cpp @@ -154,15 +154,13 @@ struct CreateRadixTree { template struct FindCollisions { + VecView queries; VecView nodeBBox_; VecView> internalChildren_; Recorder recorder; - int RecordCollision(int node, thrust::tuple& query) { - const T& queryObj = thrust::get<0>(query); - const int queryIdx = thrust::get<1>(query); - - bool overlaps = nodeBBox_[node].DoesOverlap(queryObj); + int RecordCollision(int node, const int queryIdx) { + bool overlaps = nodeBBox_[node].DoesOverlap(queries[queryIdx]); if (overlaps && IsLeaf(node)) { int leafIdx = Node2Leaf(node); if (!selfCollision || leafIdx != queryIdx) { @@ -172,14 +170,13 @@ struct FindCollisions { return overlaps && IsInternal(node); // Should traverse into node } - void operator()(thrust::tuple query) { + void operator()(const int queryIdx) { // stack cannot overflow because radix tree has max depth 30 (Morton code) + // 32 (index). int stack[64]; int top = -1; // Depth-first search int node = kRoot; - const int queryIdx = thrust::get<1>(query); // same implies that this query do not have any collision if (recorder.earlyexit(queryIdx)) return; while (1) { @@ -187,8 +184,8 @@ struct FindCollisions { int child1 = internalChildren_[internal].first; int child2 = internalChildren_[internal].second; - int traverse1 = RecordCollision(child1, query); - int traverse2 = RecordCollision(child2, query); + int traverse1 = RecordCollision(child1, queryIdx); + int traverse2 = RecordCollision(child2, queryIdx); if (!traverse1 && !traverse2) { if (top < 0) break; // done @@ -313,30 +310,30 @@ SparseIndices Collider::Collisions(const VecView& queriesIn) const { // element can store the sum when using exclusive scan if (queriesIn.size() < kSequentialThreshold) { SparseIndices queryTri; - for_each_n(ExecutionPolicy::Seq, zip(queriesIn.cbegin(), countAt(0)), - queriesIn.size(), + for_each_n(ExecutionPolicy::Seq, countAt(0), queriesIn.size(), FindCollisions>{ - nodeBBox_, internalChildren_, {queryTri}}); + queriesIn, nodeBBox_, internalChildren_, {queryTri}}); return queryTri; } else { // compute the number of collisions to determine the size for allocation and // offset, this avoids the need for atomic Vec counts(queriesIn.size() + 1, 0); Vec empty(queriesIn.size(), 0); - for_each_n(ExecutionPolicy::Par, zip(queriesIn.cbegin(), countAt(0)), - queriesIn.size(), + for_each_n(ExecutionPolicy::Par, countAt(0), queriesIn.size(), FindCollisions{ - nodeBBox_, internalChildren_, {counts, empty}}); + queriesIn, nodeBBox_, internalChildren_, {counts, empty}}); // compute start index for each query and total count exclusive_scan(ExecutionPolicy::Par, counts.begin(), counts.end(), counts.begin(), 0, std::plus()); if (counts.back() == 0) return SparseIndices(0); SparseIndices queryTri(counts.back()); // actually recording collisions - for_each_n(ExecutionPolicy::Par, zip(queriesIn.cbegin(), countAt(0)), - queriesIn.size(), + for_each_n(ExecutionPolicy::Par, countAt(0), queriesIn.size(), FindCollisions>{ - nodeBBox_, internalChildren_, {queryTri, counts, empty}}); + queriesIn, + nodeBBox_, + internalChildren_, + {queryTri, counts, empty}}); return queryTri; } } @@ -350,7 +347,7 @@ void Collider::UpdateBoxes(const VecView& leafBB) { DEBUG_ASSERT(leafBB.size() == NumLeaves(), userErr, "must have the same number of updated boxes as original"); // copy in leaf node Boxes - strided_range::Iter> leaves(nodeBBox_.begin(), nodeBBox_.end(), 2); + auto leaves = StridedRange(nodeBBox_.begin(), nodeBBox_.end(), 2); auto policy = autoPolicy(NumInternal()); copy(policy, leafBB.cbegin(), leafBB.cend(), leaves.begin()); // create global counters diff --git a/src/manifold/src/boolean3.cpp b/src/manifold/src/boolean3.cpp index 27e442950..e80cd0ec6 100644 --- a/src/manifold/src/boolean3.cpp +++ b/src/manifold/src/boolean3.cpp @@ -68,14 +68,13 @@ glm::vec4 Intersect(const glm::vec3 &pL, const glm::vec3 &pR, template struct CopyFaceEdges { const SparseIndices &p1q1; - // const int *p1q1; // x can be either vert or edge (0 or 1). SparseIndices &pXq1; VecView halfedgesQ; + const size_t offset; - void operator()(thrust::tuple in) { - int idx = 3 * thrust::get<0>(in); - size_t i = thrust::get<1>(in); + void operator()(const size_t i) { + int idx = 3 * (i + offset); int pX = p1q1.Get(i, inverted); int q2 = p1q1.Get(i, !inverted); @@ -94,10 +93,10 @@ SparseIndices Filter11(const Manifold::Impl &inP, const Manifold::Impl &inQ, const SparseIndices &p1q2, const SparseIndices &p2q1) { ZoneScoped; SparseIndices p1q1(3 * p1q2.size() + 3 * p2q1.size()); - for_each_n(autoPolicy(p1q2.size()), zip(countAt(0_z), countAt(0_z)), - p1q2.size(), CopyFaceEdges({p1q2, p1q1, inQ.halfedge_})); - for_each_n(autoPolicy(p2q1.size()), zip(countAt(p1q2.size()), countAt(0_z)), - p2q1.size(), CopyFaceEdges({p2q1, p1q1, inP.halfedge_})); + for_each_n(autoPolicy(p1q2.size()), countAt(0_z), p1q2.size(), + CopyFaceEdges({p1q2, p1q1, inQ.halfedge_, 0_z})); + for_each_n(autoPolicy(p2q1.size()), countAt(0_z), p2q1.size(), + CopyFaceEdges({p2q1, p1q1, inP.halfedge_, p1q2.size()})); p1q1.Unique(); return p1q1; } @@ -177,19 +176,21 @@ size_t monobound_quaternary_search(VecView array, int64_t key) { } struct Kernel11 { + VecView xyzz; + VecView s; VecView vertPosP; VecView vertPosQ; VecView halfedgeP; VecView halfedgeQ; - float expandP; + const float expandP; VecView normalP; const SparseIndices &p1q1; - void operator()(thrust::tuple inout) { - const int p1 = p1q1.Get(thrust::get<0>(inout), false); - const int q1 = p1q1.Get(thrust::get<0>(inout), true); - glm::vec4 &xyzz11 = thrust::get<1>(inout); - int &s11 = thrust::get<2>(inout); + void operator()(const size_t idx) { + const int p1 = p1q1.Get(idx, false); + const int q1 = p1q1.Get(idx, true); + glm::vec4 &xyzz11 = xyzz[idx]; + int &s11 = s[idx]; // For pRL[k], qRL[k], k==0 is the left and k==1 is the right. int k = 0; @@ -262,10 +263,9 @@ std::tuple, Vec> Shadow11(SparseIndices &p1q1, Vec s11(p1q1.size()); Vec xyzz11(p1q1.size()); - for_each_n(autoPolicy(p1q1.size()), - zip(countAt(0_z), xyzz11.begin(), s11.begin()), p1q1.size(), - Kernel11({inP.vertPos_, inQ.vertPos_, inP.halfedge_, inQ.halfedge_, - expandP, inP.vertNormal_, p1q1})); + for_each_n(autoPolicy(p1q1.size()), countAt(0_z), p1q1.size(), + Kernel11({xyzz11, s11, inP.vertPos_, inQ.vertPos_, inP.halfedge_, + inQ.halfedge_, expandP, inP.vertNormal_, p1q1})); p1q1.KeepFinite(xyzz11, s11); @@ -273,6 +273,8 @@ std::tuple, Vec> Shadow11(SparseIndices &p1q1, }; struct Kernel02 { + VecView s; + VecView z; VecView vertPosP; VecView halfedgeQ; VecView vertPosQ; @@ -281,11 +283,11 @@ struct Kernel02 { const SparseIndices &p0q2; const bool forward; - void operator()(thrust::tuple inout) { - const int p0 = p0q2.Get(thrust::get<0>(inout), !forward); - const int q2 = p0q2.Get(thrust::get<0>(inout), forward); - int &s02 = thrust::get<1>(inout); - float &z02 = thrust::get<2>(inout); + void operator()(const size_t idx) { + const int p0 = p0q2.Get(idx, !forward); + const int q2 = p0q2.Get(idx, forward); + int &s02 = s[idx]; + float &z02 = z[idx]; // For yzzLR[k], k==0 is the left and k==1 is the right. int k = 0; @@ -353,10 +355,9 @@ std::tuple, Vec> Shadow02(const Manifold::Impl &inP, Vec z02(p0q2.size()); auto vertNormalP = forward ? inP.vertNormal_ : inQ.vertNormal_; - for_each_n(autoPolicy(p0q2.size()), - zip(countAt(0_z), s02.begin(), z02.begin()), p0q2.size(), - Kernel02({inP.vertPos_, inQ.halfedge_, inQ.vertPos_, expandP, - vertNormalP, p0q2, forward})); + for_each_n(autoPolicy(p0q2.size()), countAt(0_z), p0q2.size(), + Kernel02({s02, z02, inP.vertPos_, inQ.halfedge_, inQ.vertPos_, + expandP, vertNormalP, p0q2, forward})); p0q2.KeepFinite(z02, s02); @@ -364,6 +365,8 @@ std::tuple, Vec> Shadow02(const Manifold::Impl &inP, }; struct Kernel12 { + VecView x; + VecView v; VecView p0q2; VecView s02; VecView z02; @@ -376,11 +379,11 @@ struct Kernel12 { const bool forward; const SparseIndices &p1q2; - void operator()(thrust::tuple inout) { - int p1 = p1q2.Get(thrust::get<0>(inout), !forward); - int q2 = p1q2.Get(thrust::get<0>(inout), forward); - int &x12 = thrust::get<1>(inout); - glm::vec3 &v12 = thrust::get<2>(inout); + void operator()(const size_t idx) { + int p1 = p1q2.Get(idx, !forward); + int q2 = p1q2.Get(idx, forward); + int &x12 = x[idx]; + glm::vec3 &v12 = v[idx]; // For xzyLR-[k], k==0 is the left and k==1 is the right. int k = 0; @@ -460,9 +463,8 @@ std::tuple, Vec> Intersect12( Vec v12(p1q2.size()); for_each_n( - autoPolicy(p1q2.size()), zip(countAt(0_z), x12.begin(), v12.begin()), - p1q2.size(), - Kernel12({p0q2.AsVec64(), s02, z02, p1q1.AsVec64(), s11, xyzz11, + autoPolicy(p1q2.size()), countAt(0_z), p1q2.size(), + Kernel12({x12, v12, p0q2.AsVec64(), s02, z02, p1q1.AsVec64(), s11, xyzz11, inP.halfedge_, inQ.halfedge_, inP.vertPos_, forward, p1q2})); p1q2.KeepFinite(v12, x12); @@ -475,26 +477,11 @@ Vec Winding03(const Manifold::Impl &inP, Vec &vertices, Vec &s02, ZoneScoped; // verts that are not shadowed (not in p0q2) have winding number zero. Vec w03(inP.NumVert(), 0); - // checking is slow, so just sort and reduce auto policy = autoPolicy(vertices.size()); - stable_sort( - policy, zip(vertices.begin(), s02.begin()), - zip(vertices.end(), s02.end()), - [](const thrust::tuple &a, const thrust::tuple &b) { - return thrust::get<0>(a) < thrust::get<0>(b); - }); - Vec w03val(w03.size()); - Vec w03vert(w03.size()); - // sum known s02 values into w03 (winding number) - auto endPair = reduce_by_key< - thrust::pair>( - policy, vertices.begin(), vertices.end(), s02.begin(), w03vert.begin(), - w03val.begin()); - scatter(policy, w03val.begin(), endPair.second, w03vert.begin(), w03.begin()); - - if (reverse) - transform(policy, w03.begin(), w03.end(), w03.begin(), - thrust::negate()); + for_each_n(policy, countAt(0), s02.size(), + [&w03, &vertices, &s02, reverse](const int i) { + AtomicAdd(w03[vertices[i]], s02[i] * (reverse ? -1 : 1)); + }); return w03; }; } // namespace diff --git a/src/manifold/src/boolean_result.cpp b/src/manifold/src/boolean_result.cpp index 9afc6150d..3f6760538 100644 --- a/src/manifold/src/boolean_result.cpp +++ b/src/manifold/src/boolean_result.cpp @@ -52,14 +52,14 @@ struct AbsSum : public thrust::binary_function { struct DuplicateVerts { VecView vertPosR; + VecView inclusion; + VecView vertR; + VecView vertPosP; - void operator()(thrust::tuple in) { - int inclusion = abs(thrust::get<0>(in)); - int vertR = thrust::get<1>(in); - glm::vec3 vertPosP = thrust::get<2>(in); - - for (int i = 0; i < inclusion; ++i) { - vertPosR[vertR + i] = vertPosP; + void operator()(const int vert) { + const int n = glm::abs(inclusion[vert]); + for (int i = 0; i < n; ++i) { + vertPosR[vertR[vert] + i] = vertPosP[vert]; } } }; @@ -77,13 +77,14 @@ template struct CountNewVerts { VecView countP; VecView countQ; + VecView i12; const SparseIndices &pq; VecView halfedges; - void operator()(thrust::tuple in) { - int edgeP = pq.Get(thrust::get<0>(in), inverted); - int faceQ = pq.Get(thrust::get<0>(in), !inverted); - int inclusion = glm::abs(thrust::get<1>(in)); + void operator()(const int idx) { + int edgeP = pq.Get(idx, inverted); + int faceQ = pq.Get(idx, !inverted); + int inclusion = glm::abs(i12[idx]); AtomicAdd(countQ[faceQ], inclusion); const Halfedge half = halfedges[edgeP]; @@ -112,16 +113,17 @@ std::tuple, Vec> SizeOutput( inP.halfedge_.end(), CountVerts({sidesPerFaceP, i03})); for_each(autoPolicy(inP.halfedge_.size()), inQ.halfedge_.begin(), inQ.halfedge_.end(), CountVerts({sidesPerFaceQ, i30})); - for_each_n(autoPolicy(i12.size()), zip(countAt(0), i12.begin()), i12.size(), + + for_each_n(autoPolicy(i12.size()), countAt(0), i12.size(), CountNewVerts( - {sidesPerFaceP, sidesPerFaceQ, p1q2, inP.halfedge_})); - for_each_n( - autoPolicy(i21.size()), zip(countAt(0), i21.begin()), i21.size(), - CountNewVerts({sidesPerFaceQ, sidesPerFaceP, p2q1, inQ.halfedge_})); + {sidesPerFaceP, sidesPerFaceQ, i12, p1q2, inP.halfedge_})); + for_each_n(autoPolicy(i21.size()), countAt(0), i21.size(), + CountNewVerts( + {sidesPerFaceQ, sidesPerFaceP, i21, p2q1, inQ.halfedge_})); Vec facePQ2R(inP.NumTri() + inQ.NumTri() + 1, 0); - auto keepFace = - thrust::make_transform_iterator(sidesPerFacePQ.begin(), NotZero()); + auto keepFace = TransformIterator(sidesPerFacePQ.begin(), NotZero()); + inclusive_scan(autoPolicy(sidesPerFacePQ.size()), keepFace, keepFace + sidesPerFacePQ.size(), facePQ2R.begin() + 1); int numFaceR = facePQ2R.back(); @@ -133,10 +135,10 @@ std::tuple, Vec> SizeOutput( inP.faceNormal_.end(), keepFace, outR.faceNormal_.begin(), thrust::identity()); if (invertQ) { - auto start = thrust::make_transform_iterator(inQ.faceNormal_.begin(), - thrust::negate()); - auto end = thrust::make_transform_iterator(inQ.faceNormal_.end(), - thrust::negate()); + auto start = + TransformIterator(inQ.faceNormal_.begin(), thrust::negate()); + auto end = + TransformIterator(inQ.faceNormal_.end(), thrust::negate()); copy_if( autoPolicy(inQ.faceNormal_.size()), start, end, keepFace + inP.NumTri(), next, thrust::identity()); @@ -394,15 +396,16 @@ struct DuplicateHalfedges { VecView halfedgesR; VecView halfedgeRef; VecView facePtr; + VecView wholeHalfedgeP; VecView halfedgesP; VecView i03; VecView vP2R; VecView faceP2R; const bool forward; - void operator()(thrust::tuple in) { - if (!thrust::get<0>(in)) return; - Halfedge halfedge = thrust::get<1>(in); + void operator()(const int idx) { + if (!wholeHalfedgeP[idx]) return; + Halfedge halfedge = halfedgesP[idx]; if (!halfedge.IsForward()) return; const int inclusion = i03[halfedge.startVert]; @@ -447,12 +450,12 @@ void AppendWholeEdges(Manifold::Impl &outR, Vec &facePtrR, const Vec &vP2R, VecView faceP2R, bool forward) { ZoneScoped; - for_each_n(ManifoldParams().deterministic ? ExecutionPolicy::Seq - : autoPolicy(inP.halfedge_.size()), - zip(wholeHalfedgeP.begin(), inP.halfedge_.begin(), countAt(0)), - inP.halfedge_.size(), - DuplicateHalfedges({outR.halfedge_, halfedgeRef, facePtrR, - inP.halfedge_, i03, vP2R, faceP2R, forward})); + for_each_n( + ManifoldParams().deterministic ? ExecutionPolicy::Seq + : autoPolicy(inP.halfedge_.size()), + countAt(0), inP.halfedge_.size(), + DuplicateHalfedges({outR.halfedge_, halfedgeRef, facePtrR, wholeHalfedgeP, + inP.halfedge_, i03, vP2R, faceP2R, forward})); } struct MapTriRef { @@ -488,6 +491,7 @@ void UpdateReference(Manifold::Impl &outR, const Manifold::Impl &inP, struct Barycentric { VecView uvw; + VecView ref; VecView vertPosP; VecView vertPosQ; VecView vertPosR; @@ -496,9 +500,8 @@ struct Barycentric { VecView halfedgeR; const float precision; - void operator()(thrust::tuple in) { - const int tri = thrust::get<0>(in); - const TriRef refPQ = thrust::get<1>(in); + void operator()(const int tri) { + const TriRef refPQ = ref[tri]; if (halfedgeR[3 * tri].startVert < 0) return; const int triPQ = refPQ.tri; @@ -530,11 +533,10 @@ void CreateProperties(Manifold::Impl &outR, const Manifold::Impl &inP, outR.meshRelation_.triProperties.resize(numTri); Vec bary(outR.halfedge_.size()); - for_each_n(autoPolicy(numTri), - zip(countAt(0), outR.meshRelation_.triRef.cbegin()), numTri, - Barycentric({bary, inP.vertPos_, inQ.vertPos_, outR.vertPos_, - inP.halfedge_, inQ.halfedge_, outR.halfedge_, - outR.precision_})); + for_each_n(autoPolicy(numTri), countAt(0), numTri, + Barycentric({bary, outR.meshRelation_.triRef, inP.vertPos_, + inQ.vertPos_, outR.vertPos_, inP.halfedge_, + inQ.halfedge_, outR.halfedge_, outR.precision_})); using Entry = std::pair; int idMissProp = outR.NumVert(); @@ -706,19 +708,15 @@ Manifold::Impl Boolean3::Result(OpType op) const { outR.vertPos_.resize(numVertR); // Add vertices, duplicating for inclusion numbers not in [-1, 1]. // Retained vertices from P and Q: - for_each_n(autoPolicy(inP_.NumVert()), - zip(i03.begin(), vP2R.begin(), inP_.vertPos_.begin()), - inP_.NumVert(), DuplicateVerts({outR.vertPos_})); - for_each_n(autoPolicy(inQ_.NumVert()), - zip(i30.begin(), vQ2R.begin(), inQ_.vertPos_.begin()), - inQ_.NumVert(), DuplicateVerts({outR.vertPos_})); + for_each_n(autoPolicy(inP_.NumVert()), countAt(0), inP_.NumVert(), + DuplicateVerts({outR.vertPos_, i03, vP2R, inP_.vertPos_})); + for_each_n(autoPolicy(inQ_.NumVert()), countAt(0), inQ_.NumVert(), + DuplicateVerts({outR.vertPos_, i30, vQ2R, inQ_.vertPos_})); // New vertices created from intersections: - for_each_n(autoPolicy(i12.size()), - zip(i12.begin(), v12R.begin(), v12_.begin()), i12.size(), - DuplicateVerts({outR.vertPos_})); - for_each_n(autoPolicy(i21.size()), - zip(i21.begin(), v21R.begin(), v21_.begin()), i21.size(), - DuplicateVerts({outR.vertPos_})); + for_each_n(autoPolicy(i12.size()), countAt(0), i12.size(), + DuplicateVerts({outR.vertPos_, i12, v12R, v12_})); + for_each_n(autoPolicy(i21.size()), countAt(0), i21.size(), + DuplicateVerts({outR.vertPos_, i21, v21R, v21_})); PRINT(nPv << " verts from inP"); PRINT(nQv << " verts from inQ"); diff --git a/src/manifold/src/constructors.cpp b/src/manifold/src/constructors.cpp index 92548060a..7d048e395 100644 --- a/src/manifold/src/constructors.cpp +++ b/src/manifold/src/constructors.cpp @@ -12,42 +12,11 @@ // See the License for the specific language governing permissions and // limitations under the License. -#include - #include "csg_tree.h" #include "impl.h" #include "par.h" #include "polygon.h" -namespace { -using namespace manifold; -using namespace thrust::placeholders; - -struct ToSphere { - float length; - void operator()(glm::vec3& v) { - v = glm::cos(glm::half_pi() * (1.0f - v)); - v = length * glm::normalize(v); - if (isnan(v.x)) v = glm::vec3(0.0); - } -}; - -struct Equals { - int val; - bool operator()(int x) { return x == val; } -}; - -struct RemoveFace { - VecView halfedge; - VecView vertLabel; - const int keepLabel; - - bool operator()(int face) { - return vertLabel[halfedge[3 * face].startVert] != keepLabel; - } -}; -} // namespace - namespace manifold { /** * Constructs a smooth version of the input mesh by creating tangents; this @@ -213,7 +182,11 @@ Manifold Manifold::Sphere(float radius, int circularSegments) { auto pImpl_ = std::make_shared(Impl::Shape::Octahedron); pImpl_->Subdivide([n](glm::vec3 edge) { return n - 1; }); for_each_n(autoPolicy(pImpl_->NumVert()), pImpl_->vertPos_.begin(), - pImpl_->NumVert(), ToSphere({radius})); + pImpl_->NumVert(), [radius](glm::vec3& v) { + v = glm::cos(glm::half_pi() * (1.0f - v)); + v = radius * glm::normalize(v); + if (isnan(v.x)) v = glm::vec3(0.0); + }); pImpl_->Finish(); // Ignore preceding octahedron. pImpl_->InitializeOriginal(); @@ -497,30 +470,34 @@ std::vector Manifold::Decompose() const { } Vec vertLabel(componentIndices); + const int numVert = NumVert(); + auto policy = autoPolicy(numVert); std::vector meshes; for (int i = 0; i < numComponents; ++i) { auto impl = std::make_shared(); // inherit original object's precision impl->precision_ = pImpl_->precision_; - impl->vertPos_.resize(NumVert()); - Vec vertNew2Old(NumVert()); - auto policy = autoPolicy(NumVert()); - auto start = zip(impl->vertPos_.begin(), vertNew2Old.begin()); - int nVert = - copy_if( - policy, zip(pImpl_->vertPos_.begin(), countAt(0)), - zip(pImpl_->vertPos_.end(), countAt(NumVert())), vertLabel.begin(), - zip(impl->vertPos_.begin(), vertNew2Old.begin()), Equals({i})) - - start; + + Vec vertNew2Old(numVert); + const int nVert = + copy_if( + policy, countAt(0), countAt(numVert), vertNew2Old.begin(), + [i, &vertLabel](int v) { return vertLabel[v] == i; }) - + vertNew2Old.begin(); impl->vertPos_.resize(nVert); + vertNew2Old.resize(nVert); + gather(policy, vertNew2Old.begin(), vertNew2Old.end(), + pImpl_->vertPos_.begin(), impl->vertPos_.begin()); Vec faceNew2Old(NumTri()); - sequence(policy, faceNew2Old.begin(), faceNew2Old.end()); - - int nFace = remove_if( - policy, faceNew2Old.begin(), faceNew2Old.end(), - RemoveFace({pImpl_->halfedge_, vertLabel, i})) - - faceNew2Old.begin(); + const auto& halfedge = pImpl_->halfedge_; + const int nFace = + copy_if( + policy, countAt(0), countAt(NumTri()), faceNew2Old.begin(), + [i, &vertLabel, &halfedge](int face) { + return vertLabel[halfedge[3 * face].startVert] == i; + }) - + faceNew2Old.begin(); faceNew2Old.resize(nFace); impl->GatherFaces(*pImpl_, faceNew2Old); diff --git a/src/manifold/src/csg_tree.cpp b/src/manifold/src/csg_tree.cpp index 781896644..a6055c9b6 100644 --- a/src/manifold/src/csg_tree.cpp +++ b/src/manifold/src/csg_tree.cpp @@ -32,9 +32,9 @@ constexpr int kParallelThreshold = 4096; namespace { using namespace manifold; struct Transform4x3 { - const glm::mat4x3 transform; + glm::mat4x3 transform; - glm::vec3 operator()(glm::vec3 position) { + glm::vec3 operator()(glm::vec3 position) const { return transform * glm::vec4(position, 1.0f); } }; @@ -253,9 +253,9 @@ Manifold::Impl CsgLeafNode::Compose( auto &oldProp = node->pImpl_->meshRelation_.properties; auto &newProp = combined.meshRelation_.properties; for (int p = 0; p < numProp; ++p) { - strided_range::IterC> oldRange(oldProp.begin() + p, - oldProp.end(), numProp); - strided_range::Iter> newRange( + auto oldRange = + StridedRange(oldProp.cbegin() + p, oldProp.cend(), numProp); + auto newRange = StridedRange( newProp.begin() + numPropOut * propVertIndices[i] + p, newProp.end(), numPropOut); copy(policy, oldRange.begin(), oldRange.end(), newRange.begin()); @@ -277,31 +277,28 @@ Manifold::Impl CsgLeafNode::Compose( } else { // no need to apply the transform to the node, just copy the vertices // and face normals and apply transform on the fly - auto vertPosBegin = thrust::make_transform_iterator( + auto vertPosBegin = TransformIterator( node->pImpl_->vertPos_.begin(), Transform4x3({node->transform_})); glm::mat3 normalTransform = glm::inverse(glm::transpose(glm::mat3(node->transform_))); - auto faceNormalBegin = thrust::make_transform_iterator( - node->pImpl_->faceNormal_.begin(), - TransformNormals({normalTransform})); + auto faceNormalBegin = + TransformIterator(node->pImpl_->faceNormal_.begin(), + TransformNormals({normalTransform})); copy_n(policy, vertPosBegin, node->pImpl_->vertPos_.size(), combined.vertPos_.begin() + vertIndices[i]); copy_n(policy, faceNormalBegin, node->pImpl_->faceNormal_.size(), combined.faceNormal_.begin() + triIndices[i]); const bool invert = glm::determinant(glm::mat3(node->transform_)) < 0; - for_each_n(policy, - zip(combined.halfedgeTangent_.begin() + edgeIndices[i], - countAt(0)), - node->pImpl_->halfedgeTangent_.size(), - TransformTangents{glm::mat3(node->transform_), invert, - node->pImpl_->halfedgeTangent_, - node->pImpl_->halfedge_}); + for_each_n( + policy, countAt(0), node->pImpl_->halfedgeTangent_.size(), + TransformTangents{combined.halfedgeTangent_, edgeIndices[i], + glm::mat3(node->transform_), invert, + node->pImpl_->halfedgeTangent_, + node->pImpl_->halfedge_}); if (invert) - for_each_n(policy, - zip(combined.meshRelation_.triRef.begin(), - countAt(triIndices[i])), - node->pImpl_->NumTri(), FlipTris({combined.halfedge_})); + for_each_n(policy, countAt(triIndices[i]), node->pImpl_->NumTri(), + FlipTris({combined.halfedge_})); } // Since the nodes may be copies containing the same meshIDs, it is // important to add an offset so that each node instance gets diff --git a/src/manifold/src/impl.cpp b/src/manifold/src/impl.cpp index 4f9fbbd83..e1e891b4d 100644 --- a/src/manifold/src/impl.cpp +++ b/src/manifold/src/impl.cpp @@ -35,15 +35,11 @@ void AtomicAddVec3(glm::vec3& target, const glm::vec3& add) { std::atomic& tar = reinterpret_cast&>(target[i]); float old_val = tar.load(std::memory_order_relaxed); while (!tar.compare_exchange_weak(old_val, old_val + add[i], - std::memory_order_relaxed)) - ; + std::memory_order_relaxed)) { + } } } -struct Normalize { - void operator()(glm::vec3& v) { v = SafeNormalize(v); } -}; - struct Transform4x3 { const glm::mat4x3 transform; @@ -53,15 +49,15 @@ struct Transform4x3 { }; struct AssignNormals { + VecView faceNormal; VecView vertNormal; VecView vertPos; VecView halfedges; const float precision; const bool calculateTriNormal; - void operator()(thrust::tuple in) { - glm::vec3& triNormal = thrust::get<0>(in); - const int face = thrust::get<1>(in); + void operator()(const int face) { + glm::vec3& triNormal = faceNormal[face]; glm::ivec3 triVerts; for (int i : {0, 1, 2}) triVerts[i] = halfedges[3 * face + i].startVert; @@ -92,26 +88,6 @@ struct AssignNormals { } }; -struct Tri2Halfedges { - VecView halfedges; - VecView edges; - - void operator()(thrust::tuple in) { - const int tri = thrust::get<0>(in); - const glm::ivec3& triVerts = thrust::get<1>(in); - for (const int i : {0, 1, 2}) { - const int j = (i + 1) % 3; - const int edge = 3 * tri + i; - halfedges[edge] = {triVerts[i], triVerts[j], -1, tri}; - // Sort the forward halfedges in front of the backward ones by setting the - // highest-order bit. - edges[edge] = glm::uint64_t(triVerts[i] < triVerts[j] ? 1 : 0) << 63 | - ((glm::uint64_t)glm::min(triVerts[i], triVerts[j])) << 32 | - glm::max(triVerts[i], triVerts[j]); - } - } -}; - struct LinkHalfedges { VecView halfedges; VecView ids; @@ -146,20 +122,6 @@ struct ReindexTriVerts { } }; -struct InitializeTriRef { - const int meshID; - VecView halfedge; - - void operator()(thrust::tuple inOut) { - TriRef& baryRef = thrust::get<0>(inOut); - int tri = thrust::get<1>(inOut); - - baryRef.meshID = meshID; - baryRef.originalID = meshID; - baryRef.tri = tri; - } -}; - struct UpdateMeshID { const HashTableD meshIDold2new; @@ -167,6 +129,8 @@ struct UpdateMeshID { }; struct CoplanarEdge { + VecView> face2face; + VecView> vert2vert; VecView triArea; VecView halfedge; VecView vertPos; @@ -178,13 +142,7 @@ struct CoplanarEdge { const float precision; // FIXME: race condition - void operator()( - thrust::tuple&, thrust::pair&, int> - inOut) { - thrust::pair& face2face = thrust::get<0>(inOut); - thrust::pair& vert2vert = thrust::get<1>(inOut); - const int edgeIdx = thrust::get<2>(inOut); - + void operator()(const int edgeIdx) { const Halfedge edge = halfedge[edgeIdx]; const Halfedge pair = halfedge[edge.pairedHalfedge]; @@ -206,8 +164,7 @@ struct CoplanarEdge { } } if (propEqual) { - vert2vert.first = prop0; - vert2vert.second = prop1; + vert2vert[edgeIdx] = thrust::make_pair(prop0, prop1); } } @@ -264,8 +221,7 @@ struct CoplanarEdge { } } - face2face.first = edge.face; - face2face.second = pair.face; + face2face[edgeIdx] = thrust::make_pair(edge.face, pair.face); } }; @@ -300,15 +256,6 @@ struct CheckCoplanarity { } }; -struct EdgeBox { - VecView vertPos; - - void operator()(thrust::tuple inout) { - const TmpEdge& edge = thrust::get<1>(inout); - thrust::get<0>(inout) = Box(vertPos[edge.first], vertPos[edge.second]); - } -}; - int GetLabels(std::vector& components, const Vec>& edges, int numNodes) { UnionFind<> uf(numNodes); @@ -609,10 +556,12 @@ void Manifold::Impl::InitializeOriginal() { const int meshID = meshRelation_.originalID; // Don't initialize if it's not an original if (meshID < 0) return; - meshRelation_.triRef.resize(NumTri()); - for_each_n(autoPolicy(NumTri()), - zip(meshRelation_.triRef.begin(), countAt(0)), NumTri(), - InitializeTriRef({meshID, halfedge_})); + auto& triRef = meshRelation_.triRef; + triRef.resize(NumTri()); + for_each_n(autoPolicy(NumTri()), countAt(0), NumTri(), + [meshID, &triRef](const int tri) { + triRef[tri] = {meshID, meshID, tri}; + }); meshRelation_.meshIDtransform.clear(); meshRelation_.meshIDtransform[meshID] = {meshID}; } @@ -626,12 +575,11 @@ void Manifold::Impl::CreateFaces(const std::vector& propertyTolerance) { Vec> face2face(halfedge_.size(), {-1, -1}); Vec> vert2vert(halfedge_.size(), {-1, -1}); Vec triArea(NumTri()); - for_each_n( - autoPolicy(halfedge_.size()), - zip(face2face.begin(), vert2vert.begin(), countAt(0)), halfedge_.size(), - CoplanarEdge({triArea, halfedge_, vertPos_, meshRelation_.triRef, - meshRelation_.triProperties, meshRelation_.properties, - propertyToleranceD, meshRelation_.numProp, precision_})); + for_each_n(autoPolicy(halfedge_.size()), countAt(0), halfedge_.size(), + CoplanarEdge({face2face, vert2vert, triArea, halfedge_, vertPos_, + meshRelation_.triRef, meshRelation_.triProperties, + meshRelation_.properties, propertyToleranceD, + meshRelation_.numProp, precision_})); if (meshRelation_.triProperties.size() > 0) { DedupePropVerts(meshRelation_.triProperties, vert2vert); @@ -677,19 +625,29 @@ void Manifold::Impl::CreateHalfedges(const Vec& triVerts) { Vec ids(numHalfedge); auto policy = autoPolicy(numTri); sequence(policy, ids.begin(), ids.end()); - for_each_n(policy, zip(countAt(0), triVerts.begin()), numTri, - Tri2Halfedges({halfedge_, edge})); + for_each_n(policy, countAt(0), numTri, + [this, &edge, &triVerts](const int tri) { + const glm::ivec3& verts = triVerts[tri]; + for (const int i : {0, 1, 2}) { + const int j = (i + 1) % 3; + const int e = 3 * tri + i; + halfedge_[e] = {verts[i], verts[j], -1, tri}; + // Sort the forward halfedges in front of the backward ones by + // setting the highest-order bit. + edge[e] = glm::uint64_t(verts[i] < verts[j] ? 1 : 0) << 63 | + ((glm::uint64_t)glm::min(verts[i], verts[j])) << 32 | + glm::max(verts[i], verts[j]); + } + }); // Stable sort is required here so that halfedges from the same face are // paired together (the triangles were created in face order). In some // degenerate situations the triangulator can add the same internal edge in // two different faces, causing this edge to not be 2-manifold. These are // fixed by duplicating verts in SimplifyTopology. - stable_sort(policy, zip(edge.begin(), ids.begin()), - zip(edge.end(), ids.end()), - [](const thrust::tuple& a, - const thrust::tuple& b) { - return thrust::get<0>(a) < thrust::get<0>(b); - }); + stable_sort( + policy, ids.begin(), ids.end(), + [&edge](const int& a, const int& b) { return edge[a] < edge[b]; }); + // Once sorted, the first half of the range is the forward halfedges, which // correspond to their backward pair at the same offset in the second half // of the range. @@ -774,15 +732,15 @@ Manifold::Impl Manifold::Impl::Transform(const glm::mat4x3& transform_) const { const bool invert = glm::determinant(glm::mat3(transform_)) < 0; if (halfedgeTangent_.size() > 0) { - for_each_n(policy, zip(result.halfedgeTangent_.begin(), countAt(0)), - halfedgeTangent_.size(), - TransformTangents({glm::mat3(transform_), invert, - halfedgeTangent_, halfedge_})); + for_each_n( + policy, countAt(0), halfedgeTangent_.size(), + TransformTangents({result.halfedgeTangent_, 0, glm::mat3(transform_), + invert, halfedgeTangent_, halfedge_})); } if (invert) { - for_each_n(policy, zip(result.meshRelation_.triRef.begin(), countAt(0)), - result.NumTri(), FlipTris({result.halfedge_})); + for_each_n(policy, countAt(0), result.NumTri(), + FlipTris({result.halfedge_})); } // This optimization does a cheap collider update if the transform is @@ -827,10 +785,11 @@ void Manifold::Impl::CalculateNormals() { faceNormal_.resize(NumTri()); calculateTriNormal = true; } - for_each_n(policy, zip(faceNormal_.begin(), countAt(0)), NumTri(), - AssignNormals({vertNormal_, vertPos_, halfedge_, precision_, - calculateTriNormal})); - for_each(policy, vertNormal_.begin(), vertNormal_.end(), Normalize()); + for_each_n(policy, countAt(0), NumTri(), + AssignNormals({faceNormal_, vertNormal_, vertPos_, halfedge_, + precision_, calculateTriNormal})); + for_each(policy, vertNormal_.begin(), vertNormal_.end(), + [](glm::vec3& v) { v = SafeNormalize(v); }); } /** @@ -865,9 +824,12 @@ SparseIndices Manifold::Impl::EdgeCollisions(const Impl& Q, Vec edges = CreateTmpEdges(Q.halfedge_); const int numEdge = edges.size(); Vec QedgeBB(numEdge); + const auto& vertPos = Q.vertPos_; auto policy = autoPolicy(numEdge); - for_each_n(policy, zip(QedgeBB.begin(), edges.cbegin()), numEdge, - EdgeBox({Q.vertPos_})); + for_each_n( + policy, countAt(0), numEdge, [&QedgeBB, &edges, &vertPos](const int e) { + QedgeBB[e] = Box(vertPos[edges[e].first], vertPos[edges[e].second]); + }); SparseIndices q1p2(0); if (inverted) diff --git a/src/manifold/src/manifold.cpp b/src/manifold/src/manifold.cpp index 26623f27f..e24aef2e2 100644 --- a/src/manifold/src/manifold.cpp +++ b/src/manifold/src/manifold.cpp @@ -29,19 +29,6 @@ using namespace thrust::placeholders; ExecutionParams manifoldParams; -struct MakeTri { - VecView halfedges; - - void operator()(thrust::tuple inOut) { - glm::ivec3& tri = thrust::get<0>(inOut); - const int face = 3 * thrust::get<1>(inOut); - - for (int i : {0, 1, 2}) { - tri[i] = halfedges[face + i].startVert; - } - } -}; - struct UpdateProperties { float* properties; const int numProp; @@ -170,9 +157,14 @@ Mesh Manifold::GetMesh() const { impl.halfedgeTangent_.end()); result.triVerts.resize(NumTri()); - // note that `triVerts` is `std::vector`, so we cannot use thrust::device - thrust::for_each_n(thrust::host, zip(result.triVerts.begin(), countAt(0)), - NumTri(), MakeTri({impl.halfedge_})); + auto& triVerts = result.triVerts; + const auto& halfedges = impl.halfedge_; + for_each_n(autoPolicy(NumTri()), countAt(0), NumTri(), + [&triVerts, &halfedges](const int tri) { + for (int i : {0, 1, 2}) { + triVerts[tri][i] = halfedges[3 * tri + i].startVert; + } + }); return result; } diff --git a/src/manifold/src/mesh_fixes.h b/src/manifold/src/mesh_fixes.h index d8cb4de89..00a2782df 100644 --- a/src/manifold/src/mesh_fixes.h +++ b/src/manifold/src/mesh_fixes.h @@ -10,9 +10,9 @@ inline int FlipHalfedge(int halfedge) { } struct TransformNormals { - const glm::mat3 transform; + glm::mat3 transform; - glm::vec3 operator()(glm::vec3 normal) { + glm::vec3 operator()(glm::vec3 normal) const { normal = glm::normalize(transform * normal); if (isnan(normal.x)) normal = glm::vec3(0.0f); return normal; @@ -20,30 +20,25 @@ struct TransformNormals { }; struct TransformTangents { + VecView tangent; + const int edgeOffset; const glm::mat3 transform; const bool invert; VecView oldTangents; VecView halfedge; - void operator()(thrust::tuple inOut) { - glm::vec4& tangent = thrust::get<0>(inOut); - int edge = thrust::get<1>(inOut); - if (invert) { - edge = halfedge[FlipHalfedge(edge)].pairedHalfedge; - } - - tangent = glm::vec4(transform * glm::vec3(oldTangents[edge]), - oldTangents[edge].w); + void operator()(const int edgeOut) { + const int edgeIn = + invert ? halfedge[FlipHalfedge(edgeOut)].pairedHalfedge : edgeOut; + tangent[edgeOut + edgeOffset] = glm::vec4( + transform * glm::vec3(oldTangents[edgeIn]), oldTangents[edgeIn].w); } }; struct FlipTris { VecView halfedge; - void operator()(thrust::tuple inOut) { - TriRef& bary = thrust::get<0>(inOut); - const int tri = thrust::get<1>(inOut); - + void operator()(const int tri) { thrust::swap(halfedge[3 * tri], halfedge[3 * tri + 2]); for (const int i : {0, 1, 2}) { diff --git a/src/manifold/src/properties.cpp b/src/manifold/src/properties.cpp index e77a2c737..6bc62fe5e 100644 --- a/src/manifold/src/properties.cpp +++ b/src/manifold/src/properties.cpp @@ -136,19 +136,8 @@ struct CurvatureAngles { } }; -struct NormalizeCurvature { - void operator()(thrust::tuple inOut) { - float& meanCurvature = thrust::get<0>(inOut); - float& gaussianCurvature = thrust::get<1>(inOut); - float area = thrust::get<2>(inOut); - float degree = thrust::get<3>(inOut); - float factor = degree / (6 * area); - meanCurvature *= factor; - gaussianCurvature *= factor; - } -}; - struct UpdateProperties { + VecView triProp; VecView properties; VecView oldProperties; @@ -161,16 +150,13 @@ struct UpdateProperties { const int meanIdx; // FIXME: race condition - void operator()(thrust::tuple inOut) { - glm::ivec3& triProp = thrust::get<0>(inOut); - const auto tri = thrust::get<1>(inOut); - + void operator()(const size_t tri) { for (const int i : {0, 1, 2}) { const int vert = halfedge[3 * tri + i].startVert; if (oldNumProp == 0) { - triProp[i] = vert; + triProp[tri][i] = vert; } - const int propVert = triProp[i]; + const int propVert = triProp[tri][i]; for (int p = 0; p < oldNumProp; ++p) { properties[numProp * propVert + p] = @@ -352,10 +338,13 @@ void Manifold::Impl::CalculateCurvature(int gaussianIdx, int meanIdx) { for_each(policy, countAt(0_z), countAt(NumTri()), CurvatureAngles({vertMeanCurvature, vertGaussianCurvature, vertArea, degree, halfedge_, vertPos_, faceNormal_})); - for_each_n(policy, - zip(vertMeanCurvature.begin(), vertGaussianCurvature.begin(), - vertArea.begin(), degree.begin()), - NumVert(), NormalizeCurvature()); + for_each_n(policy, countAt(0), NumVert(), + [&vertMeanCurvature, &vertGaussianCurvature, &vertArea, + °ree](const int vert) { + const float factor = degree[vert] / (6 * vertArea[vert]); + vertMeanCurvature[vert] *= factor; + vertGaussianCurvature[vert] *= factor; + }); const int oldNumProp = NumProp(); const int numProp = glm::max(oldNumProp, glm::max(gaussianIdx, meanIdx) + 1); @@ -367,10 +356,11 @@ void Manifold::Impl::CalculateCurvature(int gaussianIdx, int meanIdx) { } for_each_n( - policy, zip(meshRelation_.triProperties.begin(), countAt(0_z)), NumTri(), - UpdateProperties({meshRelation_.properties, oldProperties, halfedge_, - vertMeanCurvature, vertGaussianCurvature, oldNumProp, - numProp, gaussianIdx, meanIdx})); + policy, countAt(0_z), NumTri(), + UpdateProperties({meshRelation_.triProperties, meshRelation_.properties, + oldProperties, halfedge_, vertMeanCurvature, + vertGaussianCurvature, oldNumProp, numProp, gaussianIdx, + meanIdx})); CreateFaces(); Finish(); @@ -439,8 +429,7 @@ float Manifold::Impl::MinGap(const Manifold::Impl& other, SparseIndices collisions = collider_.Collisions(faceBoxOther.cview()); float minDistanceSquared = transform_reduce( - autoPolicy(collisions.size()), thrust::counting_iterator(0), - thrust::counting_iterator(collisions.size()), + autoPolicy(collisions.size()), countAt(0_z), countAt(collisions.size()), [&collisions, this, &other](int i) { const int tri = collisions.Get(i, 1); const int triOther = collisions.Get(i, 0); diff --git a/src/manifold/src/shared.h b/src/manifold/src/shared.h index 6a330bcc8..e43a79198 100644 --- a/src/manifold/src/shared.h +++ b/src/manifold/src/shared.h @@ -173,28 +173,19 @@ struct TmpEdge { }; /** @} */ -struct Halfedge2Tmp { - void operator()(thrust::tuple inout) { - const Halfedge& halfedge = thrust::get<1>(inout); - int idx = thrust::get<2>(inout); - if (!halfedge.IsForward()) idx = -1; - - thrust::get<0>(inout) = TmpEdge(halfedge.startVert, halfedge.endVert, idx); - } -}; - -struct TmpInvalid { - bool operator()(const TmpEdge& edge) { return edge.halfedgeIdx < 0; } -}; - Vec inline CreateTmpEdges(const Vec& halfedge) { Vec edges(halfedge.size()); - for_each_n(autoPolicy(edges.size()), - zip(edges.begin(), halfedge.begin(), countAt(0)), edges.size(), - Halfedge2Tmp()); + for_each_n(autoPolicy(edges.size()), countAt(0), edges.size(), + [&edges, &halfedge](const int idx) { + const Halfedge& half = halfedge[idx]; + edges[idx] = TmpEdge(half.startVert, half.endVert, + half.IsForward() ? idx : -1); + }); + size_t numEdge = remove_if( - autoPolicy(edges.size()), edges.begin(), edges.end(), TmpInvalid()) - + autoPolicy(edges.size()), edges.begin(), edges.end(), + [](const TmpEdge& edge) { return edge.halfedgeIdx < 0; }) - edges.begin(); DEBUG_ASSERT(numEdge == halfedge.size() / 2, topologyErr, "Not oriented!"); edges.resize(numEdge); diff --git a/src/manifold/src/smoothing.cpp b/src/manifold/src/smoothing.cpp index 4a3ff3263..adcd3e310 100644 --- a/src/manifold/src/smoothing.cpp +++ b/src/manifold/src/smoothing.cpp @@ -61,25 +61,9 @@ glm::vec4 CircularTangent(const glm::vec3& tangent, const glm::vec3& edgeVec) { return glm::vec4(glm::vec3(bz3) / bz3.w, bz3.w); } -struct SmoothBezier { - const Manifold::Impl* impl; - VecView vertNormal; - - void operator()(thrust::tuple inOut) { - glm::vec4& tangent = thrust::get<0>(inOut); - const Halfedge edge = thrust::get<1>(inOut); - const int edgeIdx = thrust::get<2>(inOut); - - if (impl->IsInsideQuad(edgeIdx)) { - tangent = glm::vec4(0, 0, 0, -1); - return; - } - - tangent = impl->TangentFromNormal(vertNormal[edge.startVert], edgeIdx); - } -}; - struct InterpTri { + VecView vertPos; + VecView vertBary; const Manifold::Impl* impl; static glm::vec4 Homogeneous(glm::vec4 v) { @@ -200,10 +184,10 @@ struct InterpTri { return BezierPoint(bez, y); } - void operator()(thrust::tuple inOut) const { - glm::vec3& pos = thrust::get<0>(inOut); - const int tri = thrust::get<1>(inOut).tri; - const glm::vec4 uvw = thrust::get<1>(inOut).uvw; + void operator()(const int vert) { + glm::vec3& pos = vertPos[vert]; + const int tri = vertBary[vert].tri; + const glm::vec4 uvw = vertBary[vert].uvw; const glm::ivec4 halfedges = impl->GetHalfedges(tri); const glm::mat4x3 corners = { @@ -419,13 +403,10 @@ Vec Manifold::Impl::VertFlatFace(const Vec& flatFaces) const { Vec Manifold::Impl::VertHalfedge() const { Vec vertHalfedge(NumVert()); - for_each_n(autoPolicy(halfedge_.size()), zip(countAt(0), halfedge_.begin()), - halfedge_.size(), - [&vertHalfedge](thrust::tuple in) { - const int idx = thrust::get<0>(in); - const Halfedge halfedge = thrust::get<1>(in); + for_each_n(autoPolicy(halfedge_.size()), countAt(0), halfedge_.size(), + [&vertHalfedge, this](const int idx) { // arbitrary, last one wins. - vertHalfedge[halfedge.startVert] = idx; + vertHalfedge[halfedge_[idx].startVert] = idx; }); return vertHalfedge; } @@ -674,31 +655,28 @@ void Manifold::Impl::SetNormals(int normalIdx, float minSharpAngle) { */ void Manifold::Impl::LinearizeFlatTangents() { const int n = halfedgeTangent_.size(); - for_each_n( - autoPolicy(n), zip(halfedgeTangent_.begin(), countAt(0)), n, - [this](thrust::tuple inOut) { - glm::vec4& tangent = thrust::get<0>(inOut); - const int halfedge = thrust::get<1>(inOut); - glm::vec4& otherTangent = - halfedgeTangent_[halfedge_[halfedge].pairedHalfedge]; - - const glm::bvec2 flat(tangent.w == 0, otherTangent.w == 0); - if (!halfedge_[halfedge].IsForward() || (!flat[0] && !flat[1])) { - return; - } + for_each_n(autoPolicy(n), countAt(0), n, [this](const int halfedge) { + glm::vec4& tangent = halfedgeTangent_[halfedge]; + glm::vec4& otherTangent = + halfedgeTangent_[halfedge_[halfedge].pairedHalfedge]; - const glm::vec3 edgeVec = vertPos_[halfedge_[halfedge].endVert] - - vertPos_[halfedge_[halfedge].startVert]; + const glm::bvec2 flat(tangent.w == 0, otherTangent.w == 0); + if (!halfedge_[halfedge].IsForward() || (!flat[0] && !flat[1])) { + return; + } - if (flat[0] && flat[1]) { - tangent = glm::vec4(edgeVec / 3.0f, 1); - otherTangent = glm::vec4(-edgeVec / 3.0f, 1); - } else if (flat[0]) { - tangent = glm::vec4((edgeVec + glm::vec3(otherTangent)) / 2.0f, 1); - } else { - otherTangent = glm::vec4((-edgeVec + glm::vec3(tangent)) / 2.0f, 1); - } - }); + const glm::vec3 edgeVec = vertPos_[halfedge_[halfedge].endVert] - + vertPos_[halfedge_[halfedge].startVert]; + + if (flat[0] && flat[1]) { + tangent = glm::vec4(edgeVec / 3.0f, 1); + otherTangent = glm::vec4(-edgeVec / 3.0f, 1); + } else if (flat[0]) { + tangent = glm::vec4((edgeVec + glm::vec3(otherTangent)) / 2.0f, 1); + } else { + otherTangent = glm::vec4((-edgeVec + glm::vec3(tangent)) / 2.0f, 1); + } + }); } /** @@ -905,9 +883,14 @@ void Manifold::Impl::CreateTangents(std::vector sharpenedEdges) { } } - for_each_n(autoPolicy(numHalfedge), - zip(tangent.begin(), halfedge_.cbegin(), countAt(0)), numHalfedge, - SmoothBezier({this, vertNormal})); + for_each_n(autoPolicy(numHalfedge), countAt(0), numHalfedge, + [&tangent, &vertNormal, this](const int edgeIdx) { + tangent[edgeIdx] = + IsInsideQuad(edgeIdx) + ? glm::vec4(0, 0, 0, -1) + : TangentFromNormal( + vertNormal[halfedge_[edgeIdx].startVert], edgeIdx); + }); halfedgeTangent_.swap(tangent); @@ -1023,8 +1006,8 @@ void Manifold::Impl::Refine(std::function edgeDivisions) { if (vertBary.size() == 0) return; if (old.halfedgeTangent_.size() == old.halfedge_.size()) { - for_each_n(autoPolicy(NumTri()), zip(vertPos_.begin(), vertBary.begin()), - NumVert(), InterpTri({&old})); + for_each_n(autoPolicy(NumTri()), countAt(0), NumVert(), + InterpTri({vertPos_, vertBary, &old})); // Make original since the subdivided faces have been warped into // being non-coplanar, and hence not being related to the original faces. meshRelation_.originalID = ReserveIDs(1); diff --git a/src/manifold/src/sort.cpp b/src/manifold/src/sort.cpp index 8bb91887f..2e156d178 100644 --- a/src/manifold/src/sort.cpp +++ b/src/manifold/src/sort.cpp @@ -56,46 +56,6 @@ uint32_t MortonCode(glm::vec3 position, Box bBox) { return Collider::MortonCode(position, bBox); } -struct Morton { - const Box bBox; - - void operator()(thrust::tuple inout) { - glm::vec3 position = thrust::get<1>(inout); - thrust::get<0>(inout) = MortonCode(position, bBox); - } -}; - -struct FaceMortonBox { - VecView halfedge; - VecView vertPos; - const Box bBox; - - void operator()(thrust::tuple inout) { - uint32_t& mortonCode = thrust::get<0>(inout); - Box& faceBox = thrust::get<1>(inout); - int face = thrust::get<2>(inout); - - // Removed tris are marked by all halfedges having pairedHalfedge = -1, and - // this will sort them to the end (the Morton code only uses the first 30 of - // 32 bits). - if (halfedge[3 * face].pairedHalfedge < 0) { - mortonCode = kNoCode; - return; - } - - glm::vec3 center(0.0f); - - for (const int i : {0, 1, 2}) { - const glm::vec3 pos = vertPos[halfedge[3 * face + i].startVert]; - center += pos; - faceBox.Union(pos); - } - center /= 3; - - mortonCode = MortonCode(center, bBox); - } -}; - struct Reindex { VecView indexInv; @@ -117,22 +77,6 @@ struct MarkProp { } }; -struct GatherProps { - VecView properties; - VecView oldProperties; - const int numProp; - - void operator()(thrust::tuple in) { - const int oldIdx = thrust::get<0>(in); - const int newIdx = thrust::get<1>(in); - const int keep = thrust::get<2>(in); - if (keep == 0) return; - for (int p = 0; p < numProp; ++p) { - properties[newIdx * numProp + p] = oldProperties[oldIdx * numProp + p]; - } - } -}; - struct ReindexProps { VecView old2new; @@ -143,17 +87,6 @@ struct ReindexProps { } }; -template -void Permute(Vec& inOut, const Vec& new2Old) { - Vec tmp(std::move(inOut)); - inOut.resize(new2Old.size()); - gather(autoPolicy(new2Old.size()), new2Old.begin(), new2Old.end(), - tmp.begin(), inOut.begin()); -} - -template void Permute(Vec&, const Vec&); -template void Permute(Vec&, const Vec&); - struct ReindexFace { VecView halfedge; VecView halfedgeTangent; @@ -180,28 +113,6 @@ struct ReindexFace { } }; -struct VertMortonBox { - VecView vertProperties; - const uint32_t numProp; - const float tol; - const Box bBox; - - void operator()(thrust::tuple inout) { - uint32_t& mortonCode = thrust::get<0>(inout); - Box& vertBox = thrust::get<1>(inout); - int vert = thrust::get<2>(inout); - - const glm::vec3 center(vertProperties[numProp * vert], - vertProperties[numProp * vert + 1], - vertProperties[numProp * vert + 2]); - - vertBox.min = center - tol / 2; - vertBox.max = center + tol / 2; - - mortonCode = MortonCode(center, bBox); - } -}; - struct Duplicate { thrust::pair operator()(float x) { return thrust::make_pair(x, x); @@ -293,26 +204,28 @@ void Manifold::Impl::SortVerts() { const auto numVert = NumVert(); Vec vertMorton(numVert); auto policy = autoPolicy(numVert); - for_each_n(policy, zip(vertMorton.begin(), vertPos_.cbegin()), numVert, - Morton({bBox_})); + for_each_n(policy, countAt(0), numVert, [this, &vertMorton](const int vert) { + vertMorton[vert] = MortonCode(vertPos_[vert], bBox_); + }); Vec vertNew2Old(numVert); sequence(policy, vertNew2Old.begin(), vertNew2Old.end()); - stable_sort(policy, zip(vertMorton.begin(), vertNew2Old.begin()), - zip(vertMorton.end(), vertNew2Old.end()), - [](const thrust::tuple& a, - const thrust::tuple& b) { - return thrust::get<0>(a) < thrust::get<0>(b); + stable_sort(policy, vertNew2Old.begin(), vertNew2Old.end(), + [&vertMorton](const int& a, const int& b) { + return vertMorton[a] < vertMorton[b]; }); ReindexVerts(vertNew2Old, numVert); // Verts were flagged for removal with NaNs and assigned kNoCode to sort // them to the end, which allows them to be removed. - const int newNumVert = - std::lower_bound(vertMorton.begin(), vertMorton.end(), kNoCode) - - vertMorton.begin(); + const int newNumVert = find_if( + policy, vertNew2Old.begin(), vertNew2Old.end(), + [&vertMorton](const int vert) { + return vertMorton[vert] == kNoCode; + }) - + vertNew2Old.begin(); vertNew2Old.resize(newNumVert); Permute(vertPos_, vertNew2Old); @@ -355,10 +268,18 @@ void Manifold::Impl::CompactProps() { Vec oldProp = meshRelation_.properties; const int numVertsNew = propOld2New[numVerts]; - meshRelation_.properties.resize(meshRelation_.numProp * numVertsNew); + const int numProp = meshRelation_.numProp; + auto& properties = meshRelation_.properties; + properties.resize(numProp * numVertsNew); for_each_n( - policy, zip(countAt(0), propOld2New.cbegin(), keep.cbegin()), numVerts, - GatherProps({meshRelation_.properties, oldProp, meshRelation_.numProp})); + policy, countAt(0), numVerts, + [&properties, &oldProp, &propOld2New, &keep, &numProp](const int oldIdx) { + if (keep[oldIdx] == 0) return; + for (int p = 0; p < numProp; ++p) { + properties[propOld2New[oldIdx] * numProp + p] = + oldProp[oldIdx * numProp + p]; + } + }); for_each_n(policy, meshRelation_.triProperties.begin(), NumTri(), ReindexProps({propOld2New})); } @@ -373,9 +294,28 @@ void Manifold::Impl::GetFaceBoxMorton(Vec& faceBox, ZoneScoped; faceBox.resize(NumTri()); faceMorton.resize(NumTri()); - for_each_n(autoPolicy(NumTri()), - zip(faceMorton.begin(), faceBox.begin(), countAt(0)), NumTri(), - FaceMortonBox({halfedge_, vertPos_, bBox_})); + for_each_n(autoPolicy(NumTri()), countAt(0), NumTri(), + [this, &faceBox, &faceMorton](const int face) { + // Removed tris are marked by all halfedges having pairedHalfedge + // = -1, and this will sort them to the end (the Morton code only + // uses the first 30 of 32 bits). + if (halfedge_[3 * face].pairedHalfedge < 0) { + faceMorton[face] = kNoCode; + return; + } + + glm::vec3 center(0.0f); + + for (const int i : {0, 1, 2}) { + const glm::vec3 pos = + vertPos_[halfedge_[3 * face + i].startVert]; + center += pos; + faceBox[face].Union(pos); + } + center /= 3; + + faceMorton[face] = MortonCode(center, bBox_); + }); } /** @@ -388,22 +328,22 @@ void Manifold::Impl::SortFaces(Vec& faceBox, Vec& faceMorton) { auto policy = autoPolicy(faceNew2Old.size()); sequence(policy, faceNew2Old.begin(), faceNew2Old.end()); - stable_sort(policy, zip(faceMorton.begin(), faceNew2Old.begin()), - zip(faceMorton.end(), faceNew2Old.end()), - [](const thrust::tuple& a, - const thrust::tuple& b) { - return thrust::get<0>(a) < thrust::get<0>(b); + stable_sort(policy, faceNew2Old.begin(), faceNew2Old.end(), + [&faceMorton](const int& a, const int& b) { + return faceMorton[a] < faceMorton[b]; }); // Tris were flagged for removal with pairedHalfedge = -1 and assigned kNoCode // to sort them to the end, which allows them to be removed. - const int newNumTri = - find(policy, faceMorton.begin(), - faceMorton.end(), kNoCode) - - faceMorton.begin(); - faceMorton.resize(newNumTri); + const int newNumTri = find_if( + policy, faceNew2Old.begin(), faceNew2Old.end(), + [&faceMorton](const int face) { + return faceMorton[face] == kNoCode; + }) - + faceNew2Old.begin(); faceNew2Old.resize(newNumTri); + Permute(faceMorton, faceNew2Old); Permute(faceBox, faceNew2Old); GatherFaces(faceNew2Old); } @@ -549,8 +489,7 @@ bool MeshGL::Merge() { Vec vertPropD(vertProperties); Box bBox; for (const int i : {0, 1, 2}) { - strided_range::Iter> iPos(vertPropD.begin() + i, vertPropD.end(), - numProp); + auto iPos = StridedRange(vertPropD.begin() + i, vertPropD.end(), numProp); auto minMax = transform_reduce>( autoPolicy(numVert), iPos.begin(), iPos.end(), Duplicate(), thrust::make_pair(std::numeric_limits::infinity(), @@ -566,18 +505,32 @@ bool MeshGL::Merge() { Vec vertBox(numOpenVert); Vec vertMorton(numOpenVert); - for_each_n(policy, - zip(vertMorton.begin(), vertBox.begin(), openVerts.cbegin()), - numOpenVert, VertMortonBox({vertPropD, numProp, precision, bBox})); + for_each_n(policy, countAt(0), numOpenVert, + [&vertMorton, &vertBox, &openVerts, &bBox, this](const int i) { + int vert = openVerts[i]; + + const glm::vec3 center(vertProperties[numProp * vert], + vertProperties[numProp * vert + 1], + vertProperties[numProp * vert + 2]); + + vertBox[i].min = center - precision / 2; + vertBox[i].max = center + precision / 2; - stable_sort(policy, - zip(vertMorton.begin(), vertBox.begin(), openVerts.begin()), - zip(vertMorton.end(), vertBox.end(), openVerts.end()), - [](const thrust::tuple& a, - const thrust::tuple& b) { - return thrust::get<0>(a) < thrust::get<0>(b); + vertMorton[i] = MortonCode(center, bBox); + }); + + Vec vertNew2Old(numOpenVert); + sequence(policy, vertNew2Old.begin(), vertNew2Old.end()); + + stable_sort(policy, vertNew2Old.begin(), vertNew2Old.end(), + [&vertMorton](const int& a, const int& b) { + return vertMorton[a] < vertMorton[b]; }); + Permute(vertMorton, vertNew2Old); + Permute(vertBox, vertNew2Old); + Permute(openVerts, vertNew2Old); + Collider collider(vertBox, vertMorton); SparseIndices toMerge = collider.Collisions(vertBox.cview()); diff --git a/src/manifold/src/subdivision.cpp b/src/manifold/src/subdivision.cpp index 07253e34d..2a88469f0 100644 --- a/src/manifold/src/subdivision.cpp +++ b/src/manifold/src/subdivision.cpp @@ -26,19 +26,6 @@ struct std::hash { namespace { using namespace manifold; -struct ReindexHalfedge { - VecView half2Edge; - const VecView halfedges; - - void operator()(thrust::tuple in) { - const int edge = thrust::get<0>(in); - const int halfedge = thrust::get<1>(in).halfedgeIdx; - - half2Edge[halfedge] = edge; - half2Edge[halfedges[halfedge].pairedHalfedge] = edge; - } -}; - class Partition { public: // The cached partitions don't have idx - it's added to the copy returned @@ -131,14 +118,12 @@ class Partition { const int numTri = triVert.size(); Vec newTriVert(numTri); - for_each_n( - autoPolicy(numTri), zip(newTriVert.begin(), triVert.begin()), numTri, - [&outTri, &newVerts](thrust::tuple inOut) { - for (const int j : {0, 1, 2}) { - thrust::get<0>(inOut)[outTri[j]] = - newVerts[thrust::get<1>(inOut)[j]]; - } - }); + for_each_n(autoPolicy(numTri), countAt(0), numTri, + [&newTriVert, &outTri, &newVerts, this](const int tri) { + for (const int j : {0, 1, 2}) { + newTriVert[tri][outTri[j]] = newVerts[triVert[tri][j]]; + } + }); return newTriVert; } @@ -509,29 +494,29 @@ Vec Manifold::Impl::Subdivide( const int numTri = NumTri(); Vec half2Edge(2 * numEdge); auto policy = autoPolicy(numEdge); - for_each_n(policy, zip(countAt(0), edges.begin()), numEdge, - ReindexHalfedge({half2Edge, halfedge_})); + for_each_n(policy, countAt(0), numEdge, + [&half2Edge, &edges, this](const int edge) { + const int idx = edges[edge].halfedgeIdx; + half2Edge[idx] = edge; + half2Edge[halfedge_[idx].pairedHalfedge] = edge; + }); Vec faceHalfedges(numTri); - for_each_n(policy, zip(faceHalfedges.begin(), countAt(0)), numTri, - [this](thrust::tuple inOut) { - glm::ivec4& halfedges = thrust::get<0>(inOut); - const int tri = thrust::get<1>(inOut); - halfedges = GetHalfedges(tri); - }); + for_each_n(policy, countAt(0), numTri, [&faceHalfedges, this](const int tri) { + faceHalfedges[tri] = GetHalfedges(tri); + }); Vec edgeAdded(numEdge); - for_each_n(policy, zip(edgeAdded.begin(), edges.cbegin()), numEdge, - [edgeDivisions, this](thrust::tuple inOut) { - int& divisions = thrust::get<0>(inOut); - const TmpEdge edge = thrust::get<1>(inOut); + for_each_n(policy, countAt(0), numEdge, + [&edgeAdded, &edges, edgeDivisions, this](const int i) { + const TmpEdge edge = edges[i]; if (IsMarkedInsideQuad(edge.halfedgeIdx)) { - divisions = 0; + edgeAdded[i] = 0; return; } const glm::vec3 vec = vertPos_[edge.first] - vertPos_[edge.second]; - divisions = edgeDivisions(vec); + edgeAdded[i] = edgeDivisions(vec); }); Vec edgeOffset(numEdge); @@ -541,13 +526,12 @@ Vec Manifold::Impl::Subdivide( Vec vertBary(edgeOffset.back() + edgeAdded.back()); const int totalEdgeAdded = vertBary.size() - numVert; FillRetainedVerts(vertBary); - for_each_n(policy, zip(edges.begin(), edgeAdded.begin(), edgeOffset.begin()), - numEdge, [this, &vertBary](thrust::tuple in) { - const TmpEdge edge = thrust::get<0>(in); - const int n = thrust::get<1>(in); - const int offset = thrust::get<2>(in); + for_each_n(policy, countAt(0), numEdge, + [&vertBary, &edges, &edgeAdded, &edgeOffset, this](const int i) { + const int n = edgeAdded[i]; + const int offset = edgeOffset[i]; - const BaryIndices indices = GetIndices(edge.halfedgeIdx); + const BaryIndices indices = GetIndices(edges[i].halfedgeIdx); if (indices.tri < 0) { return; // inside quad } @@ -576,14 +560,14 @@ Vec Manifold::Impl::Subdivide( }); Vec triOffset(numTri); - auto numSubTris = thrust::make_transform_iterator( - subTris.begin(), [](const Partition& part) { + auto numSubTris = + TransformIterator(subTris.begin(), [](const Partition& part) { return static_cast(part.triVert.size()); }); exclusive_scan(policy, numSubTris, numSubTris + numTri, triOffset.begin(), 0); Vec interiorOffset(numTri); - auto numInterior = thrust::make_transform_iterator( + auto numInterior = TransformIterator( subTris.begin(), [](const Partition& part) { return part.NumInterior(); }); exclusive_scan(policy, numInterior, numInterior + numTri, @@ -643,27 +627,24 @@ Vec Manifold::Impl::Subdivide( meshRelation_.triRef = triRef; Vec newVertPos(vertBary.size()); - for_each_n( - policy, zip(newVertPos.begin(), vertBary.begin()), vertBary.size(), - [this, &faceHalfedges](thrust::tuple inOut) { - glm::vec3& pos = thrust::get<0>(inOut); - const Barycentric bary = thrust::get<1>(inOut); - - const glm::ivec4 halfedges = faceHalfedges[bary.tri]; - if (halfedges[3] < 0) { - glm::mat3 triPos; - for (const int i : {0, 1, 2}) { - triPos[i] = vertPos_[halfedge_[halfedges[i]].startVert]; - } - pos = triPos * glm::vec3(bary.uvw); - } else { - glm::mat4x3 quadPos; - for (const int i : {0, 1, 2, 3}) { - quadPos[i] = vertPos_[halfedge_[halfedges[i]].startVert]; - } - pos = quadPos * bary.uvw; - } - }); + for_each_n(policy, countAt(0), vertBary.size(), + [&newVertPos, &vertBary, &faceHalfedges, this](const int vert) { + const Barycentric bary = vertBary[vert]; + const glm::ivec4 halfedges = faceHalfedges[bary.tri]; + if (halfedges[3] < 0) { + glm::mat3 triPos; + for (const int i : {0, 1, 2}) { + triPos[i] = vertPos_[halfedge_[halfedges[i]].startVert]; + } + newVertPos[vert] = triPos * glm::vec3(bary.uvw); + } else { + glm::mat4x3 quadPos; + for (const int i : {0, 1, 2, 3}) { + quadPos[i] = vertPos_[halfedge_[halfedges[i]].startVert]; + } + newVertPos[vert] = quadPos * bary.uvw; + } + }); vertPos_ = newVertPos; faceNormal_.resize(0); @@ -681,11 +662,11 @@ Vec Manifold::Impl::Subdivide( // copy interior prop verts and forward edge prop verts for_each_n( - policy, zip(countAt(numPropVert), vertBary.begin() + numVert), - addedVerts, - [this, &prop, &faceHalfedges](thrust::tuple in) { - const int vert = thrust::get<0>(in); - const Barycentric bary = thrust::get<1>(in); + policy, countAt(0), addedVerts, + [&prop, &vertBary, &faceHalfedges, numVert, numPropVert, + this](const int i) { + const int vert = numPropVert + i; + const Barycentric bary = vertBary[numVert + i]; const glm::ivec4 halfedges = faceHalfedges[bary.tri]; auto& rel = meshRelation_; @@ -713,30 +694,29 @@ Vec Manifold::Impl::Subdivide( }); // copy backward edge prop verts - for_each_n( - policy, zip(edges.begin(), edgeAdded.begin(), edgeOffset.begin()), - numEdge, - [this, &prop, propOffset, - addedVerts](thrust::tuple in) { - const TmpEdge edge = thrust::get<0>(in); - const int n = thrust::get<1>(in); - const int offset = thrust::get<2>(in) + propOffset + addedVerts; - auto& rel = meshRelation_; + for_each_n(policy, countAt(0), numEdge, + [this, &prop, &edges, &edgeAdded, &edgeOffset, propOffset, + addedVerts](const int i) { + const int n = edgeAdded[i]; + const int offset = edgeOffset[i] + propOffset + addedVerts; + auto& rel = meshRelation_; - const float frac = 1.0f / (n + 1); - const int halfedgeIdx = halfedge_[edge.halfedgeIdx].pairedHalfedge; - const int v0 = halfedgeIdx % 3; - const int tri = halfedgeIdx / 3; - const int prop0 = rel.triProperties[tri][v0]; - const int prop1 = rel.triProperties[tri][Next3(v0)]; - for (int i = 0; i < n; ++i) { - for (int p = 0; p < rel.numProp; ++p) { - prop[(offset + i) * rel.numProp + p] = glm::mix( - rel.properties[prop0 * rel.numProp + p], - rel.properties[prop1 * rel.numProp + p], (i + 1) * frac); - } - } - }); + const float frac = 1.0f / (n + 1); + const int halfedgeIdx = + halfedge_[edges[i].halfedgeIdx].pairedHalfedge; + const int v0 = halfedgeIdx % 3; + const int tri = halfedgeIdx / 3; + const int prop0 = rel.triProperties[tri][v0]; + const int prop1 = rel.triProperties[tri][Next3(v0)]; + for (int i = 0; i < n; ++i) { + for (int p = 0; p < rel.numProp; ++p) { + prop[(offset + i) * rel.numProp + p] = + glm::mix(rel.properties[prop0 * rel.numProp + p], + rel.properties[prop1 * rel.numProp + p], + (i + 1) * frac); + } + } + }); Vec triProp(triVerts.size()); for_each_n(policy, countAt(0), numTri, diff --git a/src/polygon/src/polygon.cpp b/src/polygon/src/polygon.cpp index 317f4132c..2b931f1e3 100644 --- a/src/polygon/src/polygon.cpp +++ b/src/polygon/src/polygon.cpp @@ -844,13 +844,18 @@ class EarClip { return {Collider(), itr}; } - stable_sort(autoPolicy(itr.size()), - zip(vertMorton.begin(), vertBox.begin(), itr.begin()), - zip(vertMorton.end(), vertBox.end(), itr.end()), - [](const thrust::tuple &a, - const thrust::tuple &b) { - return thrust::get<0>(a) < thrust::get<0>(b); + const int numVert = itr.size(); + auto policy = autoPolicy(numVert); + Vec vertNew2Old(numVert); + sequence(policy, vertNew2Old.begin(), vertNew2Old.end()); + + stable_sort(policy, 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}; } diff --git a/src/utilities/include/iters.h b/src/utilities/include/iters.h new file mode 100644 index 000000000..7f7408a39 --- /dev/null +++ b/src/utilities/include/iters.h @@ -0,0 +1,248 @@ +// Copyright 2024 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 +#include + +namespace manifold { + +template +struct InnerIter { + using pointer = typename std::iterator_traits::pointer; + using reference = typename std::iterator_traits::reference; + using difference_type = typename std::iterator_traits::difference_type; + using value_type = typename std::iterator_traits::value_type; + using iterator_category = + typename std::iterator_traits::iterator_category; +}; + +template +struct InnerIter>> { + using pointer = Iter; + using reference = std::remove_pointer_t&; + using difference_type = std::ptrdiff_t; + using value_type = std::remove_pointer_t; + using iterator_category = std::random_access_iterator_tag; +}; + +template +struct TransformIterator { + private: + Iter iter; + F f; + + public: + // users are not suppposed to take pointer/reference of the iterator. + using pointer = void; + using reference = void; + using difference_type = typename InnerIter::difference_type; + using value_type = + std::invoke_result_t::value_type>; + using iterator_category = typename InnerIter::iterator_category; + + TransformIterator(Iter iter, F f) : iter(iter), f(f) {} + + value_type operator*() const { return f(*iter); } + + value_type operator[](size_t i) const { return f(iter[i]); } + + // prefix increment + TransformIterator& operator++() { + iter += 1; + return *this; + } + + // postfix + TransformIterator operator++(int) { + auto old = *this; + operator++(); + return old; + } + + TransformIterator operator+(size_t n) const { + return TransformIterator(iter + n, f); + } + + TransformIterator& operator+=(size_t n) { + iter += n; + return *this; + } + + friend bool operator==(TransformIterator a, TransformIterator b) { + return a.iter == b.iter; + } + + friend bool operator!=(TransformIterator a, TransformIterator b) { + return !(a.iter == b.iter); + } + + friend bool operator<(TransformIterator a, TransformIterator b) { + return a.iter < b.iter; + } + + friend difference_type operator-(TransformIterator a, TransformIterator b) { + return a.iter - b.iter; + } + + operator TransformIterator() const { + return TransformIterator(f, iter); + } +}; + +template +struct CountingIterator { + private: + T counter; + + public: + using pointer = void; + using reference = T; + using difference_type = std::make_signed_t; + using value_type = T; + using iterator_category = std::random_access_iterator_tag; + + CountingIterator(T counter) : counter(counter) {} + + value_type operator*() const { return counter; } + value_type operator[](T i) const { return counter + i; } + + // prefix increment + CountingIterator& operator++() { + counter += 1; + return *this; + } + + // postfix + CountingIterator operator++(int) { + auto old = *this; + operator++(); + return old; + } + + CountingIterator operator+(T n) const { + return CountingIterator(counter + n); + } + + CountingIterator& operator+=(T n) { + counter += n; + return *this; + } + + friend bool operator==(CountingIterator a, CountingIterator b) { + return a.counter == b.counter; + } + + friend bool operator!=(CountingIterator a, CountingIterator b) { + return a.counter != b.counter; + } + + friend bool operator<(CountingIterator a, CountingIterator b) { + return a.counter < b.counter; + } + + friend difference_type operator-(CountingIterator a, CountingIterator b) { + return a.counter - b.counter; + } + + operator CountingIterator() const { + return CountingIterator(counter); + } +}; + +template +CountingIterator countAt(T i) { + return CountingIterator(i); +} + +template +struct StridedRange { + private: + struct StridedRangeIter { + private: + Iter iter; + size_t stride; + + public: + using pointer = typename InnerIter::pointer; + using reference = typename InnerIter::reference; + using difference_type = typename InnerIter::difference_type; + using value_type = typename InnerIter::value_type; + using iterator_category = typename InnerIter::iterator_category; + + StridedRangeIter(Iter iter, int stride) : iter(iter), stride(stride) {} + + value_type& operator*() { return *iter; } + + const value_type& operator*() const { return *iter; } + + value_type& operator[](size_t i) { return iter[i * stride]; } + + const value_type& operator[](size_t i) const { return iter[i * stride]; } + + // prefix increment + StridedRangeIter& operator++() { + iter += stride; + return *this; + } + + // postfix + StridedRangeIter operator++(int) { + auto old = *this; + operator++(); + return old; + } + + StridedRangeIter operator+(size_t n) const { + return StridedRangeIter(iter + n * stride, stride); + } + + StridedRangeIter& operator+=(size_t n) { + iter += n * stride; + return *this; + } + + friend bool operator==(StridedRangeIter a, StridedRangeIter b) { + return a.iter == b.iter; + } + + friend bool operator!=(StridedRangeIter a, StridedRangeIter b) { + return !(a.iter == b.iter); + } + + friend bool operator<(StridedRangeIter a, StridedRangeIter b) { + return a.iter < b.iter; + } + + friend difference_type operator-(StridedRangeIter a, StridedRangeIter b) { + // note that this is not well-defined if a.stride != b.stride... + return (a.iter - b.iter) / a.stride; + } + }; + Iter _start, _end; + const size_t stride; + + public: + StridedRange(Iter start, Iter end, size_t stride) + : _start(start), _end(end), stride(stride) {} + + StridedRangeIter begin() const { return StridedRangeIter(_start, stride); } + + StridedRangeIter end() const { + return StridedRangeIter(_start, stride) + + ((std::distance(_start, _end) + (stride - 1)) / stride); + } +}; + +} // namespace manifold diff --git a/src/utilities/include/par.h b/src/utilities/include/par.h index 2049a5011..fb29b48d0 100644 --- a/src/utilities/include/par.h +++ b/src/utilities/include/par.h @@ -34,12 +34,12 @@ #include #endif -#include "tbb/tbb.h" #define MANIFOLD_PAR_NS tbb #else #define MANIFOLD_PAR_NS cpp #endif +#include "iters.h" namespace manifold { enum class ExecutionPolicy { @@ -174,15 +174,34 @@ STL_DYNAMIC_BACKEND_VOID(uninitialized_fill) STL_DYNAMIC_BACKEND_VOID(uninitialized_copy) STL_DYNAMIC_BACKEND_VOID(stable_sort) STL_DYNAMIC_BACKEND_VOID(fill) -STL_DYNAMIC_BACKEND_VOID(copy) STL_DYNAMIC_BACKEND_VOID(inclusive_scan) -STL_DYNAMIC_BACKEND_VOID(copy_n) + +// there are some issues with thrust copy +template +OutputIterator copy(ExecutionPolicy policy, InputIterator first, + InputIterator last, OutputIterator result) { +#if MANIFOLD_PAR == 'T' && \ + (TBB_INTERFACE_VERSION >= 10000 && __has_include()) + if (policy == ExecutionPolicy::Par) + return std::copy(std::execution::par_unseq, first, last, result); +#endif + return std::copy(first, last, result); +} + +template +OutputIterator copy_n(ExecutionPolicy policy, InputIterator first, size_t n, + OutputIterator result) { +#if MANIFOLD_PAR == 'T' && \ + (TBB_INTERFACE_VERSION >= 10000 && __has_include()) + if (policy == ExecutionPolicy::Par) + return std::copy_n(std::execution::par_unseq, first, n, result); +#endif + return std::copy_n(first, n, result); +} // void implies that the user have to specify the return type in the template // argument, as we are unable to deduce it THRUST_DYNAMIC_BACKEND(transform_reduce, void) -THRUST_DYNAMIC_BACKEND(gather_if, void) -THRUST_DYNAMIC_BACKEND(reduce_by_key, void) STL_DYNAMIC_BACKEND(remove, void) STL_DYNAMIC_BACKEND(find, void) STL_DYNAMIC_BACKEND(find_if, void) diff --git a/src/utilities/include/sparse.h b/src/utilities/include/sparse.h index ee048adc2..9bd1fc8e6 100644 --- a/src/utilities/include/sparse.h +++ b/src/utilities/include/sparse.h @@ -135,46 +135,61 @@ class SparseIndices { size_t RemoveZeros(Vec& S) { DEBUG_ASSERT(S.size() == size(), userErr, "Different number of values than indicies!"); - VecView view = AsVec64(); - auto zBegin = zip(S.begin(), view.begin()); - auto zEnd = zip(S.end(), view.end()); - size_t size = - remove_if(autoPolicy(S.size()), zBegin, zEnd, - [](thrust::tuple x) { - return thrust::get<0>(x) == 0; - }) - - zBegin; - S.resize(size, -1); + + Vec new2Old(S.size()); + auto policy = autoPolicy(S.size()); + sequence(policy, new2Old.begin(), new2Old.end()); + + size_t size = copy_if( + policy, countAt(0_z), countAt(S.size()), new2Old.begin(), + [&S](const int i) { return S[i] != 0; }) - + new2Old.begin(); + new2Old.resize(size); + + Permute(S, new2Old); + Vec tmp(std::move(data_)); Resize(size); + gather(policy, new2Old.begin(), new2Old.end(), + reinterpret_cast(tmp.data()), + reinterpret_cast(data_.data())); + return size; } template - struct firstNonFinite { - bool NotFinite(float v) const { return !isfinite(v); } - bool NotFinite(glm::vec2 v) const { return !isfinite(v[0]); } - bool NotFinite(glm::vec3 v) const { return !isfinite(v[0]); } - bool NotFinite(glm::vec4 v) const { return !isfinite(v[0]); } - - bool operator()(thrust::tuple x) const { - bool result = NotFinite(thrust::get<0>(x)); - return result; - } + struct firstFinite { + VecView v; + + bool Finite(float v) const { return isfinite(v); } + bool Finite(glm::vec2 v) const { return isfinite(v[0]); } + bool Finite(glm::vec3 v) const { return isfinite(v[0]); } + bool Finite(glm::vec4 v) const { return isfinite(v[0]); } + + bool operator()(const int i) const { return Finite(v[i]); } }; template size_t KeepFinite(Vec& v, Vec& x) { DEBUG_ASSERT(x.size() == size(), userErr, "Different number of values than indicies!"); - VecView view = AsVec64(); - auto zBegin = zip(v.begin(), x.begin(), view.begin()); - auto zEnd = zip(v.end(), x.end(), view.end()); - size_t size = remove_if(autoPolicy(v.size()), zBegin, - zEnd, firstNonFinite()) - - zBegin; - v.resize(size); - x.resize(size); + + Vec new2Old(v.size()); + auto policy = autoPolicy(v.size()); + + size_t size = copy_if( + policy, countAt(0_z), countAt(v.size()), new2Old.begin(), + firstFinite({v})) - + new2Old.begin(); + new2Old.resize(size); + + Permute(v, new2Old); + Permute(x, new2Old); + Vec tmp(std::move(data_)); Resize(size); + gather(policy, new2Old.begin(), new2Old.end(), + reinterpret_cast(tmp.data()), + reinterpret_cast(data_.data())); + return size; } diff --git a/src/utilities/include/utils.h b/src/utilities/include/utils.h index 1c13cd0d3..e6b3c41b7 100644 --- a/src/utilities/include/utils.h +++ b/src/utilities/include/utils.h @@ -63,21 +63,6 @@ struct Timer { }; #endif -template -thrust::zip_iterator> zip(Iters... iters) { - return thrust::make_zip_iterator(thrust::make_tuple(iters...)); -} - -template -thrust::permutation_iterator perm(A a, B b) { - return thrust::make_permutation_iterator(a, b); -} - -template -thrust::counting_iterator countAt(T i) { - return thrust::make_counting_iterator(i); -} - inline int Next3(int i) { constexpr glm::ivec3 next3(1, 2, 0); return next3[i]; @@ -88,13 +73,21 @@ inline int Prev3(int i) { return prev3[i]; } +template +void Permute(Vec& inOut, const Vec& new2Old) { + Vec tmp(std::move(inOut)); + inOut.resize(new2Old.size()); + gather(autoPolicy(new2Old.size()), new2Old.begin(), new2Old.end(), + tmp.begin(), inOut.begin()); +} + template T AtomicAdd(T& target, T add) { std::atomic& tar = reinterpret_cast&>(target); T old_val = tar.load(); while (!tar.compare_exchange_weak(old_val, old_val + add, - std::memory_order_seq_cst)) - ; + std::memory_order_seq_cst)) { + } return old_val; } @@ -105,53 +98,6 @@ inline int AtomicAdd(int& target, int add) { return old_val; } -// Copied from -// https://github.com/thrust/thrust/blob/master/examples/strided_range.cu -template -class strided_range { - public: - typedef typename thrust::iterator_difference::type difference_type; - - struct stride_functor - : public thrust::unary_function { - difference_type stride; - - stride_functor(difference_type stride) : stride(stride) {} - - difference_type operator()(const difference_type& i) const { - return stride * i; - } - }; - - typedef typename thrust::counting_iterator CountingIterator; - typedef typename thrust::transform_iterator - TransformIterator; - typedef typename thrust::permutation_iterator - PermutationIterator; - - // type of the strided_range iterator - typedef PermutationIterator iterator; - - // construct strided_range for the range [first,last) - strided_range(Iterator first, Iterator last, difference_type stride) - : first(first), last(last), stride(stride) {} - strided_range() {} - - iterator begin(void) const { - return PermutationIterator( - first, TransformIterator(CountingIterator(0), stride_functor(stride))); - } - - iterator end(void) const { - return begin() + ((last - first) + (stride - 1)) / stride; - } - - protected: - Iterator first; - Iterator last; - difference_type stride; -}; - template class ConcurrentSharedPtr { public: