//----------------------------------------------------------------------------
//  EgalTech 2019-2020
//----------------------------------------------------------------------------
// File : SurfTriMeshBooleans.cpp     Data : 02.01.20     Versione : 2.2a1
// Contenuto : Implementazione delle funzioni booleane per SurfFTrimesh.
//
//
//
// Modifiche : 10.05.19 LM Creazione modulo.
//
//
//----------------------------------------------------------------------------

#include "stdafx.h"
#include "SurfTriMesh.h"
#include "CurveLine.h"
#include "CurveComposite.h"
#include "SurfFlatRegion.h"
#include "DistPointLine.h"
#include "Triangulate.h"
#include "GeoConst.h"
#include "/EgtDev/Include/EgkCurve.h"
#include "/EgtDev/Include/EgkDistPointCurve.h"
#include "/EgtDev/Include/EgkDistPointTria.h"
#include "/EgtDev/Include/EgkIntersLineTria.h"
#include "/EgtDev/Include/EgkIntersTriaTria.h"
#include "/EgtDev/Include/EGkChainCurves.h"
#include "/EgtDev/Include/EGkGeoCollection.h"
#include "/EgtDev/Include/EGkPolygon3d.h"
#include "/EgtDev/Include/EgtPerfCounter.h"
#include "/EgtDev/Include/EgnStringUtils.h"
#include <algorithm>

using namespace std ;

//----------------------------------------------------------------------------
static int
IntersRectangleTriangle( const Point3d& ptP, const Vector3d& vtL1, const Vector3d& vtL2,
                         const Triangle3d& trTria, Point3d& ptStSeg, Point3d& ptEnSeg)
{
  // Definisco i due triangoli formanti il rettangolo
   Triangle3d trTriaA ;
   trTriaA.Set( ptP, ptP + vtL1, ptP + vtL2) ;
   if ( ! trTriaA.Validate())
      return -1 ;
   Triangle3d trTriaB ;
   trTriaB.Set( ptP + vtL1, ptP + vtL1 + vtL2, ptP + vtL2) ;
   if ( ! trTriaB.Validate())
      return -1 ;
  // Interseco il triangolo con il primo dei due triangoli del rettangolo 
   int nIntA = 0 ;
   Point3d ptIntA1, ptIntA2 ;
   TRIA3DVECTOR vTriaA ;
   int nIntTypeA = IntersTriaTria( trTria, trTriaA, ptIntA1, ptIntA2, vTriaA) ;
   if ( FromSpecialToNormal( nIntTypeA) == ITTT_PNT  ||  FromSpecialToNormal( nIntTypeA) == ITTT_VERT)
      nIntA = 1 ;      
   else if ( FromSpecialToNormal( nIntTypeA) == ITTT_YES  ||  FromSpecialToNormal( nIntTypeA) == ITTT_EDGE) {
      nIntA = 2 ;
   }
  // Interseco il triangolo con il secondo dei due triangoli del rettangolo 
   int nIntB = 0 ;
   Point3d ptIntB1, ptIntB2 ;
   TRIA3DVECTOR vTriaB ;
   int nIntTypeB = IntersTriaTria( trTria, trTriaB, ptIntB1, ptIntB2, vTriaB) ;
   if ( FromSpecialToNormal( nIntTypeB) == ITTT_PNT  ||  FromSpecialToNormal( nIntTypeB) == ITTT_VERT)
      nIntB = 1 ;     
   else if ( FromSpecialToNormal( nIntTypeB) == ITTT_YES  ||  FromSpecialToNormal( nIntTypeB) == ITTT_EDGE) {
      nIntB = 2 ;
   }
  // Unisco le due intersezioni
   int nIntTot = nIntA + nIntB ;
   if ( nIntTot == 4) {
      if ( AreSamePointApprox( ptIntA2, ptIntB1)) {
         ptStSeg = ptIntA1 ;
         ptEnSeg = ptIntB2 ;
      }
      else {
         ptStSeg = ptIntB1 ;
         ptEnSeg = ptIntA2 ;
      }
      return 2 ;
   }
   else if ( nIntTot == 3) {
      if ( nIntA == 2) {
         ptStSeg = ptIntA1 ;
         ptEnSeg = ptIntA2 ;
      }
      else {
         ptStSeg = ptIntB1 ;
         ptEnSeg = ptIntB2 ;
      }
      return 2 ;
   }
   else if ( nIntTot == 2) {
      if ( nIntA == 2) {
         ptStSeg = ptIntA1 ;
         ptEnSeg = ptIntA2 ;
      }
      else if ( nIntA == 1) {
         ptStSeg = ptIntA1 ;
         ptEnSeg = ptIntB1 ;
      }
      else {
         ptStSeg = ptIntB1 ;
         ptEnSeg = ptIntB2 ;
      }
      return 2 ;
   }
   else if ( nIntTot == 1) {
      if ( nIntA == 1)
         ptStSeg = ptIntA1 ;
      else
         ptStSeg = ptIntB1 ;
      return 1 ;
   }
   else
      return 0 ;
}

