EgtGeomKernel :

- correzione a Inters3Planes (la linea deve essere infinita)
- modifiche varie e correzioni a ProjectCurveOnSurf per gestire spigoli e per far funzionare l'eliminazione dei punti superflui.

IntersPlanePlane.cpp

@@ -42,5 +42,5 @@ Inters3Planes( const Plane3d& plPlane1, const Plane3d& plPlane2, const Plane3d&
    if ( IntersPlanePlane( plPlane1, plPlane2, ptL, vtL) != IPPT_YES)
       return IPPT_NO ;
   // intersezione della linea con il terzo piano
-   return ( IntersLinePlane( ptL, vtL, 1, plPlane3, ptInt) == ILPT_YES ? IPPT_YES : IPPT_NO) ;
+   return ( IntersLinePlane( ptL, vtL, 1, plPlane3, ptInt, false) == ILPT_YES ? IPPT_YES : IPPT_NO) ;
 }

ProjectCurveSurf.cpp

@@ -19,19 +19,187 @@
 #include "/EgtDev/Include/EGkDistPointLine.h"
 #include "/EgtDev/Include/EGkDistPointCurve.h"
 #include "/EgtDev/Include/EGkDistPointSurfTm.h"
+#include "/EgtDev/Include/EGkIntersPlanePlane.h"
+#include "/EgtDev/Include/EGkIntersLinePlane.h"
 #include "/EgtDev/Include/EGkIntersLineSurfTm.h"
 #include "/EgtDev/Include/EGkProjectCurveSurf.h"
 
 using namespace std ;
 
 //----------------------------------------------------------------------------
-// Angolo limite tra normale al triangolo e direzione di proiezione 89°
+// Angolo limite tra normale al triangolo e direzione di proiezione 89°
 const double COS_ANG_LIM = 0.0175 ;
+// Massimo numero di triangoli per raffinare ricerca su spigoli
+const int MAX_FACET_FOR_CORNER = 1000 ;
+// Angolo massimo tra normali per effettuare bisezione su spigolo
+const double COS_ANG_MAX_CORNER = 0.8660 ;
+// Tipologia di punto
+const int P5AX_TO_DELETE = -1 ;        // da cancellare
+const int P5AX_OUT = 0 ;               // aggiunto prima di inizio o dopo fine
+const int P5AX_STD = 1 ;               // standard
+const int P5AX_CVEX = 2 ;              // su angolo convesso
+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 bool
+PointsInTolerance( const PNT5AXVECTOR& vPt5ax, int nPrec, int nCurr, int nNext, double dSqTol)
+{
+   for ( int i = nPrec + 1 ; i < nCurr ; ++ i) {
+      double dSqDist ;
+      if ( ! DistPointLine( vPt5ax[i].ptP, vPt5ax[nPrec].ptP, vPt5ax[nNext].ptP).GetSqDist( dSqDist)  ||  dSqDist > dSqTol)
+         return false ;
+   }
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+static bool
+RemovePointsInExcess( PNT5AXVECTOR& vPt5ax, double dLinTol, double dMaxSegmLen, bool bTestDir)
+{
+  // Parametri di riferimento
+   double dSqMaxLen = dMaxSegmLen * dMaxSegmLen ;
+   double dSqTol = dLinTol * dLinTol ;
+   const double LENREF = 200 ;
+   double dCosAngLim = 1 - dSqTol / ( 2 * LENREF * LENREF) ;
+  // Cerco gli angoli interni e marco opportunamente i punti nelle vicinanze fino ai limiti prima e dopo
+   int nInd = 0 ;
+   while ( nInd  < int( vPt5ax.size())) {
+      if ( vPt5ax[nInd].nFlag == P5AX_CONC) {
+        // analizzo i punti appena precedenti
+         int nIpv = nInd - 1 ;
+         while ( nIpv >= 0) {
+            double dSqLen = SqDist( vPt5ax[nInd].ptP, vPt5ax[nIpv].ptP) ;
+            if ( dSqLen < dSqMaxLen)
+               vPt5ax[nIpv].nFlag = P5AX_TO_DELETE ;
+            else {
+               vPt5ax[nIpv].nFlag = P5AX_BEFORE_CONC ;
+               break ;
+            }
+            -- nIpv ;
+         }
+        // analizzo i punti appena successivi
+         int nInx = nInd + 1 ;
+         while ( nInx < int( vPt5ax.size())) {
+            double dSqLen = SqDist( vPt5ax[nInd].ptP, vPt5ax[nInx].ptP) ;
+            if ( dSqLen < dSqMaxLen)
+               vPt5ax[nInx].nFlag = P5AX_TO_DELETE ;
+            else {
+               vPt5ax[nInx].nFlag = P5AX_AFTER_CONC ;
+               break ;
+            }
+            ++ nInx ;
+         }
+      }
+      ++ nInd ;
+   }
+  // Rimuovo i punti allineati entro la tolleranza e non più lontani tra loro del massimo
+   int nCnt = 0 ;
+   int nPrec = 0 ;
+   int nCurr = 1 ;
+   int nNext = 2 ;
+   while ( nNext < int( vPt5ax.size())) {
+      bool bRemove = false ;
+     // lunghezza del segmento che unisce gli adiacenti
+      double dSqLen = SqDist( vPt5ax[nPrec].ptP, vPt5ax[nNext].ptP) ;
+     // se rimovibile (Flag standard) e lunghezza inferiore al massimo, passo agli altri controlli
+      if ( vPt5ax[nCurr].nFlag == P5AX_STD  &&  dSqLen <= dSqMaxLen) {
+        // distanza del punto corrente dal segmento che unisce gli adiacenti
+         DistPointLine dPL( vPt5ax[nCurr].ptP, vPt5ax[nPrec].ptP, vPt5ax[nNext].ptP) ;
+         double dSqDist ;
+        // se distanza inferiore a tolleranza lineare
+         if ( dPL.GetSqDist( dSqDist)  &&  dSqDist < dSqTol  &&  PointsInTolerance( vPt5ax, nPrec, nCurr, nNext, dSqTol)) {
+           // verifico se errore angolare inferiore a limite
+            double dPar ; dPL.GetParamAtMinDistPoint( dPar) ;
+            if ( bTestDir) {
+               Vector3d vtNew = Media( vPt5ax[nPrec].vtDir1, vPt5ax[nNext].vtDir1, dPar) ;
+               if ( vtNew.Normalize()  &&  vtNew * vPt5ax[nCurr].vtDir1 > dCosAngLim)
+                  bRemove = true ;
+            }
+            else {
+               Vector3d vtNew = Media( vPt5ax[nPrec].vtDir2, vPt5ax[nNext].vtDir2, dPar) ;
+               if ( vtNew.Normalize()  &&  vtNew * vPt5ax[nCurr].vtDir2 > dCosAngLim)
+                  bRemove = true ;
+            }
+         }
+      }
+     // se da eliminare
+      if ( bRemove) {   
+         // dichiaro da eliminare il punto
+         vPt5ax[nCurr].nFlag = P5AX_TO_DELETE ;
+        // avanzo con corrente e successivo
+         nCurr = nNext ;
+         ++ nNext ;
+      }
+     // altrimenti da tenere
+      else {
+        // avanzo il terzetto di uno step
+         nPrec = nCurr ;
+         nCurr = nNext ;
+         ++ nNext ;
+        // incremento contatore dei punti conservati
+         ++ nCnt ;
+      }
+   }
+
+  // Copio i punti da conservare in un vettore temporaneo
+   PNT5AXVECTOR vMyPt5ax ;
+   vMyPt5ax.reserve( nCnt) ;
+   for ( const auto& Pt5ax : vPt5ax) {
+      if ( Pt5ax.nFlag != P5AX_TO_DELETE)
+         vMyPt5ax.emplace_back( Pt5ax) ;
+   }
+  // scambio i due vettori
+   vPt5ax.swap( vMyPt5ax) ;
+
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+static bool
+ProjectPointOnSurf( const Point3d& ptP, const SurfTriMesh& stmSurf, const Frame3d& frRefLine, const IntersParLinesSurfTm& intPLSTM,
+                    double dPar, Point5ax& Pt5ax)
+{
+  // intersezione retta di proiezione con superficie
+   Point3d ptL = GetToLoc( ptP, frRefLine) ;
+   ILSIVECTOR vIntRes ;
+   intPLSTM.GetInters( ptL, 1, vIntRes, false) ;
+  // cerco la prima intersezione valida a partire dall'ultima (è la più alta)
+   int nI = int( vIntRes.size()) - 1 ;
+   while ( nI >= 0  &&  abs( vIntRes[nI].dCosDN) < COS_ANG_LIM)
+      --nI ;
+  // se trovata
+   if ( nI >= 0) {
+     // calcolo il punto
+      Point3d ptInt ;
+      if ( vIntRes[nI].nILTT == ILTT_SEGM  ||  vIntRes[nI].nILTT == ILTT_SEGM_ON_EDGE)
+         ptInt = vIntRes[nI].ptI2 ;
+      else
+         ptInt = vIntRes[nI].ptI ;
+     // calcolo la normale (si calcola smooth, in caso di errore si prende quella del triangolo)
+      Triangle3dEx trTria ;
+      if ( ! stmSurf.GetTriangle( vIntRes[nI].nT, trTria))
+         return false ;
+      Vector3d vtN ;
+      if ( ! CalcNormal( ptInt, trTria, vtN))
+         vtN = trTria.GetN() ;
+     // assegno valori al punto 5assi
+      Pt5ax.ptP = ptInt ;
+      Pt5ax.vtDir1 = vtN ;
+      Pt5ax.vtDir2 = frRefLine.VersZ() ;
+      Pt5ax.dPar = dPar ;
+      Pt5ax.nFlag = P5AX_STD ;
+     // ritorno con successo
+      return true ;
+   }
+   return false ;
+}
 
 //----------------------------------------------------------------------------
 bool
 ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const Vector3d& vtDir,
-                    double dLinTol, double dMaxSegmLen, PNT5AXVECTOR& vPt5ax)
+                    double dLinTol, double dMaxSegmLen, bool bSharpEdges, PNT5AXVECTOR& vPt5ax)
 {
   // sistemazioni per tipo di superficie
    const SurfTriMesh* pSurfTm = nullptr ;
@@ -45,6 +213,8 @@ ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const Vector3d& vt
    case SRF_FLATRGN :
       pSurfTm = GetBasicSurfFlatRegion( &sfSurf)->GetAuxSurf() ;
       break ;
+   default :
+      break ;
    }
    if ( pSurfTm == nullptr)
       return false ;
