--- /dev/null
+/*
+Copyright 2012-2022 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.
+You may obtain a copy of the License at
+
+ http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+#include <map>
+#include <deque>
+#include <vector>
+#include <set>
+#include <algorithm>
+#include <cmath>
+#include <functional>
+#include <queue>
+#include <memory>
+
+#include <chrono>
+#include <numeric>
+#include <iterator>
+
+#include <vtkInformation.h>
+#include <vtkInformationVector.h>
+#include <vtkDemandDrivenPipeline.h>
+#include <vtkObjectFactory.h>
+#include <vtkPolyDataAlgorithm.h>
+#include <vtkCellData.h>
+#include <vtkPointData.h>
+#include <vtkMath.h>
+#include <vtkIdList.h>
+#include <vtkAppendPolyData.h>
+#include <vtkKdTreePointLocator.h>
+#include <vtkCleanPolyData.h>
+#include <vtkPolyDataConnectivityFilter.h>
+#include <vtkSmartPointer.h>
+#include <vtkModifiedBSPTree.h>
+#include <vtkCellArrayIterator.h>
+#include <vtkKdTree.h>
+
+#include "vtkPolyDataBooleanFilter.h"
+#include "vtkPolyDataContactFilter.h"
+
+#include "Utilities.h"
+
+vtkStandardNewMacro(vtkPolyDataBooleanFilter);
+
+vtkPolyDataBooleanFilter::vtkPolyDataBooleanFilter () {
+
+ SetNumberOfInputPorts(2);
+ SetNumberOfOutputPorts(3);
+
+ timePdA = 0;
+ timePdB = 0;
+
+ contLines = vtkSmartPointer<vtkPolyData>::New();
+
+ modPdA = vtkSmartPointer<vtkPolyData>::New();
+ modPdB = vtkSmartPointer<vtkPolyData>::New();
+
+ cellDataA = vtkSmartPointer<vtkCellData>::New();
+ cellDataB = vtkSmartPointer<vtkCellData>::New();
+
+ cellIdsA = vtkSmartPointer<vtkIdTypeArray>::New();
+ cellIdsB = vtkSmartPointer<vtkIdTypeArray>::New();
+
+ OperMode = OPER_UNION;
+
+}
+
+vtkPolyDataBooleanFilter::~vtkPolyDataBooleanFilter () {
+ // nix mehr
+}
+
+int vtkPolyDataBooleanFilter::ProcessRequest(vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector) {
+
+ if (request->Has(vtkDemandDrivenPipeline::REQUEST_DATA())) {
+
+ vtkInformation *inInfoA = inputVector[0]->GetInformationObject(0);
+ vtkInformation *inInfoB = inputVector[1]->GetInformationObject(0);
+
+ vtkPolyData *pdA = vtkPolyData::SafeDownCast(inInfoA->Get(vtkDataObject::DATA_OBJECT()));
+ vtkPolyData *pdB = vtkPolyData::SafeDownCast(inInfoB->Get(vtkDataObject::DATA_OBJECT()));
+
+ vtkInformation *outInfoA = outputVector->GetInformationObject(0);
+ vtkInformation *outInfoB = outputVector->GetInformationObject(1);
+ vtkInformation *outInfoC = outputVector->GetInformationObject(2);
+
+ resultA = vtkPolyData::SafeDownCast(outInfoA->Get(vtkDataObject::DATA_OBJECT()));
+ resultB = vtkPolyData::SafeDownCast(outInfoB->Get(vtkDataObject::DATA_OBJECT()));
+ resultC = vtkPolyData::SafeDownCast(outInfoC->Get(vtkDataObject::DATA_OBJECT()));
+
+ using clock = std::chrono::steady_clock;
+ std::vector<clock::duration> times;
+ clock::time_point start;
+
+ if (pdA->GetMTime() > timePdA || pdB->GetMTime() > timePdB) {
+ // eventuell vorhandene regionen vereinen
+
+ vtkSmartPointer<vtkCleanPolyData> cleanA = vtkSmartPointer<vtkCleanPolyData>::New();
+ cleanA->SetOutputPointsPrecision(DOUBLE_PRECISION);
+ cleanA->SetTolerance(1e-6);
+ cleanA->SetInputData(pdA);
+ cleanA->Update();
+
+ vtkSmartPointer<vtkCleanPolyData> cleanB = vtkSmartPointer<vtkCleanPolyData>::New();
+ cleanB->SetOutputPointsPrecision(DOUBLE_PRECISION);
+ cleanB->SetTolerance(1e-6);
+ cleanB->SetInputData(pdB);
+ cleanB->Update();
+
+#ifdef DEBUG
+ std::cout << "Exporting modPdA.vtk" << std::endl;
+ WriteVTK("modPdA.vtk", cleanA->GetOutput());
+
+ std::cout << "Exporting modPdB.vtk" << std::endl;
+ WriteVTK("modPdB.vtk", cleanB->GetOutput());
+#endif
+
+ // CellData sichern
+
+ cellDataA->DeepCopy(cleanA->GetOutput()->GetCellData());
+ cellDataB->DeepCopy(cleanB->GetOutput()->GetCellData());
+
+ // ermittelt kontaktstellen
+
+ start = clock::now();
+
+ vtkSmartPointer<vtkPolyDataContactFilter> cl = vtkSmartPointer<vtkPolyDataContactFilter>::New();
+ cl->SetInputConnection(0, cleanA->GetOutputPort());
+ cl->SetInputConnection(1, cleanB->GetOutputPort());
+ cl->Update();
+
+ times.push_back(clock::now()-start);
+
+ contLines->DeepCopy(cl->GetOutput());
+
+ modPdA->DeepCopy(cl->GetOutput(1));
+ modPdB->DeepCopy(cl->GetOutput(2));
+
+#ifdef DEBUG
+ std::cout << "Exporting contLines.vtk" << std::endl;
+ WriteVTK("contLines.vtk", contLines);
+
+ std::cout << "Exporting modPdA_1.vtk" << std::endl;
+ WriteVTK("modPdA_1.vtk", modPdA);
+
+ std::cout << "Exporting modPdB_1.vtk" << std::endl;
+ WriteVTK("modPdB_1.vtk", modPdB);
+#endif
+
+ if (contLines->GetNumberOfCells() == 0) {
+ vtkErrorMacro("Inputs have no contact.");
+
+ return 1;
+ }
+
+ // in den CellDatas steht drin, welche polygone einander schneiden
+
+ contsA = vtkIdTypeArray::SafeDownCast(contLines->GetCellData()->GetScalars("cA"));
+ contsB = vtkIdTypeArray::SafeDownCast(contLines->GetCellData()->GetScalars("cB"));
+
+ vtkIdTypeArray *accuracy = vtkIdTypeArray::SafeDownCast(contLines->GetPointData()->GetScalars("accuracy"));
+
+ vtkIdTypeArray *sourcesA = vtkIdTypeArray::SafeDownCast(contLines->GetCellData()->GetScalars("sourcesA"));
+ vtkIdTypeArray *sourcesB = vtkIdTypeArray::SafeDownCast(contLines->GetCellData()->GetScalars("sourcesB"));
+
+ vtkIdType i, numPts = contLines->GetNumberOfPoints();
+
+ vtkSmartPointer<vtkIdList> cells = vtkSmartPointer<vtkIdList>::New();
+
+ for (i = 0; i < numPts; i++) {
+ contLines->GetPointCells(i, cells);
+
+ if (cells->GetNumberOfIds() == 1) {
+ vtkErrorMacro("Contact ends suddenly.");
+ return 1;
+ }
+
+ if (accuracy->GetValue(i) == 1) {
+ vtkErrorMacro("Contact goes through a cell with a bad shape.");
+ return 1;
+ }
+ }
+
+ // sichert die OrigCellIds
+
+ vtkIdTypeArray *origCellIdsA = vtkIdTypeArray::SafeDownCast(modPdA->GetCellData()->GetScalars("OrigCellIds"));
+ vtkIdTypeArray *origCellIdsB = vtkIdTypeArray::SafeDownCast(modPdB->GetCellData()->GetScalars("OrigCellIds"));
+
+ cellIdsA->DeepCopy(origCellIdsA);
+ cellIdsB->DeepCopy(origCellIdsB);
+
+ vtkIdType numCellsA = modPdA->GetNumberOfCells(),
+ numCellsB = modPdB->GetNumberOfCells();
+
+ for (i = 0; i < numCellsA; i++) {
+ origCellIdsA->SetValue(i, i);
+ }
+
+ for (i = 0; i < numCellsB; i++) {
+ origCellIdsB->SetValue(i, i);
+ }
+
+ start = clock::now();
+
+ if (GetPolyStrips(modPdA, contsA, sourcesA, polyStripsA) ||
+ GetPolyStrips(modPdB, contsB, sourcesB, polyStripsB)) {
+
+ vtkErrorMacro("Strips are invalid.");
+
+ return 1;
+
+ }
+
+ times.push_back(clock::now()-start);
+
+ // trennt die polygone an den linien
+
+ start = clock::now();
+
+ if (CutCells(modPdA, polyStripsA) ||
+ CutCells(modPdB, polyStripsB)) {
+
+ vtkErrorMacro("CutCells failed.");
+
+ return 1;
+ }
+
+ times.push_back(clock::now()-start);
+
+#ifdef DEBUG
+ std::cout << "Exporting modPdA_2.vtk" << std::endl;
+ WriteVTK("modPdA_2.vtk", modPdA);
+
+ std::cout << "Exporting modPdB_2.vtk" << std::endl;
+ WriteVTK("modPdB_2.vtk", modPdB);
+#endif
+
+ start = clock::now();
+
+ RestoreOrigPoints(modPdA, polyStripsA);
+ RestoreOrigPoints(modPdB, polyStripsB);
+
+ times.push_back(clock::now()-start);
+
+#ifdef DEBUG
+ std::cout << "Exporting modPdA_3.vtk" << std::endl;
+ WriteVTK("modPdA_3.vtk", modPdA);
+
+ std::cout << "Exporting modPdB_3.vtk" << std::endl;
+ WriteVTK("modPdB_3.vtk", modPdB);
+#endif
+
+ start = clock::now();
+
+ ResolveOverlaps(modPdA, contsA, polyStripsA);
+ ResolveOverlaps(modPdB, contsB, polyStripsB);
+
+ times.push_back(clock::now()-start);
+
+#ifdef DEBUG
+ std::cout << "Exporting modPdA_4.vtk" << std::endl;
+ WriteVTK("modPdA_4.vtk", modPdA);
+
+ std::cout << "Exporting modPdB_4.vtk" << std::endl;
+ WriteVTK("modPdB_4.vtk", modPdB);
+#endif
+
+ start = clock::now();
+
+ AddAdjacentPoints(modPdA, contsA, polyStripsA);
+ AddAdjacentPoints(modPdB, contsB, polyStripsB);
+
+ times.push_back(clock::now()-start);
+
+#ifdef DEBUG
+ std::cout << "Exporting modPdA_5.vtk" << std::endl;
+ WriteVTK("modPdA_5.vtk", modPdA);
+
+ std::cout << "Exporting modPdB_5.vtk" << std::endl;
+ WriteVTK("modPdB_5.vtk", modPdB);
+#endif
+
+ start = clock::now();
+
+ DisjoinPolys(modPdA, polyStripsA);
+ DisjoinPolys(modPdB, polyStripsB);
+
+ times.push_back(clock::now()-start);
+
+#ifdef DEBUG
+ std::cout << "Exporting modPdA_6.vtk" << std::endl;
+ WriteVTK("modPdA_6.vtk", modPdA);
+
+ std::cout << "Exporting modPdB_6.vtk" << std::endl;
+ WriteVTK("modPdB_6.vtk", modPdB);
+#endif
+
+ start = clock::now();
+
+ MergePoints(modPdA, polyStripsA);
+ MergePoints(modPdB, polyStripsB);
+
+ times.push_back(clock::now()-start);
+
+#ifdef DEBUG
+ std::cout << "Exporting modPdA_7.vtk" << std::endl;
+ WriteVTK("modPdA_7.vtk", modPdA);
+
+ std::cout << "Exporting modPdB_7.vtk" << std::endl;
+ WriteVTK("modPdB_7.vtk", modPdB);
+#endif
+
+ timePdA = pdA->GetMTime();
+ timePdB = pdB->GetMTime();
+
+ }
+
+ start = clock::now();
+
+ if (CombineRegions()) {
+ vtkErrorMacro("Boolean operation failed.");
+
+ return 1;
+ }
+
+ times.push_back(clock::now()-start);
+
+// #ifdef DEBUG
+ double sum = std::chrono::duration_cast<std::chrono::duration<double>>(std::accumulate(times.begin(), times.end(), clock::duration())).count();
+
+ std::vector<clock::duration>::const_iterator itr;
+ for (itr = times.begin(); itr != times.end(); itr++) {
+ double time = std::chrono::duration_cast<std::chrono::duration<double>>(*itr).count();
+
+ std::cout << "Time " << (itr-times.begin())
+ << ": " << time << "s (" << (time/sum*100) << "%)"
+ << std::endl;
+ }
+// #endif
+
+ }
+
+ return 1;
+
+}
+
+void vtkPolyDataBooleanFilter::GetStripPoints (vtkPolyData *pd, vtkIdTypeArray *sources, PStrips &pStrips, IdsType &lines) {
+
+#ifdef DEBUG
+ std::cout << "GetStripPoints()" << std::endl;
+#endif
+
+ StripPtsType &pts = pStrips.pts;
+ const IdsType &poly = pStrips.poly;
+
+ double a[3], b[3], sA[3], sB[3], u[3], v[3], w[3], n, t, d;
+
+ std::map<vtkIdType, vtkIdType> allPts;
+
+ vtkSmartPointer<vtkIdList> line = vtkSmartPointer<vtkIdList>::New();
+
+ for (auto lineId : lines) {
+ contLines->GetCellPoints(lineId, line);
+
+ // std::cout << "? " << contsA->GetValue(lineId)
+ // << ", " << contsB->GetValue(lineId)
+ // << ", [" << line->GetId(0)
+ // << ", " << line->GetId(1)
+ // << "]"
+ // << std::endl;
+
+ allPts.emplace(line->GetId(0), sources->GetTypedComponent(lineId, 0));
+ allPts.emplace(line->GetId(1), sources->GetTypedComponent(lineId, 1));
+ }
+
+ decltype(allPts)::const_iterator itr;
+
+ for (itr = allPts.begin(); itr != allPts.end(); ++itr) {
+ StripPt sp;
+ sp.ind = itr->first;
+
+ // die koordinaten
+ double pt[3];
+ contLines->GetPoint(sp.ind, pt);
+
+ Cpy(sp.pt, pt, 3);
+
+ IdsType::const_iterator itrA, itrB;
+
+ for (itrA = poly.begin(); itrA != poly.end(); ++itrA) {
+ itrB = itrA+1;
+
+ if (itrB == poly.end()) {
+ itrB = poly.begin();
+ }
+
+ if (itr->second != NOTSET && *itrA != itr->second) {
+ continue;
+ }
+
+ pd->GetPoint(*itrA, a);
+ pd->GetPoint(*itrB, b);
+
+ vtkMath::Subtract(a, pt, sA);
+ vtkMath::Subtract(b, pt, sB);
+
+ // richtungsvektor und länge der kante
+
+ vtkMath::Subtract(b, a, u);
+ n = vtkMath::Norm(u);
+
+ // d und t zur kante
+
+ vtkMath::Subtract(pt, a, v);
+ t = vtkMath::Dot(v, u)/(n*n);
+
+ vtkMath::Cross(v, u, w);
+ d = vtkMath::Norm(w)/n;
+
+ if (d < 1e-5 && t > -1e-5 && t < 1+1e-5) {
+ sp.edge[0] = *itrA;
+ sp.edge[1] = *itrB;
+
+ sp.t = std::min(1., std::max(0., t));
+
+ if (vtkMath::Norm(sA) < 1e-5) {
+ Cpy(sp.captPt, a, 3);
+ sp.capt = Capt::A;
+
+ } else if (vtkMath::Norm(sB) < 1e-5) {
+ Cpy(sp.captPt, b, 3);
+ sp.capt = Capt::B;
+
+ } else {
+ // u ist nicht normiert
+ vtkMath::MultiplyScalar(u, t);
+
+ double x[3];
+ vtkMath::Add(a, u, x);
+
+ // projektion
+ Cpy(sp.captPt, x, 3);
+
+ sp.capt = Capt::EDGE;
+
+ }
+ }
+
+ // std::cout << "? "
+ // << sp.ind
+ // << ", " << d
+ // << ", " << t
+ // << ", " << sp.capt
+ // << std::endl;
+ }
+
+ if (itr->second != NOTSET && sp.edge[0] == NOTSET) {
+ sp.catched = false;
+ }
+
+ pts.emplace(sp.ind, std::move(sp));
+
+ }
+
+ StripPtsType::iterator itr2;
+
+ for (itr2 = pts.begin(); itr2 != pts.end(); ++itr2) {
+ StripPt &sp = itr2->second;
+
+ if (sp.capt != Capt::NOT) {
+ if (sp.capt == Capt::B) {
+ sp.t = 0;
+
+ sp.edge[0] = sp.edge[1];
+
+ auto itrA = std::find(poly.begin(), poly.end(), sp.edge[0]),
+ itrB = itrA+1;
+
+ if (itrB == poly.end()) {
+ itrB = poly.begin();
+ }
+
+ sp.edge[1] = *itrB;
+
+ sp.capt = Capt::A;
+
+ }
+
+ // für den schnitt werden die eingerasteten koordinaten verwendet
+
+ Cpy(sp.cutPt, sp.captPt, 3);
+ } else {
+
+ Cpy(sp.cutPt, sp.pt, 3);
+ }
+
+ sp.history.emplace_back(sp.edge[0], sp.edge[1]);
+
+ }
+
+#ifdef DEBUG
+ std::cout << "pts: " << std::endl;
+ for (itr2 = pts.begin(); itr2 != pts.end(); ++itr2) {
+ std::cout << itr2->first << ": " << itr2->second << std::endl;
+ }
+#endif
+
+}
+
+bool vtkPolyDataBooleanFilter::GetPolyStrips (vtkPolyData *pd, vtkIdTypeArray *conts, vtkIdTypeArray *sources, PolyStripsType &polyStrips) {
+#ifdef DEBUG
+ std::cout << "GetPolyStrips()" << std::endl;
+#endif
+
+ polyStrips.clear();
+
+ std::map<vtkIdType, IdsType> polyLines;
+
+ for (vtkIdType i = 0; i < conts->GetNumberOfTuples(); i++) {
+ vtkIdType poly = conts->GetValue(i);
+
+ /*if (poly != 1641) {
+ continue;
+ }*/
+
+ polyLines[poly].push_back(i);
+ }
+
+ std::vector<std::reference_wrapper<StripPt>> notCatched;
+
+ std::map<vtkIdType, IdsType>::iterator itr;
+
+ for (itr = polyLines.begin(); itr != polyLines.end(); ++itr) {
+
+ IdsType &lines = itr->second;
+ RemoveDuplicates(lines);
+
+ polyStrips.emplace(std::piecewise_construct,
+ std::forward_as_tuple(itr->first),
+ std::forward_as_tuple(pd, itr->first));
+
+ PStrips &pStrips = polyStrips.at(itr->first);
+
+ GetStripPoints(pd, sources, pStrips, lines);
+
+ for (auto &sp : pStrips.pts) {
+ sp.second.polyId = itr->first;
+
+ if (!sp.second.catched) {
+ notCatched.push_back(sp.second);
+ }
+ }
+
+ }
+
+ auto Next = [](const IdsType &ids, vtkIdType id) -> vtkIdType {
+ IdsType::const_iterator itr;
+
+ itr = std::find(ids.begin(), ids.end(), id);
+
+ if (++itr == ids.end()) {
+ itr = ids.begin();
+ }
+
+ return *itr;
+ };
+
+ for (StripPt &sp : notCatched) {
+ for (itr = polyLines.begin(); itr != polyLines.end(); ++itr) {
+ const PStrips &pStrips = polyStrips.at(itr->first);
+
+ try {
+ const StripPt &corr = pStrips.pts.at(sp.ind);
+
+ if (&corr != &sp) {
+ if (corr.capt == Capt::A) {
+ sp.capt = Capt::A;
+ sp.edge[0] = corr.edge[0];
+ sp.edge[1] = Next(polyStrips.at(sp.polyId).poly, sp.edge[0]);
+
+ sp.t = 0;
+
+ Cpy(sp.captPt, corr.captPt, 3);
+ Cpy(sp.cutPt, sp.captPt, 3);
+
+ sp.history.emplace_back(sp.edge[0], sp.edge[1]);
+
+ sp.catched = true;
+
+ }
+ }
+ } catch (...) {}
+
+ }
+
+ if (!sp.catched) {
+ std::cout << sp << std::endl;
+ }
+
+ assert(sp.catched);
+
+ }
+
+ vtkSmartPointer<vtkIdList> cells = vtkSmartPointer<vtkIdList>::New(),
+ line = vtkSmartPointer<vtkIdList>::New();
+
+ vtkIdType i, numCells;
+
+ StripPtsType::const_iterator itr2;
+
+ for (itr = polyLines.begin(); itr != polyLines.end(); ++itr) {
+ PStrips &pStrips = polyStrips.at(itr->first);
+
+ const IdsType &lines = itr->second;
+ const StripPtsType &pts = pStrips.pts;
+
+ for (itr2 = pts.begin(); itr2 != pts.end(); ++itr2) {
+ // ist der punkt auf einer kante, dürfen von ihm mehr als 2 linien ausgehen
+
+ const StripPt &pt = itr2->second;
+
+ if (pt.capt == Capt::NOT) {
+ contLines->GetPointCells(pt.ind, cells);
+
+ numCells = cells->GetNumberOfIds();
+
+ std::set<vtkIdType> ends;
+
+ for (i = 0; i < numCells; i++) {
+ contLines->GetCellPoints(cells->GetId(i), line);
+
+ ends.insert(pt.ind == line->GetId(0) ? line->GetId(1) : line->GetId(0));
+ }
+
+ if (ends.size() > 2) {
+ return true;
+ }
+ }
+ }
+
+ // zusammensetzen
+
+ std::deque<Pair> _lines;
+
+ vtkIdList *linePts = vtkIdList::New();
+
+ for (auto &i : lines) {
+ contLines->GetCellPoints(i, linePts);
+ _lines.emplace_back(linePts->GetId(0), linePts->GetId(1));
+ }
+
+ linePts->Delete();
+
+ decltype(_lines)::iterator _itr;
+
+ auto FindRight = [&pts, &_lines, &_itr](StripType &strip, const std::size_t &id) -> bool {
+ auto &right = strip.back();
+
+ if (pts.at(right.ind).capt == Capt::NOT) {
+ for (_itr = _lines.begin(); _itr != _lines.end(); ++_itr) {
+ if (_itr->f == right.ind) {
+ strip.emplace_back(_itr->g, id);
+
+ _lines.erase(_itr);
+ return true;
+ } else if (_itr->g == right.ind) {
+ strip.emplace_back(_itr->f, id);
+
+ _lines.erase(_itr);
+ return true;
+ }
+ }
+ }
+
+ return false;
+ };
+
+ auto FindLeft = [&pts, &_lines, &_itr](StripType &strip, const std::size_t &id) -> bool {
+ auto &left = strip.front();
+
+ if (pts.at(left.ind).capt == Capt::NOT) {
+ for (_itr = _lines.begin(); _itr != _lines.end(); ++_itr) {
+ if (_itr->f == left.ind) {
+ strip.emplace_front(_itr->g, id);
+
+ _lines.erase(_itr);
+ return true;
+ } else if (_itr->g == left.ind) {
+ strip.emplace_front(_itr->f, id);
+
+ _lines.erase(_itr);
+ return true;
+ }
+ }
+ }
+
+ return false;
+ };
+
+ std::size_t stripId {0};
+
+ while (!_lines.empty()) {
+ auto &last = _lines.back();
+
+ StripType strip {{last.f, stripId}, {last.g, stripId}};
+ _lines.pop_back();
+
+ while (FindRight(strip, stripId)) {}
+ while (FindLeft(strip, stripId)) {}
+
+ pStrips.strips.push_back(std::move(strip));
+
+ stripId++;
+
+ }
+
+ CompleteStrips(pStrips);
+
+ }
+
+ // sucht nach schnitten zw. den strips
+
+ {
+
+ PolyStripsType::const_iterator itr;
+ StripType::const_iterator itr2;
+
+ for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) {
+ const PStrips &pStrips = itr->second;
+
+ const StripsType &strips = pStrips.strips;
+ const StripPtsType &pts = pStrips.pts;
+ const Base &base = pStrips.base;
+
+ vtkSmartPointer<vtkPoints> treePts = vtkSmartPointer<vtkPoints>::New();
+
+ vtkSmartPointer<vtkPolyData> treePd = vtkSmartPointer<vtkPolyData>::New();
+ treePd->Allocate(1);
+
+ std::map<vtkIdType, vtkIdType> ptIds;
+
+ double pt[2];
+
+ for (const auto &p : pts) {
+ Transform(p.second.pt, pt, base);
+
+ ptIds.emplace(p.first, treePts->InsertNextPoint(pt[0], pt[1], 0));
+ }
+
+ for (const StripType &strip : strips) {
+ for (itr2 = strip.begin(); itr2 != strip.end()-1; ++itr2) {
+
+ vtkIdList *line = vtkIdList::New();
+ line->InsertNextId(ptIds[itr2->ind]);
+ line->InsertNextId(ptIds[(itr2+1)->ind]);
+
+ treePd->InsertNextCell(VTK_LINE, line);
+
+ line->Delete();
+ }
+ }
+
+ treePd->SetPoints(treePts);
+
+ vtkSmartPointer<vtkModifiedBSPTree> tree = vtkSmartPointer<vtkModifiedBSPTree>::New();
+ tree->SetDataSet(treePd);
+ tree->BuildLocator();
+
+ vtkIdType numA, numB;
+ const vtkIdType *lineA, *lineB;
+
+ auto lineItr = vtk::TakeSmartPointer(treePd->GetLines()->NewIterator());
+
+ for (lineItr->GoToFirstCell(); !lineItr->IsDoneWithTraversal(); lineItr->GoToNextCell()) {
+ lineItr->GetCurrentCell(numA, lineA);
+
+ double ptA[3], ptB[3];
+
+ treePts->GetPoint(lineA[0], ptA);
+ treePts->GetPoint(lineA[1], ptB);
+
+ vtkSmartPointer<vtkIdList> lineIds = vtkSmartPointer<vtkIdList>::New();
+
+ tree->IntersectWithLine(ptA, ptB, 1e-5, nullptr, lineIds);
+
+ for (vtkIdType i = 0; i < lineIds->GetNumberOfIds(); i++) {
+ treePd->GetCellPoints(lineIds->GetId(i), numB, lineB);
+
+ if (lineB[0] != lineA[0] && lineB[1] != lineA[0] && lineB[0] != lineA[1] && lineB[1] != lineA[1]) {
+ // schnitt gefunden
+
+ return true;
+ }
+ }
+ }
+
+ }
+
+ }
+
+ return false;
+
+}
+
+void vtkPolyDataBooleanFilter::RemoveDuplicates (IdsType &lines) {
+
+ IdsType unique;
+
+ // die indexe der enden auf übereinstimmung prüfen
+
+ vtkIdList *linePtsA = vtkIdList::New();
+ vtkIdList *linePtsB = vtkIdList::New();
+
+ std::size_t i, j, numLines = lines.size();
+
+ for (i = 0; i < numLines-1; i++) {
+ j = i+1;
+
+ contLines->GetCellPoints(lines[i], linePtsA);
+
+ while (j < lines.size()) {
+ contLines->GetCellPoints(lines[j], linePtsB);
+
+ if ((linePtsA->GetId(0) == linePtsB->GetId(0) && linePtsA->GetId(1) == linePtsB->GetId(1)) ||
+ (linePtsA->GetId(0) == linePtsB->GetId(1) && linePtsA->GetId(1) == linePtsB->GetId(0))) {
+ // stimmen überein
+ break;
+ }
+
+ j++;
+ }
+
+ if (j == numLines) {
+ // keine vorzeitige unterbrechung der while-schleife
+ unique.push_back(lines[i]);
+ }
+ }
+
+ unique.push_back(lines.back());
+
+ linePtsA->Delete();
+ linePtsB->Delete();
+
+ lines.swap(unique);
+
+}
+
+void vtkPolyDataBooleanFilter::CompleteStrips (PStrips &pStrips) {
+ StripsType::iterator itr;
+
+ for (itr = pStrips.strips.begin(); itr != pStrips.strips.end(); ++itr) {
+ const StripPt &start = pStrips.pts[itr->front().ind],
+ &end = pStrips.pts[itr->back().ind];
+
+ if (start.ind != end.ind) {
+ if (start.capt == Capt::NOT) {
+ StripType s(itr->rbegin(), itr->rend()-1);
+ itr->insert(itr->begin(), s.begin(), s.end());
+
+ } else if (end.capt == Capt::NOT) {
+ StripType s(itr->rbegin()+1, itr->rend());
+ itr->insert(itr->end(), s.begin(), s.end());
+
+ }
+ }
+ }
+}
+
+bool vtkPolyDataBooleanFilter::HasArea (const StripType &strip) const {
+ bool area = true;
+
+ std::size_t i, n = strip.size();
+
+ if (n%2 == 1) {
+ for (i = 0; i < (n-1)/2; i++) {
+ area = strip[i].ind != strip[n-i-1].ind;
+ }
+ }
+
+ return area;
+}
+
+void ComputeNormal (const StripPtsType &pts, const RefsType &poly, double *n) {
+ n[0] = 0; n[1] = 0; n[2] = 0;
+
+ RefsType::const_iterator itrA, itrB;
+
+ for (itrA = poly.begin(); itrA != poly.end(); ++itrA) {
+ itrB = itrA+1;
+
+ if (itrB == poly.end()) {
+ itrB = poly.begin();
+ }
+
+ const StripPtR &spA = *itrA,
+ &spB = *itrB;
+
+ auto pA = pts.find(spA.ind);
+ auto pB = pts.find(spB.ind);
+
+ const double *ptA = pA->second.cutPt,
+ *ptB = pB->second.cutPt;
+
+ n[0] += (ptA[1]-ptB[1])*(ptA[2]+ptB[2]);
+ n[1] += (ptA[2]-ptB[2])*(ptA[0]+ptB[0]);
+ n[2] += (ptA[0]-ptB[0])*(ptA[1]+ptB[1]);
+ }
+
+ vtkMath::Normalize(n);
+}
+
+bool vtkPolyDataBooleanFilter::CutCells (vtkPolyData *pd, PolyStripsType &polyStrips) {
+#ifdef DEBUG
+ std::cout << "CutCells()" << std::endl;
+#endif
+
+ vtkPoints *pdPts = pd->GetPoints();
+
+ vtkIdTypeArray *origCellIds = vtkIdTypeArray::SafeDownCast(pd->GetCellData()->GetScalars("OrigCellIds"));
+
+ PolyStripsType::iterator itr;
+
+ for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) {
+
+ vtkIdType polyInd = itr->first;
+ PStrips &pStrips = itr->second;
+
+ /*if (polyInd != 1641) {
+ continue;
+ }*/
+
+ StripsType &strips = pStrips.strips;
+ StripPtsType &pts = pStrips.pts;
+
+ IdsType &poly = pStrips.poly;
+
+ vtkIdType origId = origCellIds->GetValue(polyInd);
+
+#ifdef DEBUG
+ std::cout << "polyInd " << polyInd << ", poly [";
+ for (auto &p : poly) {
+ std::cout << p << ", ";
+ }
+ std::cout << "]" << std::endl;
+#endif
+
+ std::size_t numPts = poly.size();
+
+ std::map<vtkIdType, RefsType> edges;
+
+ StripsType::iterator itr2;
+ StripType::iterator itr3;
+
+ // holes sichern
+ StripsType holes;
+
+ auto fct = [&](const StripType &s) {
+ return pts[s.front().ind].capt == Capt::NOT && pts[s.back().ind].capt == Capt::NOT;
+ };
+
+ std::copy_if(strips.begin(), strips.end(), std::back_inserter(holes), fct);
+
+ strips.erase(std::remove_if(strips.begin(), strips.end(), fct), strips.end());
+
+ for (itr2 = strips.begin(); itr2 != strips.end(); ++itr2) {
+ StripType &strip = *itr2;
+
+#ifdef DEBUG
+ std::cout << (itr2-strips.begin()) << ". strip [";
+ for (auto &s : strip) {
+ std::cout << s.ind << ", ";
+ }
+ std::cout << "]" << std::endl;
+#endif
+
+ // init
+ if (pts[strip.front().ind].edge[0] == pts[strip.back().ind].edge[0]
+ && strip.front().ind != strip.back().ind
+ && pts[strip.front().ind].t > pts[strip.back().ind].t) {
+
+ std::reverse(strip.begin(), strip.end());
+ }
+
+ StripPt &start = pts[strip.front().ind],
+ &end = pts[strip.back().ind];
+
+ strip.front().side = Side::START;
+ strip.back().side = Side::END;
+
+ strip.front().ref = start.edge[0];
+ strip.back().ref = end.edge[0];
+
+ // nachfolgend könnte man dann anfang und ende weglassen
+
+ for (itr3 = strip.begin(); itr3 != strip.end(); ++itr3) {
+ StripPt &sp = pts[itr3->ind];
+
+ itr3->desc[0] = pdPts->InsertNextPoint(sp.cutPt);
+ itr3->desc[1] = pdPts->InsertNextPoint(sp.cutPt);
+
+#ifdef DEBUG
+ std::cout << sp << " => " << *itr3 << std::endl;
+#endif
+
+ }
+
+ // ordnet zu
+
+ edges[start.edge[0]].push_back(std::ref(strip.front()));
+ edges[end.edge[0]].push_back(std::ref(strip.back()));
+ }
+
+ // sortiert die punkte auf den kanten
+
+ std::map<vtkIdType, RefsType>::iterator itr4;
+
+ IdsType::iterator itr5;
+ StripType::reverse_iterator itr6;
+
+ RefsType::iterator itr7;
+
+ for (itr4 = edges.begin(); itr4 != edges.end(); ++itr4) {
+ RefsType &edge = itr4->second;
+
+#ifdef DEBUG
+ std::cout << "edge (" << itr4->first << ", _)" << std::endl;
+#endif
+
+ std::sort(edge.begin(), edge.end(), [&](const StripPtR &a, const StripPtR &b) {
+ const StripPt &a_ = pts[a.ind],
+ &b_ = pts[b.ind];
+
+#ifdef DEBUG
+ // std::cout << "a_: " << a_ << " -> strip " << a.strip << std::endl;
+ // std::cout << "b_: " << b_ << " -> strip " << b.strip << std::endl;
+#endif
+
+ if (a_.ind == b_.ind) {
+ // strips beginnen im gleichen punkt
+
+ if (a.strip != b.strip) {
+ // gehören nicht dem gleichen strip an
+
+ StripType &stripA = strips[a.strip],
+ &stripB = strips[b.strip];
+
+ // andere enden ermitteln
+
+ const StripPtR &eA = (&a == &(stripA.front())) ? stripA.back() : stripA.front(),
+ &eB = (&b == &(stripB.front())) ? stripB.back() : stripB.front();
+
+ const StripPt &eA_ = pts[eA.ind],
+ &eB_ = pts[eB.ind];
+
+#ifdef DEBUG
+ // std::cout << "eA_: " << eA_ << std::endl;
+ // std::cout << "eB_: " << eA_ << std::endl;
+#endif
+
+ if (eA_.ind != eB_.ind) {
+ auto i = std::find(poly.begin(), poly.end(), itr4->first)-poly.begin(),
+ iA = std::find(poly.begin(), poly.end(), eA_.edge[0])-poly.begin(),
+ iB = std::find(poly.begin(), poly.end(), eB_.edge[0])-poly.begin();
+
+ double dA = static_cast<double>(Mod(iA-i, numPts))+eA_.t,
+ dB = static_cast<double>(Mod(iB-i, numPts))+eB_.t;
+
+ if (i == iA && a_.t > eA_.t) {
+ dA += static_cast<double>(numPts);
+ }
+
+ if (i == iB && b_.t > eB_.t) {
+ dB += static_cast<double>(numPts);
+ }
+
+#ifdef DEBUG
+ // std::cout << "dA=" << dA << ", dB=" << dB << std::endl;
+#endif
+
+ return dB < dA;
+ } else {
+ RefsType poly_;
+
+ if (a.side == Side::START) {
+ poly_.insert(poly_.end(), stripA.begin(), stripA.end());
+ } else {
+ poly_.insert(poly_.end(), stripA.rbegin(), stripA.rend());
+ }
+
+ if (b.side == Side::START) {
+ poly_.insert(poly_.end(), stripB.rbegin()+1, stripB.rend()-1);
+ } else {
+ poly_.insert(poly_.end(), stripB.begin()+1, stripB.end()-1);
+ }
+
+ double n[3];
+ ComputeNormal(pts, poly_, n);
+
+ double ang = vtkMath::Dot(pStrips.n, n);
+
+#ifdef DEBUG
+ // std::cout << "ang=" << ang*180/PI << std::endl;
+#endif
+
+ return ang < .999999;
+
+ }
+ } else {
+ // gleicher strip
+
+ StripType &strip = strips[a.strip];
+
+ if (HasArea(strip)) {
+ RefsType poly_(strip.begin(), strip.end()-1);
+
+ double n[3];
+ ComputeNormal(pts, poly_, n);
+
+ double ang = vtkMath::Dot(pStrips.n, n);
+
+ if (ang > .999999) {
+ std::reverse(strip.begin(), strip.end());
+ return true;
+ } else {
+ return false;
+ }
+
+ } else {
+ // reihenfolge von a und b bereits richtig
+ return false;
+ }
+ }
+
+ } else {
+ return a_.t < b_.t;
+ }
+ });
+
+#ifdef DEBUG
+ for (auto& p : edge) {
+ std::cout << p << std::endl;
+ }
+#endif
+
+ }
+
+ // baut die strips ein
+
+ std::deque<IdsType> polys;
+ polys.push_back(pStrips.poly);
+
+ for (itr2 = strips.begin(); itr2 != strips.end(); ++itr2) {
+ StripType &strip = *itr2;
+
+ StripPtR &start = strip.front(),
+ &end = strip.back();
+
+#ifdef DEBUG
+ std::cout << "strip " << start.strip
+ << " , refs (" << start.ref << ", " << end.ref << ")"
+ << std::endl;
+#endif
+
+ std::size_t cycle = 0;
+
+ while (true) {
+ if (cycle == polys.size()) {
+ break;
+ }
+
+ IdsType next = polys.front();
+ polys.pop_front();
+
+ std::vector<IdsType> newPolys(2);
+
+ if (std::find(next.begin(), next.end(), start.ref) != next.end()) {
+ if (start.ref == end.ref) {
+ for (itr5 = next.begin(); itr5 != next.end(); ++itr5) {
+ newPolys[0].push_back(*itr5);
+
+#ifdef DEBUG
+ std::cout << "adding " << *itr5 << " to 0" << std::endl;
+#endif
+
+ if (*itr5 == start.ref) {
+ for (itr3 = strip.begin(); itr3 != strip.end(); ++itr3) {
+ newPolys[0].push_back(itr3->desc[0]);
+
+#ifdef DEBUG
+ std::cout << "adding " << itr3->desc[0] << " to 0" << std::endl;
+#endif
+
+ }
+ }
+ }
+
+ // strip selbst ist ein polygon
+
+ for (itr6 = strip.rbegin(); itr6 != strip.rend(); ++itr6) {
+ newPolys[1].push_back(itr6->desc[1]);
+
+#ifdef DEBUG
+ std::cout << "adding " << itr6->desc[1] << " to 1" << std::endl;
+#endif
+
+ }
+
+ } else {
+ std::size_t curr = 0;
+
+ for (itr5 = next.begin(); itr5 != next.end(); ++itr5) {
+ IdsType &poly = newPolys[curr];
+
+ poly.push_back(*itr5);
+
+#ifdef DEBUG
+ std::cout << "adding " << *itr5 << " to " << curr << std::endl;
+#endif
+
+ if (*itr5 == start.ref) {
+ for (itr3 = strip.begin(); itr3 != strip.end(); ++itr3) {
+ poly.push_back(itr3->desc[0]);
+
+#ifdef DEBUG
+ std::cout << "adding " << itr3->desc[0] << " to " << curr << std::endl;
+#endif
+
+ }
+
+ curr = curr == 0 ? 1 : 0;
+
+ } else if (*itr5 == end.ref) {
+ for (itr6 = strip.rbegin(); itr6 != strip.rend(); ++itr6) {
+ poly.push_back(itr6->desc[1]);
+
+#ifdef DEBUG
+ std::cout << "adding " << itr6->desc[1] << " to " << curr << std::endl;
+#endif
+
+ }
+
+ curr = curr == 0 ? 1 : 0;
+ }
+
+ }
+ }
+ }
+
+ if (!newPolys[1].empty()) {
+ // refs aktualisieren
+
+ for (itr4 = edges.begin(); itr4 != edges.end(); ++itr4) {
+ RefsType &edge = itr4->second;
+
+ RefsType::iterator itrA;
+
+ for (itrA = edge.begin()+1; itrA != edge.end(); ++itrA) {
+ StripPtR &sp = *itrA;
+
+ if (sp.strip > start.strip) {
+#ifdef DEBUG
+ std::cout << "sp: ind " << sp.ind << ", strip " << sp.strip << std::endl;
+#endif
+
+ RefsType::const_reverse_iterator itrB(itrA);
+
+ std::shared_ptr<StripPtR> _p;
+
+ for (; itrB != edge.rend(); ++itrB) {
+ const StripPtR &p = *itrB;
+
+ if (p.strip != sp.strip) {
+ if (p.strip <= start.strip) {
+#ifdef DEBUG
+ std::cout << "ref " << sp.ref;
+#endif
+
+ if (p.side == Side::END) {
+ sp.ref = p.desc[0];
+ } else {
+ sp.ref = p.desc[1];
+ }
+
+#ifdef DEBUG
+ std::cout << " -> " << sp.ref << " (from strip " << p.strip << ", ind " << p.ind << ")" << std::endl;
+#endif
+
+ _p = std::make_shared<StripPtR>(p);
+
+ break;
+
+ }
+ } else {
+#ifdef DEBUG
+ std::cout << "~1 ref " << sp.ref << " -> " << p.ref << " (from strip " << p.strip << ", ind " << p.ind << ")" << std::endl;
+#endif
+
+ sp.ref = p.ref;
+ break;
+ }
+ }
+
+ RefsType::const_iterator itrC(itrA);
+
+ ++itrC;
+
+ for (; itrC != edge.end(); ++itrC) {
+ const StripPtR &p = *itrC;
+
+ if (p.ind != sp.ind) {
+ break;
+ }
+
+ if (p.strip <= start.strip) {
+ if (_p && p.ind == _p->ind && p.strip < _p->strip) {
+ break;
+ }
+
+#ifdef DEBUG
+ std::cout << "~2 ref " << sp.ref;
+#endif
+
+ if (p.side == Side::START) {
+ sp.ref = p.desc[0];
+ } else {
+ sp.ref = p.desc[1];
+ }
+
+#ifdef DEBUG
+ std::cout << " -> " << sp.ref << " (from strip " << p.strip << ", ind " << p.ind << ")" << std::endl;
+#endif
+
+ break;
+ }
+
+ }
+ }
+ }
+
+ // erstellt die history
+
+ auto _s = std::find_if(edge.begin(), edge.end(), [&start](const StripPtR &r) {
+ return &start == &r;
+ });
+
+ if (_s != edge.end()) {
+ for (itr7 = edge.begin(); itr7 != edge.end(); ++itr7) {
+ StripPtR &sp = *itr7;
+ StripPt &_sp = pts[sp.ind];
+ auto &history = _sp.history;
+
+ if (&sp != &start) {
+ if (_sp.t < pts[start.ind].t) {
+ history.push_back({history.back().f, start.desc[0]});
+ } else {
+ history.push_back({start.desc[1], history.back().g});
+ }
+
+ }
+
+ }
+ }
+
+ auto _e = std::find_if(edge.begin(), edge.end(), [&end](const StripPtR &r) {
+ return &end == &r;
+ });
+
+ if (_e != edge.end()) {
+ for (itr7 = edge.begin(); itr7 != edge.end(); ++itr7) {
+ StripPtR &sp = *itr7;
+ StripPt &_sp = pts[sp.ind];
+ auto &history = _sp.history;
+
+ if (&sp != &end) {
+ if (_sp.t < pts[end.ind].t) {
+ history.push_back({history.back().f, end.desc[1]});
+ } else {
+ history.push_back({end.desc[0], history.back().g});
+ }
+
+ }
+
+ }
+ }
+
+ // sonderfall
+ if (edge.size() > 1) {
+ StripPtR &a = edge.front(),
+ &b = *(edge.begin()+1);
+
+ if (a.ind == b.ind
+ && b.strip == start.strip
+ && pts[a.ind].capt == Capt::A) { // sollte weg
+
+#ifdef DEBUG
+ std::cout << "~3 ref " << a.ref;
+#endif
+
+ if (b.side == Side::START) {
+ a.ref = b.desc[0];
+ } else {
+ a.ref = b.desc[1];
+ }
+
+#ifdef DEBUG
+ std::cout << " -> " << a.ref << " (from strip " << b.strip << ", ind " << b.ind << ")" << std::endl;
+#endif
+
+ }
+ }
+
+ }
+
+ // doppelte punkte entfernen
+
+ double pt[3];
+
+ for (auto &newPoly : newPolys) {
+ IdsType _newPoly;
+
+ std::map<vtkIdType, Point3d> _pts;
+
+ for (vtkIdType id : newPoly) {
+ pd->GetPoint(id, pt);
+
+ _pts.emplace(std::piecewise_construct,
+ std::forward_as_tuple(id),
+ std::forward_as_tuple(pt[0], pt[1], pt[2]));
+ }
+
+ IdsType::const_iterator itrA, itrB;
+
+ for (itrA = newPoly.begin(); itrA != newPoly.end(); ++itrA) {
+ itrB = itrA+1;
+ if (itrB == newPoly.end()) {
+ itrB = newPoly.begin();
+ }
+
+ auto _a = _pts.find(*itrA);
+ auto _b = _pts.find(*itrB);
+
+ if (_a->second == _b->second) {
+ // doppelt
+#ifdef DEBUG
+ std::cout << "removing " << *itrA << std::endl;
+#endif
+ } else {
+ _newPoly.push_back(*itrA);
+ }
+ }
+
+ newPoly.swap(_newPoly);
+
+ }
+
+ // prüft, ob die erstellten polys gültig sind
+
+ if (newPolys[0].size() > 2) {
+ polys.push_back(newPolys[0]);
+ }
+
+ if (HasArea(strip) && newPolys[1].size() > 2) {
+ polys.push_back(newPolys[1]);
+ }
+
+ cycle = 0;
+
+ break;
+
+ } else {
+ polys.push_back(next);
+
+ cycle++;
+ }
+
+ }
+
+ }
+
+ // erzeugte polys hinzufügen
+
+ IdsType descIds;
+ descIds.reserve(polys.size());
+
+ for (auto &p : polys) {
+ vtkIdList *cell = vtkIdList::New();
+
+ for (vtkIdType &id : p) {
+ cell->InsertNextId(id);
+ }
+
+ descIds.emplace_back(pd->InsertNextCell(VTK_POLYGON, cell));
+ origCellIds->InsertNextValue(origId);
+
+ cell->Delete();
+ }
+
+ pd->DeleteCell(polyInd);
+
+ // holes einbauen
+ if (!holes.empty()) {
+ try {
+ Merger(pd, pStrips, holes, descIds, origId).Run();
+ } catch (const std::runtime_error &e) {
+ return true;
+ }
+ }
+
+ }
+
+ pd->RemoveDeletedCells();
+ pd->BuildCells();
+
+ return false;
+
+}
+
+void vtkPolyDataBooleanFilter::RestoreOrigPoints (vtkPolyData *pd, PolyStripsType &polyStrips) {
+
+#ifdef DEBUG
+ std::cout << "RestoreOrigPoints()" << std::endl;
+#endif
+
+ pd->BuildLinks();
+
+ vtkKdTreePointLocator *loc = vtkKdTreePointLocator::New();
+ loc->SetDataSet(pd);
+ loc->BuildLocator();
+
+ struct Cmp {
+ bool operator() (const StripPt &l, const StripPt &r) const {
+ return l.ind < r.ind;
+ }
+ };
+
+ PolyStripsType::const_iterator itr;
+ StripPtsType::const_iterator itr2;
+
+ std::set<std::reference_wrapper<const StripPt>, Cmp> ends;
+
+ for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) {
+ const PStrips &pStrips = itr->second;
+
+ for (itr2 = pStrips.pts.begin(); itr2 != pStrips.pts.end(); ++itr2) {
+ const StripPt &sp = itr2->second;
+
+ if (sp.capt != Capt::NOT) {
+ ends.emplace(sp);
+ }
+
+ }
+ }
+
+ vtkIdList *pts = vtkIdList::New();
+
+ vtkIdType i, numPts;
+
+ for (const StripPt &sp : ends) {
+ FindPoints(loc, sp.cutPt, pts);
+ numPts = pts->GetNumberOfIds();
+
+ for (i = 0; i < numPts; i++) {
+ pd->GetPoints()->SetPoint(pts->GetId(i), sp.pt);
+ }
+ }
+
+ pts->Delete();
+
+ loc->FreeSearchStructure();
+ loc->Delete();
+
+}
+
+void vtkPolyDataBooleanFilter::DisjoinPolys (vtkPolyData *pd, PolyStripsType &polyStrips) {
+
+#ifdef DEBUG
+ std::cout << "DisjoinPolys()" << std::endl;
+#endif
+
+ pd->BuildLinks();
+
+ vtkKdTreePointLocator *loc = vtkKdTreePointLocator::New();
+ loc->SetDataSet(pd);
+ loc->BuildLocator();
+
+ struct Cmp {
+ bool operator() (const StripPt &l, const StripPt &r) const {
+ return l.ind < r.ind;
+ }
+ };
+
+ PolyStripsType::const_iterator itr;
+ StripPtsType::const_iterator itr2;
+
+ std::set<std::reference_wrapper<const StripPt>, Cmp> ends;
+
+ for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) {
+ const PStrips &pStrips = itr->second;
+
+ for (itr2 = pStrips.pts.begin(); itr2 != pStrips.pts.end(); ++itr2) {
+ const StripPt &sp = itr2->second;
+
+ if (sp.capt == Capt::A) {
+ ends.emplace(sp);
+ }
+ }
+ }
+
+ vtkIdList *pts = vtkIdList::New();
+ vtkIdList *cells = vtkIdList::New();
+
+ vtkIdType i, j, numPts, numCells;
+
+ for (const StripPt &sp : ends) {
+ FindPoints(loc, sp.pt, pts);
+ numPts = pts->GetNumberOfIds();
+
+ for (i = 0; i < numPts; i++) {
+ pd->GetPointCells(pts->GetId(i), cells);
+ numCells = cells->GetNumberOfIds();
+
+ if (numCells > 1) {
+ for (j = 0; j < numCells; j++) {
+ pd->ReplaceCellPoint(cells->GetId(j), pts->GetId(i), pd->GetPoints()->InsertNextPoint(sp.pt));
+ }
+ }
+ }
+ }
+
+ cells->Delete();
+ pts->Delete();
+
+ loc->FreeSearchStructure();
+ loc->Delete();
+
+}
+
+void vtkPolyDataBooleanFilter::ResolveOverlaps (vtkPolyData *pd, vtkIdTypeArray *conts, PolyStripsType &polyStrips) {
+
+#ifdef DEBUG
+ std::cout << "ResolveOverlaps()" << std::endl;
+#endif
+
+ pd->BuildCells();
+ pd->BuildLinks();
+
+ contLines->BuildLinks();
+
+ vtkKdTreePointLocator *loc = vtkKdTreePointLocator::New();
+ loc->SetDataSet(pd);
+ loc->BuildLocator();
+
+ struct Cmp {
+ bool operator() (const StripPt &l, const StripPt &r) const {
+ return l.ind < r.ind;
+ }
+ };
+
+ PolyStripsType::const_iterator itr;
+ StripPtsType::const_iterator itr2;
+
+ std::set<std::reference_wrapper<const StripPt>, Cmp> ends;
+
+ for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) {
+ const PStrips &pStrips = itr->second;
+
+ for (itr2 = pStrips.pts.begin(); itr2 != pStrips.pts.end(); ++itr2) {
+ const StripPt &sp = itr2->second;
+
+ if (sp.capt == Capt::EDGE) {
+ ends.emplace(sp);
+ }
+ }
+ }
+
+ vtkIdList *links = vtkIdList::New(),
+ *ptsA = vtkIdList::New(),
+ *ptsB = vtkIdList::New(),
+ *cells = vtkIdList::New(),
+ *poly = vtkIdList::New();
+
+ double ptA[3], ptB[3];
+
+ vtkIdType i, j, k, numPtsA, numPtsB, numCells, numPts, idA, idB;
+
+ std::map<Pair, IdsType> sharedPts;
+
+ for (const StripPt &sp : ends) {
+ contLines->GetPointCells(sp.ind, links);
+
+ if (links->GetNumberOfIds() == 2
+ && conts->GetValue(links->GetId(0)) != conts->GetValue(links->GetId(1))) {
+
+ std::vector<Pair> history(sp.history.rbegin(), sp.history.rend());
+
+ for (const Pair &p : history) {
+ pd->GetPoint(p.f, ptA);
+ pd->GetPoint(p.g, ptB);
+
+ FindPoints(loc, ptA, ptsA);
+ FindPoints(loc, ptB, ptsB);
+
+ numPtsA = ptsA->GetNumberOfIds();
+ numPtsB = ptsB->GetNumberOfIds();
+
+ std::vector<Pair> cellsA, cellsB;
+
+ for (i = 0; i < numPtsA; i++) {
+ pd->GetPointCells(ptsA->GetId(i), cells);
+ numCells = cells->GetNumberOfIds();
+
+ for (j = 0; j < numCells; j++) {
+ cellsA.emplace_back(ptsA->GetId(i), cells->GetId(j));
+ }
+ }
+
+ for (i = 0; i < numPtsB; i++) {
+ pd->GetPointCells(ptsB->GetId(i), cells);
+ numCells = cells->GetNumberOfIds();
+
+ for (j = 0; j < numCells; j++) {
+ cellsB.emplace_back(ptsB->GetId(i), cells->GetId(j));
+ }
+ }
+
+ for (auto &cellA : cellsA) {
+ for (auto &cellB : cellsB) {
+ if (cellA.g == cellB.g) {
+ // kante/poly existiert noch
+
+ pd->GetCellPoints(cellA.g, poly);
+ numPts = poly->GetNumberOfIds();
+
+ for (k = 0; k < numPts; k++) {
+ idA = poly->GetId(k);
+ idB = k+1 == numPts ? poly->GetId(0) : poly->GetId(k+1);
+
+ if (cellB.f == idA
+ && cellA.f == idB) {
+
+ IdsType &pts = sharedPts[Pair{sp.ind, cellA.g}];
+
+ pts.push_back(cellA.f);
+ pts.push_back(cellB.f);
+
+ break;
+ }
+ }
+
+ }
+ }
+ }
+
+
+ }
+
+ }
+ }
+
+ double pt[3];
+
+ decltype(sharedPts)::iterator itr3;
+
+ for (itr3 = sharedPts.begin(); itr3 != sharedPts.end(); ++itr3) {
+ const Pair &p = itr3->first;
+
+ IdsType &pts = itr3->second;
+
+ std::sort(pts.begin(), pts.end());
+
+ auto found = std::adjacent_find(pts.begin(), pts.end());
+
+ if (found != pts.end()) {
+ contLines->GetPoint(p.f, pt);
+
+ i = pd->GetPoints()->InsertNextPoint(pt);
+ pd->ReplaceCellPoint(p.g, *found, i);
+ }
+ }
+
+ links->Delete();
+ ptsA->Delete();
+ ptsB->Delete();
+ cells->Delete();
+ poly->Delete();
+
+ loc->FreeSearchStructure();
+ loc->Delete();
+
+}
+
+void vtkPolyDataBooleanFilter::AddAdjacentPoints (vtkPolyData *pd, vtkIdTypeArray *conts, PolyStripsType &polyStrips) {
+
+#ifdef DEBUG
+ std::cout << "AddAdjacentPoints()" << std::endl;
+#endif
+
+ pd->BuildLinks();
+
+ vtkIdTypeArray *origCellIds = vtkIdTypeArray::SafeDownCast(pd->GetCellData()->GetScalars("OrigCellIds"));
+
+ struct Cmp {
+ bool operator() (const StripPt &l, const StripPt &r) const {
+ return l.t < r.t;
+ }
+ };
+
+ vtkKdTreePointLocator *loc = vtkKdTreePointLocator::New();
+ loc->SetDataSet(pd);
+ loc->BuildLocator();
+
+ vtkIdList *cells = vtkIdList::New(),
+ *ptsA = vtkIdList::New(),
+ *ptsB = vtkIdList::New(),
+ *poly = vtkIdList::New(),
+ *newPoly = vtkIdList::New();
+
+ vtkIdType i, j, numCells, numPtsA, numPtsB, numPts, idA, idB;
+
+ PolyStripsType::const_iterator itr;
+ StripPtsType::const_iterator itr2;
+
+ for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) {
+ const PStrips &pStrips = itr->second;
+
+ std::map<Pair, std::set<std::reference_wrapper<const StripPt>, Cmp>> edgePts;
+
+ for (itr2 = pStrips.pts.begin(); itr2 != pStrips.pts.end(); ++itr2) {
+ const StripPt &sp = itr2->second;
+
+ if (sp.capt == Capt::EDGE) {
+ edgePts[{sp.edge[0], sp.edge[1]}].emplace(sp);
+ }
+ }
+
+ decltype(edgePts)::iterator itr3;
+
+ for (itr3 = edgePts.begin(); itr3 != edgePts.end(); ++itr3) {
+ const Pair &edge = itr3->first;
+ auto pts = itr3->second;
+
+ StripPt spA, spB;
+
+ pd->GetPoint(edge.f, spA.pt);
+ pd->GetPoint(edge.g, spB.pt);
+
+ spA.t = 0;
+ spB.t = 1;
+
+ pts.emplace(spA);
+ pts.emplace(spB);
+
+ std::vector<decltype(pts)::value_type> _pts(pts.rbegin(), pts.rend());
+
+ decltype(_pts)::const_iterator itrA, itrB, itrC;
+
+ itrA = _pts.begin();
+
+ while (itrA != _pts.end()-1) {
+ itrB = itrA+1;
+
+ while (itrB != _pts.end()-1) {
+ contLines->GetPointCells(itrB->get().ind, cells);
+ numCells = cells->GetNumberOfIds();
+
+ // std::cout << "[";
+ // for (i = 0; i < numCells; i++) {
+ // std::cout << conts->GetValue(cells->GetId(i)) << ", ";
+ // }
+ // std::cout << "]" << std::endl;
+
+ // break;
+
+ ++itrB;
+ }
+
+ if (std::distance(itrA, itrB) > 1) {
+ FindPoints(loc, itrA->get().pt, ptsA);
+ FindPoints(loc, itrB->get().pt, ptsB);
+
+ numPtsA = ptsA->GetNumberOfIds();
+ numPtsB = ptsB->GetNumberOfIds();
+
+ std::vector<Pair> polysA, polysB;
+
+ for (i = 0; i < numPtsA; i++) {
+ pd->GetPointCells(ptsA->GetId(i), cells);
+ numCells = cells->GetNumberOfIds();
+
+ for (j = 0; j < numCells; j++) {
+ polysA.emplace_back(cells->GetId(j), ptsA->GetId(i));
+ }
+ }
+
+ for (i = 0; i < numPtsB; i++) {
+ pd->GetPointCells(ptsB->GetId(i), cells);
+ numCells = cells->GetNumberOfIds();
+
+ for (j = 0; j < numCells; j++) {
+ polysB.emplace_back(cells->GetId(j), ptsB->GetId(i));
+ }
+ }
+
+ /*for (const Pair &pA : polysA) {
+ std::cout << "pA " << pA << std::endl;
+ }
+
+ for (const Pair &pB : polysB) {
+ std::cout << "pB " << pB << std::endl;
+ }*/
+
+ for (const Pair &pA : polysA) {
+ for (const Pair &pB : polysB) {
+ if (pA.f == pB.f && pd->GetCellType(pA.f) != VTK_EMPTY_CELL) {
+
+ pd->GetCellPoints(pA.f, poly);
+ numPts = poly->GetNumberOfIds();
+
+ newPoly->Reset();
+
+ for (j = 0; j < numPts; j++) {
+ newPoly->InsertNextId(poly->GetId(j));
+
+ idA = poly->GetId(j);
+ idB = j+1 == numPts ? poly->GetId(0) : poly->GetId(j+1);
+
+ if (pA.g == idA
+ && pB.g == idB) {
+
+ for (itrC = itrA+1; itrC != itrB; ++itrC) {
+ // newPoly->InsertNextId(pd->InsertNextLinkedPoint(itrC->get().pt, 1));
+
+ pd->InsertNextLinkedPoint(1);
+
+ newPoly->InsertNextId(pd->GetPoints()->InsertNextPoint(itrC->get().pt));
+ }
+
+ }
+
+ pd->RemoveReferenceToCell(idA, pA.f);
+ }
+
+ pd->DeleteCell(pA.f);
+
+ pd->InsertNextLinkedCell(VTK_POLYGON, newPoly->GetNumberOfIds(), newPoly->GetPointer(0));
+
+ origCellIds->InsertNextValue(origCellIds->GetValue(pA.f));
+
+ break;
+ }
+ }
+ }
+ }
+
+ itrA = itrB;
+ }
+ }
+ }
+
+ cells->Delete();
+ ptsA->Delete();
+ ptsB->Delete();
+ poly->Delete();
+ newPoly->Delete();
+
+ loc->FreeSearchStructure();
+ loc->Delete();
+
+ pd->RemoveDeletedCells();
+
+}
+
+void vtkPolyDataBooleanFilter::MergePoints (vtkPolyData *pd, PolyStripsType &polyStrips) {
+
+#ifdef DEBUG
+ std::cout << "MergePoints()" << std::endl;
+#endif
+
+ pd->BuildCells();
+ pd->BuildLinks();
+
+ contLines->BuildLinks();
+
+ vtkKdTreePointLocator *loc = vtkKdTreePointLocator::New();
+ loc->SetDataSet(pd);
+ loc->BuildLocator();
+
+ PolyStripsType::const_iterator itr;
+ StripsType::const_iterator itr2;
+
+ std::map<vtkIdType, std::set<vtkIdType>> neighPts;
+
+ vtkIdList *pts = vtkIdList::New();
+
+ vtkIdType i, j, numPts;
+
+ for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) {
+ const PStrips &pStrips = itr->second;
+
+ for (itr2 = pStrips.strips.begin(); itr2 != pStrips.strips.end(); ++itr2) {
+ const StripType &strip = *itr2;
+
+ const StripPtR &spA = strip.front(),
+ &spB = strip.back();
+
+ const auto beforeA = pStrips.pts.find((strip.begin()+1)->ind),
+ beforeB = pStrips.pts.find((strip.end()-2)->ind);
+
+ FindPoints(loc, beforeA->second.pt, pts);
+ numPts = pts->GetNumberOfIds();
+
+ for (i = 0; i < numPts; i++) {
+ neighPts[spA.ind].insert(pts->GetId(i));
+ }
+
+ FindPoints(loc, beforeB->second.pt, pts);
+ numPts = pts->GetNumberOfIds();
+
+ for (i = 0; i < numPts; i++) {
+ neighPts[spB.ind].insert(pts->GetId(i));
+ }
+ }
+ }
+
+ double pt[3];
+
+ vtkIdList *polys = vtkIdList::New(),
+ *poly = vtkIdList::New();
+
+ vtkIdType ind, polyId, _numPts, before, after;
+
+ decltype(neighPts)::const_iterator itr3;
+
+ for (itr3 = neighPts.begin(); itr3 != neighPts.end(); ++itr3) {
+ const auto &inds = itr3->second;
+
+ std::map<Point3d, std::vector<Pair>> pairs;
+
+ contLines->GetPoint(itr3->first, pt);
+
+ FindPoints(loc, pt, pts);
+ numPts = pts->GetNumberOfIds();
+
+ for (i = 0; i < numPts; i++) {
+ ind = pts->GetId(i);
+
+ pd->GetPointCells(ind, polys);
+
+ if (polys->GetNumberOfIds() > 0) {
+ polyId = polys->GetId(0);
+
+ pd->GetCellPoints(polyId, poly);
+ _numPts = poly->GetNumberOfIds();
+
+ for (j = 0; j < _numPts; j++) {
+ if (poly->GetId(j) == ind) {
+ break;
+ }
+ }
+
+ // wieder davor und danach ermitteln
+
+ before = poly->GetId(j == 0 ? _numPts-1 : j-1);
+ after = poly->GetId(j+1 == _numPts ? 0 : j+1);
+
+ if (std::find(inds.begin(), inds.end(), before) == inds.end()) {
+ pd->GetPoint(before, pt);
+ pairs[{pt[0], pt[1], pt[2]}].emplace_back(polyId, ind);
+ }
+
+ if (std::find(inds.begin(), inds.end(), after) == inds.end()) {
+ pd->GetPoint(after, pt);
+ pairs[{pt[0], pt[1], pt[2]}].emplace_back(polyId, ind);
+ }
+
+ }
+ }
+
+ std::deque<std::deque<std::reference_wrapper<const Pair>>> Pairs;
+
+ decltype(pairs)::const_iterator itr4;
+
+ for (itr4 = pairs.begin(); itr4 != pairs.end(); ++itr4) {
+ const auto &p = itr4->second;
+
+ if (p.size() == 2) {
+ auto _pts = {std::ref(p.front()), std::ref(p.back())}; // std::initializer_list
+ Pairs.emplace_back(_pts);
+ }
+ }
+
+ decltype(Pairs)::iterator itr5;
+
+ /*for (itr5 = Pairs.begin(); itr5 != Pairs.end(); ++itr5) {
+ for (auto &p : *itr5) {
+ std::cout << p.get() << ", ";
+ }
+ std::cout << std::endl;
+ }*/
+
+ decltype(Pairs)::value_type group;
+
+ decltype(group)::const_iterator itr6;
+
+ while (!Pairs.empty()) {
+ if (group.empty()) {
+ group = Pairs.front();
+ Pairs.pop_front();
+ }
+
+ itr5 = Pairs.begin();
+
+ while (itr5 != Pairs.end()) {
+ const auto &next = *itr5;
+
+ if (next.front().get() == group.front().get()) {
+ group.emplace_front(next.back());
+ Pairs.erase(itr5);
+ itr5 = Pairs.begin();
+ } else if (next.front().get() == group.back().get()) {
+ group.emplace_back(next.back());
+ Pairs.erase(itr5);
+ itr5 = Pairs.begin();
+ } else if (next.back().get() == group.front().get()) {
+ group.emplace_front(next.front());
+ Pairs.erase(itr5);
+ itr5 = Pairs.begin();
+ } else if (next.back().get() == group.back().get()) {
+ group.emplace_back(next.front());
+ Pairs.erase(itr5);
+ itr5 = Pairs.begin();
+ } else {
+ ++itr5;
+ }
+ }
+
+ if (itr5 == Pairs.end()) {
+ for (itr6 = group.begin()+1; itr6 != group.end(); ++itr6) {
+ pd->ReplaceCellPoint(itr6->get().f, itr6->get().g, group.front().get().g);
+ }
+
+ group.clear();
+ }
+ }
+
+ }
+
+ polys->Delete();
+ poly->Delete();
+ pts->Delete();
+
+ loc->FreeSearchStructure();
+ loc->Delete();
+
+}
+
+enum class Congr {
+ EQUAL,
+ OPPOSITE,
+ NOT
+};
+
+class PolyAtEdge {
+ vtkPolyData *pd;
+
+public:
+ PolyAtEdge (vtkPolyData *_pd, vtkIdType _polyId, vtkIdType _ptIdA, vtkIdType _ptIdB) : pd(_pd), polyId(_polyId), ptIdA(_ptIdA), ptIdB(_ptIdB), loc(Loc::NONE) {
+
+ double ptA[3], ptB[3];
+
+ pd->GetPoint(ptIdA, ptA);
+ pd->GetPoint(ptIdB, ptB);
+
+ vtkMath::Subtract(ptB, ptA, e);
+ vtkMath::Normalize(e);
+
+ const vtkIdType *poly;
+
+ vtkIdType numPts;
+ pd->GetCellPoints(polyId, numPts, poly);
+
+ ComputeNormal(pd->GetPoints(), n, numPts, poly);
+
+ vtkMath::Cross(e, n, r);
+
+ }
+
+ vtkIdType polyId, ptIdA, ptIdB;
+ double n[3], e[3], r[3];
+
+ Loc loc;
+
+ friend std::ostream& operator<< (std::ostream &out, const PolyAtEdge &p) {
+ out << "polyId " << p.polyId << ", ptIdA " << p.ptIdA << ", ptIdB " << p.ptIdB;
+ return out;
+ }
+
+ static constexpr double eps {.99999999}; // ~.0081deg
+
+ Congr IsCongruent (const PolyAtEdge &p) const {
+ double cong = vtkMath::Dot(n, p.n);
+
+ if (cong > eps || cong < -eps) {
+ double ang = vtkMath::Dot(r, p.r);
+
+ if (ang > eps) {
+ if (cong > eps) {
+ // normalen sind gleich ausgerichtet
+ return Congr::EQUAL;
+ } else {
+ return Congr::OPPOSITE;
+ }
+ }
+ }
+
+ return Congr::NOT;
+ }
+
+};
+
+
+class PolyPair {
+public:
+ PolyPair (PolyAtEdge _pA, PolyAtEdge _pB) : pA(_pA), pB(_pB) {}
+
+ PolyAtEdge pA, pB;
+
+ void GetLoc (PolyAtEdge &pT, vtkIdType mode) {
+
+ Congr cA = pA.IsCongruent(pT),
+ cB = pB.IsCongruent(pT);
+
+#ifdef DEBUG
+ std::cout << "GetLoc() -> polyId " << pT.polyId
+ << ", cA " << cA
+ << ", cB " << cB
+ << std::endl;
+
+ if (cA != Congr::NOT || cB != Congr::NOT) {
+ assert(cA != cB);
+ }
+
+#endif
+
+ if (cA == Congr::EQUAL || cA == Congr::OPPOSITE) {
+ if (cA == Congr::OPPOSITE) {
+ // normalen sind entgegengesetzt gerichtet
+
+ if (mode == OPER_INTERSECTION) {
+ pA.loc = Loc::OUTSIDE;
+ pT.loc = Loc::OUTSIDE;
+ } else {
+ pA.loc = Loc::INSIDE;
+ pT.loc = Loc::INSIDE;
+ }
+ } else if (mode == OPER_UNION || mode == OPER_INTERSECTION) {
+ pA.loc = Loc::INSIDE;
+ pT.loc = Loc::OUTSIDE;
+ }
+
+ } else if (cB == Congr::EQUAL || cB == Congr::OPPOSITE) {
+ if (cB == Congr::OPPOSITE) {
+ // normalen sind entgegengesetzt gerichtet
+
+ if (mode == OPER_INTERSECTION) {
+ pB.loc = Loc::OUTSIDE;
+ pT.loc = Loc::OUTSIDE;
+ } else {
+ pB.loc = Loc::INSIDE;
+ pT.loc = Loc::INSIDE;
+ }
+ } else if (mode == OPER_UNION || mode == OPER_INTERSECTION) {
+ pB.loc = Loc::INSIDE;
+ pT.loc = Loc::OUTSIDE;
+ }
+
+ } else {
+ double alpha = GetAngle(pA.r, pB.r, pA.e),
+ beta = GetAngle(pA.r, pT.r, pA.e);
+
+ if (beta > alpha) {
+ pT.loc = Loc::INSIDE;
+ } else {
+ pT.loc = Loc::OUTSIDE;
+ }
+ }
+ }
+
+};
+
+
+std::shared_ptr<PolyPair> GetEdgePolys (vtkPolyData *pd, vtkIdList *ptsA, vtkIdList *ptsB) {
+
+#ifdef DEBUG
+ std::cout << "GetEdgePolys()" << std::endl;
+#endif
+
+ std::vector<Pair> p;
+
+ vtkIdType numPtsA = ptsA->GetNumberOfIds(),
+ numPtsB = ptsB->GetNumberOfIds();
+
+ vtkIdList *polys = vtkIdList::New();
+
+ vtkIdType i, j, numCells;
+
+ for (i = 0; i < numPtsA; i++) {
+ pd->GetPointCells(ptsA->GetId(i), polys);
+ numCells = polys->GetNumberOfIds();
+
+ for (j = 0; j < numCells; j++) {
+ p.emplace_back(ptsA->GetId(i), polys->GetId(j));
+ }
+ }
+
+ for (i = 0; i < numPtsB; i++) {
+ pd->GetPointCells(ptsB->GetId(i), polys);
+ numCells = polys->GetNumberOfIds();
+
+ for (j = 0; j < numCells; j++) {
+ p.emplace_back(ptsB->GetId(i), polys->GetId(j));
+ }
+ }
+
+ polys->Delete();
+
+ std::map<vtkIdType, IdsType> pEdges;
+
+ std::vector<Pair>::const_iterator itr;
+ for (itr = p.begin(); itr != p.end(); ++itr) {
+ pEdges[itr->g].push_back(itr->f);
+ }
+
+ std::vector<PolyAtEdge> opp;
+
+ vtkIdList *poly = vtkIdList::New();
+
+ vtkIdType numPts, a, b;
+
+ std::map<vtkIdType, IdsType>::const_iterator itr2;
+
+ for (itr2 = pEdges.begin(); itr2 != pEdges.end(); ++itr2) {
+ const IdsType &pts = itr2->second;
+
+ if (pts.size() > 1) {
+ pd->GetCellPoints(itr2->first, poly);
+ numPts = poly->GetNumberOfIds();
+
+ for (i = 0; i < numPts; i++) {
+ a = poly->GetId(i);
+ b = i+1 == numPts ? poly->GetId(0) : poly->GetId(i+1);
+
+ if (std::find(pts.begin(), pts.end(), a) != pts.end()
+ && std::find(pts.begin(), pts.end(), b) != pts.end()) {
+
+ opp.emplace_back(pd, itr2->first, a, b);
+ }
+ }
+
+ }
+ }
+
+ poly->Delete();
+
+#ifdef DEBUG
+ for (auto &op : opp) {
+ std::cout << op << std::endl;
+ }
+#endif
+
+ if (opp.size() != 2) {
+ return nullptr;
+ }
+
+ return std::make_shared<PolyPair>(opp[0], opp[1]);
+
+}
+
+bool vtkPolyDataBooleanFilter::CombineRegions () {
+
+#ifdef DEBUG
+ std::cout << "CombineRegions()" << std::endl;
+#endif
+
+ vtkSmartPointer<vtkPolyData> filterdA = vtkSmartPointer<vtkPolyData>::New();
+ filterdA->DeepCopy(modPdA);
+
+ vtkSmartPointer<vtkPolyData> filterdB = vtkSmartPointer<vtkPolyData>::New();
+ filterdB->DeepCopy(modPdB);
+
+ // ungenutzte punkte löschen
+ vtkSmartPointer<vtkCleanPolyData> cleanA = vtkSmartPointer<vtkCleanPolyData>::New();
+ cleanA->PointMergingOff();
+ cleanA->SetInputData(filterdA);
+
+ vtkSmartPointer<vtkCleanPolyData> cleanB = vtkSmartPointer<vtkCleanPolyData>::New();
+ cleanB->PointMergingOff();
+ cleanB->SetInputData(filterdB);
+
+ // regionen mit skalaren ausstatten
+ vtkSmartPointer<vtkPolyDataConnectivityFilter> cfA = vtkSmartPointer<vtkPolyDataConnectivityFilter>::New();
+ cfA->SetExtractionModeToAllRegions();
+ cfA->ColorRegionsOn();
+ cfA->SetInputConnection(cleanA->GetOutputPort());
+
+ vtkSmartPointer<vtkPolyDataConnectivityFilter> cfB = vtkSmartPointer<vtkPolyDataConnectivityFilter>::New();
+ cfB->SetExtractionModeToAllRegions();
+ cfB->ColorRegionsOn();
+ cfB->SetInputConnection(cleanB->GetOutputPort());
+
+ cfA->Update();
+ cfB->Update();
+
+ vtkPolyData *pdA = cfA->GetOutput();
+ vtkPolyData *pdB = cfB->GetOutput();
+
+#ifdef DEBUG
+ std::cout << "Exporting modPdA_8.vtk" << std::endl;
+ WriteVTK("modPdA_8.vtk", pdA);
+
+ std::cout << "Exporting modPdB_8.vtk" << std::endl;
+ WriteVTK("modPdB_8.vtk", pdB);
+#endif
+
+ resultC->ShallowCopy(contLines);
+
+ if (OperMode == OPER_NONE) {
+ resultA->ShallowCopy(pdA);
+ resultB->ShallowCopy(pdB);
+
+ return false;
+ }
+
+ vtkSmartPointer<vtkKdTreePointLocator> plA = vtkSmartPointer<vtkKdTreePointLocator>::New();
+ plA->SetDataSet(pdA);
+ plA->BuildLocator();
+
+ vtkSmartPointer<vtkKdTreePointLocator> plB = vtkSmartPointer<vtkKdTreePointLocator>::New();
+ plB->SetDataSet(pdB);
+ plB->BuildLocator();
+
+ pdA->BuildLinks();
+ pdB->BuildLinks();
+
+ vtkIdTypeArray *scalarsA = vtkIdTypeArray::SafeDownCast(pdA->GetPointData()->GetScalars());
+ vtkIdTypeArray *scalarsB = vtkIdTypeArray::SafeDownCast(pdB->GetPointData()->GetScalars());
+
+ vtkSmartPointer<vtkIdList> line = vtkSmartPointer<vtkIdList>::New();
+
+ double ptA[3], ptB[3];
+
+ vtkSmartPointer<vtkIdList> fptsA = vtkSmartPointer<vtkIdList>::New();
+ vtkSmartPointer<vtkIdList> lptsA = vtkSmartPointer<vtkIdList>::New();
+
+ vtkSmartPointer<vtkIdList> fptsB = vtkSmartPointer<vtkIdList>::New();
+ vtkSmartPointer<vtkIdList> lptsB = vtkSmartPointer<vtkIdList>::New();
+
+ std::map<vtkIdType, Loc> locsA, locsB;
+
+ vtkIdType i, j, numLines = contLines->GetNumberOfCells();
+
+ IdsType _failed;
+
+ for (i = 0; i < numLines; i++) {
+
+ if (contLines->GetCellType(i) == VTK_EMPTY_CELL) {
+ continue;
+ }
+
+ contLines->GetCellPoints(i, line);
+
+ contLines->GetPoint(line->GetId(0), ptA);
+ contLines->GetPoint(line->GetId(1), ptB);
+
+ FindPoints(plA, ptA, fptsA);
+ FindPoints(plB, ptA, fptsB);
+
+#ifdef DEBUG
+ std::cout << "line " << i << std::endl;
+#else
+
+ // bereits behandelte regionen werden nicht noch einmal untersucht
+
+ vtkIdType notLocated = 0;
+
+ for (j = 0; j < fptsA->GetNumberOfIds(); j++) {
+ if (locsA.count(scalarsA->GetValue(fptsA->GetId(j))) == 0) {
+ notLocated++;
+ }
+ }
+
+ for (j = 0; j < fptsB->GetNumberOfIds(); j++) {
+ if (locsB.count(scalarsB->GetValue(fptsB->GetId(j))) == 0) {
+ notLocated++;
+ }
+ }
+
+ if (notLocated == 0) {
+ continue;
+ }
+
+#endif
+
+ FindPoints(plA, ptB, lptsA);
+ FindPoints(plB, ptB, lptsB);
+
+ auto ppA = GetEdgePolys(pdA, fptsA, lptsA);
+ auto ppB = GetEdgePolys(pdB, fptsB, lptsB);
+
+ if (ppA && ppB) {
+
+ ppB->GetLoc(ppA->pA, OperMode);
+ ppB->GetLoc(ppA->pB, OperMode);
+
+ ppA->GetLoc(ppB->pA, OperMode);
+ ppA->GetLoc(ppB->pB, OperMode);
+
+ vtkIdType fsA = scalarsA->GetValue(ppA->pA.ptIdA);
+ vtkIdType lsA = scalarsA->GetValue(ppA->pB.ptIdA);
+
+ vtkIdType fsB = scalarsB->GetValue(ppB->pA.ptIdA);
+ vtkIdType lsB = scalarsB->GetValue(ppB->pB.ptIdA);
+
+#ifdef DEBUG
+ std::cout << "polyId " << ppA->pA.polyId << ", sA " << fsA << ", loc " << ppA->pA.loc << std::endl;
+ std::cout << "polyId " << ppA->pB.polyId << ", sA " << lsA << ", loc " << ppA->pB.loc << std::endl;
+ std::cout << "polyId " << ppB->pA.polyId << ", sB " << fsB << ", loc " << ppB->pA.loc << std::endl;
+ std::cout << "polyId " << ppB->pB.polyId << ", sB " << lsB << ", loc " << ppB->pB.loc << std::endl;
+
+ if (locsA.count(fsA) > 0 && locsA[fsA] != ppA->pA.loc) {
+ std::cout << "sA " << fsA << ": " << locsA[fsA] << " -> " << ppA->pA.loc << std::endl;
+ }
+
+ if (locsA.count(lsA) > 0 && locsA[lsA] != ppA->pB.loc) {
+ std::cout << "sA " << lsA << ": " << locsA[lsA] << " -> " << ppA->pB.loc << std::endl;
+ }
+
+ if (locsB.count(fsB) > 0 && locsB[fsB] != ppB->pA.loc) {
+ std::cout << "sB " << fsB << ": " << locsB[fsB] << " -> " << ppB->pA.loc << std::endl;
+ }
+
+ if (locsB.count(lsB) > 0 && locsB[lsB] != ppB->pB.loc) {
+ std::cout << "sB " << lsB << ": " << locsB[lsB] << " -> " << ppB->pB.loc << std::endl;
+ }
+
+#endif
+
+ locsA.emplace(fsA, ppA->pA.loc);
+ locsA.emplace(lsA, ppA->pB.loc);
+
+ locsB.emplace(fsB, ppB->pA.loc);
+ locsB.emplace(lsB, ppB->pB.loc);
+
+ } else {
+ _failed.push_back(i);
+
+ // return true;
+ }
+
+ }
+
+ if (_failed.size() > 0) {
+ for (auto i : _failed) {
+ std::cout << "failed at " << i
+ << std::endl;
+ }
+
+ return true;
+ }
+
+ // reale kombination der ermittelten regionen
+
+ Loc comb[] = {Loc::OUTSIDE, Loc::OUTSIDE};
+
+ if (OperMode == OPER_INTERSECTION) {
+ comb[0] = Loc::INSIDE;
+ comb[1] = Loc::INSIDE;
+ } else if (OperMode == OPER_DIFFERENCE) {
+ comb[1] = Loc::INSIDE;
+ } else if (OperMode == OPER_DIFFERENCE2) {
+ comb[0] = Loc::INSIDE;
+ }
+
+ vtkIdType numA = cfA->GetNumberOfExtractedRegions(),
+ numB = cfB->GetNumberOfExtractedRegions();
+
+ cfA->SetExtractionModeToSpecifiedRegions();
+ cfB->SetExtractionModeToSpecifiedRegions();
+
+ std::map<vtkIdType, Loc>::const_iterator itr;
+
+ for (itr = locsA.begin(); itr != locsA.end(); itr++) {
+ if (itr->second == comb[0]) {
+ cfA->AddSpecifiedRegion(itr->first);
+ }
+ }
+
+ for (itr = locsB.begin(); itr != locsB.end(); itr++) {
+ if (itr->second == comb[1]) {
+ cfB->AddSpecifiedRegion(itr->first);
+ }
+ }
+
+ // nicht beteiligte regionen hinzufügen
+
+ if (OperMode == OPER_UNION || OperMode == OPER_DIFFERENCE) {
+ for (i = 0; i < numA; i++) {
+ if (locsA.count(i) == 0) {
+ cfA->AddSpecifiedRegion(i);
+ }
+ }
+ }
+
+ if (OperMode == OPER_UNION || OperMode == OPER_DIFFERENCE2) {
+ for (i = 0; i < numB; i++) {
+ if (locsB.count(i) == 0) {
+ cfB->AddSpecifiedRegion(i);
+ }
+ }
+ }
+
+ // nach innen zeigende normalen umkehren
+
+ cfA->Update();
+ cfB->Update();
+
+ vtkPolyData *regsA = cfA->GetOutput();
+ vtkPolyData *regsB = cfB->GetOutput();
+
+ scalarsA = vtkIdTypeArray::SafeDownCast(regsA->GetPointData()->GetScalars());
+ scalarsB = vtkIdTypeArray::SafeDownCast(regsB->GetPointData()->GetScalars());
+
+ if (OperMode != OPER_INTERSECTION) {
+ if (comb[0] == Loc::INSIDE) {
+ for (i = 0; i < regsA->GetNumberOfCells(); i++) {
+ if (locsA.count(scalarsA->GetValue(i)) == 1) {
+ regsA->ReverseCell(i);
+ }
+ }
+ }
+
+ if (comb[1] == Loc::INSIDE) {
+ for (i = 0; i < regsB->GetNumberOfCells(); i++) {
+ if (locsB.count(scalarsB->GetValue(i)) == 1) {
+ regsB->ReverseCell(i);
+ }
+ }
+ }
+ }
+
+ // OrigCellIds und CellData
+
+ vtkIdTypeArray *origCellIdsA = vtkIdTypeArray::SafeDownCast(regsA->GetCellData()->GetScalars("OrigCellIds"));
+ vtkIdTypeArray *origCellIdsB = vtkIdTypeArray::SafeDownCast(regsB->GetCellData()->GetScalars("OrigCellIds"));
+
+ vtkIdTypeArray *newOrigCellIdsA = vtkIdTypeArray::New();
+ newOrigCellIdsA->SetName("OrigCellIdsA");
+
+ vtkIdTypeArray *newOrigCellIdsB = vtkIdTypeArray::New();
+ newOrigCellIdsB->SetName("OrigCellIdsB");
+
+ vtkCellData *newCellDataA = vtkCellData::New();
+ newCellDataA->CopyAllocate(cellDataA);
+
+ vtkCellData *newCellDataB = vtkCellData::New();
+ newCellDataB->CopyAllocate(cellDataB);
+
+ vtkIdType cellId;
+
+ for (i = 0; i < regsA->GetNumberOfCells(); i++) {
+ cellId = cellIdsA->GetValue(origCellIdsA->GetValue(i));
+
+ newOrigCellIdsA->InsertNextValue(cellId);
+ newOrigCellIdsB->InsertNextValue(-1);
+
+ newCellDataA->CopyData(cellDataA, cellId, i);
+ }
+
+ for (i = 0; i < regsB->GetNumberOfCells(); i++) {
+ cellId = cellIdsB->GetValue(origCellIdsB->GetValue(i));
+
+ newOrigCellIdsB->InsertNextValue(cellId);
+ newOrigCellIdsA->InsertNextValue(-1);
+
+ newCellDataB->CopyData(cellDataB, cellId, i);
+ }
+
+ regsA->GetCellData()->Initialize();
+ regsB->GetCellData()->Initialize();
+
+ regsA->GetCellData()->ShallowCopy(newCellDataA);
+ regsB->GetCellData()->ShallowCopy(newCellDataB);
+
+ newCellDataA->Delete();
+ newCellDataB->Delete();
+
+ // zusammenführung
+
+ vtkAppendPolyData *app = vtkAppendPolyData::New();
+ app->AddInputData(regsA);
+ app->AddInputData(regsB);
+
+ // entfernt ungenutzte punkte
+ vtkCleanPolyData *cleanApp = vtkCleanPolyData::New();
+ cleanApp->PointMergingOff();
+ cleanApp->SetInputConnection(app->GetOutputPort());
+
+ // färbt die regionen nochmal neu ein, damit mehrere regionen nicht die gleiche farbe haben
+
+ vtkPolyDataConnectivityFilter *cfApp = vtkPolyDataConnectivityFilter::New();
+ cfApp->SetExtractionModeToAllRegions();
+ cfApp->ColorRegionsOn();
+ cfApp->SetInputConnection(cleanApp->GetOutputPort());
+
+ cfApp->Update();
+
+ vtkPolyData *cfPd = cfApp->GetOutput();
+
+ // resultB bleibt hier leer
+
+ resultA->ShallowCopy(cfPd);
+
+ resultA->GetCellData()->AddArray(newOrigCellIdsA);
+ resultA->GetCellData()->AddArray(newOrigCellIdsB);
+
+ // aufräumen
+
+ cfApp->Delete();
+ cleanApp->Delete();
+ app->Delete();
+
+ newOrigCellIdsB->Delete();
+ newOrigCellIdsA->Delete();
+
+ plB->FreeSearchStructure();
+ plA->FreeSearchStructure();
+
+ return false;
+
+}
+
+Merger::Merger (vtkPolyData *pd, const PStrips &pStrips, const StripsType &strips, const IdsType &descIds, vtkIdType origId) : pd(pd), pStrips(pStrips), origId(origId) {
+
+ const StripPtsType &pts = pStrips.pts;
+ const Base &base = pStrips.base;
+
+ vtkIdType i, num;
+ const vtkIdType *cell;
+
+ double pt[3];
+
+ for (auto id : descIds) {
+ pd->GetCellPoints(id, num, cell);
+
+ Poly p;
+
+ for (i = 0; i < num; i++) {
+ pd->GetPoint(cell[i], pt);
+
+ double proj[2];
+ Transform(pt, proj, base);
+
+ p.emplace_back(proj[0], proj[1], 0, cell[i]);
+
+ }
+
+ polys.push_back(p);
+
+ pd->DeleteCell(id);
+ }
+
+ for (auto &strip : strips) {
+ Poly p;
+
+ for (auto &sp : strip) {
+ const double *pt = pts.at(sp.ind).pt;
+
+ double proj[2];
+ Transform(pt, proj, base);
+
+ p.emplace_back(proj[0], proj[1], 0);
+ }
+
+ p.pop_back();
+
+ double n[3];
+ ComputeNormal(p, n);
+
+ if (n[2] < 0) {
+ Poly q(p.rbegin(), p.rend());
+ p.swap(q);
+ }
+
+ polys.push_back(p);
+ }
+
+}
+
+void Merger::Run () {
+ // mergen mit hilfe von vtkKdTree und vtkModifiedBSPTree
+
+ vtkPoints *pdPts = pd->GetPoints();
+ vtkIdTypeArray *origCellIds = vtkIdTypeArray::SafeDownCast(pd->GetCellData()->GetScalars("OrigCellIds"));
+
+ assert(origCellIds != nullptr);
+
+ const Base &base = pStrips.base;
+
+ std::vector<GroupType> groups(polys.size());
+
+ PolysType::const_iterator itrA, itrB;
+
+ for (itrA = polys.begin(); itrA != polys.end(); ++itrA) {
+ for (itrB = polys.begin(); itrB != polys.end(); ++itrB) {
+ if (itrA != itrB && PointInPoly(*itrB, *itrA->begin())) {
+ groups[static_cast<std::size_t>(itrB-polys.begin())].push_back(static_cast<std::size_t>(itrA-polys.begin()));
+ }
+ }
+ }
+
+ std::size_t parent = 0;
+
+ for (auto &group : groups) {
+ GroupType parents;
+
+ for (auto &index : group) {
+ const GroupType &_group = groups[index];
+ parents.insert(parents.end(), _group.begin(), _group.end());
+ }
+
+ std::sort(group.begin(), group.end());
+ std::sort(parents.begin(), parents.end());
+
+ GroupType _group {parent++};
+ std::set_difference(group.begin(), group.end(), parents.begin(), parents.end(), std::back_inserter(_group));
+
+ std::cout << "[";
+ for (auto &index : _group) {
+ std::cout << index << ", ";
+ }
+ std::cout << "]" << std::endl;
+
+ PolysType merged;
+
+ MergeGroup(_group, merged);
+
+ std::map<Point3d, vtkIdType> newIds;
+
+ for (auto &poly : merged) {
+ vtkSmartPointer<vtkIdList> newCell = vtkSmartPointer<vtkIdList>::New();
+
+ for (auto &p : poly) {
+ double in[] = {p.x, p.y},
+ out[3];
+
+ BackTransform(in, out, base);
+
+ vtkIdType id = p.id;
+
+ if (id == NOTSET) {
+ auto itr = newIds.find(p);
+
+ if (itr == newIds.end()) {
+ id = pdPts->InsertNextPoint(out);
+ newIds.emplace(p, id);
+ } else {
+ id = itr->second;
+ }
+ }
+
+ newCell->InsertNextId(id);
+
+ }
+
+ pd->InsertNextCell(VTK_POLYGON, newCell);
+ origCellIds->InsertNextValue(origId);
+ }
+
+ }
+
+}
+
+void Merger::MergeGroup (const GroupType &group, PolysType &merged) {
+ if (group.size() == 1) {
+ merged.push_back(polys.at(group.back()));
+
+ return;
+ }
+
+ vtkSmartPointer<vtkPoints> pts = vtkSmartPointer<vtkPoints>::New();
+
+ IndexedPolysType indexedPolys;
+
+ ReferencedPointsType refPts;
+
+ SourcesType sources;
+ std::size_t src = 0;
+
+ for (auto &index : group) {
+ const Poly &poly = polys.at(index);
+
+ IndexedPoly ids;
+
+ for (auto &p : poly) {
+ vtkIdType id = pts->InsertNextPoint(p.x, p.y, p.z);
+
+ ids.push_back(id);
+ sources.emplace(id, src);
+
+ refPts.emplace(id, p);
+ }
+
+ indexedPolys.push_back(std::move(ids));
+ src++;
+ }
+
+ vtkSmartPointer<vtkKdTree> kdTree = vtkSmartPointer<vtkKdTree>::New();
+ kdTree->OmitZPartitioning();
+ kdTree->BuildLocatorFromPoints(pts);
+
+ vtkSmartPointer<vtkPolyData> linesA = vtkSmartPointer<vtkPolyData>::New();
+ linesA->SetPoints(pts);
+ linesA->Allocate(1);
+
+ IndexedPoly::const_iterator itrA, itrB;
+
+ for (const auto &ids : indexedPolys) {
+ for (itrA = ids.begin(); itrA != ids.end(); ++itrA) {
+ itrB = itrA+1;
+ if (itrB == ids.end()) {
+ itrB = ids.begin();
+ }
+
+ vtkIdList *line = vtkIdList::New();
+ line->InsertNextId(*itrA);
+ line->InsertNextId(*itrB);
+
+ linesA->InsertNextCell(VTK_LINE, line);
+
+ line->Delete();
+ }
+ }
+
+ WriteVTK("linesA.vtk", linesA);
+
+ vtkSmartPointer<vtkModifiedBSPTree> bspTreeA = vtkSmartPointer<vtkModifiedBSPTree>::New();
+ bspTreeA->SetDataSet(linesA);
+
+ int n = 0;
+
+ PolyConnsType polyConns;
+
+ FindConns(linesA, kdTree, bspTreeA, polyConns, indexedPolys, sources, n);
+
+ PolyConnsType connected {{0, {}}};
+ std::set<vtkIdType> restricted; // keine der conns darf im gleichen punkt beginnen
+
+ vtkSmartPointer<vtkPolyData> linesB = vtkSmartPointer<vtkPolyData>::New();
+ linesB->SetPoints(pts);
+ linesB->Allocate(1);
+
+ vtkSmartPointer<vtkModifiedBSPTree> bspTreeB = vtkSmartPointer<vtkModifiedBSPTree>::New();
+ bspTreeB->SetDataSet(linesB);
+
+ ConnsType firstConns;
+
+ std::size_t i, numPolys = indexedPolys.size();
+
+ double ptA[3], ptB[3];
+
+ while (connected.size() < numPolys) {
+
+ bool foundOne = false;
+
+ for (i = 1; i < numPolys; i++) {
+ if (connected.count(i) == 0) {
+ const ConnsType &conns = polyConns[i];
+
+ for (auto &conn : conns) {
+ if (connected.count(sources.at(conn.j)) == 1
+ && restricted.count(conn.j) == 0) {
+
+ pts->GetPoint(conn.i, ptA);
+ pts->GetPoint(conn.j, ptB);
+
+ if (bspTreeB->IntersectWithLine(ptA, ptB, 1e-5, nullptr, nullptr) == 0) {
+ connected[sources.at(conn.i)].push_back(conn);
+
+ // das andere poly auch aktualisieren
+ connected[sources.at(conn.j)].emplace_back(conn.d, conn.j, conn.i);
+
+ restricted.insert(conn.i);
+ restricted.insert(conn.j);
+
+ vtkIdList *line = vtkIdList::New();
+ line->InsertNextId(conn.i);
+ line->InsertNextId(conn.j);
+
+ linesB->InsertNextCell(VTK_LINE, line);
+
+ line->Delete();
+
+ bspTreeB->Modified();
+
+ foundOne = true;
+
+ firstConns.push_back(conn);
+
+ break;
+ }
+
+ }
+ }
+ }
+ }
+
+ if (!foundOne) {
+ if (!FindConns(linesA, kdTree, bspTreeA, polyConns, indexedPolys, sources, n)) {
+ throw std::runtime_error("Merging failed.");
+ }
+ }
+ }
+
+ std::map<std::size_t, std::vector<std::size_t>> chains;
+
+ PolyConnsType::const_iterator itrC;
+
+ for (itrC = connected.begin(); itrC != connected.end(); ++itrC) {
+ auto &chain = chains[itrC->first];
+ chain.push_back(itrC->first);
+
+ while (chain.back() != 0) {
+ chain.push_back(sources.at(connected.at(chain.back()).front().j));
+ }
+ }
+
+ std::cout << connected;
+
+ decltype(chains)::const_iterator itrD;
+
+ for (itrD = chains.begin(); itrD != chains.end(); ++itrD) {
+ std::cout << itrD->first << ": [";
+ for (auto &id : itrD->second) {
+ std::cout << id << ", ";
+ }
+ std::cout << "]" << std::endl;
+ }
+
+ std::set<std::size_t> solved {0};
+
+ std::deque<std::size_t> searchInds;
+
+ for (i = 1; i < numPolys; i++) {
+ if (connected.at(i).size() == 1) {
+ searchInds.push_back(i);
+ }
+ }
+
+ while (!searchInds.empty()) {
+ PriosType prios;
+
+ for (auto ind : searchInds) {
+ PolyPriosType polyPrios;
+
+ std::cout << "ind " << ind << std::endl;
+
+ const Conn &first = connected.at(ind).back();
+
+ for (auto &conn : polyConns.at(ind)) {
+ auto &src = sources.at(conn.j);
+
+ if (polyPrios.count(src) == 1) {
+ continue;
+ }
+
+ if (conn.i != first.i
+ && conn.j != first.j
+ && restricted.count(conn.i) == 0
+ && restricted.count(conn.j) == 0) {
+
+ pts->GetPoint(conn.i, ptA);
+ pts->GetPoint(conn.j, ptB);
+
+ if (bspTreeB->IntersectWithLine(ptA, ptB, 1e-5, nullptr, nullptr) == 0) {
+ auto &chainA = chains.at(ind),
+ &chainB = chains.at(src);
+
+ std::set<std::size_t> _chainA(chainA.begin(), chainA.end()),
+ _chainB(chainB.begin(), chainB.end());
+
+ // gemeinsame eltern
+ std::set<std::size_t> shared;
+
+ std::set_intersection(_chainA.begin(), _chainA.end(), _chainB.begin(), _chainB.end(), std::inserter(shared, shared.end()));
+
+ // gemeinsame eltern müssen sich alle in solved befinden
+ if (std::includes(solved.begin(), solved.end(), shared.begin(), shared.end())) {
+ std::set<std::size_t> solvable;
+
+ std::set_difference(_chainA.begin(), _chainA.end(), solved.begin(), solved.end(), std::inserter(solvable, solvable.end()));
+ std::set_difference(_chainB.begin(), _chainB.end(), solved.begin(), solved.end(), std::inserter(solvable, solvable.end()));
+
+ polyPrios.emplace(std::piecewise_construct,
+ std::forward_as_tuple(src),
+ std::forward_as_tuple(conn, solvable, -conn.d));
+ }
+ }
+ }
+ }
+
+ PolyPriosType::const_iterator itr;
+ for (itr = polyPrios.begin(); itr != polyPrios.end(); ++itr) {
+ prios.insert(itr->second);
+ }
+ }
+
+ if (!prios.empty()) {
+ auto &prio = *prios.rbegin();
+
+ std::cout << "found " << prio << std::endl;
+
+ auto &conns = connected.at(sources.at(prio.conn.i));
+
+ conns.push_back(prio.conn);
+
+ connected.at(sources.at(prio.conn.j)).emplace_back(prio.conn.d, prio.conn.j, prio.conn.i);
+
+ restricted.insert(prio.conn.i);
+ restricted.insert(prio.conn.j);
+
+ vtkIdList *line = vtkIdList::New();
+ line->InsertNextId(prio.conn.i);
+ line->InsertNextId(prio.conn.j);
+
+ linesB->InsertNextCell(VTK_LINE, line);
+
+ line->Delete();
+
+ bspTreeB->Modified();
+
+ solved.insert(prio.solvable.begin(), prio.solvable.end());
+
+ searchInds.erase(std::find(searchInds.begin(), searchInds.end(), sources.at(prio.conn.i)));
+
+ auto itr = std::find(searchInds.begin(), searchInds.end(), sources.at(prio.conn.j));
+
+ if (itr != searchInds.end()) {
+ searchInds.erase(itr);
+ }
+ } else {
+ if (!FindConns(linesA, kdTree, bspTreeA, polyConns, indexedPolys, sources, n)) {
+ throw std::runtime_error("Merging failed.");
+ }
+ }
+ }
+
+ std::cout << connected;
+
+ WriteVTK("linesB.vtk", linesB);
+
+ // stage 1
+
+ struct Cmp {
+ bool operator() (const Conn &a, const Conn &b) const {
+ return std::tie(a.i, a.j) < std::tie(b.i, b.j);
+ }
+ };
+
+ std::set<Conn, Cmp> usedConns;
+
+ IndexedPoly polyA {indexedPolys.front()};
+
+ for (const auto &first : firstConns) {
+ auto itrA = std::find(polyA.begin(), polyA.end(), first.j);
+
+ assert(itrA != polyA.end());
+
+ auto &polyB = indexedPolys.at(sources.at(first.i));
+
+ auto itrB = std::find(polyB.begin(), polyB.end(), first.i);
+
+ assert(itrB != polyB.end());
+
+ std::rotate(polyA.begin(), itrA, polyA.end());
+ std::rotate(polyB.begin(), itrB, polyB.end());
+
+ IndexedPoly newPoly {polyA};
+ newPoly.push_back(polyA.front());
+ newPoly.push_back(polyB.front());
+
+ newPoly.insert(newPoly.end(), polyB.rbegin(), polyB.rend());
+
+ polyA.swap(newPoly);
+
+ usedConns.insert(first);
+
+ }
+
+ PolysType newPolysA;
+ GetPolys(refPts, {polyA}, newPolysA);
+
+ WritePolys("merged_stage1.vtk", newPolysA);
+
+ // stage 2
+
+ std::set<Conn, Cmp> leftConns;
+
+ for (itrC = connected.begin(); itrC != connected.end(); ++itrC) {
+ if (itrC->first == 0) {
+ continue;
+ }
+
+ auto &conns = itrC->second;
+
+ ConnsType::const_iterator itr;
+
+ for (itr = conns.begin()+1; itr != conns.end(); ++itr) {
+ Conn conn(0, itr->j, itr->i);
+
+ if (usedConns.find(conn) == usedConns.end()) {
+ if (itr->i < itr->j) {
+ leftConns.emplace(0, itr->i, itr->j);
+ } else {
+ leftConns.insert(std::move(conn));
+ }
+ }
+ }
+
+ }
+
+ std::cout << "leftConns: [";
+ for (auto &conn : leftConns) {
+ std::cout << conn << ", ";
+ }
+ std::cout << "]" << std::endl;
+
+ IndexedPolysType splitted {polyA};
+
+ for (auto &conn : leftConns) {
+ IndexedPolysType::iterator itr;
+
+ for (itr = splitted.begin(); itr != splitted.end(); ++itr) {
+ auto &poly = *itr;
+
+ auto itrA = std::find(poly.begin(), poly.end(), conn.i);
+
+ if (itrA != poly.end()) {
+ std::rotate(poly.begin(), itrA, poly.end());
+
+ auto itrB = std::find(poly.begin(), poly.end(), conn.j);
+
+ IndexedPoly newPolyA(poly.begin(), itrB+1);
+ IndexedPoly newPolyB(itrB, poly.end());
+
+ newPolyB.push_back(poly.front());
+
+ splitted.erase(itr);
+
+ splitted.push_back(std::move(newPolyA));
+ splitted.push_back(std::move(newPolyB));
+
+ break;
+ }
+ }
+ }
+
+ PolysType newPolysB;
+ GetPolys(refPts, splitted, newPolysB);
+
+ WritePolys("merged_stage2.vtk", newPolysB);
+
+ std::move(newPolysB.begin(), newPolysB.end(), std::back_inserter(merged));
+
+}
+
+bool Merger::FindConns (vtkPolyData *lines, vtkSmartPointer<vtkKdTree> kdTree, vtkSmartPointer<vtkModifiedBSPTree> bspTree, PolyConnsType &polyConns, const IndexedPolysType &indexedPolys, const SourcesType &sources, int &n) {
+
+ vtkPoints *pts = lines->GetPoints();
+
+ if (n > pts->GetNumberOfPoints()) {
+ return false;
+ }
+
+ n += 10;
+
+ vtkSmartPointer<vtkIdList> foundPts = vtkSmartPointer<vtkIdList>::New();
+
+ vtkIdType i, numPts;
+
+ vtkIdType idB;
+
+ vtkSmartPointer<vtkIdList> lineIds = vtkSmartPointer<vtkIdList>::New();
+
+ double ptA[3], ptB[3];
+
+ bool good;
+
+ vtkIdType j;
+ vtkIdType _idA, _idB;
+
+ std::map<std::size_t, std::set<Conn, ConnCmp>> _polyConns;
+
+ vtkSmartPointer<vtkIdList> line = vtkSmartPointer<vtkIdList>::New();
+
+ for (const auto &ids : indexedPolys) {
+ for (vtkIdType idA : ids) {
+ pts->GetPoint(idA, ptA);
+
+ kdTree->FindClosestNPoints(n, ptA, foundPts);
+
+ numPts = foundPts->GetNumberOfIds();
+
+ for (i = 0; i < numPts; i++) {
+ idB = foundPts->GetId(i);
+
+ auto srcA = sources.at(idA),
+ srcB = sources.at(idB);
+
+ if (srcA == srcB) {
+ continue;
+ }
+
+ pts->GetPoint(idB, ptB);
+
+ good = true;
+
+ if (bspTree->IntersectWithLine(ptA, ptB, 1e-5, nullptr, lineIds) == 1) {
+ for (j = 0; j < lineIds->GetNumberOfIds(); j++) {
+ lines->GetCellPoints(lineIds->GetId(j), line);
+
+ _idA = line->GetId(0);
+ _idB = line->GetId(1);
+
+ if (_idA != idA && _idA != idB
+ && _idB != idA && _idB != idB) {
+
+ good = false;
+ break;
+ }
+ }
+ }
+
+ if (good) {
+ double d = vtkMath::Distance2BetweenPoints(ptA, ptB);
+
+ _polyConns[srcA].emplace(d, idA, idB);
+ _polyConns[srcB].emplace(d, idB, idA);
+ }
+ }
+ }
+ }
+
+ decltype(_polyConns)::const_iterator itr;
+
+ for (itr = _polyConns.begin(); itr != _polyConns.end(); ++itr) {
+ auto &_conns = itr->second;
+
+ ConnsType conns(_conns.begin(), _conns.end());
+ std::sort(conns.begin(), conns.end());
+
+ polyConns[itr->first].swap(conns);
+
+ }
+
+ return true;
+}
--- /dev/null
+/*
+Copyright 2012-2022 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.
+You may obtain a copy of the License at
+
+ http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+#include <cmath>
+#include <vector>
+#include <map>
+#include <set>
+#include <algorithm>
+
+#include <vtkInformation.h>
+#include <vtkInformationVector.h>
+#include <vtkDemandDrivenPipeline.h>
+#include <vtkObjectFactory.h>
+#include <vtkPolyDataAlgorithm.h>
+#include <vtkPolyData.h>
+#include <vtkOBBTree.h>
+#include <vtkMatrix4x4.h>
+#include <vtkIdList.h>
+#include <vtkPoints.h>
+#include <vtkMath.h>
+#include <vtkIdTypeArray.h>
+#include <vtkCellData.h>
+#include <vtkPointData.h>
+#include <vtkCleanPolyData.h>
+#include <vtkTriangleStrip.h>
+#include <vtkSmartPointer.h>
+#include <vtkPolyDataConnectivityFilter.h>
+#include <vtkFeatureEdges.h>
+#include <vtkCellIterator.h>
+#include <vtkCellArrayIterator.h>
+
+#include <vtkCellArray.h>
+
+#include "vtkPolyDataContactFilter.h"
+#include "Utilities.h"
+
+#undef DEBUG
+
+vtkStandardNewMacro(vtkPolyDataContactFilter);
+
+vtkPolyDataContactFilter::vtkPolyDataContactFilter () {
+
+ contLines = vtkPolyData::New();
+ contLines->Allocate(1000);
+
+ contPts = vtkPoints::New();
+ contPts->SetDataTypeToDouble();
+ contLines->SetPoints(contPts);
+
+ contA = vtkIdTypeArray::New();
+ contB = vtkIdTypeArray::New();
+
+ contA->SetName("cA");
+ contB->SetName("cB");
+
+ sourcesA = vtkIdTypeArray::New();
+ sourcesA->SetNumberOfComponents(2);
+
+ sourcesB = vtkIdTypeArray::New();
+ sourcesB->SetNumberOfComponents(2);
+
+ sourcesA->SetName("sourcesA");
+ sourcesB->SetName("sourcesB");
+
+ neigsA = vtkIdTypeArray::New();
+ neigsB = vtkIdTypeArray::New();
+
+ neigsA->SetName("neigsA");
+ neigsB->SetName("neigsB");
+
+ SetNumberOfInputPorts(2);
+ SetNumberOfOutputPorts(3);
+
+ invalidA = false;
+ invalidB = false;
+
+ accuracy = vtkIdTypeArray::New();
+ accuracy->SetName("accuracy");
+}
+
+vtkPolyDataContactFilter::~vtkPolyDataContactFilter () {
+ neigsB->Delete();
+ neigsA->Delete();
+
+ sourcesB->Delete();
+ sourcesA->Delete();
+
+ contB->Delete();
+ contA->Delete();
+
+ contPts->Delete();
+ contLines->Delete();
+
+ accuracy->Delete();
+
+}
+
+int vtkPolyDataContactFilter::ProcessRequest (vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector) {
+
+ if (request->Has(vtkDemandDrivenPipeline::REQUEST_DATA())) {
+
+ vtkInformation *inInfoA = inputVector[0]->GetInformationObject(0);
+ vtkInformation *inInfoB = inputVector[1]->GetInformationObject(0);
+
+ vtkPolyData *_pdA = vtkPolyData::SafeDownCast(inInfoA->Get(vtkDataObject::DATA_OBJECT()));
+ vtkPolyData *_pdB = vtkPolyData::SafeDownCast(inInfoB->Get(vtkDataObject::DATA_OBJECT()));
+
+ vtkInformation *outInfoA = outputVector->GetInformationObject(0);
+ vtkInformation *outInfoB = outputVector->GetInformationObject(1);
+ vtkInformation *outInfoC = outputVector->GetInformationObject(2);
+
+ vtkPolyData *resultA = vtkPolyData::SafeDownCast(outInfoA->Get(vtkDataObject::DATA_OBJECT()));
+ vtkPolyData *resultB = vtkPolyData::SafeDownCast(outInfoB->Get(vtkDataObject::DATA_OBJECT()));
+ vtkPolyData *resultC = vtkPolyData::SafeDownCast(outInfoC->Get(vtkDataObject::DATA_OBJECT()));
+
+ // durchführung der aufgabe
+
+ pdA = vtkPolyData::New();
+ pdA->DeepCopy(_pdA);
+
+ pdB = vtkPolyData::New();
+ pdB->DeepCopy(_pdB);
+
+ PreparePolyData(pdA);
+ PreparePolyData(pdB);
+
+ if (pdA->GetNumberOfCells() == 0) {
+ vtkErrorMacro("First input does not contain any supported cells.");
+ return 1;
+ }
+
+ if (pdB->GetNumberOfCells() == 0) {
+ vtkErrorMacro("Second input does not contain any supported cells.");
+ return 1;
+ }
+
+ GetInvalidEdges(pdA, edgesA);
+ GetInvalidEdges(pdB, edgesB);
+
+ // anlegen der obb-trees
+
+ vtkOBBTree *obbA = vtkOBBTree::New();
+ obbA->SetDataSet(pdA);
+ obbA->SetNumberOfCellsPerNode(1);
+ obbA->BuildLocator();
+
+ vtkOBBTree *obbB = vtkOBBTree::New();
+ obbB->SetDataSet(pdB);
+ obbB->SetNumberOfCellsPerNode(1);
+ obbB->BuildLocator();
+
+ vtkMatrix4x4 *mat = vtkMatrix4x4::New();
+
+ obbA->IntersectWithOBBTree(obbB, mat, InterOBBNodes, this);
+
+ if (invalidA) {
+ vtkErrorMacro("First input has non-manifold edges.");
+ return 1;
+ }
+
+ if (invalidB) {
+ vtkErrorMacro("Second input has non-manifold edges.");
+ return 1;
+ }
+
+ contLines->GetCellData()->AddArray(contA);
+ contLines->GetCellData()->AddArray(contB);
+
+ contLines->GetCellData()->AddArray(sourcesA);
+ contLines->GetCellData()->AddArray(sourcesB);
+
+ contLines->GetPointData()->AddArray(accuracy);
+
+ contLines->RemoveDeletedCells();
+
+ vtkCleanPolyData *clean = vtkCleanPolyData::New();
+ clean->SetInputData(contLines);
+ clean->ToleranceIsAbsoluteOn();
+ clean->SetAbsoluteTolerance(1e-5);
+ clean->Update();
+
+ resultA->DeepCopy(clean->GetOutput());
+
+ vtkIdType i, numCellsA = resultA->GetNumberOfCells();
+
+ for (i = 0; i < numCellsA; i++) {
+ if (resultA->GetCellType(i) != VTK_LINE) {
+ resultA->DeleteCell(i);
+ }
+ }
+
+ resultA->RemoveDeletedCells();
+
+ clean->Delete();
+ mat->Delete();
+ obbB->Delete();
+ obbA->Delete();
+
+ resultB->DeepCopy(pdA);
+ resultC->DeepCopy(pdB);
+
+ pdB->Delete();
+ pdA->Delete();
+
+ }
+
+ return 1;
+
+}
+
+void vtkPolyDataContactFilter::PreparePolyData (vtkPolyData *pd) {
+
+ pd->GetCellData()->Initialize();
+ pd->GetPointData()->Initialize();
+
+ vtkIdTypeArray *cellIds = vtkIdTypeArray::New();
+
+ vtkCellIterator *cellItr = pd->NewCellIterator();
+ for (cellItr->InitTraversal(); !cellItr->IsDoneWithTraversal(); cellItr->GoToNextCell()) {
+ cellIds->InsertNextValue(cellItr->GetCellId());
+ }
+
+ vtkIdType cellId;
+
+ for (cellItr->InitTraversal(); !cellItr->IsDoneWithTraversal(); cellItr->GoToNextCell()) {
+ cellId = cellItr->GetCellId();
+
+ if (cellItr->GetCellType() == VTK_QUAD) {
+ vtkIdList *ptIds = cellItr->GetPointIds();
+
+ vtkPoints *pts = cellItr->GetPoints();
+
+ double n[3];
+
+ ComputeNormal(pd->GetPoints(), n, 4, ptIds->GetPointer(0));
+
+ double dA = vtkMath::Dot(n, pts->GetPoint(0)),
+ dB = vtkMath::Dot(n, pts->GetPoint(1))-dA;
+
+ if (std::abs(dB) > 1e-6) {
+ // nur wenn nicht auf einer ebene
+
+ dA = vtkMath::Distance2BetweenPoints(pts->GetPoint(0), pts->GetPoint(2));
+ dB = vtkMath::Distance2BetweenPoints(pts->GetPoint(1), pts->GetPoint(3));
+
+ vtkSmartPointer<vtkIdList> newCellA = vtkSmartPointer<vtkIdList>::New();
+ vtkSmartPointer<vtkIdList> newCellB = vtkSmartPointer<vtkIdList>::New();
+
+ newCellA->SetNumberOfIds(3);
+ newCellB->SetNumberOfIds(3);
+
+ if (dA < dB) {
+ newCellA->SetId(0, ptIds->GetId(1));
+ newCellA->SetId(1, ptIds->GetId(2));
+ newCellA->SetId(2, ptIds->GetId(3));
+
+ newCellB->SetId(0, ptIds->GetId(3));
+ newCellB->SetId(1, ptIds->GetId(0));
+ newCellB->SetId(2, ptIds->GetId(1));
+ } else {
+ newCellA->SetId(0, ptIds->GetId(0));
+ newCellA->SetId(1, ptIds->GetId(1));
+ newCellA->SetId(2, ptIds->GetId(2));
+
+ newCellB->SetId(0, ptIds->GetId(2));
+ newCellB->SetId(1, ptIds->GetId(3));
+ newCellB->SetId(2, ptIds->GetId(0));
+ }
+
+ pd->InsertNextCell(VTK_TRIANGLE, newCellA);
+ cellIds->InsertNextValue(cellId);
+
+ pd->InsertNextCell(VTK_TRIANGLE, newCellB);
+ cellIds->InsertNextValue(cellId);
+
+ pd->DeleteCell(cellId);
+ }
+
+ } else if (cellItr->GetCellType() == VTK_TRIANGLE_STRIP) {
+ vtkIdList *ptIds = cellItr->GetPointIds();
+
+ vtkCellArray *cells = vtkCellArray::New();
+
+ vtkTriangleStrip::DecomposeStrip(cellItr->GetNumberOfPoints(), ptIds->GetPointer(0), cells);
+
+ vtkIdType n;
+ const vtkIdType *pts;
+
+ for (cells->InitTraversal(); cells->GetNextCell(n, pts);) {
+ if (pts[0] != pts[1] && pts[1] != pts[2] && pts[2] != pts[0]) {
+ pd->InsertNextCell(VTK_TRIANGLE, n, pts);
+ cellIds->InsertNextValue(cellId);
+ }
+
+ }
+
+ cells->Delete();
+
+ pd->DeleteCell(cellId);
+
+ } else if (cellItr->GetCellType() != VTK_TRIANGLE && cellItr->GetCellType() != VTK_POLYGON) {
+ pd->DeleteCell(cellId);
+ }
+ }
+
+ cellItr->Delete();
+
+ cellIds->SetName("OrigCellIds");
+ pd->GetCellData()->SetScalars(cellIds);
+
+ cellIds->Delete();
+
+ pd->RemoveDeletedCells();
+ pd->BuildLinks();
+
+}
+
+void vtkPolyDataContactFilter::GetInvalidEdges (vtkPolyData *pd, InvalidEdgesType &edges) {
+ vtkSmartPointer<vtkFeatureEdges> features = vtkSmartPointer<vtkFeatureEdges>::New();
+ features->SetInputData(pd);
+
+ features->BoundaryEdgesOff();
+ features->FeatureEdgesOff();
+ features->ManifoldEdgesOff();
+ features->NonManifoldEdgesOn();
+
+ features->Update();
+
+ vtkPolyData *lines = features->GetOutput();
+
+ vtkIdType num;
+ const vtkIdType *line;
+
+ double ptA[3], ptB[3];
+
+ vtkIdType idA, idB;
+
+ auto lineItr = vtk::TakeSmartPointer(lines->GetLines()->NewIterator());
+
+ for (lineItr->GoToFirstCell(); !lineItr->IsDoneWithTraversal(); lineItr->GoToNextCell()) {
+ lineItr->GetCurrentCell(num, line);
+
+ lines->GetPoint(line[0], ptA);
+ lines->GetPoint(line[1], ptB);
+
+ idA = pd->FindPoint(ptA);
+ idB = pd->FindPoint(ptB);
+
+ edges.emplace(idA, idB);
+ edges.emplace(idB, idA);
+
+ }
+}
+
+void vtkPolyDataContactFilter::InterEdgeLine (InterPtsType &interPts, vtkPolyData *pd, vtkIdType idA, vtkIdType idB, const double *r, const double *ptA, Src src) {
+ double eA[3], eB[3];
+
+ pd->GetPoint(idA, eA);
+ pd->GetPoint(idB, eB);
+
+ double ptB[3];
+ vtkMath::Add(ptA, r, ptB);
+
+ // richtungsvektor der kante bestimmen
+
+ double e[3];
+ vtkMath::Subtract(eB, eA, e);
+ double l = vtkMath::Normalize(e);
+
+ double p[3];
+ vtkMath::Subtract(eA, ptA, p);
+
+ double w = std::abs(vtkMath::Determinant3x3(r, e, p));
+
+ if (w < 1e-4) { // ~89.995deg
+ // schnittpunkt ermitteln
+
+ double v[3];
+ vtkMath::Cross(r, e, v);
+
+ double n = vtkMath::Norm(v);
+
+ if (n > 1e-4) { // ~0.0057deg
+
+ 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;
+
+ if (s > -1e-6 && s < 1e-6) {
+ end = End::BEGIN;
+ } else if (s > l-1e-6 && s < l+1e-6) {
+ end = End::END;
+ }
+
+ interPts.emplace_back(ptA[0]+t*r[0], ptA[1]+t*r[1], ptA[2]+t*r[2], t, idA, idB, end, src);
+
+ }
+
+ } else {
+ // parallel
+
+ double vA[3], vB[3], cA[3], cB[3], dA, dB;
+
+ vtkMath::Subtract(eA, ptA, vA);
+ vtkMath::Subtract(eA, ptB, vB);
+ vtkMath::Cross(vA, vB, cA);
+
+ double dotA = vtkMath::Dot(vA, r);
+
+ vtkMath::Subtract(eB, ptA, vA);
+ vtkMath::Subtract(eB, ptB, vB);
+ vtkMath::Cross(vA, vB, cB);
+
+ double dotB = vtkMath::Dot(vA, r);
+
+ dA = vtkMath::Norm(cA);
+ dB = vtkMath::Norm(cB);
+
+ if (dA < 1e-4 || dB < 1e-4) {
+#ifdef DEBUG
+ std::cout << "congruent lines with d (" << dA << ", " << dB << ") and t (" << dotA << ", " << dotB << ") and l " << l << 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);
+
+ }
+
+ }
+
+ } else {
+ // windschief
+ }
+
+}
+
+void 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;
+#endif
+
+ std::vector<Pair> edges;
+ edges.reserve(static_cast<std::size_t>(num));
+
+ vtkIdType i, j;
+
+ for (i = 0; i < num; i++) {
+ j = i+1;
+
+ if (j == num) {
+ j = 0;
+ }
+
+ const vtkIdType &a = poly[i],
+ &b = poly[j];
+
+ vtkPolyDataContactFilter::InterEdgeLine(interPts, pd, a, b, r, pt, src);
+
+ edges.emplace_back(a, b);
+
+ }
+
+ if (interPts.empty()) {
+ return;
+ }
+
+ struct Cmp {
+ bool operator() (const double &l, const double &r) const {
+ long a = std::lround(l*1e5),
+ b = std::lround(r*1e5);
+
+ return a < b;
+ }
+ };
+
+ std::map<double, InterPtsType, Cmp> paired;
+
+ for (auto &p : interPts) {
+ paired[p.t].push_back(p);
+ }
+
+ std::vector<InterPtsType> _paired;
+
+ for (auto &p : paired) {
+ InterPtsType &pts = p.second;
+
+ if (pts.size() == 1 && pts.front().end != End::NONE) {
+ // hier fehlt der zweite punkt
+ pts.push_back(pts.back());
+ }
+
+ _paired.push_back(pts);
+ }
+
+ // trivial
+
+ if (_paired.front().size() == 2) {
+ _paired.front().pop_back();
+ }
+
+ if (_paired.back().size() == 2) {
+ _paired.back().pop_back();
+ }
+
+ // ...
+
+ std::map<vtkIdType, double> ends;
+
+ for (const auto &pts : _paired) {
+ auto &last = pts.back();
+
+ if (last.end != End::NONE) {
+ ends.emplace(last.end == End::BEGIN ? last.edge.f : last.edge.g, last.t);
+ }
+ }
+
+ double s[3], d;
+
+ vtkMath::Cross(n, r, s);
+ d = vtkMath::Dot(s, pt);
+
+ double ptA[3], ptB[3], dA, dB, vA[3], vB[3], tA, tB;
+
+ vtkIdType id, prev, next;
+
+ for (auto &pts : _paired) {
+ auto &last = pts.back();
+
+ if (last.end != End::NONE) {
+ if (last.end == End::BEGIN) {
+ id = last.edge.f;
+ next = last.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;
+ next = std::find_if(edges.begin(), edges.end(), [&id](const Pair &edge) { return edge.f == id; })->g;
+ }
+
+ if (pts.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)) {
+ // tasche
+ pts.pop_back();
+ }
+
+ continue;
+
+ } else if (ends.count(next) == 1 && ends.count(prev) == 0) {
+ 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)) {
+ // tasche
+ pts.pop_back();
+ }
+
+ continue;
+
+ }
+ }
+
+ if (ends.count(prev) == 0 && ends.count(next) == 0) {
+ pd->GetPoint(next, ptA);
+ pd->GetPoint(prev, ptB);
+
+ dA = vtkMath::Dot(s, ptA)-d;
+ dB = vtkMath::Dot(s, ptB)-d;
+
+ if (std::signbit(dA) != std::signbit(dB)) {
+ if (pts.size() == 2) {
+ pts.pop_back();
+ }
+
+ } else {
+ vtkMath::Subtract(ptA, pt, vA);
+ vtkMath::Subtract(ptB, pt, vB);
+
+ tA = vtkMath::Dot(vA, r);
+ tB = vtkMath::Dot(vB, r);
+
+ if ((tB > tA) == std::signbit(dA)) {
+ pts.clear();
+ }
+
+ }
+ }
+ }
+ }
+
+ // ...
+
+ InterPtsType _interPts;
+
+ for (const auto &pts : _paired) {
+ _interPts.insert(_interPts.end(), pts.begin(), pts.end());
+ }
+
+ interPts.swap(_interPts);
+
+}
+
+void vtkPolyDataContactFilter::InterPolys (vtkIdType idA, vtkIdType idB) {
+
+#ifdef DEBUG
+ std::cout << "InterPolys() -> idA " << idA << ", idB " << idB << std::endl;
+#endif
+
+ vtkIdType numA, numB;
+ const vtkIdType *polyA, *polyB;
+
+ pdA->GetCellPoints(idA, numA, polyA);
+ pdB->GetCellPoints(idB, numB, polyB);
+
+ // ebenen aufstellen
+
+ double nA[3], nB[3], ptA[3], ptB[3], dA, dB;
+
+ ComputeNormal(pdA->GetPoints(), nA, numA, polyA);
+ ComputeNormal(pdB->GetPoints(), nB, numB, polyB);
+
+ pdA->GetPoint(polyA[0], ptA);
+ pdB->GetPoint(polyB[0], ptB);
+
+ dA = vtkMath::Dot(nA, ptA);
+ dB = vtkMath::Dot(nB, ptB);
+
+ // richtungsvektor
+
+ double r[3];
+ 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;
+
+ for (int j = 1; j < 3; j++) {
+ if (std::abs(r[j]) > std::abs(r[i])) {
+ i = j;
+ }
+ }
+
+ int inds[] = {1, 2};
+
+ if (i == 1) {
+ inds[0] = 0; inds[1] = 2;
+ } else if (i == 2) {
+ inds[0] = 0; inds[1] = 1;
+ }
+
+ double det = nA[inds[0]]*nB[inds[1]]-nB[inds[0]]*nA[inds[1]];
+
+ if (std::abs(det) < 1e-12) {
+ return;
+ }
+
+ // ein punkt auf der schnittgeraden der beiden ebenen
+
+ double s[3];
+ s[inds[0]] = (dA*nB[inds[1]]-dB*nA[inds[1]])/det;
+ 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;
+#endif
+
+ InterPtsType intersA, intersB;
+
+ vtkPolyDataContactFilter::InterPolyLine(intersA, pdA, numA, polyA, r, s, Src::A, nA);
+ vtkPolyDataContactFilter::InterPolyLine(intersB, pdB, numB, polyB, r, s, Src::B, nB);
+
+ // probe, ob die schnittpunkte auf den kanten liegen
+ // bei ungenauen normalen ist das manchmal nicht der fall
+
+ vtkPolyDataContactFilter::CheckInters(intersA, pdA, idA, idB);
+ vtkPolyDataContactFilter::CheckInters(intersB, pdB, idA, idB);
+
+#ifdef DEBUG
+ std::cout << "intersA " << intersA.size()
+ << ", intersB " << intersB.size()
+ << std::endl;
+#endif
+
+
+ if (intersA.size() != 0 && intersB.size() != 0
+ && intersA.size()%2 == 0 && intersB.size()%2 == 0) {
+
+ AddContactLines(intersA, intersB, idA, idB);
+ }
+
+}
+
+void vtkPolyDataContactFilter::OverlapLines (OverlapsType &ols, InterPtsType &intersA, InterPtsType &intersB, vtkIdType idA, vtkIdType idB) {
+
+ auto GetNeig = [](const InterPt &pA, const InterPt &pB, vtkPolyData *pd, vtkIdType polyId) -> vtkIdType {
+ if (pA.edge == pB.edge) {
+ vtkSmartPointer<vtkIdList> neigs = vtkSmartPointer<vtkIdList>::New();
+
+ pd->GetCellEdgeNeighbors(polyId, pA.edge.f, pA.edge.g, neigs);
+
+ assert(neigs->GetNumberOfIds() == 1);
+
+ return neigs->GetId(0);
+ }
+
+ return NOTSET;
+ };
+
+ auto Add = [](InterPt &a, InterPt &b, InterPt &c, InterPt &d, vtkIdType neigA, vtkIdType neigB) {
+ a.Merge(c);
+ b.Merge(d);
+
+ return std::make_tuple(a, b, neigA, neigB);
+ };
+
+ InterPtsType::iterator itr, itr2;
+
+ vtkIdType neigA, neigB;
+
+ for (itr = intersA.begin(); itr != intersA.end(); itr += 2) {
+ neigA = GetNeig(*itr, *(itr+1), pdA, idA);
+
+ for (itr2 = intersB.begin(); itr2 != intersB.end(); itr2 += 2) {
+ neigB = GetNeig(*itr2, *(itr2+1), pdB, idB);
+
+ if (itr->t <= itr2->t && (itr+1)->t > itr2->t) {
+ if ((itr2+1)->t < (itr+1)->t) {
+ ols.push_back(Add(*itr2, *(itr2+1), *itr, *(itr+1), neigA, neigB));
+ } else {
+ ols.push_back(Add(*itr2, *(itr+1), *itr, *(itr2+1), neigA, neigB));
+ }
+ } else if (itr2->t <= itr->t && (itr2+1)->t > itr->t) {
+ if ((itr+1)->t < (itr2+1)->t) {
+ ols.push_back(Add(*itr, *(itr+1), *itr2, *(itr2+1), neigA, neigB));
+ } else {
+ ols.push_back(Add(*itr, *(itr2+1), *itr2, *(itr+1), neigA, neigB));
+ }
+ }
+ }
+ }
+
+}
+
+void vtkPolyDataContactFilter::AddContactLines (InterPtsType &intersA, InterPtsType &intersB, vtkIdType idA, vtkIdType idB) {
+
+ OverlapsType overlaps;
+ OverlapLines(overlaps, intersA, intersB, idA, idB);
+
+ OverlapsType::const_iterator itr;
+
+ for (itr = overlaps.begin(); itr != overlaps.end(); ++itr) {
+ 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;
+ }
+ }
+
+ if (s.src == Src::A) {
+ if (edgesA.count(s.edge) == 1) {
+ invalidA = true;
+ }
+ }
+
+ if (f.src == Src::B) {
+ if (edgesB.count(f.edge) == 1) {
+ invalidB = true;
+ }
+ }
+
+ if (s.src == Src::B) {
+ if (edgesB.count(s.edge) == 1) {
+ invalidB = true;
+ }
+ }
+
+ vtkIdList *linePts = vtkIdList::New();
+
+ linePts->InsertNextId(contPts->InsertNextPoint(f.pt));
+ linePts->InsertNextId(contPts->InsertNextPoint(s.pt));
+
+ contLines->InsertNextCell(VTK_LINE, linePts);
+
+ const vtkIdType tupleA[] = {f.srcA, s.srcA};
+ const vtkIdType tupleB[] = {f.srcB, s.srcB};
+
+ sourcesA->InsertNextTypedTuple(tupleA);
+ sourcesB->InsertNextTypedTuple(tupleB);
+
+ linePts->Delete();
+
+ contA->InsertNextValue(idA);
+ contB->InsertNextValue(idB);
+
+ neigsA->InsertNextValue(std::get<2>(*itr));
+ neigsB->InsertNextValue(std::get<3>(*itr));
+
+ accuracy->InsertNextValue(f.inaccurate);
+ accuracy->InsertNextValue(s.inaccurate);
+
+ }
+
+}
+
+int vtkPolyDataContactFilter::InterOBBNodes (vtkOBBNode *nodeA, vtkOBBNode *nodeB, vtkMatrix4x4 *vtkNotUsed(mat), void *caller) {
+ vtkPolyDataContactFilter *self = reinterpret_cast<vtkPolyDataContactFilter*>(caller);
+
+ vtkIdList *cellsA = nodeA->Cells;
+ vtkIdList *cellsB = nodeB->Cells;
+
+ vtkIdType numCellsA = cellsA->GetNumberOfIds(),
+ numCellsB = cellsB->GetNumberOfIds();
+
+ vtkIdType i, j, ci, cj;
+
+ for (i = 0; i < numCellsA; i++) {
+ ci = cellsA->GetId(i);
+
+ for (j = 0; j < numCellsB; j++) {
+ cj = cellsB->GetId(j);
+
+ self->InterPolys(ci, cj);
+ }
+ }
+
+ return 0;
+}
+
+void vtkPolyDataContactFilter::CheckInters (InterPtsType &interPts, vtkPolyData *pd, vtkIdType idA, vtkIdType idB) {
+ double ptA[3],
+ ptB[3],
+ v[3],
+ w[3],
+ k,
+ l,
+ alpha,
+ d;
+
+ for (auto &p : interPts) {
+ pd->GetPoint(p.edge.f, ptA);
+ pd->GetPoint(p.edge.g, ptB);
+
+ vtkMath::Subtract(ptA, ptB, v);
+ vtkMath::Normalize(v);
+ vtkMath::Subtract(ptA, p.pt, w);
+
+ k = vtkMath::Norm(w);
+ l = vtkMath::Dot(v, w);
+ alpha = std::acos(l/k);
+
+ if (std::isnan(alpha)) {
+ continue;
+ }
+
+ d = std::sin(alpha)*k;
+
+ if (d < 1e-5) {
+ continue;
+ }
+
+ // if (p.src == Src::A) {
+ // std::cout << "? A ";
+ // } else {
+ // std::cout << "? B ";
+ // }
+
+ // std::cout << idA << ", " << idB << ": "
+ // << std::fixed
+ // << d
+ // << ", [" << p.pt[0] << ", " << p.pt[1] << ", " << p.pt[2] << "], "
+ // << p.edge
+ // << std::endl;
+
+ p.inaccurate = true;
+ }
+
+}
git://git.creatis.insa-lyon.fr / creaVtk.git / commitdiff