//----------------------------------------------------------------------------
//  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 <algorithm>
#include "Tree.h"
#include "SurfBezier.h"
#include "GeoConst.h"
#include "CurveLine.h"
#include "CurveComposite.h"
#include "SurfFlatRegion.h"
#include "IntersLineLine.h"
#include "/EgtDev/Include/EGkPolyLine.h"
#include "/EgtDev/Include/EGkDistPointCurve.h"
#include "/EgtDev/Include/EGkCurve.h"

//----------------------------------------------------------------------------
Cell::Cell( void)
   : m_nId( -1),m_nTop ( -2), m_nBottom( -2), m_nLeft( -2), m_nRight ( -2), m_nParent( -2), m_nDepth( 0),
   m_nChild1( -2), m_nChild2( -2), m_nFlag( -1), m_nFlag2( 0), m_nRightEdgeIn( -1), m_ptPbl(ORIG), m_ptPtr(), m_bProcessed(false), m_bSplitVert(true)
{
   Point3d ptTr ( 1 * SBZ_TREG_COEFF, 1 * SBZ_TREG_COEFF) ;
   m_ptPtr = ptTr ;
}


//----------------------------------------------------------------------------
Cell::Cell( Point3d& ptBL, Point3d& ptTR)
   : m_nId( -1),m_nTop ( -2), m_nBottom( -2), m_nLeft( -2), m_nRight ( -2), m_nParent( -2), m_nDepth( 0),
     m_nChild1( -2), m_nChild2( -2), m_nFlag( -1), m_nFlag2( 0), m_nRightEdgeIn( -1), m_ptPbl(ptBL), m_ptPtr(ptTR), m_bProcessed(false), m_bSplitVert(true)
{}


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

//----------------------------------------------------------------------------
inline bool
Cell::IsSame( const Cell& cOtherCell) const
{
   if ( m_nId == cOtherCell.m_nId)
      return true ;
   else
      return false ;
}

//----------------------------------------------------------------------------
bool
Cell::IsLeaf ( void) const
{
   if( m_nChild1 == -2  &&  m_nChild2 == -2)
      return true ;
   else 
      return false ;
}

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

//----------------------------------------------------------------------------
Tree::Tree( const SurfBezier* pSrfBz, bool bSplitPatches)
   : m_bBilinear( false), m_bMulti( false), m_bClosed( false), m_bSplitPatches( bSplitPatches)
{
   SetSurf( pSrfBz, m_bSplitPatches) ;
}

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

//----------------------------------------------------------------------------
void
Tree::SetSurf( const SurfBezier* pSrfBz, bool bSplitPatches)
{
   m_pSrfBz = pSrfBz ;
   m_bSplitPatches = bSplitPatches ;
  // 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 ;
  // recupero i loop di trim e li divido per chunk
   if ( m_bTrimmed) {
      int nLoop = 0 ;
      INTVECTOR vChunk ;
      // dalla superficie recupero un vettore con i loop di trim
      // 
      //PtrOwner<CurveComposite> pLoop( GetBasicCurveComposite( m_pSrfBz->GetLoop( nLoop))) ;
      //while ( ! IsNull( pLoop)  &&  pLoop->IsValid()) {
      //   if ( nLoop == 0) {
      //      m_vLoop.push_back( Release( pLoop)) ;
      //      vChunck.push_back( nLoop) ;
      //      ++ nLoop ;
      //   }
      //   else {
      //      // se il loop corrente è contenuto nel precedente e sono girati in verso opposto, allora lo aggiungo allo stesso chunk
      //      if ( ) {

      //      }
      //      // salvo il precedente chunk e creo un nuovo chunk
      //      else {
      //         m_vChunk.push_back( vChunk) ;
      //         vChunk.clear() ;
      //         vChunk.push_back( nLoop) ;
      //      }
      //      ++ nLoop ;
      //      PtrOwner<CurveComposite> pLoop( GetBasicCurveComposite( m_pSrfBz->GetLoop( nLoop))) ;
      //   }
      //}

      // recupero la superficie di trim per avere accesso diretto ai loop e mantenendo le informazioni sui chunk
       
      Frame3d frSurf ;
      //const SurfFlatRegion* pTrimReg_( m_pSrfBz->GetTrimRegion()) ;
      //PtrOwner<SurfFlatRegion> pTrimReg( pTrimReg_->Clone()) ;
      PtrOwner<SurfFlatRegion> pTrimReg( m_pSrfBz->GetTrimRegion()->Clone()) ;
      for ( int i = 0 ; i < pTrimReg->GetChunkCount() ; ++ i) {
         PtrOwner<SurfFlatRegion> pChunk( pTrimReg->CloneChunk( i)) ;
         //pChunk->Scale( frSurf, 1./ SBZ_TREG_COEFF, 1./ SBZ_TREG_COEFF, 1) ;
         for ( int j = 0 ; j < pChunk->GetLoopCount( 0) ; ++ j) {
            vChunk.push_back( nLoop) ;
            // i chunk della superficie orignale non possono avere dei sub chunk
            PtrOwner<ICurve> pLoop (pChunk->GetLoop(0, j)) ;
            m_vLoop.emplace_back( Release( pLoop )) ;
            m_mChunk[nLoop] = i ;
            ++ nLoop ;
         }
         m_vChunk.push_back( vChunk) ;
         vChunk.clear() ;
      }
   }
  // salvo i vertici 3d della cella root
   Point3d ptBottom ( 0, 0) ;
   Point3d ptTop( nSpanU * SBZ_TREG_COEFF, nSpanV * SBZ_TREG_COEFF) ;
   Cell cRoot( ptBottom, ptTop) ;
   m_mTree.insert( std::pair< int, Cell>( -1, cRoot)) ;
   Point3d ptP00, ptP10, ptP11, ptP01 ;
   bool bOk = false ;
   PNTVECTOR vVert ;
   ptP00 = m_pSrfBz->GetControlPoint( 0, &bOk);
   vVert.push_back( ptP00) ;
   ptP10 = m_pSrfBz->GetControlPoint( nDegU * nSpanU, &bOk) ;
   vVert.push_back( ptP10) ;
   ptP11 = m_pSrfBz->GetControlPoint( ( nDegU * nSpanU + 1) * ( nDegV * nSpanV + 1) - 1, &bOk) ;
   vVert.push_back( ptP11) ;
   ptP01 = m_pSrfBz->GetControlPoint( ( nDegU * nSpanU + 1 ) * ( nDegV * nSpanV), &bOk) ;
   vVert.push_back( ptP01) ;
   m_mVert.insert( std::pair<int, PNTVECTOR>( -1, vVert)) ;
  // se richiesto divido preliminarmente le patches
   if ( bSplitPatches  &&  ( nSpanU > 1  ||  nSpanV > 1)) {
      int nId = -1 ;
      for ( int i = 1 ; i < nSpanU ; ++i) {
         m_mTree[nId].SetSplitDirVert( true) ;
         Split( nId, i * SBZ_TREG_COEFF) ;
         ++ nId ;
         ++ nId ;
      }
      INTVECTOR vLeaves ;
      GetHeightLeaves( -1, vLeaves) ;
      for ( int nId : vLeaves) {
         for ( int j = nSpanV - 1 ; j > 0 ; --j ) {
            m_mTree[nId].SetSplitDirVert( false) ;
            Split( nId, j * SBZ_TREG_COEFF) ;
            nId = m_mTree[nId].m_nChild2 ;
         }
      }
      vLeaves.clear() ;
      GetHeightLeaves( -1, vLeaves) ;
      m_vnParents = vLeaves ;
   }
  // se non ho già splittato le patches, controllo se la superficie è chiusa. In tal caso la splitto sul parametro su cui è chiusa
   else {
     // verifico se la superficie è chiusa ed eventualmente sistemo le adiacenze
      if ( ( AreSamePointApprox( ptP00, ptP01)  &&  AreSamePointApprox( ptP10, ptP11)) ||
           ( AreSamePointApprox( ptP00, ptP10)  &&  AreSamePointApprox( ptP01, ptP11))) {
         m_bClosed = true ;
         if ( AreSamePointApprox(ptP00, ptP01)) {
            m_mTree[-1].m_nTop = -1 ;
            m_mTree[-1].m_nBottom = -1 ;
            m_mTree[-1].SetSplitDirVert( false) ;
            Split(-1) ;
            // qui devo fare il controllo capped
           // 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 bCapped0 = true, bCapped1 = true ;
            Point3d ptV0, ptV1 ;
            // controllo se tutti i punti sull'isoparametrica sono uguali
            for ( int i = 1 ; i < nDegV * nSpanV + 1 ; ++ i) {
               ptV0 = m_pSrfBz->GetControlPoint( i * ( nDegU * nSpanU + 1), &bOk) ;
               bCapped0 = bCapped0  &&  AreSamePointApprox( ptP00, ptV0) ;
               ptV1 = m_pSrfBz->GetControlPoint( ( i + 1) * ( nDegU * nSpanU + 1) - 1, &bOk) ;
               bCapped1 = bCapped1  &&  AreSamePointApprox( ptP10, ptV1) ;
            }
            if ( bCapped0  &&  bCapped1) {
               m_mTree[0].SetSplitDirVert( true) ;
               Split( 0) ;
               m_mTree[1].SetSplitDirVert( true) ;
               Split( 1) ;
            }
         }
         if ( AreSamePointApprox(ptP00, ptP10)) {
            if( (int) m_mTree.size() == 1) {
               m_mTree[-1].m_nLeft = -1 ;
               m_mTree[-1].m_nRight = -1 ;
               m_mTree[-1].SetSplitDirVert( true) ;
               Split( -1) ;
              // 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 bCapped0 = true, bCapped1 = true ;
               Point3d ptU0, ptU1 ;
               // controllo se tutti i punti sull'isoparametrica sono uguali
               for ( int i = 1 ; i < nDegU * nSpanU + 1 ; ++ i) {
                  ptU0 = m_pSrfBz->GetControlPoint( i, &bOk) ;
                  bCapped0 = bCapped0  &&  AreSamePointApprox( ptP00, ptU0) ;
                  ptU1 = m_pSrfBz->GetControlPoint( i + ( nDegU * nSpanU + 1 ) * ( nDegV * nSpanV), &bOk) ;
                  bCapped1 = bCapped1  &&  AreSamePointApprox( ptP01, ptU1) ;
               }
               if ( bCapped0  &&  bCapped1) {
                  m_mTree[0].SetSplitDirVert( false) ;
                  Split( 0) ;
                  m_mTree[1].SetSplitDirVert( false) ;
                  Split( 1) ;
               }
            }
            else if ( (int) m_mTree.size() > 1  && (int) m_mTree.size() < 4) {
               m_mTree[0].m_nLeft = -1 ;
               m_mTree[0].m_nRight = -1 ;
               m_mTree[1].m_nLeft = -1 ;
               m_mTree[1].m_nRight = -1 ;
               m_mTree[0].SetSplitDirVert( true) ;
               Split( 0) ;
               m_mTree[1].SetSplitDirVert( true) ;
               Split( 1) ;
            }
         }
      }
   }
   
  // calcolo e salvo la distanza reale tra i vertici della cella root
   double dLen0 = Dist( ptP00, ptP10) ;
   double dLen1 = Dist( ptP10, ptP11) ;
   double dLen2 = Dist( ptP01, ptP11) ;
   double dLen3 = Dist( ptP00, ptP01) ;
   m_vDim.push_back( ( dLen0 != 0 ? dLen0 : 1)) ;
   m_vDim.push_back( ( dLen1 != 0 ? dLen1 : 1)) ;
   m_vDim.push_back( ( dLen2 != 0 ? dLen2 : 1)) ;
   m_vDim.push_back( ( dLen3 != 0 ? dLen3 : 1)) ;
}

