//----------------------------------------------------------------------------
//  EgalTech 2019-2021
//----------------------------------------------------------------------------
// File : SurfTriMeshBooleans.cpp     Data : 06.11.21     Versione : 2.3k3
// Contenuto : Implementazione delle funzioni booleane per SurfFTrimesh.
//
//
//
// Modifiche : 10.05.19 LM Creazione modulo.
//             01.10.20 LM Aggiunte scalature 1024x.
//             06.11.21 LM Rifacimento GeneralizedCut.
//
//----------------------------------------------------------------------------

#include "stdafx.h"
#include "SurfTriMesh.h"
#include "CurveLine.h"
#include "CurveComposite.h"
#include "SurfFlatRegion.h"
#include "Triangulate.h"
#include "GeoConst.h"
#include "/EgtDev/Include/EgtNumUtils.h"
#include "/EgtDev/Include/EGkCurve.h"
#include "/EgtDev/Include/EGkDistPointLine.h"
#include "/EgtDev/Include/EGkDistPointCurve.h"
#include "/EgtDev/Include/EGkDistPointTria.h"
#include "/EgtDev/Include/EGkDistPointSurfTm.h"
#include "/EgtDev/Include/EGkIntersLineTria.h"
#include "/EgtDev/Include/EGkIntersLineBox.h"
#include "/EgtDev/Include/EGkIntersPlanePlane.h"
#include "/EgtDev/Include/EGkIntersPlaneTria.h"
#include "/EgtDev/Include/EGkIntersTriaTria.h"
#include "/EgtDev/Include/EGkSfrCreate.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 ;

//----------------------------------------------------------------------------
const double BOOLEAN_SCALE = PREC_SCALE_COEFF ;

