//----------------------------------------------------------------------------
//  EgalTech 2023
//----------------------------------------------------------------------------
// File : Tree.cpp        Data : 21.04.23          Versione : 
// Contenuto : Implementazione della classe Tree.
//
//
//
// Modifiche : 21.04.23 DB Creazione modulo.
//
//
//----------------------------------------------------------------------------

//--------------------------- Include ----------------------------------------
#include "stdafx.h"
#include "Tree.h"
#include "SurfBezier.h"
#include "GeoConst.h"
#include "CurveLine.h"
#include "CurveComposite.h"
#include "SurfFlatRegion.h"
#include "IntersLineLine.h"
#include "AdjustLoops.h"
#include "/EgtDev/Include/EGkPolyLine.h"
#include "/EgtDev/Include/EGkCurve.h"
#include "/EgtDev/Include/EGkSfrCreate.h"
#include "/EgtDev/Include/EGkDistPointCurve.h"
#include "/EgtDev/Include/EGkDistLineLine.h"
#include "/EgtDev/Include/EGkDistPointTria.h"
#include "/EgtDev/Include/EGkStringUtils3d.h"
#include <algorithm>
#include <numeric>

using namespace std ;

//----------------------------------------------------------------------------
Tree::Tree( void)
   : m_pSrfBz( nullptr), m_bTrimmed( false), m_bBilinear( false), m_bMulti( false), m_bClosedU( false), m_bClosedV( false),
     m_vbPole{ false, false, false, false}
{
   Point3d ptBl( 0, 0), ptTr ( 1 * SBZ_TREG_COEFF, 1 * SBZ_TREG_COEFF) ;
   Cell cRoot( ptBl, ptTr) ;
   m_mTree.insert( pair< int, Cell>( -1, cRoot)) ;
}

//----------------------------------------------------------------------------
Tree::~Tree( void)
{
}

//----------------------------------------------------------------------------
Tree::Tree( const Point3d ptBl, const Point3d ptTr)
   : m_pSrfBz( nullptr), m_bTrimmed( false), m_bBilinear( false), m_bMulti( false), m_bClosedU( false), m_bClosedV( false),
     m_vbPole{ false, false, false, false}
{
   Cell cRoot( ptBl, ptTr) ;
   m_mTree.insert( pair< int, Cell>( -1, cRoot)) ;
}