//----------------------------------------------------------------------------
void
Tree::Split( int nId, double dSplitValue)
{
  // controllo che lo split non venga fatto sul lato della cella
   if ( ( m_mTree[nId].IsSplitVert() && dSplitValue > m_mTree[nId].GetBottomLeft().x + EPS_SMALL && dSplitValue < m_mTree[nId].GetTopRight().x - EPS_SMALL)  ||
        ( ! m_mTree[nId].IsSplitVert() && dSplitValue > m_mTree[nId].GetBottomLeft().y + EPS_SMALL && dSplitValue < m_mTree[nId].GetTopRight().y - EPS_SMALL)) {
        // per lo split a parametro libero dovrò impedire che si facciano split troppo vicini al bordo!!!!!!!!!!!!!!!!!!!
      m_mTree[nId].m_dSplit = dSplitValue ;
      Cell cChild1, cChild2 ;
      cChild1.m_nDepth = m_mTree[nId].m_nDepth + 1 ;
      cChild2.m_nDepth = m_mTree[nId].m_nDepth + 1 ;
      int nNodes = (int) m_mTree.size() ;
      cChild1.m_nId = nNodes - 1 ;
      m_mTree[nId].m_nChild1 = nNodes - 1 ;
      cChild2.m_nId = nNodes ;
      m_mTree[nId].m_nChild2 = nNodes ;
      m_mTree.insert( std::pair<int, Cell>( nNodes - 1, cChild1)) ;
      m_mTree.insert( std::pair<int, Cell>( nNodes, cChild2)) ;
      Point3d ptVert1, ptVert2 ;
      PNTVECTOR vVert ;
      m_mVert.insert( std::pair<int, PNTVECTOR>( nNodes - 1, vVert)) ;
      m_mVert.insert( std::pair<int, PNTVECTOR>( nNodes, vVert)) ;
      if ( ! m_mTree[nId].IsSplitVert())
      {
         // la cella figlio 1 è quella sopra
         Point3d ptBL( m_mTree[nId].GetBottomLeft().x, dSplitValue) ;
         m_mTree[m_mTree[nId].m_nChild1].SetBottomLeft( ptBL) ;
         m_mTree[m_mTree[nId].m_nChild1].SetTopRight( m_mTree[nId].GetTopRight()) ;
         m_mTree[m_mTree[nId].m_nChild1].m_nTop = m_mTree[nId].m_nTop ;
         m_mTree[m_mTree[nId].m_nChild1].m_nBottom = m_mTree[nId].m_nChild2 ;
         m_mTree[m_mTree[nId].m_nChild1].m_nLeft = m_mTree[nId].m_nLeft ;
         m_mTree[m_mTree[nId].m_nChild1].m_nRight = m_mTree[nId].m_nRight ;
         Point3d ptTR( m_mTree[nId].GetTopRight().x, dSplitValue) ;
         m_mTree[m_mTree[nId].m_nChild2].SetBottomLeft( m_mTree[nId].GetBottomLeft()) ;
         m_mTree[m_mTree[nId].m_nChild2].SetTopRight( ptTR) ;
         m_mTree[m_mTree[nId].m_nChild2].m_nTop = m_mTree[nId].m_nChild1 ;
         m_mTree[m_mTree[nId].m_nChild2].m_nBottom = m_mTree[nId].m_nBottom ;
         m_mTree[m_mTree[nId].m_nChild2].m_nLeft = m_mTree[nId].m_nLeft ;
         m_mTree[m_mTree[nId].m_nChild2].m_nRight = m_mTree[nId].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
         m_pSrfBz->GetPointD1D2( m_mTree[nId].GetBottomLeft().x / SBZ_TREG_COEFF, dSplitValue / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptVert1) ;
         m_pSrfBz->GetPointD1D2( m_mTree[nId].GetTopRight().x / SBZ_TREG_COEFF, dSplitValue / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptVert2) ;
         m_mVert[nNodes - 1].push_back( ptVert1) ;
         m_mVert[nNodes - 1].push_back( ptVert2) ;
         m_mVert[nNodes - 1].push_back( m_mVert[nId][2]) ;
         m_mVert[nNodes - 1].push_back( m_mVert[nId][3]) ;
         m_mVert[nNodes].push_back( m_mVert[nId][0]) ;
         m_mVert[nNodes].push_back( m_mVert[nId][1]) ;
         m_mVert[nNodes].push_back( ptVert2) ;
         m_mVert[nNodes].push_back( ptVert1) ;
      }
      else {
         // la cella figlio 1 è quella di sinistra
         Point3d ptTR( dSplitValue, m_mTree[nId].GetTopRight().y) ;
         m_mTree[m_mTree[nId].m_nChild1].SetBottomLeft( m_mTree[nId].GetBottomLeft()) ;
         m_mTree[m_mTree[nId].m_nChild1].SetTopRight( ptTR) ;
         m_mTree[m_mTree[nId].m_nChild1].m_nTop = m_mTree[nId].m_nTop ;
         m_mTree[m_mTree[nId].m_nChild1].m_nBottom = m_mTree[nId].m_nBottom ;
         m_mTree[m_mTree[nId].m_nChild1].m_nLeft = m_mTree[nId].m_nLeft ;
         m_mTree[m_mTree[nId].m_nChild1].m_nRight = m_mTree[nId].m_nChild2 ;
         Point3d ptBL( dSplitValue, m_mTree[nId].GetBottomLeft().y) ;
         m_mTree[m_mTree[nId].m_nChild2].SetBottomLeft( ptBL) ;
         m_mTree[m_mTree[nId].m_nChild2].SetTopRight( m_mTree[nId].GetTopRight()) ;
         m_mTree[m_mTree[nId].m_nChild2].m_nTop = m_mTree[nId].m_nTop ;
         m_mTree[m_mTree[nId].m_nChild2].m_nBottom = m_mTree[nId].m_nBottom ;
         m_mTree[m_mTree[nId].m_nChild2].m_nLeft = m_mTree[nId].m_nChild1 ;
         m_mTree[m_mTree[nId].m_nChild2].m_nRight = m_mTree[nId].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
         m_pSrfBz->GetPointD1D2( dSplitValue / SBZ_TREG_COEFF, m_mTree[nId].GetBottomLeft().y / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptVert2) ;
         m_pSrfBz->GetPointD1D2( dSplitValue / SBZ_TREG_COEFF, m_mTree[nId].GetTopRight().y / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptVert1) ;
         m_mVert[nNodes - 1].push_back( m_mVert[nId][0]) ;
         m_mVert[nNodes - 1].push_back( ptVert2) ;
         m_mVert[nNodes - 1].push_back( ptVert1) ;
         m_mVert[nNodes - 1].push_back( m_mVert[nId][3]) ;
         m_mVert[nNodes].push_back( ptVert2) ;
         m_mVert[nNodes].push_back( m_mVert[nId][1]) ;
         m_mVert[nNodes].push_back( m_mVert[nId][2]) ;
         m_mVert[nNodes].push_back( ptVert1) ;
      }
      m_mTree[m_mTree[nId].m_nChild1].SetParent( nId) ;
      m_mTree[m_mTree[nId].m_nChild2].SetParent( nId) ;
   }
}

//----------------------------------------------------------------------------
void
Tree::Split( int nId)
{
   double dValue ;
   if ( m_mTree[nId].IsSplitVert())
      dValue = ( m_mTree[nId].GetBottomLeft().x + m_mTree[nId].GetTopRight().x) / 2 ;
   else
      dValue = ( m_mTree[nId].GetBottomLeft().y + m_mTree[nId].GetTopRight().y) / 2 ;
   Split( nId, dValue) ;
}

//----------------------------------------------------------------------------
bool
Tree::BuildTree_test( double dLinTol_, double dSideMin, double dSideMax)
{
   //int nCToSplit = -1 ;
   //celle 0,1
   m_mTree[-1].SetSplitDirVert( true) ;
   Split( -1) ;
   //celle 2,3
   m_mTree[0].SetSplitDirVert( false) ;
   Split( 0) ;
   //celle 4,5
   m_mTree[2].SetSplitDirVert( false) ;
   Split( 2) ;
   //celle 6,7
   m_mTree[3].SetSplitDirVert( true) ;
   Split( 3) ;
   //celle 8,9
   m_mTree[1].SetSplitDirVert( false) ;
   Split( 1) ;
   //celle 10,11
   m_mTree[8].SetSplitDirVert( true) ;
   Split( 8) ;
   //celle 12,13
   m_mTree[9].SetSplitDirVert( false) ;
   Split( 9) ;
   m_vnLeaves.push_back( 4) ;
   //m_vnLeaves.push_back( 5) ;
   m_vnLeaves.push_back( 6) ;
   //m_vnLeaves.push_back( 7) ;
   //m_vnLeaves.push_back( 10) ;
   m_vnLeaves.push_back( 11) ;
   //m_vnLeaves.push_back( 12) ;
   m_vnLeaves.push_back( 13) ;

   // aggiunta di split
   //celle 14,15
   m_mTree[5].SetSplitDirVert( true) ;
   Split( 5) ;
   m_vnLeaves.push_back( 14) ;
   m_vnLeaves.push_back( 15) ;
   //celle 16,17
   m_mTree[7].SetSplitDirVert( false) ;
   Split( 7) ;
   m_vnLeaves.push_back( 16) ;
   m_vnLeaves.push_back( 17) ;
   //celle 18,19
   m_mTree[12].SetSplitDirVert( true) ;
   Split( 12) ;
   m_vnLeaves.push_back( 18) ;
   m_vnLeaves.push_back( 19) ;
   //celle 20,21
   m_mTree[10].SetSplitDirVert( false) ;
   Split( 10) ;
   m_vnLeaves.push_back( 20) ;
   m_vnLeaves.push_back( 21) ;

   return true ;
}