//----------------------------------------------------------------------------
bool
SurfTriMesh::DecomposeLoop( CHAINVECTOR& cvOpenChain, INTVECTOR& vnDegVec, PNTMATRIX& cvBoundClosedLoopVec, BOOLVECTOR& vbInOut)
{
  // Valuto se esistono loop non degeneri
   int nExistNotDeg = 0 ;
   for ( int nC = 0 ; nC < int( vnDegVec.size()) ; ++ nC) {
      if ( vnDegVec[nC] > 0)
         ++ nExistNotDeg ;
   }

  // Divido il loop di partenza in sotto-loop
   int nIterationCount = 0 ;
   while ( ! cvOpenChain.empty()) { // per ogni catena aperta...
      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 ;
           // costruisco il primo e il secondo loop che si generano
           // NB. il loop esterno è chiuso, mentre la catena aperta
           // è la successione dei tratti lineari dovuti ai triangoli dell'altra superficie
            PNTVECTOR Loop1, Loop2 ;
           // cambio il punto iniziale del loop di bordo se non coincide già con ptStart della catena aperta
            bool bChangedStart = ChangeStart( ptOpenLoopStP, cvBoundClosedLoopVec[nLoop]) ;
           // splitto
            bool bSplitted = SplitAtPoint( ptOpenLoopEnP, cvBoundClosedLoopVec[nLoop], Loop1, Loop2) ;
            if ( ! ( bChangedStart  &&  bSplitted))
               continue ; // la catena aperta non è interna al loop chiuso attuale

           // il loop 2 segue sempre la direzione della catena, il loop 1 ha dentro la catena invertita
           // ( la direzione della catena è determinata dalla normale dei triangoli che la formano;
           // avendo chiamato la chain senza ammettere inversioni, sono curve tutte concordi )
            bLoopSplitted = true ;
           // ricostrusico i due loop mediante concatenazione
            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] ;
            }
            cvBoundClosedLoopVec[nLoop] = Loop1 ;
            cvBoundClosedLoopVec[nLoop + 1] = Loop2 ;
            vbInOut[nLoop] = false ;
            vbInOut[nLoop + 1] = true ;
            ++ 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) {
           // 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  &&  ( vnDegVec[nLastOpenLoopN] == 1  ||  nExistNotDeg == 0)) {
         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 ;

  // itero sulla unorderd_map
   for ( auto it = LoopLines.begin() ; it != LoopLines.end() ; ++ it) {

     // itero sulle catene ( secondo parametro della mappa)
      for ( int nS1 = 0 ; nS1 < int( it->second.size()) - 1 ; ++ nS1) {
        // tolgo i tratti degeneri ad un punto
         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 ;
            }
         }
      }
     // se non resta nulla, allora passo alla mappa successiva
      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)
      auto itS = Ambiguos.find( it->first) ;
      if ( itS != Ambiguos.end())
         Ambiguos.erase( itS) ;

     // CREAZIONE DEI LOOP ---------------------------------------------
     // Chain
      ChainCurves LoopCreator ;
      LoopCreator.Init( false, 2 * EPS_SMALL, int( it->second.size()), 10 * BOOLEAN_SCALE) ;
     // Recupero le curve ( gli estremi dei tratti lineari) 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 = ( ! it->second.empty() ? it->second[0].ptSt : ORIG) ;
      CHAINVECTOR vChain ;
      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 ;
               }
            }
         }
      }

     // Se esistono loop divisi in catene, le unisco
      Triangle3d trTria ;
      Surf.GetTriangle( it->first, trTria) ;
      for ( int nC1 = 0 ; nC1 < nChainCnt - 1 ; ++ nC1) {
         int nFirstChainLastSegPos = int( vChain[nC1].size()) - 1 ;
         bool bFirstChainInside = nFirstChainLastSegPos >= 0  &&
                                  IsPointInsideTriangle( vChain[nC1][0].ptSt, trTria, TriangleType::OPEN)  &&
                                  IsPointInsideTriangle( vChain[nC1][nFirstChainLastSegPos].ptEn, trTria, TriangleType::OPEN) ;
         for ( int nC2 = nC1 + 1 ; nC2 < nChainCnt  &&  bFirstChainInside ; ++ nC2) {
            int nSecondChainLastSegPos = int( vChain[nC2].size()) - 1 ;
            bool bSecondChainInside = nSecondChainLastSegPos >= 0  &&
                                      IsPointInsideTriangle( vChain[nC2][0].ptSt, trTria, TriangleType::OPEN)  &&
                                      IsPointInsideTriangle( vChain[nC2][nSecondChainLastSegPos].ptEn, trTria, TriangleType::OPEN) ;
            nFirstChainLastSegPos = int( vChain[nC1].size()) - 1 ;
            bool bFisrtSecond = AreSamePointEpsilon( vChain[nC1][nFirstChainLastSegPos].ptEn, vChain[nC2][0].ptSt, 10 * EPS_SMALL) ;
            for ( int nSeg = 0 ; nSeg <= nSecondChainLastSegPos  &&  bSecondChainInside  &&  bFisrtSecond ; ++ nSeg) {
               IntSegment CurSeg ;
               CurSeg.ptSt = vChain[nC2][nSeg].ptSt ;
               CurSeg.ptEn = vChain[nC2][nSeg].ptEn ;
               CurSeg.vtOuter = vChain[nC2][nSeg].vtOuter ;
               CurSeg.bDegenerate = vChain[nC2][nSeg].bDegenerate ;
               vChain[nC1].emplace_back( CurSeg) ;
            }
            bool bSecondFirst = AreSamePointEpsilon( vChain[nC1][0].ptSt, vChain[nC2][nSecondChainLastSegPos].ptEn, 10 * EPS_SMALL) ;
            for ( int nSeg = 0 ; nSeg <= nFirstChainLastSegPos  &&  bSecondChainInside  &&  bSecondFirst  &&  ! bFisrtSecond ; ++ nSeg) {
               IntSegment CurSeg ;
               CurSeg.ptSt = vChain[nC1][nSeg].ptSt ;
               CurSeg.ptEn = vChain[nC1][nSeg].ptEn ;
               CurSeg.vtOuter = vChain[nC1][nSeg].vtOuter ;
               CurSeg.bDegenerate = vChain[nC1][nSeg].bDegenerate ;
               vChain[nC2].emplace_back( CurSeg) ;
            }
            if ( bSecondChainInside  &&  bFisrtSecond  &&  nSecondChainLastSegPos >= 0) {
               nFirstChainLastSegPos = int( vChain[nC1].size()) - 1 ;
               bFirstChainInside = nFirstChainLastSegPos >= 0  &&
                                   IsPointInsideTriangle( vChain[nC1][0].ptSt, trTria, TriangleType::OPEN)  &&
                                   IsPointInsideTriangle( vChain[nC1][nFirstChainLastSegPos].ptEn, trTria, TriangleType::OPEN) ;
               vChain.erase( vChain.begin() + nC2) ;
               -- nChainCnt ;
               -- nC2 ;
            }
            if ( bSecondChainInside  &&  bSecondFirst  &&  ! bFisrtSecond  &&  nFirstChainLastSegPos >= 0) {
               vChain.erase( vChain.begin() + nC1) ;
               -- nChainCnt ;
               -- nC2 ;
               nC2 = nC2 == nC1 ? nC2 + 1 : nC2 ;
            }
         }
      }

     // Chiudo i loop costruiti a partire dalle catene
      for ( int nC = 0 ; nC < nChainCnt ; ++ nC) {
         int nChainLastSegPos = int( vChain[nC].size()) - 1 ;
         if ( IsPointInsideTriangle( vChain[nC][0].ptSt, trTria, TriangleType::OPEN)  &&
              IsPointInsideTriangle( vChain[nC][nChainLastSegPos].ptEn, trTria, TriangleType::OPEN)  &&
              AreSamePointEpsilon( vChain[nC][0].ptSt, vChain[nC][nChainLastSegPos].ptEn, 10 * EPS_SMALL)) {
            IntSegment CurSeg ;
            CurSeg.ptSt = vChain[nC][nChainLastSegPos].ptEn ;
            CurSeg.ptEn = vChain[nC][0].ptSt ;
            CurSeg.vtOuter = ( CurSeg.ptEn - CurSeg.ptSt) ^ trTria.GetN() ;
            CurSeg.vtOuter.Normalize() ;
            CurSeg.bDegenerate = ( CurSeg.ptEn - CurSeg.ptSt).Len() > EPS_SMALL ;
            vChain[nC].emplace_back( CurSeg) ;
         }
      }

     // Fra le catene trovate separo le aperte dalle chiuse
      INTVECTOR vnDegVec ;
      CHAINVECTOR cvClosedChain ;
      CHAINVECTOR cvOpenChain ;
      for ( int nL = 0 ; nL < int( vChain.size()) ; ++ nL) {
         bool bChainDegenerate = ( vChain[nL].size() == 1  &&  AreSamePointApprox( vChain[nL][0].ptSt, vChain[nL][0].ptEn)) ;
         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 ( AreSamePointEpsilon( cvOpenChain[nCh1][nSeg1].ptSt, cvOpenChain[nCh2][nSeg2].ptSt, 100 * EPS_ZERO)  &&
                       AreSamePointEpsilon( cvOpenChain[nCh1][nSeg1].ptEn, cvOpenChain[nCh2][nSeg2].ptEn, 100 * EPS_ZERO)  &&
                       AreSameVectorEpsilon( cvOpenChain[nCh1][nSeg1].vtOuter, cvOpenChain[nCh2][nSeg2].vtOuter, 100 * EPS_ZERO)) {
                     ++ 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 ;
            }
         }
      }

     // Gestione tagli piccoli
      if ( cvOpenChain.size() == 1  &&  cvOpenChain[0].size() == 1  &&  cvClosedChain.empty()) {
         if ( AreSamePointEpsilon( cvOpenChain[0][0].ptSt, cvOpenChain[0][0].ptEn, EPS_SMALL)) {
            Plane3d plPlane ;
            plPlane.Set( cvOpenChain[0][0].ptSt, cvOpenChain[0][0].vtOuter) ;
            Point3d ptIntSt, ptIntEn ;
            int nPlaneTriaInt = IntersPlaneTria( plPlane, trTria, ptIntEn, ptIntSt) ;
            if ( nPlaneTriaInt == IPTT_YES) {
               int nOldTriaCol = Surf.m_vTria[it->first].nTFlag ;
               Surf.RemoveTriangle( it->first) ;
               int nDist[3] ;
               for ( int nV = 0 ; nV < 3 ; ++ nV) {
                  double dD = - DistPointPlane( trTria.GetP( nV), plPlane) ;
                  nDist[nV] = ( abs( dD) < EPS_SMALL ? 0 : ( dD > 0 ? 1 : - 1)) ;
               }
               int nAloneVert = ( nDist[0] == nDist[1] ? 2 : ( nDist[0] == nDist[2] ? 1 : 0)) ;
               if ( nDist[nAloneVert] == - 1) {
                  swap( ptIntSt, ptIntEn) ;
               }
               int nNewAloneId[3] = { Surf.AddVertex( trTria.GetP( nAloneVert)), Surf.AddVertex( ptIntSt), Surf.AddVertex( ptIntEn) } ;
               int nNewAloneTriaNum = Surf.AddTriangle( nNewAloneId, nOldTriaCol) ;
               if ( IsValidSvt( nNewAloneTriaNum)) {
                  Surf.m_vTria[nNewAloneTriaNum].nETempFlag[0] = 0 ;
                  Surf.m_vTria[nNewAloneTriaNum].nETempFlag[1] = 0 ;
                  Surf.m_vTria[nNewAloneTriaNum].nETempFlag[2] = 0 ;
                  Surf.m_vTria[nNewAloneTriaNum].nTempShell = nDist[nAloneVert] ;
                  bModif = true ;
               }
               int nNewCoupleId1[3] = { Surf.AddVertex( ptIntSt), Surf.AddVertex( trTria.GetP( ( nAloneVert + 1) % 3)), Surf.AddVertex( trTria.GetP( ( nAloneVert + 2) % 3)) } ;
               int nNewCoupleTriaNum1 = Surf.AddTriangle( nNewCoupleId1, nOldTriaCol) ;
               if ( IsValidSvt( nNewCoupleTriaNum1)) {
                  Surf.m_vTria[nNewCoupleTriaNum1].nETempFlag[0] = 0 ;
                  Surf.m_vTria[nNewCoupleTriaNum1].nETempFlag[1] = 0 ;
                  Surf.m_vTria[nNewCoupleTriaNum1].nETempFlag[2] = 0 ;
                  Surf.m_vTria[nNewCoupleTriaNum1].nTempShell = - nDist[nAloneVert] ;
                  bModif = true ;
               }
               int nNewCoupleId2[3] = { Surf.AddVertex( ptIntSt), Surf.AddVertex( trTria.GetP( ( nAloneVert + 2) % 3)), Surf.AddVertex( ptIntEn) } ;
               int nNewCoupleTriaNum2 = Surf.AddTriangle( nNewCoupleId2, nOldTriaCol) ;
               if ( IsValidSvt( nNewCoupleTriaNum2)) {
                  Surf.m_vTria[nNewCoupleTriaNum2].nETempFlag[0] = 0 ;
                  Surf.m_vTria[nNewCoupleTriaNum2].nETempFlag[1] = 0 ;
                  Surf.m_vTria[nNewCoupleTriaNum2].nETempFlag[2] = 0 ;
                  Surf.m_vTria[nNewCoupleTriaNum2].nTempShell = - nDist[nAloneVert] ;
                  bModif = true ;
               }
               continue ;
            }
         }
      }
     
     // 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)) ;

     // creo una matrice di punti ; ogni riga corrisponde ad una sequenza di punti che indentificano un loop
      PNTMATRIX cvBoundClosedLoopVec ;
      cvBoundClosedLoopVec.emplace_back( cvFirstLoop) ;

     // vettore per indicare se il loop è interno/esterno all'altra superficie
      BOOLVECTOR vbInOut ;
      vbInOut.push_back( true) ;

     // Divido il loop usando le catene
      if ( ! DecomposeLoop( cvOpenChain, vnDegVec, cvBoundClosedLoopVec, vbInOut)) {
         if ( cvBoundClosedLoopVec.size() == 1  &&  cvOpenChain.size() == 2) {
            Point3d ptLink0St = cvOpenChain[0][0].ptSt ;
            Point3d ptLink0En = cvOpenChain[0].back().ptEn ;
            Point3d ptLink1St = cvOpenChain.back()[0].ptSt ;
            Point3d ptLink1En = cvOpenChain.back().back().ptEn ;
            double dDist01 = sqrt( ( ptLink0En - ptLink1St) * ( ptLink0En - ptLink1St)) ;
            double dDist10 = sqrt( ( ptLink1En - ptLink0St) * ( ptLink1En - ptLink0St)) ;
            if ( dDist01 < 2 * EPS_SMALL) {
               IntSegment LinkingSeg ;
               LinkingSeg.ptSt = ptLink0En ;
               LinkingSeg.ptEn = ptLink1St ;
               LinkingSeg.vtOuter = cvOpenChain[0].back().vtOuter ;
               LinkingSeg.bDegenerate = false ;
               cvOpenChain[0].emplace_back( LinkingSeg) ;
               for ( int nLinkI = 0 ; nLinkI < int( cvOpenChain.back().size()) ; ++ nLinkI) {
                  cvOpenChain[0].emplace_back( cvOpenChain.back()[nLinkI]) ;
               }
               cvOpenChain.resize( 1) ;
               int nComplDeg = vnDegVec[0] * vnDegVec[1] ;
               vnDegVec[0] = nComplDeg ;
               vnDegVec.resize( 1) ;
            }
            else if ( dDist10 < 2 * EPS_SMALL) {
               IntSegment LinkingSeg ;
               LinkingSeg.ptSt = cvOpenChain.back().back().ptEn ;
               LinkingSeg.ptEn = cvOpenChain[0].back().ptSt ;
               LinkingSeg.vtOuter = cvOpenChain.back().back().vtOuter ;
               LinkingSeg.bDegenerate = false ;
               cvOpenChain.back().emplace_back( LinkingSeg) ;
               for ( int nLinkI = 0 ; nLinkI < int( cvOpenChain[0].size()) ; ++ nLinkI) {
                  cvOpenChain.back().emplace_back( cvOpenChain[0][nLinkI]) ;
               }
               cvOpenChain.erase( cvOpenChain.begin()) ;
               int nComplDeg = vnDegVec[0] * vnDegVec[1] ;
               vnDegVec[0] = nComplDeg ;
               vnDegVec.resize( 1) ;
            }
            else {
               Surf.m_vTria[it->first].nTempShell = 0 ;
               continue ;
            }
            vbInOut.resize( 1) ;
            vbInOut[0] = true ;
            if ( ! DecomposeLoop( cvOpenChain, vnDegVec, cvBoundClosedLoopVec, vbInOut)) {
               Surf.m_vTria[it->first].nTempShell = 0 ;
               continue ;
            }
         }
         else {
            Surf.m_vTria[it->first].nTempShell = 0 ;
            continue ;
         }
      }
     // Rimuovo il triangolo corrente
      int nOldTriaCol = Surf.m_vTria[it->first].nTFlag ;
      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 nodo corrispondente a un
        // loop è figlio del nodo corrispondente al loop che lo contiene.
         INTVECTOR vInnerLoop ;
         for ( int nCLI = 0 ; nCLI < int( cvClosedChain.size()) ; ++ nCLI) {
            for ( int nPtNum = 0 ; nPtNum < int( cvClosedChain[nCLI].size()) ; ++ nPtNum) {
               Point3d ptLoopStart = cvClosedChain[nCLI][nPtNum].ptSt ;
               double dMinDist = DBL_MAX ;
               Point3d ptMinDist ;
               bool bPointOnSt = false ;
               bool bPointOnEn = false ;
               int nSegNum = 0 ;
               int nSegMin = 0 ;
               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 ptMinDist2 ;
                  DistPointLine DistCalculator( ptLoopStart, ptS, ptE) ;
                  DistCalculator.GetMinDistPoint( ptMinDist2) ;
                  double dMinDistDot = ( ptLoopStart - ptMinDist2) * vtOut ;
                  if ( dMinDistDot < - EPS_SMALL) {
                     vInnerLoop.emplace_back( nCLI) ;
                     break ;
                  }
               }
               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 < - EPS_SMALL) {
                     vInnerLoop.emplace_back( nCLI) ;
                     break ;
                  }
               }
               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 {
                     ptNextS = ptE ;
                     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 < - EPS_SMALL) {
                     vInnerLoop.emplace_back( nCLI) ;
                     break ;
                  }
               }
            }
         }
        // Verifico se i loop interni sono validi
         bool bAllInvalid = true ;
         for ( int nInnLoop = 0 ; nInnLoop < int( vInnerLoop.size()) ; ++ nInnLoop) {
           // se chain formata da tre segmenti significa che si tratta di due linee sovrapposte ( il terzo tratto è quello aggiunto
           // per forzare la chiusura)
            if ( cvClosedChain[vInnerLoop[nInnLoop]].size() > 3) {
               bAllInvalid = false ;
               break ;
            }
         }

         if ( vInnerLoop.empty()  ||  bAllInvalid) {
           // 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, nOldTriaCol) ;
                  if ( IsValidSvt( nNewTriaNum)) {
                     Surf.m_vTria[nNewTriaNum].nETempFlag[0] = 0 ;
                     Surf.m_vTria[nNewTriaNum].nETempFlag[1] = 0 ;
                     Surf.m_vTria[nNewTriaNum].nETempFlag[2] = 0 ;
                     if ( vbInOut[nLoop])
                        Surf.m_vTria[nNewTriaNum].nTempShell = 1 ;
                     else
                        Surf.m_vTria[nNewTriaNum].nTempShell = - 1 ;
                     bModif = true ;
                  }
               }
            }
         }
         else {
           // se ho più loop, essi descrivono un poligono di n-lati
            POLYLINEVECTOR vPolygons ;
            vPolygons.emplace_back( vplPolyVec[nLoop]) ;

            for ( int nL = 0 ; nL < int( vInnerLoop.size()) ; ++ nL) {
              // per ognuno di essi, ricavo la PolyLine dai punti
               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) ;
            }

           // controllo la direzione della normale del triangolo con quella del poligono di area maggiore 
           // ( per gestire eventuali inscatolamenti)
            Vector3d vtPoly = V_NULL ;
            double dMaxArea = -1 ;
            for ( int i = 1 ; i < int( vPolygons.size()) ; i ++) {
               Plane3d plPoly ;
               double dAreaPoly ;
               vPolygons[i].IsClosedAndFlat( plPoly, dAreaPoly) ;
               if ( dAreaPoly > dMaxArea) {
                  vtPoly = plPoly.GetVersN() ;
                  dMaxArea = dAreaPoly ;
               }            
            }
            bool bCodirectedNormals = trTria.GetN() * vtPoly > 0. ;

            PNTVECTOR vPt ;
            INTVECTOR vTr ;
           // classifico i loop e ricavo la triangolazione di essi mediante classificazione
            if ( Triangulate().MakeAdvanced( 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, nOldTriaCol) ;
                  if ( IsValidSvt( nNewTriaNum)) {
                     Surf.m_vTria[nNewTriaNum].nETempFlag[0] = 0 ;
                     Surf.m_vTria[nNewTriaNum].nETempFlag[1] = 0 ;
                     Surf.m_vTria[nNewTriaNum].nETempFlag[2] = 0 ;
                     if ( bCodirectedNormals)
                        Surf.m_vTria[nNewTriaNum].nTempShell = -1 ;
                     else
                        Surf.m_vTria[nNewTriaNum].nTempShell = 1 ;
                     bModif = true ;
                  }
               }
            }

           // Divido i loop che si autointercettano
           // TO DO : da verificare. Questa porzione di codice non dovrebbe andare prima della triangolazione?
            int nInitialLoopNum = int( vPolygons.size()) ;
            for ( int nL = 1 ; nL < nInitialLoopNum ; ++ nL) {
              // Lista dei punti della PolyLine Loop corrente
               PNTULIST& LoopPointList = vPolygons[nL].GetUPointList() ;
              // Ciclo sui segmenti
               auto itSt2 = LoopPointList.begin() ;
               auto itEn2 = itSt2 ;
               ++ itEn2 ;
               for ( ; itSt2 != LoopPointList.end()  &&  itEn2 != LoopPointList.end() ; ++ itSt2, ++ itEn2) {
                 // Segmento corrente
                  Point3d ptSt = itSt2->first ;
                  Point3d ptEn = itEn2->first ;
                  Vector3d vtSeg = ptEn - ptSt ;
                  double dSegLen = vtSeg.Len() ;
                  vtSeg /= dSegLen ;
                 // Lista di punti da aggiungere al segmento corrente
                  PNTUVECTOR vAddingPointWithOrder ;
                 // Ciclo su tutti i punti non del segmento corrente
                  auto itP = LoopPointList.begin() ;
                  for ( ; itP != LoopPointList.end() ; ++ itP) {
                     if ( itP != itSt2  &&  itP != itEn2) {
                        Point3d ptP = itP->first ;
                        DistPointLine DistCalculator( ptP, ptSt, ptEn) ;
                        double dDist ;
                        DistCalculator.GetDist( dDist) ;
                        double dLongPos = ( ptP - ptSt) * vtSeg ;
                        if ( dDist < EPS_SMALL  &&  dLongPos > EPS_SMALL  &&  dLongPos < dSegLen - EPS_SMALL) {
                           POINTU NewPointU ;
                           NewPointU.first = ptP ;
                           NewPointU.second = dLongPos ;
                           vAddingPointWithOrder.emplace_back( NewPointU) ;
                        }
                     }
                  }
                 // Riordino i punti interni sul segmento esterno in funzione della distanza dall'origine di esso
                  for ( int nPi = 0 ; nPi < int( vAddingPointWithOrder.size()) - 1 ; ++ nPi) {
                     for ( int nPj = nPi + 1 ; nPj < int( vAddingPointWithOrder.size()) ; ++ nPj) {
                        if ( vAddingPointWithOrder[nPi].second > vAddingPointWithOrder[nPj].second) {
                           swap( vAddingPointWithOrder[nPi], vAddingPointWithOrder[nPj]) ;
                        }
                     }
                  }
                 // Aggiungo i punti al loop esterno
                  for ( int nPi = 0 ; nPi < int( vAddingPointWithOrder.size()) ; ++ nPi) {
                     itSt2 = LoopPointList.emplace( itEn2, vAddingPointWithOrder[nPi]) ;
                  }
               }
            
              // Spezzo i loop autointersecantesi
               POLYLINEVECTOR vAuxPolygons ;
               vAuxPolygons.emplace_back( vPolygons[nL]) ;
            
               bool bSplitted = true ;
               while ( bSplitted) {
                  bSplitted = false ;
                  for ( int nl = 0 ; nl < int( vAuxPolygons.size()) ; ++ nl) {
                     PNTULIST& PntLst = vAuxPolygons[nl].GetUPointList() ;
                     PNTVECTOR vPoint ;
                     for ( auto it2 = PntLst.begin() ; it2 != PntLst.end() ; ++ it2) {
                        vPoint.emplace_back( it2->first) ;
                     }
                     int nStartPt = -1 ;
                     int nEndPt = int( vPoint.size()) ;
                     for ( int ni = 0 ; ni < int( vPoint.size()) - 1 ; ++ ni) {
                        for ( int nj = ni + 1 ; nj < int( vPoint.size()) ; ++ nj) {
                           if ( ( ni != 0  ||  nj != int( vPoint.size()) - 1)  &&
                                AreSamePointApprox( vPoint[ni], vPoint[nj])) {
                              nStartPt = ni ;
                              nEndPt = nj ;
                              break ;
                           }
                        }
                     }
                     if ( nStartPt != -1  &&  nEndPt != int( vPoint.size())) {
                        PolyLine CurAuxLoop1, CurAuxLoop2 ;
                        for ( int ni = 0 ; ni < int( vPoint.size()) ; ++ ni) {
                           if ( ni < nStartPt  ||  ni >= nEndPt)
                              CurAuxLoop1.AddUPoint( 0., vPoint[ni]) ;
                           else
                              CurAuxLoop2.AddUPoint( 0., vPoint[ni]) ;
                        }
                        CurAuxLoop2.AddUPoint( 0., vPoint[nStartPt]) ;
                        vAuxPolygons[nl].Clear() ;
                        Point3d ptP ;
                        bool bContinue = CurAuxLoop1.GetFirstPoint( ptP) ;
                        while ( bContinue) {
                           vAuxPolygons[nl].AddUPoint( 0., ptP) ;
                           bContinue = CurAuxLoop1.GetNextPoint( ptP) ;
                        }
                        vAuxPolygons.emplace_back( CurAuxLoop2) ;
                        bSplitted = true ;
                     }
                  }
               }
               bool bReplaced = false ;
               for ( int nl = 0 ; nl < int( vAuxPolygons.size()) ; ++ nl) {
                  if ( ! bReplaced) {
                     vPolygons[nL].Clear() ;
                     Point3d ptP ;
                     bool bContinue = vAuxPolygons[nl].GetFirstPoint( ptP) ;
                     while ( bContinue) {
                        vPolygons[nL].AddUPoint( 0., ptP) ;
                        bContinue = vAuxPolygons[nl].GetNextPoint( ptP) ;
                     }
                     bReplaced = true ;
                  }
                  else {
                     vPolygons.emplace_back( vAuxPolygons[nl]) ;
                  }
               }
            }

           // tolgo il loop esterno per triangolare solo i loop interni
            vPolygons.erase( vPolygons.begin()) ;
           // tolgo le PolyLine non valide ( rejected )
            for ( int i = 0 ; i < int( vPolygons.size()) ; ) {
               if ( ! vPolygons[i].IsClosed())
                  vPolygons.erase( vPolygons.begin() + i) ;
               else
                  ++ i ;
            }
           // eventuale inversione della prima curva ( che determina il verso della triangolazione) per averla orientata 
           // come il triangolo
            if ( ! vPolygons.empty()) {
               Polygon3d pgPol ;
               pgPol.FromPolyLine( vPolygons[0]) ;
               if ( trTria.GetN() * pgPol.GetVersN() < 0.)
                  vPolygons[0].Invert() ;
            }

            if ( Triangulate().MakeAdvanced( 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, nOldTriaCol) ;
                  if ( IsValidSvt( nNewTriaNum)) {
                     Surf.m_vTria[nNewTriaNum].nETempFlag[0] = 0 ;
                     Surf.m_vTria[nNewTriaNum].nETempFlag[1] = 0 ;
                     Surf.m_vTria[nNewTriaNum].nETempFlag[2] = 0 ;
                     if ( bCodirectedNormals)
                        Surf.m_vTria[nNewTriaNum].nTempShell = 1 ;
                     else
                        Surf.m_vTria[nNewTriaNum].nTempShell = -1 ;
                     bModif = true ;
                  }
               }
            }
         }
         vInnerLoop.resize( 0) ;
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
static bool
FindTriaIncidence( const Triangle3d& trTria, const TRIA3DVECTOR& vOthTriaVec, INTMATRIX& vAdjSegToCurTria)
{
   int nNumFoundContact = 0 ;
   vAdjSegToCurTria.resize( 3) ;
   for ( int nEdge = 0 ; nEdge < 3 ; ++ nEdge) {
      for ( int nOther = 0 ; nOther < int( vOthTriaVec.size()) ; ++ nOther) {
         Triangle3d trOthTria = vOthTriaVec[nOther] ;
         Point3d ptSt = trTria.GetP( nEdge) ;
         Point3d ptEn = trTria.GetP( ( nEdge + 1) % 3) ;
         CurveLine cvEdgeLine ;
         cvEdgeLine.Set( ptSt, ptEn) ;
         for ( int nOthEdge = 0 ; nOthEdge < 3 ; ++ nOthEdge) {
            Point3d ptOthSt = trOthTria.GetP( nOthEdge) ;
            Point3d ptOthEn = trOthTria.GetP( ( nOthEdge + 1) % 3) ;
            CurveLine cvOthEdgeLine ;
            cvOthEdgeLine.Set( ptOthSt, ptOthEn) ;
            DistPointLine DistCalculatorCurOth( 0.5 * ( ptSt + ptEn), cvOthEdgeLine, false) ;
            double dSqDistCurOth ;
            DistCalculatorCurOth.GetSqDist( dSqDistCurOth) ;
            DistPointLine DistCalculatorOthCur( 0.5 * ( ptOthSt + ptOthEn), cvEdgeLine, false) ;
            double dSqDistOthCur ;
            DistCalculatorOthCur.GetSqDist( dSqDistOthCur) ;
            if ( dSqDistCurOth < EPS_SMALL * EPS_SMALL  &&  dSqDistOthCur < EPS_SMALL * EPS_SMALL) {
               vAdjSegToCurTria[nEdge].emplace_back( nOther) ;
               ++ nNumFoundContact ;
            }
         }
      }
   }
   return ( nNumFoundContact == int( vOthTriaVec.size())) ;
}

