EgtGeomKernel :

- piccoli aggiustamenti.

CAvSilhouetteSurfTm.cpp

@@ -1019,7 +1019,7 @@ CAvParSilhouettesSurfTm::Prepare( void)
    m_dDimZ = b3All.GetDimZ() ;
 
   // set utensile corrente per calcolo delle collisioni
-   double dExtraXY = m_dRad ;
+   double dExtraXY = m_dRad + m_dOffsR ;
    if ( m_dSideAng < EPS_ANG_SMALL)
       m_cavTstm.SetStdTool( m_dDimZ + m_dOffsR, m_dRad + m_dOffsR, m_dCornRad + m_dOffsR) ;
    else {
@@ -1037,7 +1037,7 @@ CAvParSilhouettesSurfTm::Prepare( void)
 
       double dStemRad = m_dRad + dDeltaRad ;
       double dTipRad = m_dRad ;
-      dExtraXY = dStemRad ;
+      dExtraXY = dStemRad + m_dOffsR ;
       m_cavTstm.SetAdvTool( m_dDimZ + m_dOffsR, dStemRad + m_dOffsR, m_dMaxMat, dTipRad + m_dOffsR, m_dCornRad + m_dOffsR) ;
    }
 

CalcPocketing.cpp

@@ -1045,14 +1045,12 @@ ExtendPath( ICurveComposite* pCompo, const ISurfFlatRegion* pSfr, const PocketPa
       vtTan.Invert() ;
    }
 
-
   // angoli di rotazione ( potrebbero servirne altri)
    double dAng45 = ANG_RIGHT / 2. ;
    static DBLVECTOR vAngles = { 0., ANG_RIGHT, 3 * ANG_RIGHT,
                                 dAng45, 3 * dAng45, 5 * dAng45, 7 * dAng45 } ;
 
   // ruoto il versore di uscita cercando un'entrata valida
