EgtGeomKernel :

- piccola correzione.

OffsetCurve.cpp

@@ -60,8 +60,8 @@ OffsetCurve::Reset( void)
       }
    }
    m_CrvLst.clear() ;
-   m_pt = P_INVALID ;
-   m_vt = V_INVALID ;
+   m_ptOffs = P_INVALID ;
+   m_vtOut = V_INVALID ;
    return true ;
 }
 
@@ -714,8 +714,8 @@ OffsetCurve::Make( const ICurve* pCrv, double dDist, int nType)
            // ( e.g. alcune curve aperte) non se ne restiusce nessuno. Da estendere, se necessario, individuando i punti dal lato 
            // di offset richiesto
             if ( vPntOffs.size() == 1) {
-               m_pt = vPntOffs[0].first ;
-               m_vt = vPntOffs[0].second ;            
+               m_ptOffs = vPntOffs[0].first ;
+               m_vtOut = vPntOffs[0].second ;            
             }        
          }
       }         
@@ -747,20 +747,20 @@ OffsetCurve::Make( const ICurve* pCrv, double dDist, int nType)
          pCrv->Scale( GLOB_FRM, 1 / dExtrOnN, 1, 1) ;
          pCrv->ToGlob( frCopy) ;
       }
-      if ( m_pt.IsValid()) {
-         m_pt.Scale( GLOB_FRM, 1 / dExtrOnN, 1, 1) ;
-         m_pt.ToGlob( frCopy) ;
-         m_vt.Scale( GLOB_FRM, 1 / dExtrOnN, 1, 1) ;
-         m_vt.ToGlob( frCopy) ;
+      if ( m_ptOffs.IsValid()) {
+         m_ptOffs.Scale( GLOB_FRM, 1 / dExtrOnN, 1, 1) ;
+         m_ptOffs.ToGlob( frCopy) ;
+         m_vtOut.Scale( GLOB_FRM, 1 / dExtrOnN, 1, 1) ;
+         m_vtOut.ToGlob( frCopy) ;
       }
 
    }
    else if ( bNeedRef) {
       for ( auto pCrv : m_CrvLst)
          pCrv->ToGlob( frCopy) ;
-      if ( m_pt.IsValid()) {       
-         m_pt.ToGlob( frCopy) ;
-         m_vt.ToGlob( frCopy) ;
+      if ( m_ptOffs.IsValid()) {       
+         m_ptOffs.ToGlob( frCopy) ;
+         m_vtOut.ToGlob( frCopy) ;
       }
    }
 
@@ -841,14 +841,14 @@ OffsetCurve::GetShorterCurve( void)
 
 //----------------------------------------------------------------------------
 bool
-OffsetCurve::GetPointOffset( Point3d& pt, Vector3d& vt)
+OffsetCurve::GetPointOffset( Point3d& ptOffs, Vector3d& vtOut)
 {
   // verifico se valori validi da restituire
-   if ( ! m_pt.IsValid()  ||  ! m_vt.IsValid())
+   if ( ! m_ptOffs.IsValid()  ||  ! m_vtOut.IsValid())
       return false ;
 
-   pt = m_pt ;
-   vt = m_vt ;
+   ptOffs = m_ptOffs ;
+   vtOut = m_vtOut ;
 
    return true ;
 }

SurfBezier.cpp

@@ -40,6 +40,7 @@
 #include "/EgtDev/Include/EGkGeoPoint3d.h"
 #include "/EgtDev/Include/EGkIntervals.h"
 #include "/EgtDev/Extern/Eigen/Dense"
+#include "/EgtDev/Include/EGkGeoObjSave.h"
 
 
 using namespace std ;
@@ -51,9 +52,11 @@ static bool ChangeStartForClosed( PolyLine& plU0, PolyLine& plU1, ICurveComposit
 static bool ParametrizeByLen( const ICurveComposite* pCurve, DBLVECTOR& vParam) ;
 static bool BuildCommonParam( const DBLMATRIX& mParam, DBLVECTOR& vCommonParam) ;
 
+static unordered_map<INTINT, DBLVECTOR, PairHashIntInt>  m_mBernCache ;  // mappa dei polinomi di bernstein
+
 //----------------------------------------------------------------------------
 SurfBezier::SurfBezier( void)
-   : m_pSTM( nullptr), m_nStatus( TO_VERIFY), m_nDegU(), m_nDegV(), m_nSpanU(), m_nSpanV(), m_bRat( false),
+   : 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}
 {
 }
@@ -92,6 +95,20 @@ SurfBezier::Init( int nDegU, int nDegV, int nSpanU, int nSpanV, bool bIsRational
    m_nStatus = TO_VERIFY ;
    m_nIsPlanar = NOT_CALCULATED ;
 
+  // setto la dimensione dei vettori
+   m_vBernU.resize( m_nDegU + 1) ;
+   m_vBernV.resize( m_nDegV + 1) ;
+   if (! m_bRat)
+      m_ptTemp.resize(m_nDegV + 1, ORIG) ;
+   else {
+      m_ptTempW.resize(m_nDegV + 1, ORIG) ;
+      m_dTempW.resize( m_nDegV + 1, 0) ;
+      m_vPtWCtrlLoc.resize( GetLocDim()) ;
+      m_vWeCtrlLoc.resize( GetLocDim()) ;
+   }
+   
+
+
   // imposto ricalcolo della grafica
    ResetAuxSurf() ;
    m_OGrMgr.Reset() ;
@@ -299,6 +316,123 @@ SurfBezier::GetCentroid( Point3d& ptCen) const
    return m_pSTM->GetCentroid( ptCen) ;
 }
 
+//----------------------------------------------------------------------------
+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()) {
+      DBLVECTOR vBern( nDegU + 1) ;
+      GetAllBernstein( dU, nDegU, vBern) ;
+      m_mBernCache[key] = vBern ;
+   }
+   vBernU = m_mBernCache[key] ;
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+bool
+SurfBezier::GetPoint( double dU, double dV, Side nUs, Side nVs,Point3d& ptPos) const
+{
+   // la curva deve essere validata
+   if ( m_nStatus != OK)
+      return false ;
+
+   // i parametri U e V devono essere compresi tra 0 e il corrispondente numero di Span
+   dU = Clamp( dU, 0., double( m_nSpanU)) ;
+   dV = Clamp( dV, 0., double( m_nSpanV)) ;
+
+   // determino gli intervalli di span e riduco i parametri in essi
+   int nBsU = min( int( dU), m_nSpanU - 1) ;
+   double dLocU = dU - nBsU ;
+   if ( abs( dLocU) < 5 * EPS_PARAM  &&  nBsU > 0  &&  nUs == ISurfBezier::FROM_MINUS) {
+      -- nBsU ;
+      dLocU = 1 ;
+   }
+   else if ( abs( dLocU) > 1 - 5 * EPS_PARAM  &&  nBsU < m_nSpanU - 1  &&  nUs == ISurfBezier::FROM_PLUS) {
+      ++ nBsU ;
+      dLocU = 0 ;
+   }
+   int nOffsU = nBsU * m_nDegU ;
+   int nBsV = min( int( dV), m_nSpanV - 1) ;
+   double dLocV = dV - nBsV ;
+   if ( abs( dLocV) < 5 * EPS_PARAM  &&  nBsV > 0  &&  nVs == ISurfBezier::FROM_MINUS) {
+      -- nBsV ;
+      dLocV = 1 ;
+   }
+   else if ( abs( dLocV) > 1 - 5 * EPS_PARAM  &&  nBsV < m_nSpanV - 1  &&  nVs == ISurfBezier::FROM_PLUS) {
+      ++ nBsV ;
+      dLocV = 0 ;
+   }
+   int nOffsV = nBsV * m_nDegV ;
+
+   // se forma polinomiale (o integrale)
+   if ( ! m_bRat) {
+
+      // calcolo dei polinomi di Bernstein per U di grado opportuno
+      GetBernstein( dLocU, m_nDegU, m_vBernU) ;
+
+      // calcolo dei punti intermedi
+      m_ptTemp.assign(GetLocDim(), ORIG) ;
+      for ( int j = 0 ; j <= m_nDegV ; ++ j) {
+         for ( int i = 0 ; i <= m_nDegU ; ++ i)
+            m_ptTemp[j] += m_vBernU[i] * m_vPtCtrl[GetInd( nOffsU + i, nOffsV + j)] ;
+      }
+
+      // calcolo dei polinomi di Bernstein per V di grado opportuno
+      GetBernstein( dLocV, m_nDegV, m_vBernV) ;
+
+      // calcolo del punto
+      ptPos = ORIG ;
+      for ( int j = 0 ; j <= m_nDegV ; ++ j)
+         ptPos += m_vBernV[j] * m_ptTemp[j] ;
+   }
+
+   // altrimenti forma razionale
+   else {
+      // recupero punti di controllo e pesi della patch in questione
+      // porto i punti in forma omogenea moltiplicandoli per i pesi
+      m_vPtWCtrlLoc.assign(GetLocDim(), ORIG) ;
+      m_vWeCtrlLoc.assign(GetLocDim(), 0) ;
+      for ( int j = 0 ; j <= m_nDegV ; ++ j) {
+         for ( int i = 0 ; i <= m_nDegU ; ++ i) {
+            m_vPtWCtrlLoc[GetLocInd( i, j)] = m_vWeCtrl[GetInd( nOffsU + i, nOffsV + j)] * m_vPtCtrl[GetInd( nOffsU + i, nOffsV + j)] ;
+            m_vWeCtrlLoc[GetLocInd( i, j)] = m_vWeCtrl[GetInd( nOffsU + i, nOffsV + j)] ;
+         }
+      }
+
+      // calcolo dei polinomi di Bernstein di grado opportuno
+      GetBernstein( dLocU, m_nDegU, m_vBernU) ;
+
+      // calcolo dei punti e pesi intermedi
+      m_ptTempW.assign(GetLocDim(), ORIG) ;
+      m_dTempW.assign(GetLocDim(), 0) ;
+      for ( int j = 0 ; j <= m_nDegV ; ++ j) {
+         for ( int i = 0 ; i <= m_nDegU ; ++ i) {
+            m_ptTempW[j] += m_vBernU[i] * m_vPtWCtrlLoc[GetLocInd( i, j)] ;
+            m_dTempW[j] += m_vBernU[i] * m_vWeCtrlLoc[GetLocInd( i, j)] ;
+         }
+      }
+
+      // calcolo dei polinomi di Bernstein per V di grado opportuno
+      GetBernstein( dLocV, m_nDegV, m_vBernV) ;
+
+      // calcolo del punto
+      double dW = 0 ;
+      ptPos = ORIG ;
+      for ( int j = 0 ; j <= m_nDegV ; ++ j) {
+         ptPos += m_vBernV[j] * m_ptTempW[j] ;
+         dW += m_vBernV[j] * m_dTempW[j] ;
+      }
+
+      // ritrasformo da forma omogenea a forma standard
+      double dInvW = 1 / ( ( dW > EPS_ZERO) ? dW : EPS_ZERO) ;
+      ptPos *= dInvW ;
+   }
+
+   return true ;
+}
+
 //----------------------------------------------------------------------------
 bool
 SurfBezier::GetPointD1D2( double dU, double dV, Side nUs, Side nVs,
@@ -1474,38 +1608,6 @@ SurfBezier::GetCurveOnVApproxLen( double dU) const
    return dLen ;
 }
 
