EgtGeomKernel :

- aggiunta dipendenza header
- aggiornata cartella header.

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

@@ -515,6 +515,10 @@ CurveToBezierCurve( const ICurve* pCrv, int nDeg, bool bMakeRatOrNot)
 ICurveBezier*
 LineToBezierCurve( const ICurveLine* pCrvLine, int nDeg, bool bMakeRatOrNot)
 {
+   // verifico che esista la linea
+   if ( pCrvLine == nullptr)
+      return nullptr ;
+
    PtrOwner<ICurveBezier> pCrvBezier( CreateCurveBezier()) ;
    // rendo tutte le curve di grado 2 e razionali così posso convertire anche archi e avere tutte curve dello stesso grado e razionali
    pCrvBezier->Init( nDeg, true) ;
@@ -526,12 +530,11 @@ LineToBezierCurve( const ICurveLine* pCrvLine, int nDeg, bool bMakeRatOrNot)
 
 //----------------------------------------------------------------------------
 ICurve*
-ArcToBezierCurve( const ICurve* pCrv, int nDeg, bool bMakeRatOrNot)
+ArcToBezierCurve( const ICurveArc* pArc, int nDeg, bool bMakeRatOrNot)
 {
    // una spirale non può essere forzata al grado 2
 
-  // verifico sia un arco
-   const CurveArc* pArc = GetBasicCurveArc( pCrv) ;
+  // verifico che esista l'arco
    if ( pArc == nullptr)
       return nullptr ;
 
@@ -541,8 +544,15 @@ ArcToBezierCurve( const ICurve* pCrv, int nDeg, bool bMakeRatOrNot)
       PtrOwner<ICurveBezier> pCrvBez( CreateBasicCurveBezier()) ;
       if ( IsNull( pCrvBez)  ||  ! pCrvBez->FromArc( *pArc))
          return nullptr ;
+      if ( ! bMakeRatOrNot) {
+          Point3d ptCen = pArc->GetCenter() ;
+          Vector3d vtN = pArc->GetNormVersor() ;
+          pCrvBez.Set( ApproxArcCurveBezierWithSingleCubic( pCrvBez, ptCen, vtN)) ;
+      }
+      if ( IsNull( pCrvBez))
+          return nullptr ;
       // aumento il grado della curva come richiesto
-      while( pCrvBez->GetDegree() < nDeg)
+      while ( pCrvBez->GetDegree() < nDeg)
          pCrvBez.Set( BezierIncreaseDegree( pCrvBez)) ;
      // restituisco la curva
       return Release( pCrvBez) ;
@@ -555,20 +565,18 @@ ArcToBezierCurve( const ICurve* pCrv, int nDeg, bool bMakeRatOrNot)
          return nullptr ;
      // inserisco nella CC le curve di Bezier equivalenti alle parti dell'arco
       int nParts = (int) ceil( abs( pArc->GetAngCenter()) / ( BEZARC_ANG_CEN_MAX + EPS_ANG_SMALL)) ;
-      //if ( ! bDeg3OrDeg2  &&  abs( pArc->GetDeltaN()) > EPS_ZERO)
-      //   nParts *= 2 ;
       nParts = max( nParts, 2) ;
       for ( int i = 0 ; i < nParts ; ++ i) {
         // copio l'arco originale
-         CurveArc cArc = *pArc ;
+         CurveArc cArc = *GetBasicCurveArc( pArc) ;
         // lo limito alla parte di interesse
-         cArc.TrimStartEndAtParam( i / double( nParts), ( i + 1) / double( nParts)) ;
+         cArc.TrimStartEndAtParam( double( i) / nParts, double( i + 1) / nParts) ;
         // creo la curva di Bezier equivalente
          PtrOwner<ICurveBezier> pCrvBez( CreateBasicCurveBezier()) ;
          if ( IsNull( pCrvBez)  ||  ! pCrvBez->FromArc( cArc))
             return nullptr ;
          // aumento il grado della curva come richiesto
-         while( pCrvBez->GetDegree() < nDeg)
+         while ( pCrvBez->GetDegree() < nDeg)
             pCrvBez.Set( BezierIncreaseDegree( pCrvBez)) ;
         // aggiungo la curva di Bezier a quella composita
          if ( ! pCrvCompo->AddCurve( Release( pCrvBez)))
@@ -586,7 +594,11 @@ ArcToBezierCurve( const ICurve* pCrv, int nDeg, bool bMakeRatOrNot)
 ICurve*
 CompositeToBezierCurve( const ICurveComposite* pCC, int nDeg, bool bMakeRatOrNot)
 {
-   // converto tutte le curve in bezier razionali di grado 2
+  // verifico che esista la curva
+   if ( pCC == nullptr)
+      return nullptr ;
+
+  // converto tutte le curve in bezier razionali di grado 2
    PtrOwner<ICurveComposite> pCCBezier( CreateCurveComposite()) ;
    for ( int i = 0 ; i < int( pCC->GetCurveCount()) ; ++i) {
       PtrOwner<ICurve> pCrvNew ;
@@ -595,7 +607,7 @@ CompositeToBezierCurve( const ICurveComposite* pCC, int nDeg, bool bMakeRatOrNot
          ICurve* pCrvBezier = ArcToBezierCurve( crArc) ;
          if ( pCrvBezier == nullptr)
             return nullptr ;
-         // se la curva è di grado superiore al secondo allora devo ricondurla al secondo grado
+        // se la curva è di grado superiore al secondo allora devo ricondurla al secondo grado
          pCrvNew.Set( pCrvBezier) ;
       }
       else if ( pCC->GetCurve(i)->GetType() == CRV_LINE) {
@@ -624,25 +636,28 @@ EditBezierCurve( const ICurveBezier* pCrvBezier, int nDeg, bool bMakeRatOrNot, d
 { 
    // se nDeg == -1 allora viene mantenuto il grado della curva originale
 
-   if( nDeg == 2 ||  nDeg == 1)
+  // verifico sia una bezier
+   if ( pCrvBezier == nullptr)
+      return nullptr ;
+
+   if ( nDeg == 2  ||  nDeg == 1)
       return nullptr ;
 
    PtrOwner<ICurveBezier> pCrvNew( pCrvBezier->Clone()) ;
    int nDegCurr = pCrvNew->GetDegree() ;
    bool bRat = pCrvNew->IsRational() ;
    int nDegWanted = nDeg == -1 ? nDegCurr : nDeg ;
-   if( ! bMakeRatOrNot) {
-      if( ! pCrvNew->MakeNonRational( dTol)) {
+   if ( ! bMakeRatOrNot) {
+      if ( ! pCrvNew->MakeNonRational( dTol)) {
          // se ho fallito la conversione diretta in curva non razionale allora la spezzo in bezier cubiche
          PtrOwner<ICurveComposite> pBezCubics( CreateCurveComposite()) ;
          pBezCubics->AddCurve( ApproxBezierWithCubics(pCrvBezier, dTol)) ;
-         if( IsNull( pBezCubics))
+         if ( IsNull( pBezCubics))
             return nullptr ;
          // adatto ogni sottocurva cubica
          PtrOwner<ICurveComposite> pCCEdited( CreateCurveComposite()) ;
-         
          for ( int i = 0 ; i < pBezCubics->GetCurveCount() ; ++i) {
-            if( ! pCCEdited->AddCurve( EditBezierCurve( GetCurveBezier( pBezCubics->GetCurve( i)), nDegWanted, bMakeRatOrNot, dTol)) )
+            if ( ! pCCEdited->AddCurve( EditBezierCurve( GetCurveBezier( pBezCubics->GetCurve( i)), nDegWanted, bMakeRatOrNot, dTol)) )
                return nullptr ;
          }
          return Release( pCCEdited) ;
@@ -661,7 +676,7 @@ EditBezierCurve( const ICurveBezier* pCrvBezier, int nDeg, bool bMakeRatOrNot, d
    else if ( nDegCurr > nDegWanted) {
       while ( nDegCurr > nDegWanted) {
          ICurveBezier* pCrvDec = BezierDecreaseDegree( pCrvNew, dTol) ;
-         if( pCrvDec == nullptr  ||  ! pCrvDec->IsValid()) {
+         if ( pCrvDec == nullptr  ||  ! pCrvDec->IsValid()) {
             // se ho fallito la riduzione di grado entro la tolleranza richiesta allora la spezzo in bezier cubiche prima di adattare
             PtrOwner<ICurveComposite> pBezCubics( CreateCurveComposite()) ;
             pBezCubics->AddCurve( ApproxBezierWithCubics(pCrvBezier, dTol)) ;
@@ -687,20 +702,22 @@ EditBezierCurve( const ICurveBezier* pCrvBezier, int nDeg, bool bMakeRatOrNot, d
 
 //----------------------------------------------------------------------------
 ICurveBezier*
-BezierIncreaseDegree(const ICurveBezier* pCrvBezier)
+BezierIncreaseDegree( const ICurveBezier* pCrvBezier)
 {
    if ( pCrvBezier == nullptr)
       return nullptr ;
-   // creo la versione con grado aumentato
+  // creo la versione con grado aumentato
    PtrOwner<ICurveBezier> pNewBezier( CreateCurveBezier()) ;
-   int nDeg = pCrvBezier->GetDegree() + 1;
-   bool bRat = pCrvBezier->IsRational();
-   pNewBezier->Init( nDeg , bRat) ;
-   // prev e curr sono riferiti alla curva di partenza
-   // salvo il primo punto
+   if ( IsNull( pNewBezier))
+      return nullptr ;
+   int nDeg = pCrvBezier->GetDegree() + 1 ;
+   bool bRat = pCrvBezier->IsRational() ;
+   pNewBezier->Init( nDeg, bRat) ;
+  // prev e curr sono riferiti alla curva di partenza
+  // salvo il primo punto
    Point3d ptCtrlPrev = pCrvBezier->GetControlPoint( 0) ;
    double dWprev = 1 ;
-   if ( bRat ) {
+   if ( bRat) {
       dWprev = pCrvBezier->GetControlWeight( 0) ;
       pNewBezier->SetControlPoint( 0, ptCtrlPrev, dWprev) ;
    }
@@ -709,26 +726,26 @@ BezierIncreaseDegree(const ICurveBezier* pCrvBezier)
    // ciclo sui punti di controllo intermedi per calcolare quelli nuovi
    Point3d ptCtrlCurr ;
    double dWcurr ;
-   for ( double i = 1 ; i < nDeg ; ++i) {
-      ptCtrlCurr = pCrvBezier->GetControlPoint( int( i)) ;
-      double dAlpha = i / nDeg ;
+   for ( int i = 1 ; i < nDeg ; ++i) {
+      ptCtrlCurr = pCrvBezier->GetControlPoint( i) ;
+      double dAlpha = double( i) / nDeg ;
       if ( bRat) {
-         dWcurr = pCrvBezier->GetControlWeight( int( i)) ;
+         dWcurr = pCrvBezier->GetControlWeight( i) ;
          double dWnew = dAlpha * dWprev + ( 1 - dAlpha) * dWcurr ;
          Point3d ptNew = dAlpha * ptCtrlPrev * dWprev + ( 1 - dAlpha) * ptCtrlCurr * dWcurr;
          ptNew /= dWnew ;
-         pNewBezier->SetControlPoint( int( i), ptNew, dWnew) ;
+         pNewBezier->SetControlPoint( i, ptNew, dWnew) ;
          dWprev = dWcurr ;
       }
       else {
          Point3d ptNew = dAlpha * ptCtrlPrev + ( 1 - dAlpha) * ptCtrlCurr ;
-         pNewBezier->SetControlPoint( int( i), ptNew) ;
+         pNewBezier->SetControlPoint( i, ptNew) ;
       }
       ptCtrlPrev = ptCtrlCurr ;
    }
    // salvo l'ultimo punto
    ptCtrlCurr = pCrvBezier->GetControlPoint( nDeg - 1) ;
-   if ( bRat ) {
+   if ( bRat) {
       dWcurr = pCrvBezier->GetControlWeight( nDeg - 1) ;
       pNewBezier->SetControlPoint( nDeg, ptCtrlCurr, dWcurr) ;
    }
@@ -889,8 +906,13 @@ BezierDecreaseDegree(const ICurveBezier* pCrvBezier, double dTol)
 
 //----------------------------------------------------------------------------
 ICurve*
-ApproxBezierWithCubics(const ICurveBezier* pCrvBezier, double dTol)
+ApproxBezierWithCubics(const ICurve* pCrv, double dTol)
 {
+   // verifico sia una bezier
+   const CurveBezier* pCrvBezier = GetBasicCurveBezier( pCrv) ;
+   if ( pCrvBezier == nullptr)
+      return nullptr ;
+
    // cerco di stimare quanti cambi di concavità ho
    // tiro una linea tra il primo punto di controllo e l'ultimo e poi scorro gli altri punti di controllo
    // controllando quante volte salto da un lato all'altro della linea
@@ -911,10 +933,8 @@ ApproxBezierWithCubics(const ICurveBezier* pCrvBezier, double dTol)
       Vector3d vtStartDir, vtEndDir ;
       pCrvBezier->GetPointD1D2( 0, ICurve::FROM_MINUS, ptStart, &vtStartDir) ;
       pCrvBezier->GetPointD1D2( 1, ICurve::FROM_PLUS, ptEnd, &vtEndDir) ;
-      Plane3d plPLane ;
-      plPLane.Set( ptStart, ptCtrl, ptEnd) ;
-      vtStartDir.Rotate( plPLane.GetVersN(), 90) ;
-      vtEndDir.Rotate( plPLane.GetVersN(), 90) ;
+      vtStartDir.Rotate( vtN, 90) ;
+      vtEndDir.Rotate( vtN, 90) ;
       Point3d ptAux1 = ptStart + vtStartDir ;
       Point3d ptAux2 = ptEnd + vtEndDir ;
       PtrOwner<CurveLine> pCL1( CreateBasicCurveLine()) ;
@@ -929,7 +949,7 @@ ApproxBezierWithCubics(const ICurveBezier* pCrvBezier, double dTol)
          double dDist2 = Dist( ptEnd, ptCen) ;
          if ( abs(dDist1 - dDist2) < EPS_SMALL ) {
             CurveArc cArc ;
-            cArc.SetC2PN( ptCen, ptStart, ptEnd, plPLane.GetVersN()) ;
+            cArc.SetC2PN( ptCen, ptStart, ptEnd, vtN) ;
             // controllo se il raggio della circonferenza risultante coincide con quello calcolato prima
             // impongo inoltre che l'angolo al centro sia minore di 90 gradi
             if ( abs( cArc.GetRadius() - dDist1) < EPS_SMALL  &&  cArc.GetAngCenter() < 90 + EPS_SMALL)
@@ -973,7 +993,7 @@ ApproxBezierWithCubics(const ICurveBezier* pCrvBezier, double dTol)
       if ( ! bIsArc)
          pCrvCubic.Set( ApproxCurveBezierWithSingleCubic( pCrvPart)) ;
       else {
-         pCrvCubic.Set( ApproxArcCurveBezierWithSingleCubic( pCrvPart, ptCen)) ;
+         pCrvCubic.Set( ApproxArcCurveBezierWithSingleCubic( pCrvPart, ptCen, vtN)) ;
          if ( IsNull( pCrvCubic)) {
             // se fallisce allora riprovo usando più di una bezier
             bIsArc = false ;
@@ -1012,8 +1032,13 @@ ApproxBezierWithCubics(const ICurveBezier* pCrvBezier, double dTol)
 
 //----------------------------------------------------------------------------
 ICurveBezier*
-ApproxCurveBezierWithSingleCubic( const ICurveBezier* pCrvBez)
+ApproxCurveBezierWithSingleCubic( const ICurve* pCrv)
 {
+   // verifico sia una bezier
+   const CurveBezier* pCrvBez = GetBasicCurveBezier( pCrv) ;
+   if ( pCrvBez == nullptr)
+      return nullptr ;
+
    Point3d ptStart, ptEnd ;
    pCrvBez->GetStartPoint( ptStart) ;
    pCrvBez->GetEndPoint( ptEnd) ;
@@ -1042,8 +1067,22 @@ ApproxCurveBezierWithSingleCubic( const ICurveBezier* pCrvBez)
 
 //----------------------------------------------------------------------------
 ICurveBezier*
-ApproxArcCurveBezierWithSingleCubic( const ICurveBezier* pCrvBez, const Point3d& ptCen)
+ApproxArcCurveBezierWithSingleCubic( const ICurve* pCrv, const Point3d& ptCen, const Vector3d& vtN)
 {
+    // verifico sia una bezier
+    const CurveBezier* pCrvBez = GetBasicCurveBezier( pCrv) ;
+    if ( pCrvBez == nullptr)
+        return nullptr ;
+
+   // mi metto nel frame della curva ( lavoro nel piano XY)
+    Point3d ptStart = pCrvBez->GetControlPoint( 0) ;
+    Point3d ptEnd = pCrvBez->GetControlPoint( 2) ;
+    Frame3d frCrv; frCrv.Set( ptStart, vtN) ;
+    ptStart = ORIG ;
+    ptEnd.ToLoc( frCrv) ;
+    Point3d ptCenLoc = ptCen ;
+    ptCenLoc.ToLoc( frCrv) ;
+
    // converto una curva di bezier che definisce un arco perfetto di circonferenza
    // dato il centro( della circonferenza), inizio e fine
    // N.B. : per archi con angolo al centro < 90 
@@ -1052,9 +1091,7 @@ ApproxArcCurveBezierWithSingleCubic( const ICurveBezier* pCrvBez, const Point3d&
    bool bRat = pCrvBez->IsRational() ;
    if ( nDeg != 2  ||  ! bRat)
       return nullptr ;
-
-   Point3d ptStart = pCrvBez->GetControlPoint( 0) ;
-   Point3d ptEnd = pCrvBez->GetControlPoint( 2) ;
+   
    if ( AreSamePointEpsilon( ptStart, ptEnd, EPS_SMALL * 50))
       return nullptr ;
    nDeg = 3 ;
@@ -1064,13 +1101,13 @@ ApproxArcCurveBezierWithSingleCubic( const ICurveBezier* pCrvBez, const Point3d&
    pCrvCubic->SetControlPoint( 3, ptEnd) ;
    // from the article of Aleksas Riškus "Approximation of a cubic bezier curve by circular arcs and vice versa"
    // with corrections from "Hans Muller's Flex Blog", page "More About Approximating Circular Arcs With a Cubic Bezier Path"
-   double ax = ptStart.x - ptCen.x ;
-   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 ax = ptStart.x - ptCenLoc.x ;
+   double ay = ptStart.y - ptCenLoc.y ;
+   double bx = ptEnd.x - ptCenLoc.x ;
+   double by = ptEnd.y - ptCenLoc.y ;
+   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 ;
@@ -1079,6 +1116,9 @@ ApproxArcCurveBezierWithSingleCubic( const ICurveBezier* pCrvBez, const Point3d&
    pCrvCubic->SetControlPoint( 1, ptCtrl1) ;
    pCrvCubic->SetControlPoint( 2, ptCtrl2) ;
 
+   // ritorno nel frame originale
+   pCrvCubic->ToGlob( frCrv) ;
+
    return Release( pCrvCubic) ;
 }
 
@@ -1086,7 +1126,9 @@ ApproxArcCurveBezierWithSingleCubic( const ICurveBezier* pCrvBez, const Point3d&
 ICurve*
 ApproxCurveWithBezier( const ICurve*, double dTol)
 {
-   // interpolazione di punti con una curva bezier
+   // campiono punti lungo la curva e poi li interpolo
+
+   // oppure faccio la fat curve e poi calcolo una bezier che stia all'interno di quella regione
 
    PtrOwner<ICurveComposite> pCC( CreateBasicCurveComposite()) ;
    return Release( pCC) ;
@@ -1262,7 +1304,7 @@ CurveToNoArcsCurve( const ICurve* pCrv)
       return nullptr ;
   // se arco, devo trasformarlo in curva di Bezier (semplice o composta)
    if ( pCrv->GetType() == CRV_ARC) {
-      return ArcToBezierCurve( pCrv) ;
+      return ArcToBezierCurve( GetCurveArc( pCrv)) ;
    }
   // se curva composita, devo trasformarla in composita senza archi
    else if ( pCrv->GetType() == CRV_COMPO) {

CurveAux.h

@@ -33,5 +33,5 @@ bool CurveGetArea( const ICurve& crvC, Plane3d& plPlane, double& dArea) ;
 bool CurveDump( const ICurve& crvC, std::string& sOut, bool bMM, const char* szNewLine) ;
 bool CopyExtrusion( const ICurve* pSouCrv, ICurve* pDestCrv) ;
 bool CopyThickness( const ICurve* pSouCrv, ICurve* pDestCrv) ;
-ICurveBezier* ApproxCurveBezierWithSingleCubic( const ICurveBezier* pCrvBez) ;
+ICurveBezier* ApproxCurveBezierWithSingleCubic( const ICurve* pCrv) ;
 Voronoi* GetCurveVoronoi( const ICurve& crvC) ;

CurveBezier.cpp

@@ -369,6 +369,7 @@ CurveBezier::CopyFrom( const CurveBezier& cbSrc)
       return true ;
    if ( ! Init( cbSrc.m_nDeg, cbSrc.m_bRat))
       return false ;
+   m_dParSing = cbSrc.m_dParSing ;
    m_vPtCtrl = cbSrc.m_vPtCtrl ;
    if ( cbSrc.m_bRat)
       m_vWeCtrl = cbSrc.m_vWeCtrl ;
@@ -2274,17 +2275,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 +2294,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 +2324,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 +2335,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 +2354,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 ;
    }
 
@@ -2383,3 +2384,11 @@ CurveBezier::IsALine( void) const
    }
    return true ;
 }
+
+//----------------------------------------------------------------------------
+PNTVECTOR
+CurveBezier::GetAllControlPoints( void) const
+{
+   PNTVECTOR vPntCtrl = m_vPtCtrl ;
+   return vPntCtrl ;
+}

CurveBezier.h

@@ -153,6 +153,7 @@ class CurveBezier : public ICurveBezier, public IGeoObjRW
       bool MakeRationalStandardForm( void) override ;
       bool MakeNonRational( double dTol) override ;
       bool IsALine( void) const override ;
+      PNTVECTOR GetAllControlPoints( void) const ;  // non aggiunta in interfaccia
 
    public :                      // IGeoObjRW
       int GetNgeId( void) const override ;

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 ;
    }
 
@@ -3041,7 +3046,7 @@ CurveComposite::ArcsToBezierCurves( void)
      // se arco, devo trasformare in una o più curve di Bezier
       if ( (*Iter)->GetType() == CRV_ARC) {
         // eseguo trasformazione
-         PtrOwner<ICurve> pNewCrv( ArcToBezierCurve( (*Iter))) ;
+         PtrOwner<ICurve> pNewCrv( ArcToBezierCurve( GetCurveArc( *Iter))) ;
          if ( IsNull( pNewCrv))
             return false ;
         // se risultato è singola curva

DistPointSurfTm.cpp

@@ -94,6 +94,8 @@ DistPointSurfTm::Calculate( const Point3d& ptP, const ISurfTriMesh& tmSurf)
 {
   // Inizializzo distanza non calcolata
    m_dDist = - 1. ;
+  // Vettore di indici dei triangoli più vicini inizialmente vuoto
+   m_vnMinDistTriaIndex.clear() ;
   // Controllo se la superficie è chiusa
    m_bIsSurfClosed = tmSurf.IsClosed() ;
 
@@ -184,6 +186,10 @@ DistPointSurfTm::Calculate( const Point3d& ptP, const ISurfTriMesh& tmSurf)
    if ( nMinDistTriaIndex == SVT_NULL)
       return ;
 
+  // Inizializzo il vettore dei triangoli a minima distanza
+   for ( auto& Tria : vTria)
+      m_vnMinDistTriaIndex.emplace_back( Tria.first) ;
+
   // salvo la distanza minima
    m_dDist = dMinDist ;
   // salvo il punto a distanza minima
@@ -279,6 +285,18 @@ DistPointSurfTm::GetMinDistTriaIndex( int& nMinDistIndex) const
    return true ;
 }
 
+//----------------------------------------------------------------------------
+bool
+DistPointSurfTm::GetMinDistTriaIndices( INTVECTOR& vMinDistTriaIndex) const
+{
+  // Distanza non valida
+   if ( m_dDist < - EPS_ZERO)
+      return false ;
+  // Distanza valida
+   vMinDistTriaIndex = m_vnMinDistTriaIndex ;
+   return true ;
+}
+
 //----------------------------------------------------------------------------
 int
 GetSurfTmNearestVertex( const Point3d& ptP, const ISurfTriMesh& tmSurf)

EGkDllMain.cpp

@@ -159,11 +159,20 @@ InitFontManager( const string& sNfeFontDir, const string& sDefaultFont)
 {
   // recupero il font manager
    FontManager& fntMgr = FontManager::GetFontManager() ;
-
   // lo inizializzo
    fntMgr.Init( sNfeFontDir, sDefaultFont) ;
 }
 
+//-----------------------------------------------------------------------------
+void
+SetDefaultFont( const string& sDefaultFont)
+{
+  // recupero il font manager
+   FontManager& fntMgr = FontManager::GetFontManager() ;
+  // imposto il dato
+   fntMgr.SetDefaultFont( sDefaultFont) ;
+}
+
 //-----------------------------------------------------------------------------
 const string&
 GetNfeFontDir( void)

EgtGeomKernel.vcxproj

@@ -117,6 +117,7 @@
       <DebugInformationFormat>ProgramDatabase</DebugInformationFormat>
       <EnablePREfast>false</EnablePREfast>
       <LanguageStandard>stdcpp17</LanguageStandard>
+      <AdditionalIncludeDirectories>C:\EgtDev\Extern\abseil\Include</AdditionalIncludeDirectories>
     </ClCompile>
     <Link>
       <SubSystem>Windows</SubSystem>

FontManager.h

@@ -29,6 +29,8 @@ class FontManager
 
    public :
       bool Init( const std::string& sNfeFontDir, const std::string& sDefaultFont) ;
+      bool SetDefaultFont( const std::string& sDefaultFont)
+               { m_sDefaultFont = sDefaultFont ; return true ; }
       bool SetCurrFont( const std::string& sFont, int nWeight, bool bItalic,
                         double dHeight, double dRatio, double dAddAdvance) ;
       const std::string& GetNfeFontDir( void) const
@@ -54,7 +56,7 @@ class FontManager
       OsFont      m_OsFont ;
 
    private :
-      FontManager( void) {}
+      FontManager( void) : m_bCurrNfeFont( false)  {}
       FontManager( FontManager const& copy) = delete ;
       FontManager& operator=( FontManager const& copy) = delete ;
 } ;

GdbGeo.cpp

@@ -336,7 +336,7 @@ GdbGeo::Scale( const Frame3d& frRef, double dCoeffX, double dCoeffY, double dCoe
         // curva originale
          ICurve* pCrv = GetCurve( m_pGeoObj) ;
         // trasformo in curva di Bezier (semplice o composta)
-         ICurve* pCrvNew = ArcToBezierCurve( pCrv) ;
+         ICurve* pCrvNew = ArcToBezierCurve( GetCurveArc( pCrv)) ;
          if ( pCrvNew == nullptr)
             return false ;
         // assegno alla nuova curva estrusione e spessore di quella originale
@@ -389,7 +389,7 @@ GdbGeo::Shear( const Point3d& ptOn, const Vector3d& vtNorm, const Vector3d& vtDi
       if ( ! pArc->IsPlane()  ||
            ! AreSameOrOppositeVectorExact( pArc->GetNormVersor(), vtNorm)) {
         // trasformo in curva di Bezier (semplice o composta)
-         ICurve* pCrvNew = ArcToBezierCurve( GetCurve( m_pGeoObj)) ;
+         ICurve* pCrvNew = ArcToBezierCurve( GetCurveArc( m_pGeoObj)) ;
          if ( pCrvNew == nullptr)
             return false ;
         // assegno alla nuova curva estrusione e spessore di quella originale

GeomDB.cpp

@@ -877,7 +877,7 @@ GeomDB::GetFirstNameInGroup( int nGroupId, const string& sName) const
      // se ha il nome o la parte iniziale di nome cercato
       string sObjName ;
       if ( pGdbO->GetName( sObjName)  &&  
-           (( ! bWild  &&  sObjName == sToFind)  ||  ( bWild  &&  sObjName.find( sToFind) == 0)))
+           (( ! bWild  &&  sObjName == sToFind)  ||  ( bWild  &&  sObjName.rfind( sToFind, 0) == 0)))
          return ( pGdbO->m_nId) ;
      // passo al successivo
       pGdbO = pGdbO->GetNext() ;
@@ -905,7 +905,7 @@ GeomDB::GetNextName( int nId, const string& sName) const
      // se ha il nome o la parte iniziale di nome cercato
       string sObjName ;
       if ( pGdbNext->GetName( sObjName)  &&  
-           (( ! bWild  &&  sObjName == sToFind)  ||  ( bWild  &&  sObjName.find( sToFind) == 0)))
+           (( ! bWild  &&  sObjName == sToFind)  ||  ( bWild  &&  sObjName.rfind( sToFind, 0) == 0)))
          return ( pGdbNext->m_nId) ;
      // passo al successivo
       pGdbNext = pGdbNext->GetNext() ;
@@ -933,7 +933,7 @@ GeomDB::GetLastNameInGroup( int nGroupId, const string& sName) const
      // se ha il nome o la parte iniziale di nome cercato
       string sObjName ;
       if ( pGdbO->GetName( sObjName)  &&  
-           (( ! bWild  &&  sObjName == sToFind)  ||  ( bWild  &&  sObjName.find( sToFind) == 0)))
+           (( ! bWild  &&  sObjName == sToFind)  ||  ( bWild  &&  sObjName.rfind( sToFind, 0) == 0)))
          return ( pGdbO->m_nId) ;
      // passo al precedente
       pGdbO = pGdbO->GetPrev() ;
@@ -961,7 +961,7 @@ GeomDB::GetPrevName( int nId, const string& sName) const
      // se ha il nome o la parte iniziale di nome cercato
       string sObjName ;
       if ( pGdbPrev->GetName( sObjName)  &&  
-           (( ! bWild  &&  sObjName == sToFind)  ||  ( bWild  &&  sObjName.find( sToFind) == 0)))
+           (( ! bWild  &&  sObjName == sToFind)  ||  ( bWild  &&  sObjName.rfind( sToFind, 0) == 0)))
          return ( pGdbPrev->m_nId) ;
      // passo al precedente
       pGdbPrev = pGdbPrev->GetPrev() ;

IntersLineSurfBez.cpp

@@ -14,11 +14,17 @@
 //--------------------------- Include ----------------------------------------
 #include "stdafx.h"
 #include "CurveLine.h"
+#include "CurveBezier.h"
+#include "CurveComposite.h"
+#include "SurfFlatRegion.h"
 #include "/EgtDev/Include/EGkDistPointLine.h"
+#include "/EgtDev/Include/EGkDistLineLine.h"
+#include "/EgtDev/Include/EGkDistPointSurfFr.h"
 #include "/EgtDev/Include/EGkIntersLineTria.h"
 #include "/EgtDev/Include/EGkIntersLineSurfTm.h"
 #include "/EgtDev/Include/EGkIntersLineSurfBez.h"
 #include "/EgtDev/Include/EGkSurfBezier.h"
+#include "/EgtDev/Include/ENkPolynomialRoots.h"
 
 using namespace std ;
 
@@ -67,11 +73,11 @@ static void
 UpdateInfoIntersLineSurfBz( const Point3d& ptL, const Vector3d& vtDir, int nILT, int nT, const Point3d& ptSP, const Point3d& ptIBz, double dCos,
                             const Point3d& ptSP2, const Point3d& ptIBz2, double dCos2, ILSBIVECTOR& vInfo)
 {
-   if ( nILT == ILTT_IN  ||  nILT == ILTT_EDGE  ||  nILT == ILTT_VERT) {
+   if ( nILT == ILTA_IN  ||  nILT == ILTA_EDGE  ||  nILT == ILTA_VERT  ||  nILT == ILTA_NO_TRIA) {
       double dU = ( ptIBz - ptL) * vtDir ;
       vInfo.emplace_back( nILT, dU, nT, dCos, ptIBz, ptSP) ;
    }
-   else if ( nILT == ILTT_SEGM  ||  nILT == ILTT_SEGM_ON_EDGE) {
+   else if ( nILT == ILTA_SEGM  ||  nILT == ILTA_SEGM_ON_EDGE) {
       double dU = ( ptIBz - ptL) * vtDir ;
       double dU2 = ( ptIBz2 - ptL) * vtDir ;
       vInfo.emplace_back( nILT, dU, dU2, nT, dCos2, ptIBz, ptIBz2, ptSP, ptSP2) ;
@@ -88,8 +94,8 @@ OrderInfoIntersLineSurfBz( ILSBIVECTOR& vInfo)
    // ordino il vettore delle intersezioni secondo il senso crescente del parametro di linea
    sort( vInfo.begin(), vInfo.end(), 
       []( const IntLinSbzInfo& a, const IntLinSbzInfo& 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) ;
+      {  double dUa = ( ( a.nILTA == ILTA_SEGM  ||  a.nILTA == ILTA_SEGM_ON_EDGE) ? ( a.dU + a.dU2) / 2 : a.dU) ;
+         double dUb = ( ( b.nILTA == ILTA_SEGM  ||  b.nILTA == ILTA_SEGM_ON_EDGE) ? ( b.dU + b.dU2) / 2 : b.dU) ;
    return ( dUa < dUb) ; }) ;
 }
 
@@ -146,7 +152,7 @@ IntersLineSurfBz( const Point3d& ptL, const Vector3d& vtL, double dLen, const IS
       double dCos = vtN * vtL ;
       double dCos2 = 0 ;
       // eventualmente ripeto tutto per ptI2 ( se ho un'intersezione con sovrapposizione)
-      if ( InfoTm.nILTT == ILTT_SEGM || InfoTm.nILTT == ILTT_SEGM_ON_EDGE ) {
+      if ( InfoTm.nILTT == ILTA_SEGM || InfoTm.nILTT == ILTA_SEGM_ON_EDGE ) {
          pSurfBz->UnprojectPointFromStm( InfoTm.nT, InfoTm.ptI2, ptSP2, InfoTm.nILTT) ;
          if ( ! RefineIntersNewton(ptL, vtL, dLen, bFinite, pSurfBz, ptSP2, ptIBz2) ) {
             int nVert[3] ;
@@ -176,7 +182,7 @@ FilterLineSurfBzInters( const ILSBIVECTOR& vInfo, INTDBLVECTOR& vInters)
    // ciclo sulle intersezioni
    for ( const auto& Info : vInfo) {
       // se intersezione puntuale
-      if ( Info.nILTT == ILTT_VERT  ||  Info.nILTT == ILTT_EDGE  ||  Info.nILTT == ILTT_IN) {
+      if ( Info.nILTA == ILTA_VERT  ||  Info.nILTA == ILTA_EDGE  ||  Info.nILTA == ILTA_IN) {
          int nFlag = LSBT_TOUCH ;
          if ( Info.dCosDN > EPS_ZERO)
             nFlag = LSBT_OUT ;
@@ -185,7 +191,7 @@ FilterLineSurfBzInters( const ILSBIVECTOR& vInfo, INTDBLVECTOR& vInters)
          vInters.emplace_back( nFlag, Info.dU) ;
       }
       // se altrimenti intersezione con coincidenza
-      else if ( Info.nILTT == ILTT_SEGM  ||  Info.nILTT == ILTT_SEGM_ON_EDGE) {
+      else if ( Info.nILTA == ILTA_SEGM  ||  Info.nILTA == ILTA_SEGM_ON_EDGE) {
          vInters.emplace_back( LSBT_TG_INI, Info.dU) ;
          vInters.emplace_back( LSBT_TG_FIN, Info.dU2) ;
       }
@@ -232,3 +238,255 @@ FilterLineSurfBzInters( const ILSBIVECTOR& vInfo, INTDBLVECTOR& vInters)
    }
    return true ;
 }
+
+//----------------------------------------------------------------------------
+//  Intersezione di una linea con una superficie di Bezier
+//----------------------------------------------------------------------------
+bool
+IntersLineSurfBzBilinear( const Point3d& ptL, const Vector3d& vtL, double dLen, const ISurfBezier* pSurfBz,
+                          ILSBIVECTOR& vInfo, bool bFinite)
+{
+   int nDegU, nDegV, nSpanU, nSpanV ;
+   bool bRat, bTrimmed ;
+   pSurfBz->GetInfo( nDegU, nDegV, nSpanU, nSpanV, bRat, bTrimmed) ;
+
+   // funzione pensata per funzionare solo con una monopatch bilineare
+   if( nDegU > 1 || nDegV > 1 ||  nSpanU > 1  || nSpanV > 1 || bRat)
+      return false ;
+
+   int nInters = int( vInfo.size()) ;
+
+   PNTVECTOR vPntCtrl ;
+   for( int p = 0 ; p < 4 ; ++p) {
+      bool bOk = false ;
+      vPntCtrl.push_back( pSurfBz->GetControlPoint( p, &bOk)) ;
+   }
+
+   Vector3d a = vPntCtrl[3] - vPntCtrl[1] + ( vPntCtrl[0] - vPntCtrl[2]) ;
+   Vector3d b = vPntCtrl[1] - vPntCtrl[0] ;
+   Vector3d c = vPntCtrl[2] - vPntCtrl[0] ;
+   Vector3d d =  vPntCtrl[0] - ORIG ;
+
+   double A1 =  a.x * vtL.z - a.z * vtL.x ;
+   double B1 =  b.x * vtL.z - b.z * vtL.x ;
+   double C1 =  c.x * vtL.z - c.z * vtL.x ;
+   double A2 =  a.y * vtL.z - a.z * vtL.y ;
+   double B2 =  b.y * vtL.z - b.z * vtL.y ;
+   double C2 =  c.y * vtL.z - c.z * vtL.y ;
+
+   double D1 = ( d.x - ptL.x) * vtL.z - ( d.z - ptL.z) * vtL.x ;
+   double D2 = ( d.y - ptL.y) * vtL.z - ( d.z - ptL.z) * vtL.y ;
+
+   DBLVECTOR vdCoeff, vdRoots ;
+   vdCoeff = { (B2 * D1 - B1 * D2), ( A2 * D1 - A1 * D2 + B2 * C1 - B1 * C2), ( A2 * C1 - A1 * C2)} ;
+   int nRoots = PolynomialRoots( 2, vdCoeff, vdRoots) ;
+   bool bFound = false ;
+   for( int w = 0 ; w < nRoots ; ++w) {
+      if ( vdRoots[w] > 0 - EPS_ZERO && vdRoots[w] < 1 + EPS_ZERO ) {
+         double dU = 0, dV = vdRoots[w] ;
+         // verifico che non sia una soluzione con molteplicità > 1
+         bool bAlreadyFound = false ;
+         for ( int k = w - 1 ; k >= 0 &&  ! bAlreadyFound ; --k)
+            bAlreadyFound = abs( dV - vdRoots[k]) < EPS_PARAM ;
+         if( ! bAlreadyFound) {
+            dU = (dV * (C1 - C2) + ( D1 - D2)) / ( dV * ( A2 - A1) + ( B2 - B1)) ;
+            if ( dU > - EPS_ZERO  &&  dU < 1 + EPS_ZERO) {
+               Point3d ptIBez, ptIBez2 ;
+               Vector3d vtN ;
+               pSurfBz->GetPointNrmD1D2(dU, dV, ISurfBezier::Side::FROM_MINUS, ISurfBezier::Side::FROM_MINUS, ptIBez, vtN) ;
+               Point3d ptSP( dU, dV, 0), ptSP2 ;
+               double dCos = vtN * vtL, dCos2 = 0 ;
+               UpdateInfoIntersLineSurfBz( ptL, vtL, ILTA_NO_TRIA, -1, ptSP, ptIBez, dCos, ptSP2, ptIBez2, dCos2, vInfo) ;
+               bFound = true ;
+            }
+         }
+      }
+   }
+
+   // se tutti i coefficienti sono zero allora potrei avere una linea che giace sulla superficie
+   // per trovare i punti di inizio e fine sovrapposizione trovo i punti a minima distanza tra la linea e gli edge della superficie
+   if( ! bFound  &&  abs( vdCoeff[0]) < EPS_ZERO  &&  abs( vdCoeff[1]) < EPS_ZERO &&  abs( vdCoeff[2]) < EPS_ZERO) {
+      ICRVCOMPOPOVECTOR vCrvEdge( 4) ;
+      vCrvEdge[0].Set(pSurfBz->GetCurveOnU( 0)) ;
+      vCrvEdge[1].Set(pSurfBz->GetCurveOnV( 1)) ;
+      vCrvEdge[2].Set(pSurfBz->GetCurveOnU( 1)) ;
+      vCrvEdge[3].Set(pSurfBz->GetCurveOnV( 0)) ;
+      double dAngTolDeg = 5 ;
+      for( int i = 0 ; i < 4 ; ++i) {
+         PolyLine plApprox ; vCrvEdge[0]->ApproxWithLines( EPS_SMALL, dAngTolDeg, ICurve::ApprLineType::APL_STD, plApprox) ;
+         //CurveComposite cCC ;
+         //cCC.FromPolyLine( plApprox) ;
+         int nClosestLine = -1 ;
+         double dMinDist = INFINITO ;
+         Point3d pt ;  plApprox.GetFirstPoint( pt) ;
+         Point3d ptClosest ;
+         int c = 0 ;
+         int nTot = plApprox.GetPointNbr() ;
+         for( int j = 0 ; j < nTot ; ++j) {
+            DistPointLine dpl( pt, ptL, vtL, dLen, bFinite) ;
+            double dDist = INFINITO ;
+            dpl.GetDist( dDist) ;
+            if ( dDist < dMinDist) {
+               nClosestLine = c ;
+               dMinDist = dDist ;
+            }
+            plApprox.GetNextPoint( pt) ;
+            ++ c ;
+         }
+
+         Point3d ptInt1, ptInt2 ;
+         if ( nClosestLine < nTot - 1  && nClosestLine > 0) {
+            // tra i due tratti dell'approssimazione che arrivano al punto selezionato come più vicino, devo trovare quale si avvicina di più
+            Point3d ptStart ; plApprox.GetFirstPoint( ptStart) ;
+            Point3d ptEnd ;
+            for( int z = 1 ; z < nClosestLine - 1 ; ++z)
+               plApprox.GetNextPoint( ptStart) ;
+            plApprox.GetNextPoint( ptEnd) ;
+            // linea precedente al punto
+            Vector3d vtLinePre = ptEnd - ptStart ;
+            double dLenPre = vtLinePre.Len() ;
+            DistLineLine dllPre( ptStart, vtLinePre, dLenPre, ptL, vtL,dLen) ;
+            double dDistPre = INFINITO ;
+            dllPre.GetDist( dDistPre) ;
+            // linea che inzia con quel punto
+            ptStart = ptEnd ;
+            plApprox.GetNextPoint( ptEnd) ;
+            Vector3d vtLineCurr = ptEnd - ptStart ;
+            double dLenCurr = vtLineCurr.Len() ;
+            DistLineLine dllCurr( ptStart, vtLineCurr, dLenCurr, ptL, vtL,dLen) ;
+            double dDistCurr = INFINITO ;
+            dllCurr.GetDist( dDistCurr) ;
+
+            if( dDistPre < dDistCurr)
+               dllPre.GetMinDistPoints( ptInt1, ptInt2) ;
+            else
+               dllCurr.GetMinDistPoints( ptInt1, ptInt2) ;
+         }
+         else if ( nClosestLine == 0){
+            // il punto più vicino è sulla prima linea
+            Point3d ptStart ; plApprox.GetFirstPoint( ptStart) ;
+            Point3d ptEnd ; plApprox.GetNextPoint( ptEnd) ;
+            Vector3d vtLineCurr = ptEnd - ptStart ;
+            double dLenCurr = vtLineCurr.Len() ;
+            DistLineLine dllCurr( ptStart, vtLineCurr, dLenCurr, ptL, vtL,dLen) ;
+            dllCurr.GetMinDistPoints( ptInt1, ptInt2) ;
+         }
+         else if ( nClosestLine == nTot- 1) {
+            // il punto più vicino è sull'ultima linea
+            Point3d ptStart ; plApprox.GetFirstPoint( ptStart) ;
+            Point3d ptEnd ;
+            for( int z = 1 ; z < nClosestLine - 1 ; ++z)
+               plApprox.GetNextPoint( ptStart) ;
+            plApprox.GetNextPoint( ptEnd) ;
+            Vector3d vtLinePre = ptEnd - ptStart ;
+            double dLenPre = vtLinePre.Len() ;
+            DistLineLine dllCurr( ptStart, vtLinePre, dLenPre, ptL, vtL,dLen) ;
+            dllCurr.GetMinDistPoints( ptInt1, ptInt2) ;
+         }
+         
+         double dU1 = 0, dV1 = 0, dU2 = 0, dV2 = 0 ;
+         // se ho trovato due punti vuol dire che la linea coincide con un edge e ho trovato tutto quello che serve
+         if( ! AreSamePointExact( ptInt2, ORIG)) {
+            if( i == 0) {
+               //dV1 = 0 ; dV2 = 0 ;
+               vCrvEdge[0]->GetParamAtPoint( ptInt1, dU1) ;
+               vCrvEdge[0]->GetParamAtPoint( ptInt2, dU2) ;
+            }
+            else if( i == 1) {
+               //dU1 = 1 ; dU2 = 1 ;
+               vCrvEdge[1]->GetParamAtPoint( ptInt1, dV1) ;
+               vCrvEdge[1]->GetParamAtPoint( ptInt2, dV2) ;
+            }
+            else if( i == 2){
+               //dV1 = 1 ; dV2 = 1 ;
+               vCrvEdge[2]->GetParamAtPoint( ptInt1, dU1) ;
+               vCrvEdge[2]->GetParamAtPoint( ptInt2, dU2) ;
+            }
+            else if( i == 3){
+               //dU1 = 0 ; dU2 = 0 ;
+               vCrvEdge[3]->GetParamAtPoint( ptInt1, dV1) ;
+               vCrvEdge[3]->GetParamAtPoint( ptInt2, dV2) ;
+            }
+            Point3d ptIBez1, ptIBez2 ;
+            Vector3d vtN1, vtN2 ;
+            pSurfBz->GetPointNrmD1D2(dU1, dV1, ISurfBezier::Side::FROM_MINUS, ISurfBezier::Side::FROM_MINUS, ptIBez1, vtN1) ;
+            pSurfBz->GetPointNrmD1D2(dU2, dV2, ISurfBezier::Side::FROM_MINUS, ISurfBezier::Side::FROM_MINUS, ptIBez2, vtN2) ;
+            Point3d ptSP1( dU1, dV1, 0) ;
+            double dCos1 = vtN1 * vtL ;
+            Point3d ptSP2( dU2, dV2, 0) ;
+            double dCos2 = vtN2 * vtL ;
+            // se avevo già trovato un punto singolo che coincide col primo punto di questa intersezione sovrapposta, allora cancello l'intersezione singola che
+            // avevo salvato e aggiungo quella sovrapposto che ho trovato ora
+            if( bFound) {
+               int nNewTot = int(vInfo.size()) ;
+               int nNewInters = nNewTot - nInters ;
+               bool bAlreadyFound = false ;
+               for( int i = 0 ; i < nNewInters ; ++i) {
+                  bAlreadyFound = AreSamePointApprox(vInfo[nNewTot - i].ptUV, ptSP1) ||  AreSamePointApprox(vInfo[nNewTot - i].ptUV, ptSP2) ;
+                  if ( bAlreadyFound) {
+                     vInfo.erase( vInfo.begin() + nNewTot - i) ;
+                     break ;
+                  }
+               }
+            }
+            UpdateInfoIntersLineSurfBz( ptL, vtL, ILTA_NO_TRIA, -1, ptSP1, ptIBez1, dCos1, ptSP2, ptIBez2, dCos2, vInfo) ;
+            bFound = true ;
+            break ;
+         }
+         // se ho trovato un punto a distanza zero dalla linea allora ho trovato l'intersezione
+         else if( dMinDist < EPS_SMALL) {
+            if( i == 0) {
+               //dV1 = 0 ;
+               vCrvEdge[0]->GetParamAtPoint( ptInt1, dU1) ;
+            }
+            else if( i == 1) {
+               //dU1 = 1 ;
+               vCrvEdge[1]->GetParamAtPoint( ptInt1, dV1) ;
+            }
+            else if( i == 2){
+               //dV1 = 1 ;
+               vCrvEdge[2]->GetParamAtPoint( ptInt1, dU1) ;
+            }
+            else if( i == 3){
+               //dU1 = 0 ;
+               vCrvEdge[3]->GetParamAtPoint( ptInt1, dV1) ;
+            }
+            Point3d ptSP1( dU1, dV1, 0), ptSP2 ;
+            // se avevo trovato già altri punti controllo di non essere esattamente su una diagonale ( e quindi avere un'intersezione con ogni edge, ma due sono doppie)
+            if( bFound) {
+               int nNewTot = int(vInfo.size()) ;
+               int nNewInters = nNewTot - nInters ;
+               bool bAlreadyFound = false ;
+               for( int i = 0 ; i < nNewInters ; ++i)
+                  bAlreadyFound = AreSamePointApprox(vInfo[nNewTot - i].ptUV, ptSP1) ;
+               if( bAlreadyFound)
+                  continue ;
+            }
+
+            Point3d ptIBez1, ptIBez2 ;
+            Vector3d vtN1, vtN2 ;
+            pSurfBz->GetPointNrmD1D2(dU1, dV1, ISurfBezier::Side::FROM_MINUS, ISurfBezier::Side::FROM_MINUS, ptIBez1, vtN1) ;
+            double dCos1 = vtN1 * vtL, dCos2 = 0 ;
+            UpdateInfoIntersLineSurfBz( ptL, vtL, ILTA_NO_TRIA, -1, ptSP1, ptIBez1, dCos1, ptSP2, ptIBez2, dCos2, vInfo) ;
+            bFound = true ;
+         }
+      }
+   }
+
+   // se la superficie è trimmed verifico che i punti trovati siano all'interno del parametrico trimmato
+   if( bTrimmed && bFound) {
+      int nNewTot = int(vInfo.size()) ;
+      int nNewInters = nNewTot - nInters ;
+      const ISurfFlatRegion* pFRTrim = pSurfBz->GetTrimRegion() ;
+      for( int i = 0 ; i < nNewInters ; ++i) {
+         Point3d ptTest = vInfo[nNewTot - i].ptUV * SBZ_TREG_COEFF ;
+         bool bInside = false ;
+         double dDist = INFINITO ;
+         IsPointInsideSurfFr( ptTest, pFRTrim, dDist, bInside) ;
+         if( ! bInside)
+            vInfo.erase( vInfo.begin() + nNewTot - i) ;
+      }
+   }
+
+   return true ;
+}

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)  ||

PolyLine.cpp

@@ -18,6 +18,8 @@
 #include "PolygonPlane.h"
 #include "PointsPCA.h"
 #include "GeoConst.h"
+#include "CurveComposite.h"
+#include "/EgtDev/Include/EGkDistPointCurve.h"
 #include "/EgtDev/Include/EGkPolyLine.h"
 #include "/EgtDev/Include/EGkPlane3d.h"
 #include "/EgtDev/Include/EGkDistPointLine.h"
@@ -796,7 +798,7 @@ PolyLine::RemoveAlignedPoints( double dToler)
       vInd.push_back( 0) ;
       if ( ! DouglasPeuckerSimplification( vPtU, dSqTol, 0, int( vPtU.size()) - 1, vInd))
          return false ;
-      vInd.push_back( vPtU.size() - 1) ;
+      vInd.push_back( int( vPtU.size()) - 1) ;
    }
   // altrimenti chiusa
    else {
@@ -817,7 +819,7 @@ PolyLine::RemoveAlignedPoints( double dToler)
       vInd.push_back( nMaxInd) ;
       if ( ! DouglasPeuckerSimplification( vPtU, dSqTol, nMaxInd, int( vPtU.size()) - 1, vInd))
          return false ;
-      vInd.push_back( vPtU.size() - 1) ;
+      vInd.push_back( int( vPtU.size()) - 1) ;
    }
 
   // ordino in senso crescente
@@ -825,7 +827,7 @@ PolyLine::RemoveAlignedPoints( double dToler)
 
   // se chiusa e almeno 4 punti rimasti, controllo allineamento dell'inizio con precedente e successivo rimasti
    if ( IsClosed()  &&  vInd.size() >= 4) {
-      if ( DistPointLine( vPtU[vInd[0]].first, vPtU[vInd[1]].first, vPtU[vInd[vInd.size()-2]].first).IsEpsilon( dToler)) {
+      if ( DistPointLine( vPtU[vInd[0]].first, vPtU[vInd[1]].first, vPtU[vInd[int(vInd.size())-2]].first).IsEpsilon( dToler)) {
          vInd.erase( vInd.begin()) ;
          vInd.back() = vInd.front() ;
       }
@@ -1445,7 +1447,7 @@ IsPointInsidePolyLine( const Point3d& ptP, const PolyLine& plPoly, double dToler
    }
   // Determino tangente di riferimento
    Vector3d vtTang ;
-   // se minima distanza nell'estremo iniziale del segmento
+  // se minima distanza nell'estremo iniziale del segmento
    if ( AreSamePointApprox( ptMinDist, prev( itMinDistEnd)->first)) {
      // direzione del segmento
       Vector3d vtCurrTg = itMinDistEnd->first - prev( itMinDistEnd)->first ;
@@ -1461,7 +1463,7 @@ IsPointInsidePolyLine( const Point3d& ptP, const PolyLine& plPoly, double dToler
       vtTang = vtPrevTg + vtCurrTg ;
       vtTang.Normalize() ;
    }
-   // se altrimenti minima distanza nell'estremo finale del segmento
+  // se altrimenti minima distanza nell'estremo finale del segmento
    else if ( AreSamePointApprox( ptMinDist, itMinDistEnd->first)) {
      // direzione del segmento
       Vector3d vtCurrTg = itMinDistEnd->first - prev( itMinDistEnd)->first ;
@@ -1477,7 +1479,7 @@ IsPointInsidePolyLine( const Point3d& ptP, const PolyLine& plPoly, double dToler
       vtTang = vtCurrTg + vtNextTg ;
       vtTang.Normalize() ;
    }
-   // altrimenti minima distanza con l'interno
+  // altrimenti minima distanza con l'interno
    else {
       vtTang = itMinDistEnd->first - prev( itMinDistEnd)->first ;
    }
@@ -1612,8 +1614,9 @@ AssociatePolyLinesMinDistPoints( const PolyLine& PL1, const PolyLine& PL2, PNTIV
    int nPnt2 = PL2.GetPointNbr() ;
    if ( nPnt1 == 0 || nPnt2 == 0)
       return false ;
-
-   bCommonInternalPoints = false ;   // indica la presenza di punti interni in comune tra le due polylines
+ 
+  // indica la presenza di punti interni in comune tra le due polylines
+   bCommonInternalPoints = false ;
 
    vPnt1.reserve( PL1.GetPointNbr()) ;
    Point3d ptP1 ;
@@ -1657,7 +1660,7 @@ AssociatePolyLinesMinDistPoints( const PolyLine& PL1, const PolyLine& PL2, PNTIV
          nMinJ = nLastJ ;
 
      // verifica se è un punto interno in comune con l'altra polyline  
-      if ( i < nTotP1 - 1  &&  dDist < EPS_SMALL  && abs( dMinDistPar - floor( dMinDistPar + 0.5)) < EPS_SMALL)
+      if ( i < nTotP1 - 1  &&  dDist < EPS_SMALL  &&  abs( dMinDistPar - floor( dMinDistPar + 0.5)) < EPS_SMALL)
          bCommonInternalPoints = true ;
 
       vPnt1[i].second = nMinJ ;
@@ -1687,12 +1690,328 @@ AssociatePolyLinesMinDistPoints( const PolyLine& PL1, const PolyLine& PL2, PNTIV
       if ( nMinI < nLastI)
          nMinI = nLastI ;
 
-      if ( j < nTotP2 - 1  &&  dDist < EPS_SMALL && abs( dMinDistPar - floor( dMinDistPar + 0.5)) < EPS_SMALL)
+      if ( j < nTotP2 - 1  &&  dDist < EPS_SMALL &&  abs( dMinDistPar - floor( dMinDistPar + 0.5)) < EPS_SMALL)
          bCommonInternalPoints = true ;
 
       vPnt2[j].second = nMinI ;
       nLastI = nMinI ;
    }
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+bool 
+MatchPolyLinesAddingPoints( const PolyLine& PL1, const PolyLine& PL2, int nType, PNTIVECTOR& vPnt1, PNTIVECTOR& vPnt2)
+{
+  // prima trovo le associazioni senza aggiunte di punti
+   Point3d ptP2 ;
+   CurveComposite cc1, cc2 ;
+   cc1.FromPolyLine( PL1) ;
+   cc2.FromPolyLine( PL2) ;
+   int nPnt1 = PL1.GetPointNbr() ;
+   int nPnt2 = PL2.GetPointNbr() ;
+   vector<POINTU> vMatch2 ;
+   PL2.GetFirstPoint( ptP2) ;
+   while ( PL2.GetNextPoint( ptP2, true)) {
+      DistPointCurve dpc( ptP2, cc1, false) ;
+      int nFlag = 0 ;
+      double dParam ; dpc.GetParamAtMinDistPoint( 0, dParam, nFlag) ;
+      Point3d ptJoint ; dpc.GetMinDistPoint( 0, ptJoint, nFlag) ;
+      vMatch2.emplace_back( ptJoint, dParam) ;
+   }
+
+   int nAtStart2 = 0 ;  // match ripetuti dalla curva U0 allo start della curva U1
+   int nAtEnd2 = 0 ;
+   Point3d ptP1 ; PL1.GetFirstPoint( ptP1) ;
+   int c = 0 ;
+   int nRep1 = 0 ;  // match interni consecutivi uguali di punti della curva U0 con punti della curva U1
+   int nRep2 = 0 ;
+   double dLastParamMatch = 0 ;
+   Point3d ptLastPointMatch = ptP1 ;
+   BOOLVECTOR vbRep1( nPnt1) ;
+   INTVECTOR vnAddedOrNextIsRep1 ;  // 0 non Rep, 1 Rep, 2 Next is Rep ;
+   DBLVECTOR vdSplit1 ;
+   DBLVECTOR vdMatch1 ;
+   fill( vbRep1.begin(), vbRep1.end(), false) ;
+   while ( PL1.GetNextPoint( ptP1, true)) {
+      // devo salvarmi se matcho più punti con lo start o l'end della curva totale
+      DistPointCurve dpc( ptP1, cc2, false) ;
+      int nFlag = 0 ;
+      double dParam ; dpc.GetParamAtMinDistPoint( 0, dParam, nFlag) ;
+      Point3d ptJoint ; dpc.GetMinDistPoint( 0, ptJoint, nFlag) ;
+      vdMatch1.push_back( dParam) ;
+      if ( dParam < EPS_SMALL ) {
+         ++nAtStart2 ;
+         vbRep1[c] = true ;
+         ++c ;
+         continue ;
+      }
+      else if ( dParam > nPnt2 - EPS_SMALL ) {
+         vbRep1[c] = nAtEnd2 == 0 ? false : true ;
+         vbRep1.back() = true ;
+         ++ nAtEnd2 ;
+         ++c ;
+         continue ;
+      }
+      if ( dParam <= dLastParamMatch  ||  AreSamePointApprox( ptJoint, ptLastPointMatch)) {
+         dParam = dLastParamMatch ;
+         vbRep1[c] = true ;
+         ++ nRep1 ;
+         ++c ;
+         continue ;
+      }
+      else {
+         dLastParamMatch = dParam ;
+         ptLastPointMatch = ptJoint ;
+         // se sono già troppo vicino ad un split esistente allora non faccio nulla
+         if ( abs(dParam - round( dParam)) < 100 * EPS_PARAM) {
+            ++c ;
+            continue ;
+         }
+         vdSplit1.push_back( dParam) ;
+         // verifico se ho un match per questo punto
+         // in tal caso vuol dire che sto creando una ripetizione nRep1
+         int nCase = 0 ;
+         for ( int j = 0 ; j < int( vMatch2.size()) ; ++j) {
+            if ( abs(vMatch2[j].second - (c + 1)) < EPS_SMALL) {
+               // devo però verificare che non ci sia un match successivo, perché in quel caso non ho una ripetizione
+               bool bFoundMatch = false ;
+               for ( int z = int( vMatch2[j].second) ; z < int( vdMatch1.size()) ; ++z) {
+                  if ( abs( vdMatch1[z] - ( j + 1 )) < EPS_SMALL ) {
+                     bFoundMatch = true ;
+                     break ;
+                  }
+               }
+
+               if ( ! bFoundMatch) {
+                  ++ nRep2 ;
+                  // capisco se il punto è rep o se lo è il suo successivo
+                  if ( j + 1 < dParam)
+                     nCase = 1 ;
+                  else
+                     nCase = 2 ;
+               }
+               break ;
+            }
+         }
+         vnAddedOrNextIsRep1.push_back( nCase) ;
+      }
+      ++c ;
+   }
+   int nAtStart1 = 0 ;
+   int nAtEnd1 = 0 ;
+   PL2.GetFirstPoint( ptP2) ;
+   c = 0 ;
+   dLastParamMatch = 0 ;
+   ptLastPointMatch = ptP2 ;
+   INTVECTOR vnAddedOrNextIsRep2 ;  // 0 non Rep, 1 Rep, 2 Next is Rep ;
+   DBLVECTOR vdSplit2 ;
+   BOOLVECTOR vbRep2( nPnt2) ;
+   fill( vbRep2.begin(), vbRep2.end(), false) ;
+   while ( PL2.GetNextPoint( ptP2, true)) {
+      DistPointCurve dpc( ptP2, cc1, false) ;
+      int nFlag = 0 ;
+      double dParam ; dpc.GetParamAtMinDistPoint( 0, dParam, nFlag) ;
+      Point3d ptJoint ; dpc.GetMinDistPoint( 0, ptJoint, nFlag) ;
+      if ( dParam < EPS_SMALL ) {
+         ++nAtStart1 ;
+         vbRep2[c] = true ;
+         ++c ;
+         continue ;
+      }
+      else if ( dParam > nPnt1 - EPS_SMALL ) {
+         vbRep2[c] = ( nAtEnd1 == 0 ? false : true) ;
+         vbRep2.back() = true ;
+         ++ nAtEnd1 ;
+         ++c ;
+         continue ;
+      }
+      if ( dParam <= dLastParamMatch  ||  AreSamePointApprox( ptJoint, ptLastPointMatch)) {
+         dParam = dLastParamMatch ;
+         vbRep2[c] = true ;
+         ++ nRep2 ;
+         ++c ;
+         continue ;
+      }
+      else {
+         dLastParamMatch = dParam ;
+         ptLastPointMatch = ptJoint ;
+        // se sono troppo vicino ad uno split esistente allora non faccio nulla
+         if ( abs( dParam - round( dParam)) < 100 * EPS_PARAM) {
+            ++c ;
+            continue ;
+         }
+         vdSplit2.push_back( dParam) ;
+        // verifico se ho un match per questo punto
+        // in tal caso vuol dire che sto creando una ripetizione nRep0
+         int nCase = 0 ;
+         for ( int j = 0 ; j < int( vdMatch1.size()) ; ++j) {
+            if ( abs( vdMatch1[j] - (c + 1)) < EPS_SMALL) {
+              // devo però verificare che non ci sia un match successivo, perché in quel caso non ho una ripetizione
+               bool bFoundMatch = false ;
+               for ( int z = int( vdMatch1[j]) ; z < int( vMatch2.size()) ; ++z) {
+                  if ( abs( vMatch2[z].second - ( j + 1 )) < EPS_SMALL) {
+                     bFoundMatch = true ;
+                     break ;
+                  }
+               }
+
+               if ( ! bFoundMatch) {
+                  ++nRep1 ;
+                 // capisco se il punto è rep o se lo è il suo successivo
+                  if ( j + 1 < dParam)
+                     nCase = 1 ;
+                  else
+                     nCase = 2 ;
+               }
+               break ;
+            }
+         }
+         vnAddedOrNextIsRep2.push_back( nCase) ;
+      }
+      ++c ;
+   }
+
+  // applico effettivamente gli split e aggiungo gli elementi ai vettori vbRep
+   int nUnit = 0 ;
+   if ( ! vdSplit1.empty())
+      nUnit = int( vdSplit1.back()) ;
+   for ( int z = int( vdSplit1.size()) - 1 ; z >= 0 ; --z) {
+      double dSplit = vdSplit1[z] ;
+      int nSplit = int( dSplit) ;
+     // se sto cercando di fare uno split sulla stessa curva che ho già splittato al passaggio precedente allora devo
+     // riscalare
+      if ( nSplit == nUnit  &&  z < int( vdSplit1.size()) - 1) {
+         dSplit = nSplit + ( dSplit - nSplit) / ( vdSplit1[z+1] - nSplit) ;
+      }
+      nUnit = nSplit ;
+      cc2.AddJoint( dSplit) ;
+      switch ( vnAddedOrNextIsRep1[z]) {
+      case 0 :
+         if ( vbRep2[nSplit])
+            ++ nRep2 ;
+        // di default aggiungerei false, ma se il successivo è già un Rep allora anche questo deve esserlo
+         vbRep2.insert( vbRep2.begin() + nSplit, vbRep2[nSplit]) ;
+         break ;
+      case 1 :
+         vbRep2.insert( vbRep2.begin() + nSplit, true) ;
+         break ;
+      case 2 :
+         if ( vbRep2[nSplit])
+            --nRep2 ;
+         else
+            vbRep2[nSplit] = true ;
+         vbRep2.insert( vbRep2.begin() + nSplit, false) ;
+         break ;
+      }
+   }
+   if ( ! vdSplit2.empty())
+      nUnit = int( vdSplit2.back()) ;
+   for ( int z = int( vdSplit2.size()) - 1 ; z >= 0 ; --z) {
+      double dSplit = vdSplit2[z] ;
+      int nSplit = int( dSplit) ;
+     // se sto cercando di fare uno split sulla stessa curva che ho già splittato al passaggio precedente allora devo
+     // riscalare
+      if ( nSplit == nUnit  &&  z < int( vdSplit2.size()) - 1) {
+         dSplit = nSplit + ( dSplit - nSplit) / (vdSplit2[z+1] - nSplit) ;
+      }
+      nUnit = nSplit ;
+      cc1.AddJoint( dSplit) ;
+      switch ( vnAddedOrNextIsRep2[z]) {
+      case 0 :
+         if ( vbRep1[nSplit])
+            ++ nRep1 ;
+        // di default aggiungerei false, ma se il successivo è già un Rep allora anche questo deve esserlo
+         vbRep1.insert( vbRep1.begin() + nSplit, vbRep1[nSplit]) ;
+         break ;
+      case 1 :
+         vbRep1.insert( vbRep1.begin() + nSplit, true) ;
+         break ;
+      case 2 :
+         if ( vbRep1[nSplit])
+            -- nRep1 ;
+         else
+            vbRep1[nSplit] = true ;
+         vbRep1.insert( vbRep1.begin() + nSplit, false) ;
+         break ;
+      }
+   }
+
+   nPnt1 = cc1.GetCurveCount() ;
+   nPnt2 = cc2.GetCurveCount() ;
+   //aggiusto i vettori delle ripetizioni in modo in modo che non arrivino mai ad essere contemporaneamente true
+   int nAddedSpan = 0 ;
+   int nCrv1 = 0 ;
+   int nCrv2 = 0 ;
+   while ( nAddedSpan < nPnt1 + nAtStart1 + nAtEnd1 + nRep2) {
+      if ( nCrv1 >= nPnt1)
+         nCrv1 = nPnt1 - 1 ;
+      if ( nCrv2 >= nPnt2)
+         nCrv2 = nPnt2 - 1 ;
+      bool bRep1 = vbRep1[nCrv1] ;
+      bool bRep2 = vbRep2[nCrv2] ;
+      if ( bRep1  &&  bRep2) {
+         vbRep1[nCrv1] = false ;
+         bRep1 = false ;
+         vbRep2[nCrv2] = false ;
+         bRep2 = false ;
+         -- nRep1 ;
+         -- nRep2 ;
+      }
+      if ( ! bRep1  ||  nCrv2 == nPnt2 - 1)
+         ++ nCrv2 ;
+      if ( ! bRep2  ||  nCrv1 == nPnt1 - 1)
+         ++ nCrv1 ;
+      ++ nAddedSpan ;
+   }
+  // se non sono arrivato all'ultima curva su U0 o U1 vuol dire che ho creato delle ripetizioni che non ho contato prima
+   if ( nCrv2 < nPnt2) {
+      nRep2 += nPnt2 - nCrv2 ;
+      for ( int z = int( vbRep2.size()) - 1 ; z >= nCrv2 ; --z)
+         vbRep2[z] = true ;
+   }
+   if( nCrv1 < nPnt1) {
+      nRep1 += nPnt1 - nCrv1 ;
+      for ( int z = int( vbRep1.size()) - 1 ; z >= nCrv1 ; --z)
+         vbRep1[z] = true ;
+   }
+
+  // trovo il numero di span che dovrà avere la superficie
+  // ( numero di sottocurve che compongono la U0 + tutte le ripetizioni dei match di punti della curva U1 con i punti di U0)
+   int nPnt = nPnt1 + nAtStart1 + nAtEnd1 + nRep2 ;
+
+   if ( nPnt != nPnt2 + nAtStart2 + nAtEnd2 + nRep1)
+      LOG_DBG_ERR( GetEGkLogger(), "There could be an error in the creation of a ruled surface in mode RLT_B_MINDIST_PLUS") ;
+
+  // aggiungo i punti di controllo scorrendo in contemporanea le due curve
+   nAddedSpan = 0 ;
+   nCrv1 = 0 ;
+   nCrv2 = 0 ;
+   bool bLast1 = false ;
+   bool bLast2 = false ;
+   while ( nAddedSpan < nPnt) {
+      if ( nCrv1 >= nPnt1) {
+         nCrv1 = nPnt1 - 1 ;
+         bLast1 = true ;
+      }
+      if ( nCrv2 >= nPnt2) {
+         nCrv2 = nPnt2 - 1 ;
+         bLast2 = true ;
+      }
+      bool bRep1 = vbRep1[nCrv1] ;
+      bool bRep2 = vbRep2[nCrv2] ;
+      const ICurve* pSubCrv1 = cc1.GetCurve( nCrv1) ;
+      Point3d ptStart1 ; pSubCrv1->GetStartPoint( ptStart1) ;
+      vPnt1.emplace_back( ptStart1, nAddedSpan) ;
+      const ICurve* pSubCrv2 = cc2.GetCurve( nCrv2) ;
+      Point3d ptStart2 ; pSubCrv2->GetStartPoint( ptStart2) ;
+      vPnt2.emplace_back( ptStart2, nAddedSpan) ;
+      if ( ! bRep2)
+         ++ nCrv1 ;
+      if ( ! bRep1)
+         ++ nCrv2 ;
+      ++ nAddedSpan ;
+   }
 
    return true ;
-}
+}

ProjectCurveSurf.cpp

@@ -1,7 +1,7 @@
 //----------------------------------------------------------------------------
-//  EgalTech 2023-2023
+//  EgalTech 2023-2025
 //----------------------------------------------------------------------------
-// File : ProjectCurveSurfTm.cpp    Data : 16.11.23        Versione : 2.5kh3
+// File : ProjectCurveSurfTm.cpp    Data : 29.08.25        Versione : 2.7h2
 // Contenuto : Implementazione funzioni proiezione curve su superficie Trimesh.
 //
 //
@@ -40,6 +40,13 @@ const int P5AX_CONC = 3 ;              // in angolo concavo
 const int P5AX_BEFORE_CONC = 4 ;       // adiacente ad angolo concavo
 const int P5AX_AFTER_CONC = 5 ;        // adiacente ad angolo concavo
 
+//----------------------------------------------------------------------------
+static double
+GetSurfBezierTol( double dLinTol)
+{
+   return max( dLinTol / 10, EPS_SMALL) ;
+}
+
 //----------------------------------------------------------------------------
 static bool
 PointsInTolerance( const PNT5AXVECTOR& vPt5ax, int nPrec, int nCurr, int nNext, double dSqTol)
@@ -154,7 +161,8 @@ RemovePointsInExcess( PNT5AXVECTOR& vPt5ax, double dLinTol, double dMaxSegmLen,
    return true ;
 }
 
-typedef std::vector<IntersParLinesSurfTm*> INTPARLINESTMPVECTOR ;   // vettore di oggetti intersezione massiva rette parallele SurfTM
+// --- vettore di oggetti intersezione massiva rette parallele SurfTM --------
+typedef std::vector<IntersParLinesSurfTm*> INTPARLINESTMPVECTOR ;
 
 //----------------------------------------------------------------------------
 static bool
@@ -230,8 +238,11 @@ ProjectCurveOnSurf( const ICurve& crCrv, const CISURFPVECTOR& vpSurf, const Vect
          pSurfTm = GetBasicSurfTriMesh( vpSurf[i]) ;
          break ;
       case SRF_BEZIER :
-         pSurfTm = GetBasicSurfBezier( vpSurf[i])->GetAuxSurf() ;
-         break ;
+       { double dOldLinTol = GetSurfBezierAuxSurfRefinedTol() ;
+         SetSurfBezierAuxSurfRefinedTol( GetSurfBezierTol( dLinTol)) ;
+         pSurfTm = GetBasicSurfBezier( vpSurf[i])->GetAuxSurfRefined() ;
+         SetSurfBezierAuxSurfRefinedTol( dOldLinTol) ;
+       } break ;
       case SRF_FLATRGN :
          pSurfTm = GetBasicSurfFlatRegion( vpSurf[i])->GetAuxSurf() ;
          break ;
@@ -440,8 +451,11 @@ ProjectCurveOnSurf( const ICurve& crCrv, const CISURFPVECTOR& vpSurf, const IGeo
          pSurfTm = GetBasicSurfTriMesh( vpSurf[i]) ;
          break ;
       case SRF_BEZIER :
-         pSurfTm = GetBasicSurfBezier( vpSurf[i])->GetAuxSurf() ;
-         break ;
+       { double dOldLinTol = GetSurfBezierAuxSurfRefinedTol() ;
+         SetSurfBezierAuxSurfRefinedTol( GetSurfBezierTol( dLinTol)) ;
+         pSurfTm = GetBasicSurfBezier( vpSurf[i])->GetAuxSurfRefined() ;
+         SetSurfBezierAuxSurfRefinedTol( dOldLinTol) ;
+       } break ;
       case SRF_FLATRGN :
          pSurfTm = GetBasicSurfFlatRegion( vpSurf[i])->GetAuxSurf() ;
          break ;
@@ -593,8 +607,11 @@ ProjectCurveOnSurf( const ICurve& crCrv, const CISURFPVECTOR& vpSurf, const ICur
          pSurfTm = GetBasicSurfTriMesh( vpSurf[i]) ;
          break ;
       case SRF_BEZIER :
-         pSurfTm = GetBasicSurfBezier( vpSurf[i])->GetAuxSurf() ;
-         break ;
+       { double dOldLinTol = GetSurfBezierAuxSurfRefinedTol() ;
+         SetSurfBezierAuxSurfRefinedTol( GetSurfBezierTol( dLinTol)) ;
+         pSurfTm = GetBasicSurfBezier( vpSurf[i])->GetAuxSurfRefined() ;
+         SetSurfBezierAuxSurfRefinedTol( dOldLinTol) ;
+       } break ;
       case SRF_FLATRGN :
          pSurfTm = GetBasicSurfFlatRegion( vpSurf[i])->GetAuxSurf() ;
          break ;
@@ -804,8 +821,11 @@ ProjectCurveOnSurf( const ICurve& crCrv, const CISURFPVECTOR& vpSurf, const ISur
          pSurfTm = GetBasicSurfTriMesh( vpSurf[i]) ;
          break ;
       case SRF_BEZIER :
-         pSurfTm = GetBasicSurfBezier( vpSurf[i])->GetAuxSurf() ;
-         break ;
+       { double dOldLinTol = GetSurfBezierAuxSurfRefinedTol() ;
+         SetSurfBezierAuxSurfRefinedTol( GetSurfBezierTol( dLinTol)) ;
+         pSurfTm = GetBasicSurfBezier( vpSurf[i])->GetAuxSurfRefined() ;
+         SetSurfBezierAuxSurfRefinedTol( dOldLinTol) ;
+       } break ;
       case SRF_FLATRGN :
          pSurfTm = GetBasicSurfFlatRegion( vpSurf[i])->GetAuxSurf() ;
          break ;
@@ -824,8 +844,11 @@ ProjectCurveOnSurf( const ICurve& crCrv, const CISURFPVECTOR& vpSurf, const ISur
       pRefTm = GetBasicSurfTriMesh( &sfRef) ;
       break ;
    case SRF_BEZIER :
-      pRefTm = GetBasicSurfBezier( &sfRef)->GetAuxSurf() ;
-      break ;
+    { double dOldLinTol = GetSurfBezierAuxSurfRefinedTol() ;
+      SetSurfBezierAuxSurfRefinedTol( GetSurfBezierTol( dLinTol)) ;
+      pRefTm = GetBasicSurfBezier(  &sfRef)->GetAuxSurfRefined() ;
+      SetSurfBezierAuxSurfRefinedTol( dOldLinTol) ;
+    } break ;
    case SRF_FLATRGN :
       pRefTm = GetBasicSurfFlatRegion( &sfRef)->GetAuxSurf() ;
       break ;

SurfBezier.cpp

@@ -41,23 +41,59 @@
 #include "/EgtDev/Include/EGkGeoPoint3d.h"
 #include "/EgtDev/Include/EGkIntervals.h"
 #include "/EgtDev/Extern/Eigen/Dense"
-
+#include "/EgtDev/Include/EgtPerfCounter.h"
 
 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) ;
+//----------------------------------------------------------------------------
+struct PairHashIntInt {
+   size_t operator()(const pair<int, int>& key) const {
+      size_t h1 = hash<int>{}( key.first) ;
+      size_t h2 = hash<int>{}( key.second) ;
+      return h1 ^ ( h2 << 1) ;  // Combine hashes
+   }
+};
+static unordered_map<INTINT, DBLVECTOR, PairHashIntInt>  s_mBernCache ;  // mappa dei polinomi di bernstein
 
-static unordered_map<INTINT, DBLVECTOR, PairHashIntInt>  m_mBernCache ;  // mappa dei polinomi di bernstein
+//----------------------------------------------------------------------------
+static double s_dAuxSurfTol = 200 * EPS_SMALL ;
+static double s_dAuxSurfRefinedTol = 50 * EPS_SMALL ;
+
+//----------------------------------------------------------------------------
+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 +401,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 +779,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 +874,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 +1724,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 +1752,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 +1779,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,30 +1791,25 @@ 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() ;
+   //debug
+   static PerformanceCounter Counter;
+   static double dTotCount = 0. ;
+   Counter.Start();
+   //debug
    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
-      }
-      if ( ! Tree.GetPolygons( vvPL, vvPL3d, m_mCCEdge, m_vCCLoop))
-         continue ;
+      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
 
      // aggiorno la chiusura della superficie
@@ -1759,6 +1819,12 @@ SurfBezier::GetApproxSurf( double dTol, double dSideMin) const
      // UpdateEdgesFromTree( Tree) ; // ora viene fatto nella GetPolygons
    }
 
+   //debug
+   dTotCount += Counter.Stop() ;
+   string sOut = "CalcTree=" + ToString( dTotCount, 3) ;
+   LOG_INFO( GetEGkLogger(), sOut.c_str())
+   //debug
+
   //// per usare i polygon basic//////////////////////
    //for (int k = 0 ; k < (int)vPL.size(); ++k) {
    //   vvPL.emplace_back() ;
@@ -1787,18 +1853,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 +1876,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 +1891,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 +1901,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 +1916,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 +1937,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 +2095,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) {
@@ -2125,7 +2183,7 @@ bool
 SurfBezier::IncreaseUV( double& dUV, double dxy, bool bUOrV, double* dUVCopy, bool bModifyOrig) const
 {
    double dUVTest ;
-   if (dUVCopy == nullptr  && !bModifyOrig)
+   if ( dUVCopy == nullptr  &&  ! bModifyOrig)
       return false ;
    if ( dUVCopy != nullptr) {
       *dUVCopy = dUV + dxy ;
@@ -2138,15 +2196,15 @@ SurfBezier::IncreaseUV( double& dUV, double dxy, bool bUOrV, double* dUVCopy, bo
    
    if ( bUOrV) {
       if ( dUVTest < 0)
-            dUVTest = 0 ;
+         dUVTest = 0 ;
       else if ( dUVTest > m_nSpanU * SBZ_TREG_COEFF )
-            dUVTest = m_nSpanU * SBZ_TREG_COEFF ;
+         dUVTest = m_nSpanU * SBZ_TREG_COEFF ;
    }
    else {
       if ( dUVTest < 0)
-            dUVTest = 0 ;
+         dUVTest = 0 ;
       else if ( dUVTest > m_nSpanV * SBZ_TREG_COEFF )
-            dUVTest = m_nSpanV * SBZ_TREG_COEFF ;
+         dUVTest = m_nSpanV * SBZ_TREG_COEFF ;
    }
    if ( bModifyOrig)
       dUV  = dUVTest ;
@@ -2812,7 +2870,6 @@ SurfBezier::UnprojectPointFromStm( int nT, const Point3d& ptI, Point3d& ptSP, in
    // se sono in polo e mi hanno passato un punto precedente allora devo prendere il triangolo di quel punto
    bool bIsPole = false ;
    bool bNearPole = false ;  // devo capire se il triangolo di riferimento ha un vertice in un polo
-   int nInters = 0 ;
    INTVECTOR vInters(4) ;
    fill( vInters.begin(), vInters.end(), 0) ;
    // se il vettore dei poli non è stato riempito vuol dire che quando è stata creata la superficie non è stata chiamata la funzione CalcPoles
@@ -2825,16 +2882,14 @@ 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)
                   bIsPole = true ;
-               nInters += vInters[c] ;
             }
             else {
                vInters[c] = m_mCCEdge[c][i]->IsPointOn(ptI) ? 1 : 0 ;
-               nInters += vInters[c] ;
             }
          }
       }
@@ -3177,8 +3232,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 +4246,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 +4274,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 +5199,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 +5366,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)
@@ -5352,3 +5407,11 @@ SurfBezier::CreateBySetOfCurves( const ICURVEPOVECTOR& vCrvBez, bool bReduceToDe
    }
    return true ;
 }
+
+//----------------------------------------------------------------------------
+PNTVECTOR
+SurfBezier::GetAllControlPoints( void) const
+{
+   PNTVECTOR vPntCtrl = m_vPtCtrl ;
+   return vPntCtrl ;
+}

SurfBezier.h

@@ -22,17 +22,6 @@
 #include "/EgtDev/Include/EGkSurfBezier.h"
 #include "/EgtDev/Include/EGkGeoCollection.h"
 
-using namespace std ;
-class Tree ;
-
-struct PairHashIntInt {
-   std::size_t operator()(const std::pair<int, int>& key) const {
-      std::size_t h1 = std::hash<int>{}(key.first);
-      std::size_t h2 = std::hash<int>{}(key.second);
-      return h1 ^ (h2 << 1);  // Combine hashes
-   }
-};
-
 //----------------------------------------------------------------------------
 class SurfBezier : public ISurfBezier, public IGeoObjRW
 {
@@ -82,10 +71,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 ;
@@ -157,6 +146,7 @@ class SurfBezier : public ISurfBezier, public IGeoObjRW
       bool CreateByPointCurve( const Point3d& pt, const ICurve* pCurve) override ;
       bool CreateByTwoCurves( const ICurve* pCurve1, const ICurve* pCurve2, int nType) override ;
       bool CreateBySetOfCurves( const ICURVEPOVECTOR& vCrvBez, bool bReduceToDeg3) override ;
+      PNTVECTOR GetAllControlPoints( void) const ;
 
    public :                      // IGeoObjRW
       int  GetNgeId( void) const override ;
@@ -168,7 +158,7 @@ class SurfBezier : public ISurfBezier, public IGeoObjRW
 
    public :
       SurfBezier( void) ;
-      SurfBezier( const SurfBezier& sbSrc)
+      SurfBezier( const SurfBezier& sbSrc) : m_pSTM( nullptr), m_pSTMRefined( nullptr), m_pTrimReg(nullptr)
                   { if ( ! CopyFrom( sbSrc))
                        LOG_ERROR( GetEGkLogger(), "SurfBezier : copy constructor error") }
       SurfBezier& operator =( const SurfBezier& sbSrc)
@@ -238,7 +228,7 @@ class SurfBezier : public ISurfBezier, public IGeoObjRW
       SurfFlatRegion* m_pTrimReg ;         // eventuale regione di trim
       int             m_nTempProp[2] ;     // vettore proprietà temporanee
       double          m_dTempParam[2] ;    // vettore parametri temporanei
-      mutable vector<ICRVCOMPOPOVECTOR> m_mCCEdge ; // vettore dei vettori che contengono le curve compo degli edge della superficie nello spazio 3D
+      mutable std::vector<ICRVCOMPOPOVECTOR> m_mCCEdge ; // vettore dei vettori che contengono le curve compo degli edge della superficie nello spazio 3D
       mutable ICRVCOMPOPOVECTOR m_vCCLoop ;         // vettore dei loop della superficie trimmata
       mutable int            m_nIsPlanar ;          // enum che indica se la superficie è piana ( -1, non è stato calcolato)
       mutable DBLVECTOR      m_vBernU ;

SurfFlatRegion.cpp

@@ -827,7 +827,7 @@ SurfFlatRegion::ConvertArcsToBezierCurves( void)
   // ciclo sui loop
    for ( auto& pLoop : m_vpLoop) {
       if ( pLoop->GetType() == CRV_ARC) {
-         ICurve* pCrvNew = ArcToBezierCurve( pLoop) ;
+         ICurve* pCrvNew = ArcToBezierCurve( GetCurveArc(pLoop)) ;
          if ( pCrvNew == nullptr)
             return false ;
          delete pLoop ;

SurfTriMesh.cpp

@@ -3637,6 +3637,18 @@ SurfTriMesh::Invert( void)
    for ( int i = 0 ; i < GetTriangleSize() ; ++ i)
       InvertTriangle( i) ;
 
+  // se bordi della sfaccettatura validi
+   if ( m_bFacEdged) {
+      for ( int nE = 0 ; nE < GetEdgeCount() ; ++ nE) {
+        // inversione dei vertici nel vettore delle sfaccettature
+         swap( m_vFacEdge[nE].nIdVert[0], m_vFacEdge[nE].nIdVert[1]) ;
+        // inversione delle adiacenze nel vettore delle sfaccettature
+         swap( m_vFacEdge[nE].nIdFacAdj[0], m_vFacEdge[nE].nIdFacAdj[1]) ;
+        // inversione dell'angolo interno
+         m_vFacEdge[nE].dIntAng *= -1. ;
+      }
+   }
+
    return true ;
 }
 

SurfTriMeshOffset.cpp

@@ -15,8 +15,21 @@
 //--------------------------- Include ----------------------------------------
 #include "stdafx.h"
 #include "VolZmap.h"
+#include "SurfTriMesh.h"
+#include "EgtDev/Include/EGkDistPointSurfTm.h"
 #include "\EgtDev\Include\EGkSurfTriMeshAux.h"
 
+#define DEBUG 0
+#if DEBUG
+   #include "/EgtDev/Include/EGkGeoObjSave.h"
+   #include "/EgtDev/Include/EGkGeoPoint3d.h"
+   #include "/EgtDev/Include/EGkGeoVector3d.h"
+   #include "/EgtDev/Include/EGkStmStandard.h"
+   #include "/EgtDev/Include/EgtPerfCounter.h"
+   std::vector<IGeoObj*> VT ;
+   std::vector<Color> VC ;
+#endif
+
 using namespace std ;
 
 //----------------------------------------------------------------------------
@@ -73,20 +86,131 @@ CreateSurfTriMeshesOffset( const CISURFTMPVECTOR& vStm, double dOffs, double dPr
    if ( vStm.empty())
       return nullptr ;
   // controllo sul valore di tolleranza lineare
-   double dMyPrec = max( dPrec, 100 * EPS_SMALL) ;
+   double dMyPrec = max( dPrec, 100. * EPS_SMALL) ;
   // --- NB. ( Il valore di Offset deve essere maggiore di 10 * EPS_SMALL in valore assoluto)
   // Nel caso sia minore, restituisco semplicemente la somma delle superfici
   // ( questo valore serve per rimanere coerente con l'Offset delle curve)
-   if ( abs( dOffs) < 10 * EPS_SMALL)
+   if ( abs( dOffs) < 10. * EPS_SMALL)
       return SumStm( vStm) ;
 
   // creo lo Zmap associato alle superfici TriMesh
-   PtrOwner<IVolZmap> pVolZmap( CreateVolZmap()) ;
-   if ( IsNull( pVolZmap)  ||  ! pVolZmap->CreateFromTriMeshOffset( vStm, dOffs, dMyPrec, nType))
+   VolZmap OneVolZmap ;
+   if ( ! OneVolZmap.CreateFromTriMeshOffset( vStm, dOffs, dMyPrec, nType))
+      return nullptr ;
+   if ( ! OneVolZmap.IsValid())
       return nullptr ;
 
-  // restituisco la superficie TriMesh
-   return ( pVolZmap->GetSurfTriMesh()) ;
+  // recupero le superfici aperte
+   CISURFTMPVECTOR vStmOpen ;
+   for ( const ISurfTriMesh* pStm : vStm) {
+      if ( pStm != nullptr  &&  pStm->IsValid() &&   ! pStm->IsClosed())
+         vStmOpen.emplace_back( pStm) ;
+   }
+
+  // --- se non ho superfici aperte
+   if ( vStmOpen.empty()) {
+     // restituisco la superficie TriMesh di Offset
+      return ( OneVolZmap.GetSurfTriMesh()) ;
+   }
+
+  // --- se ho delle superfici aperte
+  // lo Zmap creato è orientato e definisce una superficie chiusa; devo rimuovere i triangoli in eccesso
+
+  // anzitutto controllo che lo Zmap sia valido
+   if ( ! OneVolZmap.IsValid())
+      return nullptr ;
+
+  // inzializzo la superficie TriMesh da restituire
+   PtrOwner<SurfTriMesh> pStm( CreateBasicSurfTriMesh()) ;
+   if ( IsNull( pStm)  ||  ! pStm->Init( 3, 1))
+      return nullptr ;
+   PointGrid3d VertGrid ; VertGrid.Init( 50000) ;
+
+  // tolleranza di vicinanza alla superficie
+   double dTolDist = 30. * EPS_SMALL ;
+
+   #if DEBUG
+      VT.emplace_back( OneVolZmap.Clone()) ;
+      VC.emplace_back( BLACK) ;
+   #endif
+
+  // ciclo lungo i blocchi dello ZMap
+   for ( int nB = 0 ; nB < OneVolZmap.GetBlockCount() ; ++ nB) {
+
+     // recupero i triangoli
+      TRIA3DEXVECTOR vTria, vTriaSafe ;
+      OneVolZmap.GetBlockTriangles( nB, vTria) ;
+
+     // un triangolo viene ritenuto valido se è non è troppo vicino ( dOffs) alle superfici aperte
+      for ( int nT = 0 ; nT < int( vTria.size()) ; ++ nT) {
+
+        // recupero il triangolo
+         Triangle3dEx& Tria = vTria[nT] ;
+
+        // scorro le superficie aperte
+         bool bInsert = true ;
+         for ( int nS = 0 ; bInsert && nS < int( vStm.size()) ; ++ nS) {
+
+           // controllo se posso inserirlo
+            vector<DistPointSurfTm> vDistPtStm ;
+            for ( int i = 0 ; i < 3 && bInsert ; ++ i) {
+               double dDist = 0. ;
+               vDistPtStm.emplace_back( DistPointSurfTm( Tria.GetP( i), *vStm[nS])) ;
+               bInsert = ( vDistPtStm.back().GetDist( dDist) && dDist > abs( dOffs) - dTolDist) ;
+            }
+           // se il triangolo è al più a distanza di |dOffs| - dTolDist
+            if ( bInsert) {
+              // recupero i triangoli a distanza minima dai vertici del triangolo corrente
+               bool bPerp = true ;
+               for ( int i = 0 ; i < 3 && bPerp ; ++ i) {
+                  INTVECTOR vTria ;
+                  vDistPtStm[i].GetMinDistTriaIndices( vTria) ;
+                  for ( int j = 0 ; j < int( vTria.size()) && bPerp ; ++ j) {
+                     Triangle3d TriaCloser ;
+                     vStm[nS]->GetTriangle( vTria[j], TriaCloser) ;
+                     bPerp = ( abs( Tria.GetN() * TriaCloser.GetN()) < dTolDist) ;
+                  }
+               }
+              // se tutti i triangoli a distanza minima sono perpendicolari, allora non lo inserisco
+               bInsert = ( ! bPerp) ;
+            }
+         }
+
+        // se triangolo da inserire
+         if ( bInsert)
+            vTriaSafe.emplace_back( Tria) ;
+
+         #if DEBUG
+            ICurveComposite* pCompo = CreateCurveComposite() ;
+            pCompo->AddPoint( Tria.GetP( 0)) ;
+            pCompo->AddLine( Tria.GetP( 1)) ;
+            pCompo->AddLine( Tria.GetP( 2)) ;
+            pCompo->Close() ;
+            Color myCol = ( bInsert ? Color( 0., 1., 0., .5) : Color( 1., 0., 0., .5)) ;
+            VT.emplace_back( CloneCurveComposite( pCompo)) ;
+            VC.emplace_back( myCol) ;
+            ISurfFlatRegion* pSfrTria = CreateSurfFlatRegion() ;
+            pSfrTria->AddExtLoop( pCompo) ;
+            VT.emplace_back( pSfrTria) ;
+            VC.emplace_back( myCol) ;
+         #endif
+      }
+
+     // inserisco tutti i triangoli validi
+      if ( ! pStm->AddTriaFromZMap( vTriaSafe, VertGrid))
+         return nullptr ;
+   }
+
+   #if DEBUG
+      SaveGeoObj( VT, VC, "C:\\Temp\\TriangleSelection.nge") ;
+   #endif
+
+
+  // sistemo la topologia
+   if ( ! pStm->AdjustTopologyFromZMap())
+      return nullptr ;
+
+   return ( Release( pStm)) ;
 }
 
 //----------------------------------------------------------------------------
@@ -107,10 +231,12 @@ CreateSurfTriMeshesThickeningOffset( const CISURFTMPVECTOR& vStm, double dOffs,
       return SumStm( vStm) ;
 
   // creo lo Zmap associato alle superfici TriMesh
-   PtrOwner<IVolZmap> pVolZmap( CreateVolZmap()) ;
-   if ( IsNull( pVolZmap)  ||  ! pVolZmap->CreateFromTriMeshThickeningOffset( vStm, dOffs, dMyPrec, nType))
+   VolZmap OneVolZmap ;
+   if ( ! OneVolZmap.CreateFromTriMeshThickeningOffset( vStm, dOffs, dMyPrec, nType))
+      return nullptr ;
+   if ( ! OneVolZmap.IsValid())
       return nullptr ;
 
   // restituisco la superficie TriMesh
-   return ( pVolZmap->GetSurfTriMesh()) ;
+   return ( OneVolZmap.GetSurfTriMesh()) ;
 }

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, const 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,182 @@ 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 bParamDimOk = false ;
+            if( bVert)
+               bParamDimOk = ( pcToSplit->GetTopRight().x - pcToSplit->GetBottomLeft().x) / 2 > 100 * EPS_PARAM ;
+            else
+               bParamDimOk = ( pcToSplit->GetTopRight().y - pcToSplit->GetBottomLeft().y) / 2 > 100 * EPS_PARAM ;
+            bool bDimOk = ( dSideMinVal / 2 >= dSideMin  ||  ( dSideMinVal < EPS_SMALL  &&  dLengMinVal / 2 >= dSideMin))  &&  bParamDimOk ;
+            if ( dSideMaxVal > dSideMax) {
+               bSplit = true ;
+               //LOG_DBG_INFO( GetEGkLogger(), "  Split by SideMax")
             }
-            else if ( 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 +781,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 +813,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 +844,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 +888,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 +1372,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 +1768,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 +1811,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 +1821,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 +4164,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

@@ -20,15 +20,17 @@
 #include "/EgtDev/Include/EGkPolyLine.h"
 #include "/EgtDev/Include/EGkChainCurves.h"
 #include <map>
+#include "/EgtDev/Extern/abseil/Include/absl/algorithm/container.h"
+#include "/EgtDev/Extern/abseil/Include/absl/container/flat_hash_map.h"
 #include <utility>
 
 struct PairHashInt64 {
-   size_t operator()(const pair<int64_t, int64_t>& key) const {
+   size_t operator()(const std::pair<int64_t, int64_t>& key) const {
       size_t h1 = std::hash<int64_t>{}(key.first) ;
       size_t h2 = std::hash<int64_t>{}(key.second) ;
       return h1 ^ (h2 << 1);  // Combine hashes
    }
-};
+} ;
 
 //----------------------------------------------------------------------------
 struct Inters {
@@ -42,13 +44,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 +65,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 +77,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 +94,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 +127,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 +165,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 +192,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,41 +247,43 @@ 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) ;
-      bool GetPolygons( POLYLINEMATRIX& vvPolygons, POLYLINEMATRIX& vvPolygons3d, vector<ICRVCOMPOPOVECTOR>& vCCEdges3D, ICRVCOMPOPOVECTOR& vCCLoops) ;
+      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
                                                                                                                            // ad alcuni poligoni potrebbero venire tolti dei punti per evitare errori dovuti ad eventuali poli sui bordi del parametrico
-      bool GetLeaves ( std::vector<Cell>& vLeaves) const ;  // restituisce gli indici delle foglie nell'albero
-      bool GetEdges3D ( 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
+      bool GetLeaves( std::vector<Cell>& vLeaves) const ;         // restituisce gli indici delle foglie nell'albero
+      bool GetEdges3D( std::vector<ICRVCOMPOPOVECTOR>& mCCEdge, POLYLINEVECTOR& vPolygons) ;  // restituisce gli edge 3D come polyline
+      bool GetSplitLoops( ICRVCOMPOPOVECTOR& vCCLoopSplit) const  // 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 ; }
       // 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
@@ -304,31 +308,27 @@ class Tree
       bool UpdateSplitLoop( ICurveComposite* pCC, Point3d& pt) ;
       bool VerifyLoopOrientation( ICURVEPLIST& vpCrv, BOOLVECTOR& vbOrientation) const ;  // verifico l'orientazione ( CCW o CW) delle polyline in base a come sono contenute le une nelle altre
       bool AdjustLoop( PolyLine& pl, POLYLINEVECTOR& vPl, BOOLVECTOR& vbOrientation) const ;
-      bool GetPoint(double dU, double dV, Point3d& pt) const ;
-      bool SavePoint( double dU, double dV, Point3d& pt) ;
-
+      bool GetPoint(double dU, double dV, Point3d& ptP) const ;
+      bool SavePoint( double dU, double dV, const Point3d& ptP) ;
 
    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
-      unordered_map<int,int>      m_mChunk ;           // mappa in cui vengono salvati chunk di appartenza per ogni loop di trim
-      vector<tuple<PolyLine,bool>> m_vPlApprox ;       // vettore contenente le approssimazioni dei loop  // il bool indica se la curva è CCW
+      absl::flat_hash_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
       bool                        m_bBilinear ;        // superficie bilineare
       bool                        m_bMulti ;           // superficie multi-patch
       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
       int                         m_nSpanV ;           // numero di span lungo il parametro V
-      unordered_map<int,Cell>     m_mTree ;            // mappa che contiene tutti i nodi e le foglie dell'albero. -2 è puntatore Null e -1 è root
-      mutable unordered_map<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)
+      absl::flat_hash_map<int,Cell>     m_mTree ;       // mappa che contiene tutti i nodi e le foglie dell'albero. -2 è puntatore Null e -1 è root
+      mutable absl::flat_hash_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
-      vector<pair<BIPNTVECTOR, ChainCurves>>    m_vCEdge2D ;        // vettore che le chain che rappresentano ciò che resta degli edge originali, tenendo conto dei trim.
+      std::vector<std::pair<BIPNTVECTOR, ChainCurves>>    m_vCEdge2D ;        // vettore che le chain che rappresentano ciò che resta degli edge originali, tenendo conto dei trim.
 } ;

VolZmap.cpp

@@ -1929,8 +1929,6 @@ VolZmap::SetToModifyDexelBlocks( int nGrid, int nDex, int nInt)
 bool
 VolZmap::IsMapPartABox( int nMap, int nInfI, int nSupI, int nInfJ, int nSupJ, double& dMinZ, double& dMaxZ)
 {
-   if ( ! m_bIsBox)
-      return true ;
    dMinZ = m_dMaxZ[nMap] ;
    dMaxZ = m_dMinZ[nMap] ;
    for ( int i = nInfI ; i < nSupI ; ++ i) {
@@ -1960,6 +1958,7 @@ bool
 VolZmap::IsBox( void)
 {
   // Se non tridexel, non posso stabilire con il metodo seguente se è un box
+  // Verifico solo che gli spilloni di una mappa o sono nulli o hanno gli stessi estremi
    if ( m_nMapNum == 1)
       return false ;
   // Numero massimo di thread per il calcolo parallelo.
@@ -1968,9 +1967,12 @@ VolZmap::IsBox( void)
    if ( nThreadMax == 1) {
       for ( int nMap = 0 ; nMap < m_nMapNum ; ++ nMap) {
          double dMinZ, dMaxZ ;
-         if ( ! IsMapPartABox( nMap, 0, m_nNx[nMap], 0, m_nNy[nMap], dMinZ, dMaxZ))
+         if ( ! IsMapPartABox( nMap, 0, m_nNx[nMap], 0, m_nNy[nMap], dMinZ, dMaxZ)) {
+            m_bIsBox = false ;
             return false ;
+         }
       }
+      m_bIsBox = true ;
       return true ;
    }
 
@@ -2024,10 +2026,11 @@ VolZmap::IsBox( void)
          return false ;
       // Controllo che gli estremi Z siano uguali.
       for ( int nT = 1 ; nT < nThreadMax ; ++ nT) {
-         if ( abs( vMinZ[nT] - vMinZ[0]) > EPS_SMALL)
-            return false ;
-         if ( abs( vMaxZ[nT] - vMaxZ[0]) > EPS_SMALL)
+         if ( abs( vMinZ[nT] - vMinZ[0]) > EPS_SMALL  ||
+              abs( vMaxZ[nT] - vMaxZ[0]) > EPS_SMALL) {
+            m_bIsBox = false ;
             return false ;
+         }
       }
    }
 

VolZmap.h

@@ -85,9 +85,7 @@ class VolZmap : public IVolZmap, public IGeoObjRW
       bool Create( const Point3d& ptO, double dDimX, double dDimY, double dDimZ, double dStep, bool bTriDex) override ;
       bool CreateEmpty( const Point3d& ptO, double dDimX, double dDimY, double dDimZ, double dStep, bool bTriDex) override ;
       bool CreateFromFlatRegion( const ISurfFlatRegion& Surf, double dDimZ, double dStep, bool bTriDex) override ;
-      bool CreateFromTriMesh( const ISurfTriMesh& Surf, double dStep, bool bTriDex, double dExtraBox = 0.) override ;
-      bool CreateFromTriMeshOffset( const CISURFTMPVECTOR& vSurf, double dOffs, double dTol, int nType = STMOFF_FILLET) override ;
-      bool CreateFromTriMeshThickeningOffset( const CISURFTMPVECTOR& vSurf, double dOffs, double dTol, int nType = STMOFF_FILLET) override ;
+      bool CreateFromTriMesh( const ISurfTriMesh& Surf, double dStep, bool bTriDex, double dExtraBox = 0) override ;
       int  GetBlockCount( void) const override ;
       int  GetBlockUpdatingCounter( int nBlock) const override ;
       bool GetBlockTriangles( int nBlock, TRIA3DEXVECTOR& vTria) const override ;
@@ -101,7 +99,7 @@ class VolZmap : public IVolZmap, public IGeoObjRW
               { return m_nDexVoxRatio ; }   
       bool ChangeResolution( int nDexVoxRatio) override ;
       void SetShowEdges( bool bShow) override
-              { m_bShowEdges = bShow ;  // qui � necessario far ricreare la grafica
+              { m_bShowEdges = bShow ;  // qui è necessario far ricreare la grafica
                 m_OGrMgr.Clear() ; }
       bool GetShowEdges( void) const override
               { return m_bShowEdges ; }
@@ -154,6 +152,7 @@ class VolZmap : public IVolZmap, public IGeoObjRW
       bool AddSurfTm( const ISurfTriMesh* pStm) override ;
       bool SubtractSurfTm( const ISurfTriMesh* pStm) override ;
       bool MakeUniform( double dToler) override ;
+      bool Offset( double dOffs, int nType) override ;
 
    public :                      // IGeoObjRW
       int GetNgeId( void) const override ;
@@ -174,6 +173,8 @@ class VolZmap : public IVolZmap, public IGeoObjRW
                     return *this ; }
       bool GetLineIntersection( const Point3d& ptP, const Vector3d& vtD, ILZIVECTOR& vIntersInfo) const ;
       bool GetPlaneIntersection( const Plane3d& plPlane, ICURVEPOVECTOR& vpLoop) const ;
+      bool CreateFromTriMeshOffset( const CISURFTMPVECTOR& vSurf, double dOffs, double dTol, int nType = STMOFF_FILLET) ;
+      bool CreateFromTriMeshThickeningOffset( const CISURFTMPVECTOR& vSurf, double dOffs, double dTol, int nType = STMOFF_FILLET) ;
 
    private :
       enum CubeType { VOX_EXTERN = 1,
@@ -267,7 +268,7 @@ class VolZmap : public IVolZmap, public IGeoObjRW
       bool MillingTranslationStep( const Point3d& ptPs, const Point3d& ptPe, const Vector3d& vtD, const Vector3d& vtA) ;
       bool MillingGeneralMotionStep( const Point3d& ptPs, const Vector3d& vtDs, const Vector3d& vtAs,
                                      const Point3d& ptPe, const Vector3d& vtDe, const Vector3d& vtAe) ;
-      bool SelectGeneralMotion( int nGrid, const PNTVECTOR& ptPs, const PNTVECTOR& ptPe, const VCT3DVECTOR& vtLs, const VCT3DVECTOR& vtLe, int n5AxisType) ;
+      bool SelectGeneralMotion( int nGrid, const Point3d& ptPs, const Point3d& ptPe, const Vector3d& vtLs, const Vector3d& vtLe, int n5AxisType) ;
       bool SelectMotion( int nGrid, const Point3d& ptLs, const Point3d& ptLe, const Vector3d& vtL, const Vector3d& vtAL) ;
       bool InitializePointsAndVectors( const Point3d& ptPs, const Point3d& ptPe, const Vector3d& vtDs, const Vector3d& vtAs,
                                        Point3d ptLs[3], Point3d ptLe[3], Vector3d vtLs[3], Vector3d vtALs[3]) ;
@@ -306,10 +307,10 @@ class VolZmap : public IVolZmap, public IGeoObjRW
       bool GenTool_Drilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir) ;
       bool GenTool_Milling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir) ;
      // lavorazioni a 5 assi
-      bool GenTool_5AxisMilling( int nGrid, const PNTVECTOR& ptS, const PNTVECTOR& ptE, const VCT3DVECTOR& vtLs, const VCT3DVECTOR& vtLe, int nToolNum, int n5AxisType = VolZmap::Move5Axis::ACROSS) ;
-      bool Cyl_5AxisMilling( int nGrid, const PNTVECTOR& ptS, const PNTVECTOR& ptE, const VCT3DVECTOR& vtLs, const VCT3DVECTOR& vtLe, int nToolNum, double dHeightCorr = 0, int n5AxisType = VolZmap::Move5Axis::ACROSS) ;
-      bool CylBall_5AxisMilling( int nGrid, const PNTVECTOR& ptS, const PNTVECTOR& ptE, const VCT3DVECTOR& vtLs, const VCT3DVECTOR& vtLe, int nToolNum, int n5AxisType = VolZmap::Move5Axis::ACROSS) ;
-      bool Conus_5AxisMilling( int nGrid, const PNTVECTOR& ptS, const PNTVECTOR& ptE, const VCT3DVECTOR& vtLs, const VCT3DVECTOR& vtLe, int nToolNum, int n5AxisType = VolZmap::Move5Axis::ACROSS) ;
+      bool GenTool_5AxisMilling(int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtLs, const Vector3d& vtLe, int nToolNum, int n5AxisType = ACROSS) ;
+      bool Cyl_5AxisMilling(int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtLs, const Vector3d& vtLe, int nToolNum, double dHeightCorr = 0, int n5AxisType = ACROSS) ;
+      bool CylBall_5AxisMilling(int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtLs, const Vector3d& vtLe, int nToolNum, int n5AxisType = ACROSS) ;
+      bool Conus_5AxisMilling(int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtLs, const Vector3d& vtLe, int nToolNum, int n5AxisType = ACROSS) ;
 
     // COMPONENTI
      // Asse di simmetria diretto come l'asse Z 
@@ -345,13 +346,13 @@ class VolZmap : public IVolZmap, public IGeoObjRW
       bool CompPar_Milling( int nGrid, double dLenX, double dLenY, double dLenZ,
                             const Point3d& ptS, const Point3d& ptE,
                             const Vector3d& vtToolDir, const Vector3d& vtAux, int nToolNum) ;   //  E' in realtà  MillingPerp
-     // lavorazioni a 5 assi
-      bool Comp_5AxisMilling( int nGrid, const PNTVECTOR& ptS, const PNTVECTOR& ptE, const VCT3DVECTOR& vtLs, const VCT3DVECTOR& vtLe,
+      // lavorazioni a 5 assi
+      bool Comp_5AxisMilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtLs, const Vector3d& vtLe,
                               double dHeight, double dMaxRad, double dMinRad, int nToolNum, int n5AxisType) ;
-      bool CompCyl_5AxisMilling( int nGrid, const PNTVECTOR& ptS, const PNTVECTOR& ptE, const VCT3DVECTOR& vtLs, const VCT3DVECTOR& vtLe,
-                                 double dHeight, double dRadius, int nToolNum,int n5AxisType) ;
-      bool CompConus_5AxisMilling( int nGrid, const PNTVECTOR& ptS, const PNTVECTOR& ptE, const VCT3DVECTOR& vtToolDirS, const VCT3DVECTOR& vtToolDirE, double dHei, double dMaxRad, double dMinRad,
-                                   bool bTapB, bool bTapT,const Vector3d& vtArcNormMaxR, const Vector3d& vtArcNormMinR, int nToolNum, int n5AxisType) ;
+      bool CompCyl_5AxisMilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtLs, const Vector3d& vtLe,
+                                 double dHeight, double dRadius, int nToolNum, int n5AxisType) ;
+      bool CompConus_5AxisMilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDirS, const Vector3d& vtToolDirE, double dHei, double dMaxRad, double dMinRad,
+                                   bool bTapB, bool bTapT, const Vector3d& vtArcNormMaxR, const Vector3d& vtArcNormMinR, int nToolNum, int n5AxisType) ;
 
      // Generica traslazione sfera
       bool CompBall_Milling( int nGrid, const Point3d& ptS, const Point3d& ptE, double dRad, int nToolNum) ;
@@ -398,8 +399,6 @@ class VolZmap : public IVolZmap, public IGeoObjRW
                                        const Frame3d& frTruncPyramFrame, double dSegMin, double dSegMax, double dHeight,
                                        Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2) const ;
       bool TestIntersPlaneZmapBBox( const Plane3d& plPlane) const ;
-     // aggiornamento intersezioni
-      bool UpdateMaxMin( Point3d& ptBez, Vector3d& vtN, PNTVEC3DVECTOR& vInters) const ;
      // Voxel: esistenza e passaggio da N a ijk per i voxel
       bool IsValidVoxel( int nN) const ;
       bool IsValidVoxel( int nI, int nJ, int nK) const ;
@@ -462,52 +461,51 @@ class VolZmap : public IVolZmap, public IGeoObjRW
       bool SubtractMapPart( int nMap, int nInfI, int nSupI, int nInfJ, int nSupJ, const Vector3d& vtLen, const Point3d& ptMapOrig,
                             const ISurfTriMesh& Surf, IntersParLinesSurfTm& intPLSTM) ;
      // Funzioni per Offset di superfici
-      bool InitVolZMapOffs( const CISURFTMPVECTOR& vSurf, double dOffs, double dTol) ;
-      bool InitVolZMapThickeningOffs( const CISURFTMPVECTOR& vSurf, double dOffs, double dTol) ;
+      bool InitVolZMapOffset( const CISURFTMPVECTOR& vSurf, double dOffs, double dTol) ;
+      bool InitVolZMapThickeningOffset( const CISURFTMPVECTOR& vSurf, double dOffs, double dTol) ;
       bool UpdateVolZMapByOpenSurfFilletOffset( const ISurfTriMesh* Surf, double dOffs, double dTol) ;
+      bool UpdateVolZMapByOpenSurfSharpedOffset( const ISurfTriMesh* Surf, int nType, double dOffs, double dTol) ;
       bool UpdateVolZMapByClosedSurfFilletOffset( const ISurfTriMesh* Surf, double dOffs, double dTol) ;
+      bool UpdateVolZMapByClosedSurfSharpedOffset( const ISurfTriMesh* Surf, int nType, double dOffs, double dTol) ;
       bool UpdateVolZMapBySurfThickeningFilletOffset( const ISurfTriMesh* Surf, double dOffs, double dTol) ;
-      bool CreateOffsSphereOnVertex( const Point3d& ptV, double dOffs, int nGrid, int nVertexType = 0) ;
-      bool CreateOffsCylinderOnEdge( const Point3d& ptP1, const Point3d& ptP2, double dOffs, int nGrid, int nVertexType = 0) ;
+      bool UpdateVolZMapBySurfThickeningSharpedOffset( const ISurfTriMesh* Surf, int nType, double dOffs, double dTol) ;
+      bool CreateOffsetSphereOnVertex( const Point3d& ptV, double dOffs, int nGrid, int nVertexType = 0) ;
+      bool CreateOffsetCylinderOnEdge( const Point3d& ptP1, const Point3d& ptP2, double dOffs, int nGrid, int nVertexType = 0) ;
+      bool CreateFatOffsetExtrusionFace( const ISurfTriMesh* Surf, double dOffs, bool bThickle) ;
+      bool CreateOrientedOffsetExtrusionFace( const ISurfTriMesh* Surf, double dOffs) ;
       bool SubtractIntervalsForOffset( int nGrid, int nI, int nJ,
                                        double dMin, double dMax, const Vector3d& vtNMin, const Vector3d& vtNMax, 
                                        int nToolNum, bool bSkipSwap = false) ;
       bool AddIntervalsForOffset( int nGrid, int nI, int nJ,
                                   double dMin, double dMax, const Vector3d& vtNMin, const Vector3d& vtNMax,
                                   int nToolNum, bool bSkipSwap = false) ;
-   public :
-      // ------------------------- ENUM ----------------------------------------------------------------
-      enum MillingPhase {
-         COUNT_START_VOL = 0 ,
-         ONLY_LATERAL_SURF = 1 ,
-         COUNT_END_VOL = 2 ,
-         COUNT_START_END = 3
-      } ;
-
-      enum Move5Axis{
-         ALONG_CONVEX = 0 ,
-         ALONG_CONCAVE = 1 ,
-         ACROSS = 2 ,
-         NO_BASE_INTERS = 3
-      } ;
-
-      enum CuttingSurface {
-         NONE = -1,
-         TOOL = 0 ,
-         BZ = 1 ,
-         LATERAL = 2
-      } ;
+      bool CutByPlaneForOffset( const Plane3d& plCut) ;
+     // Funzioni per Offset di Zmap
+      bool OffsetFillet( double dOffs) ;
+      bool OffsetSharped( double dOffs, int nType) ;
 
    private :
-      enum Status { ERR = 0, OK = 1, TO_VERIFY = 2} ;
-      enum Shape { GENERIC = 0, BOX = 1, EXTRUSION = 2, OFFSET = 3} ;
+      enum Move5Axis {
+            ALONG_CONVEX = 0 ,
+            ALONG_CONCAVE = 1 ,
+            ACROSS = 2 ,
+            NO_BASE_INTERS = 3} ;
+      enum Status {
+            ERR = 0,
+            OK = 1,
+            TO_VERIFY = 2} ;
+      enum Shape {
+            GENERIC = 0,
+            BOX = 1,
+            EXTRUSION = 2,
+            OFFSET = 3} ;
       static const int N_MAPS = 3 ;
       static const int N_VOXBLOCK = 32 ;
 
    private :
       ObjGraphicsMgr m_OGrMgr ;             // gestore grafica dell'oggetto
       Status         m_nStatus ;            // stato
-      int            m_nTempProp[2] ;       // vettore propriet� temporanee
+      int            m_nTempProp[2] ;       // vettore proprietà temporanee
       double         m_dTempParam[2] ;      // vettore parametri temporanei
       bool           m_bShowEdges ;         // flag di visualizzazione spigoli vivi
       Frame3d        m_MapFrame ;           // riferimento intrinseco dello Zmap
@@ -541,8 +539,8 @@ class VolZmap : public IVolZmap, public IGeoObjRW
       mutable BOOLVECTOR m_BlockToUpdate ;
       mutable INTVECTOR  m_BlockUpdatingCounter ;
 
-      int m_nConnectedCompoCount ;          // Se == - 1 il numero di componenti non � noto
-                                            // Se >= 0 � il numero di componenti connesse
+      int m_nConnectedCompoCount ;          // Se == - 1 il numero di componenti non è noto
+                                            // Se >= 0 è il numero di componenti connesse
               
       mutable std::vector<VoxelContainer>  m_InterBlockVox ;
       mutable SharpTriaMatrix              m_InterBlockOriginalSharpTria ;
@@ -565,6 +563,12 @@ class VolZmap : public IVolZmap, public IGeoObjRW
       double                              m_dToolAngTolDeg ;
 } ;
 
+// Offset
+enum {
+   VOLZMAP_OFFS_FILLET = 0,
+   VOLZMAP_OFFS_CHANFER = 1,
+   VOLZMAP_OFFS_EXTENDED = 2
+} ;
 
 //-----------------------------------------------------------------------------
 inline VolZmap* CreateBasicVolZmap( void)

VolZmapCreation.cpp

@@ -505,6 +505,8 @@ VolZmap::CreateMapPart( int nMap, int nInfI, int nSupI, int nInfJ, int nSupJ, co
         nSupJ < 0  ||  nSupJ > m_nNy[nMap])
       return false ;
 
+   double dCosSmall = sin( EPS_ANG_SMALL * DEGTORAD) ;
+
    // Determinazione e ridimensionamento dei dexel interni alla trimesh
    for ( int i = nInfI ; i < nSupI ; ++ i) {
       for ( int j = nInfJ ; j < nSupJ ; ++ j) {
@@ -552,7 +554,7 @@ VolZmap::CreateMapPart( int nMap, int nInfI, int nSupI, int nInfJ, int nSupJ, co
                double dCos = IntersectionResults[k].dCosDN ;
 
                // entro nella superficie trimesh
-               if ( dCos < - EPS_SMALL) {
+               if ( dCos < - dCosSmall) {
 
                   ptIn = IntersectionResults[k].ptI ;
 
@@ -565,7 +567,7 @@ VolZmap::CreateMapPart( int nMap, int nInfI, int nSupI, int nInfJ, int nSupJ, co
                }
 
                // esco dalla superficie trimesh
-               else if ( dCos > EPS_SMALL  &&  bInside) {
+               else if ( dCos > dCosSmall  &&  bInside) {
 
                   Point3d ptOut = IntersectionResults[k].ptI ;
 
@@ -804,7 +806,9 @@ VolZmap::CreateFromTriMesh( const ISurfTriMesh& Surf, double dStep, bool bTriDex
   // quindi espandiamo il bounding box per ovviare al problema.
    if ( dExtraBox > EPS_ZERO)
       SurfBBox.Expand( dExtraBox) ;
-
+   else
+      dExtraBox = 0 ;
+    
   // Determino i punti estremi del bounding box
    Point3d ptMapOrig, ptMapEnd ; 
    SurfBBox.GetMinMax( ptMapOrig, ptMapEnd) ;
@@ -870,58 +874,69 @@ VolZmap::CreateFromTriMesh( const ISurfTriMesh& Surf, double dStep, bool bTriDex
      // Oggetto per calcolo massivo intersezioni
       IntersParLinesSurfTm intPLSTM( frMapFrame, Surf) ;
 
-      // Numero massimo di thread
-      int nThreadMax = max( 1, int( thread::hardware_concurrency()) - 1) ;
-      vector< future<bool>> vRes ;
-      vRes.resize( nThreadMax) ;
-      if ( m_nNx[nG] > m_nNy[nG]) {
-         int nDexNum = m_nNx[nG] / nThreadMax ;
-         int nRemainder = m_nNx[nG] % nThreadMax ;
-         int nInfI = 0 ;
-         int nSupI = 0 ;
-         for ( int nThread = 0 ; nThread < nThreadMax ; ++ nThread) {
-            nInfI = nSupI ;
-            nSupI = nInfI + ( nThread < nRemainder ? nDexNum + 1 : nDexNum) ;
-            vRes[nThread] = async( launch::async, &VolZmap::CreateMapPart, this, nG,
-                                   nInfI, nSupI, 0, m_nNy[nG], ref( vtLen), ref( ptMapOrig), ref( Surf), ref( intPLSTM)) ;
-         }
-      }
-      else {
-         int nDexNum = m_nNy[nG] / nThreadMax ;
-         int nRemainder = m_nNy[nG] % nThreadMax ;
-         int nInfJ = 0 ;
-         int nSupJ = 0 ;
-         for ( int nThread = 0 ; nThread < nThreadMax ; ++ nThread) {
-            nInfJ = nSupJ ;
-            nSupJ = nInfJ + ( nThread < nRemainder ? nDexNum + 1 : nDexNum) ;
-            vRes[nThread] = async( launch::async, &VolZmap::CreateMapPart, this, nG,
-                                   0, m_nNx[nG], nInfJ, nSupJ, ref( vtLen), ref( ptMapOrig), ref( Surf),ref( intPLSTM)) ;
-         }
-      }
-       
-      // Ciclo per attendere che tutti gli async abbiano terminato.
-      int nTerminated = 0 ;
-      while ( nTerminated < nThreadMax) {
-         for ( int nL = 0 ; nL < nThreadMax ; ++ nL) {
-            // Async terminato
-            if ( vRes[nL].valid()  &&  vRes[nL].wait_for( chrono::microseconds{ 1}) == future_status::ready) {
-               ++ nTerminated ;
-               bCompleted = bCompleted  &&  vRes[nL].get() ;
+     // Standarda è multithread
+      constexpr bool MULTITHREAD = true ;
+      if ( MULTITHREAD) {
+
+         // Numero massimo di thread
+         int nThreadMax = max( 1, int( thread::hardware_concurrency()) - 1) ;
+         vector< future<bool>> vRes ;
+         vRes.resize( nThreadMax) ;
+         if ( m_nNx[nG] > m_nNy[nG]) {
+            int nDexNum = m_nNx[nG] / nThreadMax ;
+            int nRemainder = m_nNx[nG] % nThreadMax ;
+            int nInfI = 0 ;
+            int nSupI = 0 ;
+            for ( int nThread = 0 ; nThread < nThreadMax ; ++ nThread) {
+               nInfI = nSupI ;
+               nSupI = nInfI + ( nThread < nRemainder ? nDexNum + 1 : nDexNum) ;
+               vRes[nThread] = async( launch::async, &VolZmap::CreateMapPart, this, nG,
+                                      nInfI, nSupI, 0, m_nNy[nG], ref( vtLen), ref( ptMapOrig), ref( Surf), ref( intPLSTM)) ;
             }
          }
+         else {
+            int nDexNum = m_nNy[nG] / nThreadMax ;
+            int nRemainder = m_nNy[nG] % nThreadMax ;
+            int nInfJ = 0 ;
+            int nSupJ = 0 ;
+            for ( int nThread = 0 ; nThread < nThreadMax ; ++ nThread) {
+               nInfJ = nSupJ ;
+               nSupJ = nInfJ + ( nThread < nRemainder ? nDexNum + 1 : nDexNum) ;
+               vRes[nThread] = async( launch::async, &VolZmap::CreateMapPart, this, nG,
+                                      0, m_nNx[nG], nInfJ, nSupJ, ref( vtLen), ref( ptMapOrig), ref( Surf),ref( intPLSTM)) ;
+            }
+         }
+       
+         // Ciclo per attendere che tutti gli async abbiano terminato.
+         int nTerminated = 0 ;
+         while ( nTerminated < nThreadMax) {
+            for ( int nL = 0 ; nL < nThreadMax ; ++ nL) {
+               // Async terminato
+               if ( vRes[nL].valid()  &&  vRes[nL].wait_for( chrono::microseconds{ 1}) == future_status::ready) {
+                  ++ nTerminated ;
+                  bCompleted = bCompleted  &&  vRes[nL].get() ;
+               }
+            }
+         }
+      }
+
+     // !!!! NON MULTITHREAD : SOLO PER DEBUG !!!!
+      else {
+         CreateMapPart( nG, 0, m_nNx[nG], 0, m_nNy[nG], vtLen, ptMapOrig, Surf, intPLSTM) ;
       }
    }
 
   // Assegno il minimo e massimo valore di Z della mappa
-   m_dMinZ[0] = 0 ;
-   m_dMaxZ[0] = vtLen.z ;
-   m_dMinZ[1] = 0 ;
-   m_dMaxZ[1] = ( bTriDex ? vtLen.x : 0) ;
-   m_dMinZ[2] = 0 ;  
-   m_dMaxZ[2] = ( bTriDex ? vtLen.y : 0) ;
+   m_dMinZ[0] = dExtraBox ;
+   m_dMaxZ[0] = vtLen.z - dExtraBox ;
+   m_dMinZ[1] = ( bTriDex ? dExtraBox : 0) ;
+   m_dMaxZ[1] = ( bTriDex ? vtLen.x - dExtraBox : 0) ;
+   m_dMinZ[2] = ( bTriDex ? dExtraBox : 0) ;  
+   m_dMaxZ[2] = ( bTriDex ? vtLen.y - dExtraBox : 0) ;
 
   // Tipologia
-   m_nShape = ( dExtraBox > EPS_ZERO && IsBox() ? BOX : GENERIC) ;
+  // Con espansione non va considerato box (calcolo trimesh va in crash)
+   m_nShape = ( dExtraBox <= EPS_ZERO && IsBox() ? BOX : GENERIC) ;
 
   // Aggiornamento dello stato
    m_nStatus = OK ;

VolZmapGraphics.cpp

@@ -950,25 +950,41 @@ VolZmap::UpdateTripleMapGraphics( void) const
       }
    }
 
-  // Calcolo i triangoli sui blocchi
-   int nBlockUpdated = 0 ;
-   vector< future<bool>> vRes ;
-   vRes.resize( m_nNumBlock) ;
-   for ( int i = 0 ; i < m_nNumBlock ; ++ i) {
-     // Se il blocco deve essere processato
-      if ( m_BlockToUpdate[i]) {
-        // processo ...
-         ++ nBlockUpdated ;
-         vRes[i] = async( launch::async, &VolZmap::ExtMarchingCubes, this, i, ref( vVoxContainerVec[i])) ;
+  // Standarda è multithread
+   constexpr bool MULTITHREAD = true ;
+   if ( MULTITHREAD) {
+
+     // Calcolo i triangoli sui blocchi
+      int nBlockUpdated = 0 ;
+      vector< future<bool>> vRes ;
+      vRes.resize( m_nNumBlock) ;
+      for ( int i = 0 ; i < m_nNumBlock ; ++ i) {
+        // Se il blocco deve essere processato
+         if ( m_BlockToUpdate[i]) {
+           // processo ...
+            ++ nBlockUpdated ;
+            vRes[i] = async( launch::async, &VolZmap::ExtMarchingCubes, this, i, ref( vVoxContainerVec[i])) ;
+         }
+      }
+      bool bOk = true ;
+      int nTerminated = 0 ;
+      while ( nTerminated < nBlockUpdated) {
+         for ( int i = 0 ; i < m_nNumBlock ; ++ i) {
+            if ( m_BlockToUpdate[i]  &&  vRes[i].valid()  &&  vRes[i].wait_for( chrono::nanoseconds{ 1}) == future_status::ready) {
+               bOk = vRes[i].get()  &&  bOk ;
+               ++ nTerminated ;
+            }
+         }
       }
    }
-   bool bOk = true ;
-   int nTerminated = 0 ;
-   while ( nTerminated < nBlockUpdated) {
+   else {
+     // Calcolo i triangoli sui blocchi
+      bool bOk = true ;
       for ( int i = 0 ; i < m_nNumBlock ; ++ i) {
-         if ( m_BlockToUpdate[i]  &&  vRes[i].valid()  &&  vRes[i].wait_for( chrono::nanoseconds{ 1}) == future_status::ready) {
-            bOk = vRes[i].get()  &&  bOk ;
-            ++ nTerminated ;
+        // Se il blocco deve essere processato
+         if ( m_BlockToUpdate[i]) {
+           // processo ...
+            bOk = ExtMarchingCubes( i, vVoxContainerVec[i]) && bOk ;
          }
       }
    }

stdafx.h

@@ -34,3 +34,4 @@
 #pragma comment(lib, EGTLIBDIR "EgtGeneral" EGTLIBVER ".lib")
 #pragma comment(lib, EGTLIBDIR "EgtNumKernel" EGTLIBVER ".lib")
 #pragma comment(lib, EGTLIBDIR "SEgtLock" EGTLIBVER ".lib")
+#pragma comment(lib, EGTLIBDIR "abseil" EGTLIBVER ".lib")