lib/creaVtk/vtkPolyDataBooleanFilter.cxx

   1 /*
   2 Copyright 2012-2022 Ronald Römer
   3
   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
   7
   8     http://www.apache.org/licenses/LICENSE-2.0
   9
  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.
  15 */
  16
  17 #include <map>
  18 #include <deque>
  19 #include <vector>
  20 #include <set>
  21 #include <algorithm>
  22 #include <cmath>
  23 #include <functional>
  24 #include <queue>
  25 #include <memory>
  26
  27 #include <chrono>
  28 #include <numeric>
  29 #include <iterator>
  30
  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>
  38 #include <vtkMath.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>
  48
  49 #include "vtkPolyDataBooleanFilter.h"
  50 #include "vtkPolyDataContactFilter.h"
  51
  52 #include "Utilities.h"
  53
  54 vtkStandardNewMacro(vtkPolyDataBooleanFilter);
  55
  56 vtkPolyDataBooleanFilter::vtkPolyDataBooleanFilter () {
  57
  58     SetNumberOfInputPorts(2);
  59     SetNumberOfOutputPorts(3);
  60
  61     timePdA = 0;
  62     timePdB = 0;
  63
  64     contLines = vtkSmartPointer<vtkPolyData>::New();
  65
  66     modPdA = vtkSmartPointer<vtkPolyData>::New();
  67     modPdB = vtkSmartPointer<vtkPolyData>::New();
  68
  69     cellDataA = vtkSmartPointer<vtkCellData>::New();
  70     cellDataB = vtkSmartPointer<vtkCellData>::New();
  71
  72     cellIdsA = vtkSmartPointer<vtkIdTypeArray>::New();
  73     cellIdsB = vtkSmartPointer<vtkIdTypeArray>::New();
  74
  75     OperMode = OPER_UNION;
  76
  77 }
  78
  79 vtkPolyDataBooleanFilter::~vtkPolyDataBooleanFilter () {
  80     // nix mehr
  81 }
  82
  83 int vtkPolyDataBooleanFilter::ProcessRequest(vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector) {
  84
  85     if (request->Has(vtkDemandDrivenPipeline::REQUEST_DATA())) {
  86
  87         vtkInformation *inInfoA = inputVector[0]->GetInformationObject(0);
  88         vtkInformation *inInfoB = inputVector[1]->GetInformationObject(0);
  89
  90         vtkPolyData *pdA = vtkPolyData::SafeDownCast(inInfoA->Get(vtkDataObject::DATA_OBJECT()));
  91         vtkPolyData *pdB = vtkPolyData::SafeDownCast(inInfoB->Get(vtkDataObject::DATA_OBJECT()));
  92
  93         vtkInformation *outInfoA = outputVector->GetInformationObject(0);
  94         vtkInformation *outInfoB = outputVector->GetInformationObject(1);
  95         vtkInformation *outInfoC = outputVector->GetInformationObject(2);
  96
  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()));
 100
 101         using clock = std::chrono::steady_clock;
 102         std::vector<clock::duration> times;
 103         clock::time_point start;
 104
 105         if (pdA->GetMTime() > timePdA || pdB->GetMTime() > timePdB) {
 106             // eventuell vorhandene regionen vereinen
 107
 108             vtkSmartPointer<vtkCleanPolyData> cleanA = vtkSmartPointer<vtkCleanPolyData>::New();
 109             cleanA->SetOutputPointsPrecision(DOUBLE_PRECISION);
 110             cleanA->SetTolerance(1e-6);
 111             cleanA->SetInputData(pdA);
 112             cleanA->Update();
 113
 114             vtkSmartPointer<vtkCleanPolyData> cleanB = vtkSmartPointer<vtkCleanPolyData>::New();
 115             cleanB->SetOutputPointsPrecision(DOUBLE_PRECISION);
 116             cleanB->SetTolerance(1e-6);
 117             cleanB->SetInputData(pdB);
 118             cleanB->Update();
 119
 120 #ifdef DEBUG
 121             std::cout << "Exporting modPdA.vtk" << std::endl;
 122             WriteVTK("modPdA.vtk", cleanA->GetOutput());
 123
 124             std::cout << "Exporting modPdB.vtk" << std::endl;
 125             WriteVTK("modPdB.vtk", cleanB->GetOutput());
 126 #endif
 127
 128             // CellData sichern
 129
 130             cellDataA->DeepCopy(cleanA->GetOutput()->GetCellData());
 131             cellDataB->DeepCopy(cleanB->GetOutput()->GetCellData());
 132
 133             // ermittelt kontaktstellen
 134
 135             start = clock::now();
 136
 137             vtkSmartPointer<vtkPolyDataContactFilter> cl = vtkSmartPointer<vtkPolyDataContactFilter>::New();
 138             cl->SetInputConnection(0, cleanA->GetOutputPort());
 139             cl->SetInputConnection(1, cleanB->GetOutputPort());
 140             cl->Update();
 141
 142             times.push_back(clock::now()-start);
 143
 144             contLines->DeepCopy(cl->GetOutput());
 145
 146             modPdA->DeepCopy(cl->GetOutput(1));
 147             modPdB->DeepCopy(cl->GetOutput(2));
 148
 149 #ifdef DEBUG
 150             std::cout << "Exporting contLines.vtk" << std::endl;
 151             WriteVTK("contLines.vtk", contLines);
 152
 153             std::cout << "Exporting modPdA_1.vtk" << std::endl;
 154             WriteVTK("modPdA_1.vtk", modPdA);
 155
 156             std::cout << "Exporting modPdB_1.vtk" << std::endl;
 157             WriteVTK("modPdB_1.vtk", modPdB);
 158 #endif
 159
 160             if (contLines->GetNumberOfCells() == 0) {
 161                 vtkErrorMacro("Inputs have no contact.");
 162
 163                 return 1;
 164             }
 165
 166             // in den CellDatas steht drin, welche polygone einander schneiden
 167
 168             contsA = vtkIdTypeArray::SafeDownCast(contLines->GetCellData()->GetScalars("cA"));
 169             contsB = vtkIdTypeArray::SafeDownCast(contLines->GetCellData()->GetScalars("cB"));
 170
 171             vtkIdTypeArray *accuracy = vtkIdTypeArray::SafeDownCast(contLines->GetPointData()->GetScalars("accuracy"));
 172
 173             vtkIdTypeArray *sourcesA = vtkIdTypeArray::SafeDownCast(contLines->GetCellData()->GetScalars("sourcesA"));
 174             vtkIdTypeArray *sourcesB = vtkIdTypeArray::SafeDownCast(contLines->GetCellData()->GetScalars("sourcesB"));
 175
 176             vtkIdType i, numPts = contLines->GetNumberOfPoints();
 177
 178             vtkSmartPointer<vtkIdList> cells = vtkSmartPointer<vtkIdList>::New();
 179
 180             for (i = 0; i < numPts; i++) {
 181                 contLines->GetPointCells(i, cells);
 182
 183                 if (cells->GetNumberOfIds() == 1) {
 184                     vtkErrorMacro("Contact ends suddenly.");
 185                     return 1;
 186                 }
 187
 188                 if (accuracy->GetValue(i) == 1) {
 189                     vtkErrorMacro("Contact goes through a cell with a bad shape.");
 190                     return 1;
 191                 }
 192             }
 193
 194             // sichert die OrigCellIds
 195
 196             vtkIdTypeArray *origCellIdsA = vtkIdTypeArray::SafeDownCast(modPdA->GetCellData()->GetScalars("OrigCellIds"));
 197             vtkIdTypeArray *origCellIdsB = vtkIdTypeArray::SafeDownCast(modPdB->GetCellData()->GetScalars("OrigCellIds"));
 198
 199             cellIdsA->DeepCopy(origCellIdsA);
 200             cellIdsB->DeepCopy(origCellIdsB);
 201
 202             vtkIdType numCellsA = modPdA->GetNumberOfCells(),
 203                 numCellsB = modPdB->GetNumberOfCells();
 204
 205             for (i = 0; i < numCellsA; i++) {
 206                 origCellIdsA->SetValue(i, i);
 207             }
 208
 209             for (i = 0; i < numCellsB; i++) {
 210                 origCellIdsB->SetValue(i, i);
 211             }
 212
 213             start = clock::now();
 214
 215             if (GetPolyStrips(modPdA, contsA, sourcesA, polyStripsA) ||
 216                 GetPolyStrips(modPdB, contsB, sourcesB, polyStripsB)) {
 217
 218                 vtkErrorMacro("Strips are invalid.");
 219
 220                 return 1;
 221
 222             }
 223
 224             times.push_back(clock::now()-start);
 225
 226             // trennt die polygone an den linien
 227
 228             start = clock::now();
 229
 230             if (CutCells(modPdA, polyStripsA) ||
 231                 CutCells(modPdB, polyStripsB)) {
 232
 233                 vtkErrorMacro("CutCells failed.");
 234
 235                 return 1;
 236             }
 237
 238             times.push_back(clock::now()-start);
 239
 240 #ifdef DEBUG
 241             std::cout << "Exporting modPdA_2.vtk" << std::endl;
 242             WriteVTK("modPdA_2.vtk", modPdA);
 243
 244             std::cout << "Exporting modPdB_2.vtk" << std::endl;
 245             WriteVTK("modPdB_2.vtk", modPdB);
 246 #endif
 247
 248             start = clock::now();
 249
 250             RestoreOrigPoints(modPdA, polyStripsA);
 251             RestoreOrigPoints(modPdB, polyStripsB);
 252
 253             times.push_back(clock::now()-start);
 254
 255 #ifdef DEBUG
 256             std::cout << "Exporting modPdA_3.vtk" << std::endl;
 257             WriteVTK("modPdA_3.vtk", modPdA);
 258
 259             std::cout << "Exporting modPdB_3.vtk" << std::endl;
 260             WriteVTK("modPdB_3.vtk", modPdB);
 261 #endif
 262
 263             start = clock::now();
 264
 265             ResolveOverlaps(modPdA, contsA, polyStripsA);
 266             ResolveOverlaps(modPdB, contsB, polyStripsB);
 267
 268             times.push_back(clock::now()-start);
 269
 270 #ifdef DEBUG
 271             std::cout << "Exporting modPdA_4.vtk" << std::endl;
 272             WriteVTK("modPdA_4.vtk", modPdA);
 273
 274             std::cout << "Exporting modPdB_4.vtk" << std::endl;
 275             WriteVTK("modPdB_4.vtk", modPdB);
 276 #endif
 277
 278             start = clock::now();
 279
 280             AddAdjacentPoints(modPdA, contsA, polyStripsA);
 281             AddAdjacentPoints(modPdB, contsB, polyStripsB);
 282
 283             times.push_back(clock::now()-start);
 284
 285 #ifdef DEBUG
 286             std::cout << "Exporting modPdA_5.vtk" << std::endl;
 287             WriteVTK("modPdA_5.vtk", modPdA);
 288
 289             std::cout << "Exporting modPdB_5.vtk" << std::endl;
 290             WriteVTK("modPdB_5.vtk", modPdB);
 291 #endif
 292
 293             start = clock::now();
 294
 295             DisjoinPolys(modPdA, polyStripsA);
 296             DisjoinPolys(modPdB, polyStripsB);
 297
 298             times.push_back(clock::now()-start);
 299
 300 #ifdef DEBUG
 301             std::cout << "Exporting modPdA_6.vtk" << std::endl;
 302             WriteVTK("modPdA_6.vtk", modPdA);
 303
 304             std::cout << "Exporting modPdB_6.vtk" << std::endl;
 305             WriteVTK("modPdB_6.vtk", modPdB);
 306 #endif
 307
 308             start = clock::now();
 309
 310             MergePoints(modPdA, polyStripsA);
 311             MergePoints(modPdB, polyStripsB);
 312
 313             times.push_back(clock::now()-start);
 314
 315 #ifdef DEBUG
 316             std::cout << "Exporting modPdA_7.vtk" << std::endl;
 317             WriteVTK("modPdA_7.vtk", modPdA);
 318
 319             std::cout << "Exporting modPdB_7.vtk" << std::endl;
 320             WriteVTK("modPdB_7.vtk", modPdB);
 321 #endif
 322
 323             timePdA = pdA->GetMTime();
 324             timePdB = pdB->GetMTime();
 325
 326         }
 327
 328         start = clock::now();
 329
 330         if (CombineRegions()) {
 331             vtkErrorMacro("Boolean operation failed.");
 332
 333             return 1;
 334         }
 335
 336         times.push_back(clock::now()-start);
 337
 338 // #ifdef DEBUG
 339         double sum = std::chrono::duration_cast<std::chrono::duration<double>>(std::accumulate(times.begin(), times.end(), clock::duration())).count();
 340
 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();
 344
 345             std::cout << "Time " << (itr-times.begin())
 346                 << ": " << time << "s (" << (time/sum*100) << "%)"
 347                 << std::endl;
 348         }
 349 // #endif
 350
 351     }
 352
 353     return 1;
 354
 355 }
 356
 357 void vtkPolyDataBooleanFilter::GetStripPoints (vtkPolyData *pd, vtkIdTypeArray *sources, PStrips &pStrips, IdsType &lines) {
 358
 359 #ifdef DEBUG
 360     std::cout << "GetStripPoints()" << std::endl;
 361 #endif
 362
 363     StripPtsType &pts = pStrips.pts;
 364     const IdsType &poly = pStrips.poly;
 365
 366     double a[3], b[3], sA[3], sB[3], u[3], v[3], w[3], n, t, d;
 367
 368     std::map<vtkIdType, vtkIdType> allPts;
 369
 370     vtkSmartPointer<vtkIdList> line = vtkSmartPointer<vtkIdList>::New();
 371
 372     for (auto lineId : lines) {
 373         contLines->GetCellPoints(lineId, line);
 374
 375         // std::cout << "? " << contsA->GetValue(lineId)
 376         //     << ", " << contsB->GetValue(lineId)
 377         //     << ", [" << line->GetId(0)
 378         //     << ", " << line->GetId(1)
 379         //     << "]"
 380         //     << std::endl;
 381
 382         allPts.emplace(line->GetId(0), sources->GetTypedComponent(lineId, 0));
 383         allPts.emplace(line->GetId(1), sources->GetTypedComponent(lineId, 1));
 384     }
 385
 386     decltype(allPts)::const_iterator itr;
 387
 388     for (itr = allPts.begin(); itr != allPts.end(); ++itr) {
 389         StripPt sp;
 390         sp.ind = itr->first;
 391
 392         // die koordinaten
 393         double pt[3];
 394         contLines->GetPoint(sp.ind, pt);
 395
 396         Cpy(sp.pt, pt, 3);
 397
 398         IdsType::const_iterator itrA, itrB;
 399
 400         for (itrA = poly.begin(); itrA != poly.end(); ++itrA) {
 401             itrB = itrA+1;
 402
 403             if (itrB == poly.end()) {
 404                 itrB = poly.begin();
 405             }
 406
 407             if (itr->second != NOTSET && *itrA != itr->second) {
 408                 continue;
 409             }
 410
 411             pd->GetPoint(*itrA, a);
 412             pd->GetPoint(*itrB, b);
 413
 414             vtkMath::Subtract(a, pt, sA);
 415             vtkMath::Subtract(b, pt, sB);
 416
 417             // richtungsvektor und länge der kante
 418
 419             vtkMath::Subtract(b, a, u);
 420             n = vtkMath::Norm(u);
 421
 422             // d und t zur kante
 423
 424             vtkMath::Subtract(pt, a, v);
 425             t = vtkMath::Dot(v, u)/(n*n);
 426
 427             vtkMath::Cross(v, u, w);
 428             d = vtkMath::Norm(w)/n;
 429
 430             if (d < 1e-5 && t > -1e-5 && t < 1+1e-5) {
 431                 sp.edge[0] = *itrA;
 432                 sp.edge[1] = *itrB;
 433
 434                 sp.t = std::min(1., std::max(0., t));
 435
 436                 if (vtkMath::Norm(sA) < 1e-5) {
 437                     Cpy(sp.captPt, a, 3);
 438                     sp.capt = Capt::A;
 439
 440                 } else if (vtkMath::Norm(sB) < 1e-5) {
 441                     Cpy(sp.captPt, b, 3);
 442                     sp.capt = Capt::B;
 443
 444                 } else {
 445                     // u ist nicht normiert
 446                     vtkMath::MultiplyScalar(u, t);
 447
 448                     double x[3];
 449                     vtkMath::Add(a, u, x);
 450
 451                     // projektion
 452                     Cpy(sp.captPt, x, 3);
 453
 454                     sp.capt = Capt::EDGE;
 455
 456                 }
 457             }
 458
 459             // std::cout << "? "
 460             //     << sp.ind
 461             //     << ", " << d
 462             //     << ", " << t
 463             //     << ", " << sp.capt
 464             //     << std::endl;
 465         }
 466
 467         if (itr->second != NOTSET && sp.edge[0] == NOTSET) {
 468             sp.catched = false;
 469         }
 470
 471         pts.emplace(sp.ind, std::move(sp));
 472
 473     }
 474
 475     StripPtsType::iterator itr2;
 476
 477     for (itr2 = pts.begin(); itr2 != pts.end(); ++itr2) {
 478         StripPt &sp = itr2->second;
 479
 480         if (sp.capt != Capt::NOT) {
 481             if (sp.capt == Capt::B) {
 482                 sp.t = 0;
 483
 484                 sp.edge[0] = sp.edge[1];
 485
 486                 auto itrA = std::find(poly.begin(), poly.end(), sp.edge[0]),
 487                     itrB = itrA+1;
 488
 489                 if (itrB == poly.end()) {
 490                     itrB = poly.begin();
 491                 }
 492
 493                 sp.edge[1] = *itrB;
 494
 495                 sp.capt = Capt::A;
 496
 497             }
 498
 499             // für den schnitt werden die eingerasteten koordinaten verwendet
 500
 501             Cpy(sp.cutPt, sp.captPt, 3);
 502         } else {
 503
 504             Cpy(sp.cutPt, sp.pt, 3);
 505         }
 506
 507         sp.history.emplace_back(sp.edge[0], sp.edge[1]);
 508
 509     }
 510
 511 #ifdef DEBUG
 512     std::cout << "pts: " << std::endl;
 513     for (itr2 = pts.begin(); itr2 != pts.end(); ++itr2) {
 514         std::cout << itr2->first << ": " << itr2->second << std::endl;
 515     }
 516 #endif
 517
 518 }
 519
 520 bool vtkPolyDataBooleanFilter::GetPolyStrips (vtkPolyData *pd, vtkIdTypeArray *conts, vtkIdTypeArray *sources, PolyStripsType &polyStrips) {
 521 #ifdef DEBUG
 522     std::cout << "GetPolyStrips()" << std::endl;
 523 #endif
 524
 525     polyStrips.clear();
 526
 527     std::map<vtkIdType, IdsType> polyLines;
 528
 529     for (vtkIdType i = 0; i < conts->GetNumberOfTuples(); i++) {
 530         vtkIdType poly = conts->GetValue(i);
 531
 532         /*if (poly != 1641) {
 533             continue;
 534         }*/
 535
 536         polyLines[poly].push_back(i);
 537     }
 538
 539     std::vector<std::reference_wrapper<StripPt>> notCatched;
 540
 541     std::map<vtkIdType, IdsType>::iterator itr;
 542
 543     for (itr = polyLines.begin(); itr != polyLines.end(); ++itr) {
 544
 545         IdsType &lines = itr->second;
 546         RemoveDuplicates(lines);
 547
 548         polyStrips.emplace(std::piecewise_construct,
 549             std::forward_as_tuple(itr->first),
 550             std::forward_as_tuple(pd, itr->first));
 551
 552         PStrips &pStrips = polyStrips.at(itr->first);
 553
 554         GetStripPoints(pd, sources, pStrips, lines);
 555
 556         for (auto &sp : pStrips.pts) {
 557             sp.second.polyId = itr->first;
 558
 559             if (!sp.second.catched) {
 560                 notCatched.push_back(sp.second);
 561             }
 562         }
 563
 564     }
 565
 566     auto Next = [](const IdsType &ids, vtkIdType id) -> vtkIdType {
 567         IdsType::const_iterator itr;
 568
 569         itr = std::find(ids.begin(), ids.end(), id);
 570
 571         if (++itr == ids.end()) {
 572             itr = ids.begin();
 573         }
 574
 575         return *itr;
 576     };
 577
 578     for (StripPt &sp : notCatched) {
 579         for (itr = polyLines.begin(); itr != polyLines.end(); ++itr) {
 580             const PStrips &pStrips = polyStrips.at(itr->first);
 581
 582             try {
 583                 const StripPt &corr = pStrips.pts.at(sp.ind);
 584
 585                 if (&corr != &sp) {
 586                     if (corr.capt == Capt::A) {
 587                         sp.capt = Capt::A;
 588                         sp.edge[0] = corr.edge[0];
 589                         sp.edge[1] = Next(polyStrips.at(sp.polyId).poly, sp.edge[0]);
 590
 591                         sp.t = 0;
 592
 593                         Cpy(sp.captPt, corr.captPt, 3);
 594                         Cpy(sp.cutPt, sp.captPt, 3);
 595
 596                         sp.history.emplace_back(sp.edge[0], sp.edge[1]);
 597
 598                         sp.catched = true;
 599
 600                     }
 601                 }
 602             } catch (...) {}
 603
 604         }
 605
 606         if (!sp.catched) {
 607             std::cout << sp << std::endl;
 608         }
 609
 610         assert(sp.catched);
 611
 612     }
 613
 614     vtkSmartPointer<vtkIdList> cells = vtkSmartPointer<vtkIdList>::New(),
 615         line = vtkSmartPointer<vtkIdList>::New();
 616
 617     vtkIdType i, numCells;
 618
 619     StripPtsType::const_iterator itr2;
 620
 621     for (itr = polyLines.begin(); itr != polyLines.end(); ++itr) {
 622         PStrips &pStrips = polyStrips.at(itr->first);
 623
 624         const IdsType &lines = itr->second;
 625         const StripPtsType &pts = pStrips.pts;
 626
 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
 629
 630             const StripPt &pt = itr2->second;
 631
 632             if (pt.capt == Capt::NOT) {
 633                 contLines->GetPointCells(pt.ind, cells);
 634
 635                 numCells = cells->GetNumberOfIds();
 636
 637                 std::set<vtkIdType> ends;
 638
 639                 for (i = 0; i < numCells; i++) {
 640                     contLines->GetCellPoints(cells->GetId(i), line);
 641
 642                     ends.insert(pt.ind == line->GetId(0) ? line->GetId(1) : line->GetId(0));
 643                 }
 644
 645                 if (ends.size() > 2) {
 646                     return true;
 647                 }
 648             }
 649         }
 650
 651         // zusammensetzen
 652
 653         std::deque<Pair> _lines;
 654
 655         vtkIdList *linePts = vtkIdList::New();
 656
 657         for (auto &i : lines) {
 658             contLines->GetCellPoints(i, linePts);
 659             _lines.emplace_back(linePts->GetId(0), linePts->GetId(1));
 660         }
 661
 662         linePts->Delete();
 663
 664         decltype(_lines)::iterator _itr;
 665
 666         auto FindRight = [&pts, &_lines, &_itr](StripType &strip, const std::size_t &id) -> bool {
 667             auto &right = strip.back();
 668
 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);
 673
 674                         _lines.erase(_itr);
 675                         return true;
 676                     } else if (_itr->g == right.ind) {
 677                         strip.emplace_back(_itr->f, id);
 678
 679                         _lines.erase(_itr);
 680                         return true;
 681                     }
 682                 }
 683             }
 684
 685             return false;
 686         };
 687
 688         auto FindLeft = [&pts, &_lines, &_itr](StripType &strip, const std::size_t &id) -> bool {
 689             auto &left = strip.front();
 690
 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);
 695
 696                         _lines.erase(_itr);
 697                         return true;
 698                     } else if (_itr->g == left.ind) {
 699                         strip.emplace_front(_itr->f, id);
 700
 701                         _lines.erase(_itr);
 702                         return true;
 703                     }
 704                 }
 705             }
 706
 707             return false;
 708         };
 709
 710         std::size_t stripId {0};
 711
 712         while (!_lines.empty()) {
 713             auto &last = _lines.back();
 714
 715             StripType strip {{last.f, stripId}, {last.g, stripId}};
 716             _lines.pop_back();
 717
 718             while (FindRight(strip, stripId)) {}
 719             while (FindLeft(strip, stripId)) {}
 720
 721             pStrips.strips.push_back(std::move(strip));
 722
 723             stripId++;
 724
 725         }
 726
 727         CompleteStrips(pStrips);
 728
 729     }
 730
 731     // sucht nach schnitten zw. den strips
 732
 733     {
 734
 735         PolyStripsType::const_iterator itr;
 736         StripType::const_iterator itr2;
 737
 738         for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) {
 739             const PStrips &pStrips = itr->second;
 740
 741             const StripsType &strips = pStrips.strips;
 742             const StripPtsType &pts = pStrips.pts;
 743             const Base &base = pStrips.base;
 744
 745             vtkSmartPointer<vtkPoints> treePts = vtkSmartPointer<vtkPoints>::New();
 746
 747             vtkSmartPointer<vtkPolyData> treePd = vtkSmartPointer<vtkPolyData>::New();
 748             treePd->Allocate(1);
 749
 750             std::map<vtkIdType, vtkIdType> ptIds;
 751
 752             double pt[2];
 753
 754             for (const auto &p : pts) {
 755                 Transform(p.second.pt, pt, base);
 756
 757                 ptIds.emplace(p.first, treePts->InsertNextPoint(pt[0], pt[1], 0));
 758             }
 759
 760             for (const StripType &strip : strips) {
 761                 for (itr2 = strip.begin(); itr2 != strip.end()-1; ++itr2) {
 762
 763                     vtkIdList *line = vtkIdList::New();
 764                     line->InsertNextId(ptIds[itr2->ind]);
 765                     line->InsertNextId(ptIds[(itr2+1)->ind]);
 766
 767                     treePd->InsertNextCell(VTK_LINE, line);
 768
 769                     line->Delete();
 770                 }
 771             }
 772
 773             treePd->SetPoints(treePts);
 774
 775             vtkSmartPointer<vtkModifiedBSPTree> tree = vtkSmartPointer<vtkModifiedBSPTree>::New();
 776             tree->SetDataSet(treePd);
 777             tree->BuildLocator();
 778
 779             vtkIdType numA, numB;
 780             const vtkIdType *lineA, *lineB;
 781
 782             auto lineItr = vtk::TakeSmartPointer(treePd->GetLines()->NewIterator());
 783
 784             for (lineItr->GoToFirstCell(); !lineItr->IsDoneWithTraversal(); lineItr->GoToNextCell()) {
 785                 lineItr->GetCurrentCell(numA, lineA);
 786
 787                 double ptA[3], ptB[3];
 788
 789                 treePts->GetPoint(lineA[0], ptA);
 790                 treePts->GetPoint(lineA[1], ptB);
 791
 792                 vtkSmartPointer<vtkIdList> lineIds = vtkSmartPointer<vtkIdList>::New();
 793
 794                 tree->IntersectWithLine(ptA, ptB, 1e-5, nullptr, lineIds);
 795
 796                 for (vtkIdType i = 0; i < lineIds->GetNumberOfIds(); i++) {
 797                     treePd->GetCellPoints(lineIds->GetId(i), numB, lineB);
 798
 799                     if (lineB[0] != lineA[0] && lineB[1] != lineA[0] && lineB[0] != lineA[1] && lineB[1] != lineA[1]) {
 800                         // schnitt gefunden
 801
 802                         return true;
 803                     }
 804                 }
 805             }
 806
 807         }
 808
 809     }
 810
 811     return false;
 812
 813 }
 814
 815 void vtkPolyDataBooleanFilter::RemoveDuplicates (IdsType &lines) {
 816
 817     IdsType unique;
 818
 819     // die indexe der enden auf übereinstimmung prüfen
 820
 821     vtkIdList *linePtsA = vtkIdList::New();
 822     vtkIdList *linePtsB = vtkIdList::New();
 823
 824     std::size_t i, j, numLines = lines.size();
 825
 826     for (i = 0; i < numLines-1; i++) {
 827         j = i+1;
 828
 829         contLines->GetCellPoints(lines[i], linePtsA);
 830
 831         while (j < lines.size()) {
 832             contLines->GetCellPoints(lines[j], linePtsB);
 833
 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))) {
 836                 // stimmen überein
 837                 break;
 838             }
 839
 840             j++;
 841         }
 842
 843         if (j == numLines) {
 844             // keine vorzeitige unterbrechung der while-schleife
 845             unique.push_back(lines[i]);
 846         }
 847     }
 848
 849     unique.push_back(lines.back());
 850
 851     linePtsA->Delete();
 852     linePtsB->Delete();
 853
 854     lines.swap(unique);
 855
 856 }
 857
 858 void vtkPolyDataBooleanFilter::CompleteStrips (PStrips &pStrips) {
 859     StripsType::iterator itr;
 860
 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];
 864
 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());
 869
 870             } else if (end.capt == Capt::NOT) {
 871                 StripType s(itr->rbegin()+1, itr->rend());
 872                 itr->insert(itr->end(), s.begin(), s.end());
 873
 874             }
 875         }
 876     }
 877 }
 878
 879 bool vtkPolyDataBooleanFilter::HasArea (const StripType &strip) const {
 880     bool area = true;
 881
 882     std::size_t i, n = strip.size();
 883
 884     if (n%2 == 1) {
 885         for (i = 0; i < (n-1)/2; i++) {
 886             area = strip[i].ind != strip[n-i-1].ind;
 887         }
 888     }
 889
 890     return area;
 891 }
 892
 893 void ComputeNormal (const StripPtsType &pts, const RefsType &poly, double *n) {
 894     n[0] = 0; n[1] = 0; n[2] = 0;
 895
 896     RefsType::const_iterator itrA, itrB;
 897
 898     for (itrA = poly.begin(); itrA != poly.end(); ++itrA) {
 899         itrB = itrA+1;
 900
 901         if (itrB == poly.end()) {
 902             itrB = poly.begin();
 903         }
 904
 905         const StripPtR &spA = *itrA,
 906             &spB = *itrB;
 907
 908         auto pA = pts.find(spA.ind);
 909         auto pB = pts.find(spB.ind);
 910
 911         const double *ptA = pA->second.cutPt,
 912             *ptB = pB->second.cutPt;
 913
 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]);
 917     }
 918
 919     vtkMath::Normalize(n);
 920 }
 921
 922 bool vtkPolyDataBooleanFilter::CutCells (vtkPolyData *pd, PolyStripsType &polyStrips) {
 923 #ifdef DEBUG
 924     std::cout << "CutCells()" << std::endl;
 925 #endif
 926
 927     vtkPoints *pdPts = pd->GetPoints();
 928
 929     vtkIdTypeArray *origCellIds = vtkIdTypeArray::SafeDownCast(pd->GetCellData()->GetScalars("OrigCellIds"));
 930
 931     PolyStripsType::iterator itr;
 932
 933     for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) {
 934
 935         vtkIdType polyInd = itr->first;
 936         PStrips &pStrips = itr->second;
 937
 938         /*if (polyInd != 1641) {
 939             continue;
 940         }*/
 941
 942         StripsType &strips = pStrips.strips;
 943         StripPtsType &pts = pStrips.pts;
 944
 945         IdsType &poly = pStrips.poly;
 946
 947         vtkIdType origId = origCellIds->GetValue(polyInd);
 948
 949 #ifdef DEBUG
 950         std::cout << "polyInd " << polyInd << ", poly [";
 951         for (auto &p : poly) {
 952             std::cout << p << ", ";
 953         }
 954         std::cout << "]" << std::endl;
 955 #endif
 956
 957         std::size_t numPts = poly.size();
 958
 959         std::map<vtkIdType, RefsType> edges;
 960
 961         StripsType::iterator itr2;
 962         StripType::iterator itr3;
 963
 964         // holes sichern
 965         StripsType holes;
 966
 967         auto fct = [&](const StripType &s) {
 968             return pts[s.front().ind].capt == Capt::NOT && pts[s.back().ind].capt == Capt::NOT;
 969         };
 970
 971         std::copy_if(strips.begin(), strips.end(), std::back_inserter(holes), fct);
 972
 973         strips.erase(std::remove_if(strips.begin(), strips.end(), fct), strips.end());
 974
 975         for (itr2 = strips.begin(); itr2 != strips.end(); ++itr2) {
 976             StripType &strip = *itr2;
 977
 978 #ifdef DEBUG
 979             std::cout << (itr2-strips.begin()) << ". strip [";
 980             for (auto &s : strip) {
 981                 std::cout << s.ind << ", ";
 982             }
 983             std::cout << "]" << std::endl;
 984 #endif
 985
 986             // init
 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) {
 990
 991                 std::reverse(strip.begin(), strip.end());
 992             }
 993
 994             StripPt &start = pts[strip.front().ind],
 995                 &end = pts[strip.back().ind];
 996
 997             strip.front().side = Side::START;
 998             strip.back().side = Side::END;
 999
1000             strip.front().ref = start.edge[0];
1001             strip.back().ref = end.edge[0];
1002
1003             // nachfolgend könnte man dann anfang und ende weglassen
1004
1005             for (itr3 = strip.begin(); itr3 != strip.end(); ++itr3) {
1006                 StripPt &sp = pts[itr3->ind];
1007
1008                 itr3->desc[0] = pdPts->InsertNextPoint(sp.cutPt);
1009                 itr3->desc[1] = pdPts->InsertNextPoint(sp.cutPt);
1010
1011 #ifdef DEBUG
1012                 std::cout << sp << " => " << *itr3 << std::endl;
1013 #endif
1014
1015             }
1016
1017             // ordnet zu
1018
1019             edges[start.edge[0]].push_back(std::ref(strip.front()));
1020             edges[end.edge[0]].push_back(std::ref(strip.back()));
1021         }
1022
1023         // sortiert die punkte auf den kanten
1024
1025         std::map<vtkIdType, RefsType>::iterator itr4;
1026
1027         IdsType::iterator itr5;
1028         StripType::reverse_iterator itr6;
1029
1030         RefsType::iterator itr7;
1031
1032         for (itr4 = edges.begin(); itr4 != edges.end(); ++itr4) {
1033             RefsType &edge = itr4->second;
1034
1035 #ifdef DEBUG
1036             std::cout << "edge (" << itr4->first << ", _)" << std::endl;
1037 #endif
1038
1039             std::sort(edge.begin(), edge.end(), [&](const StripPtR &a, const StripPtR &b) {
1040                 const StripPt &a_ = pts[a.ind],
1041                     &b_ = pts[b.ind];
1042
1043 #ifdef DEBUG
1044                 // std::cout << "a_: " << a_ << " -> strip " << a.strip << std::endl;
1045                 // std::cout << "b_: " << b_ << " -> strip " << b.strip << std::endl;
1046 #endif
1047
1048                 if (a_.ind == b_.ind) {
1049                     // strips beginnen im gleichen punkt
1050
1051                     if (a.strip != b.strip) {
1052                         // gehören nicht dem gleichen strip an
1053
1054                         StripType &stripA = strips[a.strip],
1055                             &stripB = strips[b.strip];
1056
1057                         // andere enden ermitteln
1058
1059                         const StripPtR &eA = (&a == &(stripA.front())) ? stripA.back() : stripA.front(),
1060                             &eB = (&b == &(stripB.front())) ? stripB.back() : stripB.front();
1061
1062                         const StripPt &eA_ = pts[eA.ind],
1063                             &eB_ = pts[eB.ind];
1064
1065 #ifdef DEBUG
1066                         // std::cout << "eA_: " << eA_ << std::endl;
1067                         // std::cout << "eB_: " << eA_ << std::endl;
1068 #endif
1069
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();
1074
1075                             double dA = static_cast<double>(Mod(iA-i, numPts))+eA_.t,
1076                                 dB = static_cast<double>(Mod(iB-i, numPts))+eB_.t;
1077
1078                             if (i == iA && a_.t > eA_.t) {
1079                                dA += static_cast<double>(numPts);
1080                             }
1081
1082                             if (i == iB && b_.t > eB_.t) {
1083                                dB += static_cast<double>(numPts);
1084                             }
1085
1086 #ifdef DEBUG
1087                             // std::cout << "dA=" << dA << ", dB=" << dB << std::endl;
1088 #endif
1089
1090                             return dB < dA;
1091                         } else {
1092                             RefsType poly_;
1093
1094                             if (a.side == Side::START) {
1095                                 poly_.insert(poly_.end(), stripA.begin(), stripA.end());
1096                             } else {
1097                                 poly_.insert(poly_.end(), stripA.rbegin(), stripA.rend());
1098                             }
1099
1100                             if (b.side == Side::START) {
1101                                 poly_.insert(poly_.end(), stripB.rbegin()+1, stripB.rend()-1);
1102                             } else {
1103                                 poly_.insert(poly_.end(), stripB.begin()+1, stripB.end()-1);
1104                             }
1105
1106                             double n[3];
1107                             ComputeNormal(pts, poly_, n);
1108
1109                             double ang = vtkMath::Dot(pStrips.n, n);
1110
1111 #ifdef DEBUG
1112                             // std::cout << "ang=" << ang*180/PI << std::endl;
1113 #endif
1114
1115                             return ang < .999999;
1116
1117                         }
1118                     } else {
1119                         // gleicher strip
1120
1121                         StripType &strip = strips[a.strip];
1122
1123                         if (HasArea(strip)) {
1124                             RefsType poly_(strip.begin(), strip.end()-1);
1125
1126                             double n[3];
1127                             ComputeNormal(pts, poly_, n);
1128
1129                             double ang = vtkMath::Dot(pStrips.n, n);
1130
1131                             if (ang > .999999) {
1132                                 std::reverse(strip.begin(), strip.end());
1133                                 return true;
1134                             } else {
1135                                 return false;
1136                             }
1137
1138                         } else {
1139                             // reihenfolge von a und b bereits richtig
1140                             return false;
1141                         }
1142                     }
1143
1144                 } else {
1145                     return a_.t < b_.t;
1146                 }
1147             });
1148
1149 #ifdef DEBUG
1150             for (auto& p : edge) {
1151                 std::cout << p << std::endl;
1152             }
1153 #endif
1154
1155         }
1156
1157         // baut die strips ein
1158
1159         std::deque<IdsType> polys;
1160         polys.push_back(pStrips.poly);
1161
1162         for (itr2 = strips.begin(); itr2 != strips.end(); ++itr2) {
1163             StripType &strip = *itr2;
1164
1165             StripPtR &start = strip.front(),
1166                 &end = strip.back();
1167
1168 #ifdef DEBUG
1169             std::cout << "strip " << start.strip
1170                 << " , refs (" << start.ref << ", " << end.ref << ")"
1171                 << std::endl;
1172 #endif
1173
1174             std::size_t cycle = 0;
1175
1176             while (true) {
1177                 if (cycle == polys.size()) {
1178                     break;
1179                 }
1180
1181                 IdsType next = polys.front();
1182                 polys.pop_front();
1183
1184                 std::vector<IdsType> newPolys(2);
1185
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);
1190
1191 #ifdef DEBUG
1192                             std::cout << "adding " << *itr5 << " to 0" << std::endl;
1193 #endif
1194
1195                             if (*itr5 == start.ref) {
1196                                 for (itr3 = strip.begin(); itr3 != strip.end(); ++itr3) {
1197                                     newPolys[0].push_back(itr3->desc[0]);
1198
1199 #ifdef DEBUG
1200                                     std::cout << "adding " << itr3->desc[0] << " to 0" << std::endl;
1201 #endif
1202
1203                                 }
1204                             }
1205                         }
1206
1207                         // strip selbst ist ein polygon
1208
1209                         for (itr6 = strip.rbegin(); itr6 != strip.rend(); ++itr6) {
1210                             newPolys[1].push_back(itr6->desc[1]);
1211
1212 #ifdef DEBUG
1213                             std::cout << "adding " << itr6->desc[1] << " to 1" << std::endl;
1214 #endif
1215
1216                         }
1217
1218                     } else {
1219                         std::size_t curr = 0;
1220
1221                         for (itr5 = next.begin(); itr5 != next.end(); ++itr5) {
1222                             IdsType &poly = newPolys[curr];
1223
1224                             poly.push_back(*itr5);
1225
1226 #ifdef DEBUG
1227                             std::cout << "adding " << *itr5 << " to " << curr << std::endl;
1228 #endif
1229
1230                             if (*itr5 == start.ref) {
1231                                 for (itr3 = strip.begin(); itr3 != strip.end(); ++itr3) {
1232                                     poly.push_back(itr3->desc[0]);
1233
1234 #ifdef DEBUG
1235                                     std::cout << "adding " << itr3->desc[0] << " to " << curr << std::endl;
1236 #endif
1237
1238                                 }
1239
1240                                 curr = curr == 0 ? 1 : 0;
1241
1242                             } else if (*itr5 == end.ref) {
1243                                 for (itr6 = strip.rbegin(); itr6 != strip.rend(); ++itr6) {
1244                                     poly.push_back(itr6->desc[1]);
1245
1246 #ifdef DEBUG
1247                                     std::cout << "adding " << itr6->desc[1] << " to " << curr << std::endl;
1248 #endif
1249
1250                                 }
1251
1252                                 curr = curr == 0 ? 1 : 0;
1253                             }
1254
1255                         }
1256                     }
1257                 }
1258
1259                 if (!newPolys[1].empty()) {
1260                     // refs aktualisieren
1261
1262                     for (itr4 = edges.begin(); itr4 != edges.end(); ++itr4) {
1263                         RefsType &edge = itr4->second;
1264
1265                         RefsType::iterator itrA;
1266
1267                         for (itrA = edge.begin()+1; itrA != edge.end(); ++itrA) {
1268                             StripPtR &sp = *itrA;
1269
1270                             if (sp.strip > start.strip) {
1271 #ifdef DEBUG
1272                                 std::cout << "sp: ind " << sp.ind << ", strip " << sp.strip << std::endl;
1273 #endif
1274
1275                                 RefsType::const_reverse_iterator itrB(itrA);
1276
1277                                 std::shared_ptr<StripPtR> _p;
1278
1279                                 for (; itrB != edge.rend(); ++itrB) {
1280                                     const StripPtR &p = *itrB;
1281
1282                                     if (p.strip != sp.strip) {
1283                                         if (p.strip <= start.strip) {
1284 #ifdef DEBUG
1285                                             std::cout << "ref " << sp.ref;
1286 #endif
1287
1288                                             if (p.side == Side::END) {
1289                                                 sp.ref = p.desc[0];
1290                                             } else {
1291                                                 sp.ref = p.desc[1];
1292                                             }
1293
1294 #ifdef DEBUG
1295                                             std::cout << " -> " << sp.ref << " (from strip " << p.strip << ", ind " << p.ind << ")" << std::endl;
1296 #endif
1297
1298                                             _p = std::make_shared<StripPtR>(p);
1299
1300                                             break;
1301
1302                                         }
1303                                     } else {
1304 #ifdef DEBUG
1305                                         std::cout << "~1 ref " << sp.ref << " -> " << p.ref << " (from strip " << p.strip << ", ind " << p.ind << ")" << std::endl;
1306 #endif
1307
1308                                         sp.ref = p.ref;
1309                                         break;
1310                                     }
1311                                 }
1312
1313                                 RefsType::const_iterator itrC(itrA);
1314
1315                                 ++itrC;
1316
1317                                 for (; itrC != edge.end(); ++itrC) {
1318                                     const StripPtR &p = *itrC;
1319
1320                                     if (p.ind != sp.ind) {
1321                                         break;
1322                                     }
1323
1324                                     if (p.strip <= start.strip) {
1325                                         if (_p && p.ind == _p->ind && p.strip < _p->strip) {
1326                                             break;
1327                                         }
1328
1329 #ifdef DEBUG
1330                                         std::cout << "~2 ref " << sp.ref;
1331 #endif
1332
1333                                         if (p.side == Side::START) {
1334                                             sp.ref = p.desc[0];
1335                                         } else {
1336                                             sp.ref = p.desc[1];
1337                                         }
1338
1339 #ifdef DEBUG
1340                                         std::cout << " -> " << sp.ref << " (from strip " << p.strip << ", ind " << p.ind << ")" << std::endl;
1341 #endif
1342
1343                                         break;
1344                                     }
1345
1346                                 }
1347                             }
1348                         }
1349
1350                         // erstellt die history
1351
1352                         auto _s = std::find_if(edge.begin(), edge.end(), [&start](const StripPtR &r) {
1353                             return &start == &r;
1354                         });
1355
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;
1361
1362                                 if (&sp != &start) {
1363                                     if (_sp.t < pts[start.ind].t) {
1364                                         history.push_back({history.back().f, start.desc[0]});
1365                                     } else {
1366                                         history.push_back({start.desc[1], history.back().g});
1367                                     }
1368
1369                                 }
1370
1371                             }
1372                         }
1373
1374                         auto _e = std::find_if(edge.begin(), edge.end(), [&end](const StripPtR &r) {
1375                             return &end == &r;
1376                         });
1377
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;
1383
1384                                 if (&sp != &end) {
1385                                     if (_sp.t < pts[end.ind].t) {
1386                                         history.push_back({history.back().f, end.desc[1]});
1387                                     } else {
1388                                         history.push_back({end.desc[0], history.back().g});
1389                                     }
1390
1391                                 }
1392
1393                             }
1394                         }
1395
1396                         // sonderfall
1397                         if (edge.size() > 1) {
1398                             StripPtR &a = edge.front(),
1399                                 &b = *(edge.begin()+1);
1400
1401                             if (a.ind == b.ind
1402                                 && b.strip == start.strip
1403                                 && pts[a.ind].capt == Capt::A) { // sollte weg
1404
1405 #ifdef DEBUG
1406                                 std::cout << "~3 ref " << a.ref;
1407 #endif
1408
1409                                 if (b.side == Side::START) {
1410                                     a.ref = b.desc[0];
1411                                 } else {
1412                                     a.ref = b.desc[1];
1413                                 }
1414
1415 #ifdef DEBUG
1416                                 std::cout << " -> " << a.ref << " (from strip " << b.strip << ", ind " << b.ind << ")" << std::endl;
1417 #endif
1418
1419                             }
1420                         }
1421
1422                     }
1423
1424                     // doppelte punkte entfernen
1425
1426                     double pt[3];
1427
1428                     for (auto &newPoly : newPolys) {
1429                         IdsType _newPoly;
1430
1431                         std::map<vtkIdType, Point3d> _pts;
1432
1433                         for (vtkIdType id : newPoly) {
1434                             pd->GetPoint(id, pt);
1435
1436                             _pts.emplace(std::piecewise_construct,
1437                                 std::forward_as_tuple(id),
1438                                 std::forward_as_tuple(pt[0], pt[1], pt[2]));
1439                         }
1440
1441                         IdsType::const_iterator itrA, itrB;
1442
1443                         for (itrA = newPoly.begin(); itrA != newPoly.end(); ++itrA) {
1444                             itrB = itrA+1;
1445                             if (itrB == newPoly.end()) {
1446                                 itrB = newPoly.begin();
1447                             }
1448
1449                             auto _a = _pts.find(*itrA);
1450                             auto _b = _pts.find(*itrB);
1451
1452                             if (_a->second == _b->second) {
1453                                 // doppelt
1454 #ifdef DEBUG
1455                                 std::cout << "removing " << *itrA << std::endl;
1456 #endif
1457                             } else {
1458                                 _newPoly.push_back(*itrA);
1459                             }
1460                         }
1461
1462                         newPoly.swap(_newPoly);
1463
1464                     }
1465
1466                     // prüft, ob die erstellten polys gültig sind
1467
1468                     if (newPolys[0].size() > 2) {
1469                         polys.push_back(newPolys[0]);
1470                     }
1471
1472                     if (HasArea(strip) && newPolys[1].size() > 2) {
1473                         polys.push_back(newPolys[1]);
1474                     }
1475
1476                     cycle = 0;
1477
1478                     break;
1479
1480                 } else {
1481                     polys.push_back(next);
1482
1483                     cycle++;
1484                 }
1485
1486             }
1487
1488         }
1489
1490         // erzeugte polys hinzufügen
1491
1492         IdsType descIds;
1493         descIds.reserve(polys.size());
1494
1495         for (auto &p : polys) {
1496             vtkIdList *cell = vtkIdList::New();
1497
1498             for (vtkIdType &id : p) {
1499                 cell->InsertNextId(id);
1500             }
1501
1502             descIds.emplace_back(pd->InsertNextCell(VTK_POLYGON, cell));
1503             origCellIds->InsertNextValue(origId);
1504
1505             cell->Delete();
1506         }
1507
1508         pd->DeleteCell(polyInd);
1509
1510         // holes einbauen
1511         if (!holes.empty()) {
1512             try {
1513                 Merger(pd, pStrips, holes, descIds, origId).Run();
1514             } catch (const std::runtime_error &e) {
1515                 return true;
1516             }
1517         }
1518
1519     }
1520
1521     pd->RemoveDeletedCells();
1522     pd->BuildCells();
1523
1524     return false;
1525
1526 }
1527
1528 void vtkPolyDataBooleanFilter::RestoreOrigPoints (vtkPolyData *pd, PolyStripsType &polyStrips) {
1529
1530 #ifdef DEBUG
1531     std::cout << "RestoreOrigPoints()" << std::endl;
1532 #endif
1533
1534     pd->BuildLinks();
1535
1536     vtkKdTreePointLocator *loc = vtkKdTreePointLocator::New();
1537     loc->SetDataSet(pd);
1538     loc->BuildLocator();
1539
1540     struct Cmp {
1541         bool operator() (const StripPt &l, const StripPt &r) const {
1542             return l.ind < r.ind;
1543         }
1544     };
1545
1546     PolyStripsType::const_iterator itr;
1547     StripPtsType::const_iterator itr2;
1548
1549     std::set<std::reference_wrapper<const StripPt>, Cmp> ends;
1550
1551     for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) {
1552         const PStrips &pStrips = itr->second;
1553
1554         for (itr2 = pStrips.pts.begin(); itr2 != pStrips.pts.end(); ++itr2) {
1555             const StripPt &sp = itr2->second;
1556
1557             if (sp.capt != Capt::NOT) {
1558                 ends.emplace(sp);
1559             }
1560
1561         }
1562     }
1563
1564     vtkIdList *pts = vtkIdList::New();
1565
1566     vtkIdType i, numPts;
1567
1568     for (const StripPt &sp : ends) {
1569         FindPoints(loc, sp.cutPt, pts);
1570         numPts = pts->GetNumberOfIds();
1571
1572         for (i = 0; i < numPts; i++) {
1573             pd->GetPoints()->SetPoint(pts->GetId(i), sp.pt);
1574         }
1575     }
1576
1577     pts->Delete();
1578
1579     loc->FreeSearchStructure();
1580     loc->Delete();
1581
1582 }
1583
1584 void vtkPolyDataBooleanFilter::DisjoinPolys (vtkPolyData *pd, PolyStripsType &polyStrips) {
1585
1586 #ifdef DEBUG
1587     std::cout << "DisjoinPolys()" << std::endl;
1588 #endif
1589
1590     pd->BuildLinks();
1591
1592     vtkKdTreePointLocator *loc = vtkKdTreePointLocator::New();
1593     loc->SetDataSet(pd);
1594     loc->BuildLocator();
1595
1596     struct Cmp {
1597         bool operator() (const StripPt &l, const StripPt &r) const {
1598             return l.ind < r.ind;
1599         }
1600     };
1601
1602     PolyStripsType::const_iterator itr;
1603     StripPtsType::const_iterator itr2;
1604
1605     std::set<std::reference_wrapper<const StripPt>, Cmp> ends;
1606
1607     for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) {
1608         const PStrips &pStrips = itr->second;
1609
1610         for (itr2 = pStrips.pts.begin(); itr2 != pStrips.pts.end(); ++itr2) {
1611             const StripPt &sp = itr2->second;
1612
1613             if (sp.capt == Capt::A) {
1614                 ends.emplace(sp);
1615             }
1616         }
1617     }
1618
1619     vtkIdList *pts = vtkIdList::New();
1620     vtkIdList *cells = vtkIdList::New();
1621
1622     vtkIdType i, j, numPts, numCells;
1623
1624     for (const StripPt &sp : ends) {
1625         FindPoints(loc, sp.pt, pts);
1626         numPts = pts->GetNumberOfIds();
1627
1628         for (i = 0; i < numPts; i++) {
1629             pd->GetPointCells(pts->GetId(i), cells);
1630             numCells = cells->GetNumberOfIds();
1631
1632             if (numCells > 1) {
1633                 for (j = 0; j < numCells; j++) {
1634                     pd->ReplaceCellPoint(cells->GetId(j), pts->GetId(i), pd->GetPoints()->InsertNextPoint(sp.pt));
1635                 }
1636             }
1637         }
1638     }
1639
1640     cells->Delete();
1641     pts->Delete();
1642
1643     loc->FreeSearchStructure();
1644     loc->Delete();
1645
1646 }
1647
1648 void vtkPolyDataBooleanFilter::ResolveOverlaps (vtkPolyData *pd, vtkIdTypeArray *conts, PolyStripsType &polyStrips) {
1649
1650 #ifdef DEBUG
1651     std::cout << "ResolveOverlaps()" << std::endl;
1652 #endif
1653
1654     pd->BuildCells();
1655     pd->BuildLinks();
1656
1657     contLines->BuildLinks();
1658
1659     vtkKdTreePointLocator *loc = vtkKdTreePointLocator::New();
1660     loc->SetDataSet(pd);
1661     loc->BuildLocator();
1662
1663     struct Cmp {
1664         bool operator() (const StripPt &l, const StripPt &r) const {
1665             return l.ind < r.ind;
1666         }
1667     };
1668
1669     PolyStripsType::const_iterator itr;
1670     StripPtsType::const_iterator itr2;
1671
1672     std::set<std::reference_wrapper<const StripPt>, Cmp> ends;
1673
1674     for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) {
1675         const PStrips &pStrips = itr->second;
1676
1677         for (itr2 = pStrips.pts.begin(); itr2 != pStrips.pts.end(); ++itr2) {
1678             const StripPt &sp = itr2->second;
1679
1680             if (sp.capt == Capt::EDGE) {
1681                 ends.emplace(sp);
1682             }
1683         }
1684     }
1685
1686     vtkIdList *links = vtkIdList::New(),
1687         *ptsA = vtkIdList::New(),
1688         *ptsB = vtkIdList::New(),
1689         *cells = vtkIdList::New(),
1690         *poly = vtkIdList::New();
1691
1692     double ptA[3], ptB[3];
1693
1694     vtkIdType i, j, k, numPtsA, numPtsB, numCells, numPts, idA, idB;
1695
1696     std::map<Pair, IdsType> sharedPts;
1697
1698     for (const StripPt &sp : ends) {
1699         contLines->GetPointCells(sp.ind, links);
1700
1701         if (links->GetNumberOfIds() == 2
1702             && conts->GetValue(links->GetId(0)) != conts->GetValue(links->GetId(1))) {
1703
1704             std::vector<Pair> history(sp.history.rbegin(), sp.history.rend());
1705
1706             for (const Pair &p : history) {
1707                 pd->GetPoint(p.f, ptA);
1708                 pd->GetPoint(p.g, ptB);
1709
1710                 FindPoints(loc, ptA, ptsA);
1711                 FindPoints(loc, ptB, ptsB);
1712
1713                 numPtsA = ptsA->GetNumberOfIds();
1714                 numPtsB = ptsB->GetNumberOfIds();
1715
1716                 std::vector<Pair> cellsA, cellsB;
1717
1718                 for (i = 0; i < numPtsA; i++) {
1719                     pd->GetPointCells(ptsA->GetId(i), cells);
1720                     numCells = cells->GetNumberOfIds();
1721
1722                     for (j = 0; j < numCells; j++) {
1723                         cellsA.emplace_back(ptsA->GetId(i), cells->GetId(j));
1724                     }
1725                 }
1726
1727                 for (i = 0; i < numPtsB; i++) {
1728                     pd->GetPointCells(ptsB->GetId(i), cells);
1729                     numCells = cells->GetNumberOfIds();
1730
1731                     for (j = 0; j < numCells; j++) {
1732                         cellsB.emplace_back(ptsB->GetId(i), cells->GetId(j));
1733                     }
1734                 }
1735
1736                 for (auto &cellA : cellsA) {
1737                     for (auto &cellB : cellsB) {
1738                         if (cellA.g == cellB.g) {
1739                             // kante/poly existiert noch
1740
1741                             pd->GetCellPoints(cellA.g, poly);
1742                             numPts = poly->GetNumberOfIds();
1743
1744                             for (k = 0; k < numPts; k++) {
1745                                 idA = poly->GetId(k);
1746                                 idB = k+1 == numPts ? poly->GetId(0) : poly->GetId(k+1);
1747
1748                                 if (cellB.f == idA
1749                                     && cellA.f == idB) {
1750
1751                                     IdsType &pts = sharedPts[Pair{sp.ind, cellA.g}];
1752
1753                                     pts.push_back(cellA.f);
1754                                     pts.push_back(cellB.f);
1755
1756                                     break;
1757                                 }
1758                             }
1759
1760                         }
1761                     }
1762                 }
1763
1764
1765             }
1766
1767         }
1768     }
1769
1770     double pt[3];
1771
1772     decltype(sharedPts)::iterator itr3;
1773
1774     for (itr3 = sharedPts.begin(); itr3 != sharedPts.end(); ++itr3) {
1775         const Pair &p = itr3->first;
1776
1777         IdsType &pts = itr3->second;
1778
1779         std::sort(pts.begin(), pts.end());
1780
1781         auto found = std::adjacent_find(pts.begin(), pts.end());
1782
1783         if (found != pts.end()) {
1784             contLines->GetPoint(p.f, pt);
1785
1786             i = pd->GetPoints()->InsertNextPoint(pt);
1787             pd->ReplaceCellPoint(p.g, *found, i);
1788         }
1789     }
1790
1791     links->Delete();
1792     ptsA->Delete();
1793     ptsB->Delete();
1794     cells->Delete();
1795     poly->Delete();
1796
1797     loc->FreeSearchStructure();
1798     loc->Delete();
1799
1800 }
1801
1802 void vtkPolyDataBooleanFilter::AddAdjacentPoints (vtkPolyData *pd, vtkIdTypeArray *conts, PolyStripsType &polyStrips) {
1803
1804 #ifdef DEBUG
1805     std::cout << "AddAdjacentPoints()" << std::endl;
1806 #endif
1807
1808     pd->BuildLinks();
1809
1810     vtkIdTypeArray *origCellIds = vtkIdTypeArray::SafeDownCast(pd->GetCellData()->GetScalars("OrigCellIds"));
1811
1812     struct Cmp {
1813         bool operator() (const StripPt &l, const StripPt &r) const {
1814             return l.t < r.t;
1815         }
1816     };
1817
1818     vtkKdTreePointLocator *loc = vtkKdTreePointLocator::New();
1819     loc->SetDataSet(pd);
1820     loc->BuildLocator();
1821
1822     vtkIdList *cells = vtkIdList::New(),
1823         *ptsA = vtkIdList::New(),
1824         *ptsB = vtkIdList::New(),
1825         *poly = vtkIdList::New(),
1826         *newPoly = vtkIdList::New();
1827
1828     vtkIdType i, j, numCells, numPtsA, numPtsB, numPts, idA, idB;
1829
1830     PolyStripsType::const_iterator itr;
1831     StripPtsType::const_iterator itr2;
1832
1833     for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) {
1834         const PStrips &pStrips = itr->second;
1835
1836         std::map<Pair, std::set<std::reference_wrapper<const StripPt>, Cmp>> edgePts;
1837
1838         for (itr2 = pStrips.pts.begin(); itr2 != pStrips.pts.end(); ++itr2) {
1839             const StripPt &sp = itr2->second;
1840
1841             if (sp.capt == Capt::EDGE) {
1842                 edgePts[{sp.edge[0], sp.edge[1]}].emplace(sp);
1843             }
1844         }
1845
1846         decltype(edgePts)::iterator itr3;
1847
1848         for (itr3 = edgePts.begin(); itr3 != edgePts.end(); ++itr3) {
1849             const Pair &edge = itr3->first;
1850             auto pts = itr3->second;
1851
1852             StripPt spA, spB;
1853
1854             pd->GetPoint(edge.f, spA.pt);
1855             pd->GetPoint(edge.g, spB.pt);
1856
1857             spA.t = 0;
1858             spB.t = 1;
1859
1860             pts.emplace(spA);
1861             pts.emplace(spB);
1862
1863             std::vector<decltype(pts)::value_type> _pts(pts.rbegin(), pts.rend());
1864
1865             decltype(_pts)::const_iterator itrA, itrB, itrC;
1866
1867             itrA = _pts.begin();
1868
1869             while (itrA != _pts.end()-1) {
1870                 itrB = itrA+1;
1871
1872                 while (itrB != _pts.end()-1) {
1873                     contLines->GetPointCells(itrB->get().ind, cells);
1874                     numCells = cells->GetNumberOfIds();
1875
1876                     // std::cout << "[";
1877                     // for (i = 0; i < numCells; i++) {
1878                     //     std::cout << conts->GetValue(cells->GetId(i)) << ", ";
1879                     // }
1880                     // std::cout << "]" << std::endl;
1881
1882                     // break;
1883
1884                     ++itrB;
1885                 }
1886
1887                 if (std::distance(itrA, itrB) > 1) {
1888                     FindPoints(loc, itrA->get().pt, ptsA);
1889                     FindPoints(loc, itrB->get().pt, ptsB);
1890
1891                     numPtsA = ptsA->GetNumberOfIds();
1892                     numPtsB = ptsB->GetNumberOfIds();
1893
1894                     std::vector<Pair> polysA, polysB;
1895
1896                     for (i = 0; i < numPtsA; i++) {
1897                         pd->GetPointCells(ptsA->GetId(i), cells);
1898                         numCells = cells->GetNumberOfIds();
1899
1900                         for (j = 0; j < numCells; j++) {
1901                             polysA.emplace_back(cells->GetId(j), ptsA->GetId(i));
1902                         }
1903                     }
1904
1905                     for (i = 0; i < numPtsB; i++) {
1906                         pd->GetPointCells(ptsB->GetId(i), cells);
1907                         numCells = cells->GetNumberOfIds();
1908
1909                         for (j = 0; j < numCells; j++) {
1910                             polysB.emplace_back(cells->GetId(j), ptsB->GetId(i));
1911                         }
1912                     }
1913
1914                     /*for (const Pair &pA : polysA) {
1915                         std::cout << "pA " << pA << std::endl;
1916                     }
1917
1918                     for (const Pair &pB : polysB) {
1919                         std::cout << "pB " << pB << std::endl;
1920                     }*/
1921
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) {
1925
1926                                 pd->GetCellPoints(pA.f, poly);
1927                                 numPts = poly->GetNumberOfIds();
1928
1929                                 newPoly->Reset();
1930
1931                                 for (j = 0; j < numPts; j++) {
1932                                     newPoly->InsertNextId(poly->GetId(j));
1933
1934                                     idA = poly->GetId(j);
1935                                     idB = j+1 == numPts ? poly->GetId(0) : poly->GetId(j+1);
1936
1937                                     if (pA.g == idA
1938                                         && pB.g == idB) {
1939
1940                                         for (itrC = itrA+1; itrC != itrB; ++itrC) {
1941                                             // newPoly->InsertNextId(pd->InsertNextLinkedPoint(itrC->get().pt, 1));
1942
1943                                             pd->InsertNextLinkedPoint(1);
1944
1945                                             newPoly->InsertNextId(pd->GetPoints()->InsertNextPoint(itrC->get().pt));
1946                                         }
1947
1948                                     }
1949
1950                                     pd->RemoveReferenceToCell(idA, pA.f);
1951                                 }
1952
1953                                 pd->DeleteCell(pA.f);
1954
1955                                 pd->InsertNextLinkedCell(VTK_POLYGON, newPoly->GetNumberOfIds(), newPoly->GetPointer(0));
1956
1957                                 origCellIds->InsertNextValue(origCellIds->GetValue(pA.f));
1958
1959                                 break;
1960                             }
1961                         }
1962                     }
1963                 }
1964
1965                 itrA = itrB;
1966             }
1967         }
1968     }
1969
1970     cells->Delete();
1971     ptsA->Delete();
1972     ptsB->Delete();
1973     poly->Delete();
1974     newPoly->Delete();
1975
1976     loc->FreeSearchStructure();
1977     loc->Delete();
1978
1979     pd->RemoveDeletedCells();
1980
1981 }
1982
1983 void vtkPolyDataBooleanFilter::MergePoints (vtkPolyData *pd, PolyStripsType &polyStrips) {
1984
1985 #ifdef DEBUG
1986     std::cout << "MergePoints()" << std::endl;
1987 #endif
1988
1989     pd->BuildCells();
1990     pd->BuildLinks();
1991
1992     contLines->BuildLinks();
1993
1994     vtkKdTreePointLocator *loc = vtkKdTreePointLocator::New();
1995     loc->SetDataSet(pd);
1996     loc->BuildLocator();
1997
1998     PolyStripsType::const_iterator itr;
1999     StripsType::const_iterator itr2;
2000
2001     std::map<vtkIdType, std::set<vtkIdType>> neighPts;
2002
2003     vtkIdList *pts = vtkIdList::New();
2004
2005     vtkIdType i, j, numPts;
2006
2007     for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) {
2008         const PStrips &pStrips = itr->second;
2009
2010         for (itr2 = pStrips.strips.begin(); itr2 != pStrips.strips.end(); ++itr2) {
2011             const StripType &strip = *itr2;
2012
2013             const StripPtR &spA = strip.front(),
2014                 &spB = strip.back();
2015
2016             const auto beforeA = pStrips.pts.find((strip.begin()+1)->ind),
2017                 beforeB = pStrips.pts.find((strip.end()-2)->ind);
2018
2019             FindPoints(loc, beforeA->second.pt, pts);
2020             numPts = pts->GetNumberOfIds();
2021
2022             for (i = 0; i < numPts; i++) {
2023                 neighPts[spA.ind].insert(pts->GetId(i));
2024             }
2025
2026             FindPoints(loc, beforeB->second.pt, pts);
2027             numPts = pts->GetNumberOfIds();
2028
2029             for (i = 0; i < numPts; i++) {
2030                 neighPts[spB.ind].insert(pts->GetId(i));
2031             }
2032         }
2033     }
2034
2035     double pt[3];
2036
2037     vtkIdList *polys = vtkIdList::New(),
2038         *poly = vtkIdList::New();
2039
2040     vtkIdType ind, polyId, _numPts, before, after;
2041
2042     decltype(neighPts)::const_iterator itr3;
2043
2044     for (itr3 = neighPts.begin(); itr3 != neighPts.end(); ++itr3) {
2045         const auto &inds = itr3->second;
2046
2047         std::map<Point3d, std::vector<Pair>> pairs;
2048
2049         contLines->GetPoint(itr3->first, pt);
2050
2051         FindPoints(loc, pt, pts);
2052         numPts = pts->GetNumberOfIds();
2053
2054         for (i = 0; i < numPts; i++) {
2055             ind = pts->GetId(i);
2056
2057             pd->GetPointCells(ind, polys);
2058
2059             if (polys->GetNumberOfIds() > 0) {
2060                 polyId = polys->GetId(0);
2061
2062                 pd->GetCellPoints(polyId, poly);
2063                 _numPts = poly->GetNumberOfIds();
2064
2065                 for (j = 0; j < _numPts; j++) {
2066                     if (poly->GetId(j) == ind) {
2067                         break;
2068                     }
2069                 }
2070
2071                 // wieder davor und danach ermitteln
2072
2073                 before = poly->GetId(j == 0 ? _numPts-1 : j-1);
2074                 after = poly->GetId(j+1 == _numPts ? 0 : j+1);
2075
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);
2079                 }
2080
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);
2084                 }
2085
2086             }
2087         }
2088
2089         std::deque<std::deque<std::reference_wrapper<const Pair>>> Pairs;
2090
2091         decltype(pairs)::const_iterator itr4;
2092
2093         for (itr4 = pairs.begin(); itr4 != pairs.end(); ++itr4) {
2094             const auto &p = itr4->second;
2095
2096             if (p.size() == 2) {
2097                 auto _pts = {std::ref(p.front()), std::ref(p.back())}; // std::initializer_list
2098                 Pairs.emplace_back(_pts);
2099             }
2100         }
2101
2102         decltype(Pairs)::iterator itr5;
2103
2104         /*for (itr5 = Pairs.begin(); itr5 != Pairs.end(); ++itr5) {
2105             for (auto &p : *itr5) {
2106                 std::cout << p.get() << ", ";
2107             }
2108             std::cout << std::endl;
2109         }*/
2110
2111         decltype(Pairs)::value_type group;
2112
2113         decltype(group)::const_iterator itr6;
2114
2115         while (!Pairs.empty()) {
2116             if (group.empty()) {
2117                 group = Pairs.front();
2118                 Pairs.pop_front();
2119             }
2120
2121             itr5 = Pairs.begin();
2122
2123             while (itr5 != Pairs.end()) {
2124                 const auto &next = *itr5;
2125
2126                 if (next.front().get() == group.front().get()) {
2127                     group.emplace_front(next.back());
2128                     Pairs.erase(itr5);
2129                     itr5 = Pairs.begin();
2130                 } else if (next.front().get() == group.back().get()) {
2131                     group.emplace_back(next.back());
2132                     Pairs.erase(itr5);
2133                     itr5 = Pairs.begin();
2134                 } else if (next.back().get() == group.front().get()) {
2135                     group.emplace_front(next.front());
2136                     Pairs.erase(itr5);
2137                     itr5 = Pairs.begin();
2138                 } else if (next.back().get() == group.back().get()) {
2139                     group.emplace_back(next.front());
2140                     Pairs.erase(itr5);
2141                     itr5 = Pairs.begin();
2142                 } else {
2143                     ++itr5;
2144                 }
2145             }
2146
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);
2150                 }
2151
2152                 group.clear();
2153             }
2154         }
2155
2156     }
2157
2158     polys->Delete();
2159     poly->Delete();
2160     pts->Delete();
2161
2162     loc->FreeSearchStructure();
2163     loc->Delete();
2164
2165 }
2166
2167 enum class Congr {
2168     EQUAL,
2169     OPPOSITE,
2170     NOT
2171 };
2172
2173 class PolyAtEdge {
2174     vtkPolyData *pd;
2175
2176 public:
2177     PolyAtEdge (vtkPolyData *_pd, vtkIdType _polyId, vtkIdType _ptIdA, vtkIdType _ptIdB) : pd(_pd), polyId(_polyId), ptIdA(_ptIdA), ptIdB(_ptIdB), loc(Loc::NONE) {
2178
2179         double ptA[3], ptB[3];
2180
2181         pd->GetPoint(ptIdA, ptA);
2182         pd->GetPoint(ptIdB, ptB);
2183
2184         vtkMath::Subtract(ptB, ptA, e);
2185         vtkMath::Normalize(e);
2186
2187         const vtkIdType *poly;
2188
2189         vtkIdType numPts;
2190         pd->GetCellPoints(polyId, numPts, poly);
2191
2192         ComputeNormal(pd->GetPoints(), n, numPts, poly);
2193
2194         vtkMath::Cross(e, n, r);
2195
2196     }
2197
2198     vtkIdType polyId, ptIdA, ptIdB;
2199     double n[3], e[3], r[3];
2200
2201     Loc loc;
2202
2203     friend std::ostream& operator<< (std::ostream &out, const PolyAtEdge &p) {
2204         out << "polyId " << p.polyId << ", ptIdA " << p.ptIdA << ", ptIdB " << p.ptIdB;
2205         return out;
2206     }
2207
2208     static constexpr double eps {.99999999}; // ~.0081deg
2209
2210     Congr IsCongruent (const PolyAtEdge &p) const {
2211         double cong = vtkMath::Dot(n, p.n);
2212
2213         if (cong > eps || cong < -eps) {
2214             double ang = vtkMath::Dot(r, p.r);
2215
2216             if (ang > eps) {
2217                 if (cong > eps) {
2218                     // normalen sind gleich ausgerichtet
2219                     return Congr::EQUAL;
2220                 } else {
2221                     return Congr::OPPOSITE;
2222                 }
2223             }
2224         }
2225
2226         return Congr::NOT;
2227     }
2228
2229 };
2230
2231
2232 class PolyPair {
2233 public:
2234     PolyPair (PolyAtEdge _pA, PolyAtEdge _pB) : pA(_pA), pB(_pB) {}
2235
2236     PolyAtEdge pA, pB;
2237
2238     void GetLoc (PolyAtEdge &pT, vtkIdType mode) {
2239
2240         Congr cA = pA.IsCongruent(pT),
2241             cB = pB.IsCongruent(pT);
2242
2243 #ifdef DEBUG
2244         std::cout << "GetLoc() -> polyId " << pT.polyId
2245                   << ", cA " << cA
2246                   << ", cB " << cB
2247                   << std::endl;
2248
2249         if (cA != Congr::NOT || cB != Congr::NOT) {
2250             assert(cA != cB);
2251         }
2252
2253 #endif
2254
2255         if (cA == Congr::EQUAL || cA == Congr::OPPOSITE) {
2256             if (cA == Congr::OPPOSITE) {
2257                 // normalen sind entgegengesetzt gerichtet
2258
2259                 if (mode == OPER_INTERSECTION) {
2260                     pA.loc = Loc::OUTSIDE;
2261                     pT.loc = Loc::OUTSIDE;
2262                 } else {
2263                     pA.loc = Loc::INSIDE;
2264                     pT.loc = Loc::INSIDE;
2265                 }
2266             } else if (mode == OPER_UNION || mode == OPER_INTERSECTION) {
2267                 pA.loc = Loc::INSIDE;
2268                 pT.loc = Loc::OUTSIDE;
2269             }
2270
2271         } else if (cB == Congr::EQUAL || cB == Congr::OPPOSITE) {
2272             if (cB == Congr::OPPOSITE) {
2273                 // normalen sind entgegengesetzt gerichtet
2274
2275                 if (mode == OPER_INTERSECTION) {
2276                     pB.loc = Loc::OUTSIDE;
2277                     pT.loc = Loc::OUTSIDE;
2278                 } else {
2279                     pB.loc = Loc::INSIDE;
2280                     pT.loc = Loc::INSIDE;
2281                 }
2282             } else if (mode == OPER_UNION || mode == OPER_INTERSECTION) {
2283                 pB.loc = Loc::INSIDE;
2284                 pT.loc = Loc::OUTSIDE;
2285             }
2286
2287         } else {
2288             double alpha = GetAngle(pA.r, pB.r, pA.e),
2289                 beta = GetAngle(pA.r, pT.r, pA.e);
2290
2291             if (beta > alpha) {
2292                 pT.loc = Loc::INSIDE;
2293             } else {
2294                 pT.loc = Loc::OUTSIDE;
2295             }
2296         }
2297     }
2298
2299 };
2300
2301
2302 std::shared_ptr<PolyPair> GetEdgePolys (vtkPolyData *pd, vtkIdList *ptsA, vtkIdList *ptsB) {
2303
2304 #ifdef DEBUG
2305     std::cout << "GetEdgePolys()" << std::endl;
2306 #endif
2307
2308     std::vector<Pair> p;
2309
2310     vtkIdType numPtsA = ptsA->GetNumberOfIds(),
2311         numPtsB = ptsB->GetNumberOfIds();
2312
2313     vtkIdList *polys = vtkIdList::New();
2314
2315     vtkIdType i, j, numCells;
2316
2317     for (i = 0; i < numPtsA; i++) {
2318         pd->GetPointCells(ptsA->GetId(i), polys);
2319         numCells = polys->GetNumberOfIds();
2320
2321         for (j = 0; j < numCells; j++) {
2322             p.emplace_back(ptsA->GetId(i), polys->GetId(j));
2323         }
2324     }
2325
2326     for (i = 0; i < numPtsB; i++) {
2327         pd->GetPointCells(ptsB->GetId(i), polys);
2328         numCells = polys->GetNumberOfIds();
2329
2330         for (j = 0; j < numCells; j++) {
2331             p.emplace_back(ptsB->GetId(i), polys->GetId(j));
2332         }
2333     }
2334
2335     polys->Delete();
2336
2337     std::map<vtkIdType, IdsType> pEdges;
2338
2339     std::vector<Pair>::const_iterator itr;
2340     for (itr = p.begin(); itr != p.end(); ++itr) {
2341         pEdges[itr->g].push_back(itr->f);
2342     }
2343
2344     std::vector<PolyAtEdge> opp;
2345
2346     vtkIdList *poly = vtkIdList::New();
2347
2348     vtkIdType numPts, a, b;
2349
2350     std::map<vtkIdType, IdsType>::const_iterator itr2;
2351
2352     for (itr2 = pEdges.begin(); itr2 != pEdges.end(); ++itr2) {
2353         const IdsType &pts = itr2->second;
2354
2355         if (pts.size() > 1) {
2356             pd->GetCellPoints(itr2->first, poly);
2357             numPts = poly->GetNumberOfIds();
2358
2359             for (i = 0; i < numPts; i++) {
2360                 a = poly->GetId(i);
2361                 b = i+1 == numPts ? poly->GetId(0) : poly->GetId(i+1);
2362
2363                 if (std::find(pts.begin(), pts.end(), a) != pts.end()
2364                     && std::find(pts.begin(), pts.end(), b) != pts.end()) {
2365
2366                     opp.emplace_back(pd, itr2->first, a, b);
2367                 }
2368             }
2369
2370         }
2371     }
2372
2373     poly->Delete();
2374
2375 #ifdef DEBUG
2376     for (auto &op : opp) {
2377         std::cout << op << std::endl;
2378     }
2379 #endif
2380
2381     if (opp.size() != 2) {
2382         return nullptr;
2383     }
2384
2385     return std::make_shared<PolyPair>(opp[0], opp[1]);
2386
2387 }
2388
2389 bool vtkPolyDataBooleanFilter::CombineRegions () {
2390
2391 #ifdef DEBUG
2392     std::cout << "CombineRegions()" << std::endl;
2393 #endif
2394
2395     vtkSmartPointer<vtkPolyData> filterdA = vtkSmartPointer<vtkPolyData>::New();
2396     filterdA->DeepCopy(modPdA);
2397
2398     vtkSmartPointer<vtkPolyData> filterdB = vtkSmartPointer<vtkPolyData>::New();
2399     filterdB->DeepCopy(modPdB);
2400
2401     // ungenutzte punkte löschen
2402     vtkSmartPointer<vtkCleanPolyData> cleanA = vtkSmartPointer<vtkCleanPolyData>::New();
2403     cleanA->PointMergingOff();
2404     cleanA->SetInputData(filterdA);
2405
2406     vtkSmartPointer<vtkCleanPolyData> cleanB = vtkSmartPointer<vtkCleanPolyData>::New();
2407     cleanB->PointMergingOff();
2408     cleanB->SetInputData(filterdB);
2409
2410     // regionen mit skalaren ausstatten
2411     vtkSmartPointer<vtkPolyDataConnectivityFilter> cfA = vtkSmartPointer<vtkPolyDataConnectivityFilter>::New();
2412     cfA->SetExtractionModeToAllRegions();
2413     cfA->ColorRegionsOn();
2414     cfA->SetInputConnection(cleanA->GetOutputPort());
2415
2416     vtkSmartPointer<vtkPolyDataConnectivityFilter> cfB = vtkSmartPointer<vtkPolyDataConnectivityFilter>::New();
2417     cfB->SetExtractionModeToAllRegions();
2418     cfB->ColorRegionsOn();
2419     cfB->SetInputConnection(cleanB->GetOutputPort());
2420
2421     cfA->Update();
2422     cfB->Update();
2423
2424     vtkPolyData *pdA = cfA->GetOutput();
2425     vtkPolyData *pdB = cfB->GetOutput();
2426
2427 #ifdef DEBUG
2428     std::cout << "Exporting modPdA_8.vtk" << std::endl;
2429     WriteVTK("modPdA_8.vtk", pdA);
2430
2431     std::cout << "Exporting modPdB_8.vtk" << std::endl;
2432     WriteVTK("modPdB_8.vtk", pdB);
2433 #endif
2434
2435     resultC->ShallowCopy(contLines);
2436
2437     if (OperMode == OPER_NONE) {
2438         resultA->ShallowCopy(pdA);
2439         resultB->ShallowCopy(pdB);
2440
2441         return false;
2442     }
2443
2444     vtkSmartPointer<vtkKdTreePointLocator> plA = vtkSmartPointer<vtkKdTreePointLocator>::New();
2445     plA->SetDataSet(pdA);
2446     plA->BuildLocator();
2447
2448     vtkSmartPointer<vtkKdTreePointLocator> plB = vtkSmartPointer<vtkKdTreePointLocator>::New();
2449     plB->SetDataSet(pdB);
2450     plB->BuildLocator();
2451
2452     pdA->BuildLinks();
2453     pdB->BuildLinks();
2454
2455     vtkIdTypeArray *scalarsA = vtkIdTypeArray::SafeDownCast(pdA->GetPointData()->GetScalars());
2456     vtkIdTypeArray *scalarsB = vtkIdTypeArray::SafeDownCast(pdB->GetPointData()->GetScalars());
2457
2458     vtkSmartPointer<vtkIdList> line = vtkSmartPointer<vtkIdList>::New();
2459
2460     double ptA[3], ptB[3];
2461
2462     vtkSmartPointer<vtkIdList> fptsA = vtkSmartPointer<vtkIdList>::New();
2463     vtkSmartPointer<vtkIdList> lptsA = vtkSmartPointer<vtkIdList>::New();
2464
2465     vtkSmartPointer<vtkIdList> fptsB = vtkSmartPointer<vtkIdList>::New();
2466     vtkSmartPointer<vtkIdList> lptsB = vtkSmartPointer<vtkIdList>::New();
2467
2468     std::map<vtkIdType, Loc> locsA, locsB;
2469
2470     vtkIdType i, j, numLines = contLines->GetNumberOfCells();
2471
2472     IdsType _failed;
2473
2474     for (i = 0; i < numLines; i++) {
2475
2476         if (contLines->GetCellType(i) == VTK_EMPTY_CELL) {
2477             continue;
2478         }
2479
2480         contLines->GetCellPoints(i, line);
2481
2482         contLines->GetPoint(line->GetId(0), ptA);
2483         contLines->GetPoint(line->GetId(1), ptB);
2484
2485         FindPoints(plA, ptA, fptsA);
2486         FindPoints(plB, ptA, fptsB);
2487
2488 #ifdef DEBUG
2489         std::cout << "line " << i << std::endl;
2490 #else
2491
2492         // bereits behandelte regionen werden nicht noch einmal untersucht
2493
2494         vtkIdType notLocated = 0;
2495
2496         for (j = 0; j < fptsA->GetNumberOfIds(); j++) {
2497             if (locsA.count(scalarsA->GetValue(fptsA->GetId(j))) == 0) {
2498                 notLocated++;
2499             }
2500         }
2501
2502         for (j = 0; j < fptsB->GetNumberOfIds(); j++) {
2503             if (locsB.count(scalarsB->GetValue(fptsB->GetId(j))) == 0) {
2504                 notLocated++;
2505             }
2506         }
2507
2508         if (notLocated == 0) {
2509             continue;
2510         }
2511
2512 #endif
2513
2514         FindPoints(plA, ptB, lptsA);
2515         FindPoints(plB, ptB, lptsB);
2516
2517         auto ppA = GetEdgePolys(pdA, fptsA, lptsA);
2518         auto ppB = GetEdgePolys(pdB, fptsB, lptsB);
2519
2520         if (ppA && ppB) {
2521
2522             ppB->GetLoc(ppA->pA, OperMode);
2523             ppB->GetLoc(ppA->pB, OperMode);
2524
2525             ppA->GetLoc(ppB->pA, OperMode);
2526             ppA->GetLoc(ppB->pB, OperMode);
2527
2528             vtkIdType fsA = scalarsA->GetValue(ppA->pA.ptIdA);
2529             vtkIdType lsA = scalarsA->GetValue(ppA->pB.ptIdA);
2530
2531             vtkIdType fsB = scalarsB->GetValue(ppB->pA.ptIdA);
2532             vtkIdType lsB = scalarsB->GetValue(ppB->pB.ptIdA);
2533
2534 #ifdef DEBUG
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;
2539
2540             if (locsA.count(fsA) > 0 && locsA[fsA] != ppA->pA.loc) {
2541                 std::cout << "sA " << fsA << ": " << locsA[fsA] << " -> " << ppA->pA.loc << std::endl;
2542             }
2543
2544             if (locsA.count(lsA) > 0 && locsA[lsA] != ppA->pB.loc) {
2545                 std::cout << "sA " << lsA << ": " << locsA[lsA] << " -> " << ppA->pB.loc << std::endl;
2546             }
2547
2548             if (locsB.count(fsB) > 0 && locsB[fsB] != ppB->pA.loc) {
2549                 std::cout << "sB " << fsB << ": " << locsB[fsB] << " -> " << ppB->pA.loc << std::endl;
2550             }
2551
2552             if (locsB.count(lsB) > 0 && locsB[lsB] != ppB->pB.loc) {
2553                 std::cout << "sB " << lsB << ": " << locsB[lsB] << " -> " << ppB->pB.loc << std::endl;
2554             }
2555
2556 #endif
2557
2558             locsA.emplace(fsA, ppA->pA.loc);
2559             locsA.emplace(lsA, ppA->pB.loc);
2560
2561             locsB.emplace(fsB, ppB->pA.loc);
2562             locsB.emplace(lsB, ppB->pB.loc);
2563
2564         } else {
2565             _failed.push_back(i);
2566
2567             // return true;
2568         }
2569
2570     }
2571
2572     if (_failed.size() > 0) {
2573         for (auto i : _failed) {
2574             std::cout << "failed at " << i
2575                 << std::endl;
2576         }
2577
2578         return true;
2579     }
2580
2581     // reale kombination der ermittelten regionen
2582
2583     Loc comb[] = {Loc::OUTSIDE, Loc::OUTSIDE};
2584
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;
2592     }
2593
2594     vtkIdType numA = cfA->GetNumberOfExtractedRegions(),
2595         numB = cfB->GetNumberOfExtractedRegions();
2596
2597     cfA->SetExtractionModeToSpecifiedRegions();
2598     cfB->SetExtractionModeToSpecifiedRegions();
2599
2600     std::map<vtkIdType, Loc>::const_iterator itr;
2601
2602     for (itr = locsA.begin(); itr != locsA.end(); itr++) {
2603         if (itr->second == comb[0]) {
2604             cfA->AddSpecifiedRegion(itr->first);
2605         }
2606     }
2607
2608     for (itr = locsB.begin(); itr != locsB.end(); itr++) {
2609         if (itr->second == comb[1]) {
2610             cfB->AddSpecifiedRegion(itr->first);
2611         }
2612     }
2613
2614     // nicht beteiligte regionen hinzufügen
2615
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);
2620             }
2621         }
2622     }
2623
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);
2628             }
2629         }
2630     }
2631
2632     // nach innen zeigende normalen umkehren
2633
2634     cfA->Update();
2635     cfB->Update();
2636
2637     vtkPolyData *regsA = cfA->GetOutput();
2638     vtkPolyData *regsB = cfB->GetOutput();
2639
2640     scalarsA = vtkIdTypeArray::SafeDownCast(regsA->GetPointData()->GetScalars());
2641     scalarsB = vtkIdTypeArray::SafeDownCast(regsB->GetPointData()->GetScalars());
2642
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);
2648                 }
2649             }
2650         }
2651
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);
2656                 }
2657             }
2658         }
2659     }
2660
2661     // OrigCellIds und CellData
2662
2663     vtkIdTypeArray *origCellIdsA = vtkIdTypeArray::SafeDownCast(regsA->GetCellData()->GetScalars("OrigCellIds"));
2664     vtkIdTypeArray *origCellIdsB = vtkIdTypeArray::SafeDownCast(regsB->GetCellData()->GetScalars("OrigCellIds"));
2665
2666     vtkIdTypeArray *newOrigCellIdsA = vtkIdTypeArray::New();
2667     newOrigCellIdsA->SetName("OrigCellIdsA");
2668
2669     vtkIdTypeArray *newOrigCellIdsB = vtkIdTypeArray::New();
2670     newOrigCellIdsB->SetName("OrigCellIdsB");
2671
2672     vtkCellData *newCellDataA = vtkCellData::New();
2673     newCellDataA->CopyAllocate(cellDataA);
2674
2675     vtkCellData *newCellDataB = vtkCellData::New();
2676     newCellDataB->CopyAllocate(cellDataB);
2677
2678     vtkIdType cellId;
2679
2680     for (i = 0; i < regsA->GetNumberOfCells(); i++) {
2681         cellId = cellIdsA->GetValue(origCellIdsA->GetValue(i));
2682
2683         newOrigCellIdsA->InsertNextValue(cellId);
2684         newOrigCellIdsB->InsertNextValue(-1);
2685
2686         newCellDataA->CopyData(cellDataA, cellId, i);
2687     }
2688
2689     for (i = 0; i < regsB->GetNumberOfCells(); i++) {
2690         cellId = cellIdsB->GetValue(origCellIdsB->GetValue(i));
2691
2692         newOrigCellIdsB->InsertNextValue(cellId);
2693         newOrigCellIdsA->InsertNextValue(-1);
2694
2695         newCellDataB->CopyData(cellDataB, cellId, i);
2696     }
2697
2698     regsA->GetCellData()->Initialize();
2699     regsB->GetCellData()->Initialize();
2700
2701     regsA->GetCellData()->ShallowCopy(newCellDataA);
2702     regsB->GetCellData()->ShallowCopy(newCellDataB);
2703
2704     newCellDataA->Delete();
2705     newCellDataB->Delete();
2706
2707     // zusammenführung
2708
2709     vtkAppendPolyData *app = vtkAppendPolyData::New();
2710     app->AddInputData(regsA);
2711     app->AddInputData(regsB);
2712
2713     // entfernt ungenutzte punkte
2714     vtkCleanPolyData *cleanApp = vtkCleanPolyData::New();
2715     cleanApp->PointMergingOff();
2716     cleanApp->SetInputConnection(app->GetOutputPort());
2717
2718     // färbt die regionen nochmal neu ein, damit mehrere regionen nicht die gleiche farbe haben
2719
2720     vtkPolyDataConnectivityFilter *cfApp = vtkPolyDataConnectivityFilter::New();
2721     cfApp->SetExtractionModeToAllRegions();
2722     cfApp->ColorRegionsOn();
2723     cfApp->SetInputConnection(cleanApp->GetOutputPort());
2724
2725     cfApp->Update();
2726
2727     vtkPolyData *cfPd = cfApp->GetOutput();
2728
2729     // resultB bleibt hier leer
2730
2731     resultA->ShallowCopy(cfPd);
2732
2733     resultA->GetCellData()->AddArray(newOrigCellIdsA);
2734     resultA->GetCellData()->AddArray(newOrigCellIdsB);
2735
2736     // aufräumen
2737
2738     cfApp->Delete();
2739     cleanApp->Delete();
2740     app->Delete();
2741
2742     newOrigCellIdsB->Delete();
2743     newOrigCellIdsA->Delete();
2744
2745     plB->FreeSearchStructure();
2746     plA->FreeSearchStructure();
2747
2748     return false;
2749
2750 }
2751
2752 Merger::Merger (vtkPolyData *pd, const PStrips &pStrips, const StripsType &strips, const IdsType &descIds, vtkIdType origId) : pd(pd), pStrips(pStrips), origId(origId) {
2753
2754     const StripPtsType &pts = pStrips.pts;
2755     const Base &base = pStrips.base;
2756
2757     vtkIdType i, num;
2758     const vtkIdType *cell;
2759
2760     double pt[3];
2761
2762     for (auto id : descIds) {
2763         pd->GetCellPoints(id, num, cell);
2764
2765         Poly p;
2766
2767         for (i = 0; i < num; i++) {
2768             pd->GetPoint(cell[i], pt);
2769
2770             double proj[2];
2771             Transform(pt, proj, base);
2772
2773             p.emplace_back(proj[0], proj[1], 0, cell[i]);
2774
2775         }
2776
2777         polys.push_back(p);
2778
2779         pd->DeleteCell(id);
2780     }
2781
2782     for (auto &strip : strips) {
2783         Poly p;
2784
2785         for (auto &sp : strip) {
2786             const double *pt = pts.at(sp.ind).pt;
2787
2788             double proj[2];
2789             Transform(pt, proj, base);
2790
2791             p.emplace_back(proj[0], proj[1], 0);
2792         }
2793
2794         p.pop_back();
2795
2796         double n[3];
2797         ComputeNormal(p, n);
2798
2799         if (n[2] < 0) {
2800             Poly q(p.rbegin(), p.rend());
2801             p.swap(q);
2802         }
2803
2804         polys.push_back(p);
2805     }
2806
2807 }
2808
2809 void Merger::Run () {
2810     // mergen mit hilfe von vtkKdTree und vtkModifiedBSPTree
2811
2812     vtkPoints *pdPts = pd->GetPoints();
2813     vtkIdTypeArray *origCellIds = vtkIdTypeArray::SafeDownCast(pd->GetCellData()->GetScalars("OrigCellIds"));
2814
2815     assert(origCellIds != nullptr);
2816
2817     const Base &base = pStrips.base;
2818
2819     std::vector<GroupType> groups(polys.size());
2820
2821     PolysType::const_iterator itrA, itrB;
2822
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()));
2827             }
2828         }
2829     }
2830
2831     std::size_t parent = 0;
2832
2833     for (auto &group : groups) {
2834         GroupType parents;
2835
2836         for (auto &index : group) {
2837             const GroupType &_group = groups[index];
2838             parents.insert(parents.end(), _group.begin(), _group.end());
2839         }
2840
2841         std::sort(group.begin(), group.end());
2842         std::sort(parents.begin(), parents.end());
2843
2844         GroupType _group {parent++};
2845         std::set_difference(group.begin(), group.end(), parents.begin(), parents.end(), std::back_inserter(_group));
2846
2847         std::cout << "[";
2848         for (auto &index : _group) {
2849             std::cout << index << ", ";
2850         }
2851         std::cout << "]" << std::endl;
2852
2853         PolysType merged;
2854
2855         MergeGroup(_group, merged);
2856
2857         std::map<Point3d, vtkIdType> newIds;
2858
2859         for (auto &poly : merged) {
2860             vtkSmartPointer<vtkIdList> newCell = vtkSmartPointer<vtkIdList>::New();
2861
2862             for (auto &p : poly) {
2863                 double in[] = {p.x, p.y},
2864                     out[3];
2865
2866                 BackTransform(in, out, base);
2867
2868                 vtkIdType id = p.id;
2869
2870                 if (id == NOTSET) {
2871                     auto itr = newIds.find(p);
2872
2873                     if (itr == newIds.end()) {
2874                         id = pdPts->InsertNextPoint(out);
2875                         newIds.emplace(p, id);
2876                     } else {
2877                         id = itr->second;
2878                     }
2879                 }
2880
2881                 newCell->InsertNextId(id);
2882
2883             }
2884
2885             pd->InsertNextCell(VTK_POLYGON, newCell);
2886             origCellIds->InsertNextValue(origId);
2887         }
2888
2889     }
2890
2891 }
2892
2893 void Merger::MergeGroup (const GroupType &group, PolysType &merged) {
2894     if (group.size() == 1) {
2895         merged.push_back(polys.at(group.back()));
2896
2897         return;
2898     }
2899
2900     vtkSmartPointer<vtkPoints> pts = vtkSmartPointer<vtkPoints>::New();
2901
2902     IndexedPolysType indexedPolys;
2903
2904     ReferencedPointsType refPts;
2905
2906     SourcesType sources;
2907     std::size_t src = 0;
2908
2909     for (auto &index : group) {
2910         const Poly &poly = polys.at(index);
2911
2912         IndexedPoly ids;
2913
2914         for (auto &p : poly) {
2915             vtkIdType id = pts->InsertNextPoint(p.x, p.y, p.z);
2916
2917             ids.push_back(id);
2918             sources.emplace(id, src);
2919
2920             refPts.emplace(id, p);
2921         }
2922
2923         indexedPolys.push_back(std::move(ids));
2924         src++;
2925     }
2926
2927     vtkSmartPointer<vtkKdTree> kdTree = vtkSmartPointer<vtkKdTree>::New();
2928     kdTree->OmitZPartitioning();
2929     kdTree->BuildLocatorFromPoints(pts);
2930
2931     vtkSmartPointer<vtkPolyData> linesA = vtkSmartPointer<vtkPolyData>::New();
2932     linesA->SetPoints(pts);
2933     linesA->Allocate(1);
2934
2935     IndexedPoly::const_iterator itrA, itrB;
2936
2937     for (const auto &ids : indexedPolys) {
2938         for (itrA = ids.begin(); itrA != ids.end(); ++itrA) {
2939             itrB = itrA+1;
2940             if (itrB == ids.end()) {
2941                 itrB = ids.begin();
2942             }
2943
2944             vtkIdList *line = vtkIdList::New();
2945             line->InsertNextId(*itrA);
2946             line->InsertNextId(*itrB);
2947
2948             linesA->InsertNextCell(VTK_LINE, line);
2949
2950             line->Delete();
2951         }
2952     }
2953
2954     //EED 2023-07-24
2955 //    WriteVTK("linesA.vtk", linesA);
2956
2957     vtkSmartPointer<vtkModifiedBSPTree> bspTreeA = vtkSmartPointer<vtkModifiedBSPTree>::New();
2958     bspTreeA->SetDataSet(linesA);
2959
2960     int n = 0;
2961
2962     PolyConnsType polyConns;
2963
2964     FindConns(linesA, kdTree, bspTreeA, polyConns, indexedPolys, sources, n);
2965
2966     PolyConnsType connected {{0, {}}};
2967     std::set<vtkIdType> restricted; // keine der conns darf im gleichen punkt beginnen
2968
2969     vtkSmartPointer<vtkPolyData> linesB = vtkSmartPointer<vtkPolyData>::New();
2970     linesB->SetPoints(pts);
2971     linesB->Allocate(1);
2972
2973     vtkSmartPointer<vtkModifiedBSPTree> bspTreeB = vtkSmartPointer<vtkModifiedBSPTree>::New();
2974     bspTreeB->SetDataSet(linesB);
2975
2976     ConnsType firstConns;
2977
2978     std::size_t i, numPolys = indexedPolys.size();
2979
2980     double ptA[3], ptB[3];
2981
2982     while (connected.size() < numPolys) {
2983
2984         bool foundOne = false;
2985
2986         for (i = 1; i < numPolys; i++) {
2987             if (connected.count(i) == 0) {
2988                 const ConnsType &conns = polyConns[i];
2989
2990                 for (auto &conn : conns) {
2991                     if (connected.count(sources.at(conn.j)) == 1
2992                         && restricted.count(conn.j) == 0) {
2993
2994                         pts->GetPoint(conn.i, ptA);
2995                         pts->GetPoint(conn.j, ptB);
2996
2997                         if (bspTreeB->IntersectWithLine(ptA, ptB, 1e-5, nullptr, nullptr) == 0) {
2998                             connected[sources.at(conn.i)].push_back(conn);
2999
3000                             // das andere poly auch aktualisieren
3001                             connected[sources.at(conn.j)].emplace_back(conn.d, conn.j, conn.i);
3002
3003                             restricted.insert(conn.i);
3004                             restricted.insert(conn.j);
3005
3006                             vtkIdList *line = vtkIdList::New();
3007                             line->InsertNextId(conn.i);
3008                             line->InsertNextId(conn.j);
3009
3010                             linesB->InsertNextCell(VTK_LINE, line);
3011
3012                             line->Delete();
3013
3014                             bspTreeB->Modified();
3015
3016                             foundOne = true;
3017
3018                             firstConns.push_back(conn);
3019
3020                             break;
3021                         }
3022
3023                     }
3024                 }
3025             }
3026         }
3027
3028         if (!foundOne) {
3029             if (!FindConns(linesA, kdTree, bspTreeA, polyConns, indexedPolys, sources, n)) {
3030                 throw std::runtime_error("Merging failed.");
3031             }
3032         }
3033     }
3034
3035     std::map<std::size_t, std::vector<std::size_t>> chains;
3036
3037     PolyConnsType::const_iterator itrC;
3038
3039     for (itrC = connected.begin(); itrC != connected.end(); ++itrC) {
3040         auto &chain = chains[itrC->first];
3041         chain.push_back(itrC->first);
3042
3043         while (chain.back() != 0) {
3044             chain.push_back(sources.at(connected.at(chain.back()).front().j));
3045         }
3046     }
3047
3048     std::cout << connected;
3049
3050     decltype(chains)::const_iterator itrD;
3051
3052     for (itrD = chains.begin(); itrD != chains.end(); ++itrD) {
3053         std::cout << itrD->first << ": [";
3054         for (auto &id : itrD->second) {
3055             std::cout << id << ", ";
3056         }
3057         std::cout << "]" << std::endl;
3058     }
3059
3060     std::set<std::size_t> solved {0};
3061
3062     std::deque<std::size_t> searchInds;
3063
3064     for (i = 1; i < numPolys; i++) {
3065         if (connected.at(i).size() == 1) {
3066             searchInds.push_back(i);
3067         }
3068     }
3069
3070     while (!searchInds.empty()) {
3071         PriosType prios;
3072
3073         for (auto ind : searchInds) {
3074             PolyPriosType polyPrios;
3075
3076             std::cout << "ind " << ind << std::endl;
3077
3078             const Conn &first = connected.at(ind).back();
3079
3080             for (auto &conn : polyConns.at(ind)) {
3081                 auto &src = sources.at(conn.j);
3082
3083                 if (polyPrios.count(src) == 1) {
3084                     continue;
3085                 }
3086
3087                 if (conn.i != first.i
3088                     && conn.j != first.j
3089                     && restricted.count(conn.i) == 0
3090                     && restricted.count(conn.j) == 0) {
3091
3092                     pts->GetPoint(conn.i, ptA);
3093                     pts->GetPoint(conn.j, ptB);
3094
3095                     if (bspTreeB->IntersectWithLine(ptA, ptB, 1e-5, nullptr, nullptr) == 0) {
3096                         auto &chainA = chains.at(ind),
3097                             &chainB = chains.at(src);
3098
3099                         std::set<std::size_t> _chainA(chainA.begin(), chainA.end()),
3100                             _chainB(chainB.begin(), chainB.end());
3101
3102                         // gemeinsame eltern
3103                         std::set<std::size_t> shared;
3104
3105                         std::set_intersection(_chainA.begin(), _chainA.end(), _chainB.begin(), _chainB.end(), std::inserter(shared, shared.end()));
3106
3107                         // gemeinsame eltern müssen sich alle in solved befinden
3108                         if (std::includes(solved.begin(), solved.end(), shared.begin(), shared.end())) {
3109                             std::set<std::size_t> solvable;
3110
3111                             std::set_difference(_chainA.begin(), _chainA.end(), solved.begin(), solved.end(), std::inserter(solvable, solvable.end()));
3112                             std::set_difference(_chainB.begin(), _chainB.end(), solved.begin(), solved.end(), std::inserter(solvable, solvable.end()));
3113
3114                             polyPrios.emplace(std::piecewise_construct,
3115                                 std::forward_as_tuple(src),
3116                                 std::forward_as_tuple(conn, solvable, -conn.d));
3117                         }
3118                     }
3119                 }
3120             }
3121
3122             PolyPriosType::const_iterator itr;
3123             for (itr = polyPrios.begin(); itr != polyPrios.end(); ++itr) {
3124                 prios.insert(itr->second);
3125             }
3126         }
3127
3128         if (!prios.empty()) {
3129             auto &prio = *prios.rbegin();
3130
3131             std::cout << "found " << prio << std::endl;
3132
3133             auto &conns = connected.at(sources.at(prio.conn.i));
3134
3135             conns.push_back(prio.conn);
3136
3137             connected.at(sources.at(prio.conn.j)).emplace_back(prio.conn.d, prio.conn.j, prio.conn.i);
3138
3139             restricted.insert(prio.conn.i);
3140             restricted.insert(prio.conn.j);
3141
3142             vtkIdList *line = vtkIdList::New();
3143             line->InsertNextId(prio.conn.i);
3144             line->InsertNextId(prio.conn.j);
3145
3146             linesB->InsertNextCell(VTK_LINE, line);
3147
3148             line->Delete();
3149
3150             bspTreeB->Modified();
3151
3152             solved.insert(prio.solvable.begin(), prio.solvable.end());
3153
3154             searchInds.erase(std::find(searchInds.begin(), searchInds.end(), sources.at(prio.conn.i)));
3155
3156             auto itr = std::find(searchInds.begin(), searchInds.end(), sources.at(prio.conn.j));
3157
3158             if (itr != searchInds.end()) {
3159                 searchInds.erase(itr);
3160             }
3161         } else {
3162             if (!FindConns(linesA, kdTree, bspTreeA, polyConns, indexedPolys, sources, n)) {
3163                 throw std::runtime_error("Merging failed.");
3164             }
3165         }
3166     }
3167
3168     std::cout << connected;
3169
3170 //EED 2023-07-24
3171 //    WriteVTK("linesB.vtk", linesB);
3172
3173     // stage 1
3174
3175     struct Cmp {
3176         bool operator() (const Conn &a, const Conn &b) const {
3177             return std::tie(a.i, a.j) < std::tie(b.i, b.j);
3178         }
3179     };
3180
3181     std::set<Conn, Cmp> usedConns;
3182
3183     IndexedPoly polyA {indexedPolys.front()};
3184
3185     for (const auto &first : firstConns) {
3186         auto itrA = std::find(polyA.begin(), polyA.end(), first.j);
3187
3188         assert(itrA != polyA.end());
3189
3190         auto &polyB = indexedPolys.at(sources.at(first.i));
3191
3192         auto itrB = std::find(polyB.begin(), polyB.end(), first.i);
3193
3194         assert(itrB != polyB.end());
3195
3196         std::rotate(polyA.begin(), itrA, polyA.end());
3197         std::rotate(polyB.begin(), itrB, polyB.end());
3198
3199         IndexedPoly newPoly {polyA};
3200         newPoly.push_back(polyA.front());
3201         newPoly.push_back(polyB.front());
3202
3203         newPoly.insert(newPoly.end(), polyB.rbegin(), polyB.rend());
3204
3205         polyA.swap(newPoly);
3206
3207         usedConns.insert(first);
3208
3209     }
3210
3211     PolysType newPolysA;
3212     GetPolys(refPts, {polyA}, newPolysA);
3213
3214     //EED 2023-07-24
3215 //    WritePolys("merged_stage1.vtk", newPolysA);
3216
3217     // stage 2
3218
3219     std::set<Conn, Cmp> leftConns;
3220
3221     for (itrC = connected.begin(); itrC != connected.end(); ++itrC) {
3222         if (itrC->first == 0) {
3223             continue;
3224         }
3225
3226         auto &conns = itrC->second;
3227
3228         ConnsType::const_iterator itr;
3229
3230         for (itr = conns.begin()+1; itr != conns.end(); ++itr) {
3231             Conn conn(0, itr->j, itr->i);
3232
3233             if (usedConns.find(conn) == usedConns.end()) {
3234                 if (itr->i < itr->j) {
3235                     leftConns.emplace(0, itr->i, itr->j);
3236                 } else {
3237                     leftConns.insert(std::move(conn));
3238                 }
3239             }
3240         }
3241
3242     }
3243
3244     std::cout << "leftConns: [";
3245     for (auto &conn : leftConns) {
3246         std::cout << conn << ", ";
3247     }
3248     std::cout << "]" << std::endl;
3249
3250     IndexedPolysType splitted {polyA};
3251
3252     for (auto &conn : leftConns) {
3253         IndexedPolysType::iterator itr;
3254
3255         for (itr = splitted.begin(); itr != splitted.end(); ++itr) {
3256             auto &poly = *itr;
3257
3258             auto itrA = std::find(poly.begin(), poly.end(), conn.i);
3259
3260             if (itrA != poly.end()) {
3261                 std::rotate(poly.begin(), itrA, poly.end());
3262
3263                 auto itrB = std::find(poly.begin(), poly.end(), conn.j);
3264
3265                 IndexedPoly newPolyA(poly.begin(), itrB+1);
3266                 IndexedPoly newPolyB(itrB, poly.end());
3267
3268                 newPolyB.push_back(poly.front());
3269
3270                 splitted.erase(itr);
3271
3272                 splitted.push_back(std::move(newPolyA));
3273                 splitted.push_back(std::move(newPolyB));
3274
3275                 break;
3276             }
3277         }
3278     }
3279
3280     PolysType newPolysB;
3281     GetPolys(refPts, splitted, newPolysB);
3282
3283     //EED 2023-07-24
3284 //    WritePolys("merged_stage2.vtk", newPolysB);
3285
3286     std::move(newPolysB.begin(), newPolysB.end(), std::back_inserter(merged));
3287
3288 }
3289
3290 bool Merger::FindConns (vtkPolyData *lines, vtkSmartPointer<vtkKdTree> kdTree, vtkSmartPointer<vtkModifiedBSPTree> bspTree, PolyConnsType &polyConns, const IndexedPolysType &indexedPolys, const SourcesType &sources, int &n) {
3291
3292     vtkPoints *pts = lines->GetPoints();
3293
3294     if (n > pts->GetNumberOfPoints()) {
3295         return false;
3296     }
3297
3298     n += 10;
3299
3300     vtkSmartPointer<vtkIdList> foundPts = vtkSmartPointer<vtkIdList>::New();
3301
3302     vtkIdType i, numPts;
3303
3304     vtkIdType idB;
3305
3306     vtkSmartPointer<vtkIdList> lineIds = vtkSmartPointer<vtkIdList>::New();
3307
3308     double ptA[3], ptB[3];
3309
3310     bool good;
3311
3312     vtkIdType j;
3313     vtkIdType _idA, _idB;
3314
3315     std::map<std::size_t, std::set<Conn, ConnCmp>> _polyConns;
3316
3317     vtkSmartPointer<vtkIdList> line = vtkSmartPointer<vtkIdList>::New();
3318
3319     for (const auto &ids : indexedPolys) {
3320         for (vtkIdType idA : ids) {
3321             pts->GetPoint(idA, ptA);
3322
3323             kdTree->FindClosestNPoints(n, ptA, foundPts);
3324
3325             numPts = foundPts->GetNumberOfIds();
3326
3327             for (i = 0; i < numPts; i++) {
3328                 idB = foundPts->GetId(i);
3329
3330                 auto srcA = sources.at(idA),
3331                     srcB = sources.at(idB);
3332
3333                 if (srcA == srcB) {
3334                     continue;
3335                 }
3336
3337                 pts->GetPoint(idB, ptB);
3338
3339                 good = true;
3340
3341                 if (bspTree->IntersectWithLine(ptA, ptB, 1e-5, nullptr, lineIds) == 1) {
3342                     for (j = 0; j < lineIds->GetNumberOfIds(); j++) {
3343                         lines->GetCellPoints(lineIds->GetId(j), line);
3344
3345                         _idA = line->GetId(0);
3346                         _idB = line->GetId(1);
3347
3348                         if (_idA != idA && _idA != idB
3349                             && _idB != idA && _idB != idB) {
3350
3351                             good = false;
3352                             break;
3353                         }
3354                     }
3355                 }
3356
3357                 if (good) {
3358                     double d = vtkMath::Distance2BetweenPoints(ptA, ptB);
3359
3360                     _polyConns[srcA].emplace(d, idA, idB);
3361                     _polyConns[srcB].emplace(d, idB, idA);
3362                 }
3363             }
3364         }
3365     }
3366
3367     decltype(_polyConns)::const_iterator itr;
3368
3369     for (itr = _polyConns.begin(); itr != _polyConns.end(); ++itr) {
3370         auto &_conns = itr->second;
3371
3372         ConnsType conns(_conns.begin(), _conns.end());
3373         std::sort(conns.begin(), conns.end());
3374
3375         polyConns[itr->first].swap(conns);
3376
3377     }
3378
3379     return true;
3380 }