//----------------------------------------------------------------------------
static bool
ChangeStart( const Point3d& ptNewStart, PNTVECTOR& Loop)
{
  // Cerco se esiste un tratto del loop chiuso su cui giace il punto
   int nSeg = - 1 ;  
   for ( int nPt = 0 ; nPt < int( Loop.size())  &&  nSeg == - 1 ; ++ nPt) {
     // Estremi del segmento corrente del loop
      Point3d ptSegSt = Loop[nPt] ;
      Point3d ptSegEn = Loop[( nPt + 1) % int( Loop.size())] ;
     // Vedo se il punto giace sul segmento del loop
      DistPointLine dDistCalc( ptNewStart, ptSegSt, ptSegEn) ;
      double dSqDist ;
      dDistCalc.GetSqDist( dSqDist) ;
      if ( dSqDist < SQ_EPS_SMALL) {
         nSeg = nPt ;
      }
   }
  // Se il punto non sta sul loop, errore
   if ( nSeg == - 1)
      return false ;
  // Verifico che il punto stia su un vertice, in tal caso non devo fare nulla
   bool bOnStart = AreSamePointApprox( Loop[nSeg], ptNewStart) ;
   bool bOnEnd = AreSamePointApprox( Loop[( nSeg + 1) % int( Loop.size())], ptNewStart) ;
   if ( bOnStart  ||  bOnEnd) {
      if ( bOnEnd) {
         ++ nSeg ;
         if ( nSeg % int( Loop.size()) == 0)
            return true ;
      }
      PNTVECTOR vTempVec ;
      for ( int nPt = 0 ; nPt < nSeg ; ++ nPt)
         vTempVec.emplace_back( Loop[nPt]) ;
      int nSize = int( Loop.size()) ;
      for ( int nPt = 0 ; nPt < nSize - nSeg ; ++ nPt) {
         Loop[nPt] = Loop[nPt + nSeg] ;
      }
      for ( int nPt = 0 ; nPt < int( vTempVec.size()) ; ++ nPt) {
         Loop[nPt + nSize - nSeg] = vTempVec[nPt] ;
      }
      return true ;
   }
  // Ridimensiono il loop
   Loop.resize( Loop.size() + 1) ;
  // Copio i primi punti
   PNTVECTOR LoopTemp ;
   for ( int nPt = 0 ; nPt <= nSeg ; ++ 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 - nSeg ;
   for ( int nPt = 1 ; nPt <= nLastPointNum ; ++ nPt) {
      Loop[nPt] = Loop[nPt + nSeg] ;
   }
  // Porto i primi in fondo 
   for ( int nPt = 0 ; nPt < int( LoopTemp.size()) ; ++ nPt) {
      Loop[nPt + nLastPointNum] = LoopTemp[nPt] ;
   }
   return true ;
}

//----------------------------------------------------------------------------
static bool
SplitAtPoint( const Point3d& ptStop, const PNTVECTOR& Loop, PNTVECTOR& Loop1, PNTVECTOR& Loop2)
{
  // Cerco se esiste un tratto del loop chiuso su cui giace il punto
   int nSeg = - 1 ;
   for ( int nPt = 0 ; nPt < int( Loop.size())  &&  nSeg == - 1 ; ++ nPt) {
     // Estremi del segmento corrente del loop
      Point3d ptSegSt = Loop[nPt] ;
      Point3d ptSegEn = Loop[( nPt + 1) % int( Loop.size())] ;
     // Vedo se il punto giace sul segmento del loop
      DistPointLine dDistCalc( ptStop, ptSegSt, ptSegEn) ;
      double dSqDist ;
      dDistCalc.GetSqDist( dSqDist) ;
      if ( dSqDist < SQ_EPS_SMALL) {
         nSeg = nPt ;
      }
   }
  // Se il punto non sta sul loop, errore
   if ( nSeg == - 1)
      return false ;
  // Verifico che il punto stia su un vertice, in tal caso non devo aggiungerlo
   bool bFirst = AreSamePointApprox( Loop[nSeg], ptStop) ;
   bool bLast = AreSamePointApprox( Loop[( nSeg + 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)
      ++ nSeg ;
  // Inglobo fino a nSeg nel primo loop
   for ( int nPt = 0 ; nPt <= nSeg ; ++ 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 {
      Loop2.emplace_back( Loop[nSeg]) ;
   }
  // Inglobo gli ultimi vertici in Loop2
   for ( int nPt = nSeg + 1 ; nPt < int( Loop.size()) ; ++ nPt)
      Loop2.emplace_back( Loop[nPt]) ;
   Loop2.emplace_back( Loop[0]) ;

   return true ;
}

//----------------------------------------------------------------------------
static bool
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 ;
}

//----------------------------------------------------------------------------
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 ;
  // Approssimo la curva con segmenti
   CurveComposite cvCompo ;
   PolyLine PL ;
   if ( ! cvCurve.ApproxWithLines( LIN_TOL_MIN, ANG_TOL_STD_DEG, ICurve::APL_STD, PL)  ||  ! cvCompo.FromPolyLine( PL))
      return false ;
  // 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))
      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 il numero di vertici del triangolo che cadono all'interno della curva
      int nVertInside = 0 ;
      for ( int nV = 0 ; nV < 3 ; ++ nV) {
        // Determino se il vertice cade dentro la curva
         DistPointCurve dstPC( trTria.GetP( nV), cvCompo) ;
         int nSide ;
         dstPC.GetSideAtMinDistPoint( 0, vtExtr, nSide) ;
         if ( nSide == MDS_LEFT  ||  nSide == MDS_ON)
            ++ nVertInside ;
      }
     // Vettore di catene di punti 
      CHAINVECTOR vChain ;
     // Ciclo sui segmenti
      bool bStartInside = false ;
      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) ;
        // 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 != 0) {
           // 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 ;
            if ( nInt == 2) {
               CurInters.ptSt = ptSegSt ;
               CurInters.ptEn = ptSegEn ;
               CurInters.bDegenerate = false ;
            }
            else {
               CurInters.ptSt = ptSegSt ;
               CurInters.ptEn = ptSegSt ;
               CurInters.bDegenerate = true ;
            }
            CurInters.vtOuter = ( ptEn - ptSt) ^ vtExtr ;
            CurInters.vtOuter.Normalize() ;
           // 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 {
            bChain = false ;
         }
         pCrv = cvCompo.GetNextCurve() ;
      }
     // 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 ;
         }
      }
     // semplifico catene formate da punti degeneri
      for ( int nCh = 0 ; nCh < nChainCnt ; ++ nCh) {
         if ( vChain[nCh].size() == 2  &&  ( vChain[nCh][0].bDegenerate  ||  vChain[nCh][1].bDegenerate)) {
            vChain[nCh][0].ptEn = vChain[nCh][1].ptEn ;
            vChain[nCh][0].vtOuter = ( vChain[nCh][0].bDegenerate ? vChain[nCh][1].vtOuter : vChain[nCh][0].vtOuter) ;
            vChain[nCh][0].bDegenerate = AreSamePointApprox( vChain[nCh][0].ptSt, vChain[nCh][0].ptEn) ;
            vChain[nCh].resize( 1) ;
         }
      }
     // 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
      int nDegenerateChainNum = 0 ;
      INTVECTOR vnDegVec ;
      CHAINVECTOR cvClosedChain ;
      CHAINVECTOR cvOpenChain ;
      for ( int nL = 0 ; nL < int( vChain.size()) ; ++ nL) {
         bool bChainDegenerate = false ;
         if ( vChain[nL].size() == 1  &&  AreSamePointApprox( vChain[nL][0].ptSt, vChain[nL][0].ptEn)) {
            bChainDegenerate = true ;
         }
        
         if ( bChainDegenerate)
            ++ nDegenerateChainNum ;
         int nCurLoopLast = max( int( vChain[nL].size()) - 1, 0) ;
         if ( ( ! bChainDegenerate)  &&  AreSamePointApprox( vChain[nL][0].ptSt, vChain[nL][nCurLoopLast].ptEn))
            cvClosedChain.emplace_back( vChain[nL]) ;
         else {
            cvOpenChain.emplace_back( vChain[nL]) ;
            if ( bChainDegenerate)
               vnDegVec.emplace_back( 0) ;
            else
               vnDegVec.emplace_back( 1) ;
         }
      }
      
     // Se più di una catena chiusa oppure catene chiuse e aperte, errore
      if ( cvClosedChain.size() > 1  ||
           ( cvClosedChain.size() > 0  &&  int( cvOpenChain.size()) > nDegenerateChainNum))
         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 ;
            if ( vnDegVec[nLastOpenLoopN] == 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 ;
                     CurSeg.vtOuter = - cvOpenChain[nLastOpenLoopN][nPt].vtOuter ;
                     CurSeg.bDegenerate = cvOpenChain[nLastOpenLoopN][nPt].bDegenerate ;
                     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] ;
                  }
                  int nLastPointLoop2 = int( Loop2.size()) - 1 ;
                  Vector3d vtTest = Loop1[1] - Loop1[0] ;
                  vtTest.Normalize() ;
                  bool bSecondInside = vtTest * cvOpenChain[nLastOpenLoopN][0].vtOuter < 0. ;
                  cvBoundClosedLoopVec[nLoop] = Loop1 ;
                  cvBoundClosedLoopVec[nLoop + 1] = Loop2 ;
                  vbInOut[nLoop] = bSecondInside ;
                  vbInOut[nLoop + 1] = ! bSecondInside ;
                  ++ nLoop ;
               }
            }
           // Degenere
            else {
               Point3d ptProva = 0.5 * ( cvOpenChain[nLastOpenLoopN][0].ptSt + cvOpenChain[nLastOpenLoopN][0].ptEn) ;
               Vector3d vtVecProva = cvOpenChain[nLastOpenLoopN][0].vtOuter ;
               vtVecProva.Normalize( EPS_ZERO) ;
               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][0].ptEn ;
                  // Cerco se esistono dei tratti del loop chiuso corrente che sono
                  // toccati dagli estremi del loop aperto corrente
                  int nCvFirst = -1 ;
                  int nCvSecond = -1 ;
                  for ( int nLine = 0 ; nLine < int( cvBoundClosedLoopVec[nLoop].size())  &&  nCvSecond == -1 ; ++ nLine) {
                     // Estremi del segmento corrente del loop chiuso corrente
                     Point3d ptSegSt = cvBoundClosedLoopVec[nLoop][nLine] ;
                     Point3d ptSegEn = cvBoundClosedLoopVec[nLoop][( nLine + 1) % int( cvBoundClosedLoopVec[nLoop].size())] ;
                     // Vettore congiungente i su definiti punti
                     Vector3d vtClosedLoopSeg = ptSegEn - ptSegSt ;
                     vtClosedLoopSeg.Normalize() ;
                     // Vedo se gli estremi del loop aperto stanno su un segmento del chiuso
                     DistPointLine DistCalc( ptProva, ptSegSt, ptSegEn) ;
                     double dSqDist ;
                     DistCalc.GetSqDist( dSqDist) ;
                     if ( dSqDist < 2 * SQ_EPS_SMALL) {
                        if ( nCvFirst == -1)
                           nCvFirst = nLine ;
                        else
                           nCvSecond = nLine ;
                     }
                  }
                  
                  if ( nCvFirst != nCvSecond  &&  nCvSecond != -1) {
                    // li ordino in senso crescente
                     if ( nCvFirst > nCvSecond)
                        swap( nCvFirst, nCvSecond) ;
                    // punto medio tra primo e secondo
                     int nCount = 0 ;
                     Point3d ptM12 ;
                     for ( int i = nCvFirst + 1 ; i <= nCvSecond ; ++ i) {
                        ptM12 += cvBoundClosedLoopVec[nLoop][i] ;
                        ++ nCount ;
                     }
                     ptM12 /= nCount ;
                    // Distanza quadrata media dei punti tra primo e secondo dal baricentro
                     double dVar12 = 0. ;
                     for ( int i = nCvFirst + 1 ; i <= nCvSecond ; ++ i) {
                        dVar12 += ( cvBoundClosedLoopVec[nLoop][i] - ptM12) * ( cvBoundClosedLoopVec[nLoop][i] - ptM12) ;
                     }
                     dVar12 /= nCount ;
                    // punto medio fra secondo e primo
                     nCount = 0 ;
                     Point3d ptM21 ;
                     for ( int i = nCvSecond + 1 ; i % int( cvBoundClosedLoopVec[nLoop].size()) ; ++ i) {
                        ptM21 += cvBoundClosedLoopVec[nLoop][i] ;
                        ++ nCount ;
                     }
                     for ( int i = 0 ; i <= nCvFirst ; ++ i) {
                        ptM21 += cvBoundClosedLoopVec[nLoop][i] ;
                        ++ nCount ;
                     }
                     ptM21 /= nCount ;
                    // Distanza quadrata media dei punti tra secondo e primo dal baricentro
                     double dVar21 = 0. ;
                     for ( int i = nCvSecond ; i <  i % int( cvBoundClosedLoopVec[nLoop].size()) ; ++ i) {
                        dVar21 += ( cvBoundClosedLoopVec[nLoop][i] - ptM21) * ( cvBoundClosedLoopVec[nLoop][i] - ptM21) ;
                        ++ nCount ;
                     }
                     for ( int i = 0 ; i <= nCvFirst ; ++ i) {
                        dVar21 += ( cvBoundClosedLoopVec[nLoop][i] - ptM21) * ( cvBoundClosedLoopVec[nLoop][i] - ptM21) ;
                        ++ nCount ;
                     }
                     dVar21 /= nCount ;
                    // elimino i punti dalla parte non valida
                     if ( dVar12 > dVar21) {
                       // assegno i nuovi valori
                        cvBoundClosedLoopVec[nLoop][nCvFirst] = ptProva ;
                        cvBoundClosedLoopVec[nLoop][( nCvSecond + 1) % int(cvBoundClosedLoopVec[nLoop].size())] = ptProva ;
                       // numero totale di punti
                        int nPntTot = int( cvBoundClosedLoopVec[nLoop].size());
                       // elimino i punti superflui dopo
                        for ( int i = nPntTot - 1 ; i > nCvSecond + 1 ; -- i)
                           cvBoundClosedLoopVec[nLoop].pop_back() ;
                       // elimino i punti superflui prima
                        for ( int i = 0 ; i < nCvFirst ; ++ i)
                           cvBoundClosedLoopVec[nLoop].erase( cvBoundClosedLoopVec[nLoop].begin()) ;
                       // verifico se questo punto è dalla parte valida o no
                        bool bC12 = ( ( ptM12 - ptProva) * vtVecProva < 0) ;
                        vbInOut[nLoop] = bC12 ;
                     }
                     else {
                       // assegno i nuovi valori
                        cvBoundClosedLoopVec[nLoop][nCvFirst + 1] = ptProva ;
                        cvBoundClosedLoopVec[nLoop][nCvSecond] = ptProva ;
                       // elimino i punti superflui intermedi
                        for ( int i = nCvFirst + 2 ; i < nCvSecond ; ++ i)
                           cvBoundClosedLoopVec[nLoop].erase( cvBoundClosedLoopVec[nLoop].begin() + nCvFirst + 2) ;
                       // verifico se questo punto è dalla parte valida o no
                        bool bC21 = ( ( ptM21 - ptProva) * vtVecProva < 0) ;
                        vbInOut[nLoop] = bC21 ;
                     }
                  }
               }
            }
            vnDegVec.resize( nLastOpenLoopN) ;
            cvOpenChain.resize( nLastOpenLoopN) ;
         }
        // Rimuovo il triangolo corrente
         RemoveTriangle( nT) ;
        // 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 ;
                  }
               }
            }
         }
      }
      else if ( nVertInside == 0)
         RemoveTriangle( nT) ;
   }

  // Se avvenuta modifica, aggiorno tutto
   if ( bModif)
      return ( AdjustVertices()  &&  DoCompacting()) ;

   return true ;
}