//----------------------------------------------------------------------------
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 ;
   //double dLinTol = 0.2 ;
   //double dSideMin = 1 ;
   if ( ! m_bTrimmed  ||  m_bTrimmed) {
      if ( ! m_bBilinear) {
         while ( nCToSplit != -2  &&  m_mTree[nCToSplit].IsProcessed() == false) {
           // controllo che la cella non sia già stata preliminarmente splittata
            if ( m_mTree[nCToSplit].IsLeaf()) {
              // calcolo in quale direzione ho più curvatura
              // ptP00P10 è un punto tra P00 e P10
               double dCurvU = 0, dCurvV = 0 ;
               double dLenParU = ( m_mTree[nCToSplit].GetTopRight().x - m_mTree[nCToSplit].GetBottomLeft().x) / SBZ_TREG_COEFF ;
               double dLenParV = ( m_mTree[nCToSplit].GetTopRight().y - m_mTree[nCToSplit].GetBottomLeft().y) / SBZ_TREG_COEFF;
               if ( dLenParU <= 1. / m_nDegV  ||  dLenParV <= 1. / m_nDegU) {
                  double dU = ( m_mTree[nCToSplit].GetTopRight().x + m_mTree[nCToSplit].GetBottomLeft().x) / 2 / SBZ_TREG_COEFF ;
                  double dV = ( m_mTree[nCToSplit].GetTopRight().y + m_mTree[nCToSplit].GetBottomLeft().y) / 2 / SBZ_TREG_COEFF ;
                  double dULoc = 0.5, dVLoc = 0.5 ;
                  Point3d ptPSrf, ptP00P10, ptP10P11, ptP11P01, ptP01P00 ;
                  m_pSrfBz->GetPointD1D2( dU, dV, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptPSrf) ;
                  m_pSrfBz->GetPointD1D2( dU, m_mTree[nCToSplit].GetBottomLeft().y / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptP00P10) ;
                  m_pSrfBz->GetPointD1D2( m_mTree[nCToSplit].GetTopRight().x / SBZ_TREG_COEFF, dV, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptP10P11) ;
                  m_pSrfBz->GetPointD1D2( dU, m_mTree[nCToSplit].GetTopRight().y / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptP11P01) ;
                  m_pSrfBz->GetPointD1D2( m_mTree[nCToSplit].GetBottomLeft().x / SBZ_TREG_COEFF, dV, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptP01P00) ;
                  Point3d ptV = ( 1 - dULoc) * ptP00P10 + dULoc * ptP11P01 ;
                  Point3d ptU = ( 1 - dVLoc) * ptP10P11 + dVLoc * ptP01P00 ;
                  dCurvV = Dist( ptV, ptPSrf) ;
                  dCurvU = Dist( ptU, ptPSrf) ;
               }
               // faccio un'analisi più fine della curvatura se almeno il grado di una curva di uno dei due parametri è alto e 
               // se sto ancora guardando una cella abbastanza grande
               else{
                  Point3d ptPSrf, ptP00P10, ptP10P11, ptP11P01, ptP01P00, ptPSrfMid;
                  double dStep = 1. / m_nDegU ;
                  for ( double k = dStep ; k < 1 ; k = k + dStep) {
                     double dU = ( k * m_mTree[nCToSplit].GetTopRight().x + ( 1 - k) * m_mTree[nCToSplit].GetBottomLeft().x) / SBZ_TREG_COEFF ;
                     double dV = ( m_mTree[nCToSplit].GetTopRight().y + m_mTree[nCToSplit].GetBottomLeft().y) / 2 / SBZ_TREG_COEFF;
                     m_pSrfBz->GetPointD1D2( dU, dV, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptPSrf) ;
                     if ( k == 0.5)
                        ptPSrfMid = ptPSrf ;
                     m_pSrfBz->GetPointD1D2( dU, m_mTree[nCToSplit].GetBottomLeft().y / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptP00P10) ;
                     m_pSrfBz->GetPointD1D2( dU, m_mTree[nCToSplit].GetTopRight().y / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptP11P01) ;
                     CurveLine clV ;
                     clV.Set( ptP00P10, ptP11P01) ;
                     DistPointCurve dpc( ptPSrf, clV) ;
                     double dDist ;
                     dpc.GetDist( dDist) ;
                     dCurvV = std::max( dCurvV, dDist) ;
                  }
                  dStep = 1. / m_nDegV ;
                  for ( double k = dStep ; k < 1 ; k = k + dStep) {
                     double dU = ( m_mTree[nCToSplit].GetTopRight().x + m_mTree[nCToSplit].GetBottomLeft().x) / 2 / SBZ_TREG_COEFF ;
                     double dV = ( k * m_mTree[nCToSplit].GetTopRight().y + ( 1 - k) * m_mTree[nCToSplit].GetBottomLeft().y) / SBZ_TREG_COEFF ;
                     if ( k == 0.5)
                        ptPSrf = ptPSrfMid ;
                     else
                        m_pSrfBz->GetPointD1D2( dU, dV, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptPSrf) ;
                     m_pSrfBz->GetPointD1D2( m_mTree[nCToSplit].GetTopRight().x / SBZ_TREG_COEFF, dV, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptP10P11) ;
                     m_pSrfBz->GetPointD1D2( m_mTree[nCToSplit].GetBottomLeft().x / SBZ_TREG_COEFF, dV, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptP01P00) ;
                     CurveLine clU ;
                     clU.Set( ptP01P00, ptP10P11) ;
                     DistPointCurve dpc( ptPSrf, clU) ;
                     double dDist ;
                     dpc.GetDist( dDist) ;
                     dCurvU = std::max( dCurvU, dDist) ;
                  }
               }

               // 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 ( dCurvV > dCurvU) {
                 // lungo la direzione V ho una curvatura maggiore
                  bVert = false ;
               }
               else {
                 // lungo la direzione U ho una curvatura maggiore
                  bVert = true ;
               }
               m_mTree[nCToSplit].SetSplitDirVert( bVert) ;
                Point3d ptP00, ptP10, ptP11, ptP01 ;
              // distanza reale tra i vertici della cella
               ptP00 = m_mVert[nCToSplit][0] ;
               ptP10 = m_mVert[nCToSplit][1] ;
               ptP11 = m_mVert[nCToSplit][2] ;
               ptP01 = m_mVert[nCToSplit][3] ;
               double dLen0 = Dist( ptP00, ptP10) ;
               double dLen1 = Dist( ptP10, ptP11) ;
               double dLen2 = Dist( ptP01, ptP11) ;
               double dLen3 = Dist( ptP00, ptP01) ;
              // verifico che la cella sia da splittare e che eventualmente sia abbastanza grande da poterlo fare
               double dSideMinVal = 0, dSideMaxVal = 0 ;
               if ( bVert) {
                  if ( dLen0 != 0  &&  dLen2 != 0)
                     dSideMinVal = std::max( dLen0, dLen2) ;
                  else
                     dSideMinVal = std::max( dLen0, dLen2) ;
               }
               else {
                  if ( dLen1 != 0  &&  dLen3 != 0)
                     dSideMinVal = std::max( dLen1, dLen3) ;
                  else
                     dSideMinVal = std::max( dLen1, dLen3) ;
               }
              // calcolo le diagonali per controllare la dimensione massima dei triangoli in cui dividerei la cella
               dSideMaxVal = std::max( Dist( ptP00, ptP11), Dist( ptP10, ptP01)) ;
         
              // se la cella è abbastanza grande da poter essere divisa ancora, calcolo l'errore di approssimazione
               bool bSplit = false ;
               if ( dSideMinVal / 2 >= dSideMin  &&  dSideMaxVal < dSideMax  &&  ( dCurvV > dLinTol  ||  dCurvU > dLinTol)) {
                  CurveLine cl0010, cl0001, cl1011, cl0111 ;
                 // U=0
                  cl0010.Set( ptP00, ptP10) ;
                 // U=1
                  cl0111.Set( ptP01, ptP11) ;
                  Point3d pt0010, pt0111, ptBz0, ptBz1, ptBzV ;
                  int nFlag ;
                  CurveLine clV ;
                 // determino quanti Step fare per ogni direzione parametrica
                  double dDimU = ( dLen0 >= dLen2 ? dLen0 / m_vDim[0] : dLen2 / m_vDim[2]) ;
                  dDimU = ( dDimU > 1 ? 1 : dDimU) ;
                  double dDimV = ( dLen1 >= dLen3 ? dLen1 / m_vDim[1] : dLen3 / m_vDim[3]) ;
                  dDimV = ( dDimV > 1 ? 1 : dDimV) ;
                 // numero di Step per campionare la superficie nelle due direzioni parametriche
                  int nStepsU = int( 51 * dDimU + 5 * ( 1 - dDimU)) ;
                  int nStepsV = int( 51 * dDimV + 5 * ( 1 - dDimV)) ;
                  for ( int u = 0 ; u < nStepsU  &&  ! bSplit ; ++ u) {
                     double dU = double ( u) / double ( nStepsU - 1) ;
                     double dULoc = (( 1 - dU) * m_mTree[nCToSplit].GetBottomLeft().x + dU * m_mTree[nCToSplit].GetTopRight().x) / SBZ_TREG_COEFF ;
                     if ( ! m_pSrfBz->GetPointD1D2( dULoc, m_mTree[nCToSplit].GetBottomLeft().y / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptBz0)  ||
                          ! m_pSrfBz->GetPointD1D2( dULoc, m_mTree[nCToSplit].GetTopRight().y / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptBz1))
                        return false ;
                     DistPointCurve dpc0010( ptBz0, cl0010) ;
                     DistPointCurve dpc0111( ptBz1, cl0111) ;
                     dpc0010.GetMinDistPoint( 0, pt0010, nFlag) ;
                     dpc0111.GetMinDistPoint( 0, pt0111, nFlag) ;
                     clV.Set( pt0010, pt0111) ;
                     for ( int v = 0 ; v < nStepsV ; ++ v) {
                        double dV = double ( v) / double ( nStepsV - 1) ;
                        double dVLoc = ( ( 1 - dV) * m_mTree[nCToSplit].GetBottomLeft().y + dV * m_mTree[nCToSplit].GetTopRight().y) / SBZ_TREG_COEFF ;
                        if ( ! m_pSrfBz->GetPointD1D2( dULoc, dVLoc, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptBzV))
                           return false ;
                        DistPointCurve dpc( ptBzV, clV) ;
                       // distanza di approssimazione locale
                        double dDist ; 
                        dpc.GetDist( dDist) ;
                        if ( dDist > dLinTol) {
                           bSplit = true ;
                           break ;
                        }
                     }
                  }
               }

               if ( bSplit  ||  dSideMaxVal > dSideMax) {
                  m_mTree[nCToSplit].SetSplitDirVert( bVert) ;
                 // effettuo lo split
                  Split( nCToSplit) ;

                 // procedo con lo split del Child1
                  nCToSplit = m_mTree[nCToSplit].m_nChild1 ;
               }
               else {
                 // sono arrivato ad una cella Leaf, quindi salvo la cella
                  m_vnLeaves.push_back( nCToSplit) ;
                  m_mTree[nCToSplit].SetProcessed() ;
                 // risalgo i parent finché non trovo il primo Child2 da processare
                  nCToSplit = m_mTree[nCToSplit].m_nParent ;
                  if ( nCToSplit == -2)
                     return true ;
                  if ( m_mTree[m_mTree[nCToSplit].m_nChild1].IsProcessed()  &&  m_mTree[m_mTree[nCToSplit].m_nChild2].IsProcessed())
                     m_mTree[nCToSplit].SetProcessed() ;
                  while ( m_mTree[m_mTree[nCToSplit].m_nChild2].IsProcessed()) {
                     if ( m_mTree[nCToSplit].m_nParent != -2)
                        nCToSplit = m_mTree[nCToSplit].m_nParent ;
                     if ( m_mTree[m_mTree[nCToSplit].m_nChild1].IsProcessed()  &&  m_mTree[m_mTree[nCToSplit].m_nChild2].IsProcessed())
                        m_mTree[nCToSplit].SetProcessed() ;
                     if ( nCToSplit == -1  &&  m_mTree[m_mTree[nCToSplit].m_nChild2].IsProcessed())
                           break ;
                  }
                  nCToSplit = m_mTree[nCToSplit].m_nChild2 ;
               }
            }
            else {
               nCToSplit = m_mTree[nCToSplit].m_nChild1 ;
            }
         }
         Balance() ; // da implementare quando dividerò ad un parametro a scelta e non a metà
      }
     // bilineare
      else {
         while ( nCToSplit != -2  &&  m_mTree[nCToSplit].IsProcessed() == false) {
            if ( m_mTree[nCToSplit].IsLeaf()) {
              // vertici della cella
               Point3d ptP00, ptP10, ptP11, ptP01 ;
               ptP00 = m_mVert[nCToSplit][0] ;
               ptP10 = m_mVert[nCToSplit][1] ;
               ptP11 = m_mVert[nCToSplit][2] ;
               ptP01 = m_mVert[nCToSplit][3] ;
              // 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 ;
               // calcolo in quale direzione è meglio dividere in base allo stretch
               Point3d ptPSrfU, ptPSrfV ;
               double dU = 0, dV = 0 ;
               double dDistU = 0, dDistV = 0 ;
               PNTVECTOR vPtU, vPtV ;
               if ( ! m_bMulti) {
                  if ( std::max(dLen0, dLen2) > std::max(dLen1, dLen3)) {
                     bVert = true ;
                  }
                  else {
                     bVert = false ;
                  }
               }
               else {
                  for ( double i = 0.25 ; i < 1 ; i = i + 0.25 ) {
                     dU = (( 1 - i) * m_mTree[nCToSplit].GetBottomLeft().x + i * m_mTree[nCToSplit].GetTopRight().x) / SBZ_TREG_COEFF ;
                     dV = (( 1 - i) * m_mTree[nCToSplit].GetBottomLeft().y + i * m_mTree[nCToSplit].GetTopRight().y) / SBZ_TREG_COEFF ;
                     double dVLoc = ( m_mTree[nCToSplit].GetBottomLeft().y + m_mTree[nCToSplit].GetTopRight().y) / 2 / SBZ_TREG_COEFF ;
                     double dULoc = ( m_mTree[nCToSplit].GetBottomLeft().x + m_mTree[nCToSplit].GetTopRight().x) / 2 / SBZ_TREG_COEFF ;
                     m_pSrfBz->GetPointD1D2( dU, dVLoc, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptPSrfU) ;
                     m_pSrfBz->GetPointD1D2( dULoc, dV, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptPSrfV) ;
                     vPtU.push_back( ptPSrfU) ;
                     vPtV.push_back( ptPSrfV) ;
                  }
                 // devo guardare se i tre punti in vPtU e vPtV sono allineati
                  CurveLine clU, clV;
                  clU.Set(vPtU[0], vPtU[1]) ;
                  clV.Set(vPtV[0], vPtV[1]) ;
                  DistPointCurve dpcU( vPtU[2], clU, false) ;
                  DistPointCurve dpcV( vPtV[2], clV, false) ;
                  dpcU.GetDist( dDistU) ;
                  dpcV.GetDist( dDistV) ;
                  if ( dDistU > dDistV) {
                     bVert = true ;
                  }
                  else {
                     bVert = false ;
                  }
               }
            
              // verifico che la cella sia abbastanza grande da poter essere splittata
               double dSideMinVal = 0, dSideMaxVal = 0 ;
               if ( bVert) {
                  if ( dLen0 != 0  &&  dLen2 != 0)
                     dSideMinVal = std::max( dLen0, dLen2) ;
                  else
                     dSideMinVal = std::max( dLen0, dLen2) ;
               }
               else {
                  if ( dLen1 != 0  &&  dLen3 != 0)
                     dSideMinVal = std::max( dLen1, dLen3) ;
                  else
                     dSideMinVal = std::max( dLen1, dLen3) ;
               }
              // calcolo le diagonali per controllare la dimensione massima dei triangoli in cui dividerei la cella
               dSideMaxVal = std::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) * m_mTree[nCToSplit].GetTopRight().x + i * m_mTree[nCToSplit].GetBottomLeft().x) / SBZ_TREG_COEFF ;
                        double dV =  ( ( 1 - j) * m_mTree[nCToSplit].GetTopRight().y + j * m_mTree[nCToSplit].GetBottomLeft().y) / SBZ_TREG_COEFF ;
                        m_pSrfBz->GetPointD1D2( dU, dV, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptPSrf) ;
                        dErr = std::max( abs( DistPointPlane( ptPSrf, plAppr)), dErr) ;
                     }
                  }
               }
               else {
                  dErr = 1. / 4. * ( (ptP00 - ptP01) + (ptP11 - ptP10)).Len() ;
               }
              // se la cella è abbastanza grande da poter essere divisa ancora e devo approssimare meglio, la divido
               if ( dSideMinVal / 2 >= dSideMin  &&  dSideMaxVal < dSideMax  &&  dErr > dLinTol) {
                  m_mTree[nCToSplit].SetSplitDirVert( bVert) ;
                  // effettuo lo split
                  Split( nCToSplit) ;

                  // procedo con lo split del Child1
                  nCToSplit = m_mTree[nCToSplit].m_nChild1 ;
               }
               else {
                  // sono arrivato ad una cella Leaf, quindi salvo la cella
                  m_vnLeaves.push_back( nCToSplit) ;
                  m_mTree[nCToSplit].SetProcessed() ;
                  // risalgo i parent finché non trovo il primo Child2 da processare
                  nCToSplit = m_mTree[nCToSplit].m_nParent ;
                  if ( nCToSplit == -2)
                     return true ;
                  if ( m_mTree[m_mTree[nCToSplit].m_nChild1].IsProcessed()  &&  m_mTree[m_mTree[nCToSplit].m_nChild2].IsProcessed())
                     m_mTree[nCToSplit].SetProcessed() ;
                  while ( m_mTree[m_mTree[nCToSplit].m_nChild2].IsProcessed()) {
                     if ( m_mTree[nCToSplit].m_nParent != -2)
                        nCToSplit = m_mTree[nCToSplit].m_nParent ;
                     if ( m_mTree[m_mTree[nCToSplit].m_nChild1].IsProcessed()  &&  m_mTree[m_mTree[nCToSplit].m_nChild2].IsProcessed())
                        m_mTree[nCToSplit].SetProcessed() ;
                     if ( nCToSplit == -1  &&  m_mTree[m_mTree[nCToSplit].m_nChild2].IsProcessed())
                        break ;
                  }
                  nCToSplit = m_mTree[nCToSplit].m_nChild2 ;
               }
            }
            else {
               nCToSplit = m_mTree[nCToSplit].m_nChild1 ;
            }
         }
      }
   }
  // NON IMPLEMENTATO 
  // se la superficie è trimmata creo l'albero solo all'interno dei bbox delle curve di trim ( se ho solo trim CW allora tengo il bordo)
   else {
      SurfFlatRegion sfrTrimReg ;
   }
   return true ;
}

//----------------------------------------------------------------------------
void
Tree::Balance()
{
   //for ( int i : vCheck ) {
   //   // non ancora implementato
   //   // rendo il tree balanced : ogni foglia deve avere una profondità di +- 1 rispetto alle foglie adiacenti.
   //}

   // al momento il problema viene bypassato in fase di generazione dei poligoni, considerando per ogni cella, oltre ai propri vertici
   // i vertici dei vicini che giacciono sui suoi lati
}