@@ -66,60 +236,143 @@ ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const Vector3d& vt
       return false ;
    IntersParLinesSurfTm intPLSTM( frRefLine, *pSurfTm) ;
 
-  // Vettore locale dei punti risultanti
-   PNT5AXVECTOR vMyPt5ax ;
-   vMyPt5ax.reserve( PL.GetPointNbr()) ;
+  // Pulisco e riservo spazio nel vettore dei punti risultanti
+   vPt5ax.clear() ;
+   vPt5ax.reserve( PL.GetPointNbr()) ;
 
   // proietto i punti della polilinea sulla superficie
    double dPar ;
    Point3d ptP ;
    bool bFound = PL.GetFirstUPoint( &dPar, &ptP) ;
    while ( bFound) {
-     // intersezione retta di proiezione con superficie
-      Point3d ptL = GetToLoc( ptP, frRefLine) ;
-      ILSIVECTOR vIntRes ;
-      intPLSTM.GetInters( ptL, 1, vIntRes, false) ;
-     // cerco la prima intersezione valida a partire dall'ultima (è la più alta)
-      int nI = int( vIntRes.size()) - 1 ;
-      while ( nI >= 0  &&  abs( vIntRes[nI].dCosDN) < COS_ANG_LIM)
-         --nI ;
-     // se trovata
-      if ( nI >= 0) {
-        // calcolo il punto
-         Point3d ptInt ;
-         if ( vIntRes[nI].nILTT == ILTT_SEGM  ||  vIntRes[nI].nILTT == ILTT_SEGM_ON_EDGE)
-            ptInt = vIntRes[nI].ptI2 ;
-         else
-            ptInt = vIntRes[nI].ptI ;
-        // calcolo la normale (si calcola smooth, in caso di errore si prende quella del triangolo)
-         Triangle3dEx trTria ;
-         if ( ! pSurfTm->GetTriangle( vIntRes[nI].nT, trTria))
-            return false ;
-         Vector3d vtN ;
-         if ( ! CalcNormal( ptInt, trTria, vtN))
-            vtN = trTria.GetN() ;
-        // aggiungo al vettore dei proiettati
-         vMyPt5ax.emplace_back( ptInt, vtN, vtDir, dPar, 1) ;
-      }
+     // se trovo proiezione, la salvo
+      Point5ax Pt5ax ; 
+      if ( ProjectPointOnSurf( ptP, *pSurfTm, frRefLine, intPLSTM, dPar, Pt5ax))
+         vPt5ax.emplace_back( Pt5ax) ;
+     // passo al successivo
       bFound = PL.GetNextUPoint( &dPar, &ptP) ;
    }
 
