diff --git a/CombineModels/Logic/Utilities.cxx b/CombineModels/Logic/Utilities.cxx index 22e9867..6d4b7ad 100644 --- a/CombineModels/Logic/Utilities.cxx +++ b/CombineModels/Logic/Utilities.cxx @@ -1,5 +1,5 @@ /* -Copyright 2012-2024 Ronald Römer +Copyright 2012-2025 Ronald Römer Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. @@ -14,8 +14,6 @@ See the License for the specific language governing permissions and limitations under the License. */ -#define _USE_MATH_DEFINES // Needed for making M_PI symbol available on Windows - #include "Utilities.h" #include @@ -32,33 +30,20 @@ limitations under the License. double ComputeNormal (vtkPoints *pts, double *n, vtkIdType num, const vtkIdType *poly) { n[0] = 0; n[1] = 0; n[2] = 0; - if (num == 3) { - double pA[3], pB[3], pC[3], a[3], b[3]; - - pts->GetPoint(poly[0], pA); - pts->GetPoint(poly[1], pB); - pts->GetPoint(poly[2], pC); - - vtkMath::Subtract(pB, pA, a); - vtkMath::Subtract(pC, pA, b); - - vtkMath::Cross(a, b, n); - } else { - double pA[3], pB[3]; + double pA[3], pB[3]; - vtkIdType i, a, b; + vtkIdType i, a, b; - for (i = 0; i < num; i++) { - a = poly[i]; - b = poly[i+1 == num ? 0 : i+1]; + for (i = 0; i < num; i++) { + a = poly[i]; + b = poly[i+1 == num ? 0 : i+1]; - pts->GetPoint(a, pA); - pts->GetPoint(b, pB); + pts->GetPoint(a, pA); + pts->GetPoint(b, pB); - n[0] += (pA[1]-pB[1])*(pA[2]+pB[2]); - n[1] += (pA[2]-pB[2])*(pA[0]+pB[0]); - n[2] += (pA[0]-pB[0])*(pA[1]+pB[1]); - } + n[0] += (pA[1]-pB[1])*(pA[2]+pB[2]); + n[1] += (pA[2]-pB[2])*(pA[0]+pB[0]); + n[2] += (pA[0]-pB[0])*(pA[1]+pB[1]); } return vtkMath::Normalize(n); @@ -258,3 +243,87 @@ void GetPolys (const ReferencedPointsType &pts, const IndexedPolysType &indexedP polys.push_back(std::move(newPoly)); } } + +void GetPoly (vtkPoints *pts, vtkIdType num, const vtkIdType *poly, Poly &out) { + vtkIdType i; + + double p[3]; + + for (i = 0; i < num; i++) { + pts->GetPoint(poly[i], p); + + out.emplace_back(p[0], p[1], p[2], poly[i]); + } +} + +void FlattenPoly (const Poly &poly, Poly &out, const Base &base) { + double pt[3], tr[2]; + + vtkIdType i = 0; + + for (auto &p : poly) { + pt[0] = p.x; + pt[1] = p.y; + pt[2] = p.z; + + Transform(pt, tr, base); + + out.emplace_back(tr[0], tr[1], 0, i++); + } +} + +std::shared_ptr GetEdgeProj (const Poly &poly, const Point3d &p) { + Poly::const_iterator itrA, itrB; + + double d, t; + + std::shared_ptr proj; + + for (itrA = poly.begin(); itrA != poly.end(); ++itrA) { + itrB = itrA+1; + + if (itrB == poly.end()) { + itrB = poly.begin(); + } + + ProjOnLine(*itrA, *itrB, p, &d, &t, proj); + + if (d > 0 && d < 1e-5 && t > 0 && t < 1) { + return std::make_shared(itrA->id, itrB->id, *proj, d); + } + + } + + return nullptr; +} + +void ProjOnLine (const Point3d &a, const Point3d &b, const Point3d &p, double *d, double *t, std::shared_ptr &proj) { + double v[3], w[3]; + + double v_ = Point3d::GetVec(a, b, v); + double w_ = Point3d::GetVec(a, p, w); + + double pr = std::min(std::max(vtkMath::Dot(v, w), -1.), 1.); + + *d = std::sin(std::acos(pr))*w_; + + vtkMath::MultiplyScalar(v, std::sqrt(w_*w_-*d**d)); + + proj = std::make_shared(a.x+v[0], a.y+v[1], a.z+v[2]); + + *t = pr*w_/v_; +} + +void ProjOnLine (vtkPolyData *pd, const Pair &line, const Point3d &p, std::shared_ptr &proj) { + double pA[3], pB[3]; + + pd->GetPoint(line.f, pA); + pd->GetPoint(line.g, pB); + + Point3d a(pA[0], pA[1], pA[2]); + Point3d b(pB[0], pB[1], pB[2]); + + double d, t; + + ProjOnLine(a, b, p, &d, &t, proj); +} diff --git a/CombineModels/Logic/Utilities.h b/CombineModels/Logic/Utilities.h index 79ff498..9f794b6 100644 --- a/CombineModels/Logic/Utilities.h +++ b/CombineModels/Logic/Utilities.h @@ -1,5 +1,5 @@ /* -Copyright 2012-2024 Ronald Römer +Copyright 2012-2025 Ronald Römer Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. @@ -24,6 +24,7 @@ limitations under the License. #include #include #include +#include #include #include @@ -48,7 +49,7 @@ class Point3d { vtkIdType id; Point3d () = delete; - Point3d (const double _x, const double _y, const double _z, vtkIdType _id = NOTSET) : x(_x), y(_y), z(_z), id(_id) {} + Point3d (const double x, const double y, const double z, vtkIdType id = NOTSET) : x(x), y(y), z(z), id(id) {} bool operator< (const Point3d &other) const { const long x1 = std::lround(x*1e5), y1 = std::lround(y*1e5), @@ -78,34 +79,41 @@ class Point3d { return vtkMath::Normalize(v); } + + static double GetDist (const Point3d &a, const Point3d &b) { + double dx = b.x-a.x; + double dy = b.y-a.y; + double dz = b.z-a.z; + + return dx*dx+dy*dy+dz*dz; + } }; typedef std::vector IdsType; typedef std::set _IdsType; -#ifndef __VTK_WRAP__ -template::value, bool>::type = true> -class _Pair { +class Pair { public: - T f, g; - _Pair () = delete; - _Pair (T _f, T _g) : f(_f), g(_g) {} - bool operator< (const _Pair &other) const { + vtkIdType f, g; + Pair () = delete; + Pair (vtkIdType f, vtkIdType g) : f(f), g(g) {} + bool operator< (const Pair &other) const { return std::tie(f, g) < std::tie(other.f, other.g); } - bool operator== (const _Pair &other) const { + bool operator== (const Pair &other) const { return f == other.f && g == other.g; } - friend std::ostream& operator<< (std::ostream &out, const _Pair &p) { + vtkIdType operator& (const Pair &other) const { + if (f == other.f || f == other.g) { + return f; + } + return g; + } + friend std::ostream& operator<< (std::ostream &out, const Pair &p) { out << "(" << p.f << ", " << p.g << ")"; return out; } }; -#endif // __VTK_WRAP__ - -// typedef _Pair Pair; -using Pair = _Pair; class Base { public: @@ -126,7 +134,7 @@ std::ostream& operator<< (typename std::enable_if::value, std::o return stream << static_cast::type>(e); } -typedef std::vector Poly; +typedef std::vector Poly, Points; typedef std::vector PolysType; double ComputeNormal (const Poly &poly, double *n); @@ -141,4 +149,21 @@ typedef std::map> ReferencedPoi void GetPolys (const ReferencedPointsType &pts, const IndexedPolysType &indexedPolys, PolysType &polys); -#endif +void GetPoly (vtkPoints *pts, vtkIdType num, const vtkIdType *poly, Poly &out); +void FlattenPoly (const Poly &poly, Poly &out, const Base &base); + +class Proj { +public: + Proj (vtkIdType a, vtkIdType b, const Point3d &proj, double d) : edge(a, b), proj(proj), d(d) {} + + Pair edge; + Point3d proj; + double d; +}; + +std::shared_ptr GetEdgeProj (const Poly &poly, const Point3d &p); + +void ProjOnLine (const Point3d &a, const Point3d &b, const Point3d &p, double *d, double *t, std::shared_ptr &proj); +void ProjOnLine (vtkPolyData *pd, const Pair &line, const Point3d &p, std::shared_ptr &proj); + +#endif \ No newline at end of file diff --git a/CombineModels/Logic/vtkPolyDataBooleanFilter.cxx b/CombineModels/Logic/vtkPolyDataBooleanFilter.cxx index e6e3cb8..4c4bad6 100644 --- a/CombineModels/Logic/vtkPolyDataBooleanFilter.cxx +++ b/CombineModels/Logic/vtkPolyDataBooleanFilter.cxx @@ -1,5 +1,5 @@ /* -Copyright 2012-2024 Ronald Römer +Copyright 2012-2025 Ronald Römer Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. @@ -148,6 +148,9 @@ int vtkPolyDataBooleanFilter::RequestData(vtkInformation *request, vtkInformatio modPdA->DeepCopy(cl->GetOutput(1)); modPdB->DeepCopy(cl->GetOutput(2)); + modPdA->EditableOn(); + modPdB->EditableOn(); + #ifdef DEBUG std::cout << "Exporting contLines.vtk" << std::endl; WriteVTK("contLines.vtk", contLines); @@ -175,16 +178,23 @@ int vtkPolyDataBooleanFilter::RequestData(vtkInformation *request, vtkInformatio auto cells = vtkSmartPointer::New(); + bool hasGap(false); + for (i = 0; i < numPts; i++) { contLines->GetPointCells(i, cells); if (cells->GetNumberOfIds() == 1) { - vtkErrorMacro("Contact ends suddenly."); - return 1; + vtkErrorMacro("Contact ends suddenly at point " << i << "."); + + hasGap = true; } } + if (hasGap) { + return 1; + } + // sichert die OrigCellIds vtkIdTypeArray *origCellIdsA = vtkIdTypeArray::SafeDownCast(modPdA->GetCellData()->GetScalars("OrigCellIds")); @@ -614,7 +624,9 @@ bool vtkPolyDataBooleanFilter::GetPolyStrips (vtkPolyData *pd, vtkIdTypeArray *c // sucht nach gleichen captPts { - std::map>>> collapsed; + // std::map>>> collapsed; + + std::map> collapsed; PolyStripsType::iterator itr; StripPtsType::iterator itr2; @@ -628,39 +640,24 @@ bool vtkPolyDataBooleanFilter::GetPolyStrips (vtkPolyData *pd, vtkIdTypeArray *c StripPt &sp = itr2->second; if (sp.capt & Capt::BOUNDARY) { - collapsed[{sp.cutPt[0], sp.cutPt[1], sp.cutPt[2]}][sp.ind].push_back(sp); - } - } - } + // collapsed[{sp.cutPt[0], sp.cutPt[1], sp.cutPt[2]}][sp.ind].push_back(sp); - for (auto &[pt, map] : collapsed) { - if (map.size() > 1) { - // std::cout << pt << std::endl; + auto inds = collapsed[{sp.cutPt[0], sp.cutPt[1], sp.cutPt[2]}]; - // for (auto &[ind, pts] : map) { - // for (auto &p : pts) { - // std::cout << p << std::endl; - // } - // } - - // if (map.size() == 2 && std::all_of(map.begin(), map.end(), [](const auto &m) { return m.second.size() == 1 && m.second.back().get().capt == Capt::EDGE; })) { - - // for (auto &[ind, pts] : map) { - // StripPt &sp = pts.back(); - - // sp.capt = Capt::NOT; - // sp.edge[0] = sp.edge[1] = NOTSET; - // Cpy(sp.cutPt, sp.pt, 3); - // } - - // } else { - // return true; - // } - - return true; + inds.emplace(sp.ind); + if (inds.size() > 1) { + return true; + } + } } } + + // for (auto &[pt, map] : collapsed) { + // if (map.size() > 1) { + // return true; + // } + // } } for (itr = polyLines.begin(); itr != polyLines.end(); ++itr) { diff --git a/CombineModels/Logic/vtkPolyDataBooleanFilter.h b/CombineModels/Logic/vtkPolyDataBooleanFilter.h index 9410938..b5b5158 100644 --- a/CombineModels/Logic/vtkPolyDataBooleanFilter.h +++ b/CombineModels/Logic/vtkPolyDataBooleanFilter.h @@ -1,5 +1,5 @@ /* -Copyright 2012-2024 Ronald Römer +Copyright 2012-2025 Ronald Römer Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. @@ -33,11 +33,13 @@ limitations under the License. #include "Utilities.h" -#define OPER_NONE 0 -#define OPER_UNION 1 -#define OPER_INTERSECTION 2 -#define OPER_DIFFERENCE 3 -#define OPER_DIFFERENCE2 4 +enum OperMode { + OPER_NONE = 0, + OPER_UNION, + OPER_INTERSECTION, + OPER_DIFFERENCE, + OPER_DIFFERENCE2 +}; enum class Capt { NOT = 1 << 0, diff --git a/CombineModels/Logic/vtkPolyDataContactFilter.cxx b/CombineModels/Logic/vtkPolyDataContactFilter.cxx index dbd7796..bb2d9e7 100644 --- a/CombineModels/Logic/vtkPolyDataContactFilter.cxx +++ b/CombineModels/Logic/vtkPolyDataContactFilter.cxx @@ -1,5 +1,5 @@ /* -Copyright 2012-2024 Ronald Römer +Copyright 2012-2025 Ronald Römer Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. @@ -40,11 +40,18 @@ limitations under the License. #include #include #include +#include +#include #include "vtkPolyDataContactFilter.h" #include "Utilities.h" -#undef DEBUG +// #define _debA 17920 +// #define _debB 1339 + +#if (defined(_debA) && defined(_debB)) +vtkIdType _idA, _idB; +#endif vtkStandardNewMacro(vtkPolyDataContactFilter); @@ -119,6 +126,13 @@ int vtkPolyDataContactFilter::RequestData (vtkInformation *request, vtkInformati CopyPolyData(_pdA, newPdA); CopyPolyData(_pdB, newPdB); + try { + PreventEqualCaptPoints(newPdA, newPdB).Run(); + } catch (const std::runtime_error &e) { + vtkErrorMacro("Cannot prevent equal capture points."); + return 1; + } + GetInvalidEdges(newPdA, edgesA); GetInvalidEdges(newPdB, edgesB); @@ -409,17 +423,25 @@ void vtkPolyDataContactFilter::InterEdgeLine (InterPtsType &interPts, vtkPolyDat double s = vtkMath::Determinant3x3(p, r, v)/(n*n); if (s > -1e-6 && s < l+1e-6) { + double t = vtkMath::Determinant3x3(p, e, v)/(n*n); - End end = End::NONE; + InterPt q(ptA[0]+t*r[0], ptA[1]+t*r[1], ptA[2]+t*r[2], t, idA, idB, End::NONE, src, PointSrc::CALCULATED); if (s > -1e-6 && s < 1e-6) { - end = End::BEGIN; + q.end = End::BEGIN; } else if (s > l-1e-6 && s < l+1e-6) { - end = End::END; + q.end = End::END; + } + +#if (defined(_debA) && defined(_debB)) + if (_idA == _debA && _idB == _debB) { + std::cout << "s " << s << ", " << q << std::endl; } - interPts.emplace_back(ptA[0]+t*r[0], ptA[1]+t*r[1], ptA[2]+t*r[2], t, idA, idB, end, src); +#endif + + interPts.push_back(std::move(q)); } @@ -444,12 +466,18 @@ void vtkPolyDataContactFilter::InterEdgeLine (InterPtsType &interPts, vtkPolyDat dB = std::sin(angB)*b; if ((std::isnan(dA) || dA < 1e-5) && (std::isnan(dB) || dB < 1e-5)) { -#ifdef DEBUG - std::cout << "congruent lines with d (" << dA << ", " << dB << ") and t (" << dotA << ", " << dotB << ") and l " << l << std::endl; + InterPt qA(ptA[0]+dotA*r[0], ptA[1]+dotA*r[1], ptA[2]+dotA*r[2], dotA, idA, idB, End::BEGIN, src, PointSrc::COPIED); + InterPt qB(ptA[0]+dotB*r[0], ptA[1]+dotB*r[1], ptA[2]+dotB*r[2], dotB, idA, idB, End::END, src, PointSrc::COPIED); + +#if (defined(_debA) && defined(_debB)) + if (_idA == _debA && _idB == _debB) { + std::cout << "dA " << dA << ", " << qA << std::endl; + std::cout << "dB " << dB << ", " << qB << std::endl; + } #endif - interPts.emplace_back(ptA[0]+dotA*r[0], ptA[1]+dotA*r[1], ptA[2]+dotA*r[2], dotA, idA, idB, End::BEGIN, src); - interPts.emplace_back(ptA[0]+dotB*r[0], ptA[1]+dotB*r[1], ptA[2]+dotB*r[2], dotB, idA, idB, End::END, src); + interPts.push_back(std::move(qA)); + interPts.push_back(std::move(qB)); } @@ -463,8 +491,10 @@ void vtkPolyDataContactFilter::InterEdgeLine (InterPtsType &interPts, vtkPolyDat bool vtkPolyDataContactFilter::InterPolyLine (InterPtsType &interPts, vtkPolyData *pd, vtkIdType num, const vtkIdType *poly, const double *r, const double *pt, Src src, const double *n) { -#ifdef DEBUG - std::cout << "InterPolyLine()" << std::endl; +#if (defined(_debA) && defined(_debB)) + if (_idA == _debA && _idB == _debB) { + std::cout << "InterPolyLine()" << std::endl; + } #endif std::vector edges; @@ -493,20 +523,37 @@ bool vtkPolyDataContactFilter::InterPolyLine (InterPtsType &interPts, vtkPolyDat } struct Cmp { + const double tol = 1e-5; + bool operator() (const InterPt &lhs, const InterPt &rhs) const { - if (lhs.end != End::NONE && rhs.end != End::NONE) { - const vtkIdType indA = lhs.end == End::BEGIN ? lhs.edge.f : lhs.edge.g, - indB = rhs.end == End::BEGIN ? rhs.edge.f : rhs.edge.g; + if (lhs.pointSrc == PointSrc::COPIED && rhs.pointSrc == PointSrc::COPIED) { + const vtkIdType indA = lhs.GetEnd(), + indB = rhs.GetEnd(); if (indA == indB) { - return indA < indB; + return false; + } + } else if (lhs.pointSrc == PointSrc::COPIED) { + const vtkIdType ind = lhs.GetEnd(); + + if (ind == rhs.edge.f || ind == rhs.edge.g) { + return false; + } + } else if (rhs.pointSrc == PointSrc::COPIED) { + const vtkIdType ind = rhs.GetEnd(); + + if (ind == lhs.edge.f || ind == lhs.edge.g) { + return false; } } - const long a = std::lround(lhs.t*1e5), - b = std::lround(rhs.t*1e5); + const double d = lhs.t-rhs.t; - return a < b; + if (std::abs(d) < tol) { + return false; + } else { + return d < 0; + } } }; @@ -518,57 +565,54 @@ bool vtkPolyDataContactFilter::InterPolyLine (InterPtsType &interPts, vtkPolyDat std::vector sortedPts; - for (auto &p : paired) { - InterPtsType &pts = p.second; + std::transform(paired.begin(), paired.end(), std::back_inserter(sortedPts), [](auto &kv) { + std::sort(kv.second.begin(), kv.second.end(), [](const InterPt &lhs, const InterPt &rhs) { return lhs.pointSrc > rhs.pointSrc; }); - if (pts.size() == 1 && pts.front().end != End::NONE) { - // hier fehlt der zweite punkt - pts.push_back(pts.back()); - } + return kv.second; + }); - sortedPts.push_back(pts); - } - - // trivial + std::map ends; - if (sortedPts.front().size() == 2) { - sortedPts.front().pop_back(); - } + decltype(sortedPts)::iterator itr; - if (sortedPts.back().size() == 2) { - sortedPts.back().pop_back(); - } + for (itr = sortedPts.begin(); itr != sortedPts.end(); ++itr) { + if (itr == sortedPts.begin()) { + if (itr->size() == 2) { + itr->pop_back(); + } + } else if (itr == sortedPts.end()-1) { + if (itr->size() == 2) { + itr->pop_back(); + } + } else if (itr->size() == 1 && itr->back().end != End::NONE) { + itr->push_back(itr->back()); + } - // behandelt schnitte durch kanten von sich selbst schneidenden einfachen polygonen + // validierung - for (const auto &pts : sortedPts) { - if (pts.size() > 2) { + if (itr->size() > 2) { return false; } - if (pts.size() == 2) { - if (pts.back().end == End::NONE) { + if (itr->size() == 2) { + if (itr->back().end == End::NONE) { return false; } - std::set ids {pts.front().edge.f, pts.front().edge.g, pts.back().edge.f, pts.back().edge.g}; + const Pair &edgeA = itr->front().edge, + &edgeB = itr->back().edge; - if (ids.size() == 2) { + if (edgeA.f == edgeB.g && edgeB.f == edgeA.g) { return false; } } - } - // ... + const auto &p = itr->back(); - std::map ends; - - for (const auto &pts : sortedPts) { - auto &last = pts.back(); - - if (last.end != End::NONE) { - ends.emplace(last.end == End::BEGIN ? last.edge.f : last.edge.g, last.t); + if (p.end != End::NONE) { + ends.emplace(p.end == End::BEGIN ? p.edge.f : p.edge.g, p.t); } + } double s[3], d; @@ -580,28 +624,28 @@ bool vtkPolyDataContactFilter::InterPolyLine (InterPtsType &interPts, vtkPolyDat vtkIdType id, prev, next; - for (auto &pts : sortedPts) { - auto &last = pts.back(); + for (itr = sortedPts.begin(); itr != sortedPts.end(); ++itr) { + const auto &p = itr->back(); - if (last.end != End::NONE) { - if (last.end == End::BEGIN) { - id = last.edge.f; - next = last.edge.g; + if (p.end != End::NONE) { + if (p.end == End::BEGIN) { + id = p.edge.f; + next = p.edge.g; prev = std::find_if(edges.begin(), edges.end(), [&id](const Pair &edge) { return edge.g == id; })->f; } else { - id = last.edge.g; - prev = last.edge.f; + id = p.edge.g; + prev = p.edge.f; next = std::find_if(edges.begin(), edges.end(), [&id](const Pair &edge) { return edge.f == id; })->g; } - if (pts.size() == 2) { + if (itr->size() == 2) { if (ends.count(next) == 0 && ends.count(prev) == 1) { pd->GetPoint(next, ptA); dA = vtkMath::Dot(s, ptA)-d; - if ((last.t > ends.at(prev) && dA > 0) || (last.t < ends.at(prev) && dA < 0)) { + if ((p.t > ends.at(prev) && dA > 0) || (p.t < ends.at(prev) && dA < 0)) { // tasche - pts.pop_back(); + itr->pop_back(); } continue; @@ -610,9 +654,9 @@ bool vtkPolyDataContactFilter::InterPolyLine (InterPtsType &interPts, vtkPolyDat pd->GetPoint(prev, ptB); dB = vtkMath::Dot(s, ptB)-d; - if ((last.t > ends.at(next) && dB < 0) || (last.t < ends.at(next) && dB > 0)) { + if ((p.t > ends.at(next) && dB < 0) || (p.t < ends.at(next) && dB > 0)) { // tasche - pts.pop_back(); + itr->pop_back(); } continue; @@ -628,8 +672,8 @@ bool vtkPolyDataContactFilter::InterPolyLine (InterPtsType &interPts, vtkPolyDat dB = vtkMath::Dot(s, ptB)-d; if (std::signbit(dA) != std::signbit(dB)) { - if (pts.size() == 2) { - pts.pop_back(); + if (itr->size() == 2) { + itr->pop_back(); } } else { @@ -640,7 +684,7 @@ bool vtkPolyDataContactFilter::InterPolyLine (InterPtsType &interPts, vtkPolyDat tB = vtkMath::Dot(vB, r); if ((tB > tA) == std::signbit(dA)) { - pts.clear(); + itr->clear(); } } @@ -648,8 +692,6 @@ bool vtkPolyDataContactFilter::InterPolyLine (InterPtsType &interPts, vtkPolyDat } } - // ... - InterPtsType _interPts; for (const auto &pts : sortedPts) { @@ -670,8 +712,12 @@ bool vtkPolyDataContactFilter::InterPolyLine (InterPtsType &interPts, vtkPolyDat void vtkPolyDataContactFilter::InterPolys (vtkIdType idA, vtkIdType idB) { -#ifdef DEBUG - std::cout << "InterPolys() -> idA " << idA << ", idB " << idB << std::endl; +#if (defined(_debA) && defined(_debB)) + _idA = idA; _idB = idB; + + if (_idA == _debA && _idB == _debB) { + std::cout << "InterPolys() -> idA " << idA << ", idB " << idB << std::endl; + } #endif vtkIdType numA, numB; @@ -680,6 +726,35 @@ void vtkPolyDataContactFilter::InterPolys (vtkIdType idA, vtkIdType idB) { newPdA->GetCellPoints(idA, numA, polyA); newPdB->GetCellPoints(idB, numB, polyB); +#if (defined(_debA) && defined(_debB)) + if (_idA == _debA && _idB == _debB) { + double p[3]; + + vtkIdType i; + + std::cout << _idA << " ["; + + for (i = 0; i < numA; i++) { + newPdA->GetPoint(polyA[i], p); + + std::cout << polyA[i] << ": [" << p[0] << ", " << p[1] << ", " << p[2] << "], "; + } + + std::cout << "]" << std::endl; + + std::cout << _idB << " ["; + + for (i = 0; i < numB; i++) { + newPdB->GetPoint(polyB[i], p); + + std::cout << polyB[i] << ": [" << p[0] << ", " << p[1] << ", " << p[2] << "], "; + } + + std::cout << "]" << std::endl; + } + +#endif + // ebenen aufstellen double nA[3], nB[3], ptA[3], ptB[3], dA, dB; @@ -687,6 +762,10 @@ void vtkPolyDataContactFilter::InterPolys (vtkIdType idA, vtkIdType idB) { ComputeNormal(newPdA->GetPoints(), nA, numA, polyA); ComputeNormal(newPdB->GetPoints(), nB, numB, polyB); + if (vtkMath::Dot(nA, nB) > .9999999999) { + return; + } + newPdA->GetPoint(polyA[0], ptA); newPdB->GetPoint(polyB[0], ptB); @@ -709,10 +788,6 @@ void vtkPolyDataContactFilter::InterPolys (vtkIdType idA, vtkIdType idB) { vtkMath::Cross(nA, nB, r); vtkMath::Normalize(r); - // std::cout << r[0] << ", " - // << r[1] << ", " - // << r[2] << std::endl; - // lsg. des lin. gls. mittels cramerscher regel int i = 0; @@ -744,10 +819,12 @@ void vtkPolyDataContactFilter::InterPolys (vtkIdType idA, vtkIdType idB) { s[inds[1]] = (nA[inds[0]]*dB-nB[inds[0]]*dA)/det; s[i] = 0; -#ifdef DEBUG - std::cout << "r [" << r[0] << ", " << r[1] << ", " << r[2] << "]" - << ", s [" << s[0] << ", " << s[1] << ", " << s[2] << "]" - << std::endl; +#if (defined(_debA) && defined(_debB)) + if (_idA == _debA && _idB == _debB) { + std::cout << "r [" << r[0] << ", " << r[1] << ", " << r[2] << "]" + << ", s [" << s[0] << ", " << s[1] << ", " << s[2] << "]" + << std::endl; + } #endif InterPtsType intersA, intersB; @@ -762,18 +839,23 @@ void vtkPolyDataContactFilter::InterPolys (vtkIdType idA, vtkIdType idB) { return; } -#ifdef DEBUG - std::cout << "intersA" << std::endl; +#if (defined(_debA) && defined(_debB)) + if (_idA == _debA && _idB == _debB) { - for (auto inter : intersA) { - std::cout << inter << std::endl; - } + std::cout << "intersA" << std::endl; + + for (auto inter : intersA) { + std::cout << inter << std::endl; + } + + std::cout << "intersB" << std::endl; - std::cout << "intersB" << std::endl; + for (auto inter : intersB) { + std::cout << inter << std::endl; + } - for (auto inter : intersB) { - std::cout << "B " << inter << std::endl; } + #endif if (intersA.size() != 0 && intersB.size() != 0 @@ -817,6 +899,12 @@ void vtkPolyDataContactFilter::OverlapLines (OverlapsType &ols, InterPtsType &in void vtkPolyDataContactFilter::AddContactLines (InterPtsType &intersA, InterPtsType &intersB, vtkIdType idA, vtkIdType idB) { +#if (defined(_debA) && defined(_debB)) + if (_idA == _debA && _idB == _debB) { + std::cout << "AddContactLines()" << std::endl; + } +#endif + OverlapsType overlaps; OverlapLines(overlaps, intersA, intersB, idA, idB); @@ -826,11 +914,6 @@ void vtkPolyDataContactFilter::AddContactLines (InterPtsType &intersA, InterPtsT auto &f = std::get<0>(*itr), &s = std::get<1>(*itr); -#ifdef DEBUG - std::cout << "f " << f << std::endl; - std::cout << "s " << s << std::endl; -#endif - if (f.src == Src::A) { if (edgesA.count(f.edge) == 1) { invalidA = true; @@ -855,6 +938,13 @@ void vtkPolyDataContactFilter::AddContactLines (InterPtsType &intersA, InterPtsT } } +#if (defined(_debA) && defined(_debB)) + if (_idA == _debA && _idB == _debB) { + std::cout << "line[0] " << f << std::endl; + std::cout << "line[1] " << s << std::endl; + } +#endif + vtkIdList *linePts = vtkIdList::New(); linePts->InsertNextId(contLines->GetPoints()->InsertNextPoint(f.pt)); @@ -939,3 +1029,476 @@ bool vtkPolyDataContactFilter::CheckInters (InterPtsType &interPts, vtkPolyData return true; } + +// #define _DEBUG + +PreventEqualCaptPoints::PreventEqualCaptPoints (vtkPolyData *pdA, vtkPolyData *pdB) : pdA(pdA), pdB(pdB) {} + +void PreventEqualCaptPoints::Run () { +#ifdef _DEBUG + WriteVTK("captA.vtk", pdA); + WriteVTK("captB.vtk", pdB); +#endif + + pdA->BuildLinks(); + pdB->BuildLinks(); + + Find(pdA, pdB, "A"); + +#ifdef _DEBUG + WriteVTK("modB.vtk", pdB); +#endif + + Find(pdB, pdA, "B"); + +#ifdef _DEBUG + WriteVTK("modA.vtk", pdA); +#endif +} + +void PreventEqualCaptPoints::Find (vtkPolyData *pd, vtkPolyData *other, [[maybe_unused]] const std::string &name) { +#ifdef _DEBUG + std::cout << "Find(" << name << ")" << std::endl; +#endif + + vtkIdType num; + const vtkIdType *poly; + + vtkIdType i, j; + + std::set lines; + + auto polyItr = vtk::TakeSmartPointer(pd->GetPolys()->NewIterator()); + + for (polyItr->GoToFirstCell(); !polyItr->IsDoneWithTraversal(); polyItr->GoToNextCell()) { + polyItr->GetCurrentCell(num, poly); + + for (i = 0; i < num; i++) { + j = i+1; + + if (j == num) { + j = 0; + } + + if (poly[i] < poly[j]) { + lines.emplace(poly[i], poly[j]); + } else { + lines.emplace(poly[j], poly[i]); + } + } + } + + auto tree = vtkSmartPointer::New(); + tree->SetDataSet(other); + tree->BuildLocator(); + + auto pts = vtkSmartPointer::New(); + auto cells = vtkSmartPointer::New(); + + double pA[3], pB[3]; + + vtkIdType cellId; + + double tr[2]; + +#ifdef _DEBUG + auto pdVerts = vtkSmartPointer::New(); + pdVerts->Allocate(1); + + auto ptsVerts = vtkSmartPointer::New(); +#endif + + std::map> pointSnaps; + std::map> edgeSnaps; + + for (auto &line : lines) { + pd->GetPoint(line.f, pA); + pd->GetPoint(line.g, pB); + + if (tree->IntersectWithLine(pA, pB, 1e-5, pts, cells) == 0) { + continue; + } + + for (i = 0; i < pts->GetNumberOfPoints(); i++) { + const double *pt = pts->GetPoint(i); + + Point3d sA(pt[0], pt[1], pt[2]); + + cellId = cells->GetId(i); + + other->GetCellPoints(cellId, num, poly); + + Base base(other->GetPoints(), num, poly); + + Poly polyA, polyB; + + GetPoly(other->GetPoints(), num, poly, polyA); + + FlattenPoly(polyA, polyB, base); + + Transform(pt, tr, base); + + Point3d sB(tr[0], tr[1], 0); + + if (PointInPoly(polyB, sB)) { +#ifdef _DEBUG + auto vert = vtkSmartPointer::New(); + vert->InsertNextId(ptsVerts->InsertNextPoint(pt)); + + pdVerts->InsertNextCell(VTK_VERTEX, vert); +#endif + + // snap auf ecke oder kante? + + auto snap = std::find_if(polyA.begin(), polyA.end(), [&](const Point3d &p) { return Point3d::GetDist(p, sA) < 1e-10; }); + + if (snap != polyA.end()) { + double d = Point3d::GetDist(*snap, sA); + + pointSnaps[snap->id].emplace_back(cellId, line, *snap, sA, d); + + } else { + // projektion auf kante + + auto edgeProj = GetEdgeProj(polyA, sA); + + if (edgeProj != nullptr) { + edgeSnaps[edgeProj->proj].emplace_back(cellId, line, edgeProj->edge, edgeProj->proj, sA, edgeProj->d); + } + } + } + } + } + + for (const auto& [id, snaps] : pointSnaps) { + if (snaps.size() > 1) { +#ifdef _DEBUG + std::cout << "id " << id << ", snaps" << std::endl; + + for (const auto &s : snaps) { + std::cout << s << std::endl; + } +#endif + + std::set allLines; + + std::transform(snaps.begin(), snaps.end(), std::inserter(allLines, allLines.end()), [](const SnapPoint &p) { return p.line; }); + +#ifdef _DEBUG + std::cout << "allLines" << std::endl; + + for (const auto &line : allLines) { + std::cout << line << std::endl; + } +#endif + + auto itrA = snaps.begin(); + decltype(itrA) itrB; + + double d; + + bool collapse = true; + + for (; itrA != snaps.end()-1 && collapse; ++itrA) { + for (itrB = itrA+1; itrB < snaps.end(); ++itrB) { + d = Point3d::GetDist(itrA->inter, itrB->inter); + +#ifdef _DEBUG + std::cout << d << std::endl; +#endif + + if (d > 1e-10) { + collapse = false; + break; + } + } + } + + if (collapse) { + continue; + } + + if (allLines.size() == 1) { + std::shared_ptr proj; + ProjOnLine(pd, snaps.back().line, snaps.back().point, proj); + + PreventEqualCaptPoints::MovePoint(other, id, *proj); + + } else { + // gemeinsamen punkt ermitteln + + std::set allIds; + + for (auto &line : allLines) { + allIds.insert(line.f); + allIds.insert(line.g); + } + + if (allIds.size() == allLines.size()+1) { + std::vector _allLines(allLines.begin(), allLines.end()); + + vtkIdType s = _allLines[0] & _allLines[1]; + +#ifdef _DEBUG + std::cout << "line " + << _allLines[0] + << " & " + << _allLines[1] + << " -> " + << s + << std::endl; +#endif + + double pt[3]; + pd->GetPoint(s, pt); + + Point3d p(pt[0], pt[1], pt[2]); + + PreventEqualCaptPoints::MovePoint(other, id, p); + } else { + throw std::runtime_error(""); + } + } + } + } + + using Snap = std::tuple; + + struct Cmp { + bool operator() (const Snap &l, const Snap &r) const { + return std::get<2>(l) < std::get<2>(r); + } + }; + + std::map> allEdgeSnaps; + + double pt[3], t; + + for (const auto& [proj, snaps] : edgeSnaps) { + if (snaps.size() > 1) { + if (snaps.size() > 2) { +#ifdef _DEBUG + std::cout << proj << std::endl; + + for (const auto &s : snaps) { + std::cout << s << std::endl; + + auto &line = s.line; + + pd->GetPoint(line.f, pA); + pd->GetPoint(line.g, pB); + + std::cout << Point3d(pA[0], pA[1], pA[2]) << std::endl; + std::cout << Point3d(pB[0], pB[1], pB[2]) << std::endl; + } +#endif + + continue; + } + + const auto &snapA = snaps[0]; + const auto &snapB = snaps[1]; + + if (Point3d::GetDist(snapA.inter, snapB.inter) > 1e-10) { + Pair edge(snapA.edge); + + if (edge.f > edge.g) { + std::swap(edge.f, edge.g); + } + + std::shared_ptr p; + + if (snapA.line == snapB.line) { + ProjOnLine(pd, snapA.line, snapA.proj, p); + + } else { + vtkIdType s = snapA.line & snapB.line; + + pd->GetPoint(s, pt); + + p = std::make_shared(pt[0], pt[1], pt[2]); + } + + other->GetPoint(edge.f, pt); + + Point3d q(pt[0], pt[1], pt[2]); + + t = Point3d::GetDist(q, snapA.proj); + + allEdgeSnaps[edge].emplace(snapA, *p, t); + } + } + } + + std::map newCells; + + for (const auto& [edge, data] : allEdgeSnaps) { + Points pts; + + for (const auto &d : data) { + const auto &p = std::get<1>(d); + pts.emplace_back(p.x, p.y, p.z, other->GetPoints()->InsertNextPoint(p.x, p.y, p.z)); + } + + const auto &first = std::get<0>(*(data.begin())); + + auto neigs = vtkSmartPointer::New(); + + other->GetCellEdgeNeighbors(first.cellId, first.edge.f, first.edge.g, neigs); + + if (neigs->GetNumberOfIds() != 1) { + throw std::runtime_error(""); + } + + vtkIdType neig = neigs->GetId(0); + + auto _edge = Pair(first.edge.g, first.edge.f); + + if (edge == first.edge) { + std::for_each(pts.begin(), pts.end(), [&](const Point3d &p) { newCells[first.cellId][first.edge].emplace_back(p); }); + std::for_each(pts.rbegin(), pts.rend(), [&](const Point3d &p) { newCells[neig][_edge].emplace_back(p); }); + } else { + std::for_each(pts.rbegin(), pts.rend(), [&](const Point3d &p) { newCells[first.cellId][first.edge].emplace_back(p); }); + std::for_each(pts.begin(), pts.end(), [&](const Point3d &p) { newCells[neig][_edge].emplace_back(p); }); + } + } + + for (const auto& [cellId, edges] : newCells) { + PreventEqualCaptPoints::TriangluteCell(other, cellId, edges); + } + + other->RemoveDeletedCells(); + +#ifdef _DEBUG + pdVerts->SetPoints(ptsVerts); + + auto fileName = "verts" + name + ".vtk"; + + WriteVTK(fileName.c_str(), pdVerts); +#endif + +} + +void PreventEqualCaptPoints::TriangluteCell (vtkPolyData *pd, vtkIdType cellId, const Edges &edges) { + vtkIdTypeArray *origCellIds = vtkIdTypeArray::SafeDownCast(pd->GetCellData()->GetScalars("OrigCellIds")); + + vtkIdType num; + const vtkIdType *poly; + + pd->GetCellPoints(cellId, num, poly); + + Poly _poly; + + double pt[3]; + + vtkIdType i, j; + + vtkIdType newNum = num; + + for (i = 0; i < num; i++) { + j = i+1; + + if (j == num) { + j = 0; + } + + pd->GetPoint(poly[i], pt); + + _poly.emplace_back(pt[0], pt[1], pt[2], poly[i]); + + Pair edge(poly[i], poly[j]); + + auto itr = edges.find(edge); + + if (itr != edges.end()) { + auto &pts = itr->second; + + for (const auto &p : pts) { + _poly.emplace_back(p); + + newNum++; + } + } + } + + auto vtkPoly = vtkSmartPointer::New(); + + vtkPoly->GetPointIds()->SetNumberOfIds(newNum); + vtkPoly->GetPoints()->SetNumberOfPoints(newNum); + + Base base(pd->GetPoints(), num, poly); + + Poly flattened; + + FlattenPoly(_poly, flattened, base); + + for (const auto &p : flattened) { + vtkPoly->GetPointIds()->SetId(p.id, p.id); + vtkPoly->GetPoints()->SetPoint(p.id, p.x, p.y, p.z); + } + + auto triangles = vtkSmartPointer::New(); + +#if (VTK_MAJOR_VERSION >= 9 && VTK_MINOR_VERSION > 3) + if (vtkPoly->TriangulateLocalIds(0, triangles) != 1) { + throw std::runtime_error(""); + } +#else + if (vtkPoly->Triangulate(triangles) != 1) { + throw std::runtime_error(""); + } +#endif + + auto ids = vtkSmartPointer::New(); + + for (const auto &p : _poly) { + ids->InsertNextId(p.id); + } + + vtkIdType origId = origCellIds->GetValue(cellId); + + auto triangle = vtkSmartPointer::New(); + triangle->SetNumberOfIds(3); + + for (i = 0; i < triangles->GetNumberOfIds(); i += 3) { + triangle->SetId(0, ids->GetId(triangles->GetId(i))); + triangle->SetId(1, ids->GetId(triangles->GetId(i+1))); + triangle->SetId(2, ids->GetId(triangles->GetId(i+2))); + + pd->InsertNextCell(VTK_TRIANGLE, triangle); + + origCellIds->InsertNextValue(origId); + } + + pd->DeleteCell(cellId); + +} + +void PreventEqualCaptPoints::MovePoint (vtkPolyData *pd, vtkIdType ind, const Point3d &p) { + auto cells = vtkSmartPointer::New(); + pd->GetPointCells(ind, cells); + + vtkIdType i, cellId; + + for (i = 0; i < cells->GetNumberOfIds(); i++) { + cellId = cells->GetId(i); + + if (pd->GetCellType(cellId) == VTK_POLYGON) { + PreventEqualCaptPoints::TriangluteCell(pd, cellId, {}); + } + } + +#ifdef _DEBUG + double pt[3]; + pd->GetPoints()->GetPoint(ind, pt); + + Point3d q(pt[0], pt[1], pt[2]); + + std::cout << q + << " -> " + << p + << std::endl; +#endif + + pd->GetPoints()->SetPoint(ind, p.x, p.y, p.z); +} diff --git a/CombineModels/Logic/vtkPolyDataContactFilter.h b/CombineModels/Logic/vtkPolyDataContactFilter.h index 84f5650..9e48216 100644 --- a/CombineModels/Logic/vtkPolyDataContactFilter.h +++ b/CombineModels/Logic/vtkPolyDataContactFilter.h @@ -1,5 +1,5 @@ /* -Copyright 2012-2024 Ronald Römer +Copyright 2012-2025 Ronald Römer Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. @@ -22,6 +22,7 @@ limitations under the License. #include #include #include +#include #include #include @@ -42,11 +43,16 @@ enum class End { END }; +enum class PointSrc { + CALCULATED, + COPIED +}; + class InterPt { public: InterPt () = delete; - InterPt (double x, double y, double z, double t, vtkIdType a, vtkIdType b, End end, Src src) : t(t), edge(a, b), end(end), src(src), srcA(NOTSET), srcB(NOTSET) { + InterPt (double x, double y, double z, double t, vtkIdType a, vtkIdType b, End end, Src src, PointSrc pointSrc) : t(t), edge(a, b), end(end), src(src), srcA(NOTSET), srcB(NOTSET), pointSrc(pointSrc) { pt[0] = x; pt[1] = y; pt[2] = z; @@ -58,12 +64,15 @@ class InterPt { Src src; vtkIdType srcA, srcB; + PointSrc pointSrc; + friend std::ostream& operator<< (std::ostream &out, const InterPt &s) { out << "pt [" << s.pt[0] << ", " << s.pt[1] << ", " << s.pt[2] << "]" << ", t " << s.t << ", edge " << s.edge << ", end " << s.end - << ", src " << s.src; + << ", src " << s.src + << ", pointSrc " << s.pointSrc; return out; } @@ -72,20 +81,32 @@ class InterPt { assert(src != other.src); if (src == Src::A) { - srcA = end == End::END ? edge.g : edge.f; + srcA = GetEnd(); } else { - srcB = end == End::END ? edge.g : edge.f; + srcB = GetEnd(); } if (std::abs(other.t-t) < 1e-5) { if (other.src == Src::A) { - srcA = other.end == End::END ? other.edge.g : other.edge.f; + srcA = other.GetEnd(); } else { - srcB = other.end == End::END ? other.edge.g : other.edge.f; + srcB = other.GetEnd(); } } } + inline vtkIdType GetEnd() const { + if (end == End::BEGIN) { + return edge.f; + } + + if (end == End::END) { + return edge.g; + } + + return NOTSET; + } + }; typedef std::vector InterPtsType; @@ -139,4 +160,65 @@ class VTK_SLICER_COMBINEMODELS_MODULE_LOGIC_EXPORT vtkPolyDataContactFilter : pu }; + +using Edges = std::map; + +class SnapPoint { +public: + SnapPoint () = delete; + SnapPoint (vtkIdType cellId, const Pair &line, const Point3d &point, const Point3d &inter, double d) : cellId(cellId), line(line), point(point), inter(inter), d(d) {} + + vtkIdType cellId; + Pair line; + Point3d point; + Point3d inter; + double d; + + friend std::ostream& operator<< (std::ostream &out, const SnapPoint &s) { + out << "cellId " << s.cellId + << ", line " << s.line + << ", point " << s.point + << ", inter " << s.inter + << ", d " << s.d; + + return out; + } +}; + +class SnapEdge { +public: + SnapEdge () = delete; + SnapEdge (vtkIdType cellId, const Pair &line, const Pair &edge, const Point3d &proj, const Point3d &inter, double d) : cellId(cellId), line(line), edge(edge), proj(proj), inter(inter), d(d) {} + + vtkIdType cellId; + Pair line; + Pair edge; + Point3d proj; + Point3d inter; + double d; + + friend std::ostream& operator<< (std::ostream &out, const SnapEdge &s) { + out << "cellId " << s.cellId + << ", line " << s.line + << ", edge " << s.edge + << ", proj " << s.proj + << ", inter " << s.inter + << ", d " << s.d; + + return out; + } +}; + +class PreventEqualCaptPoints { + vtkPolyData *pdA, *pdB; +public: + static void TriangluteCell (vtkPolyData *pd, vtkIdType cellId, const Edges &edges); + static void MovePoint (vtkPolyData *pd, vtkIdType ind, const Point3d &p); + + PreventEqualCaptPoints (vtkPolyData *pdA, vtkPolyData *pdB); + void Run (); +private: + void Find (vtkPolyData *pd, vtkPolyData *other, const std::string &name); +}; + #endif