//----------------------------------------------------------------------------
bool
SurfTriMesh::DecomposeLoop( CHAINVECTOR& cvOpenChain, INTVECTOR& vnDegVec, PNTMATRIX& cvBoundClosedLoopVec, BOOLVECTOR& vbInOut)
{
  // Divido il loop di partenza in sotto-loop
   int nIterationCount = 0 ;
   while ( cvOpenChain.size() > 0) {
      bool bLoopSplitted = false ;
      int nLastOpenLoopN = int( cvOpenChain.size()) - 1 ;
      if ( vnDegVec[nLastOpenLoopN] == 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 ;
            bLoopSplitted = true ;
            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 ;
               CurSeg.vtOuter = - cvOpenChain[nLastOpenLoopN][nPt].vtOuter ;
               CurSeg.bDegenerate = cvOpenChain[nLastOpenLoopN][nPt].bDegenerate ;
               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] ;
            }
            int nLastPointLoop2 = int( Loop2.size()) - 1 ;
            Vector3d vtTest = Loop1[1] - Loop1[0] ;
            vtTest.Normalize() ;
            bool bSecondInside = vtTest * cvOpenChain[nLastOpenLoopN][0].vtOuter < 0. ;
            cvBoundClosedLoopVec[nLoop] = Loop1 ;
            cvBoundClosedLoopVec[nLoop + 1] = Loop2 ;
            vbInOut[nLoop] = bSecondInside ;
            vbInOut[nLoop + 1] = ! bSecondInside ;
            ++ nLoop ;
         }
      }
     // Degenere
      else {
         Point3d ptProva = 0.5 * ( cvOpenChain[nLastOpenLoopN][0].ptSt + cvOpenChain[nLastOpenLoopN][0].ptEn) ;
         Vector3d vtVecProva = cvOpenChain[nLastOpenLoopN][0].vtOuter ;
         vtVecProva.Normalize( EPS_ZERO) ;
         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][0].ptEn ;
           // Cerco se esistono dei tratti del loop chiuso corrente che sono
           // toccati dagli estremi del loop aperto corrente
            int nCvFirst = - 1 ;
            int nCvSecond = - 1 ;
            for ( int nLine = 0 ; nLine < int( cvBoundClosedLoopVec[nLoop].size())  &&  nCvSecond == - 1 ; ++ nLine) {
              // Estremi del segmento corrente del loop chiuso corrente
               Point3d ptSegSt = cvBoundClosedLoopVec[nLoop][nLine] ;
               Point3d ptSegEn = cvBoundClosedLoopVec[nLoop][( nLine + 1) % int(cvBoundClosedLoopVec[nLoop].size())] ;
              // Vettore congiungente i su definiti punti
               Vector3d vtClosedLoopSeg = ptSegEn - ptSegSt ;
               vtClosedLoopSeg.Normalize() ;
              // Vedo se gli estremi del loop aperto stanno su un segmento del chiuso
               DistPointLine DistCalc( ptProva, ptSegSt, ptSegEn) ;
               double dSqDist ;
               DistCalc.GetSqDist( dSqDist) ;
               if ( dSqDist < 2 * SQ_EPS_SMALL) {
                  if ( nCvFirst == - 1)
                     nCvFirst = nLine ;
                  else
                     nCvSecond = nLine ;
               }
            }

            if ( nCvFirst != nCvSecond  &&  nCvSecond != - 1) {
              // li ordino in senso crescente
               if ( nCvFirst > nCvSecond)
                  swap( nCvFirst, nCvSecond) ;
              // punto medio tra primo e secondo
               int nCount = 0 ;
               Point3d ptM12 ;
               for ( int i = nCvFirst + 1 ; i <= nCvSecond ; ++ i) {
                  ptM12 += cvBoundClosedLoopVec[nLoop][i] ;
                  ++ nCount ;
               }
               ptM12 /= nCount ;
              // Distanza quadrata media dei punti tra primo e secondo dal baricentro
               double dVar12 = 0. ;
               for ( int i = nCvFirst + 1 ; i <= nCvSecond ; ++ i) {
                  dVar12 += ( cvBoundClosedLoopVec[nLoop][i] - ptM12) * ( cvBoundClosedLoopVec[nLoop][i] - ptM12) ;
               }
               dVar12 /= nCount ;
              // punto medio fra secondo e primo
               nCount = 0 ;
               Point3d ptM21 ;
               for ( int i = nCvSecond + 1 ; i % int( cvBoundClosedLoopVec[nLoop].size()) ; ++ i) {
                  ptM21 += cvBoundClosedLoopVec[nLoop][i] ;
                  ++ nCount ;
               }
               for ( int i = 0 ; i <= nCvFirst ; ++ i) {
                  ptM21 += cvBoundClosedLoopVec[nLoop][i] ;
                  ++ nCount ;
               }
               ptM21 /= nCount ;
              // Distanza quadrata media dei punti tra secondo e primo dal baricentro
               double dVar21 = 0. ;
               for ( int i = nCvSecond ; i < i % int( cvBoundClosedLoopVec[nLoop].size()) ; ++ i) {
                  dVar21 += ( cvBoundClosedLoopVec[nLoop][i] - ptM21) * ( cvBoundClosedLoopVec[nLoop][i] - ptM21) ;
                  ++ nCount ;
               }
               for ( int i = 0 ; i <= nCvFirst ; ++ i) {
                  dVar21 += ( cvBoundClosedLoopVec[nLoop][i] - ptM21) * ( cvBoundClosedLoopVec[nLoop][i] - ptM21) ;
                  ++ nCount ;
               }
               dVar21 /= nCount ;
              // elimino i punti dalla parte non valida
               if ( dVar12 > dVar21) {
                 // assegno i nuovi valori
                  cvBoundClosedLoopVec[nLoop][nCvFirst] = ptProva ;
                  cvBoundClosedLoopVec[nLoop][( nCvSecond + 1) % int( cvBoundClosedLoopVec[nLoop].size())] = ptProva ;
                 // numero totale di punti
                  int nPntTot = int( cvBoundClosedLoopVec[nLoop].size()) ;
                 // elimino i punti superflui dopo
                  for ( int i = nPntTot - 1 ; i > nCvSecond + 1 ; -- i)
                     cvBoundClosedLoopVec[nLoop].pop_back() ;
                 // elimino i punti superflui prima
                  for ( int i = 0 ; i < nCvFirst ; ++ i)
                     cvBoundClosedLoopVec[nLoop].erase( cvBoundClosedLoopVec[nLoop].begin()) ;
                 // verifico se questo punto è dalla parte valida o no
                  bool bC12 = ( ( ptM12 - ptProva) * vtVecProva < 0) ;
                  vbInOut[nLoop] = bC12 ;
               }
               else {
                 // assegno i nuovi valori
                  cvBoundClosedLoopVec[nLoop][nCvFirst + 1] = ptProva ;
                  cvBoundClosedLoopVec[nLoop][nCvSecond] = ptProva ;
                 // elimino i punti superflui intermedi
                  for ( int i = nCvFirst + 2 ; i < nCvSecond ; ++ i)
                     cvBoundClosedLoopVec[nLoop].erase( cvBoundClosedLoopVec[nLoop].begin() + nCvFirst + 2) ;
                 // verifico se questo punto è dalla parte valida o no
                  bool bC21 = ( ( ptM21 - ptProva) * vtVecProva < 0) ;
                  vbInOut[nLoop] = bC21 ;
               }
               bLoopSplitted = true ;
            }
         }
      }
      if ( ! bLoopSplitted) {
         int nCurDeg = vnDegVec[nLastOpenLoopN] ;
         vnDegVec.emplace( vnDegVec.begin(), nCurDeg) ;
         Chain CurChain ;
         for ( int nCrChSeg = 0 ; nCrChSeg < int( cvOpenChain[nLastOpenLoopN].size()) ; ++ nCrChSeg) {
            IntSegment CurChainSeg ;
            CurChainSeg.ptSt = cvOpenChain[nLastOpenLoopN][nCrChSeg].ptSt ;
            CurChainSeg.ptEn = cvOpenChain[nLastOpenLoopN][nCrChSeg].ptEn ;
            CurChainSeg.vtOuter = cvOpenChain[nLastOpenLoopN][nCrChSeg].vtOuter ;
            CurChainSeg.bDegenerate = cvOpenChain[nLastOpenLoopN][nCrChSeg].bDegenerate ;
            CurChain.emplace_back( CurChainSeg) ;
         }
         cvOpenChain.emplace( cvOpenChain.begin(), CurChain) ;
         ++ nLastOpenLoopN ;  
         ++ nIterationCount ;
      }
      else
         nIterationCount = 0 ;
      vnDegVec.resize( nLastOpenLoopN) ;
      cvOpenChain.resize( nLastOpenLoopN) ;
      if ( nIterationCount > int( cvOpenChain.size()) + 2)
         return false ;
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
SurfTriMesh::RetriangulationForBooleanOperation( CHAINMAP& LoopLines, TRIA3DVECTORMAP& Ambiguos,
                                                 SurfTriMesh& Surf, bool& bModif)
{
  // La superficie deve essere valida
   if ( ! Surf.IsValid())
      return false ;
  // Ritriangolarizzo i triangoli
   for ( auto it = LoopLines.begin() ; it != LoopLines.end() ; ++ it) {
      for ( int nS1 = 0 ; nS1 < int( it->second.size()) - 1 ; ++ nS1) {
         for ( int nS2 = nS1 + 1 ; nS2 < int( it->second.size()) ; ++ nS2) {
            if ( AreSamePointApprox( it->second[nS1].ptSt, it->second[nS2].ptEn)  &&
                 AreSamePointApprox( it->second[nS1].ptEn, it->second[nS2].ptSt)  &&
               it->second[nS1].vtOuter * it->second[nS2].vtOuter < - EPS_SMALL) {
               it->second.erase( it->second.begin() + nS2) ;
               it->second.erase( it->second.begin() + nS1) ;
               -- nS1 ;
               break ;
            }
         }
      }

      if ( int( it->second.size()) == 0)
         continue ;
     // Se il triangolo è stato sottoposto a ritriangolazione, le sue componenti sono classificabili come dentro-fuori.
     // Lo tolgo dall'insieme dei triangoli ambigui (intersezione edge-edge)
      else {
         auto itS = Ambiguos.find( it->first) ;
         if ( itS != Ambiguos.end()) {
            Ambiguos.erase( itS) ;
         }
      }

     // Recupero il triangolo
      Triangle3d trTria ;
      Surf.GetTriangle( it->first, trTria) ;
     // Lo rimuovo dalla mesh
      Surf.RemoveTriangle( it->first) ;
      bModif = true ;
      CHAINVECTOR vChain ;
     // Creo i loop
      ChainCurves LoopCreator ;
      LoopCreator.Init( false, EPS_SMALL, int( it->second.size())) ;
     // Carico le curve per concatenarle
      for ( int nCv = 0 ; nCv < int( it->second.size()); ++ nCv) {
         Point3d ptSt = it->second[nCv].ptSt ;
         Point3d ptEn = it->second[nCv].ptEn ;
         Vector3d vtDir = ptEn - ptSt ;
         vtDir.Normalize() ;
         LoopCreator.AddCurve( nCv + 1, ptSt, vtDir, ptEn, vtDir) ;
      }
     // Recupero i concatenamenti
      INTVECTOR vIds ;
      Point3d ptNearStart ;
      while ( LoopCreator.GetChainFromNear( ptNearStart, false, vIds)) {
         Chain chTemp ;
         for ( auto i : vIds) {
           // Aggiungo la linea alla curva composta.
            chTemp.emplace_back( it->second[i - 1]) ;
         }
         vChain.emplace_back( chTemp) ;
      }
     // Lavoro su loop e catene per regolarizzarle
      int 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 ;
         }
      }
     // semplifico catene formate da punti degeneri
      for ( int nCh = 0 ; nCh < nChainCnt ; ++ nCh) {
         if ( vChain[nCh].size() == 2  &&  ( vChain[nCh][0].bDegenerate  ||  vChain[nCh][1].bDegenerate)) {
            vChain[nCh][0].ptEn = vChain[nCh][1].ptEn ;
            vChain[nCh][0].vtOuter = ( vChain[nCh][0].bDegenerate ? vChain[nCh][1].vtOuter : vChain[nCh][0].vtOuter) ;
            vChain[nCh][0].bDegenerate = AreSamePointApprox( vChain[nCh][0].ptSt, vChain[nCh][0].ptEn) ;
            vChain[nCh].resize( 1) ;
         }
      }
     // 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
      int nDegenerateChainNum = 0 ;
      INTVECTOR vnDegVec ;
      CHAINVECTOR cvClosedChain ;
      CHAINVECTOR cvOpenChain ;
      for ( int nL = 0 ; nL < int( vChain.size()) ; ++ nL) {
         bool bChainDegenerate = false ;
         if ( vChain[nL].size() == 1  &&  AreSamePointApprox( vChain[nL][0].ptSt, vChain[nL][0].ptEn)) {
            bChainDegenerate = true ;
         }

         if ( bChainDegenerate)
            ++ nDegenerateChainNum ;
         int nCurLoopLast = max( int( vChain[nL].size()) - 1, 0) ;
         if ( ( ! bChainDegenerate)  &&  AreSamePointApprox( vChain[nL][0].ptSt, vChain[nL][nCurLoopLast].ptEn))
            cvClosedChain.emplace_back( vChain[nL]) ;
         else {
            cvOpenChain.emplace_back( vChain[nL]) ;
            if ( bChainDegenerate)
               vnDegVec.emplace_back( 0) ;
            else
               vnDegVec.emplace_back( 1) ;
         }
      }

      for ( int nCh1 = 0 ; nCh1 < int( cvOpenChain.size()) - 1 ; ++ nCh1) {
         for ( int nCh2 = nCh1 + 1 ; nCh2 < int( cvOpenChain.size()) ; ++ nCh2) {
            int nChainSize1 = int( cvOpenChain[nCh1].size()) ;
            int nChainSize2 = int( cvOpenChain[nCh2].size()) ;
            int nSameSeg = 0 ;
            for ( int nSeg1 = 0 ; nSeg1 < nChainSize1 ; ++ nSeg1) {
               for ( int nSeg2 = 0 ; nSeg2 < nChainSize2 ; ++ nSeg2) {
                  if ( AreSamePointExact( cvOpenChain[nCh1][nSeg1].ptSt, cvOpenChain[nCh2][nSeg2].ptSt)  &&
                       AreSamePointExact( cvOpenChain[nCh1][nSeg1].ptEn, cvOpenChain[nCh2][nSeg2].ptEn)  &&
                       AreSameVectorExact( cvOpenChain[nCh1][nSeg1].vtOuter, cvOpenChain[nCh2][nSeg2].vtOuter)) {
                     ++ nSameSeg ;
                  }
               }
            }
            if ( nChainSize1 == nSameSeg) {
               cvOpenChain.erase( cvOpenChain.begin() + nCh1) ;
               vnDegVec.erase( vnDegVec.begin() + nCh1) ;
               -- nCh1 ;
            }
            else if ( nChainSize2 == nSameSeg) {
               cvOpenChain.erase( cvOpenChain.begin() + nCh2) ;
               vnDegVec.erase( vnDegVec.begin() + nCh2) ;
               -- nCh2 ;
            }
         }
      }

     // 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 usando le catene
      bool bDecomposed = DecomposeLoop( cvOpenChain, vnDegVec, cvBoundClosedLoopVec, vbInOut) ;
     // Rimuovo il triangolo corrente
      Surf.RemoveTriangle( it->first) ;
     // 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]) ;

        // Assegno ai loop trovati i rispettivi interni
        // Assumo che i loop interni a uno dei loop creati fino ad'ora siano tutti sullo stesso livello.
        // Il caso generale si risolve con una struttura ad albero in cui il nodi corrispondente a un  
        // loop è figlio del nodo corrispondente al loop che lo contiene.
         INTVECTOR vInnerLoop ;
         for ( int nCLI = 0 ; nCLI < int( cvClosedChain.size()) ; ++ nCLI) {
            Point3d ptLoopStart = cvClosedChain[nCLI][0].ptSt ;
            double dMinDist = DBL_MAX ;
            Point3d ptMinDist ;
            bool bPointOnSt = false ;
            bool bPointOnEn = false ;
            int nSegNum = 0 ;
            int nSegMin ;
            Point3d ptS, ptE ;
            bool bContinueS = vplPolyVec[nLoop].GetFirstPoint( ptS) ;
            bool bContinueE = vplPolyVec[nLoop].GetNextPoint( ptE) ;
            while ( bContinueS  &&  bContinueE) {
               ++ nSegNum ;
               DistPointLine DistCalculator( ptLoopStart, ptS, ptE) ;
               double dDist ;
               DistCalculator.GetDist( dDist) ;
               if ( dDist < dMinDist) {
                  DistCalculator.GetMinDistPoint( ptMinDist) ;
                  bPointOnSt = AreSamePointExact( ptMinDist, ptS) ;
                  bPointOnEn = AreSamePointExact( ptMinDist, ptE) ;
                  dMinDist = dDist ;
                  nSegMin = nSegNum ;
               }
               ptS = ptE ;
               bContinueS = bContinueE ;
               bContinueE = vplPolyVec[nLoop].GetNextPoint( ptE) ;
            }
            if ( ! ( bPointOnSt  ||  bPointOnEn)) {
               vplPolyVec[nLoop].GetFirstPoint( ptS) ;
               vplPolyVec[nLoop].GetNextPoint( ptE) ;
               for ( int nSeg = 1 ; nSeg < nSegMin ; ++ nSeg) {
                  ptS = ptE ;
                  vplPolyVec[nLoop].GetNextPoint( ptE) ;
               }
               Vector3d vtTan = ptE - ptS ;
               vtTan.Normalize() ;
               Vector3d vtOut = vtTan ^ trTria.GetN() ;
               Point3d ptMinDist ;
               DistPointLine DistCalculator( ptLoopStart, ptS, ptE) ;
               DistCalculator.GetMinDistPoint( ptMinDist) ;
               double dMinDistDot = ( ptLoopStart - ptMinDist) * vtOut ;
               if ( dMinDistDot < 0.)
                  vInnerLoop.emplace_back( nCLI) ;
            }
            else if ( bPointOnSt) {
               Point3d ptPrevS, ptPrevE ;
               if ( nSegMin == 1) {
                  vplPolyVec[nLoop].GetFirstPoint( ptS) ;
                  vplPolyVec[nLoop].GetNextPoint( ptE) ;
                  vplPolyVec[nLoop].GetLastPoint( ptPrevE) ;
                  vplPolyVec[nLoop].GetPrevPoint( ptPrevS) ;
               }
               else {
                  -- nSegMin ;
                  vplPolyVec[nLoop].GetFirstPoint( ptPrevS) ;
                  vplPolyVec[nLoop].GetNextPoint( ptPrevE) ;
                  for ( int nSeg = 1 ; nSeg < nSegMin ; ++ nSeg) {
                     ptPrevS = ptPrevE ;
                     vplPolyVec[nLoop].GetNextPoint( ptPrevE) ;
                  }
                  ptS = ptPrevE ;
                  vplPolyVec[nLoop].GetNextPoint( ptE) ;
               }
               Vector3d vtTan = ptE - ptS ;
               vtTan.Normalize() ;
               Vector3d vtTanPrev = ptPrevE - ptPrevS ;
               vtTanPrev.Normalize() ;
               Vector3d vtBisector = 0.5 * ( vtTan + vtTanPrev) ^ trTria.GetN() ;
               vtBisector.Normalize() ;
               double dMinDistDot = ( ptLoopStart - ptMinDist) * vtBisector ;
               if ( dMinDistDot < 0.)
                  vInnerLoop.emplace_back( nCLI) ;
            }
            else if ( bPointOnEn) {
               Point3d ptLast ;
               vplPolyVec[nLoop].GetLastPoint( ptLast) ;
               vplPolyVec[nLoop].GetFirstPoint( ptS) ;
               vplPolyVec[nLoop].GetNextPoint( ptE) ;
               for ( int nSeg = 1 ; nSeg < nSegMin ; ++ nSeg) {
                  ptS = ptE ;
                  vplPolyVec[nLoop].GetNextPoint( ptE) ;
               }
               Point3d ptNextS, ptNextE ;
               if ( AreSamePointExact( ptE, ptLast)) {
                  vplPolyVec[nLoop].GetFirstPoint( ptNextS) ;
                  vplPolyVec[nLoop].GetNextPoint( ptNextE) ;
               }
               else {
                  vplPolyVec[nLoop].GetNextPoint( ptNextS) ;
                  vplPolyVec[nLoop].GetNextPoint( ptNextE) ;
               }
               Vector3d vtTan = ptE - ptS ;
               vtTan.Normalize() ;
               Vector3d vtTanNext = ptNextE - ptNextS ;
               vtTanNext.Normalize() ;
               Vector3d vtBisector = 0.5 * ( vtTan + vtTanNext) ^ trTria.GetN() ;
               vtBisector.Normalize() ;
               double dMinDistDot = ( ptLoopStart - ptMinDist) * vtBisector ;
               if ( dMinDistDot < 0.)
                  vInnerLoop.emplace_back( nCLI) ;
            }
         }

         if ( vInnerLoop.size() == 0) {
           // Eseguo triangolazione
            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] = { Surf.AddVertex( vPt[nNewTriaVertId[0]]),
                                    Surf.AddVertex( vPt[nNewTriaVertId[1]]),
                                    Surf.AddVertex( vPt[nNewTriaVertId[2]]) } ;
                  int nNewTriaNum = Surf.AddTriangle( nNewId) ;
                  if ( IsValidSvt( nNewTriaNum)) {
                     if ( vbInOut[nLoop])
                        Surf.m_vTria[nNewTriaNum].nTempPart = 1 ;
                     else
                        Surf.m_vTria[nNewTriaNum].nTempPart = - 1 ;
                     bModif = true ;
                  }
               }
            }
         }
         else {
            POLYLINEVECTOR vPolygons ;
            vPolygons.emplace_back( vplPolyVec[nLoop]) ;
            for ( int nL = 0 ; nL < int( vInnerLoop.size()) ; ++ nL) {
               PolyLine CurLoop ;
               for ( int nV = 0 ; nV < int( cvClosedChain[vInnerLoop[nL]].size()) ; ++ nV) {
                  CurLoop.AddUPoint( 0., cvClosedChain[vInnerLoop[nL]][nV].ptSt) ;
               }
               CurLoop.AddUPoint( 0., cvClosedChain[vInnerLoop[nL]][0].ptSt) ;
               vPolygons.emplace_back( CurLoop) ;
            }

            Polygon3d pgPol ;
            pgPol.FromPolyLine(vPolygons[1]) ;

            bool bCodirectedNormals = trTria.GetN() * pgPol.GetVersN() > 0. ;
            if ( bCodirectedNormals) {
               for ( int nL = 1 ; nL < int( vPolygons.size()) ; ++ nL) {
                  vPolygons[nL].Invert() ;
               }
            }

            PNTVECTOR vPt ;
            INTVECTOR vTr ;
            if ( Triangulate().Make( vPolygons, 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] = { Surf.AddVertex( vPt[nNewTriaVertId[0]]),
                                    Surf.AddVertex( vPt[nNewTriaVertId[1]]),
                                    Surf.AddVertex( vPt[nNewTriaVertId[2]])} ;
                  int nNewTriaNum = Surf.AddTriangle( nNewId) ;
                  if ( IsValidSvt( nNewTriaNum)) {
                     if ( bCodirectedNormals)
                        Surf.m_vTria[nNewTriaNum].nTempPart = -1 ;
                     else
                        Surf.m_vTria[nNewTriaNum].nTempPart = 1 ;
                     bModif = true ;
                  }
               }
            }

            for ( int nL = 1 ; nL < int( vPolygons.size()) ; ++ nL) {
               vPolygons[nL].Invert() ;
               if ( Triangulate().Make( vPolygons[nL], 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] = { Surf.AddVertex( vPt[nNewTriaVertId[0]]),
                                       Surf.AddVertex( vPt[nNewTriaVertId[1]]),
                                       Surf.AddVertex( vPt[nNewTriaVertId[2]])} ;
                     int nNewTriaNum = Surf.AddTriangle( nNewId) ;
                     if ( IsValidSvt( nNewTriaNum)) {
                        if ( bCodirectedNormals)
                           Surf.m_vTria[nNewTriaNum].nTempPart = 1 ;
                        else
                           Surf.m_vTria[nNewTriaNum].nTempPart = -1 ;
                        bModif = true ;
                     }
                  }
               }
            }
         }
         vInnerLoop.resize( 0) ;
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
SurfTriMesh::AmbiguosTriangleManager( TRIA3DVECTORMAP& Ambiguos, SurfTriMesh& Surf)
{ 
   for ( auto it = Ambiguos.begin() ; it != Ambiguos.end() ; ++ it) {
      Triangle3d trTria ;
      Surf.GetTriangle( it->first, trTria) ;
      trTria.Validate() ;
      Point3d ptBar = ( trTria.GetP( 0) + trTria.GetP( 1) + trTria.GetP( 2)) / 3 ;
      double dMinDist = DBL_MAX ;
      int nTriaIndex = - 1 ;
      for ( int nOthSurfT = 0 ; nOthSurfT < int( it->second.size()) ; ++ nOthSurfT) {
         Triangle3d trOthSurfTria = it->second[nOthSurfT] ;
         double dDot = ( ptBar - trOthSurfTria.GetP( 0)) * trOthSurfTria.GetN() ;
         if ( abs( dDot) > EPS_SMALL) {
            double dDist ;
            if ( DistPointTriangle( ptBar, trOthSurfTria).GetDist( dDist)  &&  dDist < dMinDist) {
               nTriaIndex = nOthSurfT ;
               dMinDist = dDist ;
            }
         }
      }
      if ( nTriaIndex != -1) {
         Triangle3d trOthSurfTria = it->second[nTriaIndex] ;
         trOthSurfTria.Validate() ;
         double dDot = ( ptBar - trOthSurfTria.GetP( 0)) * trOthSurfTria.GetN() ;
         Surf.m_vTria[it->first].nTempPart = ( dDot < 0 ? 1 : -1) ;
      }
   }   
   return true ;
}

