//----------------------------------------------------------------------------
//  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 "DistPointLine.h"
#include "Triangulate.h"
#include "GeoConst.h"
#include "/EgtDev/Include/EgtNumUtils.h"
#include "/EgtDev/Include/EGkCurve.h"
#include "/EgtDev/Include/EGkDistPointCurve.h"
#include "/EgtDev/Include/EGkDistPointTria.h"
#include "/EgtDev/Include/EGkIntersLineTria.h"
#include "/EgtDev/Include/EGkIntersLineBox.h"
#include "/EgtDev/Include/EGkIntersPlanePlane.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 = 1024 ;

//----------------------------------------------------------------------------
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.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)  ||
                 ( nLastOpenLoopPoint == 0  &&  ( Loop1.size() == 2  || Loop2.size() == 2)))
               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] ;
            }
            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 ;
  // 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) ;
         }
      }

     // Creo i loop
      ChainCurves LoopCreator ;
      LoopCreator.Init( false, 2 * EPS_SMALL, int( it->second.size()), 10 * BOOLEAN_SCALE) ;
     // 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 = ( it->second.size() > 0 ? 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
      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
      if ( ! DecomposeLoop( cvOpenChain, vnDegVec, cvBoundClosedLoopVec, vbInOut)) {
         if ( int( cvBoundClosedLoopVec.size()) == 1  &&  int( 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].nTempPart = 0 ;
               continue ;
            }
            vbInOut.resize( 1) ;
            vbInOut[0] = true ;
            if ( ! DecomposeLoop( cvOpenChain, vnDegVec, cvBoundClosedLoopVec, vbInOut)) {
               Surf.m_vTria[it->first].nTempPart = 0 ;
               continue ;
            }
         }
         else {
            Surf.m_vTria[it->first].nTempPart = 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 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) {
            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 ;
                  }
               }
            }
         }
        // Elimino loop interni non validi
         bool bDouble = true ;
         for ( int nInnLoop = 0 ; nInnLoop < int( vInnerLoop.size()) ; ++ nInnLoop) {
            if ( int( cvClosedChain[vInnerLoop[nInnLoop]].size()) > 2) {
               bDouble = false ;
               break ;
            }
         }
         
         if ( vInnerLoop.size() == 0  ||  bDouble) {
           // 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].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() ;
               }
            }

           // Aggiungo al loop esterno i punti dei loop interni che si trovano su di esso
            PNTULIST& ExternLoopList = vPolygons[0].GetUPointList() ;
           // Ciclo sui segmenti del loop esterno
            auto itSt = ExternLoopList.begin() ;
            auto itEn = itSt ;
            ++ itEn ; 
            for ( ; itSt != ExternLoopList.end()  &&  itEn != ExternLoopList.end() ; ++ itSt, ++ itEn) {
              // Estremi del segmento corrente del loop esterno e scorrispondente vettore
               Point3d ptSt = itSt->first ;
               Point3d ptEn = itEn->first ;
               Vector3d vtSeg = ptEn - ptSt ;
               double dSegLen = vtSeg.Len() ;
               vtSeg /= dSegLen ;
              // Vettore dei punti dei loop interni che stanno sul segmento del loop esterno
               PNTUVECTOR vPointWithOrder ;
              // Ciclo sui loop interni
               for ( int nInnPoly = 1 ; nInnPoly < int( vPolygons.size()) ; ++ nInnPoly) {
                 // Ciclo sui punti dei loop interni
                  Point3d ptInnPoint ;
                  bool bIsFirst = true ;
                  bool bContinue = vPolygons[nInnPoly].GetFirstPoint( ptInnPoint) ;
                  while ( bContinue) {                     
                     DistPointLine DistCalculator( ptInnPoint, ptSt, ptEn) ;
                     double dDist ;
                     DistCalculator.GetDist( dDist) ;
                     double dLongPos = ( ptInnPoint - ptSt) * vtSeg ;
                     if ( dDist < EPS_SMALL  &&  dLongPos > 0.  &&  dLongPos < dSegLen) {
                        POINTU NewPointU ;
                        NewPointU.first = ptInnPoint ;
                        NewPointU.second = dLongPos ;
                        if ( ! bIsFirst)
                           vPointWithOrder.emplace_back( NewPointU) ;
                     }
                     bIsFirst = false ;
                     bContinue = vPolygons[nInnPoly].GetNextPoint( ptInnPoint) ;
                  }
               }
              // Riordino i punti interni sul segmento esterno in funzione della distanza dall'origine di esso
               for ( int nPi = 0 ; nPi < int( vPointWithOrder.size()) - 1 ; ++ nPi) {
                  for ( int nPj = nPi + 1 ; nPj < int( vPointWithOrder.size()) ; ++ nPj) {
                     if ( vPointWithOrder[nPi].second > vPointWithOrder[nPj].second) {
                        swap( vPointWithOrder[nPi], vPointWithOrder[nPj]) ;
                     }
                  }
               }
              // Aggiungo i punti al loop esterno
               for ( int nPi = 0 ; nPi < int( vPointWithOrder.size()) ; ++ nPi) {
                  itSt = ExternLoopList.emplace( itEn, vPointWithOrder[nPi]) ;
               }                      
            }

            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, 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].nTempPart = -1 ;
                     else
                        Surf.m_vTria[nNewTriaNum].nTempPart = 1 ;
                     bModif = true ;
                  }
               }
            }
           // Divido i loop che si autointercettano
            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 > 0.  &&  dLongPos < dSegLen) {
                           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 ( true/*vAuxPolygons.GetAreaXY(double& dArea)*/) {
                     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]) ;
                     }
                  }
               }
            }
            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, 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].nTempPart = 1 ;
                        else
                           Surf.m_vTria[nNewTriaNum].nTempPart = -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].nTempPart != 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].nTempPart = nTriaClassificationByEdges[i] ;
         int j ;
         for ( j = i + 1 ; j < 3 ; ++ j) {
            if ( nTriaClassificationByEdges[j] != 0  &&  nTriaClassificationByEdges[i] != nTriaClassificationByEdges[j]) {
               Surf.m_vTria[it->first].nTempPart = 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].nTempPart = 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].nTempPart = 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() ;
   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) ;
      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) ;
         if ( nIntType == ITTT_EDGE_EDGE_SEG  ||
              nIntType == ITTT_EDGE_INT  ||
              nIntType == ITTT_INT_EDGE  ||
              nIntType == ITTT_INT_INT_SEG) {
           // Assegno i dati di intersezione
            IntSegment CurInters ;
            if ( nIntType == ITTT_EDGE_EDGE_SEG  ||  nIntType == ITTT_EDGE_INT  ||
                 nIntType == ITTT_INT_EDGE  ||  nIntType == ITTT_INT_INT_SEG) {
               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 != ITTT_EDGE_EDGE_SEG  &&  nIntType != ITTT_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 != ITTT_EDGE_EDGE_SEG  &&  nIntType != ITTT_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) {
                 // 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].nTempPart = ( ( 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].nTempPart ;
               }
            }
           // 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) {
                 // 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]].nTempPart = ( ( 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]].nTempPart ;
               }
            }
           // 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.
   bool bRetriangulated = true ;
   if ( ! bModif  &&  ( int( AmbiguosA.size()) == 0  ||  int( AmbiguosB.size()) == 0)) {
      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 ;
      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 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 i triangoli sono sovrapposti
            TRIA3DVECTOR vTriaAB ;
            Point3d ptTempA, ptTempB ;
            int nIntTypeAB = IntersTriaTria( trTriaA, trTriaB, ptTempA, ptTempB, vTriaAB) ;
            if ( nIntTypeAB == ITTT_OVERLAPS) {
               m_vTria[nTA].nTempPart = 2 ;
               SurfB.m_vTria[vNearTria[nTB]].nTempPart = 2 ;
            }
            else if ( nIntTypeAB == ITTT_COUNTER_OVERLAPS) {
               m_vTria[nTA].nTempPart = -2 ;
               SurfB.m_vTria[vNearTria[nTB]].nTempPart = -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) {
            if ( m_vTria[nT].nETempFlag[j] != 0)
               continue ;
            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 e salvo il suo flag
            int nTFlag = m_vTria[nT].nTFlag ;
            RemoveTriangle( nT) ;
           // Inserisco i nuovi triangoli che lo sostituiscono e asegno loro lo stesso flag
            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, nTFlag) ;
            }
         }
      }
   }

   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 e salvo i loro Flag
            int vnTFlag[3] ;
            for ( int nTria = 0 ; nTria < 3 ; ++ nTria) {
               vnTFlag[nTria] = m_vTria[vTriaIndex[nTria]].nTFlag ;
               RemoveTriangle( vTriaIndex[nTria]) ;
            }
           // Aggiungo il nuovo triangolo con il flag
            int nNewTFlag ;
            if ( vnTFlag[0] == vnTFlag[1])
               nNewTFlag = vnTFlag[1] ;
            else if ( vnTFlag[0] == vnTFlag[2]  ||  vnTFlag[1] == vnTFlag[2])
               nNewTFlag = vnTFlag[2] ;
            else
               nNewTFlag = 0 ;
            int nNewVert[3] = { vVertAndAdjTria[0], vVertAndAdjTria[2], vVertAndAdjTria[4] } ;                       
            int nNewTriaId = AddTriangle( nNewVert, nNewTFlag) ;
            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)
{
   m_OGrMgr.Clear() ;
   SurfTriMesh SurfB ;
   SurfB.CopyFrom( &Other) ;
   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) ;
   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 nPrevMaxTFlag = m_nMaxTFlag ;
   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) ;
         }
         if ( nPrevMaxTFlag == m_nMaxTFlag)
            ++ m_nMaxTFlag ;
         AddTriangle( nNewVert, m_nMaxTFlag) ;
      }
   }
   bool bOk = ( AdjustVertices()  &&  DoCompacting()) ;
   bOk && RemoveTripleTriangles() ;
   bOk = bOk && ( AdjustVertices()  &&  DoCompacting()) ;
   bOk && RemoveTJunctions() ;
   bOk = bOk && ( AdjustVertices()  &&  DoCompacting()) ;
   Scale( frScalingRef, 1. / BOOLEAN_SCALE, 1. / BOOLEAN_SCALE, 1. / BOOLEAN_SCALE) ;
   return bOk ;
}

//----------------------------------------------------------------------------
bool
SurfTriMesh::Intersect( const ISurfTriMesh& Other)
{
   m_OGrMgr.Clear() ;
   SurfTriMesh SurfB ;
   SurfB.CopyFrom( &Other) ;
   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) ;
   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 nPrevMaxTFlag = m_nMaxTFlag ;
   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) ;
         }
         if ( nPrevMaxTFlag == m_nMaxTFlag)
            ++ m_nMaxTFlag ;
         AddTriangle( nNewVert, m_nMaxTFlag) ;
      }
   }
   bool bOk = ( AdjustVertices()  &&  DoCompacting()) ;
   bOk && RemoveTripleTriangles() ;
   bOk = bOk && ( AdjustVertices()  &&  DoCompacting()) ;
   bOk && RemoveTJunctions() ;
   bOk = bOk && ( AdjustVertices()  &&  DoCompacting()) ;
   Scale( frScalingRef, 1. / BOOLEAN_SCALE, 1. / BOOLEAN_SCALE, 1. / BOOLEAN_SCALE) ;
   return bOk ;
}

//----------------------------------------------------------------------------
bool
SurfTriMesh::Subtract( const ISurfTriMesh& Other)
{
   m_OGrMgr.Clear() ;
   SurfTriMesh SurfB ;
   SurfB.CopyFrom( &Other) ;
   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) ;
   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 nPrevMaxTFlag = m_nMaxTFlag ;
   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]) ;
         if ( nPrevMaxTFlag == m_nMaxTFlag)
            ++ m_nMaxTFlag ;
         AddTriangle( nNewVert, m_nMaxTFlag) ;
      }
   }
   bool bOk = ( AdjustVertices()  &&  DoCompacting()) ;
   bOk && RemoveTripleTriangles() ;
   bOk = bOk && ( AdjustVertices()  &&  DoCompacting()) ;
   bOk && RemoveTJunctions() ;
   bOk = bOk && ( AdjustVertices()  &&  DoCompacting()) ;
   Scale( frScalingRef, 1. / BOOLEAN_SCALE, 1. / BOOLEAN_SCALE, 1. / BOOLEAN_SCALE) ;
   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 (  ClassifierSurf.GetVertexCount() == 0  ||  ClassifierSurf.GetTriangleCount() == 0)
      return false ;
   if ( m_vVert.empty()  ||  m_vTria.empty()) 
      return true ;
   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) ;
   IdentifyParts() ;
   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].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 ;
}

//----------------------------------------------------------------------------
bool
SurfTriMesh::CutWithOtherSurf( const ISurfTriMesh& CutterSurf, bool bInVsOut, bool bSaveOnEq)
{
   // Le superfici devono essere valide
   if ( ! IsValid()  ||  ! CutterSurf.IsValid())
      return false ;
   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) ;
   IdentifyParts() ;
   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].nTempPart == nPartToRemove  ||  m_vTria[nT].nTempPart == nCoplanarPartToRemove)
         RemoveTriangle( nT) ;
   }
   bool bOk = ( AdjustVertices()  &&  DoCompacting()) ;
   Scale( frScalingRef, 1. / BOOLEAN_SCALE, 1. / BOOLEAN_SCALE, 1. / BOOLEAN_SCALE) ;
   return bOk ;
}

//----------------------------------------------------------------------------
int
SurfTriMesh::VerifyLoopPlane( const PolyLine& ExtLoop, const Plane3d& plCutPlane)
{
  // Verifico il loop e il piano di taglio
   if ( ExtLoop.GetPointNbr() == 0  ||  ! plCutPlane.IsValid())
      return FPI_ERROR ;
  // Calcolo le distanze dei vertici della faccia dal piano
   int nExtLoopPointNum = int( ExtLoop.GetPointNbr()) ;
   DBLVECTOR vDist ;
   vDist.reserve( nExtLoopPointNum) ;
   Point3d ptP ;
   bool bContinue = ExtLoop.GetFirstPoint( ptP) ;
   while ( bContinue) {
      vDist.emplace_back( DistPointPlane( ptP, plCutPlane)) ;
      bContinue = ExtLoop.GetNextPoint( ptP) ;
   }
  // Verifico posizione della faccia rispetto al piano.
   int nVertPos = 0 ; int nVertNeg = 0 ;
   for ( const auto& dDist : vDist) {
      if ( dDist > 0.5 * EPS_SMALL)
         ++ nVertPos ;
      else if ( dDist < - 0.5 * EPS_SMALL)
         ++ nVertNeg ;
   }

  // Tutto sul piano.
   if ( nVertPos == 0  &&  nVertNeg == 0)
      return FPI_ON ;
  // Tutto dentro
   else if ( nVertPos == 0)
      return FPI_IN ;
  // Tutto fuori
   else if ( nVertNeg == 0)
      return FPI_OUT ;
  // Si intersecano
   else
      return FPI_CUT ;
}

