Tolti commenti inutili e sistemata interfaccia

- inserita opzione in compilazioni 64bit per evitare Warning su verifica riferimenti non nulli.

CurveComposite.cpp

@@ -30,6 +30,7 @@
 #include "/EgtDev/Include/EGkSfrCreate.h"
 #include "/EgtDev/Include/EGkIntervals.h"
 #include "/EgtDev/Include/EGkStringUtils3d.h"
+#include "/EgtDev/Include/EgtNumUtils.h"
 #include "/EgtDev/Include/EgtPointerOwner.h"
 #include <algorithm>
 
@@ -1321,8 +1322,8 @@ CurveComposite::ApproxWithLines( double dLinTol, double dAngTolDeg, int nType, P
   // se lineare con lato obbligato...
    if ( nType == APL_LEFT  ||  nType == APL_LEFT_CONVEX  ||
         nType == APL_RIGHT  ||  nType == APL_RIGHT_CONVEX) {
-     // prima approssimazione lineare a 10 * Epsilon
-      if ( ! ApproxWithLines( 10 * EPS_SMALL, dAngTolDeg, APL_SPECIAL, PL))
+     // prima approssimazione lineare alla tolleranza minima del programma
+      if ( ! ApproxWithLines( EPS_SMALL, dAngTolDeg, APL_SPECIAL, PL))
          return false ;
      // eliminazione dei punti in tolleranza andando solo dalla parte ammessa
       Vector3d vtExtr = ( m_VtExtr.IsSmall() ? Z_AX : m_VtExtr) ;

EgtGeomKernel.vcxproj

@@ -151,6 +151,7 @@ copy $(TargetPath) \EgtProg\DllD32</Command>
       <MultiProcessorCompilation>true</MultiProcessorCompilation>
       <DebugInformationFormat>ProgramDatabase</DebugInformationFormat>
       <LanguageStandard>stdcpp17</LanguageStandard>
+      <AdditionalOptions>-Wno-tautological-undefined-compare</AdditionalOptions>
     </ClCompile>
     <Link>
       <SubSystem>Windows</SubSystem>
@@ -243,6 +244,7 @@ copy $(TargetPath) \EgtProg\Dll32</Command>
       <WholeProgramOptimization>false</WholeProgramOptimization>
       <DebugInformationFormat>None</DebugInformationFormat>
       <LanguageStandard>stdcpp17</LanguageStandard>
+      <AdditionalOptions>-Wno-tautological-undefined-compare</AdditionalOptions>
     </ClCompile>
     <Link>
       <SubSystem>Windows</SubSystem>
@@ -392,6 +394,7 @@ copy $(TargetPath) \EgtProg\Dll64</Command>
     <ClCompile Include="SurfFlatRegionBooleans.cpp" />
     <ClCompile Include="SurfFlatRegionOffset.cpp" />
     <ClCompile Include="SurfTriMeshBooleans.cpp" />
+    <ClCompile Include="SurfTriMeshCuts.cpp" />
     <ClCompile Include="SurfTriMeshUtilities.cpp" />
     <ClCompile Include="TextureData.cpp" />
     <ClCompile Include="Tool.cpp" />
@@ -558,6 +561,7 @@ copy $(TargetPath) \EgtProg\Dll64</Command>
     <ClInclude Include="DistPointCrvComposite.h" />
     <ClInclude Include="DistPointLine.h" />
     <ClInclude Include="DllMain.h" />
+    <ClInclude Include="earcut.hpp" />
     <ClInclude Include="ExtDimension.h" />
     <ClInclude Include="ExtText.h" />
     <ClInclude Include="FontAux.h" />
@@ -586,6 +590,7 @@ copy $(TargetPath) \EgtProg\Dll64</Command>
     <ClInclude Include="IntersArcArc.h" />
     <ClInclude Include="IntersCrvCompoCrvCompo.h" />
     <ClInclude Include="IntersLineArc.h" />
+    <ClInclude Include="IntersLineBox.h" />
     <ClInclude Include="IntersLineLine.h" />
     <ClInclude Include="IntersLineSurfStd.h" />
     <ClInclude Include="IntersLineTria.h" />

EgtGeomKernel.vcxproj.filters

@@ -456,6 +456,9 @@
     <ClCompile Include="SurfTriMeshUtilities.cpp">
       <Filter>File di origine\Geo</Filter>
     </ClCompile>
+    <ClCompile Include="SurfTriMeshCuts.cpp">
+      <Filter>File di origine\Geo</Filter>
+    </ClCompile>
   </ItemGroup>
   <ItemGroup>
     <ClInclude Include="stdafx.h">
@@ -1079,6 +1082,12 @@
     <ClInclude Include="CDeRectPrismoidTria.h">
       <Filter>File di intestazione</Filter>
     </ClInclude>
+    <ClInclude Include="earcut.hpp">
+      <Filter>File di intestazione</Filter>
+    </ClInclude>
+    <ClInclude Include="IntersLineBox.h">
+      <Filter>File di intestazione</Filter>
+    </ClInclude>
   </ItemGroup>
   <ItemGroup>
     <ResourceCompile Include="EgtGeomKernel.rc">

HashGrids3d.cpp

@@ -684,7 +684,6 @@ HashGrids3d::Update( void)
                double dSize = 0 ;
                Obj.box.GetDiameter( dSize) ;
                double dCellSpan = pGrid->GetCellSpan() ;
-
                if ( dSize >= dCellSpan  ||  dSize < ( dCellSpan / hierarchyFactor)) {
                   pGrid->Remove( Obj) ;
                   addGrid( Obj) ;

IntersLineBox.cpp

@@ -1,7 +1,7 @@
 //----------------------------------------------------------------------------
-//  EgalTech 2020-2020
+//  EgalTech 2020-2022
 //----------------------------------------------------------------------------
-// File : IntersLineBox.cpp       Data : 07.05.20         Versione : 2.2e1
+// File : IntersLineBox.cpp       Data : 08.01.22         Versione : 2.4a3
 // Contenuto : Implementazione della intersezione linea/box.
 //
 //
@@ -13,13 +13,77 @@
 
 //--------------------------- Include ----------------------------------------
 #include "stdafx.h"
-#include "SurfTriMesh.h"
+#include "IntersLineBox.h"
 #include "/EgtDev/Include/EGkIntersLineBox.h"
-#include "/EgtDev/Include/EGkIntersLineSurfTm.h"
+#include "/EgtDev/Include/EgtNumUtils.h"
 #include <algorithm>
 
 using namespace std ;
 
+//----------------------------------------------------------------------------
+// Linea e box allineato assi devono essere nel medesimo sistema di riferimento.
+// In caso di intersezione viene restituito true e i parametri in dU1 e dU2.
+//----------------------------------------------------------------------------
+bool
+IntersLineBox( const Point3d& ptL, const Vector3d& vtL,
+               const Point3d& ptMin, const Point3d& ptMax,
+               double& dU1, double& dU2)
+{
+  // Verifico il versore
+   if ( vtL.IsSmall())
+      return false ;
+
+  // Verifico gli estremi del box
+   if ( ptMin.x > ptMax.x + EPS_SMALL  ||
+        ptMin.y > ptMax.y + EPS_SMALL  ||
+        ptMin.z > ptMax.z + EPS_SMALL)
+      return false ;
+
+  // Casi di intersezione impossibile
+   if ( abs( vtL.x) < EPS_ZERO  &&  ( ptL.x < ptMin.x - EPS_SMALL  ||  ptL.x > ptMax.x + EPS_SMALL))
+      return false ;
+   if ( abs( vtL.y) < EPS_ZERO  &&  ( ptL.y < ptMin.y - EPS_SMALL  ||  ptL.y > ptMax.y + EPS_SMALL))
+      return false ;
+   if ( abs( vtL.z) < EPS_ZERO  &&  ( ptL.z < ptMin.z - EPS_SMALL  ||  ptL.z > ptMax.z + EPS_SMALL))
+      return false ;
+
+  // Parametri degli estremi della retta
+   dU1 = -INFINITO ;
+   dU2 = INFINITO ;
+
+  // Confronto con piani YZ (perpendicolari ad asse X)
+   if ( vtL.x > EPS_ZERO) {
+      dU1 = max( dU1, ( ptMin.x - ptL.x) / vtL.x) ;
+      dU2 = min( dU2, ( ptMax.x - ptL.x) / vtL.x) ;
+   }
+   else if ( vtL.x < -EPS_ZERO) {
+      dU1 = max( dU1, ( ptMax.x - ptL.x) / vtL.x) ;
+      dU2 = min( dU2, ( ptMin.x - ptL.x) / vtL.x) ;
+   }
+
+  // Confronto con piani ZX (perpendicolari ad asse Y)
+   if ( vtL.y > EPS_ZERO) {
+      dU1 = max( dU1, ( ptMin.y - ptL.y) / vtL.y) ;
+      dU2 = min( dU2, ( ptMax.y - ptL.y) / vtL.y) ;
+   }
+   else if ( vtL.y < -EPS_ZERO) {
+      dU1 = max( dU1, ( ptMax.y - ptL.y) / vtL.y) ;
+      dU2 = min( dU2, ( ptMin.y - ptL.y) / vtL.y) ;
+   }
+
+  // Confronto con piani XY (perpendicolari ad asse Z)
+   if ( vtL.z > EPS_ZERO) {
+      dU1 = max( dU1, ( ptMin.z - ptL.z) / vtL.z) ;
+      dU2 = min( dU2, ( ptMax.z - ptL.z) / vtL.z) ;
+   }
+   else if ( vtL.z < -EPS_ZERO){
+      dU1 = max( dU1, ( ptMax.z - ptL.z) / vtL.z) ;
+      dU2 = min( dU2, ( ptMin.z - ptL.z) / vtL.z) ;
+   }
+
+   return ( dU2 >= dU1) ;
+}
+
 //----------------------------------------------------------------------------
 bool
 IntersLineBox( const Point3d& ptL1, const Point3d& ptL2, const BBox3d& b3Box, INTDBLVECTOR& vInters, bool bFinite)
@@ -38,114 +102,69 @@ bool
 IntersLineBox( const Point3d& ptL, const Vector3d& vtL, double dLen, const BBox3d& b3Box, INTDBLVECTOR& vInters, bool bFinite)
 {
   // Recupero i dati del Box
-   Point3d ptMin ;
-   double dDimX, dDimY, dDimZ ;
-   if ( ! b3Box.GetMinDim( ptMin, dDimX, dDimY, dDimZ))
+   Point3d ptMin, ptMax ;
+   if ( ! b3Box.GetMinMax( ptMin, ptMax))
       return false ;
-  // Contorno di base
-   PolyLine PL ;
-   PL.AddUPoint( 0, ptMin) ;
-   PL.AddUPoint( 1, ptMin + Vector3d( dDimX, 0, 0)) ;
-   PL.AddUPoint( 2, ptMin + Vector3d( dDimX, dDimY, 0)) ;
-   PL.AddUPoint( 3, ptMin + Vector3d( 0, dDimY, 0)) ;
-   PL.Close() ;
-  // Vettore altezza
-   Vector3d vtExtr( 0, 0, dDimZ) ;
-  // Creo la superficie trimesh equivalente
-   SurfTriMesh Stm ;
-   if ( ! Stm.CreateByExtrusion( PL, vtExtr))
-      return false ;
-   SurfTriMesh StmBot ;
-   if ( ! StmBot.CreateByFlatContour( PL))
-      return false ;
-   StmBot.Invert() ;
-   SurfTriMesh StmTop ;
-   if ( ! StmTop.CreateByFlatContour( PL))
-      return false ;
-   StmTop.Translate( vtExtr) ;
-   if ( ! Stm.DoSewing( StmBot)  ||  ! Stm.DoSewing( StmTop))
-      return false ;
-  // Calcolo l'intersezione
-   ILSIVECTOR vInfo ;
-   if ( ! IntersLineSurfTm( ptL, vtL, dLen, Stm, vInfo, bFinite))
-      return false ;
-  // ciclo sulle intersezioni
-   double dUcurr = -INFINITO ;
-   for ( const auto& Info : vInfo) {
-     // se intersezione puntuale
-      if ( Info.nILTT == ILTT_VERT  ||  Info.nILTT == ILTT_EDGE  ||  Info.nILTT == ILTT_IN) {
-         int nFlag = ILBT_TOUCH ;
-         if ( Info.dCosDN > EPS_ZERO)
-            nFlag = ILBT_OUT ;
-         else if ( Info.dCosDN < -EPS_ZERO)
-            nFlag = ILBT_IN ;
-         vInters.emplace_back( nFlag, Info.dU) ;
-         dUcurr = Info.dU ;
+
+  // Pulisco vettore intersezioni
+   vInters.clear() ;
+
+  // Eseguo intersezione della linea completa
+   double dU1, dU2 ;
+   bool bInters = IntersLineBox( ptL, vtL, b3Box.GetMin(), b3Box.GetMax(), dU1, dU2) ;
+
+  // Se non c'è intersezione
+   if ( ! bInters  ||  ( bFinite  &&  ( dU1 > dLen + EPS_SMALL  ||  dU2 < -EPS_SMALL)))
+      return true ;
+
+  // Se le due intersezioni coincidono
+   if ( ( dU2 - dU1) < EPS_SMALL) {
+      double dU = ( dU1 + dU2) / 2 ;
+      if ( ! bFinite) {
+         vInters.emplace_back( ILBT_TOUCH, dU) ;
       }
-     // se altrimenti intersezione con coincidenza
-      else if ( Info.nILTT == ILTT_SEGM  ||  Info.nILTT == ILTT_SEGM_ON_EDGE) {
-         if ( Info.dU > dUcurr - EPS_SMALL)
-            vInters.emplace_back( ILBT_TG_INI, Info.dU) ;
-         vInters.emplace_back( ILBT_TG_FIN, Info.dU2) ;
-         dUcurr = Info.dU2 ;
+      else {
+         if ( dU > - EPS_SMALL  &&  dU < dLen + EPS_SMALL)
+            vInters.emplace_back( ILBT_TOUCH, Clamp( dU, 0., dLen)) ;
       }
-   }
-  // elimino intersezioni ripetute
-   for ( size_t j = 1 ; j < vInters.size() ; ) {
-     // intersezione precedente
-      size_t i = j - 1 ;
-      // se sono dello stesso tipo, elimino una delle due
-      if ( vInters[i].first == vInters[j].first) {
-        // se entranti elimino la prima
-         if ( vInters[i].first == ILBT_IN  ||  vInters[i].first == ILBT_TG_INI)
-            vInters.erase( vInters.begin() + i) ;
-        // altrimenti la seconda
-         else
-            vInters.erase( vInters.begin() + j) ;
-         if ( i > 0)
-            -- j ;
-         continue ;
-      }
-     // se hanno lo stesso parametro 
-      else if ( abs( vInters[i].second - vInters[j].second) < EPS_SMALL) {
-        // se una entrante e l'altra uscente, cambio in touch ed elimino la seconda
-         if ( ( vInters[i].first == ILBT_IN  &&  vInters[j].first == ILBT_OUT)  ||
-                   ( vInters[i].first == ILBT_OUT  &&  vInters[j].first == ILBT_IN)) {
-            vInters[i].first = ILBT_TOUCH ;
-            vInters.erase( vInters.begin() + j) ;
-            if ( i > 0)
-               -- j ;
-            continue ;
-         }
-        // se prima puntuale e l'altra inizio o fine di coincidenza, elimino la prima
-         else if ( ( vInters[i].first == ILBT_IN  ||  vInters[i].first == ILBT_OUT  ||  vInters[i].first == ILBT_TOUCH)  &&
-                   ( vInters[j].first == ILBT_TG_INI  ||  vInters[j].first == ILBT_TG_FIN)) {
-            vInters.erase( vInters.begin() + i) ;
-            if ( i > 0)
-               -- j ;
-            continue ;
-         }
-        // se prima inizio o fine di coincidenza e l'altra puntuale, elimino la seconda
-         else if ( ( vInters[i].first == ILBT_TG_INI  ||  vInters[i].first == ILBT_TG_FIN)  &&
-                   ( vInters[j].first == ILBT_IN  ||  vInters[j].first == ILBT_OUT  ||  vInters[j].first == ILBT_TOUCH)) {
-            vInters.erase( vInters.begin() + j) ;
-            if ( i > 0)
-               -- j ;
-            continue ;
-         }
-        // se una fine di coincidenza e l'altra inizio di coincidenza, elimino entrambe
-         else if ( ( vInters[i].first == ILBT_TG_FIN  &&  vInters[j].first == ILBT_TG_INI)  ||
-                   ( vInters[i].first == ILBT_TG_INI  &&  vInters[j].first == ILBT_TG_FIN)) {
-            vInters.erase( vInters.begin() + j) ;
-            vInters.erase( vInters.begin() + i) ;
-            if ( i > 0)
-               -- j ;
-            continue ;
-         }
-      }
-     // passo alla successiva
-      ++ j ;
+      return true ;
    }
 
+  // Se sono due intersezioni con la linea giacente su una faccia del box
+   if ( ( abs( vtL.x) < EPS_ZERO  &&  ( abs( ptL.x - ptMin.x) < EPS_SMALL  ||  abs( ptL.x - ptMax.x) < EPS_SMALL))  ||
+        ( abs( vtL.y) < EPS_ZERO  &&  ( abs( ptL.y - ptMin.y) < EPS_SMALL  ||  abs( ptL.y - ptMax.y) < EPS_SMALL))  ||
+        ( abs( vtL.z) < EPS_ZERO  &&  ( abs( ptL.z - ptMin.z) < EPS_SMALL  ||  abs( ptL.z - ptMax.z) < EPS_SMALL))) {
+      if ( ! bFinite) {
+         vInters.emplace_back( ILBT_TG_INI, dU1) ;
+         vInters.emplace_back( ILBT_TG_FIN, dU2) ;
+      }
+      else {
+         if ( dU1 > dLen - EPS_SMALL)
+            vInters.emplace_back( ILBT_IN, dLen) ;
+         else if ( dU2 < EPS_SMALL)
+            vInters.emplace_back( ILBT_OUT, 0) ;
+         else {
+            vInters.emplace_back( ILBT_TG_INI, Clamp( dU1, 0., dLen)) ;
+            vInters.emplace_back( ILBT_TG_FIN, Clamp( dU2, 0., dLen)) ;
+         }
+      }
+      return true ;
+   }
+
+  // Altrimenti sono due intersezioni con attraversamento
+   if ( ! bFinite) {
+      vInters.emplace_back( ILBT_IN, dU1) ;
+      vInters.emplace_back( ILBT_OUT, dU2) ;
+   }
+   else {
+      if ( dU1 > dLen - EPS_SMALL)
+         vInters.emplace_back( ILBT_IN, dLen) ;
+      else if ( dU2 < EPS_SMALL)
+         vInters.emplace_back( ILBT_OUT, 0) ;
+      else {
+         vInters.emplace_back( ILBT_IN, Clamp( dU1, 0., dLen)) ;
+         vInters.emplace_back( ILBT_OUT, Clamp( dU2, 0., dLen)) ;
+      }
+   }
    return true ;
 }

IntersLineBox.h

@@ -0,0 +1,34 @@
+//----------------------------------------------------------------------------
+//  EgalTech 2022-2022
+//----------------------------------------------------------------------------
+// File : IntersLineBox.h         Data : 08.01.22         Versione : 2.4a3
+// Contenuto : Dichiarazione delle funzioni base per intersezione linea/box.
+//
+//
+//
+// Modifiche : 08.01.22 DS Creazione modulo.
+//
+//
+//----------------------------------------------------------------------------
+
+#pragma once
+
+#include "/EgtDev/Include/EGkPoint3d.h"
+
+//----------------------------------------------------------------------------
+// Linea e box allineato agli assi sono nel medesimo riferimento.
+// Con intersezione viene restituito true e i parametri in dU1 e dU2.
+//----------------------------------------------------------------------------
+bool
+IntersLineBox( const Point3d& ptL, const Vector3d& vtL,
+               const Point3d& ptMin, const Point3d& ptMax,
+               double& dU1, double& dU2) ;
+
+//----------------------------------------------------------------------------
+inline bool
+TestIntersLineBox( const Point3d& ptL, const Vector3d& vtL,
+                   const Point3d& ptMin, const Point3d& ptMax)
+{
+   double dU1, dU2 ;
+   return IntersLineBox( ptL, vtL, ptMin, ptMax, dU1, dU2) ;
+}

IntersLineSurfStd.cpp

@@ -1,7 +1,7 @@
 //----------------------------------------------------------------------------
-//  EgalTech 2017-2018
+//  EgalTech 2017-2022
 //----------------------------------------------------------------------------
-// File : IntersLineSurfStd.cpp      Data : 03.12.18        Versione : 1.9l1
+// File : IntersLineSurfStd.cpp      Data : 08.01.22        Versione : 2.4a3
 // Contenuto : Implementazione delle funzioni di intersezione
 //             componente lineare e superficie standard.
 //
@@ -13,6 +13,7 @@
 
 #include "stdafx.h"
 #include "IntersLineSurfStd.h"
+#include "IntersLineBox.h"
 #include "/EgtDev/Include/EGkIntersLinePlane.h"
 #include "/EgtDev/Include/EGkFrame3d.h"
 #include "/EgtDev/Include/ENkPolynomialRoots.h"
@@ -21,60 +22,6 @@
 
 using namespace std ;
 
-//----------------------------------------------------------------------------
-// Punti e vettore devono essere espressi nel medesimo sistema di riferimento.
-// Il box è allineato con gli assi di questo sistema di riferimento.
-// In caso di intersezione viene restituito true e i parametri in dU1 e dU2.
-//----------------------------------------------------------------------------
-bool
-IntersLineBox( const Point3d& ptP, const Vector3d& vtV,
-               const Point3d& ptMin, const Point3d& ptMax, double& dU1, double& dU2)
-{
-  // Casi di intersezione impossibile
-   if ( abs( vtV.x) < EPS_ZERO  &&  ( ptP.x < ptMin.x - EPS_SMALL  ||  ptP.x > ptMax.x + EPS_SMALL))
-      return false ;
-   if ( abs( vtV.y) < EPS_ZERO  &&  ( ptP.y < ptMin.y - EPS_SMALL  ||  ptP.y > ptMax.y + EPS_SMALL))
-      return false ;
-   if ( abs( vtV.z) < EPS_ZERO  &&  ( ptP.z < ptMin.z - EPS_SMALL  ||  ptP.z > ptMax.z + EPS_SMALL))
-      return false ;
-
-  // Inizializzo parametri intersezioni
-   dU1 = -INFINITO ;
-   dU2 = INFINITO ;
-
-  // Confronto con piani YZ (perpendicolari ad asse X)
-   if ( vtV.x > EPS_ZERO) {
-      dU1 = max( dU1, ( ptMin.x - ptP.x) / vtV.x) ;
-      dU2 = min( dU2, ( ptMax.x - ptP.x) / vtV.x) ;
-   }
-   else if ( vtV.x < -EPS_ZERO) {
-      dU1 = max( dU1, ( ptMax.x - ptP.x) / vtV.x) ;
-      dU2 = min( dU2, ( ptMin.x - ptP.x) / vtV.x) ;
-   }
-
-  // Confronto con piani ZX (perpendicolari ad asse Y)
-   if ( vtV.y > EPS_ZERO) {
-      dU1 = max( dU1, ( ptMin.y - ptP.y) / vtV.y) ;
-      dU2 = min( dU2, ( ptMax.y - ptP.y) / vtV.y) ;
-   }
-   else if ( vtV.y < -EPS_ZERO) {
-      dU1 = max( dU1, ( ptMax.y - ptP.y) / vtV.y) ;
-      dU2 = min( dU2, ( ptMin.y - ptP.y) / vtV.y) ;
-   }
-
-  // Confronto con piani XY (perpendicolari ad asse Z)
-   if ( vtV.z > EPS_ZERO) {
-      dU1 = max( dU1, ( ptMin.z - ptP.z) / vtV.z) ;
-      dU2 = min( dU2, ( ptMax.z - ptP.z) / vtV.z) ;
-   }
-   else if ( vtV.z < -EPS_ZERO) {
-      dU1 = max( dU1, ( ptMax.z - ptP.z) / vtV.z) ;
-      dU2 = min( dU2, ( ptMin.z - ptP.z) / vtV.z) ;
-   }
-
-   return ( dU2 >= dU1) ;
-}
-
 //----------------------------------------------------------------------------
 int
 LineDisc( const Point3d& ptPLine, const Vector3d& vtVLine,

IntersLineSurfStd.h

@@ -34,11 +34,6 @@
 // Costanti tipologia di componente lineare
 enum LinType { Line = 0, Ray = 1, Segment = 2} ;
 
-//----------------------------------------------------------------------------
-bool IntersLineBox( const Point3d& ptP, const Vector3d& vtV,
-                    const Point3d& ptMin, const Point3d& ptMax, double& dU1, double& dU2) ;
-
-
 // Costanti tipologia di intersezioni fra un componente lineare e un disco
 enum LinCompDiscIntersType { D_ERROR_INT = -1, D_NO_INTERS = 0, D_BOUNDARY_INT_LINE_NOT_IN_PLANE = 1,
                              D_INNER_INT_LINE_NOT_IN_PLANE = 2, D_ONE_INT_LINE_ON_PLANE = 3,

IntersLineSurfTm.cpp

@@ -1,7 +1,7 @@
 //----------------------------------------------------------------------------
-//  EgalTech 2015-2019
+//  EgalTech 2015-2022
 //----------------------------------------------------------------------------
-// File : IntersLineSurfTm.cpp      Data : 09.03.15         Versione : 2.1b6
+// File : IntersLineSurfTm.cpp      Data : 08.01.22         Versione : 2.4a3
 // Contenuto : Implementazione della intersezione linea/superficie trimesh.
 //
 //
@@ -14,8 +14,7 @@
 //--------------------------- Include ----------------------------------------
 #include "stdafx.h"
 #include "ProjPlane.h"
-#include "IntersLineSurfStd.h"
-#include "/EgtDev/Include/EGkSurfTriMesh.h"
+#include "IntersLineBox.h"
 #include "/EgtDev/Include/EGkIntersLineTria.h"
 #include "/EgtDev/Include/EGkIntersLineSurfTm.h"
 

IntersPlaneBox.cpp

@@ -1,8 +1,8 @@
 //----------------------------------------------------------------------------
 //  EgalTech 2020-2020
 //----------------------------------------------------------------------------
-// File : IntersLineBox.cpp       Data : 07.05.20         Versione : 2.2e1
-// Contenuto : Implementazione della intersezione linea/box.
+// File : IntersPlaneBox.cpp       Data : 07.05.20         Versione : 2.2e1
+// Contenuto : Implementazione della intersezione piano/box.
 //
 //
 //
@@ -20,19 +20,47 @@
 using namespace std ;
 
 //----------------------------------------------------------------------------
+// Box e piano devono essere nello stesso riferimento. Il box deve essere axis aligned.
 bool
-IntersPlaneBox( const Point3d& ptOn, const Vector3d& vtN, const BBox3d& b3Box,
+TestIntersPlaneBox( const Plane3d& plPlane, const BBox3d& b3Box) 
+{
+  // Verifica del piano
+   if ( ! plPlane.IsValid())
+      return false ;
+  // Centro e estensione del box
+   Point3d ptCen ;
+   Vector3d vtExt ;
+   if ( ! b3Box.GetCenterExtent( ptCen, vtExt))
+      return false ;
+  // Distanza del centro dal piano
+   double dDist = DistPointPlane( ptCen, plPlane) ;
+  // Proiezione dell'estensione sulla normale al piano
+   Vector3d vtN = plPlane.GetVersN() ;
+   double dProjExt = vtExt.x * abs( vtN.x) + vtExt.y * abs( vtN.y) + vtExt.z * abs( vtN.z) ;
+  // Confronto distanza del centro con estensione proiettata
+   return ( abs( dDist) < dProjExt + EPS_SMALL) ;
+}
+
+//----------------------------------------------------------------------------
+bool
+IntersPlaneBox( const Plane3d& plPlane, const BBox3d& b3Box,
                 PNTVECTOR& vPnt, BIPNTVECTOR& vBpt, TRIA3DVECTOR& vTria)
 {
-  // Creo il piano
-   Plane3d plPlane ;
-   if ( ! plPlane.Set( ptOn, vtN))
+  // Verifico il piano
+   if ( ! plPlane.IsValid())
       return false ;
   // Recupero i dati del Box
    Point3d ptMin ;
    double dDimX, dDimY, dDimZ ;
    if ( ! b3Box.GetMinDim( ptMin, dDimX, dDimY, dDimZ))
       return false ;
+  // Pulisco vettori intersezioni
+   vPnt.clear() ;
+   vBpt.clear() ;
+   vTria.clear() ;
+  // Verifico se può esistere intersezione
+   if ( ! TestIntersPlaneBox( plPlane, b3Box))
+      return true ;
   // Contorno di base
    PolyLine PL ;
    PL.AddUPoint( 0, ptMin) ;

IntersSurfTmSurfTm.cpp

@@ -116,7 +116,7 @@ IntersSurfTmSurfTm( const ISurfTriMesh& Stm1, const ISurfTriMesh& Stm2,
                   }
                }
               // se altrimenti sovrapposizione
-               else if ( nRes == ITTT_OVERLAPS) {
+               else if ( nRes == ITTT_OVERLAPS  ||  nRes == ITTT_COUNTER_OVERLAPS) {
                   for ( const auto& Tria : vIttTria) {
                     // verifico se triangolo già inserito
                      bool bFound = false ;

IntersTriaTria.cpp

@@ -83,8 +83,37 @@ IntersTriaTria( const Triangle3d& trTria1, const Triangle3d& trTria2, Point3d& p
   // se i triangoli sono complanari
    if ( nResPP == IPPT_OVERLAPS  ||
         ((( nVertPos1 == 0  &&  nVertNeg1 == 0)  ||  ( nVertPos2 == 0  &&  nVertNeg2 == 0)) &&
-         ( trTria1.GetN() ^ trTria2.GetN()).SqLen() < 25 * SIN_EPS_ANG_SMALL * SIN_EPS_ANG_SMALL))
-      return IntersCoplanarTriaTria( trTria1, trTria2, vTria) ;
+         ( trTria1.GetN() ^ trTria2.GetN()).SqLen() < 25 * SIN_EPS_ANG_SMALL * SIN_EPS_ANG_SMALL)) {
+     // verifica per triangoli con normali controverse
+      bool bCounter = false ;
+      Triangle3d trMyTria2 = trTria2 ;
+      if ( trTria1.GetN() * trTria2.GetN() < 0) {
+         Point3d ptV0 = trMyTria2.GetP( 0) ;
+         trMyTria2.SetP( 0, trMyTria2.GetP( 1)) ;
+         trMyTria2.SetP( 1, ptV0) ;
+         trMyTria2.Validate() ;
+         bCounter = true ;
+      }
+     // verifica per triangoli coincidenti
+      bool bAreSameTria = ( ( AreSamePointExact( trTria1.GetP( 0), trMyTria2.GetP( 0))  &&
+                              AreSamePointExact( trTria1.GetP( 1), trMyTria2.GetP( 1))  &&
+                              AreSamePointExact( trTria1.GetP( 2), trMyTria2.GetP( 2)))  ||
+                            ( AreSamePointExact( trTria1.GetP( 0), trMyTria2.GetP( 1))  &&
+                              AreSamePointExact( trTria1.GetP( 1), trMyTria2.GetP( 2))  &&
+                              AreSamePointExact( trTria1.GetP( 2), trMyTria2.GetP( 0)))  ||
+                            ( AreSamePointExact( trTria1.GetP( 0), trMyTria2.GetP( 2))  &&
+                              AreSamePointExact( trTria1.GetP( 1), trMyTria2.GetP( 0))  &&
+                              AreSamePointExact( trTria1.GetP( 2), trMyTria2.GetP( 1))))  &&
+                          ( AreSameVectorExact( trTria1.GetN(), trMyTria2.GetN())) ;
+      if ( bAreSameTria) {
+         vTria.emplace_back( trTria1) ;
+         return ( bCounter ? ITTT_COUNTER_OVERLAPS : ITTT_OVERLAPS) ;
+      }
+      if ( ITTT_OVERLAPS == IntersCoplanarTriaTria( trTria1, trMyTria2, vTria)) {
+         return ( bCounter ? ITTT_COUNTER_OVERLAPS : ITTT_OVERLAPS) ;
+      }
+      return ITTT_NO ;
+   }
 
   // limito la linea di intersezione con il primo triangolo
    Point3d ptSt1, ptEn1 ;
@@ -112,7 +141,7 @@ IntersTriaTria( const Triangle3d& trTria1, const Triangle3d& trTria2, Point3d& p
 int
 IntersCoplanarTriaTria( const Triangle3d& trTria1, const Triangle3d& trTria2, TRIA3DVECTOR& vTria)
 {
-  // Se i tre vertici del secondo triangolo stanno tutti a destra di almeno un lato del primo, non c'è intersezione
+  // Se i tre vertici del secondo triangolo stanno tutti a destra di almeno un lato del primo, non c'� intersezione
    for ( int i = 0 ; i < 3 ; ++ i) {
       Vector3d vtSide = trTria1.GetP( ( i + 1) % 3) - trTria1.GetP( i) ;
       Vector3d vtSN = vtSide ^ trTria1.GetN() ;
@@ -122,7 +151,7 @@ IntersCoplanarTriaTria( const Triangle3d& trTria1, const Triangle3d& trTria2, TR
            ( trTria2.GetP( 2) - trTria1.GetP( i)) * vtSN > - EPS_TRIA_H)
          return ITTT_NO ;
    }
-  // Se i tre vertici del primo triangolo stanno tutti a destra di almeno un lato del secondo, non c'è intersezione
+  // Se i tre vertici del primo triangolo stanno tutti a destra di almeno un lato del secondo, non c'� intersezione
    for ( int i = 0 ; i < 3 ; ++ i) {
       Vector3d vtSide = trTria2.GetP( ( i + 1) % 3) - trTria2.GetP( i) ;
       Vector3d vtSN = vtSide ^ trTria2.GetN() ;
@@ -205,7 +234,7 @@ int
 FindTriaTriaIntersType( const Triangle3d& trTria1, const Triangle3d& trTria2, const Point3d& ptLineSt, const Vector3d& vtLineDir,
                         double dStU1, double dEnU1, double dStU2, double dEnU2, int nRes1, int nRes2, double& dIntStU, double& dIntEnU)
 {
-  // Controllo su validità input
+  // Controllo su validit� input
    if ( ! ( trTria1.IsValid()  &&  trTria2.IsValid()  &&  vtLineDir.IsNormalized()))
       return ITTT_NO ;
   // Casi

SurfTriMesh.cpp

@@ -18,10 +18,10 @@
 #include "GeoObjFactory.h"
 #include "NgeWriter.h"
 #include "NgeReader.h"
+#include "SurfFlatRegion.h"
 #include "DistPointLine.h"
 #include "Triangulate.h"
 #include "GeoConst.h"
-#include "SurfFlatRegion.h"
 #include "/EgtDev/Include/EGkIntersLinePlane.h"
 #include "/EgtDev/Include/EGkPointGrid3d.h"
 #include "/EgtDev/Include/EGkPolyLine.h"
@@ -239,7 +239,7 @@ SurfTriMesh::RemoveTriangle( int nId)
      // indice vertice
       int nV = m_vTria[nId].nIdVert[i] ;
      // se vertice non c'è passo al prossimo
-      if ( nV < 0  ||  nV >= int( m_vTria.size()))
+      if ( nV < 0  ||  nV >= int( m_vVert.size()))
          continue ;
       if ( m_vVert[nV].nIdTria == nId) {
          int nAdjP = m_vTria[nId].nIdAdjac[i] ;
@@ -1695,8 +1695,10 @@ SurfTriMesh::TestSealing( void)
         // verifico le adiacenze
          if ( m_vTria[i].nIdAdjac[0] == SVT_NULL  ||
               m_vTria[i].nIdAdjac[1] == SVT_NULL  ||
-              m_vTria[i].nIdAdjac[2] == SVT_NULL)
+              m_vTria[i].nIdAdjac[2] == SVT_NULL) {
             bClosed = false ;
+            break ;  
+         }
       }
    }
   // aggiorno la chiusura della superficie
@@ -1716,7 +1718,7 @@ SurfTriMesh::AdjustTopology( void)
       return true ;
    }
   // dichiaro sfaccettatura da ricalcolare
-   ResetFaceting() ;
+   m_bFaceted = false ;
   // invalido calcolo connessione
    m_nParts = - 1 ;
   // verifica indici
@@ -2982,7 +2984,7 @@ SurfTriMesh::DoSewing( const ISurfTriMesh& stmOther, const Frame3d& frOther, dou
    }
    int nVIdSize = int( vVId.size()) ;
 
-  // inserisco i triangoli dell'altra trimesh
+  // aggiungo i triangoli dell'altra trimesh
    for ( int nOtId = 0 ; nOtId < pOther->GetTriangleSize() ; ++ nOtId) {
      // recupero gli indici dei vertici del triangolo
       int vOId[3] ;
@@ -3399,336 +3401,6 @@ SurfTriMesh::Invert( void)
    return true ;
 }
 
-//----------------------------------------------------------------------------
-bool
-SurfTriMesh::CutByTriangles( const Plane3d& plPlane, bool bSaveOnEq, bool& bModif)
-{
-  // la superficie deve essere validata
-   if ( m_nStatus != OK)
-      return false ;
-
-  // classifico i vertici rispetto al piano
-   double dTol = max( min( 0.5 * m_dLinTol, 10 * EPS_SMALL), EPS_SMALL) ;
-   for ( int i = 0 ; i < GetVertexSize() ; ++ i) {
-      double dDist = DistPointPlane( m_vVert[i].ptP, plPlane) ;
-      if ( abs( dDist) < dTol) {
-         m_vVert[i].nTemp = 0 ;
-         if ( abs( dDist) > EPS_SMALL)
-            m_vVert[i].ptP -= plPlane.GetVersN() * dDist ;
-      }
-      else if ( dDist > 0)
-         m_vVert[i].nTemp = +1 ;
-      else
-         m_vVert[i].nTemp = -1 ;
-   }
-
-  // aggiorno timestamp dei triangoli
-   for ( auto& Tria : m_vTria)
-      Tria.nTemp = m_nTimeStamp ;
-   ++ m_nTimeStamp ;
-
-  // sistemo i triangoli (eventualmente li elimino)
-   for ( int i = 0 ; i < GetTriangleSize() ; ++ i) {
-     // salto i triangoli cancellati e quelli aggiunti
-      if ( m_vTria[i].nIdVert[0] == SVT_DEL  ||  m_vTria[i].nTemp == m_nTimeStamp)
-         continue ;
-     // se flag abilita e giace sul piano ed equiverso, lo salvo
-      if ( bSaveOnEq  &&
-           m_vVert[ m_vTria[i].nIdVert[0]].nTemp == 0  &&
-           m_vVert[ m_vTria[i].nIdVert[1]].nTemp == 0  &&
-           m_vVert[ m_vTria[i].nIdVert[2]].nTemp == 0  &&
-           AreSameVectorApprox( m_vTria[i].vtN, plPlane.GetVersN()))
-         ;
-     // se giace sul piano o dalla parte esterna del piano, lo cancello
-      else if ( m_vVert[ m_vTria[i].nIdVert[0]].nTemp != -1  &&
-                m_vVert[ m_vTria[i].nIdVert[1]].nTemp != -1  &&
-                m_vVert[ m_vTria[i].nIdVert[2]].nTemp != -1) {
-         bModif = true ;
-         RemoveTriangle( i) ;
-      }
-     // se giace dalla parte interna con al massimo un lato sul piano, lo conservo
-      else if ( m_vVert[ m_vTria[i].nIdVert[0]].nTemp != 1  &&
-                m_vVert[ m_vTria[i].nIdVert[1]].nTemp != 1  &&
-                m_vVert[ m_vTria[i].nIdVert[2]].nTemp != 1)
-         ;
-     // altrimenti attraversa il piano, devo modificarlo
-      else {
-         bModif = true ;
-         int nDisc =  m_vVert[ m_vTria[i].nIdVert[0]].nTemp +
-                      m_vVert[ m_vTria[i].nIdVert[1]].nTemp +
-                      m_vVert[ m_vTria[i].nIdVert[2]].nTemp ;
-        // se ha un vertice all'interno e uno sul piano -> 1 nuovo vertice e 1 nuovo triangolo
-         if ( nDisc == 0) {
-           // classifico i vertici
-            int nVertOn, nVertIn, nVertOut ;
-            if ( m_vVert[ m_vTria[i].nIdVert[0]].nTemp == 0) {
-               nVertOn = 0 ;
-               nVertIn = 1 ;
-               nVertOut = 2 ;
-               if ( m_vVert[ m_vTria[i].nIdVert[1]].nTemp > 0)
-                  swap( nVertIn, nVertOut) ;
-            }
-            else if ( m_vVert[ m_vTria[i].nIdVert[1]].nTemp == 0) {
-               nVertOn = 1 ;
-               nVertIn = 2 ;
-               nVertOut = 0 ;
-               if ( m_vVert[ m_vTria[i].nIdVert[2]].nTemp > 0)
-                  swap( nVertIn, nVertOut) ;
-            }
-            else {
-               nVertOn = 2 ;
-               nVertIn = 0 ;
-               nVertOut = 1 ;
-               if ( m_vVert[ m_vTria[i].nIdVert[0]].nTemp > 0)
-                  swap( nVertIn, nVertOut) ;
-            }
-           // calcolo il punto di intersezione del lato che attraversa il piano
-            double dDistIn = - DistPointPlane( m_vVert[m_vTria[i].nIdVert[nVertIn]].ptP, plPlane) ;
-            double dDistOut = DistPointPlane( m_vVert[m_vTria[i].nIdVert[nVertOut]].ptP, plPlane) ;
-            Point3d ptInt = Media( m_vVert[m_vTria[i].nIdVert[nVertIn]].ptP,
-                                   m_vVert[m_vTria[i].nIdVert[nVertOut]].ptP,
-                                   dDistIn / ( dDistIn + dDistOut)) ;
-           // inserisco il nuovo vertice
-            int nNewV = AddVertex( ptInt) ;  
-           // inserisco il nuovo triangolo
-            int nIdVert[3] = { m_vTria[i].nIdVert[0], m_vTria[i].nIdVert[1], m_vTria[i].nIdVert[2]} ;
-            nIdVert[nVertOut] = nNewV ;
-            int nIdTria = AddTriangle( nIdVert) ;
-            if ( nIdTria == SVT_NULL)
-               return false ;
-            if ( nIdTria != SVT_DEL)
-               m_vTria[nIdTria].nTemp = m_nTimeStamp ;
-           // cancello il vecchio triangolo
-            RemoveTriangle( i) ;                      
-         }
-        // se ha un vertice all'interno e due all'esterno -> 2 nuovi vertici e 1 nuovo triangolo
-         else if ( nDisc == 1) {
-           // classifico i vertici
-            int nVertIn, nVertOut1, nVertOut2 ;
-            if ( m_vVert[ m_vTria[i].nIdVert[0]].nTemp == -1) {
-               nVertIn = 0 ;
-               nVertOut1 = 1 ;
-               nVertOut2 = 2 ;
-            }
-            else if ( m_vVert[ m_vTria[i].nIdVert[1]].nTemp == -1) {
-               nVertIn = 1 ;
-               nVertOut1 = 2 ;
-               nVertOut2 = 0 ;
-            }
-            else {
-               nVertIn = 2 ;
-               nVertOut1 = 0 ;
-               nVertOut2 = 1 ;
-            }
-           // calcolo i punti di intersezione dei lati che attraversano il piano
-            double dDistIn = - DistPointPlane( m_vVert[m_vTria[i].nIdVert[nVertIn]].ptP, plPlane) ;
-            double dDistOut1 = DistPointPlane( m_vVert[m_vTria[i].nIdVert[nVertOut1]].ptP, plPlane) ;
-            double dDistOut2 = DistPointPlane( m_vVert[m_vTria[i].nIdVert[nVertOut2]].ptP, plPlane) ;
-            Point3d ptInt1 = Media( m_vVert[m_vTria[i].nIdVert[nVertIn]].ptP,
-                                    m_vVert[m_vTria[i].nIdVert[nVertOut1]].ptP,
-                                    dDistIn / ( dDistIn + dDistOut1)) ;
-            Point3d ptInt2 = Media( m_vVert[m_vTria[i].nIdVert[nVertIn]].ptP,
-                                    m_vVert[m_vTria[i].nIdVert[nVertOut2]].ptP,
-                                    dDistIn / ( dDistIn + dDistOut2)) ;
-           // inserisco i nuovi vertici
-            int nNewV1 = AddVertex( ptInt1) ;  
-            int nNewV2 = AddVertex( ptInt2) ;  
-           // inserisco il nuovo triangolo
-            int nIdVert[3] = { m_vTria[i].nIdVert[0], m_vTria[i].nIdVert[1], m_vTria[i].nIdVert[2]} ;
-            nIdVert[nVertOut1] = nNewV1 ;
-            nIdVert[nVertOut2] = nNewV2 ;
-            int nIdTria = AddTriangle( nIdVert) ;
-            if ( nIdTria == SVT_NULL)
-               return false ;
-            if ( nIdTria != SVT_DEL)
-               m_vTria[nIdTria].nTemp = m_nTimeStamp ;
-           // cancello il vecchio triangolo
-            RemoveTriangle( i) ;                      
-         }
-        // altrimenti ha due vertici all'interno e uno all'esterno -> 2 nuovi vertici e 2 nuovi triangoli
-         else {
-           // classifico i vertici
-            int nVertIn1, nVertIn2, nVertOut ;
-            if ( m_vVert[ m_vTria[i].nIdVert[0]].nTemp == 1) {
-               nVertOut = 0 ;
-               nVertIn1 = 1 ;
-               nVertIn2 = 2 ;
-            }
-            else if ( m_vVert[ m_vTria[i].nIdVert[1]].nTemp == 1) {
-               nVertOut = 1 ;
-               nVertIn1 = 2 ;
-               nVertIn2 = 0 ;
-            }
-            else {
-               nVertOut = 2 ;
-               nVertIn1 = 0 ;
-               nVertIn2 = 1 ;
-            }
-           // calcolo i punti di intersezione dei lati che attraversano il piano
-            double dDistIn1 = - DistPointPlane( m_vVert[m_vTria[i].nIdVert[nVertIn1]].ptP, plPlane) ;
-            double dDistIn2 = - DistPointPlane( m_vVert[m_vTria[i].nIdVert[nVertIn2]].ptP, plPlane) ;
-            double dDistOut = DistPointPlane( m_vVert[m_vTria[i].nIdVert[nVertOut]].ptP, plPlane) ;
-            Point3d ptInt1 = Media( m_vVert[m_vTria[i].nIdVert[nVertIn1]].ptP,
-                                    m_vVert[m_vTria[i].nIdVert[nVertOut]].ptP,
-                                    dDistIn1 / ( dDistIn1 + dDistOut)) ;
-            Point3d ptInt2 = Media( m_vVert[m_vTria[i].nIdVert[nVertIn2]].ptP,
-                                    m_vVert[m_vTria[i].nIdVert[nVertOut]].ptP,
-                                    dDistIn2 / ( dDistIn2 + dDistOut)) ;
-           // inserisco i nuovi vertici
-            int nNewV1 = AddVertex( ptInt1) ;  
-            int nNewV2 = AddVertex( ptInt2) ;  
-           // inserisco i nuovi triangoli
-            int nIdVert1[3] = { m_vTria[i].nIdVert[0], m_vTria[i].nIdVert[1], m_vTria[i].nIdVert[2]} ;
-            int nIdVert2[3] = { m_vTria[i].nIdVert[0], m_vTria[i].nIdVert[1], m_vTria[i].nIdVert[2]} ;
-            double dSqDist1 = SqDist( m_vVert[m_vTria[i].nIdVert[nVertIn1]].ptP, ptInt2) ;
-            double dSqDist2 = SqDist( m_vVert[m_vTria[i].nIdVert[nVertIn2]].ptP, ptInt1) ;
-            if ( dSqDist1 <= dSqDist2) {
-               nIdVert1[nVertOut] = nNewV1 ;
-               nIdVert1[nVertIn2] = nNewV2 ;
-               nIdVert2[nVertOut] = nNewV2 ;
-            }
-            else {
-               nIdVert1[nVertOut] = nNewV1 ;
-               nIdVert2[nVertOut] = nNewV2 ;
-               nIdVert2[nVertIn1] = nNewV1 ;
-            }
-            int nIdTria1 = AddTriangle( nIdVert1) ;
-            int nIdTria2 = AddTriangle( nIdVert2) ;
-            if ( nIdTria1 == SVT_NULL  ||  nIdTria2 == SVT_NULL)
-               return false ;
-            if ( nIdTria1 != SVT_DEL)
-               m_vTria[nIdTria1].nTemp = m_nTimeStamp ;
-            if ( nIdTria2 != SVT_DEL)
-               m_vTria[nIdTria2].nTemp = m_nTimeStamp ;
-           // cancello il vecchio triangolo
-            RemoveTriangle( i) ;                      
-         }
-      }
-   }
-
-   return true ;
-}
-
-//----------------------------------------------------------------------------
-bool
-SurfTriMesh::CutByFacets( const Plane3d& plPlane, bool bSaveOnEq, bool& bModif)
-{
-  // Setto posizione triangoli non definita
-   for ( int nT = 0 ; nT < GetTriangleSize() ; ++ nT)
-      m_vTria[nT].nTempPart = 0 ;
-
-   INTERSCHAINMAP IntersLineMap ;
-
-  // Ciclo su tutte le facce.
-   for ( int nF = 0 ; nF < GetFacetCount() ; ++ nF) {
-     // Dati della faccia
-      POLYLINEVECTOR vLoopVec ;
-      GetFacetLoops( nF, vLoopVec) ;
-      if ( vLoopVec.empty())
-         continue ;
-     // Verifico la posizione del loop esterno
-      int nIntType = VerifyLoopPlane( vLoopVec[0], plPlane) ;
-     // Se interseca il piano di taglio, calcolo la divisione della faccia 
-      LineFacetClassVector IntersLinePart ;
-      if ( nIntType == FPI_CUT) {
-         PtrOwner<SurfFlatRegion> pReg( GetBasicSurfFlatRegion( GetSurfFlatRegionFromPolyLineVector( vLoopVec))) ;
-         if ( IsNull( pReg)  ||  ! pReg->IsValid())
-            return false ;
-         nIntType = IntersFacetPlane( *pReg, plPlane, IntersLinePart) ;
-      }
-     // Gestione dei risultati
-      if ( nIntType == FPI_CUT) {
-         for ( int nPart = 0 ; nPart < int( IntersLinePart.size()) ; ++ nPart) {
-           // Salvo intersezione per la faccia.
-            auto it = IntersLineMap.find( nF) ;
-            if ( it != IntersLineMap.end()) {
-               it->second.emplace_back( IntersInnSeg( IntersLinePart[nPart].ptSt, IntersLinePart[nPart].ptEn)) ;
-            }
-            else {
-               IntersLineMap.emplace( nF, IntersInnChain( 1, IntersInnSeg( IntersLinePart[nPart].ptSt, IntersLinePart[nPart].ptEn))) ;
-            }
-         }
-      }
-      else if ( nIntType == FPI_ON) { 
-         INTVECTOR vT ;
-         GetAllTriaInFacet( nF, vT) ;
-         Vector3d vtNf ;
-         GetFacetNormal( nF, vtNf) ;
-         for ( auto& nT : vT)
-            m_vTria[nT].nTempPart = ( vtNf * plPlane.GetVersN() > 0 ? 2 : -2) ;
-      }
-      else if ( nIntType == FPI_IN) {
-         INTVECTOR vT ;
-         GetAllTriaInFacet( nF, vT) ;
-         for ( auto& nT : vT)
-            m_vTria[nT].nTempPart = 1 ;
-      }
-      else if ( nIntType == FPI_OUT) {
-         INTVECTOR vT ;
-         GetAllTriaInFacet( nF, vT) ;
-         for ( auto& nT : vT)
-            m_vTria[nT].nTempPart = -1 ;
-      }
-      else
-         return false ;
-   }
-
-  // Divido le facce.
-   PieceMap NewFacet ;
-   SplitFacet( IntersLineMap, NewFacet) ;
-   INTERSEDGEMAP EdgeInterLineMap, EdgeEdgeLineMap ;
-   RetriangulateFacetPieces( NewFacet, EdgeInterLineMap, EdgeEdgeLineMap) ;
-
-  // Elimino i triangoli superflui
-   for ( int nT = 0 ; nT < GetTriangleSize() ; ++ nT) {
-      if (  m_vTria[nT].nTempPart == 0  ||
-            m_vTria[nT].nTempPart == -1  ||
-            m_vTria[nT].nTempPart == -2  ||
-            ( ! bSaveOnEq  &&  m_vTria[nT].nTempPart == 2)) {           
-         RemoveTriangle( nT) ;
-         bModif = true ;
-      }
-   }
-
-   return true ;
-}
-
-//----------------------------------------------------------------------------
-bool
-SurfTriMesh::Cut( const Plane3d& plPlane, bool bSaveOnEq)
-{
-  // la superficie deve essere validata
-   if ( m_nStatus != OK)
-      return false ;
-
-  // recupero il numero originale di triangoli e di facce
-   int nTriaOriCnt = GetTriangleCount() ;
-   int nFacetOriCnt = GetFacetCount() ;
-
-  // eseguo il taglio con il metodo delle faccette o dei triangoli (per ora solo coi triangoli)
-   bool bModif = false ;
-   if ( true  ||  ! CutByFacets( plPlane, bSaveOnEq, bModif)) {
-      bModif = false ;
-      if ( ! CutByTriangles( plPlane, bSaveOnEq, bModif))
-         return false ;
-   }
-
-  // se effettuate modifiche
-   if ( bModif) {
-     // aggiorno tutto
-      if ( ! AdjustVertices()  ||  ! DoCompacting())
-         return false ;
-   }
-
-  // se superficie originale a facce, cerco di semplificarle in ogni caso
-   if ( nFacetOriCnt < 200  ||  double( nTriaOriCnt) / nFacetOriCnt > 4) {
-      if ( ! SimplifyFacets( 500.0))
-         LOG_ERROR( GetEGkLogger(), "Error in SimplifyFacets of Stm::Cut")
-   }
-
-   return true ;
-}
-
 //----------------------------------------------------------------------------
 void
 SurfTriMesh::ResetHashGrids3d( void) const 

SurfTriMesh.h

@@ -94,30 +94,11 @@ typedef std::unordered_map< int, Chain> CHAINMAP ;
 typedef std::unordered_map< int, TRIA3DVECTOR> TRIA3DVECTORMAP ;
 
 //----------------------------------------------------------------------------
-// Struttura per la ricerca delle T-Junction
-struct TEdgeId {
-   int nTriaId ;
-   int nEdge ;
-} ;
-// Overloading dell'operatore < per TEdgeSet
-inline bool operator<( const TEdgeId& TEdgeA, const TEdgeId& TEdgeB) {
-   if ( TEdgeA.nTriaId != TEdgeB.nTriaId)
-      return ( TEdgeA.nTriaId < TEdgeB.nTriaId) ;
-   return ( TEdgeA.nEdge < TEdgeB.nEdge) ;
-}
-// Loop con T-Junction
-typedef std::vector<TEdgeId> TJuncLoop ;
-// Vettore di loop con T-Junction
-typedef std::vector<TJuncLoop> TJuncLoopVec ;
-// Set di TEdgeId
-typedef std::set<TEdgeId> TJEdgeSet ;
-
-// Definizione strutture intersezione linea-faccia
-//----------------------------------------------------------------------------
+// Struttura intersezione linea-faccia
 struct LineFacetClass {
    Point3d ptSt, ptEn ;
    int nTypeA, nTypeB ;
-   LineFacetClass() {
+   LineFacetClass( void) {
       nTypeA = 0 ;
       nTypeB = 0 ;  
    }
@@ -134,25 +115,24 @@ typedef std::vector<LineFacetClass> LineFacetClassVector ;
 enum FacetPlaneIntersType { FPI_ERROR = 0, FPI_CUT = 1, FPI_IN = 2, FPI_OUT = 3, FPI_ON = 4 } ;
 
 //----------------------------------------------------------------------------
-// Definizione strutture e contenitori
-// Contatti interno-interno
+// Strutture e contenitori contatti interno-interno
 struct IntersInnSeg {
    Point3d ptSt ;
    Point3d ptEn ;
    Vector3d vtOuter ;
-   IntersInnSeg(void) {
+   IntersInnSeg( void) {
       ;
    }
-   IntersInnSeg(const Point3d& ptS, const Point3d& ptE) {
-      ptSt = ptS;
-      ptEn = ptE;
+   IntersInnSeg( const Point3d& ptS, const Point3d& ptE) {
+      ptSt = ptS ;
+      ptEn = ptE ;
    }
-   IntersInnSeg(const Point3d& ptS, const Point3d& ptE, const Vector3d& vtO) {
-      ptSt = ptS;
-      ptEn = ptE;
-      vtOuter = vtO;
+   IntersInnSeg( const Point3d& ptS, const Point3d& ptE, const Vector3d& vtO) {
+      ptSt = ptS ;
+      ptEn = ptE ;
+      vtOuter = vtO ;
    }
-};
+} ;
 typedef std::vector<IntersInnSeg> IntersInnChain ;
 typedef std::vector<IntersInnChain> INNCHAINVECTOR ;
 typedef std::unordered_map<int, IntersInnChain> INTERSCHAINMAP ;
@@ -164,7 +144,7 @@ struct IntersEdge {
    IntersEdge( const Point3d& ptS, const Point3d& ptE, const INTVECTOR& vOFI) {
       ptSt = ptS;
       ptEn = ptE;
-      vOthFacetIndex = vOFI;
+      vOthFacetIndex = vOFI ;
    }
 } ;
 typedef std::vector<IntersEdge> IntersEdgeVec ;
@@ -305,6 +285,7 @@ class SurfTriMesh : public ISurfTriMesh, public IGeoObjRW
       bool GetSurfClassification( const ISurfTriMesh& ClassifierSurf,
                      INTVECTOR& vTriaIn, INTVECTOR& vTriaOut, INTVECTOR& vTriaOnP, INTVECTOR& vTriaOnM, INTVECTOR& vTriaIndef) override ;
       bool CutWithOtherSurf( const ISurfTriMesh& CutterSurf, bool bInVsOut, bool bSaveOnEq) override ;
+      bool Repair( double dMaxEdgeLen = MAX_EDGE_LEN_STD) override ;
       bool GetAllTriaOverlapBox( const BBox3d& b3Box, INTVECTOR& vT) const override ;
       const BBox3d& GetAllTriaBox( void) const override ;
       int GetPartCount( void) const override ;
@@ -370,9 +351,9 @@ class SurfTriMesh : public ISurfTriMesh, public IGeoObjRW
       bool MarchOneTria( int& nT, int& nV, int nTimeStamp, PolyLine& PL, bool& bEnd) const ;
       bool VerifyPolylinesForTwoCurves( const PolyLine&  PL1, const PolyLine&  PL2) const ;
       bool AddBiTriangle( const int nIdVert[4]) ;
-      bool ResetFaceting( void) ;
       bool VerifyFaceting( void) const ;
       bool UpdateFaceting( void) ;
+      bool UpdateOneFace( int nFacet, int nT) ;
       bool UpdateTriaFaceting( int nRefT, int nFacet, const Plane3d& plPlane, int nT) ;
       bool SetFacet( int nInd, int nT) ; 
       bool VerifyAdjacTriaFacet( int nT, INTVECTOR& vT) const ;
@@ -381,6 +362,7 @@ class SurfTriMesh : public ISurfTriMesh, public IGeoObjRW
       void ResetHashGrids3d( void) const ;
       bool VerifyHashGrids3d( void) const ;
       bool VerifyConnection( void) const ;
+      bool CutTriangleByPlane( int nTriaId, const Plane3d& plPlane, bool bSaveOnEq, bool& bModif) ;
       bool CutByTriangles( const Plane3d& plPlane, bool bSaveOnEq, bool& bModif) ;
       bool CutByFacets( const Plane3d& plPlane, bool bSaveOnEq, bool& bModif) ;
       bool DecomposeLoop( CHAINVECTOR& cvOpenChain, INTVECTOR& vnDegVec, PNTMATRIX& cvBoundClosedLoopVec, BOOLVECTOR& vbInOut) ;
@@ -397,7 +379,6 @@ class SurfTriMesh : public ISurfTriMesh, public IGeoObjRW
       bool RetriangulateFacetPieces( const PieceMap& NewFacet,
                                      const INTERSEDGEMAP& EdgeInterLineMap,
                                      const INTERSEDGEMAP& EdgeEdgeLineMap) ;
-      bool SimplifyFacets( double dMaxEdgeLen) ;
       bool EdgeInteriorContactManager( const SurfTriMesh& OthSurf,
                                        const INTERSCHAINMAP& InterInterLineMap,
                                        const INTERSEDGEMAP& EdgeInterLineMap) ;
@@ -411,6 +392,11 @@ class SurfTriMesh : public ISurfTriMesh, public IGeoObjRW
       bool RemoveTripleTriangles( void) ;
       bool ScanForTripleTriangles( bool& bModified) ;
       bool FlipTriangles( int nTA, int nTB) ;
+      bool SimplifyFacets( double dMaxEdgeLen = MAX_EDGE_LEN_STD, bool bForced = true) ;
+      bool AddChainToChain( const Chain& ChainToAdd, PNTVECTOR& OrigChain) ;
+      bool DistPointFacet( const Point3d& ptP, const POLYLINEVECTOR& vPolyVec, double& dPointFacetDist) ;
+      bool ChangeStart( const Point3d& ptNewStart, PNTVECTOR& Loop) ;
+      bool SplitAtPoint( const Point3d& ptStop, const PNTVECTOR& Loop, PNTVECTOR& Loop1, PNTVECTOR& Loop2) ;
       
    private :
       ObjGraphicsMgr m_OGrMgr ;        // gestore grafica dell'oggetto

SurfTriMeshCuts.cpp

@@ -0,0 +1,874 @@
+//----------------------------------------------------------------------------
+//  EgalTech 2019-2021
+//----------------------------------------------------------------------------
+// File : SurfTriMeshCuts.cpp     Data : 06.11.21     Versione : 2.3k3
+// Contenuto : Implementazione delle funzioni di taglio per SurfFTrimesh.
+//
+//
+//
+// Modifiche : 06.11.21 DS Creazione modulo.
+//
+//
+//----------------------------------------------------------------------------
+
+#include "stdafx.h"
+#include "CurveComposite.h"
+#include "SurfTriMesh.h"
+#include "SurfFlatRegion.h"
+#include "DistPointLine.h"
+#include "Triangulate.h"
+#include "GeoConst.h"
+#include "/EgtDev/Include/EgkSfrCreate.h"
+#include "/EgtDev/Include/EGkDistPointCurve.h"
+#include "/EgtDev/Include/EGkIntersLinePlane.h"
+#include "/EgtDev/Include/EGkIntersPlaneTria.h"
+
+using namespace std ;
+
+//----------------------------------------------------------------------------
+const double CUT_SCALE = 1024 ;
+
+//----------------------------------------------------------------------------
+bool
+SurfTriMesh::Cut( const Plane3d& plPlane, bool bSaveOnEq)
+{
+  // la superficie deve essere validata
+   if ( m_nStatus != OK)
+      return false ;
+
+  // recupero il numero originale di triangoli e di facce
+   int nTriaOriCnt = GetTriangleCount() ;
+   int nFacetOriCnt = GetFacetCount() ;
+
+  // eseguo il taglio con il metodo delle faccette o dei triangoli (per ora solo coi triangoli)
+   bool bModif = false ;
+   if ( true  ||  ! CutByFacets( plPlane, bSaveOnEq, bModif)) {
+      bModif = false ;
+      if ( ! CutByTriangles( plPlane, bSaveOnEq, bModif))
+         return false ;
+   }
+
+  // se effettuate modifiche
+   if ( bModif) {
+     // aggiorno tutto
+      if ( ! AdjustVertices()  ||  ! DoCompacting())
+         return false ;
+   }
+
+  // se superficie originale a facce, cerco di semplificarle in ogni caso
+   if ( nFacetOriCnt < 200  ||  double( nTriaOriCnt) / nFacetOriCnt > 4) {
+      if ( ! SimplifyFacets())
+         LOG_ERROR( GetEGkLogger(), "Error in SimplifyFacets of Stm::Cut")
+   }
+
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+bool
+SurfTriMesh::CutByTriangles( const Plane3d& plPlane, bool bSaveOnEq, bool& bModif)
+{
+  // la superficie deve essere validata
+   if ( m_nStatus != OK)
+      return false ;
+
+  // classifico i vertici rispetto al piano
+   double dTol = max( min( 0.5 * m_dLinTol, 10 * EPS_SMALL), EPS_SMALL) ;
+   for ( int i = 0 ; i < GetVertexSize() ; ++ i) {
+      double dDist = DistPointPlane( m_vVert[i].ptP, plPlane) ;
+      if ( abs( dDist) < dTol) {
+         m_vVert[i].nTemp = 0 ;
+         if ( abs( dDist) > EPS_SMALL)
+            m_vVert[i].ptP -= plPlane.GetVersN() * dDist ;
+      }
+      else if ( dDist > 0)
+         m_vVert[i].nTemp = +1 ;
+      else
+         m_vVert[i].nTemp = -1 ;
+   }
+
+  // aggiorno timestamp dei triangoli
+   for ( auto& Tria : m_vTria)
+      Tria.nTemp = m_nTimeStamp ;
+   ++ m_nTimeStamp ;
+
+  // sistemo i triangoli (eventualmente li elimino)
+   for ( int i = 0 ; i < GetTriangleSize() ; ++ i) {
+     // salto i triangoli cancellati e quelli aggiunti
+      if ( m_vTria[i].nIdVert[0] == SVT_DEL  ||  m_vTria[i].nTemp == m_nTimeStamp)
+         continue ;
+      CutTriangleByPlane( i, plPlane, bSaveOnEq, bModif) ;
+   }
+
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+bool
+SurfTriMesh::CutTriangleByPlane( int nT, const Plane3d& plPlane, bool bSaveOnEq, bool& bModif)
+{
+   // se flag abilita e giace sul piano ed equiverso, lo salvo
+   if ( bSaveOnEq  &&
+        m_vVert[ m_vTria[nT].nIdVert[0]].nTemp == 0  &&
+        m_vVert[ m_vTria[nT].nIdVert[1]].nTemp == 0  &&
+        m_vVert[ m_vTria[nT].nIdVert[2]].nTemp == 0  &&
+        AreSameVectorApprox( m_vTria[nT].vtN, plPlane.GetVersN()))
+      ;
+   // se giace sul piano o dalla parte esterna del piano, lo cancello
+   else if ( m_vVert[ m_vTria[nT].nIdVert[0]].nTemp != -1  &&
+             m_vVert[ m_vTria[nT].nIdVert[1]].nTemp != -1  &&
+             m_vVert[ m_vTria[nT].nIdVert[2]].nTemp != -1) {
+      bModif = true ;
+      RemoveTriangle( nT) ;
+   }
+   // se giace dalla parte interna con al massimo un lato sul piano, lo conservo
+   else if ( m_vVert[ m_vTria[nT].nIdVert[0]].nTemp != 1  &&
+             m_vVert[ m_vTria[nT].nIdVert[1]].nTemp != 1  &&
+             m_vVert[ m_vTria[nT].nIdVert[2]].nTemp != 1)
+      ;
+   // altrimenti attraversa il piano, devo modificarlo
+   else {
+      bModif = true ;
+      int nDisc =  m_vVert[ m_vTria[nT].nIdVert[0]].nTemp +
+                   m_vVert[ m_vTria[nT].nIdVert[1]].nTemp +
+                   m_vVert[ m_vTria[nT].nIdVert[2]].nTemp ;
+      // se ha un vertice all'interno e uno sul piano -> 1 nuovo vertice e 1 nuovo triangolo
+      if ( nDisc == 0) {
+         // classifico i vertici
+         int nVertOn, nVertIn, nVertOut ;
+         if ( m_vVert[ m_vTria[nT].nIdVert[0]].nTemp == 0) {
+            nVertOn = 0 ;
+            nVertIn = 1 ;
+            nVertOut = 2 ;
+            if ( m_vVert[ m_vTria[nT].nIdVert[1]].nTemp > 0)
+               swap( nVertIn, nVertOut) ;
+         }
+         else if ( m_vVert[ m_vTria[nT].nIdVert[1]].nTemp == 0) {
+            nVertOn = 1 ;
+            nVertIn = 2 ;
+            nVertOut = 0 ;
+            if ( m_vVert[ m_vTria[nT].nIdVert[2]].nTemp > 0)
+               swap( nVertIn, nVertOut) ;
+         }
+         else {
+            nVertOn = 2 ;
+            nVertIn = 0 ;
+            nVertOut = 1 ;
+            if ( m_vVert[ m_vTria[nT].nIdVert[0]].nTemp > 0)
+               swap( nVertIn, nVertOut) ;
+         }
+         // calcolo il punto di intersezione del lato che attraversa il piano
+         double dDistIn = - DistPointPlane( m_vVert[m_vTria[nT].nIdVert[nVertIn]].ptP, plPlane) ;
+         double dDistOut = DistPointPlane( m_vVert[m_vTria[nT].nIdVert[nVertOut]].ptP, plPlane) ;
+         Point3d ptInt = Media( m_vVert[m_vTria[nT].nIdVert[nVertIn]].ptP,
+                                m_vVert[m_vTria[nT].nIdVert[nVertOut]].ptP,
+                                dDistIn / ( dDistIn + dDistOut)) ;
+         // inserisco il nuovo vertice
+         int nNewV = AddVertex( ptInt) ;
+         // inserisco il nuovo triangolo
+         int nIdVert[3] = { m_vTria[nT].nIdVert[0], m_vTria[nT].nIdVert[1], m_vTria[nT].nIdVert[2]} ;
+         nIdVert[nVertOut] = nNewV ;
+         int nIdTria = AddTriangle( nIdVert) ;
+         if ( nIdTria == SVT_NULL)
+            return false ;
+         if ( nIdTria != SVT_DEL)
+            m_vTria[nIdTria].nTemp = m_nTimeStamp ;
+         // cancello il vecchio triangolo
+         RemoveTriangle( nT) ;
+      }
+      // se ha un vertice all'interno e due all'esterno -> 2 nuovi vertici e 1 nuovo triangolo
+      else if ( nDisc == 1) {
+         // classifico i vertici
+         int nVertIn, nVertOut1, nVertOut2 ;
+         if ( m_vVert[ m_vTria[nT].nIdVert[0]].nTemp == -1) {
+            nVertIn = 0 ;
+            nVertOut1 = 1 ;
+            nVertOut2 = 2 ;
+         }
+         else if ( m_vVert[ m_vTria[nT].nIdVert[1]].nTemp == -1) {
+            nVertIn = 1 ;
+            nVertOut1 = 2 ;
+            nVertOut2 = 0 ;
+         }
+         else {
+            nVertIn = 2 ;
+            nVertOut1 = 0 ;
+            nVertOut2 = 1 ;
+         }
+         // calcolo i punti di intersezione dei lati che attraversano il piano
+         double dDistIn = - DistPointPlane( m_vVert[m_vTria[nT].nIdVert[nVertIn]].ptP, plPlane) ;
+         double dDistOut1 = DistPointPlane( m_vVert[m_vTria[nT].nIdVert[nVertOut1]].ptP, plPlane) ;
+         double dDistOut2 = DistPointPlane( m_vVert[m_vTria[nT].nIdVert[nVertOut2]].ptP, plPlane) ;
+         Point3d ptInt1 = Media( m_vVert[m_vTria[nT].nIdVert[nVertIn]].ptP,
+                                 m_vVert[m_vTria[nT].nIdVert[nVertOut1]].ptP,
+                                 dDistIn / ( dDistIn + dDistOut1)) ;
+         Point3d ptInt2 = Media( m_vVert[m_vTria[nT].nIdVert[nVertIn]].ptP,
+                                 m_vVert[m_vTria[nT].nIdVert[nVertOut2]].ptP,
+                                 dDistIn / ( dDistIn + dDistOut2)) ;
+         // inserisco i nuovi vertici
+         int nNewV1 = AddVertex( ptInt1) ;
+         int nNewV2 = AddVertex( ptInt2) ;
+         // inserisco il nuovo triangolo
+         int nIdVert[3] = { m_vTria[nT].nIdVert[0], m_vTria[nT].nIdVert[1], m_vTria[nT].nIdVert[2]} ;
+         nIdVert[nVertOut1] = nNewV1 ;
+         nIdVert[nVertOut2] = nNewV2 ;
+         int nIdTria = AddTriangle( nIdVert) ;
+         if ( nIdTria == SVT_NULL)
+            return false ;
+         if ( nIdTria != SVT_DEL)
+            m_vTria[nIdTria].nTemp = m_nTimeStamp ;
+         // cancello il vecchio triangolo
+         RemoveTriangle( nT) ;
+      }
+      // altrimenti ha due vertici all'interno e uno all'esterno -> 2 nuovi vertici e 2 nuovi triangoli
+      else {
+         // classifico i vertici
+         int nVertIn1, nVertIn2, nVertOut ;
+         if ( m_vVert[ m_vTria[nT].nIdVert[0]].nTemp == 1) {
+            nVertOut = 0 ;
+            nVertIn1 = 1 ;
+            nVertIn2 = 2 ;
+         }
+         else if ( m_vVert[ m_vTria[nT].nIdVert[1]].nTemp == 1) {
+            nVertOut = 1 ;
+            nVertIn1 = 2 ;
+            nVertIn2 = 0 ;
+         }
+         else {
+            nVertOut = 2 ;
+            nVertIn1 = 0 ;
+            nVertIn2 = 1 ;
+         }
+         // calcolo i punti di intersezione dei lati che attraversano il piano
+         double dDistIn1 = - DistPointPlane( m_vVert[m_vTria[nT].nIdVert[nVertIn1]].ptP, plPlane) ;
+         double dDistIn2 = - DistPointPlane( m_vVert[m_vTria[nT].nIdVert[nVertIn2]].ptP, plPlane) ;
+         double dDistOut = DistPointPlane( m_vVert[m_vTria[nT].nIdVert[nVertOut]].ptP, plPlane) ;
+         Point3d ptInt1 = Media( m_vVert[m_vTria[nT].nIdVert[nVertIn1]].ptP,
+                                 m_vVert[m_vTria[nT].nIdVert[nVertOut]].ptP,
+                                 dDistIn1 / ( dDistIn1 + dDistOut)) ;
+         Point3d ptInt2 = Media( m_vVert[m_vTria[nT].nIdVert[nVertIn2]].ptP,
+                                 m_vVert[m_vTria[nT].nIdVert[nVertOut]].ptP,
+                                 dDistIn2 / ( dDistIn2 + dDistOut)) ;
+         // inserisco i nuovi vertici
+         int nNewV1 = AddVertex( ptInt1) ;
+         int nNewV2 = AddVertex( ptInt2) ;
+         // inserisco i nuovi triangoli
+         int nIdVert1[3] = { m_vTria[nT].nIdVert[0], m_vTria[nT].nIdVert[1], m_vTria[nT].nIdVert[2]} ;
+         int nIdVert2[3] = { m_vTria[nT].nIdVert[0], m_vTria[nT].nIdVert[1], m_vTria[nT].nIdVert[2]} ;
+         double dSqDist1 = SqDist( m_vVert[m_vTria[nT].nIdVert[nVertIn1]].ptP, ptInt2) ;
+         double dSqDist2 = SqDist( m_vVert[m_vTria[nT].nIdVert[nVertIn2]].ptP, ptInt1) ;
+         if ( dSqDist1 <= dSqDist2) {
+            nIdVert1[nVertOut] = nNewV1 ;
+            nIdVert1[nVertIn2] = nNewV2 ;
+            nIdVert2[nVertOut] = nNewV2 ;
+         }
+         else {
+            nIdVert1[nVertOut] = nNewV1 ;
+            nIdVert2[nVertOut] = nNewV2 ;
+            nIdVert2[nVertIn1] = nNewV1 ;
+         }
+         int nIdTria1 = AddTriangle( nIdVert1) ;
+         int nIdTria2 = AddTriangle( nIdVert2) ;
+         if ( nIdTria1 == SVT_NULL  ||  nIdTria2 == SVT_NULL)
+            return false ;
+         if ( nIdTria1 != SVT_DEL)
+            m_vTria[nIdTria1].nTemp = m_nTimeStamp ;
+         if ( nIdTria2 != SVT_DEL)
+            m_vTria[nIdTria2].nTemp = m_nTimeStamp ;
+         // cancello il vecchio triangolo
+         RemoveTriangle( nT) ;
+      }
+   }
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+bool
+SurfTriMesh::CutByFacets( const Plane3d& plPlane, bool bSaveOnEq, bool& bModif)
+{
+  // Setto posizione triangoli non definita
+   for ( int nT = 0 ; nT < GetTriangleSize() ; ++ nT)
+      m_vTria[nT].nTempPart = 0 ;
+
+   INTERSCHAINMAP IntersLineMap ;
+
+  // Ciclo su tutte le facce.
+   for ( int nF = 0 ; nF < GetFacetCount() ; ++ nF) {
+     // Dati della faccia
+      POLYLINEVECTOR vLoopVec ;
+      GetFacetLoops( nF, vLoopVec) ;
+      if ( vLoopVec.empty())
+         continue ;
+     // Verifico la posizione del loop esterno
+      int nIntType = VerifyLoopPlane( vLoopVec[0], plPlane) ;
+     // Se interseca il piano di taglio, calcolo la divisione della faccia 
+      LineFacetClassVector IntersLinePart ;
+      if ( nIntType == FPI_CUT) {
+         PtrOwner<SurfFlatRegion> pReg( GetBasicSurfFlatRegion( GetSurfFlatRegionFromPolyLineVector( vLoopVec))) ;
+         if ( IsNull( pReg)  ||  ! pReg->IsValid())
+            return false ;
+         nIntType = IntersFacetPlane( *pReg, plPlane, IntersLinePart) ;
+      }
+     // Gestione dei risultati
+      if ( nIntType == FPI_CUT) {
+         for ( int nPart = 0 ; nPart < int( IntersLinePart.size()) ; ++ nPart) {
+           // Salvo intersezione per la faccia.
+            auto it = IntersLineMap.find( nF) ;
+            if ( it != IntersLineMap.end()) {
+               it->second.emplace_back( IntersInnSeg( IntersLinePart[nPart].ptSt, IntersLinePart[nPart].ptEn)) ;
+            }
+            else {
+               IntersLineMap.emplace( nF, IntersInnChain( 1, IntersInnSeg( IntersLinePart[nPart].ptSt, IntersLinePart[nPart].ptEn))) ;
+            }
+         }
+      }
+      else if ( nIntType == FPI_ON) { 
+         INTVECTOR vT ;
+         GetAllTriaInFacet( nF, vT) ;
+         Vector3d vtNf ;
+         GetFacetNormal( nF, vtNf) ;
+         for ( auto& nT : vT)
+            m_vTria[nT].nTempPart = ( vtNf * plPlane.GetVersN() > 0 ? 2 : -2) ;
+      }
+      else if ( nIntType == FPI_IN) {
+         INTVECTOR vT ;
+         GetAllTriaInFacet( nF, vT) ;
+         for ( auto& nT : vT)
+            m_vTria[nT].nTempPart = 1 ;
+      }
+      else if ( nIntType == FPI_OUT) {
+         INTVECTOR vT ;
+         GetAllTriaInFacet( nF, vT) ;
+         for ( auto& nT : vT)
+            m_vTria[nT].nTempPart = -1 ;
+      }
+      else
+         return false ;
+   }
+
+  // Divido le facce.
+   PieceMap NewFacet ;
+   SplitFacet( IntersLineMap, NewFacet) ;
+   INTERSEDGEMAP EdgeInterLineMap, EdgeEdgeLineMap ;
+   RetriangulateFacetPieces( NewFacet, EdgeInterLineMap, EdgeEdgeLineMap) ;
+
+  // Elimino i triangoli superflui
+   for ( int nT = 0 ; nT < GetTriangleSize() ; ++ nT) {
+      if (  m_vTria[nT].nTempPart == 0  ||
+            m_vTria[nT].nTempPart == -1  ||
+            m_vTria[nT].nTempPart == -2  ||
+            ( ! bSaveOnEq  &&  m_vTria[nT].nTempPart == 2)) {
+         RemoveTriangle( nT) ;
+         bModif = true ;
+      }
+   }
+
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+// Risultato : 0=nessuna intersezione, 1=intersezione è un punto, 2=intersezione è un segmento
+static int
+IntersRectangleTriangle( const Point3d& ptP, const Vector3d& vtL1, const Vector3d& vtL2,
+                         const Triangle3d& trTria, Point3d& ptStSeg, Point3d& ptEnSeg)
+{
+  // Assegno tolleranza lineare
+   double dTol = EPS_SMALL ;
+  // Definisco il piano del rettangolo
+   Plane3d plRectanglePlane ;
+   if ( ! plRectanglePlane.Set( ptP, vtL1 ^ vtL2))
+      return -1 ;
+  // Interseco il piano con il triangolo (recupero estremi in ordine inverso per mantenere l'esterno a destra nei loop)
+   int nPlTrIntRes = IntersPlaneTria( plRectanglePlane, trTria, ptEnSeg, ptStSeg) ;
+   if ( nPlTrIntRes != IntPlaneTriaType::IPTT_YES)
+      return 0 ;
+  // Limito il segmento col rettangolo
+   Vector3d vtSegDir = ptEnSeg - ptStSeg ;
+   double dSegLen = vtSegDir.Len() ;
+   vtSegDir /= dSegLen ;
+   // utilizzo piani ortogonali al rettangolo passanti per i suoi estremi
+   Plane3d plTrim1 ; plTrim1.Set( ptP, -vtL1) ;
+   Plane3d plTrim2 ; plTrim2.Set( ptP + vtL1, vtL1) ;
+   if ( vtSegDir * vtL1 < 0.)
+      swap( plTrim1, plTrim2) ;
+   Point3d ptInt1 ;
+   int nLnPl1IntRes = IntersLinePlane( ptStSeg, vtSegDir, dSegLen, plTrim1, ptInt1, false) ;
+   Point3d ptInt2 ;
+   int nLnPl2IntRes = IntersLinePlane( ptStSeg, vtSegDir, dSegLen, plTrim2, ptInt2, false) ;
+  // Se il segmento giace in uno dei due piani
+   if ( ( nLnPl1IntRes == ILPT_NO  &&  nLnPl2IntRes == ILPT_INPLANE)  ||
+        ( nLnPl2IntRes == ILPT_NO  &&  nLnPl1IntRes == ILPT_INPLANE))
+      return ( AreSamePointEpsilon( ptStSeg, ptEnSeg, dTol) ? 1 : 2) ;
+  // Se non ci sono intersezioni 
+   if ( nLnPl1IntRes == ILPT_NO  &&  nLnPl2IntRes == ILPT_NO) {
+     // se è tra i due piani
+      if (( ptStSeg - plTrim1.GetPoint()) * plTrim1.GetVersN() < 0.  &&
+          ( ptStSeg - plTrim2.GetPoint()) * plTrim2.GetVersN() < 0.)
+         return ( AreSamePointEpsilon( ptStSeg, ptEnSeg, dTol) ? 1 : 2) ;
+     // altrimenti è esterno
+      else
+         return 0 ;
+   }
+  // Posizioni parametriche delle intersezioni sul segmento
+   double dLen1 = ( ptInt1 - ptStSeg) * vtSegDir ;
+   double dLen2 = ( ptInt2 - ptStSeg) * vtSegDir ;
+  // Se la prima intersezione supera la fine o la seconda viene prima dell'inizio, non ci sono intersezioni
+   if ( dLen1 > dSegLen + EPS_ZERO  ||  dLen2 < -EPS_ZERO)
+      return 0 ;
+  // Eventuale aggiustamento inizio
+   if ( dLen1 > 0.) {
+      ptStSeg = ptInt1 ;
+      dSegLen -= dLen1 ;
+      dLen2 -= dLen1 ;
+   }
+  // Eventuale aggiustamento fine
+   if ( dLen2 < dSegLen) {
+      ptEnSeg = ptInt2 ;
+      dSegLen = dLen2 ;
+   }
+  // Assegno il tipo di intersezione
+   return ( dSegLen > dTol ? 2 : 1) ;
+}
+
+//----------------------------------------------------------------------------
+bool
+SurfTriMesh::GeneralizedCut( const ICurve& cvCurve, bool bSaveOnEq)
+{
+  // La superficie deve essere valida
+   if ( m_nStatus != OK)
+      return false ;
+  // La curva deve essere valida e chiusa, il vettore estrusione deve essere non nullo
+   Vector3d vtExtr ;
+   if ( ! cvCurve.GetExtrusion( vtExtr)  ||  vtExtr.IsSmall()  ||  ! cvCurve.IsClosed())
+      return false ;
+
+  // Recupero il numero originale di triangoli e di facce
+   int nTriaOriCnt = GetTriangleCount() ;
+   int nFacetOriCnt = GetFacetCount() ;
+
+  // Approssimo la curva con segmenti
+   CurveComposite cvCompo ;
+   {
+      PolyLine PL ;
+      if ( ! cvCurve.ApproxWithLines( LIN_TOL_FINE, ANG_TOL_STD_DEG, ICurve::APL_STD, PL)  ||  ! cvCompo.FromPolyLine( PL))
+         return false ;
+   }
+
+  // Eseguo scalature
+   Frame3d frScalingRef;
+   frScalingRef.Set( m_vVert[0].ptP, X_AX, Y_AX, Z_AX) ;
+   Scale( frScalingRef, CUT_SCALE, CUT_SCALE, CUT_SCALE) ;
+   cvCompo.Scale( frScalingRef, CUT_SCALE, CUT_SCALE, CUT_SCALE) ;
+
+  // Appiattisco la polilinea nel piano perpendicolare all'estrusione
+   Frame3d frCurve ;
+   Point3d ptStart ; cvCompo.GetStartPoint( ptStart) ;
+   frCurve.Set( ptStart, vtExtr) ;
+   cvCompo.ToLoc( frCurve) ;
+   if ( ! cvCompo.Scale( GLOB_FRM, 1, 1, 0)) {
+      Scale( frScalingRef, 1. / CUT_SCALE, 1. / CUT_SCALE, 1. / CUT_SCALE) ;
+      return false ;
+   }
+   double dArea ;
+   cvCompo.GetAreaXY( dArea) ;
+   BBox3d b3Crv ;
+   cvCompo.GetLocalBBox( b3Crv) ;
+   cvCompo.ToGlob( frCurve) ;
+  // Assegno il senso di rotazione della curva (visto dalla punta del vettore estrusione)
+   bool bCCW = ( dArea > 0) ;
+
+  // Recupero Bounding-box della trimesh
+   BBox3d b3SurfBox ;
+   GetLocalBBox( b3SurfBox) ;
+  // Trovo minima e massima distanza dei vertici del bounding-box della TriMesh dal piano della curva
+   b3SurfBox.ToLoc( frCurve) ;
+   Point3d ptMin, ptMax ;
+   b3SurfBox.GetMinMax( ptMin, ptMax) ;
+   Vector3d vtMax = ( ptMax.z + 10) * vtExtr ;
+   Vector3d vtMin = ( ptMin.z - 10) * vtExtr ; 
+
+  // Ciclo sui triangoli
+   bool bModif = false ;
+   int nNumTria = GetTriangleSize() ;
+   for ( int nT = 0 ; nT < nNumTria ; ++ nT) {
+     // Recupero il triangolo
+      Triangle3d trTria ;
+      if ( ! GetTriangle( nT, trTria))
+         continue ;
+     // Box del triangolo nel riferimento locale della curva
+      BBox3d b3Tria ;
+      trTria.GetLocalBBox( b3Tria) ;
+      b3Tria.ToLoc( frCurve) ;
+     // Se il box del triangolo non interseca quello locale della curva
+      if ( ! b3Crv.OverlapsXY( b3Tria)) {
+        // Se la parte da conservare è quella all'interno della curva, elimino il triangolo 
+         if ( bCCW) {
+            RemoveTriangle( nT) ;
+            bModif = true ;
+         }
+         continue ;
+      }
+     // Determino se il centro del triangolo cade all'interno della curva
+      bool bTriaCenIn = false ;
+      {
+        // centro proiettato sul piano della curva
+         Point3d ptCen = trTria.GetCentroid() ;
+         ptCen = ptCen - ( ptCen - ptStart) * vtExtr * vtExtr ;
+        // calcolo distanza
+         DistPointCurve dstPC( ptCen, cvCompo) ;
+         double dDist ; dstPC.GetDist( dDist) ;
+        // se maggiore oltre il limite originale
+         if ( dDist > CUT_SCALE * EPS_SMALL) {
+            int nSide ;
+            if ( dstPC.GetSideAtMinDistPoint( 0, vtExtr, nSide))
+               bTriaCenIn = ( nSide == MDS_LEFT) ;
+         }
+        // altrimenti ricalcolo con centro del triangolo spostato all'interno sulla normale
+         else {
+            Point3d ptCen2 = trTria.GetCentroid() - CUT_SCALE * EPS_SMALL * trTria.GetN() ;
+            DistPointCurve dstP2C( ptCen2, cvCompo) ;
+            int nSide2 ;
+            if ( dstP2C.GetSideAtMinDistPoint( 0, vtExtr, nSide2))
+               bTriaCenIn = ( nSide2 == MDS_LEFT  &&  bSaveOnEq) ;
+         }
+      }
+
+     // Vettore di catene di punti e flag sulla complanarità del triangolo con uno dei rettangoli
+      CHAINVECTOR vChain ;
+      bool bTriaOn = false ;
+      // Segnalatore di taglio piccolo, punto e vettore per costruire il piano per eseguirlo
+      int nSmallCutState = 0 ;
+      Point3d ptSmallCutPnt ;
+      Vector3d vtSmallCutVec ;
+     // Ciclo sui segmenti
+      int nChainCnt = 0 ;
+      bool bChain = false ;
+      Point3d ptChSt, ptChEn ;
+      const ICurve* pCrv = cvCompo.GetFirstCurve() ;
+      while ( pCrv != nullptr) {
+        // estremi del segmento
+         Point3d ptSt ;  pCrv->GetStartPoint( ptSt) ;
+         Point3d ptEn ;  pCrv->GetEndPoint( ptEn) ;
+        // Se il segmento giace sul piano del triangolo, lo segnalo.
+         if ( abs( ( ptSt - trTria.GetCentroid()) * trTria.GetN()) < CUT_SCALE * EPS_SMALL  &&
+              abs( ( ptEn - trTria.GetCentroid()) * trTria.GetN()) < CUT_SCALE * EPS_SMALL  &&
+              abs( trTria.GetN() * vtExtr) < EPS_SMALL)
+            bTriaOn = true ;
+        // Intersezione fra il rettangolo (ottenuto dall'estrusione del segmento corrente) e il triangolo
+         Point3d ptSegSt, ptSegEn ;
+         int nInt = IntersRectangleTriangle( ptSt + vtMin, ptEn - ptSt, vtMax - vtMin, trTria, ptSegSt, ptSegEn) ;
+         if ( nInt == 2) {
+           // Creo nuova catena se non c'è già o se discontinuità
+            if ( ! bChain  ||  ( ! AreSamePointApprox( ptSegSt, ptChEn)  &&  ! AreSamePointApprox( ptSegEn, ptChSt))) {
+               ++ nChainCnt ;
+               vChain.resize( nChainCnt) ;
+               bChain = false ;
+            }
+           // Assegno i dati di intersezione
+            IntSegment CurInters ;
+            CurInters.ptSt = ptSegSt ;
+            CurInters.ptEn = ptSegEn ;
+           // Inserisco nella catena
+            if ( ! bChain) {
+               vChain[nChainCnt - 1].emplace_back( CurInters) ;
+               ptChSt = CurInters.ptSt ;
+               ptChEn = CurInters.ptEn ;
+            }
+            else if ( AreSamePointApprox(ptSegSt, ptChEn)) {
+               vChain[nChainCnt - 1].emplace_back(CurInters) ;
+               ptChEn = CurInters.ptEn ;
+            }
+            else {
+               vChain[nChainCnt - 1].insert(vChain[nChainCnt - 1].begin(), CurInters) ;
+               ptChSt = CurInters.ptSt ;
+            }
+            bChain = true ;
+         }
+         else {
+            // Taglio piccolo
+            if ( nInt == 1) {
+               if ( nSmallCutState == 0)
+                  nSmallCutState = 1 ;
+               else
+                  nSmallCutState = 2 ;
+               ptSmallCutPnt = ptSt ;
+               vtSmallCutVec = ptEn - ptSt ;
+            }
+            bChain = false ;
+         }
+         pCrv = cvCompo.GetNextCurve() ;
+      }
+      // Gestisco il taglio piccolo.
+      if ( nSmallCutState == 1  &&  int( vChain.size()) == 0) {
+         Plane3d plPlane ;
+         plPlane.Set( ptSmallCutPnt, vtSmallCutVec ^ ( vtMax - vtMin)) ;
+         for ( int nV = 0 ; nV < 3 ; ++ nV) {
+            double dDist = DistPointPlane( m_vVert[m_vTria[nT].nIdVert[nV]].ptP, plPlane) ;
+            if ( abs( dDist) < EPS_SMALL) {
+               m_vVert[m_vTria[nT].nIdVert[nV]].nTemp = 0 ;
+               if ( abs( dDist) > EPS_SMALL)
+                  m_vVert[m_vTria[nT].nIdVert[nV]].ptP -= plPlane.GetVersN() * dDist ;
+            }
+            else if ( dDist > 0)
+               m_vVert[m_vTria[nT].nIdVert[nV]].nTemp = +1 ;
+            else
+               m_vVert[m_vTria[nT].nIdVert[nV]].nTemp = -1 ;
+         }
+         CutTriangleByPlane( nT, plPlane, bSaveOnEq, bModif) ;
+         continue ;
+      }
+      for ( auto itI = vChain.begin() ; itI != vChain.end() ; ) {
+         bool bErased = false ;
+         auto itJ = itI ;
+         ++ itJ ;
+         for ( ; itJ != vChain.end()  &&  ! bErased ; ) {
+            if ( int( itI->size()) == 1  &&  int( itJ->size()) == 1  &&
+                 AreSamePointEpsilon( itI->back().ptSt, itJ->back().ptEn, 2 * EPS_SMALL)  &&
+                 AreSamePointEpsilon( itI->back().ptEn, itJ->back().ptSt, 2 * EPS_SMALL)) {
+               itJ = vChain.erase( itJ) ;
+               bErased = true ;
+            }
+            else
+               ++ itJ ;
+         }
+         if ( bErased) {
+            itI = vChain.erase( itI) ;
+            bErased = false ;
+         }
+         else
+            ++ itI ;
+      }
+      nChainCnt = int( vChain.size()) ;
+
+     // unisco eventuali catene estreme che sono parte di una stessa catena
+      if ( nChainCnt > 1) {
+         if ( AreSamePointApprox( vChain[0].front().ptSt, vChain[nChainCnt - 1].back().ptEn)) {
+            vChain[0].insert( vChain[0].begin(), vChain[nChainCnt - 1].begin(), vChain[nChainCnt - 1].end()) ;
+            vChain.pop_back() ;
+            -- nChainCnt ;
+         }
+         else if ( AreSamePointApprox( vChain[0].back().ptEn, vChain[nChainCnt - 1].front().ptSt)) {
+            vChain[0].insert(vChain[0].end(), vChain[nChainCnt - 1].begin(), vChain[nChainCnt - 1].end()) ;
+            vChain.pop_back() ;
+            -- nChainCnt ;
+         }
+      }
+
+     // Elimino la seconda copia di catene doppie
+      for ( int nI = 0 ; nI < nChainCnt ; ++ nI) {
+         for ( int nJ = nI + 1 ; nJ < nChainCnt ; ++ nJ) {
+            if ( vChain[nI].size() == vChain[nJ].size()) {
+               bool bSame = true ;
+               for ( int nK = 0 ; nK < int( vChain[nI].size()) ; ++ nK) {
+                  if ( ! AreSamePointApprox( vChain[nI][nK].ptSt, vChain[nJ][nK].ptSt)  ||
+                       ! AreSamePointApprox( vChain[nI][nK].ptEn, vChain[nJ][nK].ptEn)) {
+                     bSame = false ;
+                     break ;
+                  }
+               }
+               if ( bSame) {
+                  vChain.erase( vChain.begin() + nJ) ;
+                  -- nChainCnt ;
+                  -- nJ ;
+               }
+            }
+         }
+      }
+
+     // Fra le catene trovate separo le aperte dalle chiuse
+      CHAINVECTOR cvClosedChain ;
+      CHAINVECTOR cvOpenChain ;
+      for ( int nL = 0 ; nL < int( vChain.size()) ; ++ nL) {
+         int nCurLoopLast = max( int(vChain[nL].size()) - 1, 0) ;
+         if ( AreSamePointApprox( vChain[nL][0].ptSt, vChain[nL][nCurLoopLast].ptEn)  &&  nCurLoopLast > 0)
+            cvClosedChain.emplace_back( vChain[nL]) ;
+         else {
+            cvOpenChain.emplace_back( vChain[nL]) ;
+         }
+      }
+      for ( auto it = cvClosedChain.begin() ; it != cvClosedChain.end() ; ) {
+         if ( int( it->size()) < 3)
+            it = cvClosedChain.erase( it) ;
+         else
+            ++ it ;
+      }
+     // Se più di una catena chiusa oppure catene chiuse e aperte, errore
+      if ( cvClosedChain.size() > 1  ||  ( cvClosedChain.size() > 0  &&  int( cvOpenChain.size()) > 0)) {
+         Scale( frScalingRef, 1. / CUT_SCALE, 1. / CUT_SCALE, 1. / CUT_SCALE) ;
+         return false ;
+      }
+
+     // Se c'è una catena chiusa
+      if ( cvClosedChain.size() == 1) {
+        // Ne ricavo una PolyLine
+         PolyLine plInLoop ;
+         for ( int nLine = 0 ; nLine < int( cvClosedChain[0].size()) ; ++ nLine) {
+            plInLoop.AddUPoint( 0., cvClosedChain[0][nLine].ptSt) ;
+            plInLoop.AddUPoint( 0., cvClosedChain[0][nLine].ptEn) ;
+         }
+        // I tre vertici sono dalla parte interna della curva (triangolo con buco)
+         if ( ! bCCW) {
+           // Rimuovo il triangolo corrente
+            RemoveTriangle( nT) ;
+           // Definisco il loop esterno (è il triangolo)
+            PolyLine plExtLoop ;
+            plExtLoop.AddUPoint( 0., trTria.GetP( 0)) ;
+            plExtLoop.AddUPoint( 0., trTria.GetP( 1)) ;
+            plExtLoop.AddUPoint( 0., trTria.GetP( 2)) ;
+            plExtLoop.AddUPoint( 0., trTria.GetP( 0)) ;
+           // Eseguo triangolazione
+            POLYLINEVECTOR vPL ;
+            vPL.emplace_back( plExtLoop) ;
+            vPL.emplace_back( plInLoop) ;
+            PNTVECTOR vPt ;
+            INTVECTOR vTr ;
+            if ( Triangulate().Make( vPL, vPt, vTr)) {
+              // Inserisco i nuovi triangoli
+               for ( int n = 0 ; n < int( vTr.size()) - 2 ; n += 3) {
+                  int nNewTriaVertId[3] = { vTr[n], vTr[n + 1], vTr[n + 2] } ;
+                  int nNewId[3] = { AddVertex( vPt[nNewTriaVertId[0]]),
+                                    AddVertex( vPt[nNewTriaVertId[1]]),
+                                    AddVertex( vPt[nNewTriaVertId[2]]) } ;
+                  AddTriangle( nNewId) ;
+                  bModif = true ;
+               }
+            }
+         }
+        // Se nessun vertice dalla parte interna della curva (rimane solo l'area della curva)
+         else {
+           // Rimuovo il triangolo corrente
+            RemoveTriangle( nT) ;
+           // Eseguo triangolazione
+            PNTVECTOR vPt ;
+            INTVECTOR vTr ;
+            if ( Triangulate().Make( plInLoop, vPt, vTr)) {
+              // Inserisco i nuovi triangoli
+               for ( int n = 0 ; n < int( vTr.size()) - 2 ; n += 3) {
+                  int nNewTriaVertId[3] = { vTr[n], vTr[n + 1], vTr[n + 2] } ;
+                  int nNewId[3] = { AddVertex(vPt[nNewTriaVertId[0]]),
+                                    AddVertex(vPt[nNewTriaVertId[1]]),
+                                    AddVertex(vPt[nNewTriaVertId[2]]) } ;
+                  AddTriangle( nNewId) ;
+                  bModif = true ;
+               }
+            }
+         }
+      }
+
+     // Loop aperti, devo chiuderli
+      else if ( cvOpenChain.size() > 0) {
+        // Creo il loop chiuso padre di tutti, il perimetro del triangolo.
+        // Questo viene diviso in sotto-loop chiusi mediante quelli aperti.
+        // I loop chiusi trovati precedentemente sono interni a uno dei sotto-loop 
+        // chiusi di cui è formato il perimetro.
+         PNTVECTOR cvFirstLoop ;
+         cvFirstLoop.emplace_back( trTria.GetP( 0)) ;
+         cvFirstLoop.emplace_back( trTria.GetP( 1)) ;
+         cvFirstLoop.emplace_back( trTria.GetP( 2)) ;
+
+         PNTMATRIX cvBoundClosedLoopVec ;
+         cvBoundClosedLoopVec.emplace_back( cvFirstLoop) ;
+         BOOLVECTOR vbInOut ;
+         vbInOut.push_back( true) ;
+        // Divido il loop di partenza in sotto-loop
+         while ( cvOpenChain.size() > 0) {
+            int nLastOpenLoopN = int( cvOpenChain.size()) - 1 ;
+            for ( int nLoop = 0 ; nLoop < int( cvBoundClosedLoopVec.size()) ; ++ nLoop) {
+              // Estremi del loop aperto
+               int nLastOpenLoopPoint = max( int( cvOpenChain[nLastOpenLoopN].size()) - 1, 0) ;
+               Point3d ptOpenLoopStP = cvOpenChain[nLastOpenLoopN][0].ptSt ;
+               Point3d ptOpenLoopEnP = cvOpenChain[nLastOpenLoopN][nLastOpenLoopPoint].ptEn ;
+               PNTVECTOR Loop1, Loop2 ;
+               bool bChangedStart = ChangeStart( ptOpenLoopStP, cvBoundClosedLoopVec[nLoop]) ;
+               bool bSplitted = SplitAtPoint( ptOpenLoopEnP, cvBoundClosedLoopVec[nLoop], Loop1, Loop2) ;
+               if ( ! ( bChangedStart  &&  bSplitted))
+                  continue ;
+               Chain cvCounterChain ;
+               for ( int nPt = int( cvOpenChain[nLastOpenLoopN].size()) - 1 ; nPt >= 0 ; -- nPt) {
+                  IntSegment CurSeg ;
+                  CurSeg.ptSt = cvOpenChain[nLastOpenLoopN][nPt].ptEn ;
+                  CurSeg.ptEn = cvOpenChain[nLastOpenLoopN][nPt].ptSt ;
+                  cvCounterChain.emplace_back( CurSeg) ;
+               }
+               bool bAdded1 = AddChainToChain( cvCounterChain, Loop1) ;
+               bool bAdded2 = AddChainToChain( cvOpenChain[nLastOpenLoopN], Loop2) ;
+               if ( ! ( bAdded1  &&  bAdded2))
+                  continue ;
+              // Aggiungo i nuovi loop nel vettore
+               int nCurSize = int( cvBoundClosedLoopVec.size()) ;
+               cvBoundClosedLoopVec.resize( nCurSize + 1) ;
+               vbInOut.resize( nCurSize + 1) ;
+               for ( int nCL = nCurSize - 1 ; nCL > nLoop ; -- nCL) {
+                  cvBoundClosedLoopVec[nCL + 1] = cvBoundClosedLoopVec[nCL] ;
+                  vbInOut[nCL + 1] = vbInOut[nCL] ;
+               }
+               cvBoundClosedLoopVec[nLoop] = Loop1 ;
+               cvBoundClosedLoopVec[nLoop + 1] = Loop2 ;
+               vbInOut[nLoop] = false ;
+               vbInOut[nLoop + 1] = true ;
+               ++ nLoop ;
+            }
+            cvOpenChain.resize( nLastOpenLoopN) ;
+         }
+        // Rimuovo il triangolo corrente
+         RemoveTriangle( nT) ;
+        // Se il triangolo originale aveva il centro interno alla curva oppure devo tenere le parti 
+        // di superficie che hanno una parte giacente sul rettangolo della superficie di taglio, 
+        // aggiungo i nuovi triangoli.
+         if ( ! bTriaOn  ||  bSaveOnEq) {
+            // Trasformo i loop compositi in loop polyline
+            POLYLINEVECTOR vplPolyVec ;
+            vplPolyVec.resize( cvBoundClosedLoopVec.size()) ;
+            for ( int nLoop = 0 ; nLoop < int( vplPolyVec.size()) ; ++ nLoop) {
+               for ( int nLine = 0 ; nLine < int( cvBoundClosedLoopVec[nLoop].size()) ; ++ nLine) {
+                  vplPolyVec[nLoop].AddUPoint( 0., cvBoundClosedLoopVec[nLoop][nLine]) ;
+               }
+               vplPolyVec[nLoop].AddUPoint( 0., cvBoundClosedLoopVec[nLoop][0]) ;
+               if ( vbInOut[nLoop]) {
+                  // Eseguo triangolazione
+                  Triangulate CreateTriangulation ;
+                  PNTVECTOR vPt ;
+                  INTVECTOR vTr ;
+                  if ( Triangulate().Make( vplPolyVec[nLoop], vPt, vTr)) {
+                     // Inserisco i nuovi triangoli
+                     for ( int n = 0 ; n < int( vTr.size()) - 2 ; n += 3) {
+                        int nNewTriaVertId[3] = { vTr[n], vTr[n + 1], vTr[n + 2] } ;
+                        int nNewId[3] = { AddVertex(vPt[nNewTriaVertId[0]]),
+                                          AddVertex(vPt[nNewTriaVertId[1]]),
+                                          AddVertex(vPt[nNewTriaVertId[2]]) } ;
+                        AddTriangle( nNewId) ;
+                        bModif = true ;
+                     }
+                  }
+               }
+            }
+         }
+      }
+
+     // da eliminare
+      else if ( ! bTriaCenIn) {
+         RemoveTriangle( nT) ;
+         bModif = true ;
+      }
+   }
+
+  // Se avvenuta modifica, aggiorno tutto
+   if ( bModif) {
+     // aggiorno tutto
+      if ( ! AdjustVertices()  ||  ! DoCompacting()) {
+         Scale( frScalingRef, 1. / CUT_SCALE, 1. / CUT_SCALE, 1. / CUT_SCALE) ;
+         return false ;
+      }
+   }
+
+  // Ripristino scala originale
+   Scale( frScalingRef, 1. / CUT_SCALE, 1. / CUT_SCALE, 1. / CUT_SCALE) ;
+
+  // se superficie originale a facce, cerco di semplificarle in ogni caso
+   if ( nFacetOriCnt < 200  ||  double( nTriaOriCnt) / nFacetOriCnt > 4) {
+      if ( ! SimplifyFacets())
+         LOG_ERROR( GetEGkLogger(), "Error in SimplifyFacets of Stm::GeneralizedCut")
+   }
+
+   return true ;
+}

SurfTriMeshFaceting.cpp

@@ -22,15 +22,6 @@
 
 using namespace std ;
 
-//----------------------------------------------------------------------------
-bool
-SurfTriMesh::ResetFaceting( void)
-{
-   m_bFaceted = false ;
-   m_vFacet.clear() ;
-   return true ;
-}
-
 //----------------------------------------------------------------------------
 bool
 SurfTriMesh::VerifyFaceting( void) const
@@ -44,14 +35,36 @@ SurfTriMesh::VerifyFaceting( void) const
 bool
 SurfTriMesh::UpdateFaceting( void)
 {
+  // salvo vecchio vettore facce
+   INTVECTOR vOldFacet = m_vFacet ;
+
   // reset faceting
    m_bFaceted = false ;
+   m_vFacet.clear() ;
    for ( int i = 0 ; i < int( m_vTria.size()) ; ++ i)
       m_vTria[i].nIdFacet = SVT_NULL ;
 
-  // ricostruisco le sfaccettature
-   bool bOk = true ;
+  // indice faccia corrente
    int nFacet = -1 ;
+
+  // ricostruisco le sfaccettature come definite in precedenza (dove possibile)
+   bool bOk = true ;
+   for ( int j = 0 ; j < int( vOldFacet.size()) ; ++ j) {
+      int i = vOldFacet[j] ;
+     // salto triangoli inesistenti o già assegnati
+      if ( i >= int( m_vTria.size())  ||
+           m_vTria[i].nIdVert[0] == SVT_DEL  ||
+           m_vTria[i].nIdFacet != SVT_NULL)
+         continue ;  
+     // assegno indice di faccia al triangolo
+      m_vTria[i].nIdFacet = ++ nFacet ;
+      m_vFacet.push_back( i) ;
+     // aggiorno faccia
+      if ( ! UpdateOneFace( nFacet, i))
+         bOk = false ;
+   }
+
+  // ricostruisco le altre sfaccettature
    for ( int i = 0 ; i < int( m_vTria.size()) ; ++ i) {
      // salto triangoli cancellati o già assegnati
       if ( m_vTria[i].nIdVert[0] == SVT_DEL  ||
@@ -60,32 +73,9 @@ SurfTriMesh::UpdateFaceting( void)
      // assegno indice di faccia al triangolo
       m_vTria[i].nIdFacet = ++ nFacet ;
       m_vFacet.push_back( i) ;
-     // piano del triangolo
-      Plane3d plPlane ;
-      if ( ! plPlane.Set( m_vVert[m_vTria[i].nIdVert[0]].ptP, m_vTria[i].vtN)) {
-         LOG_ERROR( GetEGkLogger(), "SurfTM : UpdateFaceting error in triangle data")        
-         return false ;
-      }
-     // set di triangoli da aggiornare
-      set<int> stTria ;
-      stTria.insert( i) ;
-     // finchè set non vuoto
-      while ( ! stTria.empty()) {
-        // tolgo un triangolo dal set
-         const auto iIt = stTria.begin() ;
-         int nT = *iIt ;
-         stTria.erase( iIt) ;
-        // aggiorno i triangoli adiacenti
-         for ( int j = 0 ; j < 3 ; ++ j) {
-            int nAdjT = m_vTria[nT].nIdAdjac[j] ;
-            if ( nAdjT != SVT_NULL  &&  m_vTria[nAdjT].nIdFacet == SVT_NULL) {
-               if ( ! UpdateTriaFaceting( nT, nFacet, plPlane, nAdjT))
-                  bOk = false ;
-               if ( m_vTria[nAdjT].nIdFacet == nFacet)
-                  stTria.insert( nAdjT) ;
-            }
-         }
-      }
+     // aggiorno faccia
+      if ( ! UpdateOneFace( nFacet, i))
+         bOk = false ;
    }
 
   // se ci sono stati problemi, salvo nel log
@@ -97,6 +87,40 @@ SurfTriMesh::UpdateFaceting( void)
    return true ;
 }
 
+//----------------------------------------------------------------------------
+bool
+SurfTriMesh::UpdateOneFace( int nFacet, int nT)
+{
+  // piano del triangolo
+   Plane3d plPlane ;
+   if ( ! plPlane.Set( m_vVert[m_vTria[nT].nIdVert[0]].ptP, m_vTria[nT].vtN)) {
+      LOG_ERROR( GetEGkLogger(), "SurfTM : UpdateFaceting error in triangle data")        
+      return false ;
+   }
+  // set di triangoli da aggiornare
+   set<int> stTria ;
+   stTria.insert( nT) ;
+  // finchè set non vuoto
+   bool bOk = true ;
+   while ( ! stTria.empty()) {
+     // tolgo un triangolo dal set
+      const auto iIt = stTria.begin() ;
+      int nT = *iIt ;
+      stTria.erase( iIt) ;
+     // aggiorno i triangoli adiacenti
+      for ( int j = 0 ; j < 3 ; ++ j) {
+         int nAdjT = m_vTria[nT].nIdAdjac[j] ;
+         if ( nAdjT != SVT_NULL  &&  m_vTria[nAdjT].nIdFacet == SVT_NULL) {
+            if ( ! UpdateTriaFaceting( nT, nFacet, plPlane, nAdjT))
+               bOk = false ;
+            if ( m_vTria[nAdjT].nIdFacet == nFacet)
+               stTria.insert( nAdjT) ;
+         }
+      }
+   }
+   return bOk ;
+}
+
 //----------------------------------------------------------------------------
 bool
 SurfTriMesh::UpdateTriaFaceting( int nRefT, int nFacet, const Plane3d& plPlane, int nT)

SurfTriMeshUtilities.cpp

@@ -19,6 +19,141 @@
 
 using namespace std ;
 
+//----------------------------------------------------------------------------
+bool
+SurfTriMesh::RemoveTJunctions(void)
+{
+  // Vettore di indici dei vertici sui lati del triangolo corrente
+   unordered_map< int, INTVECTOR> TriaMap ;
+
+  // Ciclo sui triangoli della superficie per determinare gli altri vertici sul loro perimetro
+   for ( int nT = 0 ; nT < int( m_vTria.size()) ; ++ nT) {
+     // Se il triangolo non è valido, passo al successivo
+      Triangle3d trTria ;
+      if ( ! GetTriangle( nT, trTria)  ||  ! trTria.Validate( true))
+         continue ;
+     // Vettore degli altri vertici sul contorno del triangolo 
+      INTVECTOR vVertOtl ;
+     // Box del triangolo 
+      BBox3d b3Tria ;
+      trTria.GetLocalBBox( b3Tria) ;
+      INTVECTOR vNearTria ;
+      GetAllTriaOverlapBox( b3Tria, vNearTria) ;
+     // Ciclo sui lati del triangolo
+      for ( int nSeg = 0 ; nSeg < 3 ; ++ nSeg) {
+        // aggiungo al vettore il vertice iniziale del lato
+         vVertOtl.emplace_back( m_vTria[nT].nIdVert[nSeg]) ;
+        // Se in questo lato il triangolo è adiacente a un altro, lo salto.
+         if ( m_vTria[nT].nIdAdjac[nSeg] != SVT_DEL  &&  m_vTria[nT].nIdAdjac[nSeg] != SVT_NULL)
+            continue ;
+         int nPrevSize = int( vVertOtl.size()) ;
+        // recupero la geometria del lato
+         Point3d ptSegSt = trTria.GetP( nSeg) ;
+         Point3d ptSegEn = trTria.GetP( ( nSeg + 1) % 3) ;
+         Vector3d vtSeg = ptSegEn - ptSegSt ;
+         double dSegLen = vtSeg.Len() ;
+         if ( dSegLen < EPS_SMALL)
+            continue ;
+         vtSeg /= dSegLen ;
+        // Ciclo sui triangoli vicini
+         for ( int nI = 0 ; nI < int( vNearTria.size()) ; ++ nI) {
+           // Salto il triangolo se è quello di riferimento
+            if ( vNearTria[nI] == nT)
+               continue ;
+           // Cerco i vertici che stanno sul lato del triangolo
+            for ( int nVert = 0 ; nVert < 3 ; ++ nVert) {
+               Point3d ptVert ;
+               if ( ! GetVertex( m_vTria[vNearTria[nI]].nIdVert[nVert], ptVert))
+                  continue ;
+               double dProj = ( ptVert - ptSegSt) * vtSeg ;
+               double dOrt = ( ( ptVert - ptSegSt) - dProj * vtSeg).SqLen() ;               
+               if ( dProj > EPS_SMALL  &&  dProj < dSegLen - EPS_SMALL  &&  dOrt < SQ_EPS_TRIA_H)
+                  vVertOtl.emplace_back( m_vTria[vNearTria[nI]].nIdVert[nVert]) ;                
+            }
+         }
+        // Riordino i vertici sul segmento
+         auto SortVerteces = [ this, &ptSegSt, &vtSeg]( const int nV1, const int nV2)
+                                 {  Point3d ptV1, ptV2 ;
+                                    GetVertex( nV1, ptV1) ;
+                                    GetVertex( nV2, ptV2) ;
+                                    return ( ( ptV1 - ptSegSt) * vtSeg < ( ptV2 - ptSegSt) * vtSeg) ;
+                                 } ;
+         sort( vVertOtl.begin() + nPrevSize, vVertOtl.end(), SortVerteces) ;
+      }
+     // Se ci sono più di 3 vertici
+      if ( vVertOtl.size() > 3) {
+        // Elimino i vertici ripetuti
+         vVertOtl.erase( unique( vVertOtl.begin(), vVertOtl.end()), vVertOtl.end()) ;
+        // Se ci sono ancora più di 3 vertici, inserisco nel Map
+         if ( vVertOtl.size() > 3)
+            TriaMap.emplace( nT, vVertOtl) ;
+      }
+   }
+  
+  // Ciclo sui triangoli da sistemare
+   for ( auto itT = TriaMap.begin() ; itT != TriaMap.end() ; ++ itT) {
+     // Indice del triangolo
+      int nT = itT->first ;
+     // Vettore degli altri vertici sul perimetro
+      const INTVECTOR& vVertOtl = itT->second ;
+     // Se il triangolo non è valido, passo al successivo
+      Triangle3d trTria ;
+      if ( ! GetTriangle( nT, trTria)  ||  ! trTria.Validate( true))
+         continue ;
+     // Salvo eventuale indice di faccetta
+      int nFacet = m_vTria[nT].nIdFacet ;
+      if ( nFacet > int( m_vFacet.size())  ||  m_vFacet[nFacet] != nT)
+         nFacet = SVT_NULL ;
+     // Rimuovo il triangolo
+      int nTFlag = m_vTria[nT].nTFlag ;
+      RemoveTriangle( nT) ;
+     // Aggiungo i nuovi triangoli
+      int nLastNewTria = SVT_NULL ;
+     // Se ci sono 4 vertici, inserisco due triangoli
+      if ( vVertOtl.size() == 4) {
+        // se 1-2-3 è triangolo (e quindi 0-1-3)
+         int nNew1Id[3] = { vVertOtl[1], vVertOtl[2], vVertOtl[3]} ;
+         int nNew1Tria = AddTriangle( nNew1Id, nTFlag) ;
+         if ( nNew1Tria != SVT_NULL  &&  nNew1Tria != SVT_DEL) {
+            nLastNewTria = nNew1Tria ;
+            int nNew2Id[3] = { vVertOtl[0], vVertOtl[1], vVertOtl[3]} ;
+            int nNew2Tria = AddTriangle( nNew2Id, nTFlag) ;
+            if ( nNew2Tria != SVT_NULL  &&  nNew2Tria != SVT_DEL)
+               nLastNewTria = nNew2Tria ;
+         }
+        // altrimenti 0-1-2 e 2-3-0
+         else {
+            int nNew3Id[3] = { vVertOtl[0], vVertOtl[1], vVertOtl[2]} ;
+            int nNew3Tria = AddTriangle( nNew3Id, nTFlag) ;
+            if ( nNew3Tria != SVT_NULL  &&  nNew3Tria != SVT_DEL)
+               nLastNewTria = nNew3Tria ;
+            int nNew4Id[3] = { vVertOtl[2], vVertOtl[3], vVertOtl[0]} ;
+            int nNew4Tria = AddTriangle( nNew4Id, nTFlag) ;
+            if ( nNew4Tria != SVT_NULL  &&  nNew4Tria != SVT_DEL)
+               nLastNewTria = nNew4Tria ;
+         }
+      }
+     // altrimenti inserisco un ventaglio di triangoli dal centro ai vertici
+      else {
+         Point3d ptTriaCen = trTria.GetCentroid() ;
+         int nCenIndex = AddVertex( ptTriaCen) ;
+         int nVertNum = int( vVertOtl.size()) ; 
+         for ( int nStV = 0 ; nStV < nVertNum ; ++ nStV) {
+            int nEnV = ( nStV + 1) % nVertNum ;
+            int nNewId[3] = { nCenIndex, vVertOtl[nStV], vVertOtl[nEnV]} ;
+            int nNewTria = AddTriangle( nNewId, nTFlag) ;
+            if ( nNewTria != SVT_NULL  &&  nNewTria != SVT_DEL)
+               nLastNewTria = nNewTria ;
+         }
+      }
+     // Eventuale aggiustamento per indice di faccetta
+      if ( nFacet != SVT_NULL  &&  nLastNewTria != SVT_NULL)
+         m_vFacet[nFacet] = nLastNewTria ;
+   }
+   
+   return true ;
+}
+
 //----------------------------------------------------------------------------
 static bool
 IsVertex( PNTULIST& PointList, PNTULIST::const_iterator itCurr)
@@ -57,6 +192,9 @@ ChooseGoodStartPoint( PNTULIST& PointList)
   // altrimenti cerco il vertice più vicino
    for ( auto it = next( PointList.begin()) ; it != PointList.end() ; ++it) {
       if ( IsVertex( PointList, it)) {
+        // se ultimo punto non devo fare alcunché
+         if ( next( it) == PointList.end())
+            return false ;
         // cancello ultimo punto ( coincide con primo)
          PointList.pop_back() ;
         // sposto la parte iniziale dei punti alla fine
@@ -160,13 +298,13 @@ AdjustLoop( PNTULIST& PointList, double dMaxEdgeLen, bool& bModif)
 
 //----------------------------------------------------------------------------
 bool
-SurfTriMesh::SimplifyFacets( double dMaxEdgeLen)
+SurfTriMesh::SimplifyFacets( double dMaxEdgeLen, bool bForced)
 {
   // La trimesh deve essere valida
    if ( ! IsValid())
       return false ;
-  // Se la lunghezza massima del lato del triangolo sul bordo della faccia è nulla, non devo fare alcunché
-   if ( dMaxEdgeLen < EPS_SMALL)
+  // Se la lunghezza massima del lato del triangolo sul bordo della faccia è nulla e non forzata, non devo fare alcunché
+   if ( dMaxEdgeLen < EPS_SMALL  &&  ! bForced)
       return true ;
 
   // Recupero il numero delle facce (esegue anche una verifica delle stesse)
@@ -175,13 +313,13 @@ SurfTriMesh::SimplifyFacets( double dMaxEdgeLen)
   // Ciclo sulle facce della mesh per trovare quelle da ritriangolare
    unordered_map< int, pair< PNTVECTOR, INTVECTOR>> FacetMap ;
    for ( int nF = 0 ; nF < nFacetCnt ; ++ nF) {
-
+      
      // Recupero i loop della faccia (il parametro indica la faccia adiacente)
       POLYLINEVECTOR LoopVec ;
       GetFacetLoops( nF, LoopVec) ;
 
      // Ciclo sui loop della faccia
-      bool bToRetriangulate = false ;
+      bool bToRetriangulate = bForced ;
       for ( int nL = 0 ; nL < int( LoopVec.size()) ; ++ nL) {
 
         // Lista dei punti del loop
@@ -189,7 +327,7 @@ SurfTriMesh::SimplifyFacets( double dMaxEdgeLen)
 
         // Mi assicuro che il punto iniziale/finale non sia all'interno di un possibile segmento
          if ( ! ChooseGoodStartPoint( PointList))
-            return false ;
+            continue ;
 
         // Sistemo il loop
          bool bModif = false ;
@@ -199,7 +337,24 @@ SurfTriMesh::SimplifyFacets( double dMaxEdgeLen)
             bToRetriangulate = true ;
       }
 
-     // Se da ritriangolare,
+     // Se non richiesta ritriangolazione dai bordi, verifico se ci sono vertici di triangoli interni (*** disabilitato ***)
+      if ( false  &&  ! bToRetriangulate) {
+        // numero dei triangoli nella faccia
+         INTVECTOR vFacetTria ;
+         GetAllTriaInFacet( nF, vFacetTria) ;
+         int nTriaCnt = int( vFacetTria.size()) ;
+        // numero dei lati di contorno della faccia
+         int nSideCnt = 0 ;
+         for ( int nL = 0 ; nL < int( LoopVec.size()) ; ++ nL)
+            nSideCnt += LoopVec[nL].GetLineNbr() ;
+        // numero dei buchi della faccia
+         int nHoleCnt = int( LoopVec.size()) - 1 ;
+        // dalla formula di Eulero adattata al caso ( nTriaCnt = nSideCnt + 2 * ( nHoleCnt - 1))
+         if ( nTriaCnt != nSideCnt + 2 * ( nHoleCnt - 1))
+            bToRetriangulate = true ;
+      }
+
+     // Se da ritriangolare
       if ( bToRetriangulate) {
         // Eseguo la ritriangolazione della faccia
          PNTVECTOR vPt ;
@@ -214,33 +369,43 @@ SurfTriMesh::SimplifyFacets( double dMaxEdgeLen)
    }
 
   // Ciclo sulle facce da ritriangolare per eliminare i triangoli (nel contempo salvo flag colore)
+   INTVECTOR vDelTria ;
    unordered_map< int, int> ColorMap ;
    for ( auto itF = FacetMap.begin() ; itF != FacetMap.end() ; ++ itF) {
      // Recupero i triangoli della faccia
       INTVECTOR vFacetTria ;
       GetAllTriaInFacet( itF->first, vFacetTria) ;
+      vDelTria.insert( vDelTria.end(), vFacetTria.begin(), vFacetTria.end()) ;
      // Salvo il colore della faccia da flag di un suo triangolo
       ColorMap.emplace( itF->first, m_vTria[m_vFacet[itF->first]].nTFlag) ;
-     // Cancello i triangoli della faccia.
-      for ( int nT : vFacetTria)
-         RemoveTriangle( nT) ;
    }
-
+  // Cancello i triangoli
+   for ( int nT : vDelTria)
+      RemoveTriangle( nT) ;
+   
   // Applico le nuove triangolazioni delle facce
-   for ( auto itFac = FacetMap.begin() ; itFac != FacetMap.end() ; ++ itFac) {
-      const PNTVECTOR& vPt = itFac->second.first ;
-      const INTVECTOR& vTr = itFac->second.second ;
+   for ( auto itF = FacetMap.begin() ; itF != FacetMap.end() ; ++ itF) {
+      const PNTVECTOR& vPt = itF->second.first ;
+      const INTVECTOR& vTr = itF->second.second ;
      // Inserisco i nuovi triangoli
+      bool bFirstTria = true ;
       for ( int n = 0 ; n < int( vTr.size()) - 2 ; n += 3) {
          int nNewId[3] = { AddVertex( vPt[vTr[n]]),
                            AddVertex( vPt[vTr[n + 1]]),
                            AddVertex( vPt[vTr[n + 2]])} ;
-         auto itCol = ColorMap.find( itFac->first) ;
+         auto itCol = ColorMap.find( itF->first) ;
          int nTFlag = ( itCol != ColorMap.end() ? itCol->second : 0) ;
          int nNewTriaId = AddTriangle( nNewId, nTFlag) ;
+         if ( nNewTriaId != SVT_NULL  &&  nNewTriaId != SVT_DEL) {
+            m_vTria[nNewTriaId].nIdFacet = itF->first ;
+            if ( bFirstTria) {
+               m_vFacet[itF->first] = nNewTriaId ;
+               bFirstTria = false ;
+            }
+         }
       }
    }
-
+   
   // dichiaro necessità ricalcolo della grafica e di hashgrids3d
    m_OGrMgr.Reset() ;
    ResetHashGrids3d() ;
@@ -248,3 +413,192 @@ SurfTriMesh::SimplifyFacets( double dMaxEdgeLen)
   // Eseguo aggiustamenti
    return ( AdjustVertices()  &&  DoCompacting()) ;
 }
+
+//----------------------------------------------------------------------------
+bool
+SurfTriMesh::AddChainToChain( const Chain& ChainToAdd, PNTVECTOR& OrigChain)
+{
+  // Se la catena da aggiungere è vuota, non devo fare alcunchè
+   if ( ChainToAdd.size() == 0)
+      return true ;
+  // Se la catena originale è vuota, non è possibile aggiungere nulla
+   if ( OrigChain.size() == 0)
+      return false ;
+  // Se la catena originale è chiusa non posso aggiungere nulla
+   int nLastOrig = max( int( OrigChain.size()) - 1, 0) ;
+   if ( AreSamePointApprox( OrigChain[0], OrigChain[nLastOrig]))
+      return false ;
+   int nLastToAdd = max( int( ChainToAdd.size()) - 1, 0) ;  
+   if ( AreSamePointApprox( OrigChain[nLastOrig], ChainToAdd[0].ptSt)) {
+      for ( int nPt = 1 ; nPt <= nLastToAdd ; ++ nPt) {
+         if ( nPt == nLastToAdd) {
+            if ( ! AreSamePointApprox(OrigChain[0], ChainToAdd[nPt].ptSt))
+               OrigChain.emplace_back( ChainToAdd[nPt].ptSt) ;
+         }
+         else if ( nPt == 1) {
+            if ( ! AreSamePointApprox( OrigChain[nLastOrig], ChainToAdd[nPt].ptSt))
+               OrigChain.emplace_back( ChainToAdd[nPt].ptSt) ;
+         }
+         else
+            OrigChain.emplace_back( ChainToAdd[nPt].ptSt) ;
+      }
+      return true ;
+   }
+   else
+      return false ;
+}
+
+//----------------------------------------------------------------------------
+// Una faccia di una trimesh ha una sola componente connessa, con un loop esterno e possibili loop interni
+bool
+SurfTriMesh::DistPointFacet( const Point3d& ptP, const POLYLINEVECTOR& vPolyVec, double& dPointFacetDist)
+{
+  // Verifico la presenza del loop esterno
+   if ( vPolyVec.size() < 1)
+      return false ;
+  // Proietto il punto sul piano della faccia, utilizzando il loop esterno
+   Plane3d plPlane ;
+   double dArea ;
+   if ( ! vPolyVec[0].IsClosedAndFlat( plPlane, dArea))
+      return false ;
+   double dDistPtPl = DistPointPlane( ptP, plPlane) ;
+   Point3d ptProjP = ptP + dDistPtPl * plPlane.GetVersN() ;
+  // Verifico se il punto proiettato è esterno al loop esterno
+   int nPtOut = -1 ;
+   if ( ! IsPointInsidePolyLine( ptProjP, vPolyVec[0], EPS_SMALL))
+      nPtOut = 0 ;
+  // Verifico se il punto proiettato è interno ai loop interni (quindi esterno alla faccia)
+   for ( int nLoop = 1 ; nLoop < int( vPolyVec.size()) && nPtOut < 0 ; ++ nLoop) {
+      Plane3d plPlane ;
+      double dArea ;
+      if ( ! vPolyVec[nLoop].IsClosedAndFlat( plPlane, dArea))
+         return false ;
+      if ( IsPointInsidePolyLine( ptProjP, vPolyVec[nLoop], EPS_SMALL))
+         nPtOut = nLoop ;
+   }
+  // Se il punto si proietta sulla faccia, la distanza dalla faccia coincide con quella dal piano
+   if ( nPtOut < 0) {
+      dPointFacetDist = abs( dDistPtPl) ;
+      return true ;
+   }
+  // Altrimenti calcolo la minima distanza del punto dalla polilinea del contorno a cui è esterno
+   double dDist ;
+   if ( DistPointPolyLine( ptP, vPolyVec[nPtOut], dDist)) {
+      dPointFacetDist = dDist ; 
+      return true ;
+   }
+   return false ;
+}
+
+//----------------------------------------------------------------------------
+bool
+SurfTriMesh::ChangeStart( const Point3d& ptNewStart, PNTVECTOR& Loop)
+{
+  // Cerco il tratto del loop chiuso più vicino al punto
+   int nMinSeg = - 1 ;  
+   double dMinSqDinst = DBL_MAX ;
+   for ( int nPt = 0 ; nPt < int( Loop.size()) ; ++ nPt) {
+     // Estremi del segmento corrente del loop
+      Point3d ptSegSt = Loop[nPt] ;
+      Point3d ptSegEn = Loop[( nPt + 1) % int( Loop.size())] ;
+     // Distanza del punto dal segmento del loop
+      DistPointLine dDistCalc( ptNewStart, ptSegSt, ptSegEn) ;
+      double dSqDist ;
+      dDistCalc.GetSqDist( dSqDist) ;
+      if ( dSqDist < dMinSqDinst) {
+         dMinSqDinst = dSqDist ;
+         nMinSeg = nPt ;
+      }
+   }
+  // Se il punto non sta sul loop, errore
+   if ( dMinSqDinst > SQ_EPS_SMALL)
+      return false ;
+  // Verifico che il punto stia su un vertice, in tal caso non devo fare nulla
+   bool bOnStart = AreSamePointApprox( Loop[nMinSeg], ptNewStart) ;
+   bool bOnEnd = AreSamePointApprox( Loop[( nMinSeg + 1) % int( Loop.size())], ptNewStart) ;
+   if ( bOnStart  ||  bOnEnd) {
+      if ( bOnEnd) {
+         ++ nMinSeg ;
+         if ( nMinSeg % int( Loop.size()) == 0)
+            return true ;
+      }
+      PNTVECTOR vTempVec ;
+      for ( int nPt = 0 ; nPt < nMinSeg ; ++ nPt)
+         vTempVec.emplace_back( Loop[nPt]) ;
+      int nSize = int( Loop.size()) ;
+      for ( int nPt = 0 ; nPt < nSize - nMinSeg ; ++ nPt) {
+         Loop[nPt] = Loop[nPt + nMinSeg] ;
+      }
+      for ( int nPt = 0 ; nPt < int( vTempVec.size()) ; ++ nPt) {
+         Loop[nPt + nSize - nMinSeg] = vTempVec[nPt] ;
+      }
+      return true ;
+   }
+  // Ridimensiono il loop
+   Loop.resize( Loop.size() + 1) ;
+  // Copio i primi punti
+   PNTVECTOR LoopTemp ;
+   for ( int nPt = 0 ; nPt <= nMinSeg ; ++ nPt)
+      LoopTemp.emplace_back( Loop[nPt]) ;
+  // Aggiungo il nuovo punto all'inizio
+   Loop[0] = ptNewStart ;
+  // Sposto gli ultimi in testa
+  int nLastPointNum = int( Loop.size()) - 1 - nMinSeg ;
+   for ( int nPt = 1 ; nPt <= nLastPointNum ; ++ nPt) {
+      Loop[nPt] = Loop[nPt + nMinSeg] ;
+   }
+  // Porto i primi in fondo 
+   for ( int nPt = 0 ; nPt < int( LoopTemp.size()) ; ++ nPt) {
+      Loop[nPt + nLastPointNum] = LoopTemp[nPt] ;
+   }
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+bool
+SurfTriMesh::SplitAtPoint( const Point3d& ptStop, const PNTVECTOR& Loop, PNTVECTOR& Loop1, PNTVECTOR& Loop2)
+{
+  // Cerco il tratto del loop chiuso più vicino al punto
+   int nMinSeg = -1 ;
+   double dMinSqDinst = DBL_MAX ;
+   for ( int nPt = 0 ; nPt < int( Loop.size()) ; ++ nPt) {
+     // Estremi del segmento corrente del loop
+      Point3d ptSegSt = Loop[nPt] ;
+      Point3d ptSegEn = Loop[( nPt + 1) % int(Loop.size())] ;
+     // Distanza del punto dal segmento del loop
+      DistPointLine dDistCalc( ptStop, ptSegSt, ptSegEn) ;
+      double dSqDist ;
+      dDistCalc.GetSqDist( dSqDist) ;
+      if ( dSqDist < dMinSqDinst) {
+         dMinSqDinst = dSqDist ;
+         nMinSeg = nPt ;
+      }
+   }
+  // Se il punto non sta sul loop, errore
+   if ( dMinSqDinst > SQ_EPS_SMALL)
+      return false ;
+  // Verifico che il punto stia su un vertice, in tal caso non devo aggiungerlo
+   bool bFirst = AreSamePointApprox( Loop[nMinSeg], ptStop) ;
+   bool bLast = AreSamePointApprox( Loop[( nMinSeg + 1) % int( Loop.size())], ptStop) ;
+  // Se il punto è sul vertice finale del segmento, aggiungo il vertice alla lista da inglobare al primo loop
+   if ( bLast)
+      ++ nMinSeg ;
+  // Inglobo fino a nSeg nel primo loop
+   for ( int nPt = 0 ; nPt <= nMinSeg  &&  nPt < int( Loop.size()) ; ++ nPt)
+      Loop1.emplace_back( Loop[nPt]) ;
+  // Se il punto è interno al segmento, lo inglobo in entrambi i loop
+   if ( ! ( bFirst  ||  bLast)) {
+      Loop1.emplace_back( ptStop) ;
+      Loop2.emplace_back( ptStop) ;
+   }
+   else {
+      if ( nMinSeg != int( Loop.size()) )
+         Loop2.emplace_back( Loop[nMinSeg]) ;
+   }
+  // Inglobo gli ultimi vertici in Loop2
+   for ( int nPt = nMinSeg + 1 ; nPt < int( Loop.size()) ; ++ nPt)
+      Loop2.emplace_back( Loop[nPt]) ;
+   Loop2.emplace_back( Loop[0]) ;
+
+   return true ;
+}

Tool.cpp

@@ -100,9 +100,13 @@ Tool::SetStdTool( const string& sToolName, double dH, double dR, double dCornR,
       Point3d pt3( m_dRadius, - m_dHeight, 0) ;
       Point3d pt4( 0, - m_dHeight, 0) ;
       m_Outline.AddPoint( pt0) ;
-      m_Outline.AddLine( pt1);
+      m_Outline.AddLine( pt1) ;
+      m_Outline.SetCurveTempProp( 0, 1, 1) ;
       m_Outline.AddLine( pt3) ;
+      m_Outline.SetCurveTempProp( 1, 1, 1) ;
       m_Outline.AddLine( pt4) ;
+      m_Outline.SetCurveTempProp( 2, 1, 1) ;
+      m_Outline.SetTempProp( 1, 1) ;
    }
   // utensile naso di toro
    else if ( dCornR < dR - EPS_SMALL) {
@@ -120,10 +124,15 @@ Tool::SetStdTool( const string& sToolName, double dH, double dR, double dCornR,
       Point3d pt3( m_dTipRadius, - m_dHeight, 0) ;
       Point3d pt4( 0, - m_dHeight, 0) ;
       m_Outline.AddPoint( pt0) ;
-      m_Outline.AddLine( pt1);
-      m_Outline.AddLine( pt2);
+      m_Outline.AddLine( pt1) ;
+      m_Outline.SetCurveTempProp( 0, 1, 1) ;
+      m_Outline.AddLine( pt2) ;
+      m_Outline.SetCurveTempProp( 1, 1, 1) ;
       m_Outline.AddArcTg( pt3) ;
+      m_Outline.SetCurveTempProp( 2, 1, 1) ;
       m_Outline.AddLine( pt4) ;
+      m_Outline.SetCurveTempProp( 3, 1, 1) ;
+      m_Outline.SetTempProp( 1, 1) ;
      // se da approosimare
       if ( m_bApproxWithLines)
         return SetGenTool( sToolName, &m_Outline, nToolNum) ;
@@ -143,9 +152,13 @@ Tool::SetStdTool( const string& sToolName, double dH, double dR, double dCornR,
       Point3d pt2( m_dRadius, - m_dHeight + m_dTipHeight, 0) ;
       Point3d pt4( 0, - m_dHeight, 0) ;
       m_Outline.AddPoint( pt0) ;
-      m_Outline.AddLine( pt1);
-      m_Outline.AddLine( pt2);
+      m_Outline.AddLine( pt1) ;
+      m_Outline.SetCurveTempProp( 0, 1, 1) ;
+      m_Outline.AddLine( pt2) ;
+      m_Outline.SetCurveTempProp( 1, 1, 1) ;
       m_Outline.AddArcTg( pt4) ;
+      m_Outline.SetCurveTempProp( 2, 1, 1) ;
+      m_Outline.SetTempProp( 1, 1) ;
    }
   // impossibile
    else
@@ -201,9 +214,14 @@ Tool::SetAdvTool( const string& sToolName, double dH, double dR,
      // profilo
       m_Outline.AddPoint( pt0) ;
       m_Outline.AddLine( pt1) ;
+      m_Outline.SetCurveTempProp( 0, 1, 1) ;
       m_Outline.AddLine( pt2) ;
+      m_Outline.SetCurveTempProp( 1, 1, 1) ;
       m_Outline.AddLine( Point3d( m_dTipRadius, - m_dHeight, 0)) ;
+      m_Outline.SetCurveTempProp( 2, 1, 1) ;
       m_Outline.AddLine( pt5) ;
+      m_Outline.SetCurveTempProp( 3, 1, 1) ;
+      m_Outline.SetTempProp( 1, 1) ;
      // eventuali sistemazioni per altezza tagliente
       if ( ModifyForCutterHeight())
          return SetGenTool( sToolName, &m_Outline, nToolNum) ;
@@ -213,7 +231,7 @@ Tool::SetAdvTool( const string& sToolName, double dH, double dR,
 
   // Altrimenti utensile generico
 
-  // se Tip a punta ( TipRadius è raggio teorico della parte finale dell'utensile senza raccordo)
+  // se Tip a punta ( TipRadius � raggio teorico della parte finale dell'utensile senza raccordo)
    if ( m_dTipRadius < m_dRadius) {
      // se punta a sfera
       if ( m_dTipRadius < EPS_SMALL) {
@@ -242,10 +260,16 @@ Tool::SetAdvTool( const string& sToolName, double dH, double dR,
         // creazione curva composita
          m_Outline.AddPoint( pt0) ;
          m_Outline.AddLine( pt1) ;
+         m_Outline.SetCurveTempProp( 0, 1, 1) ;
          m_Outline.AddLine( pt2) ;
+         m_Outline.SetCurveTempProp( 1, 1, 1) ;
          m_Outline.AddLine( pt3) ;
+         m_Outline.SetCurveTempProp( 2, 1, 1) ;
          m_Outline.AddCurve( cvArc) ;
+         m_Outline.SetCurveTempProp( 3, 1, 1) ;
          m_Outline.AddLine( pt5) ;
+         m_Outline.SetCurveTempProp( 4, 1, 1) ;
+         m_Outline.SetTempProp( 1, 1) ;
       }
      // altrimenti punta a naso di toro
       else {
@@ -271,16 +295,22 @@ Tool::SetAdvTool( const string& sToolName, double dH, double dR,
         // creazione curva composita
          m_Outline.AddPoint( pt0) ;
          m_Outline.AddLine( pt1) ;
+         m_Outline.SetCurveTempProp( 0, 1, 1) ;
          m_Outline.AddLine( pt2) ;
+         m_Outline.SetCurveTempProp( 1, 1, 1) ;
          m_Outline.AddLine( pt3) ;
+         m_Outline.SetCurveTempProp( 2, 1, 1) ;
          m_Outline.AddCurve( Release( pArc)) ;
+         m_Outline.SetCurveTempProp( 3, 1, 1) ;
          m_Outline.AddLine( pt5) ;
+         m_Outline.SetCurveTempProp( 4, 1, 1) ;
+         m_Outline.SetTempProp( 1, 1) ;
       }
    }
 
-   // altrimenti Tip a coda di rondine ( TipRadius è raggio misurabile della parte finale dell'utensile)
+   // altrimenti Tip a coda di rondine ( TipRadius � raggio misurabile della parte finale dell'utensile)
    else {
-     // il raggio della punta non può essere inferiore al raggio corner
+     // il raggio della punta non pu� essere inferiore al raggio corner
       if ( m_dTipRadius < m_dRCorner)
          return false ;
      // Assegno il raggio di riferimento
@@ -305,10 +335,16 @@ Tool::SetAdvTool( const string& sToolName, double dH, double dR,
      // creazione curva composita
       m_Outline.AddPoint( pt0) ;
       m_Outline.AddLine( pt1) ;
+      m_Outline.SetCurveTempProp( 0, 1, 1) ;
       m_Outline.AddLine( pt2) ;
+      m_Outline.SetCurveTempProp( 1, 1, 1) ;
       m_Outline.AddLine( pt3) ;
+      m_Outline.SetCurveTempProp( 2, 1, 1) ;
       m_Outline.AddCurve( cvArc) ;
+      m_Outline.SetCurveTempProp( 3, 1, 1) ;
       m_Outline.AddLine( pt5) ;
+      m_Outline.SetCurveTempProp( 4, 1, 1) ;
+      m_Outline.SetTempProp( 1, 1) ;
    }
 
   // eventuali sistemazioni per altezza tagliente
@@ -321,7 +357,7 @@ Tool::SetAdvTool( const string& sToolName, double dH, double dR,
 bool
 Tool::ModifyForCutterHeight( void)
 {
-  // Se altezza tagliente non definita o superiore alla altezza utensile non devo fare alcunché
+  // Se altezza tagliente non definita o superiore alla altezza utensile non devo fare alcunch�
    if ( m_dCutterHeight < EPS_SMALL  ||  m_dCutterHeight > m_dHeight - EPS_SMALL)
       return false ;
   // quota di taglio
@@ -340,6 +376,8 @@ Tool::ModifyForCutterHeight( void)
       m_Outline.TrimStartAtParam( dU) ;
       m_Outline.AddLine( Point3d( 0, dYtrim, 0), false) ;
       m_Outline.AddLine( ORIG, false) ;
+      m_Outline.SetCurveTempProp( 0, 1, 1) ;
+      m_Outline.SetCurveTempProp( 1, 1, 1) ;
       return true ;
    } 
    return false ;
@@ -375,22 +413,36 @@ Tool::SetSawTool( const string& sToolName, double dH, double dR,
         // creazione profilo
          m_Outline.AddPoint( pt0) ;
          m_Outline.AddLine( pt1) ;
+         m_Outline.SetCurveTempProp( 0, 1, 1) ;
          m_Outline.AddLine( pt2) ;
+         m_Outline.SetCurveTempProp( 1, 1, 1) ;
          m_Outline.AddLine( pt3) ;
+         m_Outline.SetCurveTempProp( 2, 1, 1) ;
          m_Outline.AddLine( pt4) ;
+         m_Outline.SetCurveTempProp( 3, 1, 1) ;
          m_Outline.AddLine( pt5) ;
+         m_Outline.SetCurveTempProp( 4, 1, 1) ;
+         m_Outline.SetTempProp( 1, 1) ;
       }
      // altrimenti con raggio corner
       else {
         // creazione profilo
          m_Outline.AddPoint( pt0) ;
          m_Outline.AddLine( pt1) ;
+         m_Outline.SetCurveTempProp( 0, 1, 1) ;
          m_Outline.AddLine( pt2) ;
+         m_Outline.SetCurveTempProp( 1, 1, 1) ;
          m_Outline.AddLine( pt3 - X_AX * dCornR) ;
+         m_Outline.SetCurveTempProp( 0, 1, 1) ;
          m_Outline.AddArcTg( pt3 - Y_AX * dCornR) ;
+         m_Outline.SetCurveTempProp( 2, 1, 1) ;
          m_Outline.AddLine( pt4 + Y_AX * dCornR) ;
+         m_Outline.SetCurveTempProp( 3, 1, 1) ;
          m_Outline.AddArcTg( pt4 - X_AX * dCornR) ;
+         m_Outline.SetCurveTempProp( 4, 1, 1) ;
          m_Outline.AddLine( pt5) ;
+         m_Outline.SetCurveTempProp( 5, 1, 1) ;
+         m_Outline.SetTempProp( 1, 1) ;
       }
    }
   // altrimenti senza gambo
@@ -405,18 +457,28 @@ Tool::SetSawTool( const string& sToolName, double dH, double dR,
         // creazione profilo
          m_Outline.AddPoint( pt0) ;
          m_Outline.AddLine( pt3) ;
+         m_Outline.SetCurveTempProp( 0, 1, 1) ;
          m_Outline.AddLine( pt4) ;
+         m_Outline.SetCurveTempProp( 0, 1, 1) ;
          m_Outline.AddLine( pt5) ;
+         m_Outline.SetCurveTempProp( 0, 1, 1) ;
+         m_Outline.SetTempProp( 1, 1) ;
       }
      // altrimenti con raggio corner
       else {
         // creazione profilo
          m_Outline.AddPoint( pt0) ;
          m_Outline.AddLine( pt3 - X_AX * dCornR) ;
+         m_Outline.SetCurveTempProp( 0, 1, 1) ;
          m_Outline.AddArcTg( pt3 - Y_AX * dCornR) ;
+         m_Outline.SetCurveTempProp( 0, 1, 1) ;
          m_Outline.AddLine( pt4 + Y_AX * dCornR) ;
+         m_Outline.SetCurveTempProp( 0, 1, 1) ;
          m_Outline.AddArcTg( pt4 - X_AX * dCornR) ;
+         m_Outline.SetCurveTempProp( 0, 1, 1) ;
          m_Outline.AddLine( pt5) ;
+         m_Outline.SetCurveTempProp( 0, 1, 1) ;
+         m_Outline.SetTempProp( 1, 1) ;
       }
    }
   // Assegno il raggio di riferimento
@@ -436,17 +498,30 @@ Tool::SetGenTool( const string& sToolName, const ICurveComposite* pToolOutline,
    if ( pToolOutline != &m_Outline)
       m_Outline.Clear() ;
    m_ArcLineApprox.Clear() ;
-
+                                                           
   // Copio il profilo e garantisco sia di soli archi e rette (converto eventuali curve di Bezier)
    if ( ! m_Outline.CopyFrom( pToolOutline)  ||
         ! m_Outline.ArcsBezierCurvesToArcsPerpExtr( m_dLinTol, m_dAngTolDeg))
       return false ;
 
-  // Ciclo sulle curve componenti 
+  // Valuto se tutte le curve tagliano
+   m_ArcLineApprox.SetTempProp( 1, 1) ;
    const ICurve* pCurve = m_Outline.GetFirstCurve() ;
+   while ( pCurve != nullptr  &&  m_ArcLineApprox.GetTempProp( 1) == 1) {
+      int nCutTempProp = pCurve->GetTempProp( 1) ;
+      if ( nCutTempProp == 0) {
+         m_ArcLineApprox.SetTempProp( 0, 1) ;
+      }
+      pCurve = m_Outline.GetNextCurve() ;
+   }
+
+  // Ciclo sulle curve componenti 
+   pCurve = m_Outline.GetFirstCurve() ;
    while ( pCurve != nullptr) {
 
-     // Se la curva è un arco ed è richiesto la verifica per l'approssimazione,
+      int nCutTempProp = pCurve->GetTempProp( 1) ;
+
+     // Se la curva � un arco ed � richiesto la verifica per l'approssimazione,
      // verifico se approssimarlo
       if ( m_bApproxWithLines  &&  pCurve->GetType() == CRV_ARC) {
 
@@ -454,11 +529,11 @@ Tool::SetGenTool( const string& sToolName, const ICurveComposite* pToolOutline,
          Point3d ptO = GetBasicCurveArc( pCurve)->GetCenter() ;
          double dRadius = GetBasicCurveArc( pCurve)->GetRadius() ;
 
-        // Se il centro è fuori dall'asse devo approssimare
+        // Se il centro � fuori dall'asse devo approssimare
          bool bCurrApprox = ( abs( ptO.x) > EPS_SMALL) ;
 
         // Se una delle altre curve dista dal centro meno del raggio, devo approssimare
-         for ( int nI = 0 ; ! bCurrApprox  &&  nI < m_Outline.GetCurveCount() ; ++ nI) {
+         for ( int nI = 0 ; ! bCurrApprox  &&  nI < m_Outline.GetCurveCount()  &&  m_ArcLineApprox.GetTempProp( 1) == 1 ; ++ nI) {
             const ICurve* pOtherCrv = m_Outline.GetCurve( nI) ; 
             if ( pOtherCrv != pCurve) {
                DistPointCurve CalcDist( ptO, *pOtherCrv) ;
@@ -493,6 +568,7 @@ Tool::SetGenTool( const string& sToolName, const ICurveComposite* pToolOutline,
                m_vArcNormals.emplace_back( vtExtN) ;
                int nCvCount = m_ArcLineApprox.GetCurveCount() ;
                m_ArcLineApprox.SetCurveTempProp( nCvCount - 1, int( m_vArcNormals.size()) - 1) ;
+               m_ArcLineApprox.SetCurveTempProp( nCvCount - 1, nCutTempProp, 1) ;
             }
          }
         // altrimenti lo aggiungo semplicemente 
@@ -500,16 +576,20 @@ Tool::SetGenTool( const string& sToolName, const ICurveComposite* pToolOutline,
             m_ArcLineApprox.AddCurve( *pCurve, true) ;
             int nCvCount = m_ArcLineApprox.GetCurveCount() ;
             m_ArcLineApprox.SetCurveTempProp( nCvCount - 1, - 1) ;
+            m_ArcLineApprox.SetCurveTempProp( nCvCount - 1, nCutTempProp, 1) ;
          }
       }
-     // altrimenti è segmento e lo aggiungo semplicemente
-      else
+     // altrimenti � segmento e lo aggiungo semplicemente
+      else {
          m_ArcLineApprox.AddCurve( *pCurve, true) ;
+         int nCvCount = m_ArcLineApprox.GetCurveCount() ;
+         m_ArcLineApprox.SetCurveTempProp( nCvCount - 1, nCutTempProp, 1) ;
+      }
 
       pCurve = m_Outline.GetNextCurve() ;
    }
 
-  // Il profilo dell'utensile deve stare nel 1° e 4° quadrante del piano XY
+  // Il profilo dell'utensile deve stare nel 1� e 4� quadrante del piano XY
    BBox3d Bounding ;
    m_Outline.GetLocalBBox( Bounding) ;
    if ( Bounding.GetMin().x < - 10 * EPS_SMALL)
@@ -521,7 +601,7 @@ Tool::SetGenTool( const string& sToolName, const ICurveComposite* pToolOutline,
   // Assegno le dimensioni dell'utensile
    m_dHeight = - Bounding.GetMin().y ;
    m_dRadius = Bounding.GetMax().x ;
-  // Assegno il raggio di riferimento se non già assegnato
+  // Assegno il raggio di riferimento se non gi� assegnato
    if ( m_dRefRadius < EPS_SMALL)
       m_dRefRadius = 0.25 * m_dRadius ;
  
@@ -579,3 +659,86 @@ Tool::SetChiselTool( const string& sToolName, double dH, double dW, double dTh,
 
    return true ;
 }
+
+//----------------------------------------------------------------------------
+bool
+Tool::SetAdditiveTool( const std::string& sToolName, double dH, double dR, double dRC, int nToolNum)
+{
+   // Impostazioni generali
+   m_sName = sToolName ;
+   m_nCurrentNum = nToolNum ;
+   m_nType = UNDEF ;
+   m_Outline.Clear() ;
+   m_ArcLineApprox.Clear() ;
+
+   // verifica sulle minime dimensioni globali
+   if ( dH < EPS_SMALL  ||  dR < EPS_SMALL  ||  dRC < - EPS_SMALL)
+      return false ;
+
+   m_dHeight = dH ;
+   m_dRadius = dR ;
+   m_dRCorner = dRC ;
+   m_dTipHeight = 0 ;
+   m_dTipRadius = 0 ;
+   m_dRefRadius = 0 ;
+   m_dCutterHeight = dH ;
+
+   bool bToolDefined = true ;
+
+   double dSquareCornerRadProj = m_dRCorner * m_dRCorner - 0.25 * m_dHeight * m_dHeight ;
+   // Utensile sfiancato
+   if ( dSquareCornerRadProj > 0) {
+      double dCenX = m_dRadius - m_dRCorner ;
+      double dCylRad = dCenX + sqrt( dSquareCornerRadProj) ;
+      // Utensile mal definito
+      if ( dCylRad < EPS_SMALL)
+         return false ;
+      // Profilo
+      m_Outline.AddPoint( Point3d( 0, 0, 0)) ;
+      m_Outline.AddLine( Point3d( dCylRad, 0, 0)) ;
+      m_Outline.SetCurveTempProp( 0, 1, 1) ;
+      CurveArc cvArc ;
+      cvArc.SetC2P( Point3d( dCenX, - 0.5 * m_dHeight, 0), Point3d( dCylRad, 0, 0), Point3d( dCylRad, - m_dHeight, 0)) ;
+      m_Outline.AddCurve( cvArc) ;
+      m_Outline.SetCurveTempProp( 1, 1, 1) ;
+      m_Outline.AddLine( Point3d( 0, - m_dHeight, 0)) ;
+      m_Outline.SetCurveTempProp( 2, 1, 1) ;
+      m_Outline.SetTempProp( 1, 1) ;
+      bToolDefined = SetGenTool( sToolName, &m_Outline, nToolNum) ;
+   }
+   // Utensile cilindrico con eventuale raggio corner 
+   else {
+      // Utensile mal definito
+      if ( m_dRadius - m_dRCorner < 0)
+         return false ;
+      // Utensile sferico
+      else if ( m_dRadius - m_dRCorner < EPS_SMALL)
+         ;
+      // Raggio corner nullo: cilindro
+      else if ( m_dRCorner < EPS_SMALL) {
+         ;
+      }
+      else {
+         // Profilo
+         m_Outline.AddPoint( Point3d( 0, 0, 0)) ;
+         m_Outline.AddLine( Point3d( m_dRadius - m_dRCorner, 0, 0)) ;
+         m_Outline.SetCurveTempProp( 0, 1, 1) ;
+         CurveArc cvArc ;
+         cvArc.SetC2P( Point3d( m_dRadius - m_dRCorner, - m_dRCorner, 0), Point3d( m_dRadius - m_dRCorner, 0, 0), Point3d( m_dRadius, - m_dRCorner, 0)) ;
+         m_Outline.AddCurve( cvArc) ;
+         m_Outline.SetCurveTempProp( 1, 1, 1) ;
+         m_Outline.AddLine( Point3d( m_dRadius, - m_dHeight + m_dRCorner, 0)) ;
+         m_Outline.SetCurveTempProp( 2, 1, 1) ;
+         cvArc.SetC2P( Point3d( m_dRadius - m_dRCorner, - m_dHeight + m_dRCorner, 0), Point3d( m_dRadius, - m_dHeight + m_dRCorner, 0), Point3d( m_dRadius - m_dRCorner, - m_dHeight, 0)) ;
+         m_Outline.AddCurve( cvArc) ;
+         m_Outline.SetCurveTempProp( 3, 1, 1) ;
+         m_Outline.AddLine( Point3d( 0, - m_dHeight, 0)) ;
+         m_Outline.SetCurveTempProp( 4, 1, 1) ;
+         m_Outline.SetTempProp( 1, 1) ;
+         bToolDefined = SetGenTool( sToolName, &m_Outline, nToolNum) ;
+      }
+   } 
+
+   m_nType = ADDITIVE ;
+   return bToolDefined ;
+}

Tool.h

@@ -16,7 +16,7 @@
 #include "CurveComposite.h"
 
 //----------------------------------------------------------------------------
-class Tool
+class Tool  
 {
    public :
       Tool( bool bApproxWithLines = false) ;
@@ -33,6 +33,7 @@ class Tool
       bool SetGenTool( const std::string& sToolName, const ICurveComposite* pToolOutline, int nToolNum) ;
       bool SetMortiserTool( const std::string& sToolName, double dH, double dW, double dTh, double dRc, int nToolNum) ;
       bool SetChiselTool( const std::string& sToolName, double dH, double dW, double dTh, int nToolNum) ;
+      bool SetAdditiveTool( const std::string& sToolName, double dH, double dR, double dRC, int nToolNum) ;
       bool SetToolNum( int nToolNum)
                { m_nCurrentNum = nToolNum ; return true ; }
       int GetType() const
@@ -61,6 +62,8 @@ class Tool
                { return ( m_ArcLineApprox.GetCurveCount() == 0 ? m_Outline : m_ArcLineApprox) ; }
       const VCT3DVECTOR& GetArcNormalVec( void) const
                { return m_vArcNormals ; }
+      bool GetCuttingFlag() const
+               { return ( m_nType == GEN ? m_ArcLineApprox.GetTempProp( 1) == 1 : true ) ; }
 
    public :
       enum ToolType { UNDEF = 0,         // Utensile indefinito
@@ -70,7 +73,8 @@ class Tool
                       BULLNOSEMILL = 4,  // Naso di toro
                       CONEMILL = 5,      // Con parte terminale conica
                       MORTISER = 6,      // Mortasatrice
-                      CHISEL = 7} ;      // Scalpello
+                      CHISEL = 7,        // Scalpello
+                      ADDITIVE = 8} ;    // Additivo
 
    private :
        bool ModifyForCutterHeight( void) ;

Triangulate.cpp

@@ -1,7 +1,7 @@
 //----------------------------------------------------------------------------
-//  EgalTech 2014-2014
+//  EgalTech 2014-2021
 //----------------------------------------------------------------------------
-// File : Triangulate.cpp       Data : 23.06.14          Versione : 1.5f6
+// File : Triangulate.cpp       Data : 18.12.21          Versione : 2.3l
 // Contenuto : Implementazione della classe Triangulate.
 //
 //
@@ -16,6 +16,7 @@
 #include "DllMain.h"
 #include "Triangulate.h"
 #include "ProjPlane.h"
+#include "earcut.hpp"
 #include "/EgtDev/Include/EGkPolyLine.h"
 #include "/EgtDev/Include/EGkPlane3d.h"
 #include "/EgtDev/Include/EGkStringUtils3d.h"
@@ -29,6 +30,9 @@ enum EarStatus{ EAS_NULL = -1, EAS_NO = 0, EAS_OK = 1} ;
 //----------------------------------------------------------------------------
 static bool ChangeStartPntVector( int nNewStart, PNTVECTOR& vPi) ;
 
+//----------------------------------------------------------------------------
+static bool FORCE_EARCUT_HPP = false ;
+
 //----------------------------------------------------------------------------
 // In : PolyLine
 // Out : PNTVECTOR (Point3d Vector) : points of the polyline 
@@ -37,6 +41,9 @@ static bool ChangeStartPntVector( int nNewStart, PNTVECTOR& vPi) ;
 bool
 Triangulate::Make( const PolyLine& PL, PNTVECTOR& vPt, INTVECTOR& vTr)
 {
+  // pulisco i vettori di ritorno
+   vPt.clear() ;
+   vTr.clear() ;
   // verifico che la polilinea sia chiusa e piana e calcolo il piano medio del poligono
    double dArea ;
    Plane3d plPlane ;
@@ -46,8 +53,13 @@ Triangulate::Make( const PolyLine& PL, PNTVECTOR& vPt, INTVECTOR& vTr)
    bool bCCW ;
    if ( ! CalcProjPlane( plPlane.GetVersN(), m_nPlane, bCCW))
       return false ;
+
+  // forzo esecuzione triangolazione con earcut.hpp
+   if ( FORCE_EARCUT_HPP) {
+      return MakeByEC_HPP( PL, bCCW, vPt, vTr) ;
+   }
+
   // riempio il vettore con i vertici del poligono da triangolare
-   vPt.clear() ;
    vPt.reserve( PL.GetPointNbr() - 1) ;
    // salto il primo punto (coincide con l'ultimo)
    Point3d ptP ;
@@ -59,18 +71,12 @@ Triangulate::Make( const PolyLine& PL, PNTVECTOR& vPt, INTVECTOR& vTr)
   // se non CCW inverto il vettore dei punti
    if ( ! bCCW)
       reverse( vPt.begin(), vPt.end()) ;
-  // creo il vettore degli indici del Poligono
-   INTVECTOR vPol ;
-   int n = int( vPt.size()) ;
-   vPol.reserve( n) ;
-   for ( int i = 0 ; i < n ; ++ i)
-      vPol.push_back( i) ;
 
   // eseguo la triangolazione
-   if ( ! MakeByEC2( vPt, vPol, vTr)  &&
-        ! MakeByEC3( vPt, vPol, vTr)) {
+   if ( ! MakeByEC2( vPt, vTr)  &&
+        ! MakeByEC3( vPt, vTr)) {
       LOG_ERROR( GetEGkLogger(), "Error in MakeByEC23(1)")
-      return false ;
+      return MakeByEC_HPP( PL, bCCW, vPt, vTr) ;
    }
 
   // se era CW, devo invertire il senso dei triangoli
@@ -115,6 +121,12 @@ Triangulate::Make( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr)
            ! AreOppositeVectorApprox( plExtPlane.GetVersN(), plPlane.GetVersN()))
          return false ;
    }
+
+  // forzo esecuzione triangolazione con earcut.hpp
+   if ( FORCE_EARCUT_HPP) {
+      return MakeByEC_HPP( vPL, bCCW, vPt, vTr) ;
+   }
+
   // se non CCW inverto tutte le polilinee
    if ( ! bCCW) {
       for ( int i = 0 ; i < int( vPL.size()) ; ++i)
@@ -122,8 +134,7 @@ Triangulate::Make( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr)
    }
   // calcolo ordine decrescente su Xmax o Ymax o Zmax secondo m_nPlane per le curve interne
    INTVECTOR vOrd ;
-   if ( ! SortInternalLoops( vPL, vOrd))
-      return false ;
+   bool bOk = SortInternalLoops( vPL, vOrd) ;
   // riempio il vettore con i punti dei poligoni da triangolare
    // calcolo spazio totale e spazio massimo per un percorso interno
    int nTot = vPL[0].GetPointNbr() - 1 ;
@@ -136,22 +147,27 @@ Triangulate::Make( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr)
    // riservo spazio pari al totale dei punti (anche quelli ripetuti di chiusura che servono per i link)
    vPt.reserve( nTot) ;
    // inserisco i punti del contorno esterno (tranne il primo che coincide con l'ultimo)
-   if ( ! GetPntVectorFromPolyline( vPL[0], false, vPt))
-      return false ;
+   bOk = bOk && GetPntVectorFromPolyline( vPL[0], false, vPt) ;
    // ciclo sui percorsi interni ordinati secondo Xmax decrescente
    PNTVECTOR vPi ;
    vPi.reserve( nMax) ;
-   for ( int i = 1 ; i < int( vPL.size()) ; ++ i) {
+   for ( int i = 1 ; bOk && i < int( vPL.size()) ; ++ i) {
       // riordino il percorso per avere punto iniziale con Xmax
-      if ( ! GetPntVectorFromPolyline( vPL[vOrd[i]], true, vPi))
-         return false ;
+      if ( ! GetPntVectorFromPolyline( vPL[vOrd[i]], true, vPi)) {
+         bOk = false ;
+         break ;
+      }
       // cerco un punto del percorso esterno visibile dal punto iniziale del precedente interno
       int nI ;
-      if ( ! GetOuterPntToJoin( vPt, vPi[0], nI))
-         return false ;
+      if ( ! GetOuterPntToJoin( vPt, vPi[0], nI)) {
+         bOk = false ;
+         break ;
+      }
       // riordino percorso esterno per avere questo punto all'inizio
-      if ( ! ChangeStartPntVector( nI, vPt))
-         return false ;
+      if ( ! ChangeStartPntVector( nI, vPt)) {
+         bOk = false ;
+         break ;
+      }
       // ripeto punto iniziale
       vPt.push_back( vPt[0]) ;
       // accodo percorso interno
@@ -167,20 +183,14 @@ Triangulate::Make( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr)
       for ( int i = 0 ; i < int( vPL.size()) ; ++i)
          const_cast<PolyLine&>(vPL[i]).Invert( true) ;
    }
-
-  // creo il vettore degli indici del Poligono
-   INTVECTOR vPol ;
-   int n = int( vPt.size()) ;
-   vPol.reserve( n) ;
-   // non devo gestire separatamente CCW perchè ho già invertito i punti
-   for ( int i = 0 ; i < n ; ++ i)
-      vPol.push_back( i) ;
+   if ( ! bOk)
+      return MakeByEC_HPP( vPL, bCCW, vPt, vTr) ;
 
   // eseguo la triangolazione
-   if ( ! MakeByEC2( vPt, vPol, vTr)  &&
-        ! MakeByEC3( vPt, vPol, vTr)) {
+   if ( ! MakeByEC2( vPt, vTr)  &&
+        ! MakeByEC3( vPt, vTr)) {
       LOG_ERROR( GetEGkLogger(), "Error in MakeByEC23(N)")
-      return false ;
+      return MakeByEC_HPP( vPL, bCCW, vPt, vTr) ;
    }
 
   // se era CW, devo invertire il senso dei triangoli
@@ -190,33 +200,113 @@ Triangulate::Make( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr)
    return true ;
 }
 
+//----------------------------------------------------------------------------
+// Triangulate the CCW n-gon specified by the vertices vPt (Pt[n] != Pt[0])
+// Ear Clipping algorithm from mapbox
 //----------------------------------------------------------------------------
 bool
-Triangulate::PrepareGrid( const PNTVECTOR& vPt, const INTVECTOR& vPol,
-                          const INTVECTOR& vPrev, const INTVECTOR& vNext)
+Triangulate::MakeByEC_HPP( const PolyLine& PL, bool bCCW, PNTVECTOR& vPt, INTVECTOR& vTr)
 {
-  // points number
-   int n = int( vPol.size()) ;
-  // overall box
-   BBox3d b3All ;
-   for ( int j = 0 ; j < n ; ++ j)
-      b3All.Add( vPt[vPol[j]]) ;
-  // grid cell dimension
-   double dCellDim ;
-   b3All.GetRadius( dCellDim) ;
-   dCellDim *=  2 / sqrt( n) ;
-  // grid 
-   if ( ! m_VertGrid.Init( 2 * n, dCellDim))
+  // pulisco i vettori di ritorno
+   vPt.clear() ;
+   vTr.clear() ;
+  // riempio il vettore con i vertici del poligono da triangolare
+   vPt.reserve( PL.GetPointNbr() - 1) ;
+   // salto il primo punto (coincide con l'ultimo)
+   Point3d ptP ;
+   if ( ( bCCW  &&  ! PL.GetFirstPoint( ptP))  ||
+        ( ! bCCW  &&  ! PL.GetLastPoint( ptP)))
       return false ;
-   for ( int j = 0 ; j < n ; ++ j) {
-     // insert only reflex vertex
-      if ( ! TriangleIsCCW( vPt[vPol[vPrev[j]]], vPt[vPol[j]], vPt[vPol[vNext[j]]])) {
-         if ( ! m_VertGrid.InsertPoint( vPt[vPol[j]], j))
-            return false ;
+   // inserisco i punti
+   while ( ( bCCW  &&  PL.GetNextPoint( ptP))  ||
+            ( ! bCCW  &&   PL.GetPrevPoint( ptP)))
+      vPt.push_back( ptP) ;
+  // Create array
+   using Point = array<double, 2> ;
+   vector<vector<Point>> polygon ;
+  // Fill polygon structure with actual data. Any winding order works.
+   // The first polyline defines the main polygon.
+   polygon.resize( 1) ;
+   polygon[0].reserve( vPt.size()) ;
+   for ( int i = 0 ; i < int( vPt.size()) ;  ++i) {
+      switch ( m_nPlane) {
+      default : // PL_XY
+         polygon[0].push_back( { vPt[i].x, vPt[i].y}) ;
+         break ;
+      case PL_YZ :
+         polygon[0].push_back( { vPt[i].y, vPt[i].z}) ;
+         break ;
+      case PL_ZX :
+         polygon[0].push_back( { vPt[i].z, vPt[i].x}) ;
+         break ;
       }
    }
-   m_vVert.reserve( n / 5) ;
+  // Run tessellation
+   // Returns array of indices that refer to the vertices of the input polygon.
+   // Three subsequent indices form a triangle. Output triangles are counterclockwise.
+   vTr = mapbox::earcut<int32_t>( polygon) ;
+   if ( vTr.empty()) {
+      LOG_ERROR( GetEGkLogger(), "Error in MakeByEC_HPP(1)")
+      return false ;
+   }
+   if ( ! bCCW)
+      reverse( vTr.begin(), vTr.end()) ;
+   return true ;
+}
 
+//----------------------------------------------------------------------------
+bool
+Triangulate::MakeByEC_HPP( const POLYLINEVECTOR& vPL, bool bCCW, PNTVECTOR& vPt, INTVECTOR& vTr)
+{
+  // pulisco i vettori di ritorno
+   vPt.clear() ;
+   vTr.clear() ;
+  // Creo il vettore con i punti dei poligoni da triangolare
+   // calcolo spazio totale e spazio massimo per un percorso interno
+   int nTot = vPL[0].GetPointNbr() - 1 ;
+   for ( int i = 1 ; i < int( vPL.size()) ; ++i)
+      nTot += vPL[i].GetPointNbr() - 1 ;
+   // riservo spazio pari al totale dei punti (esclusi quelli ripetuti)
+   vPt.reserve( nTot) ;
+  // Create array
+   using Point = array<double, 2> ;
+   vector<vector<Point>> polygon ;
+   polygon.resize( vPL.size()) ;
+  // Eseguo riempimento di vettore e array
+   for ( int i = 0 ; i < int( vPL.size()) ; ++i) {
+      polygon[i].reserve( vPL[i].GetPointNbr() - 1) ;
+     // salto il primo punto (coincide con l'ultimo)
+      Point3d ptP ;
+      if ( ( bCCW  &&  ! vPL[i].GetFirstPoint( ptP))  ||
+            ( ! bCCW  &&  ! vPL[i].GetLastPoint( ptP)))
+         return false ;
+     // inserisco i punti
+      while ( ( bCCW  &&  vPL[i].GetNextPoint( ptP))  ||
+              ( ! bCCW  &&   vPL[i].GetPrevPoint( ptP))) {
+         vPt.push_back( ptP) ;
+         switch ( m_nPlane) {
+         default : // PL_XY
+            polygon[i].push_back( { ptP.x, ptP.y}) ;
+            break ;
+         case PL_YZ :
+            polygon[i].push_back( { ptP.y, ptP.z}) ;
+            break ;
+         case PL_ZX :
+            polygon[i].push_back( { ptP.z, ptP.x}) ;
+            break ;
+         }
+      }
+   }
+   // Run tessellation
+   // Returns array of indices that refer to the vertices of the input polygon.
+   // Three subsequent indices form a triangle. Output triangles are counterclockwise.
+   vTr = mapbox::earcut<int32_t>( polygon) ;
+   if ( vTr.empty()) {
+      LOG_ERROR( GetEGkLogger(), "Error in MakeByEC_HPP(N)")
+      return false ;
+   }
+   if ( ! bCCW)
+      reverse( vTr.begin(), vTr.end()) ;
    return true ;
 }
 
@@ -225,16 +315,22 @@ Triangulate::PrepareGrid( const PNTVECTOR& vPt, const INTVECTOR& vPol,
 // Ear Clipping algorithm
 //----------------------------------------------------------------------------
 bool
-Triangulate::MakeByEC( const PNTVECTOR& vPt, const INTVECTOR& vPol, INTVECTOR& vTr)
+Triangulate::MakeByEC( const PNTVECTOR& vPt, INTVECTOR& vTr)
 {
   // Clear triangle vector
    vTr.clear() ;
 
   // At least 3 points
-   int n = int( vPol.size()) ;
+   int n = int( vPt.size()) ;
    if ( n < 3)
       return false ;
 
+  // Creo il vettore degli indici del Poligono
+   INTVECTOR vPol ;
+   vPol.reserve( n) ;
+   for ( int i = 0 ; i < n ; ++ i)
+      vPol.push_back( i) ;
+
   // Preallocate triangle vector ( #triangles = n - 2)
    vTr.reserve( 3 * ( n - 2)) ;
 
@@ -297,16 +393,22 @@ Triangulate::MakeByEC( const PNTVECTOR& vPt, const INTVECTOR& vPol, INTVECTOR& v
 // Ear Clipping algorithm enhanced, choose smaller diagonal
 //----------------------------------------------------------------------------
 bool
-Triangulate::MakeByEC2( const PNTVECTOR& vPt, const INTVECTOR& vPol, INTVECTOR& vTr)
+Triangulate::MakeByEC2( const PNTVECTOR& vPt, INTVECTOR& vTr)
 {
   // Clear triangle vector
    vTr.clear() ;
 
   // At least 3 points
-   int n = int( vPol.size()) ;
+   int n = int( vPt.size()) ;
    if ( n < 3)
       return false ;
 
+  // Creo il vettore degli indici del Poligono
+   INTVECTOR vPol ;
+   vPol.reserve( n) ;
+   for ( int i = 0 ; i < n ; ++ i)
+      vPol.push_back( i) ;
+
   // Preallocate triangle vector ( #triangles = n - 2)
    vTr.reserve( 3 * ( n - 2)) ;
 
@@ -418,16 +520,22 @@ Triangulate::MakeByEC2( const PNTVECTOR& vPt, const INTVECTOR& vPol, INTVECTOR&
 // AR = ( Lmax * Lmax) / ( 2 * Area) = SqLenMax / L1 * L2
 //----------------------------------------------------------------------------
 bool
-Triangulate::MakeByEC3( const PNTVECTOR& vPt, const INTVECTOR& vPol, INTVECTOR& vTr)
+Triangulate::MakeByEC3( const PNTVECTOR& vPt, INTVECTOR& vTr)
 {
   // Clear triangle vector
    vTr.clear() ;
 
   // At least 3 points
-   int n = int( vPol.size()) ;
+   int n = int( vPt.size()) ;
    if ( n < 3)
       return false ;
 
+  // Creo il vettore degli indici del Poligono
+   INTVECTOR vPol ;
+   vPol.reserve( n) ;
+   for ( int i = 0 ; i < n ; ++ i)
+      vPol.push_back( i) ;
+
   // Preallocate triangle vector ( #triangles = n - 2)
    vTr.reserve( 3 * ( n - 2)) ;
 
@@ -525,6 +633,36 @@ Triangulate::MakeByEC3( const PNTVECTOR& vPt, const INTVECTOR& vPol, INTVECTOR&
    return true ;
 }
 
+//----------------------------------------------------------------------------
+bool
+Triangulate::PrepareGrid( const PNTVECTOR& vPt, const INTVECTOR& vPol,
+                          const INTVECTOR& vPrev, const INTVECTOR& vNext)
+{
+  // points number
+   int n = int( vPol.size()) ;
+  // overall box
+   BBox3d b3All ;
+   for ( int j = 0 ; j < n ; ++ j)
+      b3All.Add( vPt[vPol[j]]) ;
+  // grid cell dimension
+   double dCellDim ;
+   b3All.GetRadius( dCellDim) ;
+   dCellDim *=  2 / sqrt( n) ;
+  // grid 
+   if ( ! m_VertGrid.Init( 2 * n, dCellDim))
+      return false ;
+   for ( int j = 0 ; j < n ; ++ j) {
+     // insert only reflex vertex
+      if ( ! TriangleIsCCW( vPt[vPol[vPrev[j]]], vPt[vPol[j]], vPt[vPol[vNext[j]]])) {
+         if ( ! m_VertGrid.InsertPoint( vPt[vPol[j]], j))
+            return false ;
+      }
+   }
+   m_vVert.reserve( n / 5) ;
+
+   return true ;
+}
+
 //----------------------------------------------------------------------------
 bool
 Triangulate::TestTriangle( const PNTVECTOR& vPt, const INTVECTOR& vPol,
@@ -800,7 +938,7 @@ Triangulate::GetOuterPntToJoin( const PNTVECTOR& vPt, const Point3d& ptP, int& n
      // indice punto successivo
       int j = ( i + 1 < nNumPt) ? i + 1 : 0 ;
      // mi metto nel piano principale
-      switch (m_nPlane) {
+      switch ( m_nPlane) {
       default :  // PL_XY
         // se i punti coincidono
          if ( abs( vPt[i].x - ptP.x) < EPS_SMALL  &&  abs( vPt[i].y - ptP.y) < EPS_SMALL) {

Triangulate.h

@@ -24,11 +24,13 @@ class Triangulate
       bool Make( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr) ;
 
    private :
+      bool MakeByEC_HPP( const PolyLine& PL, bool bCCW, PNTVECTOR& vPt, INTVECTOR& vTr) ;
+      bool MakeByEC_HPP( const POLYLINEVECTOR& vPL, bool bCCW, PNTVECTOR& vPt, INTVECTOR& vTr) ;
+      bool MakeByEC( const PNTVECTOR& vPt, INTVECTOR& vTr) ;
+      bool MakeByEC2( const PNTVECTOR& vPt, INTVECTOR& vTr) ;
+      bool MakeByEC3( const PNTVECTOR& vPt, INTVECTOR& vTr) ;
       bool PrepareGrid( const PNTVECTOR& vPt, const INTVECTOR& vPol,
                         const INTVECTOR& vPrev, const INTVECTOR& vNext) ;
-      bool MakeByEC( const PNTVECTOR& vPt, const INTVECTOR& vPol, INTVECTOR& vTr) ;
-      bool MakeByEC2( const PNTVECTOR& vPt, const INTVECTOR& vPol, INTVECTOR& vTr) ;
-      bool MakeByEC3( const PNTVECTOR& vPt, const INTVECTOR& vPol, INTVECTOR& vTr) ;
       bool TestTriangle( const PNTVECTOR& vPt, const INTVECTOR& vPol,
                          const INTVECTOR& vPrev, INTVECTOR& vNext, int i) ;
       double CalcTriangleAspectRatio( const Point3d& ptPa, const Point3d& ptPb, const Point3d& ptPc) ;

VolZmap.cpp

@@ -35,7 +35,7 @@ GEOOBJ_REGISTER( VOL_ZMAP, NGE_V_ZMP, VolZmap) ;
 //----------------------------------------------------------------------------
 VolZmap::VolZmap(void)
    : m_nStatus( TO_VERIFY), m_dStep( 10.0), m_nMapNum( 0), m_nShape( GENERIC), m_nVoxNumPerBlock( N_VOXBLOCK),
-     m_nDexVoxRatio( 1), m_nNumBlock( 0), m_nConnectedCompoCount( 0), m_Tool( true)
+     m_nDexVoxRatio( 1), m_nNumBlock( 0), m_nConnectedCompoCount( 0), m_nCurrTool( - 1)
 {
    for ( int i = 0 ; i < N_MAPS ; ++ i) {
       m_nNx[i] = 0 ;
@@ -47,7 +47,6 @@ VolZmap::VolZmap(void)
    }
    m_nTempProp[0] = 0 ;
    m_nTempProp[1] = 0 ;
-   m_Tool.SetTolerances( LIN_TOL_STD, ANG_TOL_APPROX_DEG) ;
 }
 
 //----------------------------------------------------------------------------
@@ -76,7 +75,7 @@ VolZmap::Clear( void)
    m_dStep = EPS_SMALL ;
    m_nTempProp[0] = 0 ;
    m_nTempProp[1] = 0 ;
-   m_Tool.Clear() ;
+   ResetAllTools() ;
   // imposto ricalcolo della grafica
    m_OGrMgr.Reset() ;
 
@@ -1680,8 +1679,10 @@ VolZmap::Cut( const Plane3d& plPlane)
    Plane3d plMyPlane = plPlane ;
    plMyPlane.ToLoc( m_MapFrame) ;
   // Imposto numero fittizio di utensile per avere il colore di sezione opportuno
-   int nToolNumOld = m_Tool.GetCurrentToolNum() ;
-   m_Tool.SetToolNum( 1) ;
+   if ( m_nCurrTool < 0)
+      return false ;
+   int nToolNumOld = m_vTool[m_nCurrTool].GetCurrentToolNum() ;
+   m_vTool[m_nCurrTool].SetToolNum( 1) ;
   // Interseco lo Zmap col piano, ciclando sulle griglie
    bool bModified = false ;
    for ( int nMap = 0 ; nMap < int( m_nMapNum) ; ++ nMap) {
@@ -1735,7 +1736,7 @@ VolZmap::Cut( const Plane3d& plPlane)
    }
 
   // Ripristino numero utensile
-   m_Tool.SetToolNum( nToolNumOld) ;
+   m_vTool[m_nCurrTool].SetToolNum( nToolNumOld) ;
 
    if ( bModified == true) {
      // Imposto forma generica
@@ -1970,14 +1971,18 @@ VolZmap::CalcBlockNum( void)
 bool 
 VolZmap::SetToolTolerances( double dLinTol, double dAngTolDeg) 
 {
-   return m_Tool.SetTolerances( dLinTol, dAngTolDeg) ;
+   if ( m_nCurrTool < 0)
+      return false ;
+   return m_vTool[m_nCurrTool].SetTolerances( dLinTol, dAngTolDeg) ;
 }
 
 //----------------------------------------------------------------------------
 bool
 VolZmap::SetStdTool( const string& sToolName, double dH, double dR, double dCornR, double dCutterH, int nFlag) 
 {
-   return m_Tool.SetStdTool( sToolName, dH, dR, dCornR, dCutterH, nFlag) ;
+   if ( m_nCurrTool < 0)
+      return false ;
+   return m_vTool[m_nCurrTool].SetStdTool( sToolName, dH, dR, dCornR, dCutterH, nFlag) ;
 }
 
 //----------------------------------------------------------------------------
@@ -1985,7 +1990,9 @@ bool
 VolZmap::SetAdvTool( const string& sToolName, 
                      double dH, double dR, double dTipH, double dTipR, double dCornR, double dCutterH, int nFlag)
 {
-   return m_Tool.SetAdvTool( sToolName, dH, dR, dTipH, dTipR, dCornR, dCutterH, nFlag) ;
+   if ( m_nCurrTool < 0)
+      return false ;
+   return m_vTool[m_nCurrTool].SetAdvTool( sToolName, dH, dR, dTipH, dTipR, dCornR, dCutterH, nFlag) ;
 }
 
 //----------------------------------------------------------------------------
@@ -1993,33 +2000,101 @@ bool
 VolZmap::SetSawTool( const string& sToolName,
                      double dH, double dR, double dThick, double dStemR, double dCornR, int nFlag)
 {
-   return m_Tool.SetSawTool( sToolName, dH, dR, dThick, dStemR, dCornR, nFlag) ;
+   if ( m_nCurrTool < 0)
+      return false ;
+   return m_vTool[m_nCurrTool].SetSawTool( sToolName, dH, dR, dThick, dStemR, dCornR, nFlag) ;
 }
 
 //----------------------------------------------------------------------------
 bool 
 VolZmap::SetGenTool( const string& sToolName, const ICurveComposite* pToolOutline, int nFlag)
 {
-   return m_Tool.SetGenTool( sToolName, pToolOutline, nFlag) ; 
+   if ( m_nCurrTool < 0)
+      return false ;
+   return m_vTool[m_nCurrTool].SetGenTool( sToolName, pToolOutline, nFlag) ; 
 }
 
 //----------------------------------------------------------------------------
 bool 
 VolZmap::SetMortiserTool( const string& sToolName, double dH, double dW, double dTh, double dRc, int nFlag)
 {
-   return m_Tool.SetMortiserTool( sToolName, dH, dW, dTh, dRc, nFlag) ;
+   if ( m_nCurrTool < 0)
+      return false ;
+   return m_vTool[m_nCurrTool].SetMortiserTool( sToolName, dH, dW, dTh, dRc, nFlag) ;
 }
 
 //----------------------------------------------------------------------------
 bool 
 VolZmap::SetChiselTool( const string& sToolName, double dH, double dW, double dTh, int nFlag)
 {
-   return m_Tool.SetChiselTool( sToolName, dH, dW, dTh, nFlag) ; 
+   if ( m_nCurrTool < 0)
+      return false ;
+   return m_vTool[m_nCurrTool].SetChiselTool( sToolName, dH, dW, dTh, nFlag) ; 
+}
+
+//----------------------------------------------------------------------------
+bool
+VolZmap::SetAdditiveTool(const std::string& sToolName, double dH, double dR, double dRC, int nFlag)
+{
+   if ( m_nCurrTool < 0)
+      return false ;
+   return m_vTool[m_nCurrTool].SetAdditiveTool( sToolName, dH, dR, dRC, nFlag) ;
+}
+
+//----------------------------------------------------------------------------
+bool
+VolZmap::SetCurrTool( int nCurrTool)
+{
+   if ( nCurrTool < 0  ||  nCurrTool >= int( m_vTool.size()))
+      return false ;
+   m_nCurrTool = nCurrTool ;
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+int
+VolZmap::GetToolCount( void) const
+{
+   return int( m_vTool.size()) ;
+}
+
+//----------------------------------------------------------------------------
+bool
+VolZmap::AddTool( void)
+{
+   m_vTool.emplace_back( true) ;
+   m_vTool.back().SetTolerances( LIN_TOL_STD, ANG_TOL_APPROX_DEG) ;
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+bool
+VolZmap::ResetAllTools( void)
+{
+   bool bOk = true ;
+   for ( int n = 0 ; n < int( m_vTool.size()) ; ++ n) {
+      bOk = bOk  ||  m_vTool[n].Clear() ;
+   }
+   m_vTool.clear() ;
+   return bOk ;
 }
 
 //----------------------------------------------------------------------------
 bool 
 VolZmap::ResetTool( void)
 {
-   return m_Tool.Clear() ; 
+   if ( m_vTool.empty())
+      return true ;
+   bool bOk = m_vTool.back().Clear() ;
+   m_vTool.erase( m_vTool.begin() + int( m_vTool.size()) - 1) ; 
+   return bOk ; 
+}
+
+//----------------------------------------------------------------------------
+const ICurveComposite&
+VolZmap::GetToolOutline( bool bApprox) const
+{
+   if ( m_nCurrTool < 0)
+      return cvEmptyOutline ;
+   return ( bApprox ? m_vTool[m_nCurrTool].GetApproxOutline() : m_vTool[m_nCurrTool].GetOutline()) ;
 }

VolZmap.h

@@ -72,6 +72,7 @@ class VolZmap : public IVolZmap, public IGeoObjRW
       bool Create( const Point3d& ptO, double dDimX, double dDimY, double dDimZ, double dStep, bool bTriDex) override ;
       bool CreateFromFlatRegion( const ISurfFlatRegion& Surf, double dDimZ, double dStep, bool bTriDex) override ;
       bool CreateFromTriMesh( const ISurfTriMesh& Surf, double dStep, bool bTriDex) override ;
+      bool CreateEmptyMap( const Point3d& ptO, double dLengthX, double dLengthY, double dLengthZ, double dStep, bool bTriDex) override ;
       int  GetBlockCount( void) const override ;
       int  GetBlockUpdatingCounter( int nBlock) const override ;
       bool GetBlockTriangles( int nBlock, TRIA3DEXVECTOR& vTria) const override ;
@@ -90,9 +91,15 @@ class VolZmap : public IVolZmap, public IGeoObjRW
       bool SetGenTool( const std::string& sToolName, const ICurveComposite* pToolOutline, int nFlag) override ;
       bool SetMortiserTool( const std::string& sToolName, double dH, double dW, double dTh, double dRc, int nFlag) override ;
       bool SetChiselTool( const std::string& sToolName, double dH, double dW, double dTh, int nFlag) override ;
+      bool SetAdditiveTool( const std::string& sToolName, double dH, double dR, double dRC, int nFlag) override ;
+      bool SetCurrTool( int nCurrTool) override ;
+      int GetToolCount( void) const override ;
+      int GetCurrTool( void) const override
+             { return m_nCurrTool ; }
+      bool AddTool( void) override ;
+      bool ResetAllTools( void) override ;
       bool ResetTool( void) override ;
-      const ICurveComposite& GetToolOutline( bool bApprox = false) const override
-                  { return ( bApprox ? m_Tool.GetApproxOutline() : m_Tool.GetOutline()) ;}
+      const ICurveComposite& GetToolOutline( bool bApprox = false) const override ;
       bool MillingStep( const Point3d& ptPs, const Vector3d& vtDs, const Point3d& ptPe, const Vector3d& vtDe) override ;
       bool MillingStep( const Point3d& ptPs, const Vector3d& vtDs, const Vector3d& vtAs,
                         const Point3d& ptPe, const Vector3d& vtDe, const Vector3d& vtAe) override ;
@@ -303,10 +310,45 @@ class VolZmap : public IVolZmap, public IGeoObjRW
       inline bool TestParaBBox( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtD, const Vector3d& vtA,
                                 double dLenX, double dLenY, double dLenZ,
                                 int& nStI, int& nStJ, int& nEnI, int& nEnJ) ;
+
+    // Asportazioni superfici elementari vuote
+      bool SurfCircCrown_ZDrilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtAx, double dMaxRad, double dMinRad) ;
+      bool SurfCircCrown_Drilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtAx, double dMaxRad, double dMinRad) ;
+      bool SurfCircCrown_ZMilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtAx, double dMaxRad, double dMinRad) ;
+      bool SurfCircCrown_Milling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtAx, double dMaxRad, double dMinRad) ;
+      /*bool SurfCyl_ZDrilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir,
+                              double dHei, double dRad, bool bTapB, bool bTapT) ;
+      bool SurfCyl_Drilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir,
+                             double dHei, double dRad, bool bTapB, bool bTapT) ;*/
+      bool SurfCyl_ZMilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir,
+                             double dHei, double dRad, bool bOuterCutter) ;
+      bool SurfCyl_Milling( int nGrid, const Point3d& ptS, const Point3d& ptE,
+                            const Vector3d& vtToolDir, double dHei, double dRad, bool bOuterCutter, bool bTapB, bool bTapT) ;
+      bool SurfConus_ZDrilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, double dHei, double dMaxRad, double dMinRad,
+                                bool bOuterCutter, const Vector3d& vtArcNormMaxR, const Vector3d& vtArcNormMinR) ;
+      bool SurfConus_Drilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, 
+                               double dHei, double dMaxRad, double dMinRad, bool bOuterCutter, bool bTapB, bool bTapT,
+                               const Vector3d& vtArcNormMaxR, const Vector3d& vtArcNormMinR) ;
+      bool SurfConus_ZMilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, double dHei, double dMaxRad, double dMinRad,
+                               bool bOuterCutter, const Vector3d& vtArcNormMaxR, const Vector3d& vtArcNormMinR) ;
+      bool SurfConus_Milling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtAx, 
+                              double dHei, double dMaxRad, double dMinRad, bool bOuterCutter, bool bTapB, bool bTapT,
+                              const Vector3d& vtArcNormMaxR, const Vector3d& vtArcNormMinR) ;
+      bool SurfSphericalShellPart_Milling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtAx,
+                                           double dRad, double dInfH, double dSupH, bool bOuterCutter) ;
+      bool SurfSphericalShell_Milling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtAx,
+                                       double dRad, double dHei, bool bOuterCutter) ;
+    // Additivo
+      bool AddingMotion( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtAx) ;
+      bool AddingCylinder( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtAx, double dHei, double dRad) ;
+      bool AddingTruncatedCone( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtAx,
+                                double dMaxRad, double dMinRad, double dHei,
+                                const Vector3d& vtArcNormMaxR, const Vector3d& vtArcNormMinR) ;
+      bool AddingSphere( int nGrid, const Point3d& ptS, const Point3d& ptE, double dRad) ;
+      bool AddingGeneral( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtAx) ;
+
+
     // Intersezioni
-      bool IntersLineBox( const Point3d& ptP, const Vector3d& vtV, const Point3d& ptMin, const Point3d& ptMax) const ;
-      bool IntersLineBox( const Point3d& ptP, const Vector3d& vtV, const Point3d& ptMin, const Point3d& ptMax,
-                          double& dU1, double& dU2) const ;
       bool IntersLineZMapLattice( const Point3d& ptP, const Vector3d& vtV, double& dU1, double& dU2) const ;
       bool IntersLineZMapBBox( const Point3d& ptP, const Vector3d& vtV, double& dU1, double& dU2) const ;
       bool IntersLineDexel( const Point3d& ptP, const Vector3d& vtV, int nGrid, int nI, int nJ,
@@ -317,7 +359,7 @@ class VolZmap : public IVolZmap, public IGeoObjRW
       bool GetDepthWithVoxel( const Point3d& ptP, const Vector3d& vtDir, double& dInLength, double& dOutLength) const ;
       bool IntersLineCylinder( const Point3d& ptLineSt, const Vector3d& vtLineDir,
                                const Frame3d& CylFrame, double dH, double dRad, bool bTapLow, bool bTapUp,
-                               Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2) const ;
+                               Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2, bool bSubtracting = true) const ;
       bool IntersLineEllipticalCylinder( const Point3d& ptLineSt, const Vector3d& vtLineDir,
                                          const Frame3d& CircFrame, double dRad, double dLongMvLen, double dOrtMvLen,
                                          bool bTapLow, bool bTapUp,
@@ -327,11 +369,39 @@ class VolZmap : public IVolZmap, public IGeoObjRW
                             Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2) const ;
       bool IntersLineMyPolyhedron( const Point3d& ptLineSt, const Vector3d& vtLineDir,
                                    const Frame3d& PolyFrame, double dLenX, double dLenY, double dLenZ, double dDeltaX,
-                                   Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2) const ;
+                                   Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2, bool bSubtracting = true) const ;
       bool IntersLineTruncatedPyramid( const Point3d& ptLineSt, const Vector3d& vtLineDir,
                                        const Frame3d& frTruncPyramFrame, double dSegMin, double dSegMax, double dHeight,
                                        Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2) const ;
       bool TestIntersPlaneZmapBBox( const Plane3d& plPlane) const ;
+
+     // Intersezioni per asportazioni avanzate
+      int IntersLineCircCrown( const Point3d& ptLineP, const Vector3d& vtLineDir,
+                               const Point3d& ptCen, const Vector3d& vtAx, double dMaxRad, double dMinRad,
+                               Point3d& ptInt, Vector3d& vtN) const ;
+      int IntersLineParallelogram( const Point3d& ptLineP, const Vector3d& vtLineDir,
+                                   const Point3d& ptParOrig, const Vector3d& vtSeg1, const Vector3d& vtSeg2,
+                                   bool bExtNorm, Point3d& ptInt, Vector3d& vtN) const ;
+      int IntersLineCylinderCuttedByPlanes( const Point3d& ptLineP, const Vector3d& vtLineDir,
+                                            const Point3d& ptBaseCen, const Vector3d& vtAx, double dRad, double dH, bool bInOut,
+                                            const std::vector<Plane3d>& vPlanesVec,
+                                            Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2) const ;
+      int IntersLineCircSweptSurfCuttedByPlanes( const Point3d& ptLineP, const Vector3d& vtLineDir,
+                                                 const Point3d& ptCen, const Vector3d& vtAx, double dRad, const Vector3d& vtSweptVec, bool bInOut,
+                                                 const std::vector<Plane3d>& vPlanesVec,
+                                                 Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2) const ;
+      int IntersLineConeCuttedByPlanes( const Point3d& ptLineP, const Vector3d& vtLineDir,
+                                        const Point3d& ptVert, const Vector3d& vtAx, double dRad, double dH, bool bInOut,
+                                        const std::vector<Plane3d>& vPlanesVec,
+                                        Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2) const ;
+      int IntersLineSphereCuttedByPlanes( const Point3d& ptLineP, const Vector3d& vtLineD,
+                                          const Point3d& ptCen, double dRad, bool bInOut, 
+                                          const std::vector<Plane3d>& vPlanesVec,
+                                          Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2) const ;
+      int IntersLineCatTongue( const Point3d& ptLineP, const Vector3d& vtLineD,
+                               const Point3d& ptCenSt, const Point3d& ptCenEn, const Vector3d& vtNorm, double dRad,
+                               Point3d& ptInt, Vector3d& vtN) const ;
+     
      // Voxel: esistenza e passaggio da N a ijk per i voxel
       bool IsValidVoxel( int nN) const ;
       bool IsValidVoxel( int nI, int nJ, int nK) const ;
@@ -452,7 +522,10 @@ class VolZmap : public IVolZmap, public IGeoObjRW
       mutable std::atomic<bool>           m_bBreak ;
       std::atomic<bool>                   m_bIsBox ;
 
-      Tool m_Tool ;
+      // Utensili
+      std::vector<Tool> m_vTool ;
+      int m_nCurrTool ;
+      CurveComposite cvEmptyOutline ;
 } ;
 
 

VolZmapCalculus.cpp

@@ -16,10 +16,12 @@
 #include "CurveLine.h"
 #include "VolZmap.h"
 #include "GeoConst.h"
+#include "IntersLineBox.h"
 #include "IntersLineSurfStd.h"
 #include "/EgtDev/Include/EGkIntersLineTria.h"
 #include "/EgtDev/Include/EGkIntersLinePlane.h"
 #include "/EgtDev/Include/EGkIntersLineSphere.h"
+#include "/EgtDev/Include/EGkIntersPlaneBox.h"
 #include "/EgtDev/Include/EGkIntersPlaneTria.h"
 #include "/EgtDev/Include/EGkChainCurves.h"
 #include "/EgtDev/Include/ENkPolynomialRoots.h"
@@ -30,103 +32,6 @@
 
 using namespace std ;
 
-//----------------------------------------------------------------------------
-// Box e piano devono essere nello stesso riferimento. Il box deve essere axis aligned.
-static bool
-TestIntersPlaneBox( const BBox3d& b3Box, const Plane3d& plPlane) 
-{
-  // Calcolo le distanze con segno dei punti dal piano
-   Point3d ptE0 = b3Box.GetMin() ;
-   double dDist0 = DistPointPlane( ptE0, plPlane) ;
-   Point3d ptE1( b3Box.GetMax().x, b3Box.GetMin().y, b3Box.GetMin().z) ;
-   double dDist1 = DistPointPlane( ptE1, plPlane) ; 
-   Point3d ptE2( b3Box.GetMax().x, b3Box.GetMax().y, b3Box.GetMin().z) ;
-   double dDist2 = DistPointPlane( ptE2, plPlane) ; 
-   Point3d ptE3( b3Box.GetMin().x, b3Box.GetMax().y, b3Box.GetMin().z) ;
-   double dDist3 = DistPointPlane( ptE3, plPlane) ;  
-   Point3d ptE4( b3Box.GetMin().x, b3Box.GetMin().y, b3Box.GetMax().z) ;
-   double dDist4 = DistPointPlane( ptE4, plPlane) ; 
-   Point3d ptE5( b3Box.GetMax().x, b3Box.GetMin().y, b3Box.GetMax().z) ;
-   double dDist5 = DistPointPlane( ptE5, plPlane) ; 
-   Point3d ptE6 = b3Box.GetMax() ;
-   double dDist6 = DistPointPlane( ptE6, plPlane) ; 
-   Point3d ptE7( b3Box.GetMin().x,  b3Box.GetMax().y, b3Box.GetMax().z) ;
-   double dDist7 = DistPointPlane( ptE7, plPlane) ;
-  // Distanze tutte positive
-   if ( dDist0 > EPS_SMALL  &&  dDist1 > EPS_SMALL  &&  dDist2 > EPS_SMALL  &&  dDist3 > EPS_SMALL  &&
-        dDist4 > EPS_SMALL  &&  dDist5 > EPS_SMALL  &&  dDist6 > EPS_SMALL  &&  dDist7 > EPS_SMALL)
-      return false ;
-  // Distanze tutte negative
-   if ( dDist0 < - EPS_SMALL  &&  dDist1 < - EPS_SMALL  &&  dDist2 < - EPS_SMALL  &&  dDist3 < - EPS_SMALL  &&
-        dDist4 < - EPS_SMALL  &&  dDist5 < - EPS_SMALL  &&  dDist6 < - EPS_SMALL  &&  dDist7 < - EPS_SMALL)
-      return false ;
-  // Il piano interseca il box
-   return true ;
-}
-
-//----------------------------------------------------------------------------
-// Punti e vettore devono esse espressi nel medesimo sistema di riferimento.
-// Il box è allineato con gli assi di tale sistema di riferimento.
-// Viene restituito true in caso di intersezione, false altrimenti.
-bool
-VolZmap::IntersLineBox( const Point3d& ptP, const Vector3d& vtV, const Point3d& ptMin, const Point3d& ptMax) const 
-{
-   double dU1, dU2 ;
-   return IntersLineBox( ptP, vtV, ptMin, ptMax, dU1, dU2) ;
-}
-
-//----------------------------------------------------------------------------
-// Punti e vettore devono esse espressi nel medesimo sistema di riferimento.
-// Il box è allineato con gli assi di tale sistema di riferimento.
-// In dU1, dU2 vengono restituiti i parametri a cui si trovano le intersezioni.
-// Viene restituito true in caso di intersezione, false altrimenti.
-bool
-VolZmap::IntersLineBox( const Point3d& ptP, const Vector3d& vtV,
-                        const Point3d& ptMin, const Point3d& ptMax, double& dU1, double& dU2) const
-{
-  // Il box è allineato agli assi
-   dU1 = - INFINITO ;
-   dU2 = INFINITO ;
-
-  // confronto con piani YZ (perpendicolari ad asse X)
-   if ( vtV.x > EPS_ZERO) {
-      dU1 = max( dU1, ( ptMin.x - ptP.x) / vtV.x) ;
-      dU2 = min( dU2, ( ptMax.x - ptP.x) / vtV.x) ;
-   }
-   else if ( vtV.x < - EPS_ZERO) {
-      dU1 = max( dU1, ( ptMax.x - ptP.x) / vtV.x) ;
-      dU2 = min( dU2, ( ptMin.x - ptP.x) / vtV.x) ;
-   }
-   else if ( ptP.x < ptMin.x - EPS_SMALL  ||  ptP.x > ptMax.x + EPS_SMALL)
-      return false ;
-
-  // confronto con piani ZX (perpendicolari ad asse Y)
-   if ( vtV.y > EPS_ZERO) {
-      dU1 = max( dU1, ( ptMin.y - ptP.y) / vtV.y) ;
-      dU2 = min( dU2, ( ptMax.y - ptP.y) / vtV.y) ;
-   }
-   else if ( vtV.y < - EPS_ZERO) {
-      dU1 = max( dU1, ( ptMax.y - ptP.y) / vtV.y) ;
-      dU2 = min( dU2, ( ptMin.y - ptP.y) / vtV.y) ;
-   }
-   else if ( ptP.y < ptMin.y - EPS_SMALL  ||  ptP.y > ptMax.y + EPS_SMALL)
-      return false ;
-
-  // confronto con piani XY (perpendicolari ad asse Z)
-   if ( vtV.z > EPS_ZERO) {
-      dU1 = max( dU1, ( ptMin.z - ptP.z) / vtV.z) ;
-      dU2 = min( dU2, ( ptMax.z - ptP.z) / vtV.z) ;
-   }
-   else if ( vtV.z < - EPS_ZERO) {
-      dU1 = max( dU1, ( ptMax.z - ptP.z) / vtV.z) ;
-      dU2 = min( dU2, ( ptMin.z - ptP.z) / vtV.z) ;
-   }
-   else if ( ptP.z < ptMin.z - EPS_SMALL  ||  ptP.z > ptMax.z + EPS_SMALL)
-      return false ;
-
-   return ( dU2 >= dU1) ;
-}
-
 //----------------------------------------------------------------------------
 // La retta è rappresentata da un punto e dal versore espressi nel riferimento dello Zmap.
 // Se non vi è intersezione, viene restituito falso
@@ -2376,7 +2281,7 @@ VolZmap::AvoidSurfTm( const ISurfTriMesh& tmSurf, double dSafeDist, bool bPrecis
 bool
 VolZmap::IntersLineCylinder( const Point3d& ptLineSt, const Vector3d& vtLineDir,
                              const Frame3d& CylFrame, double dH, double dRad, bool bTapLow, bool bTapUp,
-                             Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2) const
+                             Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2, bool bSubtracting) const
 {  
   // Porto la linea nel riferimento del cilindro
    Point3d ptP = ptLineSt ;   ptP.ToLoc( CylFrame) ;
@@ -2386,14 +2291,14 @@ VolZmap::IntersLineCylinder( const Point3d& ptLineSt, const Vector3d& vtLineDir,
    int nBasInt = 0 ;
    if ( abs( vtV.z) > EPS_ZERO) {
       // le linee tangenti al cilindro non sono considerate intersecanti
-      double EpsRad = ( vtV.IsZeroXY() ? - EPS_SMALL : EPS_SMALL) ;
+      double dEpsRad = ( bSubtracting ? ( vtV.IsZeroXY() ? - EPS_SMALL : EPS_SMALL) : EPS_SMALL) ;
       ptInt1 = ptP + ( ( 0 - ptP.z) / vtV.z) * vtV ;
-      if ( ptInt1.x * ptInt1.x + ptInt1.y * ptInt1.y < dRad * dRad + 2 * dRad * EpsRad) {
+      if ( ptInt1.x * ptInt1.x + ptInt1.y * ptInt1.y < dRad * dRad + 2 * dRad * dEpsRad) {
          nBasInt += 1 ;
          vtN1 = Z_AX ;
       }
       ptInt2 = ptP + ( ( dH - ptP.z) / vtV.z) * vtV ;  
-      if ( ptInt2.x * ptInt2.x + ptInt2.y * ptInt2.y < dRad * dRad + 2 * dRad * EpsRad) {
+      if ( ptInt2.x * ptInt2.x + ptInt2.y * ptInt2.y < dRad * dRad + 2 * dRad * dEpsRad) {
          nBasInt += 2 ;
          vtN2 = - Z_AX ;
       }
@@ -2776,7 +2681,7 @@ VolZmap::IntersLineEllipticalCylinder( const Point3d& ptLineSt, const Vector3d&
 bool 
 VolZmap::IntersLineMyPolyhedron( const Point3d& ptLineSt, const Vector3d& vtLineDir,
                                  const Frame3d& PolyFrame, double dLenX, double dLenY, double dLenZ, double dDeltaZ,
-                                 Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2) const
+                                 Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2, bool bSubtracting) const
 {
   // Controllo sulle dimensioni lineari affinché sia valido il poliedro
    if ( dLenX <= EPS_SMALL  ||  dLenY <= EPS_SMALL  ||  dLenZ <= EPS_SMALL)
@@ -2800,17 +2705,19 @@ VolZmap::IntersLineMyPolyhedron( const Point3d& ptLineSt, const Vector3d& vtLine
 
   // Controlli affinché non vengano tagliati dexel a filo
   // con il passaggio dell'utensile:
-  // Controllo sulle facce 1 e 2
-   if ( abs( vtV.y) < EPS_ZERO  &&  abs( ptP.y) > ptFacet135.y)   
-      return false ;
-  // Controllo sulle facce 3 e 4
-   if ( abs( vtV.x) < EPS_ZERO  &&  ( ptP.x < ptFacet135.x  ||  ptP.x > ptFacet246.x))
-      return false ;
-  // Controllo sulle facce 5 e 6 
-   double dP1 = abs ( ( ptFacet135 - ptP) * vtOb) ;
-   double dP2 = abs ( ( ptFacet246 - ptP) * vtOb) ;
-   if ( abs( vtV * vtOb) < EPS_ZERO  &&  ( dP1 < EPS_SMALL  ||  dP2  < EPS_SMALL))
-      return false ;
+  if ( bSubtracting) {
+     // Controllo sulle facce 1 e 2
+      if ( abs( vtV.y) < EPS_ZERO  &&  abs( ptP.y) > ptFacet135.y)   
+         return false ;
+     // Controllo sulle facce 3 e 4
+      if ( abs( vtV.x) < EPS_ZERO  &&  ( ptP.x < ptFacet135.x  ||  ptP.x > ptFacet246.x))
+         return false ;
+     // Controllo sulle facce 5 e 6 
+      double dP1 = abs ( ( ptFacet135 - ptP) * vtOb) ;
+      double dP2 = abs ( ( ptFacet246 - ptP) * vtOb) ;
+      if ( abs( vtV * vtOb) < EPS_ZERO  &&  ( dP1 < EPS_SMALL  ||  dP2  < EPS_SMALL))
+         return false ;
+   }
 
   // Punti notevoli
    Point3d ptI1 = ptP + ( ( ptFacet135.y - ptP.y) / vtV.y) * vtV ;
@@ -3188,7 +3095,7 @@ VolZmap::GetLineIntersection( const Point3d& ptP, const Vector3d& vtD, ILZIVECTO
       if ( ! GetBlockBox( nBlockIJK, b3BlockBox))
          return false ;
      // Se c'è intersezione valuto tutti i voxel interni
-      if ( IntersLineBox( ptLocP, vtDirL, b3BlockBox.GetMin(), b3BlockBox.GetMax())) {
+      if ( TestIntersLineBox( ptLocP, vtDirL, b3BlockBox.GetMin(), b3BlockBox.GetMax())) {
          // Ciclo sui voxel del blocco.
          // Triangoli smooth
          for ( int nV = 0 ; nV < int( m_BlockSmoothTria[nB].size()) ; ++ nV) {
@@ -3199,7 +3106,7 @@ VolZmap::GetLineIntersection( const Point3d& ptP, const Vector3d& vtD, ILZIVECTO
                                   m_BlockSmoothTria[nB][nV].k } ;
             GetVoxelBox( nCurVoxIJK[0], nCurVoxIJK[1], nCurVoxIJK[2], b3Vox) ;
            // Se non c'è intersezione col voxel, passo al successivo.
-            if ( ! IntersLineBox( ptLocP, vtDirL, b3Vox.GetMin(), b3Vox.GetMax()))
+            if ( ! TestIntersLineBox( ptLocP, vtDirL, b3Vox.GetMin(), b3Vox.GetMax()))
                continue ;
             for ( int nT = 0 ; nT < int( m_BlockSmoothTria[nB][nV].vTria.size()) ; ++ nT) {
                Triangle3d trTria = m_BlockSmoothTria[nB][nV].vTria[nT] ;
@@ -3374,7 +3281,7 @@ VolZmap::GetPlaneIntersection( const Plane3d& plPlane, ICURVEPOVECTOR& vpLoop) c
       if ( ! GetBlockBox( nBlockIJK, b3BlockBox))
          return false ;
      // Se c'è intersezione valuto tutti i voxel interni
-      if ( TestIntersPlaneBox( b3BlockBox, plPlaneLoc)) {    
+      if ( TestIntersPlaneBox( plPlaneLoc, b3BlockBox)) {    
         // Ciclo sui voxel del blocco.    
         // Triangoli smooth
          for ( int nV = 0 ; nV < int( m_BlockSmoothTria[nB].size()) ; ++ nV) {
@@ -3382,7 +3289,7 @@ VolZmap::GetPlaneIntersection( const Plane3d& plPlane, ICURVEPOVECTOR& vpLoop) c
             BBox3d b3Vox ;
             GetVoxelBox( m_BlockSmoothTria[nB][nV].i, m_BlockSmoothTria[nB][nV].j, m_BlockSmoothTria[nB][nV].k, b3Vox) ;
             // Se non c'è intersezione col voxel, passo al successivo.
-            if ( ! TestIntersPlaneBox(b3Vox, plPlaneLoc))
+            if ( ! TestIntersPlaneBox( plPlaneLoc, b3Vox))
                continue ; 
             for ( int nT = 0 ; nT < int( m_BlockSmoothTria[nB][nV].vTria.size()) ; ++ nT) {
                Triangle3d trTria = m_BlockSmoothTria[nB][nV].vTria[nT] ;
@@ -3486,5 +3393,5 @@ bool
 VolZmap::TestIntersPlaneZmapBBox( const Plane3d& plPlane) const 
 {
    BBox3d b3Zmap( ORIG, Point3d( m_nNx[0] * m_dStep, m_nNy[0] * m_dStep, m_dMaxZ[0])) ;
-   return TestIntersPlaneBox( b3Zmap, plPlane) ;
+   return TestIntersPlaneBox( plPlane, b3Zmap) ;
 }

VolZmapCreation.cpp

@@ -665,3 +665,77 @@ VolZmap::CreateFromTriMesh( const ISurfTriMesh& Surf, double dStep, bool bTriDex
 
    return bCompleted ;
 }
+
+//----------------------------------------------------------------------------
+bool 
+VolZmap::CreateEmptyMap( const Point3d& ptO, double dLengthX, double dLengthY, double dLengthZ, double dStep, bool bTriDex)
+{
+  // Controlli sull'ammissibilità delle dimensioni lineari del grezzo e del passo
+   if ( dStep < EPS_SMALL  ||  dLengthX < EPS_SMALL  ||  dLengthY < EPS_SMALL  ||  dLengthZ < EPS_SMALL)
+      return false ;
+
+  // Il passo di discretizzazione non può essere inferiore a 100 * EPS_SMALL
+   m_dStep = max( dStep, 100 * EPS_SMALL) ;
+
+  // Aggiorno la dimensione della mappa 1 o 3
+   m_nMapNum = ( bTriDex ? 3 : 1) ;
+
+  // Disponendo i sistemi di riferimento in una successione, le coordinate x,y,z
+  // di uno si ottengono da una permutazione ciclica di quelle del precedente sistema.
+  // es: X(n) = Z(n-1), Y(n) = X(n-1), Z(n) = Y(n-1)
+
+  // Definisco il sistema di riferimento intrinseco
+   m_MapFrame.Set( ptO, X_AX, Y_AX, Z_AX) ;
+
+  // Definisco i vettori dei limiti su indici
+   m_nNx[0] = max( int( ( dLengthX + EPS_SMALL) / m_dStep + 0.5), 1) ;
+   m_nNy[0] = max( int( ( dLengthY + EPS_SMALL) / m_dStep + 0.5), 1) ;
+
+  // Numero di componenti connesse 
+   m_nConnectedCompoCount = 1 ;
+
+  // Se tridexel
+   if ( bTriDex) {
+      m_nNx[1] = m_nNy[0] ;
+      m_nNy[1] = max( int( ( dLengthZ + EPS_SMALL) / m_dStep + 0.5), 1) ;
+      m_nNx[2] = m_nNy[1] ;
+      m_nNy[2] = m_nNx[0] ;
+   }
+
+  // altrimenti mono dexel
+   else {
+      m_nNx[1] = 0 ;
+      m_nNy[1] = 0 ;
+      m_nNx[2] = 0 ;
+      m_nNy[2] = 0 ; 
+   }
+
+   // Definisco il numero di blocchi lungo x,y e z
+   if ( ! CalcBlockNum())
+      return false ;
+
+  // Creazione delle mappe
+  // Calcolo del numero di celle per ogni mappa
+   for ( int i = 0 ; i < m_nMapNum ; ++ i)
+      m_nDim[i] = m_nNx[i] * m_nNy[i] ;
+
+  // Creazione delle celle per ogni mappa
+   for ( int i = 0 ; i < m_nMapNum ; ++ i)
+      m_Values[i].resize( m_nDim[i]) ;
+
+  // Definizione delle limitazioni iniziali in Z per ogni mappa
+   m_dMinZ[0] = 0 ;
+   m_dMaxZ[0] = dLengthZ ;
+   m_dMinZ[1] = 0 ;
+   m_dMaxZ[1] = ( bTriDex ? dLengthX : 0) ;
+   m_dMinZ[2] = 0 ;  
+   m_dMaxZ[2] = ( bTriDex ? dLengthY : 0) ;
+
+   // Tipologia
+   m_nShape = GENERIC ;
+
+  // Aggiornamento dello stato
+   m_nStatus = OK ;
+
+   return true ;
+}

VolZmapGraphics.cpp

@@ -18,6 +18,7 @@
 #include "GeoConst.h"
 #include "MC_Tables.h"
 #include "PolygonPlane.h"
+#include "IntersLineBox.h"
 #include "/EgtDev/Include/EGkIntervals.h"
 #include "/EgtDev/Include/EGkStringUtils3d.h"
 #include "/EgtDev/Include/EGkChainCurves.h"
@@ -844,7 +845,7 @@ VolZmap::UpdateTripleMapGraphics( void) const
    VecTriHold.resize( m_nNumBlock) ;
 
   // Ciclo sui blocchi per eliminare le slice fra blocchi da aggiornare
-   for ( int t = 0 ; t < m_nNumBlock ; ++ t) {
+   for ( int t = 0 ; t < m_nNumBlock ; ++ t) { 
       for ( auto it = m_SliceXY[t].begin() ; it != m_SliceXY[t].end() ;) {
          int nSlIJK[3] ;
          if ( GetVoxIJKFromN( it->first, nSlIJK[0], nSlIJK[1], nSlIJK[2])) {
@@ -1083,7 +1084,7 @@ VolZmap::ExtMarchingCubes( int nBlock, VoxelContainer& vVox) const
 
            // Flag di regolarità dei campi scalare e vettoriale
             bool bReg = true ;
-
+            
            // Ciclo sui segmenti
             for ( int EdgeIndex = 0 ; EdgeIndex < 12 ; ++ EdgeIndex) {
               // Se il segmento non attraversa la superficie passo al successivo

earcut.hpp

@@ -0,0 +1,823 @@
+#pragma once
+
+#include <algorithm>
+#include <cassert>
+#include <cmath>
+#include <cstddef>
+#include <limits>
+#include <memory>
+#include <utility>
+#include <vector>
+
+namespace mapbox {
+
+namespace util {
+
+template <std::size_t I, typename T> struct nth {
+    inline static typename std::tuple_element<I, T>::type
+    get(const T& t) { return std::get<I>(t); };
+};
+
+}
+
+namespace detail {
+
+template <typename N = uint32_t>
+class Earcut {
+public:
+    std::vector<N> indices;
+    std::size_t vertices = 0;
+
+    template <typename Polygon>
+    void operator()(const Polygon& points);
+
+private:
+    struct Node {
+        Node(N index, double x_, double y_) : i(index), x(x_), y(y_) {}
+        Node(const Node&) = delete;
+        Node& operator=(const Node&) = delete;
+        Node(Node&&) = delete;
+        Node& operator=(Node&&) = delete;
+
+        const N i;
+        const double x;
+        const double y;
+
+        // previous and next vertice nodes in a polygon ring
+        Node* prev = nullptr;
+        Node* next = nullptr;
+
+        // z-order curve value
+        int32_t z = 0;
+
+        // previous and next nodes in z-order
+        Node* prevZ = nullptr;
+        Node* nextZ = nullptr;
+
+        // indicates whether this is a steiner point
+        bool steiner = false;
+    };
+
+    template <typename Ring> Node* linkedList(const Ring& points, const bool clockwise);
+    Node* filterPoints(Node* start, Node* end = nullptr);
+    void earcutLinked(Node* ear, int pass = 0);
+    bool isEar(Node* ear);
+    bool isEarHashed(Node* ear);
+    Node* cureLocalIntersections(Node* start);
+    void splitEarcut(Node* start);
+    template <typename Polygon> Node* eliminateHoles(const Polygon& points, Node* outerNode);
+    Node* eliminateHole(Node* hole, Node* outerNode);
+    Node* findHoleBridge(Node* hole, Node* outerNode);
+    bool sectorContainsSector(const Node* m, const Node* p);
+    void indexCurve(Node* start);
+    Node* sortLinked(Node* list);
+    int32_t zOrder(const double x_, const double y_);
+    Node* getLeftmost(Node* start);
+    bool pointInTriangle(double ax, double ay, double bx, double by, double cx, double cy, double px, double py) const;
+    bool isValidDiagonal(Node* a, Node* b);
+    double area(const Node* p, const Node* q, const Node* r) const;
+    bool equals(const Node* p1, const Node* p2);
+    bool intersects(const Node* p1, const Node* q1, const Node* p2, const Node* q2);
+    bool onSegment(const Node* p, const Node* q, const Node* r);
+    int sign(double val);
+    bool intersectsPolygon(const Node* a, const Node* b);
+    bool locallyInside(const Node* a, const Node* b);
+    bool middleInside(const Node* a, const Node* b);
+    Node* splitPolygon(Node* a, Node* b);
+    template <typename Point> Node* insertNode(std::size_t i, const Point& p, Node* last);
+    void removeNode(Node* p);
+
+    bool hashing;
+    double minX, maxX;
+    double minY, maxY;
+    double inv_size = 0;
+
+    template <typename T, typename Alloc = std::allocator<T>>
+    class ObjectPool {
+    public:
+        ObjectPool() { }
+        ObjectPool(std::size_t blockSize_) {
+            reset(blockSize_);
+        }
+        ~ObjectPool() {
+            clear();
+        }
+        template <typename... Args>
+        T* construct(Args&&... args) {
+            if (currentIndex >= blockSize) {
+                currentBlock = alloc_traits::allocate(alloc, blockSize);
+                allocations.emplace_back(currentBlock);
+                currentIndex = 0;
+            }
+            T* object = &currentBlock[currentIndex++];
+            alloc_traits::construct(alloc, object, std::forward<Args>(args)...);
+            return object;
+        }
+        void reset(std::size_t newBlockSize) {
+            for (auto allocation : allocations) {
+                alloc_traits::deallocate(alloc, allocation, blockSize);
+            }
+            allocations.clear();
+            blockSize = std::max<std::size_t>(1, newBlockSize);
+            currentBlock = nullptr;
+            currentIndex = blockSize;
+        }
+        void clear() { reset(blockSize); }
+    private:
+        T* currentBlock = nullptr;
+        std::size_t currentIndex = 1;
+        std::size_t blockSize = 1;
+        std::vector<T*> allocations;
+        Alloc alloc;
+        typedef typename std::allocator_traits<Alloc> alloc_traits;
+    };
+    ObjectPool<Node> nodes;
+};
+
+template <typename N> template <typename Polygon>
+void Earcut<N>::operator()(const Polygon& points) {
+    // reset
+    indices.clear();
+    vertices = 0;
+
+    if (points.empty()) return;
+
+    double x;
+    double y;
+    int threshold = 80;
+    std::size_t len = 0;
+
+    for (size_t i = 0; threshold >= 0 && i < points.size(); i++) {
+        threshold -= static_cast<int>(points[i].size());
+        len += points[i].size();
+    }
+
+    //estimate size of nodes and indices
+    nodes.reset(len * 3 / 2);
+    indices.reserve(len + points[0].size());
+
+    Node* outerNode = linkedList(points[0], true);
+    if (!outerNode || outerNode->prev == outerNode->next) return;
+
+    if (points.size() > 1) outerNode = eliminateHoles(points, outerNode);
+
+    // if the shape is not too simple, we'll use z-order curve hash later; calculate polygon bbox
+    hashing = threshold < 0;
+    if (hashing) {
+        Node* p = outerNode->next;
+        minX = maxX = outerNode->x;
+        minY = maxY = outerNode->y;
+        do {
+            x = p->x;
+            y = p->y;
+            minX = std::min<double>(minX, x);
+            minY = std::min<double>(minY, y);
+            maxX = std::max<double>(maxX, x);
+            maxY = std::max<double>(maxY, y);
+            p = p->next;
+        } while (p != outerNode);
+
+        // minX, minY and size are later used to transform coords into integers for z-order calculation
+        inv_size = std::max<double>(maxX - minX, maxY - minY);
+        inv_size = inv_size != .0 ? (1. / inv_size) : .0;
+    }
+
+    earcutLinked(outerNode);
+
+    nodes.clear();
+}
+
+// create a circular doubly linked list from polygon points in the specified winding order
+template <typename N> template <typename Ring>
+typename Earcut<N>::Node*
+Earcut<N>::linkedList(const Ring& points, const bool clockwise) {
+    using Point = typename Ring::value_type;
+    double sum = 0;
+    const std::size_t len = points.size();
+    std::size_t i, j;
+    Node* last = nullptr;
+
+    // calculate original winding order of a polygon ring
+    for (i = 0, j = len > 0 ? len - 1 : 0; i < len; j = i++) {
+        const auto& p1 = points[i];
+        const auto& p2 = points[j];
+        const double p20 = util::nth<0, Point>::get(p2);
+        const double p10 = util::nth<0, Point>::get(p1);
+        const double p11 = util::nth<1, Point>::get(p1);
+        const double p21 = util::nth<1, Point>::get(p2);
+        sum += (p20 - p10) * (p11 + p21);
+    }
+
+    // link points into circular doubly-linked list in the specified winding order
+    if (clockwise == (sum > 0)) {
+        for (i = 0; i < len; i++) last = insertNode(vertices + i, points[i], last);
+    } else {
+        for (i = len; i-- > 0;) last = insertNode(vertices + i, points[i], last);
+    }
+
+    if (last && equals(last, last->next)) {
+        removeNode(last);
+        last = last->next;
+    }
+
+    vertices += len;
+
+    return last;
+}
+
+// eliminate colinear or duplicate points
+template <typename N>
+typename Earcut<N>::Node*
+Earcut<N>::filterPoints(Node* start, Node* end) {
+    if (!end) end = start;
+
+    Node* p = start;
+    bool again;
+    do {
+        again = false;
+
+        if (!p->steiner && (equals(p, p->next) || area(p->prev, p, p->next) == 0)) {
+            removeNode(p);
+            p = end = p->prev;
+
+            if (p == p->next) break;
+            again = true;
+
+        } else {
+            p = p->next;
+        }
+    } while (again || p != end);
+
+    return end;
+}
+
+// main ear slicing loop which triangulates a polygon (given as a linked list)
+template <typename N>
+void Earcut<N>::earcutLinked(Node* ear, int pass) {
+    if (!ear) return;
+
+    // interlink polygon nodes in z-order
+    if (!pass && hashing) indexCurve(ear);
+
+    Node* stop = ear;
+    Node* prev;
+    Node* next;
+
+    int iterations = 0;
+
+    // iterate through ears, slicing them one by one
+    while (ear->prev != ear->next) {
+        iterations++;
+        prev = ear->prev;
+        next = ear->next;
+
+        if (hashing ? isEarHashed(ear) : isEar(ear)) {
+            // cut off the triangle
+            indices.emplace_back(prev->i);
+            indices.emplace_back(ear->i);
+            indices.emplace_back(next->i);
+
+            removeNode(ear);
+
+            // skipping the next vertice leads to less sliver triangles
+            ear = next->next;
+            stop = next->next;
+
+            continue;
+        }
+
+        ear = next;
+
+        // if we looped through the whole remaining polygon and can't find any more ears
+        if (ear == stop) {
+            // try filtering points and slicing again
+            if (!pass) earcutLinked(filterPoints(ear), 1);
+
+            // if this didn't work, try curing all small self-intersections locally
+            else if (pass == 1) {
+                ear = cureLocalIntersections(filterPoints(ear));
+                earcutLinked(ear, 2);
+
+            // as a last resort, try splitting the remaining polygon into two
+            } else if (pass == 2) splitEarcut(ear);
+
+            break;
+        }
+    }
+}
+
+// check whether a polygon node forms a valid ear with adjacent nodes
+template <typename N>
+bool Earcut<N>::isEar(Node* ear) {
+    const Node* a = ear->prev;
+    const Node* b = ear;
+    const Node* c = ear->next;
+
+    if (area(a, b, c) >= 0) return false; // reflex, can't be an ear
+
+    // now make sure we don't have other points inside the potential ear
+    Node* p = ear->next->next;
+
+    while (p != ear->prev) {
+        if (pointInTriangle(a->x, a->y, b->x, b->y, c->x, c->y, p->x, p->y) &&
+            area(p->prev, p, p->next) >= 0) return false;
+        p = p->next;
+    }
+
+    return true;
+}
+
+template <typename N>
+bool Earcut<N>::isEarHashed(Node* ear) {
+    const Node* a = ear->prev;
+    const Node* b = ear;
+    const Node* c = ear->next;
+
+    if (area(a, b, c) >= 0) return false; // reflex, can't be an ear
+
+    // triangle bbox; min & max are calculated like this for speed
+    const double minTX = std::min<double>(a->x, std::min<double>(b->x, c->x));
+    const double minTY = std::min<double>(a->y, std::min<double>(b->y, c->y));
+    const double maxTX = std::max<double>(a->x, std::max<double>(b->x, c->x));
+    const double maxTY = std::max<double>(a->y, std::max<double>(b->y, c->y));
+
+    // z-order range for the current triangle bbox;
+    const int32_t minZ = zOrder(minTX, minTY);
+    const int32_t maxZ = zOrder(maxTX, maxTY);
+
+    // first look for points inside the triangle in increasing z-order
+    Node* p = ear->nextZ;
+
+    while (p && p->z <= maxZ) {
+        if (p != ear->prev && p != ear->next &&
+            pointInTriangle(a->x, a->y, b->x, b->y, c->x, c->y, p->x, p->y) &&
+            area(p->prev, p, p->next) >= 0) return false;
+        p = p->nextZ;
+    }
+
+    // then look for points in decreasing z-order
+    p = ear->prevZ;
+
+    while (p && p->z >= minZ) {
+        if (p != ear->prev && p != ear->next &&
+            pointInTriangle(a->x, a->y, b->x, b->y, c->x, c->y, p->x, p->y) &&
+            area(p->prev, p, p->next) >= 0) return false;
+        p = p->prevZ;
+    }
+
+    return true;
+}
+
+// go through all polygon nodes and cure small local self-intersections
+template <typename N>
+typename Earcut<N>::Node*
+Earcut<N>::cureLocalIntersections(Node* start) {
+    Node* p = start;
+    do {
+        Node* a = p->prev;
+        Node* b = p->next->next;
+
+        // a self-intersection where edge (v[i-1],v[i]) intersects (v[i+1],v[i+2])
+        if (!equals(a, b) && intersects(a, p, p->next, b) && locallyInside(a, b) && locallyInside(b, a)) {
+            indices.emplace_back(a->i);
+            indices.emplace_back(p->i);
+            indices.emplace_back(b->i);
+
+            // remove two nodes involved
+            removeNode(p);
+            removeNode(p->next);
+
+            p = start = b;
+        }
+        p = p->next;
+    } while (p != start);
+
+    return filterPoints(p);
+}
+
+// try splitting polygon into two and triangulate them independently
+template <typename N>
+void Earcut<N>::splitEarcut(Node* start) {
+    // look for a valid diagonal that divides the polygon into two
+    Node* a = start;
+    do {
+        Node* b = a->next->next;
+        while (b != a->prev) {
+            if (a->i != b->i && isValidDiagonal(a, b)) {
+                // split the polygon in two by the diagonal
+                Node* c = splitPolygon(a, b);
+
+                // filter colinear points around the cuts
+                a = filterPoints(a, a->next);
+                c = filterPoints(c, c->next);
+
+                // run earcut on each half
+                earcutLinked(a);
+                earcutLinked(c);
+                return;
+            }
+            b = b->next;
+        }
+        a = a->next;
+    } while (a != start);
+}
+
+// link every hole into the outer loop, producing a single-ring polygon without holes
+template <typename N> template <typename Polygon>
+typename Earcut<N>::Node*
+Earcut<N>::eliminateHoles(const Polygon& points, Node* outerNode) {
+    const size_t len = points.size();
+
+    std::vector<Node*> queue;
+    for (size_t i = 1; i < len; i++) {
+        Node* list = linkedList(points[i], false);
+        if (list) {
+            if (list == list->next) list->steiner = true;
+            queue.push_back(getLeftmost(list));
+        }
+    }
+    std::sort(queue.begin(), queue.end(), [](const Node* a, const Node* b) {
+        return a->x < b->x;
+    });
+
+    // process holes from left to right
+    for (size_t i = 0; i < queue.size(); i++) {
+        outerNode = eliminateHole(queue[i], outerNode);
+        outerNode = filterPoints(outerNode, outerNode->next);
+    }
+
+    return outerNode;
+}
+
+// find a bridge between vertices that connects hole with an outer ring and and link it
+template <typename N>
+typename Earcut<N>::Node*
+Earcut<N>::eliminateHole(Node* hole, Node* outerNode) {
+    Node* bridge = findHoleBridge(hole, outerNode);
+    if (!bridge) {
+        return outerNode;
+    }
+
+    Node* bridgeReverse = splitPolygon(bridge, hole);
+
+    // filter collinear points around the cuts
+    Node* filteredBridge = filterPoints(bridge, bridge->next);
+    filterPoints(bridgeReverse, bridgeReverse->next);
+
+    // Check if input node was removed by the filtering
+    return outerNode == bridge ? filteredBridge : outerNode;
+}
+
+// David Eberly's algorithm for finding a bridge between hole and outer polygon
+template <typename N>
+typename Earcut<N>::Node*
+Earcut<N>::findHoleBridge(Node* hole, Node* outerNode) {
+    Node* p = outerNode;
+    double hx = hole->x;
+    double hy = hole->y;
+    double qx = -std::numeric_limits<double>::infinity();
+    Node* m = nullptr;
+
+    // find a segment intersected by a ray from the hole's leftmost Vertex to the left;
+    // segment's endpoint with lesser x will be potential connection Vertex
+    do {
+        if (hy <= p->y && hy >= p->next->y && p->next->y != p->y) {
+          double x = p->x + (hy - p->y) * (p->next->x - p->x) / (p->next->y - p->y);
+          if (x <= hx && x > qx) {
+            qx = x;
+            if (x == hx) {
+                if (hy == p->y) return p;
+                if (hy == p->next->y) return p->next;
+            }
+            m = p->x < p->next->x ? p : p->next;
+          }
+        }
+        p = p->next;
+    } while (p != outerNode);
+
+    if (!m) return 0;
+
+    if (hx == qx) return m; // hole touches outer segment; pick leftmost endpoint
+
+    // look for points inside the triangle of hole Vertex, segment intersection and endpoint;
+    // if there are no points found, we have a valid connection;
+    // otherwise choose the Vertex of the minimum angle with the ray as connection Vertex
+
+    const Node* stop = m;
+    double tanMin = std::numeric_limits<double>::infinity();
+    double tanCur = 0;
+
+    p = m;
+    double mx = m->x;
+    double my = m->y;
+
+    do {
+        if (hx >= p->x && p->x >= mx && hx != p->x &&
+            pointInTriangle(hy < my ? hx : qx, hy, mx, my, hy < my ? qx : hx, hy, p->x, p->y)) {
+
+            tanCur = std::abs(hy - p->y) / (hx - p->x); // tangential
+
+            if (locallyInside(p, hole) &&
+                (tanCur < tanMin || (tanCur == tanMin && (p->x > m->x || sectorContainsSector(m, p))))) {
+                m = p;
+                tanMin = tanCur;
+            }
+        }
+
+        p = p->next;
+    } while (p != stop);
+
+    return m;
+}
+
+// whether sector in vertex m contains sector in vertex p in the same coordinates
+template <typename N>
+bool Earcut<N>::sectorContainsSector(const Node* m, const Node* p) {
+    return area(m->prev, m, p->prev) < 0 && area(p->next, m, m->next) < 0;
+}
+
+// interlink polygon nodes in z-order
+template <typename N>
+void Earcut<N>::indexCurve(Node* start) {
+    assert(start);
+    Node* p = start;
+
+    do {
+        p->z = p->z ? p->z : zOrder(p->x, p->y);
+        p->prevZ = p->prev;
+        p->nextZ = p->next;
+        p = p->next;
+    } while (p != start);
+
+    p->prevZ->nextZ = nullptr;
+    p->prevZ = nullptr;
+
+    sortLinked(p);
+}
+
+// Simon Tatham's linked list merge sort algorithm
+// http://www.chiark.greenend.org.uk/~sgtatham/algorithms/listsort.html
+template <typename N>
+typename Earcut<N>::Node*
+Earcut<N>::sortLinked(Node* list) {
+    assert(list);
+    Node* p;
+    Node* q;
+    Node* e;
+    Node* tail;
+    int i, numMerges, pSize, qSize;
+    int inSize = 1;
+
+    for (;;) {
+        p = list;
+        list = nullptr;
+        tail = nullptr;
+        numMerges = 0;
+
+        while (p) {
+            numMerges++;
+            q = p;
+            pSize = 0;
+            for (i = 0; i < inSize; i++) {
+                pSize++;
+                q = q->nextZ;
+                if (!q) break;
+            }
+
+            qSize = inSize;
+
+            while (pSize > 0 || (qSize > 0 && q)) {
+
+                if (pSize == 0) {
+                    e = q;
+                    q = q->nextZ;
+                    qSize--;
+                } else if (qSize == 0 || !q) {
+                    e = p;
+                    p = p->nextZ;
+                    pSize--;
+                } else if (p->z <= q->z) {
+                    e = p;
+                    p = p->nextZ;
+                    pSize--;
+                } else {
+                    e = q;
+                    q = q->nextZ;
+                    qSize--;
+                }
+
+                if (tail) tail->nextZ = e;
+                else list = e;
+
+                e->prevZ = tail;
+                tail = e;
+            }
+
+            p = q;
+        }
+
+        tail->nextZ = nullptr;
+
+        if (numMerges <= 1) return list;
+
+        inSize *= 2;
+    }
+}
+
+// z-order of a Vertex given coords and size of the data bounding box
+template <typename N>
+int32_t Earcut<N>::zOrder(const double x_, const double y_) {
+    // coords are transformed into non-negative 15-bit integer range
+    int32_t x = static_cast<int32_t>(32767.0 * (x_ - minX) * inv_size);
+    int32_t y = static_cast<int32_t>(32767.0 * (y_ - minY) * inv_size);
+
+    x = (x | (x << 8)) & 0x00FF00FF;
+    x = (x | (x << 4)) & 0x0F0F0F0F;
+    x = (x | (x << 2)) & 0x33333333;
+    x = (x | (x << 1)) & 0x55555555;
+
+    y = (y | (y << 8)) & 0x00FF00FF;
+    y = (y | (y << 4)) & 0x0F0F0F0F;
+    y = (y | (y << 2)) & 0x33333333;
+    y = (y | (y << 1)) & 0x55555555;
+
+    return x | (y << 1);
+}
+
+// find the leftmost node of a polygon ring
+template <typename N>
+typename Earcut<N>::Node*
+Earcut<N>::getLeftmost(Node* start) {
+    Node* p = start;
+    Node* leftmost = start;
+    do {
+        if (p->x < leftmost->x || (p->x == leftmost->x && p->y < leftmost->y))
+            leftmost = p;
+        p = p->next;
+    } while (p != start);
+
+    return leftmost;
+}
+
+// check if a point lies within a convex triangle
+template <typename N>
+bool Earcut<N>::pointInTriangle(double ax, double ay, double bx, double by, double cx, double cy, double px, double py) const {
+    return (cx - px) * (ay - py) - (ax - px) * (cy - py) >= 0 &&
+           (ax - px) * (by - py) - (bx - px) * (ay - py) >= 0 &&
+           (bx - px) * (cy - py) - (cx - px) * (by - py) >= 0;
+}
+
+// check if a diagonal between two polygon nodes is valid (lies in polygon interior)
+template <typename N>
+bool Earcut<N>::isValidDiagonal(Node* a, Node* b) {
+    return a->next->i != b->i && a->prev->i != b->i && !intersectsPolygon(a, b) && // dones't intersect other edges
+           ((locallyInside(a, b) && locallyInside(b, a) && middleInside(a, b) && // locally visible
+            (area(a->prev, a, b->prev) != 0.0 || area(a, b->prev, b) != 0.0)) || // does not create opposite-facing sectors
+            (equals(a, b) && area(a->prev, a, a->next) > 0 && area(b->prev, b, b->next) > 0)); // special zero-length case
+}
+
+// signed area of a triangle
+template <typename N>
+double Earcut<N>::area(const Node* p, const Node* q, const Node* r) const {
+    return (q->y - p->y) * (r->x - q->x) - (q->x - p->x) * (r->y - q->y);
+}
+
+// check if two points are equal
+template <typename N>
+bool Earcut<N>::equals(const Node* p1, const Node* p2) {
+    return p1->x == p2->x && p1->y == p2->y;
+}
+
+// check if two segments intersect
+template <typename N>
+bool Earcut<N>::intersects(const Node* p1, const Node* q1, const Node* p2, const Node* q2) {
+    int o1 = sign(area(p1, q1, p2));
+    int o2 = sign(area(p1, q1, q2));
+    int o3 = sign(area(p2, q2, p1));
+    int o4 = sign(area(p2, q2, q1));
+
+    if (o1 != o2 && o3 != o4) return true; // general case
+
+    if (o1 == 0 && onSegment(p1, p2, q1)) return true; // p1, q1 and p2 are collinear and p2 lies on p1q1
+    if (o2 == 0 && onSegment(p1, q2, q1)) return true; // p1, q1 and q2 are collinear and q2 lies on p1q1
+    if (o3 == 0 && onSegment(p2, p1, q2)) return true; // p2, q2 and p1 are collinear and p1 lies on p2q2
+    if (o4 == 0 && onSegment(p2, q1, q2)) return true; // p2, q2 and q1 are collinear and q1 lies on p2q2
+
+    return false;
+}
+
+// for collinear points p, q, r, check if point q lies on segment pr
+template <typename N>
+bool Earcut<N>::onSegment(const Node* p, const Node* q, const Node* r) {
+    return q->x <= std::max<double>(p->x, r->x) &&
+        q->x >= std::min<double>(p->x, r->x) &&
+        q->y <= std::max<double>(p->y, r->y) &&
+        q->y >= std::min<double>(p->y, r->y);
+}
+
+template <typename N>
+int Earcut<N>::sign(double val) {
+    return (0.0 < val) - (val < 0.0);
+}
+
+// check if a polygon diagonal intersects any polygon segments
+template <typename N>
+bool Earcut<N>::intersectsPolygon(const Node* a, const Node* b) {
+    const Node* p = a;
+    do {
+        if (p->i != a->i && p->next->i != a->i && p->i != b->i && p->next->i != b->i &&
+                intersects(p, p->next, a, b)) return true;
+        p = p->next;
+    } while (p != a);
+
+    return false;
+}
+
+// check if a polygon diagonal is locally inside the polygon
+template <typename N>
+bool Earcut<N>::locallyInside(const Node* a, const Node* b) {
+    return area(a->prev, a, a->next) < 0 ?
+        area(a, b, a->next) >= 0 && area(a, a->prev, b) >= 0 :
+        area(a, b, a->prev) < 0 || area(a, a->next, b) < 0;
+}
+
+// check if the middle Vertex of a polygon diagonal is inside the polygon
+template <typename N>
+bool Earcut<N>::middleInside(const Node* a, const Node* b) {
+    const Node* p = a;
+    bool inside = false;
+    double px = (a->x + b->x) / 2;
+    double py = (a->y + b->y) / 2;
+    do {
+        if (((p->y > py) != (p->next->y > py)) && p->next->y != p->y &&
+                (px < (p->next->x - p->x) * (py - p->y) / (p->next->y - p->y) + p->x))
+            inside = !inside;
+        p = p->next;
+    } while (p != a);
+
+    return inside;
+}
+
+// link two polygon vertices with a bridge; if the vertices belong to the same ring, it splits
+// polygon into two; if one belongs to the outer ring and another to a hole, it merges it into a
+// single ring
+template <typename N>
+typename Earcut<N>::Node*
+Earcut<N>::splitPolygon(Node* a, Node* b) {
+    Node* a2 = nodes.construct(a->i, a->x, a->y);
+    Node* b2 = nodes.construct(b->i, b->x, b->y);
+    Node* an = a->next;
+    Node* bp = b->prev;
+
+    a->next = b;
+    b->prev = a;
+
+    a2->next = an;
+    an->prev = a2;
+
+    b2->next = a2;
+    a2->prev = b2;
+
+    bp->next = b2;
+    b2->prev = bp;
+
+    return b2;
+}
+
+// create a node and util::optionally link it with previous one (in a circular doubly linked list)
+template <typename N> template <typename Point>
+typename Earcut<N>::Node*
+Earcut<N>::insertNode(std::size_t i, const Point& pt, Node* last) {
+    Node* p = nodes.construct(static_cast<N>(i), util::nth<0, Point>::get(pt), util::nth<1, Point>::get(pt));
+
+    if (!last) {
+        p->prev = p;
+        p->next = p;
+
+    } else {
+        assert(last);
+        p->next = last->next;
+        p->prev = last;
+        last->next->prev = p;
+        last->next = p;
+    }
+    return p;
+}
+
+template <typename N>
+void Earcut<N>::removeNode(Node* p) {
+    p->next->prev = p->prev;
+    p->prev->next = p->next;
+
+    if (p->prevZ) p->prevZ->nextZ = p->nextZ;
+    if (p->nextZ) p->nextZ->prevZ = p->prevZ;
+}
+}
+
+template <typename N = uint32_t, typename Polygon>
+std::vector<N> earcut(const Polygon& poly) {
+    mapbox::detail::Earcut<N> earcut;
+    earcut(poly);
+    return std::move(earcut.indices);
+}
+}