//----------------------------------------------------------------------------
bool
SurfTriMesh::IntersectTriMeshTriangle( SurfTriMesh& Other)
{   
   bool bModif = false ;
   SurfTriMesh& SurfB = Other ;
  // Le superfici devono essere valide
   if ( m_nStatus != OK  ||  ! SurfB.IsValid())
      return false ;
  // Unordered map dei segmenti di intersezione
   CHAINMAP LineMapA ;
   CHAINMAP LineMapB ;
  // Unordered map dei triangoli ambigui (intersezione edge-edge)
   TRIA3DVECTORMAP AmbiguosA ;
   TRIA3DVECTORMAP AmbiguosB ;
  // Ciclo sui triangoli delle mesh
   int nTriaNumA = GetTriangleSize() ;
   int nTriaNumB = SurfB.GetTriangleSize() ;
  // Setto il triangolo come né fuori né dentro
   for ( int nTA = 0 ; nTA < nTriaNumA ; ++ nTA)
      m_vTria[nTA].nTempPart = 0 ;
   for ( int nTB = 0 ; nTB < nTriaNumB ; ++ nTB)
      SurfB.m_vTria[nTB].nTempPart = 0 ;
  // Resetto e ricalcolo la HashGrid della superficie B
   SurfB.ResetHashGrids3d() ;
   for ( int nTA = 0 ; nTA < nTriaNumA ; ++ nTA) {
     // Se il triangolo A non è valido, continuo
      Triangle3d trTriaA ;
      if ( ! GetTriangle( nTA, trTriaA)  ||  ! trTriaA.Validate( true))
         continue ;
     // Box del triangolo A
      BBox3d b3dTriaA ;
      trTriaA.GetLocalBBox( b3dTriaA) ;
     // Recupero i triangoli di B che interferiscono col box del triangolo di A
      INTVECTOR vNearTria ;
      SurfB.GetAllTriaOverlapBox( b3dTriaA, vNearTria) ;
      bool bNewTriaA = true ;
      for ( int nTB = 0 ; nTB < int( vNearTria.size()) ; ++ nTB) {
        // Se il triangolo B non è valido, continuo
         Triangle3d trTriaB ;
         if ( ! SurfB.GetTriangle( vNearTria[nTB], trTriaB)  ||  ! trTriaB.Validate( true))
            continue ;
        // Interseco i triangoli
         if ( abs( trTriaA.GetN() * trTriaB.GetN()) < 1 - EPS_ZERO) {
            Point3d ptSegSt, ptSegEn ;
            TRIA3DVECTOR vTria ;
            int nIntType = IntersTriaTria( trTriaA, trTriaB, ptSegSt, ptSegEn, vTria) ;
            if ( FromSpecialToNormal( nIntType) != ITTT_NO  &&
                 FromSpecialToNormal( nIntType) != ITTT_OVERLAPS  &&
                 FromSpecialToNormal( nIntType) != ITTTS_VERT_VERT) {
              // Assegno i dati di intersezione
               IntSegment CurInters ;
               if ( FromSpecialToNormal( nIntType) == ITTT_EDGE  ||  FromSpecialToNormal( nIntType) == ITTT_YES) {
                  CurInters.ptSt = ptSegSt ;
                  CurInters.ptEn = ptSegEn ;
                  CurInters.bDegenerate = false ;
               }
               else  {
                  CurInters.ptSt = ptSegSt ;
                  CurInters.ptEn = ptSegSt ;
                  CurInters.bDegenerate = true ;
               }
               CurInters.vtOuter = trTriaB.GetN() ;
               CurInters.vtOuter -= ( ( CurInters.vtOuter * trTriaA.GetN()) * trTriaA.GetN()) ;
               CurInters.vtOuter.Normalize() ;
              // Salvo intersezione per superficie A
               bool bIntOnEndgeA = false ;
               if ( nIntType != ITTTS_EDGE_EDGE_SEG  &&  nIntType != ITTTS_EDGE_INT) {
                  auto itA = LineMapA.find( nTA) ;
                  if ( itA != LineMapA.end()) {
                     itA->second.emplace_back( CurInters) ;
                  }
                  else {
                     Chain chTemp ;
                     chTemp.emplace_back( CurInters) ;
                     LineMapA.emplace( nTA, chTemp) ;
                  }
               }
               else
                  bIntOnEndgeA = true ;

               swap( CurInters.ptSt, CurInters.ptEn) ;
               CurInters.vtOuter = trTriaA.GetN() ;
               CurInters.vtOuter -= ( ( CurInters.vtOuter * trTriaB.GetN()) * trTriaB.GetN()) ;
               CurInters.vtOuter.Normalize() ;

              // Salvo intersezione per superficie B
               bool bIntOnEndgeB = false ;
               if ( nIntType != ITTTS_EDGE_EDGE_SEG  &&  nIntType != ITTTS_INT_EDGE) {
                  auto itB = LineMapB.find( vNearTria[nTB]) ;
                  if ( itB != LineMapB.end()) {
                     itB->second.emplace_back( CurInters) ;
                  }
                  else {
                     Chain chTemp ;
                     chTemp.emplace_back( CurInters) ;
                     LineMapB.emplace( vNearTria[nTB], chTemp) ;
                  }
               }
               else
                  bIntOnEndgeB = true ;
              // Intersezione edge-interno
               if ( bIntOnEndgeA  &&  ! bIntOnEndgeB) {
                  double dMaxDist = 0. ;
                  int nSegMaxDist = - 1 ;
                  for ( int nVA = 0 ; nVA < 3 ; ++ nVA) {
                     double dDist = abs( ( trTriaA.GetP( nVA) - trTriaB.GetP( 0)) * trTriaB.GetN()) ;
                     if ( dMaxDist < dDist) {
                        nSegMaxDist = nVA ;
                        dMaxDist = dDist ;
                     }
                  }
                  if ( nSegMaxDist >= 0) {
                     m_vTria[nTA].nTempPart = ( ( trTriaA.GetP( nSegMaxDist) - trTriaB.GetP( 0)) * trTriaB.GetN() < - EPS_SMALL ? 1 : - 1) ;
                  }
               }
              // Intersezione interno-edge
               else if ( ! bIntOnEndgeA  &&  bIntOnEndgeB) {
                  double dMaxDist = 0. ;
                  int nSegMaxDist = - 1 ;
                  for ( int nVB = 0 ; nVB < 3 ; ++ nVB) {
                     double dDist = abs( ( trTriaB.GetP( nVB) - trTriaA.GetP( 0)) * trTriaA.GetN()) ;
                     if ( dMaxDist < dDist) {
                        nSegMaxDist = nVB ;
                        dMaxDist = dDist ;
                     }
                  }
                  if ( nSegMaxDist >= 0) {
                     SurfB.m_vTria[vNearTria[nTB]].nTempPart = ( ( trTriaB.GetP( nSegMaxDist) - trTriaA.GetP( 0)) * trTriaA.GetN() < - EPS_SMALL ? 1 : - 1) ;
                  }
               }
              // Intersezione edge-edge
               else if ( bIntOnEndgeA  &&  bIntOnEndgeB) {
                  auto itA = AmbiguosA.find( nTA) ;
                  if ( itA == AmbiguosA.end()) {
                     TRIA3DVECTOR vVecTriaB ;
                     vVecTriaB.emplace_back( trTriaB) ;
                     AmbiguosA.emplace( nTA, vVecTriaB) ;
                  }
                  else {                  
                     itA->second.emplace_back( trTriaB) ;
                  }
                  auto itB = AmbiguosB.find( vNearTria[nTB]) ;
                  if ( itB == AmbiguosB.end()) {
                     TRIA3DVECTOR vVecTriaA ;
                     vVecTriaA.emplace_back( trTriaA) ;
                     AmbiguosB.emplace( vNearTria[nTB], vVecTriaA) ;
                  }
                  else {
                     itB->second.emplace_back( trTriaA) ;
                  }
               }
            }
            else {
               ;
            }
         }
      }
   }
   
  // Ritriangolarizzo i triangoli delle superfici
   RetriangulationForBooleanOperation( LineMapA, AmbiguosA, *this, bModif) ;
   RetriangulationForBooleanOperation( LineMapB, AmbiguosB, SurfB, bModif) ;
  
  // Se i triangoli delle superfici non si intersecano, una delle due è totalmente interna o esterna all'altra.
   bool bRetriangulated = true ;
   if ( ! bModif) {
      bRetriangulated = false ;
      int nVertNum = 0 ;
      Point3d ptFirstV ;
      int nCurVert = GetFirstVertex( ptFirstV) ;
      int nInOutNum = 0 ;
      while ( nInOutNum == 0  &&  nCurVert != SVT_NULL) {
         int nTriaNum = - 1 ;
         double dMinDist = DBL_MAX ;
         for ( int nTB = 0 ; nTB < nTriaNumB ; ++ nTB) {
           // Se il triangolo B non è valido, continuo
            Triangle3d trTriaB ;
            if ( ! SurfB.GetTriangle( nTB, trTriaB)  ||  ! trTriaB.Validate( true))
               continue ;
            double dDist ;
            if ( DistPointTriangle( ptFirstV, trTriaB).GetDist( dDist)  &&  dDist < dMinDist) {
               nTriaNum = nTB ;
               dMinDist = dDist ;
            }
         }
         if ( nTriaNum >= 0) {
            Triangle3d trTriaB ;
            SurfB.GetTriangle( nTriaNum, trTriaB) ;
            if ( ( ptFirstV - trTriaB.GetP(0)) * trTriaB.GetN() < - EPS_SMALL)
               nInOutNum = 1 ;
            else if ( ( ptFirstV - trTriaB.GetP(0)) * trTriaB.GetN() > EPS_SMALL)
               nInOutNum = - 1 ;
         }
         if ( nInOutNum == 0) {
            nCurVert = GetNextVertex( nVertNum, ptFirstV) ;
            ++ nVertNum ;
         }
      }
      for ( int nTA = 0 ; nTA < nTriaNumA ; ++ nTA) {
         m_vTria[nTA].nTempPart = nInOutNum ;
      }
      nVertNum = 0 ;
      ptFirstV ;
      nCurVert = SurfB.GetFirstVertex( ptFirstV) ;
      nInOutNum = 0 ;
      while ( nInOutNum == 0  &&  nCurVert != SVT_NULL) {
         int nTriaNum = - 1 ;
         double dMinDist = DBL_MAX ;
         for ( int nTA = 0 ; nTA < nTriaNumA ; ++ nTA) {
            // Se il triangolo A non è valido, continuo
            Triangle3d trTriaA ;
            if ( ! ( GetTriangle( nTA, trTriaA)  &&  trTriaA.Validate( true)))
               continue ;
            DistPointTriangle DistCalculator( ptFirstV, trTriaA) ;
            double dDist ;
            DistCalculator.GetDist( dDist) ;
            if ( dDist < dMinDist) {
               nTriaNum = nTA ;
               dMinDist = dDist ;
            }
         }
         if ( nTriaNum >= 0) {
            Triangle3d trTriaA ;
            GetTriangle( nTriaNum, trTriaA) ;
            if ( ( ptFirstV - trTriaA.GetP( 0)) * trTriaA.GetN() < - EPS_SMALL) {
               nInOutNum = 1 ;
            }
            else if ( ( ptFirstV - trTriaA.GetP(0)) * trTriaA.GetN() > EPS_SMALL) {
               nInOutNum = - 1 ;
            }
         }
         if ( nInOutNum == 0) {
            nCurVert = SurfB.GetNextVertex( nVertNum, ptFirstV) ;
            ++ nVertNum ;
         }
      }
      for ( int nTB = 0 ; nTB < nTriaNumB ; ++ nTB) {
         SurfB.m_vTria[nTB].nTempPart = nInOutNum ;
      }
   }

  // Se c'è stata una ritriangolazione di almeno un triangolo, NON siamo nel caso di tutto dentro o tutto fuori.
  // Studio i triangoli ambigui.
   if ( bRetriangulated) {
      AmbiguosTriangleManager( AmbiguosA, *this) ;
      AmbiguosTriangleManager( AmbiguosB, SurfB) ;
   }

   bool bContinue = true ;
   // Se avvenuta modifica, aggiorno tutto
   if ( bModif)
      bContinue = ( AdjustVertices()  &&  DoCompacting()  &&  SurfB.AdjustVertices()  &&  SurfB.DoCompacting()) ;
  // Triangoli sovrapposti
   if ( bContinue) {
      int nTriaNumA = GetTriangleSize() ;
      // Resetto e ricalcolo la HashGrid della superficie B
      SurfB.ResetHashGrids3d() ;
      for ( int nTA = 0 ; nTA < nTriaNumA ; ++ nTA) {
        // Se il triangolo A non è valido, continuo
         Triangle3d trTriaA ;
         if ( ! GetTriangle( nTA, trTriaA)  ||  ! trTriaA.Validate( true))
            continue ;
        // Box del triangolo A
         BBox3d b3dTriaA ;
         trTriaA.GetLocalBBox( b3dTriaA) ;
         // Recupero i triangoli di B che interferiscono col box del triangolo di A
         INTVECTOR vNearTria ;
         SurfB.GetAllTriaOverlapBox( b3dTriaA, vNearTria) ;
         bool bNewTriaA = true ;
         for ( int nTB = 0 ; nTB < int( vNearTria.size()) ; ++ nTB) {
            // Se il triangolo B non è valido, continuo
            Triangle3d trTriaB ;
            if ( ! SurfB.GetTriangle( vNearTria[nTB], trTriaB)  ||  ! trTriaB.Validate( true))
               continue ;
            // Se i triangoli sono sovrapposti
            TRIA3DVECTOR vTriaAB ;
            Point3d ptTempA, ptTempB ;
            int nIntTypeAB = IntersTriaTria( trTriaA, trTriaB, ptTempA, ptTempB, vTriaAB) ;
            if ( nIntTypeAB == ITTTS_OVERLAPS) {
               bool bInvertB = trTriaA.GetN() * trTriaB.GetN() < 0. ;
               m_vTria[nTA].nTempPart = ( bInvertB ? -2 : 2) ;
               SurfB.m_vTria[vNearTria[nTB]].nTempPart = ( bInvertB ? - 2 : 2) ;
            }
         }
      }
      return ( AdjustVertices()  &&  DoCompacting()  &&  SurfB.AdjustVertices()  &&  SurfB.DoCompacting()) ;
   }
   
   return true ;
}