-   // eventuale rimozione punti in eccesso rispetto alle tolleranze lasciata alla funzione chiamante
-
-  // copio i punti nel vettore di ritorno
-   vPt5ax.clear() ;
-   for ( const auto& Pt5ax : vMyPt5ax) {
-      if ( Pt5ax.nFlag != -1)
-         vPt5ax.emplace_back( Pt5ax) ;
+  // se richiesto, inserimento punti intermedi in presenza di spigoli
+   if ( bSharpEdges) {
+      for ( int i = 1 ; i < int( vPt5ax.size()) ; ++ i) {
+        // precedente
+         int j = i - 1 ;
+        // se normali tra corrente e precedente oltre limite e punti abbastanza lontani
+         if ( vPt5ax[i].vtDir1 * vPt5ax[j].vtDir1 < COS_ANG_MAX_CORNER  &&
+              SqDist( vPt5ax[i].ptP, vPt5ax[j].ptP) > 25 * SQ_EPS_SMALL) {
+           // intersezione tra le due facce
+            Plane3d plPlane1 ; plPlane1.Set( vPt5ax[j].ptP, vPt5ax[j].vtDir1) ;
+            Plane3d plPlane2 ; plPlane2.Set( vPt5ax[i].ptP, vPt5ax[i].vtDir1) ;
+            Point3d ptEdge ; Vector3d vtEdge ;
+            if ( IntersPlanePlane( plPlane1, plPlane2, ptEdge, vtEdge) == IPPT_YES) {
+               Plane3d plPlane3 ; plPlane3.Set( vPt5ax[i].ptP, ( vPt5ax[i].ptP - vPt5ax[j].ptP) ^ vtDir) ;
+               Point3d ptInt ;
+               if ( IntersLinePlane( ptEdge, vtEdge, 1, plPlane3, ptInt, false) == ILPT_YES) {
+                 // verifico se spigolo convesso o concavo
+                  bool bConvex ;
+                  if ( ! AreSamePointApprox( ptInt, vPt5ax[j].ptP))
+                     bConvex = ( ( vPt5ax[j].vtDir1 ^ ( ptInt - vPt5ax[j].ptP)) * vtEdge > 0) ;
+                  else
+                     bConvex = (( vPt5ax[i].vtDir1 ^ ( ptInt - vPt5ax[i].ptP)) * vtEdge < 0) ;
+                 // se convesso, metto due punti con direzione appena prima e appena dopo
+                  if ( bConvex) {
+                     Vector3d vtLine1 = ptInt - vPt5ax[j].ptP ; double dLen1 = vtLine1.Len() ;
+                     Vector3d vtLine2 = vPt5ax[i].ptP - ptInt ; double dLen2 = vtLine2.Len() ;
+                     if ( dLen1 > 2 * EPS_SMALL) {
+                        Point5ax Pt5ax ;
+                        Pt5ax.ptP = ptInt - vtLine1 / dLen1 * 2 * EPS_SMALL ;
+                        Pt5ax.vtDir1 = vPt5ax[j].vtDir1 ;
+                        Pt5ax.vtDir2 = vtDir ;
+                        Pt5ax.dPar = ( vPt5ax[i].dPar + vPt5ax[j].dPar) / 2 ;
+                        Pt5ax.nFlag = P5AX_CVEX ;
+                        vPt5ax.insert( vPt5ax.begin() + i, Pt5ax) ;
+                        ++ i ;
+                     }
+                     if ( dLen2 > 2 * EPS_SMALL) {
+                        Point5ax Pt5ax ;
+                        Pt5ax.ptP = ptInt + vtLine2 / dLen2 * 2 * EPS_SMALL ;
+                        Pt5ax.vtDir1 = vPt5ax[i].vtDir1 ;
+                        Pt5ax.vtDir2 = vtDir ;
+                        Pt5ax.dPar = ( vPt5ax[i].dPar + vPt5ax[j].dPar) / 2 ;
+                        Pt5ax.nFlag = P5AX_CVEX ;
+                        vPt5ax.insert( vPt5ax.begin() + i, Pt5ax) ;
+                        ++ i ;
+                     }
+                  }
+                 // altrimenti concavo, aggiungo un solo punto con la direzione media
+                  else {
+                     Point5ax Pt5ax ;
+                     Pt5ax.ptP = ptInt ;
+                     Pt5ax.vtDir1 = Media( vPt5ax[i].vtDir1, vPt5ax[j].vtDir1) ;
+                     Pt5ax.vtDir2 = vtDir ;
+                     Pt5ax.dPar = ( vPt5ax[i].dPar + vPt5ax[j].dPar) / 2 ;
+                     Pt5ax.nFlag = P5AX_CONC ;
+                     vPt5ax.insert( vPt5ax.begin() + i, Pt5ax) ;
+                     ++ i ;
+                  }
+               }
+            }
+         }
+      }
    }
 
+  // rimozione punti in eccesso rispetto alle tolleranze
+   RemovePointsInExcess( vPt5ax, dLinTol, dMaxSegmLen, bSharpEdges) ;
+
    return true ;
 }
 