-//----------------------------------------------------------------------------
-bool
-SurfBezier::UpdateEdgesFromTree( Tree& tr) const
-{
-   POLYLINEMATRIX mPlEdges ;
-   tr.GetEdges3D( mPlEdges) ;
-   for ( int i= 0 ; i < int( mPlEdges.size()); ++i) {
-      for ( int j = 0 ; j < int ( mPlEdges[i].size()) ; ++j) {
-         m_mCCEdge[i].emplace_back(CreateBasicCurveComposite()) ;
-         if ( ! m_mCCEdge[i].back()->FromPolyLine(mPlEdges[i][j]) ) {
-            Point3d ptStart ; 
-            if ( ! mPlEdges[i][j].GetFirstPoint( ptStart))
-               continue ;
-            m_mCCEdge[i].back()->FromPoint( ptStart) ;
-         }
-      }
-   }
-
-   if ( m_bTrimmed) {
-      POLYLINEVECTOR vPl ;
-      tr.GetSplitLoops( vPl) ;
-
-      // recupero i loop nel parametrico
-      for( int i = 0 ; i < int( vPl.size()); ++i) {
-         m_vCCLoop.emplace_back(CreateBasicCurveComposite()) ;
-         m_vCCLoop.back()->FromPolyLine(vPl[i]) ;
-      }
-   }
-
-   return true ;
-}
-
 //----------------------------------------------------------------------------
 const SurfTriMesh*
 SurfBezier::GetAuxSurf( void) const
@@ -1519,10 +1621,27 @@ SurfBezier::GetAuxSurf( void) const
    if ( m_pSTM != nullptr)
       return m_pSTM ;
   // eseguo calcolo
-   m_pSTM = GetApproxSurf( 50 * EPS_SMALL, 100 * EPS_SMALL) ;
+   m_pSTM = GetApproxSurf( 1000 * EPS_SMALL, 100 * EPS_SMALL) ;
    return m_pSTM ;
 }
 
+//----------------------------------------------------------------------------
+const SurfTriMesh*
+SurfBezier::GetAuxSurfRefined( void) const
+{
+  // la superficie deve essere validata
+   if ( m_nStatus != OK) {
+      ResetAuxSurf() ;
+      return nullptr ;
+   }
+  // se già calcolata, la restituisco
+   if ( m_pSTMRefined != nullptr)
+      return m_pSTMRefined ;
+  // eseguo calcolo
+   m_pSTMRefined = GetApproxSurf( 50 * EPS_SMALL, 100 * EPS_SMALL) ;
+   return m_pSTMRefined ;
+}
+
 //----------------------------------------------------------------------------
 SurfTriMesh*
 SurfBezier::GetApproxSurf( double dTol, double dSideMin) const
@@ -1531,10 +1650,14 @@ SurfBezier::GetApproxSurf( double dTol, double dSideMin) const
    if ( m_nStatus != OK)
       return nullptr ;
 
-  // se c'è ausiliaria e richiesta con gli stessi parametri, ne restituisco una copia
+  // se c'è ausiliaria di visualizzazione con gli stessi parametri, ne restituisco una copia
    if ( m_pSTM != nullptr  &&
-        abs( dTol - 50 * EPS_SMALL) < EPS_SMALL  &&  abs( dSideMin - 100 * EPS_SMALL) < EPS_SMALL)
+        abs( dTol - 1000 * EPS_SMALL) < 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)
+      return m_pSTMRefined->Clone() ;
   
   // costruttore della superficie
    POLYLINEMATRIX vvPL ;
@@ -1548,7 +1671,6 @@ SurfBezier::GetApproxSurf( double dTol, double dSideMin) const
    bool bTest = false ; // per debug
    // resetto il vettore degli edge
    m_mCCEdge.clear() ;
-   m_mCCEdge = vector<ICRVCOMPOPOVECTOR>(4) ;
    m_vCCLoop.clear() ;
    for ( int i = 0 ; i < (int) vTrees.size() ; ++ i) {
       Point3d ptMin = get<0>( vTrees[i]) ;
@@ -1564,15 +1686,15 @@ SurfBezier::GetApproxSurf( double dTol, double dSideMin) const
          Tree.BuildTree( dTol, dSideMin) ;
          //Tree.BuildTree( 1, 5) ;  //debug
       }
