2 Copyright 2012-2022 Ronald Römer
4 Licensed under the Apache License, Version 2.0 (the "License");
5 you may not use this file except in compliance with the License.
6 You may obtain a copy of the License at
8 http://www.apache.org/licenses/LICENSE-2.0
10 Unless required by applicable law or agreed to in writing, software
11 distributed under the License is distributed on an "AS IS" BASIS,
12 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 See the License for the specific language governing permissions and
14 limitations under the License.
31 #include <vtkInformation.h>
32 #include <vtkInformationVector.h>
33 #include <vtkDemandDrivenPipeline.h>
34 #include <vtkObjectFactory.h>
35 #include <vtkPolyDataAlgorithm.h>
36 #include <vtkCellData.h>
37 #include <vtkPointData.h>
39 #include <vtkIdList.h>
40 #include <vtkAppendPolyData.h>
41 #include <vtkKdTreePointLocator.h>
42 #include <vtkCleanPolyData.h>
43 #include <vtkPolyDataConnectivityFilter.h>
44 #include <vtkSmartPointer.h>
45 #include <vtkModifiedBSPTree.h>
46 #include <vtkCellArrayIterator.h>
47 #include <vtkKdTree.h>
49 #include "vtkPolyDataBooleanFilter.h"
50 #include "vtkPolyDataContactFilter.h"
52 #include "Utilities.h"
54 vtkStandardNewMacro(vtkPolyDataBooleanFilter);
56 vtkPolyDataBooleanFilter::vtkPolyDataBooleanFilter () {
58 SetNumberOfInputPorts(2);
59 SetNumberOfOutputPorts(3);
64 contLines = vtkSmartPointer<vtkPolyData>::New();
66 modPdA = vtkSmartPointer<vtkPolyData>::New();
67 modPdB = vtkSmartPointer<vtkPolyData>::New();
69 cellDataA = vtkSmartPointer<vtkCellData>::New();
70 cellDataB = vtkSmartPointer<vtkCellData>::New();
72 cellIdsA = vtkSmartPointer<vtkIdTypeArray>::New();
73 cellIdsB = vtkSmartPointer<vtkIdTypeArray>::New();
75 OperMode = OPER_UNION;
79 vtkPolyDataBooleanFilter::~vtkPolyDataBooleanFilter () {
83 int vtkPolyDataBooleanFilter::ProcessRequest(vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector) {
85 if (request->Has(vtkDemandDrivenPipeline::REQUEST_DATA())) {
87 vtkInformation *inInfoA = inputVector[0]->GetInformationObject(0);
88 vtkInformation *inInfoB = inputVector[1]->GetInformationObject(0);
90 vtkPolyData *pdA = vtkPolyData::SafeDownCast(inInfoA->Get(vtkDataObject::DATA_OBJECT()));
91 vtkPolyData *pdB = vtkPolyData::SafeDownCast(inInfoB->Get(vtkDataObject::DATA_OBJECT()));
93 vtkInformation *outInfoA = outputVector->GetInformationObject(0);
94 vtkInformation *outInfoB = outputVector->GetInformationObject(1);
95 vtkInformation *outInfoC = outputVector->GetInformationObject(2);
97 resultA = vtkPolyData::SafeDownCast(outInfoA->Get(vtkDataObject::DATA_OBJECT()));
98 resultB = vtkPolyData::SafeDownCast(outInfoB->Get(vtkDataObject::DATA_OBJECT()));
99 resultC = vtkPolyData::SafeDownCast(outInfoC->Get(vtkDataObject::DATA_OBJECT()));
101 using clock = std::chrono::steady_clock;
102 std::vector<clock::duration> times;
103 clock::time_point start;
105 if (pdA->GetMTime() > timePdA || pdB->GetMTime() > timePdB) {
106 // eventuell vorhandene regionen vereinen
108 vtkSmartPointer<vtkCleanPolyData> cleanA = vtkSmartPointer<vtkCleanPolyData>::New();
109 cleanA->SetOutputPointsPrecision(DOUBLE_PRECISION);
110 cleanA->SetTolerance(1e-6);
111 cleanA->SetInputData(pdA);
114 vtkSmartPointer<vtkCleanPolyData> cleanB = vtkSmartPointer<vtkCleanPolyData>::New();
115 cleanB->SetOutputPointsPrecision(DOUBLE_PRECISION);
116 cleanB->SetTolerance(1e-6);
117 cleanB->SetInputData(pdB);
121 std::cout << "Exporting modPdA.vtk" << std::endl;
122 WriteVTK("modPdA.vtk", cleanA->GetOutput());
124 std::cout << "Exporting modPdB.vtk" << std::endl;
125 WriteVTK("modPdB.vtk", cleanB->GetOutput());
130 cellDataA->DeepCopy(cleanA->GetOutput()->GetCellData());
131 cellDataB->DeepCopy(cleanB->GetOutput()->GetCellData());
133 // ermittelt kontaktstellen
135 start = clock::now();
137 vtkSmartPointer<vtkPolyDataContactFilter> cl = vtkSmartPointer<vtkPolyDataContactFilter>::New();
138 cl->SetInputConnection(0, cleanA->GetOutputPort());
139 cl->SetInputConnection(1, cleanB->GetOutputPort());
142 times.push_back(clock::now()-start);
144 contLines->DeepCopy(cl->GetOutput());
146 modPdA->DeepCopy(cl->GetOutput(1));
147 modPdB->DeepCopy(cl->GetOutput(2));
150 std::cout << "Exporting contLines.vtk" << std::endl;
151 WriteVTK("contLines.vtk", contLines);
153 std::cout << "Exporting modPdA_1.vtk" << std::endl;
154 WriteVTK("modPdA_1.vtk", modPdA);
156 std::cout << "Exporting modPdB_1.vtk" << std::endl;
157 WriteVTK("modPdB_1.vtk", modPdB);
160 if (contLines->GetNumberOfCells() == 0) {
161 vtkErrorMacro("Inputs have no contact.");
166 // in den CellDatas steht drin, welche polygone einander schneiden
168 contsA = vtkIdTypeArray::SafeDownCast(contLines->GetCellData()->GetScalars("cA"));
169 contsB = vtkIdTypeArray::SafeDownCast(contLines->GetCellData()->GetScalars("cB"));
171 vtkIdTypeArray *accuracy = vtkIdTypeArray::SafeDownCast(contLines->GetPointData()->GetScalars("accuracy"));
173 vtkIdTypeArray *sourcesA = vtkIdTypeArray::SafeDownCast(contLines->GetCellData()->GetScalars("sourcesA"));
174 vtkIdTypeArray *sourcesB = vtkIdTypeArray::SafeDownCast(contLines->GetCellData()->GetScalars("sourcesB"));
176 vtkIdType i, numPts = contLines->GetNumberOfPoints();
178 vtkSmartPointer<vtkIdList> cells = vtkSmartPointer<vtkIdList>::New();
180 for (i = 0; i < numPts; i++) {
181 contLines->GetPointCells(i, cells);
183 if (cells->GetNumberOfIds() == 1) {
184 vtkErrorMacro("Contact ends suddenly.");
188 if (accuracy->GetValue(i) == 1) {
189 vtkErrorMacro("Contact goes through a cell with a bad shape.");
194 // sichert die OrigCellIds
196 vtkIdTypeArray *origCellIdsA = vtkIdTypeArray::SafeDownCast(modPdA->GetCellData()->GetScalars("OrigCellIds"));
197 vtkIdTypeArray *origCellIdsB = vtkIdTypeArray::SafeDownCast(modPdB->GetCellData()->GetScalars("OrigCellIds"));
199 cellIdsA->DeepCopy(origCellIdsA);
200 cellIdsB->DeepCopy(origCellIdsB);
202 vtkIdType numCellsA = modPdA->GetNumberOfCells(),
203 numCellsB = modPdB->GetNumberOfCells();
205 for (i = 0; i < numCellsA; i++) {
206 origCellIdsA->SetValue(i, i);
209 for (i = 0; i < numCellsB; i++) {
210 origCellIdsB->SetValue(i, i);
213 start = clock::now();
215 if (GetPolyStrips(modPdA, contsA, sourcesA, polyStripsA) ||
216 GetPolyStrips(modPdB, contsB, sourcesB, polyStripsB)) {
218 vtkErrorMacro("Strips are invalid.");
224 times.push_back(clock::now()-start);
226 // trennt die polygone an den linien
228 start = clock::now();
230 if (CutCells(modPdA, polyStripsA) ||
231 CutCells(modPdB, polyStripsB)) {
233 vtkErrorMacro("CutCells failed.");
238 times.push_back(clock::now()-start);
241 std::cout << "Exporting modPdA_2.vtk" << std::endl;
242 WriteVTK("modPdA_2.vtk", modPdA);
244 std::cout << "Exporting modPdB_2.vtk" << std::endl;
245 WriteVTK("modPdB_2.vtk", modPdB);
248 start = clock::now();
250 RestoreOrigPoints(modPdA, polyStripsA);
251 RestoreOrigPoints(modPdB, polyStripsB);
253 times.push_back(clock::now()-start);
256 std::cout << "Exporting modPdA_3.vtk" << std::endl;
257 WriteVTK("modPdA_3.vtk", modPdA);
259 std::cout << "Exporting modPdB_3.vtk" << std::endl;
260 WriteVTK("modPdB_3.vtk", modPdB);
263 start = clock::now();
265 ResolveOverlaps(modPdA, contsA, polyStripsA);
266 ResolveOverlaps(modPdB, contsB, polyStripsB);
268 times.push_back(clock::now()-start);
271 std::cout << "Exporting modPdA_4.vtk" << std::endl;
272 WriteVTK("modPdA_4.vtk", modPdA);
274 std::cout << "Exporting modPdB_4.vtk" << std::endl;
275 WriteVTK("modPdB_4.vtk", modPdB);
278 start = clock::now();
280 AddAdjacentPoints(modPdA, contsA, polyStripsA);
281 AddAdjacentPoints(modPdB, contsB, polyStripsB);
283 times.push_back(clock::now()-start);
286 std::cout << "Exporting modPdA_5.vtk" << std::endl;
287 WriteVTK("modPdA_5.vtk", modPdA);
289 std::cout << "Exporting modPdB_5.vtk" << std::endl;
290 WriteVTK("modPdB_5.vtk", modPdB);
293 start = clock::now();
295 DisjoinPolys(modPdA, polyStripsA);
296 DisjoinPolys(modPdB, polyStripsB);
298 times.push_back(clock::now()-start);
301 std::cout << "Exporting modPdA_6.vtk" << std::endl;
302 WriteVTK("modPdA_6.vtk", modPdA);
304 std::cout << "Exporting modPdB_6.vtk" << std::endl;
305 WriteVTK("modPdB_6.vtk", modPdB);
308 start = clock::now();
310 MergePoints(modPdA, polyStripsA);
311 MergePoints(modPdB, polyStripsB);
313 times.push_back(clock::now()-start);
316 std::cout << "Exporting modPdA_7.vtk" << std::endl;
317 WriteVTK("modPdA_7.vtk", modPdA);
319 std::cout << "Exporting modPdB_7.vtk" << std::endl;
320 WriteVTK("modPdB_7.vtk", modPdB);
323 timePdA = pdA->GetMTime();
324 timePdB = pdB->GetMTime();
328 start = clock::now();
330 if (CombineRegions()) {
331 vtkErrorMacro("Boolean operation failed.");
336 times.push_back(clock::now()-start);
339 double sum = std::chrono::duration_cast<std::chrono::duration<double>>(std::accumulate(times.begin(), times.end(), clock::duration())).count();
341 std::vector<clock::duration>::const_iterator itr;
342 for (itr = times.begin(); itr != times.end(); itr++) {
343 double time = std::chrono::duration_cast<std::chrono::duration<double>>(*itr).count();
345 std::cout << "Time " << (itr-times.begin())
346 << ": " << time << "s (" << (time/sum*100) << "%)"
357 void vtkPolyDataBooleanFilter::GetStripPoints (vtkPolyData *pd, vtkIdTypeArray *sources, PStrips &pStrips, IdsType &lines) {
360 std::cout << "GetStripPoints()" << std::endl;
363 StripPtsType &pts = pStrips.pts;
364 const IdsType &poly = pStrips.poly;
366 double a[3], b[3], sA[3], sB[3], u[3], v[3], w[3], n, t, d;
368 std::map<vtkIdType, vtkIdType> allPts;
370 vtkSmartPointer<vtkIdList> line = vtkSmartPointer<vtkIdList>::New();
372 for (auto lineId : lines) {
373 contLines->GetCellPoints(lineId, line);
375 // std::cout << "? " << contsA->GetValue(lineId)
376 // << ", " << contsB->GetValue(lineId)
377 // << ", [" << line->GetId(0)
378 // << ", " << line->GetId(1)
382 allPts.emplace(line->GetId(0), sources->GetTypedComponent(lineId, 0));
383 allPts.emplace(line->GetId(1), sources->GetTypedComponent(lineId, 1));
386 decltype(allPts)::const_iterator itr;
388 for (itr = allPts.begin(); itr != allPts.end(); ++itr) {
394 contLines->GetPoint(sp.ind, pt);
398 IdsType::const_iterator itrA, itrB;
400 for (itrA = poly.begin(); itrA != poly.end(); ++itrA) {
403 if (itrB == poly.end()) {
407 if (itr->second != NOTSET && *itrA != itr->second) {
411 pd->GetPoint(*itrA, a);
412 pd->GetPoint(*itrB, b);
414 vtkMath::Subtract(a, pt, sA);
415 vtkMath::Subtract(b, pt, sB);
417 // richtungsvektor und länge der kante
419 vtkMath::Subtract(b, a, u);
420 n = vtkMath::Norm(u);
424 vtkMath::Subtract(pt, a, v);
425 t = vtkMath::Dot(v, u)/(n*n);
427 vtkMath::Cross(v, u, w);
428 d = vtkMath::Norm(w)/n;
430 if (d < 1e-5 && t > -1e-5 && t < 1+1e-5) {
434 sp.t = std::min(1., std::max(0., t));
436 if (vtkMath::Norm(sA) < 1e-5) {
437 Cpy(sp.captPt, a, 3);
440 } else if (vtkMath::Norm(sB) < 1e-5) {
441 Cpy(sp.captPt, b, 3);
445 // u ist nicht normiert
446 vtkMath::MultiplyScalar(u, t);
449 vtkMath::Add(a, u, x);
452 Cpy(sp.captPt, x, 3);
454 sp.capt = Capt::EDGE;
463 // << ", " << sp.capt
467 if (itr->second != NOTSET && sp.edge[0] == NOTSET) {
471 pts.emplace(sp.ind, std::move(sp));
475 StripPtsType::iterator itr2;
477 for (itr2 = pts.begin(); itr2 != pts.end(); ++itr2) {
478 StripPt &sp = itr2->second;
480 if (sp.capt != Capt::NOT) {
481 if (sp.capt == Capt::B) {
484 sp.edge[0] = sp.edge[1];
486 auto itrA = std::find(poly.begin(), poly.end(), sp.edge[0]),
489 if (itrB == poly.end()) {
499 // fĂ¼r den schnitt werden die eingerasteten koordinaten verwendet
501 Cpy(sp.cutPt, sp.captPt, 3);
504 Cpy(sp.cutPt, sp.pt, 3);
507 sp.history.emplace_back(sp.edge[0], sp.edge[1]);
512 std::cout << "pts: " << std::endl;
513 for (itr2 = pts.begin(); itr2 != pts.end(); ++itr2) {
514 std::cout << itr2->first << ": " << itr2->second << std::endl;
520 bool vtkPolyDataBooleanFilter::GetPolyStrips (vtkPolyData *pd, vtkIdTypeArray *conts, vtkIdTypeArray *sources, PolyStripsType &polyStrips) {
522 std::cout << "GetPolyStrips()" << std::endl;
527 std::map<vtkIdType, IdsType> polyLines;
529 for (vtkIdType i = 0; i < conts->GetNumberOfTuples(); i++) {
530 vtkIdType poly = conts->GetValue(i);
532 /*if (poly != 1641) {
536 polyLines[poly].push_back(i);
539 std::vector<std::reference_wrapper<StripPt>> notCatched;
541 std::map<vtkIdType, IdsType>::iterator itr;
543 for (itr = polyLines.begin(); itr != polyLines.end(); ++itr) {
545 IdsType &lines = itr->second;
546 RemoveDuplicates(lines);
548 polyStrips.emplace(std::piecewise_construct,
549 std::forward_as_tuple(itr->first),
550 std::forward_as_tuple(pd, itr->first));
552 PStrips &pStrips = polyStrips.at(itr->first);
554 GetStripPoints(pd, sources, pStrips, lines);
556 for (auto &sp : pStrips.pts) {
557 sp.second.polyId = itr->first;
559 if (!sp.second.catched) {
560 notCatched.push_back(sp.second);
566 auto Next = [](const IdsType &ids, vtkIdType id) -> vtkIdType {
567 IdsType::const_iterator itr;
569 itr = std::find(ids.begin(), ids.end(), id);
571 if (++itr == ids.end()) {
578 for (StripPt &sp : notCatched) {
579 for (itr = polyLines.begin(); itr != polyLines.end(); ++itr) {
580 const PStrips &pStrips = polyStrips.at(itr->first);
583 const StripPt &corr = pStrips.pts.at(sp.ind);
586 if (corr.capt == Capt::A) {
588 sp.edge[0] = corr.edge[0];
589 sp.edge[1] = Next(polyStrips.at(sp.polyId).poly, sp.edge[0]);
593 Cpy(sp.captPt, corr.captPt, 3);
594 Cpy(sp.cutPt, sp.captPt, 3);
596 sp.history.emplace_back(sp.edge[0], sp.edge[1]);
607 std::cout << sp << std::endl;
614 vtkSmartPointer<vtkIdList> cells = vtkSmartPointer<vtkIdList>::New(),
615 line = vtkSmartPointer<vtkIdList>::New();
617 vtkIdType i, numCells;
619 StripPtsType::const_iterator itr2;
621 for (itr = polyLines.begin(); itr != polyLines.end(); ++itr) {
622 PStrips &pStrips = polyStrips.at(itr->first);
624 const IdsType &lines = itr->second;
625 const StripPtsType &pts = pStrips.pts;
627 for (itr2 = pts.begin(); itr2 != pts.end(); ++itr2) {
628 // ist der punkt auf einer kante, dĂ¼rfen von ihm mehr als 2 linien ausgehen
630 const StripPt &pt = itr2->second;
632 if (pt.capt == Capt::NOT) {
633 contLines->GetPointCells(pt.ind, cells);
635 numCells = cells->GetNumberOfIds();
637 std::set<vtkIdType> ends;
639 for (i = 0; i < numCells; i++) {
640 contLines->GetCellPoints(cells->GetId(i), line);
642 ends.insert(pt.ind == line->GetId(0) ? line->GetId(1) : line->GetId(0));
645 if (ends.size() > 2) {
653 std::deque<Pair> _lines;
655 vtkIdList *linePts = vtkIdList::New();
657 for (auto &i : lines) {
658 contLines->GetCellPoints(i, linePts);
659 _lines.emplace_back(linePts->GetId(0), linePts->GetId(1));
664 decltype(_lines)::iterator _itr;
666 auto FindRight = [&pts, &_lines, &_itr](StripType &strip, const std::size_t &id) -> bool {
667 auto &right = strip.back();
669 if (pts.at(right.ind).capt == Capt::NOT) {
670 for (_itr = _lines.begin(); _itr != _lines.end(); ++_itr) {
671 if (_itr->f == right.ind) {
672 strip.emplace_back(_itr->g, id);
676 } else if (_itr->g == right.ind) {
677 strip.emplace_back(_itr->f, id);
688 auto FindLeft = [&pts, &_lines, &_itr](StripType &strip, const std::size_t &id) -> bool {
689 auto &left = strip.front();
691 if (pts.at(left.ind).capt == Capt::NOT) {
692 for (_itr = _lines.begin(); _itr != _lines.end(); ++_itr) {
693 if (_itr->f == left.ind) {
694 strip.emplace_front(_itr->g, id);
698 } else if (_itr->g == left.ind) {
699 strip.emplace_front(_itr->f, id);
710 std::size_t stripId {0};
712 while (!_lines.empty()) {
713 auto &last = _lines.back();
715 StripType strip {{last.f, stripId}, {last.g, stripId}};
718 while (FindRight(strip, stripId)) {}
719 while (FindLeft(strip, stripId)) {}
721 pStrips.strips.push_back(std::move(strip));
727 CompleteStrips(pStrips);
731 // sucht nach schnitten zw. den strips
735 PolyStripsType::const_iterator itr;
736 StripType::const_iterator itr2;
738 for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) {
739 const PStrips &pStrips = itr->second;
741 const StripsType &strips = pStrips.strips;
742 const StripPtsType &pts = pStrips.pts;
743 const Base &base = pStrips.base;
745 vtkSmartPointer<vtkPoints> treePts = vtkSmartPointer<vtkPoints>::New();
747 vtkSmartPointer<vtkPolyData> treePd = vtkSmartPointer<vtkPolyData>::New();
750 std::map<vtkIdType, vtkIdType> ptIds;
754 for (const auto &p : pts) {
755 Transform(p.second.pt, pt, base);
757 ptIds.emplace(p.first, treePts->InsertNextPoint(pt[0], pt[1], 0));
760 for (const StripType &strip : strips) {
761 for (itr2 = strip.begin(); itr2 != strip.end()-1; ++itr2) {
763 vtkIdList *line = vtkIdList::New();
764 line->InsertNextId(ptIds[itr2->ind]);
765 line->InsertNextId(ptIds[(itr2+1)->ind]);
767 treePd->InsertNextCell(VTK_LINE, line);
773 treePd->SetPoints(treePts);
775 vtkSmartPointer<vtkModifiedBSPTree> tree = vtkSmartPointer<vtkModifiedBSPTree>::New();
776 tree->SetDataSet(treePd);
777 tree->BuildLocator();
779 vtkIdType numA, numB;
780 const vtkIdType *lineA, *lineB;
782 auto lineItr = vtk::TakeSmartPointer(treePd->GetLines()->NewIterator());
784 for (lineItr->GoToFirstCell(); !lineItr->IsDoneWithTraversal(); lineItr->GoToNextCell()) {
785 lineItr->GetCurrentCell(numA, lineA);
787 double ptA[3], ptB[3];
789 treePts->GetPoint(lineA[0], ptA);
790 treePts->GetPoint(lineA[1], ptB);
792 vtkSmartPointer<vtkIdList> lineIds = vtkSmartPointer<vtkIdList>::New();
794 tree->IntersectWithLine(ptA, ptB, 1e-5, nullptr, lineIds);
796 for (vtkIdType i = 0; i < lineIds->GetNumberOfIds(); i++) {
797 treePd->GetCellPoints(lineIds->GetId(i), numB, lineB);
799 if (lineB[0] != lineA[0] && lineB[1] != lineA[0] && lineB[0] != lineA[1] && lineB[1] != lineA[1]) {
815 void vtkPolyDataBooleanFilter::RemoveDuplicates (IdsType &lines) {
819 // die indexe der enden auf Ă¼bereinstimmung prĂ¼fen
821 vtkIdList *linePtsA = vtkIdList::New();
822 vtkIdList *linePtsB = vtkIdList::New();
824 std::size_t i, j, numLines = lines.size();
826 for (i = 0; i < numLines-1; i++) {
829 contLines->GetCellPoints(lines[i], linePtsA);
831 while (j < lines.size()) {
832 contLines->GetCellPoints(lines[j], linePtsB);
834 if ((linePtsA->GetId(0) == linePtsB->GetId(0) && linePtsA->GetId(1) == linePtsB->GetId(1)) ||
835 (linePtsA->GetId(0) == linePtsB->GetId(1) && linePtsA->GetId(1) == linePtsB->GetId(0))) {
844 // keine vorzeitige unterbrechung der while-schleife
845 unique.push_back(lines[i]);
849 unique.push_back(lines.back());
858 void vtkPolyDataBooleanFilter::CompleteStrips (PStrips &pStrips) {
859 StripsType::iterator itr;
861 for (itr = pStrips.strips.begin(); itr != pStrips.strips.end(); ++itr) {
862 const StripPt &start = pStrips.pts[itr->front().ind],
863 &end = pStrips.pts[itr->back().ind];
865 if (start.ind != end.ind) {
866 if (start.capt == Capt::NOT) {
867 StripType s(itr->rbegin(), itr->rend()-1);
868 itr->insert(itr->begin(), s.begin(), s.end());
870 } else if (end.capt == Capt::NOT) {
871 StripType s(itr->rbegin()+1, itr->rend());
872 itr->insert(itr->end(), s.begin(), s.end());
879 bool vtkPolyDataBooleanFilter::HasArea (const StripType &strip) const {
882 std::size_t i, n = strip.size();
885 for (i = 0; i < (n-1)/2; i++) {
886 area = strip[i].ind != strip[n-i-1].ind;
893 void ComputeNormal (const StripPtsType &pts, const RefsType &poly, double *n) {
894 n[0] = 0; n[1] = 0; n[2] = 0;
896 RefsType::const_iterator itrA, itrB;
898 for (itrA = poly.begin(); itrA != poly.end(); ++itrA) {
901 if (itrB == poly.end()) {
905 const StripPtR &spA = *itrA,
908 auto pA = pts.find(spA.ind);
909 auto pB = pts.find(spB.ind);
911 const double *ptA = pA->second.cutPt,
912 *ptB = pB->second.cutPt;
914 n[0] += (ptA[1]-ptB[1])*(ptA[2]+ptB[2]);
915 n[1] += (ptA[2]-ptB[2])*(ptA[0]+ptB[0]);
916 n[2] += (ptA[0]-ptB[0])*(ptA[1]+ptB[1]);
919 vtkMath::Normalize(n);
922 bool vtkPolyDataBooleanFilter::CutCells (vtkPolyData *pd, PolyStripsType &polyStrips) {
924 std::cout << "CutCells()" << std::endl;
927 vtkPoints *pdPts = pd->GetPoints();
929 vtkIdTypeArray *origCellIds = vtkIdTypeArray::SafeDownCast(pd->GetCellData()->GetScalars("OrigCellIds"));
931 PolyStripsType::iterator itr;
933 for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) {
935 vtkIdType polyInd = itr->first;
936 PStrips &pStrips = itr->second;
938 /*if (polyInd != 1641) {
942 StripsType &strips = pStrips.strips;
943 StripPtsType &pts = pStrips.pts;
945 IdsType &poly = pStrips.poly;
947 vtkIdType origId = origCellIds->GetValue(polyInd);
950 std::cout << "polyInd " << polyInd << ", poly [";
951 for (auto &p : poly) {
952 std::cout << p << ", ";
954 std::cout << "]" << std::endl;
957 std::size_t numPts = poly.size();
959 std::map<vtkIdType, RefsType> edges;
961 StripsType::iterator itr2;
962 StripType::iterator itr3;
967 auto fct = [&](const StripType &s) {
968 return pts[s.front().ind].capt == Capt::NOT && pts[s.back().ind].capt == Capt::NOT;
971 std::copy_if(strips.begin(), strips.end(), std::back_inserter(holes), fct);
973 strips.erase(std::remove_if(strips.begin(), strips.end(), fct), strips.end());
975 for (itr2 = strips.begin(); itr2 != strips.end(); ++itr2) {
976 StripType &strip = *itr2;
979 std::cout << (itr2-strips.begin()) << ". strip [";
980 for (auto &s : strip) {
981 std::cout << s.ind << ", ";
983 std::cout << "]" << std::endl;
987 if (pts[strip.front().ind].edge[0] == pts[strip.back().ind].edge[0]
988 && strip.front().ind != strip.back().ind
989 && pts[strip.front().ind].t > pts[strip.back().ind].t) {
991 std::reverse(strip.begin(), strip.end());
994 StripPt &start = pts[strip.front().ind],
995 &end = pts[strip.back().ind];
997 strip.front().side = Side::START;
998 strip.back().side = Side::END;
1000 strip.front().ref = start.edge[0];
1001 strip.back().ref = end.edge[0];
1003 // nachfolgend könnte man dann anfang und ende weglassen
1005 for (itr3 = strip.begin(); itr3 != strip.end(); ++itr3) {
1006 StripPt &sp = pts[itr3->ind];
1008 itr3->desc[0] = pdPts->InsertNextPoint(sp.cutPt);
1009 itr3->desc[1] = pdPts->InsertNextPoint(sp.cutPt);
1012 std::cout << sp << " => " << *itr3 << std::endl;
1019 edges[start.edge[0]].push_back(std::ref(strip.front()));
1020 edges[end.edge[0]].push_back(std::ref(strip.back()));
1023 // sortiert die punkte auf den kanten
1025 std::map<vtkIdType, RefsType>::iterator itr4;
1027 IdsType::iterator itr5;
1028 StripType::reverse_iterator itr6;
1030 RefsType::iterator itr7;
1032 for (itr4 = edges.begin(); itr4 != edges.end(); ++itr4) {
1033 RefsType &edge = itr4->second;
1036 std::cout << "edge (" << itr4->first << ", _)" << std::endl;
1039 std::sort(edge.begin(), edge.end(), [&](const StripPtR &a, const StripPtR &b) {
1040 const StripPt &a_ = pts[a.ind],
1044 // std::cout << "a_: " << a_ << " -> strip " << a.strip << std::endl;
1045 // std::cout << "b_: " << b_ << " -> strip " << b.strip << std::endl;
1048 if (a_.ind == b_.ind) {
1049 // strips beginnen im gleichen punkt
1051 if (a.strip != b.strip) {
1052 // gehören nicht dem gleichen strip an
1054 StripType &stripA = strips[a.strip],
1055 &stripB = strips[b.strip];
1057 // andere enden ermitteln
1059 const StripPtR &eA = (&a == &(stripA.front())) ? stripA.back() : stripA.front(),
1060 &eB = (&b == &(stripB.front())) ? stripB.back() : stripB.front();
1062 const StripPt &eA_ = pts[eA.ind],
1066 // std::cout << "eA_: " << eA_ << std::endl;
1067 // std::cout << "eB_: " << eA_ << std::endl;
1070 if (eA_.ind != eB_.ind) {
1071 auto i = std::find(poly.begin(), poly.end(), itr4->first)-poly.begin(),
1072 iA = std::find(poly.begin(), poly.end(), eA_.edge[0])-poly.begin(),
1073 iB = std::find(poly.begin(), poly.end(), eB_.edge[0])-poly.begin();
1075 double dA = static_cast<double>(Mod(iA-i, numPts))+eA_.t,
1076 dB = static_cast<double>(Mod(iB-i, numPts))+eB_.t;
1078 if (i == iA && a_.t > eA_.t) {
1079 dA += static_cast<double>(numPts);
1082 if (i == iB && b_.t > eB_.t) {
1083 dB += static_cast<double>(numPts);
1087 // std::cout << "dA=" << dA << ", dB=" << dB << std::endl;
1094 if (a.side == Side::START) {
1095 poly_.insert(poly_.end(), stripA.begin(), stripA.end());
1097 poly_.insert(poly_.end(), stripA.rbegin(), stripA.rend());
1100 if (b.side == Side::START) {
1101 poly_.insert(poly_.end(), stripB.rbegin()+1, stripB.rend()-1);
1103 poly_.insert(poly_.end(), stripB.begin()+1, stripB.end()-1);
1107 ComputeNormal(pts, poly_, n);
1109 double ang = vtkMath::Dot(pStrips.n, n);
1112 // std::cout << "ang=" << ang*180/PI << std::endl;
1115 return ang < .999999;
1121 StripType &strip = strips[a.strip];
1123 if (HasArea(strip)) {
1124 RefsType poly_(strip.begin(), strip.end()-1);
1127 ComputeNormal(pts, poly_, n);
1129 double ang = vtkMath::Dot(pStrips.n, n);
1131 if (ang > .999999) {
1132 std::reverse(strip.begin(), strip.end());
1139 // reihenfolge von a und b bereits richtig
1150 for (auto& p : edge) {
1151 std::cout << p << std::endl;
1157 // baut die strips ein
1159 std::deque<IdsType> polys;
1160 polys.push_back(pStrips.poly);
1162 for (itr2 = strips.begin(); itr2 != strips.end(); ++itr2) {
1163 StripType &strip = *itr2;
1165 StripPtR &start = strip.front(),
1166 &end = strip.back();
1169 std::cout << "strip " << start.strip
1170 << " , refs (" << start.ref << ", " << end.ref << ")"
1174 std::size_t cycle = 0;
1177 if (cycle == polys.size()) {
1181 IdsType next = polys.front();
1184 std::vector<IdsType> newPolys(2);
1186 if (std::find(next.begin(), next.end(), start.ref) != next.end()) {
1187 if (start.ref == end.ref) {
1188 for (itr5 = next.begin(); itr5 != next.end(); ++itr5) {
1189 newPolys[0].push_back(*itr5);
1192 std::cout << "adding " << *itr5 << " to 0" << std::endl;
1195 if (*itr5 == start.ref) {
1196 for (itr3 = strip.begin(); itr3 != strip.end(); ++itr3) {
1197 newPolys[0].push_back(itr3->desc[0]);
1200 std::cout << "adding " << itr3->desc[0] << " to 0" << std::endl;
1207 // strip selbst ist ein polygon
1209 for (itr6 = strip.rbegin(); itr6 != strip.rend(); ++itr6) {
1210 newPolys[1].push_back(itr6->desc[1]);
1213 std::cout << "adding " << itr6->desc[1] << " to 1" << std::endl;
1219 std::size_t curr = 0;
1221 for (itr5 = next.begin(); itr5 != next.end(); ++itr5) {
1222 IdsType &poly = newPolys[curr];
1224 poly.push_back(*itr5);
1227 std::cout << "adding " << *itr5 << " to " << curr << std::endl;
1230 if (*itr5 == start.ref) {
1231 for (itr3 = strip.begin(); itr3 != strip.end(); ++itr3) {
1232 poly.push_back(itr3->desc[0]);
1235 std::cout << "adding " << itr3->desc[0] << " to " << curr << std::endl;
1240 curr = curr == 0 ? 1 : 0;
1242 } else if (*itr5 == end.ref) {
1243 for (itr6 = strip.rbegin(); itr6 != strip.rend(); ++itr6) {
1244 poly.push_back(itr6->desc[1]);
1247 std::cout << "adding " << itr6->desc[1] << " to " << curr << std::endl;
1252 curr = curr == 0 ? 1 : 0;
1259 if (!newPolys[1].empty()) {
1260 // refs aktualisieren
1262 for (itr4 = edges.begin(); itr4 != edges.end(); ++itr4) {
1263 RefsType &edge = itr4->second;
1265 RefsType::iterator itrA;
1267 for (itrA = edge.begin()+1; itrA != edge.end(); ++itrA) {
1268 StripPtR &sp = *itrA;
1270 if (sp.strip > start.strip) {
1272 std::cout << "sp: ind " << sp.ind << ", strip " << sp.strip << std::endl;
1275 RefsType::const_reverse_iterator itrB(itrA);
1277 std::shared_ptr<StripPtR> _p;
1279 for (; itrB != edge.rend(); ++itrB) {
1280 const StripPtR &p = *itrB;
1282 if (p.strip != sp.strip) {
1283 if (p.strip <= start.strip) {
1285 std::cout << "ref " << sp.ref;
1288 if (p.side == Side::END) {
1295 std::cout << " -> " << sp.ref << " (from strip " << p.strip << ", ind " << p.ind << ")" << std::endl;
1298 _p = std::make_shared<StripPtR>(p);
1305 std::cout << "~1 ref " << sp.ref << " -> " << p.ref << " (from strip " << p.strip << ", ind " << p.ind << ")" << std::endl;
1313 RefsType::const_iterator itrC(itrA);
1317 for (; itrC != edge.end(); ++itrC) {
1318 const StripPtR &p = *itrC;
1320 if (p.ind != sp.ind) {
1324 if (p.strip <= start.strip) {
1325 if (_p && p.ind == _p->ind && p.strip < _p->strip) {
1330 std::cout << "~2 ref " << sp.ref;
1333 if (p.side == Side::START) {
1340 std::cout << " -> " << sp.ref << " (from strip " << p.strip << ", ind " << p.ind << ")" << std::endl;
1350 // erstellt die history
1352 auto _s = std::find_if(edge.begin(), edge.end(), [&start](const StripPtR &r) {
1353 return &start == &r;
1356 if (_s != edge.end()) {
1357 for (itr7 = edge.begin(); itr7 != edge.end(); ++itr7) {
1358 StripPtR &sp = *itr7;
1359 StripPt &_sp = pts[sp.ind];
1360 auto &history = _sp.history;
1362 if (&sp != &start) {
1363 if (_sp.t < pts[start.ind].t) {
1364 history.push_back({history.back().f, start.desc[0]});
1366 history.push_back({start.desc[1], history.back().g});
1374 auto _e = std::find_if(edge.begin(), edge.end(), [&end](const StripPtR &r) {
1378 if (_e != edge.end()) {
1379 for (itr7 = edge.begin(); itr7 != edge.end(); ++itr7) {
1380 StripPtR &sp = *itr7;
1381 StripPt &_sp = pts[sp.ind];
1382 auto &history = _sp.history;
1385 if (_sp.t < pts[end.ind].t) {
1386 history.push_back({history.back().f, end.desc[1]});
1388 history.push_back({end.desc[0], history.back().g});
1397 if (edge.size() > 1) {
1398 StripPtR &a = edge.front(),
1399 &b = *(edge.begin()+1);
1402 && b.strip == start.strip
1403 && pts[a.ind].capt == Capt::A) { // sollte weg
1406 std::cout << "~3 ref " << a.ref;
1409 if (b.side == Side::START) {
1416 std::cout << " -> " << a.ref << " (from strip " << b.strip << ", ind " << b.ind << ")" << std::endl;
1424 // doppelte punkte entfernen
1428 for (auto &newPoly : newPolys) {
1431 std::map<vtkIdType, Point3d> _pts;
1433 for (vtkIdType id : newPoly) {
1434 pd->GetPoint(id, pt);
1436 _pts.emplace(std::piecewise_construct,
1437 std::forward_as_tuple(id),
1438 std::forward_as_tuple(pt[0], pt[1], pt[2]));
1441 IdsType::const_iterator itrA, itrB;
1443 for (itrA = newPoly.begin(); itrA != newPoly.end(); ++itrA) {
1445 if (itrB == newPoly.end()) {
1446 itrB = newPoly.begin();
1449 auto _a = _pts.find(*itrA);
1450 auto _b = _pts.find(*itrB);
1452 if (_a->second == _b->second) {
1455 std::cout << "removing " << *itrA << std::endl;
1458 _newPoly.push_back(*itrA);
1462 newPoly.swap(_newPoly);
1466 // prĂ¼ft, ob die erstellten polys gĂ¼ltig sind
1468 if (newPolys[0].size() > 2) {
1469 polys.push_back(newPolys[0]);
1472 if (HasArea(strip) && newPolys[1].size() > 2) {
1473 polys.push_back(newPolys[1]);
1481 polys.push_back(next);
1490 // erzeugte polys hinzufĂ¼gen
1493 descIds.reserve(polys.size());
1495 for (auto &p : polys) {
1496 vtkIdList *cell = vtkIdList::New();
1498 for (vtkIdType &id : p) {
1499 cell->InsertNextId(id);
1502 descIds.emplace_back(pd->InsertNextCell(VTK_POLYGON, cell));
1503 origCellIds->InsertNextValue(origId);
1508 pd->DeleteCell(polyInd);
1511 if (!holes.empty()) {
1513 Merger(pd, pStrips, holes, descIds, origId).Run();
1514 } catch (const std::runtime_error &e) {
1521 pd->RemoveDeletedCells();
1528 void vtkPolyDataBooleanFilter::RestoreOrigPoints (vtkPolyData *pd, PolyStripsType &polyStrips) {
1531 std::cout << "RestoreOrigPoints()" << std::endl;
1536 vtkKdTreePointLocator *loc = vtkKdTreePointLocator::New();
1537 loc->SetDataSet(pd);
1538 loc->BuildLocator();
1541 bool operator() (const StripPt &l, const StripPt &r) const {
1542 return l.ind < r.ind;
1546 PolyStripsType::const_iterator itr;
1547 StripPtsType::const_iterator itr2;
1549 std::set<std::reference_wrapper<const StripPt>, Cmp> ends;
1551 for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) {
1552 const PStrips &pStrips = itr->second;
1554 for (itr2 = pStrips.pts.begin(); itr2 != pStrips.pts.end(); ++itr2) {
1555 const StripPt &sp = itr2->second;
1557 if (sp.capt != Capt::NOT) {
1564 vtkIdList *pts = vtkIdList::New();
1566 vtkIdType i, numPts;
1568 for (const StripPt &sp : ends) {
1569 FindPoints(loc, sp.cutPt, pts);
1570 numPts = pts->GetNumberOfIds();
1572 for (i = 0; i < numPts; i++) {
1573 pd->GetPoints()->SetPoint(pts->GetId(i), sp.pt);
1579 loc->FreeSearchStructure();
1584 void vtkPolyDataBooleanFilter::DisjoinPolys (vtkPolyData *pd, PolyStripsType &polyStrips) {
1587 std::cout << "DisjoinPolys()" << std::endl;
1592 vtkKdTreePointLocator *loc = vtkKdTreePointLocator::New();
1593 loc->SetDataSet(pd);
1594 loc->BuildLocator();
1597 bool operator() (const StripPt &l, const StripPt &r) const {
1598 return l.ind < r.ind;
1602 PolyStripsType::const_iterator itr;
1603 StripPtsType::const_iterator itr2;
1605 std::set<std::reference_wrapper<const StripPt>, Cmp> ends;
1607 for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) {
1608 const PStrips &pStrips = itr->second;
1610 for (itr2 = pStrips.pts.begin(); itr2 != pStrips.pts.end(); ++itr2) {
1611 const StripPt &sp = itr2->second;
1613 if (sp.capt == Capt::A) {
1619 vtkIdList *pts = vtkIdList::New();
1620 vtkIdList *cells = vtkIdList::New();
1622 vtkIdType i, j, numPts, numCells;
1624 for (const StripPt &sp : ends) {
1625 FindPoints(loc, sp.pt, pts);
1626 numPts = pts->GetNumberOfIds();
1628 for (i = 0; i < numPts; i++) {
1629 pd->GetPointCells(pts->GetId(i), cells);
1630 numCells = cells->GetNumberOfIds();
1633 for (j = 0; j < numCells; j++) {
1634 pd->ReplaceCellPoint(cells->GetId(j), pts->GetId(i), pd->GetPoints()->InsertNextPoint(sp.pt));
1643 loc->FreeSearchStructure();
1648 void vtkPolyDataBooleanFilter::ResolveOverlaps (vtkPolyData *pd, vtkIdTypeArray *conts, PolyStripsType &polyStrips) {
1651 std::cout << "ResolveOverlaps()" << std::endl;
1657 contLines->BuildLinks();
1659 vtkKdTreePointLocator *loc = vtkKdTreePointLocator::New();
1660 loc->SetDataSet(pd);
1661 loc->BuildLocator();
1664 bool operator() (const StripPt &l, const StripPt &r) const {
1665 return l.ind < r.ind;
1669 PolyStripsType::const_iterator itr;
1670 StripPtsType::const_iterator itr2;
1672 std::set<std::reference_wrapper<const StripPt>, Cmp> ends;
1674 for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) {
1675 const PStrips &pStrips = itr->second;
1677 for (itr2 = pStrips.pts.begin(); itr2 != pStrips.pts.end(); ++itr2) {
1678 const StripPt &sp = itr2->second;
1680 if (sp.capt == Capt::EDGE) {
1686 vtkIdList *links = vtkIdList::New(),
1687 *ptsA = vtkIdList::New(),
1688 *ptsB = vtkIdList::New(),
1689 *cells = vtkIdList::New(),
1690 *poly = vtkIdList::New();
1692 double ptA[3], ptB[3];
1694 vtkIdType i, j, k, numPtsA, numPtsB, numCells, numPts, idA, idB;
1696 std::map<Pair, IdsType> sharedPts;
1698 for (const StripPt &sp : ends) {
1699 contLines->GetPointCells(sp.ind, links);
1701 if (links->GetNumberOfIds() == 2
1702 && conts->GetValue(links->GetId(0)) != conts->GetValue(links->GetId(1))) {
1704 std::vector<Pair> history(sp.history.rbegin(), sp.history.rend());
1706 for (const Pair &p : history) {
1707 pd->GetPoint(p.f, ptA);
1708 pd->GetPoint(p.g, ptB);
1710 FindPoints(loc, ptA, ptsA);
1711 FindPoints(loc, ptB, ptsB);
1713 numPtsA = ptsA->GetNumberOfIds();
1714 numPtsB = ptsB->GetNumberOfIds();
1716 std::vector<Pair> cellsA, cellsB;
1718 for (i = 0; i < numPtsA; i++) {
1719 pd->GetPointCells(ptsA->GetId(i), cells);
1720 numCells = cells->GetNumberOfIds();
1722 for (j = 0; j < numCells; j++) {
1723 cellsA.emplace_back(ptsA->GetId(i), cells->GetId(j));
1727 for (i = 0; i < numPtsB; i++) {
1728 pd->GetPointCells(ptsB->GetId(i), cells);
1729 numCells = cells->GetNumberOfIds();
1731 for (j = 0; j < numCells; j++) {
1732 cellsB.emplace_back(ptsB->GetId(i), cells->GetId(j));
1736 for (auto &cellA : cellsA) {
1737 for (auto &cellB : cellsB) {
1738 if (cellA.g == cellB.g) {
1739 // kante/poly existiert noch
1741 pd->GetCellPoints(cellA.g, poly);
1742 numPts = poly->GetNumberOfIds();
1744 for (k = 0; k < numPts; k++) {
1745 idA = poly->GetId(k);
1746 idB = k+1 == numPts ? poly->GetId(0) : poly->GetId(k+1);
1749 && cellA.f == idB) {
1751 IdsType &pts = sharedPts[Pair{sp.ind, cellA.g}];
1753 pts.push_back(cellA.f);
1754 pts.push_back(cellB.f);
1772 decltype(sharedPts)::iterator itr3;
1774 for (itr3 = sharedPts.begin(); itr3 != sharedPts.end(); ++itr3) {
1775 const Pair &p = itr3->first;
1777 IdsType &pts = itr3->second;
1779 std::sort(pts.begin(), pts.end());
1781 auto found = std::adjacent_find(pts.begin(), pts.end());
1783 if (found != pts.end()) {
1784 contLines->GetPoint(p.f, pt);
1786 i = pd->GetPoints()->InsertNextPoint(pt);
1787 pd->ReplaceCellPoint(p.g, *found, i);
1797 loc->FreeSearchStructure();
1802 void vtkPolyDataBooleanFilter::AddAdjacentPoints (vtkPolyData *pd, vtkIdTypeArray *conts, PolyStripsType &polyStrips) {
1805 std::cout << "AddAdjacentPoints()" << std::endl;
1810 vtkIdTypeArray *origCellIds = vtkIdTypeArray::SafeDownCast(pd->GetCellData()->GetScalars("OrigCellIds"));
1813 bool operator() (const StripPt &l, const StripPt &r) const {
1818 vtkKdTreePointLocator *loc = vtkKdTreePointLocator::New();
1819 loc->SetDataSet(pd);
1820 loc->BuildLocator();
1822 vtkIdList *cells = vtkIdList::New(),
1823 *ptsA = vtkIdList::New(),
1824 *ptsB = vtkIdList::New(),
1825 *poly = vtkIdList::New(),
1826 *newPoly = vtkIdList::New();
1828 vtkIdType i, j, numCells, numPtsA, numPtsB, numPts, idA, idB;
1830 PolyStripsType::const_iterator itr;
1831 StripPtsType::const_iterator itr2;
1833 for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) {
1834 const PStrips &pStrips = itr->second;
1836 std::map<Pair, std::set<std::reference_wrapper<const StripPt>, Cmp>> edgePts;
1838 for (itr2 = pStrips.pts.begin(); itr2 != pStrips.pts.end(); ++itr2) {
1839 const StripPt &sp = itr2->second;
1841 if (sp.capt == Capt::EDGE) {
1842 edgePts[{sp.edge[0], sp.edge[1]}].emplace(sp);
1846 decltype(edgePts)::iterator itr3;
1848 for (itr3 = edgePts.begin(); itr3 != edgePts.end(); ++itr3) {
1849 const Pair &edge = itr3->first;
1850 auto pts = itr3->second;
1854 pd->GetPoint(edge.f, spA.pt);
1855 pd->GetPoint(edge.g, spB.pt);
1863 std::vector<decltype(pts)::value_type> _pts(pts.rbegin(), pts.rend());
1865 decltype(_pts)::const_iterator itrA, itrB, itrC;
1867 itrA = _pts.begin();
1869 while (itrA != _pts.end()-1) {
1872 while (itrB != _pts.end()-1) {
1873 contLines->GetPointCells(itrB->get().ind, cells);
1874 numCells = cells->GetNumberOfIds();
1876 // std::cout << "[";
1877 // for (i = 0; i < numCells; i++) {
1878 // std::cout << conts->GetValue(cells->GetId(i)) << ", ";
1880 // std::cout << "]" << std::endl;
1887 if (std::distance(itrA, itrB) > 1) {
1888 FindPoints(loc, itrA->get().pt, ptsA);
1889 FindPoints(loc, itrB->get().pt, ptsB);
1891 numPtsA = ptsA->GetNumberOfIds();
1892 numPtsB = ptsB->GetNumberOfIds();
1894 std::vector<Pair> polysA, polysB;
1896 for (i = 0; i < numPtsA; i++) {
1897 pd->GetPointCells(ptsA->GetId(i), cells);
1898 numCells = cells->GetNumberOfIds();
1900 for (j = 0; j < numCells; j++) {
1901 polysA.emplace_back(cells->GetId(j), ptsA->GetId(i));
1905 for (i = 0; i < numPtsB; i++) {
1906 pd->GetPointCells(ptsB->GetId(i), cells);
1907 numCells = cells->GetNumberOfIds();
1909 for (j = 0; j < numCells; j++) {
1910 polysB.emplace_back(cells->GetId(j), ptsB->GetId(i));
1914 /*for (const Pair &pA : polysA) {
1915 std::cout << "pA " << pA << std::endl;
1918 for (const Pair &pB : polysB) {
1919 std::cout << "pB " << pB << std::endl;
1922 for (const Pair &pA : polysA) {
1923 for (const Pair &pB : polysB) {
1924 if (pA.f == pB.f && pd->GetCellType(pA.f) != VTK_EMPTY_CELL) {
1926 pd->GetCellPoints(pA.f, poly);
1927 numPts = poly->GetNumberOfIds();
1931 for (j = 0; j < numPts; j++) {
1932 newPoly->InsertNextId(poly->GetId(j));
1934 idA = poly->GetId(j);
1935 idB = j+1 == numPts ? poly->GetId(0) : poly->GetId(j+1);
1940 for (itrC = itrA+1; itrC != itrB; ++itrC) {
1941 // newPoly->InsertNextId(pd->InsertNextLinkedPoint(itrC->get().pt, 1));
1943 pd->InsertNextLinkedPoint(1);
1945 newPoly->InsertNextId(pd->GetPoints()->InsertNextPoint(itrC->get().pt));
1950 pd->RemoveReferenceToCell(idA, pA.f);
1953 pd->DeleteCell(pA.f);
1955 pd->InsertNextLinkedCell(VTK_POLYGON, newPoly->GetNumberOfIds(), newPoly->GetPointer(0));
1957 origCellIds->InsertNextValue(origCellIds->GetValue(pA.f));
1976 loc->FreeSearchStructure();
1979 pd->RemoveDeletedCells();
1983 void vtkPolyDataBooleanFilter::MergePoints (vtkPolyData *pd, PolyStripsType &polyStrips) {
1986 std::cout << "MergePoints()" << std::endl;
1992 contLines->BuildLinks();
1994 vtkKdTreePointLocator *loc = vtkKdTreePointLocator::New();
1995 loc->SetDataSet(pd);
1996 loc->BuildLocator();
1998 PolyStripsType::const_iterator itr;
1999 StripsType::const_iterator itr2;
2001 std::map<vtkIdType, std::set<vtkIdType>> neighPts;
2003 vtkIdList *pts = vtkIdList::New();
2005 vtkIdType i, j, numPts;
2007 for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) {
2008 const PStrips &pStrips = itr->second;
2010 for (itr2 = pStrips.strips.begin(); itr2 != pStrips.strips.end(); ++itr2) {
2011 const StripType &strip = *itr2;
2013 const StripPtR &spA = strip.front(),
2014 &spB = strip.back();
2016 const auto beforeA = pStrips.pts.find((strip.begin()+1)->ind),
2017 beforeB = pStrips.pts.find((strip.end()-2)->ind);
2019 FindPoints(loc, beforeA->second.pt, pts);
2020 numPts = pts->GetNumberOfIds();
2022 for (i = 0; i < numPts; i++) {
2023 neighPts[spA.ind].insert(pts->GetId(i));
2026 FindPoints(loc, beforeB->second.pt, pts);
2027 numPts = pts->GetNumberOfIds();
2029 for (i = 0; i < numPts; i++) {
2030 neighPts[spB.ind].insert(pts->GetId(i));
2037 vtkIdList *polys = vtkIdList::New(),
2038 *poly = vtkIdList::New();
2040 vtkIdType ind, polyId, _numPts, before, after;
2042 decltype(neighPts)::const_iterator itr3;
2044 for (itr3 = neighPts.begin(); itr3 != neighPts.end(); ++itr3) {
2045 const auto &inds = itr3->second;
2047 std::map<Point3d, std::vector<Pair>> pairs;
2049 contLines->GetPoint(itr3->first, pt);
2051 FindPoints(loc, pt, pts);
2052 numPts = pts->GetNumberOfIds();
2054 for (i = 0; i < numPts; i++) {
2055 ind = pts->GetId(i);
2057 pd->GetPointCells(ind, polys);
2059 if (polys->GetNumberOfIds() > 0) {
2060 polyId = polys->GetId(0);
2062 pd->GetCellPoints(polyId, poly);
2063 _numPts = poly->GetNumberOfIds();
2065 for (j = 0; j < _numPts; j++) {
2066 if (poly->GetId(j) == ind) {
2071 // wieder davor und danach ermitteln
2073 before = poly->GetId(j == 0 ? _numPts-1 : j-1);
2074 after = poly->GetId(j+1 == _numPts ? 0 : j+1);
2076 if (std::find(inds.begin(), inds.end(), before) == inds.end()) {
2077 pd->GetPoint(before, pt);
2078 pairs[{pt[0], pt[1], pt[2]}].emplace_back(polyId, ind);
2081 if (std::find(inds.begin(), inds.end(), after) == inds.end()) {
2082 pd->GetPoint(after, pt);
2083 pairs[{pt[0], pt[1], pt[2]}].emplace_back(polyId, ind);
2089 std::deque<std::deque<std::reference_wrapper<const Pair>>> Pairs;
2091 decltype(pairs)::const_iterator itr4;
2093 for (itr4 = pairs.begin(); itr4 != pairs.end(); ++itr4) {
2094 const auto &p = itr4->second;
2096 if (p.size() == 2) {
2097 auto _pts = {std::ref(p.front()), std::ref(p.back())}; // std::initializer_list
2098 Pairs.emplace_back(_pts);
2102 decltype(Pairs)::iterator itr5;
2104 /*for (itr5 = Pairs.begin(); itr5 != Pairs.end(); ++itr5) {
2105 for (auto &p : *itr5) {
2106 std::cout << p.get() << ", ";
2108 std::cout << std::endl;
2111 decltype(Pairs)::value_type group;
2113 decltype(group)::const_iterator itr6;
2115 while (!Pairs.empty()) {
2116 if (group.empty()) {
2117 group = Pairs.front();
2121 itr5 = Pairs.begin();
2123 while (itr5 != Pairs.end()) {
2124 const auto &next = *itr5;
2126 if (next.front().get() == group.front().get()) {
2127 group.emplace_front(next.back());
2129 itr5 = Pairs.begin();
2130 } else if (next.front().get() == group.back().get()) {
2131 group.emplace_back(next.back());
2133 itr5 = Pairs.begin();
2134 } else if (next.back().get() == group.front().get()) {
2135 group.emplace_front(next.front());
2137 itr5 = Pairs.begin();
2138 } else if (next.back().get() == group.back().get()) {
2139 group.emplace_back(next.front());
2141 itr5 = Pairs.begin();
2147 if (itr5 == Pairs.end()) {
2148 for (itr6 = group.begin()+1; itr6 != group.end(); ++itr6) {
2149 pd->ReplaceCellPoint(itr6->get().f, itr6->get().g, group.front().get().g);
2162 loc->FreeSearchStructure();
2177 PolyAtEdge (vtkPolyData *_pd, vtkIdType _polyId, vtkIdType _ptIdA, vtkIdType _ptIdB) : pd(_pd), polyId(_polyId), ptIdA(_ptIdA), ptIdB(_ptIdB), loc(Loc::NONE) {
2179 double ptA[3], ptB[3];
2181 pd->GetPoint(ptIdA, ptA);
2182 pd->GetPoint(ptIdB, ptB);
2184 vtkMath::Subtract(ptB, ptA, e);
2185 vtkMath::Normalize(e);
2187 const vtkIdType *poly;
2190 pd->GetCellPoints(polyId, numPts, poly);
2192 ComputeNormal(pd->GetPoints(), n, numPts, poly);
2194 vtkMath::Cross(e, n, r);
2198 vtkIdType polyId, ptIdA, ptIdB;
2199 double n[3], e[3], r[3];
2203 friend std::ostream& operator<< (std::ostream &out, const PolyAtEdge &p) {
2204 out << "polyId " << p.polyId << ", ptIdA " << p.ptIdA << ", ptIdB " << p.ptIdB;
2208 static constexpr double eps {.99999999}; // ~.0081deg
2210 Congr IsCongruent (const PolyAtEdge &p) const {
2211 double cong = vtkMath::Dot(n, p.n);
2213 if (cong > eps || cong < -eps) {
2214 double ang = vtkMath::Dot(r, p.r);
2218 // normalen sind gleich ausgerichtet
2219 return Congr::EQUAL;
2221 return Congr::OPPOSITE;
2234 PolyPair (PolyAtEdge _pA, PolyAtEdge _pB) : pA(_pA), pB(_pB) {}
2238 void GetLoc (PolyAtEdge &pT, vtkIdType mode) {
2240 Congr cA = pA.IsCongruent(pT),
2241 cB = pB.IsCongruent(pT);
2244 std::cout << "GetLoc() -> polyId " << pT.polyId
2249 if (cA != Congr::NOT || cB != Congr::NOT) {
2255 if (cA == Congr::EQUAL || cA == Congr::OPPOSITE) {
2256 if (cA == Congr::OPPOSITE) {
2257 // normalen sind entgegengesetzt gerichtet
2259 if (mode == OPER_INTERSECTION) {
2260 pA.loc = Loc::OUTSIDE;
2261 pT.loc = Loc::OUTSIDE;
2263 pA.loc = Loc::INSIDE;
2264 pT.loc = Loc::INSIDE;
2266 } else if (mode == OPER_UNION || mode == OPER_INTERSECTION) {
2267 pA.loc = Loc::INSIDE;
2268 pT.loc = Loc::OUTSIDE;
2271 } else if (cB == Congr::EQUAL || cB == Congr::OPPOSITE) {
2272 if (cB == Congr::OPPOSITE) {
2273 // normalen sind entgegengesetzt gerichtet
2275 if (mode == OPER_INTERSECTION) {
2276 pB.loc = Loc::OUTSIDE;
2277 pT.loc = Loc::OUTSIDE;
2279 pB.loc = Loc::INSIDE;
2280 pT.loc = Loc::INSIDE;
2282 } else if (mode == OPER_UNION || mode == OPER_INTERSECTION) {
2283 pB.loc = Loc::INSIDE;
2284 pT.loc = Loc::OUTSIDE;
2288 double alpha = GetAngle(pA.r, pB.r, pA.e),
2289 beta = GetAngle(pA.r, pT.r, pA.e);
2292 pT.loc = Loc::INSIDE;
2294 pT.loc = Loc::OUTSIDE;
2302 std::shared_ptr<PolyPair> GetEdgePolys (vtkPolyData *pd, vtkIdList *ptsA, vtkIdList *ptsB) {
2305 std::cout << "GetEdgePolys()" << std::endl;
2308 std::vector<Pair> p;
2310 vtkIdType numPtsA = ptsA->GetNumberOfIds(),
2311 numPtsB = ptsB->GetNumberOfIds();
2313 vtkIdList *polys = vtkIdList::New();
2315 vtkIdType i, j, numCells;
2317 for (i = 0; i < numPtsA; i++) {
2318 pd->GetPointCells(ptsA->GetId(i), polys);
2319 numCells = polys->GetNumberOfIds();
2321 for (j = 0; j < numCells; j++) {
2322 p.emplace_back(ptsA->GetId(i), polys->GetId(j));
2326 for (i = 0; i < numPtsB; i++) {
2327 pd->GetPointCells(ptsB->GetId(i), polys);
2328 numCells = polys->GetNumberOfIds();
2330 for (j = 0; j < numCells; j++) {
2331 p.emplace_back(ptsB->GetId(i), polys->GetId(j));
2337 std::map<vtkIdType, IdsType> pEdges;
2339 std::vector<Pair>::const_iterator itr;
2340 for (itr = p.begin(); itr != p.end(); ++itr) {
2341 pEdges[itr->g].push_back(itr->f);
2344 std::vector<PolyAtEdge> opp;
2346 vtkIdList *poly = vtkIdList::New();
2348 vtkIdType numPts, a, b;
2350 std::map<vtkIdType, IdsType>::const_iterator itr2;
2352 for (itr2 = pEdges.begin(); itr2 != pEdges.end(); ++itr2) {
2353 const IdsType &pts = itr2->second;
2355 if (pts.size() > 1) {
2356 pd->GetCellPoints(itr2->first, poly);
2357 numPts = poly->GetNumberOfIds();
2359 for (i = 0; i < numPts; i++) {
2361 b = i+1 == numPts ? poly->GetId(0) : poly->GetId(i+1);
2363 if (std::find(pts.begin(), pts.end(), a) != pts.end()
2364 && std::find(pts.begin(), pts.end(), b) != pts.end()) {
2366 opp.emplace_back(pd, itr2->first, a, b);
2376 for (auto &op : opp) {
2377 std::cout << op << std::endl;
2381 if (opp.size() != 2) {
2385 return std::make_shared<PolyPair>(opp[0], opp[1]);
2389 bool vtkPolyDataBooleanFilter::CombineRegions () {
2392 std::cout << "CombineRegions()" << std::endl;
2395 vtkSmartPointer<vtkPolyData> filterdA = vtkSmartPointer<vtkPolyData>::New();
2396 filterdA->DeepCopy(modPdA);
2398 vtkSmartPointer<vtkPolyData> filterdB = vtkSmartPointer<vtkPolyData>::New();
2399 filterdB->DeepCopy(modPdB);
2401 // ungenutzte punkte löschen
2402 vtkSmartPointer<vtkCleanPolyData> cleanA = vtkSmartPointer<vtkCleanPolyData>::New();
2403 cleanA->PointMergingOff();
2404 cleanA->SetInputData(filterdA);
2406 vtkSmartPointer<vtkCleanPolyData> cleanB = vtkSmartPointer<vtkCleanPolyData>::New();
2407 cleanB->PointMergingOff();
2408 cleanB->SetInputData(filterdB);
2410 // regionen mit skalaren ausstatten
2411 vtkSmartPointer<vtkPolyDataConnectivityFilter> cfA = vtkSmartPointer<vtkPolyDataConnectivityFilter>::New();
2412 cfA->SetExtractionModeToAllRegions();
2413 cfA->ColorRegionsOn();
2414 cfA->SetInputConnection(cleanA->GetOutputPort());
2416 vtkSmartPointer<vtkPolyDataConnectivityFilter> cfB = vtkSmartPointer<vtkPolyDataConnectivityFilter>::New();
2417 cfB->SetExtractionModeToAllRegions();
2418 cfB->ColorRegionsOn();
2419 cfB->SetInputConnection(cleanB->GetOutputPort());
2424 vtkPolyData *pdA = cfA->GetOutput();
2425 vtkPolyData *pdB = cfB->GetOutput();
2428 std::cout << "Exporting modPdA_8.vtk" << std::endl;
2429 WriteVTK("modPdA_8.vtk", pdA);
2431 std::cout << "Exporting modPdB_8.vtk" << std::endl;
2432 WriteVTK("modPdB_8.vtk", pdB);
2435 resultC->ShallowCopy(contLines);
2437 if (OperMode == OPER_NONE) {
2438 resultA->ShallowCopy(pdA);
2439 resultB->ShallowCopy(pdB);
2444 vtkSmartPointer<vtkKdTreePointLocator> plA = vtkSmartPointer<vtkKdTreePointLocator>::New();
2445 plA->SetDataSet(pdA);
2446 plA->BuildLocator();
2448 vtkSmartPointer<vtkKdTreePointLocator> plB = vtkSmartPointer<vtkKdTreePointLocator>::New();
2449 plB->SetDataSet(pdB);
2450 plB->BuildLocator();
2455 vtkIdTypeArray *scalarsA = vtkIdTypeArray::SafeDownCast(pdA->GetPointData()->GetScalars());
2456 vtkIdTypeArray *scalarsB = vtkIdTypeArray::SafeDownCast(pdB->GetPointData()->GetScalars());
2458 vtkSmartPointer<vtkIdList> line = vtkSmartPointer<vtkIdList>::New();
2460 double ptA[3], ptB[3];
2462 vtkSmartPointer<vtkIdList> fptsA = vtkSmartPointer<vtkIdList>::New();
2463 vtkSmartPointer<vtkIdList> lptsA = vtkSmartPointer<vtkIdList>::New();
2465 vtkSmartPointer<vtkIdList> fptsB = vtkSmartPointer<vtkIdList>::New();
2466 vtkSmartPointer<vtkIdList> lptsB = vtkSmartPointer<vtkIdList>::New();
2468 std::map<vtkIdType, Loc> locsA, locsB;
2470 vtkIdType i, j, numLines = contLines->GetNumberOfCells();
2474 for (i = 0; i < numLines; i++) {
2476 if (contLines->GetCellType(i) == VTK_EMPTY_CELL) {
2480 contLines->GetCellPoints(i, line);
2482 contLines->GetPoint(line->GetId(0), ptA);
2483 contLines->GetPoint(line->GetId(1), ptB);
2485 FindPoints(plA, ptA, fptsA);
2486 FindPoints(plB, ptA, fptsB);
2489 std::cout << "line " << i << std::endl;
2492 // bereits behandelte regionen werden nicht noch einmal untersucht
2494 vtkIdType notLocated = 0;
2496 for (j = 0; j < fptsA->GetNumberOfIds(); j++) {
2497 if (locsA.count(scalarsA->GetValue(fptsA->GetId(j))) == 0) {
2502 for (j = 0; j < fptsB->GetNumberOfIds(); j++) {
2503 if (locsB.count(scalarsB->GetValue(fptsB->GetId(j))) == 0) {
2508 if (notLocated == 0) {
2514 FindPoints(plA, ptB, lptsA);
2515 FindPoints(plB, ptB, lptsB);
2517 auto ppA = GetEdgePolys(pdA, fptsA, lptsA);
2518 auto ppB = GetEdgePolys(pdB, fptsB, lptsB);
2522 ppB->GetLoc(ppA->pA, OperMode);
2523 ppB->GetLoc(ppA->pB, OperMode);
2525 ppA->GetLoc(ppB->pA, OperMode);
2526 ppA->GetLoc(ppB->pB, OperMode);
2528 vtkIdType fsA = scalarsA->GetValue(ppA->pA.ptIdA);
2529 vtkIdType lsA = scalarsA->GetValue(ppA->pB.ptIdA);
2531 vtkIdType fsB = scalarsB->GetValue(ppB->pA.ptIdA);
2532 vtkIdType lsB = scalarsB->GetValue(ppB->pB.ptIdA);
2535 std::cout << "polyId " << ppA->pA.polyId << ", sA " << fsA << ", loc " << ppA->pA.loc << std::endl;
2536 std::cout << "polyId " << ppA->pB.polyId << ", sA " << lsA << ", loc " << ppA->pB.loc << std::endl;
2537 std::cout << "polyId " << ppB->pA.polyId << ", sB " << fsB << ", loc " << ppB->pA.loc << std::endl;
2538 std::cout << "polyId " << ppB->pB.polyId << ", sB " << lsB << ", loc " << ppB->pB.loc << std::endl;
2540 if (locsA.count(fsA) > 0 && locsA[fsA] != ppA->pA.loc) {
2541 std::cout << "sA " << fsA << ": " << locsA[fsA] << " -> " << ppA->pA.loc << std::endl;
2544 if (locsA.count(lsA) > 0 && locsA[lsA] != ppA->pB.loc) {
2545 std::cout << "sA " << lsA << ": " << locsA[lsA] << " -> " << ppA->pB.loc << std::endl;
2548 if (locsB.count(fsB) > 0 && locsB[fsB] != ppB->pA.loc) {
2549 std::cout << "sB " << fsB << ": " << locsB[fsB] << " -> " << ppB->pA.loc << std::endl;
2552 if (locsB.count(lsB) > 0 && locsB[lsB] != ppB->pB.loc) {
2553 std::cout << "sB " << lsB << ": " << locsB[lsB] << " -> " << ppB->pB.loc << std::endl;
2558 locsA.emplace(fsA, ppA->pA.loc);
2559 locsA.emplace(lsA, ppA->pB.loc);
2561 locsB.emplace(fsB, ppB->pA.loc);
2562 locsB.emplace(lsB, ppB->pB.loc);
2565 _failed.push_back(i);
2572 if (_failed.size() > 0) {
2573 for (auto i : _failed) {
2574 std::cout << "failed at " << i
2581 // reale kombination der ermittelten regionen
2583 Loc comb[] = {Loc::OUTSIDE, Loc::OUTSIDE};
2585 if (OperMode == OPER_INTERSECTION) {
2586 comb[0] = Loc::INSIDE;
2587 comb[1] = Loc::INSIDE;
2588 } else if (OperMode == OPER_DIFFERENCE) {
2589 comb[1] = Loc::INSIDE;
2590 } else if (OperMode == OPER_DIFFERENCE2) {
2591 comb[0] = Loc::INSIDE;
2594 vtkIdType numA = cfA->GetNumberOfExtractedRegions(),
2595 numB = cfB->GetNumberOfExtractedRegions();
2597 cfA->SetExtractionModeToSpecifiedRegions();
2598 cfB->SetExtractionModeToSpecifiedRegions();
2600 std::map<vtkIdType, Loc>::const_iterator itr;
2602 for (itr = locsA.begin(); itr != locsA.end(); itr++) {
2603 if (itr->second == comb[0]) {
2604 cfA->AddSpecifiedRegion(itr->first);
2608 for (itr = locsB.begin(); itr != locsB.end(); itr++) {
2609 if (itr->second == comb[1]) {
2610 cfB->AddSpecifiedRegion(itr->first);
2614 // nicht beteiligte regionen hinzufĂ¼gen
2616 if (OperMode == OPER_UNION || OperMode == OPER_DIFFERENCE) {
2617 for (i = 0; i < numA; i++) {
2618 if (locsA.count(i) == 0) {
2619 cfA->AddSpecifiedRegion(i);
2624 if (OperMode == OPER_UNION || OperMode == OPER_DIFFERENCE2) {
2625 for (i = 0; i < numB; i++) {
2626 if (locsB.count(i) == 0) {
2627 cfB->AddSpecifiedRegion(i);
2632 // nach innen zeigende normalen umkehren
2637 vtkPolyData *regsA = cfA->GetOutput();
2638 vtkPolyData *regsB = cfB->GetOutput();
2640 scalarsA = vtkIdTypeArray::SafeDownCast(regsA->GetPointData()->GetScalars());
2641 scalarsB = vtkIdTypeArray::SafeDownCast(regsB->GetPointData()->GetScalars());
2643 if (OperMode != OPER_INTERSECTION) {
2644 if (comb[0] == Loc::INSIDE) {
2645 for (i = 0; i < regsA->GetNumberOfCells(); i++) {
2646 if (locsA.count(scalarsA->GetValue(i)) == 1) {
2647 regsA->ReverseCell(i);
2652 if (comb[1] == Loc::INSIDE) {
2653 for (i = 0; i < regsB->GetNumberOfCells(); i++) {
2654 if (locsB.count(scalarsB->GetValue(i)) == 1) {
2655 regsB->ReverseCell(i);
2661 // OrigCellIds und CellData
2663 vtkIdTypeArray *origCellIdsA = vtkIdTypeArray::SafeDownCast(regsA->GetCellData()->GetScalars("OrigCellIds"));
2664 vtkIdTypeArray *origCellIdsB = vtkIdTypeArray::SafeDownCast(regsB->GetCellData()->GetScalars("OrigCellIds"));
2666 vtkIdTypeArray *newOrigCellIdsA = vtkIdTypeArray::New();
2667 newOrigCellIdsA->SetName("OrigCellIdsA");
2669 vtkIdTypeArray *newOrigCellIdsB = vtkIdTypeArray::New();
2670 newOrigCellIdsB->SetName("OrigCellIdsB");
2672 vtkCellData *newCellDataA = vtkCellData::New();
2673 newCellDataA->CopyAllocate(cellDataA);
2675 vtkCellData *newCellDataB = vtkCellData::New();
2676 newCellDataB->CopyAllocate(cellDataB);
2680 for (i = 0; i < regsA->GetNumberOfCells(); i++) {
2681 cellId = cellIdsA->GetValue(origCellIdsA->GetValue(i));
2683 newOrigCellIdsA->InsertNextValue(cellId);
2684 newOrigCellIdsB->InsertNextValue(-1);
2686 newCellDataA->CopyData(cellDataA, cellId, i);
2689 for (i = 0; i < regsB->GetNumberOfCells(); i++) {
2690 cellId = cellIdsB->GetValue(origCellIdsB->GetValue(i));
2692 newOrigCellIdsB->InsertNextValue(cellId);
2693 newOrigCellIdsA->InsertNextValue(-1);
2695 newCellDataB->CopyData(cellDataB, cellId, i);
2698 regsA->GetCellData()->Initialize();
2699 regsB->GetCellData()->Initialize();
2701 regsA->GetCellData()->ShallowCopy(newCellDataA);
2702 regsB->GetCellData()->ShallowCopy(newCellDataB);
2704 newCellDataA->Delete();
2705 newCellDataB->Delete();
2709 vtkAppendPolyData *app = vtkAppendPolyData::New();
2710 app->AddInputData(regsA);
2711 app->AddInputData(regsB);
2713 // entfernt ungenutzte punkte
2714 vtkCleanPolyData *cleanApp = vtkCleanPolyData::New();
2715 cleanApp->PointMergingOff();
2716 cleanApp->SetInputConnection(app->GetOutputPort());
2718 // färbt die regionen nochmal neu ein, damit mehrere regionen nicht die gleiche farbe haben
2720 vtkPolyDataConnectivityFilter *cfApp = vtkPolyDataConnectivityFilter::New();
2721 cfApp->SetExtractionModeToAllRegions();
2722 cfApp->ColorRegionsOn();
2723 cfApp->SetInputConnection(cleanApp->GetOutputPort());
2727 vtkPolyData *cfPd = cfApp->GetOutput();
2729 // resultB bleibt hier leer
2731 resultA->ShallowCopy(cfPd);
2733 resultA->GetCellData()->AddArray(newOrigCellIdsA);
2734 resultA->GetCellData()->AddArray(newOrigCellIdsB);
2742 newOrigCellIdsB->Delete();
2743 newOrigCellIdsA->Delete();
2745 plB->FreeSearchStructure();
2746 plA->FreeSearchStructure();
2752 Merger::Merger (vtkPolyData *pd, const PStrips &pStrips, const StripsType &strips, const IdsType &descIds, vtkIdType origId) : pd(pd), pStrips(pStrips), origId(origId) {
2754 const StripPtsType &pts = pStrips.pts;
2755 const Base &base = pStrips.base;
2758 const vtkIdType *cell;
2762 for (auto id : descIds) {
2763 pd->GetCellPoints(id, num, cell);
2767 for (i = 0; i < num; i++) {
2768 pd->GetPoint(cell[i], pt);
2771 Transform(pt, proj, base);
2773 p.emplace_back(proj[0], proj[1], 0, cell[i]);
2782 for (auto &strip : strips) {
2785 for (auto &sp : strip) {
2786 const double *pt = pts.at(sp.ind).pt;
2789 Transform(pt, proj, base);
2791 p.emplace_back(proj[0], proj[1], 0);
2797 ComputeNormal(p, n);
2800 Poly q(p.rbegin(), p.rend());
2809 void Merger::Run () {
2810 // mergen mit hilfe von vtkKdTree und vtkModifiedBSPTree
2812 vtkPoints *pdPts = pd->GetPoints();
2813 vtkIdTypeArray *origCellIds = vtkIdTypeArray::SafeDownCast(pd->GetCellData()->GetScalars("OrigCellIds"));
2815 assert(origCellIds != nullptr);
2817 const Base &base = pStrips.base;
2819 std::vector<GroupType> groups(polys.size());
2821 PolysType::const_iterator itrA, itrB;
2823 for (itrA = polys.begin(); itrA != polys.end(); ++itrA) {
2824 for (itrB = polys.begin(); itrB != polys.end(); ++itrB) {
2825 if (itrA != itrB && PointInPoly(*itrB, *itrA->begin())) {
2826 groups[static_cast<std::size_t>(itrB-polys.begin())].push_back(static_cast<std::size_t>(itrA-polys.begin()));
2831 std::size_t parent = 0;
2833 for (auto &group : groups) {
2836 for (auto &index : group) {
2837 const GroupType &_group = groups[index];
2838 parents.insert(parents.end(), _group.begin(), _group.end());
2841 std::sort(group.begin(), group.end());
2842 std::sort(parents.begin(), parents.end());
2844 GroupType _group {parent++};
2845 std::set_difference(group.begin(), group.end(), parents.begin(), parents.end(), std::back_inserter(_group));
2848 for (auto &index : _group) {
2849 std::cout << index << ", ";
2851 std::cout << "]" << std::endl;
2855 MergeGroup(_group, merged);
2857 std::map<Point3d, vtkIdType> newIds;
2859 for (auto &poly : merged) {
2860 vtkSmartPointer<vtkIdList> newCell = vtkSmartPointer<vtkIdList>::New();
2862 for (auto &p : poly) {
2863 double in[] = {p.x, p.y},
2866 BackTransform(in, out, base);
2868 vtkIdType id = p.id;
2871 auto itr = newIds.find(p);
2873 if (itr == newIds.end()) {
2874 id = pdPts->InsertNextPoint(out);
2875 newIds.emplace(p, id);
2881 newCell->InsertNextId(id);
2885 pd->InsertNextCell(VTK_POLYGON, newCell);
2886 origCellIds->InsertNextValue(origId);
2893 void Merger::MergeGroup (const GroupType &group, PolysType &merged) {
2894 if (group.size() == 1) {
2895 merged.push_back(polys.at(group.back()));
2900 vtkSmartPointer<vtkPoints> pts = vtkSmartPointer<vtkPoints>::New();
2902 IndexedPolysType indexedPolys;
2904 ReferencedPointsType refPts;
2906 SourcesType sources;
2907 std::size_t src = 0;
2909 for (auto &index : group) {
2910 const Poly &poly = polys.at(index);
2914 for (auto &p : poly) {
2915 vtkIdType id = pts->InsertNextPoint(p.x, p.y, p.z);
2918 sources.emplace(id, src);
2920 refPts.emplace(id, p);
2923 indexedPolys.push_back(std::move(ids));
2927 vtkSmartPointer<vtkKdTree> kdTree = vtkSmartPointer<vtkKdTree>::New();
2928 kdTree->OmitZPartitioning();
2929 kdTree->BuildLocatorFromPoints(pts);
2931 vtkSmartPointer<vtkPolyData> linesA = vtkSmartPointer<vtkPolyData>::New();
2932 linesA->SetPoints(pts);
2933 linesA->Allocate(1);
2935 IndexedPoly::const_iterator itrA, itrB;
2937 for (const auto &ids : indexedPolys) {
2938 for (itrA = ids.begin(); itrA != ids.end(); ++itrA) {
2940 if (itrB == ids.end()) {
2944 vtkIdList *line = vtkIdList::New();
2945 line->InsertNextId(*itrA);
2946 line->InsertNextId(*itrB);
2948 linesA->InsertNextCell(VTK_LINE, line);
2954 WriteVTK("linesA.vtk", linesA);
2956 vtkSmartPointer<vtkModifiedBSPTree> bspTreeA = vtkSmartPointer<vtkModifiedBSPTree>::New();
2957 bspTreeA->SetDataSet(linesA);
2961 PolyConnsType polyConns;
2963 FindConns(linesA, kdTree, bspTreeA, polyConns, indexedPolys, sources, n);
2965 PolyConnsType connected {{0, {}}};
2966 std::set<vtkIdType> restricted; // keine der conns darf im gleichen punkt beginnen
2968 vtkSmartPointer<vtkPolyData> linesB = vtkSmartPointer<vtkPolyData>::New();
2969 linesB->SetPoints(pts);
2970 linesB->Allocate(1);
2972 vtkSmartPointer<vtkModifiedBSPTree> bspTreeB = vtkSmartPointer<vtkModifiedBSPTree>::New();
2973 bspTreeB->SetDataSet(linesB);
2975 ConnsType firstConns;
2977 std::size_t i, numPolys = indexedPolys.size();
2979 double ptA[3], ptB[3];
2981 while (connected.size() < numPolys) {
2983 bool foundOne = false;
2985 for (i = 1; i < numPolys; i++) {
2986 if (connected.count(i) == 0) {
2987 const ConnsType &conns = polyConns[i];
2989 for (auto &conn : conns) {
2990 if (connected.count(sources.at(conn.j)) == 1
2991 && restricted.count(conn.j) == 0) {
2993 pts->GetPoint(conn.i, ptA);
2994 pts->GetPoint(conn.j, ptB);
2996 if (bspTreeB->IntersectWithLine(ptA, ptB, 1e-5, nullptr, nullptr) == 0) {
2997 connected[sources.at(conn.i)].push_back(conn);
2999 // das andere poly auch aktualisieren
3000 connected[sources.at(conn.j)].emplace_back(conn.d, conn.j, conn.i);
3002 restricted.insert(conn.i);
3003 restricted.insert(conn.j);
3005 vtkIdList *line = vtkIdList::New();
3006 line->InsertNextId(conn.i);
3007 line->InsertNextId(conn.j);
3009 linesB->InsertNextCell(VTK_LINE, line);
3013 bspTreeB->Modified();
3017 firstConns.push_back(conn);
3028 if (!FindConns(linesA, kdTree, bspTreeA, polyConns, indexedPolys, sources, n)) {
3029 throw std::runtime_error("Merging failed.");
3034 std::map<std::size_t, std::vector<std::size_t>> chains;
3036 PolyConnsType::const_iterator itrC;
3038 for (itrC = connected.begin(); itrC != connected.end(); ++itrC) {
3039 auto &chain = chains[itrC->first];
3040 chain.push_back(itrC->first);
3042 while (chain.back() != 0) {
3043 chain.push_back(sources.at(connected.at(chain.back()).front().j));
3047 std::cout << connected;
3049 decltype(chains)::const_iterator itrD;
3051 for (itrD = chains.begin(); itrD != chains.end(); ++itrD) {
3052 std::cout << itrD->first << ": [";
3053 for (auto &id : itrD->second) {
3054 std::cout << id << ", ";
3056 std::cout << "]" << std::endl;
3059 std::set<std::size_t> solved {0};
3061 std::deque<std::size_t> searchInds;
3063 for (i = 1; i < numPolys; i++) {
3064 if (connected.at(i).size() == 1) {
3065 searchInds.push_back(i);
3069 while (!searchInds.empty()) {
3072 for (auto ind : searchInds) {
3073 PolyPriosType polyPrios;
3075 std::cout << "ind " << ind << std::endl;
3077 const Conn &first = connected.at(ind).back();
3079 for (auto &conn : polyConns.at(ind)) {
3080 auto &src = sources.at(conn.j);
3082 if (polyPrios.count(src) == 1) {
3086 if (conn.i != first.i
3087 && conn.j != first.j
3088 && restricted.count(conn.i) == 0
3089 && restricted.count(conn.j) == 0) {
3091 pts->GetPoint(conn.i, ptA);
3092 pts->GetPoint(conn.j, ptB);
3094 if (bspTreeB->IntersectWithLine(ptA, ptB, 1e-5, nullptr, nullptr) == 0) {
3095 auto &chainA = chains.at(ind),
3096 &chainB = chains.at(src);
3098 std::set<std::size_t> _chainA(chainA.begin(), chainA.end()),
3099 _chainB(chainB.begin(), chainB.end());
3101 // gemeinsame eltern
3102 std::set<std::size_t> shared;
3104 std::set_intersection(_chainA.begin(), _chainA.end(), _chainB.begin(), _chainB.end(), std::inserter(shared, shared.end()));
3106 // gemeinsame eltern mĂ¼ssen sich alle in solved befinden
3107 if (std::includes(solved.begin(), solved.end(), shared.begin(), shared.end())) {
3108 std::set<std::size_t> solvable;
3110 std::set_difference(_chainA.begin(), _chainA.end(), solved.begin(), solved.end(), std::inserter(solvable, solvable.end()));
3111 std::set_difference(_chainB.begin(), _chainB.end(), solved.begin(), solved.end(), std::inserter(solvable, solvable.end()));
3113 polyPrios.emplace(std::piecewise_construct,
3114 std::forward_as_tuple(src),
3115 std::forward_as_tuple(conn, solvable, -conn.d));
3121 PolyPriosType::const_iterator itr;
3122 for (itr = polyPrios.begin(); itr != polyPrios.end(); ++itr) {
3123 prios.insert(itr->second);
3127 if (!prios.empty()) {
3128 auto &prio = *prios.rbegin();
3130 std::cout << "found " << prio << std::endl;
3132 auto &conns = connected.at(sources.at(prio.conn.i));
3134 conns.push_back(prio.conn);
3136 connected.at(sources.at(prio.conn.j)).emplace_back(prio.conn.d, prio.conn.j, prio.conn.i);
3138 restricted.insert(prio.conn.i);
3139 restricted.insert(prio.conn.j);
3141 vtkIdList *line = vtkIdList::New();
3142 line->InsertNextId(prio.conn.i);
3143 line->InsertNextId(prio.conn.j);
3145 linesB->InsertNextCell(VTK_LINE, line);
3149 bspTreeB->Modified();
3151 solved.insert(prio.solvable.begin(), prio.solvable.end());
3153 searchInds.erase(std::find(searchInds.begin(), searchInds.end(), sources.at(prio.conn.i)));
3155 auto itr = std::find(searchInds.begin(), searchInds.end(), sources.at(prio.conn.j));
3157 if (itr != searchInds.end()) {
3158 searchInds.erase(itr);
3161 if (!FindConns(linesA, kdTree, bspTreeA, polyConns, indexedPolys, sources, n)) {
3162 throw std::runtime_error("Merging failed.");
3167 std::cout << connected;
3169 WriteVTK("linesB.vtk", linesB);
3174 bool operator() (const Conn &a, const Conn &b) const {
3175 return std::tie(a.i, a.j) < std::tie(b.i, b.j);
3179 std::set<Conn, Cmp> usedConns;
3181 IndexedPoly polyA {indexedPolys.front()};
3183 for (const auto &first : firstConns) {
3184 auto itrA = std::find(polyA.begin(), polyA.end(), first.j);
3186 assert(itrA != polyA.end());
3188 auto &polyB = indexedPolys.at(sources.at(first.i));
3190 auto itrB = std::find(polyB.begin(), polyB.end(), first.i);
3192 assert(itrB != polyB.end());
3194 std::rotate(polyA.begin(), itrA, polyA.end());
3195 std::rotate(polyB.begin(), itrB, polyB.end());
3197 IndexedPoly newPoly {polyA};
3198 newPoly.push_back(polyA.front());
3199 newPoly.push_back(polyB.front());
3201 newPoly.insert(newPoly.end(), polyB.rbegin(), polyB.rend());
3203 polyA.swap(newPoly);
3205 usedConns.insert(first);
3209 PolysType newPolysA;
3210 GetPolys(refPts, {polyA}, newPolysA);
3212 WritePolys("merged_stage1.vtk", newPolysA);
3216 std::set<Conn, Cmp> leftConns;
3218 for (itrC = connected.begin(); itrC != connected.end(); ++itrC) {
3219 if (itrC->first == 0) {
3223 auto &conns = itrC->second;
3225 ConnsType::const_iterator itr;
3227 for (itr = conns.begin()+1; itr != conns.end(); ++itr) {
3228 Conn conn(0, itr->j, itr->i);
3230 if (usedConns.find(conn) == usedConns.end()) {
3231 if (itr->i < itr->j) {
3232 leftConns.emplace(0, itr->i, itr->j);
3234 leftConns.insert(std::move(conn));
3241 std::cout << "leftConns: [";
3242 for (auto &conn : leftConns) {
3243 std::cout << conn << ", ";
3245 std::cout << "]" << std::endl;
3247 IndexedPolysType splitted {polyA};
3249 for (auto &conn : leftConns) {
3250 IndexedPolysType::iterator itr;
3252 for (itr = splitted.begin(); itr != splitted.end(); ++itr) {
3255 auto itrA = std::find(poly.begin(), poly.end(), conn.i);
3257 if (itrA != poly.end()) {
3258 std::rotate(poly.begin(), itrA, poly.end());
3260 auto itrB = std::find(poly.begin(), poly.end(), conn.j);
3262 IndexedPoly newPolyA(poly.begin(), itrB+1);
3263 IndexedPoly newPolyB(itrB, poly.end());
3265 newPolyB.push_back(poly.front());
3267 splitted.erase(itr);
3269 splitted.push_back(std::move(newPolyA));
3270 splitted.push_back(std::move(newPolyB));
3277 PolysType newPolysB;
3278 GetPolys(refPts, splitted, newPolysB);
3280 WritePolys("merged_stage2.vtk", newPolysB);
3282 std::move(newPolysB.begin(), newPolysB.end(), std::back_inserter(merged));
3286 bool Merger::FindConns (vtkPolyData *lines, vtkSmartPointer<vtkKdTree> kdTree, vtkSmartPointer<vtkModifiedBSPTree> bspTree, PolyConnsType &polyConns, const IndexedPolysType &indexedPolys, const SourcesType &sources, int &n) {
3288 vtkPoints *pts = lines->GetPoints();
3290 if (n > pts->GetNumberOfPoints()) {
3296 vtkSmartPointer<vtkIdList> foundPts = vtkSmartPointer<vtkIdList>::New();
3298 vtkIdType i, numPts;
3302 vtkSmartPointer<vtkIdList> lineIds = vtkSmartPointer<vtkIdList>::New();
3304 double ptA[3], ptB[3];
3309 vtkIdType _idA, _idB;
3311 std::map<std::size_t, std::set<Conn, ConnCmp>> _polyConns;
3313 vtkSmartPointer<vtkIdList> line = vtkSmartPointer<vtkIdList>::New();
3315 for (const auto &ids : indexedPolys) {
3316 for (vtkIdType idA : ids) {
3317 pts->GetPoint(idA, ptA);
3319 kdTree->FindClosestNPoints(n, ptA, foundPts);
3321 numPts = foundPts->GetNumberOfIds();
3323 for (i = 0; i < numPts; i++) {
3324 idB = foundPts->GetId(i);
3326 auto srcA = sources.at(idA),
3327 srcB = sources.at(idB);
3333 pts->GetPoint(idB, ptB);
3337 if (bspTree->IntersectWithLine(ptA, ptB, 1e-5, nullptr, lineIds) == 1) {
3338 for (j = 0; j < lineIds->GetNumberOfIds(); j++) {
3339 lines->GetCellPoints(lineIds->GetId(j), line);
3341 _idA = line->GetId(0);
3342 _idB = line->GetId(1);
3344 if (_idA != idA && _idA != idB
3345 && _idB != idA && _idB != idB) {
3354 double d = vtkMath::Distance2BetweenPoints(ptA, ptB);
3356 _polyConns[srcA].emplace(d, idA, idB);
3357 _polyConns[srcB].emplace(d, idB, idA);
3363 decltype(_polyConns)::const_iterator itr;
3365 for (itr = _polyConns.begin(); itr != _polyConns.end(); ++itr) {
3366 auto &_conns = itr->second;
3368 ConnsType conns(_conns.begin(), _conns.end());
3369 std::sort(conns.begin(), conns.end());
3371 polyConns[itr->first].swap(conns);