//----------------------------------------------------------------------------
int
SurfTriMesh::IntersFacetPlane( const SurfFlatRegion& Region, const Plane3d& plCutPlane,
                               LineFacetClassVector& IntersLinePart)
{
  // Verifico la regione e il piano di taglio
   if ( ! Region.IsValid()  ||  ! plCutPlane.IsValid())
      return FPI_ERROR ;

  // Recupero il piano della regione.
   Plane3d plFacetPlane ;
   plFacetPlane.Set( Region.GetPlanePoint(), Region.GetNormVersor()) ;
   if ( ! plFacetPlane.IsValid())
      return FPI_ERROR ;

  // Pulisco il vettore delle intersezioni
   IntersLinePart.resize( 0) ;

  // Intersezione tra i piani delle due facce
   Point3d ptL ; Vector3d vtL ;
   int nResPP = IntersPlanePlane( plFacetPlane, plCutPlane, ptL, vtL) ;
  // Sono considerate non complanari, se c'è la retta d'intersezione è errore
   if ( nResPP == IPPT_NO  ||  nResPP == IPPT_OVERLAPS)
      return FPI_ERROR ;

  // Box della regione
   BBox3d b3Box ;
   Region.GetLocalBBox( b3Box) ;
   b3Box.Expand( 10) ;

  // Limito la retta nel box contenente la regione
   INTDBLVECTOR vInters ;
   if ( ! IntersLineBox( ptL, vtL, 100., b3Box, vInters, false)  ||  int( vInters.size()) < 2)
      return FPI_ERROR ;
   double dLen = vInters.back().second - vInters.front().second ;
   Point3d ptIntLineSt = ptL + vInters.front().second * vtL ;
   Point3d ptIntLineEn = ptL + vInters.back().second * vtL ;
   CurveLine cvPlaneIntersLine ;
   cvPlaneIntersLine.Set( ptIntLineSt, ptIntLineEn) ;

  // Classifico la linea di intersezione con la faccia
   CRVCVECTOR IntersectionResults ;
   bool bClassificationOk = Region.GetCurveClassification( cvPlaneIntersLine, 0.1 * EPS_SMALL, IntersectionResults) ;
  // Se non trovo le intersezioni, c'è un errore.
   if ( ! bClassificationOk)
      return FPI_ERROR ;

  // Prendo la parte interna o sul bordo.
   int nPartNum = int( IntersectionResults.size()) ;
   for ( int n = 0 ; n < nPartNum ; ++ n) {
      int nType = IntersectionResults[n].nClass ;
      if ( nType == CRVC_IN  ||  nType == CRVC_ON_P  ||  nType == CRVC_ON_M) {
         double dParS = dLen * IntersectionResults[n].dParS ;
         double dParE = dLen * IntersectionResults[n].dParE ;    
         IntersLinePart.emplace_back( LineFacetClass( ptIntLineSt + dParS * vtL, ptIntLineSt + dParE * vtL, nType, nType)) ;
      }
   }
   if ( ! IntersLinePart.empty())
      return FPI_CUT ;

  // La linea si è ridotta ad un punto, la regione è tutta dentro o tutta fuori
   Point3d ptCen ;
   if ( ! Region.GetChunkCentroid( 0, ptCen))
      return FPI_ERROR ;
   if ( DistPointPlane( ptCen, plCutPlane) > 0)
      return FPI_OUT ;
   else
      return FPI_IN ;
}

//----------------------------------------------------------------------------
bool
SurfTriMesh::IntersFacetFacet( const SurfFlatRegion& RegionA, const PolyLine& ExtLoopA,
                               const SurfFlatRegion& RegionB, const PolyLine& ExtLoopB,
                               LineFacetClassVector& IntersLinePart)
{
   IntersLinePart.resize( 0) ;

   // Determino i piani delle regioni
   Plane3d plPlaneA, plPlaneB ;
   plPlaneA.Set( RegionA.GetPlanePoint(), RegionA.GetNormVersor()) ;
   plPlaneB.Set( RegionB.GetPlanePoint(), RegionB.GetNormVersor()) ;

   // Calcolo le distanze dei vertici della faccia A dal piano della faccia B.
   int nExtLoopPointNumA = int( ExtLoopA.GetPointNbr()) ;
   vector<double> vDistA ;
   vDistA.reserve( nExtLoopPointNumA) ;
   Point3d ptP ;
   bool bContinue = ExtLoopA.GetFirstPoint( ptP) ;
   while ( bContinue) {
      vDistA.emplace_back( DistPointPlane( ptP, plPlaneB)) ;
      bContinue = ExtLoopA.GetNextPoint( ptP) ;
   }
   // Verifico posizione della faccia A rispetto al piano della faccia B.
   int nVertPosA = 0 ; int nVertNegA = 0 ;
   for ( const auto& dDist : vDistA) {
      if ( dDist > EPS_SMALL)
         ++ nVertPosA ;
      else if ( dDist < - EPS_SMALL)
         ++ nVertNegA ;
   }
   // Se la faccia A giace tutta da una parte del piano, nessuna intersezione
   if ( nVertPosA == nExtLoopPointNumA  ||  nVertNegA == nExtLoopPointNumA)
      return false ;

   // Calcolo le distanze dei vertici della faccia B dal piano della faccia A.
   int nExtLoopPointNumB = int( ExtLoopB.GetPointNbr()) ;
   vector<double> vDistB ;
   vDistB.reserve( nExtLoopPointNumB) ;
   bContinue = ExtLoopB.GetFirstPoint( ptP) ;
   while ( bContinue) {
      vDistB.emplace_back( DistPointPlane( ptP, plPlaneA)) ;
      bContinue = ExtLoopB.GetNextPoint( ptP) ;
   }
   // Verifico posizione della faccia B rispetto al piano della faccia A.
   int nVertPosB = 0 ; int nVertNegB = 0 ;
   for ( const auto& dDist : vDistB) {
      if ( dDist > EPS_SMALL)
         ++ nVertPosB ;
      else if ( dDist < - EPS_SMALL)
         ++ nVertNegB ;
   }
   // Se la faccia B giace tutta da una parte del piano, nessuna intersezione
   if ( nVertPosB == nExtLoopPointNumB  ||  nVertNegB == nExtLoopPointNumB)
      return false ;

   // Intersezione tra i piani delle due facce
   Point3d ptL ; Vector3d vtL ;
   int nResPP = IntersPlanePlane( plPlaneA, plPlaneB, ptL, vtL) ;
   if ( nResPP == IPPT_NO)
      return false ;

   // Se le facce sono complanari
   if ( nResPP == IPPT_OVERLAPS)
      return false ;
   
   // Box contenente entrambe le regioni
   BBox3d b3Box, b3BoxB;
   RegionA.GetLocalBBox( b3Box) ;
   RegionB.GetLocalBBox( b3BoxB) ;
   b3Box.Add(b3BoxB);
   b3Box.Expand( 10);

   // Limito la retta nel box contenente le regioni
   INTDBLVECTOR vInters;
   if ( ! IntersLineBox(ptL, vtL, 100., b3Box, vInters, false)  ||  int( vInters.size()) </*!=*/ 2)
      return false ;
   double dLen = vInters.back()/*[1]*/.second - vInters[0].second ;
   Point3d ptIntLineSt = ptL + vInters[0].second * vtL;
   Point3d ptIntLineEn = ptL + vInters.back()/*[1]*/.second * vtL;
   CurveLine cvPlaneIntersLine;
   cvPlaneIntersLine.Set( ptIntLineSt, ptIntLineEn);
   
   // Limito la linea di intersezione con la faccia A
   CRVCVECTOR IntersectionResultsA ;
   bool bClassificationOk = RegionA.GetCurveClassification( cvPlaneIntersLine, 0.1 * EPS_SMALL, IntersectionResultsA) ;

   // Limito la linea di intersezione con la faccia B
   CRVCVECTOR IntersectionResultsB ;
   bClassificationOk = bClassificationOk  &&  RegionB.GetCurveClassification( cvPlaneIntersLine, 0.1 * EPS_SMALL, IntersectionResultsB) ;

   if ( ! bClassificationOk)
      return false ;

   // Prendo la parte comune
   int nPartNumA = int( IntersectionResultsA.size()) ;
   int nPartNumB = int( IntersectionResultsB.size()) ;
   for ( int nA = 0 ; nA < nPartNumA ; ++ nA) {
      int nTypeA = IntersectionResultsA[nA].nClass ;
      if ( nTypeA == CRVC_IN  ||  nTypeA == CRVC_ON_P  ||  nTypeA == CRVC_ON_M) {
         for ( int nB = 0 ; nB < nPartNumB ; ++ nB) {
            int nTypeB = IntersectionResultsB[nB].nClass ;
            if ( nTypeB == CRVC_IN  ||  nTypeB == CRVC_ON_P  ||  nTypeB == CRVC_ON_M) {
               double dParSA = dLen * IntersectionResultsA[nA].dParS ;
               double dParEA = dLen * IntersectionResultsA[nA].dParE ;
               double dParSB = dLen * IntersectionResultsB[nB].dParS ;
               double dParEB = dLen * IntersectionResultsB[nB].dParE ;
               if ( dParSB < dParEA - EPS_SMALL  &&
                    dParEB > dParSA + EPS_SMALL) {
                  double dMinPar = max( dParSA, dParSB) ;
                  double dMaxPar = min( dParEA, dParEB) ;
                  IntersLinePart.emplace_back( LineFacetClass( ptIntLineSt + dMinPar * vtL, ptIntLineSt + dMaxPar * vtL, nTypeA, nTypeB)) ;
               }
            }
         }
      }
   }

   return int( IntersLinePart.size()) > 0 ;
}

//----------------------------------------------------------------------------
//bool
//CreateFlatRegionFromPolyLineVector( const POLYLINEVECTOR& PolyVec, SurfFlatRegion& Region)
//{
//   bool bOk = true ;
//   for ( int nL = 0 ; nL < int( PolyVec.size())  &&  bOk ; ++ nL) {
//      CurveComposite Loop ;
//      Point3d ptSt, ptEn ;
//      bool bContinue = PolyVec[nL].GetFirstPoint( ptSt)  &&  PolyVec[nL].GetNextPoint( ptEn) ;
//      while ( bContinue) {
//         CurveLine cvLine ;
//         cvLine.Set( ptSt, ptEn) ;
//         Loop.AddCurve( cvLine) ;
//         ptSt = ptEn ;
//         bContinue = PolyVec[nL].GetNextPoint( ptEn) ;
//      }
//      if ( nL == 0)
//         bOk = bOk  &&  Region.AddExtLoop( Loop) ;
//      else
//         bOk = bOk  &&  Region.AddIntLoop( Loop) ;
//   }
//   return bOk ;
//}

//----------------------------------------------------------------------------
struct PositionOnPolyLine {
   int nIndexInVec ;
   double dPar ;
   PositionOnPolyLine( int nIndex, double dParam) {
      nIndexInVec = nIndex ;
      dPar = dParam ;
   }
} ;