//----------------------------------------------------------------------------
void
Tree::GetTopNeigh( int nId, INTVECTOR& vTopNeighs) const
{
   if ( (int) vTopNeighs.size() == 0) {
      if ( m_mTree.at(nId).m_nTop == -2)
         return ;
      if ( m_mTree.at(m_mTree.at(nId).m_nTop).IsLeaf())
         vTopNeighs.push_back( m_mTree.at(nId).m_nTop) ;
      else {
         if ( m_mTree.at(m_mTree.at(nId).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(m_mTree.at(nId).m_nTop).GetTopRight().x - m_mTree.at(m_mTree.at(nId).m_nTop).GetBottomLeft().x <= 
               m_mTree.at(nId).GetTopRight().x - m_mTree.at(nId).GetBottomLeft().x) {
               vTopNeighs.push_back( m_mTree.at(m_mTree.at(nId).m_nTop).m_nChild1) ;
               vTopNeighs.push_back( m_mTree.at(m_mTree.at(nId).m_nTop).m_nChild2) ;
            }
           // altrimenti solo uno dei figli lo sarà
            else{
               if ( m_mTree.at(m_mTree.at(m_mTree.at(nId).m_nTop).m_nChild1).GetTopRight().x <= m_mTree.at(nId).GetBottomLeft().x ||
                  m_mTree.at(m_mTree.at(m_mTree.at(nId).m_nTop).m_nChild1).GetBottomLeft().x >= m_mTree.at(nId).GetTopRight().x  )
                  vTopNeighs.push_back( m_mTree.at(m_mTree.at(nId).m_nTop).m_nChild2) ;
               else
                  vTopNeighs.push_back( m_mTree.at(m_mTree.at(nId).m_nTop).m_nChild1) ;
            }
         }
         else {
            vTopNeighs.push_back( m_mTree.at(m_mTree.at(nId).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 <=
                  m_mTree.at(nId).GetTopRight().x - m_mTree.at(nId).GetBottomLeft().x) {
                  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 <= m_mTree.at(nId).GetBottomLeft().x ||
                     m_mTree.at(m_mTree.at(i).m_nChild1).GetBottomLeft().x >= m_mTree.at(nId).GetTopRight().x  )
                     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) ;
            }
         }
      }
   }
   std::vector<Cell> vCells ;
   for ( int k : vTopNeighs)
      vCells.push_back( m_mTree.at(k)) ;
   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) const
{
   if ( (int) vBottomNeighs.size() == 0) {
      if ( m_mTree.at(nId).m_nBottom == -2)
         return ;
      if ( m_mTree.at(m_mTree.at(nId).m_nBottom).IsLeaf())
         vBottomNeighs.push_back( m_mTree.at(nId).m_nBottom) ;
      else {
         if ( m_mTree.at(m_mTree.at(nId).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(m_mTree.at(nId).m_nBottom).GetTopRight().x - m_mTree.at(m_mTree.at(nId).m_nBottom).GetBottomLeft().x <= 
               m_mTree.at(nId).GetTopRight().x - m_mTree.at(nId).GetBottomLeft().x) {
               vBottomNeighs.push_back( m_mTree.at(m_mTree.at(nId).m_nBottom).m_nChild1) ;
               vBottomNeighs.push_back( m_mTree.at(m_mTree.at(nId).m_nBottom).m_nChild2) ;
            }
           // altrimenti solo uno dei figli lo sarà
            else{
               if ( m_mTree.at(m_mTree.at(m_mTree.at(nId).m_nBottom).m_nChild1).GetTopRight().x <= m_mTree.at(nId).GetBottomLeft().x ||
                  m_mTree.at(m_mTree.at(m_mTree.at(nId).m_nBottom).m_nChild1).GetBottomLeft().x >= m_mTree.at(nId).GetTopRight().x  )
                  vBottomNeighs.push_back( m_mTree.at(m_mTree.at(nId).m_nBottom).m_nChild2) ;
               else
                  vBottomNeighs.push_back( m_mTree.at(m_mTree.at(nId).m_nBottom).m_nChild1) ;
            }
         }
         else {
            vBottomNeighs.push_back( m_mTree.at(m_mTree.at(nId).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 <=
                  m_mTree.at(nId).GetTopRight().x - m_mTree.at(nId).GetBottomLeft().x) {
                  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 <= m_mTree.at(nId).GetBottomLeft().x  ||
                     m_mTree.at(m_mTree.at(i).m_nChild1).GetBottomLeft().x >= m_mTree.at(nId).GetTopRight().x)
                     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) ;
            }
         }
      }
   }
   std::vector<Cell> vCells ;
   for ( int k : vBottomNeighs)
      vCells.push_back( m_mTree.at(k)) ;
   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) const
{
   if ( (int) vLeftNeighs.size() == 0) {
      if ( m_mTree.at(nId).m_nLeft == -2)
         return ;
      if ( m_mTree.at(m_mTree.at(nId).m_nLeft).IsLeaf())
         vLeftNeighs.push_back( m_mTree.at(nId).m_nLeft) ;
      else {
            if ( ! m_mTree.at(m_mTree.at(nId).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(m_mTree.at(nId).m_nLeft).GetTopRight().y - m_mTree.at(m_mTree.at(nId).m_nLeft).GetBottomLeft().y <= 
                    m_mTree.at(nId).GetTopRight().y - m_mTree.at(nId).GetBottomLeft().y) {
                  vLeftNeighs.push_back( m_mTree.at(m_mTree.at(nId).m_nLeft).m_nChild1) ;
                  vLeftNeighs.push_back( m_mTree.at(m_mTree.at(nId).m_nLeft).m_nChild2) ;
               }
              // altrimenti solo uno dei figli lo sarà
               else{
                  if ( m_mTree.at(m_mTree.at(m_mTree.at(nId).m_nLeft).m_nChild1).GetTopRight().y <= m_mTree.at(nId).GetBottomLeft().y  ||
                       m_mTree.at(m_mTree.at(m_mTree.at(nId).m_nLeft).m_nChild1).GetBottomLeft().y >= m_mTree.at(nId).GetTopRight().y)
                     vLeftNeighs.push_back( m_mTree.at(m_mTree.at(nId).m_nLeft).m_nChild2) ;
                  else
                     vLeftNeighs.push_back( m_mTree.at(m_mTree.at(nId).m_nLeft).m_nChild1) ;
               }
            }
            else {
               vLeftNeighs.push_back( m_mTree.at(m_mTree.at(nId).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 <=
                    m_mTree.at(nId).GetTopRight().y - m_mTree.at(nId).GetBottomLeft().y) {
                  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 <= m_mTree.at(nId).GetBottomLeft().y  ||
                       m_mTree.at(m_mTree.at(i).m_nChild1).GetBottomLeft().y >= m_mTree.at(nId).GetTopRight().y)
                     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) ;
            }
         }
      }
   }
   std::vector<Cell> vCells ;
   for ( int k : vLeftNeighs)
      vCells.push_back( m_mTree.at(k)) ;
   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) const
{
   if ( (int) vRightNeighs.size() == 0) {
      if ( m_mTree.at(nId).m_nRight == -2)
         return ;
      if ( m_mTree.at(m_mTree.at(nId).m_nRight).IsLeaf())
         vRightNeighs.push_back( m_mTree.at(nId).m_nRight) ;
      else {
         if ( ! m_mTree.at(m_mTree.at(nId).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(m_mTree.at(nId).m_nRight).GetTopRight().y - m_mTree.at(m_mTree.at(nId).m_nRight).GetBottomLeft().y <= 
                 m_mTree.at(nId).GetTopRight().y - m_mTree.at(nId).GetBottomLeft().y) {
               vRightNeighs.push_back( m_mTree.at(m_mTree.at(nId).m_nRight).m_nChild1) ;
               vRightNeighs.push_back( m_mTree.at(m_mTree.at(nId).m_nRight).m_nChild2) ;
            }
           // altrimenti solo uno dei figli lo sarà
            else{
               if ( m_mTree.at(m_mTree.at(m_mTree.at(nId).m_nRight).m_nChild1).GetTopRight().y <= m_mTree.at(nId).GetBottomLeft().y  ||
                    m_mTree.at(m_mTree.at(m_mTree.at(nId).m_nRight).m_nChild1).GetBottomLeft().y >= m_mTree.at(nId).GetTopRight().y)
                  vRightNeighs.push_back( m_mTree.at(m_mTree.at(nId).m_nRight).m_nChild2) ;
               else
                  vRightNeighs.push_back( m_mTree.at(m_mTree.at(nId).m_nRight).m_nChild1) ;
            }
         }
         else {
            vRightNeighs.push_back( m_mTree.at(m_mTree.at(nId).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 <=
                    m_mTree.at(nId).GetTopRight().y - m_mTree.at(nId).GetBottomLeft().y) {
                  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 <= m_mTree.at(nId).GetBottomLeft().y  ||
                       m_mTree.at(m_mTree.at(i).m_nChild1).GetBottomLeft().y >= m_mTree.at(nId).GetTopRight().y)
                     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) ;
            }
         }
      }
   }
   std::vector<Cell> vCells ;
   for ( int k : vRightNeighs)
      vCells.push_back( m_mTree.at(k)) ;
   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 ( (int) vnLeaves.size() == 0) {
         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 = std::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 c = 0 ;
   while ( m_mTree.at(nId).m_nParent != nRef) {
      nId = m_mTree.at(nId).m_nParent ;
      ++ c ;
   }
   return c ;
}

struct generator
{
   int value ;
   generator() { value = -1 ;}
   int operator() () { return ++value ; }
} ;

//----------------------------------------------------------------------------
bool
Tree::GetPolygons( std::vector<POLYLINEVECTOR>& vPolygons) {
   if ( (int) m_vPolygons.size() == 0 ) {
      if ( ! m_bTrimmed) {
         vPolygons.clear() ;
         POLYLINEVECTOR vPolygonsBasic ;
         GetPolygonsBasic( vPolygonsBasic) ;
         for ( PolyLine pl : vPolygonsBasic) {
            POLYLINEVECTOR vSinglePolygon ;
            vSinglePolygon.push_back( pl) ;
            vPolygons.push_back( vSinglePolygon) ;
         }
         return  true ;
      }
      // trimmata
      else {
         if ( ! TraceLoopLabelCell())
            return false ;
         POLYLINEVECTOR vPolygonsBasic ;
         GetPolygonsBasic( vPolygonsBasic) ;
         // scorro sui poligoni delle celle non trimmate
         int nCells = (int)vPolygonsBasic.size() ;
         for ( int i = 0 ; i < nCells ; ++i ) {
            
            // costruisco i poligoni partendo dal vettore delle intersezioni, come spiegato a pag15 di Cripps
            int nId = m_vnLeaves[i] ;
            if ( m_mTree[nId].m_nFlag == 4) {
               // vettore dei poligoni ( loop) della cella nId
               POLYLINEVECTOR vCellPolygons ;
               vCellPolygons.push_back( vPolygonsBasic[i]) ;
               vPolygons.push_back( vCellPolygons) ;
            }
            else if ( m_mTree[nId].m_nFlag == 0)
               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 cotruendo nel ciclo attuale e da cui ripasserò dopo
               INTVECTOR vToCheck( (int) m_mTree[nId].m_vInters.size()) ;
               std::generate_n( vToCheck.begin(), (int) m_mTree[nId].m_vInters.size(), generator()) ;
               // 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( nPoly == 0  ||  (int)vToCheck.size() != 0) {
                  int nPolyBefore = nPoly ;
                  CreateCellPolygons( i, vPolygons, vToCheck, nPoly, vnParentChunk, vPolygonsBasic[i]) ;
                  if ( nPolyBefore == nPoly)
                     break ;
               }
               // ora analizzo anche i loop che sono contenuti nella cella
               CreateIslandAndHoles( i, vPolygons, nPoly, vnParentChunk) ;
            }
         }
         return true;
      }
   }
   else {
      vPolygons = m_vPolygons ;
      return true ;
   }
}

//----------------------------------------------------------------------------
bool
Tree::GetPolygonsBasic( POLYLINEVECTOR& vPolygons)
{
   if ( m_vPolygons.empty()) {
      PNTVECTOR vVertices ;
      INTVECTOR vNeigh ;
      bool bBottomRight , bTopLeft ;
      // scorro lungo tutte le celle leaves e oltre agli angoli della cella aggiungo alla polyline anche i vertici sui lati
      for ( int nId : m_vnLeaves) {
         vVertices.clear() ;
         vNeigh.clear() ;
         vVertices.push_back( m_mTree.at(nId).GetBottomLeft()) ;
         GetBottomNeigh( nId, vNeigh) ;
        // aggiungo i vertici che sono sul lato bottom, solo se ho più di un vicino bottom
         if ( (int) vNeigh.size() != 0  &&  (int) vNeigh.size() != 1){
            for ( int j : vNeigh )
               vVertices.push_back( m_mTree.at(j).GetTopRight()) ;
            bBottomRight = true ;
         }
         else
            bBottomRight = false ;
         vNeigh.clear() ;
         GetRightNeigh ( nId, vNeigh) ;
        // aggiungo i vertici che sono sul lato right, solo se ho più di un vicino right
         if ( (int) vNeigh.size() != 0  &&  (int) vNeigh.size() != 1){
            for ( int j : vNeigh )
               vVertices.push_back( m_mTree.at(j).GetBottomLeft()) ;
         }
        // se non l'ho già aggiunto tramite i vicini bottom aggiungo il punto bottom right
         else if ( ! bBottomRight ) {
            Point3d ptBr( m_mTree.at(nId).GetTopRight().x, m_mTree.at(nId).GetBottomLeft().y) ;
            vVertices.push_back( ptBr) ;
         }
         vNeigh.clear() ;
         vVertices.push_back( m_mTree.at(nId).GetTopRight()) ;
         GetTopNeigh ( nId, vNeigh) ;
         std::reverse( vNeigh.begin(), vNeigh.end()) ;
        // aggiungo i vertici che sono sul lato top, solo se ho più di un vicino top
         if ( (int) vNeigh.size() != 0  &&  (int) vNeigh.size() != 1) {
            for ( int j : vNeigh)
               vVertices.push_back( m_mTree.at(j).GetBottomLeft()) ;
            bTopLeft = true ;
         }
         else
            bTopLeft = false ;
         vNeigh.clear() ;
         GetLeftNeigh ( nId, vNeigh) ;
         std::reverse( vNeigh.begin(), vNeigh.end()) ;
        // aggiungo i vertici che sono sul lato left, solo se ho più di un vicino left
         if ( (int) vNeigh.size() != 0  &&  (int) vNeigh.size() != 1) {
            for ( int j : vNeigh)
               vVertices.push_back( m_mTree.at(j).GetTopRight()) ;
         }
        // se non l'ho già aggiunto tramite i vicini top aggiungo il punto top left
         else if ( ! bTopLeft) {
            Point3d ptTl( m_mTree.at(nId).GetBottomLeft().x, m_mTree.at(nId).GetTopRight().y) ;
            vVertices.push_back( ptTl) ;
         }
         vNeigh.clear() ;
         vVertices.push_back( m_mTree.at(nId).GetBottomLeft()) ;
         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, ptP00, ptP10, ptP11, ptP01;
               double dU, dV ;
               dU = ( m_mTree.at(nId).GetBottomLeft().x + m_mTree.at(nId).GetTopRight().x) / 2 / SBZ_TREG_COEFF;
               dV = ( m_mTree.at(nId).GetBottomLeft().y + m_mTree.at(nId).GetTopRight().y) / 2 / SBZ_TREG_COEFF;
               m_pSrfBz->GetPointD1D2( dU, dV, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptPSrf) ;
               ptP00 = m_mVert.at(nId).at(0) ;
               ptP10 = m_mVert.at(nId).at(1) ;
               ptP11 = m_mVert.at(nId).at(2) ;
               ptP01 = m_mVert.at(nId).at(3) ;
               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) ;
               }
            }
         }

         m_vPolygons.emplace_back() ;
         m_vPolygons.back().emplace_back() ;
         for ( int i = 0 ; i < (int) vVertices.size() ; ++i) {
            m_vPolygons.back().back().AddUPoint(i, vVertices.at(i)) ;
         }

      }
   }

   // restituisco i poligoni delle celle del tree nello spazio parametrico
   for ( int t = 0 ; t < (int)m_vPolygons.size() ; ++ t) {
      for ( int z = 0 ; z < (int)m_vPolygons[t].size(); ++z ) {
         vPolygons.push_back( m_vPolygons[t][z]) ;
      }
   }
   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, CurveLine& clTrim) const
{
   INTVECTOR nCells ;
   int nId = -1 ;
   // se fallisce ritorna un vettore vuoto
   // verifico che il punto sia all'interno dello spazio parametrico
   if ( ptToAssign.x < m_mTree.at(-1).GetBottomLeft().x  ||  ptToAssign.x > m_mTree.at( -1).GetTopRight().x  ||
         ptToAssign.y < m_mTree.at(-1).GetBottomLeft().y  ||  ptToAssign.y > m_mTree.at( -1).GetTopRight().y  ) {
      //nCells.push_back( - 2) ;
      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_bSplitPatches  &&  ( m_nSpanU > 1  ||  m_nSpanV > 1)) {
      INTVECTOR nParents = FindCell( ptToAssign, clTrim, m_vnParents) ;
      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 + EPS_SMALL) {
            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 + EPS_SMALL ) {
            nId = m_mTree.at(nId).m_nChild2 ;
         }
         else {
            nId = m_mTree.at(nId).m_nChild1 ;
         }
      }
   }

   nCells.push_back( nId) ;

   Point3d ptBr( m_mTree.at(nId).GetTopRight().x , m_mTree.at(nId).GetBottomLeft().y) ;
   Point3d ptTl( m_mTree.at(nId).GetBottomLeft().x , m_mTree.at(nId).GetTopRight().y) ;
   if ( AreSamePointApprox(ptToAssign, ptTl)  ||
      AreSamePointApprox( ptToAssign, m_mTree.at(nId).GetBottomLeft())  ||
      AreSamePointApprox(ptToAssign, ptBr)  ||
      AreSamePointApprox( ptToAssign, m_mTree.at(nId).GetTopRight()))
   {
      Point3d ptToAssignPlus ;
      double dParam ;
      Vector3d vDir ;
      clTrim.GetParamAtPoint( ptToAssign, dParam, EPS_SMALL) ;
      clTrim.GetPointTang( dParam + EPS_SMALL, ICurve::FROM_MINUS, ptToAssignPlus, vDir) ;
      if ( abs( vDir.x) > 1 - EPS_SMALL  ||  abs( vDir.y) > 1 - EPS_SMALL) {
         vDir.Rotate( Z_AX, -90) ;
         ptToAssignPlus = ptToAssignPlus + vDir * EPS_SMALL ;
      }
      nCells = FindCell( ptToAssignPlus, clTrim) ;
      if ( nCells.empty()) {
         ptToAssignPlus = ptToAssignPlus - 2 * vDir * EPS_SMALL ;
         nCells = FindCell( ptToAssignPlus, clTrim) ;
      }
   }
   return nCells ;
}

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

   // la CurveLine è il segmento di trim su cui giace ptToAssign
   // devo usare ptInters + qualcosa///////////////////////////////////////////////
   Point3d ptIntersPlus ;
   double dParam ;
   Vector3d vDir ;
   cl.GetParamAtPoint( ptToAssign, dParam, EPS_SMALL) ;
   cl.GetPointTang( dParam + EPS_SMALL, ICurve::FROM_MINUS, ptIntersPlus, vDir) ;
   INTVECTOR nCells ;
   int nId = -1 ;
   bool bFound = false ;
   for ( int nCell : vCells ) {
      if ( ptIntersPlus.x >= m_mTree.at(nCell).GetBottomLeft().x  &&  ptIntersPlus.x <= m_mTree.at(nCell).GetTopRight().x  &&
         ptIntersPlus.y >= m_mTree.at(nCell).GetBottomLeft().y  &&  ptIntersPlus.y <= m_mTree.at(nCell).GetTopRight().y) {
         nId = nCell ;
         nCells.push_back( nId) ;
      }
   }
   return nCells ;
}

