From ccd9515764b83d719668120ccf8496a8fc6e804b Mon Sep 17 00:00:00 2001 From: Eduardo DAVILA Date: Thu, 25 Aug 2022 10:46:59 +0200 Subject: [PATCH] #3493 MeshManager --- ...bcreaVtkBooleanOperationPolyDataFilter.cxx | 29 +- bbtk_creaVtk_PKG/src/bbcreaVtkMeshManager.cxx | 75 + bbtk_creaVtk_PKG/src/bbcreaVtkMeshManager.h | 56 + .../src/bbcreaVtkMeshManager_tool.cxx | 90 + .../src/bbcreaVtkMeshManager_tool.h | 54 + lib/creaVtk/MeshManagerModel.cpp | 59 + lib/creaVtk/MeshManagerModel.h | 70 + lib/creaVtk/Utilities.cxx | 239 ++ lib/creaVtk/Utilities.h | 141 + lib/creaVtk/vtkPolyDataBooleanFilter.cxx | 3376 +++++++++++++++++ lib/creaVtk/vtkPolyDataBooleanFilter.h | 298 ++ lib/creaVtk/vtkPolyDataContactFilter.cxx | 912 +++++ lib/creaVtk/vtkPolyDataContactFilter.h | 151 + 13 files changed, 5546 insertions(+), 4 deletions(-) create mode 100644 bbtk_creaVtk_PKG/src/bbcreaVtkMeshManager.cxx create mode 100644 bbtk_creaVtk_PKG/src/bbcreaVtkMeshManager.h create mode 100644 bbtk_creaVtk_PKG/src/bbcreaVtkMeshManager_tool.cxx create mode 100644 bbtk_creaVtk_PKG/src/bbcreaVtkMeshManager_tool.h create mode 100644 lib/creaVtk/MeshManagerModel.cpp create mode 100644 lib/creaVtk/MeshManagerModel.h create mode 100644 lib/creaVtk/Utilities.cxx create mode 100644 lib/creaVtk/Utilities.h create mode 100644 lib/creaVtk/vtkPolyDataBooleanFilter.cxx create mode 100644 lib/creaVtk/vtkPolyDataBooleanFilter.h create mode 100644 lib/creaVtk/vtkPolyDataContactFilter.cxx create mode 100644 lib/creaVtk/vtkPolyDataContactFilter.h diff --git a/bbtk_creaVtk_PKG/src/bbcreaVtkBooleanOperationPolyDataFilter.cxx b/bbtk_creaVtk_PKG/src/bbcreaVtkBooleanOperationPolyDataFilter.cxx index 3e8bca3..6cb69d5 100644 --- a/bbtk_creaVtk_PKG/src/bbcreaVtkBooleanOperationPolyDataFilter.cxx +++ b/bbtk_creaVtk_PKG/src/bbcreaVtkBooleanOperationPolyDataFilter.cxx @@ -4,10 +4,13 @@ #include "bbcreaVtkBooleanOperationPolyDataFilter.h" #include "bbcreaVtkPackage.h" -#include "vtkBooleanOperationPolyDataFilter.h" +// #include "vtkBooleanOperationPolyDataFilter.h" +#include "vtkPolyDataBooleanFilter.h" + #include "vtkCleanPolyData.h" #include "vtkTriangleFilter.h" + namespace bbcreaVtk { @@ -57,23 +60,41 @@ printf("EED Warnning BooleanOperationPolyDataFilter::Process Put this code at t triangle1->Update(); triangle2->Update(); + vtkPolyDataBooleanFilter *booleanOperation = vtkPolyDataBooleanFilter::New(); + booleanOperation->SetInputData(0, triangle1->GetOutput() ); + booleanOperation->SetInputData(1, triangle2->GetOutput() ); + if (bbGetInputOperation()==0 ) + { + booleanOperation->SetOperModeToUnion(); + } + if (bbGetInputOperation()==1 ) + { + booleanOperation->SetOperModeToIntersection(); + } + if (bbGetInputOperation()==2 ) + { + booleanOperation->SetOperModeToDifference(); + } + +/* vtkBooleanOperationPolyDataFilter *booleanOperation = vtkBooleanOperationPolyDataFilter::New(); booleanOperation->SetInputData(0, triangle1->GetOutput() ); booleanOperation->SetInputData(1, triangle2->GetOutput() ); - if (bbGetInputOperation()==0 ) { booleanOperation->SetOperationToUnion(); } if (bbGetInputOperation()==1 ) { - booleanOperation->SetOperationToIntersection(); + booleanOperation->SetOperationToIntersection(); } if (bbGetInputOperation()==2 ) { - booleanOperation->SetOperationToDifference(); + booleanOperation->SetOperationToDifference(); booleanOperation->SetReorientDifferenceCells( bbGetInputReorientDifferenceCells() ); } + */ + booleanOperation->Update(); bbSetOutputOut( booleanOperation->GetOutput() ); } else { diff --git a/bbtk_creaVtk_PKG/src/bbcreaVtkMeshManager.cxx b/bbtk_creaVtk_PKG/src/bbcreaVtkMeshManager.cxx new file mode 100644 index 0000000..cb4cca8 --- /dev/null +++ b/bbtk_creaVtk_PKG/src/bbcreaVtkMeshManager.cxx @@ -0,0 +1,75 @@ +//===== +// Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost) +//===== +#include "bbcreaVtkMeshManager.h" +#include "bbcreaVtkPackage.h" +namespace bbcreaVtk +{ + +BBTK_ADD_BLACK_BOX_TO_PACKAGE(creaVtk,MeshManager) +BBTK_BLACK_BOX_IMPLEMENTATION(MeshManager,bbtk::AtomicBlackBox); +//===== +// Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost) +//===== +void MeshManager::Process() +{ +// THE MAIN PROCESSING METHOD BODY +// Here we simply set the input 'In' value to the output 'Out' +// And print out the output value +// INPUT/OUTPUT ACCESSORS ARE OF THE FORM : +// void bbSet{Input|Output}NAME(const TYPE&) +// const TYPE& bbGet{Input|Output}NAME() const +// Where : +// * NAME is the name of the input/output +// (the one provided in the attribute 'name' of the tag 'input') +// * TYPE is the C++ type of the input/output +// (the one provided in the attribute 'type' of the tag 'input') +// bbSetOutputOut( bbGetInputIn() ); +// std::cout << "Output value = " <SetMeshBase( bbGetInputMesh() ); + bbSetOutputMeshBase( meshManagerModel->GetMeshBase() ); + bbSetOutputMeshTemp( meshManagerModel->GetMeshTemp() ); + bbSetOutputMeshManagerModel( meshManagerModel ); + } +} + +//===== +// Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost) +//===== +void MeshManager::bbUserSetDefaultValues() +{ +// SET HERE THE DEFAULT INPUT/OUTPUT VALUES OF YOUR BOX +// Here we initialize the input 'In' to 0 +// bbSetInputIn(0); + meshManagerModel = NULL; + bbSetInputMesh(NULL); +} + +//===== +// Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost) +//===== +void MeshManager::bbUserInitializeProcessing() +{ +// THE INITIALIZATION METHOD BODY : +// Here does nothing +// but this is where you should allocate the internal/output pointers +// if any +} + +//===== +// Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost) +//===== +void MeshManager::bbUserFinalizeProcessing() +{ +// THE FINALIZATION METHOD BODY : +// Here does nothing +// but this is where you should desallocate the internal/output pointers +// if any +} + +} // EO namespace bbcreaVtk + + diff --git a/bbtk_creaVtk_PKG/src/bbcreaVtkMeshManager.h b/bbtk_creaVtk_PKG/src/bbcreaVtkMeshManager.h new file mode 100644 index 0000000..315ee5d --- /dev/null +++ b/bbtk_creaVtk_PKG/src/bbcreaVtkMeshManager.h @@ -0,0 +1,56 @@ +//===== +// Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost) +//===== +#ifndef __bbcreaVtkMeshManager_h_INCLUDED__ +#define __bbcreaVtkMeshManager_h_INCLUDED__ + +#include "bbcreaVtk_EXPORT.h" +#include "bbtkAtomicBlackBox.h" +#include "iostream" +#include +#include + +namespace bbcreaVtk +{ + +class bbcreaVtk_EXPORT MeshManager + : + public bbtk::AtomicBlackBox +{ + BBTK_BLACK_BOX_INTERFACE(MeshManager,bbtk::AtomicBlackBox); +//===== +// Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost) +//===== + BBTK_DECLARE_INPUT(Mesh,vtkPolyData*); + BBTK_DECLARE_OUTPUT(MeshBase,vtkPolyData*); + BBTK_DECLARE_OUTPUT(MeshTemp,vtkPolyData*); + BBTK_DECLARE_OUTPUT(MeshManagerModel,MeshManagerModel*); + BBTK_PROCESS(Process); + void Process(); + + MeshManagerModel *meshManagerModel; +//===== +// Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost) +//===== +}; + +BBTK_BEGIN_DESCRIBE_BLACK_BOX(MeshManager,bbtk::AtomicBlackBox); + BBTK_NAME("MeshManager"); + BBTK_AUTHOR("InfoDev"); + BBTK_DESCRIPTION("No Description."); + BBTK_CATEGORY("empty"); + + BBTK_INPUT(MeshManager,Mesh,"Mesh",vtkPolyData*,""); + BBTK_OUTPUT(MeshManager,MeshBase,"Mesh Base",vtkPolyData*,""); + BBTK_OUTPUT(MeshManager,MeshTemp,"Mesh Temp",vtkPolyData*,""); + BBTK_OUTPUT(MeshManager,MeshManagerModel,"Mesh manager model",MeshManagerModel*,""); + +BBTK_END_DESCRIBE_BLACK_BOX(MeshManager); +//===== +// Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost) +//===== +} +// EO namespace bbcreaVtk + +#endif // __bbcreaVtkMeshManager_h_INCLUDED__ + diff --git a/bbtk_creaVtk_PKG/src/bbcreaVtkMeshManager_tool.cxx b/bbtk_creaVtk_PKG/src/bbcreaVtkMeshManager_tool.cxx new file mode 100644 index 0000000..cb12791 --- /dev/null +++ b/bbtk_creaVtk_PKG/src/bbcreaVtkMeshManager_tool.cxx @@ -0,0 +1,90 @@ +//===== +// Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost) +//===== +#include "bbcreaVtkMeshManager_tool.h" +#include "bbcreaVtkPackage.h" +namespace bbcreaVtk +{ + +BBTK_ADD_BLACK_BOX_TO_PACKAGE(creaVtk,MeshManager_tool) +BBTK_BLACK_BOX_IMPLEMENTATION(MeshManager_tool,bbtk::AtomicBlackBox); +//===== +// Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost) +//===== +void MeshManager_tool::Process() +{ +// THE MAIN PROCESSING METHOD BODY +// Here we simply set the input 'In' value to the output 'Out' +// And print out the output value +// INPUT/OUTPUT ACCESSORS ARE OF THE FORM : +// void bbSet{Input|Output}NAME(const TYPE&) +// const TYPE& bbGet{Input|Output}NAME() const +// Where : +// * NAME is the name of the input/output +// (the one provided in the attribute 'name' of the tag 'input') +// * TYPE is the C++ type of the input/output +// (the one provided in the attribute 'type' of the tag 'input') +// bbSetOutputOut( bbGetInputIn() ); +// std::cout << "Output value = " <SetMeshBase( bbGetInputMesh() ); + } + + if (bbGetInputTool()==4) // Reset + { + printf("EED Warning! MeshManager_tool Reset Not implemented.\n"); + } +} + +//===== +// Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost) +//===== +void MeshManager_tool::bbUserSetDefaultValues() +{ +// SET HERE THE DEFAULT INPUT/OUTPUT VALUES OF YOUR BOX +// Here we initialize the input 'In' to 0 + bbSetInputTool(0); + bbSetInputMesh(NULL); + bbSetInputMeshManagerModel(NULL); +} +//===== +// Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost) +//===== +void MeshManager_tool::bbUserInitializeProcessing() +{ + +// THE INITIALIZATION METHOD BODY : +// Here does nothing +// but this is where you should allocate the internal/output pointers +// if any + + +} +//===== +// Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost) +//===== +void MeshManager_tool::bbUserFinalizeProcessing() +{ + +// THE FINALIZATION METHOD BODY : +// Here does nothing +// but this is where you should desallocate the internal/output pointers +// if any + +} +} +// EO namespace bbcreaVtk + + diff --git a/bbtk_creaVtk_PKG/src/bbcreaVtkMeshManager_tool.h b/bbtk_creaVtk_PKG/src/bbcreaVtkMeshManager_tool.h new file mode 100644 index 0000000..f4ce72e --- /dev/null +++ b/bbtk_creaVtk_PKG/src/bbcreaVtkMeshManager_tool.h @@ -0,0 +1,54 @@ +//===== +// Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost) +//===== +#ifndef __bbcreaVtkMeshManager_tool_h_INCLUDED__ +#define __bbcreaVtkMeshManager_tool_h_INCLUDED__ + +#include "bbcreaVtk_EXPORT.h" +#include "bbtkAtomicBlackBox.h" +#include "iostream" +#include +#include + +namespace bbcreaVtk +{ + +class bbcreaVtk_EXPORT MeshManager_tool + : + public bbtk::AtomicBlackBox +{ + BBTK_BLACK_BOX_INTERFACE(MeshManager_tool,bbtk::AtomicBlackBox); +//===== +// Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost) +//===== + BBTK_DECLARE_INPUT(Tool,int); + BBTK_DECLARE_INPUT(Mesh,vtkPolyData*); + BBTK_DECLARE_INPUT(MeshManagerModel,MeshManagerModel*); +// BBTK_DECLARE_OUTPUT(Out,double); + BBTK_PROCESS(Process); + void Process(); +//===== +// Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost) +//===== +}; + +BBTK_BEGIN_DESCRIBE_BLACK_BOX(MeshManager_tool,bbtk::AtomicBlackBox); + BBTK_NAME("MeshManager_tool"); + BBTK_AUTHOR("InfoDev"); + BBTK_DESCRIPTION("No Description."); + BBTK_CATEGORY("empty"); + + BBTK_INPUT(MeshManager_tool,Tool,"(default 0) 0:Nothing 1:Undo 2:ReDo 3:Set 4:Reset",int,""); + BBTK_INPUT(MeshManager_tool,Mesh,"Mesh",vtkPolyData*,""); + BBTK_INPUT(MeshManager_tool,MeshManagerModel,"Mesh Manager Model",MeshManagerModel*,""); + +// BBTK_OUTPUT(MeshManager_tool,Out,"First output",double,""); +BBTK_END_DESCRIBE_BLACK_BOX(MeshManager_tool); +//===== +// Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost) +//===== +} +// EO namespace bbcreaVtk + +#endif // __bbcreaVtkMeshManager_tool_h_INCLUDED__ + diff --git a/lib/creaVtk/MeshManagerModel.cpp b/lib/creaVtk/MeshManagerModel.cpp new file mode 100644 index 0000000..bcf8c1a --- /dev/null +++ b/lib/creaVtk/MeshManagerModel.cpp @@ -0,0 +1,59 @@ +/* +# --------------------------------------------------------------------- +# +# Copyright (c) CREATIS (Centre de Recherche en Acquisition et Traitement de l'Image +# pour la Sante) +# Authors : Eduardo Davila, Frederic Cervenansky, Claire Mouton +# Previous Authors : Laurent Guigues, Jean-Pierre Roux +# CreaTools website : www.creatis.insa-lyon.fr/site/fr/creatools_accueil +# +# This software is governed by the CeCILL-B license under French law and +# abiding by the rules of distribution of free software. You can use, +# modify and/ or redistribute the software under the terms of the CeCILL-B +# license as circulated by CEA, CNRS and INRIA at the following URL +# http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html +# or in the file LICENSE.txt. +# +# As a counterpart to the access to the source code and rights to copy, +# modify and redistribute granted by the license, users are provided only +# with a limited warranty and the software's author, the holder of the +# economic rights, and the successive licensors have only limited +# liability. +# +# The fact that you are presently reading this means that you have had +# knowledge of the CeCILL-B license and that you accept its terms. +# ------------------------------------------------------------------------ +*/ + +#include "MeshManagerModel.h" + +MeshManagerModel::MeshManagerModel() +{ + _meshBase = NULL; + _meshTemp = NULL; +} + +MeshManagerModel::~MeshManagerModel() +{ +} + +void MeshManagerModel::SetMeshBase(vtkPolyData* mesh) +{ + _meshBase = mesh; + if (_meshTemp!=NULL) + { + _meshTemp->Delete(); + } // if + _meshTemp = vtkPolyData::New(); + _meshTemp->DeepCopy(_meshBase); +} + +vtkPolyData* MeshManagerModel::GetMeshBase() +{ + return _meshBase; +} + +vtkPolyData* MeshManagerModel::GetMeshTemp() +{ + return _meshTemp; +} diff --git a/lib/creaVtk/MeshManagerModel.h b/lib/creaVtk/MeshManagerModel.h new file mode 100644 index 0000000..70f33f8 --- /dev/null +++ b/lib/creaVtk/MeshManagerModel.h @@ -0,0 +1,70 @@ +/* +# --------------------------------------------------------------------- +# +# Copyright (c) CREATIS (Centre de Recherche en Acquisition et Traitement de l'Image +# pour la Sante) +# Authors : Eduardo Davila, Frederic Cervenansky, Claire Mouton +# Previous Authors : Laurent Guigues, Jean-Pierre Roux +# CreaTools website : www.creatis.insa-lyon.fr/site/fr/creatools_accueil +# +# This software is governed by the CeCILL-B license under French law and +# abiding by the rules of distribution of free software. You can use, +# modify and/ or redistribute the software under the terms of the CeCILL-B +# license as circulated by CEA, CNRS and INRIA at the following URL +# http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html +# or in the file LICENSE.txt. +# +# As a counterpart to the access to the source code and rights to copy, +# modify and redistribute granted by the license, users are provided only +# with a limited warranty and the software's author, the holder of the +# economic rights, and the successive licensors have only limited +# liability. +# +# The fact that you are presently reading this means that you have had +# knowledge of the CeCILL-B license and that you accept its terms. +# ------------------------------------------------------------------------ +*/ + +#ifndef _MESHMANAGERMODEL_H_ +#define _MESHMANAGERMODEL_H_ + +//--------------------------------------------- +// Class Name: MeshManagerModel +// [classdescription] +//--------------------------------------------- + +#include + +class MeshManagerModel +{ + +//--------------------------------------------- +//Methods and attributes exposed to other classes +//--------------------------------------------- +public : + MeshManagerModel(); + ~MeshManagerModel(); + + void SetMeshBase(vtkPolyData* mesh); + vtkPolyData* GetMeshBase(); + vtkPolyData* GetMeshTemp(); + +//--Method template---------------------------- +// void FunctionName(int& parameterA); + +//--------------------------------------------- +//Methods and attributes exposed only to classes +//that are derived from this class +//--------------------------------------------- +protected: + +//--------------------------------------------- +//Methods and attributes only visible by this class +//--------------------------------------------- +private: + vtkPolyData *_meshBase; + vtkPolyData *_meshTemp; +}; + +//-end of _MESHMANAGERMODEL_H_------------------------------------------------------ +#endif diff --git a/lib/creaVtk/Utilities.cxx b/lib/creaVtk/Utilities.cxx new file mode 100644 index 0000000..387bba0 --- /dev/null +++ b/lib/creaVtk/Utilities.cxx @@ -0,0 +1,239 @@ +/* +Copyright 2012-2022 Ronald Römer + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +#include "Utilities.h" + +#include + +#include +#include +#include +#include +#include +#include + +#include + +void ComputeNormal (vtkPoints *pts, double *n, vtkIdType num, const vtkIdType *poly) { + n[0] = 0; n[1] = 0; n[2] = 0; + + if (num == 3) { + double pA[3], pB[3], pC[3], a[3], b[3]; + + pts->GetPoint(poly[0], pA); + pts->GetPoint(poly[1], pB); + pts->GetPoint(poly[2], pC); + + vtkMath::Subtract(pB, pA, a); + vtkMath::Subtract(pC, pA, b); + + vtkMath::Cross(a, b, n); + } else { + double pA[3], pB[3]; + + vtkIdType i, a, b; + + for (i = 0; i < num; i++) { + a = poly[i]; + b = poly[i+1 == num ? 0 : i+1]; + + pts->GetPoint(a, pA); + pts->GetPoint(b, pB); + + n[0] += (pA[1]-pB[1])*(pA[2]+pB[2]); + n[1] += (pA[2]-pB[2])*(pA[0]+pB[0]); + n[2] += (pA[0]-pB[0])*(pA[1]+pB[1]); + } + } + + vtkMath::Normalize(n); +} + +void FindPoints (vtkKdTreePointLocator *pl, const double *pt, vtkIdList *pts, double tol) { + pts->Reset(); + + vtkPolyData *pd = vtkPolyData::SafeDownCast(pl->GetDataSet()); + + vtkIdList *closest = vtkIdList::New(); + + // vtkKdTree.cxx#L2505 + // arbeitet mit single-precision + pl->FindPointsWithinRadius(std::max(1e-3, tol), pt, closest); + + vtkIdType i, numPts = closest->GetNumberOfIds(); + + double c[3], v[3]; + + for (i = 0; i < numPts; i++) { + pd->GetPoint(closest->GetId(i), c); + vtkMath::Subtract(pt, c, v); + + if (vtkMath::Norm(v) < tol) { + pts->InsertNextId(closest->GetId(i)); + } + } + + closest->Delete(); +} + +void WriteVTK (const char *name, vtkPolyData *pd) { + vtkPolyDataWriter *w = vtkPolyDataWriter::New(); + w->SetInputData(pd); + w->SetFileName(name); + w->Update(); + w->Delete(); +} + +double GetAngle (double *vA, double *vB, double *n) { + // http://math.stackexchange.com/questions/878785/how-to-find-an-angle-in-range0-360-between-2-vectors + + double _vA[3]; + + vtkMath::Cross(n, vA, _vA); + double ang = std::atan2(vtkMath::Dot(_vA, vB), vtkMath::Dot(vA, vB)); + + if (ang < 0) { + ang += 2*M_PI; + } + + return ang; +} + +Base::Base (vtkPoints *pts, vtkIdType num, const vtkIdType *poly) { + ComputeNormal(pts, n, num, poly); + + double ptA[3], + ptB[3]; + + pts->GetPoint(poly[0], ptA); + pts->GetPoint(poly[1], ptB); + + ei[0] = ptB[0]-ptA[0]; + ei[1] = ptB[1]-ptA[1]; + ei[2] = ptB[2]-ptA[2]; + + vtkMath::Normalize(ei); + + ej[0] = n[1]*ei[2]-n[2]*ei[1]; + ej[1] = -n[0]*ei[2]+n[2]*ei[0]; + ej[2] = n[0]*ei[1]-n[1]*ei[0]; + + vtkMath::Normalize(ej); + + d = n[0]*ptA[0]+n[1]*ptA[1]+n[2]*ptA[2]; +} + +void Transform (const double *in, double *out, const Base &base) { + double x = base.ei[0]*in[0]+base.ei[1]*in[1]+base.ei[2]*in[2], + y = base.ej[0]*in[0]+base.ej[1]*in[1]+base.ej[2]*in[2]; + + out[0] = x; + out[1] = y; +} + +void BackTransform (const double *in, double *out, const Base &base) { + double x = in[0]*base.ei[0]+in[1]*base.ej[0]+base.d*base.n[0], + y = in[0]*base.ei[1]+in[1]*base.ej[1]+base.d*base.n[1], + z = in[0]*base.ei[2]+in[1]*base.ej[2]+base.d*base.n[2]; + + out[0] = x; + out[1] = y; + out[2] = z; +} + +void ComputeNormal (const Poly &poly, double *n) { + n[0] = 0; n[1] = 0; n[2] = 0; + + Poly::const_iterator itrA, itrB; + + for (itrA = poly.begin(); itrA != poly.end(); ++itrA) { + itrB = itrA+1; + + if (itrB == poly.end()) { + itrB = poly.begin(); + } + + const Point3d &ptA = *itrA, + &ptB = *itrB; + + n[0] += (ptA.y-ptB.y)*(ptA.z+ptB.z); + n[1] += (ptA.z-ptB.z)*(ptA.x+ptB.x); + n[2] += (ptA.x-ptB.x)*(ptA.y+ptB.y); + } + + vtkMath::Normalize(n); +} + +bool PointInPoly (const Poly &poly, const Point3d &p) { + bool in = false; + + Poly::const_iterator itrA, itrB; + + for (itrA = poly.begin(); itrA != poly.end(); ++itrA) { + itrB = itrA+1; + + if (itrB == poly.end()) { + itrB = poly.begin(); + } + + const Point3d &ptA = *itrA, + &ptB = *itrB; + + if ((ptA.x <= p.x || ptB.x <= p.x) + && ((ptA.y < p.y && ptB.y >= p.y) + || (ptB.y < p.y && ptA.y >= p.y))) { + + // schnittpunkt mit bounding box und strahlensatz + if (ptA.x+(p.y-ptA.y)*(ptB.x-ptA.x)/(ptB.y-ptA.y) < p.x) { + in = !in; + } + } + } + + return in; +} + +void WritePolys (const char *name, const PolysType &polys) { + vtkSmartPointer pts = vtkSmartPointer::New(); + + vtkSmartPointer pd = vtkSmartPointer::New(); + pd->SetPoints(pts); + pd->Allocate(1); + + for (auto &poly : polys) { + vtkSmartPointer cell = vtkSmartPointer::New(); + + for (auto &p : poly) { + cell->InsertNextId(pts->InsertNextPoint(p.x, p.y, p.z)); + } + + pd->InsertNextCell(VTK_POLYGON, cell); + } + + WriteVTK(name, pd); +} + +void GetPolys (const ReferencedPointsType &pts, const IndexedPolysType &indexedPolys, PolysType &polys) { + for (const auto &poly : indexedPolys) { + Poly newPoly; + + for (auto &id : poly) { + newPoly.push_back(pts.at(id)); + } + + polys.push_back(std::move(newPoly)); + } +} diff --git a/lib/creaVtk/Utilities.h b/lib/creaVtk/Utilities.h new file mode 100644 index 0000000..f46bfad --- /dev/null +++ b/lib/creaVtk/Utilities.h @@ -0,0 +1,141 @@ +/* +Copyright 2012-2022 Ronald Römer + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +#ifndef __Utilities_h +#define __Utilities_h + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#define NOTSET -1 + +double GetAngle (double *vA, double *vB, double *n); + +/* VTK */ +void ComputeNormal (vtkPoints *pts, double *n, vtkIdType num, const vtkIdType *poly); +void FindPoints (vtkKdTreePointLocator *pl, const double *pt, vtkIdList *pts, double tol = 1e-6); +void WriteVTK (const char *name, vtkPolyData *pd); + +class Point3d { +public: + const double x, y, z; + vtkIdType id; + + Point3d () = delete; + Point3d (const double _x, const double _y, const double _z, vtkIdType _id = NOTSET) : x(_x), y(_y), z(_z), id(_id) {} + bool operator< (const Point3d &other) const { + const long x1 = std::lround(x*1e5), + y1 = std::lround(y*1e5), + z1 = std::lround(z*1e5), + x2 = std::lround(other.x*1e5), + y2 = std::lround(other.y*1e5), + z2 = std::lround(other.z*1e5); + + return std::tie(x1, y1, z1) < std::tie(x2, y2, z2); + } + bool operator== (const Point3d &other) const { + return !(*this < other) && !(other < *this); + } + + friend std::ostream& operator<< (std::ostream &out, const Point3d &p) { + out << "Point3d(x=" << p.x + << ", y=" << p.y + << ", z=" << p.z + << ", id=" << p.id << ")"; + return out; + } +}; + +typedef std::vector IdsType; + +template::value, bool>::type = true> +class _Pair { +public: + T f, g; + _Pair () = delete; + _Pair (T _f, T _g) : f(_f), g(_g) {} + bool operator< (const _Pair &other) const { + return std::tie(f, g) < std::tie(other.f, other.g); + } + bool operator== (const _Pair &other) const { + return f == other.f && g == other.g; + } + friend std::ostream& operator<< (std::ostream &out, const _Pair &p) { + out << "(" << p.f << ", " << p.g << ")"; + return out; + } +}; + +// typedef _Pair Pair; +using Pair = _Pair; + +template::value + && std::is_integral::value + && std::is_signed::value, bool>::type = true> +T Mod (T a, U b) { + T _b = static_cast(b); + + return ((a%_b)+_b)%_b; +} + +class Base { +public: + Base () {} + Base (vtkPoints *pts, vtkIdType num, const vtkIdType *poly); + double n[3], ei[3], ej[3], d; +}; + +void Transform (const double *in, double *out, const Base &base); +void BackTransform (const double *in, double *out, const Base &base); + +inline void Cpy (double *a, const double *b, const int n = 2) { + std::copy_n(b, n, a); +} + +template +std::ostream& operator<< (typename std::enable_if::value, std::ostream>::type& stream, const T& e) { + return stream << static_cast::type>(e); +} + +typedef std::vector Poly; +typedef std::vector PolysType; + +void ComputeNormal (const Poly &poly, double *n); +bool PointInPoly (const Poly &poly, const Point3d &p); + +void WritePolys (const char *name, const PolysType &polys); + +typedef std::deque IndexedPoly; +typedef std::vector IndexedPolysType; + +typedef std::map> ReferencedPointsType; + +void GetPolys (const ReferencedPointsType &pts, const IndexedPolysType &indexedPolys, PolysType &polys); + +#endif diff --git a/lib/creaVtk/vtkPolyDataBooleanFilter.cxx b/lib/creaVtk/vtkPolyDataBooleanFilter.cxx new file mode 100644 index 0000000..ea6a9a2 --- /dev/null +++ b/lib/creaVtk/vtkPolyDataBooleanFilter.cxx @@ -0,0 +1,3376 @@ +/* +Copyright 2012-2022 Ronald Römer + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "vtkPolyDataBooleanFilter.h" +#include "vtkPolyDataContactFilter.h" + +#include "Utilities.h" + +vtkStandardNewMacro(vtkPolyDataBooleanFilter); + +vtkPolyDataBooleanFilter::vtkPolyDataBooleanFilter () { + + SetNumberOfInputPorts(2); + SetNumberOfOutputPorts(3); + + timePdA = 0; + timePdB = 0; + + contLines = vtkSmartPointer::New(); + + modPdA = vtkSmartPointer::New(); + modPdB = vtkSmartPointer::New(); + + cellDataA = vtkSmartPointer::New(); + cellDataB = vtkSmartPointer::New(); + + cellIdsA = vtkSmartPointer::New(); + cellIdsB = vtkSmartPointer::New(); + + OperMode = OPER_UNION; + +} + +vtkPolyDataBooleanFilter::~vtkPolyDataBooleanFilter () { + // nix mehr +} + +int vtkPolyDataBooleanFilter::ProcessRequest(vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector) { + + if (request->Has(vtkDemandDrivenPipeline::REQUEST_DATA())) { + + vtkInformation *inInfoA = inputVector[0]->GetInformationObject(0); + vtkInformation *inInfoB = inputVector[1]->GetInformationObject(0); + + vtkPolyData *pdA = vtkPolyData::SafeDownCast(inInfoA->Get(vtkDataObject::DATA_OBJECT())); + vtkPolyData *pdB = vtkPolyData::SafeDownCast(inInfoB->Get(vtkDataObject::DATA_OBJECT())); + + vtkInformation *outInfoA = outputVector->GetInformationObject(0); + vtkInformation *outInfoB = outputVector->GetInformationObject(1); + vtkInformation *outInfoC = outputVector->GetInformationObject(2); + + resultA = vtkPolyData::SafeDownCast(outInfoA->Get(vtkDataObject::DATA_OBJECT())); + resultB = vtkPolyData::SafeDownCast(outInfoB->Get(vtkDataObject::DATA_OBJECT())); + resultC = vtkPolyData::SafeDownCast(outInfoC->Get(vtkDataObject::DATA_OBJECT())); + + using clock = std::chrono::steady_clock; + std::vector times; + clock::time_point start; + + if (pdA->GetMTime() > timePdA || pdB->GetMTime() > timePdB) { + // eventuell vorhandene regionen vereinen + + vtkSmartPointer cleanA = vtkSmartPointer::New(); + cleanA->SetOutputPointsPrecision(DOUBLE_PRECISION); + cleanA->SetTolerance(1e-6); + cleanA->SetInputData(pdA); + cleanA->Update(); + + vtkSmartPointer cleanB = vtkSmartPointer::New(); + cleanB->SetOutputPointsPrecision(DOUBLE_PRECISION); + cleanB->SetTolerance(1e-6); + cleanB->SetInputData(pdB); + cleanB->Update(); + +#ifdef DEBUG + std::cout << "Exporting modPdA.vtk" << std::endl; + WriteVTK("modPdA.vtk", cleanA->GetOutput()); + + std::cout << "Exporting modPdB.vtk" << std::endl; + WriteVTK("modPdB.vtk", cleanB->GetOutput()); +#endif + + // CellData sichern + + cellDataA->DeepCopy(cleanA->GetOutput()->GetCellData()); + cellDataB->DeepCopy(cleanB->GetOutput()->GetCellData()); + + // ermittelt kontaktstellen + + start = clock::now(); + + vtkSmartPointer cl = vtkSmartPointer::New(); + cl->SetInputConnection(0, cleanA->GetOutputPort()); + cl->SetInputConnection(1, cleanB->GetOutputPort()); + cl->Update(); + + times.push_back(clock::now()-start); + + contLines->DeepCopy(cl->GetOutput()); + + modPdA->DeepCopy(cl->GetOutput(1)); + modPdB->DeepCopy(cl->GetOutput(2)); + +#ifdef DEBUG + std::cout << "Exporting contLines.vtk" << std::endl; + WriteVTK("contLines.vtk", contLines); + + std::cout << "Exporting modPdA_1.vtk" << std::endl; + WriteVTK("modPdA_1.vtk", modPdA); + + std::cout << "Exporting modPdB_1.vtk" << std::endl; + WriteVTK("modPdB_1.vtk", modPdB); +#endif + + if (contLines->GetNumberOfCells() == 0) { + vtkErrorMacro("Inputs have no contact."); + + return 1; + } + + // in den CellDatas steht drin, welche polygone einander schneiden + + contsA = vtkIdTypeArray::SafeDownCast(contLines->GetCellData()->GetScalars("cA")); + contsB = vtkIdTypeArray::SafeDownCast(contLines->GetCellData()->GetScalars("cB")); + + vtkIdTypeArray *accuracy = vtkIdTypeArray::SafeDownCast(contLines->GetPointData()->GetScalars("accuracy")); + + vtkIdTypeArray *sourcesA = vtkIdTypeArray::SafeDownCast(contLines->GetCellData()->GetScalars("sourcesA")); + vtkIdTypeArray *sourcesB = vtkIdTypeArray::SafeDownCast(contLines->GetCellData()->GetScalars("sourcesB")); + + vtkIdType i, numPts = contLines->GetNumberOfPoints(); + + vtkSmartPointer cells = vtkSmartPointer::New(); + + for (i = 0; i < numPts; i++) { + contLines->GetPointCells(i, cells); + + if (cells->GetNumberOfIds() == 1) { + vtkErrorMacro("Contact ends suddenly."); + return 1; + } + + if (accuracy->GetValue(i) == 1) { + vtkErrorMacro("Contact goes through a cell with a bad shape."); + return 1; + } + } + + // sichert die OrigCellIds + + vtkIdTypeArray *origCellIdsA = vtkIdTypeArray::SafeDownCast(modPdA->GetCellData()->GetScalars("OrigCellIds")); + vtkIdTypeArray *origCellIdsB = vtkIdTypeArray::SafeDownCast(modPdB->GetCellData()->GetScalars("OrigCellIds")); + + cellIdsA->DeepCopy(origCellIdsA); + cellIdsB->DeepCopy(origCellIdsB); + + vtkIdType numCellsA = modPdA->GetNumberOfCells(), + numCellsB = modPdB->GetNumberOfCells(); + + for (i = 0; i < numCellsA; i++) { + origCellIdsA->SetValue(i, i); + } + + for (i = 0; i < numCellsB; i++) { + origCellIdsB->SetValue(i, i); + } + + start = clock::now(); + + if (GetPolyStrips(modPdA, contsA, sourcesA, polyStripsA) || + GetPolyStrips(modPdB, contsB, sourcesB, polyStripsB)) { + + vtkErrorMacro("Strips are invalid."); + + return 1; + + } + + times.push_back(clock::now()-start); + + // trennt die polygone an den linien + + start = clock::now(); + + if (CutCells(modPdA, polyStripsA) || + CutCells(modPdB, polyStripsB)) { + + vtkErrorMacro("CutCells failed."); + + return 1; + } + + times.push_back(clock::now()-start); + +#ifdef DEBUG + std::cout << "Exporting modPdA_2.vtk" << std::endl; + WriteVTK("modPdA_2.vtk", modPdA); + + std::cout << "Exporting modPdB_2.vtk" << std::endl; + WriteVTK("modPdB_2.vtk", modPdB); +#endif + + start = clock::now(); + + RestoreOrigPoints(modPdA, polyStripsA); + RestoreOrigPoints(modPdB, polyStripsB); + + times.push_back(clock::now()-start); + +#ifdef DEBUG + std::cout << "Exporting modPdA_3.vtk" << std::endl; + WriteVTK("modPdA_3.vtk", modPdA); + + std::cout << "Exporting modPdB_3.vtk" << std::endl; + WriteVTK("modPdB_3.vtk", modPdB); +#endif + + start = clock::now(); + + ResolveOverlaps(modPdA, contsA, polyStripsA); + ResolveOverlaps(modPdB, contsB, polyStripsB); + + times.push_back(clock::now()-start); + +#ifdef DEBUG + std::cout << "Exporting modPdA_4.vtk" << std::endl; + WriteVTK("modPdA_4.vtk", modPdA); + + std::cout << "Exporting modPdB_4.vtk" << std::endl; + WriteVTK("modPdB_4.vtk", modPdB); +#endif + + start = clock::now(); + + AddAdjacentPoints(modPdA, contsA, polyStripsA); + AddAdjacentPoints(modPdB, contsB, polyStripsB); + + times.push_back(clock::now()-start); + +#ifdef DEBUG + std::cout << "Exporting modPdA_5.vtk" << std::endl; + WriteVTK("modPdA_5.vtk", modPdA); + + std::cout << "Exporting modPdB_5.vtk" << std::endl; + WriteVTK("modPdB_5.vtk", modPdB); +#endif + + start = clock::now(); + + DisjoinPolys(modPdA, polyStripsA); + DisjoinPolys(modPdB, polyStripsB); + + times.push_back(clock::now()-start); + +#ifdef DEBUG + std::cout << "Exporting modPdA_6.vtk" << std::endl; + WriteVTK("modPdA_6.vtk", modPdA); + + std::cout << "Exporting modPdB_6.vtk" << std::endl; + WriteVTK("modPdB_6.vtk", modPdB); +#endif + + start = clock::now(); + + MergePoints(modPdA, polyStripsA); + MergePoints(modPdB, polyStripsB); + + times.push_back(clock::now()-start); + +#ifdef DEBUG + std::cout << "Exporting modPdA_7.vtk" << std::endl; + WriteVTK("modPdA_7.vtk", modPdA); + + std::cout << "Exporting modPdB_7.vtk" << std::endl; + WriteVTK("modPdB_7.vtk", modPdB); +#endif + + timePdA = pdA->GetMTime(); + timePdB = pdB->GetMTime(); + + } + + start = clock::now(); + + if (CombineRegions()) { + vtkErrorMacro("Boolean operation failed."); + + return 1; + } + + times.push_back(clock::now()-start); + +// #ifdef DEBUG + double sum = std::chrono::duration_cast>(std::accumulate(times.begin(), times.end(), clock::duration())).count(); + + std::vector::const_iterator itr; + for (itr = times.begin(); itr != times.end(); itr++) { + double time = std::chrono::duration_cast>(*itr).count(); + + std::cout << "Time " << (itr-times.begin()) + << ": " << time << "s (" << (time/sum*100) << "%)" + << std::endl; + } +// #endif + + } + + return 1; + +} + +void vtkPolyDataBooleanFilter::GetStripPoints (vtkPolyData *pd, vtkIdTypeArray *sources, PStrips &pStrips, IdsType &lines) { + +#ifdef DEBUG + std::cout << "GetStripPoints()" << std::endl; +#endif + + StripPtsType &pts = pStrips.pts; + const IdsType &poly = pStrips.poly; + + double a[3], b[3], sA[3], sB[3], u[3], v[3], w[3], n, t, d; + + std::map allPts; + + vtkSmartPointer line = vtkSmartPointer::New(); + + for (auto lineId : lines) { + contLines->GetCellPoints(lineId, line); + + // std::cout << "? " << contsA->GetValue(lineId) + // << ", " << contsB->GetValue(lineId) + // << ", [" << line->GetId(0) + // << ", " << line->GetId(1) + // << "]" + // << std::endl; + + allPts.emplace(line->GetId(0), sources->GetTypedComponent(lineId, 0)); + allPts.emplace(line->GetId(1), sources->GetTypedComponent(lineId, 1)); + } + + decltype(allPts)::const_iterator itr; + + for (itr = allPts.begin(); itr != allPts.end(); ++itr) { + StripPt sp; + sp.ind = itr->first; + + // die koordinaten + double pt[3]; + contLines->GetPoint(sp.ind, pt); + + Cpy(sp.pt, pt, 3); + + IdsType::const_iterator itrA, itrB; + + for (itrA = poly.begin(); itrA != poly.end(); ++itrA) { + itrB = itrA+1; + + if (itrB == poly.end()) { + itrB = poly.begin(); + } + + if (itr->second != NOTSET && *itrA != itr->second) { + continue; + } + + pd->GetPoint(*itrA, a); + pd->GetPoint(*itrB, b); + + vtkMath::Subtract(a, pt, sA); + vtkMath::Subtract(b, pt, sB); + + // richtungsvektor und länge der kante + + vtkMath::Subtract(b, a, u); + n = vtkMath::Norm(u); + + // d und t zur kante + + vtkMath::Subtract(pt, a, v); + t = vtkMath::Dot(v, u)/(n*n); + + vtkMath::Cross(v, u, w); + d = vtkMath::Norm(w)/n; + + if (d < 1e-5 && t > -1e-5 && t < 1+1e-5) { + sp.edge[0] = *itrA; + sp.edge[1] = *itrB; + + sp.t = std::min(1., std::max(0., t)); + + if (vtkMath::Norm(sA) < 1e-5) { + Cpy(sp.captPt, a, 3); + sp.capt = Capt::A; + + } else if (vtkMath::Norm(sB) < 1e-5) { + Cpy(sp.captPt, b, 3); + sp.capt = Capt::B; + + } else { + // u ist nicht normiert + vtkMath::MultiplyScalar(u, t); + + double x[3]; + vtkMath::Add(a, u, x); + + // projektion + Cpy(sp.captPt, x, 3); + + sp.capt = Capt::EDGE; + + } + } + + // std::cout << "? " + // << sp.ind + // << ", " << d + // << ", " << t + // << ", " << sp.capt + // << std::endl; + } + + if (itr->second != NOTSET && sp.edge[0] == NOTSET) { + sp.catched = false; + } + + pts.emplace(sp.ind, std::move(sp)); + + } + + StripPtsType::iterator itr2; + + for (itr2 = pts.begin(); itr2 != pts.end(); ++itr2) { + StripPt &sp = itr2->second; + + if (sp.capt != Capt::NOT) { + if (sp.capt == Capt::B) { + sp.t = 0; + + sp.edge[0] = sp.edge[1]; + + auto itrA = std::find(poly.begin(), poly.end(), sp.edge[0]), + itrB = itrA+1; + + if (itrB == poly.end()) { + itrB = poly.begin(); + } + + sp.edge[1] = *itrB; + + sp.capt = Capt::A; + + } + + // für den schnitt werden die eingerasteten koordinaten verwendet + + Cpy(sp.cutPt, sp.captPt, 3); + } else { + + Cpy(sp.cutPt, sp.pt, 3); + } + + sp.history.emplace_back(sp.edge[0], sp.edge[1]); + + } + +#ifdef DEBUG + std::cout << "pts: " << std::endl; + for (itr2 = pts.begin(); itr2 != pts.end(); ++itr2) { + std::cout << itr2->first << ": " << itr2->second << std::endl; + } +#endif + +} + +bool vtkPolyDataBooleanFilter::GetPolyStrips (vtkPolyData *pd, vtkIdTypeArray *conts, vtkIdTypeArray *sources, PolyStripsType &polyStrips) { +#ifdef DEBUG + std::cout << "GetPolyStrips()" << std::endl; +#endif + + polyStrips.clear(); + + std::map polyLines; + + for (vtkIdType i = 0; i < conts->GetNumberOfTuples(); i++) { + vtkIdType poly = conts->GetValue(i); + + /*if (poly != 1641) { + continue; + }*/ + + polyLines[poly].push_back(i); + } + + std::vector> notCatched; + + std::map::iterator itr; + + for (itr = polyLines.begin(); itr != polyLines.end(); ++itr) { + + IdsType &lines = itr->second; + RemoveDuplicates(lines); + + polyStrips.emplace(std::piecewise_construct, + std::forward_as_tuple(itr->first), + std::forward_as_tuple(pd, itr->first)); + + PStrips &pStrips = polyStrips.at(itr->first); + + GetStripPoints(pd, sources, pStrips, lines); + + for (auto &sp : pStrips.pts) { + sp.second.polyId = itr->first; + + if (!sp.second.catched) { + notCatched.push_back(sp.second); + } + } + + } + + auto Next = [](const IdsType &ids, vtkIdType id) -> vtkIdType { + IdsType::const_iterator itr; + + itr = std::find(ids.begin(), ids.end(), id); + + if (++itr == ids.end()) { + itr = ids.begin(); + } + + return *itr; + }; + + for (StripPt &sp : notCatched) { + for (itr = polyLines.begin(); itr != polyLines.end(); ++itr) { + const PStrips &pStrips = polyStrips.at(itr->first); + + try { + const StripPt &corr = pStrips.pts.at(sp.ind); + + if (&corr != &sp) { + if (corr.capt == Capt::A) { + sp.capt = Capt::A; + sp.edge[0] = corr.edge[0]; + sp.edge[1] = Next(polyStrips.at(sp.polyId).poly, sp.edge[0]); + + sp.t = 0; + + Cpy(sp.captPt, corr.captPt, 3); + Cpy(sp.cutPt, sp.captPt, 3); + + sp.history.emplace_back(sp.edge[0], sp.edge[1]); + + sp.catched = true; + + } + } + } catch (...) {} + + } + + if (!sp.catched) { + std::cout << sp << std::endl; + } + + assert(sp.catched); + + } + + vtkSmartPointer cells = vtkSmartPointer::New(), + line = vtkSmartPointer::New(); + + vtkIdType i, numCells; + + StripPtsType::const_iterator itr2; + + for (itr = polyLines.begin(); itr != polyLines.end(); ++itr) { + PStrips &pStrips = polyStrips.at(itr->first); + + const IdsType &lines = itr->second; + const StripPtsType &pts = pStrips.pts; + + for (itr2 = pts.begin(); itr2 != pts.end(); ++itr2) { + // ist der punkt auf einer kante, dürfen von ihm mehr als 2 linien ausgehen + + const StripPt &pt = itr2->second; + + if (pt.capt == Capt::NOT) { + contLines->GetPointCells(pt.ind, cells); + + numCells = cells->GetNumberOfIds(); + + std::set ends; + + for (i = 0; i < numCells; i++) { + contLines->GetCellPoints(cells->GetId(i), line); + + ends.insert(pt.ind == line->GetId(0) ? line->GetId(1) : line->GetId(0)); + } + + if (ends.size() > 2) { + return true; + } + } + } + + // zusammensetzen + + std::deque _lines; + + vtkIdList *linePts = vtkIdList::New(); + + for (auto &i : lines) { + contLines->GetCellPoints(i, linePts); + _lines.emplace_back(linePts->GetId(0), linePts->GetId(1)); + } + + linePts->Delete(); + + decltype(_lines)::iterator _itr; + + auto FindRight = [&pts, &_lines, &_itr](StripType &strip, const std::size_t &id) -> bool { + auto &right = strip.back(); + + if (pts.at(right.ind).capt == Capt::NOT) { + for (_itr = _lines.begin(); _itr != _lines.end(); ++_itr) { + if (_itr->f == right.ind) { + strip.emplace_back(_itr->g, id); + + _lines.erase(_itr); + return true; + } else if (_itr->g == right.ind) { + strip.emplace_back(_itr->f, id); + + _lines.erase(_itr); + return true; + } + } + } + + return false; + }; + + auto FindLeft = [&pts, &_lines, &_itr](StripType &strip, const std::size_t &id) -> bool { + auto &left = strip.front(); + + if (pts.at(left.ind).capt == Capt::NOT) { + for (_itr = _lines.begin(); _itr != _lines.end(); ++_itr) { + if (_itr->f == left.ind) { + strip.emplace_front(_itr->g, id); + + _lines.erase(_itr); + return true; + } else if (_itr->g == left.ind) { + strip.emplace_front(_itr->f, id); + + _lines.erase(_itr); + return true; + } + } + } + + return false; + }; + + std::size_t stripId {0}; + + while (!_lines.empty()) { + auto &last = _lines.back(); + + StripType strip {{last.f, stripId}, {last.g, stripId}}; + _lines.pop_back(); + + while (FindRight(strip, stripId)) {} + while (FindLeft(strip, stripId)) {} + + pStrips.strips.push_back(std::move(strip)); + + stripId++; + + } + + CompleteStrips(pStrips); + + } + + // sucht nach schnitten zw. den strips + + { + + PolyStripsType::const_iterator itr; + StripType::const_iterator itr2; + + for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) { + const PStrips &pStrips = itr->second; + + const StripsType &strips = pStrips.strips; + const StripPtsType &pts = pStrips.pts; + const Base &base = pStrips.base; + + vtkSmartPointer treePts = vtkSmartPointer::New(); + + vtkSmartPointer treePd = vtkSmartPointer::New(); + treePd->Allocate(1); + + std::map ptIds; + + double pt[2]; + + for (const auto &p : pts) { + Transform(p.second.pt, pt, base); + + ptIds.emplace(p.first, treePts->InsertNextPoint(pt[0], pt[1], 0)); + } + + for (const StripType &strip : strips) { + for (itr2 = strip.begin(); itr2 != strip.end()-1; ++itr2) { + + vtkIdList *line = vtkIdList::New(); + line->InsertNextId(ptIds[itr2->ind]); + line->InsertNextId(ptIds[(itr2+1)->ind]); + + treePd->InsertNextCell(VTK_LINE, line); + + line->Delete(); + } + } + + treePd->SetPoints(treePts); + + vtkSmartPointer tree = vtkSmartPointer::New(); + tree->SetDataSet(treePd); + tree->BuildLocator(); + + vtkIdType numA, numB; + const vtkIdType *lineA, *lineB; + + auto lineItr = vtk::TakeSmartPointer(treePd->GetLines()->NewIterator()); + + for (lineItr->GoToFirstCell(); !lineItr->IsDoneWithTraversal(); lineItr->GoToNextCell()) { + lineItr->GetCurrentCell(numA, lineA); + + double ptA[3], ptB[3]; + + treePts->GetPoint(lineA[0], ptA); + treePts->GetPoint(lineA[1], ptB); + + vtkSmartPointer lineIds = vtkSmartPointer::New(); + + tree->IntersectWithLine(ptA, ptB, 1e-5, nullptr, lineIds); + + for (vtkIdType i = 0; i < lineIds->GetNumberOfIds(); i++) { + treePd->GetCellPoints(lineIds->GetId(i), numB, lineB); + + if (lineB[0] != lineA[0] && lineB[1] != lineA[0] && lineB[0] != lineA[1] && lineB[1] != lineA[1]) { + // schnitt gefunden + + return true; + } + } + } + + } + + } + + return false; + +} + +void vtkPolyDataBooleanFilter::RemoveDuplicates (IdsType &lines) { + + IdsType unique; + + // die indexe der enden auf übereinstimmung prüfen + + vtkIdList *linePtsA = vtkIdList::New(); + vtkIdList *linePtsB = vtkIdList::New(); + + std::size_t i, j, numLines = lines.size(); + + for (i = 0; i < numLines-1; i++) { + j = i+1; + + contLines->GetCellPoints(lines[i], linePtsA); + + while (j < lines.size()) { + contLines->GetCellPoints(lines[j], linePtsB); + + if ((linePtsA->GetId(0) == linePtsB->GetId(0) && linePtsA->GetId(1) == linePtsB->GetId(1)) || + (linePtsA->GetId(0) == linePtsB->GetId(1) && linePtsA->GetId(1) == linePtsB->GetId(0))) { + // stimmen überein + break; + } + + j++; + } + + if (j == numLines) { + // keine vorzeitige unterbrechung der while-schleife + unique.push_back(lines[i]); + } + } + + unique.push_back(lines.back()); + + linePtsA->Delete(); + linePtsB->Delete(); + + lines.swap(unique); + +} + +void vtkPolyDataBooleanFilter::CompleteStrips (PStrips &pStrips) { + StripsType::iterator itr; + + for (itr = pStrips.strips.begin(); itr != pStrips.strips.end(); ++itr) { + const StripPt &start = pStrips.pts[itr->front().ind], + &end = pStrips.pts[itr->back().ind]; + + if (start.ind != end.ind) { + if (start.capt == Capt::NOT) { + StripType s(itr->rbegin(), itr->rend()-1); + itr->insert(itr->begin(), s.begin(), s.end()); + + } else if (end.capt == Capt::NOT) { + StripType s(itr->rbegin()+1, itr->rend()); + itr->insert(itr->end(), s.begin(), s.end()); + + } + } + } +} + +bool vtkPolyDataBooleanFilter::HasArea (const StripType &strip) const { + bool area = true; + + std::size_t i, n = strip.size(); + + if (n%2 == 1) { + for (i = 0; i < (n-1)/2; i++) { + area = strip[i].ind != strip[n-i-1].ind; + } + } + + return area; +} + +void ComputeNormal (const StripPtsType &pts, const RefsType &poly, double *n) { + n[0] = 0; n[1] = 0; n[2] = 0; + + RefsType::const_iterator itrA, itrB; + + for (itrA = poly.begin(); itrA != poly.end(); ++itrA) { + itrB = itrA+1; + + if (itrB == poly.end()) { + itrB = poly.begin(); + } + + const StripPtR &spA = *itrA, + &spB = *itrB; + + auto pA = pts.find(spA.ind); + auto pB = pts.find(spB.ind); + + const double *ptA = pA->second.cutPt, + *ptB = pB->second.cutPt; + + n[0] += (ptA[1]-ptB[1])*(ptA[2]+ptB[2]); + n[1] += (ptA[2]-ptB[2])*(ptA[0]+ptB[0]); + n[2] += (ptA[0]-ptB[0])*(ptA[1]+ptB[1]); + } + + vtkMath::Normalize(n); +} + +bool vtkPolyDataBooleanFilter::CutCells (vtkPolyData *pd, PolyStripsType &polyStrips) { +#ifdef DEBUG + std::cout << "CutCells()" << std::endl; +#endif + + vtkPoints *pdPts = pd->GetPoints(); + + vtkIdTypeArray *origCellIds = vtkIdTypeArray::SafeDownCast(pd->GetCellData()->GetScalars("OrigCellIds")); + + PolyStripsType::iterator itr; + + for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) { + + vtkIdType polyInd = itr->first; + PStrips &pStrips = itr->second; + + /*if (polyInd != 1641) { + continue; + }*/ + + StripsType &strips = pStrips.strips; + StripPtsType &pts = pStrips.pts; + + IdsType &poly = pStrips.poly; + + vtkIdType origId = origCellIds->GetValue(polyInd); + +#ifdef DEBUG + std::cout << "polyInd " << polyInd << ", poly ["; + for (auto &p : poly) { + std::cout << p << ", "; + } + std::cout << "]" << std::endl; +#endif + + std::size_t numPts = poly.size(); + + std::map edges; + + StripsType::iterator itr2; + StripType::iterator itr3; + + // holes sichern + StripsType holes; + + auto fct = [&](const StripType &s) { + return pts[s.front().ind].capt == Capt::NOT && pts[s.back().ind].capt == Capt::NOT; + }; + + std::copy_if(strips.begin(), strips.end(), std::back_inserter(holes), fct); + + strips.erase(std::remove_if(strips.begin(), strips.end(), fct), strips.end()); + + for (itr2 = strips.begin(); itr2 != strips.end(); ++itr2) { + StripType &strip = *itr2; + +#ifdef DEBUG + std::cout << (itr2-strips.begin()) << ". strip ["; + for (auto &s : strip) { + std::cout << s.ind << ", "; + } + std::cout << "]" << std::endl; +#endif + + // init + if (pts[strip.front().ind].edge[0] == pts[strip.back().ind].edge[0] + && strip.front().ind != strip.back().ind + && pts[strip.front().ind].t > pts[strip.back().ind].t) { + + std::reverse(strip.begin(), strip.end()); + } + + StripPt &start = pts[strip.front().ind], + &end = pts[strip.back().ind]; + + strip.front().side = Side::START; + strip.back().side = Side::END; + + strip.front().ref = start.edge[0]; + strip.back().ref = end.edge[0]; + + // nachfolgend könnte man dann anfang und ende weglassen + + for (itr3 = strip.begin(); itr3 != strip.end(); ++itr3) { + StripPt &sp = pts[itr3->ind]; + + itr3->desc[0] = pdPts->InsertNextPoint(sp.cutPt); + itr3->desc[1] = pdPts->InsertNextPoint(sp.cutPt); + +#ifdef DEBUG + std::cout << sp << " => " << *itr3 << std::endl; +#endif + + } + + // ordnet zu + + edges[start.edge[0]].push_back(std::ref(strip.front())); + edges[end.edge[0]].push_back(std::ref(strip.back())); + } + + // sortiert die punkte auf den kanten + + std::map::iterator itr4; + + IdsType::iterator itr5; + StripType::reverse_iterator itr6; + + RefsType::iterator itr7; + + for (itr4 = edges.begin(); itr4 != edges.end(); ++itr4) { + RefsType &edge = itr4->second; + +#ifdef DEBUG + std::cout << "edge (" << itr4->first << ", _)" << std::endl; +#endif + + std::sort(edge.begin(), edge.end(), [&](const StripPtR &a, const StripPtR &b) { + const StripPt &a_ = pts[a.ind], + &b_ = pts[b.ind]; + +#ifdef DEBUG + // std::cout << "a_: " << a_ << " -> strip " << a.strip << std::endl; + // std::cout << "b_: " << b_ << " -> strip " << b.strip << std::endl; +#endif + + if (a_.ind == b_.ind) { + // strips beginnen im gleichen punkt + + if (a.strip != b.strip) { + // gehören nicht dem gleichen strip an + + StripType &stripA = strips[a.strip], + &stripB = strips[b.strip]; + + // andere enden ermitteln + + const StripPtR &eA = (&a == &(stripA.front())) ? stripA.back() : stripA.front(), + &eB = (&b == &(stripB.front())) ? stripB.back() : stripB.front(); + + const StripPt &eA_ = pts[eA.ind], + &eB_ = pts[eB.ind]; + +#ifdef DEBUG + // std::cout << "eA_: " << eA_ << std::endl; + // std::cout << "eB_: " << eA_ << std::endl; +#endif + + if (eA_.ind != eB_.ind) { + auto i = std::find(poly.begin(), poly.end(), itr4->first)-poly.begin(), + iA = std::find(poly.begin(), poly.end(), eA_.edge[0])-poly.begin(), + iB = std::find(poly.begin(), poly.end(), eB_.edge[0])-poly.begin(); + + double dA = static_cast(Mod(iA-i, numPts))+eA_.t, + dB = static_cast(Mod(iB-i, numPts))+eB_.t; + + if (i == iA && a_.t > eA_.t) { + dA += static_cast(numPts); + } + + if (i == iB && b_.t > eB_.t) { + dB += static_cast(numPts); + } + +#ifdef DEBUG + // std::cout << "dA=" << dA << ", dB=" << dB << std::endl; +#endif + + return dB < dA; + } else { + RefsType poly_; + + if (a.side == Side::START) { + poly_.insert(poly_.end(), stripA.begin(), stripA.end()); + } else { + poly_.insert(poly_.end(), stripA.rbegin(), stripA.rend()); + } + + if (b.side == Side::START) { + poly_.insert(poly_.end(), stripB.rbegin()+1, stripB.rend()-1); + } else { + poly_.insert(poly_.end(), stripB.begin()+1, stripB.end()-1); + } + + double n[3]; + ComputeNormal(pts, poly_, n); + + double ang = vtkMath::Dot(pStrips.n, n); + +#ifdef DEBUG + // std::cout << "ang=" << ang*180/PI << std::endl; +#endif + + return ang < .999999; + + } + } else { + // gleicher strip + + StripType &strip = strips[a.strip]; + + if (HasArea(strip)) { + RefsType poly_(strip.begin(), strip.end()-1); + + double n[3]; + ComputeNormal(pts, poly_, n); + + double ang = vtkMath::Dot(pStrips.n, n); + + if (ang > .999999) { + std::reverse(strip.begin(), strip.end()); + return true; + } else { + return false; + } + + } else { + // reihenfolge von a und b bereits richtig + return false; + } + } + + } else { + return a_.t < b_.t; + } + }); + +#ifdef DEBUG + for (auto& p : edge) { + std::cout << p << std::endl; + } +#endif + + } + + // baut die strips ein + + std::deque polys; + polys.push_back(pStrips.poly); + + for (itr2 = strips.begin(); itr2 != strips.end(); ++itr2) { + StripType &strip = *itr2; + + StripPtR &start = strip.front(), + &end = strip.back(); + +#ifdef DEBUG + std::cout << "strip " << start.strip + << " , refs (" << start.ref << ", " << end.ref << ")" + << std::endl; +#endif + + std::size_t cycle = 0; + + while (true) { + if (cycle == polys.size()) { + break; + } + + IdsType next = polys.front(); + polys.pop_front(); + + std::vector newPolys(2); + + if (std::find(next.begin(), next.end(), start.ref) != next.end()) { + if (start.ref == end.ref) { + for (itr5 = next.begin(); itr5 != next.end(); ++itr5) { + newPolys[0].push_back(*itr5); + +#ifdef DEBUG + std::cout << "adding " << *itr5 << " to 0" << std::endl; +#endif + + if (*itr5 == start.ref) { + for (itr3 = strip.begin(); itr3 != strip.end(); ++itr3) { + newPolys[0].push_back(itr3->desc[0]); + +#ifdef DEBUG + std::cout << "adding " << itr3->desc[0] << " to 0" << std::endl; +#endif + + } + } + } + + // strip selbst ist ein polygon + + for (itr6 = strip.rbegin(); itr6 != strip.rend(); ++itr6) { + newPolys[1].push_back(itr6->desc[1]); + +#ifdef DEBUG + std::cout << "adding " << itr6->desc[1] << " to 1" << std::endl; +#endif + + } + + } else { + std::size_t curr = 0; + + for (itr5 = next.begin(); itr5 != next.end(); ++itr5) { + IdsType &poly = newPolys[curr]; + + poly.push_back(*itr5); + +#ifdef DEBUG + std::cout << "adding " << *itr5 << " to " << curr << std::endl; +#endif + + if (*itr5 == start.ref) { + for (itr3 = strip.begin(); itr3 != strip.end(); ++itr3) { + poly.push_back(itr3->desc[0]); + +#ifdef DEBUG + std::cout << "adding " << itr3->desc[0] << " to " << curr << std::endl; +#endif + + } + + curr = curr == 0 ? 1 : 0; + + } else if (*itr5 == end.ref) { + for (itr6 = strip.rbegin(); itr6 != strip.rend(); ++itr6) { + poly.push_back(itr6->desc[1]); + +#ifdef DEBUG + std::cout << "adding " << itr6->desc[1] << " to " << curr << std::endl; +#endif + + } + + curr = curr == 0 ? 1 : 0; + } + + } + } + } + + if (!newPolys[1].empty()) { + // refs aktualisieren + + for (itr4 = edges.begin(); itr4 != edges.end(); ++itr4) { + RefsType &edge = itr4->second; + + RefsType::iterator itrA; + + for (itrA = edge.begin()+1; itrA != edge.end(); ++itrA) { + StripPtR &sp = *itrA; + + if (sp.strip > start.strip) { +#ifdef DEBUG + std::cout << "sp: ind " << sp.ind << ", strip " << sp.strip << std::endl; +#endif + + RefsType::const_reverse_iterator itrB(itrA); + + std::shared_ptr _p; + + for (; itrB != edge.rend(); ++itrB) { + const StripPtR &p = *itrB; + + if (p.strip != sp.strip) { + if (p.strip <= start.strip) { +#ifdef DEBUG + std::cout << "ref " << sp.ref; +#endif + + if (p.side == Side::END) { + sp.ref = p.desc[0]; + } else { + sp.ref = p.desc[1]; + } + +#ifdef DEBUG + std::cout << " -> " << sp.ref << " (from strip " << p.strip << ", ind " << p.ind << ")" << std::endl; +#endif + + _p = std::make_shared(p); + + break; + + } + } else { +#ifdef DEBUG + std::cout << "~1 ref " << sp.ref << " -> " << p.ref << " (from strip " << p.strip << ", ind " << p.ind << ")" << std::endl; +#endif + + sp.ref = p.ref; + break; + } + } + + RefsType::const_iterator itrC(itrA); + + ++itrC; + + for (; itrC != edge.end(); ++itrC) { + const StripPtR &p = *itrC; + + if (p.ind != sp.ind) { + break; + } + + if (p.strip <= start.strip) { + if (_p && p.ind == _p->ind && p.strip < _p->strip) { + break; + } + +#ifdef DEBUG + std::cout << "~2 ref " << sp.ref; +#endif + + if (p.side == Side::START) { + sp.ref = p.desc[0]; + } else { + sp.ref = p.desc[1]; + } + +#ifdef DEBUG + std::cout << " -> " << sp.ref << " (from strip " << p.strip << ", ind " << p.ind << ")" << std::endl; +#endif + + break; + } + + } + } + } + + // erstellt die history + + auto _s = std::find_if(edge.begin(), edge.end(), [&start](const StripPtR &r) { + return &start == &r; + }); + + if (_s != edge.end()) { + for (itr7 = edge.begin(); itr7 != edge.end(); ++itr7) { + StripPtR &sp = *itr7; + StripPt &_sp = pts[sp.ind]; + auto &history = _sp.history; + + if (&sp != &start) { + if (_sp.t < pts[start.ind].t) { + history.push_back({history.back().f, start.desc[0]}); + } else { + history.push_back({start.desc[1], history.back().g}); + } + + } + + } + } + + auto _e = std::find_if(edge.begin(), edge.end(), [&end](const StripPtR &r) { + return &end == &r; + }); + + if (_e != edge.end()) { + for (itr7 = edge.begin(); itr7 != edge.end(); ++itr7) { + StripPtR &sp = *itr7; + StripPt &_sp = pts[sp.ind]; + auto &history = _sp.history; + + if (&sp != &end) { + if (_sp.t < pts[end.ind].t) { + history.push_back({history.back().f, end.desc[1]}); + } else { + history.push_back({end.desc[0], history.back().g}); + } + + } + + } + } + + // sonderfall + if (edge.size() > 1) { + StripPtR &a = edge.front(), + &b = *(edge.begin()+1); + + if (a.ind == b.ind + && b.strip == start.strip + && pts[a.ind].capt == Capt::A) { // sollte weg + +#ifdef DEBUG + std::cout << "~3 ref " << a.ref; +#endif + + if (b.side == Side::START) { + a.ref = b.desc[0]; + } else { + a.ref = b.desc[1]; + } + +#ifdef DEBUG + std::cout << " -> " << a.ref << " (from strip " << b.strip << ", ind " << b.ind << ")" << std::endl; +#endif + + } + } + + } + + // doppelte punkte entfernen + + double pt[3]; + + for (auto &newPoly : newPolys) { + IdsType _newPoly; + + std::map _pts; + + for (vtkIdType id : newPoly) { + pd->GetPoint(id, pt); + + _pts.emplace(std::piecewise_construct, + std::forward_as_tuple(id), + std::forward_as_tuple(pt[0], pt[1], pt[2])); + } + + IdsType::const_iterator itrA, itrB; + + for (itrA = newPoly.begin(); itrA != newPoly.end(); ++itrA) { + itrB = itrA+1; + if (itrB == newPoly.end()) { + itrB = newPoly.begin(); + } + + auto _a = _pts.find(*itrA); + auto _b = _pts.find(*itrB); + + if (_a->second == _b->second) { + // doppelt +#ifdef DEBUG + std::cout << "removing " << *itrA << std::endl; +#endif + } else { + _newPoly.push_back(*itrA); + } + } + + newPoly.swap(_newPoly); + + } + + // prüft, ob die erstellten polys gültig sind + + if (newPolys[0].size() > 2) { + polys.push_back(newPolys[0]); + } + + if (HasArea(strip) && newPolys[1].size() > 2) { + polys.push_back(newPolys[1]); + } + + cycle = 0; + + break; + + } else { + polys.push_back(next); + + cycle++; + } + + } + + } + + // erzeugte polys hinzufügen + + IdsType descIds; + descIds.reserve(polys.size()); + + for (auto &p : polys) { + vtkIdList *cell = vtkIdList::New(); + + for (vtkIdType &id : p) { + cell->InsertNextId(id); + } + + descIds.emplace_back(pd->InsertNextCell(VTK_POLYGON, cell)); + origCellIds->InsertNextValue(origId); + + cell->Delete(); + } + + pd->DeleteCell(polyInd); + + // holes einbauen + if (!holes.empty()) { + try { + Merger(pd, pStrips, holes, descIds, origId).Run(); + } catch (const std::runtime_error &e) { + return true; + } + } + + } + + pd->RemoveDeletedCells(); + pd->BuildCells(); + + return false; + +} + +void vtkPolyDataBooleanFilter::RestoreOrigPoints (vtkPolyData *pd, PolyStripsType &polyStrips) { + +#ifdef DEBUG + std::cout << "RestoreOrigPoints()" << std::endl; +#endif + + pd->BuildLinks(); + + vtkKdTreePointLocator *loc = vtkKdTreePointLocator::New(); + loc->SetDataSet(pd); + loc->BuildLocator(); + + struct Cmp { + bool operator() (const StripPt &l, const StripPt &r) const { + return l.ind < r.ind; + } + }; + + PolyStripsType::const_iterator itr; + StripPtsType::const_iterator itr2; + + std::set, Cmp> ends; + + for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) { + const PStrips &pStrips = itr->second; + + for (itr2 = pStrips.pts.begin(); itr2 != pStrips.pts.end(); ++itr2) { + const StripPt &sp = itr2->second; + + if (sp.capt != Capt::NOT) { + ends.emplace(sp); + } + + } + } + + vtkIdList *pts = vtkIdList::New(); + + vtkIdType i, numPts; + + for (const StripPt &sp : ends) { + FindPoints(loc, sp.cutPt, pts); + numPts = pts->GetNumberOfIds(); + + for (i = 0; i < numPts; i++) { + pd->GetPoints()->SetPoint(pts->GetId(i), sp.pt); + } + } + + pts->Delete(); + + loc->FreeSearchStructure(); + loc->Delete(); + +} + +void vtkPolyDataBooleanFilter::DisjoinPolys (vtkPolyData *pd, PolyStripsType &polyStrips) { + +#ifdef DEBUG + std::cout << "DisjoinPolys()" << std::endl; +#endif + + pd->BuildLinks(); + + vtkKdTreePointLocator *loc = vtkKdTreePointLocator::New(); + loc->SetDataSet(pd); + loc->BuildLocator(); + + struct Cmp { + bool operator() (const StripPt &l, const StripPt &r) const { + return l.ind < r.ind; + } + }; + + PolyStripsType::const_iterator itr; + StripPtsType::const_iterator itr2; + + std::set, Cmp> ends; + + for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) { + const PStrips &pStrips = itr->second; + + for (itr2 = pStrips.pts.begin(); itr2 != pStrips.pts.end(); ++itr2) { + const StripPt &sp = itr2->second; + + if (sp.capt == Capt::A) { + ends.emplace(sp); + } + } + } + + vtkIdList *pts = vtkIdList::New(); + vtkIdList *cells = vtkIdList::New(); + + vtkIdType i, j, numPts, numCells; + + for (const StripPt &sp : ends) { + FindPoints(loc, sp.pt, pts); + numPts = pts->GetNumberOfIds(); + + for (i = 0; i < numPts; i++) { + pd->GetPointCells(pts->GetId(i), cells); + numCells = cells->GetNumberOfIds(); + + if (numCells > 1) { + for (j = 0; j < numCells; j++) { + pd->ReplaceCellPoint(cells->GetId(j), pts->GetId(i), pd->GetPoints()->InsertNextPoint(sp.pt)); + } + } + } + } + + cells->Delete(); + pts->Delete(); + + loc->FreeSearchStructure(); + loc->Delete(); + +} + +void vtkPolyDataBooleanFilter::ResolveOverlaps (vtkPolyData *pd, vtkIdTypeArray *conts, PolyStripsType &polyStrips) { + +#ifdef DEBUG + std::cout << "ResolveOverlaps()" << std::endl; +#endif + + pd->BuildCells(); + pd->BuildLinks(); + + contLines->BuildLinks(); + + vtkKdTreePointLocator *loc = vtkKdTreePointLocator::New(); + loc->SetDataSet(pd); + loc->BuildLocator(); + + struct Cmp { + bool operator() (const StripPt &l, const StripPt &r) const { + return l.ind < r.ind; + } + }; + + PolyStripsType::const_iterator itr; + StripPtsType::const_iterator itr2; + + std::set, Cmp> ends; + + for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) { + const PStrips &pStrips = itr->second; + + for (itr2 = pStrips.pts.begin(); itr2 != pStrips.pts.end(); ++itr2) { + const StripPt &sp = itr2->second; + + if (sp.capt == Capt::EDGE) { + ends.emplace(sp); + } + } + } + + vtkIdList *links = vtkIdList::New(), + *ptsA = vtkIdList::New(), + *ptsB = vtkIdList::New(), + *cells = vtkIdList::New(), + *poly = vtkIdList::New(); + + double ptA[3], ptB[3]; + + vtkIdType i, j, k, numPtsA, numPtsB, numCells, numPts, idA, idB; + + std::map sharedPts; + + for (const StripPt &sp : ends) { + contLines->GetPointCells(sp.ind, links); + + if (links->GetNumberOfIds() == 2 + && conts->GetValue(links->GetId(0)) != conts->GetValue(links->GetId(1))) { + + std::vector history(sp.history.rbegin(), sp.history.rend()); + + for (const Pair &p : history) { + pd->GetPoint(p.f, ptA); + pd->GetPoint(p.g, ptB); + + FindPoints(loc, ptA, ptsA); + FindPoints(loc, ptB, ptsB); + + numPtsA = ptsA->GetNumberOfIds(); + numPtsB = ptsB->GetNumberOfIds(); + + std::vector cellsA, cellsB; + + for (i = 0; i < numPtsA; i++) { + pd->GetPointCells(ptsA->GetId(i), cells); + numCells = cells->GetNumberOfIds(); + + for (j = 0; j < numCells; j++) { + cellsA.emplace_back(ptsA->GetId(i), cells->GetId(j)); + } + } + + for (i = 0; i < numPtsB; i++) { + pd->GetPointCells(ptsB->GetId(i), cells); + numCells = cells->GetNumberOfIds(); + + for (j = 0; j < numCells; j++) { + cellsB.emplace_back(ptsB->GetId(i), cells->GetId(j)); + } + } + + for (auto &cellA : cellsA) { + for (auto &cellB : cellsB) { + if (cellA.g == cellB.g) { + // kante/poly existiert noch + + pd->GetCellPoints(cellA.g, poly); + numPts = poly->GetNumberOfIds(); + + for (k = 0; k < numPts; k++) { + idA = poly->GetId(k); + idB = k+1 == numPts ? poly->GetId(0) : poly->GetId(k+1); + + if (cellB.f == idA + && cellA.f == idB) { + + IdsType &pts = sharedPts[Pair{sp.ind, cellA.g}]; + + pts.push_back(cellA.f); + pts.push_back(cellB.f); + + break; + } + } + + } + } + } + + + } + + } + } + + double pt[3]; + + decltype(sharedPts)::iterator itr3; + + for (itr3 = sharedPts.begin(); itr3 != sharedPts.end(); ++itr3) { + const Pair &p = itr3->first; + + IdsType &pts = itr3->second; + + std::sort(pts.begin(), pts.end()); + + auto found = std::adjacent_find(pts.begin(), pts.end()); + + if (found != pts.end()) { + contLines->GetPoint(p.f, pt); + + i = pd->GetPoints()->InsertNextPoint(pt); + pd->ReplaceCellPoint(p.g, *found, i); + } + } + + links->Delete(); + ptsA->Delete(); + ptsB->Delete(); + cells->Delete(); + poly->Delete(); + + loc->FreeSearchStructure(); + loc->Delete(); + +} + +void vtkPolyDataBooleanFilter::AddAdjacentPoints (vtkPolyData *pd, vtkIdTypeArray *conts, PolyStripsType &polyStrips) { + +#ifdef DEBUG + std::cout << "AddAdjacentPoints()" << std::endl; +#endif + + pd->BuildLinks(); + + vtkIdTypeArray *origCellIds = vtkIdTypeArray::SafeDownCast(pd->GetCellData()->GetScalars("OrigCellIds")); + + struct Cmp { + bool operator() (const StripPt &l, const StripPt &r) const { + return l.t < r.t; + } + }; + + vtkKdTreePointLocator *loc = vtkKdTreePointLocator::New(); + loc->SetDataSet(pd); + loc->BuildLocator(); + + vtkIdList *cells = vtkIdList::New(), + *ptsA = vtkIdList::New(), + *ptsB = vtkIdList::New(), + *poly = vtkIdList::New(), + *newPoly = vtkIdList::New(); + + vtkIdType i, j, numCells, numPtsA, numPtsB, numPts, idA, idB; + + PolyStripsType::const_iterator itr; + StripPtsType::const_iterator itr2; + + for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) { + const PStrips &pStrips = itr->second; + + std::map, Cmp>> edgePts; + + for (itr2 = pStrips.pts.begin(); itr2 != pStrips.pts.end(); ++itr2) { + const StripPt &sp = itr2->second; + + if (sp.capt == Capt::EDGE) { + edgePts[{sp.edge[0], sp.edge[1]}].emplace(sp); + } + } + + decltype(edgePts)::iterator itr3; + + for (itr3 = edgePts.begin(); itr3 != edgePts.end(); ++itr3) { + const Pair &edge = itr3->first; + auto pts = itr3->second; + + StripPt spA, spB; + + pd->GetPoint(edge.f, spA.pt); + pd->GetPoint(edge.g, spB.pt); + + spA.t = 0; + spB.t = 1; + + pts.emplace(spA); + pts.emplace(spB); + + std::vector _pts(pts.rbegin(), pts.rend()); + + decltype(_pts)::const_iterator itrA, itrB, itrC; + + itrA = _pts.begin(); + + while (itrA != _pts.end()-1) { + itrB = itrA+1; + + while (itrB != _pts.end()-1) { + contLines->GetPointCells(itrB->get().ind, cells); + numCells = cells->GetNumberOfIds(); + + // std::cout << "["; + // for (i = 0; i < numCells; i++) { + // std::cout << conts->GetValue(cells->GetId(i)) << ", "; + // } + // std::cout << "]" << std::endl; + + // break; + + ++itrB; + } + + if (std::distance(itrA, itrB) > 1) { + FindPoints(loc, itrA->get().pt, ptsA); + FindPoints(loc, itrB->get().pt, ptsB); + + numPtsA = ptsA->GetNumberOfIds(); + numPtsB = ptsB->GetNumberOfIds(); + + std::vector polysA, polysB; + + for (i = 0; i < numPtsA; i++) { + pd->GetPointCells(ptsA->GetId(i), cells); + numCells = cells->GetNumberOfIds(); + + for (j = 0; j < numCells; j++) { + polysA.emplace_back(cells->GetId(j), ptsA->GetId(i)); + } + } + + for (i = 0; i < numPtsB; i++) { + pd->GetPointCells(ptsB->GetId(i), cells); + numCells = cells->GetNumberOfIds(); + + for (j = 0; j < numCells; j++) { + polysB.emplace_back(cells->GetId(j), ptsB->GetId(i)); + } + } + + /*for (const Pair &pA : polysA) { + std::cout << "pA " << pA << std::endl; + } + + for (const Pair &pB : polysB) { + std::cout << "pB " << pB << std::endl; + }*/ + + for (const Pair &pA : polysA) { + for (const Pair &pB : polysB) { + if (pA.f == pB.f && pd->GetCellType(pA.f) != VTK_EMPTY_CELL) { + + pd->GetCellPoints(pA.f, poly); + numPts = poly->GetNumberOfIds(); + + newPoly->Reset(); + + for (j = 0; j < numPts; j++) { + newPoly->InsertNextId(poly->GetId(j)); + + idA = poly->GetId(j); + idB = j+1 == numPts ? poly->GetId(0) : poly->GetId(j+1); + + if (pA.g == idA + && pB.g == idB) { + + for (itrC = itrA+1; itrC != itrB; ++itrC) { + // newPoly->InsertNextId(pd->InsertNextLinkedPoint(itrC->get().pt, 1)); + + pd->InsertNextLinkedPoint(1); + + newPoly->InsertNextId(pd->GetPoints()->InsertNextPoint(itrC->get().pt)); + } + + } + + pd->RemoveReferenceToCell(idA, pA.f); + } + + pd->DeleteCell(pA.f); + + pd->InsertNextLinkedCell(VTK_POLYGON, newPoly->GetNumberOfIds(), newPoly->GetPointer(0)); + + origCellIds->InsertNextValue(origCellIds->GetValue(pA.f)); + + break; + } + } + } + } + + itrA = itrB; + } + } + } + + cells->Delete(); + ptsA->Delete(); + ptsB->Delete(); + poly->Delete(); + newPoly->Delete(); + + loc->FreeSearchStructure(); + loc->Delete(); + + pd->RemoveDeletedCells(); + +} + +void vtkPolyDataBooleanFilter::MergePoints (vtkPolyData *pd, PolyStripsType &polyStrips) { + +#ifdef DEBUG + std::cout << "MergePoints()" << std::endl; +#endif + + pd->BuildCells(); + pd->BuildLinks(); + + contLines->BuildLinks(); + + vtkKdTreePointLocator *loc = vtkKdTreePointLocator::New(); + loc->SetDataSet(pd); + loc->BuildLocator(); + + PolyStripsType::const_iterator itr; + StripsType::const_iterator itr2; + + std::map> neighPts; + + vtkIdList *pts = vtkIdList::New(); + + vtkIdType i, j, numPts; + + for (itr = polyStrips.begin(); itr != polyStrips.end(); ++itr) { + const PStrips &pStrips = itr->second; + + for (itr2 = pStrips.strips.begin(); itr2 != pStrips.strips.end(); ++itr2) { + const StripType &strip = *itr2; + + const StripPtR &spA = strip.front(), + &spB = strip.back(); + + const auto beforeA = pStrips.pts.find((strip.begin()+1)->ind), + beforeB = pStrips.pts.find((strip.end()-2)->ind); + + FindPoints(loc, beforeA->second.pt, pts); + numPts = pts->GetNumberOfIds(); + + for (i = 0; i < numPts; i++) { + neighPts[spA.ind].insert(pts->GetId(i)); + } + + FindPoints(loc, beforeB->second.pt, pts); + numPts = pts->GetNumberOfIds(); + + for (i = 0; i < numPts; i++) { + neighPts[spB.ind].insert(pts->GetId(i)); + } + } + } + + double pt[3]; + + vtkIdList *polys = vtkIdList::New(), + *poly = vtkIdList::New(); + + vtkIdType ind, polyId, _numPts, before, after; + + decltype(neighPts)::const_iterator itr3; + + for (itr3 = neighPts.begin(); itr3 != neighPts.end(); ++itr3) { + const auto &inds = itr3->second; + + std::map> pairs; + + contLines->GetPoint(itr3->first, pt); + + FindPoints(loc, pt, pts); + numPts = pts->GetNumberOfIds(); + + for (i = 0; i < numPts; i++) { + ind = pts->GetId(i); + + pd->GetPointCells(ind, polys); + + if (polys->GetNumberOfIds() > 0) { + polyId = polys->GetId(0); + + pd->GetCellPoints(polyId, poly); + _numPts = poly->GetNumberOfIds(); + + for (j = 0; j < _numPts; j++) { + if (poly->GetId(j) == ind) { + break; + } + } + + // wieder davor und danach ermitteln + + before = poly->GetId(j == 0 ? _numPts-1 : j-1); + after = poly->GetId(j+1 == _numPts ? 0 : j+1); + + if (std::find(inds.begin(), inds.end(), before) == inds.end()) { + pd->GetPoint(before, pt); + pairs[{pt[0], pt[1], pt[2]}].emplace_back(polyId, ind); + } + + if (std::find(inds.begin(), inds.end(), after) == inds.end()) { + pd->GetPoint(after, pt); + pairs[{pt[0], pt[1], pt[2]}].emplace_back(polyId, ind); + } + + } + } + + std::deque>> Pairs; + + decltype(pairs)::const_iterator itr4; + + for (itr4 = pairs.begin(); itr4 != pairs.end(); ++itr4) { + const auto &p = itr4->second; + + if (p.size() == 2) { + auto _pts = {std::ref(p.front()), std::ref(p.back())}; // std::initializer_list + Pairs.emplace_back(_pts); + } + } + + decltype(Pairs)::iterator itr5; + + /*for (itr5 = Pairs.begin(); itr5 != Pairs.end(); ++itr5) { + for (auto &p : *itr5) { + std::cout << p.get() << ", "; + } + std::cout << std::endl; + }*/ + + decltype(Pairs)::value_type group; + + decltype(group)::const_iterator itr6; + + while (!Pairs.empty()) { + if (group.empty()) { + group = Pairs.front(); + Pairs.pop_front(); + } + + itr5 = Pairs.begin(); + + while (itr5 != Pairs.end()) { + const auto &next = *itr5; + + if (next.front().get() == group.front().get()) { + group.emplace_front(next.back()); + Pairs.erase(itr5); + itr5 = Pairs.begin(); + } else if (next.front().get() == group.back().get()) { + group.emplace_back(next.back()); + Pairs.erase(itr5); + itr5 = Pairs.begin(); + } else if (next.back().get() == group.front().get()) { + group.emplace_front(next.front()); + Pairs.erase(itr5); + itr5 = Pairs.begin(); + } else if (next.back().get() == group.back().get()) { + group.emplace_back(next.front()); + Pairs.erase(itr5); + itr5 = Pairs.begin(); + } else { + ++itr5; + } + } + + if (itr5 == Pairs.end()) { + for (itr6 = group.begin()+1; itr6 != group.end(); ++itr6) { + pd->ReplaceCellPoint(itr6->get().f, itr6->get().g, group.front().get().g); + } + + group.clear(); + } + } + + } + + polys->Delete(); + poly->Delete(); + pts->Delete(); + + loc->FreeSearchStructure(); + loc->Delete(); + +} + +enum class Congr { + EQUAL, + OPPOSITE, + NOT +}; + +class PolyAtEdge { + vtkPolyData *pd; + +public: + PolyAtEdge (vtkPolyData *_pd, vtkIdType _polyId, vtkIdType _ptIdA, vtkIdType _ptIdB) : pd(_pd), polyId(_polyId), ptIdA(_ptIdA), ptIdB(_ptIdB), loc(Loc::NONE) { + + double ptA[3], ptB[3]; + + pd->GetPoint(ptIdA, ptA); + pd->GetPoint(ptIdB, ptB); + + vtkMath::Subtract(ptB, ptA, e); + vtkMath::Normalize(e); + + const vtkIdType *poly; + + vtkIdType numPts; + pd->GetCellPoints(polyId, numPts, poly); + + ComputeNormal(pd->GetPoints(), n, numPts, poly); + + vtkMath::Cross(e, n, r); + + } + + vtkIdType polyId, ptIdA, ptIdB; + double n[3], e[3], r[3]; + + Loc loc; + + friend std::ostream& operator<< (std::ostream &out, const PolyAtEdge &p) { + out << "polyId " << p.polyId << ", ptIdA " << p.ptIdA << ", ptIdB " << p.ptIdB; + return out; + } + + static constexpr double eps {.99999999}; // ~.0081deg + + Congr IsCongruent (const PolyAtEdge &p) const { + double cong = vtkMath::Dot(n, p.n); + + if (cong > eps || cong < -eps) { + double ang = vtkMath::Dot(r, p.r); + + if (ang > eps) { + if (cong > eps) { + // normalen sind gleich ausgerichtet + return Congr::EQUAL; + } else { + return Congr::OPPOSITE; + } + } + } + + return Congr::NOT; + } + +}; + + +class PolyPair { +public: + PolyPair (PolyAtEdge _pA, PolyAtEdge _pB) : pA(_pA), pB(_pB) {} + + PolyAtEdge pA, pB; + + void GetLoc (PolyAtEdge &pT, vtkIdType mode) { + + Congr cA = pA.IsCongruent(pT), + cB = pB.IsCongruent(pT); + +#ifdef DEBUG + std::cout << "GetLoc() -> polyId " << pT.polyId + << ", cA " << cA + << ", cB " << cB + << std::endl; + + if (cA != Congr::NOT || cB != Congr::NOT) { + assert(cA != cB); + } + +#endif + + if (cA == Congr::EQUAL || cA == Congr::OPPOSITE) { + if (cA == Congr::OPPOSITE) { + // normalen sind entgegengesetzt gerichtet + + if (mode == OPER_INTERSECTION) { + pA.loc = Loc::OUTSIDE; + pT.loc = Loc::OUTSIDE; + } else { + pA.loc = Loc::INSIDE; + pT.loc = Loc::INSIDE; + } + } else if (mode == OPER_UNION || mode == OPER_INTERSECTION) { + pA.loc = Loc::INSIDE; + pT.loc = Loc::OUTSIDE; + } + + } else if (cB == Congr::EQUAL || cB == Congr::OPPOSITE) { + if (cB == Congr::OPPOSITE) { + // normalen sind entgegengesetzt gerichtet + + if (mode == OPER_INTERSECTION) { + pB.loc = Loc::OUTSIDE; + pT.loc = Loc::OUTSIDE; + } else { + pB.loc = Loc::INSIDE; + pT.loc = Loc::INSIDE; + } + } else if (mode == OPER_UNION || mode == OPER_INTERSECTION) { + pB.loc = Loc::INSIDE; + pT.loc = Loc::OUTSIDE; + } + + } else { + double alpha = GetAngle(pA.r, pB.r, pA.e), + beta = GetAngle(pA.r, pT.r, pA.e); + + if (beta > alpha) { + pT.loc = Loc::INSIDE; + } else { + pT.loc = Loc::OUTSIDE; + } + } + } + +}; + + +std::shared_ptr GetEdgePolys (vtkPolyData *pd, vtkIdList *ptsA, vtkIdList *ptsB) { + +#ifdef DEBUG + std::cout << "GetEdgePolys()" << std::endl; +#endif + + std::vector p; + + vtkIdType numPtsA = ptsA->GetNumberOfIds(), + numPtsB = ptsB->GetNumberOfIds(); + + vtkIdList *polys = vtkIdList::New(); + + vtkIdType i, j, numCells; + + for (i = 0; i < numPtsA; i++) { + pd->GetPointCells(ptsA->GetId(i), polys); + numCells = polys->GetNumberOfIds(); + + for (j = 0; j < numCells; j++) { + p.emplace_back(ptsA->GetId(i), polys->GetId(j)); + } + } + + for (i = 0; i < numPtsB; i++) { + pd->GetPointCells(ptsB->GetId(i), polys); + numCells = polys->GetNumberOfIds(); + + for (j = 0; j < numCells; j++) { + p.emplace_back(ptsB->GetId(i), polys->GetId(j)); + } + } + + polys->Delete(); + + std::map pEdges; + + std::vector::const_iterator itr; + for (itr = p.begin(); itr != p.end(); ++itr) { + pEdges[itr->g].push_back(itr->f); + } + + std::vector opp; + + vtkIdList *poly = vtkIdList::New(); + + vtkIdType numPts, a, b; + + std::map::const_iterator itr2; + + for (itr2 = pEdges.begin(); itr2 != pEdges.end(); ++itr2) { + const IdsType &pts = itr2->second; + + if (pts.size() > 1) { + pd->GetCellPoints(itr2->first, poly); + numPts = poly->GetNumberOfIds(); + + for (i = 0; i < numPts; i++) { + a = poly->GetId(i); + b = i+1 == numPts ? poly->GetId(0) : poly->GetId(i+1); + + if (std::find(pts.begin(), pts.end(), a) != pts.end() + && std::find(pts.begin(), pts.end(), b) != pts.end()) { + + opp.emplace_back(pd, itr2->first, a, b); + } + } + + } + } + + poly->Delete(); + +#ifdef DEBUG + for (auto &op : opp) { + std::cout << op << std::endl; + } +#endif + + if (opp.size() != 2) { + return nullptr; + } + + return std::make_shared(opp[0], opp[1]); + +} + +bool vtkPolyDataBooleanFilter::CombineRegions () { + +#ifdef DEBUG + std::cout << "CombineRegions()" << std::endl; +#endif + + vtkSmartPointer filterdA = vtkSmartPointer::New(); + filterdA->DeepCopy(modPdA); + + vtkSmartPointer filterdB = vtkSmartPointer::New(); + filterdB->DeepCopy(modPdB); + + // ungenutzte punkte löschen + vtkSmartPointer cleanA = vtkSmartPointer::New(); + cleanA->PointMergingOff(); + cleanA->SetInputData(filterdA); + + vtkSmartPointer cleanB = vtkSmartPointer::New(); + cleanB->PointMergingOff(); + cleanB->SetInputData(filterdB); + + // regionen mit skalaren ausstatten + vtkSmartPointer cfA = vtkSmartPointer::New(); + cfA->SetExtractionModeToAllRegions(); + cfA->ColorRegionsOn(); + cfA->SetInputConnection(cleanA->GetOutputPort()); + + vtkSmartPointer cfB = vtkSmartPointer::New(); + cfB->SetExtractionModeToAllRegions(); + cfB->ColorRegionsOn(); + cfB->SetInputConnection(cleanB->GetOutputPort()); + + cfA->Update(); + cfB->Update(); + + vtkPolyData *pdA = cfA->GetOutput(); + vtkPolyData *pdB = cfB->GetOutput(); + +#ifdef DEBUG + std::cout << "Exporting modPdA_8.vtk" << std::endl; + WriteVTK("modPdA_8.vtk", pdA); + + std::cout << "Exporting modPdB_8.vtk" << std::endl; + WriteVTK("modPdB_8.vtk", pdB); +#endif + + resultC->ShallowCopy(contLines); + + if (OperMode == OPER_NONE) { + resultA->ShallowCopy(pdA); + resultB->ShallowCopy(pdB); + + return false; + } + + vtkSmartPointer plA = vtkSmartPointer::New(); + plA->SetDataSet(pdA); + plA->BuildLocator(); + + vtkSmartPointer plB = vtkSmartPointer::New(); + plB->SetDataSet(pdB); + plB->BuildLocator(); + + pdA->BuildLinks(); + pdB->BuildLinks(); + + vtkIdTypeArray *scalarsA = vtkIdTypeArray::SafeDownCast(pdA->GetPointData()->GetScalars()); + vtkIdTypeArray *scalarsB = vtkIdTypeArray::SafeDownCast(pdB->GetPointData()->GetScalars()); + + vtkSmartPointer line = vtkSmartPointer::New(); + + double ptA[3], ptB[3]; + + vtkSmartPointer fptsA = vtkSmartPointer::New(); + vtkSmartPointer lptsA = vtkSmartPointer::New(); + + vtkSmartPointer fptsB = vtkSmartPointer::New(); + vtkSmartPointer lptsB = vtkSmartPointer::New(); + + std::map locsA, locsB; + + vtkIdType i, j, numLines = contLines->GetNumberOfCells(); + + IdsType _failed; + + for (i = 0; i < numLines; i++) { + + if (contLines->GetCellType(i) == VTK_EMPTY_CELL) { + continue; + } + + contLines->GetCellPoints(i, line); + + contLines->GetPoint(line->GetId(0), ptA); + contLines->GetPoint(line->GetId(1), ptB); + + FindPoints(plA, ptA, fptsA); + FindPoints(plB, ptA, fptsB); + +#ifdef DEBUG + std::cout << "line " << i << std::endl; +#else + + // bereits behandelte regionen werden nicht noch einmal untersucht + + vtkIdType notLocated = 0; + + for (j = 0; j < fptsA->GetNumberOfIds(); j++) { + if (locsA.count(scalarsA->GetValue(fptsA->GetId(j))) == 0) { + notLocated++; + } + } + + for (j = 0; j < fptsB->GetNumberOfIds(); j++) { + if (locsB.count(scalarsB->GetValue(fptsB->GetId(j))) == 0) { + notLocated++; + } + } + + if (notLocated == 0) { + continue; + } + +#endif + + FindPoints(plA, ptB, lptsA); + FindPoints(plB, ptB, lptsB); + + auto ppA = GetEdgePolys(pdA, fptsA, lptsA); + auto ppB = GetEdgePolys(pdB, fptsB, lptsB); + + if (ppA && ppB) { + + ppB->GetLoc(ppA->pA, OperMode); + ppB->GetLoc(ppA->pB, OperMode); + + ppA->GetLoc(ppB->pA, OperMode); + ppA->GetLoc(ppB->pB, OperMode); + + vtkIdType fsA = scalarsA->GetValue(ppA->pA.ptIdA); + vtkIdType lsA = scalarsA->GetValue(ppA->pB.ptIdA); + + vtkIdType fsB = scalarsB->GetValue(ppB->pA.ptIdA); + vtkIdType lsB = scalarsB->GetValue(ppB->pB.ptIdA); + +#ifdef DEBUG + std::cout << "polyId " << ppA->pA.polyId << ", sA " << fsA << ", loc " << ppA->pA.loc << std::endl; + std::cout << "polyId " << ppA->pB.polyId << ", sA " << lsA << ", loc " << ppA->pB.loc << std::endl; + std::cout << "polyId " << ppB->pA.polyId << ", sB " << fsB << ", loc " << ppB->pA.loc << std::endl; + std::cout << "polyId " << ppB->pB.polyId << ", sB " << lsB << ", loc " << ppB->pB.loc << std::endl; + + if (locsA.count(fsA) > 0 && locsA[fsA] != ppA->pA.loc) { + std::cout << "sA " << fsA << ": " << locsA[fsA] << " -> " << ppA->pA.loc << std::endl; + } + + if (locsA.count(lsA) > 0 && locsA[lsA] != ppA->pB.loc) { + std::cout << "sA " << lsA << ": " << locsA[lsA] << " -> " << ppA->pB.loc << std::endl; + } + + if (locsB.count(fsB) > 0 && locsB[fsB] != ppB->pA.loc) { + std::cout << "sB " << fsB << ": " << locsB[fsB] << " -> " << ppB->pA.loc << std::endl; + } + + if (locsB.count(lsB) > 0 && locsB[lsB] != ppB->pB.loc) { + std::cout << "sB " << lsB << ": " << locsB[lsB] << " -> " << ppB->pB.loc << std::endl; + } + +#endif + + locsA.emplace(fsA, ppA->pA.loc); + locsA.emplace(lsA, ppA->pB.loc); + + locsB.emplace(fsB, ppB->pA.loc); + locsB.emplace(lsB, ppB->pB.loc); + + } else { + _failed.push_back(i); + + // return true; + } + + } + + if (_failed.size() > 0) { + for (auto i : _failed) { + std::cout << "failed at " << i + << std::endl; + } + + return true; + } + + // reale kombination der ermittelten regionen + + Loc comb[] = {Loc::OUTSIDE, Loc::OUTSIDE}; + + if (OperMode == OPER_INTERSECTION) { + comb[0] = Loc::INSIDE; + comb[1] = Loc::INSIDE; + } else if (OperMode == OPER_DIFFERENCE) { + comb[1] = Loc::INSIDE; + } else if (OperMode == OPER_DIFFERENCE2) { + comb[0] = Loc::INSIDE; + } + + vtkIdType numA = cfA->GetNumberOfExtractedRegions(), + numB = cfB->GetNumberOfExtractedRegions(); + + cfA->SetExtractionModeToSpecifiedRegions(); + cfB->SetExtractionModeToSpecifiedRegions(); + + std::map::const_iterator itr; + + for (itr = locsA.begin(); itr != locsA.end(); itr++) { + if (itr->second == comb[0]) { + cfA->AddSpecifiedRegion(itr->first); + } + } + + for (itr = locsB.begin(); itr != locsB.end(); itr++) { + if (itr->second == comb[1]) { + cfB->AddSpecifiedRegion(itr->first); + } + } + + // nicht beteiligte regionen hinzufügen + + if (OperMode == OPER_UNION || OperMode == OPER_DIFFERENCE) { + for (i = 0; i < numA; i++) { + if (locsA.count(i) == 0) { + cfA->AddSpecifiedRegion(i); + } + } + } + + if (OperMode == OPER_UNION || OperMode == OPER_DIFFERENCE2) { + for (i = 0; i < numB; i++) { + if (locsB.count(i) == 0) { + cfB->AddSpecifiedRegion(i); + } + } + } + + // nach innen zeigende normalen umkehren + + cfA->Update(); + cfB->Update(); + + vtkPolyData *regsA = cfA->GetOutput(); + vtkPolyData *regsB = cfB->GetOutput(); + + scalarsA = vtkIdTypeArray::SafeDownCast(regsA->GetPointData()->GetScalars()); + scalarsB = vtkIdTypeArray::SafeDownCast(regsB->GetPointData()->GetScalars()); + + if (OperMode != OPER_INTERSECTION) { + if (comb[0] == Loc::INSIDE) { + for (i = 0; i < regsA->GetNumberOfCells(); i++) { + if (locsA.count(scalarsA->GetValue(i)) == 1) { + regsA->ReverseCell(i); + } + } + } + + if (comb[1] == Loc::INSIDE) { + for (i = 0; i < regsB->GetNumberOfCells(); i++) { + if (locsB.count(scalarsB->GetValue(i)) == 1) { + regsB->ReverseCell(i); + } + } + } + } + + // OrigCellIds und CellData + + vtkIdTypeArray *origCellIdsA = vtkIdTypeArray::SafeDownCast(regsA->GetCellData()->GetScalars("OrigCellIds")); + vtkIdTypeArray *origCellIdsB = vtkIdTypeArray::SafeDownCast(regsB->GetCellData()->GetScalars("OrigCellIds")); + + vtkIdTypeArray *newOrigCellIdsA = vtkIdTypeArray::New(); + newOrigCellIdsA->SetName("OrigCellIdsA"); + + vtkIdTypeArray *newOrigCellIdsB = vtkIdTypeArray::New(); + newOrigCellIdsB->SetName("OrigCellIdsB"); + + vtkCellData *newCellDataA = vtkCellData::New(); + newCellDataA->CopyAllocate(cellDataA); + + vtkCellData *newCellDataB = vtkCellData::New(); + newCellDataB->CopyAllocate(cellDataB); + + vtkIdType cellId; + + for (i = 0; i < regsA->GetNumberOfCells(); i++) { + cellId = cellIdsA->GetValue(origCellIdsA->GetValue(i)); + + newOrigCellIdsA->InsertNextValue(cellId); + newOrigCellIdsB->InsertNextValue(-1); + + newCellDataA->CopyData(cellDataA, cellId, i); + } + + for (i = 0; i < regsB->GetNumberOfCells(); i++) { + cellId = cellIdsB->GetValue(origCellIdsB->GetValue(i)); + + newOrigCellIdsB->InsertNextValue(cellId); + newOrigCellIdsA->InsertNextValue(-1); + + newCellDataB->CopyData(cellDataB, cellId, i); + } + + regsA->GetCellData()->Initialize(); + regsB->GetCellData()->Initialize(); + + regsA->GetCellData()->ShallowCopy(newCellDataA); + regsB->GetCellData()->ShallowCopy(newCellDataB); + + newCellDataA->Delete(); + newCellDataB->Delete(); + + // zusammenführung + + vtkAppendPolyData *app = vtkAppendPolyData::New(); + app->AddInputData(regsA); + app->AddInputData(regsB); + + // entfernt ungenutzte punkte + vtkCleanPolyData *cleanApp = vtkCleanPolyData::New(); + cleanApp->PointMergingOff(); + cleanApp->SetInputConnection(app->GetOutputPort()); + + // färbt die regionen nochmal neu ein, damit mehrere regionen nicht die gleiche farbe haben + + vtkPolyDataConnectivityFilter *cfApp = vtkPolyDataConnectivityFilter::New(); + cfApp->SetExtractionModeToAllRegions(); + cfApp->ColorRegionsOn(); + cfApp->SetInputConnection(cleanApp->GetOutputPort()); + + cfApp->Update(); + + vtkPolyData *cfPd = cfApp->GetOutput(); + + // resultB bleibt hier leer + + resultA->ShallowCopy(cfPd); + + resultA->GetCellData()->AddArray(newOrigCellIdsA); + resultA->GetCellData()->AddArray(newOrigCellIdsB); + + // aufräumen + + cfApp->Delete(); + cleanApp->Delete(); + app->Delete(); + + newOrigCellIdsB->Delete(); + newOrigCellIdsA->Delete(); + + plB->FreeSearchStructure(); + plA->FreeSearchStructure(); + + return false; + +} + +Merger::Merger (vtkPolyData *pd, const PStrips &pStrips, const StripsType &strips, const IdsType &descIds, vtkIdType origId) : pd(pd), pStrips(pStrips), origId(origId) { + + const StripPtsType &pts = pStrips.pts; + const Base &base = pStrips.base; + + vtkIdType i, num; + const vtkIdType *cell; + + double pt[3]; + + for (auto id : descIds) { + pd->GetCellPoints(id, num, cell); + + Poly p; + + for (i = 0; i < num; i++) { + pd->GetPoint(cell[i], pt); + + double proj[2]; + Transform(pt, proj, base); + + p.emplace_back(proj[0], proj[1], 0, cell[i]); + + } + + polys.push_back(p); + + pd->DeleteCell(id); + } + + for (auto &strip : strips) { + Poly p; + + for (auto &sp : strip) { + const double *pt = pts.at(sp.ind).pt; + + double proj[2]; + Transform(pt, proj, base); + + p.emplace_back(proj[0], proj[1], 0); + } + + p.pop_back(); + + double n[3]; + ComputeNormal(p, n); + + if (n[2] < 0) { + Poly q(p.rbegin(), p.rend()); + p.swap(q); + } + + polys.push_back(p); + } + +} + +void Merger::Run () { + // mergen mit hilfe von vtkKdTree und vtkModifiedBSPTree + + vtkPoints *pdPts = pd->GetPoints(); + vtkIdTypeArray *origCellIds = vtkIdTypeArray::SafeDownCast(pd->GetCellData()->GetScalars("OrigCellIds")); + + assert(origCellIds != nullptr); + + const Base &base = pStrips.base; + + std::vector groups(polys.size()); + + PolysType::const_iterator itrA, itrB; + + for (itrA = polys.begin(); itrA != polys.end(); ++itrA) { + for (itrB = polys.begin(); itrB != polys.end(); ++itrB) { + if (itrA != itrB && PointInPoly(*itrB, *itrA->begin())) { + groups[static_cast(itrB-polys.begin())].push_back(static_cast(itrA-polys.begin())); + } + } + } + + std::size_t parent = 0; + + for (auto &group : groups) { + GroupType parents; + + for (auto &index : group) { + const GroupType &_group = groups[index]; + parents.insert(parents.end(), _group.begin(), _group.end()); + } + + std::sort(group.begin(), group.end()); + std::sort(parents.begin(), parents.end()); + + GroupType _group {parent++}; + std::set_difference(group.begin(), group.end(), parents.begin(), parents.end(), std::back_inserter(_group)); + + std::cout << "["; + for (auto &index : _group) { + std::cout << index << ", "; + } + std::cout << "]" << std::endl; + + PolysType merged; + + MergeGroup(_group, merged); + + std::map newIds; + + for (auto &poly : merged) { + vtkSmartPointer newCell = vtkSmartPointer::New(); + + for (auto &p : poly) { + double in[] = {p.x, p.y}, + out[3]; + + BackTransform(in, out, base); + + vtkIdType id = p.id; + + if (id == NOTSET) { + auto itr = newIds.find(p); + + if (itr == newIds.end()) { + id = pdPts->InsertNextPoint(out); + newIds.emplace(p, id); + } else { + id = itr->second; + } + } + + newCell->InsertNextId(id); + + } + + pd->InsertNextCell(VTK_POLYGON, newCell); + origCellIds->InsertNextValue(origId); + } + + } + +} + +void Merger::MergeGroup (const GroupType &group, PolysType &merged) { + if (group.size() == 1) { + merged.push_back(polys.at(group.back())); + + return; + } + + vtkSmartPointer pts = vtkSmartPointer::New(); + + IndexedPolysType indexedPolys; + + ReferencedPointsType refPts; + + SourcesType sources; + std::size_t src = 0; + + for (auto &index : group) { + const Poly &poly = polys.at(index); + + IndexedPoly ids; + + for (auto &p : poly) { + vtkIdType id = pts->InsertNextPoint(p.x, p.y, p.z); + + ids.push_back(id); + sources.emplace(id, src); + + refPts.emplace(id, p); + } + + indexedPolys.push_back(std::move(ids)); + src++; + } + + vtkSmartPointer kdTree = vtkSmartPointer::New(); + kdTree->OmitZPartitioning(); + kdTree->BuildLocatorFromPoints(pts); + + vtkSmartPointer linesA = vtkSmartPointer::New(); + linesA->SetPoints(pts); + linesA->Allocate(1); + + IndexedPoly::const_iterator itrA, itrB; + + for (const auto &ids : indexedPolys) { + for (itrA = ids.begin(); itrA != ids.end(); ++itrA) { + itrB = itrA+1; + if (itrB == ids.end()) { + itrB = ids.begin(); + } + + vtkIdList *line = vtkIdList::New(); + line->InsertNextId(*itrA); + line->InsertNextId(*itrB); + + linesA->InsertNextCell(VTK_LINE, line); + + line->Delete(); + } + } + + WriteVTK("linesA.vtk", linesA); + + vtkSmartPointer bspTreeA = vtkSmartPointer::New(); + bspTreeA->SetDataSet(linesA); + + int n = 0; + + PolyConnsType polyConns; + + FindConns(linesA, kdTree, bspTreeA, polyConns, indexedPolys, sources, n); + + PolyConnsType connected {{0, {}}}; + std::set restricted; // keine der conns darf im gleichen punkt beginnen + + vtkSmartPointer linesB = vtkSmartPointer::New(); + linesB->SetPoints(pts); + linesB->Allocate(1); + + vtkSmartPointer bspTreeB = vtkSmartPointer::New(); + bspTreeB->SetDataSet(linesB); + + ConnsType firstConns; + + std::size_t i, numPolys = indexedPolys.size(); + + double ptA[3], ptB[3]; + + while (connected.size() < numPolys) { + + bool foundOne = false; + + for (i = 1; i < numPolys; i++) { + if (connected.count(i) == 0) { + const ConnsType &conns = polyConns[i]; + + for (auto &conn : conns) { + if (connected.count(sources.at(conn.j)) == 1 + && restricted.count(conn.j) == 0) { + + pts->GetPoint(conn.i, ptA); + pts->GetPoint(conn.j, ptB); + + if (bspTreeB->IntersectWithLine(ptA, ptB, 1e-5, nullptr, nullptr) == 0) { + connected[sources.at(conn.i)].push_back(conn); + + // das andere poly auch aktualisieren + connected[sources.at(conn.j)].emplace_back(conn.d, conn.j, conn.i); + + restricted.insert(conn.i); + restricted.insert(conn.j); + + vtkIdList *line = vtkIdList::New(); + line->InsertNextId(conn.i); + line->InsertNextId(conn.j); + + linesB->InsertNextCell(VTK_LINE, line); + + line->Delete(); + + bspTreeB->Modified(); + + foundOne = true; + + firstConns.push_back(conn); + + break; + } + + } + } + } + } + + if (!foundOne) { + if (!FindConns(linesA, kdTree, bspTreeA, polyConns, indexedPolys, sources, n)) { + throw std::runtime_error("Merging failed."); + } + } + } + + std::map> chains; + + PolyConnsType::const_iterator itrC; + + for (itrC = connected.begin(); itrC != connected.end(); ++itrC) { + auto &chain = chains[itrC->first]; + chain.push_back(itrC->first); + + while (chain.back() != 0) { + chain.push_back(sources.at(connected.at(chain.back()).front().j)); + } + } + + std::cout << connected; + + decltype(chains)::const_iterator itrD; + + for (itrD = chains.begin(); itrD != chains.end(); ++itrD) { + std::cout << itrD->first << ": ["; + for (auto &id : itrD->second) { + std::cout << id << ", "; + } + std::cout << "]" << std::endl; + } + + std::set solved {0}; + + std::deque searchInds; + + for (i = 1; i < numPolys; i++) { + if (connected.at(i).size() == 1) { + searchInds.push_back(i); + } + } + + while (!searchInds.empty()) { + PriosType prios; + + for (auto ind : searchInds) { + PolyPriosType polyPrios; + + std::cout << "ind " << ind << std::endl; + + const Conn &first = connected.at(ind).back(); + + for (auto &conn : polyConns.at(ind)) { + auto &src = sources.at(conn.j); + + if (polyPrios.count(src) == 1) { + continue; + } + + if (conn.i != first.i + && conn.j != first.j + && restricted.count(conn.i) == 0 + && restricted.count(conn.j) == 0) { + + pts->GetPoint(conn.i, ptA); + pts->GetPoint(conn.j, ptB); + + if (bspTreeB->IntersectWithLine(ptA, ptB, 1e-5, nullptr, nullptr) == 0) { + auto &chainA = chains.at(ind), + &chainB = chains.at(src); + + std::set _chainA(chainA.begin(), chainA.end()), + _chainB(chainB.begin(), chainB.end()); + + // gemeinsame eltern + std::set shared; + + std::set_intersection(_chainA.begin(), _chainA.end(), _chainB.begin(), _chainB.end(), std::inserter(shared, shared.end())); + + // gemeinsame eltern müssen sich alle in solved befinden + if (std::includes(solved.begin(), solved.end(), shared.begin(), shared.end())) { + std::set solvable; + + std::set_difference(_chainA.begin(), _chainA.end(), solved.begin(), solved.end(), std::inserter(solvable, solvable.end())); + std::set_difference(_chainB.begin(), _chainB.end(), solved.begin(), solved.end(), std::inserter(solvable, solvable.end())); + + polyPrios.emplace(std::piecewise_construct, + std::forward_as_tuple(src), + std::forward_as_tuple(conn, solvable, -conn.d)); + } + } + } + } + + PolyPriosType::const_iterator itr; + for (itr = polyPrios.begin(); itr != polyPrios.end(); ++itr) { + prios.insert(itr->second); + } + } + + if (!prios.empty()) { + auto &prio = *prios.rbegin(); + + std::cout << "found " << prio << std::endl; + + auto &conns = connected.at(sources.at(prio.conn.i)); + + conns.push_back(prio.conn); + + connected.at(sources.at(prio.conn.j)).emplace_back(prio.conn.d, prio.conn.j, prio.conn.i); + + restricted.insert(prio.conn.i); + restricted.insert(prio.conn.j); + + vtkIdList *line = vtkIdList::New(); + line->InsertNextId(prio.conn.i); + line->InsertNextId(prio.conn.j); + + linesB->InsertNextCell(VTK_LINE, line); + + line->Delete(); + + bspTreeB->Modified(); + + solved.insert(prio.solvable.begin(), prio.solvable.end()); + + searchInds.erase(std::find(searchInds.begin(), searchInds.end(), sources.at(prio.conn.i))); + + auto itr = std::find(searchInds.begin(), searchInds.end(), sources.at(prio.conn.j)); + + if (itr != searchInds.end()) { + searchInds.erase(itr); + } + } else { + if (!FindConns(linesA, kdTree, bspTreeA, polyConns, indexedPolys, sources, n)) { + throw std::runtime_error("Merging failed."); + } + } + } + + std::cout << connected; + + WriteVTK("linesB.vtk", linesB); + + // stage 1 + + struct Cmp { + bool operator() (const Conn &a, const Conn &b) const { + return std::tie(a.i, a.j) < std::tie(b.i, b.j); + } + }; + + std::set usedConns; + + IndexedPoly polyA {indexedPolys.front()}; + + for (const auto &first : firstConns) { + auto itrA = std::find(polyA.begin(), polyA.end(), first.j); + + assert(itrA != polyA.end()); + + auto &polyB = indexedPolys.at(sources.at(first.i)); + + auto itrB = std::find(polyB.begin(), polyB.end(), first.i); + + assert(itrB != polyB.end()); + + std::rotate(polyA.begin(), itrA, polyA.end()); + std::rotate(polyB.begin(), itrB, polyB.end()); + + IndexedPoly newPoly {polyA}; + newPoly.push_back(polyA.front()); + newPoly.push_back(polyB.front()); + + newPoly.insert(newPoly.end(), polyB.rbegin(), polyB.rend()); + + polyA.swap(newPoly); + + usedConns.insert(first); + + } + + PolysType newPolysA; + GetPolys(refPts, {polyA}, newPolysA); + + WritePolys("merged_stage1.vtk", newPolysA); + + // stage 2 + + std::set leftConns; + + for (itrC = connected.begin(); itrC != connected.end(); ++itrC) { + if (itrC->first == 0) { + continue; + } + + auto &conns = itrC->second; + + ConnsType::const_iterator itr; + + for (itr = conns.begin()+1; itr != conns.end(); ++itr) { + Conn conn(0, itr->j, itr->i); + + if (usedConns.find(conn) == usedConns.end()) { + if (itr->i < itr->j) { + leftConns.emplace(0, itr->i, itr->j); + } else { + leftConns.insert(std::move(conn)); + } + } + } + + } + + std::cout << "leftConns: ["; + for (auto &conn : leftConns) { + std::cout << conn << ", "; + } + std::cout << "]" << std::endl; + + IndexedPolysType splitted {polyA}; + + for (auto &conn : leftConns) { + IndexedPolysType::iterator itr; + + for (itr = splitted.begin(); itr != splitted.end(); ++itr) { + auto &poly = *itr; + + auto itrA = std::find(poly.begin(), poly.end(), conn.i); + + if (itrA != poly.end()) { + std::rotate(poly.begin(), itrA, poly.end()); + + auto itrB = std::find(poly.begin(), poly.end(), conn.j); + + IndexedPoly newPolyA(poly.begin(), itrB+1); + IndexedPoly newPolyB(itrB, poly.end()); + + newPolyB.push_back(poly.front()); + + splitted.erase(itr); + + splitted.push_back(std::move(newPolyA)); + splitted.push_back(std::move(newPolyB)); + + break; + } + } + } + + PolysType newPolysB; + GetPolys(refPts, splitted, newPolysB); + + WritePolys("merged_stage2.vtk", newPolysB); + + std::move(newPolysB.begin(), newPolysB.end(), std::back_inserter(merged)); + +} + +bool Merger::FindConns (vtkPolyData *lines, vtkSmartPointer kdTree, vtkSmartPointer bspTree, PolyConnsType &polyConns, const IndexedPolysType &indexedPolys, const SourcesType &sources, int &n) { + + vtkPoints *pts = lines->GetPoints(); + + if (n > pts->GetNumberOfPoints()) { + return false; + } + + n += 10; + + vtkSmartPointer foundPts = vtkSmartPointer::New(); + + vtkIdType i, numPts; + + vtkIdType idB; + + vtkSmartPointer lineIds = vtkSmartPointer::New(); + + double ptA[3], ptB[3]; + + bool good; + + vtkIdType j; + vtkIdType _idA, _idB; + + std::map> _polyConns; + + vtkSmartPointer line = vtkSmartPointer::New(); + + for (const auto &ids : indexedPolys) { + for (vtkIdType idA : ids) { + pts->GetPoint(idA, ptA); + + kdTree->FindClosestNPoints(n, ptA, foundPts); + + numPts = foundPts->GetNumberOfIds(); + + for (i = 0; i < numPts; i++) { + idB = foundPts->GetId(i); + + auto srcA = sources.at(idA), + srcB = sources.at(idB); + + if (srcA == srcB) { + continue; + } + + pts->GetPoint(idB, ptB); + + good = true; + + if (bspTree->IntersectWithLine(ptA, ptB, 1e-5, nullptr, lineIds) == 1) { + for (j = 0; j < lineIds->GetNumberOfIds(); j++) { + lines->GetCellPoints(lineIds->GetId(j), line); + + _idA = line->GetId(0); + _idB = line->GetId(1); + + if (_idA != idA && _idA != idB + && _idB != idA && _idB != idB) { + + good = false; + break; + } + } + } + + if (good) { + double d = vtkMath::Distance2BetweenPoints(ptA, ptB); + + _polyConns[srcA].emplace(d, idA, idB); + _polyConns[srcB].emplace(d, idB, idA); + } + } + } + } + + decltype(_polyConns)::const_iterator itr; + + for (itr = _polyConns.begin(); itr != _polyConns.end(); ++itr) { + auto &_conns = itr->second; + + ConnsType conns(_conns.begin(), _conns.end()); + std::sort(conns.begin(), conns.end()); + + polyConns[itr->first].swap(conns); + + } + + return true; +} diff --git a/lib/creaVtk/vtkPolyDataBooleanFilter.h b/lib/creaVtk/vtkPolyDataBooleanFilter.h new file mode 100644 index 0000000..649e5f8 --- /dev/null +++ b/lib/creaVtk/vtkPolyDataBooleanFilter.h @@ -0,0 +1,298 @@ +/* +Copyright 2012-2022 Ronald Römer + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +#ifndef __vtkPolyDataBooleanFilter_h +#define __vtkPolyDataBooleanFilter_h + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include "Utilities.h" + +#define OPER_NONE 0 +#define OPER_UNION 1 +#define OPER_INTERSECTION 2 +#define OPER_DIFFERENCE 3 +#define OPER_DIFFERENCE2 4 + +enum class Capt { + NOT, + EDGE, + A, + B +}; + +enum class Side { + NONE, + START, + END +}; + +enum class Loc { + NONE, + INSIDE, + OUTSIDE +}; + +class StripPt { +public: + StripPt () : t(0), capt(Capt::NOT), catched(true) { + edge[0] = NOTSET; + edge[1] = NOTSET; + } + + double t; + Capt capt; + double captPt[3]; + + vtkIdType ind, edge[2]; + + double pt[3]; + double cutPt[3]; + + friend std::ostream& operator<< (std::ostream &out, const StripPt &s) { + out << "ind " << s.ind + << ", edge [" << s.edge[0] << ", " << s.edge[1] << "]" + << ", t " << s.t + << ", capt " << s.capt + << ", polyId " << s.polyId; + return out; + } + + std::vector history; + + vtkIdType polyId; + + bool catched; +}; + +class StripPtR { +public: + StripPtR () = delete; + + StripPtR (vtkIdType ind, std::size_t strip) : ind(ind), strip(strip), ref(NOTSET), side(Side::NONE) { + desc[0] = NOTSET; + desc[1] = NOTSET; + } + + vtkIdType ind; + std::size_t strip; + vtkIdType ref, desc[2]; + Side side; + + friend std::ostream& operator<< (std::ostream &out, const StripPtR &s) { + out << "ind " << s.ind + << ", desc [" << s.desc[0] << ", " << s.desc[1] << "]" + << ", strip " << s.strip + << ", side " << s.side + << ", ref " << s.ref; + return out; + } +}; + +typedef std::map StripPtsType; +typedef std::deque StripType; +typedef std::vector StripsType; + +class PStrips { +public: + PStrips (vtkPolyData *pd, vtkIdType cellId) { + const vtkIdType *cell; + vtkIdType numPts; + + pd->GetCellPoints(cellId, numPts, cell); + + for (vtkIdType i = 0; i < numPts; i++) { + poly.push_back(cell[i]); + } + + ComputeNormal(pd->GetPoints(), n, numPts, cell); + + base = Base(pd->GetPoints(), numPts, cell); + } + + StripPtsType pts; + StripsType strips; + + double n[3]; + IdsType poly; + Base base; +}; + +typedef std::map PolyStripsType; + +typedef std::vector> RefsType; + +// Merger + +typedef std::vector GroupType; + +typedef std::map SourcesType; + +class Conn { +public: + Conn () = delete; + Conn (double d, vtkIdType i, vtkIdType j) : d(d), i(i), j(j) {} + + double d; + vtkIdType i, j; + + bool operator< (const Conn &other) const { + return d < other.d; + } + + friend std::ostream& operator<< (std::ostream &out, const Conn &c) { + out << "Conn(d=" << c.d + << ", i=" << c.i + << ", j=" << c.j + << ")"; + return out; + } + +}; + +struct ConnCmp { + bool operator() (const Conn &a, const Conn &b) const { + return std::tie(a.i, a.j) < std::tie(b.i, b.j); + } +}; + +typedef std::vector ConnsType; +typedef std::map PolyConnsType; + +inline std::ostream& operator<< (std::ostream &out, const PolyConnsType& polyConns) { + PolyConnsType::const_iterator itr; + + for (itr = polyConns.begin(); itr != polyConns.end(); ++itr) { + out << itr->first << ": ["; + for (auto &conn : itr->second) { + out << conn << ", "; + } + out << "]" << std::endl; + } + + return out; +} + +class Prio { +public: + Prio () = delete; + Prio (const Conn &conn, const std::set &solvable, double d) : conn(conn), solvable(solvable), d(d) {} + + Conn conn; + std::set solvable; + double d; + + friend std::ostream& operator<< (std::ostream &out, const Prio &p) { + out << "Prio(conn=" << p.conn + << ", d=" << p.d + << ")"; + return out; + } +}; + +struct Cmp { + bool operator() (const Prio &a, const Prio &b) const { + const auto _a = a.solvable.size(), + _b = b.solvable.size(); + return std::tie(_a, a.d) < std::tie(_b, b.d); + } +}; + +typedef std::set PriosType; + +typedef std::map PolyPriosType; + +class Merger { + vtkPolyData *pd; + const PStrips &pStrips; + vtkIdType origId; + + PolysType polys; +public: + Merger (vtkPolyData *pd, const PStrips &pStrips, const StripsType &strips, const IdsType &descIds, vtkIdType origId); + void Run (); + +private: + void MergeGroup (const GroupType &group, PolysType &merged); + bool FindConns (vtkPolyData *lines, vtkSmartPointer kdTree, vtkSmartPointer bspTree, PolyConnsType &polyConns, const IndexedPolysType &indexedPolys, const SourcesType &sources, int &n); +}; + +class VTK_EXPORT vtkPolyDataBooleanFilter : public vtkPolyDataAlgorithm { + vtkPolyData *resultA, *resultB, *resultC; + + vtkSmartPointer modPdA, modPdB, contLines; + + vtkSmartPointer cellDataA, cellDataB; + vtkSmartPointer cellIdsA, cellIdsB; + + vtkIdTypeArray *contsA, *contsB; + + unsigned long timePdA, timePdB; + + PolyStripsType polyStripsA, polyStripsB; + + void GetStripPoints (vtkPolyData *pd, vtkIdTypeArray *sources, PStrips &pStrips, IdsType &lines); + bool GetPolyStrips (vtkPolyData *pd, vtkIdTypeArray *conts, vtkIdTypeArray *sources, PolyStripsType &polyStrips); + void RemoveDuplicates (IdsType &lines); + void CompleteStrips (PStrips &pStrips); + bool HasArea (const StripType &strip) const; + bool CutCells (vtkPolyData *pd, PolyStripsType &polyStrips); + void RestoreOrigPoints (vtkPolyData *pd, PolyStripsType &polyStrips); + void DisjoinPolys (vtkPolyData *pd, PolyStripsType &polyStrips); + void ResolveOverlaps (vtkPolyData *pd, vtkIdTypeArray *conts, PolyStripsType &polyStrips); + void AddAdjacentPoints (vtkPolyData *pd, vtkIdTypeArray *conts, PolyStripsType &polyStrips); + void MergePoints (vtkPolyData *pd, PolyStripsType &polyStrips); + bool CombineRegions (); + + int OperMode; + +public: + vtkTypeMacro(vtkPolyDataBooleanFilter, vtkPolyDataAlgorithm); + static vtkPolyDataBooleanFilter* New (); + + vtkSetClampMacro(OperMode, int, OPER_NONE, OPER_DIFFERENCE2); + vtkGetMacro(OperMode, int); + + void SetOperModeToNone () { OperMode = OPER_NONE; Modified(); } + void SetOperModeToUnion () { OperMode = OPER_UNION; Modified(); } + void SetOperModeToIntersection () { OperMode = OPER_INTERSECTION; Modified(); } + void SetOperModeToDifference () { OperMode = OPER_DIFFERENCE; Modified(); } + void SetOperModeToDifference2 () { OperMode = OPER_DIFFERENCE2; Modified(); } + +protected: + vtkPolyDataBooleanFilter (); + ~vtkPolyDataBooleanFilter (); + + int ProcessRequest (vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector) override; + + void PrintSelf (ostream&, vtkIndent) override {}; + +private: + vtkPolyDataBooleanFilter (const vtkPolyDataBooleanFilter&) = delete; + void operator= (const vtkPolyDataBooleanFilter&) = delete; + +}; + +#endif diff --git a/lib/creaVtk/vtkPolyDataContactFilter.cxx b/lib/creaVtk/vtkPolyDataContactFilter.cxx new file mode 100644 index 0000000..0e96c20 --- /dev/null +++ b/lib/creaVtk/vtkPolyDataContactFilter.cxx @@ -0,0 +1,912 @@ +/* +Copyright 2012-2022 Ronald Römer + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "vtkPolyDataContactFilter.h" +#include "Utilities.h" + +#undef DEBUG + +vtkStandardNewMacro(vtkPolyDataContactFilter); + +vtkPolyDataContactFilter::vtkPolyDataContactFilter () { + + contLines = vtkPolyData::New(); + contLines->Allocate(1000); + + contPts = vtkPoints::New(); + contPts->SetDataTypeToDouble(); + contLines->SetPoints(contPts); + + contA = vtkIdTypeArray::New(); + contB = vtkIdTypeArray::New(); + + contA->SetName("cA"); + contB->SetName("cB"); + + sourcesA = vtkIdTypeArray::New(); + sourcesA->SetNumberOfComponents(2); + + sourcesB = vtkIdTypeArray::New(); + sourcesB->SetNumberOfComponents(2); + + sourcesA->SetName("sourcesA"); + sourcesB->SetName("sourcesB"); + + neigsA = vtkIdTypeArray::New(); + neigsB = vtkIdTypeArray::New(); + + neigsA->SetName("neigsA"); + neigsB->SetName("neigsB"); + + SetNumberOfInputPorts(2); + SetNumberOfOutputPorts(3); + + invalidA = false; + invalidB = false; + + accuracy = vtkIdTypeArray::New(); + accuracy->SetName("accuracy"); +} + +vtkPolyDataContactFilter::~vtkPolyDataContactFilter () { + neigsB->Delete(); + neigsA->Delete(); + + sourcesB->Delete(); + sourcesA->Delete(); + + contB->Delete(); + contA->Delete(); + + contPts->Delete(); + contLines->Delete(); + + accuracy->Delete(); + +} + +int vtkPolyDataContactFilter::ProcessRequest (vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector) { + + if (request->Has(vtkDemandDrivenPipeline::REQUEST_DATA())) { + + vtkInformation *inInfoA = inputVector[0]->GetInformationObject(0); + vtkInformation *inInfoB = inputVector[1]->GetInformationObject(0); + + vtkPolyData *_pdA = vtkPolyData::SafeDownCast(inInfoA->Get(vtkDataObject::DATA_OBJECT())); + vtkPolyData *_pdB = vtkPolyData::SafeDownCast(inInfoB->Get(vtkDataObject::DATA_OBJECT())); + + vtkInformation *outInfoA = outputVector->GetInformationObject(0); + vtkInformation *outInfoB = outputVector->GetInformationObject(1); + vtkInformation *outInfoC = outputVector->GetInformationObject(2); + + vtkPolyData *resultA = vtkPolyData::SafeDownCast(outInfoA->Get(vtkDataObject::DATA_OBJECT())); + vtkPolyData *resultB = vtkPolyData::SafeDownCast(outInfoB->Get(vtkDataObject::DATA_OBJECT())); + vtkPolyData *resultC = vtkPolyData::SafeDownCast(outInfoC->Get(vtkDataObject::DATA_OBJECT())); + + // durchführung der aufgabe + + pdA = vtkPolyData::New(); + pdA->DeepCopy(_pdA); + + pdB = vtkPolyData::New(); + pdB->DeepCopy(_pdB); + + PreparePolyData(pdA); + PreparePolyData(pdB); + + if (pdA->GetNumberOfCells() == 0) { + vtkErrorMacro("First input does not contain any supported cells."); + return 1; + } + + if (pdB->GetNumberOfCells() == 0) { + vtkErrorMacro("Second input does not contain any supported cells."); + return 1; + } + + GetInvalidEdges(pdA, edgesA); + GetInvalidEdges(pdB, edgesB); + + // anlegen der obb-trees + + vtkOBBTree *obbA = vtkOBBTree::New(); + obbA->SetDataSet(pdA); + obbA->SetNumberOfCellsPerNode(1); + obbA->BuildLocator(); + + vtkOBBTree *obbB = vtkOBBTree::New(); + obbB->SetDataSet(pdB); + obbB->SetNumberOfCellsPerNode(1); + obbB->BuildLocator(); + + vtkMatrix4x4 *mat = vtkMatrix4x4::New(); + + obbA->IntersectWithOBBTree(obbB, mat, InterOBBNodes, this); + + if (invalidA) { + vtkErrorMacro("First input has non-manifold edges."); + return 1; + } + + if (invalidB) { + vtkErrorMacro("Second input has non-manifold edges."); + return 1; + } + + contLines->GetCellData()->AddArray(contA); + contLines->GetCellData()->AddArray(contB); + + contLines->GetCellData()->AddArray(sourcesA); + contLines->GetCellData()->AddArray(sourcesB); + + contLines->GetPointData()->AddArray(accuracy); + + contLines->RemoveDeletedCells(); + + vtkCleanPolyData *clean = vtkCleanPolyData::New(); + clean->SetInputData(contLines); + clean->ToleranceIsAbsoluteOn(); + clean->SetAbsoluteTolerance(1e-5); + clean->Update(); + + resultA->DeepCopy(clean->GetOutput()); + + vtkIdType i, numCellsA = resultA->GetNumberOfCells(); + + for (i = 0; i < numCellsA; i++) { + if (resultA->GetCellType(i) != VTK_LINE) { + resultA->DeleteCell(i); + } + } + + resultA->RemoveDeletedCells(); + + clean->Delete(); + mat->Delete(); + obbB->Delete(); + obbA->Delete(); + + resultB->DeepCopy(pdA); + resultC->DeepCopy(pdB); + + pdB->Delete(); + pdA->Delete(); + + } + + return 1; + +} + +void vtkPolyDataContactFilter::PreparePolyData (vtkPolyData *pd) { + + pd->GetCellData()->Initialize(); + pd->GetPointData()->Initialize(); + + vtkIdTypeArray *cellIds = vtkIdTypeArray::New(); + + vtkCellIterator *cellItr = pd->NewCellIterator(); + for (cellItr->InitTraversal(); !cellItr->IsDoneWithTraversal(); cellItr->GoToNextCell()) { + cellIds->InsertNextValue(cellItr->GetCellId()); + } + + vtkIdType cellId; + + for (cellItr->InitTraversal(); !cellItr->IsDoneWithTraversal(); cellItr->GoToNextCell()) { + cellId = cellItr->GetCellId(); + + if (cellItr->GetCellType() == VTK_QUAD) { + vtkIdList *ptIds = cellItr->GetPointIds(); + + vtkPoints *pts = cellItr->GetPoints(); + + double n[3]; + + ComputeNormal(pd->GetPoints(), n, 4, ptIds->GetPointer(0)); + + double dA = vtkMath::Dot(n, pts->GetPoint(0)), + dB = vtkMath::Dot(n, pts->GetPoint(1))-dA; + + if (std::abs(dB) > 1e-6) { + // nur wenn nicht auf einer ebene + + dA = vtkMath::Distance2BetweenPoints(pts->GetPoint(0), pts->GetPoint(2)); + dB = vtkMath::Distance2BetweenPoints(pts->GetPoint(1), pts->GetPoint(3)); + + vtkSmartPointer newCellA = vtkSmartPointer::New(); + vtkSmartPointer newCellB = vtkSmartPointer::New(); + + newCellA->SetNumberOfIds(3); + newCellB->SetNumberOfIds(3); + + if (dA < dB) { + newCellA->SetId(0, ptIds->GetId(1)); + newCellA->SetId(1, ptIds->GetId(2)); + newCellA->SetId(2, ptIds->GetId(3)); + + newCellB->SetId(0, ptIds->GetId(3)); + newCellB->SetId(1, ptIds->GetId(0)); + newCellB->SetId(2, ptIds->GetId(1)); + } else { + newCellA->SetId(0, ptIds->GetId(0)); + newCellA->SetId(1, ptIds->GetId(1)); + newCellA->SetId(2, ptIds->GetId(2)); + + newCellB->SetId(0, ptIds->GetId(2)); + newCellB->SetId(1, ptIds->GetId(3)); + newCellB->SetId(2, ptIds->GetId(0)); + } + + pd->InsertNextCell(VTK_TRIANGLE, newCellA); + cellIds->InsertNextValue(cellId); + + pd->InsertNextCell(VTK_TRIANGLE, newCellB); + cellIds->InsertNextValue(cellId); + + pd->DeleteCell(cellId); + } + + } else if (cellItr->GetCellType() == VTK_TRIANGLE_STRIP) { + vtkIdList *ptIds = cellItr->GetPointIds(); + + vtkCellArray *cells = vtkCellArray::New(); + + vtkTriangleStrip::DecomposeStrip(cellItr->GetNumberOfPoints(), ptIds->GetPointer(0), cells); + + vtkIdType n; + const vtkIdType *pts; + + for (cells->InitTraversal(); cells->GetNextCell(n, pts);) { + if (pts[0] != pts[1] && pts[1] != pts[2] && pts[2] != pts[0]) { + pd->InsertNextCell(VTK_TRIANGLE, n, pts); + cellIds->InsertNextValue(cellId); + } + + } + + cells->Delete(); + + pd->DeleteCell(cellId); + + } else if (cellItr->GetCellType() != VTK_TRIANGLE && cellItr->GetCellType() != VTK_POLYGON) { + pd->DeleteCell(cellId); + } + } + + cellItr->Delete(); + + cellIds->SetName("OrigCellIds"); + pd->GetCellData()->SetScalars(cellIds); + + cellIds->Delete(); + + pd->RemoveDeletedCells(); + pd->BuildLinks(); + +} + +void vtkPolyDataContactFilter::GetInvalidEdges (vtkPolyData *pd, InvalidEdgesType &edges) { + vtkSmartPointer features = vtkSmartPointer::New(); + features->SetInputData(pd); + + features->BoundaryEdgesOff(); + features->FeatureEdgesOff(); + features->ManifoldEdgesOff(); + features->NonManifoldEdgesOn(); + + features->Update(); + + vtkPolyData *lines = features->GetOutput(); + + vtkIdType num; + const vtkIdType *line; + + double ptA[3], ptB[3]; + + vtkIdType idA, idB; + + auto lineItr = vtk::TakeSmartPointer(lines->GetLines()->NewIterator()); + + for (lineItr->GoToFirstCell(); !lineItr->IsDoneWithTraversal(); lineItr->GoToNextCell()) { + lineItr->GetCurrentCell(num, line); + + lines->GetPoint(line[0], ptA); + lines->GetPoint(line[1], ptB); + + idA = pd->FindPoint(ptA); + idB = pd->FindPoint(ptB); + + edges.emplace(idA, idB); + edges.emplace(idB, idA); + + } +} + +void vtkPolyDataContactFilter::InterEdgeLine (InterPtsType &interPts, vtkPolyData *pd, vtkIdType idA, vtkIdType idB, const double *r, const double *ptA, Src src) { + double eA[3], eB[3]; + + pd->GetPoint(idA, eA); + pd->GetPoint(idB, eB); + + double ptB[3]; + vtkMath::Add(ptA, r, ptB); + + // richtungsvektor der kante bestimmen + + double e[3]; + vtkMath::Subtract(eB, eA, e); + double l = vtkMath::Normalize(e); + + double p[3]; + vtkMath::Subtract(eA, ptA, p); + + double w = std::abs(vtkMath::Determinant3x3(r, e, p)); + + if (w < 1e-4) { // ~89.995deg + // schnittpunkt ermitteln + + double v[3]; + vtkMath::Cross(r, e, v); + + double n = vtkMath::Norm(v); + + if (n > 1e-4) { // ~0.0057deg + + double s = vtkMath::Determinant3x3(p, r, v)/(n*n); + + if (s > -1e-6 && s < l+1e-6) { + double t = vtkMath::Determinant3x3(p, e, v)/(n*n); + + End end = End::NONE; + + if (s > -1e-6 && s < 1e-6) { + end = End::BEGIN; + } else if (s > l-1e-6 && s < l+1e-6) { + end = End::END; + } + + interPts.emplace_back(ptA[0]+t*r[0], ptA[1]+t*r[1], ptA[2]+t*r[2], t, idA, idB, end, src); + + } + + } else { + // parallel + + double vA[3], vB[3], cA[3], cB[3], dA, dB; + + vtkMath::Subtract(eA, ptA, vA); + vtkMath::Subtract(eA, ptB, vB); + vtkMath::Cross(vA, vB, cA); + + double dotA = vtkMath::Dot(vA, r); + + vtkMath::Subtract(eB, ptA, vA); + vtkMath::Subtract(eB, ptB, vB); + vtkMath::Cross(vA, vB, cB); + + double dotB = vtkMath::Dot(vA, r); + + dA = vtkMath::Norm(cA); + dB = vtkMath::Norm(cB); + + if (dA < 1e-4 || dB < 1e-4) { +#ifdef DEBUG + std::cout << "congruent lines with d (" << dA << ", " << dB << ") and t (" << dotA << ", " << dotB << ") and l " << l << std::endl; +#endif + + interPts.emplace_back(ptA[0]+dotA*r[0], ptA[1]+dotA*r[1], ptA[2]+dotA*r[2], dotA, idA, idB, End::BEGIN, src); + interPts.emplace_back(ptA[0]+dotB*r[0], ptA[1]+dotB*r[1], ptA[2]+dotB*r[2], dotB, idA, idB, End::END, src); + + } + + } + + } else { + // windschief + } + +} + +void vtkPolyDataContactFilter::InterPolyLine (InterPtsType &interPts, vtkPolyData *pd, vtkIdType num, const vtkIdType *poly, const double *r, const double *pt, Src src, const double *n) { + +#ifdef DEBUG + std::cout << "InterPolyLine()" << std::endl; +#endif + + std::vector edges; + edges.reserve(static_cast(num)); + + vtkIdType i, j; + + for (i = 0; i < num; i++) { + j = i+1; + + if (j == num) { + j = 0; + } + + const vtkIdType &a = poly[i], + &b = poly[j]; + + vtkPolyDataContactFilter::InterEdgeLine(interPts, pd, a, b, r, pt, src); + + edges.emplace_back(a, b); + + } + + if (interPts.empty()) { + return; + } + + struct Cmp { + bool operator() (const double &l, const double &r) const { + long a = std::lround(l*1e5), + b = std::lround(r*1e5); + + return a < b; + } + }; + + std::map paired; + + for (auto &p : interPts) { + paired[p.t].push_back(p); + } + + std::vector _paired; + + for (auto &p : paired) { + InterPtsType &pts = p.second; + + if (pts.size() == 1 && pts.front().end != End::NONE) { + // hier fehlt der zweite punkt + pts.push_back(pts.back()); + } + + _paired.push_back(pts); + } + + // trivial + + if (_paired.front().size() == 2) { + _paired.front().pop_back(); + } + + if (_paired.back().size() == 2) { + _paired.back().pop_back(); + } + + // ... + + std::map ends; + + for (const auto &pts : _paired) { + auto &last = pts.back(); + + if (last.end != End::NONE) { + ends.emplace(last.end == End::BEGIN ? last.edge.f : last.edge.g, last.t); + } + } + + double s[3], d; + + vtkMath::Cross(n, r, s); + d = vtkMath::Dot(s, pt); + + double ptA[3], ptB[3], dA, dB, vA[3], vB[3], tA, tB; + + vtkIdType id, prev, next; + + for (auto &pts : _paired) { + auto &last = pts.back(); + + if (last.end != End::NONE) { + if (last.end == End::BEGIN) { + id = last.edge.f; + next = last.edge.g; + prev = std::find_if(edges.begin(), edges.end(), [&id](const Pair &edge) { return edge.g == id; })->f; + } else { + id = last.edge.g; + prev = last.edge.f; + next = std::find_if(edges.begin(), edges.end(), [&id](const Pair &edge) { return edge.f == id; })->g; + } + + if (pts.size() == 2) { + if (ends.count(next) == 0 && ends.count(prev) == 1) { + pd->GetPoint(next, ptA); + dA = vtkMath::Dot(s, ptA)-d; + + if ((last.t > ends.at(prev) && dA > 0) || (last.t < ends.at(prev) && dA < 0)) { + // tasche + pts.pop_back(); + } + + continue; + + } else if (ends.count(next) == 1 && ends.count(prev) == 0) { + pd->GetPoint(prev, ptB); + dB = vtkMath::Dot(s, ptB)-d; + + if ((last.t > ends.at(next) && dB < 0) || (last.t < ends.at(next) && dB > 0)) { + // tasche + pts.pop_back(); + } + + continue; + + } + } + + if (ends.count(prev) == 0 && ends.count(next) == 0) { + pd->GetPoint(next, ptA); + pd->GetPoint(prev, ptB); + + dA = vtkMath::Dot(s, ptA)-d; + dB = vtkMath::Dot(s, ptB)-d; + + if (std::signbit(dA) != std::signbit(dB)) { + if (pts.size() == 2) { + pts.pop_back(); + } + + } else { + vtkMath::Subtract(ptA, pt, vA); + vtkMath::Subtract(ptB, pt, vB); + + tA = vtkMath::Dot(vA, r); + tB = vtkMath::Dot(vB, r); + + if ((tB > tA) == std::signbit(dA)) { + pts.clear(); + } + + } + } + } + } + + // ... + + InterPtsType _interPts; + + for (const auto &pts : _paired) { + _interPts.insert(_interPts.end(), pts.begin(), pts.end()); + } + + interPts.swap(_interPts); + +} + +void vtkPolyDataContactFilter::InterPolys (vtkIdType idA, vtkIdType idB) { + +#ifdef DEBUG + std::cout << "InterPolys() -> idA " << idA << ", idB " << idB << std::endl; +#endif + + vtkIdType numA, numB; + const vtkIdType *polyA, *polyB; + + pdA->GetCellPoints(idA, numA, polyA); + pdB->GetCellPoints(idB, numB, polyB); + + // ebenen aufstellen + + double nA[3], nB[3], ptA[3], ptB[3], dA, dB; + + ComputeNormal(pdA->GetPoints(), nA, numA, polyA); + ComputeNormal(pdB->GetPoints(), nB, numB, polyB); + + pdA->GetPoint(polyA[0], ptA); + pdB->GetPoint(polyB[0], ptB); + + dA = vtkMath::Dot(nA, ptA); + dB = vtkMath::Dot(nB, ptB); + + // richtungsvektor + + double r[3]; + vtkMath::Cross(nA, nB, r); + vtkMath::Normalize(r); + + // std::cout << r[0] << ", " + // << r[1] << ", " + // << r[2] << std::endl; + + // lsg. des lin. gls. mittels cramerscher regel + + int i = 0; + + for (int j = 1; j < 3; j++) { + if (std::abs(r[j]) > std::abs(r[i])) { + i = j; + } + } + + int inds[] = {1, 2}; + + if (i == 1) { + inds[0] = 0; inds[1] = 2; + } else if (i == 2) { + inds[0] = 0; inds[1] = 1; + } + + double det = nA[inds[0]]*nB[inds[1]]-nB[inds[0]]*nA[inds[1]]; + + if (std::abs(det) < 1e-12) { + return; + } + + // ein punkt auf der schnittgeraden der beiden ebenen + + double s[3]; + s[inds[0]] = (dA*nB[inds[1]]-dB*nA[inds[1]])/det; + s[inds[1]] = (nA[inds[0]]*dB-nB[inds[0]]*dA)/det; + s[i] = 0; + +#ifdef DEBUG + std::cout << "r [" << r[0] << ", " << r[1] << ", " << r[2] << "]" + << ", s [" << s[0] << ", " << s[1] << ", " << s[2] << "]" + << std::endl; +#endif + + InterPtsType intersA, intersB; + + vtkPolyDataContactFilter::InterPolyLine(intersA, pdA, numA, polyA, r, s, Src::A, nA); + vtkPolyDataContactFilter::InterPolyLine(intersB, pdB, numB, polyB, r, s, Src::B, nB); + + // probe, ob die schnittpunkte auf den kanten liegen + // bei ungenauen normalen ist das manchmal nicht der fall + + vtkPolyDataContactFilter::CheckInters(intersA, pdA, idA, idB); + vtkPolyDataContactFilter::CheckInters(intersB, pdB, idA, idB); + +#ifdef DEBUG + std::cout << "intersA " << intersA.size() + << ", intersB " << intersB.size() + << std::endl; +#endif + + + if (intersA.size() != 0 && intersB.size() != 0 + && intersA.size()%2 == 0 && intersB.size()%2 == 0) { + + AddContactLines(intersA, intersB, idA, idB); + } + +} + +void vtkPolyDataContactFilter::OverlapLines (OverlapsType &ols, InterPtsType &intersA, InterPtsType &intersB, vtkIdType idA, vtkIdType idB) { + + auto GetNeig = [](const InterPt &pA, const InterPt &pB, vtkPolyData *pd, vtkIdType polyId) -> vtkIdType { + if (pA.edge == pB.edge) { + vtkSmartPointer neigs = vtkSmartPointer::New(); + + pd->GetCellEdgeNeighbors(polyId, pA.edge.f, pA.edge.g, neigs); + + assert(neigs->GetNumberOfIds() == 1); + + return neigs->GetId(0); + } + + return NOTSET; + }; + + auto Add = [](InterPt &a, InterPt &b, InterPt &c, InterPt &d, vtkIdType neigA, vtkIdType neigB) { + a.Merge(c); + b.Merge(d); + + return std::make_tuple(a, b, neigA, neigB); + }; + + InterPtsType::iterator itr, itr2; + + vtkIdType neigA, neigB; + + for (itr = intersA.begin(); itr != intersA.end(); itr += 2) { + neigA = GetNeig(*itr, *(itr+1), pdA, idA); + + for (itr2 = intersB.begin(); itr2 != intersB.end(); itr2 += 2) { + neigB = GetNeig(*itr2, *(itr2+1), pdB, idB); + + if (itr->t <= itr2->t && (itr+1)->t > itr2->t) { + if ((itr2+1)->t < (itr+1)->t) { + ols.push_back(Add(*itr2, *(itr2+1), *itr, *(itr+1), neigA, neigB)); + } else { + ols.push_back(Add(*itr2, *(itr+1), *itr, *(itr2+1), neigA, neigB)); + } + } else if (itr2->t <= itr->t && (itr2+1)->t > itr->t) { + if ((itr+1)->t < (itr2+1)->t) { + ols.push_back(Add(*itr, *(itr+1), *itr2, *(itr2+1), neigA, neigB)); + } else { + ols.push_back(Add(*itr, *(itr2+1), *itr2, *(itr+1), neigA, neigB)); + } + } + } + } + +} + +void vtkPolyDataContactFilter::AddContactLines (InterPtsType &intersA, InterPtsType &intersB, vtkIdType idA, vtkIdType idB) { + + OverlapsType overlaps; + OverlapLines(overlaps, intersA, intersB, idA, idB); + + OverlapsType::const_iterator itr; + + for (itr = overlaps.begin(); itr != overlaps.end(); ++itr) { + auto &f = std::get<0>(*itr), + &s = std::get<1>(*itr); + +#ifdef DEBUG + std::cout << "f " << f << std::endl; + std::cout << "s " << s << std::endl; +#endif + + if (f.src == Src::A) { + if (edgesA.count(f.edge) == 1) { + invalidA = true; + } + } + + if (s.src == Src::A) { + if (edgesA.count(s.edge) == 1) { + invalidA = true; + } + } + + if (f.src == Src::B) { + if (edgesB.count(f.edge) == 1) { + invalidB = true; + } + } + + if (s.src == Src::B) { + if (edgesB.count(s.edge) == 1) { + invalidB = true; + } + } + + vtkIdList *linePts = vtkIdList::New(); + + linePts->InsertNextId(contPts->InsertNextPoint(f.pt)); + linePts->InsertNextId(contPts->InsertNextPoint(s.pt)); + + contLines->InsertNextCell(VTK_LINE, linePts); + + const vtkIdType tupleA[] = {f.srcA, s.srcA}; + const vtkIdType tupleB[] = {f.srcB, s.srcB}; + + sourcesA->InsertNextTypedTuple(tupleA); + sourcesB->InsertNextTypedTuple(tupleB); + + linePts->Delete(); + + contA->InsertNextValue(idA); + contB->InsertNextValue(idB); + + neigsA->InsertNextValue(std::get<2>(*itr)); + neigsB->InsertNextValue(std::get<3>(*itr)); + + accuracy->InsertNextValue(f.inaccurate); + accuracy->InsertNextValue(s.inaccurate); + + } + +} + +int vtkPolyDataContactFilter::InterOBBNodes (vtkOBBNode *nodeA, vtkOBBNode *nodeB, vtkMatrix4x4 *vtkNotUsed(mat), void *caller) { + vtkPolyDataContactFilter *self = reinterpret_cast(caller); + + vtkIdList *cellsA = nodeA->Cells; + vtkIdList *cellsB = nodeB->Cells; + + vtkIdType numCellsA = cellsA->GetNumberOfIds(), + numCellsB = cellsB->GetNumberOfIds(); + + vtkIdType i, j, ci, cj; + + for (i = 0; i < numCellsA; i++) { + ci = cellsA->GetId(i); + + for (j = 0; j < numCellsB; j++) { + cj = cellsB->GetId(j); + + self->InterPolys(ci, cj); + } + } + + return 0; +} + +void vtkPolyDataContactFilter::CheckInters (InterPtsType &interPts, vtkPolyData *pd, vtkIdType idA, vtkIdType idB) { + double ptA[3], + ptB[3], + v[3], + w[3], + k, + l, + alpha, + d; + + for (auto &p : interPts) { + pd->GetPoint(p.edge.f, ptA); + pd->GetPoint(p.edge.g, ptB); + + vtkMath::Subtract(ptA, ptB, v); + vtkMath::Normalize(v); + vtkMath::Subtract(ptA, p.pt, w); + + k = vtkMath::Norm(w); + l = vtkMath::Dot(v, w); + alpha = std::acos(l/k); + + if (std::isnan(alpha)) { + continue; + } + + d = std::sin(alpha)*k; + + if (d < 1e-5) { + continue; + } + + // if (p.src == Src::A) { + // std::cout << "? A "; + // } else { + // std::cout << "? B "; + // } + + // std::cout << idA << ", " << idB << ": " + // << std::fixed + // << d + // << ", [" << p.pt[0] << ", " << p.pt[1] << ", " << p.pt[2] << "], " + // << p.edge + // << std::endl; + + p.inaccurate = true; + } + +} diff --git a/lib/creaVtk/vtkPolyDataContactFilter.h b/lib/creaVtk/vtkPolyDataContactFilter.h new file mode 100644 index 0000000..544938e --- /dev/null +++ b/lib/creaVtk/vtkPolyDataContactFilter.h @@ -0,0 +1,151 @@ +/* +Copyright 2012-2022 Ronald Römer + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +*/ + +#ifndef __vtkPolyDataContactFilter_h +#define __vtkPolyDataContactFilter_h + +#include +#include +#include + +#include +#include + +#include "Utilities.h" + +class vtkOBBNode; +class vtkMatrix4x4; + +enum class Src { + A, + B +}; + +enum class End { + NONE, + BEGIN, + END +}; + +class InterPt { +public: + InterPt () = delete; + + InterPt (double x, double y, double z, double t, vtkIdType a, vtkIdType b, End end, Src src) : t(t), edge(a, b), end(end), src(src), srcA(NOTSET), srcB(NOTSET), inaccurate(false) { + pt[0] = x; + pt[1] = y; + pt[2] = z; + } + + double pt[3], t; + Pair edge; + End end; + Src src; + vtkIdType srcA, srcB; + bool inaccurate; + + friend std::ostream& operator<< (std::ostream &out, const InterPt &s) { + out << "pt [" << s.pt[0] << ", " << s.pt[1] << ", " << s.pt[2] << "]" + << ", t " << s.t + << ", edge " << s.edge + << ", end " << s.end + << ", src " << s.src; + + return out; + } + + void Merge (const InterPt &other) { + assert(src != other.src); + + if (src == Src::A) { + srcA = end == End::END ? edge.g : edge.f; + } else { + srcB = end == End::END ? edge.g : edge.f; + } + + if (std::abs(other.t-t) < 1e-5) { + if (other.src == Src::A) { + srcA = other.end == End::END ? other.edge.g : other.edge.f; + } else { + srcB = other.end == End::END ? other.edge.g : other.edge.f; + } + } + + if (other.inaccurate) { + inaccurate = true; + } + } + +}; + +typedef std::vector InterPtsType; + +typedef std::vector> OverlapsType; + +typedef std::set InvalidEdgesType; + +class VTK_EXPORT vtkPolyDataContactFilter : public vtkPolyDataAlgorithm { + + void PreparePolyData (vtkPolyData *pd); + + static void InterEdgeLine (InterPtsType &interPts, vtkPolyData *pd, vtkIdType idA, vtkIdType idB, const double *r, const double *pt, Src src); + static void InterPolyLine (InterPtsType &interPts, vtkPolyData *pd, vtkIdType num, const vtkIdType *poly, const double *r, const double *pt, Src src, const double *n); + void InterPolys (vtkIdType idA, vtkIdType idB); + void OverlapLines (OverlapsType &ols, InterPtsType &intersA, InterPtsType &intersB, vtkIdType idA, vtkIdType idB); + void AddContactLines (InterPtsType &intersA, InterPtsType &intersB, vtkIdType idA, vtkIdType idB); + + static void CheckInters (InterPtsType &interPts, vtkPolyData *pd, vtkIdType idA, vtkIdType idB); + + vtkIdTypeArray *contA, *contB; + + vtkPolyData *contLines; + vtkPoints *contPts; + + vtkPolyData *pdA, *pdB; + + vtkIdTypeArray *sourcesA, *sourcesB; + + vtkIdTypeArray *neigsA, *neigsB; + + bool invalidA, invalidB; + InvalidEdgesType edgesA, edgesB; + + void GetInvalidEdges (vtkPolyData *pd, InvalidEdgesType &edges); + + vtkIdTypeArray *accuracy; + +public: + vtkTypeMacro(vtkPolyDataContactFilter, vtkPolyDataAlgorithm); + + static vtkPolyDataContactFilter* New(); + + static int InterOBBNodes (vtkOBBNode *nodeA, vtkOBBNode *nodeB, vtkMatrix4x4 *mat, void *caller); + +protected: + vtkPolyDataContactFilter (); + ~vtkPolyDataContactFilter (); + + int ProcessRequest (vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector) override; + + void PrintSelf (ostream&, vtkIndent) override {}; + +private: + vtkPolyDataContactFilter (const vtkPolyDataContactFilter&) = delete; + void operator= (const vtkPolyDataContactFilter&) = delete; + +}; + +#endif -- 2.47.1