+//----------------------------------------------------------------------------
+static bool
+ProjectPointOnSurf( const Point3d& ptP, const SurfTriMesh& stmSurf, const IGeoPoint3d& gpRef, double dPar, Point5ax& Pt5ax)
+{
+  // punto di riferimento
+   Point3d ptMin = gpRef.GetPoint() ;
+  // intersezione della retta di minima distanza con la superficie
+   Vector3d vtLine = ptP - ptMin ;
+   double dLineLen = vtLine.Len() ;
+   if ( dLineLen > EPS_SMALL) {        
+      vtLine /= dLineLen ;
+      ILSIVECTOR vIntRes ;
+      if ( IntersLineSurfTm( ptP, vtLine, dLineLen, stmSurf, vIntRes, false)) {
+        // cerco la prima intersezione valida a partire dall'ultima (è la più alta)
+         int nI = int( vIntRes.size()) - 1 ;
+         while ( nI >= 0  &&  abs( vIntRes[nI].dCosDN) < COS_ANG_LIM)
+            --nI ;
+        // se trovata
+         if ( nI >= 0) {
+           // calcolo il punto
+            Point3d ptInt ;
+            if ( vIntRes[nI].nILTT == ILTT_SEGM  ||  vIntRes[nI].nILTT == ILTT_SEGM_ON_EDGE)
+               ptInt = vIntRes[nI].ptI2 ;
+            else
+               ptInt = vIntRes[nI].ptI ;
+           // calcolo la normale (si calcola smooth, in caso di errore si prende quella del triangolo)
+            Triangle3dEx trTria ;
+            if ( ! stmSurf.GetTriangle( vIntRes[nI].nT, trTria))
+               return false ;
+            Vector3d vtN ;
+            if ( ! CalcNormal( ptInt, trTria, vtN))
+               vtN = trTria.GetN() ;
+           // assegno valori al punto 5assi
+            Pt5ax.ptP = ptInt ;
+            Pt5ax.vtDir1 = vtN ;
+            Pt5ax.vtDir2 = vtLine ;
+            Pt5ax.dPar = dPar ;
+            Pt5ax.nFlag = P5AX_STD ;
+           // ritorno con successo
+            return true ;
+         }
+      }
+   }
+   return false ;
+}
+
 //----------------------------------------------------------------------------
 bool
 ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const IGeoPoint3d& gpRef,
-                    double dLinTol, double dMaxSegmLen, PNT5AXVECTOR& vPt5ax)
+                    double dLinTol, double dMaxSegmLen, bool bSharpEdges, PNT5AXVECTOR& vPt5ax)
 {
   // sistemazioni per tipo di superficie
    const SurfTriMesh* pSurfTm = nullptr ;
@@ -133,6 +386,8 @@ ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const IGeoPoint3d&
    case SRF_FLATRGN :
       pSurfTm = GetBasicSurfFlatRegion( &sfSurf)->GetAuxSurf() ;
       break ;
+   default :
+      break ;
    }
    if ( pSurfTm == nullptr)
       return false ;
@@ -141,7 +396,7 @@ ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const IGeoPoint3d&
    dLinTol = max( dLinTol, LIN_TOL_MIN) ;
    dMaxSegmLen = max( dMaxSegmLen, 10 * EPS_SMALL) ;
 
-  // approssimo la curva con una polilinea entro la metà della tolleranza
+  // approssimo la curva con una polilinea entro la metà della tolleranza
    PolyLine PL ;
    if ( ! crCrv.ApproxWithLines( dLinTol, ANG_TOL_STD_DEG, ICurve::APL_STD, PL))
       return false ;
@@ -149,25 +404,67 @@ ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const IGeoPoint3d&
    if ( ! PL.AdjustForMaxSegmentLen( MAX_SEG_LEN))
       return false ;
 
-  // Vettore locale dei punti risultanti
-   PNT5AXVECTOR vMyPt5ax ;
-   vMyPt5ax.reserve( PL.GetPointNbr()) ;
+  // Pulisco e riservo spazio nel vettore dei punti risultanti
+   vPt5ax.clear() ;
+   vPt5ax.reserve( PL.GetPointNbr()) ;
 
   // proietto i punti della polilinea sulla superficie con direzione data dal punto di riferimento
    double dPar ;
    Point3d ptP ;
    bool bFound = PL.GetFirstUPoint( &dPar, &ptP) ;
    while ( bFound) {
-     // punto di riferimento
-      Point3d ptMin = gpRef.GetPoint() ;
+     // se trovo proiezione, la salvo
+      Point5ax Pt5ax ; 
+      if ( ProjectPointOnSurf( ptP, *pSurfTm, gpRef, dPar, Pt5ax))
+         vPt5ax.emplace_back( Pt5ax) ;
+     // passo al successivo
+      bFound = PL.GetNextUPoint( &dPar, &ptP) ;
+   }
+
+  // se superfici con non troppi triangoli, inserimento punti intermedi in presenza di spigoli
+   if ( pSurfTm->GetFacetCount() < MAX_FACET_FOR_CORNER) {
+      for ( int i = 1 ; i < int( vPt5ax.size()) ; ++ i) {
+        // precedente
+         int j = i - 1 ;
+        // se normali tra corrente e precedente oltre limite e punti abbastanza lontani
+         if ( vPt5ax[i].vtDir1 * vPt5ax[j].vtDir1 < COS_ANG_MAX_CORNER  &&
+              SqDist( vPt5ax[i].ptP, vPt5ax[j].ptP) > 25 * SQ_EPS_SMALL) {
+           // punto medio
+            Point3d ptMid = Media( vPt5ax[i].ptP, vPt5ax[j].ptP) ;
+            double dMid = ( vPt5ax[i].dPar + vPt5ax[j].dPar) / 2 ;
+           // se trovo proiezione, la salvo
+            Point5ax Pt5ax ; 
+            if ( ProjectPointOnSurf( ptMid, *pSurfTm, gpRef, dMid, Pt5ax)) {
+               vPt5ax.insert( vPt5ax.begin() + i, Pt5ax) ;
+               -- i ;
+            }
+         }
+      }
+   }
+
+  // rimozione punti in eccesso rispetto alle tolleranze
+   RemovePointsInExcess( vPt5ax, dLinTol, dMaxSegmLen, bSharpEdges) ;
+
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+static bool
+ProjectPointOnSurf( const Point3d& ptP, const SurfTriMesh& stmSurf, const ICurve& crRef, double dPar, Point5ax& Pt5ax)
+{
+  // punto a minima distanza
+   DistPointCurve dPC( ptP, crRef) ;
+   Point3d ptMin ;
+   int nFlag ;
+   if ( dPC.GetMinDistPoint( 0, ptMin, nFlag)) {
      // intersezione della retta di minima distanza con la superficie
       Vector3d vtLine = ptP - ptMin ;
       double dLineLen = vtLine.Len() ;
       if ( dLineLen > EPS_SMALL) {        
          vtLine /= dLineLen ;
          ILSIVECTOR vIntRes ;
-         if ( IntersLineSurfTm( ptP, vtLine, dLineLen, *pSurfTm, vIntRes, false)) {
-           // cerco la prima intersezione valida a partire dall'ultima (è la più alta)
+         if ( IntersLineSurfTm( ptP, vtLine, dLineLen, stmSurf, vIntRes, false)) {
+           // cerco la prima intersezione valida a partire dall'ultima (è la più alta)
             int nI = int( vIntRes.size()) - 1 ;
             while ( nI >= 0  &&  abs( vIntRes[nI].dCosDN) < COS_ANG_LIM)
                --nI ;
@@ -181,37 +478,32 @@ ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const IGeoPoint3d&
                   ptInt = vIntRes[nI].ptI ;
               // calcolo la normale (si calcola smooth, in caso di errore si prende quella del triangolo)
                Triangle3dEx trTria ;
-               if ( ! pSurfTm->GetTriangle( vIntRes[nI].nT, trTria))
+               if ( ! stmSurf.GetTriangle( vIntRes[nI].nT, trTria))
                   return false ;
                Vector3d vtN ;
                if ( ! CalcNormal( ptInt, trTria, vtN))
                   vtN = trTria.GetN() ;
-              // aggiungo al vettore dei proiettati
-               vMyPt5ax.emplace_back( ptInt, vtN, vtLine, dPar, 1) ;
+              // assegno valori al punto 5assi
+               Pt5ax.ptP = ptInt ;
+               Pt5ax.vtDir1 = vtN ;
+               Pt5ax.vtDir2 = vtLine ;
+               Pt5ax.dPar = dPar ;
+               Pt5ax.nFlag = P5AX_STD ;
+              // ritorno con successo
+               return true ;
             }
          }
       }
-      bFound = PL.GetNextUPoint( &dPar, &ptP) ;
    }
-
-  // eventuale rimozione punti in eccesso rispetto alle tolleranze lasciata alla funzione chiamante
-
-  // copio i punti nel vettore di ritorno
-   vPt5ax.clear() ;
-   for ( const auto& Pt5ax : vMyPt5ax) {
-      if ( Pt5ax.nFlag != -1)
-         vPt5ax.emplace_back( Pt5ax) ;
-   }
-
-   return true ;
+   return false ;
 }
 
 //----------------------------------------------------------------------------
 bool
 ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const ICurve& crRef,