//----------------------------------------------------------------------------
bool
SurfTriMesh::AmbiguosTriangleManager( TRIA3DVECTORMAP& Ambiguos, SurfTriMesh& Surf)
{   
   for ( auto it = Ambiguos.begin() ; it != Ambiguos.end() ; ++ it) {
     // Se il triangolo ha l'indice diverso da zero vuol dire che oltre a un
     // contatto edge-edge ha avuto dei contatti che lo hanno già classificato.
      if ( Surf.m_vTria[it->first].nTempShell != 0)
         continue ;
     // Recupero il triangolo corrente
      Triangle3d trTria ;
      Surf.GetTriangle( it->first, trTria) ;
      trTria.Validate() ;
     // Vettore dei triangoli i cui edge incidono su quelli del triangolo corrente
      TRIA3DVECTOR& vOthTriaVec = it->second ;
     // Vettore degli indici dei segmenti del triangolo adiacenti agli altri triangoli
      INTMATRIX vAdjSegToCurTria ;
      if ( ! FindTriaIncidence( trTria, vOthTriaVec, vAdjSegToCurTria))
         return false ;
     // Classifico il triangolo in base ai triangoli che incidono sui suoi edge
      int nTriaClassificationByEdges[3] = { 0, 0, 0 } ;
      for ( int nEdge = 0 ; nEdge < 3 ; ++ nEdge) {
         if ( int( vAdjSegToCurTria[nEdge].size()) == 0)
            continue ;
        // Trovo due triangoli che incidono sull'edge corrente e che non sono sulla stessa faccia.
        // Si assume che ci siano solo due facce dell'altra superficie per edge del triangolo
         Triangle3d trOthTria1 = vOthTriaVec[vAdjSegToCurTria[nEdge][0]] ;
         Triangle3d trOthTria2 ;
         bool bFound = false ;
         for ( int nTr = 1 ; nTr < int( vAdjSegToCurTria[nEdge].size()) ; ++ nTr) {
            if ( ! AreSameVectorApprox( trOthTria1.GetN(), vOthTriaVec[vAdjSegToCurTria[nEdge][nTr]].GetN())) {
               trOthTria2 = vOthTriaVec[vAdjSegToCurTria[nEdge][nTr]] ;
               bFound = true ;
               break ;
            }
         }
        // Calcolo il vettore ortogonale all'edge corrente che punta al baricentro del triangolo corrente
         Point3d ptBar = ( trTria.GetP( 0) + trTria.GetP( 1) + trTria.GetP( 2)) / 3. ;
         Vector3d vtBarVec = ptBar - trTria.GetP( nEdge) ;
         Vector3d vtEdgeDir = trTria.GetP( ( nEdge + 1) % 3) - trTria.GetP( nEdge) ;
         vtEdgeDir.Normalize() ;
         vtBarVec -= ( vtBarVec * vtEdgeDir) * vtEdgeDir ;
        // Caso con due facce
         if ( bFound) {
            Point3d ptOthBar1 = ( trOthTria1.GetP( 0) + trOthTria1.GetP( 1) + trOthTria1.GetP( 2)) / 3. ;
            Point3d ptOthBar2 = ( trOthTria2.GetP( 0) + trOthTria2.GetP( 1) + trOthTria2.GetP( 2)) / 3. ;
            Vector3d vtBarBar12 = ptOthBar2 - ptOthBar1 ;
            vtBarBar12.Normalize() ;
           // Caso convesso
            if ( vtBarBar12 * trOthTria1.GetN() < EPS_ZERO) {
               double dDot1 = vtBarVec * trOthTria1.GetN() ;
               double dDot2 = vtBarVec * trOthTria2.GetN() ;
               nTriaClassificationByEdges[nEdge] = dDot1 < 0.  &&  dDot2 < 0. ? 1 : -1 ;             
            }
           // Caso concavo
            else {
               double dDot1 = vtBarVec * trOthTria1.GetN() ;
               double dDot2 = vtBarVec * trOthTria2.GetN() ;
               nTriaClassificationByEdges[nEdge] = dDot1 > 0.  &&  dDot2 > 0. ? -1 : 1 ;
            }
         }
        // Caso con una faccia
         else {
            double dDot1 = vtBarVec * trOthTria1.GetN() ;
            nTriaClassificationByEdges[nEdge] = dDot1 < 0 ? 1 : -1 ;
         }
      }
     // Verifico che le classificazioni siano coerenti
      for ( int i = 0 ; i < 3 ; ++ i) {
         if ( nTriaClassificationByEdges[i] == 0)
            continue ;
         Surf.m_vTria[it->first].nTempShell = nTriaClassificationByEdges[i] ;
         int j ;
         for ( j = i + 1 ; j < 3 ; ++ j) {
            if ( nTriaClassificationByEdges[j] != 0  &&  nTriaClassificationByEdges[i] != nTriaClassificationByEdges[j]) {
               Surf.m_vTria[it->first].nTempShell = 0 ;
               break ;
            }
         }
         if ( j < 3)
            break ;
      }
     // Se la classificazione è coerente segno gli edge di contatto come invalicabili
      Surf.m_vTria[it->first].nETempFlag[0] = int( vAdjSegToCurTria[0].size()) > 0 ? 1 : 0 ;
      Surf.m_vTria[it->first].nETempFlag[1] = int( vAdjSegToCurTria[1].size()) > 0 ? 1 : 0 ;
      Surf.m_vTria[it->first].nETempFlag[2] = int( vAdjSegToCurTria[2].size()) > 0 ? 1 : 0 ;
   }
   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].nTempShell = 0 ;
      m_vTria[nTA].nETempFlag[0] = 0 ;
      m_vTria[nTA].nETempFlag[1] = 0 ;
      m_vTria[nTA].nETempFlag[2] = 0 ;
   }
   for ( int nTB = 0 ; nTB < nTriaNumB ; ++ nTB) {
      SurfB.m_vTria[nTB].nTempShell = 0 ;
      SurfB.m_vTria[nTB].nETempFlag[0] = 0 ;
      SurfB.m_vTria[nTB].nETempFlag[1] = 0 ;
      SurfB.m_vTria[nTB].nETempFlag[2] = 0 ;
   }

  // Resetto e ricalcolo la HashGrid della superficie B
   SurfB.ResetHashGrids3d() ;

  // Scorro i triangoli della superficie A
   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) ;

     // Scorro tutti i triangoli di B che intersecano il box di A
      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
         Point3d ptSegSt, ptSegEn ;
         TRIA3DVECTOR vTria ;
         int nIntType = IntersTriaTria( trTriaA, trTriaB, ptSegSt, ptSegEn, vTria) ;

        // se l'intersezione è identificabile con un segmento...
         if ( nIntType == ITTT_EDGE_EDGE_SEG  ||  nIntType == ITTT_EDGE_INT  ||
              nIntType == ITTT_INT_EDGE  ||  nIntType == ITTT_INT_INT_SEG) {

           // Assegno i dati di intersezione per la superficie A
            IntSegment CurInters ;
            CurInters.ptSt = ptSegSt ;             // punto iniziale del segmento
            CurInters.ptEn = ptSegEn ;             // punto finale del segmento
            CurInters.bDegenerate = false ;        // segmento non degenere
            CurInters.vtOuter = trTriaB.GetN() ;   // perpendicolare alla tangente ( oritentata )
            CurInters.vtOuter -= ( ( CurInters.vtOuter * trTriaA.GetN()) * trTriaA.GetN()) ;
            CurInters.vtOuter.Normalize() ;
           // Salvo intersezione per superficie A
            bool bIntOnEndgeA = false ;
            if ( nIntType != ITTT_EDGE_EDGE_SEG  &&  nIntType != ITTT_EDGE_INT) {
              // controllo se tale segmento non è già contenuto
               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 ;

           // Inverto i dati del segmento per intersezione con superficie B
            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 != ITTT_EDGE_EDGE_SEG  &&  nIntType != ITTT_INT_EDGE) {
              // controllo se tale segmento non è già contenuto
               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 ;

           // caso di Intersezione edge-interno ( per superificie A con B)
            if ( bIntOnEndgeA  &&  ! bIntOnEndgeB) {
              // calcolo la massima distanza e il tratto ad essa associato
               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) {
                 // Cerco qual è il segmento di contatto per dichiararlo come invalicabile
                  int nVA ;
                  for ( nVA = 0 ; nVA < 3 ; ++ nVA) {
                     if ( abs( ( trTriaA.GetP( nVA) - trTriaB.GetP( 0)) * trTriaB.GetN()) < EPS_SMALL  &&
                          abs( ( trTriaA.GetP( ( nVA + 1) % 3) - trTriaB.GetP( 0)) * trTriaB.GetN()) < EPS_SMALL)
                        break ;
                  }
                  m_vTria[nTA].nTempShell = ( ( trTriaA.GetP( nSegMaxDist) - trTriaB.GetP( 0)) * trTriaB.GetN() < - EPS_SMALL ? 1 : - 1) ;
                  if ( nVA >= 0  && nVA <= 2)
                     m_vTria[nTA].nETempFlag[nVA] = m_vTria[nTA].nTempShell ;
               }
            }
           // caso di Intersezione interno-edge ( per superficie A con B)
            else if ( ! bIntOnEndgeA  &&  bIntOnEndgeB) {
              // calcolo la massima distanza e il tratto ad essa associato
               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) {
                 // Cerco qual è il segmento di contatto per dichiararlo come invalicabile
                  int nVB ;
                  for ( nVB = 0 ; nVB < 3 ; ++ nVB) {
                     if ( abs( ( trTriaB.GetP( nVB) - trTriaA.GetP(0)) * trTriaA.GetN()) < EPS_SMALL  &&
                          abs( ( trTriaB.GetP( ( nVB + 1) % 3) - trTriaA.GetP( 0)) * trTriaA.GetN()) < EPS_SMALL)
                        break ;
                  }
                  SurfB.m_vTria[vNearTria[nTB]].nTempShell = ( ( trTriaB.GetP( nSegMaxDist) - trTriaA.GetP( 0)) * trTriaA.GetN() < - EPS_SMALL ? 1 : - 1) ;
                  if ( nVB >= 0  &&  nVB <= 2)
                     SurfB.m_vTria[vNearTria[nTB]].nETempFlag[nVB] = SurfB.m_vTria[vNearTria[nTB]].nTempShell ;
               }
            }
           // Intersezione edge-edge: salvo indice e vettore triangoli
           // Uso i triangoli perchè, se un triangolo fosse cancellato, non potrei accedervi poi usando l'indice.
           // Salvando i triangoli risolvo il problema perchè ai fini dello studio di questi contatti, triangolo
           // e sua ritriangolazione portano al medesimo risultato.
            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) ;
               }
            }
         }
      }
   }

  // 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.
  // non mi basta fare un controllo sulle bbox perché non so come sono orientate le superfici e che potrebbero anche non essere chiuse
   bool bRetriangulated = true ;
   if ( ! bModif  &&  ( int( AmbiguosA.size()) == 0  ||  int( AmbiguosB.size()) == 0)) {
      bRetriangulated = false ;
     // devo assegnare a tutti i triangoli della superficie la medesima proprietà ( definita da nInOutNum)
     // ( -1 -> esterno | 0 -> indefinito | +1 -> interno )
     // devo farlo sia per la SurfA( *this) che per la SurfB
      int nInOutNum = 0 ;
      for ( int v = 0 ; v < int( m_vVert.size())  &&  nInOutNum == 0 ; ++ v) {
         double dDist = 0. ;
         DistPointSurfTm distCalculator( m_vVert[v].ptP, SurfB) ;
         if ( distCalculator.GetDist( dDist)  &&  dDist > EPS_SMALL)
            nInOutNum = ( distCalculator.IsPointOnLeftSide() ? 1 : -1) ;
      }
      for ( int nTA = 0 ; nTA < nTriaNumA ; ++ nTA)
         m_vTria[nTA].nTempShell = nInOutNum ;

      nInOutNum = 0 ;
      for ( int v = 0 ; v < int( SurfB.m_vVert.size())  &&  nInOutNum == 0 ; ++ v) {
         double dDist = 0. ;
         DistPointSurfTm distCalculator( SurfB.m_vVert[v].ptP, *this) ;
         if ( distCalculator.GetDist( dDist)  &&  dDist > EPS_SMALL)
            nInOutNum = ( distCalculator.IsPointOnLeftSide() ? 1 : -1) ;
      }
      for ( int nTB = 0 ; nTB < nTriaNumB ; ++ nTB)
         SurfB.m_vTria[nTB].nTempShell = 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 nTriaNum2A = GetTriangleSize() ;
     // Resetto e ricalcolo la HashGrid della superficie B
      SurfB.ResetHashGrids3d() ;
      for ( int nTA = 0 ; nTA < nTriaNum2A ; ++ 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) ;
         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 sono già stati classificati entrambi come sovrapposti continuo
            if ( abs( m_vTria[nTA].nTempShell) == 2  &&  abs( SurfB.m_vTria[vNearTria[nTB]].nTempShell) == 2)
               continue ;

           // Se i triangoli sono sovrapposti
            TRIA3DVECTOR vTriaAB ;
            Point3d ptTempA, ptTempB ;
            int nIntTypeAB = IntersTriaTria( trTriaA, trTriaB, ptTempA, ptTempB, vTriaAB) ;
            if ( nIntTypeAB == ITTT_OVERLAPS) {
               m_vTria[nTA].nTempShell = 2 ;
               SurfB.m_vTria[vNearTria[nTB]].nTempShell = 2 ;
            }
            else if ( nIntTypeAB == ITTT_COUNTER_OVERLAPS) {
               m_vTria[nTA].nTempShell = -2 ;
               SurfB.m_vTria[vNearTria[nTB]].nTempShell = -2 ;  
            }
         }
      }
      return true ;
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
SurfTriMesh::IdentifyShells( 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].nTempShell) != 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) {
            if ( m_vTria[nT].nETempFlag[j] != 0)
               continue ;
            int nAdjT = m_vTria[nT].nIdAdjac[j] ;
            if ( nAdjT != SVT_NULL  &&  m_vTria[nAdjT].nTempShell == 0) {
               m_vTria[nAdjT].nTempShell = m_vTria[nT].nTempShell ;
               stTria.insert( nAdjT) ;
            }
         }
      }
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
SurfTriMesh::Add( const ISurfTriMesh& Other)
{

  // le superfici devono essere valide
   if ( ! IsValid()  ||  ! Other.IsValid())
      return false ;

  // se la seconda è vuota non devo fare alcunchè
   if ( Other.IsEmpty())
      return true ;

   m_OGrMgr.Clear() ;

  // se la prima è vuota, copio la seconda nella prima
   if ( IsEmpty()) {
      CopyFrom( &Other) ;
      return true ;
   }

  // clono la superficie B
   SurfTriMesh SurfB ;
   SurfB.CopyFrom( &Other) ;

  // creazione del frame per scalare A e B
   Frame3d frScalingRef ;
   frScalingRef.Set( m_vVert[0].ptP, X_AX, Y_AX, Z_AX) ;
  // scalo A e B
   Scale( frScalingRef, BOOLEAN_SCALE, BOOLEAN_SCALE, BOOLEAN_SCALE) ;
   SurfB.Scale( frScalingRef, BOOLEAN_SCALE, BOOLEAN_SCALE, BOOLEAN_SCALE) ;

  // tengo una copia di B ( la superficie B viene modificata durante la ritriangolazione )
   SurfTriMesh SurfB_cl ;
   SurfB_cl.CopyFrom( &SurfB) ;

  // tengo una copia di B ( la superficie B viene modificata durante la ritriangolazione )
   SurfTriMesh SurfA_cl ;
   SurfA_cl.CopyFrom( this) ;

  // ritriangolo le due superfici mediante ogni intersezione Triangolo-Triangolo
   IntersectTriMeshTriangle( SurfB) ;

  // assegno un medesimo indice ai triangoli che non interferiscono con altri
   IdentifyShells() ;
   SurfB.IdentifyShells() ;

  // rimozione dei triangoli di A con proprietà 1 e -2 ( e gestione dei triangoli ambigui, 0)
   int nTriaNumA = GetTriangleSize() ;
   for ( int nTA = 0 ; nTA < nTriaNumA ; ++ nTA) {
      if ( m_vTria[nTA].nTempShell == 1  ||  m_vTria[nTA].nTempShell == - 2)
         RemoveTriangle( nTA) ;
     // se triangolo ambiguo
      if ( m_vTria[nTA].nTempShell == 0) {
         Triangle3d TriaA ;
         GetTriangle( nTA, TriaA) ;
        // rimuovo il triangolo se interno a surfB ( basta controllare un solo vertice(?) )
         DistPointSurfTm distCalculator( TriaA.GetP( 0), SurfB_cl) ;
         if ( distCalculator.IsPointOnLeftSide())
            RemoveTriangle( nTA) ;
      }
   }

  // aggiunta di tutti i triangoli di B con proprietà -1 ( e gestione dei triangoli ambigui, 0)
   int nPrevMaxTFlag = m_nMaxTFlag ;
   int nTriaNumB = SurfB.GetTriangleSize() ;
   for ( int nTB = 0 ; nTB < nTriaNumB ; ++ nTB) {
      bool bAdd = ( SurfB.m_vTria[nTB].nTempShell == - 1) ;
      if ( ! bAdd) {
        // se ambiguo
         if ( SurfB.m_vTria[nTB].nTempShell == 0) {
           // aggiungo se il triangolo è a surfA ( basta controllare un solo vertice(?) )
            DistPointSurfTm distCalculator( SurfB.m_vVert[SurfB.m_vTria[nTB].nIdVert[0]].ptP, SurfA_cl) ;
            bAdd = ! distCalculator.IsPointOnLeftSide() ;
         }
      }
      if ( bAdd) {
         int nNewVert[3] ;
         for ( int nV = 0 ; nV < 3 ; ++ nV)
            nNewVert[nV] = AddVertex( SurfB.m_vVert[SurfB.m_vTria[nTB].nIdVert[nV]].ptP) ;
         if ( nPrevMaxTFlag == m_nMaxTFlag)
            ++ m_nMaxTFlag ;
         AddTriangle( nNewVert, m_nMaxTFlag) ;
      }
   }

  // sistemazioni varie
   bool bOk = ( AdjustVertices()  &&  DoCompacting()) ;
   bool bModified = false ;
   bOk = bOk && RemoveDoubleTriangles( bModified) ;
   if ( bModified)
      bOk = bOk && ( AdjustVertices()  &&  DoCompacting()) ;
   bOk = bOk && RemoveTJunctions( bModified) ;
   if ( bModified)
      bOk = bOk && ( AdjustVertices()  &&  DoCompacting()) ;

  // scalo alla dimensioni originali
   Scale( frScalingRef, 1. / BOOLEAN_SCALE, 1. / BOOLEAN_SCALE, 1. / BOOLEAN_SCALE) ;

  // semplifico eventuale geometria delle facce
   if ( ! SimplifyFacets())
      LOG_ERROR( GetEGkLogger(), "Error in SimplifyFacets of Stm::Add")

   return bOk ;
}