//----------------------------------------------------------------------------
bool
Tree::VerifyLoopOrientation( ICURVEPLIST& vpCrv, BOOLVECTOR& vbOrientation) const
{
   // verifico che il verso dei loop sia corretto controllando le relazioni tra loop
   vpCrv.sort( []( ICurve* a, ICurve* b) { double AreaA, AreaB ; a->GetAreaXY( AreaA) ; b->GetAreaXY( AreaB) ; return abs(AreaA) > abs(AreaB) ;}) ;
   auto it = vpCrv.begin() ;
   for ( ++it ; it != vpCrv.end() ; ++it) {
      bool bIdentified = false ;
      for ( auto k = it ; k != vpCrv.begin() ;) {
         --k ;
         IntersCurveCurve icc( **k, **it) ;
         int nRes = icc.GetRegionCurveClassification() ;
         if ( nRes == CCREGC_IN1) {
            bIdentified = true ;
            vbOrientation.push_back( ! vbOrientation[ distance( vpCrv.begin(), k)]) ; // l'orientazione deve essere opposta alla prima curva che contiene la corrente
            break ;
         }
      }
      if ( ! bIdentified)
         vbOrientation.push_back( vbOrientation[0]) ;
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
Tree::AdjustLoop( PolyLine& pl, POLYLINEVECTOR& vPl, BOOLVECTOR& vbOrientation) const
{
   PtrOwner<ICurveComposite> pCC( CreateCurveComposite()) ;
   pCC->FromPolyLine( pl) ;
   ICURVEPLIST vCrv ;
   AdjustLoops( Release( pCC), vCrv, false) ;

   if ( vCrv.size() > 1)
      VerifyLoopOrientation( vCrv, vbOrientation) ;
   for ( auto itCrv = vCrv.begin() ; itCrv != vCrv.end() ; ++itCrv) {
      PolyLine plApprox ; (*itCrv)->ApproxWithLines( 0.01, 15, 0, plApprox) ;
      vPl.push_back( plApprox) ;
      delete (*itCrv) ;
   }
   return true ;
}

//----------------------------------------------------------------------------
static double
GetHalfKey( double dPar)
{
   const int TWO_TO_15 = 32768 ;
   return ( dPar * TWO_TO_15) ;
}

//----------------------------------------------------------------------------
bool
Tree::GetPoint( double dU, double dV, Point3d& ptP) const
{
   pair<int64_t, int64_t> key( int64_t( GetHalfKey( dU)), int64_t( GetHalfKey( dV))) ;
   bool bOk = true ;
   if ( m_mPt3d.find( key) == m_mPt3d.end()) {
      bOk = bOk &&  m_pSrfBz->GetPoint( dU / SBZ_TREG_COEFF, dV / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptP) ;
      m_mPt3d[key] = ptP ;
   }
   ptP = m_mPt3d[key] ;
   return bOk ;
}

//----------------------------------------------------------------------------
bool
Tree::SavePoint( double dU, double dV, const Point3d& ptP)
{
   pair<int64_t, int64_t> key( int64_t( GetHalfKey( dU)), int64_t( GetHalfKey( dV))) ;
   if ( m_mPt3d.find( key) == m_mPt3d.end())
      m_mPt3d[key] = ptP ;
   return true ;
}

//----------------------------------------------------------------------------
bool
Tree::SetSurf( const SurfBezier* pSrfBz, const Point3d& ptMin, const Point3d& ptMax)
{
   if ( pSrfBz == nullptr  ||  ! pSrfBz->IsValid())
      return false ;
  // pulisco i vettori membri
   m_mTree.clear() ;
   m_vnLeaves.clear() ;
   m_vnParents.clear() ;
   m_mPt3d.clear() ;
   m_mChunk.clear() ;
   m_vPlApprox.clear() ;
   m_vCCLoop2D.clear() ;

   m_pSrfBz = pSrfBz ;
  // le coordinate delle celle sono nello spazio parametrico
   int nDegU, nDegV, nSpanU, nSpanV ;
   bool bIsRat, bTrimmed ;
   m_pSrfBz->GetInfo( nDegU, nDegV, nSpanU, nSpanV, bIsRat, bTrimmed) ;
   m_bTrimmed = bTrimmed ;
   m_nDegU = nDegU ;
   m_nDegV = nDegV ;
   m_nSpanU = nSpanU ;
   m_nSpanV = nSpanV ;
   if ( nDegU == 1  &&  nDegV == 1) 
      m_bBilinear = true ;
   if ( nSpanU * nSpanV != 1)
      m_bMulti = true ;
  // creo la cella Root
   Point3d ptTop( nSpanU * SBZ_TREG_COEFF, nSpanV * SBZ_TREG_COEFF) ;
   bool bLimited = false ;
   if ( ! AreSamePointXYExact( ptMax,ORIG)) {
      if ( ! AreSamePointXYExact( ptMax,ptTop))
         ptTop = ptMax ;
      bLimited = true ;
   }
   Cell cRoot( ptMin, ptTop) ;
   m_mTree.insert( pair< int, Cell>( -1, cRoot)) ;
  // recupero i loop di trim e li divido per chunk
   if ( m_bTrimmed) {
      int nLoop = 0 ;
     // recupero la superficie di trim per avere accesso diretto ai loop e mantenendo le informazioni sui chunk
      PtrOwner<SurfFlatRegion> pTrimReg( m_pSrfBz->GetTrimRegion()->Clone()) ;
      double dLinTol = 0.005 ; // questo è riferito allo spazio parametrico
      double dAngTolDeg = 5 ;
      for ( int i = 0 ; i < pTrimReg->GetChunkCount() ; ++ i) {
         PtrOwner<SurfFlatRegion> pChunk( pTrimReg->CloneChunk( i)) ;
         for ( int j = 0 ; j < pChunk->GetLoopCount( 0) ; ++ j) {
           // i chunk della falt region sono ancora flat region composte da 1 chunk
            // rimuovo i difetti dei loop prima di salvarli
            PtrOwner<ICurveComposite> pLoop( GetBasicCurveComposite( pChunk->GetLoop( 0, j))) ;
            if ( IsNull( pLoop))
               return false ;
            pLoop->MergeCurves( dLinTol, dAngTolDeg) ;
            pLoop->RemoveSmallDefects( dLinTol, dAngTolDeg, true) ;
            pLoop->RemoveSmallParts( dLinTol, dAngTolDeg) ;
           // approssimo i loop di trim con delle spezzate
            PolyLine plApprox ;
            int nType = 0 ;
            pLoop->ApproxWithLines( dLinTol,dAngTolDeg, nType, plApprox) ;
           // calcolo se il loop è CCW o CW
            double dArea ;
            Plane3d plExtPlane ;
            bool bCCW ;
            plApprox.IsClosedAndFlat( plExtPlane, dArea, 50 * EPS_SMALL) ;
            if ( plExtPlane.GetVersN().z > 0)
               bCCW = true ;
            else
               bCCW = false ;

            POLYLINEVECTOR vPlAdjusted ;
            BOOLVECTOR vbOrientation ;
            vbOrientation.push_back( bCCW) ;
            AdjustLoop( plApprox, vPlAdjusted, vbOrientation) ;

            nLoop = int ( m_vPlApprox.size()) ;
            for ( int k = 0 ; k < int( vPlAdjusted.size()) ; ++k )
               m_vPlApprox.push_back( tuple<PolyLine,bool>(vPlAdjusted[k], vbOrientation[k])) ;
           // aggiorno la mappa del chunk di appartenenza per tutti i loop aggiunti  // do per scontato che siano tutti dello stesso chunk anche se li ho separati
            for ( int k = nLoop ; k < int( m_vPlApprox.size()); ++k)
               m_mChunk[k] = i ;
         }
      }
   }
  // salvo i vertici 3d della cella root
   Point3d ptP00, ptP10, ptP11, ptP01 ;
   bool bOk = false ;
   if ( ! bLimited) {
      ptP00 = m_pSrfBz->GetControlPoint( 0, &bOk) ;
      SavePoint( 0, 0, ptP00) ;
      ptP10 = m_pSrfBz->GetControlPoint( nDegU * nSpanU, &bOk) ;
      SavePoint( nSpanU * SBZ_TREG_COEFF, 0, ptP10) ;
      ptP11 = m_pSrfBz->GetControlPoint( ( nDegU * nSpanU + 1) * ( nDegV * nSpanV + 1) - 1, &bOk) ;
      SavePoint( nSpanU * SBZ_TREG_COEFF, nSpanV * SBZ_TREG_COEFF, ptP11) ;
      ptP01 = m_pSrfBz->GetControlPoint( ( nDegU * nSpanU + 1) * ( nDegV * nSpanV), &bOk) ;
      SavePoint( 0, nSpanV * SBZ_TREG_COEFF, ptP01) ;
   }
   else {
      GetPoint( ptMin.x, ptMin.y, ptP00) ;
      GetPoint( ptMax.x, ptMin.y, ptP10) ;
      GetPoint( ptMax.x, ptMax.y, ptP11) ;
      GetPoint( ptMin.x, ptMax.y, ptP01) ;
   }
  // se richiesto divido preliminarmente le patches
   m_vnParents.clear() ;
   bool bIsPlanar = m_pSrfBz->IsPlanar() ;
   if ( ! bIsPlanar  ||  m_bMulti) {
     // se la superficie è chiusa lungo il parametro U, sistemo le adiacenze al bordo
      if ( AreSamePointApprox( ptP00, ptP10) && AreSamePointApprox( ptP01, ptP11) ) {
         Point3d ptP0001 ;
         GetPoint( ptMin.x, ( ptMin.y + ptTop.y) / 2, ptP0001) ;
         Point3d ptP1011 ;
         GetPoint( ptTop.x, ( ptMin.y + ptTop.y) / 2, ptP1011) ;
         if ( AreSamePointApprox( ptP0001, ptP1011)) {
            m_mTree[-1].m_nLeft = -1 ;
            m_mTree[-1].m_nRight = -1 ;
            m_bClosedU = true ;
         }
      }
     // se la superficie è chiusa lungo il parametro V, sistemo le adiacenze al bordo
      if ( ( AreSamePointApprox( ptP00, ptP01) && AreSamePointApprox( ptP10, ptP11) ) ) {
         Point3d ptP0010 ;
         GetPoint( ( ptMin.x + ptTop.x) / 2, ptMin.y, ptP0010) ;
         Point3d ptP0111 ;
         GetPoint( (ptMin.x + ptTop.x) / 2, ptTop.y, ptP0111) ;
         if ( AreSamePointApprox( ptP0010, ptP0111)) {
            m_mTree[-1].m_nTop = -1 ;
            m_mTree[-1].m_nBottom = -1 ;
            m_bClosedV = true ;
         }
      }
      if ( nSpanU > 1  ||  nSpanV > 1  ||  m_bClosedU  ||  m_bClosedV) {
         int nId = -1 ;
         int nSplit = nSpanU ;
         bool bSplitInHalf = false ;
        // se la superficie è chiusa lungo un parametro la splitto lungo quel parametro almeno una volta
         if ( m_bClosedU  &&  nSpanU == 1) {
            nSplit = 2 ;
            bSplitInHalf = true ;
         }
         for ( int i = 1 ; i < nSplit ; ++i) {
            m_mTree[nId].SetSplitDirVert( true) ;
            double dSplit = ( bSplitInHalf ? i * SBZ_TREG_COEFF / 2 : i * SBZ_TREG_COEFF) ;
            if ( Split( nId, dSplit))
               nId = m_mTree[nId].m_nChild2 ;
         }
         INTVECTOR vLeaves ;
         GetHeightLeaves( -1, vLeaves) ;
         nSplit = nSpanV - 1 ;
         bSplitInHalf = false ;
        // se la superficie è chiusa lungo un parametro la splitto lungo quel parametro almeno una volta
         if ( m_bClosedV &&  nSpanV == 1) {
            nSplit = 1 ;
            bSplitInHalf = true ;
         }
         for ( int nId : vLeaves) {
            for ( int j = nSplit ; j > 0 ; --j) {
               m_mTree[nId].SetSplitDirVert( false) ;
               double dSplit = bSplitInHalf ? j * SBZ_TREG_COEFF / 2 : j * SBZ_TREG_COEFF ;
               if ( Split( nId, dSplit))
                  nId = m_mTree[nId].m_nChild2 ;
            }
         }
         vLeaves.clear() ;
      }
     // controllo se la superficie è chiusa.
     // se è chiusa e non ho già fatto split preliminare, splitto sul parametro su cui è chiusa
     // e sistemo le adiacenze
      if ( ( AreSamePointApprox( ptP00, ptP01)  ||  AreSamePointApprox( ptP10, ptP11)) ||
           ( AreSamePointApprox( ptP00, ptP10)  ||  AreSamePointApprox( ptP01, ptP11))) {
         if ( ( AreSamePointApprox( ptP00, ptP01)  ||  AreSamePointApprox( ptP10, ptP11))) {

           // qui devo fare il controllo capped ( chiusura a semisfera)
           // devo controllare se i punti ai parametri U=0 e U=1 sono tutti coincidenti
           // in caso devo fare uno split nell'altra direzione
            bool bOk = false ;
            bool bPole0 = true, bPole1 = true ;
            Point3d ptU0, ptU1 ;
           // controllo se tutti i punti di controllo sull'isoparametrica sono uguali
            for ( int i = 1 ; i < nDegV * nSpanV + 1 ; ++ i) {
               ptU0 = m_pSrfBz->GetControlPoint( i * ( nDegU * nSpanU + 1), &bOk) ;
               bPole0 = bPole0  &&  AreSamePointApprox( ptP00, ptU0) ;
               ptU1 = m_pSrfBz->GetControlPoint( ( i + 1) * ( nDegU * nSpanU + 1) - 1, &bOk) ;
               bPole1 = bPole1  &&  AreSamePointApprox( ptP10, ptU1) ;
            }
            m_vbPole[1] = bPole0 ;
            m_vbPole[3] = bPole1 ;
            if ( bPole0  &&  bPole1  &&  m_mTree.size() == 3) {
               m_mTree[0].SetSplitDirVert( true) ;
               Split( 0) ;
               m_mTree[1].SetSplitDirVert( true) ;
               Split( 1) ;
            }
         }
        // nella condizione di questo if non controllo eventuali divisioni preliminari, perché ne tengo conto dopo
         if ( AreSamePointApprox( ptP00, ptP10)  ||  AreSamePointApprox( ptP01, ptP11)) {

           // devo controllare se i punti ai parametri V=0 e V=1 sono tutti coincidenti
           // in caso devo fare uno split nell'altra direzione
            bool bOk = false ;
            bool bPole0 = true, bPole1 = true ;
           // controllo se tutti i punti sull'isoparametrica sono uguali
            for ( int i = 1 ; i < nDegU * nSpanU + 1 ; ++ i) {
               Point3d ptV0 = m_pSrfBz->GetControlPoint( i, &bOk) ;
               bPole0 = bPole0  &&  AreSamePointApprox( ptP00, ptV0) ;
               Point3d ptV1 = m_pSrfBz->GetControlPoint( i + ( nDegU * nSpanU + 1) * ( nDegV * nSpanV), &bOk) ;
               bPole1 = bPole1  &&  AreSamePointApprox( ptP01, ptV1) ;
            }
            m_vbPole[0] = bPole1 ;
            m_vbPole[2] = bPole0 ;
            if ( bPole0  &&  bPole1  &&  int( m_mTree.size()) == 3) {
               m_mTree[0].SetSplitDirVert( false) ;
               Split( 0) ;
               m_mTree[1].SetSplitDirVert( false) ;
               Split( 1) ;
            }
           // se ho fatto solo 1 split orizzontale e ho due celle foglie nId = 0 e nId = 1
            if ( m_mTree.size() == 3  &&  ! m_mTree.at(-1).IsSplitVert()) {
               m_mTree[0].SetSplitDirVert( true) ;
               Split( 0) ;
               m_mTree[1].SetSplitDirVert( true) ;
               Split( 1) ;
            }
         }
      }
   }
  // calcolo i parent che ho creato con le eventuali divisioni preliminari
   INTVECTOR vLeaves ;
   GetHeightLeaves( -1, vLeaves) ;
   m_vnParents = vLeaves ;

   return true ;
}

//----------------------------------------------------------------------------
static bool
AddOrMergeBBox( const BBox3d& bBox3dA, vector<BBox3d>& vBBox, bool bAdd = true, int nInd = 0)
{
   Point3d ptMin = bBox3dA.GetMin() ;
   Point3d ptMax = bBox3dA.GetMax() ;
   if ( vBBox.empty()) {
      vBBox.push_back( bBox3dA) ;
      return true ;
   }
   bool bAdded = false ;
   for ( int b = 0 ; b < int( vBBox.size()) ; ++b) {
      BBox3d bBox3dB = vBBox[b] ;
      BBox3d b3Int ;
     // se sono celle diverse e ho un'intersezione faccio il merge
      if ( ! ( AreSamePointXYExact( ptMin, bBox3dB.GetMin())  &&  AreSamePointXYExact( ptMax, bBox3dB.GetMax()))  &&
           bBox3dA.FindIntersectionXY( bBox3dB, b3Int)) {
         vBBox[b].Add( bBox3dA) ;
         if ( ! bAdd ) {
            vBBox.erase( vBBox.begin() + nInd) ;
            -- b ;
         }
        // se ho fatto un merge devo controllare se ora la nuova bbox ha delle intersezioni
         AddOrMergeBBox( vBBox[b], vBBox, false, b) ;
         bAdded = true ;
         break ;
      }
   }
   if ( ! bAdded  &&  bAdd)
      vBBox.push_back( bBox3dA) ;
   return true ;
}

//----------------------------------------------------------------------------
bool
Tree::GetIndependentTrees( BIPNTVECTOR& vTrees)
{
   if ( ! m_bTrimmed) {
      vTrees.emplace_back( ORIG, ORIG) ;
   }
   else {
     // se ho dei loop di trim trovo le loro BBox3d per costruire l'albero solo all'interno di queste BBox
      BOXVECTOR vBBox ;
      for ( int i = 0 ; i < int( m_vPlApprox.size()) ; ++ i) {
         PolyLine& plLoop = get<0>( m_vPlApprox[i]) ;
        // calcolo la BBox3d
         Point3d ptP ;
         plLoop.GetFirstPoint( ptP) ;
         BBox3d bBox3dA( ptP) ;
         while ( plLoop.GetNextPoint( ptP))
            bBox3dA.Add( ptP) ;
        // controllo se ho intersezioni con altre bbox
        // se ho intersezioni unisco le bbox, altrimenti le lascio indipendenti
         AddOrMergeBBox( bBox3dA, vBBox) ;
      }
     // controllo se dopo aver unito le bbox ho ottenuto Root
      bool bIsRoot = false ;
      Point3d ptTR( m_nSpanU * SBZ_TREG_COEFF, m_nSpanV * SBZ_TREG_COEFF) ;
      for ( int i = 0 ; i < int( vBBox.size()) ; ++ i) {
         BBox3d& bBox3d = vBBox[i] ;
         if ( AreSamePointXYEpsilon( bBox3d.GetMin(), ORIG, 10)  &&  AreSamePointXYEpsilon( bBox3d.GetMax(), ptTR, 10)) {
            bIsRoot = true ;
            break ;
         }
      }
     // restituisco le celle parent di partenza a partire dalle bbox che ho ottenuto
      if ( ! bIsRoot) {
         for ( int i = 0 ; i < int( vBBox.size()) ; ++ i) {
            Point3d ptMin = vBBox[i].GetMin() ;
            Point3d ptMax = vBBox[i].GetMax() ;
            vTrees.emplace_back( ptMin, ptMax) ;
         }
      }
      else
         vTrees.emplace_back( ORIG, ORIG);
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
Tree::Split( int nId, double dSplitValue)
{
   Cell& cToSplit = m_mTree.at(nId) ;
  // controllo che lo split non venga fatto sul lato della cella
   bool bGoodSplitVert = cToSplit.IsSplitVert()  &&  dSplitValue > cToSplit.GetBottomLeft().x + 10 * EPS_SMALL  &&
                         dSplitValue < cToSplit.GetTopRight().x - 10 * EPS_SMALL ;
   bool bGoodSplitHoriz = ! cToSplit.IsSplitVert()  &&  dSplitValue > cToSplit.GetBottomLeft().y + 10 * EPS_SMALL  &&
                          dSplitValue < cToSplit.GetTopRight().y - 10 * EPS_SMALL ;
   Point3d ptP00, ptP01, ptP10, ptP11 ;
   
   if ( bGoodSplitVert) {
      if ( cToSplit.GetBottomRight().x - dSplitValue > dSplitValue - cToSplit.GetBottomLeft().x) {
         GetPoint( cToSplit.GetBottomLeft().x, cToSplit.GetBottomLeft().y, ptP00) ;
         GetPoint( dSplitValue, cToSplit.GetBottomRight().y, ptP10) ;
         GetPoint( cToSplit.GetTopLeft().x, cToSplit.GetTopLeft().y, ptP01) ;
         GetPoint( dSplitValue, cToSplit.GetTopRight().y, ptP11) ;
      }
      else {
         GetPoint( dSplitValue, cToSplit.GetBottomLeft().y, ptP00) ;
         GetPoint( cToSplit.GetBottomRight().x, cToSplit.GetBottomRight().y, ptP10) ;
         GetPoint( dSplitValue, cToSplit.GetTopLeft().y, ptP01) ;
         GetPoint( cToSplit.GetTopRight().x, cToSplit.GetTopRight().y, ptP11) ;
      }
      if ( AreSamePointApprox( ptP00, ptP10)  &&  AreSamePointApprox( ptP01, ptP11)  &&
          ( cToSplit.GetBottomRight().x - dSplitValue < SBZ_TREG_COEFF - EPS_SMALL  ||
            dSplitValue - cToSplit.GetBottomLeft().x < SBZ_TREG_COEFF - EPS_SMALL))
         bGoodSplitVert = false ;
   }
   if ( bGoodSplitHoriz) {
      if ( cToSplit.GetTopLeft().y - dSplitValue > dSplitValue - cToSplit.GetBottomLeft().y) {
         GetPoint( cToSplit.GetBottomLeft().x, cToSplit.GetBottomLeft().y, ptP00) ;
         GetPoint( cToSplit.GetBottomRight().x, cToSplit.GetBottomRight().y, ptP10) ;
         GetPoint( cToSplit.GetTopLeft().x, dSplitValue, ptP01) ;
         GetPoint( cToSplit.GetTopRight().x, dSplitValue, ptP11) ;
      }
      else {
         GetPoint( cToSplit.GetBottomLeft().x, dSplitValue, ptP00) ;
         GetPoint( cToSplit.GetBottomRight().x, dSplitValue, ptP10) ;
         GetPoint( cToSplit.GetTopLeft().x, cToSplit.GetTopLeft().y, ptP01) ;
         GetPoint( cToSplit.GetTopRight().x, cToSplit.GetTopRight().y, ptP11) ;
      }
      if ( AreSamePointApprox( ptP00, ptP01)  &&  AreSamePointApprox( ptP10, ptP11)  &&
          ( cToSplit.GetTopLeft().y - dSplitValue < SBZ_TREG_COEFF - EPS_SMALL  ||
            dSplitValue - cToSplit.GetBottomLeft().y < SBZ_TREG_COEFF - EPS_SMALL))
         bGoodSplitHoriz = false ;
   }
   
   if ( bGoodSplitVert  ||  bGoodSplitHoriz) {
      cToSplit.m_dSplit = dSplitValue ;
      Cell cChild1, cChild2 ;
      cChild1.m_nDepth = cToSplit.m_nDepth + 1 ;
      cChild2.m_nDepth = cToSplit.m_nDepth + 1 ;
      int nNodes = (int) m_mTree.size() ;
      cChild1.m_nId = nNodes - 1 ;
      cToSplit.m_nChild1 = nNodes - 1 ;
      cChild2.m_nId = nNodes ;
      cToSplit.m_nChild2 = nNodes ;
      Point3d ptVert1, ptVert2 ;
      if ( ! cToSplit.IsSplitVert()) {
         // la cella figlio 1 è quella sopra
         Point3d ptBL( cToSplit.GetBottomLeft().x, dSplitValue) ;
         cChild1.SetBottomLeft( ptBL) ;
         cChild1.SetTopRight( cToSplit.GetTopRight()) ;
         cChild1.m_nTop = cToSplit.m_nTop ;
         cChild1.m_nBottom = cToSplit.m_nChild2 ;
         cChild1.m_nLeft = cToSplit.m_nLeft ;
         cChild1.m_nRight = cToSplit.m_nRight ;
         Point3d ptTR( cToSplit.GetTopRight().x, dSplitValue) ;
         cChild2.SetBottomLeft( cToSplit.GetBottomLeft()) ;
         cChild2.SetTopRight( ptTR) ;
         cChild2.m_nTop = cToSplit.m_nChild1 ;
         cChild2.m_nBottom = cToSplit.m_nBottom ;
         cChild2.m_nLeft = cToSplit.m_nLeft ;
         cChild2.m_nRight = cToSplit.m_nRight ;
         // metto i corrispondenti 3d dei punti dello split nella mappa m_mVert
         // per ogni cella i punti devono essere nell'ordine ptP00, ptP10, ptP11, ptP01
         GetPoint( cToSplit.GetBottomLeft().x, dSplitValue, ptVert1) ;
         GetPoint( cToSplit.GetTopRight().x, dSplitValue, ptVert2) ;
      }
      else {
         // la cella figlio 1 è quella di sinistra
         Point3d ptTR( dSplitValue, cToSplit.GetTopRight().y) ;
         cChild1.SetBottomLeft( cToSplit.GetBottomLeft()) ;
         cChild1.SetTopRight( ptTR) ;
         cChild1.m_nTop = cToSplit.m_nTop ;
         cChild1.m_nBottom = cToSplit.m_nBottom ;
         cChild1.m_nLeft = cToSplit.m_nLeft ;
         cChild1.m_nRight = cToSplit.m_nChild2 ;
         Point3d ptBL( dSplitValue, cToSplit.GetBottomLeft().y) ;
         cChild2.SetBottomLeft( ptBL) ;
         cChild2.SetTopRight( cToSplit.GetTopRight()) ;
         cChild2.m_nTop = cToSplit.m_nTop ;
         cChild2.m_nBottom = cToSplit.m_nBottom ;
         cChild2.m_nLeft = cToSplit.m_nChild1 ;
         cChild2.m_nRight = cToSplit.m_nRight ;
         // metto i corrispondenti 3d dei punti dello split nella mappa m_mVert
         // per ogni cella i punti devono essere nell'ordine ptP00, ptP10, ptP11, ptP01
         GetPoint( dSplitValue, cToSplit.GetBottomLeft().y, ptVert2) ;
         GetPoint( dSplitValue, cToSplit.GetTopRight().y, ptVert1) ;
      }
      cChild1.SetParent( nId) ;
      cChild2.SetParent( nId) ;
      // inserisco nell'albero
      m_mTree.insert( pair<int, Cell>( nNodes - 1, cChild1)) ;
      m_mTree.insert( pair<int, Cell>( nNodes, cChild2)) ;
      return true ;
   }
   return false ;
}

//----------------------------------------------------------------------------
bool
Tree::Split( int nId)
{
   double dValue ;
   Cell& cCell =  m_mTree.at( nId) ;
   if ( cCell.IsSplitVert())
      dValue = ( cCell.GetBottomLeft().x + cCell.GetTopRight().x) / 2 ;
   else
      dValue = ( cCell.GetBottomLeft().y + cCell.GetTopRight().y) / 2 ;
   return Split( nId, dValue) ;
}

//----------------------------------------------------------------------------
bool
Tree::BuildTree( double dLinTol, double dSideMin, double dSideMax)
{
  // suddivido lo spazio parametrico con divisioni a metà su uno dei due parametri
   int nCToSplit = -1 ;
   Cell* pcToSplit = &m_mTree[nCToSplit] ;
   bool bIsPlanar = m_pSrfBz->IsPlanar() ;
   if ( ! m_bBilinear) {
      while ( nCToSplit != -2  &&  ! pcToSplit->IsProcessed()) {

        // se la pezza non è già stata preliminarmente splittata
         if ( pcToSplit->IsLeaf()) {

           // dimensioni parametriche della pezza
            double dLenParU = ( pcToSplit->GetTopRight().x - pcToSplit->GetBottomLeft().x) / SBZ_TREG_COEFF ;
            double dLenParV = ( pcToSplit->GetTopRight().y - pcToSplit->GetBottomLeft().y) / SBZ_TREG_COEFF ;

           // vertici e centro della pezza
            Point3d ptP00, ptP10, ptP11, ptP01, ptCen ;
            GetPoint( pcToSplit->GetBottomLeft().x, pcToSplit->GetBottomLeft().y, ptP00) ;
            GetPoint( pcToSplit->GetTopRight().x, pcToSplit->GetBottomLeft().y, ptP10) ;
            GetPoint( pcToSplit->GetTopRight().x, pcToSplit->GetTopRight().y, ptP11) ;
            GetPoint( pcToSplit->GetBottomLeft().x, pcToSplit->GetTopRight().y, ptP01) ;
            GetPoint( pcToSplit->GetCenter().x, pcToSplit->GetCenter().y, ptCen) ;
            //{ string sLog = "P00=" + ToString( ptP00, 3) + " P10=" + ToString( ptP10, 3) +
            //                " P11=" + ToString( ptP11, 3) + " P01=" + ToString( ptP01, 3) +
            //                " Cen=" + ToString( ptCen, 3);
            //  LOG_DBG_INFO( GetEGkLogger(), sLog.c_str())}

           // distanze sul contorno tra i quattro vertici
            double dLen0 = Dist( ptP00, ptP10) ;
            double dLen1 = Dist( ptP10, ptP11) ;
            double dLen2 = Dist( ptP01, ptP11) ;
            double dLen3 = Dist( ptP00, ptP01) ;
           // distanza reale tra i vertici della pezza
            if ( dLen0 < EPS_ZERO  &&  dLen2 < EPS_ZERO ) {
               double dV = pcToSplit->GetCenter().y / SBZ_TREG_COEFF ;
               PtrOwner<CurveComposite> pCrvV( m_pSrfBz->GetCurveOnU( dV)) ;
               double dLenU0, dLenU1 ;
               pCrvV->GetLengthAtParam( pcToSplit->GetBottomLeft().x / SBZ_TREG_COEFF, dLenU0) ;
               pCrvV->GetLengthAtParam( pcToSplit->GetTopRight().x / SBZ_TREG_COEFF, dLenU1) ;
               dLen0 = abs( dLenU1 - dLenU0) ;
            }
            if ( dLen1 < EPS_ZERO  &&  dLen3 < EPS_ZERO ) {
               double dU = pcToSplit->GetCenter().x / SBZ_TREG_COEFF ;
               PtrOwner<CurveComposite> pCrvU( m_pSrfBz->GetCurveOnV( dU)) ;
               double dLenV0, dLenV1 ;
               pCrvU->GetLengthAtParam( pcToSplit->GetBottomLeft().y / SBZ_TREG_COEFF, dLenV0) ;
               pCrvU->GetLengthAtParam( pcToSplit->GetTopRight().y / SBZ_TREG_COEFF, dLenV1) ;
               dLen1 = abs( dLenV1 - dLenV0) ;
            }

           // calcolo massimo della freccia (sagitta)
            double dSagU = 0, dSagV = 0 ;

           // da distanza tra centro e triangoli di approssimazione
            if ( dSagU < dLinTol  &&  dSagV < dLinTol) {
              // distanza del centro dai due triangoli con la prima diagonale
               Triangle3d Tria1A, Tria1B ;
               Tria1A.Set( ptP00, ptP10, ptP11) ; Tria1A.Validate( true) ;
               Tria1B.Set( ptP00, ptP11, ptP01) ; Tria1B.Validate( true) ;
               double dDist1A = NAN ;
               DistPointTriangle( ptCen, Tria1A).GetDist( dDist1A) ;
               double dDist1B = NAN ;
               DistPointTriangle( ptCen, Tria1B).GetDist( dDist1B) ;
               double dDist1 = min( dDist1A, dDist1B) ;
              // distanza del centro dai due triangoli con la seconda diagonale
               Triangle3d Tria2A, Tria2B ;
               Tria2A.Set( ptP10, ptP11, ptP01) ; Tria2A.Validate( true) ;
               Tria2B.Set( ptP10, ptP01, ptP00) ; Tria2B.Validate( true) ;
               double dDist2A = NAN ;
               DistPointTriangle( ptCen, Tria2A).GetDist( dDist2A) ;
               double dDist2B = NAN ;
               DistPointTriangle( ptCen, Tria2B).GetDist( dDist2B) ;
               double dDist2 = min( dDist2A, dDist2B) ;
              // prendo la minore delle due distanze
               double dDist = min( dDist1, dDist2) ;
              // se maggiore della tolleranza, imposto come freccia
               if ( isfinite( dDist)  &&  dDist > dLinTol) {
                 // divido in base alle distanze tra i vertici dei lati
                  if ( max( dLen0, dLen2) >= max( dLen1, dLen3))
                     dSagU = dDist ;
                  else
                     dSagV = dDist ;
                  //{ string sLog = " Da triangoli : FrecciaU=" + ToString( dSagU, 3) + " FrecciaV=" + ToString( dSagV, 3) ;
                  //  LOG_DBG_INFO( GetEGkLogger(), sLog.c_str())}
               }
            }

            // calcolo se la parte di superficie nella cella è piatta
            PolyLine PL ;
            PL.AddUPoint( 0, ptP00) ;
            PL.AddUPoint( 1, ptP10) ;
            PL.AddUPoint( 2, ptP11) ;
            PL.AddUPoint( 3, ptP01) ;
            PL.AddUPoint( 4, ptCen) ;
            Plane3d plPlane ;
            bool bIsFlat = PL.IsFlat( plPlane, dLinTol) ;

           // su isoparametriche in U e V
            if ( dSagU < dLinTol  &&  dSagV < dLinTol  &&  ! bIsFlat) {
              // step di verifica in U e in V
               int nStepU = ( dLenParU > 1. / m_nDegU ? m_nDegU + 1 : 2) ;
               int nStepV = ( dLenParV > 1. / m_nDegV ? m_nDegV + 1 : 2) ;
              // verifico lungo curve in U
               for ( int i = 0 ; i <= nStepU ; ++ i) {
                 // parametro U in analisi
                  double dCoeffU = double( i) / nStepU ;
                  double dU = ( 1 - dCoeffU) * pcToSplit->GetBottomLeft().x + dCoeffU * pcToSplit->GetTopRight().x ;
                 // linea tra gli estremi a questa U (ptP00P10 è un punto tra P00 e P10 e così via)
                  Point3d ptP00P10 ;
                  GetPoint( dU, pcToSplit->GetBottomLeft().y, ptP00P10) ;
                  Point3d ptP11P01 ;
                  GetPoint( dU, pcToSplit->GetTopRight().y, ptP11P01) ;
                  CurveLine clV ;
                  clV.Set( ptP00P10, ptP11P01) ;
                 // verifico alcuni punti intermedi in V
                  for ( int j = 1 ; j < nStepV ; ++ j) {
                    // parametro in V
                     double dCoeffV = double( j) / nStepV ;
                     double dV = ( 1 - dCoeffV) * pcToSplit->GetBottomLeft().y + dCoeffV * pcToSplit->GetTopRight().y ;
                    // punto a U, V
                     Point3d ptPSrf ;
                     GetPoint( dU, dV, ptPSrf) ;
                    // distanza del punto dalla precedente linea
                     DistPointCurve dpc( ptPSrf, clV) ;
                     double dDist ;
                     dpc.GetDist( dDist) ;
                     dSagV = max( dSagV, dDist) ;
                  }
               }
              // verifico lungo curve in V
               for ( int i = 0 ; i <= nStepV ; ++ i) {
                 // parametro in V in analisi
                  double dCoeffV = double( i) / nStepV ;
                  double dV = ( 1 - dCoeffV) * pcToSplit->GetBottomLeft().y + dCoeffV * pcToSplit->GetTopRight().y ;
                 // linea tra gli estremi a questa V
                  Point3d ptP10P11 ;
                  GetPoint( pcToSplit->GetTopRight().x, dV, ptP10P11) ;
                  Point3d ptP01P00 ;
                  GetPoint( pcToSplit->GetBottomLeft().x, dV, ptP01P00) ;
                  CurveLine clU ;
                  clU.Set( ptP01P00, ptP10P11) ;
                 // verifico alcuni punti intermedi in U
                  for ( int j = 1 ; j < nStepU ; ++ j) {
                    // parametro U
                     double dCoeffU = double( j) / nStepU ;
                     double dU = ( 1 - dCoeffU) * pcToSplit->GetBottomLeft().x + dCoeffU * pcToSplit->GetTopRight().x ;
                    // punto a U, V
                     Point3d ptPSrf ;
                     GetPoint( dU, dV, ptPSrf) ;
                    // distanza del punto dalla precedente linea
                     DistPointCurve dpc( ptPSrf, clU) ;
                     double dDist ;
                     dpc.GetDist( dDist) ;
                     dSagU = max( dSagU, dDist) ;
                  }
               }
               //{ string sLog = " Da Isoparam : FrecciaU=" + ToString( dSagU, 3) + " FrecciaV=" + ToString( dSagV, 3) ;
               //  LOG_DBG_INFO( GetEGkLogger(), sLog.c_str())}
            }
            else if ( dSagU < dLinTol  &&  dSagV < dLinTol  &&  bIsFlat) {
               // se la cella è piatta devo verificare che i bordi siano dei tratti retti, altrimenti potrei commettere un errore di approssimazione
               
               // bordo inferiore e superiore
               double dMaxDist = 0 ;
               for ( int i = 0 ; i < 2 ; ++i) {
                  CurveLine clU ;
                  if ( i == 0)
                     clU.Set( ptP00, ptP10) ;
                  else if ( i == 1)
                     clU.Set( ptP01, ptP11) ;
                  double dV = 0 ;
                  if ( i == 0)
                     dV = pcToSplit->GetBottomLeft().y ;
                  else if ( i == 1)
                     dV = pcToSplit->GetTopRight().y ;


                  int nStepU = 4 ;
                  for ( int j = 1 ; j < nStepU ; ++ j) {
                    // parametro U
                     double dCoeffU = double( j) / nStepU ;
                     double dU = ( 1 - dCoeffU) * pcToSplit->GetBottomLeft().x + dCoeffU * pcToSplit->GetTopRight().x ;
                     Point3d ptBez ;
                     GetPoint( dU, dV, ptBez) ;
                     DistPointCurve dpc( ptBez, clU) ;
                     double dDist = 0 ;
                     dpc.GetDist( dDist) ;
                     if ( dDist > dMaxDist)
                        dMaxDist = dDist ;
                  }
               }
               if ( dMaxDist > dLinTol)
                  dSagU = dMaxDist ;

               // bordo sinistro e destro
               dMaxDist = 0 ;
               for ( int i = 0 ; i < 2 ; ++i) {
                  CurveLine clV ;
                  if ( i == 0)
                     clV.Set( ptP00, ptP01) ;
                  else if ( i == 1)
                     clV.Set( ptP10, ptP11) ;
                  double dU = 0 ;
                  if ( i == 0)
                     dU = pcToSplit->GetBottomLeft().x ;
                  else if ( i == 1)
                     dU = pcToSplit->GetTopRight().x ;

                  int nStepV = 4 ;
                  for ( int j = 1 ; j < nStepV ; ++ j) {
                    // parametro in V
                     double dCoeffV = double( j) / nStepV ;
                     double dV = ( 1 - dCoeffV) * pcToSplit->GetBottomLeft().y + dCoeffV * pcToSplit->GetTopRight().y ;
                     Point3d ptBez ;
                     GetPoint( dU, dV, ptBez) ;
                     DistPointCurve dpc( ptBez, clV) ;
                     double dDist = 0 ;
                     dpc.GetDist( dDist) ;
                     if ( dDist > dMaxDist)
                        dMaxDist = dDist ;
                  }
               }
               if ( dMaxDist > dLinTol)
                  dSagV = dMaxDist ;
            }


           // per lo split scelgo la direzione che è più vicina alla superficie originale nel punto di maggior distanza
           // misura approssimativa della curvatura in una direzione
            bool bVert ;
            if ( max( dLen0, dLen2) > 5 * max( dLen1, dLen3))
               bVert = true ;
            else if ( max( dLen1, dLen3) > 5 * max( dLen0, dLen2))
               bVert = false ;
            else
               bVert = ( dSagV <= dSagU) ;
            bool bFirstTry = true ;
            retry :
            // verifico che la cella sia da splittare e che eventualmente sia abbastanza grande da poterlo fare
            double dSideMinVal = 0 ;
            double dLengMinVal = 0 ;
            if ( bVert) {
               if ( dLen0 > EPS_ZERO  &&  dLen2 > EPS_ZERO)
                  dSideMinVal = min( dLen0, dLen2) ;
               else
                  dSideMinVal = max( dLen0, dLen2) ;
               if ( dLen1 > EPS_ZERO  &&  dLen3 > EPS_ZERO)
                  dLengMinVal = min( dLen1, dLen3) ;
               else
                  dLengMinVal = max( dLen1, dLen3) ;
            }
            else {
               if ( dLen1 > EPS_ZERO  &&  dLen3 > EPS_ZERO)
                  dSideMinVal = min( dLen1, dLen3) ;
               else
                  dSideMinVal = max( dLen1, dLen3) ;
               if ( dLen0 > EPS_ZERO  &&  dLen2 > EPS_ZERO)
                  dLengMinVal = min( dLen0, dLen2) ;
               else
                  dLengMinVal = max( dLen0, dLen2) ;
            }
           // calcolo le diagonali per controllare la dimensione massima dei triangoli in cui dividerei la cella
            double dSideMaxVal = max( Dist( ptP00, ptP11), Dist( ptP10, ptP01)) ;
         
           // se la cella è abbastanza grande da poter essere divisa ancora, calcolo l'errore di approssimazione
           // ( dSideMinVal è zero se entrambi i lati da splittare sono collassati in un punto, controllo dLengMinVal)
            bool bSplit = false ;
            bool bParamDimOk = false ;
            if ( bVert)
               bParamDimOk = ( pcToSplit->GetTopRight().x - pcToSplit->GetBottomLeft().x) / 2 > 100 * EPS_PARAM ;
            else
               bParamDimOk = ( pcToSplit->GetTopRight().y - pcToSplit->GetBottomLeft().y) / 2 > 100 * EPS_PARAM ;
            bool bDimOk = ( dSideMinVal / 2 >= dSideMin  ||  ( dSideMinVal < EPS_SMALL  &&  dLengMinVal / 2 >= dSideMin))  &&  bParamDimOk ;
            if ( dSideMaxVal > dSideMax) {
               bSplit = true ; 
               //LOG_DBG_INFO( GetEGkLogger(), "  Split by SideMax")
            }
            else if ( dSagV > dLinTol  ||  dSagU > dLinTol) {
               bSplit = bDimOk ;
               if ( ! bSplit  &&  bFirstTry) {
                  bFirstTry = false ;
                  bVert = ! bVert ;
                  goto retry ;
               }
               //if ( bSplit)
               //   LOG_DBG_INFO( GetEGkLogger(), "  Split by SagittaUV")
            }

            if ( bSplit) {
               pcToSplit->SetSplitDirVert( bVert) ;
              // effettuo lo split
               if ( Split( nCToSplit)) {
                 // procedo con lo split del Child1
                  nCToSplit = pcToSplit->m_nChild1 ;
                  pcToSplit = &m_mTree[nCToSplit] ;
               }
               else
                  return false ;
            }
            else {
              // sono arrivato ad una cella Leaf, quindi salvo la cella
               m_vnLeaves.push_back( nCToSplit) ;
               pcToSplit->SetProcessed() ;
               if ( AreSamePointApprox( ptP00, ptP01)  &&  AreSamePointApprox( ptP10, ptP11))
                  pcToSplit->m_nCollapsed = Cell::Collapsed::VERT_EDGES ;
               else if ( AreSamePointApprox( ptP00, ptP10)  &&  AreSamePointApprox( ptP01, ptP11))
                  pcToSplit->m_nCollapsed = Cell::Collapsed::HORIZ_EDGES ;
               else
                  pcToSplit->m_nCollapsed = Cell::Collapsed::NO_COLLAPSE ;

              // risalgo i parent finché non trovo il primo Child2 da processare
               nCToSplit = pcToSplit->m_nParent ;
               pcToSplit = &m_mTree[nCToSplit] ;
               if ( nCToSplit == -2)
                  return true ;
               if ( m_mTree[pcToSplit->m_nChild1].IsProcessed()  &&  m_mTree[pcToSplit->m_nChild2].IsProcessed())
                  pcToSplit->SetProcessed() ;
               while ( m_mTree[pcToSplit->m_nChild2].IsProcessed()) {
                  if ( pcToSplit->m_nParent != -2) {
                     nCToSplit = pcToSplit->m_nParent ;
                     pcToSplit = &m_mTree[nCToSplit] ;
                  }
                  if ( m_mTree[pcToSplit->m_nChild1].IsProcessed()  &&  m_mTree[pcToSplit->m_nChild2].IsProcessed())
                     pcToSplit->SetProcessed() ;
                  if ( nCToSplit == -1  &&  m_mTree[pcToSplit->m_nChild2].IsProcessed())
                        break ;
               }
               nCToSplit = pcToSplit->m_nChild2 ;
               pcToSplit = &m_mTree[nCToSplit] ;
            }
         }
         else {
            nCToSplit = pcToSplit->m_nChild1 ;
            pcToSplit = &m_mTree[nCToSplit] ;
         }
      }
   }

  // bilineare
   else {
      while ( nCToSplit != -2  &&  pcToSplit->IsProcessed() == false) {
         if ( pcToSplit->IsLeaf()) {
           // vertici della cella
            Point3d ptP00, ptP10, ptP11, ptP01 ;
            GetPoint(pcToSplit->GetBottomLeft().x, pcToSplit->GetBottomLeft().y, ptP00) ;
            GetPoint(pcToSplit->GetTopRight().x, pcToSplit->GetBottomLeft().y, ptP10) ;
            GetPoint(pcToSplit->GetTopRight().x, pcToSplit->GetTopRight().y, ptP11) ;
            GetPoint(pcToSplit->GetBottomLeft().x, pcToSplit->GetTopRight().y, ptP01) ;
           // distanza reale tra i vertici della cella
            double dLen0 = Dist( ptP00, ptP10) ;
            double dLen1 = Dist( ptP10, ptP11) ;
            double dLen2 = Dist( ptP01, ptP11) ;
            double dLen3 = Dist( ptP00, ptP01) ;
               
            bool bVert = false ;
           // per capire in quale direzione splittare devo guardare quale coppia di lati opposti è più sghemba
            Vector3d vtU0 = ptP01 - ptP00 ;
            Vector3d vtU1 = ptP11 - ptP10 ;
            double dLU0, dFU0, dTU0, dLU1, dFU1, dTU1 ;
            vtU0.ToSpherical( &dLU0, &dFU0, &dTU0) ;
            vtU1.ToSpherical( &dLU1, &dFU1, &dTU1) ;
            double dSkewnessV = abs( dFU1 - dFU0) + abs( dTU1 - dTU0) ;
            Vector3d vtV0 = ptP10 - ptP00 ;
            Vector3d vtV1 = ptP11 - ptP01 ;
            double dLV0, dFV0, dTV0, dLV1, dFV1, dTV1 ;
            vtV0.ToSpherical( &dLV0, &dFV0, &dTV0) ;
            vtV1.ToSpherical( &dLV1, &dFV1, &dTV1) ;
            double dSkewnessU = abs( dFV1 - dFV0) + abs( dTV1 - dTV0) ;
            if ( dSkewnessU > dSkewnessV)
               bVert = false ;
            else
               bVert = true ;
            
           // verifico che la cella sia abbastanza grande da poter essere splittata
            double dSideMinVal = ( bVert ? max( dLen0, dLen2) : max( dLen1, dLen3)) ;
           // calcolo le diagonali per controllare la dimensione massima dei triangoli in cui dividerei la cella
            double dSideMaxVal = max( Dist( ptP00, ptP11), Dist( ptP10, ptP01)) ;
           
            double dErr = 0 ;
            if ( m_bMulti) {
               Point3d ptPSrf ;
               Plane3d plAppr ;
               if ( ! AreSamePointApprox( ptP00, ptP10)  &&  ! AreSamePointApprox( ptP00, ptP01))
                  plAppr.Set( ptP00, ( ptP00 - ptP01) ^ ( ptP00 - ptP10)) ;
               else if ( AreSamePointApprox( ptP00, ptP10)) {
                  plAppr.Set( ptP01, ( ptP00 - ptP01) ^ ( ptP01 - ptP11)) ;
               }
               else if ( AreSamePointApprox( ptP00, ptP01)) {
                  plAppr.Set( ptP10, ( ptP10 - ptP11) ^ ( ptP00 - ptP10)) ;
               }
               for ( double i = 0.25 ; i < 1 ; i = i + 0.25) {
                  for ( double j = 0.25 ; j < 1 ; j = j + 0.25) {
                     double dU =  ( ( 1 - i) * pcToSplit->GetTopRight().x + i * pcToSplit->GetBottomLeft().x) ;
                     double dV =  ( ( 1 - j) * pcToSplit->GetTopRight().y + j * pcToSplit->GetBottomLeft().y) ;
                     GetPoint( dU, dV, ptPSrf) ;
                     dErr = max( abs( DistPointPlane( ptPSrf, plAppr)), dErr) ;
                  }
               }
            }
            else if ( ! bIsPlanar) {
               dErr = ( ( ptP00 - ptP01) + ( ptP11 - ptP10)).Len() / 4 ;
            }
           // se la cella è abbastanza grande da poter essere divisa ancora e devo approssimare meglio, la divido
            if ( dSideMinVal / 2 >= dSideMin  &&  dSideMaxVal < dSideMax  &&  dErr > dLinTol) {
               pcToSplit->SetSplitDirVert( bVert) ;
              // effettuo lo split
               if ( Split( nCToSplit)){
                  // procedo con lo split del Child1
                  nCToSplit = pcToSplit->m_nChild1 ;
                  pcToSplit = &m_mTree[nCToSplit] ;
               }
               else
                  return false ;
            }
            else {
              // sono arrivato ad una cella Leaf, quindi salvo la cella
               pcToSplit->SetProcessed() ;
               if ( AreSamePointApprox( ptP00, ptP01)  &&  AreSamePointApprox( ptP10, ptP11))
                  pcToSplit->m_nCollapsed = Cell::Collapsed::VERT_EDGES ;
               else if ( AreSamePointApprox( ptP00, ptP10)  &&  AreSamePointApprox( ptP01, ptP11))
                  pcToSplit->m_nCollapsed = Cell::Collapsed::HORIZ_EDGES ;
               else {
                  pcToSplit->m_nCollapsed = Cell::Collapsed::NO_COLLAPSE ;
                  m_vnLeaves.push_back( nCToSplit) ;
               }

              // risalgo i parent finché non trovo il primo Child2 da processare
               nCToSplit = pcToSplit->m_nParent ;
               pcToSplit = &m_mTree[nCToSplit] ;
               if ( nCToSplit == -2)
                  return true ;
               if ( m_mTree[pcToSplit->m_nChild1].IsProcessed()  &&  m_mTree[pcToSplit->m_nChild2].IsProcessed())
                  pcToSplit->SetProcessed() ;
               while ( m_mTree[pcToSplit->m_nChild2].IsProcessed()) {
                  if ( pcToSplit->m_nParent != -2) {
                     nCToSplit = pcToSplit->m_nParent ;
                     pcToSplit = &m_mTree[nCToSplit] ;
                  }
                  if ( m_mTree[pcToSplit->m_nChild1].IsProcessed()  &&  m_mTree[pcToSplit->m_nChild2].IsProcessed())
                     pcToSplit->SetProcessed() ;
                  if ( nCToSplit == -1  &&  m_mTree[pcToSplit->m_nChild2].IsProcessed())
                     break ;
               }
               nCToSplit = pcToSplit->m_nChild2 ;
               pcToSplit = &m_mTree[nCToSplit] ;
            }
         }
         else {
            nCToSplit = pcToSplit->m_nChild1 ;
            pcToSplit = &m_mTree[nCToSplit] ;
         }
      }
   }

   return true ;
}

//----------------------------------------------------------------------------
void
Tree::GetTopNeigh( int nId, INTVECTOR& vTopNeighs, DBLDBL ddInt) const
{
   const Cell& cell = m_mTree.at( nId) ;
   double dMax = cell.GetTopRight().x ;
   double dMin = cell.GetBottomLeft().x ;
   if ( ddInt.second > 0) {
      dMax = ddInt.second ;
      dMin = ddInt.first ;
   }

   // le celle restituite sono ordinate per x crescente
   if ( vTopNeighs.empty()) {
      if ( cell.m_nTop == -2)
         return ;
      if ( m_mTree.at( cell.m_nTop).IsLeaf())
         vTopNeighs.push_back( cell.m_nTop) ;
      else {
         if ( m_mTree.at( cell.m_nTop).IsSplitVert()) {
           // se la cella vicina è più piccola della cella indagata, allora entrambi i figli saranno vicini di quest'ultima
            if ( m_mTree.at( cell.m_nTop).GetTopRight().x - m_mTree.at( cell.m_nTop).GetBottomLeft().x <= 
               dMax - dMin) {
               vTopNeighs.push_back( m_mTree.at( cell.m_nTop).m_nChild1) ;
               vTopNeighs.push_back( m_mTree.at( cell.m_nTop).m_nChild2) ;
            }
           // altrimenti solo uno dei figli lo sarà
            else {
               if ( m_mTree.at( m_mTree.at( cell.m_nTop).m_nChild1).GetTopRight().x <= dMin  ||
                    m_mTree.at( m_mTree.at( cell.m_nTop).m_nChild1).GetBottomLeft().x >= dMax )
                  vTopNeighs.push_back( m_mTree.at( cell.m_nTop).m_nChild2) ;
               else
                  vTopNeighs.push_back( m_mTree.at( cell.m_nTop).m_nChild1) ;
            }
         }
         else {
            vTopNeighs.push_back( m_mTree.at( cell.m_nTop).m_nChild2) ;
         }
      }
      bool bAllLeaves = true ;
      for ( int i : vTopNeighs) {
         if ( ! m_mTree.at( i).IsLeaf())
            bAllLeaves = false ;
      }
      if ( ! bAllLeaves)
        // almeno una cella tra i vicini trovati non è leaf quindi devo richiamare ricorsivamente questa funzione per trovare i suoi child
         GetTopNeigh( nId, vTopNeighs) ;
   }
   else {
      for ( int j = 0 ; j != (int) vTopNeighs.size() ; ++ j) {
         int i = vTopNeighs.at( j) ;
         if ( m_mTree.at( i).IsLeaf())
            continue ;
         else {
           // se la cella non è leaf la tolgo dal vettore delle foglie e aggiungo invece i suoi child
            vTopNeighs.erase( remove( vTopNeighs.begin(),vTopNeighs.end(),i)) ;
            -- j ;
            if ( m_mTree.at( i).IsSplitVert()) {
              // se la cella è più piccola della cella indagata, allora entrambi i figli saranno vicini di quest'ultima
               if ( m_mTree.at( i).GetTopRight().x - m_mTree.at( i).GetBottomLeft().x <=
                  dMax - dMin) {
                  vTopNeighs.push_back( m_mTree.at( i).m_nChild1) ;
                  vTopNeighs.push_back( m_mTree.at( i).m_nChild2) ;
               }
              // altrimenti solo uno dei figli lo sarà
               else {
                  if ( m_mTree.at( m_mTree.at( i).m_nChild1).GetTopRight().x <= dMin ||
                     m_mTree.at( m_mTree.at( i).m_nChild1).GetBottomLeft().x >= dMax)
                     vTopNeighs.push_back( m_mTree.at( i).m_nChild2) ;
                  else
                     vTopNeighs.push_back( m_mTree.at( i).m_nChild1) ;
               }
            }
            else {
               vTopNeighs.push_back( m_mTree.at( i).m_nChild2) ;
            }
         }
      }
   }
   vector<Cell> vCells ;
   for ( int k : vTopNeighs)
      vCells.push_back( m_mTree.at( k)) ;
   // le celle restituite sono ordinate per x crescente
   std::sort( vCells.begin(), vCells.end(), Cell::minorX) ;
   vTopNeighs.clear() ;
   for ( Cell& c : vCells)
      vTopNeighs.push_back( c.m_nId) ;
}


//----------------------------------------------------------------------------
void
Tree::GetBottomNeigh( int nId, INTVECTOR& vBottomNeighs, DBLDBL ddInt) const
{
   const Cell& cell = m_mTree.at( nId) ;
   double dMax = cell.GetTopRight().x ;
   double dMin = cell.GetBottomLeft().x ;
   if ( ddInt.second > 0) {
      dMax = ddInt.second ;
      dMin = ddInt.first ;
   }

   // le celle restituite sono ordinate per x crescente
   if ( vBottomNeighs.empty()) {
      if ( cell.m_nBottom == -2)
         return ;
      if ( m_mTree.at( cell.m_nBottom).IsLeaf())
         vBottomNeighs.push_back( cell.m_nBottom) ;
      else {
         if ( m_mTree.at( cell.m_nBottom).IsSplitVert()) {
           // se la cella vicina è più piccola della cella indagata, allora entrambi i figli saranno vicini di quest'ultima
            if ( m_mTree.at( cell.m_nBottom).GetTopRight().x - m_mTree.at( cell.m_nBottom).GetBottomLeft().x <= 
               dMax - dMin) {
               vBottomNeighs.push_back( m_mTree.at( cell.m_nBottom).m_nChild1) ;
               vBottomNeighs.push_back( m_mTree.at( cell.m_nBottom).m_nChild2) ;
            }
           // altrimenti solo uno dei figli lo sarà
            else{
               if ( m_mTree.at( m_mTree.at( cell.m_nBottom).m_nChild1).GetTopRight().x <= dMin ||
                  m_mTree.at( m_mTree.at( cell.m_nBottom).m_nChild1).GetBottomLeft().x >= dMax )
                  vBottomNeighs.push_back( m_mTree.at( cell.m_nBottom).m_nChild2) ;
               else
                  vBottomNeighs.push_back( m_mTree.at( cell.m_nBottom).m_nChild1) ;
            }
         }
         else {
            vBottomNeighs.push_back( m_mTree.at( cell.m_nBottom).m_nChild1) ;
         }
      }
      bool bAllLeaves = true ;
      for ( int i : vBottomNeighs) {
         if ( ! m_mTree.at( i).IsLeaf())
            bAllLeaves = false ;
      }
      if ( ! bAllLeaves)
        // almeno una cella tra i vicini trovati non è leaf quindi devo richiamare ricorsivamente questa funzione per trovare i suoi child
         GetBottomNeigh( nId, vBottomNeighs) ;
   }
   else {
      for ( int j = 0 ; j < int( vBottomNeighs.size()) ; ++ j) {
         int i = vBottomNeighs.at( j) ;
         if ( m_mTree.at( i).IsLeaf())
            continue ;
         else {
           // se la cella non è leaf la tolgo dal vettore delle foglie e aggiungo invece i suoi child
            vBottomNeighs.erase( remove( vBottomNeighs.begin(),vBottomNeighs.end(),i)) ;
            -- j ;
            if ( m_mTree.at( i).IsSplitVert()) {
              // se la cella è più piccola della cella indagata, allora entrambi i figli saranno vicini di quest'ultima
               if ( m_mTree.at( i).GetTopRight().x - m_mTree.at( i).GetBottomLeft().x <=
                  dMax - dMin) {
                  vBottomNeighs.push_back( m_mTree.at( i).m_nChild1) ;
                  vBottomNeighs.push_back( m_mTree.at( i).m_nChild2) ;
               }
              // altrimenti solo uno dei figli lo sarà
               else {
                  if ( m_mTree.at( m_mTree.at( i).m_nChild1).GetTopRight().x <= dMin  ||
                     m_mTree.at( m_mTree.at( i).m_nChild1).GetBottomLeft().x >= dMax)
                     vBottomNeighs.push_back( m_mTree.at( i).m_nChild2) ;
                  else
                     vBottomNeighs.push_back( m_mTree.at( i).m_nChild1) ;
               }
            }
            else {
               vBottomNeighs.push_back( m_mTree.at( i).m_nChild1) ;
            }
         }
      }
   }
   vector<Cell> vCells ;
   for ( int k : vBottomNeighs)
      vCells.push_back( m_mTree.at( k)) ;
   // le celle restituite sono ordinate per x crescente
   std::sort( vCells.begin(), vCells.end(), Cell::minorX) ;
   vBottomNeighs.clear() ;
   for ( Cell& c : vCells)
      vBottomNeighs.push_back( c.m_nId) ;
}

//----------------------------------------------------------------------------
void
Tree::GetLeftNeigh( int nId, INTVECTOR& vLeftNeighs, DBLDBL ddInt) const
{
   const Cell& cell = m_mTree.at( nId) ;
   double dMax = cell.GetTopRight().y ;
   double dMin = cell.GetBottomLeft().y ;
   if ( ddInt.second > 0) {
      dMax = ddInt.second ;
      dMin = ddInt.first ;
   }

   // le celle restituite sono ordinate per y crescente
   if ( vLeftNeighs.empty()) {
      if ( cell.m_nLeft == -2)
         return ;
      if ( m_mTree.at( cell.m_nLeft).IsLeaf())
         vLeftNeighs.push_back( cell.m_nLeft) ;
      else {
            if ( ! m_mTree.at( cell.m_nLeft).IsSplitVert()) {
              // se la cella vicina è più piccola della cella indagata, allora entrambi i figli saranno vicini di quest'ultima
               if ( m_mTree.at( cell.m_nLeft).GetTopRight().y - m_mTree.at( cell.m_nLeft).GetBottomLeft().y <= 
                    dMax - dMin) {
                  vLeftNeighs.push_back( m_mTree.at( cell.m_nLeft).m_nChild1) ;
                  vLeftNeighs.push_back( m_mTree.at( cell.m_nLeft).m_nChild2) ;
               }
              // altrimenti solo uno dei figli lo sarà
               else{
                  if ( m_mTree.at( m_mTree.at( cell.m_nLeft).m_nChild1).GetTopRight().y <= dMin  ||
                       m_mTree.at( m_mTree.at( cell.m_nLeft).m_nChild1).GetBottomLeft().y >= dMax)
                     vLeftNeighs.push_back( m_mTree.at( cell.m_nLeft).m_nChild2) ;
                  else
                     vLeftNeighs.push_back( m_mTree.at( cell.m_nLeft).m_nChild1) ;
               }
            }
            else {
               vLeftNeighs.push_back( m_mTree.at( cell.m_nLeft).m_nChild2) ;
            }
      }
      bool bAllLeaves = true ;
      for ( int i : vLeftNeighs) {
         if ( ! m_mTree.at( i).IsLeaf())
            bAllLeaves = false ;
      }
      if ( ! bAllLeaves)
        // almeno una cella tra i vicini trovati non è leaf quindi devo richiamare ricorsivamente questa funzione per trovare i suoi child
         GetLeftNeigh( nId, vLeftNeighs) ;
   }
   else {
      for (int j = 0 ; j < int( vLeftNeighs.size()) ; ++ j) {
         int i = vLeftNeighs.at( j) ;
         if ( m_mTree.at( i).IsLeaf())
            continue ;
         else {
           // se la cella non è leaf la tolgo dal vettore delle foglie e aggiungo invece i suoi child
            vLeftNeighs.erase( remove( vLeftNeighs.begin(),vLeftNeighs.end(),i)) ;
            -- j ;
            if ( ! m_mTree.at( i).IsSplitVert()) {
              // se la cella è più piccola della cella indagata, allora entrambi i figli saranno vicini di quest'ultima
               if ( m_mTree.at( i).GetTopRight().y - m_mTree.at( i).GetBottomLeft().y <=
                    dMax - dMin) {
                  vLeftNeighs.push_back( m_mTree.at( i).m_nChild1) ;
                  vLeftNeighs.push_back( m_mTree.at( i).m_nChild2) ;
               }
              // altrimenti solo uno dei figli lo sarà
               else {
                  if ( m_mTree.at( m_mTree.at( i).m_nChild1).GetTopRight().y <= dMin  ||
                       m_mTree.at( m_mTree.at( i).m_nChild1).GetBottomLeft().y >= dMax)
                     vLeftNeighs.push_back( m_mTree.at( i).m_nChild2) ;
                  else
                     vLeftNeighs.push_back( m_mTree.at( i).m_nChild1) ;
               }
            }
            else {
               vLeftNeighs.push_back( m_mTree.at( i).m_nChild2) ;
            }
         }
      }
   }
   vector<Cell> vCells ;
   for ( int k : vLeftNeighs)
      vCells.push_back( m_mTree.at( k)) ;
   // le celle restituite sono ordinate per y crescente
   std::sort( vCells.begin(), vCells.end(), Cell::minorY) ;
   vLeftNeighs.clear() ;
   for ( Cell& c : vCells)
      vLeftNeighs.push_back( c.m_nId) ;
}

//----------------------------------------------------------------------------
void
Tree::GetRightNeigh( int nId, INTVECTOR& vRightNeighs, DBLDBL ddInt) const
{
   const Cell& cell = m_mTree.at( nId) ;
   double dMax = cell.GetTopRight().y ;
   double dMin = cell.GetBottomLeft().y ;
   if ( ddInt.second > 0) {
      dMax = ddInt.second ;
      dMin = ddInt.first ;
   }

   // le celle restituite sono ordinate per y crescente
   if ( vRightNeighs.empty()) {
      if ( cell.m_nRight == -2)
         return ;
      if ( m_mTree.at( cell.m_nRight).IsLeaf())
         vRightNeighs.push_back( cell.m_nRight) ;
      else {
         if ( ! m_mTree.at( cell.m_nRight).IsSplitVert()) {
           // se la cella vicina è più piccola della cella indagata, allora entrambi i figli saranno vicini di quest'ultima
            if ( m_mTree.at( cell.m_nRight).GetTopRight().y - m_mTree.at( cell.m_nRight).GetBottomLeft().y <= 
                 dMax - dMin) {
               vRightNeighs.push_back( m_mTree.at( cell.m_nRight).m_nChild1) ;
               vRightNeighs.push_back( m_mTree.at( cell.m_nRight).m_nChild2) ;
            }
           // altrimenti solo uno dei figli lo sarà
            else{
               if ( m_mTree.at( m_mTree.at( cell.m_nRight).m_nChild1).GetTopRight().y <= dMin  ||
                    m_mTree.at( m_mTree.at( cell.m_nRight).m_nChild1).GetBottomLeft().y >= dMax)
                  vRightNeighs.push_back( m_mTree.at( cell.m_nRight).m_nChild2) ;
               else
                  vRightNeighs.push_back( m_mTree.at( cell.m_nRight).m_nChild1) ;
            }
         }
         else {
            vRightNeighs.push_back( m_mTree.at( cell.m_nRight).m_nChild1) ;
         }
      }
      bool bAllLeaves = true ;
      for ( int i : vRightNeighs) {
         if ( ! m_mTree.at( i).IsLeaf())
            bAllLeaves = false ;
      }
      if ( ! bAllLeaves)
        // almeno una cella tra i vicini trovati non è leaf quindi devo richiamare ricorsivamente questa funzione per trovare i suoi child
         GetRightNeigh( nId, vRightNeighs) ;
   }
   else {
      for (int j = 0 ; j < int( vRightNeighs.size()) ; ++ j) {
         int i = vRightNeighs.at( j) ;
         if ( m_mTree.at( i).IsLeaf())
            continue ;
         else {
           // se la cella non è leaf la tolgo dal vettore delle foglie e aggiungo invece i suoi child
            vRightNeighs.erase( remove( vRightNeighs.begin(),vRightNeighs.end(), i)) ;
            -- j ;
            if ( ! m_mTree.at( i).IsSplitVert()) {
              // se la cella è più piccola della cella indagata, allora entrambi i figli saranno vicini di quest'ultima
               if ( m_mTree.at( i).GetTopRight().y - m_mTree.at( i).GetBottomLeft().y <=
                    dMax - dMin) {
                  vRightNeighs.push_back( m_mTree.at( i).m_nChild1) ;
                  vRightNeighs.push_back( m_mTree.at( i).m_nChild2) ;
               }
              // altrimenti solo uno dei figli lo sarà
               else {
                  if ( m_mTree.at( m_mTree.at( i).m_nChild1).GetTopRight().y <= dMin  ||
                       m_mTree.at( m_mTree.at( i).m_nChild1).GetBottomLeft().y >= dMax)
                     vRightNeighs.push_back( m_mTree.at( i).m_nChild2) ;
                  else
                     vRightNeighs.push_back( m_mTree.at( i).m_nChild1) ;
               }
            }
            else {
               vRightNeighs.push_back( m_mTree.at( i).m_nChild1) ;
            }
         }
      }
   }
   vector<Cell> vCells ;
   for ( int k : vRightNeighs)
      vCells.push_back( m_mTree.at( k)) ;
   // le celle restituite sono ordinate per y crescente
   std::sort( vCells.begin(), vCells.end(), Cell::minorY) ;
   vRightNeighs.clear() ;
   for ( Cell& c : vCells)
      vRightNeighs.push_back( c.m_nId) ;
}

//----------------------------------------------------------------------------
void
Tree::GetRootNeigh( int nEdge, INTVECTOR& vNeigh)
{
   int nId = -1 ;
   bool bMod = false ;
   if ( nEdge == 0) {
      if ( m_mTree[nId].m_nBottom == -2) {
         m_mTree[nId].m_nBottom = -1 ;
         bMod = true ;
      }
      GetBottomNeigh( nId, vNeigh) ;
      if ( bMod)
         m_mTree[nId].m_nBottom = -2 ;
   }
   else if ( nEdge == 1) {
      if ( m_mTree[nId].m_nRight == -2) {
         m_mTree[nId].m_nRight = -1 ;
         bMod = true ;
      }
      GetRightNeigh( nId, vNeigh) ;
      if ( bMod)
         m_mTree[nId].m_nRight = -2 ;
   }
   else if ( nEdge == 2) {
      if ( m_mTree[nId].m_nTop == -2) {
         m_mTree[nId].m_nTop = -1 ;
         bMod = true ;
      }
      GetTopNeigh( nId, vNeigh) ;
      if ( bMod)
         m_mTree[nId].m_nTop = -2 ;
   }
   else if ( nEdge == 3) {
      if ( m_mTree[nId].m_nLeft == -2) {
         m_mTree[nId].m_nLeft = -1 ;
         bMod = true ;
      }
      GetLeftNeigh( nId, vNeigh) ;
      if ( bMod)
         m_mTree[nId].m_nLeft = -2 ;
   } 
}

//----------------------------------------------------------------------------
int
Tree::GetHeightLeaves( int nId, INTVECTOR& vnLeaves, int d) const
{
   if ( nId == -1  &&  m_mTree.at( -1).IsLeaf()) {
      vnLeaves.push_back( -1) ;
      return 0 ;
   }
   else {
      if ( vnLeaves.empty()) {
         if ( m_mTree.at( nId).IsLeaf())
            return d ;
         else {
            vnLeaves.push_back( m_mTree.at( nId).m_nChild1) ;
            vnLeaves.push_back( m_mTree.at( nId).m_nChild2) ;
            if ( ! m_mTree.at( m_mTree.at( nId).m_nChild1).IsLeaf()  ||  ! m_mTree.at( m_mTree.at( nId).m_nChild2).IsLeaf())
              // almeno un child non è leaf quindi devo richiamare ricorsivamente questa funzione sui child in questione
               d = GetHeightLeaves( nId, vnLeaves, m_mTree.at( m_mTree.at( nId).m_nChild1).m_nDepth) ;
         }
      }
      else {
         for ( int j = 0 ; j < int( vnLeaves.size()) ; ++ j) {
            int i = vnLeaves.at( j) ;
            if ( m_mTree.at( i).IsLeaf()) {
               continue ;
            }
            else {
              // se la cella non è leaf la tolgo dal vettore delle foglie e aggiungo invece i suoi child
               vnLeaves.erase( remove( vnLeaves.begin(),vnLeaves.end(),i)) ;
               -- j ;
               vnLeaves.push_back( m_mTree.at( i).m_nChild1) ;
               vnLeaves.push_back( m_mTree.at( i).m_nChild2) ;
               d = max( d, m_mTree.at( m_mTree.at( i).m_nChild1).m_nDepth) ;
            }
         }
         return d ;
      }
      return ( d - m_mTree.at( nId).m_nDepth) ;
   }
}

//----------------------------------------------------------------------------
int
Tree::GetDepth( int nId, int nRef = -2) const 
{
   int i = 0 ;
   int nCell = nId ;
   while ( m_mTree.at( nCell).m_nParent != nRef) {
      nCell = m_mTree.at( nCell).m_nParent ;
      ++ i ;
   }
   return i ;
}

//----------------------------------------------------------------------------
bool
Tree::GetPolygons( POLYLINEMATRIX& vvPolygons, POLYLINEMATRIX& vvPolygons3d, vector<ICRVCOMPOPOVECTOR>& vvCCEdges3D, ICRVCOMPOPOVECTOR& vCCLoops, bool bUpdateEdges)
{
   POLYLINEVECTOR vPolygonsBasic ;
   POLYLINEVECTOR vPolygonsCorrected ;
   POLYLINEVECTOR vPolygonsBasic3d ;
   if ( ! m_bTrimmed) {
      vvPolygons.clear() ;
         
      GetPolygonsBasic( vPolygonsBasic, vPolygonsCorrected, vPolygonsBasic3d) ;
      if (vPolygonsBasic.empty())
         return false ;
      int c = 0;
      for ( PolyLine& pl : vPolygonsCorrected) {
         vvPolygons.emplace_back() ;
         vvPolygons.back().emplace_back(std::move(pl)) ;
         vvPolygons3d.emplace_back() ;
         vvPolygons3d.back().emplace_back(std::move(vPolygonsBasic3d[c])) ;
         ++c ;
      }
   }
   // trimmata
   else {
      GetPolygonsBasic( vPolygonsBasic, vPolygonsCorrected, vPolygonsBasic3d) ;
      if (vPolygonsBasic.empty())
         return false ;
      // aggiungo 4 elementi al vettore che contiene ciò che resta degli edge dopo il trim
      m_vCEdge2D.clear() ;
      for ( int i = 0 ; i < 4  ; ++i) {
         m_vCEdge2D.emplace_back() ;
         m_vCEdge2D.back().second.Init( false, EPS_SMALL, 1) ;
      }
      // percorro i loop, trovo le intersezioni con le celle e le categorizzo
      if ( ! TraceLoopLabelCell( vPolygonsBasic))
         return false ;
      // scorro sulle celle e costruisco i poligoni
      int nPolyInd = 0 ; // indice del poligono corrispondente alla cella, nel vettore dei poligoni
      for ( int nId: m_vnLeaves) {
         // costruisco i poligoni partendo dal vettore delle intersezioni, come spiegato a pag15 di Cripps
         if ( m_mTree[nId].m_nFlag == 4) {
            // vettore dei poligoni ( loop) della cella nId
            vvPolygons.emplace_back() ;
            vvPolygons.back().push_back(std::move(vPolygonsCorrected[nPolyInd])) ;
            vvPolygons3d.emplace_back() ;
            vvPolygons3d.back().push_back(std::move(vPolygonsBasic3d[nPolyInd])) ;
         }
         else if ( m_mTree[nId].m_nFlag == 0) {
            ++ nPolyInd ;
            continue ;
         }
         else if ( m_mTree[nId].m_nCollapsed != Cell::Collapsed::NO_COLLAPSE)
            continue ;
         else {
            // vettore in cui salvo il chunk di appartenenza di ogni loop che attraversa la cella
            INTVECTOR vnParentChunk ;
            // vettore in cui salvo i loop che non appartengono al poligono che sto costruendo nel ciclo attuale e da cui ripasserò dopo
            INTVECTOR vToCheck( m_mTree[nId].m_vInters.size()) ;
            iota( vToCheck.begin(), vToCheck.end(), 0) ;
            // numero di poligoni aggiunti
            int nPoly = 0 ;
            // scorro sui vettori intersezione della cella nId e sui suoi vertici
            // in questo for analizzo solo i loop che tagliano la cella
            while ( ! vToCheck.empty()) {
               int nPolyBefore = nPoly ;
               PolyLine pl3d ;
               CreateCellPolygons( nId, vvPolygons, vvPolygons3d, vToCheck, nPoly, vnParentChunk, vPolygonsCorrected[nPolyInd], vPolygonsBasic3d[nPolyInd]) ;
               if ( nPolyBefore == nPoly)
                  break ;
            }
            // ora analizzo anche i loop che sono contenuti nella cella
            CreateIslandAndHoles( nId, vvPolygons, vvPolygons3d, nPoly, vnParentChunk, vPolygonsCorrected[nPolyInd], vPolygonsBasic3d[nPolyInd]) ;
         }
         ++ nPolyInd ;
      }
   }

   // se richiesto aggiorno gli edge 3d e i loop della superficie
   if ( bUpdateEdges  ||  vvCCEdges3D.empty()) {
      GetEdges3D( vvCCEdges3D, vPolygonsBasic) ;
      GetSplitLoops( vCCLoops) ;
   }
   return  true ;
}

//----------------------------------------------------------------------------
bool
Tree::GetPolygonsBasic( POLYLINEVECTOR& vPolygonsBasic, POLYLINEVECTOR& vPolygonsCorrected, POLYLINEVECTOR& vPolygons3d)
{
   // condizioni per il calcolo dei poligoni di base
   PNTVECTOR vVertices ;
   PNTVECTOR vVerticesCorr ;
   PNTVECTOR vVertices3d ;

   INTVECTOR vNeigh ;
   // setto le celle che sono sul LeftEdge e sul TopEdge
   if ( m_bClosedU) {
      GetRootNeigh( 1, vNeigh) ;
      for ( int k : vNeigh) {
         m_mTree[k].m_bOnLeftEdge = true ;
      }
      vNeigh.clear() ;
   }
   if ( m_bClosedV) {
      GetRootNeigh( 0, vNeigh) ;
      for ( int k : vNeigh) {
         m_mTree[k].m_bOnTopEdge = true ;
      }
   }
   bool bBottomRight , bTopLeft ;
   // scorro lungo tutte le celle leaves e oltre agli angoli della cella aggiungo alla polyline della cella anche i vertici, delle celle adiacenti, 
   // che sono sui lati della cella corrente
   // N.B. :i poligoni sono costruiti a partire dal ptBL !!!
   int c = 0 ;
   for ( int nId : m_vnLeaves) {
      Cell& cell = m_mTree.at( nId) ;
      if ( cell.m_nCollapsed != Cell::Collapsed::NO_COLLAPSE)
         continue ;
      vVertices.clear() ;
      vVertices3d.clear() ;
      vNeigh.clear() ;
      vVertices.push_back( cell.GetBottomLeft()) ;
      Point3d pt3d ; GetPoint( cell.GetBottomLeft().x, cell.GetBottomLeft().y, pt3d) ;
      vVertices3d.push_back( pt3d) ;
      INTVECTOR vnVert ;
      BOOLVECTOR vbBonusVert(4) ;
      std::fill( vbBonusVert.begin(), vbBonusVert.end(), false) ;
      vnVert.push_back( int(vVertices.size()) - 1) ;
      GetBottomNeigh( nId, vNeigh) ;
      if ( ! vNeigh.empty()  &&  m_mTree.at(vNeigh[0]).m_nCollapsed == Cell::Collapsed::VERT_EDGES)
         GetBottomNeigh( vNeigh[0], vNeigh, DBLDBL( cell.GetBottomLeft().x, cell.GetTopRight().x)) ;
      Point3d ptP00, ptP10, ptP11, ptP01 ;
      GetPoint( cell.GetBottomLeft().x, cell.GetBottomLeft().y, ptP00) ;
      GetPoint( cell.GetTopRight().x, cell.GetBottomLeft().y, ptP10) ;
      GetPoint( cell.GetTopRight().x, cell.GetTopRight().y, ptP11) ;
      GetPoint( cell.GetBottomLeft().x, cell.GetTopRight().y, ptP01) ;

      // aggiungo i vertici che sono sul lato bottom, solo se ho più di un vicino bottom
      if ( vNeigh.size() > 1){
         // se la superficie è chiusa lungo il parametro V e le celle vicine bottom sono sul lato Top
         // devo aggiungere i vertici tenendo conto della periodicità dello spazio parametrico.
         if ( m_bClosedV  &&  m_mTree.at(vNeigh[0]).m_bOnTopEdge) {
            for ( int j : vNeigh) {
               Cell& cNeigh = m_mTree.at(j) ;
               Point3d pt( cNeigh.GetTopRight().x, cell.GetBottomLeft().y) ;
               vVertices.push_back( pt) ;
               Point3d pt3d ; GetPoint( cNeigh.GetTopRight().x, cNeigh.GetTopRight().y, pt3d) ;
               vVertices3d.push_back( pt3d) ;
            }
         }
         else {
            for ( int j : vNeigh) {
               Cell& cNeigh = m_mTree.at(j) ;
               vVertices.push_back( cNeigh.GetTopRight()) ;
               Point3d pt3d ; GetPoint( cNeigh.GetTopRight().x, cNeigh.GetTopRight().y, pt3d) ;
               vVertices3d.push_back( pt3d) ;
            }
         }
         bBottomRight = true ;
         vbBonusVert[2] = true ;
      }
      else
         bBottomRight = false ;
      vNeigh.clear() ;
      GetRightNeigh ( nId, vNeigh) ;
      if ( ! vNeigh.empty()  &&  m_mTree.at(vNeigh[0]).m_nCollapsed == Cell::Collapsed::HORIZ_EDGES)
         GetRightNeigh( vNeigh[0], vNeigh, DBLDBL( cell.GetBottomLeft().y, cell.GetTopRight().y)) ;
      // aggiungo i vertici che sono sul lato right, solo se ho più di un vicino right
      if ( vNeigh.size() > 1){
         // se la superficie è chiusa lungo il parametro U e le celle vicine right sono sul lato Left
         // devo aggiungere i vertici tenendo conto della periodicità dello spazio parametrico.
         vnVert.push_back( int(vVertices.size())) ;
         if ( m_bClosedU  &&  m_mTree.at( vNeigh[0]).m_bOnLeftEdge ) {
            for ( int j : vNeigh) {
               Cell& cNeigh = m_mTree.at(j) ;
               Point3d pt( cell.GetTopRight().x, cNeigh.GetBottomLeft().y) ;
               vVertices.push_back( pt) ;
               Point3d pt3d ; GetPoint( cNeigh.GetBottomLeft().x, cNeigh.GetBottomLeft().y, pt3d) ;
               vVertices3d.push_back( pt3d) ;
            }
         }
         else {
            for ( int j : vNeigh) {
               Cell& cNeigh = m_mTree.at(j) ;
               vVertices.push_back( cNeigh.GetBottomLeft()) ;
               Point3d pt3d ; GetPoint( cNeigh.GetBottomLeft().x, cNeigh.GetBottomLeft().y, pt3d) ;
               vVertices3d.push_back( pt3d) ;
            }
         }
         vbBonusVert[3] = true ;
      }
      // se non l'ho già aggiunto tramite i vicini bottom aggiungo il punto bottom right
      else if ( ! bBottomRight) {  
         Point3d ptBr( cell.GetTopRight().x, cell.GetBottomLeft().y) ;
         vVertices.push_back( ptBr) ;
         Point3d pt3d ; GetPoint( cell.GetTopRight().x, cell.GetBottomLeft().y, pt3d) ;
         vVertices3d.push_back( pt3d) ;
         vnVert.push_back( int(vVertices.size()) - 1) ;
      }
      else
         vnVert.push_back( int(vVertices.size()) - 1) ;  // ho già aggiunto il punto bottom right, ne salvo la posizione
      vNeigh.clear() ;
      vVertices.push_back( cell.GetTopRight()) ;
      pt3d = ORIG ; GetPoint( cell.GetTopRight().x, cell.GetTopRight().y, pt3d) ;
      vVertices3d.push_back( pt3d) ;
      vnVert.push_back( int(vVertices.size()) - 1) ;
      GetTopNeigh ( nId, vNeigh) ;
      if ( ! vNeigh.empty()  &&  m_mTree.at(vNeigh[0]).m_nCollapsed == Cell::Collapsed::VERT_EDGES)
         GetTopNeigh( vNeigh[0], vNeigh, DBLDBL( cell.GetBottomLeft().x, cell.GetTopRight().x)) ;
      std::reverse( vNeigh.begin(), vNeigh.end()) ;
      // aggiungo i vertici che sono sul lato top, solo se ho più di un vicino top
      if ( vNeigh.size() > 1) {
         // se la superficie è chiusa lungo il parametro V e la cella è sul lato top
         // devo aggiungere i vertici tenendo conto della periodicità dello spazio parametrico.
         vnVert.push_back( int(vVertices.size())) ;
         if ( m_bClosedV  &&  cell.m_bOnTopEdge) {
            for ( int j : vNeigh) {
               Point3d pt( m_mTree.at( j).GetBottomLeft().x, cell.GetTopRight().y) ;
               vVertices.push_back( pt) ;
               Cell& cNeigh = m_mTree.at(j) ;
               Point3d pt3d ; GetPoint( cNeigh.GetBottomLeft().x, cNeigh.GetBottomLeft().y, pt3d) ;
               vVertices3d.push_back( pt3d) ;
            }
         }
         else {
            for ( int j : vNeigh) {
               vVertices.push_back( m_mTree.at( j).GetBottomLeft()) ;
               Cell& cNeigh = m_mTree.at(j) ;
               Point3d pt3d ; GetPoint( cNeigh.GetBottomLeft().x, cNeigh.GetBottomLeft().y, pt3d) ;
               vVertices3d.push_back( pt3d) ;
            }
         }
         bTopLeft = true ;
         vbBonusVert[0] = true ;
      }
      else
         bTopLeft = false ;
      vNeigh.clear() ;
      GetLeftNeigh ( nId, vNeigh) ;
      if ( ! vNeigh.empty()  &&  m_mTree.at(vNeigh[0]).m_nCollapsed == Cell::Collapsed::HORIZ_EDGES)
         GetLeftNeigh( vNeigh[0], vNeigh, DBLDBL( cell.GetBottomLeft().y, cell.GetTopRight().y)) ;
      std::reverse( vNeigh.begin(), vNeigh.end()) ;
      // aggiungo i vertici che sono sul lato left, solo se ho più di un vicino left
      if ( vNeigh.size() > 1) {
         // se la superficie è chiusa lungo il parametro U e la cella è sul lato left
         // devo aggiungere i vertici tenendo conto della periodicità dello spazio parametrico.
         vnVert.push_back( int(vVertices.size())) ;
         if ( m_bClosedU  &&  cell.m_bOnLeftEdge) {
            for ( int j : vNeigh) {
               Cell& cNeigh = m_mTree.at(j) ;
               Point3d pt( cell.GetBottomLeft().x, cNeigh.GetTopRight().y) ;
               vVertices.push_back( pt) ;
               Point3d pt3d ; GetPoint( cNeigh.GetTopRight().x, cNeigh.GetTopRight().y, pt3d) ;
               vVertices3d.push_back( pt3d) ;
            }
         }
         else {
            for ( int j : vNeigh) {
               Cell& cNeigh = m_mTree.at(j) ;
               vVertices.push_back( cNeigh.GetTopRight()) ;
               Point3d pt3d ; GetPoint( cNeigh.GetTopRight().x, cNeigh.GetTopRight().y, pt3d) ;
               vVertices3d.push_back( pt3d) ;
            }
         }
         vbBonusVert[1] = true ;
      }
      // se non l'ho già aggiunto tramite i vicini top aggiungo il punto top left
      else if ( ! bTopLeft) {
         Point3d ptTl( cell.GetBottomLeft().x, cell.GetTopRight().y) ;
         vVertices.push_back( ptTl) ;
         Point3d pt3d ; GetPoint( cell.GetBottomLeft().x, cell.GetTopRight().y, pt3d) ;
         vVertices3d.push_back( pt3d) ;
         vnVert.push_back( int(vVertices.size()) - 1) ;
      }
      else
         vnVert.push_back( int(vVertices.size()) - 1) ;
      vNeigh.clear() ;
      vVertices.push_back( cell.GetBottomLeft()) ;
      pt3d = ORIG ; GetPoint( cell.GetBottomLeft().x, cell.GetBottomLeft().y, pt3d) ;
      vVertices3d.push_back( pt3d) ;
      vnVert.push_back( int(vVertices.size()) - 1) ;

      BOOLVECTOR vbKeepPoint( vVertices.size()) ;
      std::fill( vbKeepPoint.begin(), vbKeepPoint.end(), true) ;

      // se non è trimmata aggiusto subito il vettore dei vertici, sennò i salvo le informazioni per quando creerò i poligoni trimmato
      vVerticesCorr = vVertices ;
      // ora devo controllare se uno dei lati della cella è collassato in un punto.
      // se così, devo guardare se sui due lati adiacenti a quello di polo ho messo dei punti extra
      // se ne ho messi solo su uno dei due allora devo togliere il vertice dell'altro lato
      bool bAlreadyCropped = false ;
      if ( AreSamePointApprox(ptP00, ptP10)) {
         // sennò devo togliere l'estremo del lato su cui NON ho punti bonus
         if ( vbBonusVert[1]  &&  ! vbBonusVert[3])  {// lati contati a partire da quello sopra in senso CCW
            if ( ! m_bTrimmed) {
               vVerticesCorr.erase(vVerticesCorr.begin() + vnVert[1]) ;  // vertici della cella contati a partire da ptBL in senso CCW
               vVertices3d.erase(vVertices3d.begin() + vnVert[1]) ;
            }
            vbKeepPoint[vnVert[1]] = false ;
            cell.m_nVertToErase = 1 ; // ptBr
         }
         else if ( ! vbBonusVert[1]  &&  vbBonusVert[3] ) {
            if ( ! m_bTrimmed) {
               // dovrei eliminare ptBL, quindi devo eliminare il primo e l'ultimo punto della polyline e poi chiuderla con quello che era il penultimo punto
               vVerticesCorr.pop_back() ;
               vVerticesCorr[0] = vVerticesCorr.end()[-1] ;
               vVertices3d.pop_back() ;
               vVertices3d[0] = vVertices3d.end()[-1] ;
            }
            vbKeepPoint[0] = false ;
            vbKeepPoint.back() = false ;
            cell.m_nVertToErase = 0 ; // ptBL
         }
         // se non ho bonus su nessuno dei due
         else {
            // ne scelgo uno dei due da togliere
            vbKeepPoint[vnVert[1]] = false ;
            cell.m_nVertToErase = 1 ; // ptBr
         }
         bAlreadyCropped = true ;
      }
      if ( AreSamePointApprox(ptP10, ptP11)) {
         if ( bAlreadyCropped)
            continue ;
         // sennò devo togliere l'estremo del lato su cui NON ho punti bonus
         if ( vbBonusVert[0]  &&  ! vbBonusVert[2]){  // lati contati a partire da quello sopra in senso CCW
            if ( ! m_bTrimmed) {
               vVerticesCorr.erase(vVerticesCorr.begin() + vnVert[1]) ;  // vertici della cella contati a partire da ptBL in senso CCW
               vVertices3d.erase(vVertices3d.begin() + vnVert[1]) ;
            }
            vbKeepPoint[vnVert[1]] = false ;
            cell.m_nVertToErase = 1 ; // ptBr
         }
         else if ( ! vbBonusVert[0]  &&  vbBonusVert[2] ){
            if ( ! m_bTrimmed) {
               vVerticesCorr.erase(vVerticesCorr.begin() + vnVert[2]) ;
               vVertices3d.erase(vVertices3d.begin() + vnVert[2]) ;
            }
            vbKeepPoint[vnVert[2]] = false ;
            cell.m_nVertToErase = 2 ; // ptTR
         }
         // se non ho bonus su nessuno dei due
         else {
            // ne scelgo uno dei due da togliere
            vbKeepPoint[vnVert[2]] = false ;
            cell.m_nVertToErase = 2 ; // ptTR
         }
         bAlreadyCropped = true ;
      }
      if ( AreSamePointApprox(ptP11, ptP01)) {
         if ( bAlreadyCropped)
            continue ;
         // sennò devo togliere l'estremo del lato su cui NON ho punti bonus
         if ( vbBonusVert[1]  &&  ! vbBonusVert[3]){  // lati contati a partire da quello sopra in senso CCW
            if ( ! m_bTrimmed) {
               vVerticesCorr.erase(vVerticesCorr.begin() + vnVert[2]) ;  // vertici della cella contati a partire da ptBL in senso CCW
               vVertices3d.erase(vVertices3d.begin() + vnVert[2]) ;
            }
            vbKeepPoint[vnVert[2]] = false ;
            cell.m_nVertToErase = 2 ; // ptTR
         }
         else if ( ! vbBonusVert[1]  &&  vbBonusVert[3]) {
            if ( ! m_bTrimmed) {
               vVerticesCorr.erase(vVerticesCorr.begin() + vnVert[3]) ;
               vVertices3d.erase(vVertices3d.begin() + vnVert[3]) ;
            }
            vbKeepPoint[vnVert[3]] = false ;
            cell.m_nVertToErase = 3 ; // ptTl
         }
         // se non ho bonus su nessuno dei due
         else {
            // ne scelgo uno dei due da togliere
            vbKeepPoint[vnVert[3]] = false ;
            cell.m_nVertToErase = 3 ; // ptTl
         }
         bAlreadyCropped = true ;
      }
      if ( AreSamePointApprox(ptP01, ptP00)) {
         if ( bAlreadyCropped)
            continue ;
         // sennò devo togliere l'estremo del lato su cui NON ho punti bonus
         if ( vbBonusVert[0]  &&  ! vbBonusVert[2]) {  // lati contati a partire da quello sopra in senso CCW
            if ( ! m_bTrimmed) {
               // dovrei eliminare ptBL, quindi devo eliminare il primo e l'ultimo punto della polyline e poi chiuderla con quello che era il penultimo punto
               vVerticesCorr.pop_back() ;
               vVerticesCorr[0] = vVerticesCorr.end()[-1] ;
               vVertices3d.pop_back() ;
               vVertices3d[0] = vVertices3d.end()[-1] ;
            }
            vbKeepPoint[0] = false ;
            vbKeepPoint.back() = false ;
            cell.m_nVertToErase = 0 ; // ptBL
         }
         else if ( ! vbBonusVert[0]  &&  vbBonusVert[2]) {
            if ( ! m_bTrimmed) {
               vVerticesCorr.erase(vVerticesCorr.begin() + vnVert[3]) ;
               vVertices3d.erase(vVertices3d.begin() + vnVert[3]) ;
            }
            vbKeepPoint[vnVert[3]] = false ;
            cell.m_nVertToErase = 3 ; // ptTl
         }
         // se non ho bonus su nessuno dei due
         else {
            // ne scelgo uno dei due da togliere
            vbKeepPoint[vnVert[3]] = false ;
            cell.m_nVertToErase = 3 ; // ptTl
         }
         bAlreadyCropped = true ;
      }

      if ( ! m_bTrimmed) {
         // se ho una cella con vicino dello stesso grado ( quindi il poligono ha solo 5 punti) controllo la curvatura nella cella e
         // se necessario cambio l'ordine dei vertici per scegliere la diagonale di split migliore
         if ( vVertices.size() == 5) {
            Point3d ptPSrf ;
            double dU, dV ;
            dU = ( cell.GetBottomLeft().x + cell.GetTopRight().x) / 2 / SBZ_TREG_COEFF ;
            dV = ( cell.GetBottomLeft().y + cell.GetTopRight().y) / 2 / SBZ_TREG_COEFF ;
            GetPoint( dU, dV, ptPSrf) ;
            Point3d ptP00P11 = ( ptP00 + ptP11) / 2 ;
            Point3d ptP10P01 = ( ptP10 + ptP01) / 2 ;
            // ho la curvatura maggiore sulla diagonale tra P10 e P01, ruoto l'ordine dei vertici, in modo che triangulate prenda la diagonale giusta
            if ( Dist( ptP00P11, ptPSrf) + EPS_SMALL > Dist( ptP10P01, ptPSrf)) {
               rotate( vVertices.begin(), vVertices.begin() + 1,vVertices.end()) ;
               vVertices.back() = vVertices.at( 0) ;
               rotate( vVertices3d.begin(), vVertices3d.begin() + 1,vVertices3d.end()) ;
               vVertices3d.back() = vVertices3d.at( 0) ;
               rotate( vbKeepPoint.begin(), vbKeepPoint.begin() + 1,vbKeepPoint.end()) ;
               vbKeepPoint.back() = vbKeepPoint.at( 0) ;
               rotate( vVerticesCorr.begin(), vVerticesCorr.begin() + 1,vVerticesCorr.end()) ;
               vVerticesCorr.back() = vVerticesCorr.at( 0) ;
            }
         }
      }

      cell.AddPoly( c) ;
      // aggiorno il contatore dei poligoni
      c += 1 ;
      vPolygonsBasic.emplace_back() ;
      vPolygonsCorrected.emplace_back() ;
      vPolygons3d.emplace_back() ;
      int nVertCorr = int(vVerticesCorr.size());
      for ( int i = 0 ; i < (int) vVertices.size() ; ++i) {
         int dPar = i ;
         if ( ! vbKeepPoint[i])
            dPar = -1 ;
         vPolygonsBasic.back().AddUPoint( dPar, vVertices.at( i)) ;
         if ( i < nVertCorr){
            vPolygonsCorrected.back().AddUPoint( dPar, vVerticesCorr.at( i)) ;
            vPolygons3d.back().AddUPoint( dPar, vVertices3d.at( i)) ;
         }
      }
   }

   return true ;
}

//----------------------------------------------------------------------------
void
Tree::ResetTree( void)
{  
   // setto tutte le foglie a Processed = false
   for ( int nC : m_vnLeaves) {
      m_mTree[nC].SetProcessed( false) ;
   }
}

//----------------------------------------------------------------------------
INTVECTOR
Tree::FindCell( const Point3d& ptToAssign, const CurveLine& clTrim, bool bRecurs) const
{
   INTVECTOR nCells ;
   int nId = -1 ;
   // se fallisce ritorna un vettore vuoto
   // verifico che il punto sia all'interno dello spazio parametrico
   // allargo i bordi in modo da tenere anche i punti sul bordo dello spazio parametrico
   if ( ptToAssign.x < m_mTree.at( -1).GetBottomLeft().x - EPS_SMALL  ||
        ptToAssign.x > m_mTree.at( -1).GetTopRight().x  + EPS_SMALL  ||
        ptToAssign.y < m_mTree.at( -1).GetBottomLeft().y - EPS_SMALL  ||
        ptToAssign.y > m_mTree.at( -1).GetTopRight().y + EPS_SMALL) {
      return nCells ;
   }
   
   // se ho diviso preliminarmente le patches e in uno dei due parametri ho un numero dispari di patches devo individuare a mano la cella parent
   // in cui individuare la foglia giusta
   if ( m_nSpanU > 1  ||  m_nSpanV > 1) {
      INTVECTOR nParents = FindCell( ptToAssign, clTrim, m_vnParents) ;
      if ( nParents.empty())
         return nCells ;
      nId = nParents.back() ;
   }
   // individuo la foglia in cui ho lo start del loop
   while ( ! m_mTree.at( nId).IsLeaf()) {
      if ( m_mTree.at( nId).IsSplitVert()) {
         double dMid = ( m_mTree.at( nId).GetBottomLeft().x + m_mTree.at( nId).GetTopRight().x) / 2 ;
         if ( ptToAssign.x < dMid) {
            nId = m_mTree.at( nId).m_nChild1 ;
         }
         else {
            nId = m_mTree.at( nId).m_nChild2 ;
         }
      }
      else {
         double dMid = ( m_mTree.at( nId).GetBottomLeft().y + m_mTree.at( nId).GetTopRight().y) / 2 ;
         if ( ptToAssign.y < dMid) {
            nId = m_mTree.at( nId).m_nChild2 ;
         }
         else {
            nId = m_mTree.at( nId).m_nChild1 ;
         }
      }
   }

   nCells.push_back( nId) ;
   if ( ! bRecurs) {
      // se sono in un vertice o su un lato devo controllare di aver trovato la cella giusta
      int nEdge = - 2 ;
      OnWhichEdge( nId, ptToAssign, nEdge) ;
      if ( nEdge != - 2) {
         Vector3d vtDir ;
         clTrim.GetStartDir( vtDir) ;
         // proseguo lungo la curva di trim di EPS_SMALL
         Point3d ptToAssignPlus = ptToAssign + vtDir * EPS_SMALL ;
         // se la curva di trim è praticamente parallela ad un lato allora giro a destra, perché voglio considerare intersecate le celle esterne ( sul bordo) al loop
         if ( abs( vtDir.x) > 1 - EPS_SMALL  ||  abs( vtDir.y) > 1 - EPS_SMALL) {
            Vector3d vtDirDX = vtDir ; vtDirDX.Rotate( Z_AX, 90) ;
            ptToAssignPlus = ptToAssignPlus + vtDir * EPS_SMALL ;
            // controllo di non essere uscito dallo spazio parametrico ed eventualmente giro a sinistra
            if ( ptToAssignPlus.x < m_mTree.at( -1).GetBottomLeft().x - EPS_SMALL  ||
                 ptToAssignPlus.x > m_mTree.at( -1).GetTopRight().x  + EPS_SMALL  ||
                 ptToAssignPlus.y < m_mTree.at( -1).GetBottomLeft().y - EPS_SMALL  ||
                 ptToAssignPlus.y > m_mTree.at( -1).GetTopRight().y + EPS_SMALL) {
               // rispetto al punto di partenza avanzo lungo la curva di trim
               ptToAssignPlus = ptToAssign + vtDir * EPS_SMALL ;
               Vector3d vtDirSX = vtDir ; vtDirSX.Rotate( Z_AX, -90) ;
               // e poi giro a sinistra
               ptToAssignPlus = ptToAssign + vtDir * EPS_SMALL ;
            }
         }
         INTVECTOR nCellsSave = nCells ;
         nCells = FindCell( ptToAssignPlus, clTrim, true) ;
         if ( nCells.empty()) {
            // se nonostante tutto non trovo nulla allora tengo il risultato ottenuto scorrendo l'albero
            nCells = nCellsSave ;
         }
      }
   }
   return nCells ;
}

//----------------------------------------------------------------------------
INTVECTOR
Tree::FindCell( const Point3d& ptToAssign, const CurveLine& cl, INTVECTOR vCells, bool bRecurs) const
{
   // se non trova nulla restituisce un vettore vuoto
   // restituisce sempre solo una cella

   // questo punto va trovato con precisione, perché se seleziono la cella sbagliata si incasina tutto
   INTVECTOR nCells ;
   int nId = -1 ;
   for ( int nCell : vCells) {
      if ( ptToAssign.x > m_mTree.at( nCell).GetBottomLeft().x  &&  ptToAssign.x < m_mTree.at( nCell).GetTopRight().x  &&
         ptToAssign.y > m_mTree.at( nCell).GetBottomLeft().y  &&  ptToAssign.y < m_mTree.at( nCell).GetTopRight().y) {
         nId = nCell ;
         nCells.push_back( nId) ;
      }
   }
   // se non ho trovato nulla vuol dire che sono su un vertice o su un lato
   if ( nCells.empty()  &&  ! bRecurs) {
      int nEdge = -2 ;
      for ( int nCell : vCells) {
         OnWhichEdge( nCell, ptToAssign, nEdge) ;
         if ( nEdge != -2  &&  (m_mTree.at(nCell).m_nCollapsed == Cell::Collapsed::NO_COLLAPSE  ||  m_mTree.at(nCell).m_nCollapsed == Cell::Collapsed::TO_VERIFY))
            nCells.push_back( nCell) ;
         nEdge = -2 ;
      }
      if ( ssize( nCells) == 1)
         return nCells ;

      Vector3d vtDir ;
      cl.GetStartDir( vtDir) ;
      Point3d ptIntersPlus = ptToAssign ;
      // N.B. DESTRA e SINISTRA SONO SEMPRE RIFERITE al guardare nella direzione del vettore di partenza della curva di trim
      // potrei essere su un lato
      if ( nEdge == 0  ||  nEdge == 1  ||  nEdge == 2  ||  nEdge == 3) {
         // se è orientato con il lato mi sposto a destra
         if ( ( nEdge == 0 || nEdge == 2 ) && abs(vtDir.x) > 1 - EPS_SMALL/100) {
            Vector3d vtDirDX = vtDir ; vtDirDX.Rotate( Z_AX, -90) ;
            ptIntersPlus = ptIntersPlus + vtDirDX * EPS_SMALL ;
         }
         else if ( (nEdge == 0  ||  nEdge == 2)  &&  abs(vtDir.y) > EPS_SMALL/100)
            ptIntersPlus.y +=  (vtDir.y > 0 ? 1 : -1) * EPS_SMALL ;
         else if ( ( nEdge == 1 || nEdge == 3 ) && abs(vtDir.y) > 1 - EPS_SMALL/100 ) {
            Vector3d vtDirDX = vtDir ; vtDirDX.Rotate( Z_AX, -90) ;
            ptIntersPlus = ptIntersPlus + vtDirDX * EPS_SMALL ;
         }
         else if ( (nEdge == 1  ||  nEdge == 3)  &&  abs(vtDir.x) > EPS_SMALL/100)
            ptIntersPlus.x += (vtDir.x > 0 ? 1 : -1) * EPS_SMALL ;
      }
      // o in un vertice
      else if ( nEdge == 4  ||  nEdge == 5  ||  nEdge == 6  ||  nEdge == 7) {
         // se non è orientato con gli assi procedo con la direzione del trim
         if ( abs(vtDir.x) < 1 - EPS_SMALL/100  &&  abs(vtDir.y) < 1 - EPS_SMALL/100 )
            ptIntersPlus = ptIntersPlus + vtDir * EPS_SMALL ;
         // altrimenti ruoto a destra
         else if (( nEdge == 4 && vtDir.x > 1 - EPS_SMALL / 100)  ||  ( nEdge == 6 && vtDir.x < - 1 + EPS_SMALL / 100) ||
                  ( nEdge == 5 && vtDir.y > 1 - EPS_SMALL / 100)  ||  ( nEdge == 7 && vtDir.y < - 1 + EPS_SMALL / 100)) {
            Vector3d vtDirDX = vtDir ; vtDirDX.Rotate( Z_AX, -45) ;
            ptIntersPlus = ptIntersPlus + vtDirDX * EPS_SMALL ;
         }
         // altrimenti ruoto a sinistra
         else /*if (( nEdge == 4 && vtDir.y < - 1 + EPS_SMALL / 100)  ||  ( nEdge == 6 && vtDir.y < 1 - EPS_SMALL / 100) ||
                    ( nEdge == 5 && vtDir.x > 1 - EPS_SMALL / 100)  ||  ( nEdge == 7 && vtDir.x < - 1 + EPS_SMALL / 100))    // + tutti gli altri casi */ {
            Vector3d vtDirDX = vtDir ; vtDirDX.Rotate( Z_AX, 45) ;
            ptIntersPlus = ptIntersPlus + vtDirDX * EPS_SMALL ;
         }
      }
      //
      else
         return nCells ;
      // controllo di non essere uscito dallo spazio parametrico (in tal caso riparto da ptToAssign e mi muovo a sinistra)
      if ( ptIntersPlus.x < m_mTree.at( -1).GetBottomLeft().x  ||  ptIntersPlus.x > m_mTree.at( -1).GetTopRight().x  ||
         ptIntersPlus.y < m_mTree.at( -1).GetBottomLeft().y  ||  ptIntersPlus.y > m_mTree.at( -1).GetTopRight().y) {
         Vector3d vtDirSX = vtDir ; vtDirSX.Rotate(Z_AX, 45) ;
         ptIntersPlus = ptToAssign + vtDirSX * EPS_SMALL ;
      }
      // rilancio la rigida ricerca
      nCells = FindCell( ptIntersPlus, cl, nCells, true) ;
   }

   return nCells ;
}

//----------------------------------------------------------------------------
bool
Tree::UpdateSplitLoop( ICurveComposite* pCC, Point3d& pt)
{
   Point3d ptLast ;
   bool bAdded = false ;
   if ( ! pCC->GetEndPoint( ptLast)) {
      // potrei avere una compo vuota o con solo un punto
      // se non riesco ad aggiungere una linea allora era una compo vuota
      if ( ! pCC->GetOnlyPoint( ptLast)) {
         pCC->AddPoint( pt) ;
         ptLast = pt ;
      }
      else {
         pCC->Clear() ;
         pCC->AddPoint( ptLast) ;
         pCC->AddLine( pt) ;
      }
      bAdded = true ;
   }
   if ( pCC->IsValid()) {
      if ( ! bAdded)
         pCC->AddLine( pt) ;
      Vector3d vtDir = pt - ptLast ;
      if ( ! vtDir.Normalize())
         return false ;
      // se sono allineato con gli edge della cella ROOT
      int nEdge = -1 ;
      if ( abs( vtDir.x) > 1 - EPS_SMALL) {
         if ( pt.y < EPS_SMALL)
            nEdge = 2 ;
         else if ( m_nSpanV * SBZ_TREG_COEFF - pt.y < EPS_SMALL)
            nEdge = 0 ;
      }
      else if ( abs( vtDir.y) > 1 - EPS_SMALL) {
         if ( pt.x < EPS_SMALL )
            nEdge = 1 ;
         else if ( m_nSpanU * SBZ_TREG_COEFF - pt.x < EPS_SMALL )
            nEdge = 3 ;
      }
      if ( nEdge != -1) {
         m_vCEdge2D[nEdge].second.AddCurve( m_vCEdge2D[nEdge].first.size() + 1, ptLast, vtDir, pt, vtDir) ;
         m_vCEdge2D[nEdge].first.emplace_back( BIPOINT( ptLast, pt)) ;
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
Tree::TraceLoopLabelCell( const POLYLINEVECTOR& vplPolygons)
{
   // l'idea è di seguire il loop di trim, passando di cella in cella, restando, quando possibile, a destra del loop.
   // ( destra e sinistra sono determinate dal verso del loop).

   // il valore è negativo perché voglio considerare contenuto anche un punto che sta su un lato
   double dLinTol = - EPS_SMALL ;
   // percorro i loop trovando le interezioni con le celle e riempiendo i vettori m_vInters delle varie celle
   for ( int i = 0 ; i < (int) m_vPlApprox.size() ; ++ i) {
      PolyLine plLoop = get<0>( m_vPlApprox[i]) ;
      // creo la compo che aggiunge alla polyline originale degli split dove interseca le celle ( serve per ricostruire gli edge aperti della superficie)
      PtrOwner<ICurveComposite> pCCLoopSplit( CreateCurveComposite()) ;
      // controllo se il loop è CCW o CW
      bool bCCW = get<1>( m_vPlApprox[i]) ;
      // trovo in quale cella è il ptStart
      Point3d ptStart ;
      plLoop.GetFirstPoint( ptStart) ;
      PNTULIST lPt = plLoop.GetUPointList() ;
      PNTULIST:: iterator ptFirst = lPt.begin() ;
      PNTULIST:: iterator ptSecond = ptFirst ;
      advance( ptSecond, 1) ;
      CurveLine clFirst ;
      clFirst.Set( ptFirst->first, ptSecond->first) ;
      // aggiorno la compo splittata
      UpdateSplitLoop( pCCLoopSplit, ptFirst->first) ;
      // individuo la cella da cui parte il loop
      INTVECTOR nCells = FindCell( ptStart, clFirst) ;
      int nId ;
      if ( ! nCells.empty())
         nId = nCells.back() ;
      else
         // il loop è risultato fuori dallo spazio parametrico // può capitare se ho una superficie trimmata con più chunk con bbox separate.
         // ma potrebbe anche essere indicatore di un errore
         continue ; 
      int nFirstCell = nId ;
      // trovo quali punti della polyline sono nella cella e l'intersezione
      PNTVECTOR vptInters ;
      vptInters.push_back( ptStart) ;
      // qui mi devo salvare quanti elementi ho già nel vettore m_vInters della cella
      // per poter fare il merge con l'ultimo vptInters della curva di trim, che sarà ancora in questa cella
      // e terminerà nel punto di start
      Cell* pCell = &m_mTree[nId] ;
      int nPass = (int) pCell->m_vInters.size() ;
      pCell->m_vInters.emplace_back() ;
      // salvo il verso del loop
      pCell->m_vInters.back().bCCW = bCCW ;
      // salvo il chunk del loop
      pCell->m_vInters.back().nChunk = m_mChunk[i] ;
      bool bLoopInside = true ;
      Point3d ptCurr ;
      int nIdPolygon = - 1;
      if ( ! pCell->m_vnPolyId.empty())
         nIdPolygon = pCell->m_vnPolyId[0] ;
      else
         return false ;
      bool bEraseNextPoint = false ;
      while ( plLoop.GetNextPoint( ptCurr)) {
         Point3d ptTStart, ptTEnd ;
         plLoop.GetPrevPoint( ptTStart) ;
         plLoop.GetNextPoint( ptTEnd) ;
         CurveLine clTrim ;
         clTrim.Set( ptTStart, ptTEnd) ;
         Point3d ptLastInters = P_INVALID ;
         int nInfiniteLoopCount = 0 ;
         // qui devo mettere una tolleranza negativa per poter tener conto anche dei punti che sono SULLA curva
         while ( ! IsPointInsidePolyLine( ptCurr, vplPolygons[nIdPolygon], dLinTol)) {
            // sto uscendo dalla cella, quindi cerco l'intersezione
            if ( bEraseNextPoint  &&  ! vptInters.empty()) {
               vptInters.pop_back() ;
               bEraseNextPoint = false ;
            }
            bLoopInside = false ;
            // trovo l'intersezione e passo alla cella successiva. nId viene aggiornato dalla funzione FindInters
            // se non trovo l'intersezione vuol dire che non sono nella cella giusta!
            // al precedente FindInters avrei dovuto passare di cella
            if ( ! FindInters( nId, clTrim, vplPolygons[nIdPolygon], vptInters, true)) {
               // scarterò il punto molto vicino al lato e tengo solo l'intersezione del trim col lato
               if ( ! vptInters.empty())
                  vptInters.pop_back() ;
               plLoop.GetPrevPoint( ptTEnd) ;
               plLoop.GetPrevPoint( ptTStart) ;
               plLoop.GetNextPoint( ptCurr) ;
               clTrim.Set( ptTStart, ptTEnd) ;
               clTrim.ExtendEndByLen( EPS_SMALL * 2) ;
               clTrim.ExtendStartByLen( EPS_SMALL * 2) ;
               if ( ! FindInters( nId, clTrim, vplPolygons[nIdPolygon], vptInters, true))
                  return false ;
               bEraseNextPoint = true ;
            }
            // se l'intersezione e la stessa della precedente allora potrei essere entrato in un loop infinito
            // se per più volte il punto di intersezione resta più o meno lo stesso allora blocco tutto
            if ( AreSamePointEpsilon( vptInters.back(), ptLastInters, 10 * EPS_SMALL)) {
               ++ nInfiniteLoopCount ;
               if ( nInfiniteLoopCount == 4) {
                  LOG_ERROR( GetEGkLogger(), "Error Triangulating SurfBezier: infinte while loop occured in Tree::TraceLoopLabelCell") 
                  return false ;
               }
            }
            ptLastInters = vptInters.back() ;
            // aggiorno il puntatore alla cella
            pCell = &m_mTree[nId] ;
            // recupero l'indice del poligono base associato alla cella
            if ( ! pCell->m_vnPolyId.empty())
               nIdPolygon = pCell->m_vnPolyId[0] ;
            else
               return false ;
            // salvo il verso del loop
            pCell->m_vInters.back().bCCW = bCCW ;
            // salvo il chunk del loop
            pCell->m_vInters.back().nChunk = m_mChunk[i] ;
            // aggiorno la compo splittata
            UpdateSplitLoop( pCCLoopSplit, vptInters.back()) ;
         }
         // controllo di aggiungere un punto abbastanza distante dal precedente
         if ( ! AreSamePointXYEpsilon( vptInters.back(), ptCurr, 10 * EPS_SMALL)) {
            // aggiungo la fine del segmento nel vettore delle intersezioni
            vptInters.push_back( ptCurr) ;
            // aggiorno la compo splittata
            UpdateSplitLoop( pCCLoopSplit, ptCurr) ;
         }
      }
      if ( nId == nFirstCell)
         vptInters.pop_back() ;
      pCell->m_vInters.back().vpt = vptInters ;
      if ( bLoopInside) {
         // setto la categoria della cella
         if ( pCell->m_nFlag == -1)
            pCell->m_nFlag = 2 ;
         else if ( pCell->m_nFlag == 1)
            pCell->m_nFlag = 3 ;
         // setto i lati di ingresso e uscita a -1 per indicare che ho un loop interno alla cella
         pCell->m_vInters.back().nIn = -1 ;
         pCell->m_vInters.back().nOut = -1 ;
      }
      // sono tornato alla cella di partenza, quindi devo fare il merge dei due vettori di intersezione che ho creato per questa cella
      // per lo stesso loop
      else {
         // verifico se sono effettivamente nella cella di partenza o in una cella adiacente
         if ( nId != nFirstCell) {
            Point3d ptFirst =  m_mTree[nFirstCell].m_vInters[nPass].vpt[0] ;
            Point3d ptLast = pCell->m_vInters.back().vpt.back() ;
            //sistemo l'ingresso della prima cella
            int nEdge ;
            if ( ! OnWhichEdge( nFirstCell, ptFirst, nEdge))
               return false ;
            m_mTree[nFirstCell].m_vInters[nPass].nIn = nEdge ;
            // sistemo l'uscita dell'ultima cella
            if ( ! OnWhichEdge( nId, ptLast, nEdge))
               return false ;
            pCell->m_vInters.back().nOut = nEdge ;
            // sistemo il flag dell'ultima cella
            if ( pCell->m_nFlag == -1)
               pCell->m_nFlag = 1 ;
            else if ( pCell->m_nFlag == 2)
               pCell->m_nFlag = 3 ;
         }
         // sono tornato nella cella iniziale, quindi giunto i due vettori intersezione
         else if ( nId == nFirstCell) {
            int nOut = pCell->m_vInters[nPass].nOut ;
            pCell->m_vInters.back().vpt.insert( pCell->m_vInters.back().vpt.end(),
                                                      pCell->m_vInters[nPass].vpt.begin(),
                                                      pCell->m_vInters[nPass].vpt.end()) ;
            pCell->m_vInters[nPass] = pCell->m_vInters.back() ;
            pCell->m_vInters.pop_back() ;
            // sistemo il lato d'uscita
            pCell->m_vInters[nPass].nOut = nOut ;
         }
      }
      // salvo la compo splittata
      m_vCCLoop2D.emplace_back( Release( pCCLoopSplit)) ;
   }

   // riordino i vettori di intersezione per ogni cella e setto il flag RightEdgeIn
   for ( int nId : m_vnLeaves) {
      std::sort( m_mTree[nId].m_vInters.begin(), m_mTree[nId].m_vInters.end()) ;
      SetRightEdgeIn( nId) ;
   }

   // devo riconoscere le celle dentro i loop che sono ancora con label m_bLabelled = false
   ResetTree() ;
   INTVECTOR vNeigh, vFirst ;
   GetRootNeigh( 1, vFirst) ;
   for ( int k : vFirst) {
      m_mTree[k].m_bOnLeftEdge = true ;
   }
   int nLastLeft = vFirst.back() ;
   int nCell = vFirst[0] ;
   GetRightNeigh( nCell, vNeigh) ;
   // proseguo finché non sono sull'elemento più alto di vFirst e tutti i suoi vicini sono processati/categorizzati
   bool bAllDone = false ;
   // mentre scorro a destra mi basta trovare un vicino da cui non sono passato
   // mentre torno indietro a sinistra, mi fermo quando trovo una cella non processata ( con vicini a destra da cui non sono passato)
   while ( ! bAllDone) {
      // categorizzo la cella
      m_mTree[nCell].m_bLabelled = true ;
      CategorizeCell( nCell) ;
      bool bDone = false ;
      // fintanto che la cella ha tra i vicini a destra una cella non elaborata mi sposto a destra
      // definisco una cella Processed se tutto il ramo a destra è categorizzato
      while ( ! m_mTree[nCell].IsProcessed()) {
         // verso la cella a destra più in basso da cui non sono ancora passato
         bool bProceeded = false ;
         for ( int i = 0 ; i < (int)vNeigh.size() ; ++ i) {
            if ( m_mTree[vNeigh[i]].m_bLabelled == false) {
               nCell = vNeigh[i] ;
               bProceeded = true ;
               break ;
            }
         }
         if ( ! bProceeded) {
            m_mTree[nCell].SetProcessed() ;
         }
         else {
            //categorizzo la cella
            m_mTree[nCell].m_bLabelled = true ;
            CategorizeCell( nCell) ;
         }
         
         if ( ! m_mTree[nCell].IsProcessed()) {
            // guardo i vicini a destra per passare alla prossima cella
            vNeigh.clear() ;
            GetRightNeigh( nCell, vNeigh) ;
            bDone = true ;
            //controllo che tra i vicini di destra ce ne sia almeno uno da cui non sono passato
            // ( e controllo che non sia una cella sul lato sinistro ( adiacenza in caso di superficie chiusa))
            for ( int t: vNeigh) {
               if ( m_mTree[t].m_bLabelled == false ) {
                  if ( ! m_bClosedU) {
                     bDone = false ;
                     break ;
                  }
                  else {
                     // controllo che non sia sul lato sinistro
                     if ( ! m_mTree[t].m_bOnLeftEdge) {
                        bDone = false ;
                        break ;
                     }
                  }
               }
            }
            m_mTree[nCell].SetProcessed( bDone) ;
         }
      }
      vNeigh.clear() ;
      GetRightNeigh( nCell, vNeigh) ;
      // se non ho vicini a destra o se i vicini sono già tutti categorizzati
      // torno indietro a sinistra alla cella già categorizzata più bassa
      while ( m_mTree[nCell].IsProcessed()) {
         // trovo il vicino a sinistra, già categorizzato, più basso
         vNeigh.clear() ;
         GetLeftNeigh( nCell, vNeigh) ;
         // però se sono già sul lato sinistro non cerco di spostarmi ancora più a sinistra.
         if ( ! m_mTree[nCell].m_bOnLeftEdge) {
            for ( int p = 0 ; p < (int)vNeigh.size() ; ++ p) {
               if ( m_mTree[vNeigh[p]].m_bLabelled != false) {
                  nCell = vNeigh[p] ;
                  break ;
               }
            }
         }
         if ( ! m_bClosedU) {
            // se non ho vicini a sinistra sono tornato sul lato sinistro
            // se questa cella non è processata devo procedere alla cella più bassa non processata sul lato sinistro
            if ( vNeigh.empty()  &&  m_mTree[nCell].IsProcessed()) {
               for ( int p = 0 ; p < (int)vFirst.size() ; ++ p) {
                  if ( m_mTree[vFirst[p]].m_bLabelled == false) {
                     nCell = vFirst[p] ;
                     break ;
                  }
               }
               if ( nCell == nLastLeft  &&  m_mTree[nCell].IsProcessed()) {
                  bAllDone = true ;
                  break ;
               }
            }
         }
         else {
            // verifico se sono tornato sul lato sinistro, in questo caso procedo con la prima cella di vFirst non processata
            if ( m_mTree[nCell].m_bOnLeftEdge  &&  m_mTree[nCell].IsProcessed()) {
               for ( int p = 0 ; p < (int) vFirst.size(); ++ p) {
                  if ( nCell == vFirst[p]  &&  nCell != nLastLeft) {
                     nCell = vFirst[p + 1] ;
                     break ;
                  }
                  else if ( nCell == nLastLeft  &&  m_mTree[nCell].IsProcessed()){
                     bAllDone = true ;
                  }
               }
            }
            if ( bAllDone)
               break ;
         }
         // controllo se tutti i vicini di destra sono categorizzati
         vNeigh.clear() ;
         GetRightNeigh( nCell, vNeigh) ;
         bool bDone = true ;
         for ( int k = 0 ; k < (int)vNeigh.size() ; ++ k) {
            if ( ! m_mTree[vNeigh[k]].IsProcessed()) {
               bDone = false ;
               break ;
            }
         }
         if ( bDone) {
            m_mTree[nCell].SetProcessed( bDone) ;
            if ( m_mTree[nCell].m_bLabelled == false) {
               m_mTree[nCell].m_bLabelled = true ;
               CategorizeCell( nCell) ;
            }
         }
      }
      vNeigh.clear() ;
      GetRightNeigh( nCell, vNeigh) ;
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
Tree::FindInters( int& nId, const CurveLine& clTrim, const PolyLine& plPolygon, PNTVECTOR& vptInters, bool bFirstInters)
{
   Cell* pCell = &m_mTree.at(nId) ;
   Point3d ptTR = pCell->GetTopRight() ;
   Point3d ptBL = pCell->GetBottomLeft() ;
   Point3d ptTl( ptBL.x , ptTR.y) ;
   Point3d ptBr( ptTR.x , ptBL.y) ;
   int nEdge ; // flag che indica il lato su cui ho l'intersezione a partire dal lato top in senso antiorario
   Point3d ptInters ;

   PtrOwner<ICurveComposite> pCC( CreateCurveComposite()) ;
   PolyLine plSimplePolygon ;
   plSimplePolygon = plPolygon ;
   plSimplePolygon.RemoveAlignedPoints() ;
   pCC->FromPolyLine( plSimplePolygon) ;
   IntersCurveCurve icc( clTrim, *pCC) ;
   if ( icc.GetCrossOrOverlapIntersCount() < 1) {
      CurveLine clTrimExtended = clTrim ;
      clTrimExtended.ExtendEndByLen( 5 * EPS_SMALL) ;
      clTrimExtended.ExtendStartByLen( 5 * EPS_SMALL) ;
      IntersCurveCurve iccRetry( clTrimExtended, *pCC) ;
      icc = iccRetry ;
   }
   if ( icc.GetCrossOrOverlapIntersCount() < 1) {
      CurveLine clTrimExtended = clTrim ;
      clTrimExtended.ExtendEndByLen( 5 * EPS_SMALL) ;
      clTrimExtended.ExtendStartByLen( 5 * EPS_SMALL) ;
      IntersCurveCurve iccRetry( clTrimExtended, *pCC) ;
      icc = iccRetry ;
   }
   IntCrvCrvInfo iccInfo ;
   if ( icc.GetOverlaps()) {
      for ( int i = 0 ; i < icc.GetIntersCount() ; ++i) {
         icc.GetIntCrvCrvInfo( i, iccInfo) ;
         if ( iccInfo.bOverlap){
            ptInters = iccInfo.IciA[1].ptI ;
            break ;
         }
      }
   }
   else {
      int nLastInters = icc.GetIntersCount() -1 ;
      if ( nLastInters < 0)
         return false ;
      icc.GetIntCrvCrvInfo( nLastInters, iccInfo) ;
      ptInters = iccInfo.IciA[0].ptI ;
   }

   // determino il lato/vertice di uscita
   if ( ! OnWhichEdge( nId, ptInters, nEdge))
      return false ;
   pCell->m_vInters.back().nOut = nEdge ;
   // se è risultato che sono abbastanza vicino( di EPS_SMALL) ad un vertice allora il punto aggiunto sarà il vertice vero e proprio
   // se è abbastanza vicino ad un lato allora modifico le sue coordinate in modo che sia esattamente sul lato
   if ( nEdge == 0)
      ptInters.y = ptTl.y ;
   else if ( nEdge == 1)
      ptInters.x = ptTl.x ;
   else if ( nEdge == 2)
      ptInters.y = ptBL.y ;
   else if ( nEdge == 3)
      ptInters.x = ptTR.x ;
   else if ( nEdge == 4)
      ptInters = ptTl ;
   else if ( nEdge == 5)
      ptInters = ptBL ;
   else if ( nEdge == 6)
      ptInters = ptBr ;
   else if ( nEdge == 7)
      ptInters = ptTR ;
   // aggiungo il nuovo punto al vettore delle intersezioni
   if ( vptInters.empty()  ||  ! AreSamePointXYEpsilon( ptInters , vptInters.back(), 10 * EPS_SMALL))
      vptInters.push_back( ptInters) ;
   else {
      // se l'ultimo punto del vettore delle intersezioni è quasi uguale al punto che devo aggiungere allora lo sostituisco con quest'ultimo
      vptInters.pop_back() ;
      vptInters.push_back( ptInters) ;
   }

   // salvo il vettore intersezione per la cella e capisco in quale altra cella passare
   if ( (int)vptInters.size() == 1)
      pCell->m_vInters.back().vpt.push_back( vptInters[0]) ;
   else
      pCell->m_vInters.back().vpt = vptInters ;
   vptInters.clear() ;
   // setto la categoria della cella
   if ( pCell->m_nFlag == -1) 
      pCell->m_nFlag = 1 ;
   else if ( pCell->m_nFlag == 2)
      pCell->m_nFlag = 3 ;

   // seleziono la cella successiva da analizzare
   // 
   // se la superficie è chiusa su un parametro devo stare attento a considerare lo spazio parametrico come se non fosse periodico
   // es: non posso saltare dal lato destro al lato sinistro anche se una cella sul lato sinistro risulta la RightNeigh di una cella sul lato destro!
   INTVECTOR vNeigh, vNeigh1 ;
   if ( nEdge == 0) {
      GetTopNeigh( nId, vNeigh) ;
      std::reverse( vNeigh.begin(), vNeigh.end()) ;
      for ( int j : vNeigh) {
         if ( ptInters.x >= m_mTree[j].GetBottomLeft().x) {
            nId = j ;
            break ;
         }
      }
      // aggiorno il puntatore alla cella
      pCell = &m_mTree.at(nId) ;
      pCell->m_vInters.emplace_back() ;
      Point3d ptBr( pCell->GetTopRight().x, pCell->GetBottomLeft().y) ;
      if ( ptInters.x > pCell->GetBottomLeft().x + EPS_SMALL  &&  ptInters.x < pCell->GetTopRight().x - EPS_SMALL)
         pCell->m_vInters.back().nIn = 2 ;
      else if ( AreSamePointXYApprox( ptInters, pCell->GetBottomLeft()))
         pCell->m_vInters.back().nIn = 5 ;
      else if ( AreSamePointXYApprox( ptInters, ptBr))
         pCell->m_vInters.back().nIn = 6 ;
      else
         pCell->m_vInters.back().nIn = 2 ;
   }
   else if ( nEdge == 1) {
      GetLeftNeigh( nId, vNeigh) ;
      std::reverse( vNeigh.begin(), vNeigh.end()) ;
      for ( int j : vNeigh) {
         if ( ptInters.y >= m_mTree[j].GetBottomLeft().y) {
            nId = j ;
            break ;
         }
      }
      // aggiorno il puntatore alla cella
      pCell = &m_mTree.at(nId) ;
      pCell->m_vInters.emplace_back() ;
      Point3d ptBr( pCell->GetTopRight().x, pCell->GetBottomLeft().y) ;
      if ( ptInters.y > pCell->GetBottomLeft().y + EPS_SMALL  &&  ptInters.y < pCell->GetTopRight().y - EPS_SMALL)
         pCell->m_vInters.back().nIn = 3 ;
      else if ( AreSamePointXYApprox( ptInters, ptBr))
         pCell->m_vInters.back().nIn = 6 ;
      else if ( AreSamePointXYApprox( ptInters, pCell->GetTopRight()))
         pCell->m_vInters.back().nIn = 7 ;
      else
         pCell->m_vInters.back().nIn = 3 ;
   }
   else if ( nEdge == 2) {
      GetBottomNeigh( nId, vNeigh) ;
      for ( int j : vNeigh) {
         if ( ptInters.x <= m_mTree[j].GetTopRight().x) {
            nId = j ;
            break ;
         }
      }
      // aggiorno il puntatore alla cella
      pCell = &m_mTree.at(nId) ;
      pCell->m_vInters.emplace_back() ;
      Point3d ptTl( pCell->GetBottomLeft().x, pCell->GetTopRight().y) ;
      if ( ptInters.x > pCell->GetBottomLeft().x + EPS_SMALL  &&  ptInters.x < pCell->GetTopRight().x - EPS_SMALL)
         pCell->m_vInters.back().nIn = 0 ;
      else if ( AreSamePointXYApprox( ptInters, ptTl))
         pCell->m_vInters.back().nIn = 4 ;
      else if ( AreSamePointXYApprox( ptInters, pCell->GetTopRight()))
         pCell->m_vInters.back().nIn = 7 ;
      else
         pCell->m_vInters.back().nIn = 0 ;
   }
   else if ( nEdge == 3) {
      GetRightNeigh( nId, vNeigh) ;
      for ( int j : vNeigh) {
         if ( ptInters.y <= m_mTree[j].GetTopRight().y) {
            nId = j ;
            break ;
         }
      }
      // aggiorno il puntatore alla cella
      pCell = &m_mTree.at(nId) ;
      pCell->m_vInters.emplace_back() ;
      Point3d ptTl( pCell->GetBottomLeft().x, pCell->GetTopRight().y) ;
      if ( ptInters.y > pCell->GetBottomLeft().y + EPS_SMALL  &&  ptInters.y < pCell->GetTopRight().y - EPS_SMALL)
         pCell->m_vInters.back().nIn = 1 ;
      else if ( AreSamePointXYApprox( ptInters, pCell->GetBottomLeft()))
         pCell->m_vInters.back().nIn = 5 ;
      else if ( AreSamePointXYApprox( ptInters, ptTl))
         pCell->m_vInters.back().nIn = 4 ;
      else
         pCell->m_vInters.back().nIn = 1 ;
   }
   // esco da uno dei vertici
   else if ( nEdge == 4) {
      GetTopNeigh( nId, vNeigh) ;
      GetLeftNeigh( nId, vNeigh1) ;
      INTVECTOR nPossible, nPossible1 ;
      nPossible = FindCell( ptInters, clTrim, vNeigh) ;
      nPossible1 = FindCell( ptInters, clTrim, vNeigh1) ;
      // ingresso dal basso
      if ( ! nPossible.empty()) {
         nId = nPossible[0] ;
         // aggiorno il puntatore alla cella
         pCell = &m_mTree.at(nId) ;
         pCell->m_vInters.emplace_back() ;
         // controllo se entro in un vertice o a metà lato
         Point3d ptBr( pCell->GetTopRight().x, pCell->GetBottomLeft().y) ;
         if ( AreSamePointXYExact( ptInters, ptBr))
            pCell->m_vInters.back().nIn = 6 ;
         else if ( AreSamePointXYExact( ptInters, pCell->GetBottomLeft()))
            pCell->m_vInters.back().nIn = 5 ;
         else
            pCell->m_vInters.back().nIn = 2 ;
      }
      // ingresso da destra
      else if ( ! nPossible1.empty()) {
         nId = nPossible1.back() ;
         // aggiorno il puntatore alla cella
         pCell = &m_mTree.at(nId) ;
         pCell->m_vInters.emplace_back() ;
         // controllo se entro in un vertice o a metà lato
         Point3d ptBr( pCell->GetTopRight().x, pCell->GetBottomLeft().y) ;
         if ( AreSamePointXYExact( ptInters, ptBr))
            pCell->m_vInters.back().nIn = 6 ;
         else if ( AreSamePointXYExact( ptInters, pCell->GetTopRight()))
            pCell->m_vInters.back().nIn = 7 ;
         else
            pCell->m_vInters.back().nIn = 3 ;
      }
      // ingresso in diagonale
      else {
         if ( ! vNeigh.empty()) {
            int nIdTemp = vNeigh[0] ;
            vNeigh.clear() ;
            GetLeftNeigh( nIdTemp, vNeigh) ;
            nId = vNeigh[0] ;
            m_mTree[nId].m_vInters.emplace_back() ;
            m_mTree[nId].m_vInters.back().nIn = 6 ;
         }
         else if (! vNeigh1.empty()) {
            nId = vNeigh1.back() ;
            m_mTree[nId].m_vInters.emplace_back() ;
            m_mTree[nId].m_vInters.back().nIn = 7 ;
         }
         else
            return false ;
      }
   }
   else if ( nEdge == 5) {
      GetLeftNeigh( nId, vNeigh) ;
      GetBottomNeigh( nId, vNeigh1) ;
      INTVECTOR nPossible, nPossible1 ;
      nPossible = FindCell( ptInters, clTrim, vNeigh) ;
      nPossible1 = FindCell( ptInters, clTrim, vNeigh1) ;
      // ingresso dal destra
      if ( ! nPossible.empty()) {
         nId = nPossible[0] ;
         // aggiorno il puntatore alla cella
         pCell = &m_mTree.at(nId) ;
         pCell->m_vInters.emplace_back() ;
         // controllo se entro in un vertice o a metà lato
         Point3d ptBr( pCell->GetTopRight().x, pCell->GetBottomLeft().y) ;
         if ( AreSamePointXYExact( ptInters, ptBr))
            pCell->m_vInters.back().nIn = 6 ;
         else if ( AreSamePointXYExact( ptInters, pCell->GetTopRight()))
            pCell->m_vInters.back().nIn = 7 ;
         else
            pCell->m_vInters.back().nIn = 3 ;
      }
      // ingresso dall'alto
      else if ( ! nPossible1.empty()) {
         nId = nPossible1[0] ;
         // aggiorno il puntatore alla cella
         pCell = &m_mTree.at(nId) ;
         pCell->m_vInters.emplace_back() ;
         // controllo se entro in un vertice o a metà lato
         Point3d ptTl( pCell->GetBottomLeft().x, pCell->GetTopRight().y) ;
         if ( AreSamePointXYExact( ptInters, ptTl))
            pCell->m_vInters.back().nIn = 4 ;
         else if ( AreSamePointXYExact( ptInters, pCell->GetTopRight()))
            pCell->m_vInters.back().nIn = 7 ;
         else
            pCell->m_vInters.back().nIn = 0 ;
      }
      // ingresso in diagonale
      else {
         if ( ! vNeigh.empty()) {
            int nIdTemp = vNeigh[0] ;
            vNeigh.clear() ;
            GetBottomNeigh( nIdTemp, vNeigh) ;
            nId = vNeigh.back() ;
            m_mTree[nId].m_vInters.emplace_back() ;
            m_mTree[nId].m_vInters.back().nIn = 7 ;
         }
         else if (! vNeigh1.empty()) {
            nId = vNeigh1[0] ;
            m_mTree[nId].m_vInters.emplace_back() ;
            m_mTree[nId].m_vInters.back().nIn = 4 ;
         }
         else
            return false ;
      }
   }
   else if ( nEdge == 6) {
      GetBottomNeigh( nId, vNeigh) ;
      GetRightNeigh( nId, vNeigh1) ;
      INTVECTOR nPossible, nPossible1 ;
      nPossible = FindCell( ptInters, clTrim, vNeigh) ;
      nPossible1 = FindCell( ptInters, clTrim, vNeigh1) ;
      // ingresso dall'alto
      if ( ! nPossible.empty()) {
         nId = nPossible.back() ;
         // aggiorno il puntatore alla cella
         pCell = &m_mTree.at(nId) ;
         pCell->m_vInters.emplace_back() ;
         // controllo se entro in un vertice o a metà lato
         Point3d ptTl( pCell->GetBottomLeft().x, pCell->GetTopRight().y) ;
         if ( AreSamePointXYExact( ptInters, ptTl))
            pCell->m_vInters.back().nIn = 4 ;
         else if ( AreSamePointXYExact( ptInters, pCell->GetTopRight()))
            pCell->m_vInters.back().nIn = 7 ;
         else
            pCell->m_vInters.back().nIn = 0 ;
      }
      // ingresso da sinistra
      else if ( ! nPossible1.empty()) {
         nId = nPossible1[0] ;
         // aggiorno il puntatore alla cella
         pCell = &m_mTree.at(nId) ;
         pCell->m_vInters.emplace_back() ;
         // controllo se entro in un vertice o a metà lato
         Point3d ptTl( pCell->GetBottomLeft().x, pCell->GetTopRight().y) ;
         if ( AreSamePointXYExact( ptInters, ptTl))
            pCell->m_vInters.back().nIn = 4 ;
         else if ( AreSamePointXYExact( ptInters, pCell->GetBottomLeft()))
            pCell->m_vInters.back().nIn = 5 ;
         else
            pCell->m_vInters.back().nIn = 1 ;
      }
      // ingresso in diagonale
      else {
         if ( ! vNeigh.empty()){
            int nIdTemp = vNeigh.back() ;
            vNeigh.clear() ;
            GetRightNeigh( nIdTemp, vNeigh) ;
            nId = vNeigh.back() ;
            m_mTree[nId].m_vInters.emplace_back() ;
            m_mTree[nId].m_vInters.back().nIn = 4 ;
         }
         else if (! vNeigh1.empty()) {
            nId = vNeigh1[0] ;
            m_mTree[nId].m_vInters.emplace_back() ;
            m_mTree[nId].m_vInters.back().nIn = 5 ;
         }
         else
            return false ;
      }
   }
   else if ( nEdge == 7) {
      GetRightNeigh( nId, vNeigh) ;
      GetTopNeigh( nId, vNeigh1) ;
      INTVECTOR nPossible, nPossible1 ;
      nPossible = FindCell( ptInters, clTrim, vNeigh) ;
      nPossible1 = FindCell( ptInters, clTrim, vNeigh1) ;
      // ingresso da sinistra
      if ( ! nPossible.empty()) {
         nId = nPossible.back() ;
         // aggiorno il puntatore alla cella
         pCell = &m_mTree.at(nId) ;
         pCell->m_vInters.emplace_back() ;
         // controllo se entro in un vertice o a metà lato
         Point3d ptTl( pCell->GetBottomLeft().x, pCell->GetTopRight().y) ;
         if ( AreSamePointXYExact( ptInters, ptTl))
            pCell->m_vInters.back().nIn = 4 ;
         else if ( AreSamePointXYExact( ptInters, pCell->GetBottomLeft()))
            pCell->m_vInters.back().nIn = 5 ;
         else
            pCell->m_vInters.back().nIn = 1 ;
      }
      // ingresso dal basso
      else if ( ! nPossible1.empty()) {
         nId = nPossible1.back() ;
         // aggiorno il puntatore alla cella
         pCell = &m_mTree.at(nId) ;
         pCell->m_vInters.emplace_back() ;
         // controllo se entro in un vertice o a metà lato
         Point3d ptBr( pCell->GetTopRight().x, pCell->GetBottomLeft().y) ;
         if ( AreSamePointXYExact( ptInters, ptBr))
            pCell->m_vInters.back().nIn = 6 ;
         else if ( AreSamePointXYExact( ptInters, pCell->GetBottomLeft()))
            pCell->m_vInters.back().nIn = 5 ;
         else
            pCell->m_vInters.back().nIn = 2 ;
      }
      // ingresso in diagonale
      else {
         if ( ! vNeigh.empty()  /*&&  ! m_mTree[vNeigh[0]].m_bOnLeftEdge */) {  // questa aggiunta in teoria non serve
            int nIdTemp = vNeigh.back() ;
            vNeigh.clear() ;
            GetTopNeigh( nIdTemp, vNeigh) ;
            nId = vNeigh[0] ;
            m_mTree[nId].m_vInters.emplace_back() ;
            m_mTree[nId].m_vInters.back().nIn = 5 ;
         }
         else if (! vNeigh1.empty()) {
            nId = vNeigh1.back() ;
            m_mTree[nId].m_vInters.emplace_back() ;
            m_mTree[nId].m_vInters.back().nIn = 6 ;
         }
         else
            return false ;
      }
   }
   // aggiungo l'intersezione al vettore delle intersezioni della prossima cella
   vptInters.push_back( ptInters) ;
   return true ;
}

//----------------------------------------------------------------------------
bool
Tree::CreateCellPolygons( int nId, POLYLINEMATRIX& vPolygons, POLYLINEMATRIX& vPolygons3d, INTVECTOR& vToCheck, int& nPoly,
                          INTVECTOR& vnParentChunk, const PolyLine& plCell, const PolyLine& plCell3d)
{
  // conto quanti vertici in più ho per lato e creo un vettore dei vertici per lato
   Cell& cCell = m_mTree[nId] ;
   PNTMATRIX vEdgeVertex(4) ;
   PNTMATRIX vEdgeVertex3d(4) ;
   Point3d ptBL = cCell.GetBottomLeft() ;
   Point3d ptTR = cCell.GetTopRight() ;
   Point3d ptTl = cCell.GetTopLeft() ;
   Point3d ptBr = cCell.GetBottomRight() ;
   int nVertToErase = cCell.m_nVertToErase ;
   if ( nVertToErase != 2)
      vEdgeVertex[0].push_back( ptTR) ;
   if ( nVertToErase != 3)
      vEdgeVertex[1].push_back( ptTl) ;
   if ( nVertToErase != 0)
      vEdgeVertex[2].push_back( ptBL) ;
   if ( nVertToErase != 1)
      vEdgeVertex[3].push_back( ptBr) ;
   // capisco se sono in modalità ForTriangulation
   bool bForTriangulation = plCell3d.GetPointNbr() != 0 ;
   if ( bForTriangulation) {
      if ( nVertToErase != 2) {
         Point3d pt3d ; GetPoint( cCell.GetTopRight().x, cCell.GetTopRight().y, pt3d) ;
         vEdgeVertex3d[0].push_back( pt3d) ;
      }
      if ( nVertToErase != 3) {
         Point3d pt3d ; GetPoint( cCell.GetBottomLeft().x, cCell.GetTopRight().y, pt3d) ;
         vEdgeVertex3d[1].push_back( pt3d) ;
      }
      if ( nVertToErase != 0) {
         Point3d pt3d ; GetPoint( cCell.GetBottomLeft().x, cCell.GetBottomLeft().y, pt3d) ;
         vEdgeVertex3d[2].push_back( pt3d) ;
      }
      if ( nVertToErase != 1) {
         Point3d pt3d ; GetPoint( cCell.GetTopRight().x, cCell.GetBottomLeft().y, pt3d) ;
         vEdgeVertex3d[3].push_back( pt3d) ;
      }
   }

  // la PolyLine è riempita a partire dal lato bottom
   Point3d ptStart, pt3d ;
   plCell.GetFirstPoint( ptStart) ;
   if ( bForTriangulation)
      plCell3d.GetFirstPoint( pt3d) ;
   INTVECTOR vEdge = { 2, 3, 0, 1} ;
   for ( int p = 0 ; p < 4 ; ++ p) {
      int j = vEdge[p] ;
      int next = j +1 ;
      if ( j == 3)
         next = 0 ;
      Point3d ptToAdd ;
      Point3d ptNextVert ;
      switch ( next ) {
      case 0 : ptNextVert = ptTR ; break ;
      case 1 : ptNextVert = ptTl ; break ;
      case 2 : ptNextVert = ptBL ; break ;
      case 3 : ptNextVert = ptBr ; break ;
      }
      if ( ! bForTriangulation) {
         while ( plCell.GetNextPoint( ptToAdd)  &&  ! AreSamePointXYExact( ptToAdd, ptNextVert))
            vEdgeVertex[j].push_back( ptToAdd) ;
      }
      else {
         double dPar = 0 ;
         while ( plCell.GetNextUPoint( &dPar, &ptToAdd)  &&  ! AreSamePointXYExact( ptToAdd, ptNextVert)) {
            vEdgeVertex[j].push_back( ptToAdd) ;
            if ( dPar > 0)
               plCell3d.GetNextPoint( pt3d) ;
            vEdgeVertex3d[j].push_back( pt3d) ;
         }
         if ( dPar > 0)
            plCell3d.GetNextPoint( pt3d) ;
      }
   }
   
   // comincio a costruire il poligono
   INTVECTOR vToCheckNow = vToCheck ;
   // costruisco i poligoni partendo dal vettore delle intersezioni, come spiegato a pag15 di Cripps
   PolyLine plTrimmedPoly ;
   PolyLine plTrimmedPoly3d ;
   // numero di volte che la cella è stata attraversata da una curva di trim
   int nPassToCheck = (int) vToCheckNow.size() ;
   // numero di vertici aggiunti al nuovo poligono
   int c = 0 ;
   // scorro sui vettori intersezione della cella nId e sui suoi vertici
   // in questo for analizzo solo i loop che tagliano la cella 
   int nEdgeIn = -1 ;
   int nFirstLoopInPoly = -1 ;
   INTVECTOR vAddedLoops ;
   Point3d ptLastAdded( -1, -1, 0) ;
   for ( int w = 0 ; w < (int)vToCheckNow.size() ; ++ w) {
      if ( cCell.m_vInters[w].vpt.size() < 2) {
         continue ;
      }
      // indice del loop in m_vInters
      int j = vToCheckNow[w] ;
      Inters inA = cCell.m_vInters[j] ;
      if ( inA.nIn != -1) {
         int nEdge ;
         if ( nEdgeIn == -1) {
            // salvo il lato di ingresso del primo lato del poligono che sto costruendo
            nEdgeIn = inA.nIn ;
            nFirstLoopInPoly = j ;
         }
         for ( Point3d& ptInt : inA.vpt) {
            AddVertex( nId, vEdgeVertex, vEdgeVertex3d, plTrimmedPoly, c, ptInt, plTrimmedPoly3d, bForTriangulation, ptLastAdded) ;
         }
         vAddedLoops.push_back( j) ;
         nEdge = inA.nOut ;
         // devo verificare di non essere uscito in un vertice, con un tratto sovrapposto al lato che ripercorrerò tra poco
         Point3d ptLast, ptSecondToLast ;
         plTrimmedPoly.GetLastPoint( ptLast) ;
         plTrimmedPoly.GetPrevPoint( ptSecondToLast) ;
         Vector3d vLast = ptLast - ptSecondToLast ;
         vLast.Normalize() ;
         Vector3d vEdge ;
         // estendo: se l'ultimo tratto è sovrapposto e controverso allora elimino l'ultimo punto
         bool bNotEquiverseOverlap = false ;
         int nEdgeLast = -1 ;
         OnWhichEdge(nId,ptLast, nEdgeLast) ;
         if ( (nEdge == 0 || nEdge == 7)  &&  nEdge == nEdgeLast) {
            vEdge.Set( 1,0,0) ;
            if ( AreOppositeVectorApprox( vLast, vEdge)) {
                  plTrimmedPoly.EraseLastUPoint() ;
                  plTrimmedPoly3d.EraseLastUPoint() ;
                  nEdge = 0 ;
                  bNotEquiverseOverlap = true ;
            }
         }
         else if ( ( nEdge == 1 || nEdge == 4 )  &&  nEdge == nEdgeLast) {
            vEdge.Set( 0,-1,0) ;
            if ( AreOppositeVectorApprox( vLast, vEdge)) {
                  plTrimmedPoly.EraseLastUPoint() ;
                  plTrimmedPoly3d.EraseLastUPoint() ;
                  nEdge = 1 ;
                  bNotEquiverseOverlap = true ;
            }
         }
         else if ( ( nEdge == 2 || nEdge == 5)  &&  nEdge == nEdgeLast) {
            vEdge.Set( -1,0,0) ;
            if ( AreOppositeVectorApprox( vLast, vEdge)) {
                  plTrimmedPoly.EraseLastUPoint() ;
                  plTrimmedPoly3d.EraseLastUPoint() ;
                  nEdge = 2 ;
                  bNotEquiverseOverlap = true ;
            }
         }
         else if ( ( nEdge == 3 || nEdge == 6)  &&  nEdge == nEdgeLast) {
            vEdge.Set( 0,1,0) ;
            if ( AreOppositeVectorApprox( vLast, vEdge)) {
                  plTrimmedPoly.EraseLastUPoint() ;
                  plTrimmedPoly3d.EraseLastUPoint() ;
                  nEdge = 3 ;
                  bNotEquiverseOverlap = true ;
            }
         }

         // se mi è rimasto solo un punto sulla polyline vuol dire che avevo solo un tratto parallelo ad lato
         // quindi salto al prossimo loop
         if ( plTrimmedPoly.GetPointNbr() == 1  &&  bNotEquiverseOverlap) {
            plTrimmedPoly.Clear() ;
            plTrimmedPoly3d.Clear() ;
            if ( j == nFirstLoopInPoly)
               nEdgeIn = -1 ;
            continue ;
         }
         if ( nEdge > 3 && nEdge != 7) {
            nEdge = nEdge - 3 ;
         }
         else if ( nEdge == 7) {
            nEdge = 0 ;
         }
         // se ho altri Pass vado avanti ad aggiungere vertici finché trovo il prossimo o finché non sono tornato sul lato di partenza
         bool bNotCameBack = true ;
         bool bValidNextStart = false ;
         bool bAtNextStart = false ;
         int nEdgeWithVertexSkipped ;
         switch ( nVertToErase) {
         case 0 : nEdgeWithVertexSkipped = 2 ; break ;
         case 1 : nEdgeWithVertexSkipped = 3 ; break ;
         case 2 : nEdgeWithVertexSkipped = 0 ; break ;
         case 3 : nEdgeWithVertexSkipped = 1 ; break ;
         default : nEdgeWithVertexSkipped = -1 ; break ;
         }
         if ( w < nPassToCheck - 1) {
            int nSecondCheck = 0 ;
            int nNext ;
            // ciclo sui loop successivi per vedere se ne ho uno con un valid start
            for ( int t = w + 1 ; t < nPassToCheck ; ++ t) {
               bValidNextStart = CheckIfBetween( cCell.m_vInters[j], cCell.m_vInters[vToCheckNow[t]]) ;
               if ( bValidNextStart) {
                  nNext = t ;
                  bAtNextStart = AreSameEdge( cCell.m_vInters[j].nOut, cCell.m_vInters[vToCheckNow[t]].nIn) ;
                  break ;
               }
               else {
                  ++ nSecondCheck ;
               }
            }

            bNotCameBack =  ! ( AreSameEdge( nEdge, nEdgeIn)  &&  CheckIfBefore( plTrimmedPoly, nEdge)) ;

            while ( ! ( bValidNextStart  &&  bAtNextStart)  &&  bNotCameBack) {
               Point3d ptVert ;
               if ( nEdge == 0)
                  ptVert = ptTl ;
               else if ( nEdge == 1)
                  ptVert = ptBL ;
               else if ( nEdge == 2)
                  ptVert = ptBr ;
               else if ( nEdge == 3)
                  ptVert = ptTR ;
               AddVertex( nId, vEdgeVertex, vEdgeVertex3d, plTrimmedPoly, c, ptVert, plTrimmedPoly3d, bForTriangulation, ptLastAdded) ;
               if ( nEdge > 3  &&  nEdge != 7)
                  nEdge = nEdge - 4 ;
               else if ( nEdge < 3)
                  ++ nEdge ;
               else
                  nEdge = 0 ;

               // aggiorno le condizioni per il while
               if ( bValidNextStart)
                  bAtNextStart = AreSameEdge( nEdge, cCell.m_vInters[vToCheckNow[nNext]].nIn) ;

               bNotCameBack =  ! ( AreSameEdge( nEdge, nEdgeIn)  &&  CheckIfBefore( plTrimmedPoly, nEdge)) ;
            }

            // se ho trovato un altro loop salto all'inizio del for, dopo aver aggiunto eventuali punti intermedi
            if ( bValidNextStart) {
               Point3d ptLast ; plTrimmedPoly.GetLastPoint( ptLast) ;
               for ( int p = nEdge == nEdgeWithVertexSkipped ? 0 : 1 ; p < (int) vEdgeVertex[nEdge].size() ; ++ p) {
                  if ( CheckIfBefore( nEdge, vEdgeVertex[nEdge][p], cCell.m_vInters[vToCheckNow[nNext]].vpt[0])  &&
                       CheckIfBefore( nEdge, ptLast, vEdgeVertex[nEdge][p])) {
                     plTrimmedPoly.AddUPoint( c, vEdgeVertex[nEdge][p]) ;
                     plTrimmedPoly3d.AddUPoint( c, vEdgeVertex3d[nEdge][p]) ;
                     ++ c ;
                  }
               }
               w = w + nSecondCheck ;
               continue ;
            }
            // sono tornato indietro
            else if ( ! bNotCameBack){
               Point3d ptStart ;
               plTrimmedPoly.GetFirstPoint( ptStart) ;
               Point3d ptLast ; plTrimmedPoly.GetLastPoint( ptLast) ;
               for ( int p = nEdge == nEdgeWithVertexSkipped ? 0 : 1 ; p < (int) vEdgeVertex[nEdge].size() ; ++ p) {
                  if ( CheckIfBefore( nEdge, vEdgeVertex[nEdge][p], ptStart)  &&  CheckIfBefore( nEdge, ptLast, vEdgeVertex[nEdge][p])) {
                     plTrimmedPoly.AddUPoint( c, vEdgeVertex[nEdge][p]) ;
                     plTrimmedPoly3d.AddUPoint( c, vEdgeVertex3d[nEdge][p]) ;
                     ++ c ;
                  }
               }
            }
         }
         // non ho altri loop quindi aggiungo vertici finché torno al punto di partenza
         else {
            bNotCameBack =  ! ( AreSameEdge( nEdge, nEdgeIn)  &&  CheckIfBefore( plTrimmedPoly, nEdge)) ;

            Point3d ptStart ;
            plTrimmedPoly.GetFirstPoint( ptStart) ;
            while ( bNotCameBack) {
               Point3d ptVert ;
               if ( nEdge == 0)
                  ptVert = ptTl ;
               else if ( nEdge == 1)
                  ptVert = ptBL ;
               else if ( nEdge == 2)
                  ptVert = ptBr ;
               else if ( nEdge == 3)
                  ptVert = ptTR ;
               AddVertex( nId, vEdgeVertex, vEdgeVertex3d, plTrimmedPoly, c, ptVert, plTrimmedPoly3d, bForTriangulation, ptLastAdded) ;
               if ( nEdge > 3  &&  nEdge != 7)
                  nEdge = nEdge - 4 ;
               else if ( nEdge < 3)
                  ++ nEdge ;
               else
                  nEdge = 0 ;
               // aggiorno le condizioni per il while
               bNotCameBack = ! ( AreSameEdge( nEdge, nEdgeIn)  &&  CheckIfBefore( plTrimmedPoly, nEdge)) ;
            }
            Point3d ptLast ; plTrimmedPoly.GetLastPoint( ptLast) ;
            for ( int p = nEdge == nEdgeWithVertexSkipped ? 0 : 1 ; p < (int) vEdgeVertex[nEdge].size() ; ++ p) {
               if ( CheckIfBefore( nEdge, vEdgeVertex[nEdge][p], ptStart)  &&  CheckIfBefore( nEdge, ptLast, vEdgeVertex[nEdge][p])) {
                  plTrimmedPoly.AddUPoint( c, vEdgeVertex[nEdge][p]) ;
                  plTrimmedPoly3d.AddUPoint( c, vEdgeVertex3d[nEdge][p]) ;
                  ++ c ;
               }
            }

         }
         plTrimmedPoly.Close() ;
         plTrimmedPoly3d.Close() ;
         ptLast.Set( -1,-1,0) ;
         // controllo sull'area del poligono, se è 0 ( quindi un segmento), non lo aggiungo
         double dArea ;
         plTrimmedPoly.GetAreaXY( dArea) ;
         //if ( dArea > SQ_EPS_SMALL) {
         if ( dArea > 0) {
            vPolygons.emplace_back() ;
            vPolygons.back().push_back(std::move(plTrimmedPoly)) ;
            if ( bForTriangulation) {
               vPolygons3d.emplace_back() ;
               vPolygons3d.back().push_back(std::move(plTrimmedPoly3d)) ;
            }
            ++ nPoly ;
            vnParentChunk.push_back( inA.nChunk) ;
         }

         c = 0 ;
         plTrimmedPoly.Clear() ;
         plTrimmedPoly3d.Clear() ;
         nEdgeIn = -1 ;
         // devo verificare se tra i loop che sono finiti in vToCheck in realtà qualcuno l'ho usato per fare un poligono
         for ( int k = 0 ; k < (int)vToCheck.size() ; ++ k) {
            for ( int i = 0 ; i < (int)vAddedLoops.size() ; ++ i) {
               if ( vToCheck[k] == vAddedLoops[i]) {
                  vToCheck.erase( vToCheck.begin() + k) ;
               }
            }
         }
      }
      else
         continue ;
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
Tree::CreateIslandAndHoles( int nId, POLYLINEMATRIX& vPolygons, POLYLINEMATRIX& vPolygons3d, int& nPoly, INTVECTOR& vnParentChunk,
                            const PolyLine& plPolygonsBasic, const PolyLine& plPolygonsBasic3d)
{
   // costruisco i poligoni partendo dal vettore delle intersezioni
   Cell& cCell = m_mTree[nId] ;
   //capisco se sono in modalità ForTriangulation

   // loop interni in una cella intersecata
   int nChunkBiggestCW = -1 ;
   if ( cCell.m_nFlag == 3  ||  cCell.m_nFlag == 2) {
      PolyLine plInLoop ;
      PolyLine plInLoop3d ;
      Inters inA ;
      // se ho almeno un loop CW che non è contenuto in un altro poligono o in un loop interno CCW devo aggiungere il bordo
      bool bAllContained = true ;
      bool bContained = false ;
      int nInters = (int) cCell.m_vInters.size() ;
      for ( int n = 0 ; n < nInters ; ++ n) {
         inA = cCell.m_vInters[n] ;
         if ( inA.nIn == -1) {
            // per ogni loop CW verifico che ci sia un loop CCW dello stesso chunk ( che quindi lo contiene)
            if ( ! inA.bCCW) {
               if ( nChunkBiggestCW == -1)
                  nChunkBiggestCW = inA.nChunk ;
               bContained = false ;
               Inters inB = cCell.m_vInters[0] ;
               for ( int c = 0 ; c < nInters ; ++ c){
                  inB = cCell.m_vInters[c] ;
                  if ( inB.nIn == -1) {
                     if ( inB != inA  &&  inB.nChunk == inA.nChunk  &&  inB.bCCW) {
                        bContained = true ;
                        break ;
                     }
                  }
                  else
                     break ;
               }
               bAllContained = bAllContained  &&  bContained ;
            }
         }
         else
            break ;
      }

      if ( cCell.m_nFlag == 2  &&  ! bAllContained) {
         // i loop esterni sono CW, quindi prima dei loop di trim aggiungo il bordo cella
         Point3d ptVert = cCell.GetTopRight() ;
         plInLoop.AddUPoint( 0, ptVert) ;
         ptVert = cCell.GetTopLeft() ;
         plInLoop.AddUPoint( 1, ptVert) ;
         ptVert = cCell.GetBottomLeft() ;
         plInLoop.AddUPoint( 2, ptVert) ;
         ptVert = cCell.GetBottomRight() ;
         plInLoop.AddUPoint( 3, ptVert) ;
         plInLoop.Close() ;
         vPolygons.emplace_back() ;
         vPolygons.back().push_back( std::move(plInLoop)) ;
         ++ nPoly ;
         Point3d ptP00, ptP10, ptP11, ptP01 ;
         GetPoint( cCell.GetBottomLeft().x, cCell.GetBottomLeft().y, ptP00) ;
         GetPoint( cCell.GetTopRight().x, cCell.GetBottomLeft().y, ptP10) ;
         GetPoint( cCell.GetTopRight().x, cCell.GetTopRight().y, ptP11) ;
         GetPoint( cCell.GetBottomLeft().x, cCell.GetTopRight().y, ptP01) ;
         plInLoop3d.AddUPoint( 0, ptP11) ;
         plInLoop3d.AddUPoint( 1, ptP01) ;
         plInLoop3d.AddUPoint( 2, ptP00) ;
         plInLoop3d.AddUPoint( 3, ptP10) ;
         plInLoop3d.Close() ;
         vPolygons3d.emplace_back() ;
         vPolygons3d.back().push_back( std::move(plInLoop3d)) ;
         plInLoop3d.Clear() ;
         // imposto il chunk del loop CW più grande ( il primo che ho incontrato)
         vnParentChunk.push_back( nChunkBiggestCW) ;
         plInLoop.Clear() ;
      }
      for ( int nLoop = 0 ; nLoop < nInters ; ++ nLoop) {
         inA = cCell.m_vInters[nLoop] ;
         if ( inA.nIn == -1) {
            // numero di vertici aggiunti al nuovo poligono
            int k = 0 ;
            for ( Point3d& ptInt : inA.vpt) {
               plInLoop.AddUPoint( k, ptInt) ;
               Point3d pt3d ; GetPoint( ptInt.x, ptInt.y, pt3d) ;
               plInLoop3d.AddUPoint( k, pt3d) ;
               ++ k ;
            }
            plInLoop.Close() ;
            plInLoop3d.Close() ;
            bool bAdded = false ;
            // se il loop è CW devo controllare in quale altro dei poligoni che ho già aggiunto è contenuto
            if ( ! inA.bCCW) {
               Point3d ptStart ;
               plInLoop.GetFirstPoint( ptStart) ;
               int nOtherPoly = (int)vPolygons.size() ;
               for ( int r = 0 ; r < nPoly ; ++r) {
                  if ( IsPointInsidePolyLine( ptStart, vPolygons[nOtherPoly - r - 1][0], -0.01)  &&  vnParentChunk[nPoly - r - 1] == inA.nChunk) {
                     vPolygons[nOtherPoly - r - 1].push_back( plInLoop) ;
                     plInLoop.Clear() ;
                     vPolygons3d[nOtherPoly - r - 1].push_back( std::move(plInLoop3d)) ;
                     plInLoop3d.Clear() ;
                     bAdded = true ;
                     break ;
                  }
               }
            }
            if ( ! bAdded) {
               vPolygons.emplace_back() ;
               vPolygons.back().push_back( std::move(plInLoop)) ;
               plInLoop.Clear() ;
               vPolygons3d.emplace_back( ) ;
               vPolygons3d.back().push_back( std::move(plInLoop3d)) ;
               plInLoop3d.Clear() ;
               ++ nPoly ;
               vnParentChunk.push_back( inA.nChunk) ;
            }
            plInLoop.Clear() ;
            plInLoop3d.Clear() ;
         }
         else
            continue ;
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
Tree::CheckIfBefore( const PolyLine& pl, int nEdge) const
{ 
   // controllo se ptEnd è prima di ptStart sul lato nEdge rispetto al senso antiorario ( quindi se è dopo in senso orario)
   Point3d ptStart, ptEnd ;
   pl.GetFirstPoint( ptStart) ;
   pl.GetLastPoint( ptEnd) ;
   if ( AreSameEdge( nEdge, 0)) {
      return ptEnd.x > ptStart.x ;
   }
   else if ( AreSameEdge( nEdge, 1)) {
      return ptEnd.y > ptStart.y ;
   }
   else if ( AreSameEdge( nEdge, 2)) {
      return ptEnd.x < ptStart.x ;
   }
   else if ( AreSameEdge( nEdge, 3)) {
      return ptEnd.y < ptStart.y ;
   }
   return false ;
}

//----------------------------------------------------------------------------
bool
Tree::CheckIfBefore( const Inters& inA) const
{ 
   // questa funzione è pensata in riferimento al lato 3, quindi nessuno dei due punti può stare su Edge = 3 ( funzione usata per definire RightEdgeIn)
   // controllo se l'ingresso è prima dell'uscita
   int nEdge1 = inA.nIn ;
   int nEdge2 = inA.nOut ;
   if ( nEdge1 == -1)
      return false ;  // questo return è usato in modo improprio perché non fa capire che l'intersezione analizzata non era corretta
   PolyLine pl ;
   pl.AddUPoint( 0, inA.vpt.back()) ;
   pl.AddUPoint( 1, inA.vpt[0]) ;
   INTVECTOR vEdges = { 7, 0, 4, 1, 5, 2, 6} ;
   // controllo se nEdge1 viene prima di  nEdge2. la partenza è da ptTR e l'arrivo è ptBr

   auto iter1 = find( vEdges.begin(), vEdges.end(), nEdge1) ;
   int nPos1 = distance( vEdges.begin(), iter1) ;
   auto iter2 = find( vEdges.begin(), vEdges.end(), nEdge2) ;
   int nPos2 = distance( vEdges.begin(), iter2) ;
   if ( nPos1 < nPos2)
      return true ;
   else if ( nPos1 > nPos2)
      return false ;
   // nPos1 == nPos2
   else {
      if ( CheckIfBefore( pl, vEdges[nPos1])) {
         return true ;
      }
      else
         return false ;
   }
}

//----------------------------------------------------------------------------
bool
Tree::CheckIfBefore( int nEdge1, const Point3d& ptP1, int nEdge2, const Point3d& ptP2) const
{
   if ( nEdge1 == -1  ||  nEdge2 == -1)
      return false ;
   // questa funzione è pensata in riferimento al lato 3, quindi nessuno dei due punti può stare su Edge = 3
   INTVECTOR vEdges = { 7, 0, 4, 1, 5, 2, 6} ;
   // controllo se ptP1, che è su nEdge1, viene prima di ptP2, che è su nEdge2. la partenza è da ptTR e l'arrivo è ptBr
   auto iter1 = find( vEdges.begin(), vEdges.end(), nEdge1) ;
   int nPos1 = distance( vEdges.begin(), iter1) ;
   auto iter2 = find( vEdges.begin(), vEdges.end(), nEdge2) ;
   int nPos2 = distance( vEdges.begin(), iter2) ;
   if ( nPos1 < nPos2)
      return true ;
   else if ( nPos1 > nPos2)
      return false ;
   // ( nPos1 == nPos2)
   else {
      PolyLine pl ;
      pl.AddUPoint( 0, ptP2) ;
      pl.AddUPoint( 1, ptP1) ;
      if ( CheckIfBefore( pl, vEdges[nPos1]))
         return true ;
      else
         return false ;
   }
}

//----------------------------------------------------------------------------
bool
Tree::CheckIfBefore( int nEdge, const Point3d& ptP1, const Point3d& ptP2, int nEdge2) const
{
  // i punti devono essere sullo stesso lato nEdge
  // nEdge2 è di backup, in caso nEdge sia un vertice, per capire di quale lato si tratta
   int nEdgeRef = -1 ;
   if ( nEdge != nEdge2  &&  nEdge2 != -1) {
      if ( min( nEdge, nEdge2) < 4)
         nEdgeRef = min( nEdge, nEdge2) ;
      else if ( AreSameEdge( nEdge, 0)  &&  AreSameEdge( nEdge2, 0))
         nEdgeRef = 0 ;
      else if ( AreSameEdge( nEdge, 1)  &&  AreSameEdge( nEdge2, 1))
         nEdgeRef = 1 ;
      else if ( AreSameEdge( nEdge, 2)  &&  AreSameEdge( nEdge2, 2))
         nEdgeRef = 2 ;
      else if ( AreSameEdge( nEdge, 3)  &&  AreSameEdge( nEdge2, 3))
         nEdgeRef = 3 ;
   }
   else
      nEdgeRef = nEdge ;
  // sul lato nEdge controllo se ptP1 viene prima di ptP2.
  // i lati vengono percorsi in senso antiorario
   if ( AreSameEdge( nEdgeRef, 0)) {
      return ( ptP1.x > ptP2.x) ;
   }
   else if ( AreSameEdge( nEdgeRef, 1)) {
      return ( ptP1.y > ptP2.y) ;
   }
   else if ( AreSameEdge( nEdgeRef, 2)) {
      return ( ptP1.x < ptP2.x) ;
   }
   else if ( AreSameEdge( nEdgeRef, 3)) {
      return ( ptP1.y < ptP2.y) ;
   }
   return false ;
}

//----------------------------------------------------------------------------
bool
Tree::AreSameEdge( int nEdge1, int nEdge2) const
{ 
   if ( nEdge1 == 0)
      return ( nEdge2 == 4  ||  nEdge2 == 0  ||  nEdge2 == 7) ;
   else if ( nEdge1 == 1)
      return ( nEdge2 == 4  ||  nEdge2 == 1  ||  nEdge2 == 5) ;
   else if ( nEdge1 == 2)
      return ( nEdge2 == 6  ||  nEdge2 == 2  ||  nEdge2 == 5) ;
   else if ( nEdge1 == 3)
      return ( nEdge2 == 6  ||  nEdge2 == 3  ||  nEdge2 == 7) ;
   else if ( nEdge1 == 4)
      return ( nEdge2 == 0  ||  nEdge2 == 1  ||  nEdge2 == 7  || nEdge2 == 5) ;
   else if ( nEdge1 == 5)
      return ( nEdge2 == 2  ||  nEdge2 == 1  ||  nEdge2 == 4  || nEdge2 == 6) ;
   else if ( nEdge1 == 6)
      return ( nEdge2 == 2  ||  nEdge2 == 3  ||  nEdge2 == 5  || nEdge2 == 7) ;
   else if ( nEdge1 == 7)
      return ( nEdge2 == 0  ||  nEdge2 == 3  ||  nEdge2 == 4  || nEdge2 == 6) ;
   else
      return false ;
}

//----------------------------------------------------------------------------
bool
Tree::AddVertex( int nId, const PNTMATRIX& vEdgeVertex, const PNTMATRIX& vEdgeVertex3d, PolyLine& plTrimmedPoly, int& c,
                 const Point3d& ptToAdd, PolyLine& plTrimmedPoly3d, bool bForTriangulation, Point3d& ptLast) const
{
   Point3d ptBr = m_mTree.at(nId).GetBottomRight() ;
   Point3d ptTR = m_mTree.at(nId).GetTopRight() ;
   Point3d ptTl = m_mTree.at(nId).GetTopLeft() ;
   Point3d ptBL = m_mTree.at(nId).GetBottomLeft() ;
   int nVertToSkip = m_mTree.at(nId).m_nVertToErase ;
  // se è il primo punto della PolyLine lo aggiungo
   if ( plTrimmedPoly.GetPointNbr() == 0) {
     // se cerco di aggiungere il vertice che devo saltare, non faccio nulla
      if (( nVertToSkip == 0 && AreSamePointXYApprox( ptToAdd, ptBL))  ||
          ( nVertToSkip == 1 && AreSamePointXYApprox( ptToAdd, ptBr))  ||
          ( nVertToSkip == 2 && AreSamePointXYApprox( ptToAdd, ptTR))  ||
          ( nVertToSkip == 3 && AreSamePointXYApprox( ptToAdd, ptTl))) {
        // aggiorno l'ultimo punto aggiunto
         ptLast.Set( ptToAdd.x, ptToAdd.y, 0) ;
         return true ;
      }

      plTrimmedPoly.AddUPoint( c, ptToAdd) ;
      if ( bForTriangulation) {
         Point3d pt3d ;
         GetPoint( ptToAdd.x, ptToAdd.y, pt3d) ;
         plTrimmedPoly3d.AddUPoint( c, pt3d) ;
      }
      ++ c ;
      // aggiorno l'ultimo punto aggiunto
      ptLast.Set( ptToAdd.x, ptToAdd.y, 0) ;
      return true ;
   }

  // verifico di essere allineato con un lato, sennò aggiungo e basta
   Vector3d vDir ;
   if ( ! AreSamePointXYApprox( ptToAdd, ptLast))
      vDir = ptToAdd - ptLast ;
   else
      return false ;
  // se non riesco a normalizzare perché sono troppo vicino ad un vertice allora aggiungo direttamente il vertice
   if ( ! vDir.Normalize()) {
      plTrimmedPoly.EraseLastUPoint() ;
      plTrimmedPoly3d.EraseLastUPoint() ;
      Point3d ptVert ;
      int nVert = -1 ;
      if ( AreSamePointXYApprox( ptToAdd, ptBr)){
         ptVert = ptBr ;
         nVert = 1 ;
      }
      else if ( AreSamePointXYApprox( ptToAdd, ptTR)) {
         ptVert = ptTR ;
         nVert = 2 ;
      }
      else if ( AreSamePointXYApprox( ptToAdd, ptTl)) {
         ptVert = ptTl ;
         nVert = 3 ;
      }
      else if ( AreSamePointXYApprox( ptToAdd, ptBL)) {
         ptVert = ptBL ;
         nVert = 0 ;
      }
      else
         ptVert = ptToAdd ;
      if ( nVert != nVertToSkip) {
         plTrimmedPoly.AddUPoint( c, ptVert) ;
         if ( bForTriangulation){
            Point3d pt3d ; GetPoint( ptVert.x, ptVert.y, pt3d) ;
            plTrimmedPoly3d.AddUPoint( c, pt3d) ;
         }
         ++ c ;
      }
      // aggiorno l'ultimo punto aggiunto
      ptLast.Set( ptToAdd.x, ptToAdd.y, 0) ;
      return true ;
   }
   if ( abs( vDir.x) > 1 - EPS_SMALL  ||  abs( vDir.y) > 1 - EPS_SMALL) {
     // se su un edge devo fare dei controlli
     // edge 0
      if ( ptToAdd.x >= ptBL.x  &&  ptToAdd.x <= ptTR.x  &&  ptToAdd.y == ptTR.y  &&  abs( vDir.x) > 1 - EPS_SMALL) {
         for ( int t = nVertToSkip == 2 ? 0 : 1 ; t < (int)vEdgeVertex[0].size() ; ++ t) {
            Point3d ptIntermed = vEdgeVertex[0][t] ;
            if ( ptIntermed.x > ptToAdd.x  &&  ptIntermed.x < ptLast.x) {
               plTrimmedPoly.AddUPoint( c, ptIntermed) ;
               if ( bForTriangulation)
                  plTrimmedPoly3d.AddUPoint( c, vEdgeVertex3d[0][t]) ;
               ++ c ;
            }
         }
         bool bVert = AreSamePointXYApprox( ptToAdd, ptTl) ;
         if ( ! ( nVertToSkip == 3 && bVert)) {
            plTrimmedPoly.AddUPoint( c, ptToAdd) ;
            if ( bForTriangulation){
               Point3d pt3d ; GetPoint( ptToAdd.x, ptToAdd.y, pt3d) ;
               plTrimmedPoly3d.AddUPoint( c, pt3d) ;
            }
            ++ c ;
         }
      }
     // edge 1
      else if ( ptToAdd.y >= ptBL.y  &&  ptToAdd.y <= ptTR.y  &&  ptToAdd.x == ptBL.x  &&  abs( vDir.y) > 1 - EPS_SMALL) {
         for ( int t = nVertToSkip == 3 ? 0 : 1 ; t < (int)vEdgeVertex[1].size() ; ++ t) {
            Point3d ptIntermed = vEdgeVertex[1][t] ;
            if ( ptIntermed.y > ptToAdd.y  &&  ptIntermed.y < ptLast.y) {
               plTrimmedPoly.AddUPoint( c, ptIntermed) ;
               if ( bForTriangulation)
                  plTrimmedPoly3d.AddUPoint( c, vEdgeVertex3d[1][t]) ;
               ++ c ;
            }
         }
         bool bVert = AreSamePointXYApprox( ptToAdd, ptBL) ;
         if ( ! ( nVertToSkip == 0  &&  bVert)) {
            plTrimmedPoly.AddUPoint( c, ptToAdd) ;
            if ( bForTriangulation) {
               Point3d pt3d ; GetPoint( ptToAdd.x, ptToAdd.y, pt3d) ;
               plTrimmedPoly3d.AddUPoint( c, pt3d) ;
            }
            ++ c ;
         }
      }
     // edge 2
      else if ( ptToAdd.x >= ptBL.x  &&  ptToAdd.x <= ptTR.x  &&  ptToAdd.y == ptBL.y  &&  abs( vDir.x) > 1 - EPS_SMALL) {
         for ( int t = nVertToSkip == 0 ? 0 : 1 ; t < (int)vEdgeVertex[2].size() ; ++ t) {
            Point3d ptIntermed = vEdgeVertex[2][t] ;
            if ( ptIntermed.x < ptToAdd.x  &&  ptIntermed.x > ptLast.x) {
               plTrimmedPoly.AddUPoint( c, ptIntermed) ;
               if ( bForTriangulation)
                  plTrimmedPoly3d.AddUPoint( c, vEdgeVertex3d[2][t]) ;
               ++ c ;
            }
         }
         bool bVert = AreSamePointXYApprox( ptToAdd, ptBr) ;
         if ( ! ( nVertToSkip == 1  &&  bVert)) {
            plTrimmedPoly.AddUPoint( c, ptToAdd) ;
            if ( bForTriangulation){
               Point3d pt3d ; GetPoint( ptToAdd.x, ptToAdd.y, pt3d) ;
               plTrimmedPoly3d.AddUPoint( c, pt3d) ;
            }
            ++ c ;
         }
      }
     // edge 3
      else if ( ptToAdd.y >= ptBL.y  &&  ptToAdd.y <= ptTR.y  &&  ptToAdd.x == ptTR.x  &&  abs( vDir.y) > 1 - EPS_SMALL) {
         for ( int t = nVertToSkip == 1 ? 0 : 1 ; t < (int)vEdgeVertex[3].size() ; ++ t) {
            Point3d ptIntermed = vEdgeVertex[3][t] ;
            if ( ptIntermed.y < ptToAdd.y  &&  ptIntermed.y > ptLast.y) {
               plTrimmedPoly.AddUPoint( c, ptIntermed) ;
               if ( bForTriangulation)
                  plTrimmedPoly3d.AddUPoint( c, vEdgeVertex3d[3][t]) ;
               ++ c ;
            }
         }
         bool bVert = AreSamePointXYApprox( ptToAdd, ptTR) ;
         if ( ! ( nVertToSkip == 2  &&  bVert)) {
            plTrimmedPoly.AddUPoint( c, ptToAdd) ;
            if ( bForTriangulation) {
               Point3d pt3d ; GetPoint( ptToAdd.x, ptToAdd.y, pt3d) ;
               plTrimmedPoly3d.AddUPoint( c, pt3d) ;
            }
            ++ c ;
         }
      }
     // sono allineato con un lato, ma NON sono su un lato
     // aggiungo e basta
      else {
         plTrimmedPoly.AddUPoint( c, ptToAdd) ;
         if ( bForTriangulation) {
            Point3d pt3d ; GetPoint( ptToAdd.x, ptToAdd.y, pt3d) ;
            plTrimmedPoly3d.AddUPoint( c, pt3d) ;
         }
         ++ c ;
      }
   }
  // non su un edge, quindi aggiungo e basta
   else {
      plTrimmedPoly.AddUPoint( c, ptToAdd) ;
      if ( bForTriangulation) {
         Point3d pt3d ; GetPoint( ptToAdd.x, ptToAdd.y, pt3d) ;
         plTrimmedPoly3d.AddUPoint( c, pt3d) ;
      }
      ++ c ;
   }
   // aggiorno l'ultimo punto aggiunto
   ptLast.Set( ptToAdd.x, ptToAdd.y, 0) ;
   return true ;
}

//usando le intersezioni con le celle
//----------------------------------------------------------------------------
bool
Tree::SetRightEdgeIn( int nId)
{
   Cell& cCell = m_mTree[nId] ;
  // categorizzo la cella in base a quanta parte del lato destro è contenuta all'interno delle curve di trim
  // RightEdgeIn -> 0 non contenuto ; 1 contenuto ; 2 in parte contenuto
   int nPass = (int) cCell.m_vInters.size() ;
   if ( nPass == 0) {
      cCell.m_nRightEdgeIn = -1 ;
      return true ;
   }
   bool bDone = false ;
  // se ho solo loop interni devo controllare se il più esterno è CCW ( lato destro esterno) o CW ( lato destro interno)
   if ( cCell.m_nFlag == 2) {
      bool bAllContained = true ;
      bool bContained = false ;
      Inters inA ;
      int nInters = (int) cCell.m_vInters.size() ;
      for ( int n = 0 ; n < nInters ; ++ n) {
         inA = cCell.m_vInters[n] ;
         if ( inA.nIn == -1) {
            if ( ! inA.bCCW) {
               bContained = false ;
               Inters inB = cCell.m_vInters[0] ;
               for ( int c = 0 ; c < nInters ; ++ c) {
                  inB = cCell.m_vInters[c] ;
                  if ( inB.nIn == -1) {
                     if ( inB != inA  &&  inB.nChunk == inA.nChunk  &&  inB.bCCW) {
                        bContained = true ;
                        break ;
                     }
                  }
                  else
                     break ;
               }
               bAllContained = bAllContained && bContained ;
            }
         }
         else
            break ;
      }
      // se i loop CW interni sono tutti contenuti in un loop CW allora il right edge è esterno
      bool bRightEdgeIn = ! bAllContained ;
      if ( bRightEdgeIn)
         cCell.m_nRightEdgeIn = 1 ;
      else
         cCell.m_nRightEdgeIn = 0 ;
      return true ;
   }
   // se ho inters sul lato destro ( nEdge == 3) allora in parte è dentro
   for ( int k = 0 ; k < nPass ; ++ k) {
      if ( cCell.m_vInters[k].nIn == 3  ||  cCell.m_vInters[k].nOut == 3) {
         cCell.m_nRightEdgeIn = 2 ;
         bDone = true ;
         break ;
      }
      // considero anche ingressi/ uscite dai vertici 6 e 7
      // controllo nei vertici
      else if ( cCell.m_vInters[k].nOut == 6  &&  cCell.m_vInters[k].nIn == 7) {
         cCell.m_nRightEdgeIn = 1 ;
         bDone = true ;
         break ;
      }
      else if ( cCell.m_vInters[k].nOut == 7  &&  cCell.m_vInters[k].nIn == 6) {
         cCell.m_nRightEdgeIn = 0 ;
         bDone = true ;
         break ;
      }
   }
   // se non ho inters sul lato destro devo verificare se è tutto dentro o tutto fuori
   if ( ! bDone) {
      // ciclo sulle inters e cerco quella più a destra, se ha un'uscita più bassa di un'entrata allora il right edge è compreso!
      bool bFound = false ;
      bool bRightMost = false ;
      int nEdgeUp = 6, nEdgeDown = 7, nLoop ;
      Point3d ptDown( cCell.GetTopRight().x, cCell.GetBottomLeft().y) ;
      Point3d ptUp = cCell.GetTopRight() ;
      for ( int k = 0 ; k < nPass ; ++ k) {
         if ( cCell.m_vInters[k].nIn != -1) {
            // trovo il loop che ha l'ingresso o l'uscita più a destra
            if (  CheckIfBefore( cCell.m_vInters[k].nIn, cCell.m_vInters[k].vpt[0], nEdgeUp, ptUp) ||
                  CheckIfBefore( cCell.m_vInters[k].nOut, cCell.m_vInters[k].vpt.back(), nEdgeUp, ptUp)  ||
                  CheckIfBefore( nEdgeDown, ptDown, cCell.m_vInters[k].nIn, cCell.m_vInters[k].vpt[0]) ||
                  CheckIfBefore( nEdgeDown, ptDown, cCell.m_vInters[k].nOut, cCell.m_vInters[k].vpt.back())) {
               nLoop = k ;
               if ( CheckIfBefore( cCell.m_vInters[k])) {
                  ptUp = cCell.m_vInters[k].vpt[0] ;
                  nEdgeUp = cCell.m_vInters[k].nIn ;
                  ptDown = cCell.m_vInters[k].vpt.back() ;
                  nEdgeDown = cCell.m_vInters[k].nOut ;
               }
               else {
                  ptUp = cCell.m_vInters[k].vpt.back() ;
                  nEdgeUp = cCell.m_vInters[k].nOut ;
                  ptDown = cCell.m_vInters[k].vpt[0] ;
                  nEdgeDown = cCell.m_vInters[k].nIn ;
               }
               bFound = true ;
            }
         }
         else
            continue ;
      }
      if ( bFound  &&  CheckIfBefore( cCell.m_vInters[nLoop])) {
         bRightMost = true ;
      }
      // se il mio campione attraversa dall'alto al basso allora il lato destro è dentro!!!
      if ( bFound && bRightMost)
         cCell.m_nRightEdgeIn = 1 ;
      else
         cCell.m_nRightEdgeIn = 0 ;
   }

   return ( cCell.m_nRightEdgeIn != -1) ;
}

//----------------------------------------------------------------------------
bool
Tree::CategorizeCell( int nId)
{
   Cell& cCell = m_mTree[nId] ;
   if ( cCell.m_nFlag != -1) {
      return true ;
   }
   INTVECTOR vNeigh, vNeigh1 ;
   GetLeftNeigh( nId, vNeigh) ;
   // mi servono i vicini di sinistra per capire se sono dentro il loop o fuori
   // nRightEdgeIn // 0 right edge fuori, 1 right edge dentro, 2 metà e metà
   // nFlag // 0 fuori, 1 intersecata, 2 contiene loop, 3 = 1 & 2, 4 dentro
   if ( vNeigh.empty()) {
      cCell.m_nFlag = 0 ;
   }
   else if ( vNeigh.size() == 1) {
      if ( m_mTree[vNeigh[0]].m_nRightEdgeIn == 1)
         cCell.m_nFlag = 4 ;
      else if ( m_mTree[vNeigh[0]].m_nRightEdgeIn == 0  ||  m_mTree[vNeigh[0]].m_nRightEdgeIn == -1) {
         // devo verificare se la cella è intersecata
         if ( m_mTree[vNeigh[0]].m_nFlag == 1  ||  m_mTree[vNeigh[0]].m_nFlag == 3)
            cCell.m_nFlag = 0 ;
         else {
            if ( m_mTree[vNeigh[0]].m_nFlag == 4)
               cCell.m_nFlag = 4 ;
            else if ( m_mTree[vNeigh[0]].m_nFlag == 0)
               cCell.m_nFlag = 0 ;
         }
      }
      // se solo parte del right edge del vicino è compreso, allora devo verificare se la cella è contenuta o no
      // guardando nFlag del vicino bottom, che è già categorizzato!
      else if ( m_mTree[vNeigh[0]].m_nRightEdgeIn == 2) {
         GetBottomNeigh( nId, vNeigh1) ;
         if ( ! vNeigh1.empty()) {
            int nNeigh = vNeigh1[0] ;
            if ( m_mTree[nNeigh].m_nFlag == 0  ||  m_mTree[nNeigh].m_nFlag == 2) {
               cCell.m_nFlag = 0 ;
            }
            else if ( m_mTree[nNeigh].m_nFlag == 4) {
               cCell.m_nFlag = 4 ;
            }
            else if ( m_mTree[nNeigh].m_nFlag == 1  ||  m_mTree[nNeigh].m_nFlag == 3) {
               if ( m_mTree[nNeigh].m_nRightEdgeIn == 0  ||  m_mTree[nNeigh].m_nRightEdgeIn == 1){
                  // devo capire dove viene tagliato il vicino bottom per poter trarre qualche informazione
                  // se sul lato top del vicino bottom l'intersezione più vicina alla cella in esame (al suo ptBr) è un'uscita, allora sono dentro
                  // in tutti gli altri casi sono fuori
                  int nPass = int( m_mTree[nNeigh].m_vInters.size()) ;
                  bool bIntersTop = false ;
                  for ( int r = 0 ; r < nPass ; ++ r) {
                     if ( m_mTree[nNeigh].m_vInters[r].nOut == 0  ||  m_mTree[nNeigh].m_vInters[r].nOut == 7) {
                        bIntersTop = true ;
                        break ;
                     }
                  }
                  // se nella cella bottom non ho intersezioni con il lato top allora la cella in esame è fuori
                  if ( ! bIntersTop  &&  m_mTree[nNeigh].m_nRightEdgeIn == 0)
                     cCell.m_nFlag = 0 ;
                  else if ( ! bIntersTop  &&  m_mTree[nNeigh].m_nRightEdgeIn == 1)
                     cCell.m_nFlag = 4 ;
                  else {
                     Point3d ptBr = cCell.GetBottomRight() ;
                     bool bFound = false ;
                     Point3d ptLeftMostInters = m_mTree[nNeigh].GetTopRight() ;
                     for ( int r = 0 ; r < nPass ; ++ r) {
                        Point3d ptIntersOut = m_mTree[nNeigh].m_vInters[r].vpt.back() ;
                        Point3d ptIntersIn = m_mTree[nNeigh].m_vInters[r].vpt[0] ;
                        if ( (m_mTree[nNeigh].m_vInters[r].nOut == 0  ||  m_mTree[nNeigh].m_vInters[r].nOut == 7) &&
                             ( CheckIfBefore( 2, ptBr, ptIntersOut)  ||  AreSamePointApprox( ptBr, ptIntersOut))) {
                           if (ptIntersOut.x <= ptLeftMostInters.x) {
                              ptLeftMostInters = ptIntersOut ;
                              bFound = true ;
                           }
                        }
                        if ( (m_mTree[nNeigh].m_vInters[r].nIn == 0  ||  m_mTree[nNeigh].m_vInters[r].nIn == 7) &&
                              ( CheckIfBefore( 2, ptBr, ptIntersIn)  ||  AreSamePointApprox( ptBr, ptIntersIn))) {
                           if (ptIntersIn.x <= ptLeftMostInters.x) {
                              ptLeftMostInters = ptIntersIn ;
                              bFound = false ;
                           }
                        }
                     }
                     if ( bFound)
                        cCell.m_nFlag = 4 ;
                     else
                        cCell.m_nFlag = 0 ;
                  }
               }
               else if ( m_mTree[nNeigh].m_nRightEdgeIn == 2) {
                  // in questo caso devo verificare le intersezioni sul lato destro e top del vicino bottom
                  // scorro e cerco l'intersezione più vicina al ptBr della cella in esame( percorrendo il bordo della cella bottom in senso orario)
                  // se è un nIn allora la cella in esame è dentro, altrimenti è fuori
                  int nPass = (int) m_mTree[nNeigh].m_vInters.size() ;
                  bool bFound = false ;
                  Point3d ptTopLeftMostInters = m_mTree[nNeigh].GetBottomRight() ;
                  for ( int r = 0 ; r < nPass ; ++ r) {
                     int nOut = m_mTree[nNeigh].m_vInters[r].nOut ;
                     int nIn = m_mTree[nNeigh].m_vInters[r].nIn ;
                     if ( (AreSameEdge( nOut, 3)  || nOut == 0)) {
                        Point3d ptInters = m_mTree[nNeigh].m_vInters[r].vpt.back() ;
                        if ( ( abs( ptInters.x - ptTopLeftMostInters.x) > EPS_SMALL  &&  ptInters.x < ptTopLeftMostInters.x)  ||
                           (abs( ptInters.y - ptTopLeftMostInters.y) > EPS_SMALL  &&  ptInters.y > ptTopLeftMostInters.y)) {
                           ptTopLeftMostInters = ptInters ;
                           bFound = true ;
                        }
                     }
                     if ( (AreSameEdge( nIn, 3)  || nIn == 0)) {
                        Point3d ptInters = m_mTree[nNeigh].m_vInters[r].vpt[0] ;
                        if ( ( abs( ptInters.x - ptTopLeftMostInters.x) > EPS_SMALL  &&  ptInters.x < ptTopLeftMostInters.x)  ||
                           (abs( ptInters.y - ptTopLeftMostInters.y) > EPS_SMALL  &&  ptInters.y > ptTopLeftMostInters.y)) {
                           ptTopLeftMostInters = ptInters ;
                           bFound = false ;
                        }
                     }
                  }
                  if ( bFound)
                     cCell.m_nFlag = 4 ;
                  else
                     cCell.m_nFlag = 0 ;
               }
            }
         }
        // se non ho vicini bottom devo per forza guardare il vicino a sinistra
         else {
            Point3d ptInters = m_mTree[vNeigh[0]].GetTopRight() ;
            int nEdge = 3 ;
            bool bIntersIn = false ;
            for ( int k = 0 ; k < (int) m_mTree[vNeigh[0]].m_vInters.size() ; ++ k) {
               if ( AreSameEdge( m_mTree[vNeigh[0]].m_vInters[k].nIn, 3)) {
                  if ( CheckIfBefore( nEdge, m_mTree[vNeigh[0]].m_vInters[k].vpt[0], ptInters, 3)) {
                     ptInters = m_mTree[vNeigh[0]].m_vInters[k].vpt[0] ;
                     bIntersIn = true ;
                  }
               }
               if ( AreSameEdge( m_mTree[vNeigh[0]].m_vInters[k].nOut, 3)) {
                  if ( CheckIfBefore( nEdge, m_mTree[vNeigh[0]].m_vInters[k].vpt.back(), ptInters, 3)) {
                     ptInters = m_mTree[vNeigh[0]].m_vInters[k].vpt.back() ;
                     bIntersIn = false ;
                  }
               }
            }
            if ( bIntersIn)
               cCell.m_nFlag = 4 ;
            else
               cCell.m_nFlag = 0 ;
         }
      }
   }
   else {
      for ( int r = 0 ; r < (int)vNeigh.size() ; ++ r) {
         if ( m_mTree[vNeigh[r]].m_nRightEdgeIn == 1) {
            cCell.m_nFlag = 4 ;
            break ;
         }
         else if ( m_mTree[vNeigh[r]].m_nRightEdgeIn == 0  ||   m_mTree[vNeigh[r]].m_nRightEdgeIn == -1) {
            if ( m_mTree[vNeigh[r]].m_nFlag == 4) {
               cCell.m_nFlag = 4 ;
               break ;
            }
            else if ( m_mTree[vNeigh[r]].m_nFlag == 0) {
               cCell.m_nFlag = 0 ;
               break ;
            }
         }
      }
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
Tree::CheckIfBetween( const Inters& inA, const Inters& inB) const
{
   // controllo se inB è compreso tra l'end e lo start di inA
   // ( dall'end di A percorro i bordi della cella fino a tornare allo start e devo incontrare In e Out di B)
   INTVECTOR vEdges ;
   int nEdge = inA.nOut ;
   if ( nEdge > 3  &&  nEdge != 7) {
      nEdge = nEdge - 4 ; 
   }
   else if ( nEdge  == 7)
      nEdge = 0 ;
   if ( AreSameEdge( inA.nIn, inA.nOut)  &&  ! CheckIfBefore( inA.nIn, inA.vpt[0], inA.vpt.back(), nEdge))
      vEdges.push_back( nEdge) ;
   int nEdgeIn = inA.nIn ;
   if ( nEdgeIn > 3  &&  nEdgeIn != 7) {
      nEdgeIn = nEdgeIn - 4 ; 
   }
   else if ( nEdgeIn  == 7)
      nEdgeIn = 0 ;
   // creo la sequenza di Edges da scorrere per trovare i possibili validNextStart
   while ( ! AreSameEdge( nEdge, nEdgeIn)  ||  vEdges.empty()) {
      vEdges.push_back( nEdge) ;
      if ( nEdge == 3)
         nEdge = 0 ;
      else 
         ++ nEdge ;
   }
   if ( ! AreSameEdge( inA.nIn, inA.nOut))
      vEdges.push_back( nEdge) ;
   bool bFound = false ;
   for ( int i : vEdges) {
      if ( AreSameEdge( inB.nIn, i)) {
         if ( AreSameEdge( inB.nIn, inA.nIn)  &&  AreSameEdge( inA.nIn, inA.nOut)  &&  AreSameEdge( inB.nIn, inA.nOut)) {
            nEdge = inA.nIn ;
            //se l'inizio di A è prima della fine, allora devo controllare che B sia compreso tra Out e In ( esterno)
            if ( CheckIfBefore( nEdge, inA.vpt[0], inA.vpt.back(), inA.nOut)) {
               if ( CheckIfBefore( nEdge, inA.vpt.back(), inB.vpt[0], inA.nOut )  ||  CheckIfBefore( nEdge, inB.vpt[0], inA.vpt[0], inA.nOut))
                  bFound = true ;
            }
            // se l'inizio di A è dopo la fine, allora devo controllare che B sia compreso tra In e Out ( interno)
            else {
               if ( CheckIfBefore( nEdge, inA.vpt.back(), inB.vpt[0], inA.nOut)  &&  CheckIfBefore(  nEdge, inB.vpt[0], inA.vpt[0], inA.nOut))
                  bFound = true ;
            }
         }
         else if ( AreSameEdge( inB.nIn, inA.nOut)) {
            PolyLine pl ;
            pl.AddUPoint( 0, inA.vpt[0]) ;
            pl.AddUPoint( 1, inA.vpt.back()) ;
            Point3d ptEnd ; pl.GetLastPoint( ptEnd) ;
            if ( CheckIfBefore( i, ptEnd, inB.vpt[0]))
               bFound = true ;
         }
         else if ( AreSameEdge( inB.nIn, inA.nIn)) {
            //devo controllare il loop b sia prima dell'inizio di A
            if ( CheckIfBefore( inA.nIn, inB.vpt[0], inA.vpt[0], inB.nIn))
               bFound = true ;
         }
         else
            // devo controllare che inB sia prima di OutB
            if ( AreSameEdge( inB.nOut, inB.nIn)  &&  CheckIfBefore( inB.nOut, inB.vpt[0], inB.vpt.back(), inB.nIn)) {
               bFound = true ;
            }
            else if ( ! AreSameEdge( inB.nOut,inB.nIn))
               bFound = true ;
      }
      // se incontro prima l'uscita dell'ingresso, sicuramente InB NON è un validNextStart
      if ( AreSameEdge( inB.nOut, i)  && ! bFound  &&  CheckIfBefore( i, inA.vpt[0], inB.vpt.back())  &&  CheckIfBefore( i, inA.vpt.back(), inB.vpt.back()))
         break ;
   }
   return bFound ;
}

//----------------------------------------------------------------------------
bool
Tree::CheckIfBefore( int nEdgeA, int nEdgeB) const
{
   // questa funzione è pensata per edge che sono uguali per la funzione AreSameEdge
   // su quell'edge verifico se EdgeA viene prima di EdgeB ( il senso di precorrenza è sempre antiorario)
   switch ( nEdgeA ) {
   case 0 : return nEdgeB == 4 ? true : false ; break ;
   case 1 : return nEdgeB == 5 ? true : false ; break ;
   case 2 : return nEdgeB == 6 ? true : false ; break ;
   case 3 : return nEdgeB == 7 ? true : false ; break ;
   case 4 :
      if ( nEdgeB == 1  ||  nEdgeB == 5)
         return true ;
      break ;
   case 5 :
      if ( nEdgeB == 2  ||  nEdgeB == 6)
         return true ;
      break ;
   case 6 :
      if ( nEdgeB == 3  ||  nEdgeB == 7)
         return true ;
      break ;
   case 7 :
      if ( nEdgeB == 0  ||  nEdgeB == 4)
         return true ;
      break ;
   default : return false ;
   }
   return false ;
}

//----------------------------------------------------------------------------
bool
Tree::GetLeaves( vector<Cell>& vLeaves) const
{
   vLeaves.clear() ;
   for ( int k : m_vnLeaves) {
      Cell cToAdd = m_mTree.at( k) ;
      vLeaves.push_back( cToAdd) ;
   }
   
   return true ;
}

//----------------------------------------------------------------------------
bool
Tree::OnWhichEdge( int nId, const Point3d& ptToAssign, int& nEdge) const
{
   Point3d ptTR = m_mTree.at( nId).GetTopRight() ;
   Point3d ptBL = m_mTree.at( nId).GetBottomLeft() ;
   Point3d ptTl ( ptBL.x, ptTR.y) ;
   Point3d ptBr ( ptTR.x, ptBL.y) ;

   if ( ptToAssign.y < ptBL.y - EPS_SMALL  ||  ptToAssign.y > ptTR.y + EPS_SMALL  || 
        ptToAssign.x < ptBL.x - EPS_SMALL  ||  ptToAssign.x > ptTR.x + EPS_SMALL)
      return false ;
   else if ( AreSamePointXYApprox( ptToAssign, ptTR))
      nEdge = 7 ;
   else if ( AreSamePointXYApprox( ptToAssign, ptTl))
      nEdge = 4 ;
   else if ( AreSamePointXYApprox( ptToAssign, ptBL))
      nEdge = 5 ;
   else if ( AreSamePointXYApprox( ptToAssign, ptBr))
      nEdge = 6 ;
   else if ( ptToAssign.x > ptBL.x  &&  ptToAssign.x < ptTR.x  &&  abs( ptToAssign.y - ptTR.y) < EPS_SMALL)
      nEdge = 0 ;
   else if ( ptToAssign.y > ptBL.y  &&  ptToAssign.y < ptTR.y  &&  abs( ptToAssign.x - ptBL.x) < EPS_SMALL)
      nEdge = 1 ;
   else if ( ptToAssign.x > ptBL.x  &&  ptToAssign.x < ptTR.x  &&  abs( ptToAssign.y - ptBL.y) < EPS_SMALL)
      nEdge = 2 ;
   else if ( ptToAssign.y > ptBL.y  &&  ptToAssign.y < ptTR.y  &&  abs( ptToAssign.x - ptTR.x) < EPS_SMALL)
      nEdge = 3 ;
   else if ( ptToAssign.y > ptBL.y  &&  ptToAssign.y < ptTR.y  &&  ptToAssign.x > ptBL.x  &&  ptToAssign.x < ptTR.x)
      nEdge = -1 ;
   else
      return false ;
   return true ;
}

//----------------------------------------------------------------------------
bool
Tree::GetEdges3D( vector<ICRVCOMPOPOVECTOR>& mCCEdges, POLYLINEVECTOR& vPolygons)
{  
   // se la superficie non è trimmata ricostruisco dalle celle al bordo
   if ( ! m_bTrimmed) {
      INTMATRIX vEdges ; // le righe sono gli edge a partire dallo 0, le colonne sono le celle che compongono quell'edge
      // recupero le celle sui quattro bordi
      vEdges.emplace_back() ;
      GetRootNeigh( 0, vEdges[0]) ;
      // le celle sui bordi orizzontali sono ordinate per x o y crescente, ma i io voglio costruire gli edge in senso antiorario a partire dal ptTR,
      // quindi devo invertire gli Edge 0 e 1
      std::reverse( vEdges[0].begin(), vEdges[0].end()) ;
      vEdges.emplace_back() ;
      GetRootNeigh( 1, vEdges[1]) ;
      std::reverse( vEdges[1].begin(), vEdges[1].end()) ;
      vEdges.emplace_back() ;
      GetRootNeigh( 2, vEdges[2]) ;
      vEdges.emplace_back() ;
      GetRootNeigh( 3, vEdges[3]) ;
   
      // scorro sui gruppi di polyline che rappresentano i poligoni delle celle lungo un lato
      for ( int i = 0 ; i < int( vEdges.size()) ; ++i) {
         mCCEdges.emplace_back() ;
         mCCEdges.back().emplace_back(CreateBasicCurveComposite()) ;
         // scorro sui poligoni delle celle di un lato
         for ( int c = 0 ; c < int( vEdges[i].size()) ; ++c) {
            Cell& cNeigh = m_mTree.at(vEdges[i][c]) ;
            if ( cNeigh.m_vnPolyId.empty())
               continue ;
            PolyLine& plCell = vPolygons[cNeigh.m_vnPolyId[0]] ;
            Point3d pt ; plCell.GetFirstPoint( pt) ;
            Point3d pt3d ;
            // a seconda del lato controllo di stare scorrendo il poligono prendendo solo i punti su quel lato
            if ( i == 0) {
               while ( ! AreSamePointXYApprox(pt, cNeigh.GetTopRight())  &&  plCell.GetNextPoint(pt)) {
                  continue ;
               }
               Point3d pt3d ; GetPoint( cNeigh.GetTopRight().x, cNeigh.GetTopRight().y, pt3d) ;
               mCCEdges.back().back()->AddPoint( pt3d) ;
               // scorro fino alla fine di quel lato
               while ( plCell.GetNextPoint( pt)  &&  ! AreSamePointXYApprox(pt, cNeigh.GetTopLeft())) {
                  GetPoint( pt.x, pt.y, pt3d) ;
                  mCCEdges.back().back()->AddLine( pt3d) ;
               }
               pt3d = ORIG ; GetPoint( cNeigh.GetBottomLeft().x, cNeigh.GetTopRight().y, pt3d) ;
               mCCEdges.back().back()->AddLine( pt3d) ;
               if ( ! mCCEdges.back().back()->IsValid())
                  mCCEdges.back().back()->AddPoint( pt3d) ;
            }
            else if ( i == 1 ) {
               while ( ! AreSamePointXYApprox(pt, cNeigh.GetTopLeft())  &&  plCell.GetNextPoint( pt)) {
                  continue ;
               }
               Point3d pt3d ; GetPoint( cNeigh.GetBottomLeft().x, cNeigh.GetTopRight().y, pt3d) ;
               mCCEdges.back().back()->AddPoint( pt3d) ;
               // scorro fino alla fine di quel lato
               while ( plCell.GetNextPoint( pt)  &&  ! AreSamePointXYApprox(pt, cNeigh.GetBottomLeft())) {
                  GetPoint( pt.x, pt.y, pt3d) ;
                  mCCEdges.back().back()->AddLine( pt3d) ;
               }
               pt3d = ORIG ; GetPoint( cNeigh.GetBottomLeft().x, cNeigh.GetBottomLeft().y, pt3d) ;
               mCCEdges.back().back()->AddLine( pt3d) ;
               if ( ! mCCEdges.back().back()->IsValid())
                  mCCEdges.back().back()->AddPoint( pt3d) ;
            }
            else if ( i == 2) {
               while ( ! AreSamePointXYApprox(pt, cNeigh.GetBottomLeft())  &&  plCell.GetNextPoint( pt)) {
                  continue ;
               }
               Point3d pt3d ; GetPoint( cNeigh.GetBottomLeft().x, cNeigh.GetBottomLeft().y, pt3d) ;
               mCCEdges.back().back()->AddPoint( pt3d) ;
               // scorro fino alla fine di quel lato
               while ( plCell.GetNextPoint( pt)  &&  ! AreSamePointXYApprox(pt, cNeigh.GetBottomRight())) {
                  GetPoint( pt.x, pt.y, pt3d) ;
                  mCCEdges.back().back()->AddLine( pt3d) ;
               }
               pt3d = ORIG ; GetPoint( cNeigh.GetTopRight().x, cNeigh.GetBottomLeft().y, pt3d) ;
               mCCEdges.back().back()->AddLine( pt3d) ;
               if ( ! mCCEdges.back().back()->IsValid())
                  mCCEdges.back().back()->AddPoint( pt3d) ;
            }
            else if ( i == 3) {
               while ( ! AreSamePointXYApprox(pt, cNeigh.GetBottomRight())  &&  plCell.GetNextPoint( pt)) {
                  continue ;
               }
               Point3d pt3d ; GetPoint( cNeigh.GetTopRight().x, cNeigh.GetBottomLeft().y, pt3d) ;
               mCCEdges.back().back()->AddPoint( pt3d) ;
               // scorro fino alla fine di quel lato
               while ( plCell.GetNextPoint( pt)  &&  ! AreSamePointXYApprox(pt, cNeigh.GetTopRight())) {
                  GetPoint( pt.x, pt.y, pt3d) ;
                  mCCEdges.back().back()->AddLine( pt3d) ;
               }
               pt3d = ORIG ; GetPoint( cNeigh.GetTopRight().x, cNeigh.GetTopRight().y, pt3d) ;
               mCCEdges.back().back()->AddLine( pt3d) ;
               if ( ! mCCEdges.back().back()->IsValid())
                  mCCEdges.back().back()->AddPoint( pt3d) ;
            }
         }
      }
   }
   // se la superficie è trimmata ricostruisco dai loop splittati ricostruiti durante il TraceLoop
   else {
      // per ogni edge creo le compo che compongono l'edge dopo i trim ( possono essere più compo separate tra loro)
      for ( int i = 0 ; i < 4 ; ++i) {
         mCCEdges.emplace_back() ;
         INTVECTOR vId ;
         Point3d ptNear = m_mTree.at(-1).GetBottomLeft() ;
         while ( m_vCEdge2D[i].second.GetChainFromNear(ptNear, false, vId)) {
            PtrOwner<ICurveComposite> pCC3D( CreateCurveComposite()) ;
            int nInd = abs( vId[0]) - 1 ;
            Point3d pt2D = m_vCEdge2D[i].first[nInd].first ;
            Point3d pt3D ; GetPoint( pt2D.x, pt2D.y, pt3D) ;
            pCC3D->AddPoint( pt3D) ;
            for ( int j = 1 ; j < ssize( vId) ; ++j) {
               nInd = abs( vId[j]) - 1 ;
               pt2D = m_vCEdge2D[i].first[nInd].second ;
               GetPoint( pt2D.x, pt2D.y, pt3D) ;
               pCC3D->AddLine( pt3D) ;
            }
            if ( ! pCC3D->IsValid())
               pCC3D->AddPoint( pt3D) ;
            // qui devo fare dei controlli prima di aggiungere questa polyline?
            mCCEdges[i].emplace_back( Release( pCC3D)) ; 
         }
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
Tree::AddCutsToRoot( POLYLINEVECTOR& vCuts)
{
   // questa funzione dà per scontato di stare lavorando sulla root di uno spazio parametrico di cui si sta solamente calcolando 
   // il Contour. La funzione da chiamare dopo questa è la CreateCellContour.
   // i tagli sono sempre delle linee che tagliano da parte a parte la cella, quindi sono curve aperte
   if ( m_mTree.size() > 1)
      return false ;
   int nRoot = -1 ;
   for ( int i = 0 ; i < int( vCuts.size()); ++i) {
      PolyLine pl = vCuts.at( i) ;
      m_mTree.at( nRoot).m_vInters.emplace_back() ;
      Point3d pt ; pl.GetFirstPoint( pt) ; 
      int nEdgeIn ;
      // se non trovo il primo punto esattamente su un lato, estendo il primo tratto della polyline all'inditro finchè trovo un'intersezione con un lato
      // questa intersezione diventa il nuovo primo punto del vettore intersezione
      if ( ! OnWhichEdge(nRoot, pt, nEdgeIn) ) {
         Point3d ptSecond ; pl.GetNextPoint( ptSecond) ;
         PtrOwner<ICurveLine> pCL( CreateCurveLine()) ; pCL->SetPVL( ptSecond, pt - ptSecond, 1e6) ;
         PtrOwner<ICurveLine> pCL0( CreateCurveLine()) ; pCL0->Set( m_mTree.at(-1).GetTopRight(), m_mTree.at(-1).GetTopLeft()) ;
         PtrOwner<ICurveLine> pCL1( CreateCurveLine()) ; pCL1->Set( m_mTree.at(-1).GetTopLeft(), m_mTree.at(-1).GetBottomLeft()) ;
         PtrOwner<ICurveLine> pCL2( CreateCurveLine()) ; pCL2->Set( m_mTree.at(-1).GetBottomLeft(), m_mTree.at(-1).GetBottomRight()) ;
         PtrOwner<ICurveLine> pCL3( CreateCurveLine()) ; pCL3->Set( m_mTree.at(-1).GetBottomRight(), m_mTree.at(-1).GetTopRight()) ;
         IntersCurveCurve icc0( *pCL, *pCL0) ;
         IntersCurveCurve icc1( *pCL, *pCL1) ;
         IntersCurveCurve icc2( *pCL, *pCL2) ;
         IntersCurveCurve icc3( *pCL, *pCL3) ;
         IntCrvCrvInfo iccInfo ;
         if (icc0.GetIntersCount() != 0)
            icc0.GetIntCrvCrvInfo( 0, iccInfo) ;
         else if (icc1.GetIntersCount() != 0)
            icc1.GetIntCrvCrvInfo( 0, iccInfo) ;
         else if (icc2.GetIntersCount() != 0)
            icc2.GetIntCrvCrvInfo( 0, iccInfo) ;
         else if (icc3.GetIntersCount() != 0)
            icc3.GetIntCrvCrvInfo( 0, iccInfo) ;
         pt = iccInfo.IciA[0].ptI ;
         if ( ! OnWhichEdge(nRoot, pt, nEdgeIn))
            return false ;
      }
      m_mTree.at( nRoot).m_vInters.back().nIn = nEdgeIn ;
      PNTVECTOR vInters ;
      vInters.emplace_back( pt) ;
      while ( pl.GetNextPoint( pt))
         vInters.emplace_back( pt) ;
      pl.GetLastPoint( pt) ;
      int nEdgeOut ;
      // se non trovo l'ultimo punto esattamente su un lato, estendo l'ultimo tratto della polyline in avnti finchè trovo un'intersezione con un lato
      // questa intersezione diventa il nuovo ultimo punto del vettore intersezione
      if ( ! OnWhichEdge(nRoot, pt, nEdgeOut) ) {
         Point3d ptSecondToLast ; pl.GetNextPoint( ptSecondToLast) ;
         PtrOwner<ICurveLine> pCL( CreateCurveLine()) ; pCL->SetPVL( ptSecondToLast, pt - ptSecondToLast, 1e6) ;
         PtrOwner<ICurveLine> pCL0( CreateCurveLine()) ; pCL0->Set( m_mTree.at(-1).GetTopRight(), m_mTree.at(-1).GetTopLeft()) ;
         PtrOwner<ICurveLine> pCL1( CreateCurveLine()) ; pCL1->Set( m_mTree.at(-1).GetTopLeft(), m_mTree.at(-1).GetBottomLeft()) ;
         PtrOwner<ICurveLine> pCL2( CreateCurveLine()) ; pCL2->Set( m_mTree.at(-1).GetBottomLeft(), m_mTree.at(-1).GetBottomRight()) ;
         PtrOwner<ICurveLine> pCL3( CreateCurveLine()) ; pCL3->Set( m_mTree.at(-1).GetBottomRight(), m_mTree.at(-1).GetTopRight()) ;
         IntersCurveCurve icc0( *pCL, *pCL0) ;
         IntersCurveCurve icc1( *pCL, *pCL1) ;
         IntersCurveCurve icc2( *pCL, *pCL2) ;
         IntersCurveCurve icc3( *pCL, *pCL3) ;
         IntCrvCrvInfo iccInfo ;
         if (icc0.GetIntersCount() != 0)
            icc0.GetIntCrvCrvInfo( 0, iccInfo) ;
         else if (icc1.GetIntersCount() != 0)
            icc1.GetIntCrvCrvInfo( 0, iccInfo) ;
         else if (icc2.GetIntersCount() != 0)
            icc2.GetIntCrvCrvInfo( 0, iccInfo) ;
         else if (icc3.GetIntersCount() != 0)
            icc3.GetIntCrvCrvInfo( 0, iccInfo) ;
         pt = iccInfo.IciA[0].ptI ;
         vInters.pop_back() ;
         vInters.emplace_back( pt) ;
         if ( ! OnWhichEdge(nRoot, pt, nEdgeOut))
            return false ;
      }
      m_mTree.at( nRoot).m_vInters.back().vpt = vInters ;
      m_mTree.at( nRoot).m_vInters.back().nOut = nEdgeOut ;
   }
   // chiamo una funzione per renderli coerenti
   AdjustCuts() ;
   return true ;
}

//----------------------------------------------------------------------------
bool
Tree::AdjustCuts( void)
{
   // correzione del verso dei tagli aperti
   if ( int( m_mTree.at( -1).m_vInters.size()) == 1)
      return true ;
   // li riordino per ordine di quali taglio incontrerei percorrendo il bordo della cella a partire da ptTR
   std::sort( m_mTree.at( -1).m_vInters.begin(), m_mTree.at( -1).m_vInters.end(), [](Inters& a, Inters& b){ return Inters::FirstEncounter(a,b) ;}) ;
   // ora controllo che le intersezioni che trovo siano ingressi alternati ad uscite, sennò inverto l'intersezione
   bool bPreviousWasStart =  m_mTree.at( -1).m_vInters.at(0).bSortedbyStart ;
   for ( int i = 0 ; i < int(  m_mTree.at( -1).m_vInters.size()); ++i) {
      if ( m_mTree.at( -1).m_vInters.at(i).bSortedbyStart == bPreviousWasStart) {
         std::reverse( m_mTree.at( -1).m_vInters.at(i).vpt.begin(), m_mTree.at( -1).m_vInters.at(i).vpt.end()) ;
         int nEdgeOutNew = m_mTree.at( -1).m_vInters.at(i).nIn ;
         m_mTree.at( -1).m_vInters.at(i).nIn = m_mTree.at( -1).m_vInters.at(i).nOut ;
         m_mTree.at( -1).m_vInters.at(i).nOut = nEdgeOutNew ;
         bPreviousWasStart = m_mTree.at( -1).m_vInters.at(i).bSortedbyStart ;
      }
      else
         bPreviousWasStart = ! m_mTree.at( -1).m_vInters.at(i).bSortedbyStart ;
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
Tree::CreateCellContour( POLYLINEMATRIX& vPolygons)
{
  // questa funzione è pensata per essere chiamata dopo la AddCutsToRoot, per creare il poligono di un'unica cella a cui sono stati aggiunti dei tagli
   if ( m_mTree.size() > 1)
      return false ;
   int nRoot = -1 ;
  // preparo tutto per poter chiamare la createCellPolygon
   m_vnLeaves.push_back( nRoot) ;
   INTVECTOR vToCheck( m_mTree.at(nRoot).m_vInters.size()) ;
   iota( vToCheck.begin(),  vToCheck.end(), 0) ;
   int nPoly = 0 ;
   INTVECTOR vnParentChunk ;
   PolyLine pl ;
   pl.AddUPoint(0, m_mTree.at(nRoot).GetTopRight()) ;
   pl.AddUPoint(1, m_mTree.at(nRoot).GetTopLeft()) ;
   pl.AddUPoint(2, m_mTree.at(nRoot).GetBottomLeft()) ;
   pl.AddUPoint(3, m_mTree.at(nRoot).GetBottomRight()) ;
   pl.Close() ;
   PolyLine pl3d ;
   Cell& cRoot = m_mTree.at(nRoot) ;
   Point3d pt3d11 ; GetPoint( cRoot.GetTopRight().x, cRoot.GetTopRight().y, pt3d11) ;
   pl3d.AddUPoint(0, pt3d11) ;
   Point3d pt3d01 ; GetPoint( cRoot.GetBottomLeft().x, cRoot.GetTopRight().y, pt3d01) ;
   pl3d.AddUPoint(1, pt3d01) ;
   Point3d pt3d00 ; GetPoint( cRoot.GetBottomLeft().x, cRoot.GetBottomLeft().y, pt3d00) ;
   pl3d.AddUPoint(2, pt3d00) ;
   Point3d pt3d10 ; GetPoint( cRoot.GetTopRight().x, cRoot.GetBottomLeft().y, pt3d10) ;
   pl3d.AddUPoint(3, pt3d10) ;
   pl3d.Close() ;
   // ora posso creare il poligono della cella con i tagli
   POLYLINEMATRIX vPolygons3d ;
   while ( ! vToCheck.empty()) {
      int nPolyBefore = nPoly ;
      CreateCellPolygons( 0, vPolygons, vPolygons3d, vToCheck, nPoly, vnParentChunk, pl, pl3d) ;
      // l'aggiunta di vPolygons3d forse è un problema in questo punto.
      // sarà da valutare quando si aggiungerà la funzione per aggiungere tagli aperti
      // e se si vuole usare la funzione su uno spazio 2D ideale non collegato ad una sup di Bezier
      if ( nPolyBefore == nPoly)
         break ;
   }
   return true ;
}