//----------------------------------------------------------------------------
bool
SurfTriMesh::SplitFacet( const INTERSCHAINMAP& IntersLineMap, PieceMap& NewFacet)
{
   for ( auto it = IntersLineMap.begin() ; it != IntersLineMap.end() ; ++ it) {
     // Normale alla faccia
      Vector3d vtFacetNorm ;
      GetFacetNormal( it->first, vtFacetNorm) ;
     // Creo i loop
      ChainCurves LoopCreator ;
      LoopCreator.Init( false, 10 * EPS_SMALL, int( it->second.size()), 10) ;
     // 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 e li divido fra chiusi e aperti.
      INNCHAINVECTOR cvClosedChain ;
      INNCHAINVECTOR cvOpenChain ;
      INTVECTOR vIds ;
      Point3d ptNearStart = ( it->second.size() > 0 ? it->second[0].ptSt : ORIG) ;
      while ( LoopCreator.GetChainFromNear( ptNearStart, false, vIds)) {
         IntersInnChain chTemp ;
        // Aggiungo la linea alla curva composta.
         for ( auto i : vIds)
            chTemp.emplace_back( it->second[i - 1]) ;
        // Inserisco opportunamente negli aperti o chiusi
         int nCurLoopLast = int( chTemp.size()) - 1 ;
         if ( nCurLoopLast > 0  &&  AreSamePointEpsilon( chTemp[0].ptSt, chTemp[nCurLoopLast].ptEn, 10 * EPS_SMALL))
            cvClosedChain.emplace_back( chTemp) ;
         else
            cvOpenChain.emplace_back( chTemp) ;
      }
     // Elimino la seconda copia di catene aperte doppie
      for ( int nI = 0 ; nI < int( cvOpenChain.size()) - 1 ; ++ nI) {
         for ( int nJ = nI + 1 ; nJ < int( cvOpenChain.size()) ; ++ nJ) {
            if ( cvOpenChain[nI].size() == cvOpenChain[nJ].size()) {
               bool bSame = true ;
               for ( int nK = 0 ; nK < int( cvOpenChain[nI].size()) ; ++ nK) {
                  if ( ! AreSamePointEpsilon( cvOpenChain[nI][nK].ptSt, cvOpenChain[nJ][nK].ptSt, EPS_SMALL)  ||
                       ! AreSamePointEpsilon( cvOpenChain[nI][nK].ptEn, cvOpenChain[nJ][nK].ptEn, EPS_SMALL)) {
                     bSame = false ;
                     break ;
                  }
               }
               if ( bSame) {
                  cvOpenChain.erase( cvOpenChain.begin() + nJ) ;
                  -- nJ ;
               }
            }
         }
      }
     // Elimino la seconda copia di catene chiuse doppie
      for ( int nI = 0 ; nI < int( cvClosedChain.size()) - 1 ; ++ nI) {
         for ( int nJ = nI + 1 ; nJ < int( cvClosedChain.size()) ; ++ nJ) {
            if ( cvClosedChain[nI].size() == cvClosedChain[nJ].size()) {
               bool bSame = true ;
               for ( int nK = 0 ; nK < int( cvClosedChain[nI].size()) ; ++ nK) {
                  if ( ! AreSamePointEpsilon( cvClosedChain[nI][nK].ptSt, cvClosedChain[nJ][nK].ptSt, EPS_SMALL)  ||
                       ! AreSamePointEpsilon( cvClosedChain[nI][nK].ptEn, cvClosedChain[nJ][nK].ptEn, EPS_SMALL)) {
                     bSame = false ;
                     break ;
                  }
               }
               if ( bSame) {
                  cvClosedChain.erase( cvClosedChain.begin() + nJ) ;
                  -- nJ ;
               }
            }
         }
      }
     // Recupero i loop della faccia da dividere
      vector<FacetPiece>& vNewPieces = ( NewFacet.emplace( it->first, vector<FacetPiece>()).first)->second ;
      vNewPieces.emplace_back();
      vNewPieces.back().nPiecePart = 1 ;
      GetFacetLoops( it->first, vNewPieces.back().vPieceLoop) ;
     // Divido la faccia in parti.
     // Ciclo sulle catene finché non esauriscono.
      int nIterNumber = 0 ;
      int nPrevLastChainNum = int( cvOpenChain.size()) ;
      while ( int( cvOpenChain.size()) > 0) {
         // Controllo di non essere in un loop infinito.
         int nLastChainNum = int( cvOpenChain.size()) - 1 ;
         if ( nPrevLastChainNum == nLastChainNum) {
            cvOpenChain.insert( cvOpenChain.begin(), cvOpenChain.back()) ;
            cvOpenChain.pop_back() ;
            ++ nIterNumber ;
         }
         else {
            nIterNumber = 0 ;
            nPrevLastChainNum = nLastChainNum;
         }
         if ( nIterNumber == int( cvOpenChain.size())) {
            cvOpenChain.resize( int( cvOpenChain.size()) - 1) ;
            -- nLastChainNum ;
         }
         if ( int( cvOpenChain.size()) == 0)
            break ;
        // Ciclo su tutte le parti della faccia.
         int nPartLoopNum = int( vNewPieces.size()) ;
         for ( int nPart = 0 ; nPart < nPartLoopNum ; ++ nPart) {
           // Cerco i loop su cui la catena corrente inizia e finisce.
            INTVECTOR vLoopIndexes ;
           // Ciclo sui loop della faccia.
            for ( int nLoop = 0 ; nLoop < int( vNewPieces[nPart].vPieceLoop.size())  &&  int( vLoopIndexes.size()) < 2 ; ++ nLoop) {
               double dPar ;
               if ( GetPointParamOnPolyLine( cvOpenChain[nLastChainNum][0].ptSt, vNewPieces[nPart].vPieceLoop[nLoop], 2 * EPS_SMALL, dPar)) {
                  vLoopIndexes.emplace_back( nLoop) ;
                  if ( vLoopIndexes.size() > 1) {
                     swap( vLoopIndexes[0], vLoopIndexes[1]) ;
                  }
               }  
               if ( GetPointParamOnPolyLine( cvOpenChain[nLastChainNum][cvOpenChain[nLastChainNum].size() - 1].ptEn,
                                                  vNewPieces[nPart].vPieceLoop[nLoop], 2 * EPS_SMALL, dPar)) {
                  vLoopIndexes.emplace_back( nLoop) ;
               }
            }
           // Se non ho trovato almeno un taglio completo, vado al prossimo.
            if ( int( vLoopIndexes.size()) < 2)
               continue ;
           // La catena finisce sul loop ove inizia. Divido la nuova parte.
            if ( vLoopIndexes[0] == vLoopIndexes[1]) {
              // Cambio inizio al loop iniziale.
               ChangePolyLineStart( vNewPieces[nPart].vPieceLoop[vLoopIndexes[0]], cvOpenChain[nLastChainNum].back().ptEn, 2 * EPS_SMALL) ;
              // Divido il loop della parte
              // Loop1
               PolyLine NewLoop1 ;
              // Inserisco i punti della catena nella PolyLine del nuovo loop.
               for ( int m = 0 ; m < int( cvOpenChain[nLastChainNum].size()) ; ++ m) {
                  NewLoop1.AddUPoint( 0., cvOpenChain[nLastChainNum][m].ptSt) ;
                  if ( m == int( cvOpenChain[nLastChainNum].size()) - 1)
                     NewLoop1.AddUPoint( 0., cvOpenChain[nLastChainNum][m].ptEn) ;
               }
              // Spezzo il loop successivo alla catena nel punto in cui comincia la catena successiva.
               PolyLine SplitLoop1, SplitLoop2 ;
               SplitPolyLineAtPoint( vNewPieces[nPart].vPieceLoop[vLoopIndexes[0]],
                                     cvOpenChain[nLastChainNum][0].ptSt, 2 * EPS_SMALL,
                                     SplitLoop1, SplitLoop2) ;
              // Aggiungo i punti precedenti il punto di frattura nella NewLoop1.
               NewLoop1.Join( SplitLoop1) ;
              // Loop2
               PolyLine NewLoop2 ;
               NewLoop2.Join( SplitLoop2) ;
              // Inserisco i punti della catena nella PolyLine del nuovo loop.
               for ( int m = int( cvOpenChain[nLastChainNum].size()) - 1 ; m >= 0 ; -- m) {
                  NewLoop2.AddUPoint( 0., cvOpenChain[nLastChainNum][m].ptEn) ;
                  if ( m == 0)
                     NewLoop2.AddUPoint( 0., cvOpenChain[nLastChainNum][m].ptSt) ;
               }
              // Loop esterno coinvolto
               if ( vLoopIndexes[0] == 0) {
                 // Creo i pezzi nuovi, i cui loop esterni sono quelli appena definiti.
                  FacetPiece PieceInn, PieceOut ;
                  PieceInn.vPieceLoop.emplace_back( NewLoop1) ;
                  PieceInn.nPiecePart = 1 ;
                  PieceOut.vPieceLoop.emplace_back( NewLoop2) ;
                  PieceOut.nPiecePart = - 1 ;
                 // Assegno i loop interni del vecchio pezzo non tagliati da catene ai nuovi.
                  for ( int nIL = 1 ; nIL < int( vNewPieces[nPart].vPieceLoop.size()) ; ++ nIL) {
                     if ( nIL == vLoopIndexes[0])
                        continue ;   
                     Point3d ptPointInnerLoop ;
                     vNewPieces[nPart].vPieceLoop[nIL].GetFirstPoint( ptPointInnerLoop) ;
                     if ( PieceInn.vPieceLoop[0].GetPointNbr() < PieceOut.vPieceLoop[0].GetPointNbr()) {
                        if ( IsPointInsidePolyLine( ptPointInnerLoop, PieceInn.vPieceLoop[0], EPS_SMALL)) {
                           PieceInn.vPieceLoop.emplace_back( vNewPieces[nPart].vPieceLoop[nIL]) ;
                        }
                        else {
                           PieceOut.vPieceLoop.emplace_back( vNewPieces[nPart].vPieceLoop[nIL]) ;
                        }
                     }
                     else {
                        if ( IsPointInsidePolyLine( ptPointInnerLoop, PieceOut.vPieceLoop[0], EPS_SMALL)) {
                           PieceOut.vPieceLoop.emplace_back( vNewPieces[nPart].vPieceLoop[nIL]) ;
                        }
                        else {
                           PieceInn.vPieceLoop.emplace_back( vNewPieces[nPart].vPieceLoop[nIL]) ;
                        }
                     }
                  }
                  vNewPieces.erase( vNewPieces.begin() + nPart) ;
                  vNewPieces.emplace_back( PieceInn) ;
                  vNewPieces.emplace_back( PieceOut) ;
               }
              // Loop esterno non coinvolto
               else {
                  Plane3d plLoopPlane1, plLoopPlane2, plContLoopPlane;
                  double dArea1, dArea2, dAreaCont;
                  NewLoop1.IsClosedAndFlat(plLoopPlane1, dArea1);
                  NewLoop2.IsClosedAndFlat(plLoopPlane2, dArea2);
                  vNewPieces[nPart].vPieceLoop[0].IsClosedAndFlat(plContLoopPlane, dAreaCont);
                  Vector3d vtLoopPlaneNorm1 = plLoopPlane1.GetVersN();
                  Vector3d vtContLoopPlaneNorm = plContLoopPlane.GetVersN();
                  PolyLine NewLoopCCW, NewLoopCW;
                  bool bFirstLoopIsCounter = vtLoopPlaneNorm1 * vtContLoopPlaneNorm > EPS_SMALL;
                  if (bFirstLoopIsCounter) {
                     NewLoopCCW = NewLoop1;
                     NewLoopCW = NewLoop2;
                  }
                  else {
                     NewLoopCCW = NewLoop2;
                     NewLoopCW = NewLoop1;
                  }
                 // Creo il pezzo nuovo, il cui loop esterno è stato appena definito.
                  FacetPiece DetachedPiece;
                  DetachedPiece.vPieceLoop.emplace_back(NewLoopCCW);
                  DetachedPiece.nPiecePart = bFirstLoopIsCounter ? 1 : -1;
                 // Al pezzo staccato assegno i loop interni del vecchio pezzo non tagliati da catene.
                  for ( int nIL = 1 ; nIL < int( vNewPieces[nPart].vPieceLoop.size()) ; ++ nIL) {
                     if ( nIL == vLoopIndexes[0])
                        continue ;
                     Point3d ptPointInnerLoop ;
                     vNewPieces[nPart].vPieceLoop[nIL].GetFirstPoint( ptPointInnerLoop) ;
                     if ( IsPointInsidePolyLine( ptPointInnerLoop, DetachedPiece.vPieceLoop[0], EPS_SMALL)) {
                        DetachedPiece.vPieceLoop.emplace_back( vNewPieces[nPart].vPieceLoop[nIL]) ;
                        vNewPieces[nPart].vPieceLoop.erase( vNewPieces[nPart].vPieceLoop.begin() + nIL) ;
                       // Cambio il numero dei loop tagliati
                        for ( int nIntersLoop = 0 ; nIntersLoop < int( vLoopIndexes.size()) ; ++ nIntersLoop) {
                           if ( vLoopIndexes[nIntersLoop] > nIL)
                              -- vLoopIndexes[nIntersLoop] ;
                        }
                        -- nIL ;
                     }
                  }
                 // Aggiungo al pezzo, da cui si stacca quello appena creato, il loop complementare a quello che appartiene al pezzo che si stacca.
                  /*if ( ! bFirstLoopIsCounter)*/
                     vNewPieces[nPart].vPieceLoop.erase(vNewPieces[nPart].vPieceLoop.begin() + vLoopIndexes[0]);
                  vNewPieces[nPart].vPieceLoop.emplace_back( NewLoopCW) ;
                  vNewPieces[nPart].nPiecePart = bFirstLoopIsCounter ? - 1 : 1 ;
                  vNewPieces.emplace_back( DetachedPiece) ; 
               }
              // Elimino la catena usata.
               cvOpenChain.erase( cvOpenChain.begin() + nLastChainNum) ;
              // Interrompo il ciclo sulle parti.
               break ;
            }
           // Cerco le catene per dividere la nuova parte.
            else {
               INTVECTOR vChainIndex ;
               vChainIndex.emplace_back( nLastChainNum) ;
               while ( vLoopIndexes.back() != vLoopIndexes[0]) {
                 // Cambio inizio del loop corrente in modo che il punto di inizio sia il punto finale dell'ultima catena trovata.
                  ChangePolyLineStart( vNewPieces[nPart].vPieceLoop[vLoopIndexes.back()], cvOpenChain[vChainIndex.back()].back().ptEn, 2 * EPS_SMALL) ;
                 // Cerco catene che iniziano sul loop coorente
                  vector<PositionOnPolyLine> vChainStartingOnLoop ;
                  for ( int nCh = 0 ; nCh < nLastChainNum ; ++ nCh) {
                     double dPar ;
                     if ( GetPointParamOnPolyLine( cvOpenChain[nCh][0].ptSt, vNewPieces[nPart].vPieceLoop[vLoopIndexes.back()], 2 * EPS_SMALL, dPar)) {
                        vChainStartingOnLoop.emplace_back( PositionOnPolyLine(nCh, dPar)) ;
                     }
                  }
                 // Ordino le catene secondo la vicinanza lungo il loop del loro punto d'inizio al punto d'inizio del loop stesso.
                  sort( vChainStartingOnLoop.begin(), vChainStartingOnLoop.end(), [] ( PositionOnPolyLine Ch1, PositionOnPolyLine Ch2) {
                                                                                          return Ch1.dPar < Ch2.dPar ; } ) ;
                 // Cerco la prima catena che non termina sul loop corrente.
                  for ( int n = 0 ; n < int( vChainStartingOnLoop.size()) ; ++ n) {
                    // Indice della catena e indice del segmento finale
                     int nCh = vChainStartingOnLoop[n].nIndexInVec ;
                     int nLastLineIndex = int( cvOpenChain[nCh].size()) - 1 ;
                    // Ciclo sui loop
                     int nLoopNum = int( vNewPieces[nPart].vPieceLoop.size()) ;
                     bool bFound = false ;
                     for ( int nLoop = 0 ; nLoop < nLoopNum ; ++ nLoop) {
                       // Salto il loop corrente
                        if ( nLoop == vLoopIndexes.back())
                           continue ;
                        double dPar ;
                       // La catena termina su questo loop diverso da quello corrente.
                        if ( GetPointParamOnPolyLine( cvOpenChain[nCh][nLastLineIndex].ptEn, vNewPieces[nPart].vPieceLoop[nLoop], 2 * EPS_SMALL, dPar)) {
                          // Salvo l'indice della catena e quello del loop.
                           vChainIndex.emplace_back( nCh) ;
                           vLoopIndexes.emplace_back( nLoop) ;
                           bFound = true ;
                           break ;
                        }
                     }
                     if ( bFound)
                        break ;
                  }
               }
              // Cambio inizio al loop iniziale.
               ChangePolyLineStart( vNewPieces[nPart].vPieceLoop[vLoopIndexes[0]], cvOpenChain[vChainIndex.back()].back().ptEn, 2 * EPS_SMALL) ;
              // Divido i loop della parte.
               POLYLINEVECTOR vPolySecondPartVec ;
               PolyLine NewLoop1 ;
               for ( int n = 0 ; n < int( vChainIndex.size()) ; ++ n) {
                 // Inserisco i punti della catena nella PolyLine del nuovo loop.
                  for ( int m = 0 ; m < int( cvOpenChain[vChainIndex[n]].size()) ; ++ m) {
                     NewLoop1.AddUPoint( 0., cvOpenChain[vChainIndex[n]][m].ptSt) ;
                     if ( m == int( cvOpenChain[vChainIndex[n]].size()) - 1)
                        NewLoop1.AddUPoint( 0., cvOpenChain[vChainIndex[n]][m].ptEn) ;
                  }
                 // Spezzo il loop successivo alla catena nel punto in cui comincia la catena successiva.
                  PolyLine SplitLoop1, SplitLoop2 ;
                  SplitPolyLineAtPoint( vNewPieces[nPart].vPieceLoop[vLoopIndexes[n + 1]],
                                        cvOpenChain[vChainIndex[( n + 1) % int( vChainIndex.size())]][0].ptSt, 2 * EPS_SMALL,
                                        SplitLoop1, SplitLoop2) ;
                 // Aggiungo i punti precedenti il punto di frattura in NewLoop1.
                  NewLoop1.Join( SplitLoop1) ;
                 // Salvo la parte dopo il punto di frattura in un apposito vettore.
                  vPolySecondPartVec.emplace_back( SplitLoop2) ;
               }
               PolyLine NewLoop2 ;
               for ( int n = int( vChainIndex.size()) - 1 ; n >= 0 ; -- n) {
                 // Aggiungo i punti successivi il punto di frattura in NewLoop2.
                  NewLoop2.Join( vPolySecondPartVec[n]) ;
                 // Inserisco i punti della catena nella PolyLine del nuovo loop.
                  for ( int m = int( cvOpenChain[vChainIndex[n]].size()) - 1 ; m >= 0 ; -- m) {
                     NewLoop2.AddUPoint( 0., cvOpenChain[vChainIndex[n]][m].ptEn) ;
                     if ( m == 0)
                        NewLoop2.AddUPoint( 0., cvOpenChain[vChainIndex[n]][m].ptSt) ;
                  }
               }
              // Controllo che il loop esterno sia interessato.
               bool bExtCutted = false ;
               for ( int n = 0 ; n < int( vLoopIndexes.size())  &&  ! bExtCutted ; ++ n) {
                  if ( vLoopIndexes[n] == 0)
                     bExtCutted = true ;
               }
              // Loop esterno interessato
               if ( bExtCutted) {
                 // Creo i pezzi nuovi, i cui loop esterni sono quelli appena definiti.
                  FacetPiece PieceInn, PieceOut ;
                  PieceInn.vPieceLoop.emplace_back( NewLoop1) ;
                  PieceInn.nPiecePart = 1 ;
                  PieceOut.vPieceLoop.emplace_back( NewLoop2) ;
                  PieceOut.nPiecePart = - 1 ;
                 // Assegno i loop interni del vecchio pezzo non tagliati da catene ai nuovi.
                  for ( int nIL = 1 ; nIL < int( vNewPieces[nPart].vPieceLoop.size()) ; ++ nIL) {
                     bool bUsed = false ;
                     for ( int nU = 0 ; nU < int( vLoopIndexes.size()) - 1 ; ++ nU) {
                        if ( nIL == vLoopIndexes[nU]) {
                           bUsed = true ;
                           break ;
                        }
                     }
                     if ( bUsed)
                        continue ;
                     Point3d ptPointInnerLoop ;
                     vNewPieces[nPart].vPieceLoop[nIL].GetFirstPoint( ptPointInnerLoop) ;
                     if ( PieceInn.vPieceLoop[0].GetPointNbr() < PieceOut.vPieceLoop[0].GetPointNbr()) {
                        if ( IsPointInsidePolyLine( ptPointInnerLoop, PieceInn.vPieceLoop[0], EPS_SMALL)) {
                           PieceInn.vPieceLoop.emplace_back( vNewPieces[nPart].vPieceLoop[nIL]) ;
                        }
                        else {
                           PieceOut.vPieceLoop.emplace_back( vNewPieces[nPart].vPieceLoop[nIL]) ;
                        }
                     }
                     else {
                        if ( IsPointInsidePolyLine( ptPointInnerLoop, PieceOut.vPieceLoop[0], EPS_SMALL)) {
                           PieceOut.vPieceLoop.emplace_back( vNewPieces[nPart].vPieceLoop[nIL]) ;
                        }
                        else {
                           PieceInn.vPieceLoop.emplace_back( vNewPieces[nPart].vPieceLoop[nIL]) ;
                        }
                     }
                  }
                 // Aggiungo i due nuovi pezzi ed elimino quello da cui sono nati.
                  vNewPieces.erase( vNewPieces.begin() + nPart) ;
                  vNewPieces.emplace_back( PieceInn) ;
                  vNewPieces.emplace_back( PieceOut) ;
               }
              // Loop esterno non interessato
               else {
                  Plane3d plLoopPlane1, plLoopPlane2, plContLoopPlane ;
                  double dArea1, dArea2, dAreaCont ;
                  NewLoop1.IsClosedAndFlat( plLoopPlane1, dArea1) ;
                  NewLoop2.IsClosedAndFlat( plLoopPlane2, dArea2) ;
                  vNewPieces[nPart].vPieceLoop[0].IsClosedAndFlat( plContLoopPlane, dAreaCont) ;
                  Vector3d vtLoopPlaneNorm1 = plLoopPlane1.GetVersN();
                  Vector3d vtContLoopPlaneNorm = plContLoopPlane.GetVersN();
                  PolyLine NewLoopCCW, NewLoopCW;
                  bool bFirstLoopIsCounter = vtLoopPlaneNorm1 * vtContLoopPlaneNorm > EPS_SMALL;
                  if ( bFirstLoopIsCounter) {
                     NewLoopCCW = NewLoop1 ;
                     NewLoopCW = NewLoop2 ;
                  }
                  else {
                     NewLoopCCW = NewLoop2 ;
                     NewLoopCW = NewLoop1 ;
                  }
                 // // Creo il pezzo nuovo, il cui loop esterno è stato appena definito.
                  // FacetPiece PieceInn ;
                  // PieceInn.vPieceLoop.emplace_back( NewLoop1) ;
                  // PieceInn.nPiecePart = 1 ;
                  // // Assegno i loop interni del vecchio pezzo non tagliati da catene ai nuovi.
                  // for ( int nIL = 1 ; nIL < int( vNewPieces[nPart].vPieceLoop.size()) ; ++ nIL) {
                  //    bool bUsed = false ;
                  //    for ( int nU = 1 ; nU < int( vLoopIndexes.size()) - 1 ; ++ nU) {
                  //       if ( nIL == vLoopIndexes[nU]) {
                  //          bUsed = true ;
                  //          break ;
                  //       }
                  //    }
                  //    if ( bUsed)
                  //       continue ;
                  //    Point3d ptPointInnerLoop ;
                  //    vNewPieces[nPart].vPieceLoop[nIL].GetFirstPoint( ptPointInnerLoop) ;  
                  //    if ( IsPointInsidePolyLine( ptPointInnerLoop, PieceInn.vPieceLoop[0])) {
                  //       PieceInn.vPieceLoop.emplace_back( vNewPieces[nPart].vPieceLoop[nIL]) ;
                  //       vNewPieces[nPart].vPieceLoop.erase( vNewPieces[nPart].vPieceLoop.begin() + nIL) ;
                  //       -- nIL ;
                  //    }
                  // }
                  // // Aggiungo al pezzo, da cui si stacca quello appena creato, il loop complementare a quello che appartiene al nuovo pezzo.
                  // vNewPieces[nPart].vPieceLoop.emplace_back( NewLoop2) ;
                  // vNewPieces[nPart].nPiecePart = - 1 ;
                  // vNewPieces.emplace_back( PieceInn) ;
                 // Creo il pezzo nuovo, il cui loop esterno è stato appena definito.
                  FacetPiece DetachedPiece ;
                  DetachedPiece.vPieceLoop.emplace_back( NewLoopCCW) ;
                  DetachedPiece.nPiecePart = bFirstLoopIsCounter ? 1 : - 1 ;
                 // Al pezzo che si stacca assegno i loop interni del vecchio pezzo non tagliati da catene.
                  for ( int nIL = 1 ; nIL < int( vNewPieces[nPart].vPieceLoop.size()) ; ++ nIL) {
                     bool bUsed = false ;
                     for ( int nU = 1 ; nU < int( vLoopIndexes.size()) - 1 ; ++ nU) {
                        if ( nIL == vLoopIndexes[nU]) {
                           bUsed = true ;
                           break ;
                        }
                     }
                     if ( bUsed)
                        continue ;
                     Point3d ptPointInnerLoop ;
                     vNewPieces[nPart].vPieceLoop[nIL].GetFirstPoint( ptPointInnerLoop) ;  
                     if ( IsPointInsidePolyLine( ptPointInnerLoop, DetachedPiece.vPieceLoop[0], EPS_SMALL)) {
                        DetachedPiece.vPieceLoop.emplace_back( vNewPieces[nPart].vPieceLoop[nIL]) ;
                        vNewPieces[nPart].vPieceLoop.erase( vNewPieces[nPart].vPieceLoop.begin() + nIL) ;
                       // Cambio il numero dei loop talgiati
                        for ( int nIntersLoop = 0 ; nIntersLoop < int( vLoopIndexes.size()) ; ++ nIntersLoop) {
                           if ( vLoopIndexes[nIntersLoop] > nIL)
                              -- vLoopIndexes[nIntersLoop] ;
                        }
                        -- nIL ;
                     }
                  }
                 // Aggiungo al pezzo, da cui si stacca quello appena creato, il loop complementare a quello che appartiene al pezzo che si stacca.
                  vLoopIndexes.resize(int(vLoopIndexes.size()) - 1);
                  sort( vLoopIndexes.begin(), vLoopIndexes.end(), [] ( int nIndex1, int nIndex2) { return nIndex1 > nIndex2 ; }) ;
                  for (int nLoopInd = 0; nLoopInd < int(vLoopIndexes.size()); ++nLoopInd) {
                     vNewPieces[nPart].vPieceLoop.erase(vNewPieces[nPart].vPieceLoop.begin() + vLoopIndexes[nLoopInd]);
                  } 
                  vNewPieces[nPart].vPieceLoop.emplace_back( NewLoopCW) ;
                  vNewPieces[nPart].nPiecePart = bFirstLoopIsCounter ? - 1 : 1 ;
                  vNewPieces.emplace_back( DetachedPiece) ;
               }
              // Elimino catene usate.
               sort( vChainIndex.begin(), vChainIndex.end(), [] ( int nIndex1, int nIndex2) { return nIndex1 > nIndex2 ; }) ;
               for ( int nChInd = 0 ; nChInd < int( vChainIndex.size()) ; ++ nChInd) {
                  cvOpenChain.erase( cvOpenChain.begin() + vChainIndex[nChInd]) ;
               }
              // Interrompo il ciclo sulle parti.
               break ;
            }
         }
      }
     // Assegno i loop chiusi interni alle rispettive parti di faccia.
     // Ciclo sui pezzi di facet.
      int nPieceNum = int( vNewPieces.size()) ;
      for ( int nPieceN = 0 ; nPieceN < nPieceNum ; ++ nPieceN) {
        // Ciclo sui loop interni
         for ( int nInnL = 0 ; nInnL < int( cvClosedChain.size()) ; ++ nInnL) {
           // Trasformo il loop interno in PolyLine.
            PolyLine CurInnerLoop ;
            CurInnerLoop.AddUPoint( 0., cvClosedChain[nInnL][0].ptSt) ;
            for ( int nSeg = 0 ; nSeg < int( cvClosedChain[nInnL].size()) ; ++ nSeg) {
               CurInnerLoop.AddUPoint( 0., cvClosedChain[nInnL][nSeg].ptEn) ;
            }
            Plane3d plPlane ;
            double dArea ;
            if ( ! CurInnerLoop.IsClosedAndFlat( plPlane, dArea)  ||  dArea < EPS_SMALL)
               continue ;
           // Ciclo sui punti del loop interno.
            PNTULIST& LoopList = CurInnerLoop.GetUPointList() ;
            for ( auto itInn = LoopList.begin() ; itInn != LoopList.end() ; ++ itInn) {
               Point3d ptInnP = itInn->first ;
               Point3d ptInnNextP ;
               auto itInnNext = itInn ;
               ++ itInnNext ;
               if ( itInnNext != LoopList.end())
                  ptInnNextP = itInnNext->first ;
               if ( IsPointInsidePolyLine( ptInnP, vNewPieces[nPieceN].vPieceLoop[0], EPS_SMALL)  ||
                    ( itInnNext != LoopList.end()  &&
                      IsPointInsidePolyLine( 0.5 * ( ptInnP + ptInnNextP), vNewPieces[nPieceN].vPieceLoop[0], EPS_SMALL))) {
                  Polygon3d AuxPolygon ;
                  AuxPolygon.FromPolyLine( CurInnerLoop) ;
                  Vector3d vtInnLoopNorm =  AuxPolygon.GetVersN() ;
                  if ( vtFacetNorm * vtInnLoopNorm < 0.) {
                    // Aggiungo loop al pezzo.
                     vNewPieces[nPieceN].nPiecePart = 1 ;
                     vNewPieces[nPieceN].vPieceLoop.emplace_back( CurInnerLoop) ;
                    // Aggiungo nuovo pezzo.
                     vNewPieces.emplace_back() ;
                     vNewPieces.back().vPieceLoop.emplace_back( CurInnerLoop) ;
                     vNewPieces.back().vPieceLoop.back().Invert() ;
                     vNewPieces.back().nPiecePart = - 1 ;
                    // Cerco loop interni a quello appena aggiunto
                     INTVECTOR vSecondLevel ;
                     for ( int nSI = 1 ; nSI < int( vNewPieces[nPieceN].vPieceLoop.size()) ; ++ nSI) {
                        PNTULIST& SecLevLoopList = vNewPieces[nPieceN].vPieceLoop[nSI].GetUPointList() ;
                        for ( auto itSL = SecLevLoopList.begin() ; itSL != SecLevLoopList.end() ; ++ itSL) {
                           Point3d ptPntSL = itSL->first ;
                           if ( IsPointInsidePolyLine( ptPntSL, vNewPieces.back().vPieceLoop.back(), EPS_SMALL)) {
                              vSecondLevel.emplace_back( nSI) ;
                              break ;
                           }
                        }
                     }
                    // Aggiungo i loop al nuovo pezzo e li tolgo dal precedente.
                     for ( int nSI = 0 ; nSI < int( vSecondLevel.size()) ; ++ nSI)
                        vNewPieces.back().vPieceLoop.emplace_back( vNewPieces[nPieceN].vPieceLoop[vSecondLevel[nSI]]) ;
                     for ( int nSI = 0 ; nSI < int( vSecondLevel.size()) ; ++ nSI)
                        vNewPieces[nPieceN].vPieceLoop.erase( vNewPieces[nPieceN].vPieceLoop.begin() + vSecondLevel[nSI]) ;
                     break ;
                  }
                  else {
                    // Aggiungo loop al pezzo.
                     vNewPieces[nPieceN].nPiecePart = - 1 ;
                     vNewPieces[nPieceN].vPieceLoop.emplace_back( CurInnerLoop) ;
                     vNewPieces[nPieceN].vPieceLoop.back().Invert() ;
                    // Aggiungo nuovo pezzo.
                     vNewPieces.emplace_back() ;
                     vNewPieces.back().vPieceLoop.emplace_back( CurInnerLoop) ;
                     vNewPieces.back().nPiecePart = 1 ;
                    // Cerco loop interni a quello appena aggiunto
                     INTVECTOR vSecondLevel ;
                     for ( int nSI = 1 ; nSI < int( vNewPieces[nPieceN].vPieceLoop.size()) ; ++ nSI) {
                        PNTULIST& SecLevLoopList = vNewPieces[nPieceN].vPieceLoop[nSI].GetUPointList() ;
                        for ( auto itSL = SecLevLoopList.begin() ; itSL != SecLevLoopList.end() ; ++ itSL) {
                           Point3d ptPntSL = itSL->first ;
                           if ( IsPointInsidePolyLine( ptPntSL, vNewPieces.back().vPieceLoop.back(), EPS_SMALL)) {
                              vSecondLevel.emplace_back( nSI) ;
                              break ;
                           }
                        }
                     }
                    // Aggiungo i loop al nuovo pezzo e li tolgo dal precedente.
                     for ( int nSI = 0 ; nSI < int( vSecondLevel.size()) ; ++ nSI)
                        vNewPieces.back().vPieceLoop.emplace_back( vNewPieces[nPieceN].vPieceLoop[vSecondLevel[nSI]]) ;
                     for ( int nSI = 0 ; nSI < int( vSecondLevel.size()) ; ++ nSI)
                        vNewPieces[nPieceN].vPieceLoop.erase( vNewPieces[nPieceN].vPieceLoop.begin() + vSecondLevel[nSI]) ;
                     break ;
                  }
               }
            }
         }
      }
     // Aggiungo al loop esterno i punti dei loop interni che lo toccano.
     // Ciclo sui pezzi della faccia.
      for ( int nPieceN = 0 ; nPieceN < nPieceNum ; ++ nPieceN) {
        // Ciclo sui segmenti del loop esterno.
         PNTULIST& ExtLoopList = vNewPieces[nPieceN].vPieceLoop[0].GetUPointList() ;
         auto itSt = ExtLoopList.begin() ; auto itEn = itSt ; ++ itEn ;
         for ( ; itEn != ExtLoopList.end(); ++itSt, ++itEn) {
            Point3d ptSt = itSt->first ;
            Point3d ptEn = itEn->first ;
            Vector3d vtSeg = ptEn - ptSt ;
            double dSegLen = vtSeg.Len() ;
            vtSeg /= dSegLen ;
           // Vettore dei punti dei loop interni che stanno sul segmento del loop esterno
            PNTUVECTOR vPointWithOrder ;
           // Ciclo sui loop interni del pezzo
            int nInnerLoopNum = int( vNewPieces[nPieceN].vPieceLoop.size()) ;
            for ( int nInnLoop = 1 ; nInnLoop < nInnerLoopNum ; ++ nInnLoop) {
              // Ciclo sui punti del loop interno.
               Point3d ptInnPoint ;
               bool bIsFirst = true ;
               bool bContinue =  vNewPieces[nPieceN].vPieceLoop[nInnLoop].GetFirstPoint( ptInnPoint) ;
               while ( bContinue) {                     
                  DistPointLine DistCalculator( ptInnPoint, ptSt, ptEn) ;
                  double dDist ;
                  DistCalculator.GetDist( dDist) ;
                  double dLongPos = ( ptInnPoint - ptSt) * vtSeg ;
                  if ( dDist < 10 * EPS_SMALL  &&  dLongPos > 0.  &&  dLongPos < dSegLen) {
                     POINTU NewPointU ;
                     NewPointU.first = ptInnPoint ;
                     NewPointU.second = dLongPos ;
                     if ( ! bIsFirst)
                        vPointWithOrder.emplace_back( NewPointU) ;
                  }
                  bIsFirst = false ;
                  bContinue = vNewPieces[nPieceN].vPieceLoop[nInnLoop].GetNextPoint( ptInnPoint) ;
               }  
            }
           // Riordino i punti interni sul segmento esterno in funzione della distanza dall'origine di esso
            for ( int nPi = 0 ; nPi < int( vPointWithOrder.size()) - 1 ; ++ nPi) {
               for ( int nPj = nPi + 1 ; nPj < int( vPointWithOrder.size()) ; ++ nPj) {
                  if ( vPointWithOrder[nPi].second > vPointWithOrder[nPj].second) {
                     swap( vPointWithOrder[nPi], vPointWithOrder[nPj]) ;
                  }
               }
            }
           // Aggiungo i punti al loop esterno
            for ( int nPi = 0 ; nPi < int( vPointWithOrder.size()) ; ++ nPi) {
               itSt = ExtLoopList.emplace( itEn, vPointWithOrder[nPi]) ;
            }
         }       
      }
     // Se non ho diviso la faccia, la elimino
      if ( int( vNewPieces.size()) == 1)
         NewFacet.erase( it->first) ;
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
SurfTriMesh::EdgeInteriorContactManager( const SurfTriMesh& OthSurf,
                                         const INTERSCHAINMAP& InterInterLineMap,
                                         const INTERSEDGEMAP& EdgeInterLineMap)
{
   // Ciclo su tutte le facce.
   for ( auto it = EdgeInterLineMap.begin() ; it != EdgeInterLineMap.end() ; ++ it) {
      int nFacet = it->first ;
      const vector<IntersEdge>& vEdgeContact = it->second ;
      // Triangoli della faccia
      INTVECTOR vFacetTria ;
      GetAllTriaInFacet( nFacet, vFacetTria) ;
      // Versore normale alla faccia.
      Vector3d vtN = m_vTria[m_vFacet[nFacet]].vtN ;
      // Se la faccia è stata tagliata, salto alla prossima faccia.
      if ( InterInterLineMap.find( nFacet) != InterInterLineMap.end())
         continue ;
      int nInOutPart = 0 ;
      // Ciclo sui contatti edge-interno della faccia.
      for ( const IntersEdge& Edge : vEdgeContact) {
         // Estremi del segmento di contatto
         Point3d ptSegSt = Edge.ptSt ;
         Point3d ptSegEn = Edge.ptEn ;
         Vector3d vtSegDir = ptSegEn - ptSegSt ;
         double dLen = vtSegDir.Len() ;
         if ( dLen < EPS_SMALL)
            continue ;
         vtSegDir /= dLen ;
         // Cerco il tratto della frontiera della faccia su cui giace il segmento di contatto.
         Vector3d vtLoopSeg ;
         POLYLINEVECTOR vFacetLoop ;
         GetFacetLoops( nFacet, vFacetLoop) ;
         for ( int nL = 0 ; nL < int( vFacetLoop.size())  &&  nInOutPart == 0 ; ++ nL) {
            Point3d ptLs, ptLe ;
            bool bContinue = vFacetLoop[nL].GetFirstPoint( ptLs)  &&  vFacetLoop[nL].GetNextPoint( ptLe) ;
            while ( bContinue) {
               DistPointLine StartPointLineDistCalc( ptLs, ptSegSt, ptSegEn, false) ;
               double dDistStart ;
               StartPointLineDistCalc.GetDist( dDistStart) ;
               DistPointLine EndPointLineDistCalc( ptLe, ptSegSt, ptSegEn, false) ;
               double dDistEnd ;
               EndPointLineDistCalc.GetDist( dDistEnd) ;
               double dUS = ( ptLs - ptSegSt) * vtSegDir ;
               double dUE = ( ptLe - ptSegSt) * vtSegDir ;
               // Segmento sovrapposto
               if ( dDistStart < EPS_SMALL  &&  dDistEnd < EPS_SMALL  &&  dUS < dLen + EPS_SMALL &&  dUE > - EPS_SMALL) {
                  vtLoopSeg = ptLe - ptLs ;
                  vtLoopSeg.Normalize() ;
                  bContinue = false ;
               }
               else {
                  ptLs = ptLe ;
                  bContinue = vFacetLoop[nL].GetNextPoint( ptLe) ;
               }
            }
         }
         // Determino il valore dell'indice interno/esterno.
         Vector3d vtIn = vtN ^ vtLoopSeg ;
         // Ciclo sulle facce di contatto dell'altra superficie.
         for ( int nOthFacet = 0 ; nOthFacet < int( Edge.vOthFacetIndex.size())  &&  nInOutPart == 0 ; ++ nOthFacet) {
            int nOthFacetIndex = Edge.vOthFacetIndex[nOthFacet] ;
            // Versore normale all'altra faccia.
            Vector3d vtOthN = OthSurf.m_vTria[OthSurf.m_vFacet[nOthFacetIndex]].vtN ;
            // Valore dell'indice interno/esterno.
            nInOutPart = ( vtIn * vtOthN < - EPS_SMALL ? 1 : ( vtIn * vtOthN > EPS_SMALL ? - 1 : 0)) ;
         }
         // Ciclo sui triangoli della faccia
         for ( int& nT : vFacetTria) {
            // Assegno indice interno/esterno al triangolo
            m_vTria[nT].nTempPart = nInOutPart ;
            Triangle3d trTria ;
            GetTriangle( nT, trTria) ;
            trTria.Validate() ;
            // Se un segmento del triangolo si sovrappone a quello di contatto, lo marco come invalicabile.
            for ( int n = 0 ; n < 3 ; ++ n) {
               Point3d ptTS = trTria.GetP( n) ;
               Point3d ptTE = trTria.GetP( ( n + 1) % 3) ;
               DistPointLine StartPointLineDistCalc( ptTS, ptSegSt, ptSegEn, false) ;
               double dDistStart ;
               StartPointLineDistCalc.GetDist( dDistStart) ;
               DistPointLine EndPointLineDistCalc( ptTE, ptSegSt, ptSegEn, false) ;
               double dDistEnd ;
               EndPointLineDistCalc.GetDist( dDistEnd) ;
               double dUS = ( ptTS - ptSegSt) * vtSegDir ;
               double dUE = ( ptTE - ptSegSt) * vtSegDir ;
               // Segmento sovrapposto
               if ( dDistStart < EPS_SMALL  &&  dDistEnd < EPS_SMALL  &&  dUS < dLen + EPS_SMALL  &&  dUE > - EPS_SMALL) {
                  m_vTria[nT].nETempFlag[n] = 1 ;   
               }
            }
         }
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
SurfTriMesh::EdgeEdgeContactManager( const SurfTriMesh& OthSurf,
                                     const INTERSCHAINMAP& InterInterLineMap,
                                     const INTERSEDGEMAP& EdgeEdgeLineMap)
{
   // Ciclo su tutte le facce.
   for ( auto it = EdgeEdgeLineMap.begin() ; it != EdgeEdgeLineMap.end() ; ++ it) {
      int nFacet = it->first ;
      const vector<IntersEdge>& vEdgeContact = it->second ;
      // Triangoli della faccia
      INTVECTOR vFacetTria ;
      GetAllTriaInFacet( nFacet, vFacetTria) ;
      // Versore normale alla faccia.
      Vector3d vtN = m_vTria[m_vFacet[nFacet]].vtN ;
      // Se la faccia è stata tagliata, salto alla prossima faccia.
      if ( InterInterLineMap.find( nFacet) != InterInterLineMap.end())
         continue ;
      int nInOutPart = 0 ;
      // Ciclo sui contatti edge-interno della faccia.
      for ( const IntersEdge& Edge : vEdgeContact) {
         // Estremi del segmento di contatto
         Point3d ptSegSt = Edge.ptSt ;
         Point3d ptSegEn = Edge.ptEn ;
         Vector3d vtSegDir = ptSegEn - ptSegSt ;
         double dLen = vtSegDir.Len() ;
         if ( dLen < EPS_SMALL)
            continue ;
         vtSegDir /= dLen ;
         // Cerco il tratto della frontiera della faccia su cui giace il segmento di contatto.
         Vector3d vtLoopSeg ;
         POLYLINEVECTOR vFacetLoop ;
         GetFacetLoops( nFacet, vFacetLoop) ;
         for ( int nL = 0 ; nL < int( vFacetLoop.size())  &&  nInOutPart == 0 ; ++ nL) {
            Point3d ptLs, ptLe ;
            bool bContinue = vFacetLoop[nL].GetFirstPoint( ptLs)  &&  vFacetLoop[nL].GetNextPoint( ptLe) ;
            while ( bContinue) {
               DistPointLine StartPointLineDistCalc( ptLs, ptSegSt, ptSegEn, false) ;
               double dDistStart ;
               StartPointLineDistCalc.GetDist( dDistStart) ;
               DistPointLine EndPointLineDistCalc( ptLe, ptSegSt, ptSegEn, false) ;
               double dDistEnd ;
               EndPointLineDistCalc.GetDist( dDistEnd) ;
               double dUS = ( ptLs - ptSegSt) * vtSegDir ;
               double dUE = ( ptLe - ptSegSt) * vtSegDir ;
               // Segmento sovrapposto
               if ( dDistStart < EPS_SMALL  &&  dDistEnd < EPS_SMALL  &&  dUS < dLen + EPS_SMALL  &&  dUE > - EPS_SMALL) {
                  vtLoopSeg = ptLe - ptLs ;
                  vtLoopSeg.Normalize() ;
                  bContinue = false ;
               }
               else {
                  ptLs = ptLe ;
                  bContinue = vFacetLoop[nL].GetNextPoint( ptLe) ;
               }
            }
         }
         // Determino il valore dell'indice interno/esterno.
         // Vettore diretto verso l'interno della faccia.
         Vector3d vtIn = vtN ^ vtLoopSeg ;
         // Ciclo sulle facce di contatto dell'altra superficie.
         int nOthFacet1 ;
         int nOthFacet2 = - 1 ;
         for ( nOthFacet1 = 0 ; nOthFacet1 < int( Edge.vOthFacetIndex.size()) - 1  &&  nInOutPart == 0 ; ++ nOthFacet1) {
            int nOthFacetIndex1 = Edge.vOthFacetIndex[nOthFacet1] ;
            Vector3d vtOthN1 = OthSurf.m_vTria[OthSurf.m_vFacet[nOthFacetIndex1]].vtN ;
            for ( nOthFacet2 = nOthFacet1 + 1 ; nOthFacet2 < int( Edge.vOthFacetIndex.size())  &&  nInOutPart == 0 ; ++ nOthFacet2) {
               int nOthFacetIndex2 = Edge.vOthFacetIndex[nOthFacet2] ;
               Vector3d vtOthN2 = OthSurf.m_vTria[OthSurf.m_vFacet[nOthFacetIndex2]].vtN ;
               if ( ! AreSameVectorApprox( vtOthN1, vtOthN2)) {
                  break ;
               }
            }
            if ( nOthFacet2 != int( Edge.vOthFacetIndex.size()))
               break ;
         }
         // Se l'indice interno/esterno non è stato assegnato
         if ( nInOutPart == 0) {
            // Se non ho trovato una coppia di facce, uso la prima.
            if ( ( nOthFacet1 == int( Edge.vOthFacetIndex.size()) - 1  &&  nOthFacet2 == int( Edge.vOthFacetIndex.size()))  ||  nOthFacet2 == - 1) {
               ;
               //int nOthFacetIndex1 = Edge.vOthFacetIndex[0] ;
               //Vector3d vtOthN1 = OthSurf.m_vTria[OthSurf.m_vFacet[nOthFacetIndex1]].vtN ;
               //// Valore dell'indice interno/esterno.
               //nInOutPart = ( vtIn * vtOthN1 < - EPS_SMALL ? 1 : ( vtIn * vtOthN1 > EPS_SMALL ? - 1 : 0)) ;
            }
            // Due facce
            else {
               int nOthFacetIndex1 = Edge.vOthFacetIndex[nOthFacet1] ;
               Vector3d vtOthN1 = OthSurf.m_vTria[OthSurf.m_vFacet[nOthFacetIndex1]].vtN ;
               int nOthFacetIndex2 = Edge.vOthFacetIndex[nOthFacet2] ;
               Vector3d vtOthN2 =  OthSurf.m_vTria[OthSurf.m_vFacet[nOthFacetIndex2]].vtN ;
               POLYLINEVECTOR vOthFacetLoop2 ;
               OthSurf.GetFacetLoops( nOthFacetIndex2, vOthFacetLoop2) ;
               Vector3d vtContactEdgeLoop2 ;
               Point3d ptLs, ptLe ;
               bool bContinue = true ;
               for ( int nL = 0 ; nL < int( vOthFacetLoop2.size())  &&  bContinue ; ++ nL) {
                  bContinue = vOthFacetLoop2[nL].GetFirstPoint( ptLs)  &&  vOthFacetLoop2[nL].GetNextPoint( ptLe) ;
                  while ( bContinue) {
                     DistPointLine StartPointLineDistCalc( ptLs, ptSegSt, ptSegEn, false) ;
                     double dDistStart ;
                     StartPointLineDistCalc.GetDist( dDistStart) ;
                     DistPointLine EndPointLineDistCalc( ptLe, ptSegSt, ptSegEn, false) ;
                     double dDistEnd ;
                     EndPointLineDistCalc.GetDist( dDistEnd) ;
                     double dUS = ( ptLs - ptSegSt) * vtSegDir ;
                     double dUE = ( ptLe - ptSegSt) * vtSegDir ;
                     // Segmento sovrapposto
                     if ( dDistStart < EPS_SMALL  &&  dDistEnd < EPS_SMALL  &&  dUS < dLen + EPS_SMALL  &&  dUE > - EPS_SMALL) {
                        vtContactEdgeLoop2 = ptLe - ptLs ;
                        vtContactEdgeLoop2.Normalize() ;
                        bContinue = false ;
                     }
                     else {
                        ptLs = ptLe ;
                        bContinue = vOthFacetLoop2[nL].GetNextPoint( ptLe) ;
                     }
                  }
               }
               Vector3d vtInFacet2 = vtN ^ vtContactEdgeLoop2 ;
               // Convesso
               if ( vtInFacet2 * vtOthN1 < 0.) {
                  nInOutPart = vtIn * vtOthN1 < 0  &&  vtIn * vtOthN2 < 0 ? 1 : - 1 ; 
               }
               // Concavo
               else {
                  nInOutPart = vtIn * vtOthN1 > 0  &&  vtIn * vtOthN2 > 0 ? - 1 : 1 ;
               }
            }
         }
         // Ciclo sui triangoli della faccia
         for ( int& nT : vFacetTria) {
            // Assegno indice interno/esterno al triangolo
            m_vTria[nT].nTempPart = nInOutPart ;
            Triangle3d trTria ;
            GetTriangle( nT, trTria) ;
            trTria.Validate() ;
            // Se un segmento del triangolo si sovrappone a quello di contatto, lo marco come invalicabile.
            for ( int n = 0 ; n < 3 ; ++ n) {
               Point3d ptTS = trTria.GetP( n) ;
               Point3d ptTE = trTria.GetP( ( n + 1) % 3) ;
               DistPointLine StartPointLineDistCalc( ptTS, ptSegSt, ptSegEn, false) ;
               double dDistStart ;
               StartPointLineDistCalc.GetDist( dDistStart) ;
               DistPointLine EndPointLineDistCalc( ptTE, ptSegSt, ptSegEn, false) ;
               double dDistEnd ;
               EndPointLineDistCalc.GetDist( dDistEnd) ;
               double dUS = ( ptTS - ptSegSt) * vtSegDir ;
               double dUE = ( ptTE - ptSegSt) * vtSegDir ;
               // Segmento sovrapposto
               if ( dDistStart < EPS_SMALL  &&  dDistEnd < EPS_SMALL  &&  dUS < dLen + EPS_SMALL  &&  dUE > - EPS_SMALL) {
                  m_vTria[nT].nETempFlag[n] = 1 ;   
               }
            }
         }
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
SurfTriMesh::RetriangulateFacetPieces( const PieceMap& NewFacet,
                                       const INTERSEDGEMAP& EdgeInterLineMap,
                                       const INTERSEDGEMAP& EdgeEdgeLineMap)
{
   // Vettore degli indici delle facce
   INTVECTOR vFacetIndexes ;
   // Ciclo sulle facce.
   for ( auto it = NewFacet.begin() ; it != NewFacet.end() ; ++ it) {
      // Accedo al colore di un triangolo della faccia e cancello i triangoli della faccia.
      INTVECTOR vFacetTria ;
      GetAllTriaInFacet( it->first, vFacetTria) ;
      for ( int& nT : vFacetTria)
         RemoveTriangle( nT) ;
   }
   for ( auto it = NewFacet.begin() ; it != NewFacet.end() ; ++ it) {
      const vector<FacetPiece>& PiecesVector = it->second ;
      // Ciclo sui nuovi pezzi di faccia.
      int nPartNum = int( PiecesVector.size()) ;
      for ( int nPart = 0 ; nPart < nPartNum ; ++ nPart) {
         PNTVECTOR vPt ;
         INTVECTOR vTr ;
         if ( Triangulate().Make( PiecesVector[nPart].vPieceLoop, 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]])} ;
               // Il colore passato nel secondo parametro è definito perché il vettore delle facce ha ancora salvato l'indice di un suo triangolo e
               // il triangolo è cancellato semplicente assegnadno la costante apposita VT_DEL (-2) a nIdVert[0]; ma il suo colore resta definito. 
               int nNewTriaNum = AddTriangle( nNewId, m_vTria[m_vFacet[it->first]].nTFlag) ;
               if ( IsValidSvt( nNewTriaNum)) {
                  // Assegno l'indice della parte interno/esterno.
                  m_vTria[nNewTriaNum].nTempPart = PiecesVector[nPart].nPiecePart ;
                  // Indici parte interno/esterno invalicabili
                  m_vTria[nNewTriaNum].nETempFlag[0] = 0 ;
                  m_vTria[nNewTriaNum].nETempFlag[1] = 0 ;
                  m_vTria[nNewTriaNum].nETempFlag[2] = 0 ;
                  // Se la faccia ha avuto un contatto edge-interior, cerco edge invalicabili.
                  auto itEdgeInterEdge = EdgeInterLineMap.find( it->first) ;
                  auto itEdgeEdgeEdge = EdgeEdgeLineMap.find( it->first) ;
                  if ( itEdgeInterEdge != EdgeInterLineMap.end()) {
                     const vector<IntersEdge>& vEdgeContact = itEdgeInterEdge->second ;
                     Triangle3d trTria ;
                     GetTriangle( nNewTriaNum, trTria) ;
                     // Ciclo sui segmenti del triangolo.
                     for ( int n = 0 ; n < 3 ; ++ n) {
                        Point3d ptTS = trTria.GetP( n) ;
                        Point3d ptTE = trTria.GetP( ( n + 1) % 3) ;
                        for (int m = 0; m < int(vEdgeContact.size()) ; ++ m) {
                           Point3d ptEdgeSt = vEdgeContact[m].ptSt ;
                           Point3d ptEdgeEn = vEdgeContact[m].ptEn ;
                           Vector3d vtEdgeDir = ptEdgeEn - ptEdgeSt ;
                           double dLen = vtEdgeDir.Len() ;
                           if ( dLen < EPS_SMALL)
                              continue ; 
                           DistPointLine StartPointLineDistCalc( ptTS, ptEdgeSt, ptEdgeEn, false) ;
                           double dDistStart ;
                           StartPointLineDistCalc.GetDist( dDistStart) ;
                           DistPointLine EndPointLineDistCalc( ptTE, ptEdgeSt, ptEdgeEn, false) ;
                           double dDistEnd ;
                           EndPointLineDistCalc.GetDist( dDistEnd) ;
                           double dUS = ( ptTS - ptEdgeSt) * vtEdgeDir ;
                           double dUE = ( ptTE - ptEdgeSt) * vtEdgeDir ;
                           // Segmento sovrapposto
                           if ( dDistStart < EPS_SMALL  &&  dDistEnd < EPS_SMALL  &&  dUS < dLen  &&  dUE > 0) {
                              m_vTria[nNewTriaNum].nETempFlag[n] = 1 ;
                              break ;
                           }
                        }
                     }
                  }
                  else if ( itEdgeEdgeEdge != EdgeEdgeLineMap.end()) {
                     const vector<IntersEdge>& vEdgeContact = itEdgeEdgeEdge->second ;
                     Triangle3d trTria ;
                     GetTriangle( nNewTriaNum, trTria) ;
                     // Ciclo sui segmenti del triangolo.
                     for ( int n = 0 ; n < 3 ; ++ n) {
                        Point3d ptTS = trTria.GetP( n) ;
                        Point3d ptTE = trTria.GetP( ( n + 1) % 3) ;
                        for (int m = 0; m < int(vEdgeContact.size()) ; ++ m) {
                           Point3d ptEdgeSt = vEdgeContact[m].ptSt ;
                           Point3d ptEdgeEn = vEdgeContact[m].ptEn ;
                           Vector3d vtEdgeDir = ptEdgeEn - ptEdgeSt ;
                           double dLen = vtEdgeDir.Len() ;
                           if ( dLen < EPS_SMALL)
                              continue ; 
                           DistPointLine StartPointLineDistCalc( ptTS, ptEdgeSt, ptEdgeEn, false) ;
                           double dDistStart ;
                           StartPointLineDistCalc.GetDist( dDistStart) ;
                           DistPointLine EndPointLineDistCalc( ptTE, ptEdgeSt, ptEdgeEn, false) ;
                           double dDistEnd ;
                           EndPointLineDistCalc.GetDist( dDistEnd) ;
                           double dUS = ( ptTS - ptEdgeSt) * vtEdgeDir ;
                           double dUE = ( ptTE - ptEdgeSt) * vtEdgeDir ;
                           // Segmento sovrapposto
                           if ( dDistStart < EPS_SMALL  &&  dDistEnd < EPS_SMALL  &&  dUS < dLen  &&  dUE > 0) {
                              m_vTria[nNewTriaNum].nETempFlag[n] = 1 ;
                              break ;
                           }
                        }
                     }
                  }
               }
            }
         }
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
SurfTriMesh::IntersectTriMeshFacets( SurfTriMesh& Other)
{
   ////////////////////////////////////////////////////////////////////////////////////////////////////////
   static int nTime = 0 ;
   nTime ++ ;
   ////////////////////////////////////////////////////////////////////////////////////////////////////////
   SurfTriMesh& SurfB = Other ;
   // Le superfici devono essere valide
   if ( m_nStatus != OK  ||  ! SurfB.IsValid())
      return false ;
   // Unordered map dei segmenti di intersezione (contatti interno-interno) e dei contatti con edge
   INTERSCHAINMAP IntersLineMapA, IntersLineMapB ;
   INTERSEDGEMAP EdgeInnLineMapA, EdgeInnLineMapB ;
   INTERSEDGEMAP EdgeEdgeLineMapA, EdgeEdgeLineMapB ;
   // Setto il triangolo come né fuori né dentro
   int nTriaNumA = GetTriangleSize();
   int nTriaNumB = SurfB.GetTriangleSize();
   for ( int nTA = 0 ; nTA < nTriaNumA ; ++ nTA) {
      m_vTria[nTA].nTempPart = 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].nTempPart = 0 ;
      SurfB.m_vTria[nTB].nETempFlag[0] = 0 ;
      SurfB.m_vTria[nTB].nETempFlag[1] = 0 ;
      SurfB.m_vTria[nTB].nETempFlag[2] = 0 ;
   }
   // Ciclo sulle facce delle mesh
   int nFacetNumA = GetFacetCount() ;
   for ( int nFA = 0 ; nFA < nFacetNumA ; ++ nFA) {
      // Dati della faccia
      POLYLINEVECTOR LoopVecA ;
      GetFacetLoops( nFA, LoopVecA) ;
      PtrOwner<SurfFlatRegion> pRegA( GetBasicSurfFlatRegion( GetSurfFlatRegionFromPolyLineVector( LoopVecA))) ;
      if ( pRegA == nullptr)
         return false ;
      SurfFlatRegion RegionA = *pRegA ;
      // Recupero tutti i triangoli della superficie B che cadono nel box della faccia di A.
      BBox3d b3BoxA ;
      RegionA.GetLocalBBox( b3BoxA) ;
      INTVECTOR vNearTria ;
      SurfB.GetAllTriaOverlapBox( b3BoxA, vNearTria) ;
      // Determino le facce a cui appartengono i tirangoli che cadono nel box
      INTVECTOR vFacetIndexesB ;
      for ( int& nT : vNearTria) {
         int nF = SurfB.GetFacetFromTria(nT) ;
         int n ;
         for ( n = 0 ; n < int( vFacetIndexesB.size()) ; ++ n) {
            if ( nF == vFacetIndexesB[n])
               break ;
         }
         if ( n == int(vFacetIndexesB.size()))
            vFacetIndexesB.emplace_back( nF) ;
      }
      for ( int& nFB : vFacetIndexesB) {
         // Dati della faccia
         POLYLINEVECTOR LoopVecB ;
         SurfB.GetFacetLoops( nFB, LoopVecB) ;
         PtrOwner<SurfFlatRegion> pRegB( GetBasicSurfFlatRegion( GetSurfFlatRegionFromPolyLineVector( LoopVecB))) ;
         if ( pRegB == nullptr)
            return false ;
         SurfFlatRegion RegionB = *pRegB ;
         // Interseco le due facce
         LineFacetClassVector CommonIntersection ;
         // Intersezione fra le facce
         if ( IntersFacetFacet( RegionA, LoopVecA[0], RegionB, LoopVecB[0], CommonIntersection)) {
            // Ciclo sulle parti dell'intersezione
            for ( int nPart = 0 ; nPart < int( CommonIntersection.size()) ; ++ nPart) {
               // Intersezione nell'interno di entrambe le facce
               if ( CommonIntersection[nPart].nTypeA == CRVC_IN  &&
                    CommonIntersection[nPart].nTypeB == CRVC_IN) {
                  // Salvo intersezione per la faccia A
                  auto itA = IntersLineMapA.find( nFA) ;
                  if ( itA != IntersLineMapA.end()) {
                     itA->second.emplace_back( IntersInnSeg( CommonIntersection[nPart].ptSt, CommonIntersection[nPart].ptEn)) ;
                  }
                  else {
                     IntersLineMapA.emplace( nFA, IntersInnChain( 1, IntersInnSeg( CommonIntersection[nPart].ptSt, CommonIntersection[nPart].ptEn))) ;
                  }
                  // Salvo intersezione per la faccia B
                  auto itB = IntersLineMapB.find( nFB) ;
                  if ( itB != IntersLineMapB.end()) {
                     itB->second.emplace_back( IntersInnSeg( CommonIntersection[nPart].ptEn, CommonIntersection[nPart].ptSt)) ;
                  }
                  else {
                     IntersLineMapB.emplace( nFB, IntersInnChain( 1, IntersInnSeg( CommonIntersection[nPart].ptEn, CommonIntersection[nPart].ptSt))) ;
                  }
               }
               // Intersezione all'interno della faccia A
               else if ( CommonIntersection[nPart].nTypeA == CRVC_IN) {
                  // Salvo intersezione per la faccia A
                  auto itA = IntersLineMapA.find( nFA) ;
                  if ( itA != IntersLineMapA.end()) {
                     itA->second.emplace_back( IntersInnSeg( CommonIntersection[nPart].ptSt, CommonIntersection[nPart].ptEn)) ;
                  }
                  else {
                     IntersLineMapA.emplace( nFA, IntersInnChain( 1, IntersInnSeg( CommonIntersection[nPart].ptSt, CommonIntersection[nPart].ptEn))) ;
                  }
                  // Salvo intersezione per la faccia B
                  auto itB = EdgeInnLineMapB.find( nFB) ;
                  if ( itB != EdgeInnLineMapB.end()) {
                     int nE ;
                     for ( nE = 0 ; nE < int( itB->second.size()) ; ++ nE) {
                        if ( AreSamePointExact( itB->second[nE].ptSt, CommonIntersection[nPart].ptEn)  &&
                             AreSamePointExact( itB->second[nE].ptEn, CommonIntersection[nPart].ptSt)) {
                           itB->second[nE].vOthFacetIndex.emplace_back( nFA) ;
                           break ;
                        }
                     }
                     if ( nE == int( itB->second.size())) {
                        itB->second.emplace_back( IntersEdge( CommonIntersection[nPart].ptEn, CommonIntersection[nPart].ptSt, INTVECTOR( 1, nFA))) ;
                     }
                  }
                  else {
                     EdgeInnLineMapB.emplace( nFB, IntersEdgeVec( 1, IntersEdge( CommonIntersection[nPart].ptEn, CommonIntersection[nPart].ptSt,
                                                                                 INTVECTOR( 1, nFA)))) ;
                  }
               }
               // Intersezione all'interno della faccia B
               else  if ( CommonIntersection[nPart].nTypeB == CRVC_IN) {
                  // Salvo intersezione per la faccia A
                  auto itA = EdgeInnLineMapA.find( nFA) ;
                  if ( itA != EdgeInnLineMapA.end()) {
                     int nE ;
                     for ( nE = 0 ; nE < int( itA->second.size()) ; ++ nE) {
                        if ( AreSamePointExact( itA->second[nE].ptSt, CommonIntersection[nPart].ptSt)  &&
                             AreSamePointExact( itA->second[nE].ptEn, CommonIntersection[nPart].ptEn)) {
                           itA->second[nE].vOthFacetIndex.emplace_back( nFB) ;
                           break ;
                        }
                     }
                     if ( nE == int( itA->second.size())) {
                        itA->second.emplace_back( IntersEdge( CommonIntersection[nPart].ptSt, CommonIntersection[nPart].ptEn, INTVECTOR( 1, nFB))) ;  
                     }
                  }
                  else {
                     EdgeInnLineMapA.emplace( nFA, IntersEdgeVec( 1, IntersEdge( CommonIntersection[nPart].ptSt, CommonIntersection[nPart].ptEn,
                                                                                 INTVECTOR( 1, nFB)))) ;
                  }
                  // Salvo intersezione per la faccia B
                  auto itB = IntersLineMapB.find( nFB) ;
                  if ( itB != IntersLineMapB.end()) {
                     itB->second.emplace_back( IntersInnSeg( CommonIntersection[nPart].ptEn, CommonIntersection[nPart].ptSt)) ;
                  }
                  else {
                     IntersLineMapB.emplace( nFB, IntersInnChain( 1, IntersInnSeg( CommonIntersection[nPart].ptEn, CommonIntersection[nPart].ptSt))) ;
                  }
               }
               // Intersezione sulla frontiera di entrambe le facce
               else {
                  // Salvo intersezione per la faccia A
                  auto itA = EdgeEdgeLineMapA.find( nFA) ;
                  if ( itA != EdgeEdgeLineMapA.end()) {
                     int nE ;
                     for ( nE = 0 ; nE < int( itA->second.size()) ; ++ nE) {
                        if ( ( AreSamePointExact( itA->second[nE].ptSt, CommonIntersection[nPart].ptSt)  &&
                               AreSamePointExact( itA->second[nE].ptEn, CommonIntersection[nPart].ptEn))  ||
                             ( AreSamePointExact( itA->second[nE].ptSt, CommonIntersection[nPart].ptEn)  &&
                               AreSamePointExact( itA->second[nE].ptEn, CommonIntersection[nPart].ptSt))) {
                           itA->second[nE].vOthFacetIndex.emplace_back( nFB) ;
                           break ;
                        }
                     }
                     if ( nE == int( itA->second.size())) {
                        itA->second.emplace_back( IntersEdge( CommonIntersection[nPart].ptSt, CommonIntersection[nPart].ptEn, INTVECTOR( 1, nFB)));
                     }
                  }
                  else {
                     EdgeEdgeLineMapA.emplace( nFA, IntersEdgeVec( 1, IntersEdge( CommonIntersection[nPart].ptSt, CommonIntersection[nPart].ptEn,
                                                                                 INTVECTOR( 1, nFB)))) ;
                  }
                  // Salvo intersezione per la faccia B
                  auto itB = EdgeEdgeLineMapB.find( nFB) ;
                  if ( itB != EdgeEdgeLineMapB.end()) {
                     int nE ;
                     for ( nE = 0 ; nE < int( itB->second.size()) ; ++ nE) {
                        if ( ( AreSamePointExact( itB->second[nE].ptSt, CommonIntersection[nPart].ptEn)  &&
                               AreSamePointExact( itB->second[nE].ptEn, CommonIntersection[nPart].ptSt))  ||
                             ( AreSamePointExact( itB->second[nE].ptSt, CommonIntersection[nPart].ptSt)  &&
                               AreSamePointExact( itB->second[nE].ptEn, CommonIntersection[nPart].ptEn))) {
                           itB->second[nE].vOthFacetIndex.emplace_back( nFA) ;
                           break ;
                        }
                     }
                     if ( nE == int( itB->second.size())) {
                        itB->second.emplace_back( IntersEdge( CommonIntersection[nPart].ptEn, CommonIntersection[nPart].ptSt, INTVECTOR( 1, nFA))) ;
                     }
                  }
                  else {
                     EdgeEdgeLineMapB.emplace( nFB, IntersEdgeVec( 1, IntersEdge( CommonIntersection[nPart].ptEn, CommonIntersection[nPart].ptSt,
                                                                                 INTVECTOR( 1, nFA)))) ;
                  }
               }
            }
         }
      }  
   }

   // Divido le facce
   PieceMap NewFacetA, NewFacetB ;
   SplitFacet( IntersLineMapA, NewFacetA) ;
   SurfB.SplitFacet( IntersLineMapB, NewFacetB) ;  

   // Facce EdgeInn
   EdgeInteriorContactManager( SurfB, IntersLineMapA, EdgeInnLineMapA) ;
   SurfB.EdgeInteriorContactManager( *this, IntersLineMapB, EdgeInnLineMapB) ;

   // Facce Edge Edge
   EdgeEdgeContactManager( SurfB, IntersLineMapA, EdgeEdgeLineMapA) ;
   SurfB.EdgeEdgeContactManager( *this, IntersLineMapB, EdgeEdgeLineMapB) ;

   // Se non ci sono stati tagli, valuto se una è tutta interna all'altra.
   if ( int( NewFacetA.size() + NewFacetB.size()) == 0) {
      int nVertNum = 0 ;
      Point3d ptVert ;
      int nCurVert = GetFirstVertex( ptVert) ;
      int nInOutNumA = 0 ;
      while ( nInOutNumA == 0  &&  nCurVert != SVT_NULL) {
         int nFacetNum = - 1 ;
         double dMinDist = DBL_MAX ;
         for ( int nFB = 0 ; nFB < SurfB.GetFacetCount() ; ++ nFB) {
            // Loop della faccia della superficie B
            POLYLINEVECTOR vFacetLoopVec ;
            SurfB.GetFacetLoops( nFB, vFacetLoopVec) ;
            double dDist ;
            if ( DistPointFacet( ptVert, vFacetLoopVec, dDist)  &&  dDist < dMinDist) {
               dMinDist = dDist ;
               nFacetNum = nFB ;
            }
         }
         if ( nFacetNum >= 0) {
            Triangle3d trTriaB ;
            SurfB.GetTriangle( SurfB.m_vFacet[nFacetNum], trTriaB) ;
            trTriaB.Validate() ;
            if ( ( ptVert - trTriaB.GetP( 0)) * trTriaB.GetN() < - EPS_SMALL)
               nInOutNumA = 1 ;
            else if ( ( ptVert - trTriaB.GetP( 0)) * trTriaB.GetN() > EPS_SMALL)
               nInOutNumA = - 1 ;
         }
         if ( nInOutNumA == 0) {
            nCurVert = GetNextVertex( nVertNum, ptVert) ;
            ++ nVertNum ;
         }
      }
      // Se la superficie A non è interna valuto la posizione della superficie B.
      int nInOutNumB = nInOutNumA == 1 ? - 1 : 0 ;
      if ( nInOutNumB == 0) {
         nVertNum = 0 ;
         nCurVert = SurfB.GetFirstVertex( ptVert) ;
         while ( nInOutNumB == 0  &&  nCurVert != SVT_NULL) {
            int nFacetNum = - 1 ;
            double dMinDist = DBL_MAX ;
            for ( int nFA = 0 ; nFA < GetFacetCount() ; ++ nFA) {
               // Loop della faccia della superficie A
               POLYLINEVECTOR vFacetLoopVec ;
               GetFacetLoops( nFA, vFacetLoopVec) ;
               double dDist ;
               if ( DistPointFacet( ptVert, vFacetLoopVec, dDist)  &&  dDist < dMinDist) {
                  dMinDist = dDist ;
                  nFacetNum = nFA ;
               }
            }
            if ( nFacetNum >= 0) {
               Triangle3d trTriaA ;
               GetTriangle( m_vFacet[nFacetNum], trTriaA) ;
               trTriaA.Validate() ;
               if ( ( ptVert - trTriaA.GetP( 0)) * trTriaA.GetN() < - EPS_SMALL)
                  nInOutNumB = 1 ;
               else if ( ( ptVert - trTriaA.GetP( 0)) * trTriaA.GetN() > EPS_SMALL)
                  nInOutNumB = - 1 ;

            }
            if ( nInOutNumB == 0) {
               nCurVert = SurfB.GetNextVertex( nVertNum, ptVert) ;
               ++ nVertNum ;
            }
         }
      }
      // Se la posizione della superficie A rispetto alla B non è definita e
      // B è esterna ad A, allora assumiamo che siano reciprocamente esterne.
      if ( nInOutNumA == 0  &&  nInOutNumB == - 1)
         nInOutNumA = - 1 ;
      // Assegno gli indici interni/esterni.
      for ( int nTA = 0 ; nTA < nTriaNumA ; ++ nTA) {
         m_vTria[nTA].nTempPart = nInOutNumA ;
      }
      for ( int nTB = 0 ; nTB < nTriaNumB ; ++ nTB) {
         SurfB.m_vTria[nTB].nTempPart = nInOutNumB ;
      }
   }

#define UseTria 0
#if ! UseTria
   // Gestisco facce sovrapposte.
   // Ciclo sulle facce di A.
   for ( int nFA = 0 ; nFA < nFacetNumA ; ++ nFA) {
      auto itFacetA = NewFacetA.find( nFA) ;
      // Faccia A divisa in pezzi.
      if ( itFacetA != NewFacetA.end()) {
         // Ciclo sui pezzi della faccia A.
         for ( int nPieceNumA = 0 ; nPieceNumA < int( itFacetA->second.size()) ; ++ nPieceNumA) {
            // Dati del pezzo di faccia di A
            POLYLINEVECTOR& vPieceLoopsVecA = itFacetA->second[nPieceNumA].vPieceLoop ;
            PtrOwner<SurfFlatRegion> pRegA( GetBasicSurfFlatRegion( GetSurfFlatRegionFromPolyLineVector( vPieceLoopsVecA))) ;
            if ( pRegA == nullptr)
               return false ;
            SurfFlatRegion RegionA = *pRegA ;
            Plane3d plPlaneA ;
            plPlaneA.Set( RegionA.GetPlanePoint(), RegionA.GetNormVersor()) ;
            // Recupero tutti i triangoli della superficie B che cadono nel box del pezzo di faccia di A.
            BBox3d b3BoxA ;
            RegionA.GetLocalBBox( b3BoxA) ;
            INTVECTOR vNearTria ;
            SurfB.GetAllTriaOverlapBox( b3BoxA, vNearTria) ;
            // Determino le facce a cui appartengono i tirangoli che cadono nel box
            INTVECTOR vFacetIndexesB ;
            for ( int& nT : vNearTria) {
               int nF = SurfB.GetFacetFromTria( nT) ;
               int n ;
               for ( n = 0 ; n < int( vFacetIndexesB.size()) ; ++ n) {
                  if ( nF == vFacetIndexesB[n])
                     break ;
               }
               if ( n == int( vFacetIndexesB.size()))
                  vFacetIndexesB.emplace_back( nF) ;
            }
            // Ciclo sulle facce di B
            for ( int& nFB : vFacetIndexesB) {
               auto itFacetB = NewFacetB.find( nFB) ;
               // Faccia B divisa in pezzi
               if ( itFacetB != NewFacetB.end()) {
                  // Ciclo sui pezzi della faccia B.
                  for ( int nPieceNumB = 0 ; nPieceNumB < int( itFacetB->second.size()) ; ++ nPieceNumB) {
                     // Dati del pezzo di faccia di B
                     POLYLINEVECTOR& vPieceLoopsVecB = itFacetB->second[nPieceNumB].vPieceLoop ;
                     PtrOwner<SurfFlatRegion> pRegB( GetBasicSurfFlatRegion( GetSurfFlatRegionFromPolyLineVector( vPieceLoopsVecB))) ;
                     if ( pRegB == nullptr)
                        return false ;
                     SurfFlatRegion RegionB = *pRegB ;
                     Plane3d plPlaneB ;
                     plPlaneB.Set( RegionB.GetPlanePoint(), RegionB.GetNormVersor()) ;
                     // Intersezione tra i piani delle due facce: se non esiste esse sono complanari.
                     Point3d ptL ; Vector3d vtL ;
                     int nResPP = IntersPlanePlane( plPlaneA, plPlaneB, ptL, vtL) ;
                     if ( /*nResPP == IPPT_NO ||*/ nResPP == IPPT_OVERLAPS) {
                        SurfFlatRegion frIntersRegion = RegionA ;
                        frIntersRegion.Intersect( RegionB) ;
                        double dIntersRegArea ;
                        frIntersRegion.GetArea( dIntersRegArea) ;
                        // Assegno indice di parte ai pezzi.
                        if ( dIntersRegArea > EPS_SMALL) {
                           bool bCodirectedNorm = RegionA.GetNormVersor() * RegionB.GetNormVersor() > 0;
                           int nPartIndex = bCodirectedNorm ? 2 : - 2 ;
                           itFacetA->second[nPieceNumA].nPiecePart = nPartIndex;
                           itFacetB->second[nPieceNumB].nPiecePart = nPartIndex;
                        }
                     }
                  }   
               }
               // Faccia B non divisa in pezzi
               else {
                  // Dati della faccia
                  POLYLINEVECTOR vFacetLoopVecB ;
                  SurfB.GetFacetLoops( nFB, vFacetLoopVecB) ;
                  PtrOwner<SurfFlatRegion> pRegB( GetBasicSurfFlatRegion( GetSurfFlatRegionFromPolyLineVector( vFacetLoopVecB))) ;
                  if ( pRegB == nullptr)
                     return false ;
                  SurfFlatRegion RegionB = *pRegB ;
                  Plane3d plPlaneB ;
                  plPlaneB.Set( RegionB.GetPlanePoint(), RegionB.GetNormVersor()) ;
                  // Intersezione tra i piani delle due facce: se non esiste esse sono complanari.
                  Point3d ptL ; Vector3d vtL ;
                  int nResPP = IntersPlanePlane( plPlaneA, plPlaneB, ptL, vtL) ;
                  if ( /*nResPP == IPPT_NO  ||*/  nResPP == IPPT_OVERLAPS) {
                     SurfFlatRegion frIntersRegion = RegionA ;
                     frIntersRegion.Intersect( RegionB) ;
                     double dIntersRegArea ;
                     frIntersRegion.GetArea( dIntersRegArea) ;
                     // Assegno indice di parte al pezzo di faccia A e ai triangoli della faccia B.
                     if ( dIntersRegArea > EPS_SMALL) {
                        bool bCodirectedNorm = RegionA.GetNormVersor() * RegionB.GetNormVersor() > 0 ;
                        int nPartIndex = bCodirectedNorm ? 2 : - 2 ;
                        itFacetA->second[nPieceNumA].nPiecePart = nPartIndex ;
                        INTVECTOR vTB ;
                        SurfB.GetAllTriaInFacet( nFB, vTB) ;
                        for ( int& nTB : vTB)
                           SurfB.m_vTria[nTB].nTempPart = nPartIndex ;
                     }
                  }
               }
            }
         }  
      }
      // Faccia A non divisa in pezzi
      else {
         // Dati della faccia
         POLYLINEVECTOR vFacetLoopsVecA ;
         GetFacetLoops( nFA, vFacetLoopsVecA) ;
         PtrOwner<SurfFlatRegion> pRegA( GetBasicSurfFlatRegion( GetSurfFlatRegionFromPolyLineVector( vFacetLoopsVecA))) ;
         if ( pRegA == nullptr)
            return false ;
         SurfFlatRegion RegionA = *pRegA ;
         Plane3d plPlaneA ;
         plPlaneA.Set( RegionA.GetPlanePoint(), RegionA.GetNormVersor()) ;
         // Recupero tutti i triangoli della superficie B che cadono nel box del pezzo di faccia di A.
         BBox3d b3BoxA ;
         RegionA.GetLocalBBox( b3BoxA) ;
         INTVECTOR vNearTria ;
         SurfB.GetAllTriaOverlapBox( b3BoxA, vNearTria) ;
         // Determino le facce a cui appartengono i tirangoli che cadono nel box
         INTVECTOR vFacetIndexesB ;
         for ( int& nT : vNearTria) {
            int nF = SurfB.GetFacetFromTria( nT) ;
            int n ;
            for ( n = 0 ; n < int( vFacetIndexesB.size()) ; ++ n) {
               if ( nF == vFacetIndexesB[n])
                  break ;
            }
            if ( n == int( vFacetIndexesB.size()))
               vFacetIndexesB.emplace_back( nF) ;
         }
         // Ciclo sulle facce di B
         for ( int& nFB : vFacetIndexesB) {
            auto itFacetB = NewFacetB.find( nFB) ;
            // Faccia B divisa in pezzi
            if ( itFacetB != NewFacetB.end()) {
               for ( int nPieceNumB = 0 ; nPieceNumB < int( itFacetB->second.size()) ; ++ nPieceNumB) {
                  POLYLINEVECTOR& vPieceLoopsVecB = itFacetB->second[nPieceNumB].vPieceLoop ;
                  PtrOwner<SurfFlatRegion> pRegB( GetBasicSurfFlatRegion( GetSurfFlatRegionFromPolyLineVector( vPieceLoopsVecB))) ;
                  if ( pRegB == nullptr)
                     return false ;
                  SurfFlatRegion RegionB = *pRegB ;
                  Plane3d plPlaneB ;
                  plPlaneB.Set( RegionB.GetPlanePoint(), RegionB.GetNormVersor()) ;
                  // Intersezione tra i piani delle due facce: se non esiste esse sono complanari.
                  Point3d ptL ; Vector3d vtL ;
                  int nResPP = IntersPlanePlane( plPlaneA, plPlaneB, ptL, vtL) ;
                  if ( /*nResPP == IPPT_NO  ||*/  nResPP == IPPT_OVERLAPS) {
                     SurfFlatRegion frIntersRegion = RegionA ;
                     frIntersRegion.Intersect( RegionB) ;
                     double dIntersRegArea ;
                     frIntersRegion.GetArea( dIntersRegArea) ;
                     // Assegno indice di parte ai triangoli di A e ai pezzi di B.
                     if ( dIntersRegArea > EPS_SMALL) {
                        bool bCodirectedNorm = RegionA.GetNormVersor() * RegionB.GetNormVersor() > 0 ;
                        int nPartIndex = bCodirectedNorm ? 2 : - 2 ;
                        INTVECTOR vTA ;
                        GetAllTriaInFacet( nFA, vTA) ;
                        for ( int& nTA : vTA)
                           m_vTria[nTA].nTempPart = nPartIndex ;
                        itFacetB->second[nPieceNumB].nPiecePart = nPartIndex ;
                     }
                  }
               }   
            }
            // Faccia B non divisa in pezzi
            else {
               // Dati della faccia
               POLYLINEVECTOR vFacetLoopVecB ;
               SurfB.GetFacetLoops( nFB, vFacetLoopVecB) ;
               PtrOwner<SurfFlatRegion> pRegB( GetBasicSurfFlatRegion( GetSurfFlatRegionFromPolyLineVector( vFacetLoopVecB))) ;
               if ( pRegB == nullptr)
                  return false ;
               SurfFlatRegion RegionB = *pRegB ;
               Plane3d plPlaneB ;
               plPlaneB.Set( RegionB.GetPlanePoint(), RegionB.GetNormVersor()) ;
               // Intersezione tra i piani delle due facce: se non esiste esse sono complanari.
               Point3d ptL ; Vector3d vtL ;
               int nResPP = IntersPlanePlane( plPlaneA, plPlaneB, ptL, vtL) ;
               if ( /*nResPP == IPPT_NO  ||*/  nResPP == IPPT_OVERLAPS) {
                  SurfFlatRegion frIntersRegion = RegionA ;
                  frIntersRegion.Intersect( RegionB) ;
                  double dIntersRegArea ;
                  frIntersRegion.GetArea( dIntersRegArea) ;
                  // Assegno indice di parte al pezzo di faccia A e ai triangoli della faccia B.
                  if ( dIntersRegArea > EPS_SMALL) {
                     bool bCodirectedNorm = RegionA.GetNormVersor() * RegionB.GetNormVersor() > 0 ;
                     int nPartIndex = bCodirectedNorm ? 2 : - 2 ;
                     INTVECTOR vTA ;
                     GetAllTriaInFacet( nFA, vTA) ;
                     for ( int& nTA : vTA)
                        m_vTria[nTA].nTempPart = nPartIndex ;
                     INTVECTOR vTB ;
                     SurfB.GetAllTriaInFacet( nFB, vTB) ;
                     for ( int& nTB : vTB)
                        SurfB.m_vTria[nTB].nTempPart = nPartIndex ;
                  }
               }
            }
         }
      }
   }
#endif
   
   // Elimino i triangoli delle facce e li sostituisco con i nuovi con i loro indici
   RetriangulateFacetPieces( NewFacetA, EdgeInnLineMapA, EdgeEdgeLineMapB) ;
   SurfB.RetriangulateFacetPieces( NewFacetB, EdgeInnLineMapB, EdgeEdgeLineMapB) ;

   bool bOk = ( AdjustVertices()  &&  DoCompacting()  &&  SurfB.AdjustVertices()  &&  SurfB.DoCompacting()) ;
   
#if UseTria
   // Gestione triangoli sovrapposti
   if ( bOk) {
      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 i triangoli sono sovrapposti
            TRIA3DVECTOR vTriaAB ;
            Point3d ptTempA, ptTempB ;
            int nIntTypeAB = IntersTriaTria( trTriaA, trTriaB, ptTempA, ptTempB, vTriaAB) ;
            // Verifica della correttezza del'intersezione
            bool bSuccesfullInters = true ;
            if ( int( vTriaAB.size()) == 1) {
               Point3d ptIntTriaCentroid = vTriaAB[0].GetCentroid() ;
               bSuccesfullInters = IsPointInsideTriangle( ptIntTriaCentroid, trTriaA, TriangleType::EXACT)  &&
                                   IsPointInsideTriangle( ptIntTriaCentroid, trTriaB, TriangleType::EXACT) ;
            }
            // Se l'intersezeione è corretta e i triangoli sono sovrapposti aggiorno gli indici.
            if ( nIntTypeAB == ITTT_OVERLAPS  &&  bSuccesfullInters) {
               m_vTria[nTA].nTempPart = 2 ;
               SurfB.m_vTria[vNearTria[nTB]].nTempPart = 2 ;
            }
            else if ( nIntTypeAB == ITTT_COUNTER_OVERLAPS  &&  bSuccesfullInters) {
               m_vTria[nTA].nTempPart = - 2 ;
               SurfB.m_vTria[vNearTria[nTB]].nTempPart = - 2 ;  
            }
         }
      }
      return ( AdjustVertices()  &&  DoCompacting()  &&  SurfB.AdjustVertices()  &&  SurfB.DoCompacting()) ;
   }
#endif
   return bOk ;
}