//----------------------------------------------------------------------------
bool
Tree::TraceLoopLabelCell( void)
{
   // approssimo i loop di trim con delle spezzate
   POLYLINEVECTOR vPlApprox ;
   //for ( PtrOwner<ICurve> pCrv : m_vLoop){
   for ( int p = 0 ; p < (int) m_vLoop.size();  ++ p){
      PtrOwner<ICurve> pCrv( m_vLoop[p]->Clone()) ;
      double dLinTol = 0.01 ; // questo è riferito allo spazio parametrico perché le curve di loop sono già state riscalate!!!!!! 
      double dAngTolDeg = 5 ;
      PolyLine plApprox ;
      int nType = 0 ;
      pCrv->ApproxWithLines( dLinTol,dAngTolDeg, nType, plApprox) ;
      vPlApprox.push_back( plApprox) ;
   }
   // ora le curve di trim hanno lo stesso indice delle loro corrispondenti spezzate nel vettore vPlApprox
   
   double dLinTol = - EPS_SMALL ;  // questo è riferito allo spazio parametrico, quando è già stato riporatato al range 0..1 !!!!!!
                              // il valore è negativo perché voglio considerare contenuto anche un punto che sta su un lato
   POLYLINEVECTOR vplPolygons ;
   GetPolygonsBasic( vplPolygons) ;
   // percorro i loop trovando le interezioni con le celle e riempiendo i vettori m_vInters delle varie celle
   for ( int i = 0 ; i < (int) vPlApprox.size(); ++ i) {
      PolyLine plLoop = vPlApprox[i] ;
     // calcolo se il loop è CCW o Cw
      double dArea ;
      Plane3d plExtPlane ;
      bool bCCW ;
      plLoop.IsClosedAndFlat( plExtPlane, dArea, 50 * EPS_SMALL) ;
      if ( plExtPlane.GetVersN().z > 0 )
         bCCW = true ;
      else
         bCCW = false ;
     // trovo in quale cella è il ptStart
      Point3d ptStart ;
      plLoop.GetFirstPoint( ptStart) ;
      PNTULIST lPt = plLoop.GetUPointList() ;
      PNTULIST:: iterator ptFirst = lPt.begin();
      PNTULIST:: iterator ptSecond = ptFirst ;
      std::advance( ptSecond, 1) ;
      CurveLine clFirst ;
      clFirst.Set( ptFirst->first, ptSecond->first) ;
      // individuo la cella da cui parte il loop
      INTVECTOR nCells = FindCell( ptStart, clFirst) ;
      int nId = nCells.back() ;
      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
      int nPass = (int) m_mTree[nId].m_vInters.size() ;
      m_mTree[nId].m_vInters.emplace_back() ;
      // salvo il verso del loop
      m_mTree[nId].m_vInters.back().bCCW = bCCW ;
      // salvo il chunk del loop
      m_mTree[nId].m_vInters.back().nChunk = m_mChunk[i] ;
      bool bLoopInside = true ;
      Point3d ptCurr ;
      INTVECTOR :: iterator iter = find( m_vnLeaves.begin(), m_vnLeaves.end(), nId) ;
      int nIdPolygon = std::distance( m_vnLeaves.begin(), iter) ;
      bool bEraseNextPoint = false ;
      while ( plLoop.GetNextPoint( ptCurr)) {
         // sto uscendo dalla cella, quindi cerco l'intersezione
         Point3d ptTStart, ptTEnd ;
         plLoop.GetPrevPoint( ptTStart) ;
         plLoop.GetNextPoint( ptTEnd) ;
         CurveLine clTrim ;
         clTrim.Set( ptTStart, ptTEnd) ;
         while( ! IsPointInsidePolyLine( ptCurr, vplPolygons[nIdPolygon], dLinTol)) {  /// qui devo mettere una tolleranza negativa per poter tener conto anche dei punti che sono SULLA curva
            if ( bEraseNextPoint ) {
               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, vptInters, true)) {
               // scarterò il punto molto vicino al lato e tengo solo l'intersezione del trim col lato
               m_mTree[nId].m_vInters.back().vpt.pop_back() ;
               plLoop.GetPrevPoint( ptTEnd) ;
               plLoop.GetPrevPoint( ptTStart) ;
               plLoop.GetNextPoint( ptCurr) ;
               clTrim.Set( ptTStart, ptTEnd) ;
               //clTrim.ExtendEndByLen( EPS_SMALL * 2) ;
               vptInters.clear() ;
               if ( FindInters( nId, clTrim, vptInters, true))
                  return false ;
               bEraseNextPoint = true ;
            }
            // ricalcolo la posizione di nId nel vettore delle foglie
            iter = find( m_vnLeaves.begin(), m_vnLeaves.end(), nId) ;
            nIdPolygon = std::distance( m_vnLeaves.begin(), iter) ;
            // salvo il verso del loop
            m_mTree[nId].m_vInters.back().bCCW = bCCW ;
            // salvo il chunk del loop
            m_mTree[nId].m_vInters.back().nChunk = m_mChunk[i] ;
         }
         // aggiungo la fine del segmento nel vettore delle intersezioni
         vptInters.push_back( ptCurr) ;
      }
      if ( nId == nFirstCell)
         vptInters.pop_back() ;
      m_mTree[nId].m_vInters.back().vpt = vptInters ;
      if ( bLoopInside) {
         //m_mTree[nId].m_vnLoop.push_back( i) ;
         // setto la categoria della cella
         if ( m_mTree[nId].m_nFlag == -1)
            m_mTree[nId].m_nFlag = 2 ;
         else if ( m_mTree[nId].m_nFlag == 1)
            m_mTree[nId].m_nFlag = 3 ;
         // setto i lati di ingresso e uscita a -1 per indicare che ho un loop interno alla cella
         m_mTree[nId].m_vInters.back().nIn = -1 ;
         m_mTree[nId].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 = m_mTree[nId].m_vInters.back().vpt.back() ;
            //sistemo l'ingresso della prima cella
            int nEdge ;
            OnWhichEdge( nFirstCell, ptFirst, nEdge) ;
            m_mTree[nFirstCell].m_vInters[nPass].nIn = nEdge ;
            // sistemo l'uscita dell'ultima cella
            OnWhichEdge( nId, ptLast, nEdge) ;
            m_mTree[nId].m_vInters.back().nOut = nEdge ;
            // sistemo il flag dell'ultima cella
            if ( m_mTree[nId].m_nFlag == -1)
               m_mTree[nId].m_nFlag = 1 ;
            else if ( m_mTree[nId].m_nFlag == 2)
               m_mTree[nId].m_nFlag = 3 ;
         }
         // sono tornato nella cella iniziale, quindi giunto i due vettori intersezione
         else if ( nId == nFirstCell) {
            int nOut = m_mTree[nId].m_vInters[nPass].nOut ;
            m_mTree[nId].m_vInters.back().vpt.insert( m_mTree[nId].m_vInters.back().vpt.end(),
                                                         m_mTree[nId].m_vInters[nPass].vpt.begin(),
                                                         m_mTree[nId].m_vInters[nPass].vpt.end()) ;
            m_mTree[nId].m_vInters[nPass] = m_mTree[nId].m_vInters.back() ;
            m_mTree[nId].m_vInters.pop_back() ;
            // sistemo il lato d'uscita
            m_mTree[nId].m_vInters[nPass].nOut = nOut ;
         }
      }
   }

   // 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 nFlag2 = 0
   ResetTree() ;
   INTVECTOR vNeigh, vFirst ;
   GetRootNeigh( 1, vFirst) ;
   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 ;
   while ( nCell != nLastLeft  &&  ! bAllDone) {
      // categorizzo la cella
      m_mTree[nCell].m_nFlag2 = 1 ;
      CategorizeCell( nCell) ;
      bool bDone = false ;
      int nInProcessing = -1 ;
      int c = 0 ;
      int nDisplay = -1 ;
      // 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 ( ((int)vNeigh.size() > 0  &&  ! bDone)  ||  ! m_mTree[nCell].IsProcessed()) {
         // per debug
         if ( nInProcessing != nCell) {
            nInProcessing = nCell ;
            c = 0 ;
         }
         else
             ++ c ;
         // 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_nFlag2 == 0 ) {
               nCell = vNeigh[i] ;
               nDisplay = m_mTree[nCell].m_nFlag ;
               bProceeded = true ;
               break ;
            }
         }
         if ( ! bProceeded ) {
            m_mTree[nCell].SetProcessed() ;
         }
         else { //categorizzo la cella
            m_mTree[nCell].m_nFlag2 = 1 ;
            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 non processato
            for ( int t: vNeigh) {
               if ( ! m_mTree[t].IsProcessed()) {
                  bDone = false ;
                  break ;
               }
            }
         }
         // per debug
         if ( c > 3)
            break ;
      }
      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
      //bDone = true ;
      c = 0 ;
      while ( (int) vNeigh.size() == 0  ||  m_mTree[nCell].IsProcessed()) {
         // per debug
         if ( nInProcessing != nCell) {
            nInProcessing = nCell ;
            c = 0 ;
         }
         else
            ++ c ;
         // trovo il vicino a sinistra, già categorizzato, più basso
         vNeigh.clear() ;
         GetLeftNeigh( nCell, vNeigh) ;
         for ( int p = 0; p < (int)vNeigh.size() ; ++ p) {
            if ( m_mTree[vNeigh[p]].m_nFlag2 != 0 ) {
               nCell = vNeigh[p] ;
               break ;
            }
         }
         // se non ho vicini a sinistra sono tornato sul lato sinistro e quindi devo procedere alla cella più bassa non processata sul lato sinistro
         if ( vNeigh.empty() ) {
            for ( int p = 0 ; p < (int)vFirst.size() ; ++ p ) {
               if ( m_mTree[vFirst[p]].m_nFlag2 == 0 ) {
                  nCell = vFirst[p] ;
                  break ;
               }
            }
            if ( nCell == nLastLeft  &&  m_mTree[nCell].IsProcessed() ) {
               // categorizzo
               m_mTree[nCell].m_nFlag2 = 1 ;
               CategorizeCell( nCell) ;
               bAllDone = true ;
               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() ) {
               //nCell = vNeigh[k] ;
               bDone = false ;
               break ;
            }
         }
         if ( bDone) {
            m_mTree[nCell].SetProcessed( bDone) ;
            if ( m_mTree[nCell].m_nFlag2 == 0) {
               m_mTree[nCell].m_nFlag2 = 1 ;
               CategorizeCell( nCell) ;
            }
         }
         // per debug
         if ( c > 3)
            break ;
      }
      vNeigh.clear() ;
      GetRightNeigh( nCell, vNeigh) ;
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
Tree::FindInters( int& nId, CurveLine& clTrim, PNTVECTOR& vptInters, bool bFirstInters)
{
   CurveLine clEdge , clEdge2 ;
   Point3d ptStart , ptEnd ;
   clTrim.GetStartPoint( ptStart) ;
   clTrim.GetEndPoint( ptEnd) ;
   // trovo da quale lato sto uscendo
   Point3d ptTR = m_mTree[nId].GetTopRight() ;
   Point3d ptBL = m_mTree[nId].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
   // oltre il 3 sono le celle adiacenti in diagonale al vertice-> 4 corrisponde al ptTl e da lì in senso antiorario
   // -1 se la curva è sempre dentro la cella
   bool bIntersFound = false ;
   Point3d ptInters ;
   if ( ptEnd.y >= ptTR.y  &&  ptEnd.x <= ptTR.x) {
      nEdge = 0 ;
      // lato sopra
      clEdge.Set( ptTR, ptTl) ;
      // lato sinistro
      if ( ptEnd.x < ptBL.x)
         clEdge2.Set( ptTl, ptBL) ;
      else if ( AreSamePointExact( ptEnd, ptTl)  &&  bFirstInters) {
         nEdge = 4 ;
         ptInters = ptTl ;
         bIntersFound = true ;
      }
   }
   else if ( ptEnd.x <= ptBL.x  &&  ptEnd.y <= ptTR.y) {
      nEdge = 1 ;
      // lato sinistro
      clEdge.Set( ptTl, ptBL) ;
      // lato sotto
      if ( ptEnd.y < ptBL.y)
         clEdge2.Set( ptBL, ptBr) ;
      else if ( AreSamePointExact( ptEnd, ptBL)  &&  bFirstInters) {
         nEdge = 5 ;
         ptInters = ptBL ;
         bIntersFound = true ;
      }
   }
   else if ( ptEnd.y <= ptBL.y  &&  ptEnd.x >= ptBL.x) {
      nEdge = 2 ;
      // lato sotto
      clEdge.Set( ptBL, ptBr) ;
      // lato destro
      if ( ptEnd.x > ptTR.x)
         clEdge2.Set( ptBr, ptTR) ;
      else if ( AreSamePointExact( ptEnd, ptBr)  &&  bFirstInters) {
         nEdge = 6 ;
         ptInters = ptBr ;
         bIntersFound = true ;
      }
   }
   else if ( ptEnd.x >= ptTR.x  &&  ptEnd.y >= ptBL.y) {
      nEdge = 3 ;
      // lato desto
      clEdge.Set( ptBr, ptTR) ;
      // lato sopra
      if ( ptEnd.y > ptTR.y)
         clEdge2.Set( ptTR, ptTl) ;
      else if ( AreSamePointExact( ptEnd, ptTR)  &&  bFirstInters) {
         nEdge = 7 ;
         ptInters = ptTR ;
         bIntersFound = true ;
      }
   }
   else
      return false ;


   if ( ! bIntersFound) {
      // intersezione e controlli
      IntersLineLine illExit( clTrim, clEdge, true) ;
      IntCrvCrvInfo aInfo, aInfo2 ;
      bool bIntersOn1Found = true ;
      if ( ! illExit.GetIntCrvCrvInfo( aInfo)) {
         bIntersOn1Found = false ;
         if ( ! clEdge2.IsValid())
            return false ;
      }
      else if ( aInfo.bOverlap  &&  ! bFirstInters) {
         ptInters = aInfo.IciA[1].ptI ;
         bIntersFound = true ;
      }
      if ( clEdge2.IsValid()  &&  ! bIntersFound){
         IntersLineLine illExit2( clTrim, clEdge2, true) ;
         // verifico su quale dei due lati ho l'intersezione
         if ( ! illExit2.GetIntCrvCrvInfo( aInfo2)){
            if ( bIntersOn1Found) {
               // se ho intersezione su Edge1 con sovrapposizione, seleziono o il primo punto o il secondo.
               if ( aInfo.bOverlap  &&  ! bFirstInters)
                  ptInters = aInfo.IciA[1].ptI ;
               else
                  ptInters = aInfo.IciA[0].ptI ;
            }
            else
               return false ;
         }
         else {
            //// solo intersezione sul lato 2
            if ( aInfo2.bOverlap  &&  ! bFirstInters)
               ptInters = aInfo2.IciA[1].ptI ;
            else
               ptInters = aInfo2.IciA[0].ptI ;
         }
      }
      else
         ptInters = aInfo.IciA[0].ptI ;
   }
   OnWhichEdge( nId, ptInters, nEdge) ;
   m_mTree[nId].m_vInters.back().nOut = nEdge ;
   if ( (int)vptInters.size() == 0)
      vptInters.push_back( ptInters) ;
   else if ( ! AreSamePointExact(ptInters , vptInters.back()))
      vptInters.push_back( ptInters) ;
   // salvo il vettore intersezione per la cella e capisco in quale altra cella passare
   if ( (int)vptInters.size() == 1)
      m_mTree[nId].m_vInters.back().vpt.push_back( vptInters[0]) ;
   else
      m_mTree[nId].m_vInters.back().vpt = vptInters ;
   vptInters.clear() ;
   // setto la categoria della cella
   if ( m_mTree[nId].m_nFlag == -1) 
      m_mTree[nId].m_nFlag = 1 ;
   else if ( m_mTree[nId].m_nFlag == 2)
      m_mTree[nId].m_nFlag = 3 ;
   // seleziono la cella successiva da analizzare
   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 ;
         }
      }
      m_mTree[nId].m_vInters.emplace_back() ;
      Point3d ptBr( m_mTree[nId].GetTopRight().x, m_mTree[nId].GetBottomLeft().y) ;
      if ( ptInters.x > m_mTree[nId].GetBottomLeft().x + EPS_SMALL  &&  ptInters.x < m_mTree[nId].GetTopRight().x - EPS_SMALL)
         m_mTree[nId].m_vInters.back().nIn = 2 ;
      else if ( AreSamePointApprox( ptInters, m_mTree[nId].GetBottomLeft()))
         m_mTree[nId].m_vInters.back().nIn = 5 ;
      else if ( AreSamePointApprox( ptInters, ptBr))
         m_mTree[nId].m_vInters.back().nIn = 6 ;
      else
         m_mTree[nId].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 ;
         }
      }
      m_mTree[nId].m_vInters.emplace_back() ;
      Point3d ptBr( m_mTree[nId].GetTopRight().x, m_mTree[nId].GetBottomLeft().y) ;
      if ( ptInters.y > m_mTree[nId].GetBottomLeft().y + EPS_SMALL  &&  ptInters.y < m_mTree[nId].GetTopRight().y - EPS_SMALL)
         m_mTree[nId].m_vInters.back().nIn = 3 ;
      else if ( AreSamePointApprox( ptInters, ptBr))
         m_mTree[nId].m_vInters.back().nIn = 6 ;
      else if ( AreSamePointApprox( ptInters, m_mTree[nId].GetTopRight()))
         m_mTree[nId].m_vInters.back().nIn = 7 ;
      else
         m_mTree[nId].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 ;
         }
      }
      m_mTree[nId].m_vInters.emplace_back() ;
      Point3d ptTl( m_mTree[nId].GetBottomLeft().x, m_mTree[nId].GetTopRight().y) ;
      if ( ptInters.x > m_mTree[nId].GetBottomLeft().x + EPS_SMALL  &&  ptInters.x < m_mTree[nId].GetTopRight().x - EPS_SMALL)
         m_mTree[nId].m_vInters.back().nIn = 0 ;
      else if ( AreSamePointApprox( ptInters, ptTl))
         m_mTree[nId].m_vInters.back().nIn = 4 ;
      else if ( AreSamePointApprox( ptInters, m_mTree[nId].GetTopRight()))
         m_mTree[nId].m_vInters.back().nIn = 7 ;
      else
         m_mTree[nId].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 ;
         }
      }
      m_mTree[nId].m_vInters.emplace_back() ;
      Point3d ptTl( m_mTree[nId].GetBottomLeft().x, m_mTree[nId].GetTopRight().y) ;
      if ( ptInters.y > m_mTree[nId].GetBottomLeft().y + EPS_SMALL  &&  ptInters.y < m_mTree[nId].GetTopRight().y - EPS_SMALL)
         m_mTree[nId].m_vInters.back().nIn = 1 ;
      else if ( AreSamePointApprox( ptInters, m_mTree[nId].GetBottomLeft()))
         m_mTree[nId].m_vInters.back().nIn = 5 ;
      else if ( AreSamePointApprox( ptInters, ptTl))
         m_mTree[nId].m_vInters.back().nIn = 4 ;
      else
         m_mTree[nId].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] ;
         m_mTree[nId].m_vInters.emplace_back() ;
         // controllo se entro in un vertice o a metà lato
         Point3d ptBr( m_mTree[nId].GetTopRight().x, m_mTree[nId].GetBottomLeft().y) ;
         if ( AreSamePointExact( ptInters, ptBr))
            m_mTree[nId].m_vInters.back().nIn = 6 ;
         else if ( AreSamePointExact( ptInters, m_mTree[nId].GetBottomLeft()))
            m_mTree[nId].m_vInters.back().nIn = 5 ;
         else
            m_mTree[nId].m_vInters.back().nIn = 2 ;
      }
      // ingresso da destra
      else if ( ! nPossible1.empty()) {
         nId = nPossible1.back() ;
         m_mTree[nId].m_vInters.emplace_back() ;
         // controllo se entro in un vertice o a metà lato
         Point3d ptBr( m_mTree[nId].GetTopRight().x, m_mTree[nId].GetBottomLeft().y) ;
         if ( AreSamePointExact( ptInters, ptBr))
            m_mTree[nId].m_vInters.back().nIn = 6 ;
         else if ( AreSamePointExact( ptInters, m_mTree[nId].GetTopRight()))
            m_mTree[nId].m_vInters.back().nIn = 7 ;
         else
            m_mTree[nId].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 {
            nId = vNeigh1.back() ;
            m_mTree[nId].m_vInters.emplace_back() ;
            m_mTree[nId].m_vInters.back().nIn = 7 ;
         }
      }
   }
   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] ;
         m_mTree[nId].m_vInters.emplace_back() ;
         // controllo se entro in un vertice o a metà lato
         Point3d ptBr( m_mTree[nId].GetTopRight().x, m_mTree[nId].GetBottomLeft().y) ;
         if ( AreSamePointExact( ptInters, ptBr))
            m_mTree[nId].m_vInters.back().nIn = 6 ;
         else if ( AreSamePointExact( ptInters, m_mTree[nId].GetTopRight()))
            m_mTree[nId].m_vInters.back().nIn = 7 ;
         else
            m_mTree[nId].m_vInters.back().nIn = 3 ;
      }
      // ingresso dall'alto
      else if ( ! nPossible1.empty()) {
         nId = nPossible1[0] ;
         m_mTree[nId].m_vInters.emplace_back() ;
         // controllo se entro in un vertice o a metà lato
         Point3d ptTl( m_mTree[nId].GetBottomLeft().x, m_mTree[nId].GetTopRight().y) ;
         if ( AreSamePointExact( ptInters, ptTl))
            m_mTree[nId].m_vInters.back().nIn = 4 ;
         else if ( AreSamePointExact( ptInters, m_mTree[nId].GetTopRight()))
            m_mTree[nId].m_vInters.back().nIn = 7 ;
         else
            m_mTree[nId].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 {
            nId = vNeigh[0] ;
            m_mTree[nId].m_vInters.emplace_back() ;
            m_mTree[nId].m_vInters.back().nIn = 4 ;
         }
      }
   }
   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() ;
         m_mTree[nId].m_vInters.emplace_back() ;
         // controllo se entro in un vertice o a metà lato
         Point3d ptTl( m_mTree[nId].GetBottomLeft().x, m_mTree[nId].GetTopRight().y) ;
         if ( AreSamePointExact( ptInters, ptTl))
            m_mTree[nId].m_vInters.back().nIn = 4 ;
         else if ( AreSamePointExact( ptInters, m_mTree[nId].GetTopRight()))
            m_mTree[nId].m_vInters.back().nIn = 7 ;
         else
            m_mTree[nId].m_vInters.back().nIn = 0 ;
      }
      // ingresso da sinistra
      else if ( ! nPossible1.empty()) {
         nId = nPossible1[0] ;
         m_mTree[nId].m_vInters.emplace_back() ;
         // controllo se entro in un vertice o a metà lato
         Point3d ptTl( m_mTree[nId].GetBottomLeft().x, m_mTree[nId].GetTopRight().y) ;
         if ( AreSamePointExact( ptInters, ptTl))
            m_mTree[nId].m_vInters.back().nIn = 4 ;
         else if ( AreSamePointExact( ptInters, m_mTree[nId].GetBottomLeft()))
            m_mTree[nId].m_vInters.back().nIn = 5 ;
         else
            m_mTree[nId].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 {
            GetRightNeigh( nId, vNeigh) ;
            nId = vNeigh[0] ;
            m_mTree[nId].m_vInters.emplace_back() ;
            m_mTree[nId].m_vInters.back().nIn = 5 ;
         }
      }
   }
   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() ;
         m_mTree[nId].m_vInters.emplace_back() ;
         // controllo se entro in un vertice o a metà lato
         Point3d ptTl( m_mTree[nId].GetBottomLeft().x, m_mTree[nId].GetTopRight().y) ;
         if ( AreSamePointExact( ptInters, ptTl))
            m_mTree[nId].m_vInters.back().nIn = 4 ;
         else if ( AreSamePointExact( ptInters, m_mTree[nId].GetBottomLeft()))
            m_mTree[nId].m_vInters.back().nIn = 5 ;
         else
            m_mTree[nId].m_vInters.back().nIn = 1 ;
      }
      // ingresso dal basso
      else if ( ! nPossible1.empty()) {
         nId = nPossible1.back() ;
         m_mTree[nId].m_vInters.emplace_back() ;
         // controllo se entro in un vertice o a metà lato
         Point3d ptBr( m_mTree[nId].GetTopRight().x, m_mTree[nId].GetBottomLeft().y) ;
         if ( AreSamePointExact( ptInters, ptBr))
            m_mTree[nId].m_vInters.back().nIn = 6 ;
         else if ( AreSamePointExact( ptInters, m_mTree[nId].GetBottomLeft()))
            m_mTree[nId].m_vInters.back().nIn = 5 ;
         else
            m_mTree[nId].m_vInters.back().nIn = 2 ;
      }
      // ingresso in diagonale
      else {
         if ( ! vNeigh.empty()) {
            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 {
            nId = vNeigh.back() ;
            m_mTree[nId].m_vInters.emplace_back() ;
            m_mTree[nId].m_vInters.back().nIn = 6 ;
         }
      }
   }
   // aggiungo l'intersezione al vettore delle intersezioni della prossima cella
   vptInters.push_back( ptInters) ;
   return true ;
}