-      if ( ! Tree.GetPolygons( vvPL, vvPL3d))
+      if ( ! Tree.GetPolygons( vvPL, vvPL3d, m_mCCEdge, m_vCCLoop))
          continue ;
       //Tree.GetPolygonsBasic( vPL, true) ;  // per usare i polygon basic
 
      // aggiorno la chiusura della superficie
       m_bClosedU = m_bClosedU  ||  Tree.IsClosedU() ;
       m_bClosedV = m_bClosedV  ||  Tree.IsClosedV() ;
-     // salvo i bordi in 3d, che servono in caso si voglia trimmare la superficie DOPO aver costruito la trimesh ausiliaria
-      UpdateEdgesFromTree( Tree) ;
+     //// salvo i bordi in 3d, che servono in caso si voglia trimmare la superficie DOPO aver costruito la trimesh ausiliaria
+     // UpdateEdgesFromTree( Tree) ; // ora viene fatto nella GetPolygons
    }
 
   //// per usare i polygon basic//////////////////////
@@ -1859,6 +1981,7 @@ SurfBezier::GetLeaves( vector<tuple<int, Point3d, Point3d>>& vLeaves) const
          //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) ;
       }
       vector<Cell> vCells ;
@@ -1895,6 +2018,9 @@ SurfBezier::ResetAuxSurf( void) const
    if ( m_pSTM != nullptr)
       delete( m_pSTM) ;
    m_pSTM = nullptr ;
+   if ( m_pSTMRefined != nullptr)
+      delete( m_pSTMRefined) ;
+   m_pSTMRefined = nullptr ;
 
    m_mCCEdge.clear() ;
    m_vCCLoop.clear() ;