-                    double dLinTol, double dMaxSegmLen, PNT5AXVECTOR& vPt5ax)
+                    double dLinTol, double dMaxSegmLen, bool bSharpEdges, PNT5AXVECTOR& vPt5ax)
 {
-  // sistemazioni per tipo di superficie
+  // Sistemazioni per tipo di superficie
    const SurfTriMesh* pSurfTm = nullptr ;
    switch ( sfSurf.GetType()) {
    case SRF_TRIMESH :
@@ -223,15 +515,17 @@ ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const ICurve& crRe
    case SRF_FLATRGN :
       pSurfTm = GetBasicSurfFlatRegion( &sfSurf)->GetAuxSurf() ;
       break ;
+   default :
+      break ;
    }
    if ( pSurfTm == nullptr)
       return false ;
 
-  // controllo le tolleranze
+  // Controllo le tolleranze
    dLinTol = max( dLinTol, LIN_TOL_MIN) ;
    dMaxSegmLen = max( dMaxSegmLen, 10 * EPS_SMALL) ;
 
-  // approssimo la curva con una polilinea alla massima risoluzione
+  // Approssimo la curva con una polilinea alla massima risoluzione
    PolyLine PL ;
    if ( ! crCrv.ApproxWithLines( EPS_SMALL, ANG_TOL_STD_DEG, ICurve::APL_STD, PL))
       return false ;
@@ -239,71 +533,163 @@ ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const ICurve& crRe
    if ( ! PL.AdjustForMaxSegmentLen( MAX_SEG_LEN))
       return false ;
 
-  // Vettore locale dei punti risultanti
-   PNT5AXVECTOR vMyPt5ax ;
-   vMyPt5ax.reserve( PL.GetPointNbr()) ;
+  // Pulisco e riservo spazio nel vettore dei punti risultanti
+   vPt5ax.clear() ;
+   vPt5ax.reserve( PL.GetPointNbr()) ;
 
   // proietto i punti della polilinea sulla superficie con direzione normale alla curva di riferimento
    double dPar ;
    Point3d ptP ;
    bool bFound = PL.GetFirstUPoint( &dPar, &ptP) ;
    while ( bFound) {
-     // punto sulla curva a minima distanza
-      DistPointCurve dPC( ptP, crRef) ;
-      Point3d ptMin ;
-      int nFlag ;
-      if ( dPC.GetMinDistPoint( 0, ptMin, nFlag)) {
-        // intersezione della retta di minima distanza con la superficie
-         Vector3d vtLine = ptP - ptMin ;
-         double dLineLen = vtLine.Len() ;
-         if ( dLineLen > EPS_SMALL) {        
-            vtLine /= dLineLen ;
-            ILSIVECTOR vIntRes ;
-            if ( IntersLineSurfTm( ptP, vtLine, dLineLen, *pSurfTm, vIntRes, false)) {
-              // cerco la prima intersezione valida a partire dall'ultima (è la più alta)
-               int nI = int( vIntRes.size()) - 1 ;
-               while ( nI >= 0  &&  abs( vIntRes[nI].dCosDN) < COS_ANG_LIM)
-                  --nI ;
-              // se trovata
-               if ( nI >= 0) {
-                 // calcolo il punto
-                  Point3d ptInt ;
-                  if ( vIntRes[nI].nILTT == ILTT_SEGM  ||  vIntRes[nI].nILTT == ILTT_SEGM_ON_EDGE)
-                     ptInt = vIntRes[nI].ptI2 ;
+     // se trovo proiezione, la salvo
+      Point5ax Pt5ax ; 
+      if ( ProjectPointOnSurf( ptP, *pSurfTm, crRef, dPar, Pt5ax))
+         vPt5ax.emplace_back( Pt5ax) ;
+     // passo al successivo
+      bFound = PL.GetNextUPoint( &dPar, &ptP) ;
+   }
+
+  // se richiesto, inserimento punti intermedi in presenza di spigoli
+   if ( bSharpEdges) {
+      for ( int i = 1 ; i < int( vPt5ax.size()) ; ++ i) {
+        // precedente
+         int j = i - 1 ;
+        // se normali tra corrente e precedente oltre limite e punti abbastanza lontani
+         if ( vPt5ax[i].vtDir1 * vPt5ax[j].vtDir1 < COS_ANG_MAX_CORNER  &&
+              SqDist( vPt5ax[i].ptP, vPt5ax[j].ptP) > 25 * SQ_EPS_SMALL) {
+           // intersezione tra le due facce
+            Plane3d plPlane1 ; plPlane1.Set( vPt5ax[j].ptP, vPt5ax[j].vtDir1) ;
+            Plane3d plPlane2 ; plPlane2.Set( vPt5ax[i].ptP, vPt5ax[i].vtDir1) ;
+            Point3d ptEdge ; Vector3d vtEdge ;
+            if ( IntersPlanePlane( plPlane1, plPlane2, ptEdge, vtEdge) == IPPT_YES) {
+               Plane3d plPlane3 ; plPlane3.Set( vPt5ax[i].ptP, vPt5ax[i].vtDir2 ^ vPt5ax[j].vtDir2) ;
+               Point3d ptInt ;
+               if ( IntersLinePlane( ptEdge, vtEdge, 1, plPlane3, ptInt, false) == ILPT_YES) {
+                 // verifico se spigolo convesso o concavo
+                  bool bConvex ;
+                  if ( ! AreSamePointApprox( ptInt, vPt5ax[j].ptP))
+                     bConvex = ( ( vPt5ax[j].vtDir1 ^ ( ptInt - vPt5ax[j].ptP)) * vtEdge > 0) ;
                   else
-                     ptInt = vIntRes[nI].ptI ;
-                 // calcolo la normale (si calcola smooth, in caso di errore si prende quella del triangolo)
-                  Triangle3dEx trTria ;
-                  if ( ! pSurfTm->GetTriangle( vIntRes[nI].nT, trTria))
-                     return false ;
-                  Vector3d vtN ;
-                  if ( ! CalcNormal( ptInt, trTria, vtN))
-                     vtN = trTria.GetN() ;
-                 // aggiungo al vettore dei proiettati
-                  vMyPt5ax.emplace_back( ptInt, vtN, vtLine, dPar, 1) ;
+                     bConvex = (( vPt5ax[i].vtDir1 ^ ( ptInt - vPt5ax[i].ptP)) * vtEdge < 0) ;
+                 // se convesso, metto due punti con direzione appena prima e appena dopo
+                  if ( bConvex) {
+                     Vector3d vtLine1 = ptInt - vPt5ax[j].ptP ; double dLen1 = vtLine1.Len() ;
+                     Vector3d vtLine2 = vPt5ax[i].ptP - ptInt ; double dLen2 = vtLine2.Len() ;
+                     if ( dLen1 > 2 * EPS_SMALL) {
+                        Point5ax Pt5ax ;
+                        Pt5ax.ptP = ptInt - vtLine1 / dLen1 * 2 * EPS_SMALL ;
+                        Pt5ax.vtDir1 = vPt5ax[j].vtDir1 ;
+                        Pt5ax.vtDir2 = Media( vPt5ax[i].vtDir2, vPt5ax[j].vtDir2) ;
+                        Pt5ax.dPar = ( vPt5ax[i].dPar + vPt5ax[j].dPar) / 2 ;
+                        Pt5ax.nFlag = P5AX_CVEX ;
+                        vPt5ax.insert( vPt5ax.begin() + i, Pt5ax) ;
+                        ++ i ;
+                     }
+                     if ( dLen2 > 2 * EPS_SMALL) {
+                        Point5ax Pt5ax ;
+                        Pt5ax.ptP = ptInt + vtLine2 / dLen2 * 2 * EPS_SMALL ;
+                        Pt5ax.vtDir1 = vPt5ax[i].vtDir1 ;
+                        Pt5ax.vtDir2 = Media( vPt5ax[i].vtDir2, vPt5ax[j].vtDir2) ;
+                        Pt5ax.dPar = ( vPt5ax[i].dPar + vPt5ax[j].dPar) / 2 ;
+                        Pt5ax.nFlag = P5AX_CVEX ;
+                        vPt5ax.insert( vPt5ax.begin() + i, Pt5ax) ;
+                        ++ i ;
+                     }
+                  }
+                 // altrimenti concavo, aggiungo un solo punto con la direzione media
+                  else {
+                     Point5ax Pt5ax ;
+                     Pt5ax.ptP = ptInt ;
+                     Pt5ax.vtDir1 = Media( vPt5ax[i].vtDir1, vPt5ax[j].vtDir1) ;
+                     Pt5ax.vtDir2 = Media( vPt5ax[i].vtDir2, vPt5ax[j].vtDir2) ;
+                     Pt5ax.dPar = ( vPt5ax[i].dPar + vPt5ax[j].dPar) / 2 ;
+                     Pt5ax.nFlag = P5AX_CONC ;
+                     vPt5ax.insert( vPt5ax.begin() + i, Pt5ax) ;
+                     ++ i ;
+                  }
                }
             }
          }
       }
-      bFound = PL.GetNextUPoint( &dPar, &ptP) ;
    }
 
-  // eventuale rimozione punti in eccesso rispetto alle tolleranze lasciata alla funzione chiamante
-
-  // copio i punti nel vettore di ritorno
-   vPt5ax.clear() ;
-   for ( const auto& Pt5ax : vMyPt5ax) {
-      if ( Pt5ax.nFlag != -1)
-         vPt5ax.emplace_back( Pt5ax) ;
-   }
+  // rimozione punti in eccesso rispetto alle tolleranze
+   RemovePointsInExcess( vPt5ax, dLinTol, dMaxSegmLen, bSharpEdges) ;
 
    return true ;
 }
 
+//----------------------------------------------------------------------------
+static bool
+ProjectPointOnSurf( const Point3d& ptP, const SurfTriMesh& stmSurf, const SurfTriMesh& stmRef, double dPar, Point5ax& Pt5ax)
+{
+  // punto sulla superficie guida a minima distanza
+   DistPointSurfTm dPS( ptP, stmRef) ;
+   Point3d ptMin ;
+   int nTriaMin ;
+   if ( dPS.GetMinDistPoint( ptMin)  &&  dPS.GetMinDistTriaIndex ( nTriaMin)) {
+     // recupero direzione della retta di minima distanza, altrimenti normale alla superficie
+      Vector3d vtLine = ptP - ptMin ;
+      double dLineLen = vtLine.Len() ;
+      if ( dLineLen > EPS_SMALL)       
+         vtLine /= dLineLen ;
+      else {
+        // calcolo la normale della superficie guida
+         Triangle3dEx trGuide ;
+         if ( ! stmRef.GetTriangle( nTriaMin, trGuide))
+            return false ;
+         if ( ! CalcNormal( ptMin, trGuide, vtLine))
+            vtLine = trGuide.GetN() ;
+         dLineLen = 100 ;
+      }
+     // intersezione della retta con la superficie
+      ILSIVECTOR vIntRes ;
+      if ( IntersLineSurfTm( ptP, vtLine, dLineLen, stmSurf, vIntRes, false)) {
+        // cerco la prima intersezione valida a partire dall'ultima (è la più alta)
+         int nI = int( vIntRes.size()) - 1 ;
+         while ( nI >= 0  &&  abs( vIntRes[nI].dCosDN) < COS_ANG_LIM)
+            --nI ;
+        // se trovata
+         if ( nI >= 0) {
+           // calcolo il punto
+            Point3d ptInt ;
+            if ( vIntRes[nI].nILTT == ILTT_SEGM  ||  vIntRes[nI].nILTT == ILTT_SEGM_ON_EDGE)
+               ptInt = vIntRes[nI].ptI2 ;
+            else
+               ptInt = vIntRes[nI].ptI ;
+           // calcolo la normale (si calcola smooth, in caso di errore si prende quella del triangolo)
+            Triangle3dEx trTria ;
+            if ( ! stmSurf.GetTriangle( vIntRes[nI].nT, trTria))
+               return false ;
+            Vector3d vtN ;
+            if ( ! CalcNormal( ptMin, trTria, vtN))
+               vtN = trTria.GetN() ;
+           // calcolo la normale della superficie guida
+            Triangle3dEx trGuide ;
+            if ( ! stmRef.GetTriangle( nTriaMin, trGuide))
+               return false ;
+            Vector3d vtN2 ;
+            if ( ! CalcNormal( ptMin, trGuide, vtN2))
+               vtN2 = trGuide.GetN() ;
+           // assegno valori al punto 5assi
+            Pt5ax.ptP = ptInt ;
+            Pt5ax.vtDir1 = vtN ;
+            Pt5ax.vtDir2 = vtN2 ;
+            Pt5ax.dPar = dPar ;
+            Pt5ax.nFlag = P5AX_STD ;
+           // ritorno con successo
+            return true ;
+         }
+      }
+   }
+   return false ;
+}
+
 //----------------------------------------------------------------------------
 bool
 ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const ISurf& sfRef,
-                    double dLinTol, double dMaxSegmLen, PNT5AXVECTOR& vPt5ax)
+                    double dLinTol, double dMaxSegmLen, bool bSharpEdges, PNT5AXVECTOR& vPt5ax)
 {
   // sistemazioni per tipo di superficie
    const SurfTriMesh* pSurfTm = nullptr ;
@@ -317,6 +703,8 @@ ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const ISurf& sfRef
    case SRF_FLATRGN :
       pSurfTm = GetBasicSurfFlatRegion( &sfSurf)->GetAuxSurf() ;
       break ;
+   default :
+      break ;
    }
    if ( pSurfTm == nullptr)
       return false ;
@@ -333,6 +721,8 @@ ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const ISurf& sfRef
    case SRF_FLATRGN :
       pRefTm = GetBasicSurfFlatRegion( &sfRef)->GetAuxSurf() ;
       break ;
+   default :
+      break ;
    }
    if ( pRefTm == nullptr)
       return false ;
@@ -341,7 +731,7 @@ ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const ISurf& sfRef
    dLinTol = max( dLinTol, LIN_TOL_MIN) ;
    dMaxSegmLen = max( dMaxSegmLen, 10 * EPS_SMALL) ;
 
-  // approssimo la curva con una polilinea entro la metà della tolleranza
+  // approssimo la curva con una polilinea entro la metà della tolleranza
    PolyLine PL ;
    if ( ! crCrv.ApproxWithLines( dLinTol, ANG_TOL_STD_DEG, ICurve::APL_STD, PL))
       return false ;
@@ -349,147 +739,46 @@ ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const ISurf& sfRef
    if ( ! PL.AdjustForMaxSegmentLen( MAX_SEG_LEN))
       return false ;
 
-  // Vettore locale dei punti risultanti
-   PNT5AXVECTOR vMyPt5ax ;
-   vMyPt5ax.reserve( PL.GetPointNbr()) ;
+  // Pulisco e riservo spazio nel vettore dei punti risultanti
+   vPt5ax.clear() ;
+   vPt5ax.reserve( PL.GetPointNbr()) ;
 
   // proietto i punti della polilinea sulla superficie con direzione normale alla curva di riferimento
    double dPar ;
    Point3d ptP ;
    bool bFound = PL.GetFirstUPoint( &dPar, &ptP) ;
    while ( bFound) {
-     // punto sulla superficie guida a minima distanza
-      DistPointSurfTm dPS( ptP, *pRefTm) ;
-      Point3d ptMin ;
-      int nTriaMin ;
-      if ( dPS.GetMinDistPoint( ptMin)  &&  dPS.GetMinDistTriaIndex ( nTriaMin)) {
-        // recupero direzione della retta di minima distanza, altrimenti normale alla superficie
-         Vector3d vtLine = ptP - ptMin ;
-         double dLineLen = vtLine.Len() ;
-         if ( dLineLen > EPS_SMALL)       
-            vtLine /= dLineLen ;
-         else {
-           // calcolo la normale della superficie guida
-            Triangle3dEx trGuide ;
-            if ( ! pRefTm->GetTriangle( nTriaMin, trGuide))
-               return false ;
-            if ( ! CalcNormal( ptMin, trGuide, vtLine))
-               vtLine = trGuide.GetN() ;
-            dLineLen = 100 ;
-         }
-        // intersezione della retta con la superficie
-         ILSIVECTOR vIntRes ;
-         if ( IntersLineSurfTm( ptP, vtLine, dLineLen, *pSurfTm, vIntRes, false)) {
-           // cerco la prima intersezione valida a partire dall'ultima (è la più alta)
-            int nI = int( vIntRes.size()) - 1 ;
-            while ( nI >= 0  &&  abs( vIntRes[nI].dCosDN) < COS_ANG_LIM)
-               --nI ;
-           // se trovata
-            if ( nI >= 0) {
-              // calcolo il punto
-               Point3d ptInt ;
-               if ( vIntRes[nI].nILTT == ILTT_SEGM  ||  vIntRes[nI].nILTT == ILTT_SEGM_ON_EDGE)
-                  ptInt = vIntRes[nI].ptI2 ;
-               else
-                  ptInt = vIntRes[nI].ptI ;
-              // calcolo la normale (si calcola smooth, in caso di errore si prende quella del triangolo)
-               Triangle3dEx trTria ;
-               if ( ! pSurfTm->GetTriangle( vIntRes[nI].nT, trTria))
-                  return false ;
-               Vector3d vtN ;
-               if ( ! CalcNormal( ptMin, trTria, vtN))
-                  vtN = trTria.GetN() ;
-              // calcolo la normale della superficie guida
-               Triangle3dEx trGuide ;
-               if ( ! pRefTm->GetTriangle( nTriaMin, trGuide))
-                  return false ;
-               Vector3d vtN2 ;
-               if ( ! CalcNormal( ptMin, trGuide, vtN2))
-                  vtN2 = trGuide.GetN() ;
-              // aggiungo al vettore dei proiettati
-               vMyPt5ax.emplace_back( ptInt, vtN, vtN2, dPar, 1) ;
-            }
-         }
-      }
+     // se trovo proiezione, la salvo
+      Point5ax Pt5ax ; 
+      if ( ProjectPointOnSurf( ptP, *pSurfTm, *pRefTm, dPar, Pt5ax))
+         vPt5ax.emplace_back( Pt5ax) ;
+     // passo al successivo
       bFound = PL.GetNextUPoint( &dPar, &ptP) ;
    }
 
-  // eventuale rimozione punti in eccesso rispetto alle tolleranze lasciata alla funzione chiamante
-
-  // copio i punti nel vettore di ritorno
-   vPt5ax.clear() ;
-   for ( const auto& Pt5ax : vMyPt5ax) {
-      if ( Pt5ax.nFlag != -1)
-         vPt5ax.emplace_back( Pt5ax) ;
-   }
-
-   return true ;
-}
-
-//----------------------------------------------------------------------------
-static bool
-PointsInTolerance( const PNT5AXVECTOR& vPt5ax, int nPrec, int nCurr, int nNext, double dSqTol)
-{
-   for ( int i = nPrec + 1 ; i < nCurr ; ++ i) {
-      double dSqDist ;
-      if ( ! DistPointLine( vPt5ax[i].ptP, vPt5ax[nPrec].ptP, vPt5ax[nNext].ptP).GetSqDist( dSqDist)  ||  dSqDist > dSqTol)
-         return false ;
-   }
-   return true ;
-}
-
-//----------------------------------------------------------------------------
-bool
-RemovePointsInExcess( PNT5AXVECTOR& vMyPt5ax, double dLinTol, double dMaxSegmLen, bool bTestDir)
-{
-  // rimuovo i punti allineati entro la tolleranza e non più lontani tra loro del massimo
-   double dSqMaxLen = dMaxSegmLen * dMaxSegmLen ;
-   double dSqTol = dLinTol * dLinTol ;
-   const double LENREF = 100 ;
-   double dCosAngLim = 1 - dSqTol / ( 2 * LENREF * LENREF) ;
-   int nPrec = 0 ;
-   int nCurr = 1 ;
-   int nNext = 2 ;
-   while ( nNext < int( vMyPt5ax.size())) {
-      bool bRemove = false ;
-     // lunghezza del segmento che unisce gli adiacenti
-      double dSqLen = SqDist( vMyPt5ax[nPrec].ptP, vMyPt5ax[nNext].ptP) ;
-     // se lunghezza inferiore al massimo, passo agli altri controlli
-      if ( dSqLen <= dSqMaxLen) {
-        // distanza del punto corrente dal segmento che unisce gli adiacenti
-         DistPointLine dPL( vMyPt5ax[nCurr].ptP, vMyPt5ax[nPrec].ptP, vMyPt5ax[nNext].ptP) ;
-         double dSqDist ;
-        // se distanza inferiore a tolleranza lineare
-         if ( dPL.GetSqDist( dSqDist)  &&  dSqDist < dSqTol  &&  PointsInTolerance( vMyPt5ax, nPrec, nCurr, nNext, dSqTol)) {
-           // verifico se errore angolare inferiore a limite
-            double dPar ; dPL.GetParamAtMinDistPoint( dPar) ;
-            if ( bTestDir) {
-               Vector3d vtNew = Media( vMyPt5ax[nPrec].vtDir, vMyPt5ax[nNext].vtDir, dPar) ;
-               if ( vtNew.Normalize()  &&  vtNew * vMyPt5ax[nCurr].vtDir > dCosAngLim)
-                  bRemove = true ;
-            }
-            else {
-               Vector3d vtNew = Media( vMyPt5ax[nPrec].vtDir2, vMyPt5ax[nNext].vtDir2, dPar) ;
-               if ( vtNew.Normalize()  &&  vtNew * vMyPt5ax[nCurr].vtDir2 > dCosAngLim)
-                  bRemove = true ;
+  // se superfici con non troppi triangoli, inserimento punti intermedi in presenza di spigoli
+   if ( pSurfTm->GetFacetCount() < MAX_FACET_FOR_CORNER) {
+      for ( int i = 1 ; i < int( vPt5ax.size()) ; ++ i) {
+        // precedente
+         int j = i - 1 ;
+        // se normali tra corrente e precedente oltre limite e punti abbastanza lontani
+         if ( vPt5ax[i].vtDir1 * vPt5ax[j].vtDir1 < COS_ANG_MAX_CORNER  &&
+              SqDist( vPt5ax[i].ptP, vPt5ax[j].ptP) > 25 * SQ_EPS_SMALL) {
+           // punto medio
+            Point3d ptMid = Media( vPt5ax[i].ptP, vPt5ax[j].ptP) ;
+            double dMid = ( vPt5ax[i].dPar + vPt5ax[j].dPar) / 2 ;
+           // se trovo proiezione, la salvo
+            Point5ax Pt5ax ; 
+            if ( ProjectPointOnSurf( ptMid, *pSurfTm, *pRefTm, dMid, Pt5ax)) {
+               vPt5ax.insert( vPt5ax.begin() + i, Pt5ax) ;
+               -- i ;
             }
          }
       }
-     // se da eliminare
-      if ( bRemove) {   
-         // dichiaro da eliminare il punto
-         vMyPt5ax[nCurr].nFlag = -1 ;
-        // avanzo con corrente e successivo
-         nCurr = nNext ;
-         ++ nNext ;
-      }
-     // altrimenti da tenere
-      else {
-        // avanzo il terzetto di uno step
-         nPrec = nCurr ;
-         nCurr = nNext ;
-         ++ nNext ;
-      }
    }
+
+  // rimozione punti in eccesso rispetto alle tolleranze
+   RemovePointsInExcess( vPt5ax, dLinTol, dMaxSegmLen, bSharpEdges) ;
+
    return true ;
 }