//----------------------------------------------------------------------------
bool
Tree::CreateCellPolygons( int nLeafId, std::vector<POLYLINEVECTOR>& vPolygons, INTVECTOR& vToCheck, int& nPoly, INTVECTOR& vnParentChunk, PolyLine plCell)
{
   // conto quanti vertici in più ho per lato e creo un vettore dei vertici per lato
   int nId = m_vnLeaves[nLeafId] ;
   std::vector<PNTVECTOR> vEdgeVertex ;
   Point3d ptTl(m_mTree[nId].GetBottomLeft().x, m_mTree[nId].GetTopRight().y) ;
   Point3d ptBr(m_mTree[nId].GetTopRight().x, m_mTree[nId].GetBottomLeft().y) ;
   vEdgeVertex.emplace_back() ;
   vEdgeVertex.back().push_back( m_mTree[nId].GetTopRight()) ;
   vEdgeVertex.emplace_back() ;
   vEdgeVertex.back().push_back( ptTl) ;
   vEdgeVertex.emplace_back() ;
   vEdgeVertex.back().push_back( m_mTree[nId].GetBottomLeft()) ;
   vEdgeVertex.emplace_back() ;
   vEdgeVertex.back().push_back( ptBr) ;
   // la PolyLine è riempita a partire dal lato bottom
   Point3d ptStart ;
   plCell.GetFirstPoint( ptStart) ;
   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 ;
      while ( plCell.GetNextPoint( ptToAdd)  &&  ! AreSamePointExact( ptToAdd, vEdgeVertex[next][0])) {
         vEdgeVertex[j].push_back( ptToAdd) ;
      }
      //vEdgeVertex[next].push_back( ptToAdd) ;
   }
   
   // comincio a costruire il poligono
   INTVECTOR vToCheckNow = vToCheck ;
   //vToCheck.clear() ;
   // vettore dei poligoni ( loop) della cella nId
   POLYLINEVECTOR vCellPolygons ;
   // costruisco i poligoni partendo dal vettore delle intersezioni, come spiegato a pag15 di Cripps
   PolyLine plTrimmedPoly ;
   // 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 ;
   for ( int w = 0 ; w < (int)vToCheckNow.size() ; ++ w) {
      if ( m_mTree[nId].m_vInters[w].vpt.size() < 2) {
         continue ;
      }
      // indice del loop in m_vInters
      int j = vToCheckNow[w] ;
      Inters inA = m_mTree[nId].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, plTrimmedPoly, c, ptInt) ;
         }
         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 ;
         if ( AreSameEdge( nEdge, 0)) {
            if ( AreSamePointExact(ptLast, ptTl) ) {
               vEdge = m_mTree[nId].GetBottomLeft() - ptTl ;
               vEdge.Normalize() ;
               if ( AreOppositeVectorApprox(vLast, vEdge) ) {
                  plTrimmedPoly.EraseLastUPoint() ;
                  nEdge = 1 ;
               }
            }
         }
         else if ( AreSameEdge( nEdge, 1) && AreSamePointExact(ptLast, m_mTree[nId].GetBottomLeft()) ) {
            vEdge = ptBr - m_mTree[nId].GetBottomLeft() ;
            vEdge.Normalize() ;
            if ( AreOppositeVectorApprox(vLast, vEdge) ) {
               plTrimmedPoly.EraseLastUPoint() ;
               nEdge = 2 ;
            }
         }
         else if ( AreSameEdge( nEdge, 2) ) {
            if ( AreSamePointExact(ptLast, ptBr) ) {
               vEdge = m_mTree[nId].GetTopRight() - ptBr ;
               vEdge.Normalize() ;
               if ( AreOppositeVectorApprox(vLast, vEdge) ) {
                  plTrimmedPoly.EraseLastUPoint() ;
                  nEdge = 3 ;
               }
            }
         }
         else if ( AreSameEdge( nEdge, 3) && AreSamePointExact(ptLast, m_mTree[nId].GetTopRight()) ) {
            vEdge = ptTl - m_mTree[nId].GetTopRight() ;
            vEdge.Normalize() ;
            if ( AreOppositeVectorApprox(vLast, vEdge) ) {
               plTrimmedPoly.EraseLastUPoint() ;
               nEdge = 0 ;
            }
         }

         // 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 ) {
            plTrimmedPoly.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 ;
         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( m_mTree[nId].m_vInters[j], m_mTree[nId].m_vInters[vToCheckNow[t]]) ;
               if ( bValidNextStart) {
                  nNext = t ;
                  bAtNextStart = AreSameEdge( m_mTree[nId].m_vInters[j].nOut, m_mTree[nId].m_vInters[vToCheckNow[t]].nIn) ;
                  break ;
               }
               else {
                  ++ nSecondCheck ;
                  //vToCheck.push_back( vToCheckNow[t]) ;
               }
            }

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

            while ( ! ( bValidNextStart  &&  bAtNextStart)  &&  bNotCameBack) {
               Point3d ptVert ;
               if ( nEdge == 0)
                  ptVert = ptTl ;
               else if ( nEdge == 1)
                  ptVert = m_mTree[nId].GetBottomLeft() ;
               else if ( nEdge == 2)
                  ptVert = ptBr ;
               else if ( nEdge == 3)
                  ptVert = m_mTree[nId].GetTopRight() ;
               AddVertex( nId, vEdgeVertex, plTrimmedPoly, c, ptVert) ;
               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, m_mTree[nId].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) {
               for ( int p = 1 ; p < (int) vEdgeVertex[nEdge].size() ; ++ p) {
                  if ( CheckIfBefore( nEdge, vEdgeVertex[nEdge][p], m_mTree[nId].m_vInters[vToCheckNow[nNext]].vpt[0]) ) {
                     plTrimmedPoly.AddUPoint( c, vEdgeVertex[nEdge][p]) ;
                     ++ c ;
                  }
               }
               w = w + nSecondCheck ;
               continue ;
            }
            // sono tornato indietro
            else if ( ! bNotCameBack){
               Point3d ptStart ;
               plTrimmedPoly.GetFirstPoint( ptStart) ;
               for ( int p = 1 ; p < (int) vEdgeVertex[nEdge].size() ; ++ p) {
                  if ( CheckIfBefore( nEdge, vEdgeVertex[nEdge][p], ptStart) ) {
                     plTrimmedPoly.AddUPoint( c, vEdgeVertex[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 = m_mTree[nId].GetBottomLeft() ;
               else if ( nEdge == 2)
                  ptVert = ptBr ;
               else if ( nEdge == 3)
                  ptVert = m_mTree[nId].GetTopRight() ;
               AddVertex( nId, vEdgeVertex, plTrimmedPoly, c, ptVert) ;
               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)) ;
            }
            for ( int p = 1 ; p < (int) vEdgeVertex[nEdge].size() ; ++ p) {
               if ( CheckIfBefore( nEdge, vEdgeVertex[nEdge][p], ptStart)) {
                  plTrimmedPoly.AddUPoint( c, vEdgeVertex[nEdge][p]) ;
                  ++ c ;
               }
            }

         }
         plTrimmedPoly.Close() ;
         // 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) {
            vCellPolygons.push_back( plTrimmedPoly) ;
            vPolygons.push_back( vCellPolygons) ;
            ++ nPoly ;
            vnParentChunk.push_back( inA.nChunk) ;
            vCellPolygons.clear() ;
         }

         c = 0 ;
         plTrimmedPoly.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 nLeafId, std::vector<POLYLINEVECTOR>& vPolygons, int& nPoly, INTVECTOR& vnParentChunk)
{
   // vettore dei poligoni ( loop) della cella nId
   POLYLINEVECTOR vCellPolygons ;
   // costruisco i poligoni partendo dal vettore delle intersezioni, come spiegato a pag15 di Cripps
   int nId = m_vnLeaves[nLeafId] ;
   PolyLine plTrimmedPoly ;
   // numero di volte che la cella è stata attraversata da una curva di trim
   int nPass = (int) m_mTree[nId].m_vInters.size() ;
   // numero di vertici aggiunti al nuovo poligono
   int c = 0 ;

   // loop interni in una cella intersecata
   if ( m_mTree[nId].m_nFlag == 3  ||  m_mTree[nId].m_nFlag == 2) {
      PolyLine plInLoop ;
      // numero dei loop interni passati
      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 bEdgeToAdd = true ;
      bool bAllContained = true ;
      bool bContained = false ;
      int nInters = (int) m_mTree[nId].m_vInters.size() ;
      for ( int n = 0 ; n < nInters ; ++ n) {
         inA = m_mTree[nId].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) {
               bContained = false ;
               Inters inB = m_mTree[nId].m_vInters[0] ;
               for( int c = 0 ; c < nInters ; ++ c){
                  inB = m_mTree[nId].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 ( m_mTree[nId].m_nFlag == 2  &&  ( bEdgeToAdd  ||  ! bAllContained)) {
      if ( m_mTree[nId].m_nFlag == 2  &&  ! bAllContained) {
         // i loop esterni sono CW, quindi prima dei loop di trim aggiungo il bordo cella
         Point3d ptVert = m_mTree[nId].GetTopRight() ;
         plInLoop.AddUPoint( 0, ptVert) ;
         ptVert.x = m_mTree[nId].GetBottomLeft().x ;
         ptVert.y = m_mTree[nId].GetTopRight().y ;
         plInLoop.AddUPoint( 1, ptVert) ;
         ptVert = m_mTree[nId].GetBottomLeft() ;
         plInLoop.AddUPoint( 2, ptVert) ;
         ptVert.x = m_mTree[nId].GetTopRight().x ;
         ptVert.y = m_mTree[nId].GetBottomLeft().y ;
         plInLoop.AddUPoint( 3, ptVert) ;
         plInLoop.Close();
         vCellPolygons.push_back( plInLoop) ;
         vPolygons.push_back( vCellPolygons) ;
         ++ nPoly ;
         vnParentChunk.push_back( inA.nChunk) ;
         vCellPolygons.clear() ;
         plInLoop.Clear() ;
      }
      for ( int i = 0 ; i < nInters ; ++ i) {
         inA = m_mTree[nId].m_vInters[i] ;
         if ( inA.nIn == -1) {
            // numero di vertici aggiunti al nuovo poligono
            int k = 0 ;
            for ( Point3d ptInt : inA.vpt) {
               plInLoop.AddUPoint( k, ptInt) ;
               ++ k ;
            }
            plInLoop.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() ;
                     vnParentChunk.push_back( inA.nChunk) ;
                     bAdded = true ;
                     break ;
                  }
               }
            }
            if ( ! bAdded) {
               vCellPolygons.push_back( plInLoop) ;
               vPolygons.push_back( vCellPolygons) ;
               ++ nPoly ;
               vnParentChunk.push_back( inA.nChunk) ;
               plInLoop.Clear() ;
               vCellPolygons.clear() ;
            }
            plInLoop.Clear() ;
         }
         else
            break ;
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
Tree::CheckIfBefore( 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( Inters& inA) const
{ 
   // questa funzione è pensata in riferimento al lato 3, quindi nessuno dei due punti può stare su Edge = 3
   // controllo se l'ingresso è prima dell'uscita
   int nEdge1 = inA.nIn ;
   int nEdge2 = inA.nOut ;
   if ( nEdge1 == -1)
      return false ;
   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

   INTVECTOR :: iterator iter1 = find( vEdges.begin(), vEdges.end(), nEdge1) ;
   int nPos1 = std::distance( vEdges.begin(), iter1) ;
   INTVECTOR :: iterator iter2 = find( vEdges.begin(), vEdges.end(), nEdge2) ;
   int nPos2 = std::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, Point3d& ptP1, int& nEdge2, 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
   INTVECTOR :: iterator iter1 = find( vEdges.begin(), vEdges.end(), nEdge1) ;
   int nPos1 = std::distance( vEdges.begin(), iter1) ;
   INTVECTOR :: iterator iter2 = find( vEdges.begin(), vEdges.end(), nEdge2) ;
   int nPos2 = std::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, Point3d& ptP1, Point3d& ptP2) const
{
   // sul lato nEdge controllo se ptP1 viene prima di ptP2.
   // i lati vengono percorsi in senso antiorario
   if ( AreSameEdge(nEdge, 0) ) {
      return ptP1.x > ptP2.x ;
   }
   else if ( AreSameEdge(nEdge, 1) ) {
      return ptP1.y > ptP2.y ;
   }
   else if ( AreSameEdge(nEdge, 2) ) {
      return ptP1.x < ptP2.x ;
   }
   else if ( AreSameEdge(nEdge, 3) ) {
      return ptP1.y < ptP2.y ;
   }
   return false ;
}

//----------------------------------------------------------------------------
bool
Tree::CheckIfAfter( PolyLine& pl, Point3d& ptNextStart, int nEdge) const
{ 
   // controllo se ptNextStart è più avanti rispetto a ptEnd in senso antiorario
   // ATTENZIONE i due punti devono essere sullo stesso lato!

   // se i punti non sono sullo stesso lato ritorno false
   Point3d ptEnd ;
   pl.GetLastPoint( ptEnd) ;
   if ( AreSameEdge(nEdge, 0) ) {
      if ( ptNextStart.y > ptEnd.y - EPS_SMALL  &&  ptNextStart.y < ptEnd.y + EPS_SMALL)
         return ptEnd.x > ptNextStart.x ;
      else
         return false ;
   }
   else if ( AreSameEdge(nEdge, 1) ) {
      if ( ptNextStart.x > ptEnd.x - EPS_SMALL  &&  ptNextStart.x < ptEnd.x + EPS_SMALL)
         return ptEnd.y > ptNextStart.y ;
      else
         return false ;
   }
   else if ( AreSameEdge(nEdge, 2) ) {
      if ( ptNextStart.y > ptEnd.y - EPS_SMALL  &&  ptNextStart.y < ptEnd.y + EPS_SMALL)
         return ptEnd.x < ptNextStart.x ;
      else
         return false ;
   }
   else if ( AreSameEdge(nEdge, 3) ) {
      if ( ptNextStart.x > ptEnd.x - EPS_SMALL  &&  ptNextStart.x < ptEnd.x + EPS_SMALL)
         return ptEnd.y < ptNextStart.y ;
      else
         return false ;
   }
   return false ;
}

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

//----------------------------------------------------------------------------
bool
Tree::AddVertex( int nId, std::vector<PNTVECTOR>& vEdgeVertex, PolyLine& plTrimmedPoly, int& c, Point3d& ptToAdd) const
{
   // se è il primo punto della PolyLine lo aggiungo
   if ( plTrimmedPoly.GetPointNbr() == 0) {
      plTrimmedPoly.AddUPoint( c, ptToAdd) ;
      ++ c ;
      return true ;
   }
   Point3d ptBr = vEdgeVertex[3][0] ;
   Point3d ptTR = vEdgeVertex[0][0] ;
   Point3d ptTl = vEdgeVertex[1][0] ;
   Point3d ptBL = vEdgeVertex[2][0] ;
   Point3d ptLast ;
   plTrimmedPoly.GetLastPoint( ptLast) ;
   // verifico di essere allineato con un lato, sennò aggiungo e basta
   Vector3d vDir;
   if ( ! AreSamePointExact(ptToAdd, ptLast))
      vDir = ptToAdd - ptLast ;
   else
      return true ;
   // se non riesco a normalizzare perché sono troppo vicino ad un vertice allora aggiungo direttamente il vertice
   if ( ! vDir.Normalize()) {
      plTrimmedPoly.EraseLastUPoint() ;
      if ( AreSamePointApprox( ptToAdd, ptBr))
         plTrimmedPoly.AddUPoint( c, ptBr) ;
      else if ( AreSamePointApprox( ptToAdd, ptTR))
         plTrimmedPoly.AddUPoint( c, ptTR) ;
      else if ( AreSamePointApprox( ptToAdd, ptTl))
         plTrimmedPoly.AddUPoint( c, ptTl) ;
      else if ( AreSamePointApprox( ptToAdd, ptBL))
         plTrimmedPoly.AddUPoint( c, ptBL) ;
      ++ c ;
      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 = 0 ; t < (int)vEdgeVertex[0].size() ; ++ t) {
            Point3d ptIntermed = vEdgeVertex[0][t] ;
            if ( ptIntermed.x > ptToAdd.x  &&  ptIntermed.x < ptLast.x) {
               plTrimmedPoly.AddUPoint( c, ptIntermed) ;
               ++ c ;
            }
         }
         plTrimmedPoly.AddUPoint( c, ptToAdd) ;
         ++ 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 = 0 ; t < (int)vEdgeVertex[1].size() ; ++ t) {
            Point3d ptIntermed = vEdgeVertex[1][t] ;
            if ( ptIntermed.y > ptToAdd.y  &&  ptIntermed.y < ptLast.y) {
               plTrimmedPoly.AddUPoint( c, ptIntermed) ;
               ++ c ;
            }
         }
         plTrimmedPoly.AddUPoint( c, ptToAdd) ;
         ++ 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 = 0 ; t < (int)vEdgeVertex[2].size() ; ++ t) {
            Point3d ptIntermed = vEdgeVertex[2][t] ;
            if ( ptIntermed.x < ptToAdd.x  &&  ptIntermed.x > ptLast.x) {
               plTrimmedPoly.AddUPoint( c, ptIntermed) ;
               ++ c ;
            }
         }
         plTrimmedPoly.AddUPoint( c, ptToAdd) ;
         ++ 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 = 0 ; t < (int)vEdgeVertex[3].size() ; ++ t) {
            Point3d ptIntermed = vEdgeVertex[3][t] ;
            if ( ptIntermed.y < ptToAdd.y  &&  ptIntermed.y > ptLast.y) {
               plTrimmedPoly.AddUPoint( c, ptIntermed) ;
               ++ c ;
            }
         }
         plTrimmedPoly.AddUPoint( c, ptToAdd) ;
         ++ c ;
      }
      // sono allineato con un lato, ma NON sono su un lato
      // aggiungo e basta
      else {
         plTrimmedPoly.AddUPoint( c, ptToAdd);
         ++ c ;
      }
   }
   // non su un edge, quindi aggiungo e basta
   else {
      plTrimmedPoly.AddUPoint( c, ptToAdd);
      ++ c ;
   }
   return true ;
}