@@ -1937,6 +2063,8 @@ bool
 SurfBezier::IncreaseUV( double& dUV, double dxy, bool bUOrV, double* dUVCopy, bool bModifyOrig) const
 {
    double dUVTest ;
+   if (dUVCopy == nullptr  && !bModifyOrig)
+      return false ;
    if ( dUVCopy != nullptr) {
       *dUVCopy = dUV + dxy ;
       dUVTest = *dUVCopy ;
@@ -1993,6 +2121,7 @@ SurfBezier::UnprojectCurveFromStm( const ICurveComposite* pCC, ICRVCOMPOPVECTOR&
    // aggiungo tutti i successivi
    BIPNTVECTOR vBPnt ;
 
+   // per il primo punto non mi interessa sapere se è passato attraverso un edge, tanto ha già le coordinate parametriche giuste
    bThroughEdge = false ;
    int nRejected = 0 ;
    for ( int i = 0 ; i < int( pCC->GetCurveCount()) ; ++i) {
@@ -2158,7 +2287,7 @@ SurfBezier::CreateTrimRegionFromCuts( ICRVCOMPOPOVECTOR& vpCCOpen, ICRVCOMPOPOVE
       else
          pNewTrim.Set( GetSurfFlatRegionRectangle( SBZ_TREG_COEFF * m_nSpanU, SBZ_TREG_COEFF * m_nSpanV)) ;
 
-     // costruisco la mappa delle intersezioni, trovando tutte le intersezioni tra i trim e i loop dei vari chunk della falr region
+     // costruisco la mappa delle intersezioni, trovando tutte le intersezioni tra i trim e i loop dei vari chunk della flat region
       unordered_map<TRINT,ICCIVECTOR> mInters ;
       int nInters = 0 ;
       bool bStartFound = false ;
@@ -2298,28 +2427,34 @@ SurfBezier::CreateTrimRegionFromCuts( ICRVCOMPOPOVECTOR& vpCCOpen, ICRVCOMPOPOVE
          vbAdded[nNewToAdd] = true ;
         // di questo taglio mi salvo il chunk e loop di start e end
          TRINT tiStart, tiEnd ;
+         bool bStartFound = false ;
+         bool bEndFound = false ;
          for (const auto& pair : mInters) {
             if ( get<2>(pair.first) == nNewToAdd ) {
                for (int p = 0 ; p < int(pair.second.size()) ; ++p) {
                   if ( pair.second[p].IciB->dU < EPS_SMALL) {
                      tiStart = pair.first ;
+                     bStartFound = true ;
                      if ( bFirstCurveOfEdge){
                        // salvo l'inizio del taglio che è la prima curva di questa curva compo
                         tiFirstInters = pair.first ;
                         bFirstCurveOfEdge = false ;
                      }
                   }
-                  else
+                  else {
                      tiEnd = pair.first ;
+                     bEndFound = true ;
+                  }
                }
             }
          }
-
+         if ( ! bEndFound || ! bStartFound)
+            return nullptr ;
         // devo trovare fino a che punto seguire il loop che ho trovato come prosecuzione del taglio corrente
         // devo quindi trovare la prossima intersezione con un taglio
          int nInters = -1 ;
          double dNextCut = INFINITO ;
-         double dEndCurrentCut ;
+         double dEndCurrentCut = 0 ;
          for ( int i = 0 ; i < int( mInters[tiEnd].size()); ++i) {
            // se ho trovato l'intersezione con la fine del taglio corrente, salvo il parametro sul loop
             if ( mInters[tiEnd][i].IciB->dU > EPS_SMALL)
@@ -2422,7 +2557,7 @@ SurfBezier::Cut( const Plane3d& plPlane, bool bSaveOnEq)
    PNTVECTOR vPnt ;
    BIPNTVECTOR vBPnt ;
    TRIA3DVECTOR vTria ;
-   IntersPlaneSurfTm( plPlane, *GetAuxSurf(), vPnt, vBPnt, vTria) ;
+   IntersPlaneSurfTm( plPlane, *GetAuxSurfRefined(), vPnt, vBPnt, vTria) ;
 
    // concateno le curve 3d
    ChainCurves chainC ;
@@ -2493,8 +2628,7 @@ SurfBezier::Cut( const Plane3d& plPlane, bool bSaveOnEq)
       vpCCClosed.back()->FromPolyLine( vPLTria[i]) ;
    }
 
-   // ora posso chiamare la costruzione dello spazio parametrico trimmato
-   
+  // ora posso chiamare la costruzione dello spazio parametrico trimmato
    PtrOwner<ISurfFlatRegion> pSFR( CreateTrimRegionFromCuts( vpCCOpen, vpCCClosed)) ;
    if ( IsNull( pSFR)  ||  ! pSFR->IsValid())
       return false ;
@@ -2510,9 +2644,8 @@ SurfBezier::Cut( const Plane3d& plPlane, bool bSaveOnEq)
    double dAreaMin = 1e30 ;
    bool bPos = false ;
    for ( int c = 0 ; c < int( pSFR->GetChunkCount()); ++c) {
-      PtrOwner<ISurfFlatRegion> pSurf( pSFR->CloneChunk( c)) ;
-         for ( int l = 0 ; l < pSurf->GetLoopCount( 0); ++l) {
-            PtrOwner<ICurve> pCrv( pSurf->GetLoop( 0, l)) ;
+         for ( int l = 0 ; l < pSFR->GetLoopCount( c); ++l) {
+            PtrOwner<ICurve> pCrv( pSFR->GetLoop( c, l)) ;
             double dArea ; pCrv->GetAreaXY( dArea) ;
             if ( abs( dArea) < dAreaMin) {
                nChunkMin = c ;
@@ -2574,6 +2707,9 @@ SurfBezier::Cut( const Plane3d& plPlane, bool bSaveOnEq)
          return false ;
    }
 
+   //// taglio la trimesh col piano  // questa operazione richiederebbe di ricavare anche i nuovi edge 2D e edge3D, che normalmente sarebbero calcolati durante la GetAuxSurf dal tree
+   //m_pSTM->Cut( plPlane, bSaveOnEq) ;
+
    // imposto ricalcolo della grafica
    m_OGrMgr.Reset() ;
 
@@ -2595,7 +2731,7 @@ SurfBezier::UnprojectPointFromStm( int nT, const Point3d& ptI, Point3d& ptSP, in
    if ( bThroughEdge != nullptr)
       *bThroughEdge = false ;
    // dato un punto sulla trimesh ausiliaria, ne ricavo le coordinate parametriche
-   const ISurfTriMesh* pSurfTm = GetAuxSurf() ;
+   const ISurfTriMesh* pSurfTm = GetAuxSurfRefined() ;
    int nTriaIndex = nT ;
    if ( nT == -1) {
       DistPointSurfTm distPtStm0( ptI, *pSurfTm) ;
@@ -2622,6 +2758,8 @@ SurfBezier::UnprojectPointFromStm( int nT, const Point3d& ptI, Point3d& ptSP, in
                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 {
@@ -2633,9 +2771,7 @@ SurfBezier::UnprojectPointFromStm( int nT, const Point3d& ptI, Point3d& ptSP, in
 
       // se ho tre intersezioni vuol dire che un lato è collassato in un punto e il punto di cui voglio la controimmagine è esattamente nel polo
       // oppure sono su un lato di chiusura
-      if ( nInters == 3  ||  ( m_bClosedU &&  ( vInters[1] == 1 || vInters[3] == 1))  ||  ( m_bClosedV  &&  ( vInters[0] == 1  ||  vInters[2] == 1))) {
-         if ( nInters == 3)
-            bIsPole = true ;
+      if ( bIsPole  ||  ( m_bClosedU &&  ( vInters[1] == 1 || vInters[3] == 1))  ||  ( m_bClosedV  &&  ( vInters[0] == 1  ||  vInters[2] == 1))) {
          // visto che sono in un polo o su un lato di chiusura devo verificare di aver ricevuto il triangolo giusto
          // se è stato passato il punto successivo o precedente mi sposto verso quello e ricalcolo il triangolo di appartenenza
          if ( ! ptIPrevOrNext.IsValid())
@@ -2660,7 +2796,7 @@ SurfBezier::UnprojectPointFromStm( int nT, const Point3d& ptI, Point3d& ptSP, in
    pSurfTm->GetVertex( nVert[0], vPT[0]) ;
    pSurfTm->GetVertex( nVert[1], vPT[1]) ;
    pSurfTm->GetVertex( nVert[2], vPT[2]) ;
-   // devo capire se il trinagolo ha un vertice su un polo
+   // devo capire se il triangolo ha un vertice su un polo
    int nVertOnPole = -1 ;
    int nEdge = -1 ;
    // do per scontato che al più un vertice possa essere in un polo
@@ -2694,13 +2830,16 @@ SurfBezier::UnprojectPointFromStm( int nT, const Point3d& ptI, Point3d& ptSP, in
    }
 
    // se la superficie è chiusa controllo se devo tenere conto della periodicità nel prendere le coordinate parametriche dei vertici
-   double dParamH, dParamL ;
+   double dParamH = 0, dParamL = 0 ;
+   INTVECTOR vOn(3) ;
+   fill( vOn.begin(), vOn.end(), -1) ;
+   bool bOneVertexOnClosureButNotPole = false ;
    if ( m_bClosedU || m_bClosedV) {
       dParamH = m_nSpanV * SBZ_TREG_COEFF ;
       dParamL = m_nSpanU * SBZ_TREG_COEFF ;
       // devo trovare il lato più lungo e confrontarlo con le dimensioni dello spazio parametrico
       Vector3d vtDir ;
-      double dDist ;
+      double dDist = 0 ;
       if ( DistXY( vPtPa[0], vPtPa[1]) > DistXY( vPtPa[1], vPtPa[2])  &&  DistXY( vPtPa[0], vPtPa[1]) > Dist( vPtPa[0], vPtPa[2])){
          vtDir = vPtPa[1] - vPtPa[0] ;
          dDist = DistXY( vPtPa[0], vPtPa[1]) ;
@@ -2719,8 +2858,6 @@ SurfBezier::UnprojectPointFromStm( int nT, const Point3d& ptI, Point3d& ptSP, in
       // do per scontato che se la superficie è chiusa lungo un parametro i lati di chiusura non siano dei poli
       if ( m_bClosedU  &&  abs(vtDir.x) > abs( vtDir.y) && dDist > dParamL * 0.5 ) {
          // trovo se dei vertici del triangolo sono sul bordo dello spazio parametrico
-         INTVECTOR vOn(3) ;
-         fill( vOn.begin(), vOn.end(), -1) ;
          INTVECTOR vEdgesClosed = { 1, 3} ;
          // scorro sui vertici
          for ( int p = 0 ; p < 3; ++p ) {
@@ -2735,12 +2872,12 @@ SurfBezier::UnprojectPointFromStm( int nT, const Point3d& ptI, Point3d& ptSP, in
          }
          // controllo che almeno un vertice sia su un edge e se è l'unico vertice, che non sia su un polo
          if ( vOn[0] > 0  ||  vOn[1] > 0  ||  vOn[2] > 0) {
-            // se ho più un vertice sul lato oppure se ne ho solo uno ma non è sul polo allora procedo alla correzione delle coordinate
-            if ( vOn[0] * vOn[1] * vOn[2] < 0  ||
-               (vOn[0] > 0  &&  0 != nVertOnPole)  ||
-               (vOn[1] > 0  &&  1 != nVertOnPole)  ||
-               (vOn[2] > 0  &&  2 != nVertOnPole)) {
-               double dRightX ;
+            // se ho più di un vertice sul lato oppure se ne ho solo uno ma non è sul polo allora procedo alla correzione delle coordinate
+            bOneVertexOnClosureButNotPole = (vOn[0] > 0  &&  0 != nVertOnPole)  ||
+                                          (vOn[1] > 0  &&  1 != nVertOnPole)  ||
+                                          (vOn[2] > 0  &&  2 != nVertOnPole) ;
+            if ( vOn[0] * vOn[1] * vOn[2] < 0  || bOneVertexOnClosureButNotPole) {
+               double dRightX = 0 ;
                // tengo per buone le coordinate dei vertici che NON sono sul bordo dello spazio parametrico
                for ( int p = 0 ; p < 3; ++p) {
                   if ( vOn[p] == -1) {
@@ -2749,7 +2886,8 @@ SurfBezier::UnprojectPointFromStm( int nT, const Point3d& ptI, Point3d& ptSP, in
                   }
                }
                for ( int p = 0 ; p < 3; ++p) {
-                  if ( abs(vPtPa[p].x - dRightX) > EPS_SMALL ) {
+                  //if ( abs(vPtPa[p].x - dRightX) > EPS_SMALL ) {
+                  if ( abs(vPtPa[p].x - dRightX) > SBZ_TREG_COEFF / 2) {
                      if ( vPtPa[p].x < EPS_SMALL)
                         vPtPa[p].x = dParamL ;
                      else
@@ -2760,8 +2898,6 @@ SurfBezier::UnprojectPointFromStm( int nT, const Point3d& ptI, Point3d& ptSP, in
          }
       }
       else if ( m_bClosedV  &&  abs(vtDir.y) > abs(vtDir.x)  &&  dDist > dParamH * 0.5) {
-         INTVECTOR vOn(3) ;
-         fill( vOn.begin(), vOn.end(), -1) ;
          int nVertOnPole = -1 ;
          INTVECTOR vEdgesClosed = { 0, 2} ;
          // scorro sui vertici
@@ -2778,11 +2914,11 @@ SurfBezier::UnprojectPointFromStm( int nT, const Point3d& ptI, Point3d& ptSP, in
          // controllo che almeno un vertice sia su un edge
          if ( vOn[0] > 0  ||  vOn[1] > 0  ||  vOn[2] > 0) {
             // se ho più un vertice sul lato oppure se ne ho solo uno ma non è sul polo allora procedo alla correzione delle coordinate
-            if ( vOn[0] * vOn[1] * vOn[2] < 0  ||
-               (vOn[0] > 0  &&  0 != nVertOnPole)  ||
-               (vOn[1] > 0  &&  1 != nVertOnPole)  ||
-               (vOn[2] > 0  &&  2 != nVertOnPole)) {
-               double dRightY ;
+            bOneVertexOnClosureButNotPole = (vOn[0] > 0  &&  0 != nVertOnPole)  ||
+                                             (vOn[1] > 0  &&  1 != nVertOnPole)  ||
+                                             (vOn[2] > 0  &&  2 != nVertOnPole) ;
+            if ( vOn[0] * vOn[1] * vOn[2] < 0  ||  bOneVertexOnClosureButNotPole) {
+               double dRightY = 0 ;
                // tengo per buone le coordinate dei vertici che NON sono sul bordo dello spazio parametrico
                for ( int p = 0 ; p < 3; ++p) {
                   if ( vOn[p] == -1) {
@@ -2791,7 +2927,8 @@ SurfBezier::UnprojectPointFromStm( int nT, const Point3d& ptI, Point3d& ptSP, in
                   }
                }
                for ( int p = 0 ; p < 3; ++p) {
-                  if ( abs(vPtPa[p].y - dRightY) > EPS_SMALL) {
+                  //if ( abs(vPtPa[p].y - dRightY) > EPS_SMALL) {
+                  if ( abs(vPtPa[p].y - dRightY) > SBZ_TREG_COEFF / 2) {
                      if ( vPtPa[p].y < EPS_SMALL)
                         vPtPa[p].y = dParamH ;
                      else
@@ -2801,6 +2938,25 @@ SurfBezier::UnprojectPointFromStm( int nT, const Point3d& ptI, Point3d& ptSP, in
             }
          }
       }
+      else {
+         for ( int p = 0 ; p < 3; ++p ) {
+            //scorro sugli edge
+            for (int ed = 0 ; ed < 4 ; ++ed) {
+               // scorro sui tratti che compongono l'edge
+               for ( int i = 0 ; i < int( m_mCCEdge[ed].size()) ; ++i) {
+                  if( m_mCCEdge[ed][i]->IsPointOn( vPT[p])  &&  vOn[p] == -1)
+                     vOn[p] = ed ;
+               }
+            }
+         }
+         // anche se non ho un triangolo molto stretchato potrei avere il punto di polo con coordinate parametriche sballate
+         bOneVertexOnClosureButNotPole =  ((m_bClosedU && (vOn[0] == 1  ||  vOn[0] == 3)) && 0 != nVertOnPole) ||
+                                          ((m_bClosedU && (vOn[1] == 1  ||  vOn[1] == 3)) && 1 != nVertOnPole) ||
+                                          ((m_bClosedU && (vOn[2] == 1  ||  vOn[2] == 3)) && 2 != nVertOnPole) ||
+                                          ((m_bClosedV && (vOn[0] == 0  ||  vOn[0] == 2)) && 0 != nVertOnPole) ||
+                                          ((m_bClosedV && (vOn[1] == 0  ||  vOn[1] == 2)) && 1 != nVertOnPole) ||
+                                          ((m_bClosedV && (vOn[2] == 0  ||  vOn[2] == 2)) && 2 != nVertOnPole) ;
+      }
    }
    // devo anche tener conto della possibilità che i lati siano collassati in poli
    if ( bIsPole  ||  bNearPole) {
@@ -2808,16 +2964,28 @@ SurfBezier::UnprojectPointFromStm( int nT, const Point3d& ptI, Point3d& ptSP, in
       if ( nEdge == -1)
          return false ;
       // trovo la coordinata giusta da tenere ( x o y a seconda dell'edge)
-      double dRightX, dRightY ;
+      double dRightX = 0, dRightY = 0 ;
       for ( int p = 0 ; p < 3; ++p) {
          if ( p != nVertOnPole) {
             if ( nEdge == 0 || nEdge == 2) {
-               dRightX = vPtPa[p].x ;
                dRightY = nEdge == 0 ? dParamH : 0 ;
+               if ( vPtPa[p].x > dRightX )
+                  dRightX= vPtPa[p].x ;
+               // se un vertice è su un lato di chiusura, visto che sono vicino ad un polo, mi accerto che un lato del triangolo sia SUL lato di chiusura
+               if( m_bClosedU  &&  bOneVertexOnClosureButNotPole) {
+                  dRightX = ( dParamL - vPtPa[p].x > vPtPa[p].x ? 0 : dParamL) ;
+                  break ;
+               }
             }
             else if ( nEdge == 1 || nEdge == 3) {
+               // se un vertice è su un lato di chiusura, visto che sono vicino ad un polo, mi accerto che un lato del triangolo sia SUL lato di chiusura
                dRightX = nEdge == 1 ? 0 : dParamL ;
-               dRightY = vPtPa[p].y ;
+               if ( vPtPa[p].y > dRightY )
+                  dRightY = vPtPa[p].y ;
+               if( m_bClosedV  &&  bOneVertexOnClosureButNotPole) {
+                  dRightY = ( dParamH - vPtPa[p].y > vPtPa[p].y ? 0 : dParamH) ;
+                  break ;
+               }
             }
          }
       }
@@ -2891,7 +3059,7 @@ SurfBezier::UnprojectPoint( const Point3d& pt3D, Point3d& ptParam, const Point3d
    // trovato un primo candidato ptParam, ne calcolo l'immagine sulla superficie ( ptBez) e ne calcolo la distanza con il punto pt3D
    // ripeto cercando di avvicinarmi il più possibile
    // per trovare il primo punto trovo il triangolo della trimesh ausiliaria più vicino e il punto più vicino
-   DistPointSurfTm dptSurfTm( pt3D, *GetAuxSurf()) ;
+   DistPointSurfTm dptSurfTm( pt3D, *GetAuxSurfRefined()) ;
    Point3d ptI ; dptSurfTm.GetMinDistPoint( ptI) ;
    if ( ! UnprojectPointFromStm( -1, ptI, ptParam, 5, ptIPrev, bThroughEdge))
       return false ;
@@ -3056,7 +3224,7 @@ SurfBezier::GetLoops( ICRVCOMPOPOVECTOR& vCC, bool bLineOrBezier, int nEdge) con
    // se nEdge non è definito ( == -1) allora restituisco tutti gli edge in ordine
    // se bOpenOrAll è true allora restituisco solo gli edge aperti e che non sono di polo
    if ( m_pSTM == nullptr)
-      GetAuxSurf() ;
+      GetAuxSurfRefined() ;
    
    // se la superficie non è trimmata mi basta recuperare gli edge della superficie
    if ( ! m_bTrimmed) {
@@ -3119,10 +3287,10 @@ SurfBezier::GetLoops( ICRVCOMPOPOVECTOR& vCC, bool bLineOrBezier, int nEdge) con
          bool bInters = false ;
          // devo controllare se i loop si appoggiano a bordi chiusi
          // se ne ho più di uno che si appoggia allora devo vedere se trovo dei gruppi di loop che in realtà in 3d sono un edge unico, ma nel 2d risultano a cavallo dell'edge di chiusura
-         PtrOwner<ICurveComposite> pLoop( m_vCCLoop[l]->Clone()) ;
+         //PtrOwner<ICurveComposite> pLoop( m_vCCLoop[l]->Clone()) ;
          for ( int j = 0 ; j < 4 ; ++j) {
             // il primo è un loop di trim, il secondo è un edge dello spazio parametrico
-            IntersCurveCurve icc( *pLoop, *vEdge[j]) ;
+            IntersCurveCurve icc( *m_vCCLoop[l], *vEdge[j]) ;
             ICCIVECTOR vICC ;
             for ( int i = 0 ; i < int(icc.GetIntersCount()) ; ++i ) {
                IntCrvCrvInfo iccInfo ;
@@ -3142,7 +3310,7 @@ SurfBezier::GetLoops( ICRVCOMPOPOVECTOR& vCC, bool bLineOrBezier, int nEdge) con
          // scorro i loop che fanno intersezioni con l'edge 0 ( V=1) e le riporto sull'edge opposto
          for ( pair<int, ICCIVECTOR> vInters0 : vmInters[0]) {
             // scorro le intersezioni del loop
-            for ( IntCrvCrvInfo inters0 : vInters0.second) {  
+            for ( IntCrvCrvInfo& inters0 : vInters0.second) {  
                if ( ! inters0.bOverlap)
                   continue ;
                Point3d ptStart0 = inters0.IciA[0].ptI ;
@@ -3152,9 +3320,9 @@ SurfBezier::GetLoops( ICRVCOMPOPOVECTOR& vCC, bool bLineOrBezier, int nEdge) con
                ptEnd0.y = 0 ;
                // scorro le intersezioni sull'edge opposto e quando trovo un loop che tra le sue intersezioni contiene uno dei punti che ho appena portato
                // su qusto edge allora cancello la parte comune
-               for ( pair<int, ICCIVECTOR> vInters2 : vmInters[2] ) {
+               for ( auto& vInters2 : vmInters[2] ) {
                   // scorro le intersezioni di questo loop
-                  for ( IntCrvCrvInfo inters2 : vInters2.second) {
+                  for ( IntCrvCrvInfo& inters2 : vInters2.second) {
                      // se lo start o end del loop corrente è compreso tra lo start e l'end di un loop che fa intesezione sull'edge opposto allora 
                      // devo cancellare la parte comune
                      Point3d ptStart2 = inters2.IciA[0].ptI ;
@@ -3245,9 +3413,9 @@ SurfBezier::GetLoops( ICRVCOMPOPOVECTOR& vCC, bool bLineOrBezier, int nEdge) con
 
       if ( m_bClosedU) {
          // scorro i loop che fanno intersezioni con l'edge 1 ( U=0) e le riporto sull'edge opposto
-         for ( pair<int, ICCIVECTOR> vInters1 : vmInters[1]) {
+         for ( auto& vInters1 : vmInters[1]) {
             // scorro le intersezioni del loop
-            for ( IntCrvCrvInfo inters1 : vInters1.second) {
+            for ( IntCrvCrvInfo& inters1 : vInters1.second) {
                if ( ! inters1.bOverlap)
                   continue ;
                Point3d ptStart1 = inters1.IciA[0].ptI ;
@@ -3259,7 +3427,7 @@ SurfBezier::GetLoops( ICRVCOMPOPOVECTOR& vCC, bool bLineOrBezier, int nEdge) con
                // su qusto edge allora cancello la parte comune
                for ( pair<int, ICCIVECTOR> vInters3 : vmInters[3] ) {
                   // scorro le intersezioni di questo loop
-                  for ( IntCrvCrvInfo inters3 : vInters3.second) {
+                  for ( IntCrvCrvInfo& inters3 : vInters3.second) {
                      // se lo start o end del loop corrente è compreso tra lo start e l'end di un loop che fa intesezione sull'edge opposto allora 
                      // devo cancellare la parte comune
                      Point3d ptStart3 = inters3.IciA[0].ptI ;
@@ -3402,7 +3570,7 @@ ICurveComposite*
 SurfBezier::GetSingleEdge3D( bool bLineOrBezier, int nEdge) const
 {
    if ( m_mCCEdge.empty()  &&  bLineOrBezier)
-      GetAuxSurf() ;
+      GetAuxSurfRefined() ;
    // questa funzione dà per scontato che la superficie NON sia trimmata
    if ( nEdge < 0  || nEdge > 3  ||  m_bTrimmed)
       return nullptr ;

SurfBezier.h

@@ -25,6 +25,14 @@
 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
 {
@@ -103,6 +111,7 @@ class SurfBezier : public ISurfBezier, public IGeoObjRW
                { return GetControlWeight( GetInd( nIndU, nIndV), pbOk) ; }
       double GetControlWeight( int nInd, bool* pbOk) const override ;
       bool IsAPoint( void) const override ;
+      bool GetPoint( double dU, double dV, Side nUs, Side nVs, Point3d& ptPos) const override ;
       bool GetPointD1D2( double dU, double dV, Side nUs, Side nVs,
                          Point3d& ptPos,
                          Vector3d* pvtDerU = nullptr, Vector3d* pvtDerV = nullptr,
@@ -113,23 +122,24 @@ class SurfBezier : public ISurfBezier, public IGeoObjRW
                             Vector3d* pvtDerUU = nullptr, Vector3d* pvtDerVV = nullptr, Vector3d* pvtDerUV = nullptr) const override ;
       CurveComposite* GetCurveOnU( double dV) const override ;
       CurveComposite* GetCurveOnV( double dU) const override ;
-      CurveComposite* GetLoop( int nLoop) const override ;                            // nLoop 0-based (1°esterno, successivi interni) 
+      CurveComposite* GetLoop( int nLoop) const override ;                            // nLoop 0-based (1°esterno, successivi interni) 
       bool GetControlCurveOnU( int nIndV, PolyLine& plCtrlU) const override ;
       bool GetControlCurveOnV( int nIndU, PolyLine& plCtrlV) const override ;
       const SurfTriMesh* GetAuxSurf( void) const override ;
+      const SurfTriMesh* GetAuxSurfRefined( void) const override ;
       SurfTriMesh* GetApproxSurf( double dTol, double dSideMin = 100 * EPS_SMALL) const override ;
       // funzione per ottenere la suddivisione dello spazio parametrico nelle celle utilizzate per la triangolazione.
       bool GetLeaves( std::vector<std::tuple<int, Point3d, Point3d>>& vLeaves) const override ;
       bool GetTriangles2D( std::vector<std::tuple<int,Point3d, Point3d, Point3d>>& vTria2D) const override ;
       // funzioni che servono per ricavare l'immagine nel parametrico di un punto appartenente alla trimesh ausiliaria della superficie di Bezier
-      // a nIL si può passare 5 come valore di default
+      // a nIL si può passare 5 come valore di default
       bool UnprojectPointFromStm( int nT, const Point3d& ptI, Point3d& ptSP, int nIL = 5) const override ;
       bool UnprojectPointFromStm( int nT, const Point3d& ptI, Point3d& ptSP, int nIL, const Point3d& ptIPrev, bool* bTroughEdge = nullptr) const override ;
       // restituisce il corrispettivo parametrico di un punto qualunque della trimesh associata alla superficie
-      // ptIPrev è un punto addizionale che precede o segue il punto pt3D nel caso in cui il punto faccia parte di una curva 3d sulla superficie
-      // pPlCut è il piano di taglio su cui dovrebbe giacere il punto raffinato
+      // ptIPrev è un punto addizionale che precede o segue il punto pt3D nel caso in cui il punto faccia parte di una curva 3d sulla superficie
+      // pPlCut è il piano di taglio su cui dovrebbe giacere il punto raffinato
       bool UnprojectPoint( const Point3d& pt3D, Point3d& ptParam, const Point3d& ptIPrev, bool* bTroughEdge = nullptr, const Plane3d* plCut = nullptr) const override ;
-      // pPlCut è il piano di taglio su cui giace la curva
+      // pPlCut è il piano di taglio su cui giace la curva
       bool UnprojectCurveFromStm( const ICurveComposite* pCC, ICRVCOMPOPVECTOR& vpCC, const Plane3d* pPlCut) const override ;
       // funzione per tagliare una superficie di bezier con un piano ( cancello la parte dal lato positivo della normale del piano).
       // bSaveOnEq indica se tenere i triangoli (della trimesh associata) che sono sul piano
@@ -138,8 +148,8 @@ class SurfBezier : public ISurfBezier, public IGeoObjRW
       bool IncreaseUV( double& dU, double dx, bool bUOrV, double* dUVCopy = nullptr, bool bModifyOrig = true) const override ;
       bool IncreaseUV( Point3d& ptUV, Vector3d vtH , Point3d* ptUVCopy, bool bModifyOrig) const override ;
       // funzione che restituisce gli edge della superficie o in forma di linea spezzata o in forma di curva di Bezier
-      // se la superficie è trimmata restituisce i loop dello spazio parametrico in forma di linee spezzate
-      bool GetLoops( ICRVCOMPOPOVECTOR& vCC, bool bLineOrBezier, int nEdge = -1) const override ;  // se la superficie non è trimmata restituisce un vettore di 4 elementi. Se la superficie è chiusa lungo un parametro i lati algi estremi di quel parametro saranno null.
+      // se la superficie è trimmata restituisce i loop dello spazio parametrico in forma di linee spezzate
+      bool GetLoops( ICRVCOMPOPOVECTOR& vCC, bool bLineOrBezier, int nEdge = -1) const override ;  // se la superficie non è trimmata restituisce un vettore di 4 elementi. Se la superficie è chiusa lungo un parametro i lati algi estremi di quel parametro saranno null.
       bool IsPlanar( void) const override ;
       bool CreateByFlatContour( const PolyLine& PL) override ;
       bool CreateByRegion( const POLYLINEVECTOR& vPL) override ;
@@ -203,16 +213,17 @@ class SurfBezier : public ISurfBezier, public IGeoObjRW
       ISurfFlatRegion* CreateTrimRegionFromCuts( ICRVCOMPOPOVECTOR& vpCCOpen, ICRVCOMPOPOVECTOR& vpCCClosed) const ;
       // restituisce il singolo edge della superficie non trimmata
       ICurveComposite* GetSingleEdge3D( bool bLineOrBezier, int nEdge) const ;
-      bool UpdateEdgesFromTree( Tree& tr) const ;
       // funzione che calcola se gli edge sono collassati in poli
       bool CalcPoles( void) const ;
       bool FindMatchByParam( const PolyLine& pl0, const PolyLine& pl1, INTVECTOR& vMatch, int& nLong) const ;
       bool ReorderPntVector( const POLYLINEVECTOR& vPL, bool bTriangulatedIn3D, const PNTVECTOR& vPnt, const POLYLINEVECTOR& vPLToOrd, PNTVECTOR& vPntOrd) const ;
       bool ReorderPntEnhancedVector( const POLYLINEVECTOR& vPL, bool bTriangulatedIn3D, const PNTVECTOR& vPnt, const POLYLINEVECTOR& vPLToOrd, PNTVECTOR& vPntOrd) const ;
+      bool GetBernstein( double dU, int nDegU, DBLVECTOR& vBernU) const ;
 
    private :
       ObjGraphicsMgr  m_OGrMgr ;        // gestore grafica dell'oggetto
-      mutable SurfTriMesh* m_pSTM ;     // superficie trimesh ausiliaria
+      mutable SurfTriMesh* m_pSTM ;     // superficie trimesh ausiliaria per la visualizzazione
+      mutable SurfTriMesh* m_pSTMRefined ; // superficie trimesh ausiliaria raffinata per i calcoli
       Status          m_nStatus ;       // stato
       int             m_nDegU ;         // grado in U
       int             m_nDegV ;         // grado in V
@@ -220,17 +231,24 @@ class SurfBezier : public ISurfBezier, public IGeoObjRW
       int             m_nSpanV ;        // numero di pezze in V
       bool            m_bRat ;          // flag di razionale/polinomiale
       bool            m_bTrimmed ;      // flag per presenza regione di trim
-      mutable bool    m_bClosedU ;      // flag che indica se la superficie è chiusa lungo il parametro U
-      mutable bool    m_bClosedV ;      // flag che indica se la superficie è chiusa lungo il parametro V
+      mutable bool    m_bClosedU ;      // flag che indica se la superficie è chiusa lungo il parametro U
+      mutable bool    m_bClosedV ;      // flag che indica se la superficie è chiusa lungo il parametro V
       mutable BOOLVECTOR m_vbPole ;  // vettore di flag che indicano se i lati sono collassati in dei poli
       PNTVECTOR       m_vPtCtrl ;       // vettore dei punti di controllo
       DBLVECTOR       m_vWeCtrl ;       // vettore dei pesi di controllo
       SurfFlatRegion* m_pTrimReg ;      // eventuale regione di trim
-      int             m_nTempProp[2] ;  // vettore proprietà temporanee
+      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 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 int    m_nIsPlanar ;     // enum che indica se la superficie è piana ( -1, non è stato calcolato)
+      mutable DBLVECTOR      m_vBernU ;
+      mutable PNTVECTOR      m_ptTemp ;
+      mutable DBLVECTOR      m_vBernV ;
+      mutable PNTVECTOR      m_ptTempW ;
+      mutable DBLVECTOR      m_dTempW ;
+      mutable PNTVECTOR      m_vPtWCtrlLoc ;
+      mutable DBLVECTOR      m_vWeCtrlLoc ;
 } ;
 
 //-----------------------------------------------------------------------------

Tree.h

@@ -20,6 +20,15 @@
 #include "/EgtDev/Include/EGkPolyLine.h"
 #include "/EgtDev/Include/EGkChainCurves.h"
 #include <map>
+#include <utility>
+
+struct PairHashInt64 {
+   size_t operator()(const 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 {
@@ -155,7 +164,7 @@ class Cell
       ~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_nFlag2( 0), m_nRightEdgeIn( -1), m_bOnLeftEdge( false), m_bOnTopEdge( false), m_nVertToErase( -1),
+           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) ;
@@ -163,7 +172,7 @@ class Cell
          }
       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_nFlag2( 0), m_nRightEdgeIn( -1), m_bOnLeftEdge( false), m_bOnTopEdge( false), m_nVertToErase( -1),
+           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) ; }
@@ -194,9 +203,11 @@ class Cell
       void SetProcessed( bool bProcessed = true)
               { m_bProcessed = bProcessed ; }
       static bool minorX( const Cell& c1, const Cell& c2)
-                     { return c1.m_ptPbl.x < c2.m_ptPbl.x ; }
+              { return c1.m_ptPbl.x < c2.m_ptPbl.x ; }
       static bool minorY( const Cell& c1, const Cell& c2)
-                     { return c1.m_ptPbl.y < c2.m_ptPbl.y ; }
+              { return c1.m_ptPbl.y < c2.m_ptPbl.y ; }
+      void AddPoly( int nPolyId)
+              { m_vnPolyId.push_back( nPolyId) ;}
 
    public :
       int                  m_nId ;              // Id della cella
@@ -209,9 +220,9 @@ class Cell
       double               m_dSplit ;           // parametro a cui è stata splittata la cella
       int                  m_nChild1 ;          // prima cella figlio
       int                  m_nChild2 ;          // seconda cella figlio
-      int                  m_nFlag ;            // falg che indica la caratterizzazione della cella rispetto ai loop di trim
+      int                  m_nFlag ;            // flag che indica la caratterizzazione della cella rispetto ai loop di trim
                                                 // 0 esterna, 1 intersecata, 2 contiene un loop, 3 intersecata e contenente un loop, 4 contenuta in un loop
-      int                  m_nFlag2 ;           // falg che indica se la cella è stata attraversata durante l'ultima fase del labelling
+      bool                 m_bLabelled ;        // flag che indica se la cella è stata attraversata durante l'ultima fase del labelling
       int                  m_nRightEdgeIn ;     // 0 right edge fuori, 1 right edge dentro, 2 metà e metà
       bool                 m_bOnLeftEdge ;      // flag che indica se la cella è sul lato sinistro ( per superfici chiuse sul parametro U)
       bool                 m_bOnTopEdge ;       // flag che indica se la cella è sul lato top ( per superfici chiuse sul parametro V)
@@ -219,6 +230,7 @@ class Cell
                                                 // ogni elemento del vettore è l'insieme dei punti che caratterizza un attraversamento della cella
       int                  m_nVertToErase ;     // vertice da eliminare dal poligono della cella, in caso di lato sovrapposto ad un lato di polo
                                                 // contati in senso CCW a partire dal bottom left
+      INTVECTOR            m_vnPolyId ;         // indici dei poligoni associati a questa cella nel vettore m_vPolygons del Tree
 
    private :
       Point3d          m_ptPbl ;            // punto bottom left
@@ -240,18 +252,14 @@ class Tree
       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) ;
-      bool GetPolygons( POLYLINEMATRIX& vvPolygons, POLYLINEMATRIX& vvPolygons3d) ;
-      bool GetPolygons( POLYLINEMATRIX& vPolygons, bool bForTriangulation, POLYLINEMATRIX& vvPolygons3d) ;
-      bool GetPolygonsBasic( POLYLINEVECTOR& vPolygons, POLYLINEVECTOR& vPolygonsCorrected,                 // restituisce il poligono corrispondente ad ogni cella foglia dell'albero
-                              bool bForTriangulation, POLYLINEVECTOR& vPolygons3d, INTVECTOR vCells = {}) ; // ad ogni poligono sono stati aggiunti tutti i vertici dei vicini posizionati sui suoi lati
-                                                                                                            // se richiesti per la triangolazione ad alcuni poligoni potrebbero venire tolti dei punti per evitare errori dovuti ad eventuali poli sui bordi del parametrico
-      bool GetPolygonsBasic( POLYLINEVECTOR& vPolygons, POLYLINEVECTOR& vPolygonsCorrected, POLYLINEVECTOR& vPolygons3d) ;
-      bool GetPolygonsBasic( POLYLINEVECTOR& vPolygons, INTVECTOR vCells = {}) ;
+      bool GetPolygons( POLYLINEMATRIX& vvPolygons, POLYLINEMATRIX& vvPolygons3d, 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 ( POLYLINEMATRIX& mPLEdges) ;  // restituisce gli edge 3D come polyline
-      bool GetSplitLoops( POLYLINEVECTOR& vPl) const // funzione che restituisce i loop splitatti ai confini delle celle
-            { for ( int i = 0 ; i < int( m_vPlLoop2D.size()); ++i) vPl.emplace_back( m_vPlLoop2D[i]) ; return true ; };
+      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
       // funzioni da usare per ricostruire tagli che vanno aggiunti allo spazio parametrico
       bool AddCutsToRoot( POLYLINEVECTOR& vCuts) ;          // aggiunge i tagli al tree
@@ -261,7 +269,6 @@ class Tree
       std::vector<bool> GetPoles( void) { return m_vbPole ;} ; // funzione che restituisce i flag che indicano se i lati sono collassati in dei poli
 
    private :
-      bool LimitLoop( PolyLine& pl, POLYLINEVECTOR& vPl, BOOLVECTOR& vbOrientation) const ;              // funzione che limita i loop di trim allo spazio parametrico
       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
@@ -279,7 +286,7 @@ class Tree
       bool FindInters( int& nId, const CurveLine& clTrim, const PolyLine& plPolygon, PNTVECTOR& vptInters, bool bFirstInters = true) ;   // trova le intersezioni tra una cella e una linea di trim
                                                                               // resituisce l'id della cella verso cui la curva di trim esce e il vettore delle intersezioni per la cella successiva con il primo punto
       bool CreateCellPolygons( int nLeafId, POLYLINEMATRIX& vPolygons, POLYLINEMATRIX& vPolygons3d, INTVECTOR& vToCheck, int& nPoly, INTVECTOR& vnParentChunk, const PolyLine& plCell, const PolyLine& plCell3d) ;  // crea i poligoni della cella passata. richiede anche la funzione CreateIslandAndHoles per completare i poligoni.
-      bool CreateIslandAndHoles( int nLeafId, POLYLINEMATRIX& vPolygons, POLYLINEMATRIX& vvPolygons3d, int& nPoly, INTVECTOR& vnParentChunk, bool bForTriangulation, 
+      bool CreateIslandAndHoles( int nLeafId, POLYLINEMATRIX& vPolygons, POLYLINEMATRIX& vvPolygons3d, int& nPoly, INTVECTOR& vnParentChunk, 
                                  const PolyLine& plPolygonsBasic, const PolyLine& plPolygonsBasic3d) ;  // ai poligoni generati da CreatePolygonsCell aggiunge i loop che creano isole o buchi all'interno della singola cella
       bool CheckIfBefore( const PolyLine& pl, int nEdge) const ;  // controllo se ptEnd è prima di ptStart sul lato nEdge rispetto al senso antiorario
       bool CheckIfBefore( const Inters& inA) const ;  // controlla se l'ingresso è prima dell'uscita in senso antiorario a partire da ptTR.
@@ -294,21 +301,19 @@ class Tree
       bool CheckIfBetween( const Inters& inA, const Inters& inB) const ; // / controllo se inB è compreso tra l'end e lo start di inA (in senso CCW)
       bool OnWhichEdge( int nId, const Point3d& ptToAssign, int& nEdge) const ;  // indica a quale edge o vertice il punto è vicino entro EPS_SMALL
       bool AdjustCuts( void) ;
-      bool UpdateSplitLoop( PolyLine& pl, int& nCount, Point3d& pt) ;
-      bool CloseOpenCuts( void) ;
-      bool CloseOpenCuts( POLYLINEVECTOR& vPL, PolyLine& pl) const ;
+      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) ;
 
 
    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
-      //INTMATRIX                   m_vChunk ;           // elenco dei loop divisi per chunk
-      std::unordered_map<int,int>   m_mChunk ;           // mappa in cui vengono salvati chunk di appartenza per ogni loop di trim
-      //ICURVEPOVECTOR              m_vLoop ;            // curve di loop
-      std::vector<std::tuple<PolyLine,bool>> m_vPlApprox ;  // vettore contenente le approssimazioni dei loop  // il bool indica se la curva è CCW
+      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
       bool                        m_bBilinear ;        // superficie bilineare
       bool                        m_bMulti ;           // superficie multi-patch
       bool                        m_bClosedU ;         // superficie chiusa lungo il parametro U
@@ -319,14 +324,11 @@ class Tree
       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
-      POLYLINEMATRIX              m_vPolygons ;        // matrice dei poligoni del tree
-      POLYLINEMATRIX              m_vPolygonsCorr ;    // matrice dei poligoni del tree, corretti per i punti che sono nei poli
-      POLYLINEMATRIX              m_vPolygons3d ;        // matrice dei poligoni3d del tree
-      std::unordered_map<int,Cell>          m_mTree ;            // mappa che contiene tutti i nodi e le foglie dell'albero. -2 è puntatore Null e -1 è root
-      std::unordered_map<int,PNTVECTOR>     m_mVert ;            // mappa che contiene tutti i vertici 3d delle celle del tree. L'Id è lo stesso che la cella ha in m_mTree. I punti sono nell'ordine P00, P10, P11, P01
+      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)
       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
-      POLYLINEVECTOR              m_vPlLoop2D ;        // vettore che contiene le polyline che rappresentano i loop di trim tenendo conto della divisione in celle
-      std::vector<std::pair<BIPNTVECTOR, ChainCurves>>    m_vCEdge2D ;        // vettore che le chain che rappresentano ciò che resta degli edge originali, tenendo conto dei trim.
+      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.
 } ;

Voronoi.cpp

@@ -810,6 +810,7 @@ Voronoi::CalcSpecialPointOffset( PNTVECTVECTOR& vResult, double dOffs)
            // calcolo il punto sul bisettore in corrispondenza dell'offset
             Point3d pt ;
             m_vroni->GetBisectorPointAtParam( i, abs( dOffs), pt.v) ;
+            pt.ToGlob( m_Frame) ;
 
            // calcolo il vettore tangente
             PtrOwner<ICurve> pCrv( GetBisectorCurve( i)) ;