//----------------------------------------------------------------------------
bool
SurfTriMesh::Intersect( const ISurfTriMesh& Other)
{

  // Le superfici devono essere valide
   if ( ! IsValid()  ||  ! Other.IsValid())
      return false ;

  // Se la prima è vuota non devo fare alcunchè
   if ( IsEmpty())
      return true ;

   m_OGrMgr.Clear() ;

  // Se la seconda è vuota, basta vuotare la prima
   if ( Other.IsEmpty()) {
      Clear() ;
      m_bOriented = true ;
      m_bClosed = true ;
      return true ;
   }

  // Surf A è la superficie attuale ( this), SurfB è l'altra
   SurfTriMesh SurfB ;
   SurfB.CopyFrom( &Other) ;

  // scalo la superficie
   Frame3d frScalingRef ;
   frScalingRef.Set( m_vVert[0].ptP, X_AX, Y_AX, Z_AX) ;
   Scale( frScalingRef, BOOLEAN_SCALE, BOOLEAN_SCALE, BOOLEAN_SCALE) ;
   SurfB.Scale( frScalingRef, BOOLEAN_SCALE, BOOLEAN_SCALE, BOOLEAN_SCALE) ;

  // tengo una copia di B ( la superficie B viene modificata durante la ritriangolazione )
   SurfTriMesh SurfB_cl ;
   SurfB_cl.CopyFrom( &SurfB) ;

  // tengo una copia di B ( la superficie B viene modificata durante la ritriangolazione )
   SurfTriMesh SurfA_cl ;
   SurfA_cl.CopyFrom( this) ;

  // ritriangolo le due superfici mediante ogni intersezione Triangolo-Triangolo
   IntersectTriMeshTriangle( SurfB) ;

  // assegno un medesimo indice ai triangoli che non interferiscono con altri
   IdentifyShells() ;
   SurfB.IdentifyShells() ;

  // rimozione dei triangoli di A con proprietà -1 e -2 (e gestione dei triangoli ambigui, 0)
   int nTriaNumA = GetTriangleSize() ;
   for ( int nTA = 0 ; nTA < nTriaNumA ; ++ nTA) {
      if ( m_vTria[nTA].nTempShell == - 1  ||  m_vTria[nTA].nTempShell == - 2)
         RemoveTriangle( nTA) ;
     // se triangolo ambiguo
      else if ( m_vTria[nTA].nTempShell == 0) {
         Triangle3d TriaA ;
         GetTriangle( nTA, TriaA) ;
        // rimuovo il triangolo se fuori a surfB ( basta controllare un solo vertice(?) )
         DistPointSurfTm distCalculator( TriaA.GetP( 0), SurfB_cl) ;
         if ( ! distCalculator.IsPointOnLeftSide())
            RemoveTriangle( nTA) ;
      }
   }

  // aggiunta dei triangoli di B con proprietà +1 ( e gestione dei triangoli ambigui, 0)
   int nPrevMaxTFlag = m_nMaxTFlag ;
   int nTriaNumB = SurfB.GetTriangleSize() ;
   for ( int nTB = 0 ; nTB < nTriaNumB ; ++ nTB) {
      bool bAdd = ( SurfB.m_vTria[nTB].nTempShell == 1) ;
      if ( ! bAdd  &&  SurfB.m_vTria[nTB].nTempShell == 0) {
        // aggiungo se il triangolo è interno a surfA ( basta controllare un solo vertice(?) )
         DistPointSurfTm distCalculator( SurfB.m_vVert[SurfB.m_vTria[nTB].nIdVert[0]].ptP, SurfA_cl) ;
         bAdd = distCalculator.IsPointOnLeftSide() ;
      }
      if ( bAdd) {
         int nNewVert[3] ;
         for ( int nV = 0 ; nV < 3 ; ++ nV)
            nNewVert[nV] = AddVertex( SurfB.m_vVert[SurfB.m_vTria[nTB].nIdVert[nV]].ptP) ;
         if ( nPrevMaxTFlag == m_nMaxTFlag)
            ++ m_nMaxTFlag ;
         AddTriangle( nNewVert, m_nMaxTFlag) ;
      }
   }
  
  // sistemazioni varie
   bool bOk = ( AdjustVertices()  &&  DoCompacting()) ;
   bool bModified = false ;
   bOk = bOk && RemoveDoubleTriangles( bModified) ; 
   if ( bModified)   
      bOk = bOk && ( AdjustVertices()  &&  DoCompacting()) ;
   bOk = bOk && RemoveTJunctions( bModified) ;
   if ( bModified)
      bOk = bOk && ( AdjustVertices()  &&  DoCompacting()) ;

  // scalo alle dimensioni originali
   Scale( frScalingRef, 1. / BOOLEAN_SCALE, 1. / BOOLEAN_SCALE, 1. / BOOLEAN_SCALE) ;

   if ( ! SimplifyFacets())
      LOG_ERROR( GetEGkLogger(), "Error in SimplifyFacets of Stm::Intersect")

   return bOk ;
}