//usando i poligoni
////----------------------------------------------------------------------------
//bool
//Tree::SetRightEdgeIn( int nId, std::vector<PNTVECTOR>& vEdgeVertex, PolyLine& plTrimmedPoly)
//{
//   Point3d ptBr( m_mTree.at(nId).GetTopRight().x , m_mTree.at(nId).GetBottomLeft().y) ;
//   Point3d ptTR = m_mTree.at(nId).GetTopRight() ;
//   int nPos0 = -1 , nPos1 = -1 ;
//   bool bSave = false ;
//   PNTLIST lPtInt ;
//   for ( int k = 0 ; k < plTrimmedPoly.GetPointNbr() ; ++ k) {
//      Point3d ptIter ;
//      plTrimmedPoly. GetNextPoint( ptIter) ;
//      if ( AreSamePointExact( ptIter, ptBr) ) {
//         nPos0 = k ;
//         bSave = bSave ^ true ;
//      }
//      if ( AreSamePointExact( ptIter, ptTR) ) {
//         nPos1 = k ;
//         bSave = bSave ^ true ;
//      }
//      if ( bSave) {
//         lPtInt.push_back( ptIter) ;
//      }
//   }
//   lPtInt.pop_front() ;
//   // verifico se ci sono tutti e due i punti
//   if ( nPos0 != -1 && nPos1 != -1 ) {
//      // se sono di seguito tutto il lato destro è dentro
//      if ( nPos1 == nPos0 + 1  ||  nPos0 == nPos1 + 1) {
//         m_mTree[nId].m_nRightEdgeIn = 1 ;
//      }
//      // ci sono dei punti in mezzo, devo capire se sono allineati o no
//      else {
//         // verifico se tutti i punti intermedi sono allineati o no
//         CurveLine clRightEdge ;
//         clRightEdge.Set( ptBr, ptTR) ;
//         bool bAllOn = true ;
//
//         PNTLIST::iterator lIter = lPtInt.begin() ;
//         std::advance(lIter, 0) ;
//         for ( int p = 0 ; p < (int)lPtInt.size() ; ++ p) {
//            if ( ! clRightEdge.IsPointOn( *lIter)) {
//               bAllOn = false ;
//               break ;
//            }
//         }
//      // divisi da punti allineati e quindi lato dentro
//         if ( bAllOn) {
//            m_mTree[nId].m_nRightEdgeIn = 1 ;
//         }
//      // se sono dentro, ma non allineati, allora solo in parte
//         else {
//            m_mTree[nId].m_nRightEdgeIn = 2 ;
//         }
//      }
//   }
//   // il lato destro è sicuramente esterno
//   else {
//      m_mTree[nId].m_nRightEdgeIn = 0 ;
//   }
//   return true ;
//}

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

   if ( m_mTree[nId].m_nRightEdgeIn != -1) 
      return true ;
   else
      return false ;
}