-   
    double dMinDist = PockParams.dRad + PockParams.dRadialOffset ;
    for ( int i = 0 ; i < int( vAngles.size()) ; ++ i) {
      // ruoto il versore d'uscita
@@ -2420,6 +2418,13 @@ CalcCircleSpiral( const Point3d& ptCen, const Vector3d& vtN, double dOutRad, dou
    PtrOwner<ICurveArc> pArc( CreateCurveArc()) ;
    if ( IsNull( pArc)  ||  ! pArc->Set( ptCen, vtN, dOutRad))
       return false ;
+  // cambio il punto inizile se necessario
+   if ( PockParams.ptStart.IsValid()) {
+      double dUStart = 0. ;
+      int nFlag = 0 ;
+      DistPointCurve( PockParams.ptStart, *pArc).GetParamAtMinDistPoint( 0., dUStart, nFlag) ;
+      pArc->ChangeStartPoint( dUStart) ;
+   }
 
   // creo la superificie per la Feed
    PtrOwner<ISurfFlatRegion> pSrfRemoved( CreateSurfFlatRegion()) ;
@@ -3231,10 +3236,12 @@ CalcTrapezoidSpiralXCoord( const ICurveComposite* pCrvPocket, int nBase, int nSe
                for ( int i = 0 ; i < intCC2.GetIntersCount() ; ++ i) {
                   IntCrvCrvInfo ccClass2 ;
                   if ( intCC2.GetIntCrvCrvInfo( i, ccClass2)) {
+                     double dDiffY = ccClass2.IciA[0].ptI.y - dYCoord ;
+                     double dOffs = sqrt( dRad * dRad - dDiffY * dDiffY) ;
                      if ( bStart)
-                        dXCoord = max( dXCoord, Box3d.GetMax().x + sqrt( dRad * dRad - pow(( ccClass2.IciA[0].ptI.y - dYCoord), 2))) ;
+                        dXCoord = max( dXCoord, Box3d.GetMax().x + dOffs) ;
                      else
-                        dXCoord = min( dXCoord, Box3d.GetMin().x - sqrt( dRad * dRad - pow(( ccClass2.IciA[0].ptI.y - dYCoord), 2))) ;
+                        dXCoord = min( dXCoord, Box3d.GetMin().x - dOffs) ;
                   }
                }
             }
@@ -6338,7 +6345,9 @@ CalcZigZag( const ISurfFlatRegion* pSrfZigZag, const PocketParams& PockParams, I
       CRVCVECTOR ccClass ;
       pSrfZigZag->GetCurveClassification( *pLine, EPS_SMALL, ccClass) ;
       for ( int j = 0 ; j < int( ccClass.size()) ; ++ j) {
-         if ( ccClass[j].nClass == CRVC_IN) { // memorizzo il segmento
+        // se interno
+         if ( ccClass[j].nClass == CRVC_IN) {
+           // recupero il tratto
             Section currSec ;
             currSec.bActive = true ;
             PtrOwner<ICurveLine> pSeg( GetCurveLine( pLine->CopyParamRange( ccClass[j].dParS, ccClass[j].dParE))) ;
@@ -6346,6 +6355,9 @@ CalcZigZag( const ISurfFlatRegion* pSrfZigZag, const PocketParams& PockParams, I
                continue ;
             Point3d ptS, ptE ;
             pSeg->GetStartPoint( ptS) ; pSeg->GetEndPoint( ptE) ;
+           // se troppo piccolo, lo scarto
+            if ( SqDist( ptS, ptE) < 250 * 250 * SQ_EPS_SMALL)
+               continue ;
             currSec.ptS = ptS ; currSec.ptE = ptE ;
             for ( int k = 0 ; k < int( vFirstOff.size()) ; ++ k) {
                if ( vFirstOff[k]->IsPointOn( ptS)) {
@@ -7361,10 +7373,10 @@ CalcConformalLink( ICurveComposite* pCrvOffs0, ICurveComposite* pCrvOffs1, Pocke
          vCrvChunk.emplace_back( vCrvClassBorder[vIndOffs0[i]]->Clone()) ;
      // se la curva successiva è chiusa, aggiorno il suo punto iniziale a quello a minima distanza
       if ( ! bOpen1) {
-         Point3d ptS ; pCrvOffs0->GetStartPoint( ptS) ;
+         Point3d ptE ; pCrvOffs0->GetEndPoint( ptE) ;
          double dPar ;
          int nFlag ;
-         if ( DistPointCurve( ptS, *pCrvOffs1).GetParamAtMinDistPoint( 0., dPar, nFlag))
+         if ( DistPointCurve( ptE, *pCrvOffs1).GetParamAtMinDistPoint( 0., dPar, nFlag))
             pCrvOffs1->ChangeStartPoint( dPar) ;
       }
       if ( ! CutCurveToConnect( pCrvOffs0, pCrvOffs1, vCrvChunk, PockParams,
@@ -7810,13 +7822,27 @@ CalcConformalOffsAndLinks( VICRVCOMPOPOVECTOR& vvCrvOffs, const ISurfFlatRegion*
      I link che collegano due curve di Offset aperte devono partire dalla prima e, seguendo il
      bordo della regione di classificazione ( pSfrClass) arrivare alla seconda.
   */
+  /* NB.
+     La regione di classificazione ( pSfrClass) deve tenere conto della regione limite; nel caso di
+     bordi chiusi, bisogna stare a distanza opportuna da essi, non raccordandoci i Links
+  */
   /* NB.
      I link che collegano due curve di Offset chiuse tra di loro vengono calcolati esattamente
      come per i percorsi SPIRAL, pertanto devo definire un insieme di curve di primo Offset sui
      quali smussare i Link per evitare di uscire dalla regione di svuotatura
   */
    ICRVCOMPOPOVECTOR vCrvSfrClass ;
-   if ( ! GetSfrCrvCompoLoops( pSfrClass, vCrvSfrClass))
+   PtrOwner<ISurfFlatRegion> pSfrClassReal( CloneSurfFlatRegion( pSfrClass)) ;
+   if ( IsNull( pSfrClassReal)  ||  ! pSfrClassReal->IsValid())
+      return false ;
+   if ( PockParams.SfrLimit.IsValid()) {
+      PtrOwner<ISurfFlatRegion> pSfrLimit( CloneSurfFlatRegion( &PockParams.SfrLimit)) ;
+      if ( IsNull( pSfrLimit)  ||  ! pSfrLimit->IsValid()  ||
+           ! pSfrLimit->Offset( PockParams.dRad + PockParams.dRadialOffset, ICurve::OFF_FILLET)  ||
+           ! pSfrClassReal->Subtract( *pSfrLimit))
+         return false ;
+   }
+   if ( ! GetSfrCrvCompoLoops( pSfrClassReal, vCrvSfrClass))
         return false ;
 
    /* NB.
@@ -8049,8 +8075,12 @@ OrderAndExtendConformalPaths( ICRVCOMPOPOVECTOR& vCrvPaths, const ISurfFlatRegio
   // superficie di controllo per parti isolate
    PtrOwner<ISurfFlatRegion> pSfrToRemove( CloneSurfFlatRegion( pSfrChunk)) ;
    if ( IsNull( pSfrToRemove)  ||  ! pSfrToRemove->IsValid()  ||
-        ! pSfrToRemove->Intersect( *pSfrOrig))
-      return false ;
+        ! pSfrToRemove->Intersect( *pSfrOrig)) {
+     // riprovo con una leggera tolleranza
+      pSfrToRemove->Offset( - 10 * EPS_SMALL, ICurve::OFF_FILLET) ;
+      if ( ! pSfrToRemove->Intersect( *pSfrOrig))
+         return false ;
+   }
 
   // creo le regioni piane di svuotatura dei percorsi
    ISURFFRPOVECTOR vSfrRemoved ; vSfrRemoved.resize( vCrvPaths.size()) ;
@@ -8586,7 +8616,7 @@ GetConformalOffsets( const ISurfFlatRegion* pSfrChunk, const ISurfFlatRegion* pS
    if ( PockParam.SfrLimit.IsValid()) {
       if ( ! pSfrLimit.Set( PockParam.SfrLimit.Clone())  ||
            ! pSfrLimit->IsValid()  ||
-           ! pSfrLimit->Offset( PockParam.dRad + PockParam.dRadialOffset + 5 * EPS_SMALL, ICurve::OFF_FILLET))
+           ! pSfrLimit->Offset( PockParam.dRad + PockParam.dRadialOffset - 5 * EPS_SMALL, ICurve::OFF_FILLET))
          return false ;
    }
 
@@ -8644,8 +8674,8 @@ GetConformalOffsets( const ISurfFlatRegion* pSfrChunk, const ISurfFlatRegion* pS
                      }
                   }
                   else {
-                     // almeno una curva in trovata
-                      bStop = false ;
+                    // almeno una curva in trovata
+                     bStop = false ;
                     // inverto se necessario
                      if ( bInvert)
                         pMyCrv->Invert() ;

CurveAux.cpp

@@ -1098,9 +1098,9 @@ ApproxArcCurveBezierWithSingleCubic( const ICurve* pCrv, const Point3d& ptCen)
    double ay = ptStart.y - ptCen.y ;
    double bx = ptEnd.x - ptCen.x ;
    double by = ptEnd.y - ptCen.y ;
-   double q1 = pow( ax, 2) + pow( ay, 2) ;
+   double q1 = ax * ax + ay * ay ;
    double q2 = q1 + ax * bx + ay * by ;
-   double k2 = 4./3 * ( sqrt( 2 * q1 * q2) - q2) / ( ax * by  - ay * bx) ;
+   double k2 = 4. / 3 * ( sqrt( 2 * q1 * q2) - q2) / ( ax * by  - ay * bx) ;
    Point3d ptCtrl1, ptCtrl2 ;
    ptCtrl1.x = ptStart.x - k2 * ay ;
    ptCtrl1.y = ptStart.y + k2 * ax ;

CurveBezier.cpp

@@ -2274,17 +2274,17 @@ CurveBezier::MakeRationalStandardForm( void)
 { 
    if ( ! m_bRat)
       return false ;
-   double dW0 = m_vWeCtrl[0] ;
+   double dW0 = m_vWeCtrl.front() ;
    double dWn = m_vWeCtrl.back() ;
-   if( dW0 > 1- EPS_ZERO  &&  dWn > 1 - EPS_ZERO)
+   if ( dW0 > 1 - EPS_ZERO  &&  dWn > 1 - EPS_ZERO)
       return true ;
-   if( dW0 < EPS_ZERO  ||  dWn < EPS_ZERO)
+   if ( dW0 < EPS_ZERO  ||  dWn < EPS_ZERO)
       return false ;
 
-   // formula del Farin
+  // formula del Farin
    double dCoeff = pow( dW0 / dWn, 1. / m_nDeg) ;
    for ( int i = 0 ; i < m_nDeg + 1 ; ++i)
-      m_vWeCtrl[i] *= pow( dCoeff, i) / dW0 ;
+      m_vWeCtrl[i] *= Pow( dCoeff, i) / dW0 ;
 
    return true ;
 }
@@ -2293,13 +2293,13 @@ CurveBezier::MakeRationalStandardForm( void)
 bool
 CurveBezier::MakeNonRational( double dTol)
 {
-   if( ! m_bRat)
+   if ( ! m_bRat)
       return true ;
 
    // controllo se i pesi sono tutti == 1 allora è una finta razionale e mi basta fare una copia dei punti di controllo
    bool bIsActualRat = false ;
    for ( int i = 0 ; i < m_nDeg ; ++i) {
-      if ( abs(m_vWeCtrl[i] - 1) > EPS_SMALL) {
+      if ( abs( m_vWeCtrl[i] - 1) > EPS_SMALL) {
          bIsActualRat = true ;
          break ;
       }
@@ -2323,8 +2323,8 @@ CurveBezier::MakeNonRational( double dTol)
          pNewBez->Init( nDeg, false) ;
          PNTVECTOR vPntCtrl ;
          PNTVECTOR vPntSampling ;
-         for ( int p = 0 ; p < nDeg + 1; ++p) {
-            Point3d pt ; GetPointD1D2( double(p) / nDeg, pt) ;
+         for ( int p = 0 ; p < nDeg + 1 ; ++p) {
+            Point3d pt ; GetPointD1D2( double( p) / nDeg, pt) ;
             pNewBez->SetControlPoint( p, pt) ;
             vPntCtrl.push_back( pt) ;
          }
@@ -2334,8 +2334,8 @@ CurveBezier::MakeNonRational( double dTol)
          while ( dErr > dTol  &&  c < 100) {
             double dErrMax = 0 ;
             // calcolo le differenze tra i punti di sampling sulla nuova curva e quelli sulla curva originale
-            for ( int p = 0 ; p < nDeg + 1; ++p) {
-               Point3d pt ; pNewBez->GetPointD1D2( double(p) / nDeg, pt) ;
+            for ( int p = 0 ; p < nDeg + 1 ; ++p) {
+               Point3d pt ; pNewBez->GetPointD1D2( double( p) / nDeg, pt) ;
                Vector3d vDiff = vPntSampling[p] - pt ;
                double dErrLoc = vDiff.Len() ;
                if( dErrLoc > dErrMax)
@@ -2353,15 +2353,15 @@ CurveBezier::MakeNonRational( double dTol)
          // calcolo l'errore di approssimazione sulla curva
          CalcBezierApproxError( this, pNewBez, dErr) ;
          bOk = dErr < dTol ;
-         if( bOk) {
+         if ( bOk) {
             // aggiorno la curva di bezier originale con quella approssimata
             Init( nDeg, false) ;
-            for( int i = 0 ; i < nDeg + 1 ; ++i) {
+            for ( int i = 0 ; i < nDeg + 1 ; ++i) {
                SetControlPoint( i, pNewBez->GetControlPoint( i)) ;
                SetControlWeight( i, pNewBez->GetControlWeight( i)) ;
             }
          }
-         else if( nDeg < m_nDeg + 4)
+         else if ( nDeg < m_nDeg + 4)
             goto retry ;
    }
 

CurveComposite.cpp

@@ -381,9 +381,10 @@ CurveComposite::FromPolyLine( const PolyLine& PL)
       return false ;
 
   // ciclo di inserimento dei segmenti che uniscono i punti
+   double dParIni, dParFin ;
    Point3d ptIni, ptFin ;
-   PL.GetFirstPoint( ptIni) ;
-   while ( PL.GetNextPoint( ptFin)) {
+   PL.GetFirstUPoint( &dParIni, &ptIni) ;
+   while ( PL.GetNextUPoint( &dParFin, &ptFin)) {
      // se i punti della coppia coincidono, passo alla coppia successiva
       if ( AreSamePointApprox( ptIni, ptFin))
          continue ;
@@ -394,10 +395,14 @@ CurveComposite::FromPolyLine( const PolyLine& PL)
      // assegno i punti estremi
       if ( ! pCrvLine->Set( ptIni, ptFin))
          return false ;
+     // assegno i parametri degli estremi
+      pCrvLine->SetTempParam( dParIni, 0) ;
+      pCrvLine->SetTempParam( dParFin, 1) ;
      // aggiungo la retta alla curva composita
       if ( ! AddSimpleCurve( Release( pCrvLine)))
          return false ;
      // aggiorno dati prossimo punto iniziale
+      dParIni = dParFin ;
       ptIni = ptFin ;
    }
 

IntersLineSurfTm.cpp

@@ -23,20 +23,20 @@ using namespace std ;
 //----------------------------------------------------------------------------
 static void
 UpdateInfoIntersLineSurfTm( const Point3d& ptL, const Vector3d& vtDir, double dLen,
-                            int nT, const Triangle3d& Tria, ILSIVECTOR& vInfo, bool bFinite)
+                            int nStm, int nT, const Triangle3d& Tria, ILSIVECTOR& vInfo, bool bFinite)
 {
    Point3d ptInt, ptInt2 ;
    int nRes = IntersLineTria( ptL, vtDir, dLen, Tria, ptInt, ptInt2, bFinite) ;
    if ( nRes == ILTT_IN  ||  nRes == ILTT_EDGE  ||  nRes == ILTT_VERT) {
       double dU = ( ptInt - ptL) * vtDir ;
       double dCosDN = vtDir * Tria.GetN() ;
-      vInfo.emplace_back( nRes, dU, nT, dCosDN, ptInt) ;
+      vInfo.emplace_back( nRes, dU, nStm, nT, dCosDN, ptInt) ;
    }
    else if ( nRes == ILTT_SEGM  ||  nRes == ILTT_SEGM_ON_EDGE) {
       double dU = ( ptInt - ptL) * vtDir ;
       double dU2 = ( ptInt2 - ptL) * vtDir ;
       double dCosDN = vtDir * Tria.GetN() ;
-      vInfo.emplace_back( nRes, dU, dU2, nT, dCosDN, ptInt, ptInt2) ;
+      vInfo.emplace_back( nRes, dU, dU2, nStm, nT, dCosDN, ptInt, ptInt2) ;
    }
 }
 
@@ -48,7 +48,7 @@ OrderInfoIntersLineSurfTm( ILSIVECTOR& vInfo)
    if ( vInfo.empty())
       return ;
   // ordino il vettore delle intersezioni secondo il senso crescente del parametro di linea
-   sort( vInfo.begin(), vInfo.end(), 
+   sort( vInfo.begin(), vInfo.end(),
             []( const IntLinStmInfo& a, const IntLinStmInfo& b)
                {  double dUa = ( ( a.nILTT == ILTT_SEGM  ||  a.nILTT == ILTT_SEGM_ON_EDGE) ? ( a.dU + a.dU2) / 2 : a.dU) ;
                   double dUb = ( ( b.nILTT == ILTT_SEGM  ||  b.nILTT == ILTT_SEGM_ON_EDGE) ? ( b.dU + b.dU2) / 2 : b.dU) ;
@@ -108,7 +108,7 @@ IntersLineSurfTm( const Point3d& ptL, const Vector3d& vtL, double dLen, const IS
                Triangle3d Tria ;
                Stm.GetTriangle( nT, Tria) ;
               // aggiorno info con intersezione
-               UpdateInfoIntersLineSurfTm( ptL, vtDir, dLen, nT, Tria, vInfo, bFinite) ;
+               UpdateInfoIntersLineSurfTm( ptL, vtDir, dLen, 0, nT, Tria, vInfo, bFinite) ;
             }
          }
       }
@@ -124,19 +124,23 @@ IntersLineSurfTm( const Point3d& ptL, const Vector3d& vtL, double dLen, const IS
 //  Intersezione di molte linee parallele con una superficie TriMesh
 //----------------------------------------------------------------------------
 IntersParLinesSurfTm::IntersParLinesSurfTm( const Frame3d& frLines, const ISurfTriMesh& Stm)
-   : m_bOk( false), m_frLines( frLines), m_pSTm( &Stm)
+   : m_bOk( false), m_frLines( frLines), m_vpSTm( {&Stm})
 {
   // verifico esistenza superficie
-   if ( m_pSTm == nullptr  ||  ! m_pSTm->IsValid())
+   if ( m_vpSTm[0] == nullptr  ||  ! m_vpSTm[0]->IsValid())
       return ;
 
-  // creo HashGrid 2d
-   const int LIM_HG_TRIA = 127 ;
-   m_HGrids.SetActivationGrid( m_pSTm->GetTriangleCount() > LIM_HG_TRIA) ;
+  // aggiorno il vettore degli indici di base per mappare i triangoli con le rispettivi superfici
+  // ( in questo caso la superficie è unica, quindi ho solo due elementi)
+   m_vBaseInd = { 0, m_vpSTm[0]->GetTriangleCount()} ;
 
-  // riempio HashGrid 
+  // creo HashGrid 2d ed eventualmente attivo la griglia
+   const int LIM_HG_TRIA = 127 ;
+   m_HGrids.SetActivationGrid( m_vBaseInd.back() > LIM_HG_TRIA) ;
+
+  // riempio HashGrid
    Triangle3d Tria ;
-   int nT = Stm.GetFirstTriangle( Tria) ;
+   int nT = m_vpSTm[0]->GetFirstTriangle( Tria) ;
    while ( nT != SVT_NULL) {
      // calcolo il BBox del triangolo nel riferimento scelto
       Tria.ToLoc( m_frLines) ;
@@ -145,10 +149,55 @@ IntersParLinesSurfTm::IntersParLinesSurfTm( const Frame3d& frLines, const ISurfT
       b3Tria.Add( Tria.GetP( 1)) ;
       b3Tria.Add( Tria.GetP( 2)) ;
      // inserisco nella griglia
-      if ( ! m_HGrids.Add( nT, b3Tria))
+      if ( ! m_HGrids.Add( nT, b3Tria))  // ( 0 + nT, Tria)
          return ;
-     // passo al prossimo triangolo 
-      nT = Stm.GetNextTriangle( nT, Tria) ;
+     // passo al prossimo triangolo
+      nT = m_vpSTm[0]->GetNextTriangle( nT, Tria) ;
+   }
+  // aggiorno
+   m_bOk = m_HGrids.Update() ;
+}
+
+//----------------------------------------------------------------------------
+//  Intersezione di molte linee parallele con un vettore di superfici TriMesh
+//----------------------------------------------------------------------------
+IntersParLinesSurfTm::IntersParLinesSurfTm( const Frame3d& frLines, const CISURFTMPVECTOR& vStm)
+   : m_bOk( false), m_frLines( frLines), m_vpSTm( vStm), m_vBaseInd( {0})
+{
+  // verifico esistenza superfici
+   if ( m_vpSTm.empty())
+      return ;
+   for ( const ISurfTriMesh* pStm : m_vpSTm) {
+      if ( pStm == nullptr  ||  ! pStm->IsValid())
+         return ;
+   }
+
+  // aggiorno il vettore degli indici di base per mappare i triangoli con le rispettivi superfici
+  // NB. dal costruttore è già inizializzato a {0}
+   for ( int i = 0 ; i < int( m_vpSTm.size()) ; ++ i)
+      m_vBaseInd.emplace_back( m_vBaseInd.back() + m_vpSTm[i]->GetTriangleCount()) ;
+
+  // creo HashGrid 2d ed eventualmente attivo la griglia
+   const int LIM_HG_TRIA = 256 ;
+   m_HGrids.SetActivationGrid( m_vBaseInd.back() > LIM_HG_TRIA) ;
+
+  // riempio HashGrid
+   for ( int i = 0 ; i < int( m_vpSTm.size()) ; ++ i) {
+      Triangle3d Tria ;
+      int nT = m_vpSTm[i]->GetFirstTriangle( Tria) ;
+      while ( nT != SVT_NULL) {
+        // calcolo il BBox del triangolo nel riferimento scelto
+         Tria.ToLoc( m_frLines) ;
+         BBox3d b3Tria ;
+         b3Tria.Add( Tria.GetP( 0)) ;
+         b3Tria.Add( Tria.GetP( 1)) ;
+         b3Tria.Add( Tria.GetP( 2)) ;
+        // inserisco nella griglia ( aggiungo shift per indice del triangolo)
+         if ( ! m_HGrids.Add( m_vBaseInd[i] + nT, b3Tria))
+            return ;
+        // passo al prossimo triangolo
+         nT = m_vpSTm[i]->GetNextTriangle( nT, Tria) ;
+      }
    }
   // aggiorno
    m_bOk = m_HGrids.Update() ;
@@ -162,7 +211,7 @@ IntersParLinesSurfTm::GetInters( const Point3d& ptL, double dLen, ILSIVECTOR& vI
    if ( &vInfo == nullptr)
       return false ;
    vInfo.clear() ;
-  // verifico validità
+  // verifico validità
    if ( ! m_bOk)
       return false ;
 
@@ -174,15 +223,22 @@ IntersParLinesSurfTm::GetInters( const Point3d& ptL, double dLen, ILSIVECTOR& vI
    Point3d ptLL = ptL ;
    ptLL.ToGlob( m_frLines) ;
 
-  // recupero indici triangoli che intersecano box in 2d 
+  // recupero indici triangoli che intersecano box in 2d
    INTVECTOR vnIds ;
    if ( m_HGrids.Find( b3Line, vnIds)) {
       for ( int i = 0 ; i < int( vnIds.size()) ; ++ i) {
-         int nT = vnIds[i] ;
+        // recupero la superficie
+         int nInd = vnIds[i] ;
+         int nSurf = GetSurfInd( nInd) ;
+         if ( nSurf == -1)
+            return false ;
+        // recupero il triangolo
+         int nT = nInd - m_vBaseInd[nSurf] ;
          Triangle3d Tria ;
-         m_pSTm->GetTriangle( nT, Tria) ;
+         if ( ! m_vpSTm[nSurf]->GetTriangle( nT, Tria))
+            return false ;
         // aggiorno info con intersezione
-         UpdateInfoIntersLineSurfTm( ptLL, m_frLines.VersZ(), dLen, nT, Tria, vInfo, bFinite) ;
+         UpdateInfoIntersLineSurfTm( ptLL, m_frLines.VersZ(), dLen, nSurf, nT, Tria, vInfo, bFinite) ;
       }
    }
 
@@ -192,6 +248,35 @@ IntersParLinesSurfTm::GetInters( const Point3d& ptL, double dLen, ILSIVECTOR& vI
    return true ;
 }
 
+//----------------------------------------------------------------------------
+int
+IntersParLinesSurfTm::GetSurfInd( int nT) const
+{
+  // verifico la presenza di almeno un intervallo
+   if ( m_vBaseInd.size() < 2)
+      return -1 ;
+  // se la superficie è unica, allora non devo cercarla
+   if ( int( m_vBaseInd.size()) == 2)
+      return 0 ;
+  // ricerca binaria dell'intervallo contenente la posizione del triangolo
+   int nS = 0 ;
+   int nE = int( m_vBaseInd.size()) - 1 ;
+   while ( true) {
+      if ( nT < m_vBaseInd[nS]  ||  nT >= m_vBaseInd[nE])
+         return -1 ;
+      if ( nE - nS == 1)
+         return nS ;
+      int nM = ( nS + nE) / 2 ;
+      if ( nT == m_vBaseInd[nM])
+         return nM ;
+      if ( nT < m_vBaseInd[nM])
+         nE = nM ;
+      else
+         nS = nM ;
+   }
+   return -1 ;
+}
+
 //----------------------------------------------------------------------------
 bool
 FilterLineSurfTmInters( const ILSIVECTOR& vInfo, INTDBLVECTOR& vInters)
@@ -217,7 +302,7 @@ FilterLineSurfTmInters( const ILSIVECTOR& vInfo, INTDBLVECTOR& vInters)
    for ( size_t j = 1 ; j < vInters.size() ; ) {
      // intersezione precedente
       size_t i = j - 1 ;
-     // se hanno lo stesso parametro 
+     // se hanno lo stesso parametro
       if ( abs( vInters[i].second - vInters[j].second) < EPS_SMALL) {
         // se sono entrambe entranti o uscenti, elimino la seconda
          if ( ( vInters[i].first == LST_IN  &&  vInters[j].first == LST_IN)  ||

SurfBezier.cpp

@@ -42,22 +42,48 @@
 #include "/EgtDev/Include/EGkIntervals.h"
 #include "/EgtDev/Extern/Eigen/Dense"
 
-
 using namespace std ;
 
 //----------------------------------------------------------------------------
 GEOOBJ_REGISTER( SRF_BEZIER, NGE_S_BEZ, SurfBezier) ;
 
-static bool ChangeStartForClosed( PolyLine& plU0, PolyLine& plU1, ICurveComposite* pCrvU0, ICurveComposite* pCrvU1) ;
-static bool ParametrizeByLen( const ICurveComposite* pCurve, DBLVECTOR& vParam) ;
-static bool BuildCommonParam( const DBLMATRIX& mParam, DBLVECTOR& vCommonParam) ;
+//----------------------------------------------------------------------------
+static unordered_map<INTINT, DBLVECTOR, PairHashIntInt>  s_mBernCache ;  // mappa dei polinomi di bernstein
+static double s_dAuxSurfTol = 200 * EPS_SMALL ;
+static double s_dAuxSurfRefinedTol = 50 * EPS_SMALL ;
 
-static unordered_map<INTINT, DBLVECTOR, PairHashIntInt>  m_mBernCache ;  // mappa dei polinomi di bernstein
+//----------------------------------------------------------------------------
+void
+SetSurfBezierAuxSurfTol( double dTol)
+{
+   s_dAuxSurfTol = max( dTol, EPS_SMALL) ;
+}
+
+//----------------------------------------------------------------------------
+double
+GetSurfBezierAuxSurfTol( void)
+{
+   return s_dAuxSurfTol ;
+}
+
+//----------------------------------------------------------------------------
+void
+SetSurfBezierAuxSurfRefinedTol( double dTol)
+{
+   s_dAuxSurfRefinedTol = max( dTol, EPS_SMALL) ;
+}
+
+//----------------------------------------------------------------------------
+double
+GetSurfBezierAuxSurfRefinedTol( void)
+{
+   return s_dAuxSurfRefinedTol ;
+}
 
 //----------------------------------------------------------------------------
 SurfBezier::SurfBezier( void)
    : m_pSTM( nullptr), m_pSTMRefined( nullptr), m_nStatus( TO_VERIFY), m_nDegU(), m_nDegV(), m_nSpanU(), m_nSpanV(), m_bRat( false),
-     m_bTrimmed( false), m_bClosedU( false), m_bClosedV( false), m_pTrimReg(nullptr), m_nTempProp{0,0}, m_dTempParam{0.0,0.0}
+     m_bTrimmed( false), m_bClosedU( false), m_bClosedV( false), m_pTrimReg( nullptr), m_nTempProp{0,0}, m_dTempParam{0.0,0.0}, m_nIsPlanar( -1)
 {
 }
 
@@ -365,24 +391,26 @@ SurfBezier::GetCentroid( Point3d& ptCen) const
       return false ;
   // inizio con centro nell'origine
    ptCen = ORIG ;
-  // calcolo il baricentro
-   if ( m_pSTM == nullptr)
-      if ( ! GetAuxSurf())
-         return false ;
-   return m_pSTM->GetCentroid( ptCen) ;
+  // calcolo il baricentro della superficie ausiliaria
+   GetAuxSurf() ;
+   if ( m_pSTM != nullptr)
+      return m_pSTM->GetCentroid( ptCen) ;
+   else
+      return false ;
 }
 
 //----------------------------------------------------------------------------
 bool
 SurfBezier::GetBernstein( double dU, int nDegU, DBLVECTOR& vBernU) const
 {
-   INTINT key( nDegU, int( dU * pow(2,24))) ;
-   if ( m_mBernCache.find( key) == m_mBernCache.end()) {
+   const int TWO_TO_24 = 16777216 ;
+   INTINT key( nDegU, int( dU * TWO_TO_24)) ;
+   if ( s_mBernCache.find( key) == s_mBernCache.end()) {
       DBLVECTOR vBern( nDegU + 1) ;
       GetAllBernstein( dU, nDegU, vBern) ;
-      m_mBernCache[key] = vBern ;
+      s_mBernCache[key] = vBern ;
    }
-   vBernU = m_mBernCache[key] ;
+   vBernU = s_mBernCache[key] ;
    return true ;
 }
 
@@ -741,42 +769,49 @@ SurfBezier::GetPointNrmD1D2( double dU, double dV, Side nUs, Side nVs,
    if ( vtN.Normalize())
       return true ;
 
-  // se solo una delle due derivate è piccola, mi sposto lungo il relativo parametro e uso le tangenti
+  // se solo una sola delle due derivate è piccola, mi sposto lungo l'altro parametro
    if ( pvtDerU->Len() < EPS_SMALL  &&  pvtDerV->Len() > 10 * EPS_SMALL) {
-      double dCoeff = ( dU - 1000 * EPS_PARAM < 0. ? 1 : -1) ;
-      double dUm = dU + 1000 * EPS_PARAM * dCoeff ;
+      double dVm = dV ;
+      if ( nVs == ISurfBezier::FROM_MINUS)
+         dVm = ( dV - 100 * EPS_PARAM < 0 ? dV + 100 * EPS_PARAM : dV - 100 * EPS_PARAM) ;
+      else
+         dVm = ( dV + 100 * EPS_PARAM > m_nSpanV ? dV - 100 * EPS_PARAM : dV + 100 * EPS_PARAM) ;
       Point3d ptTmp ;
       Vector3d vtTmpU, vtTmpV ;
-      GetPointD1D2( dUm, dV, nUs, nVs, ptTmp, &vtTmpU, &vtTmpV) ;
-      vtN = ( *pvtDerV ^ vtTmpV) * dCoeff ;
+      GetPointD1D2( dU, dVm, nUs, nVs, ptTmp, &vtTmpU, &vtTmpV) ;
+      vtN = vtTmpU ^ vtTmpV ;
       if ( vtN.Normalize())
          return true ;
    }
    if ( pvtDerU->Len() > 10 * EPS_SMALL  &&  pvtDerV->Len() < EPS_SMALL) {
-      double dCoeff = ( dV - 1000 * EPS_PARAM < 0. ? 1 : -1) ;
-      double dVm = dV + 1000 * EPS_PARAM * dCoeff ;
+      double dUm = dU ;
+      if ( nUs == ISurfBezier::FROM_MINUS)
+         dUm = ( dU - 100 * EPS_PARAM < 0 ? dU + 100 * EPS_PARAM : dU - 100 * EPS_PARAM) ;
+      else
+         dUm = ( dU + 100 * EPS_PARAM > m_nSpanU ? dU - 100 * EPS_PARAM : dU + 100 * EPS_PARAM) ;
       Point3d ptTmp ;
       Vector3d vtTmpU, vtTmpV ;
-      GetPointD1D2( dU, dVm, nUs, nVs, ptTmp, &vtTmpU, &vtTmpV) ;
-      vtN = ( vtTmpU ^ *pvtDerU) * dCoeff ;
+      GetPointD1D2( dUm, dV, nUs, nVs, ptTmp, &vtTmpU, &vtTmpV) ;
+      vtN = vtTmpU ^ vtTmpV ;
       if ( vtN.Normalize())
          return true ;
    }
 
   // ricalcolo con una piccola variazione di entrambi i parametri 
-   double dUm ;
-   if ( dU - 100 * EPS_PARAM < 0.)
-      dUm = dU + 100 * EPS_PARAM ;
+   double dUm = dU ;
+   if ( nUs == ISurfBezier::FROM_MINUS)
+      dUm = ( dU - 100 * EPS_PARAM < 0 ? dU + 100 * EPS_PARAM : dU - 100 * EPS_PARAM) ;
    else
-      dUm = dU - 100 * EPS_PARAM ;
-   double dVm ;
-   if ( dV - 100 * EPS_PARAM < 0.)
-      dVm = dV + 100 * EPS_PARAM ;
+      dUm = ( dU + 100 * EPS_PARAM > m_nSpanU ? dU - 100 * EPS_PARAM : dU + 100 * EPS_PARAM) ;
+   double dVm = dV ;
+   if ( nVs == ISurfBezier::FROM_MINUS)
+      dVm = ( dV - 100 * EPS_PARAM < 0 ? dV + 100 * EPS_PARAM : dV - 100 * EPS_PARAM) ;
    else
-      dVm = dV - 100 * EPS_PARAM ;
+      dVm = ( dV + 100 * EPS_PARAM > m_nSpanV ? dV - 100 * EPS_PARAM : dV + 100 * EPS_PARAM) ;
    Point3d ptTmp ;
-   GetPointD1D2( dUm, dVm, nUs, nVs, ptTmp, pvtDerU, pvtDerV) ;
-   vtN = *pvtDerU ^ *pvtDerV ;
+   Vector3d vtTmpU, vtTmpV ;
+   GetPointD1D2( dUm, dVm, nUs, nVs, ptTmp, &vtTmpU, &vtTmpV) ;
+   vtN = vtTmpU ^ vtTmpV ;
    return vtN.Normalize() ;
 }
 
@@ -829,7 +864,7 @@ SurfBezier::CopyFrom( const SurfBezier& sbSrc)
          m_mCCEdge.back().emplace_back( sbSrc.m_mCCEdge[i][j]->Clone()) ;
       }
    }
-   if( sbSrc.m_bTrimmed) {
+   if ( sbSrc.m_bTrimmed) {
       for ( int i = 0 ; i < int( sbSrc.m_vCCLoop.size()) ; ++i)
          m_vCCLoop.emplace_back( sbSrc.m_vCCLoop[i]->Clone());
    }
@@ -1679,11 +1714,21 @@ SurfBezier::GetAuxSurf( void) const
       ResetAuxSurf() ;
       return nullptr ;
    }
-  // se già calcolata, la restituisco
-   if ( m_pSTM != nullptr)
-      return m_pSTM ;
+  // se già calcolata
+   if ( m_pSTM != nullptr) {
+     // se con la stessa tolleranza, la restituisco
+      if ( abs( s_dAuxSurfTol - m_pSTM->GetTempParam()) < EPS_SMALL)
+         return m_pSTM ;
+     // altrimenti la cancello prima di ricalcolarla
+      else {
+         delete( m_pSTM) ;
+         m_pSTM = nullptr ;
+      }
+   }
   // eseguo calcolo
-   m_pSTM = GetApproxSurf( 1000 * EPS_SMALL, 100 * EPS_SMALL) ;
+   m_pSTM = GetApproxSurf( s_dAuxSurfTol, 100 * EPS_SMALL) ;
+   if ( m_pSTM != nullptr)
+      m_pSTM->SetTempParam( s_dAuxSurfTol) ;
    return m_pSTM ;
 }
 
@@ -1697,10 +1742,20 @@ SurfBezier::GetAuxSurfRefined( void) const
       return nullptr ;
    }
   // se già calcolata, la restituisco
-   if ( m_pSTMRefined != nullptr)
-      return m_pSTMRefined ;
+   if ( m_pSTMRefined != nullptr) {
+     // se con la stessa tolleranza, la restituisco
+      if ( abs( s_dAuxSurfRefinedTol - m_pSTMRefined->GetTempParam()) < EPS_SMALL)
+         return m_pSTMRefined ;
+     // altrimenti la cancello prima di ricalcolarla
+      else {
+         delete( m_pSTMRefined) ;
+         m_pSTMRefined = nullptr ;
+      }
+   }
   // eseguo calcolo
-   m_pSTMRefined = GetApproxSurf( 50 * EPS_SMALL, 100 * EPS_SMALL) ;
+   m_pSTMRefined = GetApproxSurf( s_dAuxSurfRefinedTol, 100 * EPS_SMALL) ;
+   if ( m_pSTMRefined != nullptr)
+      m_pSTMRefined->SetTempParam( s_dAuxSurfRefinedTol) ;
    return m_pSTMRefined ;
 }
 
@@ -1714,11 +1769,11 @@ SurfBezier::GetApproxSurf( double dTol, double dSideMin) const
 
   // se c'è ausiliaria di visualizzazione con gli stessi parametri, ne restituisco una copia
    if ( m_pSTM != nullptr  &&
-        abs( dTol - 1000 * EPS_SMALL) < EPS_SMALL  &&  abs( dSideMin - 100 * EPS_SMALL) < EPS_SMALL)
+        abs( dTol - s_dAuxSurfTol) < EPS_SMALL  &&  abs( dSideMin - 100 * EPS_SMALL) < EPS_SMALL)
       return m_pSTM->Clone() ;
   // se c'è ausiliaria e richiesta con gli stessi parametri, ne restituisco una copia
    if ( m_pSTMRefined != nullptr  &&
-        abs( dTol - 50 * EPS_SMALL) < EPS_SMALL  &&  abs( dSideMin - 100 * EPS_SMALL) < EPS_SMALL)
+        abs( dTol - s_dAuxSurfRefinedTol) < EPS_SMALL  &&  abs( dSideMin - 100 * EPS_SMALL) < EPS_SMALL)
       return m_pSTMRefined->Clone() ;
   
   // costruttore della superficie
@@ -1726,28 +1781,18 @@ SurfBezier::GetApproxSurf( double dTol, double dSideMin) const
    POLYLINEMATRIX vvPL3d ;
    //POLYLINEVECTOR vPL ;    // per usare i polygon basic
    Tree Tree ;
-   if ( ! Tree.SetSurf( this, true))
+   if ( ! Tree.SetSurf( this))
       return nullptr ;
    BIPNTVECTOR vTrees ;
    Tree.GetIndependentTrees( vTrees) ;
-   bool bTest = false ; // per debug
    // resetto il vettore degli edge
    m_mCCEdge.clear() ;
    m_vCCLoop.clear() ;
    for ( int i = 0 ; i < (int) vTrees.size() ; ++ i) {
       Point3d ptMin = get<0>( vTrees[i]) ;
       Point3d ptMax = get<1>( vTrees[i]) ;
-      Tree.SetSurf( this, true, ptMin, ptMax) ;
-      if ( bTest) {
-         Tree.BuildTree_test() ;        // per debug
-         //Tree.BuildTree( 5 * LIN_TOL_FINE, 1) ; // per debug
-         Tree.SetTestMode() ;
-      }
-      else {
-         //Tree.BuildTree( 5 * LIN_TOL_FINE, 0.1) ;
-         Tree.BuildTree( dTol, dSideMin) ;
-         //Tree.BuildTree( 1, 5) ;  //debug
-      }
+      Tree.SetSurf( this, ptMin, ptMax) ;
+      Tree.BuildTree( dTol, dSideMin) ;
       if ( ! Tree.GetPolygons( vvPL, vvPL3d, m_mCCEdge, m_vCCLoop))
          continue ;
       //Tree.GetPolygonsBasic( vPL, true) ;  // per usare i polygon basic
@@ -1787,18 +1832,18 @@ SurfBezier::GetApproxSurf( double dTol, double dSideMin) const
       POLYLINEVECTOR vPL3d = vvPL3d[c] ;
 
       PNTVECTOR vPnt3d ;
-      for( int i = 0 ; i < int( vPL3d.size()) ; ++i) {
+      for ( int i = 0 ; i < int( vPL3d.size()) ; ++i) {
          PolyLine& pl3d = vPL3d[i] ;
          Point3d pt3d ; pl3d.GetFirstPoint( pt3d) ;
-         if( vPL3d.size() > 1)
+         if ( vPL3d.size() > 1)
             vPnt3d.push_back( pt3d) ;
          while ( pl3d.GetNextPoint( pt3d)) {
             vPnt3d.push_back( pt3d) ;
          }
       }
 
-      // se ho ottenuto meno triangoli di quelli che avrei dovuto avere allora potrei aver avuto un problema in prossimità di un polo
-      // se comunque ho corrispondenza tra vPnt e vPnt3d allora eseguo la trinagolazione in 3d
+     // se ho ottenuto meno triangoli di quelli che avrei dovuto avere allora potrei aver avuto un problema in prossimità di un polo
+     // se comunque ho corrispondenza tra vPnt e vPnt3d allora eseguo la trinagolazione in 3d
       int nTriaNumber = ( vPnt.size() - 2) + (2 * (vPL.size() - 1)) - (vPL.size() > 2 ? vPL.size() - 2 : 0)  - (vPL.size() != 1 ? vPL.size() : 0) ;
       bool bTriangulationFailedIn2D = nTriaNumber != ( vTria.size()) / 3 ;
       //bool bMismatch2D3D = vPnt.size() != vPnt3d.size() ;
@@ -1810,13 +1855,13 @@ SurfBezier::GetApproxSurf( double dTol, double dSideMin) const
          if ( ! Tri.Make( vPL3d, vPnt, vTria))
             bTriangulationSucceded = false ;
          bTriangulationSucceded = bTriangulationSucceded  &&  ( vPnt.size() - 2) + (2 * (vPL.size() - 1)) - (vPL.size() > 2 ? vPL.size() - 2 : 0) - (vPL.size() != 1 ? vPL.size() : 0) == ( vTria.size() ) / 3 ;
-         if( bTriangulationFailedIn2D) {
+         if ( bTriangulationFailedIn2D) {
             if ( bTriangulationSucceded)
                LOG_INFO( GetEGkLogger(), "Info : Last problem in MakeByEC23(1) RESOLVED")
             else
                LOG_INFO( GetEGkLogger(), "Info : This problem in MakeByEC23(1) is the same as the previous one")
          }
-         if( bTriangulationSucceded) {
+         if ( bTriangulationSucceded) {
             bTriangulatedIn3D = true ;
             // calcolo la normale della polyline per flippare eventualmente i triangoli se hanno la normale sbagliata
             PolyLine& pl3d = vPL3d[0] ;
@@ -1825,7 +1870,7 @@ SurfBezier::GetApproxSurf( double dTol, double dSideMin) const
             Vector3d vtN = plPlane.GetVersN() ;
             Triangle3d tria ;
             tria.Set( vPnt[vTria[0]], vPnt[vTria[1]], vPnt[vTria[2]]) ;
-            if( tria.GetN() * vtN < 0)
+            if ( tria.GetN() * vtN < 0)
                reverse( vTria.begin(), vTria.end()) ;
          }
          else {
@@ -1835,11 +1880,11 @@ SurfBezier::GetApproxSurf( double dTol, double dSideMin) const
       }
 
       // riordino il vettore dei punti su cui non ho fatto la triangolazione
-      if( int(vPL.size()) == 2) {
+      if ( vPL.size() == 2) {
          //if( vPnt.size() != vPnt3d.size())
          //   return nullptr ;
          PNTVECTOR vPntOrd ;
-         if( bTriangulatedIn3D) {
+         if ( bTriangulatedIn3D) {
             ReorderPntVector( vPL3d, true, vPnt, vPL, vPntOrd) ;
             vPnt3d = vPnt ;
             vPnt = vPntOrd ;
@@ -1850,7 +1895,7 @@ SurfBezier::GetApproxSurf( double dTol, double dSideMin) const
          }
       }
       else if ( bTriangulatedIn3D) {
-         if( vPL.size() == 1) {
+         if ( vPL.size() == 1) {
             vPnt3d = vPnt ;
             PNTVECTOR vPnt2d ;
             for( int i = 0 ; i < int( vPL.size()) ; ++i) {
@@ -1871,9 +1916,9 @@ SurfBezier::GetApproxSurf( double dTol, double dSideMin) const
          }
       }
 
-      //controllo che i due vettori vPnt e vPnt3d abbiano la stessa lunghezza, sennò vuol dire che nel vettore vPnt3d ho avuto dei Rejected e devo ricalcolarli
-      // i punti in eccesso verranno poi scartati dalla trimesh
-      if( vPnt.size() != vPnt3d.size()) {
+     // controllo che i due vettori vPnt e vPnt3d abbiano la stessa lunghezza, sennò vuol dire che nel vettore vPnt3d ho avuto dei Rejected e devo ricalcolarli
+     // i punti in eccesso verranno poi scartati dalla trimesh
+      if ( vPnt.size() != vPnt3d.size()) {
          vPnt3d.clear() ;
          for ( int i = 0 ; i < int( vPnt.size()) ; ++ i) {
             Point3d pt3d ;
@@ -2029,23 +2074,15 @@ bool
 SurfBezier::GetLeaves( vector<tuple<int, Point3d, Point3d>>& vLeaves) const
 {
    Tree Tree ;
-   if ( ! Tree.SetSurf( this, true))
+   if ( ! Tree.SetSurf( this))
       return false ;
    BIPNTVECTOR vTrees ;
    Tree.GetIndependentTrees( vTrees) ;
    for ( int i = 0 ; i < int( vTrees.size()) ; ++ i) {
       Point3d ptMin = get<0>( vTrees[i]) ;
       Point3d ptMax = get<1>( vTrees[i]) ;
-      Tree.SetSurf( this, true, ptMin, ptMax) ;
-      bool bTest = false ; // per debug
-      if ( bTest) {
-         Tree.BuildTree_test() ;                    // per debug
-         //Tree.BuildTree( 5 * LIN_TOL_FINE, 1) ;     // per debug
-      }
-      else {
-         //Tree.BuildTree( 100 * LIN_TOL_FINE, 0.1) ;
-         Tree.BuildTree( 5 * LIN_TOL_FINE, 0.1) ;
-      }
+      Tree.SetSurf( this, ptMin, ptMax) ;
+      Tree.BuildTree( s_dAuxSurfTol, 100 * EPS_SMALL) ;
       vector<Cell> vCells ;
       Tree.GetLeaves( vCells) ;
       for ( int k = 0 ; k < int( vCells.size()) ; ++ k) {
@@ -2825,7 +2862,7 @@ SurfBezier::UnprojectPointFromStm( int nT, const Point3d& ptI, Point3d& ptSP, in
          for ( int i = 0 ; i < int( m_mCCEdge[c].size()) ; ++i) {
             if ( ! m_mCCEdge[c][i]->IsValid()) {
                Point3d pt ;
-               if ( ! m_mCCEdge[c][i]->GetOnlyPoint(pt)) 
+               if ( ! m_mCCEdge[c][i]->GetOnlyPoint( pt)) 
                   return false ;
                vInters[c] = AreSamePointApprox( pt, ptI) ? 1 : 0 ;
                if ( vInters[c] == 1)
@@ -3177,8 +3214,8 @@ SurfBezier::UnprojectPoint( const Point3d& pt3D, Point3d& ptParam, const Point3d
       }
       // calcolo le nuove coordinate
       Vector3d vtDir ;
-      double dASum = abs(dfdU) + abs(dfdV) ;
-      double dSSum = sqrt( pow(dfdU,2) + pow( dfdV,2)) ;
+      double dASum = abs( dfdU) + abs( dfdV) ;
+      double dSSum = sqrt( dfdU * dfdU + dfdV * dfdV) ;
       vtDir.Set( - dfdU, - dfdV, 0) ;
       if ( ! vtDir.Normalize() )
          vtDir.Set( - dfdU / dSSum, - dfdV / dSSum, 0) ;
@@ -4191,8 +4228,8 @@ SurfBezier::CreateByPointCurve( const Point3d& pt, const ICurve* pCurve)
       return false ;
    // recupero span e grado nel parametro U
    int nSpanU = int( CrvU.GetCurveCount()) ;
-   int nDegU = ( GetCurveBezier( CrvU.GetCurve(0)))->GetDegree() ;
-   bool bRat = (GetCurveBezier( CrvU.GetCurve(0)))->IsRational() ;
+   int nDegU = GetCurveBezier( CrvU.GetCurve(0))->GetDegree() ;
+   bool bRat = GetCurveBezier( CrvU.GetCurve(0))->IsRational() ;
    // in V decido che la funzione è una patch di grado 1 
    int nSpanV = 1 ;
    int nDegV = 1 ;
@@ -4219,6 +4256,52 @@ SurfBezier::CreateByPointCurve( const Point3d& pt, const ICurve* pCurve)
    return true ;
 }
 
+//----------------------------------------------------------------------------
+static bool
+ChangeStartForClosed( PolyLine& plU0, PolyLine& plU1, ICurveComposite* pCrvU0, ICurveComposite* pCrvU1) {
+   // se sono chiuse devo controllare che gli start siano il più allineati possibile, se non lo sono cambio gli start
+   if ( plU0.IsClosed() && plU1.IsClosed()) {
+      vector<tuple<double,int,Point3d,int,Point3d>> vDistVert ;
+      tuple<double,int,Point3d,int,Point3d> tMatch ;
+      Point3d pt0 ;
+      bool bOk0 = plU0.GetFirstPoint( pt0) ;
+      int c0 = 0 ;
+      double dMinDist = INFINITO ;
+      while ( bOk0) {
+         Point3d pt1 ;
+         bool bOk1 = plU1.GetFirstPoint( pt1) ;
+         int c1 = 0 ;
+         while ( bOk1) {
+            double dDist = Dist( pt0, pt1) ;
+            if ( dDist < dMinDist) {
+               tMatch = make_tuple( dDist, c0, pt0, c1, pt1) ;
+               dMinDist = dDist ;
+            }
+            ++c1 ;
+            bOk1 = plU1.GetNextPoint( pt1) ;
+         }
+         vDistVert.push_back( tMatch) ;
+         dMinDist = INFINITO ;
+         c1 = 0 ;
+         ++c0 ;
+         bOk0 = plU0.GetNextPoint( pt0) ;
+      }
+      int nMin = 0 ;
+      dMinDist = INFINITO ;
+      for ( int i = 0 ; i < int( vDistVert.size()) ; ++i) {
+         if( get<0>(vDistVert[i]) < dMinDist) {
+            dMinDist = get<0>(vDistVert[i]) ;
+            nMin = i ;
+         }
+      }
+      ChangePolyLineStart(plU0, get<2>(vDistVert[nMin]), EPS_SMALL) ;
+      ChangePolyLineStart(plU1, get<4>(vDistVert[nMin]), EPS_SMALL) ;
+      pCrvU0->ChangeStartPoint( double( get<1>(vDistVert[nMin]))) ;
+      pCrvU1->ChangeStartPoint( double( get<3>(vDistVert[nMin]))) ;
+   }
+   return true ;
+}
+
 //----------------------------------------------------------------------------
 bool
 SurfBezier::CreateByTwoCurves( const ICurve* pCurve0, const ICurve* pCurve1, int nRuledType)
@@ -5098,52 +5181,6 @@ SplitByCommonParam( ICURVEPOVECTOR& vCrvBezUnif, DBLVECTOR& vCommonParam, DBLMAT
    return true ;
 }
 
-static bool
-ChangeStartForClosed( PolyLine& plU0, PolyLine& plU1, ICurveComposite* pCrvU0, ICurveComposite* pCrvU1) {
-   // se sono chiuse devo controllare che gli start siano il più allineati possibile, se non lo sono cambio gli start
-   if ( plU0.IsClosed() && plU1.IsClosed()) {
-      vector<tuple<double,int,Point3d,int,Point3d>> vDistVert ;
-      tuple<double,int,Point3d,int,Point3d> tMatch ;
-      Point3d pt0 ;
-      bool bOk0 = plU0.GetFirstPoint( pt0) ;
-      int c0 = 0 ;
-      double dMinDist = INFINITO ;
-      while ( bOk0) {
-         Point3d pt1 ;
-         bool bOk1 = plU1.GetFirstPoint( pt1) ;
-         int c1 = 0 ;
-         while ( bOk1) {
-            double dDist = Dist( pt0, pt1) ;
-            if ( dDist < dMinDist) {
-               tMatch = make_tuple( dDist, c0, pt0, c1, pt1) ;
-               dMinDist = dDist ;
-            }
-            ++c1 ;
-            bOk1 = plU1.GetNextPoint( pt1) ;
-         }
-         vDistVert.push_back( tMatch) ;
-         dMinDist = INFINITO ;
-         c1 = 0 ;
-         ++c0 ;
-         bOk0 = plU0.GetNextPoint( pt0) ;
-      }
-      int nMin = 0 ;
-      dMinDist = INFINITO ;
-      for ( int i = 0 ; i < int( vDistVert.size()) ; ++i) {
-         if( get<0>(vDistVert[i]) < dMinDist) {
-            dMinDist = get<0>(vDistVert[i]) ;
-            nMin = i ;
-         }
-      }
-      ChangePolyLineStart(plU0, get<2>(vDistVert[nMin]), EPS_SMALL) ;
-      ChangePolyLineStart(plU1, get<4>(vDistVert[nMin]), EPS_SMALL) ;
-      pCrvU0->ChangeStartPoint( double( get<1>(vDistVert[nMin]))) ;
-      pCrvU1->ChangeStartPoint( double( get<3>(vDistVert[nMin]))) ;
-   }
-   return true ;
-}
-
-
 //----------------------------------------------------------------------------
 bool
 SurfBezier::CreateBySetOfCurves( const ICURVEPOVECTOR& vCrvBez, bool bReduceToDeg3)
@@ -5311,7 +5348,7 @@ SurfBezier::CreateBySetOfCurves( const ICURVEPOVECTOR& vCrvBez, bool bReduceToDe
                Point3d pt2 = ptP2 ; pt2.ToLoc( frParab) ;
 
                Eigen::Matrix3d mA ;
-               mA.col(0) << pow(pt0.x,2), pow(pt1.x,2) , pow(pt2.x,2) ;
+               mA.col(0) << pt0.x * pt0.x, pt1.x * pt1.x, pt2.x * pt2.x ;
                mA.col(1) << pt0.x, pt1.x , pt2.x ;
                mA.col(2) << 1, 1, 1 ;
                if( abs( mA.determinant()) < EPS_SMALL)

SurfBezier.h

@@ -81,10 +81,10 @@ class SurfBezier : public ISurfBezier, public IGeoObjRW
       bool GetVolume( double& dVolume) const override
                { if ( &dVolume == nullptr)
                     return false ;
-                 if( m_pSTM == nullptr)
+                 if ( m_pSTM == nullptr)
                     GetAuxSurf() ;
                  dVolume = 0 ;
-                 if( m_pSTM != nullptr)
+                 if ( m_pSTM != nullptr)
                     m_pSTM->GetVolume( dVolume) ;
                  return true ; }
       bool GetCentroid( Point3d& ptCen) const override ;

Tree.cpp

@@ -21,12 +21,13 @@
 #include "SurfFlatRegion.h"
 #include "IntersLineLine.h"
 #include "AdjustLoops.h"
-#include "/EgtDev/Include/EGkDistLineLine.h"
 #include "/EgtDev/Include/EGkPolyLine.h"
-#include "/EgtDev/Include/EGkDistPointCurve.h"
 #include "/EgtDev/Include/EGkCurve.h"
-#include "DistPointCrvComposite.h"
 #include "/EgtDev/Include/EGkSfrCreate.h"
+#include "/EgtDev/Include/EGkDistPointCurve.h"
+#include "/EgtDev/Include/EGkDistLineLine.h"
+#include "/EgtDev/Include/EGkDistPointTria.h"
+#include "/EgtDev/Include/EGkStringUtils3d.h"
 #include <algorithm>
 #include <numeric>
 
@@ -34,8 +35,8 @@ using namespace std ;
 
 //----------------------------------------------------------------------------
 Tree::Tree( void)
-   : m_pSrfBz( nullptr), m_bTrimmed( false), m_bBilinear( false), m_bMulti( false), m_bClosedU( false), m_bClosedV( false), m_vbPole( { false, false, false, false}),
-     m_bSplitPatches( true), m_bTestMode( false)
+   : m_pSrfBz( nullptr), m_bTrimmed( false), m_bBilinear( false), m_bMulti( false), m_bClosedU( false), m_bClosedV( false),
+     m_vbPole{ false, false, false, false}
 {
    Point3d ptBl( 0, 0), ptTr ( 1 * SBZ_TREG_COEFF, 1 * SBZ_TREG_COEFF) ;
    Cell cRoot( ptBl, ptTr) ;
@@ -49,8 +50,8 @@ Tree::~Tree( void)
 
 //----------------------------------------------------------------------------
 Tree::Tree( const Point3d ptBl, const Point3d ptTr)
-   : m_pSrfBz( nullptr), m_bTrimmed( false), m_bBilinear( false), m_bMulti( false), m_bClosedU( false), m_bClosedV( false), m_vbPole( { false, false, false, false}),
-     m_bSplitPatches( true), m_bTestMode( false)
+   : m_pSrfBz( nullptr), m_bTrimmed( false), m_bBilinear( false), m_bMulti( false), m_bClosedU( false), m_bClosedV( false),
+     m_vbPole{ false, false, false, false}
 {
    Cell cRoot( ptBl, ptTr) ;
    m_mTree.insert( pair< int, Cell>( -1, cRoot)) ;
@@ -101,32 +102,40 @@ Tree::AdjustLoop( PolyLine& pl, POLYLINEVECTOR& vPl, BOOLVECTOR& vbOrientation)
 }
 
 //----------------------------------------------------------------------------
-bool
-Tree::GetPoint( double dU, double dV, Point3d& pt) const
+static double
+GetHalfKey( double dPar)
 {
-   pair<int64_t, int64_t> key (static_cast<int64_t>(dU * pow(2,15)), static_cast<int64_t>(dV * pow(2,15))) ;
+   const int TWO_TO_15 = 32768 ;
+   return ( dPar * TWO_TO_15) ;
+}
+
+//----------------------------------------------------------------------------
+bool
+Tree::GetPoint( double dU, double dV, Point3d& ptP) const
+{
+   pair<int64_t, int64_t> key( int64_t( GetHalfKey( dU)), int64_t( GetHalfKey( dV))) ;
    bool bOk = true ;
    if ( m_mPt3d.find( key) == m_mPt3d.end()) {
-      bOk = bOk &&  m_pSrfBz->GetPoint( dU / SBZ_TREG_COEFF, dV / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, pt) ;
-      m_mPt3d[key] = pt ;
+      bOk = bOk &&  m_pSrfBz->GetPoint( dU / SBZ_TREG_COEFF, dV / SBZ_TREG_COEFF, ISurfBezier::FROM_MINUS, ISurfBezier::FROM_MINUS, ptP) ;
+      m_mPt3d[key] = ptP ;
    }
-   pt = m_mPt3d[key] ;
+   ptP = m_mPt3d[key] ;
    return bOk ;
 }
 
 //----------------------------------------------------------------------------
 bool
-Tree::SavePoint( double dU, double dV, Point3d& pt)
+Tree::SavePoint( double dU, double dV, Point3d& ptP)
 {
-   pair<int64_t, int64_t> key (static_cast<int64_t>(dU * pow(2,15)), static_cast<int64_t>(dV * pow(2,15))) ;
+   pair<int64_t, int64_t> key( int64_t( GetHalfKey( dU)), int64_t( GetHalfKey( dV))) ;
    if ( m_mPt3d.find( key) == m_mPt3d.end())
-      m_mPt3d[key] = pt ;
+      m_mPt3d[key] = ptP ;
    return true ;
 }
 
 //----------------------------------------------------------------------------
 bool
-Tree::SetSurf( const SurfBezier* pSrfBz, bool bSplitPatches, const Point3d& ptMin, const Point3d& ptMax)
+Tree::SetSurf( const SurfBezier* pSrfBz, const Point3d& ptMin, const Point3d& ptMax)
 {
    if ( pSrfBz == nullptr  ||  ! pSrfBz->IsValid())
       return false ;
@@ -140,7 +149,6 @@ Tree::SetSurf( const SurfBezier* pSrfBz, bool bSplitPatches, const Point3d& ptMi
    m_vCCLoop2D.clear() ;
 
    m_pSrfBz = pSrfBz ;
-   m_bSplitPatches = bSplitPatches ;
   // le coordinate delle celle sono nello spazio parametrico
    int nDegU, nDegV, nSpanU, nSpanV ;
    bool bIsRat, bTrimmed ;
@@ -230,43 +238,52 @@ Tree::SetSurf( const SurfBezier* pSrfBz, bool bSplitPatches, const Point3d& ptMi
   // se richiesto divido preliminarmente le patches
    m_vnParents.clear() ;
    bool bIsPlanar = m_pSrfBz->IsPlanar() ;
-   if( ! bIsPlanar  ||  m_bMulti) {
-      // se la superficie è chiusa lungo il parametro U, sistemo le adiacenze al bordo
-      if ( AreSamePointApprox(ptP00, ptP10) && AreSamePointApprox(ptP01, ptP11) ) {
-         m_mTree[-1].m_nLeft = -1 ;
-         m_mTree[-1].m_nRight = -1 ;
-         m_bClosedU = true ;
+   if ( ! bIsPlanar  ||  m_bMulti) {
+     // se la superficie è chiusa lungo il parametro U, sistemo le adiacenze al bordo
+      if ( AreSamePointApprox( ptP00, ptP10) && AreSamePointApprox( ptP01, ptP11) ) {
+         Point3d ptP0001 ;
+         GetPoint( ptMin.x, ( ptMin.y + ptTop.y) / 2, ptP0001) ;
+         Point3d ptP1011 ;
+         GetPoint( ptTop.x, ( ptMin.y + ptTop.y) / 2, ptP1011) ;
+         if ( AreSamePointApprox( ptP0001, ptP1011)) {
+            m_mTree[-1].m_nLeft = -1 ;
+            m_mTree[-1].m_nRight = -1 ;
+            m_bClosedU = true ;
+         }
       }
-      // se la superficie è chiusa lungo il parametro V, sistemo le adiacenze al bordo
-      if ( ( AreSamePointApprox(ptP00, ptP01) && AreSamePointApprox(ptP10, ptP11) ) ) {
-         m_mTree[-1].m_nTop = -1 ;
-         m_mTree[-1].m_nBottom = -1 ;
-         m_bClosedV = true ;
+     // se la superficie è chiusa lungo il parametro V, sistemo le adiacenze al bordo
+      if ( ( AreSamePointApprox( ptP00, ptP01) && AreSamePointApprox( ptP10, ptP11) ) ) {
+         Point3d ptP0010 ;
+         GetPoint( ( ptMin.x + ptTop.x) / 2, ptMin.y, ptP0010) ;
+         Point3d ptP0111 ;
+         GetPoint( (ptMin.x + ptTop.x) / 2, ptTop.y, ptP0111) ;
+         if ( AreSamePointApprox( ptP0010, ptP0111)) {
+            m_mTree[-1].m_nTop = -1 ;
+            m_mTree[-1].m_nBottom = -1 ;
+            m_bClosedV = true ;
+         }
       }
-      if ( (m_bSplitPatches  &&  ( nSpanU > 1  ||  nSpanV > 1))  ||  m_bClosedU  ||  m_bClosedV) {
+      if ( nSpanU > 1  ||  nSpanV > 1  ||  m_bClosedU  ||  m_bClosedV) {
          int nId = -1 ;
          int nSplit = nSpanU ;
          bool bSplitInHalf = false ;
-         // se la superficie è chiusa lungo un parametro la splitto lungo quel parametro almeno una volta
-         if( m_bClosedU  &&  nSpanU == 1) {
+        // se la superficie è chiusa lungo un parametro la splitto lungo quel parametro almeno una volta
+         if ( m_bClosedU  &&  nSpanU == 1) {
             nSplit = 2 ;
             bSplitInHalf = true ;
          }
          for ( int i = 1 ; i < nSplit ; ++i) {
             m_mTree[nId].SetSplitDirVert( true) ;
-            double dSplit = bSplitInHalf ? i * SBZ_TREG_COEFF / 2 : i * SBZ_TREG_COEFF ;
-            if ( Split( nId, dSplit)) {
-               //++ nId ;
-               //++ nId ;
+            double dSplit = ( bSplitInHalf ? i * SBZ_TREG_COEFF / 2 : i * SBZ_TREG_COEFF) ;
+            if ( Split( nId, dSplit))
                nId = m_mTree[nId].m_nChild2 ;
-            }
          }
          INTVECTOR vLeaves ;
          GetHeightLeaves( -1, vLeaves) ;
          nSplit = nSpanV - 1 ;
          bSplitInHalf = false ;
-         // se la superficie è chiusa lungo un parametro la splitto lungo quel parametro almeno una volta
-         if( m_bClosedV &&  nSpanV == 1) {
+        // se la superficie è chiusa lungo un parametro la splitto lungo quel parametro almeno una volta
+         if ( m_bClosedV &&  nSpanV == 1) {
             nSplit = 1 ;
             bSplitInHalf = true ;
          }
@@ -284,20 +301,8 @@ Tree::SetSurf( const SurfBezier* pSrfBz, bool bSplitPatches, const Point3d& ptMi
      // se è chiusa e non ho già fatto split preliminare, splitto sul parametro su cui è chiusa
      // e sistemo le adiacenze
       if ( ( AreSamePointApprox( ptP00, ptP01)  ||  AreSamePointApprox( ptP10, ptP11)) ||
-            ( AreSamePointApprox( ptP00, ptP10)  ||  AreSamePointApprox( ptP01, ptP11))) {
+           ( AreSamePointApprox( ptP00, ptP10)  ||  AreSamePointApprox( ptP01, ptP11))) {
          if ( ( AreSamePointApprox( ptP00, ptP01)  ||  AreSamePointApprox( ptP10, ptP11))) {
-            
-            ////questo in teoria non serve più perché lo faccio appena sopra
-            if( int( m_mTree.size()) == 1) { 
-               if ( AreSamePointApprox( ptP00, ptP01)  &&  AreSamePointApprox( ptP10, ptP11)) {
-                  m_mTree[-1].m_nTop = -1 ;
-                  m_mTree[-1].m_nBottom = -1 ;
-                  m_bClosedV = true ;
-               }
-               m_mTree[-1].SetSplitDirVert( false) ;
-               Split( -1) ;
-            }
-            ////////
 
            // qui devo fare il controllo capped ( chiusura a semisfera)
            // devo controllare se i punti ai parametri U=0 e U=1 sono tutti coincidenti
@@ -314,7 +319,7 @@ Tree::SetSurf( const SurfBezier* pSrfBz, bool bSplitPatches, const Point3d& ptMi
             }
             m_vbPole[1] = bPole0 ;
             m_vbPole[3] = bPole1 ;
-            if ( bPole0  &&  bPole1  &&  int( m_mTree.size() == 3)) {
+            if ( bPole0  &&  bPole1  &&  m_mTree.size() == 3) {
                m_mTree[0].SetSplitDirVert( true) ;
                Split( 0) ;
                m_mTree[1].SetSplitDirVert( true) ;
@@ -323,29 +328,16 @@ Tree::SetSurf( const SurfBezier* pSrfBz, bool bSplitPatches, const Point3d& ptMi
          }
         // nella condizione di questo if non controllo eventuali divisioni preliminari, perché ne tengo conto dopo
          if ( AreSamePointApprox( ptP00, ptP10)  ||  AreSamePointApprox( ptP01, ptP11)) {
-            
-            //// questo in teoria non serve più perché lo faccio sopra
-            if ( m_mTree.size() == 1) {
-               if ( AreSamePointApprox( ptP00, ptP10)  &&  AreSamePointApprox( ptP01, ptP11)) {
-                  m_mTree[-1].m_nLeft = -1 ;
-                  m_mTree[-1].m_nRight = -1 ;
-                  m_bClosedU = true ;
-               }
-               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
+           // devo controllare se i punti ai parametri V=0 e V=1 sono tutti coincidenti
+           // in caso devo fare uno split nell'altra direzione
             bool bOk = false ;
             bool bPole0 = true, bPole1 = true ;
-            Point3d ptV0, ptV1 ;
-            // controllo se tutti i punti sull'isoparametrica sono uguali
+           // controllo se tutti i punti sull'isoparametrica sono uguali
             for ( int i = 1 ; i < nDegU * nSpanU + 1 ; ++ i) {
-               ptV0 = m_pSrfBz->GetControlPoint( i, &bOk) ;
+               Point3d ptV0 = m_pSrfBz->GetControlPoint( i, &bOk) ;
                bPole0 = bPole0  &&  AreSamePointApprox( ptP00, ptV0) ;
-               ptV1 = m_pSrfBz->GetControlPoint( i + ( nDegU * nSpanU + 1) * ( nDegV * nSpanV), &bOk) ;
+               Point3d ptV1 = m_pSrfBz->GetControlPoint( i + ( nDegU * nSpanU + 1) * ( nDegV * nSpanV), &bOk) ;
                bPole1 = bPole1  &&  AreSamePointApprox( ptP01, ptV1) ;
             }
             m_vbPole[0] = bPole1 ;
@@ -357,7 +349,7 @@ Tree::SetSurf( const SurfBezier* pSrfBz, bool bSplitPatches, const Point3d& ptMi
                Split( 1) ;
             }
            // se ho fatto solo 1 split orizzontale e ho due celle foglie nId = 0 e nId = 1
-            if ( int( m_mTree.size() == 3)  &&  ! m_mTree.at(-1).IsSplitVert()) {
+            if ( m_mTree.size() == 3  &&  ! m_mTree.at(-1).IsSplitVert()) {
                m_mTree[0].m_nLeft = -1 ;
                m_mTree[0].m_nRight = -1 ;
                m_mTree[1].m_nLeft = -1 ;
@@ -370,41 +362,11 @@ Tree::SetSurf( const SurfBezier* pSrfBz, bool bSplitPatches, const Point3d& ptMi
          }
       }
    }
-   // calcolo i parent che ho creato con le eventuali divisioni preliminari
+  // calcolo i parent che ho creato con le eventuali divisioni preliminari
    INTVECTOR vLeaves ;
    GetHeightLeaves( -1, vLeaves) ;
    m_vnParents = vLeaves ;
-   
-  // calcolo e salvo la lunghezza reale delle curve di bezier di bordo
-   PtrOwner<CurveComposite> pCrvV0( m_pSrfBz->GetCurveOnU( 0)) ;
-   PtrOwner<CurveComposite> pCrvV1( m_pSrfBz->GetCurveOnU( double(nSpanV))) ;
-   PtrOwner<CurveComposite> pCrvU0( m_pSrfBz->GetCurveOnV( 0)) ;
-   PtrOwner<CurveComposite> pCrvU1( m_pSrfBz->GetCurveOnV( double(nSpanU))) ;
-   double dLen0 ; pCrvV0->GetApproxLength( dLen0) ;
-   double dLen1 ; pCrvU1->GetApproxLength( dLen1) ;
-   double dLen2 ; pCrvV1->GetApproxLength( dLen2) ;
-   double dLen3 ; pCrvU0->GetApproxLength( dLen3) ;
-   if ( dLen0 < EPS_ZERO  &&  dLen2 < EPS_ZERO ) {
-      PtrOwner<CurveComposite> pCrvV( m_pSrfBz->GetCurveOnU( double(nSpanV) / 2)) ;
-      pCrvV->GetApproxLength( dLen0) ;
-      if ( dLen0 < EPS_ZERO ) {
-         pCrvV.Set( m_pSrfBz->GetCurveOnU( double(nSpanV) / 4)) ;
-         pCrvV->GetApproxLength( dLen0) ;
-      }
-   }
-   if ( dLen1 < EPS_ZERO  &&  dLen3 < EPS_ZERO ) {
-      PtrOwner<CurveComposite> pCrvU( m_pSrfBz->GetCurveOnV( double(nSpanU) / 2)) ;
-      pCrvU->GetApproxLength( dLen1) ;
-      if ( dLen1 < EPS_ZERO ) {
-         pCrvU.Set( m_pSrfBz->GetCurveOnV( double(nSpanU) / 4)) ;
-         pCrvU->GetApproxLength( dLen1) ;
-      }
-   }
-   m_vDim.clear() ;
-   m_vDim.push_back( ( dLen0 > EPS_ZERO ? dLen0 : 1)) ;
-   m_vDim.push_back( ( dLen1 > EPS_ZERO ? dLen1 : 1)) ;
-   m_vDim.push_back( ( dLen2 > EPS_ZERO ? dLen2 : 1)) ;
-   m_vDim.push_back( ( dLen3 > EPS_ZERO ? dLen3 : 1)) ;
+
    return true ;
 }
 
@@ -419,7 +381,7 @@ AddOrMergeBBox( const BBox3d& bBox3dA, vector<BBox3d>& vBBox, bool bAdd = true,
       return true ;
    }
    bool bAdded = false ;
-   for ( int b = 0 ; b < (int)vBBox.size() ; ++b) {
+   for ( int b = 0 ; b < int( vBBox.size()) ; ++b) {
       BBox3d bBox3dB = vBBox[b] ;
       BBox3d b3Int ;
      // se sono celle diverse e ho un'intersezione faccio il merge
@@ -565,7 +527,7 @@ bool
 Tree::Split( int nId)
 {
    double dValue ;
-   Cell& cCell =  m_mTree.at(nId) ;
+   Cell& cCell =  m_mTree.at( nId) ;
    if ( cCell.IsSplitVert())
       dValue = ( cCell.GetBottomLeft().x + cCell.GetTopRight().x) / 2 ;
    else
@@ -573,223 +535,44 @@ Tree::Split( int nId)
    return Split( nId, dValue) ;
 }
 
-//----------------------------------------------------------------------------
-bool
-Tree::BuildTree_test( double dLinTol, double dSideMin, double dSideMax)
-{
-   // per poter usare questa funzione, anziché quella normale, bisogna:
-   // - commentare la parte di funzione di SetSurf dove si fanno gli split preliminare
-   // - se si usa anche la funzione GetLeaves, bisogna anche lì usare BuildTree_test al posto di BuildTree
-
-   //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) ;
-   // riempio anche la lista dei parent delle celle
-   m_vnParents = m_vnLeaves ;
-
-   return true ;
-}
-
 //----------------------------------------------------------------------------
 bool
 Tree::BuildTree( double dLinTol, double dSideMin, double dSideMax)
 {
-   // suddivido lo spazio parametrico con divisioni a metà su uno dei due parametri
+  // suddivido lo spazio parametrico con divisioni a metà su uno dei due parametri
    int nCToSplit = -1 ;
    Cell* pcToSplit = &m_mTree[nCToSplit] ;
    bool bIsPlanar = m_pSrfBz->IsPlanar() ;
    if ( ! m_bBilinear) {
-      while ( nCToSplit != -2  &&  pcToSplit->IsProcessed() == false) {
-         // controllo che la cella non sia già stata preliminarmente splittata
+      while ( nCToSplit != -2  &&  ! pcToSplit->IsProcessed()) {
+
+        // se la pezza non è già stata preliminarmente splittata
          if ( pcToSplit->IsLeaf()) {
-            // calcolo in quale direzione ho più curvatura
-            // ptP00P10 è un punto tra P00 e P10
-            double dCurvU = 0, dCurvV = 0 ;
+
+           // dimensioni parametriche della pezza
             double dLenParU = ( pcToSplit->GetTopRight().x - pcToSplit->GetBottomLeft().x) / SBZ_TREG_COEFF ;
             double dLenParV = ( pcToSplit->GetTopRight().y - pcToSplit->GetBottomLeft().y) / SBZ_TREG_COEFF ;
-            Point3d ptP00, ptP10, ptP11, ptP01 ;
-            // i vertici della cella
-            GetPoint(pcToSplit->GetBottomLeft().x, pcToSplit->GetBottomLeft().y, ptP00) ;
-            GetPoint(pcToSplit->GetTopRight().x, pcToSplit->GetBottomLeft().y, ptP10) ;
-            GetPoint(pcToSplit->GetTopRight().x, pcToSplit->GetTopRight().y, ptP11) ;
-            GetPoint(pcToSplit->GetBottomLeft().x, pcToSplit->GetTopRight().y, ptP01) ;
-            if ( dLenParU <= 1. / m_nDegV  ||  dLenParV <= 1. / m_nDegU  ||  Dist(ptP00, ptP11) < dSideMin * 2
-                 ||  Dist(ptP10, ptP01) < dSideMin * 2) {
-               double dU = ( pcToSplit->GetTopRight().x + pcToSplit->GetBottomLeft().x) / 2 ;
-               double dV = ( pcToSplit->GetTopRight().y + pcToSplit->GetBottomLeft().y) / 2 ;
-               double dULoc = 0.5, dVLoc = 0.5 ;
-               Point3d ptPSrf, ptP00P10, ptP10P11, ptP11P01, ptP01P00 ;
-               GetPoint( dU, dV, ptPSrf) ;
-               GetPoint( dU, pcToSplit->GetBottomLeft().y, ptP00P10) ;
-               GetPoint( pcToSplit->GetTopRight().x, dV, ptP10P11) ;
-               GetPoint( dU, pcToSplit->GetTopRight().y, ptP11P01) ;
-               GetPoint( pcToSplit->GetBottomLeft().x, dV, 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 * 2) ;
-               for ( double k = dStep ; k < 1 + EPS_SMALL ; k = k + dStep) {
-                  double dU = ( k * pcToSplit->GetTopRight().x + ( 1 - k) * pcToSplit->GetBottomLeft().x) ;
-                  double dV = ( pcToSplit->GetTopRight().y + pcToSplit->GetBottomLeft().y) / 2 ;
-                  GetPoint( dU, dV, ptPSrf) ;
-                  if ( k == 0.5)
-                     ptPSrfMid = ptPSrf ;
-                  GetPoint( dU, pcToSplit->GetBottomLeft().y, ptP00P10) ;
-                  GetPoint( dU, pcToSplit->GetTopRight().y, ptP11P01) ;
-                  CurveLine clV ;
-                  clV.Set( ptP00P10, ptP11P01) ;
-                  DistPointCurve dpc( ptPSrf, clV) ;
-                  double dDist ;
-                  dpc.GetDist( dDist) ;
-                  dCurvV = max( dCurvV, dDist) ;
-               }
-               dStep = 1. / ( m_nDegV * 2) ;
-               for ( double k = dStep ; k < 1 + EPS_SMALL ; k = k + dStep) {
-                  double dU = ( pcToSplit->GetTopRight().x + pcToSplit->GetBottomLeft().x) / 2 ;
-                  double dV = ( k * pcToSplit->GetTopRight().y + ( 1 - k) * pcToSplit->GetBottomLeft().y) ;
-                  if ( k == 0.5  &&  ! AreSamePointApprox( ORIG, ptPSrfMid))
-                     ptPSrf = ptPSrfMid ;
-                  else
-                     GetPoint( dU, dV, ptPSrf) ;
-                  GetPoint( pcToSplit->GetTopRight().x, dV, ptP10P11) ;
-                  GetPoint( pcToSplit->GetBottomLeft().x, dV, ptP01P00) ;
-                  CurveLine clU ;
-                  clU.Set( ptP01P00, ptP10P11) ;
-                  DistPointCurve dpc( ptPSrf, clU) ;
-                  double dDist ;
-                  dpc.GetDist( dDist) ;
-                  dCurvU = max( dCurvU, dDist) ;
-               }
-            }
 
-            // devo calcolare anche il twist, in caso di bordi rettilinei ( superficie di grado maggiore di 1, ma che in realtà è una bilineare)
-            // NON posso guardare la distanza tra il punto medio delle diagonali e il punto centrale della cella ( uLoc = 0.5, vLoc = 0.5)
-            // posso guardare la distanza tra le due diagonali
-            bool bTwist = false ;
-            // serve una valutazione più fine, sennò approssimo la superficie in modo troppo grossolano
-            DistLineLine dll( ptP00, ptP11, ptP10, ptP01, true, true) ;
-            double dDist = 0 ; dll.GetDist( dDist) ;
-            if ( dDist > max(dCurvU, dCurvV)  ||  dDist < 5 * EPS_SMALL) {
-               bool bFlat = false ;
-               // controllo se la cella è twistata di 180 gradi e quindi piatta
-               Triangle3d tria1, tria2 ;
-               tria1.Set( ptP00, ptP10, ptP11) ; tria1.Validate( true) ;
-               tria2.Set( ptP00, ptP11, ptP01) ; tria2.Validate( true) ;
-               if( AreOppositeVectorEpsilon(tria1.GetN(), tria2.GetN(), 5 * EPS_SMALL)) {
-                  bTwist = true ;
-                  bFlat = true ;
-               }
-               // controllo che la cella non sia piatta
-               if( ! bTwist) {
-                  PolyLine plCell ;
-                  plCell.AddUPoint(0,ptP00) ;
-                  plCell.AddUPoint(1,ptP10) ;
-                  plCell.AddUPoint(2,ptP11) ;
-                  double dU = (pcToSplit->GetTopRight().x + pcToSplit->GetBottomLeft().x) / 2 ;
-                  double dV = (pcToSplit->GetTopRight().y + pcToSplit->GetBottomLeft().y) / 2 ;
-                  Point3d ptCen ; GetPoint( dU, dV, ptCen) ;
-                  plCell.AddUPoint(3,ptCen) ;
-                  plCell.AddUPoint(4,ptP01) ;
-                  plCell.Close() ;
-                  Plane3d plPlane ; double dArea = 0 ;
-                  bFlat = plCell.IsClosedAndFlat( plPlane, dArea, 10 * EPS_SMALL) ;
-               }
-               if( ! bFlat  &&  dDist > 5 * EPS_SMALL) {
-                  bTwist = true ;
-                  // devo decidere in quale direzione splittare
-                  // dovrei capire in quale delle due direzioni è più torta la superficie
-                  Vector3d vtU0 = ptP10 - ptP00 ;
-                  Vector3d vtU1 = ptP11 - ptP01 ;
-                  double dAngU ;
-                  bool bDetU = false ;
-                  vtU0.GetRotation( vtU1, vtU0 ^ vtU1, dAngU, bDetU) ;
-                  Vector3d vtV0 = ptP01 - ptP00 ;
-                  Vector3d vtV1 = ptP11 - ptP10 ;
-                  double dAngV ;
-                  bool bDetV = false ;
-                  // faccio la get rotation tra le coppie vettori, usando come asse il loro prodotto vettoriale per ottenere la rotazione tra i due lati
-                  // splitto nella direzione perpendicolare alla coppia di vettori più torti tra loro.
-                  vtV0.GetRotation( vtV1, vtV0 ^ vtV1, dAngV, bDetV) ;
-                  if ( dAngU > dAngV)
-                     dCurvV = dDist ;
-                  else
-                     dCurvU = dDist ;
-               }
-            }
+           // vertici e centro della pezza
+            Point3d ptP00, ptP10, ptP11, ptP01, ptCen ;
+            GetPoint( pcToSplit->GetBottomLeft().x, pcToSplit->GetBottomLeft().y, ptP00) ;
+            GetPoint( pcToSplit->GetTopRight().x, pcToSplit->GetBottomLeft().y, ptP10) ;
+            GetPoint( pcToSplit->GetTopRight().x, pcToSplit->GetTopRight().y, ptP11) ;
+            GetPoint( pcToSplit->GetBottomLeft().x, pcToSplit->GetTopRight().y, ptP01) ;
+            GetPoint( pcToSplit->GetCenter().x, pcToSplit->GetCenter().y, ptCen) ;
+            //{ string sLog = "P00=" + ToString( ptP00, 3) + " P10=" + ToString( ptP10, 3) +
+            //                " P11=" + ToString( ptP11, 3) + " P01=" + ToString( ptP01, 3) +
+            //                " Cen=" + ToString( ptCen, 3);
+            //  LOG_DBG_INFO( GetEGkLogger(), sLog.c_str())}
 
-            // 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 ;
-            }
-            pcToSplit->SetSplitDirVert( bVert) ;
-            // distanza reale tra i vertici della cella
+           // distanze sul contorno tra i quattro vertici
             double dLen0 = Dist( ptP00, ptP10) ;
             double dLen1 = Dist( ptP10, ptP11) ;
             double dLen2 = Dist( ptP01, ptP11) ;
             double dLen3 = Dist( ptP00, ptP01) ;
+           // distanza reale tra i vertici della pezza
             if ( dLen0 < EPS_ZERO  &&  dLen2 < EPS_ZERO ) {
-               double dV = ( pcToSplit->GetBottomLeft().y + pcToSplit->GetTopRight().y) / 2 / SBZ_TREG_COEFF ;
+               double dV = pcToSplit->GetCenter().y / SBZ_TREG_COEFF ;
                PtrOwner<CurveComposite> pCrvV( m_pSrfBz->GetCurveOnU( dV)) ;
                double dLenU0, dLenU1 ;
                pCrvV->GetLengthAtParam( pcToSplit->GetBottomLeft().x / SBZ_TREG_COEFF, dLenU0) ;
@@ -797,126 +580,177 @@ Tree::BuildTree( double dLinTol, double dSideMin, double dSideMax)
                dLen0 = abs( dLenU1 - dLenU0) ;
             }
             if ( dLen1 < EPS_ZERO  &&  dLen3 < EPS_ZERO ) {
-               double dU = ( pcToSplit->GetBottomLeft().x + pcToSplit->GetTopRight().x) / 2 / SBZ_TREG_COEFF ;
+               double dU = pcToSplit->GetCenter().x / SBZ_TREG_COEFF ;
                PtrOwner<CurveComposite> pCrvU( m_pSrfBz->GetCurveOnV( dU)) ;
                double dLenV0, dLenV1 ;
                pCrvU->GetLengthAtParam( pcToSplit->GetBottomLeft().y / SBZ_TREG_COEFF, dLenV0) ;
                pCrvU->GetLengthAtParam( pcToSplit->GetTopRight().y / SBZ_TREG_COEFF, dLenV1) ;
                dLen1 = abs( dLenV1 - dLenV0) ;
             }
+
+           // calcolo massimo della freccia (sagitta)
+            double dSagU = 0, dSagV = 0 ;
+
+           // da distanza tra centro e triangoli di approssimazione
+            if ( dSagU < dLinTol  &&  dSagV < dLinTol) {
+              // distanza del centro dai due triangoli con la prima diagonale
+               Triangle3d Tria1A, Tria1B ;
+               Tria1A.Set( ptP00, ptP10, ptP11) ; Tria1A.Validate( true) ;
+               Tria1B.Set( ptP00, ptP11, ptP01) ; Tria1B.Validate( true) ;
+               double dDist1A = NAN ;
+               DistPointTriangle( ptCen, Tria1A).GetDist( dDist1A) ;
+               double dDist1B = NAN ;
+               DistPointTriangle( ptCen, Tria1B).GetDist( dDist1B) ;
+               double dDist1 = min( dDist1A, dDist1B) ;
+              // distanza del centro dai due triangoli con la seconda diagonale
+               Triangle3d Tria2A, Tria2B ;
+               Tria2A.Set( ptP10, ptP11, ptP01) ; Tria2A.Validate( true) ;
+               Tria2B.Set( ptP10, ptP01, ptP00) ; Tria2B.Validate( true) ;
+               double dDist2A = NAN ;
+               DistPointTriangle( ptCen, Tria2A).GetDist( dDist2A) ;
+               double dDist2B = NAN ;
+               DistPointTriangle( ptCen, Tria2B).GetDist( dDist2B) ;
+               double dDist2 = min( dDist2A, dDist2B) ;
+              // prendo la minore delle due distanze
+               double dDist = min( dDist1, dDist2) ;
+              // se maggiore della tolleranza, imposto come freccia
+               if ( isfinite( dDist)  &&  dDist > dLinTol) {
+                 // divido in base alle distanze tra i vertici dei lati
+                  if ( max( dLen0, dLen2) >= max( dLen1, dLen3))
+                     dSagU = dDist ;
+                  else
+                     dSagV = dDist ;
+                  //{ string sLog = " Da triangoli : FrecciaU=" + ToString( dSagU, 3) + " FrecciaV=" + ToString( dSagV, 3) ;
+                  //  LOG_DBG_INFO( GetEGkLogger(), sLog.c_str())}
+               }
+            }
+
+           // su isoparametriche in U e V
+            if ( dSagU < dLinTol  &&  dSagV < dLinTol) {
+              // step di verifica in U e in V
+               int nStepU = ( dLenParU > 1. / m_nDegU ? m_nDegU + 1 : 2) ;
+               int nStepV = ( dLenParV > 1. / m_nDegV ? m_nDegV + 1 : 2) ;
+              // verifico lungo curve in U
+               for ( int i = 0 ; i <= nStepU ; ++ i) {
+                 // parametro U in analisi
+                  double dCoeffU = double( i) / nStepU ;
+                  double dU = ( 1 - dCoeffU) * pcToSplit->GetBottomLeft().x + dCoeffU * pcToSplit->GetTopRight().x ;
+                 // linea tra gli estremi a questa U (ptP00P10 è un punto tra P00 e P10 e così via)
+                  Point3d ptP00P10 ;
+                  GetPoint( dU, pcToSplit->GetBottomLeft().y, ptP00P10) ;
+                  Point3d ptP11P01 ;
+                  GetPoint( dU, pcToSplit->GetTopRight().y, ptP11P01) ;
+                  CurveLine clV ;
+                  clV.Set( ptP00P10, ptP11P01) ;
+                 // verifico alcuni punti intermedi in V
+                  for ( int j = 1 ; j < nStepV ; ++ j) {
+                    // parametro in V
+                     double dCoeffV = double( j) / nStepV ;
+                     double dV = ( 1 - dCoeffV) * pcToSplit->GetBottomLeft().y + dCoeffV * pcToSplit->GetTopRight().y ;
+                    // punto a U, V
+                     Point3d ptPSrf ;
+                     GetPoint( dU, dV, ptPSrf) ;
+                    // distanza del punto dalla precedente linea
+                     DistPointCurve dpc( ptPSrf, clV) ;
+                     double dDist ;
+                     dpc.GetDist( dDist) ;
+                     dSagV = max( dSagV, dDist) ;
+                  }
+               }
+              // verifico lungo curve in V
+               for ( int i = 0 ; i <= nStepV ; ++ i) {
+                 // parametro in V in analisi
+                  double dCoeffV = double( i) / nStepV ;
+                  double dV = ( 1 - dCoeffV) * pcToSplit->GetBottomLeft().y + dCoeffV * pcToSplit->GetTopRight().y ;
+                 // linea tra gli estremi a questa V
+                  Point3d ptP10P11 ;
+                  GetPoint( pcToSplit->GetTopRight().x, dV, ptP10P11) ;
+                  Point3d ptP01P00 ;
+                  GetPoint( pcToSplit->GetBottomLeft().x, dV, ptP01P00) ;
+                  CurveLine clU ;
+                  clU.Set( ptP01P00, ptP10P11) ;
+                 // verifico alcuni punti intermedi in U
+                  for ( int j = 1 ; j < nStepU ; ++ j) {
+                    // parametro U
+                     double dCoeffU = double( j) / nStepU ;
+                     double dU = ( 1 - dCoeffU) * pcToSplit->GetBottomLeft().x + dCoeffU * pcToSplit->GetTopRight().x ;
+                    // punto a U, V
+                     Point3d ptPSrf ;
+                     GetPoint( dU, dV, ptPSrf) ;
+                    // distanza del punto dalla precedente linea
+                     DistPointCurve dpc( ptPSrf, clU) ;
+                     double dDist ;
+                     dpc.GetDist( dDist) ;
+                     dSagU = max( dSagU, dDist) ;
+                  }
+               }
+               //{ string sLog = " Da Isoparam : FrecciaU=" + ToString( dSagU, 3) + " FrecciaV=" + ToString( dSagV, 3) ;
+               //  LOG_DBG_INFO( GetEGkLogger(), sLog.c_str())}
+            }
+
+           // per lo split scelgo la direzione che è più vicina alla superficie originale nel punto di maggior distanza
+           // misura approssimativa della curvatura in una direzione
+            bool bVert ;
+            if ( max( dLen0, dLen2) > 5 * max( dLen1, dLen3))
+               bVert = true ;
+            else if ( max( dLen1, dLen3) > 5 * max( dLen0, dLen2))
+               bVert = false ;
+            else
+               bVert = ( dSagV <= dSagU) ;
+            pcToSplit->SetSplitDirVert( bVert) ;
             // verifico che la cella sia da splittare e che eventualmente sia abbastanza grande da poterlo fare
-            double dSideMinVal = 0, dSideMaxVal = 0 ;
+            double dSideMinVal = 0 ;
+            double dLengMinVal = 0 ;
             if ( bVert) {
                if ( dLen0 > EPS_ZERO  &&  dLen2 > EPS_ZERO)
                   dSideMinVal = min( dLen0, dLen2) ;
                else
                   dSideMinVal = max( dLen0, dLen2) ;
+               if ( dLen1 > EPS_ZERO  &&  dLen3 > EPS_ZERO)
+                  dLengMinVal = min( dLen1, dLen3) ;
+               else
+                  dLengMinVal = max( dLen1, dLen3) ;
             }
             else {
                if ( dLen1 > EPS_ZERO  &&  dLen3 > EPS_ZERO)
                   dSideMinVal = min( dLen1, dLen3) ;
                else
                   dSideMinVal = max( dLen1, dLen3) ;
+               if ( dLen0 > EPS_ZERO  &&  dLen2 > EPS_ZERO)
+                  dLengMinVal = min( dLen0, dLen2) ;
+               else
+                  dLengMinVal = max( dLen0, dLen2) ;
             }
-            // calcolo le diagonali per controllare la dimensione massima dei triangoli in cui dividerei la cella
-            dSideMaxVal = max( Dist( ptP00, ptP11), Dist( ptP10, ptP01)) ;
+           // calcolo le diagonali per controllare la dimensione massima dei triangoli in cui dividerei la cella
+            double dSideMaxVal = max( Dist( ptP00, ptP11), Dist( ptP10, ptP01)) ;
          
-            // se la cella è abbastanza grande da poter essere divisa ancora, calcolo l'errore di approssimazione
+           // se la cella è abbastanza grande da poter essere divisa ancora, calcolo l'errore di approssimazione
+           // ( dSideMinVal è zero se entrambi i lati da splittare sono collassati in un punto, controllo dLengMinVal)
             bool bSplit = false ;
-            // dSideMinVal potrebbe essere zero se entrambi i lati che dovrei splittare sono collassati in un punto, ma questo non vuol
-            // dire che non dovrei eseguire lo split
-            if ( ! bTwist  &&  (dSideMinVal / 2 >= dSideMin  ||  dSideMinVal < EPS_SMALL)  &&  dSideMaxVal < dSideMax  &&  ( dCurvV > dLinTol  ||  dCurvU > dLinTol)) {
-               CurveLine cl0010, cl0001, cl1011, cl0111 ;
-               // V=0
-               cl0010.Set( ptP00, ptP10) ;
-               // V=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) * pcToSplit->GetBottomLeft().x + dU * pcToSplit->GetTopRight().x) ;
-                  if ( ! GetPoint( dULoc, pcToSplit->GetBottomLeft().y, ptBz0)  ||
-                        ! GetPoint( dULoc, pcToSplit->GetTopRight().y, ptBz1))
-                     return false ;
-                  // verifico che la cella non sia uno spicchio in verticale, cioè con ptP00 == ptP01 && ptP10 == ptP11
-                  // ( vedi disegno sotto per uno spicchio verticale)
-                  // sennò i punti che cerco sono semplicemente i vertici
-                  if ( cl0010.IsValid()) {
-                     DistPointCurve dpc0010( ptBz0, cl0010) ;
-                     dpc0010.GetMinDistPoint( 0, pt0010, nFlag) ;
-                  }
-                  else
-                     pt0010 = ptP00 ;
-                  if ( cl0111.IsValid()) {
-                     DistPointCurve dpc0111( ptBz1, cl0111) ;
-                     dpc0111.GetMinDistPoint( 0, pt0111, nFlag) ;
-                  }
-                  else
-                     pt0111 = ptP01 ;
-                  // curva a parametro U fisso, con V che scorre
-                  clV.Set( pt0010, pt0111) ;
-                  for ( int v = 0 ; v < nStepsV ; ++ v) {
-                     double dV = double ( v) / double ( nStepsV - 1) ;
-                     double dVLoc = ( ( 1 - dV) * pcToSplit->GetBottomLeft().y + dV * pcToSplit->GetTopRight().y) ;
-                     if ( ! GetPoint( dULoc, dVLoc, ptBzV))
-                        return false ;
-                     DistPointCurve dpc( ptBzV, clV) ;
-                     // distanza di approssimazione locale
-                     double dDist ; 
-                     dpc.GetDist( dDist) ;
-                     // se la cella è uno spicchio, quindi con due lati collassati, devo calcolare in modo diverso dist
-                    //     ptP00 == ptP01
-                    //         / \
-                    //        /   \
-                    //       /     \
-                    //      (       )
-                    //       \     /
-                    //        \   /
-                    //         \ /
-                    //     ptP10 == ptP11
-                     if ( ! clV.IsValid()  &&  AreSamePointApprox( ptP00, ptP01)  &&  AreSamePointApprox( ptP10, ptP11)) {
-                        DistPointCurve dpcSlice( ptBzV, cl0010) ;
-                        dpcSlice.GetDist( dDist) ;
-                     }
-                     if ( dDist > dLinTol) {
-                        bSplit = true ;
-                        break ;
-                     }
-                  }
-               }
+            bool bDimOk = ( dSideMinVal / 2 >= dSideMin  ||  ( dSideMinVal < EPS_SMALL  &&  dLengMinVal / 2 >= dSideMin)) ;
+            if ( dSideMaxVal > dSideMax) {
+               bSplit = true ;
+               //LOG_DBG_INFO( GetEGkLogger(), "  Split by SideMax")
             }
-            else if ( bTwist ) {
-               // se la cella è twistata allora l'errore lo calcolo come nelle bilineari
-               //double dErr = 1. / 4. * ( ( ptP00 - ptP01) + ( ptP11 - ptP10)).Len() ;
-               double dErr = ( ( ptP00 - ptP01) + ( ptP11 - ptP10)).Len() / 20 ;
-               if ( dErr > dLinTol)
-                  bSplit = true ;
+            else if ( dSagV > dLinTol  ||  dSagU > dLinTol) {
+               bSplit = bDimOk ;
+               //if ( bSplit)
+               //   LOG_DBG_INFO( GetEGkLogger(), "  Split by SagittaUV")
             }
 
-            if ( bSplit  ||  dSideMaxVal > dSideMax) {
+            if ( bSplit) {
                pcToSplit->SetSplitDirVert( bVert) ;
-               // effettuo lo split
+              // effettuo lo split
                Split( nCToSplit) ;
-
-               // procedo con lo split del Child1
+              // procedo con lo split del Child1
                nCToSplit = pcToSplit->m_nChild1 ;
                pcToSplit = &m_mTree[nCToSplit] ;
             }
             else {
-               // sono arrivato ad una cella Leaf, quindi salvo la cella
+              // sono arrivato ad una cella Leaf, quindi salvo la cella
                m_vnLeaves.push_back( nCToSplit) ;
                pcToSplit->SetProcessed() ;
-               // risalgo i parent finché non trovo il primo Child2 da processare
+              // risalgo i parent finché non trovo il primo Child2 da processare
                nCToSplit = pcToSplit->m_nParent ;
                pcToSplit = &m_mTree[nCToSplit] ;
                if ( nCToSplit == -2)
@@ -942,26 +776,26 @@ Tree::BuildTree( double dLinTol, double dSideMin, double dSideMax)
             pcToSplit = &m_mTree[nCToSplit] ;
          }
       }
-      Balance() ; // da implementare quando dividerò ad un parametro a scelta e non a metà  // probabilmente mi servirà salvare nella cella il livello di profondità
    }
-   // bilineare
+
+  // bilineare
    else {
       while ( nCToSplit != -2  &&  pcToSplit->IsProcessed() == false) {
          if ( pcToSplit->IsLeaf()) {
-            // vertici della cella
+           // vertici della cella
             Point3d ptP00, ptP10, ptP11, ptP01 ;
             GetPoint(pcToSplit->GetBottomLeft().x, pcToSplit->GetBottomLeft().y, ptP00) ;
             GetPoint(pcToSplit->GetTopRight().x, pcToSplit->GetBottomLeft().y, ptP10) ;
             GetPoint(pcToSplit->GetTopRight().x, pcToSplit->GetTopRight().y, ptP11) ;
             GetPoint(pcToSplit->GetBottomLeft().x, pcToSplit->GetTopRight().y, ptP01) ;
-            // distanza reale tra i vertici della cella
+           // distanza reale tra i vertici della cella
             double dLen0 = Dist( ptP00, ptP10) ;
             double dLen1 = Dist( ptP10, ptP11) ;
             double dLen2 = Dist( ptP01, ptP11) ;
             double dLen3 = Dist( ptP00, ptP01) ;
                
             bool bVert = false ;
-            ///per capire in quale direzione splittare devo guardare quale coppia di lati opposti è più sghemba
+           // per capire in quale direzione splittare devo guardare quale coppia di lati opposti è più sghemba
             Vector3d vtU0 = ptP01 - ptP00 ;
             Vector3d vtU1 = ptP11 - ptP10 ;
             double dLU0, dFU0, dTU0, dLU1, dFU1, dTU1 ;
@@ -974,7 +808,7 @@ Tree::BuildTree( double dLinTol, double dSideMin, double dSideMax)
             vtV0.ToSpherical( &dLV0, &dFV0, &dTV0) ;
             vtV1.ToSpherical( &dLV1, &dFV1, &dTV1) ;
             double dSkewnessU = abs( dFV1 - dFV0) + abs( dTV1 - dTV0) ;
-            if( dSkewnessU > dSkewnessV)
+            if ( dSkewnessU > dSkewnessV)
                bVert = false ;
             else
                bVert = true ;
@@ -1005,25 +839,24 @@ Tree::BuildTree( double dLinTol, double dSideMin, double dSideMax)
                   }
                }
             }
-            else if ( ! bIsPlanar){
-               dErr = 1. / 4. * ( ( ptP00 - ptP01) + ( ptP11 - ptP10)).Len() ;
-               //int dErr2 = 1. / 4. * ( ( ptP10 - ptP01) + ( ptP11 - ptP00)).Len() ;  //correzione che mi verrebbe intuitiva, ma che fa dividere la superficie molto di più ( probabilmente troppo)!! quindi probabilmente sbagliata
+            else if ( ! bIsPlanar) {
+               dErr = ( ( ptP00 - ptP01) + ( ptP11 - ptP10)).Len() / 4 ;
             }
-            // se la cella è abbastanza grande da poter essere divisa ancora e devo approssimare meglio, la divido
+           // se la cella è abbastanza grande da poter essere divisa ancora e devo approssimare meglio, la divido
             if ( dSideMinVal / 2 >= dSideMin  &&  dSideMaxVal < dSideMax  &&  dErr > dLinTol) {
                pcToSplit->SetSplitDirVert( bVert) ;
-               // effettuo lo split
+              // effettuo lo split
                Split( nCToSplit) ;
 
-               // procedo con lo split del Child1
+              // procedo con lo split del Child1
                nCToSplit = pcToSplit->m_nChild1 ;
                pcToSplit = &m_mTree[nCToSplit] ;
             }
             else {
-               // sono arrivato ad una cella Leaf, quindi salvo la cella
+              // sono arrivato ad una cella Leaf, quindi salvo la cella
                m_vnLeaves.push_back( nCToSplit) ;
                pcToSplit->SetProcessed() ;
-               // risalgo i parent finché non trovo il primo Child2 da processare
+              // risalgo i parent finché non trovo il primo Child2 da processare
                nCToSplit = pcToSplit->m_nParent ;
                pcToSplit = &m_mTree[nCToSplit] ;
                if ( nCToSplit == -2)
@@ -1050,22 +883,10 @@ Tree::BuildTree( double dLinTol, double dSideMin, double dSideMax)
          }
       }
    }
+
    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
@@ -1546,8 +1367,8 @@ Tree::GetPolygons( POLYLINEMATRIX& vvPolygons, POLYLINEMATRIX& vvPolygons3d, vec
          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()) ;
+            // vettore in cui salvo i loop che non appartengono al poligono che sto costruendo nel ciclo attuale e da cui ripasserò dopo
+            INTVECTOR vToCheck( m_mTree[nId].m_vInters.size()) ;
             iota( vToCheck.begin(), vToCheck.end(), 0) ;
             // numero di poligoni aggiunti
             int nPoly = 0 ;
@@ -1942,23 +1763,20 @@ Tree::FindCell( const Point3d& ptToAssign, const CurveLine& clTrim, bool bRecurs
    // se fallisce ritorna un vettore vuoto
    // verifico che il punto sia all'interno dello spazio parametrico
    // allargo i bordi in modo da tenere anche i punti sul bordo dello spazio parametrico
-   if ( ptToAssign.x < m_mTree.at( -1).GetBottomLeft().x - EPS_SMALL  ||  ptToAssign.x > m_mTree.at( -1).GetTopRight().x  + EPS_SMALL||
-         ptToAssign.y < m_mTree.at( -1).GetBottomLeft().y - EPS_SMALL  ||  ptToAssign.y > m_mTree.at( -1).GetTopRight().y + EPS_SMALL) {
-      //nCells.push_back( - 2) ;
+   if ( ptToAssign.x < m_mTree.at( -1).GetBottomLeft().x - EPS_SMALL  ||
+        ptToAssign.x > m_mTree.at( -1).GetTopRight().x  + EPS_SMALL  ||
+        ptToAssign.y < m_mTree.at( -1).GetBottomLeft().y - EPS_SMALL  ||
+        ptToAssign.y > m_mTree.at( -1).GetTopRight().y + EPS_SMALL) {
       return nCells ;
    }
    
    // se ho diviso preliminarmente le patches e in uno dei due parametri ho un numero dispari di patches devo individuare a mano la cella parent
    // in cui individuare la foglia giusta
-   if (  (m_bSplitPatches  &&  ( m_nSpanU > 1  ||  m_nSpanV > 1))  || m_bTestMode) {
+   if ( m_nSpanU > 1  ||  m_nSpanV > 1) {
       INTVECTOR nParents = FindCell( ptToAssign, clTrim, m_vnParents) ;
       if ( nParents.empty())
          return nCells ;
       nId = nParents.back() ;
-      if ( m_bTestMode ) {
-         nCells.push_back(nId) ;
-         return nCells ;
-      }
    }
    // individuo la foglia in cui ho lo start del loop
    while ( ! m_mTree.at( nId).IsLeaf()) {
@@ -1988,8 +1806,7 @@ Tree::FindCell( const Point3d& ptToAssign, const CurveLine& clTrim, bool bRecurs
       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 ( abs( ptToAssign.x - ptTl.x) < EPS_SMALL  ||  abs( ptToAssign.x - ptBr.x) < EPS_SMALL  ||
-           abs( ptToAssign.y - ptTl.y) < EPS_SMALL  ||  abs( ptToAssign.y - ptBr.y) < EPS_SMALL)
-      {
+           abs( ptToAssign.y - ptTl.y) < EPS_SMALL  ||  abs( ptToAssign.y - ptBr.y) < EPS_SMALL) {
          Vector3d vtDir ;
          clTrim.GetStartDir( vtDir) ;
          // proseguo lungo la curva di trim di EPS_SMALL
@@ -1999,8 +1816,10 @@ Tree::FindCell( const Point3d& ptToAssign, const CurveLine& clTrim, bool bRecurs
             Vector3d vtDirDX = vtDir ; vtDirDX.Rotate( Z_AX, 90) ;
             ptToAssignPlus = ptToAssignPlus + vtDir * EPS_SMALL ;
             // controllo di non essere uscito dallo spazio parametrico ed eventualmente giro a sinistra
-            if ( ptToAssignPlus.x < m_mTree.at( -1).GetBottomLeft().x - EPS_SMALL  ||  ptToAssignPlus.x > m_mTree.at( -1).GetTopRight().x  + EPS_SMALL||
-               ptToAssignPlus.y < m_mTree.at( -1).GetBottomLeft().y - EPS_SMALL  ||  ptToAssignPlus.y > m_mTree.at( -1).GetTopRight().y + EPS_SMALL) {
+            if ( ptToAssignPlus.x < m_mTree.at( -1).GetBottomLeft().x - EPS_SMALL  ||
+                 ptToAssignPlus.x > m_mTree.at( -1).GetTopRight().x  + EPS_SMALL  ||
+                 ptToAssignPlus.y < m_mTree.at( -1).GetBottomLeft().y - EPS_SMALL  ||
+                 ptToAssignPlus.y > m_mTree.at( -1).GetTopRight().y + EPS_SMALL) {
                // rispetto al punto di partenza avanzo lungo la curva di trim
                ptToAssignPlus = ptToAssign + vtDir * EPS_SMALL ;
                Vector3d vtDirSX = vtDir ; vtDirSX.Rotate( Z_AX, -90) ;
@@ -4340,7 +4159,7 @@ Tree::CreateCellContour( POLYLINEMATRIX& vPolygons)
    int nRoot = -1 ;
    // preparo tutto per poter chiamare la createCellPolygon
    m_vnLeaves.push_back( nRoot) ;
-   INTVECTOR vToCheck( (int) m_mTree.at(nRoot).m_vInters.size()) ;
+   INTVECTOR vToCheck( m_mTree.at(nRoot).m_vInters.size()) ;
    iota( vToCheck.begin(),  vToCheck.end(), 0) ;
    int nPoly = 0 ;
    INTVECTOR vnParentChunk ;

Tree.h

@@ -28,7 +28,7 @@ struct PairHashInt64 {
       size_t h2 = std::hash<int64_t>{}(key.second) ;
       return h1 ^ (h2 << 1);  // Combine hashes
    }
-};
+} ;
 
 //----------------------------------------------------------------------------
 struct Inters {
@@ -42,13 +42,13 @@ struct Inters {
    // se ho più intersezioni che entrano in un lato le riordino considerando che percorro i lati in senso antiorario a partire da ptTR
    bool operator < ( Inters& b)
    {
-      // trovo in che ordine stanno i due start, tenendo conto anche della possibilità che siano vertici
+     // trovo in che ordine stanno i due start, tenendo conto anche della possibilità che siano vertici
       INTVECTOR vEdges = { 7, 0, 4, 1, 5, 2, 6, 3} ;
       const auto iter1 = find( vEdges.begin(), vEdges.end(), nIn) ;
       int nPos1 = std::distance( vEdges.begin(), iter1) ;
       const auto iter2 = find( vEdges.begin(), vEdges.end(), b.nIn) ;
       int nPos2 = std::distance( vEdges.begin(), iter2) ;
-      // se sono loop interni li ordino in modo decrescente rispetto all'area
+     // se sono loop interni li ordino in modo decrescente rispetto all'area
       bool bEqIn = ( nIn == b.nIn) ;
       double dAreaA = 0 , dAreaB = 0 ;
       if ( bEqIn  &&  nIn == -1) {
@@ -63,7 +63,7 @@ struct Inters {
          pl.Close() ;
          pl.GetAreaXY( dAreaB) ;
       }
-      // se nIn è un vertice sistemo il valore
+     // se nIn è un vertice sistemo il valore
       int nEdgeIn = nIn ;
       if ( nIn > 3)
          nEdgeIn = nIn - 4 ;
@@ -75,13 +75,13 @@ struct Inters {
                ( bEqIn  &&  nEdgeIn == 3  &&  vpt[0].y < b.vpt[0].y)) ;
    }
 
-   static bool FirstEncounter ( Inters& a, Inters& b)
+   static bool FirstEncounter( Inters& a, Inters& b)
    {
-      // riordino in base al lato toccato, o dall'uscita o  dall'ingresso, che viene prima.
-      // ottengo l'ordine che avrei percorrendo il bordo da ptTR e considerando i loop che incontro, indipendentemente se li incontro nel punto di uscita o ingresso
-      // nell'intersezione salvo se il taglio è stato ordinato guardando l'ingresso o l'uscita
+     // riordino in base al lato toccato, o dall'uscita o  dall'ingresso, che viene prima.
+     // ottengo l'ordine che avrei percorrendo il bordo da ptTR e considerando i loop che incontro, indipendentemente se li incontro nel punto di uscita o ingresso
+     // nell'intersezione salvo se il taglio è stato ordinato guardando l'ingresso o l'uscita
       INTVECTOR vEdges = { 7, 0, 4, 1, 5, 2, 6, 3} ;
-      // trovo i lati di ingresso e uscita
+     // trovo i lati di ingresso e uscita
       const auto iter1 = find( vEdges.begin(), vEdges.end(), a.nIn) ;
       int nPos1 = std::distance( vEdges.begin(), iter1) ;
       const auto iter2 = find( vEdges.begin(), vEdges.end(), a.nOut) ;
@@ -92,13 +92,13 @@ struct Inters {
       int nPos4 = std::distance( vEdges.begin(), iter4) ;
       int nFirstA = 0 ;
       int nFirstB = 0 ;
-      // salvo l'indice del primo punto dell'intersezione che ho incontrato scorrendo il bordo da ptTR
-      // salvo il lato che viene prima confrontando ingresso e uscita
+     // salvo l'indice del primo punto dell'intersezione che ho incontrato scorrendo il bordo da ptTR
+     // salvo il lato che viene prima confrontando ingresso e uscita
       if ( nPos2 < nPos1) {
          nPos1 = nPos2 ;
          nFirstA = int( a.vpt.size()) - 1 ;
       }
-      // se ingresso e uscita sono sullo stesso lato allora confronto le coordinate per capire se viene prima l'ingresso o l'uscita
+     // se ingresso e uscita sono sullo stesso lato allora confronto le coordinate per capire se viene prima l'ingresso o l'uscita
       else if ( nPos2 == nPos1 ) {
          if ( nPos1 == 0 )
             nFirstA = a.vpt[0].x > a.vpt.back().x ? 0 : ( int( a.vpt.size()) - 1) ;
@@ -125,11 +125,11 @@ struct Inters {
       }
       a.bSortedbyStart = nFirstA == 0 ;
       b.bSortedbyStart = nFirstB == 0 ;
-      // se sono diversi ritorno il confronto
+     // se sono diversi ritorno il confronto
       if ( nPos1 != nPos3)
          return nPos1 < nPos3 ;
-      // se sono uguali devo valutare il punto di intersezione
-       return ( nPos1 == 0  &&  a.vpt[nFirstA].x > b.vpt[nFirstB].x)  ||
+     // se sono uguali devo valutare il punto di intersezione
+      return ( nPos1 == 0  &&  a.vpt[nFirstA].x > b.vpt[nFirstB].x)  ||
               ( nPos1 == 1  &&  a.vpt[nFirstA].y > b.vpt[nFirstB].y)  ||
               ( nPos1 == 2  &&  a.vpt[nFirstA].x < b.vpt[nFirstB].x)  ||
               ( nPos1 == 3  &&  a.vpt[nFirstA].y < b.vpt[nFirstB].y) ;
@@ -163,17 +163,15 @@ class Cell
    public :
       ~Cell( void) {}
       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_bLabelled( false), m_nRightEdgeIn( -1), m_bOnLeftEdge( false), m_bOnTopEdge( false), m_nVertToErase( -1),
-           m_ptPbl( ORIG), m_ptPtr( SBZ_TREG_COEFF, SBZ_TREG_COEFF, 0), m_bProcessed( false), m_bSplitVert( true)
-         {
-            Point3d ptTr ( 1 * SBZ_TREG_COEFF, 1 * SBZ_TREG_COEFF) ;
-            m_ptPtr = ptTr ;
-         }
+         : m_nId( -1), m_nTop( -2), m_nBottom( -2), m_nLeft( -2), m_nRight( -2), m_nParent( -2), m_nDepth( 0),
+           m_dSplit( 0), m_nChild1( -2), m_nChild2( -2), m_nFlag( -1), m_bLabelled( false), m_nRightEdgeIn( -1),
+           m_bOnLeftEdge( false), m_bOnTopEdge( false), m_nVertToErase( -1), m_ptPbl( ORIG),
+           m_ptPtr( SBZ_TREG_COEFF, SBZ_TREG_COEFF, 0), m_bProcessed( false), m_bSplitVert( true) {}
       Cell( const Point3d& ptBL, const 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_bLabelled( 0), m_nRightEdgeIn( -1), m_bOnLeftEdge( false), m_bOnTopEdge( false), m_nVertToErase( -1),
-           m_ptPbl( ptBL), m_ptPtr( ptTR), m_bProcessed( false), m_bSplitVert( true) {}
+         : m_nId( -1), m_nTop( -2), m_nBottom( -2), m_nLeft( -2), m_nRight( -2), m_nParent( -2), m_nDepth( 0),
+           m_dSplit( 0), m_nChild1( -2), m_nChild2( -2), m_nFlag( -1), m_bLabelled( 0), m_nRightEdgeIn( -1),
+           m_bOnLeftEdge( false), m_bOnTopEdge( false), m_nVertToErase( -1), m_ptPbl( ptBL),
+           m_ptPtr( ptTR), m_bProcessed( false), m_bSplitVert( true) {}
       bool IsSame( const Cell& cOtherCell) const
               { return ( m_nId == cOtherCell.m_nId) ; }
       void SetBottomLeft( const Point3d& ptBL)
@@ -192,9 +190,11 @@ class Cell
               { return Point3d( m_ptPbl.x, m_ptPtr.y) ; }
       Point3d GetBottomRight( void) const
               { return Point3d( m_ptPtr.x, m_ptPbl.y); }
+      Point3d GetCenter( void) const
+              { return ( m_ptPbl + m_ptPtr) / 2 ; }
       double GetSplitValue( void) const
               { return m_dSplit ; }
-      bool IsSplitVert( void) const                 // se true la cella verrebbe splittata verticalmente, sennò orizzontalmente
+      bool IsSplitVert( void) const                 // se true la cella verrebbe splittata verticalmente, altrimenti orizzontalmente
               { return m_bSplitVert ; }
       bool IsLeaf( void) const                      // flag che indica se la cella ha figli o se è una foglia
               { return ( m_nChild1 == -2  &&  m_nChild2 == -2) ; }
@@ -245,13 +245,12 @@ class Tree
    public :
       ~Tree( void) ;
       Tree( void) ;
-      //Tree ( const SurfBezier* pSrfBz, bool bSplitPatches = true, const Point3d& ptMin = ORIG, const Point3d& ptMax = ORIG) ;
-      Tree( const Point3d ptBl, const Point3d ptTr) ; // creatore da usare solo nel caso in cui si voglia aggiungere tagli ad un'unica cella e del risultato ottenere il contorno 
-      bool SetSurf( const SurfBezier* pSrfBz, bool bSplitPatches = true, const Point3d& ptMin = ORIG, const Point3d& ptMax = ORIG) ;
-      bool GetIndependentTrees( BIPNTVECTOR& vTrees) ;                             // calcolo la suddivisione della superficie solo sulle singole bbox dei loop di trim ( unendo quelli vicini)
-      bool BuildTree( double dLinTol = LIN_TOL_STD, double dSideMin = 1, double dSideMax = INFINITO) ; // dSideMax è il massimo per la dimensione maggiore di un triangolo della trimesh
-                                                                                                       // dSideMin è lunghezza minima del lato di una cella nello spazio reale
-      bool BuildTree_test( double dLinTol = LIN_TOL_STD, double dSideMin = 1, double dSideMax = INFINITO) ;
+      Tree( const Point3d ptBl, const Point3d ptTr) ;    // da usare solo nel caso in cui si voglia aggiungere tagli ad un'unica cella e del risultato ottenere il contorno 
+      bool SetSurf( const SurfBezier* pSrfBz, const Point3d& ptMin = ORIG, const Point3d& ptMax = ORIG) ;
+      bool GetIndependentTrees( BIPNTVECTOR& vTrees) ;       // calcolo la suddivisione della superficie solo sulle singole bbox dei loop di trim ( unendo quelli vicini)
+      bool BuildTree( double dLinTol = LIN_TOL_STD,
+                      double dSideMin = 1,                   // è la minima lunghezza del lato di una cella
+                      double dSideMax = INFINITO) ;          // è la massima dimensione di un triangolo della trimesh
       bool GetPolygons( POLYLINEMATRIX& vvPolygons, POLYLINEMATRIX& vvPolygons3d, std::vector<ICRVCOMPOPOVECTOR>& vCCEdges3D, ICRVCOMPOPOVECTOR& vCCLoops) ;
       bool GetPolygonsBasic( POLYLINEVECTOR& vPolygons, POLYLINEVECTOR& vPolygonsCorrected, POLYLINEVECTOR& vPolygons3d) ; // restituisce il poligono corrispondente ad ogni cella foglia dell'albero
                                                                                                                            // ad ogni poligono sono stati aggiunti tutti i vertici dei vicini posizionati sui suoi lati
@@ -260,26 +259,27 @@ class Tree
       bool GetEdges3D ( std::vector<ICRVCOMPOPOVECTOR>& mCCEdge, POLYLINEVECTOR& vPolygons) ;  // restituisce gli edge 3D come polyline
       bool GetSplitLoops( ICRVCOMPOPOVECTOR& vCCLoopSplit) const // funzione che restituisce i loop splitatti ai confini delle celle
              { for ( int i = 0 ; i < int( m_vCCLoop2D.size()); ++i) vCCLoopSplit.emplace_back( m_vCCLoop2D[i]->Clone()) ; return true ; };
-      void SetTestMode( void) { m_bTestMode = true ;} ;        // attivando la test mode, per la costruzione dell'albero viene usata la funzione BuiltTree_test e viene corretta di conseguenza la FindCell
       // funzioni da usare per ricostruire tagli che vanno aggiunti allo spazio parametrico
       bool AddCutsToRoot( POLYLINEVECTOR& vCuts) ;          // aggiunge i tagli al tree
       bool CreateCellContour( POLYLINEMATRIX& vPolygons) ;  // crea il nuovo contorno esterno, tenendo conto dei tagli
-      bool IsClosedU( void) const { return m_bClosedU ;} ;        // funzione che riferisce se la superficie è chiusa lungo il parametro U
-      bool IsClosedV( void) const { return m_bClosedV ;} ;        // funzione che riferisce se la superficie è chiusa lungo il parametro V
-      std::vector<bool> GetPoles( void) { return m_vbPole ;} ; // funzione che restituisce i flag che indicano se i lati sono collassati in dei poli
+      bool IsClosedU( void) const                                 // restituisce flag di chiusara in U
+             { return m_bClosedU ; }
+      bool IsClosedV( void) const                                 // restituisce flag di chiusara in V
+             { return m_bClosedV ; }
+      BOOLVECTOR GetPoles( void)                                  // restituisce i flag che indicano se i lati sono collassati in dei poli
+             { return m_vbPole ; }
 
    private :
       bool Split( int nId, double dSplitValue) ;                              // funzione di split di una cella al parametro indicato nella direzione data da bVert
       bool Split( int nId) ;                                                  // funzione di split di una cella dell'albero a metà nella direzione data da bVert
-      void Balance( void) ;                                                       // creo rami in modo che tutte tutte le foglie abbiano come adiacenti foglie ad una profondità di +- 1
-      int GetHeightLeaves( int nId, INTVECTOR& vnLeaves, int d = 0) const ;  // altezza del subtree a partire dal nodo nId
-      int GetDepth( int nId, int nRef) const ;                               // livello del nodo nId
+      int  GetHeightLeaves( int nId, INTVECTOR& vnLeaves, int d = 0) const ;  // altezza del subtree a partire dal nodo nId
+      int  GetDepth( int nId, int nRef) const ;                               // livello del nodo nId
       void GetTopNeigh( int nId, INTVECTOR& vTopNeighs) const ;               // restituisce le celle foglie che sono adiacenti al lato top
       void GetBottomNeigh( int nId, INTVECTOR& vBottomNeighs) const ;         // restituisce le celle foglie che sono adiacenti al lato bottom
       void GetLeftNeigh( int nId, INTVECTOR& vLeftNeighs) const ;             // restituisce le celle foglie che sono adiacenti al lato left
       void GetRightNeigh( int nId, INTVECTOR& vRightNeighs) const ;           // restituisce le celle foglie che sono adiacenti al lato right
       void GetRootNeigh( int nEdge, INTVECTOR& vNeigh) ;                      // restituisce le foglie dell'albero che sono adiacenti al lato nEdge, numerato a partire dal top ( 0) in senso antiorario 
-      void ResetTree( void) ;                                                // resetto m_bProcessed a false per tutti i nodi dell'albero
+      void ResetTree( void) ;                                                 // resetto m_bProcessed a false per tutti i nodi dell'albero
       INTVECTOR FindCell( const Point3d& ptToAssign, const CurveLine& cl, bool bRecurs = false) const ;   // dato un punto, trova la cella foglia a cui appartiene
       INTVECTOR FindCell( const Point3d& ptToAssign, const CurveLine& cl, INTVECTOR vCells, bool bRecurs = false) const ;   // dato un punto, trova la cella foglia a cui appartiene
       bool TraceLoopLabelCell( const POLYLINEVECTOR& vplPolygons) ;               // tracing dei loop e labelling delle celle
@@ -307,10 +307,8 @@ class Tree
       bool GetPoint(double dU, double dV, Point3d& pt) const ;
       bool SavePoint( double dU, double dV, Point3d& pt) ;
 
-
    private :
       const SurfBezier*           m_pSrfBz ;           // superficie di bezier
-      DBLVECTOR                   m_vDim ;             // distanze tra i vertici della superficie di bezier in 3d in ordine antiorario a partire da ptP00
       bool                        m_bTrimmed ;         // superficie trimmata
       std::unordered_map<int,int>      m_mChunk ;           // mappa in cui vengono salvati chunk di appartenza per ogni loop di trim
       std::vector<std::tuple<PolyLine,bool>> m_vPlApprox ;       // vettore contenente le approssimazioni dei loop  // il bool indica se la curva è CCW
@@ -319,7 +317,6 @@ class Tree
       bool                        m_bClosedU ;         // superficie chiusa lungo il parametro U
       bool                        m_bClosedV ;         // superficie chiusa lungo il parametro V
       BOOLVECTOR                  m_vbPole ;           // vettore che indica se i vari lati sono collassati in poli ( indici riferiti all'ordine degli edge)
-      bool                        m_bSplitPatches ;    // flag che indica se le patches sono state divise prima della creazione dell'albero
       int                         m_nDegU ;            // grado della superficie nel parametro U
       int                         m_nDegV ;            // grado della superficie nel parametro V
       int                         m_nSpanU ;           // numero di span lungo il parametro U
@@ -328,7 +325,6 @@ class Tree
       mutable std::unordered_map<std::pair<int64_t, int64_t>,Point3d, PairHashInt64>     m_mPt3d ;    // mappa che contiene tutti i punti 3d della superficie calcolati (la chiave sono le coordinate, moltiplicate per 2^24 e trasformate in int)
       INTVECTOR                   m_vnLeaves ;         // vettore delle foglie
       INTVECTOR                   m_vnParents ;        // vettore delle celle ottenute dalla divisione preliminare in singole patch
-      bool                        m_bTestMode ;        // bool che indica se la test mode è attiva
       ICRVCOMPOPOVECTOR           m_vCCLoop2D ;        // vettore che contiene le CurveCompo che rappresentano i loop di trim tenendo conto della divisione in celle
       std::vector<std::pair<BIPNTVECTOR, ChainCurves>>    m_vCEdge2D ;        // vettore che le chain che rappresentano ciò che resta degli edge originali, tenendo conto dei trim.
 } ;

VolZmapVolume.cpp

@@ -1080,12 +1080,12 @@ VolZmap::MillingStep( int nCurrTool,
    Vector3d vtALe = GetToLoc( vtAe, m_MapFrame) ;
    vtALe.Normalize() ;
 
-   //static PerformanceCounter Counter ;
-   //{
-   //string sOut = "Draw=" + ToString( Counter.Stop(), 3) ;
-   //LOG_INFO( GetEGkLogger(), sOut.c_str())
-   //Counter.Start() ;
-   //}
+   static PerformanceCounter Counter ;
+   {
+   string sOut = "Draw=" + ToString( Counter.Stop(), 3) ;
+   LOG_INFO( GetEGkLogger(), sOut.c_str())
+   Counter.Start() ;
+   }
 
   // Se pura traslazione
    bool bOk ;
@@ -1095,11 +1095,20 @@ VolZmap::MillingStep( int nCurrTool,
    else
       bOk = MillingGeneralMotionStep( ptPLs, vtDLs, vtALs, ptPLe, vtDLe, vtALe) ;
 
-   //{
-   //string sOut = "Calc=" + ToString( Counter.Stop(), 3) ;
-   //LOG_INFO( GetEGkLogger(), sOut.c_str())
-   //Counter.Start() ;
-   //}
+   {
+   // debug per 5 assi
+   double dPartCount = Counter.Stop() ;
+   if ( ! AreSameVectorApprox( vtDLs, vtDLe)) {
+      dTotCount += dPartCount ;
+      string sOut = "Calc5Axis=" + ToString( dTotCount, 3) ;
+      LOG_INFO( GetEGkLogger(), sOut.c_str())
+   }
+   //debug per 5 assi
+
+   string sOut = "Calc=" + ToString( Counter.Stop(), 3) ;
+   LOG_INFO( GetEGkLogger(), sOut.c_str())
+   Counter.Start() ;
+   }
 
    return bOk ;
 }
@@ -2107,7 +2116,7 @@ VolZmap::Comp_5AxisMilling( int nGrid, const PNTVECTOR& ptS, const PNTVECTOR& pt
                   vvPtCtrl.emplace_back( PNTVECTOR({ vPntTipStartBack[i+1], vPntTipStartBack[i], vPntTipEndBack[i+1], vPntTipEndBack[i]})) ;
                // superfici superiori
                for ( int i = 0 ; i < nSub ; ++i)
-                  vvPtCtrl.emplace_back( PNTVECTOR({ vPntTopStartFront[i], vPntTopStartFront[i+1], vPntTopEndFront[i], vPntTopEndFront[i+1]})) ;
+                  vvPtCtrl.emplace_back( PNTVECTOR({ vPntTopStartFront[i+1], vPntTopStartFront[i], vPntTopEndFront[i+1], vPntTopEndFront[i]})) ;
             }
          }
          // aggiungo le superfici bilineari anche nell'altro verso sia per il top che per il tip