//----------------------------------------------------------------------------
bool
SurfTriMesh::Subtract( const ISurfTriMesh& Other)
{

  // le superfici devono essere valide
   if ( ! IsValid()  ||  ! Other.IsValid())
      return false ;

  // se una delle due è vuota non devo fare alcunchè
   if ( IsEmpty()  ||  Other.IsEmpty())
      return true ;

  // pulisco grafica
   m_OGrMgr.Clear() ;

  // clono superficie B
   SurfTriMesh SurfB ;
   SurfB.CopyFrom( &Other) ;

  // creazione del frame per scalare A e B
   Frame3d frScalingRef ;
   frScalingRef.Set( m_vVert[0].ptP, X_AX, Y_AX, Z_AX) ;
  // scalo A e B
   Scale( frScalingRef, BOOLEAN_SCALE, BOOLEAN_SCALE, BOOLEAN_SCALE) ;
   SurfB.Scale( frScalingRef, BOOLEAN_SCALE, BOOLEAN_SCALE, BOOLEAN_SCALE) ;

  // tengo una copia di B ( la superficie B viene modificata durante la ritriangolazione )
   SurfTriMesh SurfB_cl ;
   SurfB_cl.CopyFrom( &SurfB) ;

  // tengo una copia di B ( la superficie B viene modificata durante la ritriangolazione )
   SurfTriMesh SurfA_cl ;
   SurfA_cl.CopyFrom( this) ;
 
  // ritriangolo le due superfici mediante ogni intersezione Triangolo-Triangolo
   IntersectTriMeshTriangle( SurfB) ;

  // assegno un medesimo indice ai triangoli che non interferiscono con altri
   IdentifyShells() ;
   SurfB.IdentifyShells() ;

  // rimozione dei triangoli di A con proprietà 1 e 2 ( e gestione triangoli ambigui, 0)
   int nTriaNumA = GetTriangleSize() ;
   for ( int nTA = 0 ; nTA < nTriaNumA ; ++ nTA) { 
      if  ( m_vTria[nTA].nTempShell == 1  ||  m_vTria[nTA].nTempShell == 2)
         RemoveTriangle( nTA) ;
     // se triangolo ambiguo...
      if ( m_vTria[nTA].nTempShell == 0) {
         Triangle3d TriaA ;
         GetTriangle( nTA, TriaA) ;
        // rimuovo il triangolo se interno a SurfB ( basta controllare un solo vertice(?) )
         DistPointSurfTm distCalculator( TriaA.GetP( 0), SurfB_cl) ;
         if ( distCalculator.IsPointOnLeftSide())
            RemoveTriangle( nTA) ;
      }
   }

  // aggiunta dei triangoli di B con proprietà +1 ( e gestione triangoli ambigui, 0)
   int nPrevMaxTFlag = m_nMaxTFlag ;
   int nTriaNumB = SurfB.GetTriangleSize() ;
   for ( int nTB = 0 ; nTB < nTriaNumB ; ++ nTB) {
      bool bAdd = SurfB.m_vTria[nTB].nTempShell == 1 ;
      if ( ! bAdd) {
        // se ambiguo
         if ( SurfB.m_vTria[nTB].nTempShell == 0) {
           // aggiungo il triangolo se interno alla SurfA ( basta controllare un solo vertice(?) )
            DistPointSurfTm distCalculator( SurfB.m_vVert[SurfB.m_vTria[nTB].nIdVert[0]].ptP, SurfA_cl) ;
            bAdd = distCalculator.IsPointOnLeftSide() ;
         }
      }
      if ( bAdd) {
         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]) ;
         if ( nPrevMaxTFlag == m_nMaxTFlag)
            ++ m_nMaxTFlag ;
         AddTriangle( nNewVert, m_nMaxTFlag) ;
      }
   }

  // sistemazioni varie
   bool bOk = ( AdjustVertices()  &&  DoCompacting()) ;
   bool bModified = false ;
   bOk = bOk && RemoveDoubleTriangles( bModified) ;
   if ( bModified)
      bOk = bOk && ( AdjustVertices()  &&  DoCompacting()) ;
   bOk = bOk && RemoveTJunctions( bModified) ;
   if ( bModified)
      bOk = bOk && ( AdjustVertices()  &&  DoCompacting()) ;

  // scalo alle dimensioni originali
   Scale( frScalingRef, 1. / BOOLEAN_SCALE, 1. / BOOLEAN_SCALE, 1. / BOOLEAN_SCALE) ;

  // semplifico le facce
   if ( ! SimplifyFacets())
      LOG_ERROR( GetEGkLogger(), "Error in SimplifyFacets of Stm::Subtract")

   return bOk ;
}

