//----------------------------------------------------------------------------
//  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 "/EgtDev/Include/EGkPolyLine.h"
#include "/EgtDev/Include/EGkDistPointCurve.h"

using namespace std ;

//----------------------------------------------------------------------------
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_ptPbl( ORIG), m_ptPtr(), m_bProcessed ( false) , m_bSplitVert ( true)
{
   Point3d ptTr ( 1, 1) ;
   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_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)
{
   Point3d ptBl( 0, 0), ptTr ( 1, 1) ;
   Cell cRoot( ptBl, ptTr) ;
   m_mTree.insert( pair< int, Cell>( -1, cRoot)) ;
}

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

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

//----------------------------------------------------------------------------
void Tree::SetSurf( const SurfBezier* pSrfBz, bool bSplitPatches)
{
   m_pSrfBz = pSrfBz ;
  // le coordinate delle celle sono nello spazio parametrico
   int nDegU, nDegV, nSpanU, nSpanV ;
   bool bIsRat, bTrimmed ;
   m_pSrfBz->GetInfo( nDegU, nDegV, nSpanU, nSpanV, bIsRat, bTrimmed) ;
   m_bTrimmed = bTrimmed ;
   m_nDegU = nDegU ;
   m_nDegV = nDegV ;
   if ( nDegU == 1  &&  nDegV == 1) 
      m_bBilinear = true ;
   if ( nSpanU * nSpanV != 1)
      m_bMulti = true ;
   
  // salvo i vertici 3d della cella root
   Point3d ptBottom ( 0, 0) ;
   Point3d ptTop( nSpanU, nSpanV) ;
   Cell cRoot( ptBottom, ptTop) ;
   m_mTree.insert( 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( 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) ;
         ++ 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) ;
            nId = m_mTree[nId].m_nChild2 ;
         }
      }
      // split preliminari per dividere le patch in modo da triangolarle indipendentemente////////////////////////////////////////////////////////
   }
  // 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( pair<int, Cell>( nNodes - 1, cChild1)) ;
      m_mTree.insert( pair<int, Cell>( nNodes, cChild2)) ;
      Point3d ptVert1, ptVert2 ;
      PNTVECTOR vVert ;
      m_mVert.insert( pair<int, PNTVECTOR>( nNodes - 1, vVert)) ;
      m_mVert.insert( 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, dSplitValue, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptVert1) ;
         m_pSrfBz->GetPointD1D2( m_mTree[nId].GetTopRight().x, dSplitValue, 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, m_mTree[nId].GetBottomLeft().y, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptVert2) ;
         m_pSrfBz->GetPointD1D2( dSplitValue, m_mTree[nId].GetTopRight().y, 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( 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) {
      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 ;
               double dLenParV = m_mTree[nCToSplit].GetTopRight().y - m_mTree[nCToSplit].GetBottomLeft().y ;
               if ( dLenParU <= 1. / m_nDegV  ||  dLenParV <= 1. / m_nDegU) {
                  double dU = ( m_mTree[nCToSplit].GetTopRight().x + m_mTree[nCToSplit].GetBottomLeft().x) / 2 ;
                  double dV = ( m_mTree[nCToSplit].GetTopRight().y + m_mTree[nCToSplit].GetBottomLeft().y) / 2 ;
                  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, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptP00P10) ;
                  m_pSrfBz->GetPointD1D2( m_mTree[nCToSplit].GetTopRight().x, dV, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptP10P11) ;
                  m_pSrfBz->GetPointD1D2( dU, m_mTree[nCToSplit].GetTopRight().y, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptP11P01) ;
                  m_pSrfBz->GetPointD1D2( m_mTree[nCToSplit].GetBottomLeft().x, 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 ;
                     double dV = ( m_mTree[nCToSplit].GetTopRight().y + m_mTree[nCToSplit].GetBottomLeft().y) / 2 ;
                     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, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptP00P10) ;
                     m_pSrfBz->GetPointD1D2( dU, m_mTree[nCToSplit].GetTopRight().y, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptP11P01) ;
                     CurveLine clV ;
                     clV.Set( ptP00P10, ptP11P01) ;
                     DistPointCurve dpc( ptPSrf, clV) ;
                     double dDist ;
                     dpc.GetDist( dDist) ;
                     dCurvV = 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 ;
                     double dV = k * m_mTree[nCToSplit].GetTopRight().y + ( 1 - k) * m_mTree[nCToSplit].GetBottomLeft().y ;
                     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, dV, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptP10P11) ;
                     m_pSrfBz->GetPointD1D2( m_mTree[nCToSplit].GetBottomLeft().x, dV, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptP01P00) ;
                     CurveLine clU ;
                     clU.Set( ptP01P00, ptP10P11) ;
                     DistPointCurve dpc( ptPSrf, clU) ;
                     double dDist ;
                     dpc.GetDist( dDist) ;
                     dCurvU = 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 = min( dLen0, dLen2) ;
                  else
                     dSideMinVal = max( dLen0, dLen2) ;
               }
               else {
                  if ( dLen1 != 0  &&  dLen3 != 0)
                     dSideMinVal = min( dLen1, dLen3) ;
                  else
                     dSideMinVal = max( dLen1, dLen3) ;
               }
              // calcolo le diagonali per controllare la dimensione massima dei triangoli in cui dividerei la cella
               dSideMaxVal = 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 ;
                     if ( ! m_pSrfBz->GetPointD1D2( dULoc, m_mTree[nCToSplit].GetBottomLeft().y, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptBz0)  ||
                          ! m_pSrfBz->GetPointD1D2( dULoc, m_mTree[nCToSplit].GetTopRight().y, 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 ;
                        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].Processed() ;
                 // risalgo i parent finché non trovo il primo Child2 da processare
                  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].Processed() ;
                  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].Processed() ;
                     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 ( max(dLen0, dLen2) > 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 ;
                     dV = ( 1 - i) * m_mTree[nCToSplit].GetBottomLeft().y + i * m_mTree[nCToSplit].GetTopRight().y ;
                     double dVLoc = ( m_mTree[nCToSplit].GetBottomLeft().y + m_mTree[nCToSplit].GetTopRight().y) / 2 ;
                     double dULoc = ( m_mTree[nCToSplit].GetBottomLeft().x + m_mTree[nCToSplit].GetTopRight().x) / 2 ;
                     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 = min( dLen0, dLen2) ;
                  else
                     dSideMinVal = max( dLen0, dLen2) ;
               }
               else {
                  if ( dLen1 != 0  &&  dLen3 != 0)
                     dSideMinVal = min( dLen1, dLen3) ;
                  else
                     dSideMinVal = max( dLen1, dLen3) ;
               }
              // calcolo le diagonali per controllare la dimensione massima dei triangoli in cui dividerei la cella
               dSideMaxVal = max( Dist( ptP00, ptP11), Dist( ptP10, ptP01)) ;

           
               double dErr = 0 ;
               if ( m_bMulti) {
                  Point3d ptPSrf ;
                  Plane3d plAppr ;
                  if ( ! AreSamePointApprox( ptP00, ptP10)  &&  ! AreSamePointApprox( ptP00, ptP01))
                     plAppr.Set( ptP00, ( ptP00 - ptP01) ^ ( ptP00 - ptP10)) ;
                  else if ( AreSamePointApprox( ptP00, ptP10)) {
                     plAppr.Set( ptP01, ( ptP00 - ptP01) ^ ( ptP01 - ptP11)) ;
                  }
                  else if ( AreSamePointApprox( ptP00, ptP01)) {
                     plAppr.Set( ptP10, ( ptP10 - ptP11) ^ ( ptP00 - ptP10)) ;
                  }
                  for ( double i = 0.25 ; i < 1 ; i = i + 0.25) {
                     for ( double j = 0.25 ; j < 1 ; j = j + 0.25) {
                        double dU =  ( 1 - i) * m_mTree[nCToSplit].GetTopRight().x + i * m_mTree[nCToSplit].GetBottomLeft().x ;
                        double dV =  ( 1 - j) * m_mTree[nCToSplit].GetTopRight().y + j * m_mTree[nCToSplit].GetBottomLeft().y ;
                        m_pSrfBz->GetPointD1D2( dU, dV, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptPSrf) ;
                        dErr = 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].Processed() ;
                  // risalgo i parent finché non trovo il primo Child2 da processare
                  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].Processed() ;
                  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].Processed() ;
                     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 ;
            }
         }
      }
   }
  // se la superficie è trimmata
   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.
   //}
}

//----------------------------------------------------------------------------
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) ;
            }
         }
      }
   }
   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) ;
            }
         }
      }
   }
   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) ;
            }
         }
      }
   }
   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) ;
            }
         }
      }
   }
   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) ;
}

//----------------------------------------------------------------------------
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 = 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 ;
}

//----------------------------------------------------------------------------
bool Tree::GetPolygons( POLYLINEVECTOR& vPolygons)
{
   if ( m_vPolygons.empty()) {
      PNTVECTOR vVertices ;
      INTVECTOR vNeigh ;
      bool bBottomRight , bTopLeft ;
      // scorro lungo tutte le celle leaves ( dell'albero bilanciato) 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()) ;
         // 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 ;
            dV = ( m_mTree.at(nId).GetBottomLeft().y + m_mTree.at(nId).GetTopRight().y) / 2 ;
            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() ;
         for ( int i = 0 ; i < (int) vVertices.size() ; ++i) {
            m_vPolygons.back().AddUPoint( i, vVertices.at(i)) ;
         }
      }
   }

   // restituisco i poligoni delle celle del tree nello spazio parametrico
   vPolygons = m_vPolygons ;
   return true ;
}