//----------------------------------------------------------------------------
bool
SurfTriMesh::IdentifyParts( void) const
{
   for ( int i = 0 ; i < int( m_vTria.size()) ; ++ i) {
     // salto triangoli cancellati o già assegnati
      if ( m_vTria[i].nIdVert[0] == SVT_DEL  ||
           abs( m_vTria[i].nTempPart) != 1)
         continue ;
     // set di triangoli da aggiornare
      set<int> stTria ;
      stTria.insert( i) ;
      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].nTempPart == 0) {
               m_vTria[nAdjT].nTempPart = m_vTria[nT].nTempPart ;
               stTria.insert( nAdjT) ;
            }
         }
      }
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
SurfTriMesh::RemoveTJunctions(void)
{
   //PerformanceCounter Counter ;

  // Vettore di indici dei vertici sui lati del triangolo corrente
   INTMATRIX vvIndexMatrix( m_vTria.size()) ;
  // Ciclo sui triangoli della superficie
   for ( int nT = 0 ; nT < int( m_vTria.size()) ; ++ nT) {
     // Riga di matrice per il triangolo 
      INTVECTOR& vIndexRow = vvIndexMatrix[nT] ;
     // Se il triangolo non è valido, passo al successivo
      Triangle3d trTria ;
      if ( ! GetTriangle( nT, trTria)  ||  ! trTria.Validate( true))
         continue ;
     // 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 alla riga della matrice il vertice iniziale del lato
         vIndexRow.emplace_back( m_vTria[nT].nIdVert[nSeg]) ;
         int nPrevSize = int( vIndexRow.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 ;
               GetVertex( m_vTria[vNearTria[nI]].nIdVert[nVert], ptVert) ;
               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) {
                  vIndexRow.emplace_back( m_vTria[vNearTria[nI]].nIdVert[nVert]) ;                
               }
            }
         }
        // Riordino i vertici sul segmento
         auto SortVerteces = [ this, &ptSegSt, &vtSeg]( const int nV, const int nVV)
                                 {  Point3d ptPV, ptPVV ;
                                    GetVertex( nV, ptPV) ;
                                    GetVertex( nVV, ptPVV) ;
                                    return ( ( ptPV - ptSegSt) * vtSeg < ( ptPVV - ptSegSt) * vtSeg) ;
                                 } ;
         sort( vIndexRow.begin() + nPrevSize, vIndexRow.end(), SortVerteces) ;
      }
   }

   //double dTime1 = Counter.Stop() ; Counter.Start() ;  

   for ( int nT = 0 ; nT < int( vvIndexMatrix.size()) ; ++ nT) {
      if ( vvIndexMatrix[nT].size() > 3) {
         PolyLine Polygon ;
         for ( int nV = 0 ; nV < int( vvIndexMatrix[nT].size()) ; ++ nV) {
            Point3d ptPt;
            GetVertex( vvIndexMatrix[nT][nV], ptPt) ;
            Polygon.AddUPoint( 0., ptPt) ;
         }
         Polygon.Close() ;
         PNTVECTOR vPt ;
         INTVECTOR vTr ;
         if ( Triangulate().Make( Polygon, vPt, vTr)) {
           // Rimuovo il triangolo originale
            RemoveTriangle( nT) ;
           // Inserisco i nuovi triangoli che lo sostituiscono
            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) ;
            }
         }
      }
   }

   //double dTime2 = Counter.Stop() ; Counter.Start() ;  

   //string sOut = "T1=" + ToString( dTime1, 1) +
   //              " T2=" + ToString( dTime2, 1) ;
   //LOG_INFO( GetEGkLogger(), sOut.c_str()) ;

   return true ;
}