//----------------------------------------------------------------------------
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 ;
   if ( IsEmpty()  ||  ClassifierSurf.IsEmpty())
      return false ;

   SurfTriMesh SurfC ;
   SurfC.CopyFrom( &ClassifierSurf) ;

   Frame3d frScalingRef ;
   frScalingRef.Set( m_vVert[0].ptP, X_AX, Y_AX, Z_AX) ;
   Scale( frScalingRef, BOOLEAN_SCALE, BOOLEAN_SCALE, BOOLEAN_SCALE) ;
   SurfC.Scale( frScalingRef, BOOLEAN_SCALE, BOOLEAN_SCALE, BOOLEAN_SCALE) ;

   IntersectTriMeshTriangle( SurfC) ;
   IdentifyShells() ;
   Scale( frScalingRef, 1. / BOOLEAN_SCALE, 1. / BOOLEAN_SCALE, 1. / BOOLEAN_SCALE) ;
  
   int nTriaNum = GetTriangleSize() ;
   for ( int nT = 0 ; nT < nTriaNum ; ++ nT) {
      if ( m_vTria[nT].nIdVert[0] == SVT_DEL)
         continue ;
      switch ( m_vTria[nT].nTempShell) {
      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 ;
}

//----------------------------------------------------------------------------
bool
SurfTriMesh::CutWithOtherSurf( const ISurfTriMesh& CutterSurf, bool bInVsOut, bool bSaveOnEq)
{
  // Le superfici devono essere valide
   if ( ! IsValid()  ||  ! CutterSurf.IsValid())
      return false ;

  // Se una delle due è vuota non devo fare alcunchè
   if ( IsEmpty()  ||  CutterSurf.IsEmpty())
      return true ;

   m_OGrMgr.Clear() ;

   SurfTriMesh SurfC ;
   SurfC.CopyFrom( &CutterSurf) ;

   Frame3d frScalingRef ;
   frScalingRef.Set( m_vVert[0].ptP, X_AX, Y_AX, Z_AX) ;
   Scale( frScalingRef, BOOLEAN_SCALE, BOOLEAN_SCALE, BOOLEAN_SCALE) ;
   SurfC.Scale( frScalingRef, BOOLEAN_SCALE, BOOLEAN_SCALE, BOOLEAN_SCALE) ;

   IntersectTriMeshTriangle( SurfC) ;
   IdentifyShells() ;

   int nPartToRemove = ( bInVsOut ? -1 : 1) ;
   int nCoplanarPartToRemove = ( bSaveOnEq ? ( nPartToRemove ? -2 : 2) : 5) ;
   int nTriaNum = GetTriangleSize() ;
   for ( int nT = 0 ; nT < nTriaNum ; ++ nT) {
      if ( m_vTria[nT].nTempShell == nPartToRemove  ||  m_vTria[nT].nTempShell == nCoplanarPartToRemove)
         RemoveTriangle( nT) ;
   }

   bool bOk = ( AdjustVertices()  &&  DoCompacting()) ;

   bool bModified = false ;
   bOk = bOk && RemoveDoubleTriangles( bModified) ;
   if ( bModified)
      bOk = bOk && ( AdjustVertices()  &&  DoCompacting()) ;

   bOk = bOk && RemoveTJunctions( bModified) ;
   if ( bModified)
      bOk = bOk && ( AdjustVertices()  &&  DoCompacting()) ;

   Scale( frScalingRef, 1. / BOOLEAN_SCALE, 1. / BOOLEAN_SCALE, 1. / BOOLEAN_SCALE) ;

   if ( ! SimplifyFacets())
      LOG_ERROR( GetEGkLogger(), "Error in SimplifyFacets of Stm::CutWithOtherSurf")

   return bOk ;
}

//----------------------------------------------------------------------------
bool
SurfTriMesh::Repair( double dMaxEdgeLen)
{
  // La superficie deve essere valida
   if ( ! IsValid())
      return false ;

  // Forzo aggiornamento grafica
   m_OGrMgr.Clear() ;

  // Aggiustamento generale
   if ( ! AdjustVertices()  ||  ! DoCompacting())
      return false ;

  // Elimino triangoli coincidenti
   bool bModified = false ;
   RemoveDoubleTriangles( bModified) ;
   if ( bModified) {
      if ( ! AdjustVertices()  ||  ! DoCompacting())
         return false ;
   }

  // Elimino giunzioni a T, con opportuno inserimento di triangoli
   RemoveTJunctions( bModified) ;
   if ( bModified) {
      if ( ! AdjustVertices()  ||  ! DoCompacting())
         return false ;
   }

  // Ritriangolo le facce
   if ( ! SimplifyFacets( dMaxEdgeLen, true))
      LOG_ERROR( GetEGkLogger(), "Error in SimplifyFacets of Stm::Repair")

   return true ;
}