diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPolyhedronGeneric.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPolyhedronGeneric.h
new file mode 100644
index 0000000000000000000000000000000000000000..f367de4e8bf2bdb5e1df85dedb828c1cf8b0abfa
--- /dev/null
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPolyhedronGeneric.h
@@ -0,0 +1,147 @@
+/*
+ Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef GEOMODELKERNEL_GEOPOLYHEDRONGENERIC_H
+#define GEOMODELKERNEL_GEOPOLYHEDRONGENERIC_H
+
+/**
+ * @class: GeoPolyhedronGeneric
+ *
+ * @brief This shape represents a generic polyhedron, described by an array of vertices
+ * and faces defined by a set of ordered vertices
+ */
+
+#include <vector>
+#include "GeoModelKernel/GeoShape.h"
+
+class GeoPolyhedronGeneric : public GeoShape {
+ public:
+ // Default constructor
+ GeoPolyhedronGeneric(){};
+
+ // Constructor
+ GeoPolyhedronGeneric(const GeoShape*);
+
+ // Copy constructor
+ GeoPolyhedronGeneric(const GeoPolyhedronGeneric& gpg);
+
+ // Returns the volume of the shape, for mass inventory
+ virtual double volume() const;
+
+ // Returns the bonding box of the shape
+ virtual void extent (double& xmin, double& ymin, double& zmin,
+ double& xmax, double& ymax, double& zmax) const;
+
+ // Returns true if the shape contains the point, false otherwise.
+ virtual bool contains (double x, double y, double z) const;
+
+ // Returns true if the shape contains the point, false otherwise.
+ virtual bool contains (GeoTrf::Vector3D point, double precision = 1.e-6, bool debugMode = false) const;
+
+ // Returns the shape type, as a string
+ virtual const std::string& type() const {
+ return getClassType();
+ }
+
+ // Returns the shape type, as a coded integer
+ virtual ShapeType typeID() const{
+ return getClassTypeID();
+ }
+
+ // For type identification
+ static const std::string& getClassType() {
+ return s_classType;
+ }
+
+ // For type identification
+ static ShapeType getClassTypeID() {
+ return s_classTypeID;
+ }
+
+ // Executes a GeoShapeAction
+ virtual void exec(GeoShapeAction *action) const;
+
+ // Define vertices
+ void addVertices(std::vector<GeoTrf::Vector3D> vertices);
+
+ // Add another facet to the polyhedron. A minimum of four
+ // faces are required to create a valid solid.
+ void addFace(std::vector<size_t> face);
+
+ // normals can be uploaded or calculated
+ void addNormal(GeoTrf::Vector3D n);
+
+ // finalize GPG
+ bool construct(bool debugMode = false);
+
+ // Returns the number of vertices of the polyhedron
+ unsigned int getNVertices() const { return m_Vertices.size(); }
+
+ // Returns the number of (boundary) faces of the polyhedron
+ unsigned int getNFaces() const { return m_Faces.size(); }
+
+ // check of the hermeticity of the polyhedron (all edge need to be mirrored), plus at least four faces required, false otherwise
+ bool isValid ( bool debugMode = false) const;
+
+ // Returns coordinates of the specified vertex
+ GeoTrf::Vector3D getVertex(unsigned int i) const {
+ return m_Vertices.at(i);
+ }
+
+ // Returns definition of the specified face
+ std::vector<size_t> getFace(unsigned int i) const {
+ return m_Faces.at(i);
+ }
+
+ // return normal ( pointing outwards )
+ GeoTrf::Vector3D faceNormal(unsigned int i) const;
+
+ // true for point laying on the face
+ bool insideFace(unsigned int i, GeoTrf::Vector3D point, double precision = 1.e-6, bool debugMode = false) const;
+
+ // method to intersect ( extra-polyhedron) edge with face
+ std::vector< std::pair< double, std::pair<size_t,size_t> > > edgeFaceIntersect( int i, GeoTrf::Vector3D eInit, GeoTrf::Vector3D edge ) const;
+
+ GeoTrf::Vector3D getRefPoint(int face = -1) const;
+
+ int getOrigin() const;
+ void setOrigin( int );
+
+ void cleanAfterBooleanOperation();
+
+ // Returns true
+ virtual bool isPolyhedron () const {
+ return true;
+ }
+
+
+ virtual ~GeoPolyhedronGeneric() = default;
+protected:
+
+ private:
+ static const std::string s_classType;
+ static const ShapeType s_classTypeID;
+
+ std::vector<GeoTrf::Vector3D> m_Vertices{};
+ std::vector< std::vector<size_t> > m_Faces{};
+ GeoTrf::Vector3D m_cog{GeoTrf::Vector3D(0.,0.,0.)};
+ double m_span{0.};
+ std::vector<GeoTrf::Vector3D> m_face_cog;
+ std::vector<double> m_face_span;
+ mutable std::vector<GeoTrf::Vector3D> m_Normals{};
+
+ bool intersectLines(GeoTrf::Vector3D A, GeoTrf::Vector3D a, GeoTrf::Vector3D B, GeoTrf::Vector3D b, double& da, double& db, double& dx) const;
+
+ // fill reference span for volume and face : run automatically for transcript from GeoShape
+ void fillSpans();
+
+ // retrieve cached span value
+ double getRefSpan(int face = -1) const;
+
+ // for debugging purposes
+ int m_origin{0};
+
+};
+
+#endif
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPolyhedronManipulator.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPolyhedronManipulator.h
new file mode 100644
index 0000000000000000000000000000000000000000..c0b1fe90fcf5c84fd07f86e5f1e1361421fa3332
--- /dev/null
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPolyhedronManipulator.h
@@ -0,0 +1,211 @@
+/*
+ Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef GEOMODELKERNEL_GEOPOLYHEDRONMANIPULATOR_H
+#define GEOMODELKERNEL_GEOPOLYHEDRONMANIPULATOR_H
+
+/**
+ * @class GeoPolyhedronManipulator
+ *
+ * @brief This helper class collects methods for analytical and numerical processing of boolean shapes involving polyhedra.
+ *
+ */
+
+#include "GeoModelKernel/GeoShape.h"
+#include <GeoModelKernel/GeoDefinitions.h>
+#include "GeoModelKernel/GeoShapeShift.h"
+#include "GeoModelKernel/GeoPolyhedronGeneric.h"
+#include <string>
+
+class GeoShapeIntersection;
+class GeoShapeUnion;
+class GeoShapeSubtraction;
+class GeoPolyhedronGeneric;
+
+struct ShiftedGeoShape {
+ const GeoShape* shape{nullptr};
+ const GeoPolyhedronGeneric* phed{nullptr};
+ GeoTrf::Transform3D transform;
+ double volume{-1}; // -1 not calculated/not analytical
+ // constructor
+ ShiftedGeoShape(const GeoShape* inputShape) {
+ shape = inputShape;
+ transform = GeoTrf::Transform3D::Identity();
+ // resolve (nested) shift transform
+ while ( shape->type()=="Shift" ) {
+ const GeoShapeShift* shift = dynamic_cast<const GeoShapeShift*> (shape);
+ transform = shift->getX()*transform;
+ shape = shift->getOp();
+ }
+ if (!(shape->type()=="Union" || shape->type()=="Subtraction" || shape->type()=="Intersection") ) {
+ volume = shape->volume();
+ phed = new GeoPolyhedronGeneric(shape);
+ if (!phed->isValid()) phed=nullptr; // temporary / to deal with non-polyhedrized shapes
+ }
+ }
+};
+
+struct NPL_vertex{
+ int index; bool shared{false}; int overlapping_face{-1}; bool used{false}; int onEdge{-1}; bool onEdgeDummy{false};
+ int shared_adjacent{-1}; bool oriented{false}; // int alternative_face;
+ bool onEdgeProtected{false}; // vertex with protected line connection
+ std::vector<int> active_lines;
+ std::pair<int,int> faces; // oriented : (incoming adjacent, outgoing hit face : edge->edge, edge->line, line->edge)
+ std::pair<int,int> faces_shared_loop; // non-oriented : line->line
+ NPL_vertex() { overlapping_face = -1; used = false; shared_adjacent = 0; onEdgeDummy = false; oriented = false;}
+};
+
+struct NPL_edge{
+ std::pair<int,int> endpoints{std::make_pair<int>(-1,-1)};
+ std::pair<int,int> endpoint_faces{std::make_pair<int>(-1,-1)};
+ std::pair<int,int> endpoint_hitvx{std::make_pair<int>(-1,-1)};
+ int adjacent_face{-1}; int overlap{-1}; int generating_edge{-1};
+ bool used{false}; bool shared{false}; int symmetrized{-1};
+ bool protectedLine{false}; int dummy_overlap{-1};
+ // defoult constructor
+ NPL_edge(){}
+ // constructor
+ NPL_edge(std::pair<int,int> ep, int adjacent){ endpoints = ep; adjacent_face = adjacent; }
+};
+
+struct NPL_ordered{
+ std::vector<size_t> vtx;
+ std::vector<int> edges;
+ int top_overlap{-1}; bool shared{false}; bool oriented{false};
+ void reorder(size_t a) {
+ while ( vtx.front()!= a) {
+ vtx.push_back(vtx.front());
+ vtx.erase(vtx.begin());
+ edges.push_back(edges.front());
+ edges.erase(edges.begin());
+ }
+ }
+};
+
+struct NPH_plane{
+ int index_ph;
+ std::vector<size_t> vertices;
+ std::vector<size_t> edges;
+ std::vector<NPL_vertex> hit_vertices;
+ std::vector<std::pair< int, std::vector< NPL_vertex > > > active_planes;
+ std::vector<NPL_vertex> vlog;
+ std::vector<NPL_vertex> vlog_ori;
+ //std::vector<NPL_ordered> shared_subface;
+ std::vector<NPL_ordered> subfaces;
+ std::vector<NPL_edge> plane_edges;
+ GeoTrf::Vector3D normal;
+ std::vector<int> overlapping;
+ NPH_plane(int index) { index_ph = index;}
+};
+
+struct NPH_vertex{
+ int index{-1};
+ int index_ph{-1};
+ bool shared{false};
+ int onDoublEdge{-1}; int hitDoublEdge{-1};
+ std::pair<int,int> hit_edge_edge;
+ std::vector<std::pair<int,int> > edges;
+ NPH_vertex() {hit_edge_edge = std::pair<int,int>(-1,-1);}
+};
+
+struct NPH_vtx{ // vertex created by intersection of edge with face or edge or vertex of another polyhedron
+ int index;
+ int edge; // -1 for original vertices, generating edge otherwise
+ double distance; // along generating edge
+ int hitPlane; // index of intersected plane
+ int hitEdge; // index of intersected edge
+ int hitVertex; // index of intersected vertex
+ std::pair<int,int> faces;
+ NPH_vtx(int vtx_index) { index = vtx_index; hitPlane = -1; hitEdge = -1; hitVertex = -1; }
+};
+
+struct NPH_edge{
+ int index_ph; int edge_inverse;
+ std::pair<int,int> endpoints;
+ std::pair<int,int> endpoint_faces;
+ std::pair<int,int> faces;
+ std::vector<int> inPlane;
+ std::vector<int> doubleEdge;
+ std::vector<NPH_vtx> evtx; // vector of vertices along edge, index+relative distance ( within (0,1) interval )
+ // defoult constructor
+ NPH_edge() { endpoints = std::make_pair<int>(-1,-1); faces = std::make_pair<int>(-1,-1); endpoint_faces = std::make_pair<int>(-1,-1); edge_inverse = -1; }
+ // constructor
+ NPH_edge(std::pair<int,int> ep){ endpoints = ep; faces = std::make_pair<int>(-1,-1); endpoint_faces = std::make_pair<int>(-1,-1); edge_inverse = -1; }
+};
+
+struct NPH_line{
+ std::pair<int,int> faces;
+ std::pair<int,int> endpoints;
+ std::vector<int> evtx; // vector of vertices along edge, index+relative distance ( within (0,1) interval )
+ // constructor
+ NPH_line(std::pair<int,int> ef){ faces = ef; }
+};
+
+class GeoPolyhedronManipulator
+{
+ public:
+
+ GeoPolyhedronGeneric* intersection( const GeoPolyhedronGeneric* phA, const GeoPolyhedronGeneric* phB, bool debugMode);
+ GeoPolyhedronGeneric* subtraction( const GeoPolyhedronGeneric* phA, const GeoPolyhedronGeneric* phB, bool debugMode);
+ GeoPolyhedronGeneric* uni( const GeoPolyhedronGeneric* phA, const GeoPolyhedronGeneric* phB, bool debugMode);
+ GeoPolyhedronGeneric* resolveBoolean(const GeoPolyhedronGeneric* phA, const GeoPolyhedronGeneric* phB, int modeBoolean, bool debugMode);
+
+ private:
+
+ // mutable unsigned int m_counterNumerical = 0;
+ double m_tolerance = 1.e-5;
+ // cache input polyhedra
+ const GeoPolyhedronGeneric* m_phA{nullptr};
+ const GeoPolyhedronGeneric* m_phB{nullptr};
+
+ // method to find (distance to) the intersection of a (polyhedron) edge with (another) polyhedron face plane
+ // arguments : edge: (1) starting point, (2) endpoint-starting point, plane: (3) point , (4) normal
+ double edgeFaceIntersection( GeoTrf::Vector3D eInit, GeoTrf::Vector3D edge, GeoTrf::Vector3D facePoint, GeoTrf::Vector3D faceNormal ) const;
+
+ // auxiliary method : input edges
+ std::vector<NPH_edge> collectEdges(std::vector<NPH_plane>& planes, std::vector<NPH_vertex>& vlog, bool debugMode = false ) const;
+
+ // auxiliary method : intersections
+ void collectLines( std::vector<NPH_plane>& planes, std::vector<GeoTrf::Vector3D>& vtx, std::vector<NPH_vertex>& vlog,
+ std::vector<NPH_edge>& edges, bool debugMode ) const;
+
+ std::vector<NPL_edge> evaluateVertexLogic(int plane_index, std::vector<NPH_vertex>& vlog, std::vector<NPH_edge> edges,
+ std::vector<NPH_plane>& planes, std::vector<GeoTrf::Vector3D> vtx) const;
+
+ bool insideFace(unsigned int plane_index, GeoTrf::Vector3D point, double precision = 1.e-5, bool debugMode = false) const {
+ unsigned int nPA = m_phA->getNFaces();
+ return ( plane_index < m_phA->getNFaces() ? m_phA->insideFace(plane_index,point, precision, debugMode) :
+ m_phB->insideFace(plane_index-nPA, point, precision, debugMode) ) ;
+ }
+
+ void intersectFace(unsigned int in_edge, std::vector<NPH_edge>& edges, unsigned int plane_index, std::vector<NPH_plane>& planes, std::vector<GeoTrf::Vector3D>& vtx, std::vector<NPH_vertex>& vlog) const;
+
+ bool intersectLines(GeoTrf::Vector3D A, GeoTrf::Vector3D a, GeoTrf::Vector3D B, GeoTrf::Vector3D b, double& da, double& db, double& d) const;
+
+ bool inPlane(NPH_vertex vlog, std::vector<NPH_edge> edges, int ip) const;
+
+ NPL_vertex filterLog(NPH_vertex vlog, std::vector<NPH_edge> edges, int ip, std::vector<NPH_plane> planes) const;
+ void filterLogPlane(NPL_vertex& plv, NPH_vertex& vlog, std::vector<NPH_edge> edges, int ip, std::vector<NPH_plane> planes) const;
+
+ bool matchPlanes( int , int, std::vector<NPH_plane> planes ) const;
+
+ std::vector<int> findHitPlane( NPH_vertex vlog, int ph_index, std::vector<NPH_edge> edges, std::vector<NPH_plane> planes) const;
+
+ std::pair<int,int> hitDoubleEdge(int vx, int plane_index, int ph_index, std::vector<NPH_vertex> vlog, std::vector<NPH_edge> edges) const;
+
+ void setDoubleEdge(int ie,int ia, std::vector<NPH_edge>& edges) const;
+
+ NPL_edge symmetrize( NPL_edge line ) const;
+
+ void saveHitOnEdge(NPH_vtx nvtx, std::vector<NPH_edge>& edges) const;
+
+ NPL_edge createLine(int iv1, int iv2, int adj, NPH_plane plane) const;
+
+ NPL_ordered mergeLoops( NPL_ordered loop1, NPL_ordered loop2, std::vector< GeoTrf::Vector3D > vtx);
+
+ GeoTrf::Vector3D areaVec( std::vector<size_t> loop, std::vector< GeoTrf::Vector3D > vtx) const;
+
+};
+
+#endif
diff --git a/GeoModelCore/GeoModelKernel/src/GeoPolyhedronGeneric.cxx b/GeoModelCore/GeoModelKernel/src/GeoPolyhedronGeneric.cxx
new file mode 100755
index 0000000000000000000000000000000000000000..28e9b1825ca4e8bf3b5ed0dd510b8c1d460d0141
--- /dev/null
+++ b/GeoModelCore/GeoModelKernel/src/GeoPolyhedronGeneric.cxx
@@ -0,0 +1,648 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "GeoModelKernel/GeoPolyhedronGeneric.h"
+#include "GeoModelKernel/GeoShapeAction.h"
+#include "GeoModelKernel/GeoBox.h"
+#include "GeoModelKernel/GeoTrd.h"
+#include "GeoModelKernel/GeoTube.h"
+#include "GeoModelKernel/GeoSimplePolygonBrep.h"
+#include <cmath>
+#include <stdexcept>
+#include <iostream>
+#include <map>
+
+
+const std::string GeoPolyhedronGeneric::s_classType = "PolyhedronGeneric";
+const ShapeType GeoPolyhedronGeneric::s_classTypeID = 0x25;
+
+GeoPolyhedronGeneric::GeoPolyhedronGeneric(const GeoShape* shape) {
+
+ std::vector<GeoTrf::Vector3D> vtx;
+
+ if (shape->type() == "Box") {
+ const GeoBox* box = dynamic_cast<const GeoBox*> (shape);
+ double hx = box->getXHalfLength();
+ double hy = box->getYHalfLength();
+ double hz = box->getZHalfLength();
+
+ vtx.push_back(GeoTrf::Vector3D(-hx,-hy, hz));
+ vtx.push_back(GeoTrf::Vector3D(-hx,+hy, hz));
+ vtx.push_back(GeoTrf::Vector3D(+hx,+hy, hz));
+ vtx.push_back(GeoTrf::Vector3D(+hx,-hy, hz));
+ vtx.push_back(GeoTrf::Vector3D(+hx,-hy,-hz));
+ vtx.push_back(GeoTrf::Vector3D(+hx,+hy,-hz));
+ vtx.push_back(GeoTrf::Vector3D(-hx,+hy,-hz));
+ vtx.push_back(GeoTrf::Vector3D(-hx,-hy,-hz));
+
+ m_origin = 1;
+ }
+
+ if (shape->type() == "Trd") {
+ const GeoTrd* trd = dynamic_cast<const GeoTrd*> (shape);
+ double hx1 = trd->getXHalfLength1();
+ double hx2 = trd->getXHalfLength2();
+ double hy1 = trd->getYHalfLength1();
+ double hy2 = trd->getYHalfLength2();
+ double hz = trd->getZHalfLength();
+
+ vtx.push_back(GeoTrf::Vector3D(-hx2,-hy2, hz));
+ vtx.push_back(GeoTrf::Vector3D(-hx2,+hy2, hz));
+ vtx.push_back(GeoTrf::Vector3D(+hx2,+hy2, hz));
+ vtx.push_back(GeoTrf::Vector3D(+hx2,-hy2, hz));
+ vtx.push_back(GeoTrf::Vector3D(+hx1,-hy1,-hz));
+ vtx.push_back(GeoTrf::Vector3D(+hx1,+hy1,-hz));
+ vtx.push_back(GeoTrf::Vector3D(-hx1,+hy1,-hz));
+ vtx.push_back(GeoTrf::Vector3D(-hx1,-hy1,-hz));
+
+ m_origin = 2;
+ }
+
+ if (shape->type() == "Box" || shape->type() == "Trd") {
+
+ addVertices(vtx);
+ std::vector<size_t> f0 = {0,1,2,3};
+ addFace(f0);
+ std::vector<size_t> f1 = {4,5,6,7};
+ addFace(f1);
+ std::vector<size_t> f2 = {0,7,6,1};
+ addFace(f2);
+ std::vector<size_t> f3 = {3,4,7,0};
+ addFace(f3);
+ std::vector<size_t> f4 = {2,5,4,3};
+ addFace(f4);
+ std::vector<size_t> f5 = {1,6,5,2};
+ addFace(f5);
+
+ for (unsigned int i=0; i<getNFaces(); i++)
+ m_Normals.push_back(
+ ((m_Vertices[m_Faces[i][0]]-m_Vertices[m_Faces[i][1]]).cross(m_Vertices[m_Faces[i][2]]-m_Vertices[m_Faces[i][1]])).normalized() );
+ }
+
+ if (shape->type() == "SimplePolygonBrep") {
+
+ const GeoSimplePolygonBrep* spb = dynamic_cast<const GeoSimplePolygonBrep*> (shape);
+ double hz = spb->getDZ();
+ unsigned int nvtx = spb->getNVertices();
+
+ for (unsigned int i=0; i < nvtx; i++) vtx.push_back(GeoTrf::Vector3D(spb->getXVertex(i),spb->getYVertex(i),hz));
+ for (unsigned int i=0; i < nvtx; i++) vtx.push_back(GeoTrf::Vector3D(spb->getXVertex(i),spb->getYVertex(i),-hz));
+
+ // consistent coding requires points to be ordered clockwise
+ double clockwise = 0.;
+ for (unsigned int iv = 0; iv < nvtx; iv++) {
+ unsigned int iv1 = iv < nvtx-1 ? iv + 1 : iv + 1 - nvtx;
+ unsigned int iv2 = iv < nvtx-2 ? iv + 2 : iv + 2 - nvtx;
+ double pn = ( (vtx[iv]-vtx[iv1]).cross(vtx[iv2]-vtx[iv1]).normalized() ).z();
+ if ( std::abs(pn) > 1.e-3 ) clockwise += pn;
+ }
+
+ addVertices(vtx);
+
+ std::vector<size_t> f0; std::vector<size_t> f1;
+ if (clockwise>0) {
+ for (unsigned int i=0; i < nvtx; i++) f0.push_back(i);
+ for (unsigned int i=2*nvtx-1; i >= nvtx; i--) f1.push_back(i);
+ } else {
+ for (int i=nvtx-1; i >= 0; i--) f0.push_back(i);
+ for (int i=nvtx; i < 2*nvtx; i++) f1.push_back(i);
+ }
+ addFace(f0); addNormal(GeoTrf::Vector3D(0.,0.,1.));
+ addFace(f1); addNormal(GeoTrf::Vector3D(0.,0.,-1.));
+
+ for (unsigned int i=0; i<nvtx; i++) {
+ std::vector<size_t> f; f.clear();
+ if (clockwise>0) {
+ f.push_back( i == nvtx-1 ? 0 : i+1); f.push_back(i); f.push_back(nvtx+i); f.push_back( i == nvtx-1 ? nvtx : nvtx+i+1);
+ } else {
+ f.push_back( i == nvtx-1 ? nvtx : nvtx+i+1); f.push_back(nvtx+i); f.push_back(i); f.push_back( i == nvtx-1 ? 0 : i+1);
+ }
+ addFace(f); addNormal( (m_Vertices[f[0]]-m_Vertices[f[1]]).cross(m_Vertices[f[2]]-m_Vertices[f[1]]).normalized() );
+ }
+
+ fillSpans();
+
+ // double-check normals and consistent ordering of faces
+ for (unsigned int ip=2; ip<getNFaces(); ip++) {
+ if (m_Normals[ip].norm()!=1) {
+ for (unsigned int iv=1; iv<m_Faces[ip].size()-2; iv++) {
+ // recalculate normal using other vertices
+ GeoTrf::Vector3D tn = ((m_Vertices[m_Faces[ip][iv]]-m_Vertices[m_Faces[ip][iv+1]]).cross(m_Vertices[m_Faces[ip][iv+2]]-m_Vertices[m_Faces[ip][iv+1]])).normalized();
+ if (tn.norm()>0) { m_Normals[ip] = tn; break; }
+ }
+ }
+ }
+
+ if (!isValid() ) std::cout << " GeoPolyhedronGeneric constructor from SimplePolygonBrep not valid" << std::endl;
+
+ m_origin = 3;
+ }
+
+ if (shape->type() == "Tube") { // polyhedrize : TODO : make the precision tunable
+
+ const GeoTube* tube = dynamic_cast<const GeoTube*> (shape);
+ double hz = tube->getZHalfLength();
+ double rmin = tube->getRMin();
+ double rmax = tube->getRMax();
+
+ int nsplit = 16;
+
+ // build outer shell
+ double ff = acos(-1.)/nsplit;
+ double rx = rmax * sqrt( ff / sin(ff) / cos(ff) );
+ for (unsigned int i=0; i < nsplit; i++) vtx.push_back(GeoTrf::Vector3D(rx*cos(-i*2*ff),rx*sin(-i*2*ff),hz));
+ for (unsigned int i=0; i < nsplit; i++) vtx.push_back(GeoTrf::Vector3D(rx*cos(-i*2*ff),rx*sin(-i*2*ff),-hz));
+
+ std::vector<size_t> f0; std::vector<size_t> f1;
+
+ for (unsigned int i=0; i < nsplit; i++) f0.push_back(i);
+ for (unsigned int i=2*nsplit-1; i >= nsplit; i--) f1.push_back(i);
+
+ if (rmin>0) {
+
+ f0.push_back(0); f1.push_back(2*nsplit-1); // needed for combined loop
+
+ double rn = rmin * sqrt( ff / sin(ff) / cos(ff) );
+ for (unsigned int i=0; i < nsplit; i++) vtx.push_back(GeoTrf::Vector3D(rn*cos(-i*2*ff),rn*sin(-i*2*ff),-hz));
+ for (unsigned int i=0; i < nsplit; i++) vtx.push_back(GeoTrf::Vector3D(rn*cos(-i*2*ff),rn*sin(-i*2*ff),hz));
+
+ for (unsigned int i=2*nsplit; i < 3*nsplit; i++) f1.push_back(i);
+ for (unsigned int i=4*nsplit-1; i >= 3*nsplit; i--) f0.push_back(i);
+
+ f0.push_back(4*nsplit-1); f1.push_back(2*nsplit); // needed for combined loop
+ }
+
+ addFace(f0); addNormal(GeoTrf::Vector3D(0.,0.,1.));
+ addFace(f1); addNormal(GeoTrf::Vector3D(0.,0.,-1.));
+
+ addVertices(vtx);
+
+ for (unsigned int i=0; i<nsplit; i++) {
+ std::vector<size_t> f; f.clear();
+ f.push_back( i == nsplit-1 ? 0 : i+1); f.push_back(i); f.push_back(nsplit+i); f.push_back( i == nsplit-1 ? nsplit : nsplit+i+1);
+ addFace(f); addNormal( (m_Vertices[f[0]]-m_Vertices[f[1]]).cross(m_Vertices[f[2]]-m_Vertices[f[1]]).normalized() );
+ }
+
+ if (rmin>0) {
+ std::vector<size_t> f; f.clear();
+ f.push_back(3*nsplit); f.push_back(0); f.push_back(nsplit); f.push_back(2*nsplit);
+ addFace(f); addNormal( (m_Vertices[f[0]]-m_Vertices[f[1]]).cross(m_Vertices[f[2]]-m_Vertices[f[1]]).normalized() );
+
+ for (unsigned int i=0; i<nsplit; i++) {
+ f.clear();
+ f.push_back( i == 3*nsplit-1 ? 2*nsplit : 2*nsplit+i+1); f.push_back(2*nsplit+i); f.push_back(3*nsplit+i); f.push_back( i == 3*nsplit-1 ? 3*nsplit : 3*nsplit+i+1);
+ addFace(f); addNormal( (m_Vertices[f[0]]-m_Vertices[f[1]]).cross(m_Vertices[f[2]]-m_Vertices[f[1]]).normalized() );
+ }
+
+ f.clear();
+ f.push_back(2*nsplit); f.push_back(nsplit); f.push_back(0); f.push_back(3*nsplit);
+ addFace(f); addNormal( (m_Vertices[f[0]]-m_Vertices[f[1]]).cross(m_Vertices[f[2]]-m_Vertices[f[1]]).normalized() );
+ }
+
+ fillSpans();
+
+ std::cout << "tube->polyhedron:rmin,rmax,hz:"<<rmin <<","<<rmax<<","<<hz<<":valid:" << this->isValid() <<":volume:" << this->volume() << std::endl;
+
+ m_origin = 4;
+
+ }
+
+ if (m_face_span.size()<m_Faces.size()) fillSpans();
+
+}
+
+GeoPolyhedronGeneric::GeoPolyhedronGeneric(const GeoPolyhedronGeneric& gpg) {
+
+ m_Vertices = gpg.m_Vertices;
+ m_Faces = gpg.m_Faces;
+ m_Normals = gpg.m_Normals;
+ m_origin = gpg.m_origin; // only used for debugging of initial conversion from GeoShape
+
+ if (isValid()) fillSpans();
+
+}
+
+double GeoPolyhedronGeneric::volume () const {
+ if (!isValid())
+ throw std::runtime_error ("Volume requested for incomplete polyhedron");
+ //
+ GeoTrf::Vector3D refVtx = m_Vertices[0]; // reference point
+ double vol = 0.;
+ std::vector<GeoTrf::Vector3D> vtx;
+
+ for (unsigned int i=0; i< getNFaces(); i++) {
+
+ // retrieve vertices
+ std::vector<size_t> faceDef = getFace(i);
+ vtx.clear();
+ for (auto vi : faceDef ) vtx.push_back(m_Vertices[vi]);
+ // normal
+ GeoTrf::Vector3D normal = faceNormal(i);
+
+ // area
+ int n = vtx.size();
+ GeoTrf::Vector3D area = vtx[0].cross(vtx[n - 1]);
+ for (int k = 1; k < n; ++k) area += vtx[k].cross(vtx[k-1]);
+ vol += 0.5 * area.norm() * (vtx[1]-refVtx).dot(normal) / 3.;
+ }
+
+ return std::abs(vol);
+}
+
+void GeoPolyhedronGeneric::extent (double& xmin, double& ymin, double& zmin,
+ double& xmax, double& ymax, double& zmax) const
+{
+ if (!isValid())
+ throw std::runtime_error ("Extent requested for incomplete polyhedron");
+ xmin = xmax = m_Vertices[0].x();
+ ymin = ymax = m_Vertices[0].y();
+ zmin = zmax = m_Vertices[0].z();
+ for (size_t k = 1; k < m_Vertices.size(); ++k)
+ {
+ double x = m_Vertices[k].x();
+ double y = m_Vertices[k].y();
+ double z = m_Vertices[k].z();
+ xmin = std::min(xmin, x);
+ xmax = std::max(xmax, x);
+ ymin = std::min(ymin, y);
+ ymax = std::max(ymax, y);
+ zmin = std::min(zmin, z);
+ zmax = std::max(zmax, z);
+ }
+}
+
+bool GeoPolyhedronGeneric::contains (double x, double y, double z) const
+{
+ return this->contains( GeoTrf::Vector3D(x,y,z) );
+}
+
+bool GeoPolyhedronGeneric::contains ( GeoTrf::Vector3D point, double precision, bool debugMode ) const
+{
+ // TODO : make the method independent from the direction of normals ( so that it can be used for their validation )
+ if (debugMode) std::cout <<"test inside for point:" << point.x()<<"," <<point.y()<<"," << point.z()
+ <<":distance from ref:"<< ( point - getRefPoint() ).norm() <<" compare with span:" << getRefSpan() << std::endl;
+
+ // compare with object extent
+ if ( ( point - getRefPoint() ).norm() > getRefSpan() ) return false;
+
+ // algorithm : count intersection with faces along probe line : point+(0,0,1)*t
+ // odd/even number of intersections : inside/outside
+
+ bool ineg = false; bool ipos = false;
+ GeoTrf::Vector3D probeDir(0.,0.,1.);
+ std::vector<GeoTrf::Vector3D> intersections; // keep trace of intersections to intervene when an edge is hit
+
+ for (auto vi : m_Vertices) if ( (vi - point).norm() < precision) return true;
+
+ if (debugMode) for (auto vi : m_Vertices) std::cout <<"input vertices:" << vi.x() <<"," << vi.y() <<"," << vi.z() << std::endl;
+
+ double h,pn,dist;
+ for (size_t i=0; i< getNFaces() ; i++) {
+ std::vector<size_t> face= getFace(i);
+ if (face.size() < 3) continue;
+ // check if point lays within the face
+ GeoTrf::Vector3D n = this->faceNormal(i);
+
+ h = ( m_Vertices[m_Faces[i][0]] - point).dot(n); // nearest distance point - plane
+
+ if ( std::abs(h) < precision ) {
+ if ( insideFace(i, point, precision) ) {
+ if (debugMode) std::cout << "point at face:"<< i << std::endl;
+ return true;
+ }
+ }
+
+ pn = probeDir.dot(n); // dot product normal*probe direction
+ //if (debugMode) std::cout <<" face normal:" << i <<":" << n.x() <<"," << n.y()<<"," << n.z() << std::endl;
+ if ( std::abs(pn) > 1.e-9 ) {
+ dist = h / pn;
+ GeoTrf::Vector3D ix = point + dist*probeDir ;
+
+ if ( insideFace(i, GeoTrf::Vector3D( ix ), precision, debugMode ) ) {
+ //if (debugMode)
+ // std::cout << "plane intersection:" << (point+dist*probeDir).x() <<"," <<(point+dist*probeDir).y() <<"," <<(point+dist*probeDir).z()
+ // <<": at distance:" << dist<<":valid:" << validIntersection << std::endl;
+ if ( std::abs(dist) < precision ) {
+ if (debugMode) std::cout << "point at face:"<< i << std::endl;
+ return true;
+ }
+ bool edgeHit = false;
+ for ( auto px : intersections ) if ((ix - px).norm() < precision ) edgeHit = true;
+ if (!edgeHit) {
+ if (dist > precision) ipos ^= true;
+ else if (dist < -precision) ineg ^= true;
+ if (debugMode) std::cout << "valid intersection with face:" << i << ":at distance:" << dist <<":pos:" << ipos <<":neg:" << ineg << std::endl;
+ intersections.push_back(ix);
+ }
+ }
+ }
+ }
+
+ return ipos && ineg; // TODO check if not compatible
+}
+
+void GeoPolyhedronGeneric::addVertices(std::vector<GeoTrf::Vector3D> vertices)
+{
+ m_Vertices = vertices;
+}
+
+void GeoPolyhedronGeneric::addFace(std::vector<size_t> faceVtx)
+{
+ m_Faces.push_back(faceVtx);
+}
+
+void GeoPolyhedronGeneric::addNormal(GeoTrf::Vector3D n)
+{
+ m_Normals.push_back(n);
+}
+
+
+GeoTrf::Vector3D GeoPolyhedronGeneric::faceNormal(unsigned int i) const
+{
+
+ if (!m_Normals.size() ) { // recalculate normals
+ for (unsigned int i=0; i<getNFaces(); i++) {
+ m_Normals.push_back(
+ ( (m_Vertices[m_Faces[i][0]]-m_Vertices[m_Faces[i][1]]).cross(m_Vertices[m_Faces[i][2]]-m_Vertices[m_Faces[i][1]])).normalized() ) ;
+ // check orientation
+ for (unsigned int iv=0; iv<m_Faces[i].size()-2; iv++) {
+ GeoTrf::Vector3D ref = 0.5*(m_Vertices[m_Faces[i][iv]] + 0.5*(m_Vertices[m_Faces[i][iv+1]]+m_Vertices[m_Faces[i][iv+2]]));
+ if (insideFace(i,ref)) {
+ bool cn = contains( ref+1.e-3*m_Normals.back());
+ if ( cn != contains( ref-1.e-3*m_Normals.back()) ) { // valid test point
+ if (cn) m_Normals.back() = -1.* m_Normals.back();
+ break;
+ }
+ }
+ }
+ }
+ }
+
+ return m_Normals[i];
+}
+
+bool GeoPolyhedronGeneric::insideFace(unsigned int i, GeoTrf::Vector3D point, double precision, bool debugMode) const
+{
+ std::vector<GeoTrf::Vector3D> vtx;
+ std::vector<size_t> face= getFace(i);
+ for (auto vi : face ) vtx.push_back(m_Vertices[vi]);
+
+ if ( std::abs(( point - vtx[0]).dot(faceNormal(i) ) ) > precision ) {
+ if (debugMode) std::cout << "insideFace called for point not on plane" << std::endl;
+ return false;
+ }
+
+ GeoTrf::Vector3D probe = getRefPoint(i) - point ;
+
+ // check span
+ if ( probe.norm() > getRefSpan() ) {
+ if (debugMode) std::cout << "insideFace called for point outside span" << std::endl;
+ return false;
+ }
+
+ // vertex
+ for (auto vx : vtx) if ( (vx - point).norm()< precision ) {
+ if (debugMode) std::cout << "insideFace vertex hit" << std::endl;
+ return true;
+ }
+
+ // save valid intersections in order to handle case when probe hits a vertex ( double intersection )
+ std::vector<GeoTrf::Vector3D> valid_intersections;
+
+ if ( probe.norm() < precision) return true; // on edge
+
+ bool in = false;
+ double d,dist,t;
+
+ if (debugMode) {
+ std::cout << "insideFace:" << i <<":test point:" << point.x() <<"," << point.y() <<"," << point.z() << std::endl;
+ std::cout << "insideFace:probe:" << probe.x() <<"," << probe.y() <<"," << probe.z() <<":norm:" << probe.norm()<< std::endl;
+ for ( auto iv : vtx ) std::cout << "face vertices:" << iv.x() <<"," << iv.y() <<"," << iv.z() << std::endl;
+ }
+
+ for (unsigned int iv = 0; iv < vtx.size(); iv++) {
+
+ GeoTrf::Vector3D edge = (iv<vtx.size()-1) ? vtx[iv+1] - vtx[iv] : vtx[0] - vtx[iv];
+ double ee = edge.norm();
+ if (debugMode) std::cout << "insideFace:edge:" << edge.x() <<"," << edge.y() <<"," << edge.z() <<":norm:" << ee<< std::endl;
+
+ // dist and t in units (probe, edge)
+ bool lineIntersect = intersectLines( point, probe, vtx[iv], edge, dist, t, d);
+ if (!lineIntersect && std::abs(d) > precision ) continue;
+
+ if (!lineIntersect) { // point on edge line
+ double td = ( point - vtx[iv] ).dot(edge)/ee/ee;
+ if (td*ee > -precision && (td-1.)*ee < precision) return true; // does not count as intersection otherwise
+ }
+ if (!lineIntersect) continue;
+ //
+ if (debugMode) std::cout << "insideFace:intersecting line:" <<iv <<":d:"<< d << ":dist:" << dist <<":t:" << t << std::endl;
+ if ( std::abs(d)< 1.e-3 && t*ee >- precision && (t-1.)*ee < precision ) { // TODO calculate error on d to remove hardcoded tolerance
+ //
+ if (fabs(dist*probe.norm()) < precision ) return true;
+ // check if not found already
+ GeoTrf::Vector3D xn = vtx[iv]+t*edge;
+ bool known = false; for ( auto x : valid_intersections ) if ( (xn -x).norm() < precision ) known = true;
+ if (!known) {
+ in ^= dist>0;
+ valid_intersections.push_back(xn);
+ if (debugMode) std::cout <<"insideFace:collecting valid intersection:" << xn.x() <<"," << xn.y() <<"," << xn.z()<< ":insideFace current:" << in << std::endl;
+ }
+ if ( (xn-point).norm() < precision ) return true; // point on edge
+ }
+
+ }
+
+ if (debugMode) std::cout << "insideFace returns:" << in << std::endl;
+
+ return in;
+
+}
+
+void GeoPolyhedronGeneric::exec(GeoShapeAction *action) const
+{
+ action->handleShape(this);
+}
+
+
+bool GeoPolyhedronGeneric::intersectLines(GeoTrf::Vector3D A, GeoTrf::Vector3D a, GeoTrf::Vector3D B, GeoTrf::Vector3D b, double& da, double& db, double& dx) const
+{
+ // lines parametrization : A + a*da , B + b*db ; a,b define "units" of distances da, db
+ // return value : closest distance between lines in native [mm] units
+
+ // scalars
+ double ab = a.dot(b); double a2 = a.dot(a); double b2 = b.dot(b);
+ double ABa = (A-B).dot(a);
+ double ABb = (A-B).dot(b);
+
+
+ // special cases
+ if ( std::abs(ab)/sqrt(a2*b2)<1.e-9 ) { // to ensure sufficient precision
+ da = -ABa / a2;
+ db = ABb / b2;
+ dx = (A + da*a -B - db*b).norm();
+ return true;
+ }
+
+ if ( a2*b2 == ab*ab ) { // parallel lines
+ da = 0; db = 0;
+ dx = ((A-B).cross(b)).norm()/b.norm();
+ return false;
+ }
+
+ da = (ABb*ab -b2*ABa) / (a2*b2 - ab*ab);
+ db = (ABa + a2*da) / ab;
+ dx = (A + da*a -B - db*b).norm();
+ return true;
+
+}
+
+// this is obsolete
+std::vector< std::pair< double, std::pair<size_t,size_t> > > GeoPolyhedronGeneric::edgeFaceIntersect( int i, GeoTrf::Vector3D eInit, GeoTrf::Vector3D edge ) const
+{
+ std::vector< std::pair< double, std::pair<size_t,size_t> > > output;
+
+ std::vector<GeoTrf::Vector3D> vtx;
+ std::vector<size_t> face= getFace(i);
+ for (auto vi : face ) vtx.push_back(m_Vertices[vi]);
+
+ double d,dist,t;
+ for (unsigned int iv = 0; iv < vtx.size(); iv++) {
+
+ std::pair<size_t,size_t> fvx(face[iv], (iv<face.size()-1) ? face[iv+1] : face[0]);
+
+ GeoTrf::Vector3D fi = (iv<vtx.size()-1) ? vtx[iv+1] - vtx[iv] : vtx[0] - vtx[iv];
+
+ // dist and t in units (probe, edge)
+ bool lineIntersect = intersectLines( vtx[iv], fi, eInit, edge, dist, t, d);
+
+ if (lineIntersect) {
+ if ( std::abs(d)<1.e-6 && dist>=0 && dist <= 1 && t>=0. && t<=1.) {
+ std::pair< double, std::pair<size_t,size_t> > ix(t, fvx);
+ output.push_back(ix);
+ }
+ }
+ }
+ return output;
+}
+
+
+bool GeoPolyhedronGeneric::isValid( bool debugMode ) const
+{
+ if (getNFaces()<4) return false;
+
+ // loop over edges, check that each appears with its mirrored version ( attached faces )
+ for ( unsigned int ic=0; ic < getNFaces(); ic++) {
+ std::vector<size_t> vtx = getFace(ic);
+ if (vtx.size()<2) continue; // should not happen though
+ for (unsigned int iv=0; iv<vtx.size(); iv++) {
+ std::pair<size_t,size_t> edge(vtx[iv], iv<vtx.size()-1 ? vtx[iv+1] : vtx[0]); // probe
+ bool mirrored = false; // search mirrored edge
+ for ( unsigned int icm=0; icm < getNFaces(); icm++) {
+ std::vector<size_t> vtm = getFace(icm);
+ if (vtm.size()<2) continue; // should not happen though
+ for (unsigned int ivm=0; ivm<vtm.size(); ivm++) {
+ std::pair<size_t,size_t> edgm(vtm[ivm], ivm<vtm.size()-1 ? vtm[ivm+1] : vtm[0]); // test
+ if ( edgm.first == edge.second && edgm.second == edge.first) mirrored = true;
+ if (mirrored) break;
+ } // end loop over other edges
+ if (mirrored) break;
+ } // end loop over other faces
+ if (!mirrored) {
+ if (debugMode) std::cout << "problem in polyhedron, edge not mirrored:face:edge:" << ic<<":" <<iv << std::endl;
+ return false;
+ }
+ }
+ }
+
+ return true;
+}
+
+void GeoPolyhedronGeneric::fillSpans() {
+
+ // calculate reference points and span for volumes
+ m_cog = GeoTrf::Vector3D(0.,0.,0.); m_span = 0.;
+ if (m_Vertices.size()) {
+ for ( auto vx : m_Vertices ) m_cog = m_cog + vx; m_cog = 1./m_Vertices.size()*m_cog;
+ for ( auto vx : m_Vertices ) m_span = std::max( m_span, (vx-m_cog).norm() );
+ }
+ // and faces
+ m_face_cog.clear(); m_face_span.clear();
+ for ( auto fi : m_Faces ) { GeoTrf::Vector3D fc(0.,0.,0.); double fs=0.;
+ if (fi.size()) {
+ for ( auto vi : fi ) fc = fc + m_Vertices[vi]; fc = 1./fi.size() * fc;
+ for ( auto vi : fi) fs = std::max( fs, (m_Vertices[vi] - fc).norm() );
+ }
+ m_face_cog.push_back(fc); m_face_span.push_back(fs);
+ }
+}
+
+GeoTrf::Vector3D GeoPolyhedronGeneric::getRefPoint(int face) const {
+
+ // first check that info available
+ if ( m_face_cog.size() < m_Faces.size() ) {
+ std::cout << "WARNING: GeoPolyhedronGeneric object not constructed: run ::construct() first"<< std::endl;
+ return GeoTrf::Vector3D(1.e5,1.e5,1.e5); // dummy return value
+ }
+
+ if (face<0) return m_cog;
+ else if ( face < m_face_cog.size() ) return m_face_cog[face];
+ else return GeoTrf::Vector3D(1.e5,1.e5,1.e5); // return dummy
+
+}
+
+double GeoPolyhedronGeneric::getRefSpan(int face) const {
+
+ // first check that info available
+ if ( m_face_span.size() < m_Faces.size() ) {
+ std::cout << "WARNING: GeoPolyhedronGeneric object not constructed: run ::construct() first"<< std::endl;
+ return 1.e5; // dummy return value
+ }
+
+ if (face<0) return m_span;
+ else if ( face < m_face_span.size() ) return m_face_span[face];
+ else return 1.e5; // return dummy
+
+}
+
+int GeoPolyhedronGeneric::getOrigin() const { return m_origin; }
+
+void GeoPolyhedronGeneric::setOrigin( int ori ) { m_origin = ori; }
+
+void GeoPolyhedronGeneric::cleanAfterBooleanOperation() {
+
+ // remove dummy faces ( less than 3 vertices )
+ std::vector<std::vector<size_t> > nfaces;
+ std::vector< GeoTrf::Vector3D > normals;
+ for (unsigned ip = 0; ip < m_Faces.size(); ip++) {
+ if (m_Faces[ip].size()>2) {
+ nfaces.push_back(m_Faces[ip]);
+ normals.push_back(m_Normals[ip]);
+ }
+ }
+
+ // removing unused vertices
+ std::vector<GeoTrf::Vector3D> nvtx;
+ std::map<size_t,size_t> m;
+ for (unsigned ip = 0; ip < nfaces.size(); ip++) {
+ for (unsigned iv = 0; iv < nfaces[ip].size(); iv++ ) {
+ auto const& pair = m.try_emplace(nfaces[ip][iv], nvtx.size());
+ if (pair.second) nvtx.push_back(m_Vertices[nfaces[ip][iv]]);
+ nfaces[ip][iv] = (*pair.first).second;
+ }
+ }
+
+ m_Vertices = nvtx;
+ m_Faces = nfaces;
+ m_Normals = normals;
+
+ fillSpans();
+}
+
+bool GeoPolyhedronGeneric::construct(bool debugMode) {
+
+ fillSpans();
+ return isValid(debugMode);
+
+}
diff --git a/GeoModelCore/GeoModelKernel/src/GeoPolyhedronManipulator.cxx b/GeoModelCore/GeoModelKernel/src/GeoPolyhedronManipulator.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..b1632ff9d48df1bb948d1a2579db2b30203a5792
--- /dev/null
+++ b/GeoModelCore/GeoModelKernel/src/GeoPolyhedronManipulator.cxx
@@ -0,0 +1,1564 @@
+/*
+ Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "GeoModelKernel/GeoPolyhedronManipulator.h"
+#include "GeoModelKernel/GeoShape.h"
+#include "GeoModelKernel/GeoShapeIntersection.h"
+#include "GeoModelKernel/GeoShapeUnion.h"
+#include "GeoModelKernel/GeoShapeSubtraction.h"
+#include "GeoModelKernel/GeoBox.h"
+#include <cmath>
+#include <cstdint>
+#include <iostream>
+#include <utility>
+
+GeoPolyhedronGeneric* GeoPolyhedronManipulator::intersection( const GeoPolyhedronGeneric* phA, const GeoPolyhedronGeneric* phB, bool debugMode)
+{ return resolveBoolean( phA, phB, 0, debugMode); }
+
+GeoPolyhedronGeneric* GeoPolyhedronManipulator::subtraction( const GeoPolyhedronGeneric* phA, const GeoPolyhedronGeneric* phB, bool debugMode)
+{ return resolveBoolean( phA, phB, 1, debugMode); }
+
+GeoPolyhedronGeneric* GeoPolyhedronManipulator::uni( const GeoPolyhedronGeneric* phA, const GeoPolyhedronGeneric* phB, bool debugMode)
+{ return resolveBoolean( phA, phB, 2, debugMode); }
+
+GeoPolyhedronGeneric* GeoPolyhedronManipulator::resolveBoolean(const GeoPolyhedronGeneric* phA, const GeoPolyhedronGeneric* phB, int modeBoolean, bool debugMode) {
+ // modeBoolean : 0 intersection (A and B): shared(A,B) non-overlapping + shared A overlapping
+ // 1 subtraction (A and not B): non-shared(A) + inverted shared (B)
+ // 2 union (A or B) : non-shared(A,B) + shared A overlapping
+
+ // check inputs
+ if ( !phA->isValid() || !phB->isValid()) {
+ return nullptr;
+ } else {
+ m_phA = phA; m_phB = phB;
+ }
+ // collect all planes / vertices
+ size_t nVA = phA->getNVertices(); size_t nPA = phA->getNFaces();
+ std::vector<GeoTrf::Vector3D> vtx; std::vector< NPH_vertex > vlog; std::vector< NPH_plane > planes;
+ for (unsigned int iv=0; iv < nVA; iv++) {
+ vtx.push_back( phA->getVertex(iv) );
+ NPH_vertex phv; phv.index=vlog.size(); phv.index_ph=0; phv.shared = phB->contains( vtx.back(), 1.e-3 ); vlog.push_back( phv );
+ }
+ for (unsigned int ip=0; ip < phA->getNFaces(); ip++) {
+ planes.push_back(NPH_plane(0));
+ planes.back().vertices = phA->getFace(ip);
+ planes.back().normal = phA->faceNormal(ip);
+ if (debugMode) std::cout <<"phA:face:" << ip <<":vtx:"; for ( auto jv : planes.back().vertices ) std::cout << jv<<","; std::cout <<"" << std::endl;
+ }
+ // phB + overlapping vertices + overlapping planes
+ std::vector<std::pair<int,int>> overlapping_vertices;
+ for (unsigned int iv=0; iv < phB->getNVertices(); iv++) {
+ GeoTrf::Vector3D new_vtx = phB->getVertex(iv);
+ for (unsigned int jv=0; jv<vtx.size(); jv++) if ( (vtx[jv]-new_vtx).norm() < m_tolerance ) overlapping_vertices.push_back(std::pair<int,int>(iv+nVA,jv));
+ vtx.push_back( new_vtx );
+ NPH_vertex phv; phv.index = vlog.size(); phv.index_ph = 1; phv.shared = phA->contains( vtx.back(), 1.e-3 ); vlog.push_back( phv );
+ }
+ //
+ for (unsigned int ip=0; ip < phB->getNFaces(); ip++) {
+ planes.push_back(NPH_plane(1));
+ std::vector<size_t> vtx_ph = phB->getFace(ip); // auxiliary
+ for (auto vh : vtx_ph ) {
+ int overlap = -1; for (auto ov : overlapping_vertices) if ( ov.first == vh+nVA ) overlap = ov.second;
+ planes.back().vertices.push_back( overlap>=0 ? overlap : vh+nVA);
+ }
+ planes.back().normal = phB->faceNormal(ip);
+ if (debugMode) std::cout <<"phB:face:" << ip <<":vtx:"; for ( auto jv : planes.back().vertices ) std::cout << jv<<","; std::cout <<"" << std::endl;
+ // overlapping planes
+ for (unsigned int ipA=0; ipA < nPA; ipA++) { NPH_plane plA = planes[ipA];
+ double nn = plA.normal.dot(phB->faceNormal(ip));
+ if (std::abs(std::abs( nn ) - 1.) < m_tolerance ) {
+ if ( std::abs( (vtx[plA.vertices[0]]-vtx[planes.back().vertices[0]]).dot(plA.normal)) < m_tolerance ) {
+ planes[ipA].overlapping.push_back( nn>0 ? nPA+ip : nPA+ip+1000);
+ planes[nPA+ip].overlapping.push_back( nn>0 ? ipA : ipA+1000);
+ }
+ }
+ }
+ }
+
+ if (debugMode) {
+ for (unsigned int iv=0; iv<vlog.size(); iv++) if (vlog[iv].shared)
+ std::cout << "original shared vertices:" << iv <<":origin:" << vlog[iv].index_ph << std::endl;
+ }
+
+ // collect edges
+ std::vector<NPH_edge> edges = collectEdges(planes, vlog, debugMode);
+
+ // collect new vertices ( edge-plane intersections )
+ collectLines(planes, vtx, vlog, edges, debugMode);
+
+ // ===================== loop over faces ==================================================
+
+ bool facesOK = true; unsigned int fail = 0;
+ for (unsigned int ip = 0; ip< planes.size(); ip++) {
+ NPH_plane plane = planes[ip];
+ // collect shared vertices
+ unsigned int ns = 0; for (auto iv : plane.vertices ) {
+ if (vlog[iv].shared) {
+ bool foundVtx = false; for ( auto ihv : plane.hit_vertices ) if (ihv.index == iv) {foundVtx = true; break;}
+ if (!foundVtx) { NPL_vertex plv = filterLog(vlog[iv],edges, ip, planes);
+ planes[ip].hit_vertices.push_back(plv);
+ }
+ }
+ }
+ if (debugMode) {
+ std::cout << "processing plane:" << ip <<":origin:" << plane.index_ph <<":overlapping planes:";
+ for (auto op : plane.overlapping) std::cout << op <<","; std::cout <<" " << std::endl;
+ }
+ for ( auto ov : plane.overlapping ) { // add shared vertices TO DO : this loop seem obsolete / superflu
+ for ( auto iv : planes[ov%1000].vertices ) { if (!vlog[iv].shared) continue;
+ if (!insideFace( ip, vtx[iv]) ) continue;
+ bool foundVtx = false; for ( auto ihv : plane.hit_vertices ) if (ihv.index == iv) {foundVtx = true; break;}
+ if (!foundVtx) { NPL_vertex plv = filterLog(vlog[iv],edges, ov%1000, planes);
+ planes[ip].hit_vertices.push_back(plv);
+ }
+ }
+ }
+
+ plane = planes[ip];
+ if (debugMode) std::cout << "number of shared vertices:" << plane.hit_vertices.size() << std::endl;
+
+ // =================== all vertices collected ==========================================
+
+ // run first loop over edges to (1) collect endpoint basic info (2) mark "onEdge" ( on "native" edge) vertices
+ // plane edges
+ std::vector<NPL_edge> plane_edges;
+ for (unsigned int ie = 0; ie<plane.edges.size(); ie++) {
+ NPH_edge edge = edges[plane.edges[ie]];
+ int ihit = -1; int ihit_start = -1; double dist = 0.;
+ int vin = edge.endpoints.first; int fin = plane.vlog_ori[ie].faces.first;
+ if (vlog[vin].shared) {
+ for (unsigned int ivl = 0; ivl < planes[ip].hit_vertices.size(); ivl++) {
+ if ( vin == plane.hit_vertices[ivl].index ) { planes[ip].hit_vertices[ivl].onEdge = ie; break; }
+ }
+ }
+ for (unsigned int iv = 0; iv < edge.evtx.size() ; iv++) {
+ // mark vertex on edge
+ int ivx_hit = -1;
+ for (unsigned int ivl = 0; ivl < planes[ip].hit_vertices.size(); ivl++) {
+ if ( edge.evtx[iv].index == planes[ip].hit_vertices[ivl].index ) { planes[ip].hit_vertices[ivl].onEdge = ie; ivx_hit = ivl; break; }
+ }
+ // collect edge piece
+ NPL_edge ple( std::pair<int,int>(vin,edge.evtx[iv].index),edge.faces.second); ple.generating_edge = plane.edges[ie];
+ // check mid-edge before shared classification
+ if (vlog[ple.endpoints.first].shared && vlog[ple.endpoints.second].shared ) {
+ GeoTrf::Vector3D vxmid = 0.5*(vtx[ple.endpoints.first] + vtx[ple.endpoints.second]);
+ if (phA->contains(vxmid,1.e-3) && phB->contains(vxmid,1.e-3) ) ple.shared = true;
+ if (!ple.shared ) {
+ if (debugMode) { std::cout <<"mid-edge point cancels shared edge:"<< ple.endpoints.first<<"," << ple.endpoints.second << std::endl;
+ if (!phA->contains(vxmid,1.e-3)) std::cout << "repeat in debug mode A:" << phA->contains(vxmid,1.e-3,true) << std::endl;
+ if (!phB->contains(vxmid,1.e-3)) std::cout << "repeat in debug mode B:" << phB->contains(vxmid,1.e-3,true) << std::endl;
+ }
+ }
+ }
+ if (ivx_hit>-1) ple.endpoint_hitvx.second = ivx_hit;
+ plane_edges.push_back(ple);
+ // prepare next piece
+ vin = edge.evtx[iv].index;
+ }
+ // collect last
+ if (vin != edge.endpoints.second) {
+ NPL_edge ple( std::pair<int,int>(vin,edge.endpoints.second),edge.faces.second);
+ if (vlog[ple.endpoints.first].shared && vlog[ple.endpoints.second].shared ) {
+ GeoTrf::Vector3D vxmid = 0.5*(vtx[ple.endpoints.first] + vtx[ple.endpoints.second]);
+ if (plane.index_ph == 0 && phB->contains(vxmid,1.e-3) ) ple.shared = true;
+ else if (plane.index_ph == 1 && phA->contains(vxmid,1.e-3) ) ple.shared = true;
+ if (!ple.shared) {
+ if (debugMode) {
+ std::cout <<"mid-edge point cancels last shared edge:"<< ple.endpoints.first<<"," << ple.endpoints.second << std::endl;
+ if (!phA->contains(vxmid)) std::cout << "repeat in debug mode A:" << phA->contains(vxmid,1.e-3,true) << std::endl;
+ if (!phB->contains(vxmid)) std::cout << "repeat in debug mode B:" << phB->contains(vxmid,1.e-3,true) << std::endl;
+ }
+ }
+
+ }
+ if ( vlog[ple.endpoints.second].shared ) {
+ int ivx_hit = -1;
+ for (unsigned int ivl = 0; ivl < plane.hit_vertices.size(); ivl++) {
+ if ( edge.endpoints.second == plane.hit_vertices[ivl].index ) { plane.hit_vertices[ivl].onEdge = ie; ivx_hit = ivl; break; }
+ }
+ ple.endpoint_hitvx.second = ivx_hit;
+ }
+ ple.generating_edge = plane.edges[ie];
+ plane_edges.push_back(ple);
+ }
+ }
+ // loop over edges to identify vertices with protected line connection
+ for (unsigned int ie = 0; ie< plane_edges.size()-1; ie++) {
+ if ( plane_edges[ie].shared != plane_edges[ie+1].shared ) planes[ip].hit_vertices[plane_edges[ie].endpoint_hitvx.second].onEdgeProtected = true;
+ }
+ if ( plane_edges.back().shared != plane_edges.front().shared ) planes[ip].hit_vertices[plane_edges.back().endpoint_hitvx.second].onEdgeProtected = true;
+ // loop over edges to identify dummy overlaps
+ for (unsigned int ie = 0; ie< plane_edges.size()-1; ie++) { NPL_edge ple = plane_edges[ie];
+ for ( auto ov : planes[ple.adjacent_face].overlapping ) { if ( ov < 1000 ) continue;
+ if ( insideFace( ov%1000, vtx[ple.endpoints.first]) && insideFace( ov%1000, vtx[ple.endpoints.second]) ) {
+ plane_edges[ie].dummy_overlap = ov;
+ }
+ }
+ }
+ //
+ // ======================== prepare active lines ( mixed edges ) ====================================
+
+ std::vector<NPL_edge> resolved_lines = evaluateVertexLogic(ip, vlog, edges, planes, vtx);
+ if (debugMode) {
+ std::cout << "resolved lines:" << resolved_lines.size() <<":"; for (auto il : resolved_lines) std::cout << il.adjacent_face <<",";
+ std::cout <<""<<std::endl;
+ }
+ plane = planes[ip];
+
+ // ====================== second loop over edges and lines =========================================
+
+ // resolve endpoint faces
+ unsigned int plane_edge_non_shared=0;
+ for (unsigned int ie = 0; ie<plane_edges.size(); ie++) {
+ NPH_edge edge = edges[plane_edges[ie].generating_edge];
+ int ie_next = ie<plane_edges.size()-1 ? ie+1 : 0;
+ int ie_prev = ie>0 ? ie-1 : plane_edges.size()-1;
+ NPL_edge ple = plane_edges[ie];
+ NPL_edge ple_next = plane_edges[ie_next];
+ NPL_edge ple_prev = plane_edges[ie_prev];
+ plane_edges[ie].endpoint_hitvx.first = ple_prev.endpoint_hitvx.second;
+
+ if (!ple.shared && !ple_next.shared) { // non-shared edges
+ if (!vlog[ple.endpoints.second].shared) {
+ plane_edges[ie].endpoint_faces.second = ple_next.adjacent_face;
+ plane_edges[ie_next].endpoint_faces.first = ple.adjacent_face;
+ } else {
+ if (plane.hit_vertices[ple.endpoint_hitvx.second].active_lines.size() == 2) {
+ plane_edges[ie].endpoint_faces.second = ple_next.adjacent_face;
+ plane.hit_vertices[ple.endpoint_hitvx.second].faces.first = ple.adjacent_face;
+ plane.hit_vertices[ple.endpoint_hitvx.second].faces.second = ple_next.adjacent_face;
+ plane.hit_vertices[ple.endpoint_hitvx.second].faces_shared_loop.first = plane.hit_vertices[ple.endpoint_hitvx.second].active_lines[0];
+ plane.hit_vertices[ple.endpoint_hitvx.second].faces_shared_loop.second = plane.hit_vertices[ple.endpoint_hitvx.second].active_lines[1];
+ } else if (plane.hit_vertices[ple.endpoint_hitvx.second].active_lines.size() == 0) {
+ plane_edges[ie].endpoint_faces.second = ple_next.adjacent_face;
+ plane_edges[ie_next].endpoint_faces.first = ple.adjacent_face;
+ plane.hit_vertices[ple.endpoint_hitvx.second].faces.first = ple.adjacent_face;
+ plane.hit_vertices[ple.endpoint_hitvx.second].faces.second = ple_next.adjacent_face;
+ plane.hit_vertices[ple.endpoint_hitvx.second].faces_shared_loop.first = ple.adjacent_face;
+ plane.hit_vertices[ple.endpoint_hitvx.second].faces_shared_loop.second = ple_next.adjacent_face;
+ plane.hit_vertices[ple.endpoint_hitvx.second].oriented = true;
+ }
+ }
+ } else if (ple.shared) {
+ // single valid line + edge not in plane
+ std::vector<int> pruned_lines; bool edgeInPlane;
+ for ( auto ap : plane.hit_vertices[ple.endpoint_hitvx.second].active_lines ) {
+ edgeInPlane = false; for ( auto ep : edge.inPlane ) if (ap == ep) edgeInPlane = true;
+ if (!edgeInPlane && inPlane(vlog[ple_next.endpoints.second],edges,ap) ) edgeInPlane = true;
+ if (!edgeInPlane) pruned_lines.push_back(ap);
+ if (edgeInPlane && ( plane.hit_vertices[ple.endpoint_hitvx.first].onEdgeProtected || plane.hit_vertices[ple.endpoint_hitvx.second].onEdgeProtected) )
+ std::cout << "pruning removed protected line:" << ap <<":endpoints:" << ple.endpoints.first <<"," << ple.endpoints.second << std::endl;
+ }
+ int nlines = pruned_lines.size();
+ if ( nlines == 1 ) {
+ int al = pruned_lines[0];
+ bool activeLine = false; for ( auto jl : resolved_lines ) if ( jl.adjacent_face == al ) activeLine = true;
+
+ if (activeLine) {
+ plane.hit_vertices[ple.endpoint_hitvx.second].faces.first = ple.adjacent_face ;
+ plane.hit_vertices[ple.endpoint_hitvx.second].faces.second = al;
+ plane.hit_vertices[ple.endpoint_hitvx.second].oriented = true;
+ plane.hit_vertices[ple.endpoint_hitvx.second].faces_shared_loop.first = ple.adjacent_face;
+ plane.hit_vertices[ple.endpoint_hitvx.second].faces_shared_loop.second = al;
+ plane_edges[ie].endpoint_faces.second = al;
+ plane_edges[ie_next].endpoint_faces.first = al;
+ for (unsigned int ia=0; ia<resolved_lines.size(); ia++) {if (resolved_lines[ia].adjacent_face!=al) continue;
+ if (!ple.shared) resolved_lines[ia].endpoint_hitvx.second = ple.endpoint_hitvx.second;
+ else resolved_lines[ia].endpoint_hitvx.first = ple.endpoint_hitvx.second;
+ }
+ }
+ } else if (nlines==0) {
+ plane_edges[ie].endpoint_faces.second = ple_next.adjacent_face;
+ plane_edges[ie_next].endpoint_faces.first = ple.adjacent_face;
+ plane.hit_vertices[ple.endpoint_hitvx.second].faces.first = ple.adjacent_face ;
+ plane.hit_vertices[ple.endpoint_hitvx.second].faces.second = ple_next.adjacent_face;
+ } else { // implies shared vertex adjacent to double line
+ // there will be an ( unoriented back and forth loop ) + two active overlapping lines
+ // initial choice random but needs to be consistent for endpoints
+ NPL_edge l0; for (auto il : resolved_lines) if ( pruned_lines[0] == il.adjacent_face ) l0 = il;
+ NPL_edge l1; for (auto il : resolved_lines) if ( pruned_lines[1] == il.adjacent_face ) l1 = il;
+ if (l0.endpoints.first == l1.endpoints.first && l0.endpoints.second == l1.endpoints.second ) {
+ l1.endpoints.first = l0.endpoints.second; l1.endpoints.second = l0.endpoints.first;
+ if (debugMode) std::cout << "double line, hit vertices?" << l0.endpoint_hitvx.first << "," << l0.endpoint_hitvx.second << ":" <<
+ l1.endpoint_hitvx.first << "," << l1.endpoint_hitvx.second<<std::endl;
+ }
+ // anticipate link to sym lines
+ if ( plane_edges[ie].endpoints.second == l0.endpoints.first ) {
+ plane_edges[ie].endpoint_faces.second = pruned_lines[1];
+ plane_edges[ie_next].endpoint_faces.first = pruned_lines[0];
+ } else {
+ plane_edges[ie].endpoint_faces.second = pruned_lines[0];
+ plane_edges[ie_next].endpoint_faces.first = pruned_lines[1];
+ }
+ plane.hit_vertices[ple.endpoint_hitvx.second].faces.first = ple.adjacent_face ;
+ plane.hit_vertices[ple.endpoint_hitvx.second].faces.second = plane_edges[ie].endpoint_faces.second;
+ plane.hit_vertices[ple_next.endpoint_hitvx.second].faces.second = ple_next.adjacent_face ;
+ plane.hit_vertices[ple_next.endpoint_hitvx.second].faces.first = plane_edges[ie].endpoint_faces.second;
+ }
+ } // end non/shared vertices
+ } // end second loop over edges
+
+ if (debugMode) { // print overview of shared edges
+ for (unsigned int ie = 0; ie<plane_edges.size(); ie++) {
+ if (!plane_edges[ie].shared ) continue;
+ NPL_edge ple = plane_edges[ie];
+ std::cout <<"edge resolved:" << ple.endpoints.first <<":" << ple.endpoints.second<<"("<<ple.endpoint_faces.first <<"," << ple.adjacent_face
+ <<"," << ple.endpoint_faces.second <<")"<<":faces:"
+ << plane.hit_vertices[ple.endpoint_hitvx.first].faces.second<< "," << plane.hit_vertices[ple.endpoint_hitvx.first].faces.first
+ << ":" << plane.hit_vertices[ple.endpoint_hitvx.second].faces.first <<","
+ << plane.hit_vertices[ple.endpoint_hitvx.second].faces.second <<":faces shared loop:"
+ << plane.hit_vertices[ple.endpoint_hitvx.first].faces_shared_loop.second <<","
+ << plane.hit_vertices[ple.endpoint_hitvx.first].faces_shared_loop.first <<":"
+ << plane.hit_vertices[ple.endpoint_hitvx.second].faces_shared_loop.first<<","
+ << plane.hit_vertices[ple.endpoint_hitvx.second].faces_shared_loop.second << std::endl;
+ }
+ }
+ // ============== merge edges and lines ================================================
+ // at this point, lines are not oriented and have not been symmetrized yet
+ // the ordering ensures oriented edges are recovered first
+
+ // lines duplicating shared edges are maintained but pushed to the end of the buffer
+
+ std::vector<NPL_edge> duplicated_edges;
+
+ for ( unsigned int il=0; il<resolved_lines.size(); il++) {
+ // sanity check : copy problems TODO use map index(global)->hit_vx(in plane)
+ if (plane.hit_vertices[resolved_lines[il].endpoint_hitvx.first].index != resolved_lines[il].endpoints.first) {
+ for (unsigned int ih=0; ih<plane.hit_vertices.size(); ih++) {
+ if (plane.hit_vertices[ih].index == resolved_lines[il].endpoints.first ) { resolved_lines[il].endpoint_hitvx.first = ih; break; }
+ }
+ }
+ if (plane.hit_vertices[resolved_lines[il].endpoint_hitvx.second].index != resolved_lines[il].endpoints.second) {
+ for (unsigned int ih=0; ih<plane.hit_vertices.size(); ih++) {
+ if (plane.hit_vertices[ih].index == resolved_lines[il].endpoints.second ) { resolved_lines[il].endpoint_hitvx.second = ih; break; }
+ }
+ }
+
+ // catch duplicated edges
+ bool foundDuplicate = false; int ei = resolved_lines[il].endpoints.first; int eo = resolved_lines[il].endpoints.second;
+ for (unsigned int ie=0; ie<plane_edges.size(); ie++) {
+ if ( ( ei == plane_edges[ie].endpoints.first && eo == plane_edges[ie].endpoints.second )
+ || ( ei == plane_edges[ie].endpoints.second && eo == plane_edges[ie].endpoints.first ) ) {
+ foundDuplicate = true; // plane_edges[ie].shared = false;
+ }
+ if (foundDuplicate) break;
+ }
+
+ if (!foundDuplicate || resolved_lines[il].protectedLine) plane_edges.push_back(resolved_lines[il]);
+ else duplicated_edges.push_back(resolved_lines[il]);
+ }
+ if (debugMode && duplicated_edges.size() ) std::cout << "number of duplicated edges:" << duplicated_edges.size() << std::endl;
+
+ if ( plane.hit_vertices.size()>2 ) {
+ for (auto ld : duplicated_edges) { bool protectedLine = false;
+ if ( plane.hit_vertices[ld.endpoint_hitvx.first].onEdgeProtected && plane.hit_vertices[ld.endpoint_hitvx.first].active_lines.size()==1)
+ protectedLine = true;
+ if ( plane.hit_vertices[ld.endpoint_hitvx.second].onEdgeProtected && plane.hit_vertices[ld.endpoint_hitvx.second].active_lines.size()==1)
+ protectedLine = true;
+ if (protectedLine) {
+ if (debugMode) std::cout << "protected line added:" << ld.endpoints.first <<":" << ld.endpoints.second <<":active lines:"
+ << plane.hit_vertices[ld.endpoint_hitvx.first].active_lines.size()<<","
+ << plane.hit_vertices[ld.endpoint_hitvx.second].active_lines.size() <<std::endl;
+ plane_edges.push_back(ld);
+ }
+ }
+ }
+
+ if (debugMode) {
+ std::cout << "combined edges:" << plane_edges.size() <<":(lines)"<< resolved_lines.size() << std::endl;
+ for ( auto ie : plane_edges) std::cout << ie.endpoints.first <<":" << ie.endpoints.second <<":adj:" << ie.adjacent_face<<":gen.edge:" << ie.generating_edge <<":sym:" << ie.symmetrized <<":shared:" << ie.shared <<":hitvx:" << ie.endpoint_hitvx.first <<"," << ie.endpoint_hitvx.second <<":dummy overlap:" << ie.dummy_overlap << std::endl;
+ }
+
+ // ============== shared loops : overlapping faces =========================================================
+
+ std::vector<NPL_ordered> shared_loops;
+ for ( auto ov : planes[ip].overlapping ) {
+ std::vector<int> ovtx;
+ for (unsigned int iv=0; iv<plane.hit_vertices.size(); iv++ ) {
+ if (debugMode) {
+ if ( !insideFace( ov%1000, vtx[plane.hit_vertices[iv].index]) &&
+ inPlane(vlog[plane.hit_vertices[iv].index],edges,ov%1000) )
+ std::cout << "fake overlap vtx?" << plane.hit_vertices[iv].index << std::endl;
+ }
+
+ if ( insideFace( ov%1000, vtx[plane.hit_vertices[iv].index])
+ || inPlane(vlog[plane.hit_vertices[iv].index],edges,ov%1000)) {
+ if (!ovtx.size() || plane.hit_vertices[iv].onEdge>=0 || plane.hit_vertices[ovtx.front()].oriented) ovtx.push_back(iv);
+ else if ( plane.hit_vertices[iv].oriented ) ovtx.insert(ovtx.begin(),iv);
+ else if ( plane.hit_vertices[iv].onEdge>=0 ) ovtx.insert(ovtx.begin(),iv);
+ else ovtx.push_back(iv);
+ }
+ }
+ if (ovtx.size()<2) continue;
+
+ int nedge = 0; std::vector<int> adjDiff; std::vector<int> edgeSelect;
+ for (unsigned int ie=0; ie<plane_edges.size(); ie++ )
+ if ( (insideFace( ov%1000, vtx[plane_edges[ie].endpoints.first]) || inPlane(vlog[plane_edges[ie].endpoints.first],edges,ov%1000)) &&
+ (insideFace( ov%1000, vtx[plane_edges[ie].endpoints.second]) || inPlane(vlog[plane_edges[ie].endpoints.second],edges,ov%1000) )
+ && !plane_edges[ie].used) {
+ GeoTrf::Vector3D vxmid = 0.5*(vtx[plane_edges[ie].endpoints.first] + vtx[plane_edges[ie].endpoints.second]);
+ bool midIn = insideFace( ov%1000, vxmid);
+ if (midIn) {
+ nedge++; plane_edges[ie].overlap = ov;
+ if (!adjDiff.size()) adjDiff.push_back(plane_edges[ie].adjacent_face);
+ else {
+ bool adjFound = false; for (auto ia : adjDiff) if (ia == plane_edges[ie].adjacent_face ) adjFound= true;
+ if (!adjFound ) adjDiff.push_back(plane_edges[ie].adjacent_face);
+ }
+ edgeSelect.push_back(ie);
+ }
+ }
+
+ if ( nedge<2 ) continue;
+ if ( adjDiff.size()<3 ) continue; // loop needs at least 3 different edges
+
+ NPL_ordered loop; loop.edges.push_back(edgeSelect[0]); plane_edges[edgeSelect[0]].used = true;
+ loop.vtx.push_back(plane_edges[edgeSelect[0]].endpoints.first);
+ loop.vtx.push_back(plane_edges[edgeSelect[0]].endpoints.second);
+ loop.shared=true; loop.top_overlap=ov; loop.oriented = plane_edges[edgeSelect[0]].generating_edge>=0 ? true : false;
+ std::vector<int>::iterator it = edgeSelect.begin();
+ while ( loop.edges.size() < ovtx.size() && loop.vtx.front() != loop.vtx.back() && it != edgeSelect.end()) {
+ if (!plane_edges[*it].used && plane_edges[*it].endpoints.first == loop.vtx.back()) {
+ loop.edges.push_back(*it); loop.vtx.push_back(plane_edges[*it].endpoints.second);
+ plane_edges[*it].used = true; it = edgeSelect.begin();
+ } else if (!plane_edges[*it].used && plane_edges[*it].generating_edge<0 && plane_edges[*it].symmetrized<0
+ && plane_edges[*it].endpoints.second == loop.vtx.back()) { // need to symmetrize the line now
+ NPL_edge symLine = symmetrize(plane_edges[*it]); symLine.used=true; plane_edges[*it].symmetrized = plane_edges.size(); symLine.symmetrized=(*it);
+ (*it) = plane_edges.size(); plane_edges.push_back(symLine);
+ loop.edges.push_back(*it); loop.vtx.push_back(symLine.endpoints.second); it = edgeSelect.begin();
+ } else it++;
+ }
+
+ if (debugMode) {
+ std::cout << "overlapping shared loop done:" << loop.vtx.size() << ":" << loop.edges.size() <<std::endl;
+ for ( auto ie : loop.edges) std::cout << "loop content:" << plane_edges[ie].endpoints.first <<"," << plane_edges[ie].endpoints.second <<":" <<
+ plane_edges[ie].used << std::endl;
+ }
+
+ shared_loops.push_back(loop);
+
+ }
+
+ // =================== shared loops : base ===========================================
+ // select shared non-overlapping first, then add all unused shared
+
+ int nedge = 0; std::vector<int> edgeSelect;
+ for (unsigned int ie=0; ie<plane_edges.size(); ie++ ) {
+ if ( !plane_edges[ie].used && plane_edges[ie].shared && plane_edges[ie].overlap<0 ) {
+ edgeSelect.push_back(ie); nedge++;
+ }
+ }
+
+ if (nedge>1) {
+ // add overlaps: symmetrize if not done already
+ for (unsigned int ie=0; ie<plane_edges.size(); ie++ ) {
+ if (plane_edges[ie].shared && plane_edges[ie].generating_edge<0 && plane_edges[ie].overlap>=0 && plane_edges[ie].symmetrized<0 ) {
+ NPL_edge symLine = symmetrize(plane_edges[ie]); plane_edges[ie].symmetrized=plane_edges.size(); symLine.symmetrized = ie;
+ plane_edges.push_back(symLine);
+ }
+ }
+ for (unsigned int ie=0; ie<plane_edges.size(); ie++ ) {
+ if (!plane_edges[ie].used && plane_edges[ie].shared && plane_edges[ie].overlap>=0 ) {
+ edgeSelect.push_back(ie); nedge++;
+ }
+ }
+
+ if (edgeSelect.size()>1) {
+ std::vector<int>::iterator it = edgeSelect.begin();
+ NPL_ordered loop; loop.edges.push_back(*it); plane_edges[*it].used = true;
+ loop.vtx.push_back(plane_edges[*it].endpoints.first);
+ loop.vtx.push_back(plane_edges[*it].endpoints.second);
+ loop.shared=true; loop.oriented = plane_edges[*it].generating_edge>=0 ? true : false; // after including check on dummy shared edges, oriented edges may appear in the middle of the loop
+ it = edgeSelect.begin();
+ while ( loop.vtx.front() != loop.vtx.back() && it != edgeSelect.end()) {
+ if ( loop.edges.size() && plane_edges[loop.edges.back()].symmetrized == (*it) ) {
+ it++;
+ } else if (!plane_edges[*it].used && plane_edges[*it].endpoints.first == loop.vtx.back() ) {
+ //
+ loop.edges.push_back(*it); loop.vtx.push_back(plane_edges[*it].endpoints.second);
+ plane_edges[*it].used = true; it = edgeSelect.begin(); // std::cout << ":adding:" << (*it) << std::endl;
+
+ } else if (!plane_edges[*it].used && plane_edges[*it].generating_edge<0 && plane_edges[*it].symmetrized<0
+ && plane_edges[*it].endpoints.second == loop.vtx.back()) { // need to symmetrize the line now
+ NPL_edge symLine = symmetrize(plane_edges[*it]); plane_edges[*it].symmetrized=plane_edges.size(); symLine.used=true; symLine.symmetrized = (*it);
+ (*it) = plane_edges.size(); plane_edges.push_back(symLine); // std::cout << ":adding symmetrized:" << (*it) << std::endl;
+ loop.edges.push_back(*it); loop.vtx.push_back(symLine.endpoints.second); it = edgeSelect.begin();
+
+ } else it++;
+ }
+
+ if (loop.vtx.front() != loop.vtx.back() || loop.edges.size()<2 ) {
+ for ( auto ie : loop.edges ) plane_edges[ie].used = false; // undo loop if not closed
+ } else
+ shared_loops.push_back(loop);
+
+ int nseed = 0; // keep trace of non-symmetrized unused shared edges
+ for (auto ie : edgeSelect ) if (!plane_edges[ie].used) {
+ if (plane_edges[ie].symmetrized<0 ) nseed++;
+ else if (!plane_edges[plane_edges[ie].symmetrized].used) nseed++;
+ }
+
+ if (debugMode) {
+ std::cout << "base shared loop:selected:" << edgeSelect.size() <<":loop size:vtx:"<< loop.vtx.size()
+ <<":loop size:edge:" << loop.edges.size() << std::endl;
+ if (nseed>0) {
+ std::cout <<"WARNING:remaining seeds:" << nseed << std::endl;
+ for (auto iv : planes[ip].vertices) std::cout << "vtx ori:" <<iv<<":" << vtx[iv].x() <<"," << vtx[iv].y() <<"," << vtx[iv].z() << std::endl;
+ for (auto iv : planes[ip].hit_vertices) std::cout << "vtx hit:" <<iv.index<<":" << vtx[iv.index].x() <<"," << vtx[iv.index].y() <<"," << vtx[iv.index].z() << std::endl;
+ }
+ }
+ }
+ }
+
+ //================ symmetrized remaining active lines =====================================
+
+ for (unsigned int ie=0; ie<plane_edges.size(); ie++ ) {
+ if (plane_edges[ie].shared && plane_edges[ie].generating_edge<0 && plane_edges[ie].symmetrized < 0 ) {
+ NPL_edge symLine = symmetrize(plane_edges[ie]); plane_edges[ie].symmetrized=plane_edges.size(); symLine.symmetrized = ie;
+ plane_edges.push_back(symLine);
+ }
+ }
+
+ //============== symmetrized inverted non-oriented shared loops ( -> non-shared ) ==========
+ //============== orientation resolved with help of area calculation ========================
+
+ for (unsigned int il = 0; il< shared_loops.size(); il++) shared_loops[il].vtx.pop_back();
+ std::vector< NPL_ordered> shared_inverted_subface;
+
+ for (auto il : shared_loops) {
+ if ( il.edges.size()>2 ) {
+ if (il.oriented) planes[ip].subfaces.push_back(il);
+ else { // prepare inverted loop
+ NPL_ordered loopi; loopi.vtx.insert(loopi.vtx.begin(),il.vtx.rbegin(),il.vtx.rend());
+ for (unsigned int iv = 0; iv<loopi.vtx.size(); iv++) {
+ int ivn = iv<loopi.vtx.size()-1 ? iv+1 : 0;
+ for (unsigned int ie=0; ie<plane_edges.size(); ie++ ) { if (plane_edges[ie].used) continue;
+ if (plane_edges[ie].endpoints.first == loopi.vtx[iv] && plane_edges[ie].endpoints.second == loopi.vtx[ivn]) {
+ loopi.edges.push_back(ie); plane_edges[ie].used = true; break;
+ }
+ }
+ }
+ loopi.top_overlap = il.top_overlap;
+ loopi.shared = il.shared;
+ // resolve orientation
+ if ( areaVec(il.vtx,vtx).dot(plane.normal)>=0 ) {
+ il.oriented = true; loopi.oriented = true;
+ planes[ip].subfaces.push_back(il);
+ shared_inverted_subface.push_back(loopi);
+ } else {
+ il.oriented = true; loopi.oriented = true;
+ planes[ip].subfaces.push_back(loopi);
+ shared_inverted_subface.push_back(il);
+ }
+ } // end orientation
+ } // valid shared loops
+ } // shared loops
+
+
+ //=============== non-shared loop =========================================================
+ // vertices may be unshared ( no endpoint_hitvx info )
+ std::vector<NPL_ordered> non_shared_loops;
+
+ nedge = 0; std::vector<int> edgeNS; int istart =-1;
+ for (unsigned int ie=0; ie<plane_edges.size(); ie++ ) {
+ if ( !plane_edges[ie].used ) {
+ edgeNS.push_back(ie); nedge++;
+ if (istart<0) istart = edgeNS.size()-1;
+ }
+ }
+
+ if (debugMode) {
+ std::cout <<"non-shared loop(s):" <<edgeNS.size() <<":" << istart << std::endl;
+ for ( auto ie : edgeNS) std::cout << plane_edges[ie].endpoints.first <<":" << plane_edges[ie].endpoints.second <<":adj:" << plane_edges[ie].adjacent_face<<":gen.edge:" << plane_edges[ie].generating_edge <<":sym:" << plane_edges[ie].symmetrized << std::endl;
+ }
+
+ int istart_last=istart; int nfail = 0;
+
+ while ( nedge > 1 && istart>=0 ) {
+
+ NPL_ordered loop; loop.edges.push_back(edgeNS[istart]); plane_edges[edgeNS[istart]].used = true;
+ loop.vtx.push_back(plane_edges[edgeNS[istart]].endpoints.first);
+ loop.vtx.push_back(plane_edges[edgeNS[istart]].endpoints.second);
+ loop.oriented = plane_edges[edgeNS[istart]].generating_edge>=0 ? true : false;
+ std::vector<int>::iterator it = edgeNS.begin();
+ while ( loop.vtx.front() != loop.vtx.back() && it != edgeNS.end()) {
+ if ( loop.edges.size() && plane_edges[loop.edges.back()].symmetrized == (*it) ) {
+ it++;
+ } else if (!plane_edges[*it].used && plane_edges[*it].endpoints.first == loop.vtx.back() ) {
+ loop.edges.push_back(*it); loop.vtx.push_back(plane_edges[*it].endpoints.second);
+ plane_edges[*it].used = true; it = edgeNS.begin();
+ } else it++;
+ }
+
+ istart_last = istart;
+ istart = -1; for (unsigned int ie=0; ie<edgeNS.size(); ie++ ) if (!plane_edges[edgeNS[ie]].used) { istart=ie; break; }
+ if (debugMode) std::cout << "non-shared loop done:" << loop.edges.size() << ":start next:" << istart << std::endl;
+
+ if (loop.vtx.front() != loop.vtx.back() ) {
+ // find alternative loop before releasing loop
+ if (debugMode) std::cout << "non-shared loop not closed, undo" << std::endl;
+ for ( auto ie : loop.edges ) plane_edges[ie].used = false;
+ if (debugMode) {
+ for ( auto ie : loop.edges) std::cout <<"failed loop:" << plane_edges[ie].endpoints.first <<":" << plane_edges[ie].endpoints.second <<":adj:" << plane_edges[ie].adjacent_face<<":gen.edge:" << plane_edges[ie].generating_edge <<":sym:" << plane_edges[ie].symmetrized << std::endl;
+ }
+ nfail++; if (istart == istart_last || nfail>10) break;
+ } else {
+ // shared loop left-over ?
+ bool nshared = false; for (auto iv : loop.vtx ) if (!vlog[iv].shared) nshared = true;
+ if (nshared) non_shared_loops.push_back(loop);
+ else { loop.shared = true; loop.vtx.pop_back(); shared_loops.push_back(loop); }
+ nedge -= loop.edges.size();
+ }
+ }
+
+ if (debugMode) {
+ std::cout << "non-shared loops:" << non_shared_loops.size() <<":remaining:"<< nedge << std::endl;
+ if (nedge>0 && plane.hit_vertices.size() > 2 ) {
+ std::cout <<"WARNING:edges left over:" << nedge <<":hit vertices:" << plane.hit_vertices.size() << std::endl;
+ for (auto ie : edgeNS) if (!plane_edges[ie].used) std::cout <<"not used:" << plane_edges[ie].endpoints.first <<","
+ << plane_edges[ie].endpoints.second <<":shared:"
+ << plane_edges[ie].shared <<":gen.edge:"
+ << plane_edges[ie].generating_edge <<":adj:"
+ <<plane_edges[ie].adjacent_face << std::endl;
+ }
+ }
+
+ // =============================== loop overview ===============================================
+
+ for (unsigned int il = 0; il< non_shared_loops.size(); il++) non_shared_loops[il].vtx.pop_back();
+
+ //std::cout << "original face:" <<ip<<":areaVec.dot(normal):" << areaVec(plane.vertices,vtx).dot(plane.normal) << std::endl;
+ if (debugMode) {
+ std::cout <<"loop overview for face:" << ip << ":original vertices:"<< plane.vertices.size() << std::endl;
+ for (auto iv : plane.vertices) std::cout << iv <<"," ; std::cout <<"" << std::endl;
+ std::cout <<"loops:shared:" << shared_loops.size() <<":inverted:" << shared_inverted_subface.size() <<":non-shared:" << non_shared_loops.size() << std::endl;
+ }
+
+ for (auto il : plane.subfaces) {
+ if (debugMode) {
+ std::cout <<"to face:" << ip <<":shared:" << il.shared <<":oriented:" << il.oriented <<":size:" << il.vtx.size() <<"," << il.edges.size() << std::endl;
+ for (auto iv : il.vtx) std::cout << iv <<"," ; std::cout <<"" << std::endl;
+ }
+ }
+
+ /* this part removes the self-eliminating parts of loop, simplifying the description
+ for (auto il : shared_loops) {
+ if (debugMode) {
+ std::cout <<"to face:" << ip <<":shared:" << il.shared <<":oriented:" << il.oriented <<":overlap:" << il.top_overlap << ":size:" << il.vtx.size() <<"," << il.edges.size() << std::endl;
+ std::cout << "shared face area:" << areaVec(il.vtx,vtx).dot(plane.normal) <<":oriented:" << il.oriented << std::endl;
+ for (auto iv : il.vtx) std::cout << iv <<"," ; std::cout <<"" << std::endl;
+ for (auto ie : il.edges) std::cout << plane_edges[ie].dummy_overlap <<"," ; std::cout <<"" << std::endl;
+ }
+ if (il.vtx.size()>2) {
+ std::vector<size_t> vx_lobe; std::vector<int> ee_lobe; bool lobe = false; unsigned int vmax = il.vtx.size()-1;
+ for (unsigned int iv=0; iv< il.vtx.size(); iv++) {
+ if (iv>vmax) {
+ vx_lobe.push_back(il.vtx[iv]);
+ ee_lobe.push_back(il.edges[iv]);
+ continue;
+ }
+ unsigned int ivm2 = iv<il.vtx.size()-2 ? iv+2 : iv+2-il.vtx.size();
+ unsigned int ivm1 = iv<il.vtx.size()-1 ? iv+1 : iv+1-il.vtx.size();
+ unsigned int ivn1 = iv>0 ? iv-1 : il.vtx.size()-1;
+ if ( il.vtx[iv] == il.vtx[ivm2] && ( plane_edges[il.edges[iv]].dummy_overlap>=1000 || plane_edges[il.edges[ivm1]].dummy_overlap>=1000) ) {
+ if (debugMode) std::cout << "lobe found at:" << iv <<":" << il.vtx[iv] << std::endl;
+ lobe = true; iv++;
+ } else if (il.vtx[ivm1] == il.vtx[ivn1] && ( plane_edges[il.edges[ivn1]].dummy_overlap>=1000 || plane_edges[il.edges[iv]].dummy_overlap>=1000) ) {
+ lobe = true; iv++; vmax = il.vtx.size()-2;
+ } else {
+ vx_lobe.push_back(il.vtx[iv]);
+ ee_lobe.push_back(il.edges[iv]);
+ }
+ }
+ if (lobe) { il.vtx = vx_lobe; il.edges = ee_lobe; }
+ }
+ planes[ip].shared_subface.push_back(il);
+ }
+ */
+
+ if ( non_shared_loops.size() ) {
+ // merge inverted shared with non-shared loops
+ for (unsigned int il=0; il<shared_inverted_subface.size(); il++ ) {
+ non_shared_loops[0] = mergeLoops(non_shared_loops[0], shared_inverted_subface[il], vtx);
+ }
+ }
+
+ for (auto il : non_shared_loops) {
+ if (debugMode) {
+ std::cout <<"to face:" << ip <<":shared:" << il.shared <<":oriented:" << il.oriented <<":size:" << il.vtx.size() <<"," << il.edges.size() << std::endl;
+ for (auto iv : il.vtx) std::cout << iv <<"," ; std::cout <<"" << std::endl;
+ }
+ planes[ip].subfaces.push_back(il);
+ }
+
+ // save edges for further processing
+ planes[ip].plane_edges = plane_edges;
+
+ /*
+ // sanity check : matching areas
+ double sa = 0.;
+ for ( auto il : planes[ip].subfaces ) {
+ sa += areaVec(il.vtx,vtx).dot(planes[ip].normal);
+ }
+ double aori = areaVec(planes[ip].vertices,vtx).dot(planes[ip].normal);
+ if ( std::abs(aori - sa) > 1 ) std::cout << "sanity check plane:area:"<< ip << ":" << sa << ":original:" << aori << std::endl;
+ */
+
+ } // end loop over planes
+
+ GeoPolyhedronGeneric* gpg = nullptr;
+
+ // construct intersection
+ GeoPolyhedronGeneric intersection;
+ intersection.addVertices(vtx);
+ for ( unsigned int ip=0; ip<planes.size(); ip++ ) { NPH_plane plane = planes[ip];
+ for ( unsigned int il=0; il<plane.subfaces.size(); il++) {
+ if ( !plane.subfaces[il].shared ) continue;
+ NPL_ordered loop = plane.subfaces[il];
+ if ( loop.top_overlap<0 || (plane.index_ph==0 && loop.top_overlap<1000) ) {
+ intersection.addFace(loop.vtx);
+ intersection.addNormal( plane.index_ph==0 ? phA->faceNormal(ip) : phB->faceNormal(ip-nPA) );
+ }
+ }
+ }
+ if (debugMode && intersection.getNFaces()>0 && !intersection.construct(true)) {
+ std::cout <<"VERIFY intersection:GPG origin:" << phA->getOrigin() <<":"<<phB->getOrigin() << std::endl;
+ for (unsigned int ii=0; ii<intersection.getNFaces(); ii++) {
+ std::vector<size_t> vtx = intersection.getFace(ii);
+ std::cout <<"face:"<<ii<<":"; for (auto vx : vtx) std::cout <<vx <<","; std::cout <<"" << std::endl;
+ }
+ for (unsigned int ip=0; ip<planes.size(); ip++) { NPH_plane plane = planes[ip];
+ int nShared = 0; for ( auto il : plane.subfaces ) if (il.shared) nShared++;
+ if ( plane.hit_vertices.size()>2 && nShared<1 ) std::cout <<"PROBLEM in plane:" << ip << std::endl;
+ }
+ }
+
+ if ( modeBoolean == 0 && intersection.getNFaces() > 3 ) gpg = new GeoPolyhedronGeneric(intersection);
+
+ // construct subtraction
+ GeoPolyhedronGeneric subtraction;
+ subtraction.addVertices(vtx);
+ for ( unsigned int ip=0; ip<planes.size(); ip++ ) { NPH_plane plane = planes[ip];
+ if ( plane.index_ph == 0 ) { // collect non-shared && inverted
+ for ( unsigned int il=0; il<plane.subfaces.size(); il++) { if (plane.subfaces[il].shared) continue;
+ NPL_ordered loop = plane.subfaces[il];
+ subtraction.addFace(loop.vtx);
+ subtraction.addNormal( phA->faceNormal(ip));
+ }
+ for ( unsigned int il=0; il<plane.subfaces.size(); il++) { NPL_ordered loop = plane.subfaces[il];
+ if ( !loop.shared ) continue;
+ if ( loop.top_overlap<1000) continue;
+ subtraction.addFace(loop.vtx);
+ subtraction.addNormal( phA->faceNormal(ip));
+ }
+ } else { // invert shared non-overlapping
+ for ( unsigned int il=0; il<plane.subfaces.size(); il++) {
+ if (!plane.subfaces[il].shared) continue;
+ NPL_ordered loop = plane.subfaces[il];
+ if ( loop.top_overlap >= 0 ) continue;
+ std::vector<size_t> lvtx; for (auto iv : loop.vtx) lvtx.insert(lvtx.begin(),iv);
+ subtraction.addFace(lvtx);
+ subtraction.addNormal( (-1.)* phB->faceNormal(ip-nPA) );
+ }
+ }
+ }
+ if (debugMode && subtraction.getNFaces()>0 && !subtraction.construct(true)) {
+ std::cout <<"VERIFY subtraction:GPG origin:" << phA->getOrigin() <<":"<<phB->getOrigin() << std::endl;
+ for (unsigned int ii=0; ii<subtraction.getNFaces(); ii++) {
+ std::vector<size_t> vtx = subtraction.getFace(ii);
+ std::cout <<"face:"<<ii<<":"; for (auto vx : vtx) std::cout <<vx <<","; std::cout <<"" << std::endl;
+ }
+ }
+
+ if ( modeBoolean == 1 && subtraction.getNFaces()>3 ) gpg = new GeoPolyhedronGeneric(subtraction);
+
+ // construct union : non-shared + half of overlaps
+ GeoPolyhedronGeneric uni;
+ uni.addVertices(vtx);
+ for ( unsigned int ip=0; ip<planes.size(); ip++ ) { NPH_plane plane = planes[ip];
+ for ( unsigned int il=0; il<plane.subfaces.size(); il++) { if (plane.subfaces[il].shared) continue;
+ NPL_ordered loop = plane.subfaces[il];
+ if (loop.shared) continue;
+ uni.addFace(loop.vtx);
+ uni.addNormal( plane.index_ph==0 ? phA->faceNormal(ip) : phB->faceNormal(ip-nPA) );
+ }
+
+ if ( plane.index_ph == 0 ) {
+ for ( unsigned int il=0; il<plane.subfaces.size(); il++) { NPL_ordered loop = plane.subfaces[il];
+ if ( loop.shared && loop.top_overlap>=0 && loop.top_overlap<1000 ) {
+ uni.addFace(loop.vtx);
+ uni.addNormal( plane.index_ph==0 ? phA->faceNormal(ip) : phB->faceNormal(ip-nPA) );
+ }
+ }
+ }
+ }
+ if (debugMode && !uni.construct(true)) {
+ std::cout <<"VERIFY union:GPG origin:" << phA->getOrigin() <<":"<<phB->getOrigin()<<":valid? " << uni.isValid() << std::endl;
+ std::cout << "review input polyhedrons:A:vtx:faces:" << phA->getNVertices() <<":" << phA->getNFaces() << std::endl;
+ for (unsigned int iv=0 ; iv<phA->getNVertices(); iv++) { GeoTrf::Vector3D xx = phA->getVertex(iv);
+ std::cout << iv <<":" << xx.x() << "," << xx.y() <<"," << xx.z() << std::endl;
+ }
+ for (unsigned int ip=0; ip<phA->getNFaces(); ip++) {
+ std::vector<size_t> vtx = phA->getFace(ip);
+ GeoTrf::Vector3D nA = phA->faceNormal(ip);
+ std::cout << "face A:" << ip <<":normal:" << nA.x() <<"," << nA.y() <<"," << nA.z() << std::endl;
+ for ( auto iv : vtx ) std::cout << iv<<","; std::cout<<""<<std::endl;
+ }
+ std::cout << "review input polyhedrons:B:vtx:faces:" << phB->getNVertices() <<":" << phB->getNFaces() << std::endl;
+ for (unsigned int iv=0 ; iv<phB->getNVertices(); iv++) { GeoTrf::Vector3D xx = phB->getVertex(iv);
+ std::cout << iv <<":" << xx.x() << "," << xx.y() <<"," << xx.z() << std::endl;
+ }
+ for (unsigned int ip=0; ip<phB->getNFaces(); ip++) {
+ std::vector<size_t> vtx = phB->getFace(ip);
+ GeoTrf::Vector3D nB = phB->faceNormal(ip);
+ std::cout << "face B:" << ip <<":normal:" << nB.x() <<"," << nB.y() <<"," << nB.z() << std::endl;
+ for ( auto iv : vtx ) std::cout << iv<<","; std::cout<<""<<std::endl;
+ }
+ //
+ for (unsigned int ii=0; ii<uni.getNFaces(); ii++) {
+ std::vector<size_t> vtx = uni.getFace(ii);
+ std::cout <<"face:"<<ii<<":"; for (auto vx : vtx) std::cout <<vx <<","; std::cout <<"" << std::endl;
+ }
+ }
+
+ if ( modeBoolean == 2 ) gpg = new GeoPolyhedronGeneric(uni);
+
+ if (debugMode) {
+ if (uni.isValid() && intersection.isValid()) std::cout <<"VOLUME check uni:" << phA->volume()+phB->volume()-intersection.volume()
+ <<"=?" << uni.volume() << std::endl;
+ if (subtraction.isValid() && intersection.isValid()) std::cout <<"VOLUME check subtr:" << phA->volume() - intersection.volume()
+ <<"=?" << subtraction.volume() << std::endl;
+ }
+
+ double v0 = (gpg && gpg->isValid()) ? gpg->volume() : 0.;
+
+ if (debugMode && gpg) std::cout << "test gpg before cleanup:" << gpg->isValid() <<":volume:" << v0 << std::endl;
+
+ if (gpg && gpg->isValid()) gpg->cleanAfterBooleanOperation();
+
+ if (debugMode) {
+ if (gpg) { std::cout << "gpg conversion test:valid:" << gpg->isValid();
+ if (gpg->isValid()) std::cout <<":volume:" << gpg->volume() <<":v0:"<<v0;
+ std::cout<<""<< std::endl;
+ }
+ }
+
+ return gpg;
+}
+
+std::vector<NPH_edge> GeoPolyhedronManipulator::collectEdges(std::vector<NPH_plane>& planes, std::vector<NPH_vertex>& vlog, bool debugMode ) const {
+
+ std::vector<NPH_edge> edges;
+ // each edge belongs to 2 planes, attach defining vertices
+ unsigned int nf = planes.size();
+ for (unsigned int index = 0; index<2; index++) {
+ for (unsigned int iplane = 0; iplane < planes.size(); iplane++) {
+ NPH_plane plane = planes[iplane];
+ if ( plane.index_ph != index ) continue;
+ std::vector<size_t> vs = plane.vertices;
+ for (unsigned int iv = 0; iv < vs.size(); iv++) {
+ int evi = vs[iv]; int evf = iv < vs.size()-1 ? vs[iv+1] : vs[0];
+ NPH_edge new_edge(std::pair<int,int>(evi,evf));
+ new_edge.faces.first = iplane; new_edge.index_ph = index;
+ // find the associated plane - full loop needed before put in use
+ for (unsigned int ie=0; ie<edges.size(); ie++) {
+ NPH_edge edge = edges[ie]; if (edge.index_ph != index ) continue;
+ if (edge.endpoints.first == new_edge.endpoints.second && edge.endpoints.second == new_edge.endpoints.first )
+ { edges[ie].faces.second = iplane; new_edge.faces.second = edge.faces.first; new_edge.edge_inverse=ie; edges[ie].edge_inverse=edges.size(); break; }
+ }
+ edges.push_back(new_edge);
+ vlog[evf].edges.push_back(std::make_pair<int>(edges.size()-1,-1));
+ planes[iplane].edges.push_back(edges.size()-1);
+ }
+ }
+ // second loop to add corresponding faces
+ for (unsigned int iplane = 0; iplane < nf; iplane++) {
+ NPH_plane plane = planes[iplane];
+ if ( plane.index_ph != index ) continue;
+ std::vector<size_t> vs = plane.vertices;
+ for (unsigned int iv = 0; iv < vs.size(); iv++) {
+ int evi = vs[iv]; int evf = iv < vs.size()-1 ? vs[iv+1] : vs[0];
+ // search edge with first face iplane and starting vertex at vs[iv+1]
+ for (unsigned int ie = 0; ie<edges.size(); ie++) {
+ if (edges[ie].index_ph != index) continue;
+ if (edges[ie].faces.first==iplane && edges[ie].endpoints.first==evf) {
+ for (unsigned int ii=0; ii<vlog[evf].edges.size(); ii++) {
+ if (edges[vlog[evf].edges[ii].first].faces.first==iplane
+ && edges[vlog[evf].edges[ii].first].endpoints.second==evf) { vlog[evf].edges[ii].second = edges[ie].faces.second; break; }
+ }
+ }
+ } // end search edge
+ } // end loop vertices
+ } // end loop planes
+
+ // third loop to retrieve vertex logic
+ for (unsigned int iplane = 0; iplane < nf; iplane++) {
+ NPH_plane plane = planes[iplane];
+ if ( plane.index_ph != index ) continue;
+ std::vector<size_t> vs = plane.vertices;
+ for (unsigned int iv = 0; iv < vs.size(); iv++) {
+ // collect vertex ordering logic
+ NPL_vertex plv = filterLog(vlog[vs[iv]],edges,iplane,planes);
+ planes[iplane].vlog_ori.push_back(plv);
+ } // end loop vertices
+ } // end loop planes
+
+ } // end loop ph
+
+ if (debugMode) {
+ for (unsigned int iv=0; iv<vlog.size(); iv++) {
+ std::cout << "vertex logic:" << iv<< std::endl;
+ for (unsigned int ie = 0; ie<vlog[iv].edges.size(); ie++) std::cout<< ie <<":" << edges[vlog[iv].edges[ie].first].faces.first<<"," <<
+ edges[vlog[iv].edges[ie].first].faces.second<< ":" << vlog[iv].edges[ie].second << std::endl;
+ }
+ }
+
+ return edges;
+}
+
+void GeoPolyhedronManipulator::collectLines(std::vector<NPH_plane>& planes, std::vector<GeoTrf::Vector3D>& vtx, std::vector<NPH_vertex>& vlog,
+ std::vector<NPH_edge>& edges, bool debugMode ) const {
+ GeoTrf::Vector3D eInit; GeoTrf::Vector3D edge; GeoTrf::Vector3D point; GeoTrf::Vector3D normal;
+ for (unsigned int ie = 0; ie < edges.size(); ie++ ) {
+ eInit = vtx[edges[ie].endpoints.first];
+ edge = vtx[edges[ie].endpoints.second] - eInit;
+ if (debugMode) std::cout << "collectLines:edge:" << ie <<":" << edges[ie].endpoints.first <<"," << edges[ie].endpoints.second
+ <<":ori:" <<edges[ie].index_ph<<":faces:" << edges[ie].faces.first <<":" << edges[ie].faces.second << std::endl;
+ for (unsigned int ip = 0; ip < planes.size() ; ip++) {
+ if (edges[ie].index_ph == planes[ip].index_ph ) continue;
+ std::vector<size_t> pvx = planes[ip].vertices;
+ for (auto vx : pvx) if ( std::abs(edge.cross(vtx[vx] - eInit).norm()) > m_tolerance ) {
+ point = vtx[vx]; break; //if (debugMode) std::cout << "setting ref point:" << point.x()<<"," <<point.y()<<"," << point.z() << std::endl; break;
+ }
+ normal = planes[ip].normal;
+ if (std::abs(edge.dot(normal)) > m_tolerance) {
+ // intersection of the edge line with plane
+ double d = edgeFaceIntersection(eInit, edge, point, normal);
+ if (d*edge.norm()<-m_tolerance || (d-1.)*edge.norm()>m_tolerance) continue;
+ GeoTrf::Vector3D xyz = eInit + edge * d;
+ if ( insideFace(ip,xyz,1.e-5) ) {
+ //if (debugMode) std::cout <<"intersecting plane:" << ip <<":ori:" << planes[ip].index_ph <<":d:" << d
+ // <<":intersection:" <<xyz.x()<<","<<xyz.y()<<","<<xyz.z() <<std::endl;
+ int index = -1; // valid vertex : already found ?
+ double closest = 1.e10;
+ if (std::abs(d*edge.norm())<m_tolerance) { index = edges[ie].endpoints.first;
+ if (!vlog[index].shared) std::cout << "ERROR:i/ endpoints crossing found for unshared vertex:" << index << std::endl;
+ }
+ else if (std::abs((d-1.)*edge.norm())<m_tolerance) { index = edges[ie].endpoints.second;
+ if (!vlog[index].shared) std::cout << "ERROR:ii/ endpoints crossing found for unshared vertex:" << index << std::endl;
+ }
+ else {
+ for (unsigned int iv = 0; iv<vtx.size(); iv++ ) {
+ if ( (vtx[iv] - xyz).norm() < m_tolerance ) { index = iv; break; }
+ if ( (vtx[iv] - xyz).norm() < closest ) { closest = (vtx[iv] - xyz).norm(); }
+ }
+ if (index<0) { // new vertex
+ vtx.push_back(xyz); index = vtx.size()-1;
+ NPH_vertex vc; vc.shared = true;
+ vlog.push_back(vc);
+ }
+ }
+ if (debugMode) std::cout <<"vertex:" << index << ":merge: valid intersection with plane:" << ip << ":normal:"<<normal.x()<<","<<normal.y()<<","<<normal.z()<<"),"<< d << std::endl;
+ bool foundOnEdge = false;
+ if (!foundOnEdge) for (auto iv : edges[ie].evtx) if ( iv.index == index ) { foundOnEdge = true; break; }
+ if (!foundOnEdge) {
+ NPH_vtx nvtx(index); nvtx.edge=ie ; nvtx.hitPlane = ip; nvtx.distance = d; nvtx.faces = std::pair<int,int>(edges[ie].faces.second,ip);
+ std::vector<NPH_vtx>::iterator vIt = edges[ie].evtx.begin();
+ while ( vIt< edges[ie].evtx.end()) {
+ if ( (vtx[(*vIt).index]-eInit).dot(edge)/pow(edge.norm(),2) < d) vIt++;
+ else break;
+ }
+ edges[ie].evtx.insert(vIt,nvtx);
+ if (debugMode) std::cout <<"0:inserting evtx:"<< nvtx.index <<":at:" << (vtx[nvtx.index]-eInit).dot(edge)/pow(edge.norm(),2)
+ <<":endpoints:" << edges[ie].endpoints.first << "," << edges[ie].endpoints.second<<":hit plane:" << nvtx.hitPlane
+ <<":edge faces:" << edges[ie].faces.first << "," << edges[ie].faces.second <<":closest:" << closest << std::endl;
+ }
+ // add the hit vertex to the list of new vertices on plane
+ bool foundOnFace = false; for (auto iv : planes[ip].hit_vertices ) if ( iv.index == index ) { foundOnFace = true; break; }
+ if (!foundOnFace) {
+ NPL_vertex plv; plv.index = index; plv.shared = true; plv.faces = edges[ie].faces;
+ planes[ip].hit_vertices.push_back(plv);
+ }
+ bool foundOnFace1 = false; for (auto iv : planes[edges[ie].faces.first].hit_vertices ) if ( iv.index == index ) { foundOnFace1 = true; break; }
+ if (!foundOnFace1) {
+ NPL_vertex plv1; plv1.index = index; plv1.shared = true; plv1.faces.first = edges[ie].faces.second; plv1.faces.second = ip;
+ planes[edges[ie].faces.first].hit_vertices.push_back(plv1);
+ }
+ bool foundOnFace2 = false; for (auto iv : planes[edges[ie].faces.second].hit_vertices ) if ( iv.index == index ) { foundOnFace2 = true; break; }
+ if (!foundOnFace2) {
+ NPL_vertex plv2; plv2.index = index; plv2.shared = true; plv2.faces.first = ip; plv2.faces.second = edges[ie].faces.first;
+ planes[edges[ie].faces.second].hit_vertices.push_back(plv2);
+ }
+ //
+ vlog[index].edges.push_back(std::make_pair<int>(ie,ip));
+ } // valid edge intersection with face
+ } else {
+ // special configuration: edge laying on face plane, need to calculate edge intersections
+ if ( std::abs((eInit-point).dot(normal))< m_tolerance) {
+ edges[ie].inPlane.push_back(ip);
+ intersectFace(ie, edges, ip, planes, vtx, vlog);
+ }
+ }
+ } // end loop over faces
+ } // end loop over edges
+
+ // loop over intersection vertices to resolve vertex logic
+ for (unsigned int ie = 0; ie < edges.size(); ie++ ) {
+ for (unsigned int iv = 0; iv < edges[ie].evtx.size(); iv++ ) { NPH_vtx vx = edges[ie].evtx[iv];
+ if ( vlog[vx.index].index_ph >= 0 && vlog[vx.index].index_ph != edges[ie].index_ph ) {
+ vx.hitVertex = vx.index;
+ } else if (vlog[vx.index].index_ph<0) {
+ for (auto vl : vlog[vx.index].edges) {
+ if ( vl.first != ie && vl.first != edges[ie].edge_inverse ) {
+ if ( vx.hitEdge < 0 ) { vx.hitEdge = vl.first;
+ vlog[vx.index].hit_edge_edge.first = edges[ie].index_ph == 0 ? ie : vl.first;
+ vlog[vx.index].hit_edge_edge.second = edges[ie].index_ph == 1 ? ie : vl.first;
+ }
+ }
+ }
+ } else if ( vlog[vx.index].index_ph == edges[ie].index_ph ) {
+ for (auto vl : vlog[vx.index].edges) {
+ if ( edges[vl.first].index_ph == edges[ie].index_ph ) continue;
+ if ( vx.hitEdge < 0 ) vx.hitEdge = vl.first;
+ }
+ }
+ }
+ }
+
+ // loop over intersections to remove double hits
+ std::vector<NPH_vtx>::iterator vIt; std::vector<NPH_vtx>::iterator wIt;
+ for (unsigned int ie = 0; ie < edges.size(); ie++ ) {
+ vIt = edges[ie].evtx.begin(); int vin = edges[ie].endpoints.first;
+ while ( vIt != edges[ie].evtx.end() ) {
+ if ( (*vIt).index == vin ) edges[ie].evtx.erase(vIt);
+ else { vin = (*vIt).index; vIt++;}
+ }
+ if (vin == edges[ie].endpoints.second && edges[ie].evtx.size()>0 ) edges[ie].evtx.pop_back();
+ }
+
+ // check symmetrization
+ for (unsigned int ie = 0; ie < edges.size(); ie++ ) {
+ for ( auto iv : edges[ie].evtx) {
+ bool foundHit = false;
+ for (auto jv : edges[edges[ie].edge_inverse].evtx) if ( jv.index == iv.index) foundHit = true;
+ if (!foundHit) {
+ NPH_vtx nvtx(iv.index); nvtx.edge = edges[ie].edge_inverse ; nvtx.hitPlane = iv.hitPlane; nvtx.distance = 1.-iv.distance;
+ nvtx.faces = std::pair<int,int>(iv.faces.second,iv.faces.first);
+ saveHitOnEdge( nvtx, edges);
+ }
+ }
+ }
+
+}
+
+
+void GeoPolyhedronManipulator::intersectFace(unsigned int in_edge, std::vector<NPH_edge>& edges, unsigned int plane_index, std::vector<NPH_plane>& planes, std::vector<GeoTrf::Vector3D>& vtx, std::vector<NPH_vertex>& vlog) const
+{
+ // define probe
+ GeoTrf::Vector3D point = vtx[edges[in_edge].endpoints.first];
+ GeoTrf::Vector3D probe = vtx[edges[in_edge].endpoints.second] - point;
+ double pn = probe.norm();
+
+ double d,dist,t; double en;
+
+ for (unsigned int ie = 0; ie < edges.size(); ie++) {
+
+ if (edges[ie].faces.first != plane_index ) continue;
+
+ GeoTrf::Vector3D eInit = vtx[edges[ie].endpoints.first];
+ GeoTrf::Vector3D edge = vtx[edges[ie].endpoints.second] - eInit;
+ en = edge.norm();
+
+ // dist and t in units (probe, edge)
+ bool lineIntersect = intersectLines( point, probe, eInit, edge, dist, t, d);
+ if (!lineIntersect) {
+ if (std::abs(d)< m_tolerance) { // overlapping edges ? TODO : figure out the actual overlap
+ setDoubleEdge(ie,in_edge,edges);
+ }
+ continue;
+ }
+ if ( std::abs(d) > m_tolerance || t*en < -m_tolerance || t*en > en+m_tolerance || dist*pn < -m_tolerance || dist*pn > pn+m_tolerance ) continue;
+ // valid intersection - may be known already
+ int index = -1; bool foundHit;
+ if (std::abs(t*en)<m_tolerance) index = edges[ie].endpoints.first;
+ else if (std::abs(dist*pn)<m_tolerance) index = edges[in_edge].endpoints.first;
+ else if ( std::abs(t*en -en)<m_tolerance) index = edges[ie].endpoints.second;
+ else if ( std::abs(dist*pn -pn)<m_tolerance) index = edges[in_edge].endpoints.second;
+ else {
+ GeoTrf::Vector3D hit = eInit + t*edge;
+ for (unsigned int jv = 0; jv<vtx.size(); jv++) if ( (hit - vtx[jv]).norm() < m_tolerance ) { index = jv ; break ;}
+ if (index < 0) { vtx.push_back(hit); index = vtx.size()-1; NPH_vertex pv; pv.index = index; pv.shared = true; vlog.push_back(pv); }
+ }
+ // save on probe edge
+ NPH_vtx nvtx(index); nvtx.edge=in_edge ; nvtx.hitPlane = edges[in_edge].faces.second ; nvtx.distance = dist; nvtx.faces = std::pair<int,int>(edges[ie].faces);
+ saveHitOnEdge(nvtx,edges);
+ // save on plane edge
+ NPH_vtx nptx(index); nptx.edge=ie ; nptx.hitPlane = edges[ie].faces.second ; nptx.distance = t; nptx.faces = std::pair<int,int>(edges[in_edge].faces);
+ saveHitOnEdge(nptx,edges);
+ //
+ vlog[index].edges.push_back(std::pair<int,int>(in_edge,edges[ie].faces.second));
+ vlog[index].edges.push_back(std::pair<int,int>(ie,edges[in_edge].faces.second));
+ vlog[index].hit_edge_edge.first = edges[in_edge].index_ph == 0 ? in_edge : ie;
+ vlog[index].hit_edge_edge.second = edges[in_edge].index_ph == 1 ? ie : in_edge;
+ // add the hit vertex to the list of new vertices on plane
+ bool foundOnFace = false; for (auto iv : planes[plane_index].hit_vertices ) if ( iv.index == index ) { foundOnFace = true; break; }
+ if (!foundOnFace) {
+ NPL_vertex plv; plv.index = index; plv.shared = true; plv.faces = edges[in_edge].faces;
+ planes[plane_index].hit_vertices.push_back(plv);
+ }
+ bool foundOnFaca = false; for (auto iv : planes[edges[ie].faces.second].hit_vertices ) if ( iv.index == index ) { foundOnFaca = true; break; }
+ if (!foundOnFaca) { // vertex saved on adjacent face
+ NPL_vertex plva; plva.index = index; plva.shared = true; plva.faces = edges[in_edge].faces; // plva.faces = edges[in_edge].faces;
+ planes[edges[ie].faces.second].hit_vertices.push_back(plva);
+ }
+ bool foundOnFace1 = false; for (auto iv : planes[edges[in_edge].faces.first].hit_vertices ) if ( iv.index == index ) { foundOnFace1 = true; break; }
+ if (!foundOnFace1) { // vertex to be saved on probe edge faces
+ NPL_vertex plv1; plv1.index = index; plv1.shared = true; plv1.faces = edges[ie].faces;
+ planes[edges[in_edge].faces.first].hit_vertices.push_back(plv1);
+ }
+ bool foundOnFace2 = false; for (auto iv : planes[edges[in_edge].faces.second].hit_vertices ) if ( iv.index == index ) { foundOnFace2 = true; break; }
+ if (!foundOnFace2) { // vertex to be saved on probe edge faces
+ NPL_vertex plv2; plv2.index = index; plv2.shared = true; plv2.faces = edges[ie].faces;
+ planes[edges[in_edge].faces.second].hit_vertices.push_back(plv2);
+ }
+ } // end loop over face edges
+
+ return;
+}
+
+NPL_vertex GeoPolyhedronManipulator::filterLog(NPH_vertex vlog, std::vector<NPH_edge> edges, int ip, std::vector<NPH_plane> planes) const {
+
+ NPL_vertex plv;
+ plv.index = vlog.index;
+ plv.shared = vlog.shared;
+ plv.overlapping_face = -1;
+ plv.oriented = false;
+
+ int edgeIn = -1;
+ // original vertices define orientation
+ if ( vlog.index_ph == planes[ip].index_ph ) {
+ for ( auto ee : vlog.edges ) {
+ if ( ee.second == ip && edges[ee.first].index_ph == vlog.index_ph ) { plv.faces = edges[ee.first].faces; edgeIn = ee.first; }
+ }
+ return plv;
+ }
+
+ for ( auto ee : vlog.edges ) {
+ if ( ee.second == ip) { plv.faces = edges[ee.first].faces; edgeIn = ee.first; }
+ }
+
+ // check overlaps
+ for ( auto ee : vlog.edges ) {
+ if ( ee.first == edgeIn && ee.second!=ip ) { plv.overlapping_face = ee.second ; plv.shared = true; } // shared used to be 2 here, verify
+ }
+ return plv;
+}
+
+void GeoPolyhedronManipulator::filterLogPlane(NPL_vertex& plv, NPH_vertex& vlog, std::vector<NPH_edge> edges, int ip, std::vector<NPH_plane> planes) const {
+
+ for ( auto ee : vlog.edges ) {
+ if ( ee.second == ip ) {
+ plv.faces = edges[ee.first].faces;
+ if (edges[ee.first].index_ph == planes[ip].index_ph) plv.oriented = true;
+ }
+ }
+ for (auto ov : planes[ip].overlapping) {
+ for ( auto ee : vlog.edges ) {
+ if ( ee.second == ov%1000 ) {
+ plv.faces_shared_loop = edges[ee.first].faces;
+ }
+ }
+ }
+
+}
+
+bool GeoPolyhedronManipulator::intersectLines(GeoTrf::Vector3D A, GeoTrf::Vector3D a, GeoTrf::Vector3D B, GeoTrf::Vector3D b, double& da, double& db, double& d) const
+{
+ // lines parametrization : A + a*da , B + b*db ; a,b define "units" of distances da, db
+ // return value : closest distance between lines in native [mm] units
+
+ // scalars
+ double ab = a.dot(b); double a2 = a.dot(a); double b2 = b.dot(b);
+ double ABa = (A-B).dot(a);
+ double ABb = (A-B).dot(b);
+
+ // special cases
+ //if ( !ab ) {
+ if (std::abs(ab)/sqrt(a2*b2) < 1.e-9 ) { // to ensure sufficient precision
+ da = -ABa / a2;
+ db = ABb / b2;
+ d = (A + da*a -B - db*b).norm() ;
+ return true;
+ }
+
+ if ( a2*b2 == ab*ab ) { // parallel lines
+ da = 0; db = 0;
+ d = ((A-B).cross(b)).norm()/b.norm();
+ return false;
+ }
+
+ da = (ABb*ab -b2*ABa) / (a2*b2 - ab*ab);
+ db = (ABa + a2*da) / ab;
+ d = (A + da*a -B - db*b).norm() ;
+ return true;
+
+}
+
+// method to find intersection of a (polyhedron) edge with (another) polyhedron face plane
+// arguments : edge: (1) starting point, (2) endpoint-starting point, plane: (3) point , (4) normal
+double GeoPolyhedronManipulator::edgeFaceIntersection( GeoTrf::Vector3D eInit, GeoTrf::Vector3D edge, GeoTrf::Vector3D facePoint, GeoTrf::Vector3D faceNormal ) const
+{
+ if (!edge.dot(faceNormal)) {
+ std::cout << "ERROR, edgeFace intersection called for edge parallel to plane" << std::endl;
+ return 1.e10;
+ }
+
+ // scalars
+ double en = edge.dot(faceNormal);
+ double pen = ( facePoint - eInit ).dot(faceNormal);
+
+ return pen / en ;
+}
+
+std::vector<NPL_edge> GeoPolyhedronManipulator::evaluateVertexLogic(int plane_index, std::vector<NPH_vertex>& vlog, std::vector<NPH_edge> edges, std::vector<NPH_plane>& planes, std::vector<GeoTrf::Vector3D> vtx) const {
+
+ // loop over vertices to determine the level of ambiguity : results save as NPL_vertex properties
+ // default mode processes hit vertices ( intersections )
+
+ NPH_plane plane = planes[plane_index];
+
+ std::vector<NPL_edge> lines; std::vector< std::pair< int, std::vector< int > > > active_intersected;
+
+ for (unsigned int ih = 0; ih< plane.hit_vertices.size(); ih++) { NPL_vertex iplv = plane.hit_vertices[ih]; NPH_vertex iv = vlog[iplv.index];
+ // overview of adjacent planes
+ std::vector<int> np; std::vector<int> native_planes; std::vector<int> intersected_planes; int pl; bool found;
+ for (auto ie : iv.edges) {
+ if (edges[ie.first].index_ph==0) iv.hit_edge_edge.first = ie.first;
+ else if ( edges[ie.first].index_ph==1) iv.hit_edge_edge.second = ie.first;
+ //
+ pl = edges[ie.first].faces.first;
+ if (pl>=0) { found = false; for (auto pp : np ) if (pl==pp) {found=true; break; } if (!found) np.push_back(pl);}
+ pl = edges[ie.first].faces.second;
+ if (pl>=0) { found = false; for (auto pp : np ) if (pl==pp) {found=true; break; } if (!found) np.push_back(pl);}
+ pl = ie.second;
+ if (pl>=0) { found = false; for (auto pp : np ) if (pl==pp) {found=true; break; } if (!found) np.push_back(pl);}
+ }
+ //
+ for (auto pp : np ) if (planes[pp].index_ph==plane.index_ph) native_planes.push_back(pp); else intersected_planes.push_back(pp);
+ for (auto ip : intersected_planes) { bool foundIntersected = false ;
+ for (unsigned int ii=0; ii<active_intersected.size(); ii++) {
+ if ( active_intersected[ii].first == ip ) { foundIntersected = true;
+ active_intersected[ii].second.push_back(iplv.index);
+ }
+ }
+ if (!foundIntersected) active_intersected.push_back(std::pair<int, std::vector<int> > (ip, std::vector<int> (1,iplv.index)) );
+ }
+ //
+ if (iv.index_ph < 0) { // new vertex / mixed planes : "native" fit the top plane index
+ //
+ if ( native_planes.size()==1 && intersected_planes.size() == 2 ) { filterLogPlane(iplv, iv, edges, plane_index, planes); iplv.oriented = false; }
+ else if ( native_planes.size()==2 && intersected_planes.size() == 1 ) {
+ std::vector<int> hitFaces = findHitPlane( iv, plane_index, edges, planes);
+ if ( hitFaces.size() >1 ) { iplv.faces.first = hitFaces[0]; iplv.faces.second = hitFaces[1]; }
+ } else { // see if number of hit faces reasonable when taking into account overlaps
+ std::vector<int> hitFaces = findHitPlane( iv, plane_index, edges, planes);
+ std::vector<int>::iterator hIt = hitFaces.begin();
+ while ( hIt!= hitFaces.end() ) {
+ if ( matchPlanes(plane_index, (*hIt), planes ) ) hitFaces.erase(hIt);
+ else hIt++;
+ }
+ if (hitFaces.size() == 2) { iplv.faces.first = hitFaces[0]; iplv.faces.second = hitFaces[1]; }
+ else if (intersected_planes.size()>=2) {
+ std::vector<int> nov; for ( auto ii : intersected_planes ) if ( !matchPlanes(plane_index,ii,planes) ) nov.push_back(ii);
+ if (nov.size()==2) { iplv.faces.first = nov[0]; iplv.faces.second = nov[1]; }
+ else { iplv.faces.first = -1; iplv.faces.second = -1; }
+ } else { iplv.faces.first = -1; iplv.faces.second = -1; }
+ }
+
+ } else {
+
+ filterLogPlane( planes[plane_index].hit_vertices[ih],iv,edges,plane_index, planes);
+ }
+
+ plane.hit_vertices[ih].overlapping_face = -1;
+ for (auto ov : planes[plane_index].overlapping) if (inPlane(iv,edges,ov%1000) ) plane.hit_vertices[ih].overlapping_face=ov%1000;
+
+ if ( np.size() == 5) {
+ std::pair<int,int> doublEdg = hitDoubleEdge(iplv.index, plane_index, plane.index_ph, vlog, edges); // indicates adjacent face
+ if (doublEdg.first>=0 && iv.index_ph==plane.index_ph) {
+ vlog[iplv.index].hitDoublEdge = doublEdg.first;
+ vlog[iplv.index].onDoublEdge = (planes[doublEdg.first].index_ph == plane.index_ph) ? 0 : 1;
+ }
+ if (doublEdg.first>=0 && iv.index_ph!=plane.index_ph) {
+ planes[plane_index].hit_vertices[ih].faces.first = doublEdg.second;
+ // check if valid intersecting plane
+ bool validPlane = false;
+ for ( unsigned int jh=0; jh<plane.hit_vertices.size(); jh++) { if (jh==ih) continue;
+ if ( inPlane(vlog[plane.hit_vertices[jh].index],edges,doublEdg.first) ) validPlane = true;
+ }
+ planes[plane_index].hit_vertices[ih].faces.second = validPlane ? doublEdg.first : doublEdg.second;
+ if (!validPlane) iplv.onEdgeDummy = true;
+ vlog[iplv.index].hitDoublEdge = validPlane ? doublEdg.first : doublEdg.second ;
+ vlog[iplv.index].onDoublEdge = validPlane ? 1 : 0;
+ }
+ }
+ }
+
+ for (auto ii : active_intersected) { // std::cout << "protoline?" << ii.first <<":" << ii.second.size() ;
+ bool ol = false; for (auto ov : plane.overlapping) if ( ov%1000 == ii.first) ol = true;
+ if (ol) continue; // do not consider overlapping faces
+ if ( ii.second.size()<2) continue; // skip incomplete connections
+ // order vertices along line ( needs acces to vertex position )
+ std::vector<std::pair<int,double>> vxord;
+ vxord.push_back(std::pair<int,double>(ii.second[0],0.) );
+ vxord.push_back(std::pair<int,double>(ii.second[1],1.) );
+ double al = (vtx[ii.second[1]] -vtx[ii.second[0]]).norm();
+ std::vector<std::pair<int,double> >::iterator ita;
+ for (unsigned int iv = 2; iv<ii.second.size(); iv++) {
+ double da = (vtx[ii.second[iv]] -vtx[ii.second[0]]).dot(vtx[ii.second[1]] -vtx[ii.second[0]])/al/al;
+ ita = vxord.begin();
+ while ( ita!= vxord.end() ) {
+ if (da < (*ita).second ) { vxord.insert(ita,std::pair<int,double>(ii.second[iv],da) ); break; }
+ else ita++;
+ }
+ if (ita==vxord.end()) vxord.insert(ita,std::pair<int,double>(ii.second[iv],da) );
+ }
+ // count number of hit vtx
+ int nhit = 0; int native = 0; std::vector<int> hits;
+ for (auto vv : vxord) {
+ if (vlog[vv.first].index_ph < 0) { nhit++;
+ for (unsigned int ih=0; ih<plane.hit_vertices.size(); ih++) {
+ if (vv.first == plane.hit_vertices[ih].index) hits.push_back(ih);
+ }
+ }
+ if (vlog[vv.first].index_ph == planes[plane_index].index_ph ) native++;
+ }
+
+ int noe = 0; // not-on-edge
+ for (auto vv : vxord) {
+ for (unsigned int ih=0; ih<plane.hit_vertices.size(); ih++) {
+ if (vv.first == plane.hit_vertices[ih].index) {
+ if ( plane.hit_vertices[ih].onEdge < 0 ) noe++;
+ break;
+ }
+ }
+ }
+
+ // examine different configurations
+ bool splitLine = true;
+
+ if (vxord.size()>2 && vlog[vxord.front().first].index_ph < 0 && vlog[vxord.back().first].index_ph < 0 ) splitLine = false;
+ if (noe>0) splitLine = true;
+
+ // create line edge
+ for (unsigned int jl =0; jl<vxord.size()-1; jl++) {
+ lines.push_back(createLine(vxord[jl].first,vxord[jl+1].first,ii.first, plane) );
+ if (!splitLine) lines.back().protectedLine = true;
+ }
+
+ for (auto iv : ii.second) { // std::cout <<":index:"<< iv;
+ for (unsigned int ih=0; ih<plane.hit_vertices.size(); ih++) {
+ if (iv == plane.hit_vertices[ih].index) {
+ //std::cout <<":faces:" << plane.hit_vertices[ih].faces.first <<"," << plane.hit_vertices[ih].faces.second ;
+ planes[plane_index].hit_vertices[ih].active_lines.push_back(ii.first);
+ }
+ }
+ }
+ //std::cout <<""<< std::endl;
+ }
+
+ return lines;
+}
+
+bool GeoPolyhedronManipulator::inPlane(NPH_vertex vlog, std::vector<NPH_edge> edges, int ip) const {
+
+ bool planeVtx = false;
+
+ int f1,f2,f3; bool vin = false;
+ for ( auto ef : vlog.edges ) {
+ f1 = edges[ef.first].faces.first;
+ f2 = edges[ef.first].faces.second;
+ f3 = ef.second;
+ vin = ( (ip==f1) || (ip==f2) || (ip==f3) ) ;
+ if (vin) planeVtx = true;
+ }
+
+ return planeVtx;
+
+}
+
+bool GeoPolyhedronManipulator::matchPlanes( int p1, int p2, std::vector<NPH_plane> planes) const {
+
+ bool match = false;
+
+ if (p1 == p2) return match=true;
+
+ for ( auto ov : planes[p1].overlapping ) if ( ov%1000 == p2 ) match = true;
+
+ return match;
+
+}
+
+std::vector<int> GeoPolyhedronManipulator::findHitPlane( NPH_vertex vlog, int plane_index, std::vector<NPH_edge> edges, std::vector<NPH_plane> planes) const {
+
+ // ph_index stands for the polyhedron index of the currently processed plane
+ std::vector<int> faceOut; bool foundFace;
+ for ( auto ee : vlog.edges ) {
+ if ( !matchPlanes(edges[ee.first].faces.first,plane_index,planes) ) {
+ foundFace = false ; for (auto iv : faceOut ) if ( iv == edges[ee.first].faces.first ) foundFace = true;
+ if (!foundFace ) faceOut.push_back(edges[ee.first].faces.first);
+ }
+ if ( !matchPlanes(edges[ee.first].faces.second,plane_index,planes) ) {
+ foundFace = false ; for (auto iv : faceOut ) if ( iv == edges[ee.first].faces.second ) foundFace = true;
+ if (!foundFace) faceOut.push_back(edges[ee.first].faces.second);
+ }
+ if ( !matchPlanes(ee.second,plane_index,planes) ) {
+ foundFace = false ; for (auto iv : faceOut ) if ( iv == ee.second ) foundFace = true;
+ if (!foundFace) faceOut.push_back(ee.second);
+ }
+ }
+
+ return faceOut;
+}
+
+std::pair<int,int> GeoPolyhedronManipulator::hitDoubleEdge(int index, int plane_index, int ph_index, std::vector<NPH_vertex> vlog, std::vector<NPH_edge> edges) const {
+
+ int hitFace = -1; int adjFace = -1;
+
+ int edgeExt = vlog[index].index_ph==0 ? vlog[index].hit_edge_edge.second : vlog[index].hit_edge_edge.first ;
+ int edgeInt = ph_index==0 ? vlog[index].hit_edge_edge.first : vlog[index].hit_edge_edge.second ;
+ if ( edgeExt < 0 ) return std::pair<int,int>(hitFace,adjFace);
+
+ if (edges[edgeInt].faces.first == plane_index) adjFace = edges[edgeInt].faces.second;
+ else adjFace = edges[edgeInt].faces.first;
+
+ for (auto ie : vlog[index].edges) { if ( edges[ie.first].index_ph != vlog[index].index_ph ) continue;
+ bool inPlaneFirst = false; for ( auto ii : edges[ie.first].inPlane ) if ( ii== edges[edgeExt].faces.first ) inPlaneFirst = true;
+ bool inPlaneSecond = false; for ( auto ii : edges[ie.first].inPlane ) if ( ii== edges[edgeExt].faces.second ) inPlaneSecond = true;
+ if (inPlaneFirst && inPlaneSecond ) hitFace = ie.second;
+ }
+
+ return std::pair<int,int>(hitFace,adjFace);
+
+}
+
+void GeoPolyhedronManipulator::setDoubleEdge(int ie,int ia, std::vector<NPH_edge>& edges) const {
+
+ bool foundEdge = false; for (auto je : edges[ie].doubleEdge) if (je == ia) foundEdge = true;
+ if (!foundEdge) edges[ie].doubleEdge.push_back(ia);
+
+ foundEdge = false; for (auto je : edges[ia].doubleEdge) if (je == ie) foundEdge = true;
+ if (!foundEdge) edges[ia].doubleEdge.push_back(ie);
+
+}
+
+NPL_edge GeoPolyhedronManipulator::symmetrize( NPL_edge line ) const {
+
+ NPL_edge sym(std::pair<int,int>(line.endpoints.second,line.endpoints.first), line.adjacent_face);
+ sym.endpoint_faces.first = line.endpoint_faces.second;
+ sym.endpoint_faces.second = line.endpoint_faces.first;
+ sym.endpoint_hitvx.first = line.endpoint_hitvx.second;
+ sym.endpoint_hitvx.second = line.endpoint_hitvx.first;
+ sym.overlap = line.overlap;
+ sym.shared = line.shared;
+ sym.protectedLine = line.protectedLine;
+ sym.dummy_overlap = line.dummy_overlap;
+ return sym;
+}
+
+
+void GeoPolyhedronManipulator::saveHitOnEdge(NPH_vtx nvtx, std::vector<NPH_edge>& edges) const{
+
+ if (edges[nvtx.edge].endpoints.first == nvtx.index ) return;
+ if (edges[nvtx.edge].endpoints.second == nvtx.index ) return;
+ for (auto iv : edges[nvtx.edge].evtx ) if (iv.index == nvtx.index) return;
+
+ std::vector<NPH_vtx>::iterator vIt = edges[nvtx.edge].evtx.begin();
+ while ( vIt< edges[nvtx.edge].evtx.end()) {
+ if ( (*vIt).distance < nvtx.distance) vIt++;
+ else break;
+ }
+ edges[nvtx.edge].evtx.insert(vIt,nvtx);
+ //std::cout <<"inserting evtx(method):" << nvtx.index <<":at:" << nvtx.distance <<":endpoints:" <<
+ // edges[nvtx.edge].endpoints.first << "," << edges[nvtx.edge].endpoints.second<<":hit plane:" << nvtx.hitPlane<< std::endl;
+}
+
+NPL_edge GeoPolyhedronManipulator::createLine(int iv1, int iv2, int adj, NPH_plane plane) const {
+
+ NPL_edge line( std::pair<int,int>(iv1,iv2),adj); line.generating_edge = -1; line.shared = true;
+
+ int h1=-1; int h2=-1;
+ for (int ih=0; ih<plane.hit_vertices.size(); ih++) {
+ if ( iv1 == plane.hit_vertices[ih].index) h1 = ih;
+ if ( iv2 == plane.hit_vertices[ih].index) h2 = ih;
+ }
+
+ if (h1>=0) {
+ line.endpoint_hitvx.first = h1;
+ line.endpoint_faces.first = adj == plane.hit_vertices[h1].faces.first ?
+ plane.hit_vertices[h1].faces.second : plane.hit_vertices[h1].faces.first ;
+ }
+ if (h2>=0) {
+ line.endpoint_hitvx.second = h2;
+ line.endpoint_faces.second = adj == plane.hit_vertices[h2].faces.first ?
+ plane.hit_vertices[h2].faces.second : plane.hit_vertices[h2].faces.first ;
+ }
+
+ return line;
+}
+
+NPL_ordered GeoPolyhedronManipulator::mergeLoops( NPL_ordered loop1, NPL_ordered loop2, std::vector< GeoTrf::Vector3D > vtx) {
+
+ NPL_ordered merged = loop1;
+
+ if ( !loop2.vtx.size() ) return merged;
+
+ std::vector<size_t> shared;
+
+ for (auto iv1 : loop1.vtx) {
+ for ( auto iv2 : loop2.vtx) {
+ if ( iv1 == iv2 ) shared.push_back(iv1);
+ }
+ }
+
+ if (shared.size()) {
+
+ loop2.reorder(shared[0]);
+ merged.reorder(shared[0]);
+ merged.vtx.insert(merged.vtx.begin(),loop2.vtx.begin(),loop2.vtx.end());
+ merged.edges.insert(merged.edges.begin(),loop2.edges.begin(),loop2.edges.end());
+
+ } else { // there is no common vertex : new dummy symmetrized edge created
+
+ // find closest
+ double dist = 1.e10; int im1 = -1; int im2 = -1;
+ for (unsigned int iv1=0; iv1<loop1.vtx.size(); iv1++) {
+ for (unsigned int iv2=0; iv2<loop2.vtx.size(); iv2++) {
+ double d = ( vtx[loop1.vtx[iv1]] - vtx[loop2.vtx[iv2]] ).norm();
+ if ( d < dist ) {
+ im1 = loop1.vtx[iv1] ; im2 = loop2.vtx[iv2] ; dist = d;
+ }
+ }}
+ merged.reorder(im1);
+ loop2.reorder(im2);
+
+ merged.vtx.push_back(merged.vtx.front());
+ loop2.vtx.push_back(loop2.vtx.front());
+
+ merged.vtx.insert(merged.vtx.begin(),loop2.vtx.begin(),loop2.vtx.end());
+ }
+
+ return merged;
+}
+
+
+GeoTrf::Vector3D GeoPolyhedronManipulator::areaVec( std::vector<size_t> loop, std::vector< GeoTrf::Vector3D > vtx) const {
+
+ int n = loop.size();
+ GeoTrf::Vector3D area = vtx[loop[0]].cross(vtx[loop[n - 1]]);
+ for (int k = 1; k < n; ++k) area += vtx[loop[k]].cross(vtx[loop[k-1]]);
+
+ return area;
+
+}