//----------------------------------------------------------------------------
bool
SurfTriMesh::RemoveCaps(void)
{
   INTVECTOR vIndex ;
   int nTriaNum = GetTriangleSize() ;
   for ( int nT = 0 ; nT < nTriaNum ; ++ nT) {
      Triangle3d trTria ;
      GetTriangle( nT, trTria) ;
      int nMinIndex = -1 ;
      double dMinSqDist = DBL_MAX ;
      for ( int nV = 0 ; nV < 3 ; ++ nV) {         
         Point3d ptVert, ptSegSt, ptSegEn ;
         ptVert = trTria.GetP( nV) ;
         ptSegSt = trTria.GetP( ( nV + 1) % 3) ;
         ptSegEn = trTria.GetP( ( nV + 2) % 3) ;
         CurveLine cvSeg ;
         cvSeg.Set( ptSegSt, ptSegEn) ;
         DistPointLine DistCalculator( ptVert, cvSeg, false) ;
         double dSqDist ;
         DistCalculator.GetSqDist( dSqDist) ;
         Vector3d vtSeg = ptSegEn - ptSegSt ;
         double dSegLen = vtSeg.Len() ;
         vtSeg /= dSegLen ;
         double dProj = ( ptVert - ptSegSt) * vtSeg ;
         if ( dSqDist < dMinSqDist  &&  dProj > EPS_SMALL  &&  dProj < dSegLen - EPS_SMALL) {
            dMinSqDist = dSqDist ;
            nMinIndex = nV ;
         }
      }
      if ( dMinSqDist < EPS_SMALL) {
         vIndex.emplace_back( nT) ;
         vIndex.emplace_back( nMinIndex) ;
      }
   }
   for ( int nPos = 0 ; nPos < int( vIndex.size()) ; nPos += 2) {
      int nT = vIndex[nPos] ;
      int nV = vIndex[nPos + 1] ;
      int nAdjT = m_vTria[nT].nIdAdjac[( nV + 1) % 3] ;
      if ( nAdjT < 0)
         continue ;
      int nAdjV ; 
      for ( nAdjV = 0 ; nAdjV < 3 ; ++ nAdjV) {
         if ( m_vTria[nAdjT].nIdAdjac[nAdjV] == nT)
            break ;
      }
      Triangle3d trTriaAdj ;
      GetTriangle( nAdjT, trTriaAdj) ;
      Vector3d vtLen = trTriaAdj.GetP( ( nAdjV + 1) % 3) - trTriaAdj.GetP( nAdjV) ;
      if ( vtLen * vtLen > 4 * EPS_SMALL * EPS_SMALL) {
         FlipTriangles( nT, nAdjT) ;        
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
SurfTriMesh::RemoveTripleTriangles()
{
   int nCounter = 1 ;
   const int N_STOP = 100 ;
   bool bContinue = true ;
   while ( bContinue  &&  nCounter <= N_STOP) {
      ScanForTripleTriangles( bContinue) ;
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
SurfTriMesh::ScanForTripleTriangles( bool& bModified)
{
  // Setto il flag come nessuna modifica eseguita
   bModified = false ;
  // Ciclo sui vertici
   int nNumVert = GetVertexCount() ;
   for ( int nV = 0 ; nV < nNumVert ; ++ nV) {
      bool bCirc ;
      INTVECTOR vTriaIndex ;
      int nTriaNum = GetAllTriaAroundVertex( nV, vTriaIndex, bCirc) ;
     // Se il vertice è condiviso da tre triangoli tutti adiacenti fra loro a due a due
      if ( nTriaNum == 3  &&  bCirc) {
        // Valuto parallelismo fra le normali
         int nParallel = 0 ;
         for ( int nT1 = 0 ; nT1 < 2 ; ++ nT1) {
            for ( int nT2 = nT1 + 1 ; nT2 < 3 ; ++ nT2) {
               if ( AreSameVectorExact( m_vTria[vTriaIndex[nT1]].vtN, m_vTria[vTriaIndex[nT2]].vtN)) {
                  ++ nParallel ;
               }
            }
         }
        // Se sono tutte parallele fra loro unisco i tre triangoli in uno
         if ( nParallel == 3) {
           // Cerco gli altri triangoli adiacenti a quelli che sostituirò con il triangolo grande.
           // Cerco anche quali lati dei tre triangoli che saranno eliminati sono adiacenti a quelli esterni
            INTVECTOR vVertAndAdjTria ;
            for ( int nTria = 0 ; nTria < 3 ; ++ nTria) {
               for ( int nVert = 0 ; nVert < 3 ; ++ nVert) {
                  if ( m_vTria[vTriaIndex[nTria]].nIdVert[nVert] == nV) {
                     vVertAndAdjTria.emplace_back( m_vTria[vTriaIndex[nTria]].nIdVert[(nVert + 1) % 3]) ;
                     vVertAndAdjTria.emplace_back( m_vTria[vTriaIndex[nTria]].nIdAdjac[(nVert + 1) % 3]) ;
                     break ;
                  }
               }
            }
           // Cerco fra questi altri triangoli esterni sono adiacenti ai triangoli che eliminerò
            int nAdjTriaContactEdge[3] ;
            for ( int nAdjTria = 1 ; nAdjTria < 6 ; nAdjTria += 2) {
               for ( int nAdjEdge = 0 ; nAdjEdge < 3 ; ++ nAdjEdge) {
                  if ( vVertAndAdjTria[nAdjTria] >= 0  &&  m_vTria[vVertAndAdjTria[nAdjTria]].nIdAdjac[nAdjEdge] == vTriaIndex[( nAdjTria - 1) / 2])
                     nAdjTriaContactEdge[(nAdjTria - 1) / 2] = nAdjEdge ;
               }
            }
           // Elimino i triangoli
            for ( int nTria = 0 ; nTria < 3 ; ++ nTria) {
               RemoveTriangle( vTriaIndex[nTria]) ;
            }
           // Aggiungo il novo triangolo
            int nNewVert[3] = { vVertAndAdjTria[0], vVertAndAdjTria[2], vVertAndAdjTria[4] } ;                       
            int nNewTriaId = AddTriangle( nNewVert) ;
            if ( nNewTriaId != SVT_NULL) {
              // Sistemo le adiacenze
               m_vTria[nNewTriaId].nIdAdjac[0] = vVertAndAdjTria[1] ;
               m_vTria[nNewTriaId].nIdAdjac[1] = vVertAndAdjTria[3] ;
               m_vTria[nNewTriaId].nIdAdjac[2] = vVertAndAdjTria[5] ;
               if ( vVertAndAdjTria[1] >= 0)
                  m_vTria[vVertAndAdjTria[1]].nIdAdjac[nAdjTriaContactEdge[0]] = nNewTriaId ;
               if ( vVertAndAdjTria[3] >= 0)
                  m_vTria[vVertAndAdjTria[3]].nIdAdjac[nAdjTriaContactEdge[1]] = nNewTriaId ;
               if ( vVertAndAdjTria[5] >= 0)
                  m_vTria[vVertAndAdjTria[5]].nIdAdjac[nAdjTriaContactEdge[2]] = nNewTriaId ;
              // Setto il flag come modifica eseguita
               bModified = true ;
            }
         }
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
SurfTriMesh::FlipTriangles( int nTA, int nTB)
{
  // Verifico esistenza triangoli
   int nNumTria = GetTriangleCount() ;
   if ( nTA < 0  ||  nTA >= nNumTria  ||  nTB < 0  ||  nTB >= nNumTria)
      return false ;
  // Verifico adiacenza triangoli
   int nEdgeA ;
   for ( nEdgeA = 0 ; nEdgeA < 3 ; ++ nEdgeA) {
      if ( m_vTria[nTA].nIdAdjac[nEdgeA] == nTB)
         break ;
   }
   int nEdgeB ;
   for ( nEdgeB = 0 ; nEdgeB < 3 ; ++ nEdgeB) {
      if ( m_vTria[nTB].nIdAdjac[nEdgeB] == nTA)
         break ;
   }
  // Se uno è adiacente all'altro ma non viceversa, c'è un errore
   if ( ( nEdgeA == 3  &&  nEdgeB < 3)  ||  ( nEdgeA < 3  &&  nEdgeB == 3))
      return false ;
  // Se non sono adiacenti, ho finito
   else if ( nEdgeA == 3  &&  nEdgeB == 3)
      return true ;
  // Se non trovo i vertici del triangolo A, c'è un errorre
   Point3d ptSegSt, ptSegEn, ptVertA ;  
   if ( ! GetVertex( m_vTria[nTA].nIdVert[nEdgeA], ptSegSt)  ||
        ! GetVertex( m_vTria[nTA].nIdVert[( nEdgeA + 1) % 3], ptSegEn)  ||
        ! GetVertex( m_vTria[nTA].nIdVert[( nEdgeA + 2) % 3], ptVertA))
      return false ;
  // Se non trovo i vertici del triangolo B, c'è un errorre
   Point3d ptVertB ;
   if ( ! GetVertex( m_vTria[nTB].nIdVert[( nEdgeB + 2) % 3], ptVertB))
      return false ;
  // Calcolo le proiezioni dei Vertici fuori dal segmento su quest'ultimo
   Vector3d vtVec = ptSegEn - ptSegSt ;
   double dLen = vtVec.Len() ;
   vtVec /= dLen ; 
   double dProjA = ( ptVertA - ptSegSt) * vtVec ;
   double dProjB = ( ptVertB - ptSegSt) * vtVec ;
  // Se le proiezioni sono interne al segmento comune eseguo il flipping
   if ( dProjA > EPS_SMALL  &&  dProjA < dLen - EPS_SMALL  &&
        dProjB > EPS_SMALL  &&  dProjB < dLen - EPS_SMALL) {
      m_vTria[nTA].nIdVert[nEdgeA] = m_vTria[nTB].nIdVert[( nEdgeB + 2) % 3] ;
      m_vTria[nTB].nIdVert[nEdgeB] = m_vTria[nTA].nIdVert[( nEdgeA + 2) % 3] ;
      m_vTria[nTA].nIdAdjac[nEdgeA] = m_vTria[nTB].nIdAdjac[( nEdgeB + 2) % 3] ;
      m_vTria[nTA].nIdAdjac[( nEdgeA + 2) % 3] = nTB ;
      m_vTria[nTB].nIdAdjac[nEdgeB] = m_vTria[nTA].nIdAdjac[( nEdgeA + 2) % 3] ;
      m_vTria[nTB].nIdAdjac[( nEdgeB + 2) % 3] = nTA ;
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
SurfTriMesh::Add( const ISurfTriMesh& Other)
{
   SurfTriMesh SurfB ;
   SurfB.CopyFrom( &Other) ;
   IntersectTriMeshTriangle( SurfB) ;
   IdentifyParts() ;
   SurfB.IdentifyParts() ;
   int nTriaNumA = GetTriangleSize() ;
   for ( int nTA = 0 ; nTA < nTriaNumA ; ++ nTA) {
      if ( m_vTria[nTA].nTempPart == 1  ||  m_vTria[nTA].nTempPart == - 2)
         RemoveTriangle( nTA) ;
   }
   int nTriaNumB = SurfB.GetTriangleSize() ;
   for ( int nTB = 0 ; nTB < nTriaNumB ; ++ nTB) {
      if ( SurfB.m_vTria[nTB].nTempPart == - 1) {
         int nNewVert[3] ;
         for ( int nV = 0 ; nV < 3 ; ++ nV) {
            nNewVert[nV] = AddVertex( SurfB.m_vVert[SurfB.m_vTria[nTB].nIdVert[nV]].ptP) ;
         }
         AddTriangle( nNewVert) ;
      }
   }
   if ( ! AdjustVertices()  ||  ! DoCompacting())
      return false ;
   RemoveTripleTriangles() ;
   if ( ! AdjustVertices()  ||  ! DoCompacting())
      return false ;
   RemoveTJunctions() ;
   return ( AdjustVertices()  &&  DoCompacting()) ;
}

//----------------------------------------------------------------------------
bool
SurfTriMesh::Intersect( const ISurfTriMesh& Other)
{
   SurfTriMesh SurfB ;
   SurfB.CopyFrom( &Other) ;
   IntersectTriMeshTriangle( SurfB) ;
   IdentifyParts() ;
   SurfB.IdentifyParts() ;
   int nTriaNumA = GetTriangleSize() ;
   for ( int nTA = 0 ; nTA < nTriaNumA ; ++ nTA) {
      if ( m_vTria[nTA].nTempPart == - 1  ||  m_vTria[nTA].nTempPart == - 2)
         RemoveTriangle( nTA) ;
   }
   int nTriaNumB = SurfB.GetTriangleSize() ;
   for ( int nTB = 0 ; nTB < nTriaNumB ; ++ nTB) {
      if ( SurfB.m_vTria[nTB].nTempPart == 1) {
         int nNewVert[3] ;
         for ( int nV = 0 ; nV < 3 ; ++ nV) {
            nNewVert[nV] = AddVertex( SurfB.m_vVert[SurfB.m_vTria[nTB].nIdVert[nV]].ptP) ;
         }
         AddTriangle( nNewVert) ;
      }
   }
   if ( ! AdjustVertices()  ||  ! DoCompacting())
      return false ;
   RemoveTripleTriangles() ;
   if ( ! AdjustVertices()  ||  ! DoCompacting())
      return false ;
   RemoveTJunctions() ;
   return ( AdjustVertices()  &&  DoCompacting()) ;
}

//----------------------------------------------------------------------------
bool
SurfTriMesh::Subtract( const ISurfTriMesh& Other)
{
   SurfTriMesh SurfB ;
   SurfB.CopyFrom( &Other) ;
   IntersectTriMeshTriangle( SurfB) ;
   IdentifyParts() ;
   SurfB.IdentifyParts() ;
   int nTriaNumA = GetTriangleSize() ;
   for ( int nTA = 0 ; nTA < nTriaNumA ; ++ nTA) { 
      if  ( m_vTria[nTA].nTempPart == 1  ||  m_vTria[nTA].nTempPart == 2)
         RemoveTriangle( nTA) ;
   }
   int nTriaNumB = SurfB.GetTriangleSize() ;
   for ( int nTB = 0 ; nTB < nTriaNumB ; ++ nTB) {
      if ( SurfB.m_vTria[nTB].nTempPart == 1) {
         int nNewVert[3] ;
         for ( int nV = 0 ; nV < 3 ; ++ nV) {
            nNewVert[nV] = AddVertex( SurfB.m_vVert[SurfB.m_vTria[nTB].nIdVert[nV]].ptP) ;
         }
         swap( nNewVert[1], nNewVert[2]) ;
         AddTriangle( nNewVert) ;
      }
   }
   if ( ! AdjustVertices()  ||  ! DoCompacting())
      return false ;
   RemoveTripleTriangles() ;
   if ( ! AdjustVertices()  ||  ! DoCompacting())
      return false ;
   RemoveTJunctions() ;
   return ( AdjustVertices()  &&  DoCompacting()) ;
}

//----------------------------------------------------------------------------
bool
SurfTriMesh::GetSurfClassification( const ISurfTriMesh& ClassifierSurf,
                                    INTVECTOR& vTriaIn, INTVECTOR& vTriaOut, INTVECTOR& vTriaOnP, INTVECTOR& vTriaOnM, INTVECTOR& vTriaIndef)
{
   // Le superfici devono essere valide
   if ( ! IsValid()  ||  ! ClassifierSurf.IsValid())
      return false ;
   SurfTriMesh SurfC ;
   SurfC.CopyFrom( &ClassifierSurf) ;
   IntersectTriMeshTriangle( SurfC) ;
   IdentifyParts() ;
  
   int nTriaNum = GetTriangleSize() ;
   for ( int nT = 0 ; nT < nTriaNum ; ++ nT) {
      if ( m_vTria[nT].nIdVert[0] == SVT_DEL)
         continue ;
      switch ( m_vTria[nT].nTempPart) {
      case -2 : 
         vTriaOnM.push_back( nT) ;
         break ;
      case -1 :
         vTriaOut.push_back( nT) ;
         break ;
      case 0 :
         vTriaIndef.push_back( nT) ;
         break ;
      case 1 :
         vTriaIn.push_back( nT) ;
         break ;
      case 2 :
         vTriaOnP.push_back( nT) ;
         break ;
      }
   }
   return true ;
}

///////////////////////// DEBUG /////////////////////////////////////////////////////////////////////////////////////////////////
//if ( ( AreSamePointApprox(trTriaA.GetP(0), Point3d(20, 40, 46)) &&
//       AreSamePointApprox(trTriaA.GetP(1), Point3d(100, 40, 30)) &&
//       AreSamePointApprox(trTriaA.GetP(2), Point3d(20, 40, 30))) ||
//     ( AreSamePointApprox(trTriaA.GetP(2), Point3d(20, 40, 46)) &&
//       AreSamePointApprox(trTriaA.GetP(0), Point3d(100, 40, 30)) &&
//       AreSamePointApprox(trTriaA.GetP(1), Point3d(20, 40, 30))) ||
//     ( AreSamePointApprox(trTriaA.GetP(1), Point3d(20, 40, 46)) &&
//       AreSamePointApprox(trTriaA.GetP(2), Point3d(100, 40, 30)) &&
//       AreSamePointApprox(trTriaA.GetP(0), Point3d(20, 40, 30)))) {
//   int a = 0 ;        
//}
//
//if ( ( AreSamePointApprox(trTriaB.GetP(0), Point3d(76, 40, 30)) &&
//       AreSamePointApprox(trTriaB.GetP(1), Point3d(60, 40, 50)) &&
//       AreSamePointApprox(trTriaB.GetP(2), Point3d(20, 40, 30))) ||
//     ( AreSamePointApprox(trTriaB.GetP(2), Point3d(76, 40, 30)) &&
//       AreSamePointApprox(trTriaB.GetP(0), Point3d(60, 40, 50)) &&
//       AreSamePointApprox(trTriaB.GetP(1), Point3d(20, 40, 30))) ||
//     ( AreSamePointApprox(trTriaB.GetP(1), Point3d(76, 40, 30)) &&
//       AreSamePointApprox(trTriaB.GetP(2), Point3d(60, 40, 50)) &&
//       AreSamePointApprox(trTriaB.GetP(0), Point3d(20, 40, 30)))) {
//   int b = 0 ;
//}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////
//if ((AreSamePointApprox(trTriaA.GetP(0), Point3d(20, 40, 46)) &&
//   AreSamePointApprox(trTriaA.GetP(1), Point3d(100, 40, 30)) &&
//   AreSamePointApprox(trTriaA.GetP(2), Point3d(20, 40, 30))) ||
//   (AreSamePointApprox(trTriaA.GetP(2), Point3d(20, 40, 46)) &&
//      AreSamePointApprox(trTriaA.GetP(0), Point3d(100, 40, 30)) &&
//      AreSamePointApprox(trTriaA.GetP(1), Point3d(20, 40, 30))) ||
//      (AreSamePointApprox(trTriaA.GetP(1), Point3d(20, 40, 46)) &&
//         AreSamePointApprox(trTriaA.GetP(2), Point3d(100, 40, 30)) &&
//         AreSamePointApprox(trTriaA.GetP(0), Point3d(20, 40, 30)))) {
//
//   if ((AreSamePointApprox(trTriaB.GetP(0), Point3d(76, 40, 30)) &&
//      AreSamePointApprox(trTriaB.GetP(1), Point3d(60, 40, 50)) &&
//      AreSamePointApprox(trTriaB.GetP(2), Point3d(20, 40, 30))) ||
//      (AreSamePointApprox(trTriaB.GetP(2), Point3d(76, 40, 30)) &&
//         AreSamePointApprox(trTriaB.GetP(0), Point3d(60, 40, 50)) &&
//         AreSamePointApprox(trTriaB.GetP(1), Point3d(20, 40, 30))) ||
//         (AreSamePointApprox(trTriaB.GetP(1), Point3d(76, 40, 30)) &&
//            AreSamePointApprox(trTriaB.GetP(2), Point3d(60, 40, 50)) &&
//            AreSamePointApprox(trTriaB.GetP(0), Point3d(20, 40, 30)))) {
//      int b = 0;
//   }
//}
///////////////////////////////////////////////////////////////////////////////////////////////////
//bool
//FacesRetriangulation()
//{
//  // Mi assicuro che la trimesh sia ok
//  // Assicurarsi che la trimesh sia ok ......
//  // Pongo tutti i triangoli come non visitati
//   for ( int nT = 0 ; nT < GetTriangleSize() ; ++ nT) {
//      m_vTria[nT].nTemp = 0 ;
//   }
//  // Ciclo sui triangoli
//   for ( int nT = 0 ; nT < GetTriangleSize() ; ++ nT) {
//     // Se triangolo non visitato
//      if ( m_vTria[nT].nTemp == 0) {
//         m_vTria[nT].nTemp = 1 ;
//         unordered_set<int> TriaIndexSet ;// devi usare stack
//         TriaIndexSet.emplace( nT) ;
//         while ( ! TriaIndexSet.empty()) {
//            
//         }
//      }
//   }
//   return true ;
//}