//----------------------------------------------------------------------------
bool
Tree::CategorizeCell( int& nId)
{
   if ( m_mTree[nId].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 ( (int)vNeigh.size() == 0 ) {
      m_mTree[nId].m_nFlag = 0 ;
   }
   else if ( (int)vNeigh.size() == 1) {
      if ( m_mTree[vNeigh[0]].m_nRightEdgeIn == 1)
         m_mTree[nId].m_nFlag = 4 ;
      else if ( m_mTree[vNeigh[0]].m_nRightEdgeIn == 0)
         // devo verificare se la cella è intersecata
         if ( m_mTree[vNeigh[0]].m_nFlag == 1  ||  m_mTree[vNeigh[0]].m_nFlag == 3)
            m_mTree[nId].m_nFlag = 0 ;
         else {
            if ( m_mTree[vNeigh[0]].m_nFlag == 4) {
               m_mTree[nId].m_nFlag = 4 ;
            }
            else if ( m_mTree[vNeigh[0]].m_nFlag == 0) {
               m_mTree[nId].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 ( (int) vNeigh1.size() > 0) {
            int nNeigh = vNeigh1[0] ;
            if ( m_mTree[nNeigh].m_nFlag == 0  ||  m_mTree[nNeigh].m_nFlag == 2) {
               m_mTree[nId].m_nFlag = 0 ;
            }
            else if ( m_mTree[nNeigh].m_nFlag == 4) {
               m_mTree[nId].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[nId].m_nFlag = 0 ;
               else if ( m_mTree[nNeigh].m_nRightEdgeIn == 1)
                  m_mTree[nId].m_nFlag = 4 ;
               else if ( m_mTree[nNeigh].m_nRightEdgeIn == 2 ) {
                  // in questo caso devo verificare le intersezioni sul lato destro del vicino bottom
                  // scorro e cerco l'intersezione più alta, se è out la cella è dentro, se è in devo verificare l'out se è più a sinistra o a destra.
                  int nPass = (int) m_mTree[nNeigh].m_vInters.size() ;
                  bool bTopMost = false ;
                  bool bFound = false;
                  Point3d ptInters, ptBr( m_mTree[nNeigh].GetTopRight().x, m_mTree[nNeigh].GetBottomLeft().y) ;
                  ptInters = ptBr ;  // inizializzato a vertice bottom right
                  int nEdgeIn = 3, nEdgeOut = 3, nLoop ;
                  for ( int r = 0 ; r < nPass; ++ r ) {
                     // trovo il loop che ha l'ingresso o l'uscita più in alto sul lato destro
                     // verifico che o l'ingresso o l'uscita siano sul lato destro e che sia l'intersezione più alta su quel lato
                     if ( m_mTree[nNeigh].m_vInters[r].nIn == 3  &&  ! CheckIfBefore( m_mTree[nNeigh].m_vInters[r].nIn, m_mTree[nNeigh].m_vInters[r].vpt[0], ptInters)) {
                        nLoop = r ;
                        ptInters = m_mTree[nNeigh].m_vInters[r].vpt[0] ;
                        bFound = true ;
                     }
                     if ( m_mTree[nNeigh].m_vInters[r].nOut == 3  &&  ! CheckIfBefore( m_mTree[nNeigh].m_vInters[r].nOut, m_mTree[nNeigh].m_vInters[r].vpt.back(), ptInters)) {
                        nLoop = r ;
                        ptInters = m_mTree[nNeigh].m_vInters[r].vpt.back() ;
                        bFound = true ;
                     }
                     // controllo nei vertici
                     if ( m_mTree[nNeigh].m_vInters[r].nOut == 7  ||  m_mTree[nNeigh].m_vInters[r].nIn == 7) {
                        nLoop = r ;
                        ptInters = m_mTree[nNeigh].GetTopRight() ;
                        bFound = true ;
                        break ;
                     }
                     if ( m_mTree[nNeigh].m_vInters[r].nOut == 6 || m_mTree[nNeigh].m_vInters[r].nIn == 6 ) {
                        nLoop = r ;
                        ptInters = ptBr ;
                        bFound = true ;
                     }
                  }
                  if ( bFound  &&  CheckIfBefore( m_mTree[nNeigh].m_vInters[nLoop])) {
                     bTopMost = true ;
                  }
                  if ( bFound  &&  bTopMost) {
                     if ( AreSameEdge(m_mTree[nNeigh].m_vInters[nLoop].nIn, 3) && AreSameEdge(m_mTree[nNeigh].m_vInters[nLoop].nOut, 3) ) {
                        if ( CheckIfBefore( m_mTree[nNeigh].m_vInters[nLoop]))
                           m_mTree[nId].m_nFlag = 4 ;
                        else
                           m_mTree[nId].m_nFlag = 0 ;
                     }
                     else if ( AreSameEdge(m_mTree[nNeigh].m_vInters[nLoop].nOut, 3))
                        m_mTree[nId].m_nFlag = 4 ;
                     else if ( AreSameEdge(m_mTree[nNeigh].m_vInters[nLoop].nIn, 3) ) {
                        // devo verificare se l'uscita è più a sinistra o più a destra della cella
                        Point3d ptOut = m_mTree[nNeigh].m_vInters[nLoop].vpt.back() ;
                        Point3d ptBr( m_mTree[nId].GetTopRight().x, m_mTree[nId].GetBottomLeft().y) ;
                        if ( ptOut.x >= ptBr.x)
                           m_mTree[nId].m_nFlag = 4 ;
                        else
                           m_mTree[nId].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 nLoop ;
            int nEdge = 3 ;
            bool bIntersIn ;
            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)) {
                     ptInters = m_mTree[vNeigh[0]].m_vInters[k].vpt[0] ;
                     nLoop = k ;
                     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) ) {
                     ptInters = m_mTree[vNeigh[0]].m_vInters[k].vpt.back() ;
                     nLoop = k ;
                     bIntersIn = false ;
                  }
               }
            }
            if ( bIntersIn)
               m_mTree[nId].m_nFlag = 4 ;
            else
               m_mTree[nId].m_nFlag = 0 ;
         }
      }
   }
   else {
      for ( int r = 0 ; r < (int)vNeigh.size() ; ++ r) {
         if ( m_mTree[vNeigh[r]].m_nRightEdgeIn == 1) {
            m_mTree[nId].m_nFlag = 4 ;
            break ;
         }
         else if ( m_mTree[vNeigh[r]].m_nRightEdgeIn == 0) {
            if ( m_mTree[vNeigh[r]].m_nFlag == 4) {
               m_mTree[nId].m_nFlag = 4 ;
               break ;
            }
            else if ( m_mTree[vNeigh[r]].m_nFlag == 0) {
               m_mTree[nId].m_nFlag = 0 ;
               break ;
            }
         }
      }
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
Tree::CheckIfBetween( Inters& inA, 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 ;
   while ( ! AreSameEdge( nEdge, inA.nIn) || (int) vEdges.size() == 0) {
      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)) {
            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()) ) {
               if ( CheckIfBefore( nEdge, inA.vpt.back(), inB.vpt[0])  ||  CheckIfBefore(  nEdge, inB.vpt[0], inA.vpt[0]))
                  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])  &&  CheckIfBefore(  nEdge, inB.vpt[0], inA.vpt[0]))
                  bFound = true ;
            }
         }
         else if ( AreSameEdge(inB.nIn, inA.nOut)) {
            PolyLine pl ;
            pl.AddUPoint( 0, inA.vpt[0]) ;
            pl.AddUPoint( 1, inA.vpt.back()) ;
            if ( CheckIfAfter( pl, inB.vpt[0], i))
               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]))
               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())) {
               bFound = true ;
            }
            else if ( ! AreSameEdge(inB.nOut,inB.nIn))
               bFound = true ;
      }
      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::GetLeaves( std::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, 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.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 ;
   if ( AreSamePointApprox( ptToAssign, ptTR))
      nEdge = 7 ;
   else if ( AreSamePointApprox( ptToAssign, ptTl))
      nEdge = 4 ;
   else if ( AreSamePointApprox( ptToAssign, ptBL))
      nEdge = 5 ;
   else if ( AreSamePointApprox( ptToAssign, ptBr))
      nEdge = 6 ;
   else
      return false ;
   return true ;
}