<commit_before>/*
Copyright 2012-2019 <NAME>

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

#include <map>
#include <deque>
#include <vector>
#include <set>
#include <algorithm>
#include <cmath>
#include <functional>
#include <queue>
#include <sstream>

#include <vtkInformation.h>
#include <vtkInformationVector.h>
#include <vtkDemandDrivenPipeline.h>
#include <vtkObjectFactory.h>
#include <vtkPolyDataAlgorithm.h>
#include <vtkCellData.h>
#include <vtkPointData.h>
#include <vtkMath.h>
#include <vtkIdList.h>
#include <vtkCell.h>
#include <vtkAppendPolyData.h>
#include <vtkKdTreePointLocator.h>
#include <vtkCleanPolyData.h>
#include <vtkPolyDataConnectivityFilter.h>

#include "vtkPolyDataBooleanFilter.h"
#include "vtkPolyDataContactFilter.h"

#include "Utilities.h"

#include "Merger.h"
#include "Decomposer.h"

#ifdef DEBUG
#include <chrono>
#include <numeric>
#include <iterator>
#endif

// #include <csignal>

vtkStandardNewMacro(vtkPolyDataBooleanFilter);

vtkPolyDataBooleanFilter::vtkPolyDataBooleanFilter () {

    SetNumberOfInputPorts(2);
    SetNumberOfOutputPorts(2);

    timePdA = 0;
    timePdB = 0;

    contLines = vtkPolyData::New();

    modPdA = vtkPolyData::New();
    modPdB = vtkPolyData::New();

    cellDataA = vtkCellData::New();
    cellDataB = vtkCellData::New();

    cellIdsA = vtkIntArray::New();
    cellIdsB = vtkIntArray::New();

    OperMode = OPER_UNION;

    MergeRegs = false;
    DecPolys = true;

}

vtkPolyDataBooleanFilter::~vtkPolyDataBooleanFilter () {

    cellIdsA->Delete();
    cellIdsB->Delete();

    cellDataB->Delete();
    cellDataA->Delete();

    modPdB->Delete();
    modPdA->Delete();

    contLines->Delete();

}

int vtkPolyDataBooleanFilter::ProcessRequest(vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector) {

    if (request->Has(vtkDemandDrivenPipeline::REQUEST_DATA())) {

        vtkInformation *inInfoA = inputVector[0]->GetInformationObject(0);
        vtkInformation *inInfoB = inputVector[1]->GetInformationObject(0);

        vtkPolyData *pdA = vtkPolyData::SafeDownCast(inInfoA->Get(vtkDataObject::DATA_OBJECT()));
        vtkPolyData *pdB = vtkPolyData::SafeDownCast(inInfoB->Get(vtkDataObject::DATA_OBJECT()));

        vtkInformation *outInfoA = outputVector->GetInformationObject(0);
        vtkInformation *outInfoB = outputVector->GetInformationObject(1);

        resultA = vtkPolyData::SafeDownCast(outInfoA->Get(vtkDataObject::DATA_OBJECT()));
        resultB = vtkPolyData::SafeDownCast(outInfoB->Get(vtkDataObject::DATA_OBJECT()));

#ifdef DEBUG
        using clock = std::chrono::steady_clock;
        std::vector<clock::duration> times;
        clock::time_point start;
#endif

        if (pdA->GetMTime() > timePdA || pdB->GetMTime() > timePdB) {

            // eventuell vorhandene regionen vereinen

            vtkCleanPolyData *cleanA = vtkCleanPolyData::New();
            cleanA->SetOutputPointsPrecision(DOUBLE_PRECISION);
            cleanA->SetTolerance(1e-6);
            cleanA->SetInputData(pdA);
            cleanA->Update();

            vtkCleanPolyData *cleanB = vtkCleanPolyData::New();
            cleanB->SetOutputPointsPrecision(DOUBLE_PRECISION);
            cleanB->SetTolerance(1e-6);
            cleanB->SetInputData(pdB);
            cleanB->Update();

#ifdef DEBUG
            std::cout << "Exporting modPdA.vtk" << std::endl;
            WriteVTK("modPdA.vtk", cleanA->GetOutput());

            std::cout << "Exporting modPdB.vtk" << std::endl;
            WriteVTK("modPdB.vtk", cleanB->GetOutput());
#endif

            // CellData sichern

            cellDataA->DeepCopy(cleanA->GetOutput()->GetCellData());
            cellDataB->DeepCopy(cleanB->GetOutput()->GetCellData());

            // ermittelt kontaktstellen

#ifdef DEBUG

            start = clock::now();
#endif

            vtkPolyDataContactFilter *cl = vtkPolyDataContactFilter::New();
            cl->SetInputConnection(0, cleanA->GetOutputPort());
            cl->SetInputConnection(1, cleanB->GetOutputPort());
            cl->Update();

#ifdef DEBUG
            times.push_back(clock::now()-start);
#endif

            contLines->DeepCopy(cl->GetOutput());

#ifdef DEBUG
            std::cout << "Exporting contLines.vtk" << std::endl;
            WriteVTK("contLines.vtk", contLines);

            std::cout << "Exporting modPdA_1.vtk" << std::endl;
            WriteVTK("modPdA_1.vtk", cl->GetOutput(1));

            std::cout << "Exporting modPdB_1.vtk" << std::endl;
            WriteVTK("modPdB_1.vtk", cl->GetOutput(2));
#endif

            modPdA->DeepCopy(cl->GetOutput(1));
            modPdB->DeepCopy(cl->GetOutput(2));

            if (contLines->GetNumberOfCells() == 0) {
                vtkErrorMacro("Inputs have no contact.");

                return 1;
            }

            // in den CellDatas steht drin, welche polygone einander schneiden

            vtkIntArray *contsA = vtkIntArray::SafeDownCast(contLines->GetCellData()->GetScalars("cA"));
            vtkIntArray *contsB = vtkIntArray::SafeDownCast(contLines->GetCellData()->GetScalars("cB"));

            vtkIntArray *sourcesA = vtkIntArray::SafeDownCast(contLines->GetCellData()->GetScalars("sourcesA"));
            vtkIntArray *sourcesB = vtkIntArray::SafeDownCast(contLines->GetCellData()->GetScalars("sourcesB"));

            // gültigkeit des schnitts prüfen

            bool valid = true;

            vtkIdList *cells = vtkIdList::New();

            for (int i = 0; i < contLines->GetNumberOfPoints() && valid; i++) {
                cells->Reset();
                contLines->GetPointCells(i, cells);

                // der schnitt endet abrupt
                if (cells->GetNumberOfIds() == 1) {
                    valid = false;
                    break;
                }

                vtkIdList *line = vtkIdList::New();

                std::map<int, std::set<int> > countsA, countsB;

                for (int j = 0; j < cells->GetNumberOfIds(); j++) {
                    contLines->GetCellPoints(cells->GetId(j), line);

                    int target = line->GetId(0) == i ? line->GetId(1) : line->GetId(0);

                    int indA = contsA->GetValue(cells->GetId(j));
                    int indB = contsB->GetValue(cells->GetId(j));

                    // sind am punkt sind mehr als zwei linien beteiligt?

                    countsA[indA].insert(target);

                    if (countsA[indA].size() > 2) {
                        valid = false;
                        break;
                    }

                    countsB[indB].insert(target);

                    if (countsB[indB].size() > 2) {
                        valid = false;
                        break;
                    }

                }

                line->Delete();
            }

            cells->Delete();

            if (!valid) {
                vtkErrorMacro("Contact is ambiguous or incomplete.");

                return 1;
            }

            // sichert die OrigCellIds

            vtkIntArray *origCellIdsA = vtkIntArray::SafeDownCast(modPdA->GetCellData()->GetScalars("OrigCellIds"));
            vtkIntArray *origCellIdsB = vtkIntArray::SafeDownCast(modPdB->GetCellData()->GetScalars("OrigCellIds"));

            cellIdsA->DeepCopy(origCellIdsA);
            cellIdsB->DeepCopy(origCellIdsB);

            for (int i = 0; i < modPdA->GetNumberOfCells(); i++) {
                origCellIdsA->SetValue(i, i);
            }

            for (int i = 0; i < modPdB->GetNumberOfCells(); i++) {
                origCellIdsB->SetValue(i, i);
            }


#ifdef DEBUG
            start = clock::now();
#endif

            GetPolyStrips(modPdA, contsA, sourcesA, polyStripsA);
            GetPolyStrips(modPdB, contsB, sourcesB, polyStripsB);

#ifdef DEBUG
            times.push_back(clock::now()-start);
#endif

            // löst ein sehr spezielles problem

#ifdef DEBUG
            start = clock::now();
#endif

            _Test(modPdA, polyStripsA);
            _Test(modPdB, polyStripsB);

#ifdef DEBUG
            times.push_back(clock::now()-start);
#endif

            // trennt die polygone an den linien

#ifdef DEBUG
            start = clock::now();
#endif

            CutCells(modPdA, polyStripsA);
            CutCells(modPdB, polyStripsB);

#ifdef DEBUG
            times.push_back(clock::now()-start);
#endif

#ifdef DEBUG
            std::cout << "Exporting modPdA_2.vtk" << std::endl;
            WriteVTK("modPdA_2.vtk", modPdA);

            std::cout << "Exporting modPdB_2.vtk" << std::endl;
            WriteVTK("modPdB_2.vtk", modPdB);
#endif

#ifdef DEBUG
            start = clock::now();
#endif

            RestoreOrigPoints(modPdA, polyStripsA);
            RestoreOrigPoints(modPdB, polyStripsB);

#ifdef DEBUG
            times.push_back(clock::now()-start);
#endif

#ifdef DEBUG
            std::cout << "Exporting modPdA_3.vtk" << std::endl;
            WriteVTK("modPdA_3.vtk", modPdA);

            std::cout << "Exporting modPdB_3.vtk" << std::endl;
            WriteVTK("modPdB_3.vtk", modPdB);
#endif

#ifdef DEBUG
            start = clock::now();
#endif

            ResolveOverlaps(modPdA, contsA, polyStripsA);
            ResolveOverlaps(modPdB, contsB, polyStripsB);

#ifdef DEBUG
            times.push_back(clock::now()-start);
#endif

#ifdef DEBUG
            std::cout << "Exporting modPdA_4.vtk" << std::endl;
            WriteVTK("modPdA_4.vtk", modPdA);

            std::cout << "Exporting modPdB_4.vtk" << std::endl;
            WriteVTK("modPdB_4.vtk", modPdB);
#endif

#ifdef DEBUG
            start = clock::now();
#endif

            AddAdjacentPoints(modPdA, contsA, polyStripsA);
            AddAdjacentPoints(modPdB, contsB, polyStripsB);

#ifdef DEBUG
            times.push_back(clock::now()-start);
#endif

#ifdef DEBUG
            std::cout << "Exporting modPdA_5.vtk" << std::endl;
            WriteVTK("modPdA_5.vtk", modPdA);

            std::cout << "Exporting modPdB_5.vtk" << std::endl;
            WriteVTK("modPdB_5.vtk", modPdB);
#endif

#ifdef DEBUG
            start = clock::now();
#endif

            DisjoinPolys(modPdA, polyStripsA);
            DisjoinPolys(modPdB, polyStripsB);

#ifdef DEBUG
            times.push_back(clock::now()-start);
#endif

#ifdef DEBUG
            std::cout << "Exporting modPdA_6.vtk" << std::endl;
            WriteVTK("modPdA_6.vtk", modPdA);

            std::cout << "Exporting modPdB_6.vtk" << std::endl;
            WriteVTK("modPdB_6.vtk", modPdB);
#endif

#ifdef DEBUG
            start = clock::now();
#endif

            MergePoints(modPdA, polyStripsA);
            MergePoints(modPdB, polyStripsB);

#ifdef DEBUG
            times.push_back(clock::now()-start);
#endif

#ifdef DEBUG
            std::cout << "Exporting modPdA_7.vtk" << std::endl;
            WriteVTK("modPdA_7.vtk", modPdA);

            std::cout << "Exporting modPdB_7.vtk" << std::endl;
            WriteVTK("modPdB_7.vtk", modPdB);
#endif

            involvedA.clear();
            involvedB.clear();

            int numLines = contLines->GetNumberOfCells();
            for (int i = 0; i < numLines; i++) {
                involvedA.insert(contsA->GetValue(i));
                involvedB.insert(contsB->GetValue(i));
            }

            relsA.clear();
            relsB.clear();

            // aufräumen

            cl->Delete();
            cleanB->Delete();
            cleanA->Delete();

            timePdA = pdA->GetMTime();
            timePdB = pdB->GetMTime();

        }

#ifdef DEBUG
        times.push_back(clock::now()-start);
#endif

        DecPolys_(modPdA, involvedA, relsA);
        DecPolys_(modPdB, involvedB, relsB);

#ifdef DEBUG
        times.push_back(clock::now()-start);
#endif

#ifdef DEBUG
        std::cout << "Exporting modPdA_8.vtk" << std::endl;
        WriteVTK("modPdA_8.vtk", modPdA);

        std::cout << "Exporting modPdB_8.vtk" << std::endl;
        WriteVTK("modPdB_8.vtk", modPdB);
#endif

#ifdef DEBUG
        start = clock::now();
#endif

        if (MergeRegs) {
            MergeRegions();
        } else {
            CombineRegions();
        }

#ifdef DEBUG
        times.push_back(clock::now()-start);
#endif


#ifdef DEBUG
        double sum = std::chrono::duration_cast<std::chrono::duration<double>>(std::accumulate(times.begin(), times.end(), clock::duration())).count();

        std::vector<clock::duration>::const_iterator itr;
        for (itr = times.begin(); itr != times.end(); itr++) {
            double time = std::chrono::duration_cast<std::chrono::duration<double>>(*itr).count();

            std::cout << "Time " << (itr-times.begin())
                << ": " << time << "s (" << (time/sum) << "%)"
                << std::endl;
        }
#endif

    }

    return 1;

}

void vtkPolyDataBooleanFilter::GetStripPoints (vtkPolyData *pd, vtkIntArray *sources, PStrips &pStrips, IdsType &lines) {

#ifdef DEBUG
    std::cout << "GetStripPoints()" << std::endl;
#endif

    StripPtsType &pts = pStrips.pts;
    IdsType &poly = pStrips.poly;

    int numPts = poly.size();

    double tol = 1e-5;

    IdsType::const_iterator itr;

    for (itr = lines.begin(); itr != lines.end(); itr++) {

        vtkIdList *linePts = vtkIdList::New();
        contLines->GetCellPoints(*itr, linePts);

        // diese punkte durchlaufen

        for (int i = 0; i < 2; i++) {

            int realInd = linePts->GetId(i);

            if (pts.count(realInd) == 0) {
                // lage analysieren

                pts[realInd].ind = realInd;

                // die koordinaten
                double pt[3];
                contLines->GetPoint(realInd, pt);

                Cpy(pts[realInd].pt, pt, 3);

                int src = sources->GetComponent(*itr, i);

                pts[realInd].src = src;

                double lastD = DBL_MAX;

                // jetzt muss man die kanten durchlaufen
                for (int j = 0; j < numPts; j++) {

                    if (src > -1 && !(j == src || (j+1)%numPts == src)) {
                        continue;
                    }

                    int indA, indB;
                    double a[3], b[3];

                    indA = poly[j];
                    indB = poly[(j+1)%numPts];

                    pd->GetPoint(indA, a);
                    pd->GetPoint(indB, b);

                    double sA[3], sB[3];
                    vtkMath::Subtract(a, pt, sA);
                    vtkMath::Subtract(b, pt, sB);

                    // richtungsvektor und länge der kante

                    double u[3];
                    vtkMath::Subtract(b, a, u);

                    double n = vtkMath::Norm(u);

                    // d und t zur kante

                    double v[3];
                    vtkMath::Subtract(pt, a, v);

                    double t = vtkMath::Dot(v, u)/(n*n);

                    double w[3];
                    vtkMath::Cross(v, u, w);

                    double d = vtkMath::Norm(w)/n;

                    if (d < tol && d < lastD && t > -tol && t < 1+tol) {
                        // liegt im toleranzbereich der kante

                        pts[realInd].edge[0] = indA;
                        pts[realInd].edge[1] = indB;

                        pts[realInd].t = std::min(1., std::max(0., t));

                        if (vtkMath::Norm(sA) < tol) {
                            Cpy(pts[realInd].captPt, a, 3);
                            pts[realInd].capt = CAPT_A;

                            break;

                        } else if (vtkMath::Norm(sB) < tol) {
                            Cpy(pts[realInd].captPt, b, 3);
                            pts[realInd].capt = CAPT_B;

                            break;

                        } else {
                            // u ist nicht normiert
                            vtkMath::MultiplyScalar(u, t);

                            double x[3];
                            vtkMath::Add(a, u, x);

                            // projektion
                            Cpy(pts[realInd].captPt, x, 3);

                            pts[realInd].capt = CAPT_EDGE;

                            lastD = d;
                        }

                    }

                }

                if (src != NO_USE && pts[realInd].edge[0] == NO_USE) {
                    pts[realInd].catched = false;
                }

            }
        }

        linePts->Delete();

    }

    StripPtsType::iterator itr2;

    for (itr2 = pts.begin(); itr2 != pts.end(); ++itr2) {
        StripPt &sp = itr2->second;

        if (sp.capt != CAPT_NOT) {
            if (sp.capt == CAPT_B) {
                sp.t = 0;

                sp.edge[0] = sp.edge[1];

                int i = std::find(poly.begin(), poly.end(), sp.edge[0])-poly.begin();

                sp.edge[1] = poly[(i+1)%numPts];

                sp.capt = CAPT_A;

            }

            // für den schnitt werden die eingerasteten koordinaten verwendet

            Cpy(sp.cutPt, sp.captPt, 3);
        } else {

            Cpy(sp.cutPt, sp.pt, 3);
        }

        sp.history.push_back({sp.edge[0], sp.edge[1]});

    }

#ifdef DEBUG
    std::cout << "pts: " << std::endl;
    for (itr2 = pts.begin(); itr2 != pts.end(); ++itr2) {
        std::cout << itr2->first << ": " << itr2->second << std::endl;
    }
#endif

}

void vtkPolyDataBooleanFilter::GetPolyStrips (vtkPolyData *pd, vtkIntArray *conts, vtkIntArray *sources, PolyStripsType &polyStrips) {
#ifdef DEBUG
    std::cout << "GetPolyStrips()" << std::endl;
#endif

    polyStrips.clear();

    std::map<int, IdsType> polyLines;

    for (int i = 0; i < conts->GetNumberOfTuples(); i++) {
        int poly = conts->GetValue(i);
        polyLines[poly].push_back(i);
    }

    std::vector<std::reference_wrapper<StripPt>> notCatched;

    std::map<int, IdsType>::iterator itr;

    for (itr = polyLines.begin(); itr != polyLines.end(); ++itr) {
        IdsType &lines = itr->second;
        RemoveDuplicates(lines);

        PStrips &pStrips = polyStrips[itr->first];

        vtkIdList *polyPts = vtkIdList::New();
        pd->GetCellPoints(itr->first, polyPts);

        int numPts = polyPts->GetNumberOfIds();

        for (int i = 0; i < numPts; i++) {
            pStrips.poly.push_back(polyPts->GetId(i));
        }

        ComputeNormal(pd->GetPoints(), pStrips.n, polyPts);

        GetStripPoints(pd, sources, pStrips, lines);

        for (auto &sp : pStrips.pts) {
            sp.second.polyId = itr->first;

            if (!sp.second.catched) {
                notCatched.push_back(sp.second);
            }
        }

        polyPts->Delete();
    }

    auto Next = [](const IdsType &ids, int id) -> int {
        IdsType::const_iterator itr;

        itr = std::find(ids.begin(), ids.end(), id);

        if (++itr == ids.end()) {
            itr = ids.begin();
        }

        return *itr;
    };

    for (StripPt &sp : notCatched) {
        for (itr = polyLines.begin(); itr != polyLines.end(); ++itr) {
            const PStrips &pStrips = polyStrips[itr->first];

            try {
                const StripPt &corr = pStrips.pts.at(sp.ind);

                if (&corr != &sp) {
                    if (corr.capt == CAPT_A) {
                        sp.capt = CAPT_A;
                        sp.edge[0] = corr.edge[0];
                        sp.edge[1] = Next(polyStrips[sp.polyId].poly, sp.edge[0]);

                        sp.t = 0;

                        Cpy(sp.captPt, corr.captPt, 3);
                        Cpy(sp.cutPt, sp.captPt, 3);

                        sp.history.push_back({sp.edge[0], sp.edge[1]});

                        sp.catched = true;

                    }
                }
            } catch (...) {}

        }

        assert(sp.catched);

    }

    for (itr = polyLines.begin(); itr != polyLines.end(); ++itr) {
        IdsType &lines = itr->second;

        PStrips &pStrips = polyStrips[itr->first];

        StripType strip;

        // zusammensetzen

        std::deque<int> _lines(lines.begin(), lines.end());

        int i = 0;

        while (_lines.size() > 0) {
            vtkIdList *linePts = vtkIdList::New();
            contLines->GetCellPoints(_lines[i], linePts);

            int indA = linePts->GetId(0);
            int indB = linePts->GetId(1);

            if (strip.empty()) {
                strip.push_back(StripPtR(indA));
                strip.push_back(StripPtR(indB));

                _lines.erase(_lines.begin());
            } else {
                StripPt &start = pStrips.pts[strip.front().ind],
                    &end = pStrips.pts[strip.back().ind];

                if (end.capt == CAPT_NOT && end.ind == indA) {
                    strip.push_back(StripPtR(indB));
                    _lines.erase(_lines.begin()+i);
                    i = 0;

                } else if (end.capt == CAPT_NOT && end.ind == indB) {
                    strip.push_back(StripPtR(indA));
                    _lines.erase(_lines.begin()+i);
                    i = 0;

                } else if (start.capt == CAPT_NOT && start.ind == indA) {
                    strip.push_front(StripPtR(indB));
                    _lines.erase(_lines.begin()+i);
                    i = 0;

                } else if (start.capt == CAPT_NOT && start.ind == indB) {
                    strip.push_front(StripPtR(indA));
                    _lines.erase(_lines.begin()+i);
                    i = 0;

                } else {
                    i++;

                }
            }

            StripPt &_start = pStrips.pts[strip.front().ind],
                &_end = pStrips.pts[strip.back().ind];

            if ((_start.capt != CAPT_NOT && _end.capt != CAPT_NOT)
                || (_start.ind == _end.ind)
                || (i == _lines.size())) {

                // einen neuen strip anlegen
                pStrips.strips.push_back(strip);
                strip.clear();

                i = 0;
            }

            linePts->Delete();

        }

        CompleteStrips(pStrips);

    }

}

void vtkPolyDataBooleanFilter::RemoveDuplicates (IdsType &lines) {

    IdsType unique;

    int i, j;

    // die indexe der enden auf übereinstimmung prüfen

    vtkIdList *linePtsA = vtkIdList::New();
    vtkIdList *linePtsB = vtkIdList::New();

    int numLines = lines.size();

    for (i = 0; i < numLines-1; i++) {
        j = i+1;

        contLines->GetCellPoints(lines[i], linePtsA);

        while (j < lines.size()) {
            contLines->GetCellPoints(lines[j], linePtsB);

            if ((linePtsA->GetId(0) == linePtsB->GetId(0) && linePtsA->GetId(1) == linePtsB->GetId(1)) ||
                (linePtsA->GetId(0) == linePtsB->GetId(1) && linePtsA->GetId(1) == linePtsB->GetId(0))) {
                // stimmen überein
                break;
            }

            j++;
        }

        if (j == numLines) {
            // keine vorzeitige unterbrechung der while-schleife
            unique.push_back(lines[i]);
        }
    }

    unique.push_back(lines.back());

    linePtsA->Delete();
    linePtsB->Delete();

    lines.swap(unique);

}

void vtkPolyDataBooleanFilter::CompleteStrips (PStrips &pStrips) {
    StripsType::iterator itr;

    for (itr = pStrips.strips.begin(); itr != pStrips.strips.end(); ++itr) {
        StripPt &start = pStrips.pts[itr->front().ind],
            &end = pStrips.pts[itr->back().ind];

        if (start.ind != end.ind) {
            if (start.capt == CAPT_NOT) {
                StripType s(itr->rbegin(), itr->rend()-1);
                itr->insert(itr->begin(), s.begin(), s.end());

            } else if (end.capt == CAPT_NOT) {
                StripType s(itr->rbegin()+1, itr->rend());
                itr->insert(itr->end(), s.begin(), s.end());

            }
        }
    }
}

bool vtkPolyDataBooleanFilter::HasArea (StripType &strip) {
    bool area = true;

    int n = strip.size();
    if (n%2 == 1) {
        for (int i = 0; i < (n-1)/2; i++) {
            area = strip[i].ind != strip[n-i-1].ind;
        }
    }

    return area;
}

void vtkPolyDataBooleanFilter::CutCells (vtkPolyData *pd, PolyStripsType &polyStrips) {
#ifdef DEBUG
    std::cout << "CutCells()" << std::endl;
#endif

    vtkPoints *pdPts = pd->GetPoints();

    vtkIntArray *origCellIds = vtkIntArray::SafeDownCast(pd->GetCellData()->GetScalars("OrigCellIds"));

    PolyStripsType::iterator itr;

    for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) {

        int polyInd = itr->first;
        PStrips &pStrips = itr->second;

#ifdef DEBUG
        if (polyInd != 204) {
            //continue;
        }
#endif

        StripsType &strips = pStrips.strips;
        StripPtsType &pts = pStrips.pts;

        IdsType &poly = pStrips.poly;

        int origId = origCellIds->GetValue(polyInd);

#ifdef DEBUG
        IdsType::iterator itr_;
        std::cout << "polyInd=" << polyInd << ", poly=[";
        for (itr_ = poly.begin(); itr_ != poly.end(); ++itr_) {
            std::cout << *itr_ << ", ";
        }
        std::cout << "]" << std::endl;
#endif

        int numPts = poly.size();

        std::map<int, RefsType> edges;

        StripsType::iterator itr2;
        StripType::iterator itr3;

        // holes sammeln

        HolesType holes;

        for (itr2 = strips.begin(); itr2 != strips.end(); ++itr2) {
            StripType &s = *itr2;
            if (pts[s.front().ind].capt == CAPT_NOT && pts[s.back().ind].capt == CAPT_NOT) {
                IdsType hole;

                for (auto& sp : s) {
                    hole.push_back(pdPts->InsertNextPoint(pts[sp.ind].pt));
                }

                // anfang und ende sind ja gleich
                hole.pop_back();

                holes.push_back(std::move(hole));

            }
        }

        // holes löschen
        strips.erase(std::remove_if(strips.begin(), strips.end(), [&](const StripType &s) {
            return pts[s.front().ind].capt == CAPT_NOT && pts[s.back().ind].capt == CAPT_NOT; }), strips.end());

        for (itr2 = strips.begin(); itr2 != strips.end(); ++itr2) {
            StripType &strip = *itr2;

#ifdef DEBUG
            std::cout << "strip [";
            for (auto &s : strip) {
                std::cout << s.ind << ", ";
            }
            std::cout << "]" << std::endl;
#endif

            // init
            if (pts[strip.front().ind].edge[0] == pts[strip.back().ind].edge[0]
                && strip.front().ind != strip.back().ind
                && pts[strip.front().ind].t > pts[strip.back().ind].t) {

                std::reverse(strip.begin(), strip.end());
            }

            StripPt &start = pts[strip.front().ind],
                &end = pts[strip.back().ind];

            strip.front().side = SIDE_START;
            strip.back().side = SIDE_END;

            strip.front().ref = start.edge[0];
            strip.back().ref = end.edge[0];

            int ind = itr2-strips.begin();

            strip.front().strip = strip.back().strip = ind;

            for (itr3 = strip.begin(); itr3 != strip.end(); ++itr3) {
                StripPt &sp = pts[itr3->ind];

                itr3->desc[0] = pdPts->InsertNextPoint(sp.cutPt);
                itr3->desc[1] = pdPts->InsertNextPoint(sp.cutPt);

#ifdef DEBUG
                std::cout << sp << " => " << *itr3 << std::endl;
#endif

            }

            // ordnet zu

            edges[start.edge[0]].push_back(std::ref(strip.front()));
            edges[end.edge[0]].push_back(std::ref(strip.back()));
        }

        // sortiert die punkte auf den kanten

        std::map<int, RefsType>::iterator itr4;

        IdsType::iterator itr5;
        StripType::reverse_iterator itr6;

        RefsType::iterator itr7;

        std::vector<IdsType>::iterator itr8;


        for (itr4 = edges.begin(); itr4 != edges.end(); ++itr4) {
            RefsType &edge = itr4->second;

#ifdef DEBUG
            std::cout << "edge (" << itr4->first << ", _)" << std::endl;
#endif

            std::sort(edge.begin(), edge.end(), [&](const StripPtR &a, const StripPtR &b) {
                StripPt &a_ = pts[a.ind],
                    &b_ = pts[b.ind];

#ifdef DEBUG
                std::cout << "a_: " << a_ << " -> strip " << a.strip << std::endl;
                std::cout << "b_: " << b_ << " -> strip " << b.strip << std::endl;
#endif

                if (a_.ind == b_.ind) {
                    // strips beginnen im gleichen punkt

                    if (a.strip != b.strip) {
                        // gehören nicht dem gleichen strip an

                        StripType &stripA = strips[a.strip],
                            &stripB = strips[b.strip];

                        // andere enden ermitteln

                        StripPtR &eA = (&a == &(stripA.front())) ? stripA.back() : stripA.front(),
                            &eB = (&b == &(stripB.front())) ? stripB.back() : stripB.front();

                        StripPt &eA_ = pts[eA.ind],
                            &eB_ = pts[eB.ind];

#ifdef DEBUG
                        std::cout << "eA_: " << eA_ << std::endl;
                        std::cout << "eB_: " << eA_ << std::endl;
#endif

                        if (eA_.ind != eB_.ind) {
                            int i = std::find(poly.begin(), poly.end(), itr4->first)-poly.begin();

                            int iA = std::find(poly.begin(), poly.end(), eA_.edge[0])-poly.begin(),
                                iB = std::find(poly.begin(), poly.end(), eB_.edge[0])-poly.begin();

                            double dA = Mod(iA-i, numPts)+eA_.t,
                                dB = Mod(iB-i, numPts)+eB_.t;

                            if (i == iA && a_.t > eA_.t) {
                               dA += numPts;
                            }

                            if (i == iB && b_.t > eB_.t) {
                               dB += numPts;
                            }

#ifdef DEBUG
                            std::cout << "dA=" << dA << ", dB=" << dB << std::endl;
#endif

                            return dB < dA;
                        } else {
                            RefsType poly_;

                            if (a.side == SIDE_START) {
                                poly_.insert(poly_.end(), stripA.begin(), stripA.end());
                            } else {
                                poly_.insert(poly_.end(), stripA.rbegin(), stripA.rend());
                            }

                            if (b.side == SIDE_START) {
                                poly_.insert(poly_.end(), stripB.rbegin()+1, stripB.rend()-1);
                            } else {
                                poly_.insert(poly_.end(), stripB.begin()+1, stripB.end()-1);
                            }

                            int num = poly_.size();

                            vtkPoints *pts_ = vtkPoints::New();
                            pts_->SetNumberOfPoints(num);

                            for (itr7 = poly_.begin(); itr7 != poly_.end(); ++itr7) {
                                StripPtR &sp = *itr7;
                                int i = itr7-poly_.begin();

                                pts_->SetPoint(i, pts[sp.ind].cutPt);
                            }

                            double n[3];
                            ComputeNormal(pts_, n);

                            pts_->Delete();

                            double ang = vtkMath::Dot(pStrips.n, n);

#ifdef DEBUG
                            std::cout << "ang=" << ang*180/PI << std::endl;
#endif

                            return ang < .999999;

                        }
                    } else {
                        // gleicher strip

                        StripType &strip = strips[a.strip];

                        if (HasArea(strip)) {
                            RefsType poly_(strip.begin(), strip.end()-1);

                            int num = poly_.size();

                            vtkPoints *pts_ = vtkPoints::New();
                            pts_->SetNumberOfPoints(num);

                            for (itr7 = poly_.begin(); itr7 != poly_.end(); ++itr7) {
                                StripPtR &sp = *itr7;
                                int i = itr7-poly_.begin();

                                pts_->SetPoint(i, pts[sp.ind].cutPt);
                            }

                            double n[3];
                            ComputeNormal(pts_, n);

                            pts_->Delete();

                            double ang = vtkMath::Dot(pStrips.n, n);

                            if (ang > .999999) {
                                std::reverse(strip.begin(), strip.end());
                                return true;
                            } else {
                                return false;
                            }

                        } else {
                            // reihenfolge von a und b bereits richtig
                            return false;
                        }
                    }

                } else {
                    return a_.t < b_.t;
                }
            });

#ifdef DEBUG
            std::cout << "after sort" << std::endl;
            for (itr7 = edge.begin(); itr7 != edge.end(); ++itr7) {
                std::cout << itr7->get().strip << ": " << pts[itr7->get().ind] << std::endl;
            }
#endif

        }

        // baut die strips ein

        std::deque<IdsType> polys;
        polys.push_back(pStrips.poly);

        int i = 0;

        std::deque<IdsType>::iterator itr9;

        for (itr2 = strips.begin(); itr2 != strips.end(); ++itr2) {
            StripType &strip = *itr2;

            StripPtR &start = strip.front(),
                &end = strip.back();

#ifdef DEBUG
            std::cout << "strip " << start.strip
                << " refs=[" << start.ref << ", " << end.ref << "]"
                << std::endl;
#endif

            int cycle = 0;

            while (true) {

                if (cycle == polys.size()) {
                    break;
                }

                IdsType p = polys.front();
                polys.pop_front();

                std::vector<IdsType> divided(2);

                if (std::find(p.begin(), p.end(), start.ref) != p.end()) {
                    if (start.ref == end.ref) {
                        for (itr5 = p.begin(); itr5 != p.end(); ++itr5) {
                            divided[0].push_back(*itr5);

#ifdef DEBUG
                            std::cout << "add d[0] " << *itr5 << std::endl;
#endif

                            if (*itr5 == start.ref) {
                                for (itr3 = strip.begin(); itr3 != strip.end(); ++itr3) {
                                    divided[0].push_back(itr3->desc[0]);

#ifdef DEBUG
                                    std::cout << "add d[0] " << itr3->desc[0] << std::endl;
#endif

                                }
                            }
                        }

                        // strip selbst ist ein polygon

                        for (itr6 = strip.rbegin(); itr6 != strip.rend(); ++itr6) {
                            divided[1].push_back(itr6->desc[1]);

#ifdef DEBUG
                            std::cout << "add d[1] " << itr6->desc[1] << std::endl;
#endif

                        }

                    } else {
                        int d = 0;
                        for (itr5 = p.begin(); itr5 != p.end(); ++itr5) {
                            divided[d].push_back(*itr5);

#ifdef DEBUG
                            std::cout << "add d[" << d << "] " << *itr5 << std::endl;
#endif

                            if (*itr5 == start.ref) {
                                for (itr3 = strip.begin(); itr3 != strip.end(); ++itr3) {
                                    divided[d].push_back(itr3->desc[0]);

#ifdef DEBUG
                                    std::cout << "add d[" << d << "] " << itr3->desc[0] << std::endl;
#endif

                                }
                                d++;
                            } else if (*itr5 == end.ref) {
                                for (itr6 = strip.rbegin(); itr6 != strip.rend(); ++itr6) {
                                    divided[d].push_back(itr6->desc[1]);

#ifdef DEBUG
                                    std::cout << "add d[" << d << "] " << itr6->desc[1] << std::endl;
#endif

                                }
                                d++;
                            }

                            d %= 2;
                        }
                    }
                }

                if (divided[1].size() > 0) {

                    // refs aktualisieren

                    int stripInd = itr2-strips.begin();

                    for (itr4 = edges.begin(); itr4 != edges.end(); ++itr4) {
                        RefsType &edge = itr4->second;

                        for (itr7 = edge.begin()+1; itr7 != edge.end(); ++itr7) {
                            StripPtR &sp = *itr7;

                            RefsType::reverse_iterator itr7_(itr7);

                            for (; itr7_ != edge.rend(); ++itr7_) {
                                // pt davor
                                StripPtR &pre = *itr7_;

                                if (pre.strip != sp.strip) {
                                    if (pre.strip <= stripInd) {
                                        // der davor wurde bereits verwendet

                                        if (pre.side == SIDE_END) {
                                            sp.ref = pre.desc[0];
                                        } else {
                                            sp.ref = pre.desc[1];
                                        }

                                        break;
                                    }
                                } else {
                                    // gehört dem gleichen strip an
                                    sp.ref = pre.ref;

                                    break;
                                }
                            }

                        }

                        // erstellt die history

                        auto _s = std::find_if(edge.begin(), edge.end(), [&start](const StripPtR &r) {
                            return &start == &r;
                        });

                        if (_s != edge.end()) {
                            for (itr7 = edge.begin(); itr7 != edge.end(); ++itr7) {
                                StripPtR &sp = *itr7;
                                StripPt &_sp = pts[sp.ind];
                                auto &history = _sp.history;

                                if (&sp != &start) {
                                    if (_sp.t < pts[start.ind].t) {
                                        history.push_back({history.back().f, start.desc[0]});
                                    } else {
                                        history.push_back({start.desc[1], history.back().g});
                                    }

                                }

                            }
                        }

                        auto _e = std::find_if(edge.begin(), edge.end(), [&end](const StripPtR &r) {
                            return &end == &r;
                        });

                        if (_e != edge.end()) {
                            for (itr7 = edge.begin(); itr7 != edge.end(); ++itr7) {
                                StripPtR &sp = *itr7;
                                StripPt &_sp = pts[sp.ind];
                                auto &history = _sp.history;

                                if (&sp != &end) {
                                    if (_sp.t < pts[end.ind].t) {
                                        history.push_back({history.back().f, end.desc[1]});
                                    } else {
                                        history.push_back({end.desc[0], history.back().g});
                                    }

                                }

                            }
                        }

                        // sonderfall
                        if (edge.size() > 1) {
                            StripPtR &a = edge.front(),
                                &b = *(edge.begin()+1);

                            if (a.ind == b.ind
                                && b.strip == stripInd
                                && pts[a.ind].capt == CAPT_A) {

                                if (b.side == SIDE_START) {
                                    a.ref = b.desc[0];
                                } else {
                                    a.ref = b.desc[1];
                                }

                            }
                        }
                    }

                    // doppelte punkte entfernen

                    for (itr8 = divided.begin(); itr8 != divided.end(); ++itr8) {
#ifdef DEBUG
                        std::cout << "d[" << itr8-divided.begin() << "]" << std::endl;
#endif

                        IdsType &div = *itr8;

                        IdsType _div;

                        int num = div.size();

                        for (int i = 0; i < num; i++) {
                            double ptA[3], ptB[3];
                            pd->GetPoint(div[i], ptA);
                            pd->GetPoint(div[(i+1)%num], ptB);

                            double d = GetD(ptA, ptB);

#ifdef DEBUG
                            std::cout << "edge (" << div[i] << ", "
                                      << div[(i+1)%num] << "), d=" << d
                                      << std::endl;
#endif

                            if (d > 1e-6) {
                                _div.push_back(div[i]);
                            } else {
#ifdef DEBUG
                                std::cout << "rm " << div[i] << std::endl;
#endif
                            }
                        }

                        div = std::move(_div);
                    }

                    // prüfen ob erstellte polygone gültig

                    if (divided[0].size() > 2) {
                        polys.push_back(divided[0]);
                    }

                    if (HasArea(strip) && divided[1].size() > 2) {
                        polys.push_back(divided[1]);
                    }

                    cycle = 0;

                    break;

                } else {
                    polys.push_back(p);

                    cycle++;
                }

            }

            i++;

            if (i == 4) {
                //break;
            }

        }

        // erzeugte polys hinzufügen

        IdsType descIds;

        for (itr9 = polys.begin(); itr9 != polys.end(); ++itr9) {
            IdsType &p = *itr9;

            int num = p.size();

            vtkIdList *cell = vtkIdList::New();
            cell->SetNumberOfIds(num);

            for (int i = 0; i < num; i++) {
                cell->SetId(i, p[i]);
            }

            descIds.push_back(pd->InsertNextCell(VTK_POLYGON, cell));

            origCellIds->InsertNextValue(origId);

            cell->Delete();
        }

        pd->DeleteCell(polyInd);

        // holes verarbeiten

        if (!holes.empty()) {
            _Wrapper w(pd, descIds, origId);

            for (auto& hole : holes) {
                w.Add(hole);
            }

            w.MergeAll();
        }

    }

    pd->RemoveDeletedCells();
    pd->BuildCells();

}

void vtkPolyDataBooleanFilter::_Test (vtkPolyData *pd, PolyStripsType &polyStrips) {

#ifdef DEBUG
    std::cout << "_Test()" << std::endl;
#endif

    struct Cmp {
        bool operator() (const StripPt &l, const StripPt &r) const {
            const int x1 = static_cast<int>(l.cutPt[0]*1e5),
                y1 = static_cast<int>(l.cutPt[1]*1e5),
                z1 = static_cast<int>(l.cutPt[2]*1e5),
                x2 = static_cast<int>(r.cutPt[0]*1e5),
                y2 = static_cast<int>(r.cutPt[1]*1e5),
                z2 = static_cast<int>(r.cutPt[2]*1e5);

            return std::tie(x1, y1, z1) < std::tie(x2, y2, z2);
        }
    };

    struct Cmp2 {
        bool operator() (const StripPt &l, const StripPt &r) const {
            return l.ind < r.ind;
        }
    };

    std::map<StripPt, std::set<std::reference_wrapper<StripPt>, Cmp2>, Cmp> test;

    PolyStripsType::iterator itr;
    StripPtsType::iterator itr2;

    for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) {
        PStrips &pStrips = itr->second;

        for (itr2 = pStrips.pts.begin(); itr2 != pStrips.pts.end(); ++itr2) {
            StripPt &sp = itr2->second;

            test[sp].insert(sp);
        }
    }

    vtkIdList *cells = vtkIdList::New();

    for (auto &s : test) {
        auto &pts = s.second;

        if (pts.size() > 1) {
            StripPt &a = *(pts.begin()),
                &b = *(std::next(pts.begin()));

            int indA = a.ind,
                indB = b.ind;

#ifdef DEBUG
            std::cout << "collapsing " << b.ind << " -> " << a.ind << std::endl;
#endif

            contLines->GetPointCells(b.ind, cells);

            for (int i = 0; i < cells->GetNumberOfIds(); i++) {
                contLines->ReplaceCellPoint(cells->GetId(i), indB, indA);

                // aktualisiert die links
                contLines->RemoveReferenceToCell(indB, cells->GetId(i));
                contLines->ResizeCellList(indA, 1);
                contLines->AddReferenceToCell(indA, cells->GetId(i));
            }

            std::set<Pair> pairs;

            if (a.capt == CAPT_A && b.capt == CAPT_A) {
                std::map<int, IdsType> shared;

                vtkIdList *lines = vtkIdList::New(),
                    *cell = vtkIdList::New();
                contLines->GetPointCells(indA, lines);

                for (int i = 0; i < lines->GetNumberOfIds(); i++) {
                    contLines->GetCellPoints(lines->GetId(i), cell);

                    int pA = cell->GetId(0),
                        pB = cell->GetId(1);

                    shared[pA == indA ? pB : pA].push_back(lines->GetId(i));
                }

                cell->Delete();
                lines->Delete();

                for (auto &s : shared) {
                    if (s.second.size() > 1) {
                        assert(s.second.size() == 2);

                        for (int l : s.second) {
                            contLines->DeleteCell(l);
                            contLines->RemoveCellReference(l);

                            pairs.insert({indA, s.first});
                            pairs.insert({s.first, indA});
                        }
                    }
                }

            }

            if (a.capt == CAPT_A && b.capt == CAPT_EDGE) {
                vtkIdList *lines = vtkIdList::New(),
                    *cell = vtkIdList::New();
                contLines->GetPointCells(indA, lines);

                for (int i = 0; i < lines->GetNumberOfIds(); i++) {
                    if (contLines->GetCellType(lines->GetId(i)) != VTK_EMPTY_CELL) {
                        contLines->GetCellPoints(lines->GetId(i), cell);

                        int pA = cell->GetId(0),
                            pB = cell->GetId(1);

                        if (pA == indA && pB == indA) {
                            contLines->DeleteCell(lines->GetId(i));
                            contLines->RemoveCellReference(lines->GetId(i));
                        }
                    }

                }

                cell->Delete();
                lines->Delete();
            }

            auto Fct = [&](PolyStripsType &polyStrips_) -> void {
                for (auto &ps : polyStrips_) {
                    StripsType &strips = ps.second.strips;
                    StripPtsType &pts = ps.second.pts;

                    if (pts.count(indB) == 1) {
                        if (pts.count(indA) == 0) {
                            if (a.capt == pts[indB].capt) {
                                pts[indA] = pts[indB];
                                pts[indA].ind = indA;
                                Cpy(pts[indA].pt, a.pt, 3);
                            } else if (a.capt == CAPT_A && pts[indB].capt == CAPT_EDGE) {
                                pts[indA] = a;
                            } else {
                                // hier brauche ich ein beispiel
                                assert(false);<commit_msg>Another bugfix in _Test(2)<commit_after>/*
Copyright 2012-2019 <NAME>

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

#include <map>
#include <deque>
#include <vector>
#include <set>
#include <algorithm>
#include <cmath>
#include <functional>
#include <queue>
#include <sstream>

#include <vtkInformation.h>
#include <vtkInformationVector.h>
#include <vtkDemandDrivenPipeline.h>
#include <vtkObjectFactory.h>
#include <vtkPolyDataAlgorithm.h>
#include <vtkCellData.h>
#include <vtkPointData.h>
#include <vtkMath.h>
#include <vtkIdList.h>
#include <vtkCell.h>
#include <vtkAppendPolyData.h>
#include <vtkKdTreePointLocator.h>
#include <vtkCleanPolyData.h>
#include <vtkPolyDataConnectivityFilter.h>

#include "vtkPolyDataBooleanFilter.h"
#include "vtkPolyDataContactFilter.h"

#include "Utilities.h"

#include "Merger.h"
#include "Decomposer.h"

#ifdef DEBUG
#include <chrono>
#include <numeric>
#include <iterator>
#endif

// #include <csignal>

vtkStandardNewMacro(vtkPolyDataBooleanFilter);

vtkPolyDataBooleanFilter::vtkPolyDataBooleanFilter () {

    SetNumberOfInputPorts(2);
    SetNumberOfOutputPorts(2);

    timePdA = 0;
    timePdB = 0;

    contLines = vtkPolyData::New();

    modPdA = vtkPolyData::New();
    modPdB = vtkPolyData::New();

    cellDataA = vtkCellData::New();
    cellDataB = vtkCellData::New();

    cellIdsA = vtkIntArray::New();
    cellIdsB = vtkIntArray::New();

    OperMode = OPER_UNION;

    MergeRegs = false;
    DecPolys = true;

}

vtkPolyDataBooleanFilter::~vtkPolyDataBooleanFilter () {

    cellIdsA->Delete();
    cellIdsB->Delete();

    cellDataB->Delete();
    cellDataA->Delete();

    modPdB->Delete();
    modPdA->Delete();

    contLines->Delete();

}

int vtkPolyDataBooleanFilter::ProcessRequest(vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector) {

    if (request->Has(vtkDemandDrivenPipeline::REQUEST_DATA())) {

        vtkInformation *inInfoA = inputVector[0]->GetInformationObject(0);
        vtkInformation *inInfoB = inputVector[1]->GetInformationObject(0);

        vtkPolyData *pdA = vtkPolyData::SafeDownCast(inInfoA->Get(vtkDataObject::DATA_OBJECT()));
        vtkPolyData *pdB = vtkPolyData::SafeDownCast(inInfoB->Get(vtkDataObject::DATA_OBJECT()));

        vtkInformation *outInfoA = outputVector->GetInformationObject(0);
        vtkInformation *outInfoB = outputVector->GetInformationObject(1);

        resultA = vtkPolyData::SafeDownCast(outInfoA->Get(vtkDataObject::DATA_OBJECT()));
        resultB = vtkPolyData::SafeDownCast(outInfoB->Get(vtkDataObject::DATA_OBJECT()));

#ifdef DEBUG
        using clock = std::chrono::steady_clock;
        std::vector<clock::duration> times;
        clock::time_point start;
#endif

        if (pdA->GetMTime() > timePdA || pdB->GetMTime() > timePdB) {

            // eventuell vorhandene regionen vereinen

            vtkCleanPolyData *cleanA = vtkCleanPolyData::New();
            cleanA->SetOutputPointsPrecision(DOUBLE_PRECISION);
            cleanA->SetTolerance(1e-6);
            cleanA->SetInputData(pdA);
            cleanA->Update();

            vtkCleanPolyData *cleanB = vtkCleanPolyData::New();
            cleanB->SetOutputPointsPrecision(DOUBLE_PRECISION);
            cleanB->SetTolerance(1e-6);
            cleanB->SetInputData(pdB);
            cleanB->Update();

#ifdef DEBUG
            std::cout << "Exporting modPdA.vtk" << std::endl;
            WriteVTK("modPdA.vtk", cleanA->GetOutput());

            std::cout << "Exporting modPdB.vtk" << std::endl;
            WriteVTK("modPdB.vtk", cleanB->GetOutput());
#endif

            // CellData sichern

            cellDataA->DeepCopy(cleanA->GetOutput()->GetCellData());
            cellDataB->DeepCopy(cleanB->GetOutput()->GetCellData());

            // ermittelt kontaktstellen

#ifdef DEBUG

            start = clock::now();
#endif

            vtkPolyDataContactFilter *cl = vtkPolyDataContactFilter::New();
            cl->SetInputConnection(0, cleanA->GetOutputPort());
            cl->SetInputConnection(1, cleanB->GetOutputPort());
            cl->Update();

#ifdef DEBUG
            times.push_back(clock::now()-start);
#endif

            contLines->DeepCopy(cl->GetOutput());

#ifdef DEBUG
            std::cout << "Exporting contLines.vtk" << std::endl;
            WriteVTK("contLines.vtk", contLines);

            std::cout << "Exporting modPdA_1.vtk" << std::endl;
            WriteVTK("modPdA_1.vtk", cl->GetOutput(1));

            std::cout << "Exporting modPdB_1.vtk" << std::endl;
            WriteVTK("modPdB_1.vtk", cl->GetOutput(2));
#endif

            modPdA->DeepCopy(cl->GetOutput(1));
            modPdB->DeepCopy(cl->GetOutput(2));

            if (contLines->GetNumberOfCells() == 0) {
                vtkErrorMacro("Inputs have no contact.");

                return 1;
            }

            // in den CellDatas steht drin, welche polygone einander schneiden

            vtkIntArray *contsA = vtkIntArray::SafeDownCast(contLines->GetCellData()->GetScalars("cA"));
            vtkIntArray *contsB = vtkIntArray::SafeDownCast(contLines->GetCellData()->GetScalars("cB"));

            vtkIntArray *sourcesA = vtkIntArray::SafeDownCast(contLines->GetCellData()->GetScalars("sourcesA"));
            vtkIntArray *sourcesB = vtkIntArray::SafeDownCast(contLines->GetCellData()->GetScalars("sourcesB"));

            // gültigkeit des schnitts prüfen

            bool valid = true;

            vtkIdList *cells = vtkIdList::New();

            for (int i = 0; i < contLines->GetNumberOfPoints() && valid; i++) {
                cells->Reset();
                contLines->GetPointCells(i, cells);

                // der schnitt endet abrupt
                if (cells->GetNumberOfIds() == 1) {
                    valid = false;
                    break;
                }

                vtkIdList *line = vtkIdList::New();

                std::map<int, std::set<int> > countsA, countsB;

                for (int j = 0; j < cells->GetNumberOfIds(); j++) {
                    contLines->GetCellPoints(cells->GetId(j), line);

                    int target = line->GetId(0) == i ? line->GetId(1) : line->GetId(0);

                    int indA = contsA->GetValue(cells->GetId(j));
                    int indB = contsB->GetValue(cells->GetId(j));

                    // sind am punkt sind mehr als zwei linien beteiligt?

                    countsA[indA].insert(target);

                    if (countsA[indA].size() > 2) {
                        valid = false;
                        break;
                    }

                    countsB[indB].insert(target);

                    if (countsB[indB].size() > 2) {
                        valid = false;
                        break;
                    }

                }

                line->Delete();
            }

            cells->Delete();

            if (!valid) {
                vtkErrorMacro("Contact is ambiguous or incomplete.");

                return 1;
            }

            // sichert die OrigCellIds

            vtkIntArray *origCellIdsA = vtkIntArray::SafeDownCast(modPdA->GetCellData()->GetScalars("OrigCellIds"));
            vtkIntArray *origCellIdsB = vtkIntArray::SafeDownCast(modPdB->GetCellData()->GetScalars("OrigCellIds"));

            cellIdsA->DeepCopy(origCellIdsA);
            cellIdsB->DeepCopy(origCellIdsB);

            for (int i = 0; i < modPdA->GetNumberOfCells(); i++) {
                origCellIdsA->SetValue(i, i);
            }

            for (int i = 0; i < modPdB->GetNumberOfCells(); i++) {
                origCellIdsB->SetValue(i, i);
            }


#ifdef DEBUG
            start = clock::now();
#endif

            GetPolyStrips(modPdA, contsA, sourcesA, polyStripsA);
            GetPolyStrips(modPdB, contsB, sourcesB, polyStripsB);

#ifdef DEBUG
            times.push_back(clock::now()-start);
#endif

            // löst ein sehr spezielles problem

#ifdef DEBUG
            start = clock::now();
#endif

            _Test(modPdA, polyStripsA);
            _Test(modPdB, polyStripsB);

#ifdef DEBUG
            times.push_back(clock::now()-start);
#endif

            // trennt die polygone an den linien

#ifdef DEBUG
            start = clock::now();
#endif

            CutCells(modPdA, polyStripsA);
            CutCells(modPdB, polyStripsB);

#ifdef DEBUG
            times.push_back(clock::now()-start);
#endif

#ifdef DEBUG
            std::cout << "Exporting modPdA_2.vtk" << std::endl;
            WriteVTK("modPdA_2.vtk", modPdA);

            std::cout << "Exporting modPdB_2.vtk" << std::endl;
            WriteVTK("modPdB_2.vtk", modPdB);
#endif

#ifdef DEBUG
            start = clock::now();
#endif

            RestoreOrigPoints(modPdA, polyStripsA);
            RestoreOrigPoints(modPdB, polyStripsB);

#ifdef DEBUG
            times.push_back(clock::now()-start);
#endif

#ifdef DEBUG
            std::cout << "Exporting modPdA_3.vtk" << std::endl;
            WriteVTK("modPdA_3.vtk", modPdA);

            std::cout << "Exporting modPdB_3.vtk" << std::endl;
            WriteVTK("modPdB_3.vtk", modPdB);
#endif

#ifdef DEBUG
            start = clock::now();
#endif

            ResolveOverlaps(modPdA, contsA, polyStripsA);
            ResolveOverlaps(modPdB, contsB, polyStripsB);

#ifdef DEBUG
            times.push_back(clock::now()-start);
#endif

#ifdef DEBUG
            std::cout << "Exporting modPdA_4.vtk" << std::endl;
            WriteVTK("modPdA_4.vtk", modPdA);

            std::cout << "Exporting modPdB_4.vtk" << std::endl;
            WriteVTK("modPdB_4.vtk", modPdB);
#endif

#ifdef DEBUG
            start = clock::now();
#endif

            AddAdjacentPoints(modPdA, contsA, polyStripsA);
            AddAdjacentPoints(modPdB, contsB, polyStripsB);

#ifdef DEBUG
            times.push_back(clock::now()-start);
#endif

#ifdef DEBUG
            std::cout << "Exporting modPdA_5.vtk" << std::endl;
            WriteVTK("modPdA_5.vtk", modPdA);

            std::cout << "Exporting modPdB_5.vtk" << std::endl;
            WriteVTK("modPdB_5.vtk", modPdB);
#endif

#ifdef DEBUG
            start = clock::now();
#endif

            DisjoinPolys(modPdA, polyStripsA);
            DisjoinPolys(modPdB, polyStripsB);

#ifdef DEBUG
            times.push_back(clock::now()-start);
#endif

#ifdef DEBUG
            std::cout << "Exporting modPdA_6.vtk" << std::endl;
            WriteVTK("modPdA_6.vtk", modPdA);

            std::cout << "Exporting modPdB_6.vtk" << std::endl;
            WriteVTK("modPdB_6.vtk", modPdB);
#endif

#ifdef DEBUG
            start = clock::now();
#endif

            MergePoints(modPdA, polyStripsA);
            MergePoints(modPdB, polyStripsB);

#ifdef DEBUG
            times.push_back(clock::now()-start);
#endif

#ifdef DEBUG
            std::cout << "Exporting modPdA_7.vtk" << std::endl;
            WriteVTK("modPdA_7.vtk", modPdA);

            std::cout << "Exporting modPdB_7.vtk" << std::endl;
            WriteVTK("modPdB_7.vtk", modPdB);
#endif

            involvedA.clear();
            involvedB.clear();

            int numLines = contLines->GetNumberOfCells();
            for (int i = 0; i < numLines; i++) {
                involvedA.insert(contsA->GetValue(i));
                involvedB.insert(contsB->GetValue(i));
            }

            relsA.clear();
            relsB.clear();

            // aufräumen

            cl->Delete();
            cleanB->Delete();
            cleanA->Delete();

            timePdA = pdA->GetMTime();
            timePdB = pdB->GetMTime();

        }

#ifdef DEBUG
        times.push_back(clock::now()-start);
#endif

        DecPolys_(modPdA, involvedA, relsA);
        DecPolys_(modPdB, involvedB, relsB);

#ifdef DEBUG
        times.push_back(clock::now()-start);
#endif

#ifdef DEBUG
        std::cout << "Exporting modPdA_8.vtk" << std::endl;
        WriteVTK("modPdA_8.vtk", modPdA);

        std::cout << "Exporting modPdB_8.vtk" << std::endl;
        WriteVTK("modPdB_8.vtk", modPdB);
#endif

#ifdef DEBUG
        start = clock::now();
#endif

        if (MergeRegs) {
            MergeRegions();
        } else {
            CombineRegions();
        }

#ifdef DEBUG
        times.push_back(clock::now()-start);
#endif


#ifdef DEBUG
        double sum = std::chrono::duration_cast<std::chrono::duration<double>>(std::accumulate(times.begin(), times.end(), clock::duration())).count();

        std::vector<clock::duration>::const_iterator itr;
        for (itr = times.begin(); itr != times.end(); itr++) {
            double time = std::chrono::duration_cast<std::chrono::duration<double>>(*itr).count();

            std::cout << "Time " << (itr-times.begin())
                << ": " << time << "s (" << (time/sum) << "%)"
                << std::endl;
        }
#endif

    }

    return 1;

}

void vtkPolyDataBooleanFilter::GetStripPoints (vtkPolyData *pd, vtkIntArray *sources, PStrips &pStrips, IdsType &lines) {

#ifdef DEBUG
    std::cout << "GetStripPoints()" << std::endl;
#endif

    StripPtsType &pts = pStrips.pts;
    IdsType &poly = pStrips.poly;

    int numPts = poly.size();

    double tol = 1e-5;

    IdsType::const_iterator itr;

    for (itr = lines.begin(); itr != lines.end(); itr++) {

        vtkIdList *linePts = vtkIdList::New();
        contLines->GetCellPoints(*itr, linePts);

        // diese punkte durchlaufen

        for (int i = 0; i < 2; i++) {

            int realInd = linePts->GetId(i);

            if (pts.count(realInd) == 0) {
                // lage analysieren

                pts[realInd].ind = realInd;

                // die koordinaten
                double pt[3];
                contLines->GetPoint(realInd, pt);

                Cpy(pts[realInd].pt, pt, 3);

                int src = sources->GetComponent(*itr, i);

                pts[realInd].src = src;

                double lastD = DBL_MAX;

                // jetzt muss man die kanten durchlaufen
                for (int j = 0; j < numPts; j++) {

                    if (src > -1 && !(j == src || (j+1)%numPts == src)) {
                        continue;
                    }

                    int indA, indB;
                    double a[3], b[3];

                    indA = poly[j];
                    indB = poly[(j+1)%numPts];

                    pd->GetPoint(indA, a);
                    pd->GetPoint(indB, b);

                    double sA[3], sB[3];
                    vtkMath::Subtract(a, pt, sA);
                    vtkMath::Subtract(b, pt, sB);

                    // richtungsvektor und länge der kante

                    double u[3];
                    vtkMath::Subtract(b, a, u);

                    double n = vtkMath::Norm(u);

                    // d und t zur kante

                    double v[3];
                    vtkMath::Subtract(pt, a, v);

                    double t = vtkMath::Dot(v, u)/(n*n);

                    double w[3];
                    vtkMath::Cross(v, u, w);

                    double d = vtkMath::Norm(w)/n;

                    if (d < tol && d < lastD && t > -tol && t < 1+tol) {
                        // liegt im toleranzbereich der kante

                        pts[realInd].edge[0] = indA;
                        pts[realInd].edge[1] = indB;

                        pts[realInd].t = std::min(1., std::max(0., t));

                        if (vtkMath::Norm(sA) < tol) {
                            Cpy(pts[realInd].captPt, a, 3);
                            pts[realInd].capt = CAPT_A;

                            break;

                        } else if (vtkMath::Norm(sB) < tol) {
                            Cpy(pts[realInd].captPt, b, 3);
                            pts[realInd].capt = CAPT_B;

                            break;

                        } else {
                            // u ist nicht normiert
                            vtkMath::MultiplyScalar(u, t);

                            double x[3];
                            vtkMath::Add(a, u, x);

                            // projektion
                            Cpy(pts[realInd].captPt, x, 3);

                            pts[realInd].capt = CAPT_EDGE;

                            lastD = d;
                        }

                    }

                }

                if (src != NO_USE && pts[realInd].edge[0] == NO_USE) {
                    pts[realInd].catched = false;
                }

            }
        }

        linePts->Delete();

    }

    StripPtsType::iterator itr2;

    for (itr2 = pts.begin(); itr2 != pts.end(); ++itr2) {
        StripPt &sp = itr2->second;

        if (sp.capt != CAPT_NOT) {
            if (sp.capt == CAPT_B) {
                sp.t = 0;

                sp.edge[0] = sp.edge[1];

                int i = std::find(poly.begin(), poly.end(), sp.edge[0])-poly.begin();

                sp.edge[1] = poly[(i+1)%numPts];

                sp.capt = CAPT_A;

            }

            // für den schnitt werden die eingerasteten koordinaten verwendet

            Cpy(sp.cutPt, sp.captPt, 3);
        } else {

            Cpy(sp.cutPt, sp.pt, 3);
        }

        sp.history.push_back({sp.edge[0], sp.edge[1]});

    }

#ifdef DEBUG
    std::cout << "pts: " << std::endl;
    for (itr2 = pts.begin(); itr2 != pts.end(); ++itr2) {
        std::cout << itr2->first << ": " << itr2->second << std::endl;
    }
#endif

}

void vtkPolyDataBooleanFilter::GetPolyStrips (vtkPolyData *pd, vtkIntArray *conts, vtkIntArray *sources, PolyStripsType &polyStrips) {
#ifdef DEBUG
    std::cout << "GetPolyStrips()" << std::endl;
#endif

    polyStrips.clear();

    std::map<int, IdsType> polyLines;

    for (int i = 0; i < conts->GetNumberOfTuples(); i++) {
        int poly = conts->GetValue(i);
        polyLines[poly].push_back(i);
    }

    std::vector<std::reference_wrapper<StripPt>> notCatched;

    std::map<int, IdsType>::iterator itr;

    for (itr = polyLines.begin(); itr != polyLines.end(); ++itr) {
        IdsType &lines = itr->second;
        RemoveDuplicates(lines);

        PStrips &pStrips = polyStrips[itr->first];

        vtkIdList *polyPts = vtkIdList::New();
        pd->GetCellPoints(itr->first, polyPts);

        int numPts = polyPts->GetNumberOfIds();

        for (int i = 0; i < numPts; i++) {
            pStrips.poly.push_back(polyPts->GetId(i));
        }

        ComputeNormal(pd->GetPoints(), pStrips.n, polyPts);

        GetStripPoints(pd, sources, pStrips, lines);

        for (auto &sp : pStrips.pts) {
            sp.second.polyId = itr->first;

            if (!sp.second.catched) {
                notCatched.push_back(sp.second);
            }
        }

        polyPts->Delete();
    }

    auto Next = [](const IdsType &ids, int id) -> int {
        IdsType::const_iterator itr;

        itr = std::find(ids.begin(), ids.end(), id);

        if (++itr == ids.end()) {
            itr = ids.begin();
        }

        return *itr;
    };

    for (StripPt &sp : notCatched) {
        for (itr = polyLines.begin(); itr != polyLines.end(); ++itr) {
            const PStrips &pStrips = polyStrips[itr->first];

            try {
                const StripPt &corr = pStrips.pts.at(sp.ind);

                if (&corr != &sp) {
                    if (corr.capt == CAPT_A) {
                        sp.capt = CAPT_A;
                        sp.edge[0] = corr.edge[0];
                        sp.edge[1] = Next(polyStrips[sp.polyId].poly, sp.edge[0]);

                        sp.t = 0;

                        Cpy(sp.captPt, corr.captPt, 3);
                        Cpy(sp.cutPt, sp.captPt, 3);

                        sp.history.push_back({sp.edge[0], sp.edge[1]});

                        sp.catched = true;

                    }
                }
            } catch (...) {}

        }

        assert(sp.catched);

    }

    for (itr = polyLines.begin(); itr != polyLines.end(); ++itr) {
        IdsType &lines = itr->second;

        PStrips &pStrips = polyStrips[itr->first];

        StripType strip;

        // zusammensetzen

        std::deque<int> _lines(lines.begin(), lines.end());

        int i = 0;

        while (_lines.size() > 0) {
            vtkIdList *linePts = vtkIdList::New();
            contLines->GetCellPoints(_lines[i], linePts);

            int indA = linePts->GetId(0);
            int indB = linePts->GetId(1);

            if (strip.empty()) {
                strip.push_back(StripPtR(indA));
                strip.push_back(StripPtR(indB));

                _lines.erase(_lines.begin());
            } else {
                StripPt &start = pStrips.pts[strip.front().ind],
                    &end = pStrips.pts[strip.back().ind];

                if (end.capt == CAPT_NOT && end.ind == indA) {
                    strip.push_back(StripPtR(indB));
                    _lines.erase(_lines.begin()+i);
                    i = 0;

                } else if (end.capt == CAPT_NOT && end.ind == indB) {
                    strip.push_back(StripPtR(indA));
                    _lines.erase(_lines.begin()+i);
                    i = 0;

                } else if (start.capt == CAPT_NOT && start.ind == indA) {
                    strip.push_front(StripPtR(indB));
                    _lines.erase(_lines.begin()+i);
                    i = 0;

                } else if (start.capt == CAPT_NOT && start.ind == indB) {
                    strip.push_front(StripPtR(indA));
                    _lines.erase(_lines.begin()+i);
                    i = 0;

                } else {
                    i++;

                }
            }

            StripPt &_start = pStrips.pts[strip.front().ind],
                &_end = pStrips.pts[strip.back().ind];

            if ((_start.capt != CAPT_NOT && _end.capt != CAPT_NOT)
                || (_start.ind == _end.ind)
                || (i == _lines.size())) {

                // einen neuen strip anlegen
                pStrips.strips.push_back(strip);
                strip.clear();

                i = 0;
            }

            linePts->Delete();

        }

        CompleteStrips(pStrips);

    }

}

void vtkPolyDataBooleanFilter::RemoveDuplicates (IdsType &lines) {

    IdsType unique;

    int i, j;

    // die indexe der enden auf übereinstimmung prüfen

    vtkIdList *linePtsA = vtkIdList::New();
    vtkIdList *linePtsB = vtkIdList::New();

    int numLines = lines.size();

    for (i = 0; i < numLines-1; i++) {
        j = i+1;

        contLines->GetCellPoints(lines[i], linePtsA);

        while (j < lines.size()) {
            contLines->GetCellPoints(lines[j], linePtsB);

            if ((linePtsA->GetId(0) == linePtsB->GetId(0) && linePtsA->GetId(1) == linePtsB->GetId(1)) ||
                (linePtsA->GetId(0) == linePtsB->GetId(1) && linePtsA->GetId(1) == linePtsB->GetId(0))) {
                // stimmen überein
                break;
            }

            j++;
        }

        if (j == numLines) {
            // keine vorzeitige unterbrechung der while-schleife
            unique.push_back(lines[i]);
        }
    }

    unique.push_back(lines.back());

    linePtsA->Delete();
    linePtsB->Delete();

    lines.swap(unique);

}

void vtkPolyDataBooleanFilter::CompleteStrips (PStrips &pStrips) {
    StripsType::iterator itr;

    for (itr = pStrips.strips.begin(); itr != pStrips.strips.end(); ++itr) {
        StripPt &start = pStrips.pts[itr->front().ind],
            &end = pStrips.pts[itr->back().ind];

        if (start.ind != end.ind) {
            if (start.capt == CAPT_NOT) {
                StripType s(itr->rbegin(), itr->rend()-1);
                itr->insert(itr->begin(), s.begin(), s.end());

            } else if (end.capt == CAPT_NOT) {
                StripType s(itr->rbegin()+1, itr->rend());
                itr->insert(itr->end(), s.begin(), s.end());

            }
        }
    }
}

bool vtkPolyDataBooleanFilter::HasArea (StripType &strip) {
    bool area = true;

    int n = strip.size();
    if (n%2 == 1) {
        for (int i = 0; i < (n-1)/2; i++) {
            area = strip[i].ind != strip[n-i-1].ind;
        }
    }

    return area;
}

void vtkPolyDataBooleanFilter::CutCells (vtkPolyData *pd, PolyStripsType &polyStrips) {
#ifdef DEBUG
    std::cout << "CutCells()" << std::endl;
#endif

    vtkPoints *pdPts = pd->GetPoints();

    vtkIntArray *origCellIds = vtkIntArray::SafeDownCast(pd->GetCellData()->GetScalars("OrigCellIds"));

    PolyStripsType::iterator itr;

    for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) {

        int polyInd = itr->first;
        PStrips &pStrips = itr->second;

#ifdef DEBUG
        if (polyInd != 204) {
            //continue;
        }
#endif

        StripsType &strips = pStrips.strips;
        StripPtsType &pts = pStrips.pts;

        IdsType &poly = pStrips.poly;

        int origId = origCellIds->GetValue(polyInd);

#ifdef DEBUG
        IdsType::iterator itr_;
        std::cout << "polyInd=" << polyInd << ", poly=[";
        for (itr_ = poly.begin(); itr_ != poly.end(); ++itr_) {
            std::cout << *itr_ << ", ";
        }
        std::cout << "]" << std::endl;
#endif

        int numPts = poly.size();

        std::map<int, RefsType> edges;

        StripsType::iterator itr2;
        StripType::iterator itr3;

        // holes sammeln

        HolesType holes;

        for (itr2 = strips.begin(); itr2 != strips.end(); ++itr2) {
            StripType &s = *itr2;
            if (pts[s.front().ind].capt == CAPT_NOT && pts[s.back().ind].capt == CAPT_NOT) {
                IdsType hole;

                for (auto& sp : s) {
                    hole.push_back(pdPts->InsertNextPoint(pts[sp.ind].pt));
                }

                // anfang und ende sind ja gleich
                hole.pop_back();

                holes.push_back(std::move(hole));

            }
        }

        // holes löschen
        strips.erase(std::remove_if(strips.begin(), strips.end(), [&](const StripType &s) {
            return pts[s.front().ind].capt == CAPT_NOT && pts[s.back().ind].capt == CAPT_NOT; }), strips.end());

        for (itr2 = strips.begin(); itr2 != strips.end(); ++itr2) {
            StripType &strip = *itr2;

#ifdef DEBUG
            std::cout << "strip [";
            for (auto &s : strip) {
                std::cout << s.ind << ", ";
            }
            std::cout << "]" << std::endl;
#endif

            // init
            if (pts[strip.front().ind].edge[0] == pts[strip.back().ind].edge[0]
                && strip.front().ind != strip.back().ind
                && pts[strip.front().ind].t > pts[strip.back().ind].t) {

                std::reverse(strip.begin(), strip.end());
            }

            StripPt &start = pts[strip.front().ind],
                &end = pts[strip.back().ind];

            strip.front().side = SIDE_START;
            strip.back().side = SIDE_END;

            strip.front().ref = start.edge[0];
            strip.back().ref = end.edge[0];

            int ind = itr2-strips.begin();

            strip.front().strip = strip.back().strip = ind;

            for (itr3 = strip.begin(); itr3 != strip.end(); ++itr3) {
                StripPt &sp = pts[itr3->ind];

                itr3->desc[0] = pdPts->InsertNextPoint(sp.cutPt);
                itr3->desc[1] = pdPts->InsertNextPoint(sp.cutPt);

#ifdef DEBUG
                std::cout << sp << " => " << *itr3 << std::endl;
#endif

            }

            // ordnet zu

            edges[start.edge[0]].push_back(std::ref(strip.front()));
            edges[end.edge[0]].push_back(std::ref(strip.back()));
        }

        // sortiert die punkte auf den kanten

        std::map<int, RefsType>::iterator itr4;

        IdsType::iterator itr5;
        StripType::reverse_iterator itr6;

        RefsType::iterator itr7;

        std::vector<IdsType>::iterator itr8;


        for (itr4 = edges.begin(); itr4 != edges.end(); ++itr4) {
            RefsType &edge = itr4->second;

#ifdef DEBUG
            std::cout << "edge (" << itr4->first << ", _)" << std::endl;
#endif

            std::sort(edge.begin(), edge.end(), [&](const StripPtR &a, const StripPtR &b) {
                StripPt &a_ = pts[a.ind],
                    &b_ = pts[b.ind];

#ifdef DEBUG
                std::cout << "a_: " << a_ << " -> strip " << a.strip << std::endl;
                std::cout << "b_: " << b_ << " -> strip " << b.strip << std::endl;
#endif

                if (a_.ind == b_.ind) {
                    // strips beginnen im gleichen punkt

                    if (a.strip != b.strip) {
                        // gehören nicht dem gleichen strip an

                        StripType &stripA = strips[a.strip],
                            &stripB = strips[b.strip];

                        // andere enden ermitteln

                        StripPtR &eA = (&a == &(stripA.front())) ? stripA.back() : stripA.front(),
                            &eB = (&b == &(stripB.front())) ? stripB.back() : stripB.front();

                        StripPt &eA_ = pts[eA.ind],
                            &eB_ = pts[eB.ind];

#ifdef DEBUG
                        std::cout << "eA_: " << eA_ << std::endl;
                        std::cout << "eB_: " << eA_ << std::endl;
#endif

                        if (eA_.ind != eB_.ind) {
                            int i = std::find(poly.begin(), poly.end(), itr4->first)-poly.begin();

                            int iA = std::find(poly.begin(), poly.end(), eA_.edge[0])-poly.begin(),
                                iB = std::find(poly.begin(), poly.end(), eB_.edge[0])-poly.begin();

                            double dA = Mod(iA-i, numPts)+eA_.t,
                                dB = Mod(iB-i, numPts)+eB_.t;

                            if (i == iA && a_.t > eA_.t) {
                               dA += numPts;
                            }

                            if (i == iB && b_.t > eB_.t) {
                               dB += numPts;
                            }

#ifdef DEBUG
                            std::cout << "dA=" << dA << ", dB=" << dB << std::endl;
#endif

                            return dB < dA;
                        } else {
                            RefsType poly_;

                            if (a.side == SIDE_START) {
                                poly_.insert(poly_.end(), stripA.begin(), stripA.end());
                            } else {
                                poly_.insert(poly_.end(), stripA.rbegin(), stripA.rend());
                            }

                            if (b.side == SIDE_START) {
                                poly_.insert(poly_.end(), stripB.rbegin()+1, stripB.rend()-1);
                            } else {
                                poly_.insert(poly_.end(), stripB.begin()+1, stripB.end()-1);
                            }

                            int num = poly_.size();

                            vtkPoints *pts_ = vtkPoints::New();
                            pts_->SetNumberOfPoints(num);

                            for (itr7 = poly_.begin(); itr7 != poly_.end(); ++itr7) {
                                StripPtR &sp = *itr7;
                                int i = itr7-poly_.begin();

                                pts_->SetPoint(i, pts[sp.ind].cutPt);
                            }

                            double n[3];
                            ComputeNormal(pts_, n);

                            pts_->Delete();

                            double ang = vtkMath::Dot(pStrips.n, n);

#ifdef DEBUG
                            std::cout << "ang=" << ang*180/PI << std::endl;
#endif

                            return ang < .999999;

                        }
                    } else {
                        // gleicher strip

                        StripType &strip = strips[a.strip];

                        if (HasArea(strip)) {
                            RefsType poly_(strip.begin(), strip.end()-1);

                            int num = poly_.size();

                            vtkPoints *pts_ = vtkPoints::New();
                            pts_->SetNumberOfPoints(num);

                            for (itr7 = poly_.begin(); itr7 != poly_.end(); ++itr7) {
                                StripPtR &sp = *itr7;
                                int i = itr7-poly_.begin();

                                pts_->SetPoint(i, pts[sp.ind].cutPt);
                            }

                            double n[3];
                            ComputeNormal(pts_, n);

                            pts_->Delete();

                            double ang = vtkMath::Dot(pStrips.n, n);

                            if (ang > .999999) {
                                std::reverse(strip.begin(), strip.end());
                                return true;
                            } else {
                                return false;
                            }

                        } else {
                            // reihenfolge von a und b bereits richtig
                            return false;
                        }
                    }

                } else {
                    return a_.t < b_.t;
                }
            });

#ifdef DEBUG
            std::cout << "after sort" << std::endl;
            for (itr7 = edge.begin(); itr7 != edge.end(); ++itr7) {
                std::cout << itr7->get().strip << ": " << pts[itr7->get().ind] << std::endl;
            }
#endif

        }

        // baut die strips ein

        std::deque<IdsType> polys;
        polys.push_back(pStrips.poly);

        int i = 0;

        std::deque<IdsType>::iterator itr9;

        for (itr2 = strips.begin(); itr2 != strips.end(); ++itr2) {
            StripType &strip = *itr2;

            StripPtR &start = strip.front(),
                &end = strip.back();

#ifdef DEBUG
            std::cout << "strip " << start.strip
                << " refs=[" << start.ref << ", " << end.ref << "]"
                << std::endl;
#endif

            int cycle = 0;

            while (true) {

                if (cycle == polys.size()) {
                    break;
                }

                IdsType p = polys.front();
                polys.pop_front();

                std::vector<IdsType> divided(2);

                if (std::find(p.begin(), p.end(), start.ref) != p.end()) {
                    if (start.ref == end.ref) {
                        for (itr5 = p.begin(); itr5 != p.end(); ++itr5) {
                            divided[0].push_back(*itr5);

#ifdef DEBUG
                            std::cout << "add d[0] " << *itr5 << std::endl;
#endif

                            if (*itr5 == start.ref) {
                                for (itr3 = strip.begin(); itr3 != strip.end(); ++itr3) {
                                    divided[0].push_back(itr3->desc[0]);

#ifdef DEBUG
                                    std::cout << "add d[0] " << itr3->desc[0] << std::endl;
#endif

                                }
                            }
                        }

                        // strip selbst ist ein polygon

                        for (itr6 = strip.rbegin(); itr6 != strip.rend(); ++itr6) {
                            divided[1].push_back(itr6->desc[1]);

#ifdef DEBUG
                            std::cout << "add d[1] " << itr6->desc[1] << std::endl;
#endif

                        }

                    } else {
                        int d = 0;
                        for (itr5 = p.begin(); itr5 != p.end(); ++itr5) {
                            divided[d].push_back(*itr5);

#ifdef DEBUG
                            std::cout << "add d[" << d << "] " << *itr5 << std::endl;
#endif

                            if (*itr5 == start.ref) {
                                for (itr3 = strip.begin(); itr3 != strip.end(); ++itr3) {
                                    divided[d].push_back(itr3->desc[0]);

#ifdef DEBUG
                                    std::cout << "add d[" << d << "] " << itr3->desc[0] << std::endl;
#endif

                                }
                                d++;
                            } else if (*itr5 == end.ref) {
                                for (itr6 = strip.rbegin(); itr6 != strip.rend(); ++itr6) {
                                    divided[d].push_back(itr6->desc[1]);

#ifdef DEBUG
                                    std::cout << "add d[" << d << "] " << itr6->desc[1] << std::endl;
#endif

                                }
                                d++;
                            }

                            d %= 2;
                        }
                    }
                }

                if (divided[1].size() > 0) {

                    // refs aktualisieren

                    int stripInd = itr2-strips.begin();

                    for (itr4 = edges.begin(); itr4 != edges.end(); ++itr4) {
                        RefsType &edge = itr4->second;

                        for (itr7 = edge.begin()+1; itr7 != edge.end(); ++itr7) {
                            StripPtR &sp = *itr7;

                            RefsType::reverse_iterator itr7_(itr7);

                            for (; itr7_ != edge.rend(); ++itr7_) {
                                // pt davor
                                StripPtR &pre = *itr7_;

                                if (pre.strip != sp.strip) {
                                    if (pre.strip <= stripInd) {
                                        // der davor wurde bereits verwendet

                                        if (pre.side == SIDE_END) {
                                            sp.ref = pre.desc[0];
                                        } else {
                                            sp.ref = pre.desc[1];
                                        }

                                        break;
                                    }
                                } else {
                                    // gehört dem gleichen strip an
                                    sp.ref = pre.ref;

                                    break;
                                }
                            }

                        }

                        // erstellt die history

                        auto _s = std::find_if(edge.begin(), edge.end(), [&start](const StripPtR &r) {
                            return &start == &r;
                        });

                        if (_s != edge.end()) {
                            for (itr7 = edge.begin(); itr7 != edge.end(); ++itr7) {
                                StripPtR &sp = *itr7;
                                StripPt &_sp = pts[sp.ind];
                                auto &history = _sp.history;

                                if (&sp != &start) {
                                    if (_sp.t < pts[start.ind].t) {
                                        history.push_back({history.back().f, start.desc[0]});
                                    } else {
                                        history.push_back({start.desc[1], history.back().g});
                                    }

                                }

                            }
                        }

                        auto _e = std::find_if(edge.begin(), edge.end(), [&end](const StripPtR &r) {
                            return &end == &r;
                        });

                        if (_e != edge.end()) {
                            for (itr7 = edge.begin(); itr7 != edge.end(); ++itr7) {
                                StripPtR &sp = *itr7;
                                StripPt &_sp = pts[sp.ind];
                                auto &history = _sp.history;

                                if (&sp != &end) {
                                    if (_sp.t < pts[end.ind].t) {
                                        history.push_back({history.back().f, end.desc[1]});
                                    } else {
                                        history.push_back({end.desc[0], history.back().g});
                                    }

                                }

                            }
                        }

                        // sonderfall
                        if (edge.size() > 1) {
                            StripPtR &a = edge.front(),
                                &b = *(edge.begin()+1);

                            if (a.ind == b.ind
                                && b.strip == stripInd
                                && pts[a.ind].capt == CAPT_A) {

                                if (b.side == SIDE_START) {
                                    a.ref = b.desc[0];
                                } else {
                                    a.ref = b.desc[1];
                                }

                            }
                        }
                    }

                    // doppelte punkte entfernen

                    for (itr8 = divided.begin(); itr8 != divided.end(); ++itr8) {
#ifdef DEBUG
                        std::cout << "d[" << itr8-divided.begin() << "]" << std::endl;
#endif

                        IdsType &div = *itr8;

                        IdsType _div;

                        int num = div.size();

                        for (int i = 0; i < num; i++) {
                            double ptA[3], ptB[3];
                            pd->GetPoint(div[i], ptA);
                            pd->GetPoint(div[(i+1)%num], ptB);

                            double d = GetD(ptA, ptB);

#ifdef DEBUG
                            std::cout << "edge (" << div[i] << ", "
                                      << div[(i+1)%num] << "), d=" << d
                                      << std::endl;
#endif

                            if (d > 1e-6) {
                                _div.push_back(div[i]);
                            } else {
#ifdef DEBUG
                                std::cout << "rm " << div[i] << std::endl;
#endif
                            }
                        }

                        div = std::move(_div);
                    }

                    // prüfen ob erstellte polygone gültig

                    if (divided[0].size() > 2) {
                        polys.push_back(divided[0]);
                    }

                    if (HasArea(strip) && divided[1].size() > 2) {
                        polys.push_back(divided[1]);
                    }

                    cycle = 0;

                    break;

                } else {
                    polys.push_back(p);

                    cycle++;
                }

            }

            i++;

            if (i == 4) {
                //break;
            }

        }

        // erzeugte polys hinzufügen

        IdsType descIds;

        for (itr9 = polys.begin(); itr9 != polys.end(); ++itr9) {
            IdsType &p = *itr9;

            int num = p.size();

            vtkIdList *cell = vtkIdList::New();
            cell->SetNumberOfIds(num);

            for (int i = 0; i < num; i++) {
                cell->SetId(i, p[i]);
            }

            descIds.push_back(pd->InsertNextCell(VTK_POLYGON, cell));

            origCellIds->InsertNextValue(origId);

            cell->Delete();
        }

        pd->DeleteCell(polyInd);

        // holes verarbeiten

        if (!holes.empty()) {
            _Wrapper w(pd, descIds, origId);

            for (auto& hole : holes) {
                w.Add(hole);
            }

            w.MergeAll();
        }

    }

    pd->RemoveDeletedCells();
    pd->BuildCells();

}

void vtkPolyDataBooleanFilter::_Test (vtkPolyData *pd, PolyStripsType &polyStrips) {

#ifdef DEBUG
    std::cout << "_Test()" << std::endl;
#endif

    struct Cmp {
        bool operator() (const StripPt &l, const StripPt &r) const {
            const int x1 = static_cast<int>(l.cutPt[0]*1e5),
                y1 = static_cast<int>(l.cutPt[1]*1e5),
                z1 = static_cast<int>(l.cutPt[2]*1e5),
                x2 = static_cast<int>(r.cutPt[0]*1e5),
                y2 = static_cast<int>(r.cutPt[1]*1e5),
                z2 = static_cast<int>(r.cutPt[2]*1e5);

            return std::tie(x1, y1, z1) < std::tie(x2, y2, z2);
        }
    };

    struct Cmp2 {
        bool operator() (const StripPt &l, const StripPt &r) const {
            return l.ind < r.ind;
        }
    };

    std::map<StripPt, std::set<std::reference_wrapper<StripPt>, Cmp2>, Cmp> test;

    PolyStripsType::iterator itr;
    StripPtsType::iterator itr2;

    for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) {
        PStrips &pStrips = itr->second;

        for (itr2 = pStrips.pts.begin(); itr2 != pStrips.pts.end(); ++itr2) {
            StripPt &sp = itr2->second;

            test[sp].insert(sp);
        }
    }

    vtkIdList *cells = vtkIdList::New();

    for (auto &s : test) {
        auto &pts = s.second;

        if (pts.size() > 1) {
            StripPt &_a = *(pts.begin()),
                &_b = *(std::next(pts.begin()));

            std::reference_wrapper<StripPt> __a(_a), __b(_b);

            if (_a.capt == CAPT_EDGE && _b.capt == CAPT_A) {
                // so ist a niemals CAPT_EDGE, wenn b CAPT_A ist
                __a = _b;
                __b = _a;
            }

            StripPt &a = __a,
                &b = __b;

            int indA = a.ind,
                indB = b.ind;

#ifdef DEBUG
            std::cout << "collapsing " << b.ind << " -> " << a.ind << std::endl;
#endif

            contLines->GetPointCells(b.ind, cells);

            for (int i = 0; i < cells->GetNumberOfIds(); i++) {
                contLines->ReplaceCellPoint(cells->GetId(i), indB, indA);

                // aktualisiert die links
                contLines->RemoveReferenceToCell(indB, cells->GetId(i));
                contLines->ResizeCellList(indA, 1);
                contLines->AddReferenceToCell(indA, cells->GetId(i));
            }

            std::set<Pair> pairs;

            if (a.capt == CAPT_A && b.capt == CAPT_A) {
                std::map<int, IdsType> shared;

                vtkIdList *lines = vtkIdList::New(),
                    *cell = vtkIdList::New();
                contLines->GetPointCells(indA, lines);

                for (int i = 0; i < lines->GetNumberOfIds(); i++) {
                    contLines->GetCellPoints(lines->GetId(i), cell);

                    int pA = cell->GetId(0),
                        pB = cell->GetId(1);

                    shared[pA == indA ? pB : pA].push_back(lines->GetId(i));
                }

                cell->Delete();
                lines->Delete();

                for (auto &s : shared) {
                    if (s.second.size() > 1) {
                        assert(s.second.size() == 2);

                        for (int l : s.second) {
                            contLines->DeleteCell(l);
                            contLines->RemoveCellReference(l);

                            pairs.insert({indA, s.first});
                            pairs.insert({s.first, indA});
                        }
                    }
                }

            }

            if (a.capt == CAPT_A && b.capt == CAPT_EDGE) {
                vtkIdList *lines = vtkIdList::New(),
                    *cell = vtkIdList::New();
                contLines->GetPointCells(indA, lines);

                for (int i = 0; i < lines->GetNumberOfIds(); i++) {
                    if (contLines->GetCellType(lines->GetId(i)) != VTK_EMPTY_CELL) {
                        contLines->GetCellPoints(lines->GetId(i), cell);

                        int pA = cell->GetId(0),
                            pB = cell->GetId(1);

                        if (pA == indA && pB == indA) {
                            contLines->DeleteCell(lines->GetId(i));
                            contLines->RemoveCellReference(lines->GetId(i));
                        }
                    }

                }

                cell->Delete();
                lines->Delete();
            }

            auto Fct = [&](PolyStripsType &polyStrips_) -> void {
                for (auto &ps : polyStrips_) {
                    StripsType &strips = ps.second.strips;
                    StripPtsType &pts = ps.second.pts;

                    if (pts.count(indB) == 1) {
                        if (pts.count(indA) == 0) {
                            if (a.capt == pts[indB].capt) {
                                pts[indA] = pts[indB];
                                pts[indA].ind = indA;
                                Cpy(pts[indA].pt, a.pt, 3);
                            } else /*if (a.capt == CAPT_A && pts[indB].capt == CAPT_EDGE)*/ {
                                pts[indA] = a;
{'fraction_non_alphanumeric': 0.1403004838298956, 'fraction_numerical': 0.008423733129615482, 'mean_word_length': 1.2907809692698975, 'pattern_counts': {'":': 6, '<': 824, '<?xml version=': 0, '>': 604, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 123, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3194444', 'n_tokens_mistral': 31075, 'n_tokens_neox': 29168, 'n_words': 5473}