EgtGeomKernel :

- Aggiunte funzioni per identificazione di Parts e Shells nelle TriMesh.

ArcPntDirTgCurve.cpp

@@ -16,9 +16,9 @@
 #include "CurveBezier.h"
 #include "CurveComposite.h"
 #include "CreateCurveAux.h"
-#include "DistPointLine.h"
 #include "GeoConst.h"
 #include "/EgtDev/Include/EGkArcPntDirTgCurve.h"
+#include "/EgtDev/Include/EGkDistPointLine.h"
 #include "/EgtDev/Include/EGkDistPointCurve.h"
 #include "/EgtDev/Include/EGkArcSpecial.h"
 #include "/EgtDev/Include/EgtPointerOwner.h"

CAvToolSurfTm.cpp

@@ -14,8 +14,8 @@
 #include "stdafx.h"
 #include "CAvToolTriangle.h"
 #include "CAvToolSurfTm.h"
-#include "DistPointLine.h"
 #include "DllMain.h"
+#include "/EgtDev/Include/EGkDistPointLine.h"
 #include "/EgtDev/Include/EGnStringUtils.h"
 #include <thread>
 #include <future>
@@ -212,8 +212,9 @@ CAvToolSurfTm::TestSubPath( int nId, PNTULIST& lPntM, const Vector3d& vtDir, dou
    while ( itPntMCurr != lPntM.end()) {
      // verifico il punto
       ptCurr = itPntMCurr->first ;
-      itPntMCurr->second = MyTestPositionHG( itPntMCurr->first, vtDir) ;
-      if ( itPntMCurr->second < - EPS_SMALL)
+      double dMove = MyTestPositionHG( itPntMCurr->first, vtDir) ;
+      itPntMCurr->second = - dMove ;
+      if ( dMove < - EPS_SMALL)
          return false ;
      // se esiste il punto precedente devo verificare il medio
       if ( itPntMPrev != lPntM.end()) {

CAvToolTriangle.cpp

@@ -17,8 +17,8 @@
 #include "CAvToolTriangle.h"
 #include "IntersLineSurfStd.h"
 #include "IntersLineTria.h"
-#include "DistPointLine.h"
 #include "CDeUtility.h"
+#include "/EgtDev/Include/EGkDistPointLine.h"
 #include "/EgtDev/Include/EGkIntervals.h"
 #include "/EgtDev/Include/ENkPolynomialRoots.h"
 #include "/EgtDev/Include/EgtNumUtils.h"

CDeCylTria.cpp

@@ -13,9 +13,9 @@
 
 //--------------------------- Include ----------------------------------------
 #include "stdafx.h"
-#include "DistPointLine.h"
 #include "CDeCylTria.h"
 #include "CDeConvexTorusTria.h"
+#include "/EgtDev/Include/EGkDistPointLine.h"
 #include "/EgtDev/Include/EGkPolygon3d.h"
 
 using namespace std ;

CDeTriaTria.cpp

@@ -12,11 +12,10 @@
 //----------------------------------------------------------------------------
 
 #include "stdafx.h"
-#include "DistLineLine.h"
 #include "CDeTriaTria.h"
+#include "/EgtDev/Include/EGkDistLineLine.h"
 #include "/EgtDev/Include/EGkIntersPlanePlane.h"
 #include <array>
-#include <algorithm>
 
 using namespace std ;
 

CDeUtility.cpp

@@ -1,6 +1,6 @@
 #include "stdafx.h"
 #include "CDeUtility.h"
-#include "DistPointLine.h"
+#include "/EgtDev/Include/EGkDistPointLine.h"
 #include "/EgtDev/Include/ENkPolynomialRoots.h"
 
 using namespace std ;

CalcPocketing.cpp

@@ -759,7 +759,7 @@ GetUnclearedRegion( ICRVCOMPOPOVECTOR& vFirstOffs, ICRVCOMPOPOVECTOR& vCrvs, ICU
          PtrOwner<ISurfFlatRegion> pSrfToolRegLinki( GetSurfFlatRegionFromFatCurve( Release( pCrvLink_i), s_dRad + 5 * EPS_SMALL, false, false)) ;
          if ( ! IsNull( pSrfToolRegLinki)) {
             if ( ! pSrfTool_Links->IsValid()  ||  pSrfTool_Links->GetChunkCount() == 0)
-               pSrfTool_Links.Set( GetBasicSurfFlatRegion( pSrfToolRegLinki)) ;
+               pSrfTool_Links.Set( GetBasicSurfFlatRegion( Release( pSrfToolRegLinki))) ;
             else
                pSrfTool_Links->Add( *pSrfToolRegLinki) ;
          }
@@ -2536,7 +2536,7 @@ CalcZigZag( const ISurfFlatRegion* pSrfZigZag, ICRVCOMPOPOVECTOR& vpCrvs, double
              vAddedLinks.emplace_back( pCrvLink->Clone()) ;
 
          // pSeg2 ora diventa pSeg1 e il procedimento si itera per tutto il percorso
-          pLastSeg.Set( Release( pSeg2)) ;
+          pLastSeg.Set( pSeg2) ;
       }
      // cerco nella stessa fila o in quella successiva sezione successiva raccordabile tramite il contorno
       double dUstart = Sec.dOe ;
@@ -3164,7 +3164,7 @@ AdjustLinkDifferentY( const ICurve* pCrvSegPrec, const ICurve* pCrvSegSucc, cons
    if ( IsNull( pSfrTrap_c))
       return false ;
    if ( ! pSfrTrap->Offset( - s_dRad / 5, ICurve::OFF_CHAMFER))
-      pSfrTrap.Set( Release( pSfrTrap_c)) ;
+      pSfrTrap.Set( pSfrTrap_c) ;
 
   // se il Link interseca questa regione, allora va eliminato
    CRVCVECTOR ccClass ;

CircleCenTgCurve.cpp

@@ -16,8 +16,8 @@
 #include "CurveBezier.h"
 #include "CurveComposite.h"
 #include "CreateCurveAux.h"
-#include "DistPointLine.h"
 #include "GeoConst.h"
+#include "/EgtDev/Include/EGkDistPointLine.h"
 #include "/EgtDev/Include/EGkCircleCenTgCurve.h"
 #include "/EgtDev/Include/EGkDistPointCurve.h"
 #include "/EgtDev/Include/EgtPointerOwner.h"

CurveAux.cpp

@@ -403,9 +403,20 @@ CopyThickness( const ICurve* pSouCrv, ICurve* pDestCrv)
    return true ;
 }
 
+//----------------------------------------------------------------------------
+ICurveBezier*
+LineToBezierCurve( const ICurveLine* pCrvLine)
+{
+   PtrOwner<ICurveBezier> pCrvBezier( CreateCurveBezier()) ;
+   // rendo tutte le curve di grado 2 e razionali così posso convertire anche archi e avere tutte curve dello stesso grado e razionali
+   pCrvBezier->Init( 2, true) ;
+   pCrvBezier->FromLine( *pCrvLine) ;
+   return Release( pCrvBezier) ;
+}
+
 //----------------------------------------------------------------------------
 ICurve*
-ArcToBezierCurve( const ICurve* pCrv)
+ArcToBezierCurve( const ICurve* pCrv, bool bDeg3OrDeg2)
 {
   // verifico sia un arco
    const CurveArc* pArc = GetBasicCurveArc( pCrv) ;
@@ -413,7 +424,8 @@ ArcToBezierCurve( const ICurve* pCrv)
       return nullptr ;
 
   // se angolo al centro sotto il limite, basta una curva
-   if ( abs( pArc->GetAngCenter()) < BEZARC_ANG_CEN_MAX + EPS_ANG_SMALL) {
+   if ( ( abs( pArc->GetAngCenter()) < BEZARC_ANG_CEN_MAX + EPS_ANG_SMALL  &&  abs( pArc->GetDeltaN()) < EPS_ZERO  &&  ! bDeg3OrDeg2)  ||
+        ( abs( pArc->GetAngCenter()) < BEZARC_ANG_CEN_MAX + EPS_ANG_SMALL  &&  bDeg3OrDeg2)) {
      // creo la curva di Bezier equivalente all'arco
       PtrOwner<CurveBezier> pCrvBez( CreateBasicCurveBezier()) ;
       if ( IsNull( pCrvBez)  ||  ! pCrvBez->FromArc( *pArc))
@@ -429,6 +441,8 @@ ArcToBezierCurve( const ICurve* pCrv)
          return nullptr ;
      // inserisco nella CC le curve di Bezier equivalenti alle parti dell'arco
       int nParts = (int) ceil( abs( pArc->GetAngCenter()) / ( BEZARC_ANG_CEN_MAX + EPS_ANG_SMALL)) ;
+      if ( ! bDeg3OrDeg2  &&  abs( pArc->GetDeltaN()) > EPS_ZERO)
+         nParts *= 2 ;
       nParts = max( nParts, 2) ;
       for ( int i = 0 ; i < nParts ; ++ i) {
         // copio l'arco originale
@@ -451,6 +465,260 @@ ArcToBezierCurve( const ICurve* pCrv)
    }
 }
 
+//----------------------------------------------------------------------------
+ICurve*
+CompositeToBezierCurve( const ICurveComposite* pCC)
+{
+   // converto tutte le curve in bezier razionali di grado 2
+   PtrOwner<ICurveComposite> pCCBezier( CreateCurveComposite()) ;
+   for ( int i = 0 ; i < int( pCC->GetCurveCount()) ; ++i) {
+      PtrOwner<ICurve> pCrvNew ;
+      if ( pCC->GetCurve(i)->GetType() == CRV_ARC) {
+         const CurveArc* crArc = GetBasicCurveArc( pCC->GetCurve(i)) ;
+         ICurve* pCrvBezier = ArcToBezierCurve( crArc) ;
+         if ( pCrvBezier == nullptr)
+            return nullptr ;
+         // se la curva è di grado superiore al secondo allora devo ricondurla al secondo grado
+         pCrvNew.Set( pCrvBezier) ;
+      }
+      else if ( pCC->GetCurve(i)->GetType() == CRV_LINE) {
+         const CurveLine* crLine = GetBasicCurveLine( pCC->GetCurve(i)) ;
+         ICurve* pCrvBezier = LineToBezierCurve( crLine) ;
+         if ( pCrvBezier == nullptr)
+            return nullptr ;
+         pCrvNew.Set( pCrvBezier) ;
+      }
+      else if ( pCC->GetCurve(i)->GetType() == CRV_BEZIER ) {
+         const CurveBezier* crvBezier = GetBasicCurveBezier( pCC->GetCurve(i)) ;
+         ICurve* pCrvBezier = BezierToBasicBezierCurve( crvBezier) ;
+         if ( pCrvBezier == nullptr)
+            return nullptr ;
+         pCrvNew.Set( pCrvBezier) ;
+      }
+
+      pCCBezier->AddCurve( Release( pCrvNew)) ;
+   }
+   return Release( pCCBezier) ;
+}
+
+//----------------------------------------------------------------------------
+ICurve*
+BezierToBasicBezierCurve( const ICurveBezier* pCrvBezier)
+{ 
+   // resta da calcolare un errore sull'approssimazione oppure usare la tecnica di spezzare la curva originale in sottocurve e approssimarle con bezier cubiche
+   // per ridurre molto l'errore
+   
+   // dovrei restituire una bezier di grado 2, razionale per poter essere uniforme con le altre curve trasmorate in bezier
+   PtrOwner<ICurveBezier> pCrvNew( pCrvBezier->Clone()) ;
+   int nDeg = pCrvBezier->GetDegree() ;
+   bool bRat = pCrvBezier->IsRational() ;
+   // se la curva è già nella forma giusta la restituisco
+   if ( nDeg == 2  &&  bRat)
+      return Release( pCrvNew) ;
+   // sennò la riduco di grado fino al 2
+   PtrOwner<ICurveBezier> pBasicBezier ( CreateCurveBezier()) ;
+   pBasicBezier->Init( 2, true) ;
+   double dW = 1 ;
+   if ( nDeg == 1) {
+      Point3d ptStart = pCrvBezier->GetControlPoint( 0) ;
+      Point3d ptEnd = pCrvBezier->GetControlPoint( 1) ;
+      pBasicBezier->SetControlPoint( 0, ptStart, dW) ;
+      pBasicBezier->SetControlPoint( 1, 0.5 * (ptStart + ptEnd), dW) ;
+      pBasicBezier->SetControlPoint( 2, ptEnd, dW) ;
+      return Release( pBasicBezier) ;
+   }
+   while ( nDeg > 2 ) {
+      pCrvNew.Set( BezierDecreaseDegree( pCrvNew, 1)) ;
+      if ( IsNull( pCrvNew)  ||  ! pCrvNew->IsValid())
+         return nullptr ;
+      -- nDeg ;
+   }
+   for ( int i = 0 ; i < 3 ; ++i) {
+      Point3d ptCopy = pCrvNew->GetControlPoint( i) ;
+      if ( bRat)
+         dW = pCrvNew->GetControlWeight( i) ;
+      pBasicBezier->SetControlPoint( i, ptCopy, dW) ;
+   }
+
+   return Release( pBasicBezier) ;
+}
+
+ICurveBezier*
+BezierIncreaseDegree(const ICurveBezier* pCrvBezier)
+{
+   if ( pCrvBezier == nullptr)
+      return nullptr ;
+   // creo la versione con grado aumentato
+   PtrOwner<ICurveBezier> pNewBezier( CreateCurveBezier()) ;
+   int nDeg = pCrvBezier->GetDegree() + 1;
+   bool bRat = pCrvBezier->IsRational();
+   pNewBezier->Init( nDeg , bRat) ;
+   // prev e curr sono riferiti alla curva di partenza
+   // salvo il primo punto
+   Point3d ptCtrlPrev = pCrvBezier->GetControlPoint( 0) ;
+   double dWprev = 1 ;
+   if ( bRat ) {
+      dWprev = pCrvBezier->GetControlWeight( 0) ;
+      pNewBezier->SetControlPoint( 0, ptCtrlPrev, dWprev) ;
+   }
+   else
+      pNewBezier->SetControlPoint( 0, ptCtrlPrev) ;
+   // ciclo sui punti di controllo intermedi per calcolare quelli nuovi
+   Point3d ptCtrlCurr ;
+   double dWcurr ;
+   for ( double i = 1 ; i < nDeg ; ++i) {
+      ptCtrlCurr = pCrvBezier->GetControlPoint( int( i)) ;
+      Point3d ptNew = ( i / nDeg) * ptCtrlPrev + ( 1 - i / ( nDeg)) * ptCtrlCurr ;
+      if ( bRat) {
+         dWcurr = pCrvBezier->GetControlWeight( int( i)) ;
+         double dWnew = ( i / nDeg) * dWprev + ( 1 - i / ( nDeg)) * dWcurr ;
+         pNewBezier->SetControlPoint( int( i), ptNew, dWnew) ;
+         dWprev = dWnew ;
+      }
+      else
+         pNewBezier->SetControlPoint( int( i), ptNew) ;
+      ptCtrlPrev = ptCtrlCurr ;
+   }
+   // salvo l'ultimo punto
+   ptCtrlCurr = pCrvBezier->GetControlPoint( nDeg - 1) ;
+   if ( bRat ) {
+      dWcurr = pCrvBezier->GetControlWeight( nDeg - 1) ;
+      pNewBezier->SetControlPoint( nDeg, ptCtrlCurr, dWcurr) ;
+   }
+   else
+      pNewBezier->SetControlPoint( nDeg, ptCtrlCurr) ;
+   return Release( pNewBezier) ;
+}
+
+ICurveBezier*
+BezierDecreaseDegree(const ICurveBezier* pCrvBezier, double dTol)
+{
+   if ( pCrvBezier == nullptr)
+      return nullptr ;
+   PtrOwner<ICurveBezier> pNewBezier( CreateCurveBezier()) ;
+   int nDeg = pCrvBezier->GetDegree() - 1 ;
+   if ( nDeg == 0)
+      return nullptr ;
+   bool bRat = pCrvBezier->IsRational() ;
+   pNewBezier->Init( nDeg, bRat) ;
+   // prev e curr sono riferiti alla nuova curva
+   //salvo il primo punto
+   Point3d ptCtrlPrev = pCrvBezier->GetControlPoint( 0) ;
+   double dWprev = 1 ;
+   if ( bRat ) {
+      dWprev = pCrvBezier->GetControlWeight( 0) ;
+      pNewBezier->SetControlPoint( 0, ptCtrlPrev, dWprev) ;
+   }
+   else
+      pNewBezier->SetControlPoint( 0, ptCtrlPrev) ;
+   // ciclo sui punti intermedi per trovare i nuovi punti di controllo  // equazioni da NURBS book
+   Point3d ptCtrlCurr ;
+   double dWcurr ;
+   int r = int (nDeg / 2) ;
+   // prima metà
+   for ( double i = 1 ; i < r ; ++i) {
+      double dAlpha = i / ( nDeg + 1) ;
+      // old è riferito alla curva originale
+      Point3d ptOld = pCrvBezier->GetControlPoint( int( i)) ;
+      ptCtrlCurr = ( ptOld + (- dAlpha * ptCtrlPrev)) / ( 1 - dAlpha) ;
+      if ( bRat) {
+         double dWold = pCrvBezier->GetControlWeight( int( i)) ;
+         dWcurr = ( dWold + (- dAlpha * dWprev)) / ( 1 - dAlpha) ;
+         pNewBezier->SetControlPoint( int( i), ptCtrlCurr, dWcurr) ;
+         dWprev = dWcurr ;
+      }
+      else
+         pNewBezier->SetControlPoint( int( i), ptCtrlCurr) ;
+      ptCtrlPrev = ptCtrlCurr ;
+   }
+   // risolvo il punto di mezzo per il caso nDeg pari
+   if ( ( nDeg + 1) % 2 == 0) {  // pari
+      double dAlpha = r / ( nDeg + 1.) ;
+      Point3d ptOld = pCrvBezier->GetControlPoint( r) ;
+      ptCtrlCurr = ( ptOld + (- dAlpha * ptCtrlPrev)) / ( 1 - dAlpha) ;
+      if ( bRat) {
+         double dWold = pCrvBezier->GetControlWeight( r) ;
+         dWcurr = ( dWold + (- dAlpha * dWprev)) / ( 1 - dAlpha) ;
+         pNewBezier->SetControlPoint( r, ptCtrlCurr, dWcurr) ;
+         dWprev = dWcurr ;
+      }
+      else
+         pNewBezier->SetControlPoint( r, ptCtrlCurr) ;
+   }
+   // salvo l'ultimo punto
+   ptCtrlCurr = pCrvBezier->GetControlPoint( nDeg + 1) ;
+   if ( bRat ) {
+      dWcurr = pCrvBezier->GetControlWeight( nDeg + 1) ;
+      pNewBezier->SetControlPoint( nDeg, ptCtrlCurr, dWcurr) ;
+   }
+   else
+      pNewBezier->SetControlPoint( nDeg, ptCtrlCurr) ;
+   ptCtrlPrev = ptCtrlCurr ;
+   if ( bRat)
+      dWprev = dWcurr ;
+   // seconda metà
+   for ( double i = nDeg - 1 ; i > r ; --i) {
+      double dAlpha = ( i + 1) / ( nDeg + 1) ;
+      Point3d ptOld = pCrvBezier->GetControlPoint( int( i) + 1) ;
+      ptCtrlCurr = ( ptOld + ( - ( 1 - dAlpha) * ptCtrlPrev)) / dAlpha ;
+      if ( bRat) {
+         double dWold = pCrvBezier->GetControlWeight( int( i) + 1) ;
+         dWcurr = ( dWold - ( ( 1 - dAlpha) * dWprev)) / dAlpha ;
+         pNewBezier->SetControlPoint( int( i), ptCtrlCurr, dWcurr) ;
+         dWprev = dWcurr ;
+      }
+      else
+         pNewBezier->SetControlPoint( int( i), ptCtrlCurr) ;
+      ptCtrlPrev = ptCtrlCurr ;
+   }
+   // risolvo il punto di mezzo per il caso nDeg dispari
+   // calcolo anche l'errore di approssimazione
+   double dErr = 0 ;
+   if ( (nDeg + 1) % 2 != 0) { // dispari
+      // puntdo di sinistra
+      double dAlpha = r / ( nDeg + 1.) ;
+      ptCtrlPrev = pNewBezier->GetControlPoint( r - 1) ;
+      Point3d ptOld = pCrvBezier->GetControlPoint( r) ;
+      Point3d ptL = ( ptOld + ( - dAlpha * ptCtrlPrev)) / ( 1 - dAlpha) ;
+      double dWL = 1 ;
+      if ( bRat) {
+         dWprev = pNewBezier->GetControlWeight( r - 1) ;
+         double dWold =pCrvBezier->GetControlWeight( r) ;
+         dWL = ( dWold - ( dAlpha * dWprev)) / ( 1 - dAlpha) ;
+      }
+      // punto di destra
+      dAlpha = ( r + 1.) / ( nDeg + 1.) ;
+      ptCtrlPrev = pNewBezier->GetControlPoint( r + 1) ;
+      ptOld = pCrvBezier->GetControlPoint( r + 1) ;
+      double dWR = 1 ;
+      if ( bRat) {
+         dWprev = pNewBezier->GetControlWeight( r + 1) ;
+         double dWold =pCrvBezier->GetControlWeight( r + 1) ;
+         dWR = ( dWold - ( ( 1 - dAlpha) * dWprev)) / dAlpha ;
+      }
+      Point3d ptR = ( ptOld + ( - ( 1 - dAlpha) * ptCtrlPrev)) / dAlpha ;
+      // calcolo il punto di mezzo e lo aggiungo
+      Point3d ptNew = ( ptL + ptR) / 2 ;
+      double dWnew = ( dWL + dWR) / 2 ;
+      if ( bRat)
+         pNewBezier->SetControlPoint( r, ptNew, dWnew) ;
+      else
+         pNewBezier->SetControlPoint( r, ptNew) ;
+      dErr = Dist( ptL, ptR) ;
+   }
+   else {
+      ptCtrlCurr = pNewBezier->GetControlPoint( r + 1) ;
+      ptCtrlPrev = pNewBezier->GetControlPoint( r) ;
+      Point3d ptOld = pCrvBezier->GetControlPoint( r + 1) ;
+      dErr = Dist( ptOld, 0.5 * ( ptCtrlPrev + ptCtrlCurr)) ;
+   }
+   //// se l'approssimazione dà un errore troppo alto allora annullo tutto  // da controllare
+   //if ( dErr > dTol)
+   //   return nullptr ;
+
+   return Release( pNewBezier) ;
+}
+
 //----------------------------------------------------------------------------
 ICurve*
 CurveToNoArcsCurve( const ICurve* pCrv)

CurveBezier.cpp

@@ -19,7 +19,6 @@
 #include "BiArcs.h"
 #include "GeoConst.h"
 #include "PolygonPlane.h"
-#include "DistPointLine.h"
 #include "GeoObjFactory.h"
 #include "NgeWriter.h"
 #include "NgeReader.h"
@@ -27,12 +26,16 @@
 #include "Bernstein.h"
 #include "deCasteljau.h"
 #include "Voronoi.h"
+#include "IntersLineLine.h"
 #include "/EgtDev/Include/EGkCurveArc.h"
+#include "/EgtDev/Include/EGkDistPointLine.h"
 #include "/EgtDev/Include/EGkStringUtils3d.h"
 #include "/EgtDev/Include/EGkUiUnits.h"
 #include "/EgtDev/Include/ENkPolynomial.h"
 #include "/EgtDev/Include/EgtNumUtils.h"
 #include "/EgtDev/Include/EgtPointerOwner.h"
+#include "/EgtDev/Include/EGkGeoPoint3d.h"
+#include "/EgtDev/Include/EGkCurveLine.h"
 #include <array>
 
 using namespace std ;
@@ -195,6 +198,35 @@ CurveBezier::FromArc( const ICurveArc& crArc)
       SetControlPoint( 1, ptNew1, dW) ;
       SetControlPoint( 2, ptNew2, dW) ;
       SetControlPoint( 3, ptEnd, 1) ;
+
+   //// anziché usare una bezier cubica approssimo le eliche con più bezier quadratiche razionali.
+   //// più è grande l'angolo al centro dell'arco e maggiore sarà l'errore di approssimazione
+   //
+   //  // quadratica razionale
+   //  // peso del punto di controllo intermedio
+   //   double dW = dCosAhalf ;
+   //  // calcolo dei punti di controllo
+   //   Point3d  ptStart ;
+   //   crArc.GetStartPoint( ptStart) ;
+   //   Point3d  ptEnd ;
+   //   crArc.GetEndPoint( ptEnd) ;
+   //  // Point3d ptMed = Media( ptStart, ptEnd, 0.5) ;
+   //  //Vector3d vtDir = ptMed - crArc.GetCenter() ;
+   //  // Point3d ptNew = crArc.GetCenter() + vtDir / ( dCosAhalf * dCosAhalf) ;
+   //   Vector3d vtStart ; crArc.GetStartDir( vtStart) ;
+   //   Vector3d vtEnd ; crArc.GetEndDir( vtEnd) ;
+   //   PtrOwner<CurveLine> pCrvLine1( CreateBasicCurveLine()) ;
+   //   pCrvLine1->SetPVL( ptStart, vtStart, 10) ;
+   //   PtrOwner<CurveLine> pCrvLine2( CreateBasicCurveLine()) ;
+   //   pCrvLine2->SetPVL( ptEnd, vtEnd, 10) ;
+   //   IntersLineLine ill( *pCrvLine1, *pCrvLine2, false) ;
+   //   IntCrvCrvInfo iccInfo ; ill.GetIntCrvCrvInfo( iccInfo) ;
+   //   Point3d ptNew = 0.5 * (iccInfo.IciA->ptI + iccInfo.IciB->ptI) ;
+   //  // inserimento nella curva dei punti di controllo con i pesi
+   //   Init( 2, true) ;
+   //   SetControlPoint( 0, ptStart, 1) ;
+   //   SetControlPoint( 1, ptNew, dW) ;
+   //   SetControlPoint( 2, ptEnd, 1) ;
    }
 
   // copio estrusione e spessore
@@ -204,6 +236,30 @@ CurveBezier::FromArc( const ICurveArc& crArc)
    return true ;
 }
 
+//----------------------------------------------------------------------------
+bool
+CurveBezier::FromLine( const ICurveLine& crLine)
+{
+   if ( ! crLine.IsValid())
+      return false ;
+   double dWeight = 1 ;
+   int nCount = 0 ;
+   Point3d ptStart ; crLine.GetStartPoint( ptStart) ;
+   SetControlPoint( nCount, ptStart, dWeight) ;
+   ++nCount ;
+   double dPart = 1. / m_nDeg ;
+   for ( int i = 1 ; i < m_nDeg ; ++i) {
+      double dU = i * dPart ;
+      Point3d ptMid ; crLine.GetPointD1D2( dU, ICurve::FROM_MINUS, ptMid) ;
+      SetControlPoint( nCount, ptMid, dWeight) ;
+      ++nCount ;
+   }
+   Point3d ptEnd ; crLine.GetEndPoint( ptEnd) ;
+   SetControlPoint( nCount, ptEnd, dWeight) ;
+   ++nCount ;
+   return true ;
+}
+
 //----------------------------------------------------------------------------
 bool
 CurveBezier::IsAPoint( void) const

CurveBezier.h

@@ -138,6 +138,7 @@ class CurveBezier : public ICurveBezier, public IGeoObjRW
       bool SetControlPoint( int nInd, const Point3d& ptCtrl) override ;
       bool SetControlPoint( int nInd, const Point3d& ptCtrl, double dW) override ;
       bool FromArc( const ICurveArc& crArc) override ;
+      bool FromLine( const ICurveLine& crLine) override ;
       int  GetDegree( void) const override
                 { return m_nDeg ; }
       bool IsRational( void) const override
@@ -190,7 +191,7 @@ class CurveBezier : public ICurveBezier, public IGeoObjRW
 
    private :
       enum Status { ERR = 0, OK = 1, TO_VERIFY = 2} ;
-      static const int MAXDEG = 11 ;
+      static const int MAXDEG = 21 ;
 
    private :
       ObjGraphicsMgr m_OGrMgr ;         // gestore grafica dell'oggetto

CurveByApprox.cpp

@@ -16,8 +16,8 @@
 #include "CurveComposite.h"
 #include "CalcDerivate.h"
 #include "BiArcs.h"
-#include "DistPointLine.h"
 #include "RemoveCurveDefects.h"
+#include "/EgtDev/Include/EGkDistPointLine.h"
 #include "/EgtDev/Include/EGkCurveByApprox.h"
 #include "/EgtDev/Include/EGkPolyLine.h"
 #include "/EgtDev/Include/EGkPolyArc.h"

CurveComposite.cpp

@@ -14,7 +14,6 @@
 //--------------------------- Include ----------------------------------------
 #include "stdafx.h"
 #include "CurveComposite.h"
-#include "DistPointLine.h"
 #include "DistPointCrvComposite.h"
 #include "CurveLine.h"
 #include "CurveArc.h"
@@ -27,6 +26,7 @@
 #include "NgeWriter.h"
 #include "NgeReader.h"
 #include "Voronoi.h"
+#include "/EgtDev/Include/EGkDistPointLine.h"
 #include "/EgtDev/Include/EGkCurveByApprox.h"
 #include "/EgtDev/Include/EGkArcSpecial.h"
 #include "/EgtDev/Include/EGkSfrCreate.h"
@@ -745,7 +745,7 @@ CurveComposite::Load( NgeReader& ngeIn)
       ICurve* pCrv = ::GetCurve( pGeoO) ;
       bOk = bOk && ( pCrv != nullptr  &&  pCrv->IsSimple()) ;
      // aggiungo questa curva (sicuramente semplice)
-      bOk = bOk && AddSimpleCurve( pCrv) ;
+      bOk = bOk && AddSimpleCurve( pCrv, true, 10 * EPS_SMALL) ;
      // se errore 
      if ( ! bOk)
          return false ;
@@ -878,10 +878,13 @@ CurveComposite::TestClosure( void)
    Point3d ptEnd ; m_CrvSmplS.back()->GetEndPoint( ptEnd) ;
    // se distanza superiore al limite ridotto forzo i punti a coincidere
    if ( ! AreSamePointEpsilon( ptStart, ptEnd, EPS_CONNECT)) {
+     // se un solo arco
+      if ( m_CrvSmplS.size() == 1  &&  m_CrvSmplS.front()->GetType() == CRV_ARC)
+         return GetBasicCurveArc( m_CrvSmplS.front())->ChangeAngCenter( ANG_FULL) ;
+     // caso generale
       Point3d ptM = Media( ptStart, ptEnd) ;
-      if ( ! m_CrvSmplS.front()->ModifyStart( ptM)  ||
-           ! m_CrvSmplS.back()->ModifyEnd( ptM))
-         return false ;
+      return ( m_CrvSmplS.front()->ModifyStart( ptM)  &&
+               m_CrvSmplS.back()->ModifyEnd( ptM)) ;
    }
    return true ;
 }
@@ -3222,61 +3225,64 @@ MergeTwoCurves( ICurve* pCrvP, ICurve* pCrvC, double& dCurrLinTol, double dCosAn
      // verifico di non superare l'angolo giro al centro
       if ( abs( pArcP->GetAngCenter() + pArcC->GetAngCenter()) > ANG_FULL + EPS_ANG_SMALL)
          return 0 ;
-     // se archi piatti
-      if ( pArcP->IsPlane()  &&  pArcC->IsPlane()) {
-        // se calcolo nuovo arco ok, procedo con l'unione
-         Point3d ptP1 ;
-         pArcP->GetStartPoint( ptP1) ;
-         Point3d ptP2 ;
-         pArcP->GetEndPoint( ptP2) ;
-         Point3d ptP3 ;
-         pArcC->GetEndPoint( ptP3) ;
-         // verifico se circonferenza completa
-         bool bCirc = ( AreSamePointApprox( ptP1, ptP3)) ;
-         if ( bCirc)
-            pArcC->GetMidPoint( ptP3) ;
-         CurveArc NewArc ;
-         if ( NewArc.Set3P( ptP1, ptP2, ptP3, bCirc)) {
-           // verifico normale al piano dell'arco
-            if ( NewArc.GetNormVersor() * pArcC->GetNormVersor() < 0)
-               NewArc.InvertN() ;
-           // se curve originali con la stessa proprietà, la riporto
-            if ( nTpr0P == nTpr0C)
-               NewArc.SetTempProp( nTpr0C, 0) ;
-            if ( nTpr1P == nTpr1C)
-               NewArc.SetTempProp( nTpr1C, 1) ;
-           // aggiorno l'arco corrente e torno flag modifica
-            *pArcC = NewArc ;
-            return -1 ;
-         }
-         else
+     // verifico se archi piatti
+      bool bPlaneArcs = pArcP->IsPlane()  &&  pArcC->IsPlane() ;
+     // se archi non piatti verifico coincidenza pendenza sulla normale
+      if ( ! bPlaneArcs) {        
+         double dN = pArcP->GetNormVersor() * pArcC->GetNormVersor() ;
+         if ( abs(( pArcC->GetDeltaN() * pArcP->GetAngCenter() - dN * pArcP->GetDeltaN() * pArcC->GetAngCenter()) /
+            ( pArcP->GetAngCenter() + pArcC->GetAngCenter())) > dCurrLinTol)
             return 0 ;
       }
-     // verifico coincidenza pendenza sulla normale
-      double dN = pArcP->GetNormVersor() * pArcC->GetNormVersor() ;
-      if ( abs(( pArcC->GetDeltaN() * pArcP->GetAngCenter() - dN * pArcP->GetDeltaN() * pArcC->GetAngCenter()) /
-                ( pArcP->GetAngCenter() + pArcC->GetAngCenter())) < dCurrLinTol) {
-        // se calcolo nuovo arco ok, procedo con l'unione
-         Point3d ptP1 ;
-         pArcP->GetStartPoint( ptP1) ;
-         Vector3d vtDir1 ;
-         pArcP->GetStartDir( vtDir1) ;
-         Point3d ptP3 ;
-         pArcC->GetEndPoint( ptP3) ;
-         CurveArc NewArc ;
-         if ( NewArc.Set2PVN( ptP1, ptP3, vtDir1, pArcC->GetNormVersor())) {
-           // se curve originali con la stessa proprietà, la riporto
-            if ( nTpr0P == nTpr0C)
-               NewArc.SetTempProp( nTpr0C, 0) ;
-            if ( nTpr1P == nTpr1C)
-               NewArc.SetTempProp( nTpr1C, 1) ;
-           // aggiorno l'arco corrente e torno flag modifica
-            *pArcC = NewArc ;
-            return -1 ;
-         }
-         else
-            return 0 ;
+
+     // se calcolo nuovo arco ok, procedo con l'unione
+      Point3d ptP1 ;
+      pArcP->GetStartPoint( ptP1) ;
+      Point3d ptP2 ;
+      pArcP->GetEndPoint( ptP2) ;
+      Point3d ptP3 ;
+      pArcC->GetEndPoint( ptP3) ;
+
+     // se archi non piani costruisco arco sul piano definito dalla normale e dal punto di partenza del primo arco
+      Frame3d frRef ;
+      if ( ! frRef.Set( ptP1, pArcP->GetNormVersor()))
+         return 0 ;
+      if ( ! bPlaneArcs) {
+         ptP1.Scale( frRef, 1, 1, 0) ;
+         ptP2.Scale( frRef, 1, 1, 0) ;
+         ptP3.Scale( frRef, 1, 1, 0) ;      
       }
+
+     // verifico se circonferenza completa
+      bool bCirc = ( AreSamePointApprox( ptP1, ptP3)) ;
+      if ( bCirc) {
+         pArcC->GetMidPoint( ptP3) ;
+         if ( ! bPlaneArcs)
+            ptP3.Scale( frRef, 1, 1, 0) ;
+      }
+         
+      CurveArc NewArc ;
+      if ( NewArc.Set3P( ptP1, ptP2, ptP3, bCirc)) {
+        // verifico normale al piano dell'arco
+         if ( NewArc.GetNormVersor() * pArcC->GetNormVersor() < 0)
+            NewArc.InvertN() ;
+        // se archi non piani ripristino il deltaN
+         if ( ! bPlaneArcs) {
+            double dDeltaN1 = pArcP->GetDeltaN() ;
+            double dDeltaN2 = pArcC->GetDeltaN() ;
+            NewArc.ChangeDeltaN( dDeltaN1 + dDeltaN2) ;
+         }
+        // se curve originali con la stessa proprietà, la riporto
+         if ( nTpr0P == nTpr0C)
+            NewArc.SetTempProp( nTpr0C, 0) ;
+         if ( nTpr1P == nTpr1C)
+            NewArc.SetTempProp( nTpr1C, 1) ;
+        // aggiorno l'arco corrente e torno flag modifica
+         *pArcC = NewArc ;
+         return -1 ;
+      }
+      else
+         return 0 ;
    }
 
   // nessuna fusione

CurveLine.cpp

@@ -14,11 +14,11 @@
 //--------------------------- Include ----------------------------------------
 #include "stdafx.h"
 #include "CurveLine.h"
-#include "DistPointLine.h"
 #include "GeoObjFactory.h"
 #include "NgeWriter.h"
 #include "NgeReader.h"
 #include "Voronoi.h"
+#include "/EgtDev/Include/EGkDistPointLine.h"
 #include "/EgtDev/Include/EGkStringUtils3d.h"
 #include "/EgtDev/Include/EgtNumUtils.h"
 #include "/EgtDev/Include/EgtPointerOwner.h"

DistLineLine.cpp

@@ -1,7 +1,7 @@
 //----------------------------------------------------------------------------
-//  EgalTech 2020-2022
+//  EgalTech 2020-2024
 //----------------------------------------------------------------------------
-// File : DistLineLine.h      Data : 12.08.22         Versione : 2.4h1
+// File : DistLineLine.cpp      Data : 10.05.24         Versione : 2.6e3
 // Contenuto : Implementazione della classe distanza fra elementi lineari.
 //
 //
@@ -12,11 +12,10 @@
 //----------------------------------------------------------------------------
 
 #include "stdafx.h"
-#include "DistLineLine.h"
+#include "/EgtDev/Include/EGkDistLineLine.h"
 #include "/EgtDev/Include/EgtNumUtils.h"
 #include "/EgtDev/Include/EGkGeoCollection.h"
 #include "/EgtDev/Include/EGkGeoConst.h"
-#include <algorithm>
 
 using namespace std ;
 

DistLineLine.h

@@ -1,53 +0,0 @@
-//----------------------------------------------------------------------------
-//  EgalTech 2020-2020
-//----------------------------------------------------------------------------
-// File : DistLineLine.h      Data : 06.11.20         Versione : 2.2k1
-// Contenuto : Dichiarazione della classe distanza fra elementi lineari.
-//
-//
-//
-// Modifiche : 06.11.20 LM Creazione modulo.
-//
-//
-//----------------------------------------------------------------------------
-
-#pragma once
-
-#include "/EgtDev/Include/EGkVector3d.h"
-#include "/EgtDev/Include/EGkPoint3d.h"
-
-//----------------------------------------------------------------------------
-class DistLineLine
-{
-   public :
-      DistLineLine( const Point3d& ptSt1, const Point3d& ptEn1,
-                    const Point3d& ptSt2, const Point3d& ptEn2,
-                    bool bIsSegment1 = true, bool bIsSegment2 = true) ;
-      DistLineLine( const Point3d& ptSt1, const Vector3d& vtD1, double dLen1,
-                    const Point3d& ptSt2, const Vector3d& vtD2, double dLen2,
-                    bool bIsSegment1 = true, bool bIsSegment2 = true) ;
-
-   public :
-      bool GetSqDist( double& dSqDist) const ;
-      bool GetDist( double& dDist) const ;
-      bool IsEpsilon( double dTol) const
-              { double dSqDist ; return ( GetSqDist( dSqDist)  &&  ( dSqDist < SQ_EPS_ZERO  ||  dSqDist < dTol * dTol)) ; }
-      bool IsSmall( void) const
-              { return IsEpsilon( EPS_SMALL) ; }
-      bool IsZero( void) const
-              { return IsEpsilon( EPS_ZERO) ;  }
-      bool GetMinDistPoints( Point3d& ptMinDist1, Point3d& ptMinDist2) const ;
-      bool GetPositionsAtMinDistPoints( double& dPos1, double& dPos2) const ;
-
-   private :
-      void Calculate( const Point3d& ptSt1, const Vector3d& vtD1, double dLen1,
-                      const Point3d& ptSt2, const Vector3d& vtD2, double dLen2,
-                      bool bIsSegment1, bool bIsSegment2) ;
-   private:
-      double  m_dSqDist ;
-      mutable double  m_dDist ;
-      double  m_dPos1 ;
-      double  m_dPos2 ;
-      Point3d m_ptMinDist1 ;
-      Point3d m_ptMinDist2 ;
-} ;

DistPointCrvAux.cpp

@@ -14,9 +14,9 @@
 //--------------------------- Include ----------------------------------------
 #include "stdafx.h"
 #include "DllMain.h"
-#include "DistPointCrvAux.h"
-#include "DistPointLine.h"
 #include "GeoConst.h"
+#include "DistPointCrvAux.h"
+#include "/EgtDev/Include/EGkDistPointLine.h"
 
 
 //----------------------------------------------------------------------------

DistPointCurve.cpp

@@ -13,10 +13,10 @@
 
 //--------------------------- Include ----------------------------------------
 #include "stdafx.h"
-#include "DistPointLine.h"
 #include "DistPointArc.h"
 #include "DistPointCrvBezier.h"
 #include "DistPointCrvComposite.h"
+#include "/EgtDev/Include/EGkDistPointLine.h"
 #include "/EgtDev/Include/EGkDistPointCurve.h"
 
 

DistPointLine.cpp

@@ -1,19 +1,19 @@
 //----------------------------------------------------------------------------
-//  EgalTech 2013-2013
+//  EgalTech 2013-2024
 //----------------------------------------------------------------------------
-// File : DistPointLine.cpp       Data : 17.12.13          Versione : 1.4l1
+// File : DistPointLine.cpp       Data : 20.05.24          Versione : 2.6e5
 // Contenuto : Implementazione della classe distanza punto da linea/segmento.
 //
 //
 //
 // Modifiche : 17.12.13 DS Creazione modulo.
-//
+//             20.05.24 DS Reso pubblico in Include.
 //
 //----------------------------------------------------------------------------
 
 //--------------------------- Include ----------------------------------------
 #include "stdafx.h"
-#include "DistPointLine.h"
+#include "/EgtDev/Include/EGkDistPointLine.h"
 
 
 //----------------------------------------------------------------------------

DistPointLine.h

@@ -1,58 +0,0 @@
-//----------------------------------------------------------------------------
-//  EgalTech 2013-2014
-//----------------------------------------------------------------------------
-// File : DistPointLine.h      Data : 02.01.14          Versione : 1.5a1
-// Contenuto : Dichiarazione della classe distanza punto da linea/segmento.
-//
-//
-//
-// Modifiche : 30.12.12 DS Creazione modulo.
-//
-//
-//----------------------------------------------------------------------------
-
-#pragma once
-
-#include "/EgtDev/Include/EGkPoint3d.h"
-#include "/EgtDev/Include/EGkCurveLine.h"
-
-
-//-----------------------------------------------------------------------------
-class DistPointLine
-{
-   friend class DistPointCurve ;
-
-   public :
-      DistPointLine( const Point3d& ptP,
-                     const ICurveLine& crvLine, bool bIsSegment = true) ;
-      DistPointLine( const Point3d& ptP,
-                     const Point3d& ptIni, const Point3d& ptFin, bool bIsSegment = true) ;
-      DistPointLine( const Point3d& ptP,
-                     const Point3d& ptIni, const Vector3d& vtDir, double dLen, bool bIsSegment = true) ;
-
-   public :
-      bool GetSqDist( double& dSqDist) const ;
-      bool GetDist( double& dDist) const ;
-      bool IsEpsilon( double dTol) const
-              { double dSqDist ; return ( GetSqDist( dSqDist) && ( dSqDist < SQ_EPS_ZERO || dSqDist < dTol * dTol)) ; }
-      bool IsSmall( void) const
-              { return IsEpsilon( EPS_SMALL) ; }
-      bool IsZero( void) const
-              { return IsEpsilon( EPS_ZERO) ; }
-      int  GetNbrMinDist( void) const
-              { return (( m_dSqDist < 0) ? 0 : 1) ; }
-      bool GetMinDistPoint( Point3d& ptMinDist) const ;
-      bool GetParamAtMinDistPoint( double& dParam) const ;
-
-   private :
-      DistPointLine( void) ;
-      void Calculate( const Point3d& ptP,
-                      const Point3d& ptIni, const Vector3d& vtDir, double dLen, bool bIsSegment) ;
-
-   private :
-      double  m_dSqDist ;
-      mutable double  m_dDist ;
-      double  m_dParam ;
-      Point3d m_ptMinDist ;
-} ;
-

DistPointSurfTm.cpp

@@ -15,6 +15,7 @@
 #include "SurfTriMesh.h"
 #include "/EgtDev/Include/EGkDistPointTria.h"
 #include "/EgtDev/Include/EGkDistPointSurfTm.h"
+#include "/EgtDev/Include/EGkIntersLineTria.h"
 
 using namespace std ;
 
@@ -93,6 +94,8 @@ DistPointSurfTm::Calculate( const Point3d& ptP, const ISurfTriMesh& tmSurf)
 {
   // Inizializzo distanza non calcolata
    m_dDist = - 1. ;
+  // Controllo se la superficie è chiusa
+   m_bIsSurfClosed = tmSurf.IsClosed() ;
 
   // Lavoro con l'oggetto superficie trimesh di base
    const SurfTriMesh* pStm = GetBasicSurfTriMesh( &tmSurf) ;
@@ -121,6 +124,9 @@ DistPointSurfTm::Calculate( const Point3d& ptP, const ISurfTriMesh& tmSurf)
   // Finché non si verifica la condizione di terminazione ingrandisco il box.
    pStm->ResetTempInts() ;
    bool bContinue = true ;
+
+  // creazione del vettore dei triangoli più vicini a ptP
+   vector<pair<int, Triangle3d>> vTria ; // <indice triangolo, Triangolo>
    while ( bContinue) {
      // Calcolo il box differenza con il precedente per non esplorare parti già considerate
       BOXVECTOR vBox ;
@@ -134,8 +140,8 @@ DistPointSurfTm::Calculate( const Point3d& ptP, const ISurfTriMesh& tmSurf)
             continue ;
         // ricerca sui triangoli nel box
          bCollide = true ;
-         INTVECTOR vnIds ; 
-         if ( pStm->GetAllTriaOverlapBox( b3Int, vnIds)) { 
+         INTVECTOR vnIds ;
+         if ( pStm->GetAllTriaOverlapBox( b3Int, vnIds)) {
            // Ciclo sui triangoli del sotto-box corrente
             for ( auto nT : vnIds) {
                int nTriaTemp ;
@@ -144,18 +150,24 @@ DistPointSurfTm::Calculate( const Point3d& ptP, const ISurfTriMesh& tmSurf)
                   pStm->SetTempInt( nT, 1) ;
                   DistPointTriangle distPT( ptP, trCurTria) ;
                   double dCurSqDist ;
-                 // Se la distanza del triangolo č valida e minore di quella attuale aggiorno 
-                  if ( distPT.GetSqDist( dCurSqDist)  &&  dCurSqDist < dMinSqDist) {
-                     dMinSqDist = dCurSqDist ;
-                     nMinDistTriaIndex = nT ;
-                     distPT.GetMinDistPoint( ptMinDistPoint) ;
+                 // Se la distanza del triangolo è valida e minore di quella attuale aggiorno
+                  if ( distPT.GetSqDist( dCurSqDist)) {
+                     if ( abs( dCurSqDist - dMinSqDist) < EPS_SMALL) // se distanze uguali...
+                        vTria.emplace_back( make_pair( nT, trCurTria)) ; // aggiungo il triangolo
+                     else if ( dCurSqDist < dMinSqDist) { // se minore...
+                        vTria.clear() ; // pulisco il vettore
+                        dMinSqDist = dCurSqDist ;
+                        nMinDistTriaIndex = nT ;
+                        distPT.GetMinDistPoint( ptMinDistPoint) ;
+                        vTria.emplace_back( make_pair( nT, trCurTria)) ; // aggiungo il triangolo
+                     }
                   }
                }
             }
          }
       }
      // Se si verifica la condizione di terminazione arresto il ciclo altrimenti aggiorno i box 
-      if ( ! bCollide  ||  dMinSqDist < EPS_SMALL * EPS_SMALL)
+      if ( ! bCollide  ||  dMinSqDist < SQ_EPS_SMALL)
          bContinue = false ;
       else {
         boxPPrev = boxP ;
@@ -163,15 +175,70 @@ DistPointSurfTm::Calculate( const Point3d& ptP, const ISurfTriMesh& tmSurf)
       }
    }
 
-   if ( nMinDistTriaIndex != SVT_NULL) {
-      m_dDist = sqrt( max( dMinSqDist, 0.)) ;
+  // se non ho trovato nessun triangolo, esco
+   if ( nMinDistTriaIndex == SVT_NULL)
+      return ;
+
+  // salvo la distanza minima
+   m_dDist = sqrt( max( dMinSqDist, 0.)) ;
+  // salvo il punto a distanza minima
+   m_ptMinDistPoint = ptMinDistPoint ;
+  // se il punto è sulla TriMesh...
+   if ( m_dDist < EPS_SMALL) {
       m_nMinDistTriaIndex = nMinDistTriaIndex ;
-      m_ptMinDistPoint = ptMinDistPoint ;
-      Triangle3d trMinDistTria ;
-      pStm->GetTriangle( m_nMinDistTriaIndex, trMinDistTria) ;
-      trMinDistTria.Validate() ;
-      m_bIsInside = ( ( ptP - m_ptMinDistPoint) * trMinDistTria.GetN() < - EPS_SMALL)  &&  pStm->IsClosed() ;
+      m_bIsInside = false ;
+      return ;
    }
+  // se ho un solo triangolo, allora deduco le informazioni da lui
+   else if ( int( vTria.size()) == 1) {
+      m_nMinDistTriaIndex = vTria.back().first ;
+      m_bIsInside = ( ( ptP - m_ptMinDistPoint) * vTria.back().second.GetN() < - EPS_SMALL) ;
+      return ;
+   }
+  
+  // controllo se tutti i triangoli a minima distanza forniscono la stessa informazione
+  // ( il punto potrebbe essere esterno a tutti, interno a tutti o indefinito )
+   bool bInside = false ;
+   bool bOutside = false ;
+   for ( int i = 0 ; i < int( vTria.size()) ; ++ i) { // scorro i triangoli a minima distanza
+      if ( ( ptP - vTria[i].second.GetP( 0)) * vTria[i].second.GetN() < - EPS_SMALL)
+         bInside = true ;
+      else
+         bOutside = true ;
+   }
+
+  // inizializzo le variabili membro
+   m_nMinDistTriaIndex = nMinDistTriaIndex ;
+   m_bIsInside = false ;
+
+  // se le informazioni non sono coerenti, allora :
+  // 1) calcolo i centroidi dei triangoli in questione
+  // 2) ottengo il punto medio di questi centroidi
+  // 3) controllo quale triangolo interseca il segmento che parte da ptP e arriva a tale punto
+  // 4) userò questo triangolo per classificare ptP
+   if ( bOutside == bInside) {
+     // calcolo il baricentro complessivo
+      Point3d ptBar_tot ;
+      for ( auto& Tria : vTria)
+         ptBar_tot += Tria.second.GetCentroid() ;
+      ptBar_tot /= int( vTria.size()) ;
+     // per ogni triangolo, cerco quello che interseca il segmento
+      for ( auto& Tria : vTria) {
+         Point3d ptInters1, ptInters2 ;
+         int nType = IntersLineTria( ptP, ptBar_tot, Tria.second, ptInters1, ptInters2) ;
+         if ( nType == ILTT_IN) { // se intersezione ho finito
+            DistPointTriangle( ptP, Tria.second).GetMinDistPoint( m_ptMinDistPoint) ;
+            m_bIsInside = ( ( ptP - m_ptMinDistPoint) * Tria.second.GetN() < - EPS_SMALL) ;
+            m_nMinDistTriaIndex = Tria.first ;
+            break ;
+         }
+      }
+   }
+   else // se informazioni coerenti
+      m_bIsInside = bInside ;
+
+   return ;
+
 }
 
 //----------------------------------------------------------------------------

DistPointTria.cpp

@@ -14,7 +14,7 @@
 //--------------------------- Include ----------------------------------------
 #include "stdafx.h"
 #include "ProjPlane.h"
-#include "DistPointLine.h"
+#include "/EgtDev/Include/EGkDistPointLine.h"
 #include "/EgtDev/Include/EGkDistPointTria.h"
 
 

EgtGeomKernel.vcxproj

@@ -328,6 +328,8 @@ copy $(TargetPath) \EgtProg\Dll64</Command>
     <ClInclude Include="..\Include\EGkCDeConeFrustumClosedSurfTm.h" />
     <ClInclude Include="..\Include\EGkCDeConvexTorusClosedSurfTm.h" />
     <ClInclude Include="..\Include\EGkCDeRectPrismoidClosedSurfTm.h" />
+    <ClInclude Include="..\Include\EGkDistLineLine.h" />
+    <ClInclude Include="..\Include\EGkDistPointLine.h" />
     <ClInclude Include="..\Include\EGkIntersCurveSurfTm.h" />
     <ClInclude Include="..\Include\EGkIntersLineBox.h" />
     <ClInclude Include="..\Include\EGkIntersLineVolZmap.h" />
@@ -595,12 +597,10 @@ copy $(TargetPath) \EgtProg\Dll64</Command>
     <ClInclude Include="CAvSimpleSurfFrMove.h" />
     <ClInclude Include="CAvToolSurfTm.h" />
     <ClInclude Include="CreateCurveAux.h" />
-    <ClInclude Include="DistLineLine.h" />
     <ClInclude Include="DistPointArc.h" />
     <ClInclude Include="DistPointCrvAux.h" />
     <ClInclude Include="DistPointCrvBezier.h" />
     <ClInclude Include="DistPointCrvComposite.h" />
-    <ClInclude Include="DistPointLine.h" />
     <ClInclude Include="DllMain.h" />
     <ClInclude Include="earcut.hpp" />
     <ClInclude Include="ExtDimension.h" />

EgtGeomKernel.vcxproj.filters

@@ -593,9 +593,6 @@
     <ClInclude Include="DistPointArc.h">
       <Filter>File di intestazione</Filter>
     </ClInclude>
-    <ClInclude Include="DistPointLine.h">
-      <Filter>File di intestazione</Filter>
-    </ClInclude>
     <ClInclude Include="DistPointCrvBezier.h">
       <Filter>File di intestazione</Filter>
     </ClInclude>
@@ -1118,9 +1115,6 @@
     <ClInclude Include="..\Include\EGkIntersPlaneBox.h">
       <Filter>File di intestazione\Include</Filter>
     </ClInclude>
-    <ClInclude Include="DistLineLine.h">
-      <Filter>File di intestazione</Filter>
-    </ClInclude>
     <ClInclude Include="CDeUtility.h">
       <Filter>File di intestazione</Filter>
     </ClInclude>
@@ -1214,6 +1208,12 @@
     <ClInclude Include="..\Include\EGkQuaternion.h">
       <Filter>File di intestazione\Include</Filter>
     </ClInclude>
+    <ClInclude Include="..\Include\EGkDistLineLine.h">
+      <Filter>File di intestazione\Include</Filter>
+    </ClInclude>
+    <ClInclude Include="..\Include\EGkDistPointLine.h">
+      <Filter>File di intestazione\Include</Filter>
+    </ClInclude>
   </ItemGroup>
   <ItemGroup>
     <ResourceCompile Include="EgtGeomKernel.rc">

Frame3d.cpp

@@ -37,6 +37,23 @@ Frame3d::Set( const Point3d& ptOrig, const Vector3d& vtDirX,
         ! m_vtVersZ.Normalize())
       return false ;
 
+  // se ci sono errori molto piccoli di ortogonalità, li correggo
+   double dOrtXZ = m_vtVersX * m_vtVersZ ;
+   if ( dOrtXZ > EPS_ZERO  &&  dOrtXZ < 10 * EPS_ZERO) {
+      m_vtVersX = OrthoCompo( m_vtVersX, m_vtVersZ) ;
+      m_vtVersX.Normalize() ;
+   }
+   double dOrtYX = m_vtVersY * m_vtVersX ;
+   if ( dOrtYX > EPS_ZERO  &&  dOrtYX < 10 * EPS_ZERO) {
+      m_vtVersY = OrthoCompo( m_vtVersY, m_vtVersX) ;
+      m_vtVersY.Normalize() ;
+   }
+   double dOrtYZ = m_vtVersY * m_vtVersZ ;
+   if ( dOrtYZ > EPS_ZERO  &&  dOrtYZ < 10 * EPS_ZERO) {
+      m_vtVersY = OrthoCompo( m_vtVersY, m_vtVersZ) ;
+      m_vtVersY.Normalize() ;
+   }
+
   // verifica della ortogonalità dei versori e del senso destrorso
    if ( ! Verify())
       return false ;

IntersLineCaps.cpp

@@ -14,7 +14,7 @@
 //--------------------------- Include ----------------------------------------
 #include "stdafx.h"
 #include "IntersLineCaps.h"
-#include "DistLineLine.h"
+#include "/EgtDev/Include/EGkDistLineLine.h"
 #include "/EgtDev/Include/EGkIntersLineSphere.h"
 
 using namespace std ;

IntersLineSurfBez.cpp

@@ -13,12 +13,12 @@
 
 //--------------------------- Include ----------------------------------------
 #include "stdafx.h"
+#include "CurveLine.h"
+#include "/EgtDev/Include/EGkDistPointLine.h"
 #include "/EgtDev/Include/EGkIntersLineTria.h"
 #include "/EgtDev/Include/EGkIntersLineSurfTm.h"
 #include "/EgtDev/Include/EGkIntersLineSurfBez.h"
 #include "/EgtDev/Include/EGkSurfBezier.h"
-#include "DistPointLine.h"
-#include "CurveLine.h"
 
 using namespace std ;
 

IntersLineSurfTm.cpp

@@ -52,6 +52,8 @@ OrderInfoIntersLineSurfTm( ILSIVECTOR& vInfo)
             []( const IntLinStmInfo& a, const IntLinStmInfo& b)
                {  double dUa = ( ( a.nILTT == ILTT_SEGM  ||  a.nILTT == ILTT_SEGM_ON_EDGE) ? ( a.dU + a.dU2) / 2 : a.dU) ;
                   double dUb = ( ( b.nILTT == ILTT_SEGM  ||  b.nILTT == ILTT_SEGM_ON_EDGE) ? ( b.dU + b.dU2) / 2 : b.dU) ;
+                  if ( abs( dUa - dUb) < EPS_SMALL)
+                     return ( a.dCosDN < b.dCosDN) ;
                   return ( dUa < dUb) ; }) ;
 }
 

IntersLineTria.cpp

@@ -17,8 +17,8 @@
 #include "CurveLine.h"
 #include "IntersLineLine.h"
 #include "IntersLineTria.h"
-#include "DistPointLine.h"
-#include "DistLineLine.h"
+#include "/EgtDev/Include/EGkDistPointLine.h"
+#include "/EgtDev/Include/EGkDistLineLine.h"
 #include "/EgtDev/Include/EGkIntersLinePlane.h"
 #include "/EgtDev/Include/EGkFrame3d.h"
 #include <array>

IntersPlaneSurfTm.cpp

@@ -15,7 +15,7 @@
 #include "stdafx.h"
 #include "ProjPlane.h"
 #include "CurveLine.h"
-#include "DistPointLine.h"
+#include "/EgtDev/Include/EGkDistPointLine.h"
 #include "/EgtDev/Include/EGkIntersPlaneSurfTm.h"
 #include "/EgtDev/Include/EGkIntersPlaneTria.h"
 #include "/EgtDev/Include/EGkIntersLineSurfTm.h"

IntersSurfTmSurfTm.cpp

@@ -13,9 +13,9 @@
 
 //--------------------------- Include ----------------------------------------
 #include "stdafx.h"
-#include "DistPointLine.h"
 #include "IntersLineTria.h"
 #include "DllMain.h"
+#include "/EgtDev/Include/EGkDistPointLine.h"
 #include "/EgtDev/Include/EGkIntersSurfTmSurfTm.h"
 #include "/EgtDev/Include/EGkIntersTriaTria.h"
 #include "/EgtDev/Include/EGkPointGrid3d.h"

LineTgCurvePerpCurve.cpp

@@ -14,7 +14,7 @@
 //--------------------------- Include ----------------------------------------
 #include "stdafx.h"
 #include "CreateCurveAux.h"
-#include "DistPointLine.h"
+#include "/EgtDev/Include/EGkDistPointLine.h"
 #include "/EgtDev/Include/EGkLineTgCurvePerpCurve.h"
 #include "/EgtDev/Include/EGkLinePntTgCurve.h"
 #include "/EgtDev/Include/EgtPointerOwner.h"

LineTgTwoCurves.cpp

@@ -14,7 +14,7 @@
 //--------------------------- Include ----------------------------------------
 #include "stdafx.h"
 #include "CreateCurveAux.h"
-#include "DistPointLine.h"
+#include "/EgtDev/Include/EGkDistPointLine.h"
 #include "/EgtDev/Include/EGkLineTgTwoCurves.h"
 #include "/EgtDev/Include/EgtPointerOwner.h"
 

PolyArc.cpp

@@ -13,9 +13,9 @@
 
 //--------------------------- Include ----------------------------------------
 #include "stdafx.h"
-#include "DistPointLine.h"
 #include "GeoConst.h"
 #include "CurveArc.h"
+#include "/EgtDev/Include/EGkDistPointLine.h"
 #include "/EgtDev/Include/EGkPolyArc.h"
 #include "/EgtDev/Include/EGkFrame3d.h"
 #include <algorithm>

PolyLine.cpp

@@ -1,4 +1,4 @@
-//----------------------------------------------------------------------------
+//----------------------------------------------------------------------------
 //  EgalTech 2013-2013
 //----------------------------------------------------------------------------
 // File : PolyLine.cpp          Data : 22.12.13          Versione : 1.4l3
@@ -14,13 +14,13 @@
 //--------------------------- Include ----------------------------------------
 #include "stdafx.h"
 #include "CurveLine.h"
-#include "DistPointLine.h"
 #include "IntersLineLine.h"
 #include "PolygonPlane.h"
 #include "PointsPCA.h"
 #include "GeoConst.h"
 #include "/EgtDev/Include/EGkPolyLine.h"
 #include "/EgtDev/Include/EGkPlane3d.h"
+#include "/EgtDev/Include/EGkDistPointLine.h"
 #include "/EgtDev/Include/EGnStringUtils.h"
 #include "/EgtDev/Include/EgtNumUtils.h"
 
@@ -53,7 +53,7 @@ PolyLine::AddUPoint( double dPar, const Point3d& ptP, bool bEndOrStart)
 {
   // se da aggiungere in coda
    if ( bEndOrStart) {
-     // se il punto è uguale all'ultimo (ignoro parametro), non lo inserisco ma ok
+     // se il punto è uguale all'ultimo (ignoro parametro), non lo inserisco ma ok
       if ( m_lUPoints.size() > 0  &&  AreSamePointApprox( ptP, m_lUPoints.back().first)) {
          ++ m_nRejected ;
          return true ;
@@ -68,7 +68,7 @@ PolyLine::AddUPoint( double dPar, const Point3d& ptP, bool bEndOrStart)
    }
   // altrimenti si aggiunge in testa
    else {
-     // se il punto è uguale al primo (ignoro parametro), non lo inserisco ma ok
+     // se il punto è uguale al primo (ignoro parametro), non lo inserisco ma ok
       if ( m_lUPoints.size() > 0  &&  AreSamePointApprox( ptP, m_lUPoints.front().first)) {
          ++ m_nRejected ;
          return true ;
@@ -92,7 +92,7 @@ PolyLine::Close( void)
   // ci devono essere almeno 2 punti
    if ( m_lUPoints.size() < 2)
       return false ;
-  // verifico non sia già chiuso
+  // verifico non sia già chiuso
    if ( AreSamePointApprox( m_lUPoints.front().first, m_lUPoints.back().first))
       return false ;
   // aggiungo un punto uguale al primo in coda
@@ -219,7 +219,7 @@ PolyLine::ToLoc( const Frame3d& frRef)
 bool
 PolyLine::LocToLoc( const Frame3d& frOri, const Frame3d& frDest)
 {
-  // se i due riferimenti coincidono, non devo fare alcunché
+  // se i due riferimenti coincidono, non devo fare alcunché
    if ( AreSameFrame( frOri, frDest))
       return true ;
   // ciclo sui punti
@@ -233,7 +233,7 @@ PolyLine::LocToLoc( const Frame3d& frOri, const Frame3d& frDest)
 bool
 PolyLine::Join( PolyLine& PL, double dOffsetPar)
 {
-  // se l'altra polilinea non contiene alcunchè, esco con ok
+  // se l'altra polilinea non contiene alcunchè, esco con ok
    if ( PL.m_lUPoints.size() == 0)
       return true ;
   // verifico che l'ultimo punto di questa polilinea coincida con il primo dell'altra
@@ -385,7 +385,7 @@ PolyLine::GetPrevUPoint( double* pdPar, Point3d* pptP, bool bNotFirst) const
 bool
 PolyLine::GetCurrUPoint( double* pdPar, Point3d* pptP) const
 {
-  // verifico validità punto corrente
+  // verifico validità punto corrente
    if ( m_iter == m_lUPoints.end())
       return false ;
 
@@ -426,7 +426,7 @@ PolyLine::GetFirstULine( double* pdIni, Point3d* pptIni, double* pdFin, Point3d*
 bool
 PolyLine::GetNextULine( double* pdIni, Point3d* pptIni, double* pdFin, Point3d* pptFin) const
 {
-  // parametro e punto iniziali (è il precedente finale)
+  // parametro e punto iniziali (è il precedente finale)
    if ( m_iter == m_lUPoints.end())
       return false ;
    if ( pdIni != nullptr)
@@ -510,19 +510,19 @@ PolyLine::IsFlat( int& nRank, Point3d& ptCen, Vector3d& vtDir, double dToler) co
       ptsPCA.AddPoint( Media( ptP1, ptP2, 0.25), dLen / 2) ;
       ptsPCA.AddPoint( Media( ptP1, ptP2, 0.75), dLen / 2) ;
    }
-  // recupero il rango, ovvero la dimensionalità dell'insieme di punti
+  // recupero il rango, ovvero la dimensionalità dell'insieme di punti
    nRank = ptsPCA.GetRank() ;
   // se dimensione nulla, o non ci sono punti o sono tutti praticamente coincidenti
    if ( nRank == 0)
       return ptsPCA.GetCenter( ptCen) ;
   // se dimensione 1, allora i punti sono distribuiti su una linea
    if ( nRank == 1) {
-     // assegno il centro e la direzione della linea (il verso è indifferente)
+     // assegno il centro e la direzione della linea (il verso è indifferente)
       ptsPCA.GetCenter( ptCen) ;
       ptsPCA.GetPrincipalComponent( 0, vtDir) ;
       return true ;
    }
-  // altrimenti dimensione 2 o 3, allora è determinato un piano principale, verifico se tutti i punti vi giacciono
+  // altrimenti dimensione 2 o 3, allora è determinato un piano principale, verifico se tutti i punti vi giacciono
    // Center and normal vector
    ptsPCA.GetCenter( ptCen) ;
    Vector3d vtX, vtY ;
@@ -530,9 +530,9 @@ PolyLine::IsFlat( int& nRank, Point3d& ptCen, Vector3d& vtDir, double dToler) co
    ptsPCA.GetPrincipalComponent( 1, vtY) ;
    vtDir = vtX ^ vtY ;
    if ( ! vtDir.Normalize()) {
-     // riduco la dimensionalità a lineare
+     // riduco la dimensionalità a lineare
       nRank = 1 ;
-     // assegno il centro e la direzione della linea (il verso è indifferente)
+     // assegno il centro e la direzione della linea (il verso è indifferente)
       ptsPCA.GetCenter( ptCen) ;
       vtDir = vtX ;
       return true ;
@@ -561,12 +561,12 @@ PolyLine::IsFlat( Plane3d& plPlane, double dToler) const
       plPlane.Reset() ;
       return false ;
    }
-  // recupero dati sulla planarità della polilinea
+  // recupero dati sulla planarità della polilinea
    int nRank ;
    Point3d ptCen ;
    Vector3d vtDir ;
    bool bFlat = IsFlat( nRank, ptCen, vtDir, dToler) ;
-   // imposto il piano a seconda della dimensionalità
+   // imposto il piano a seconda della dimensionalità
    switch ( nRank) {
    case 0 :      // punto
       plPlane.Set( ptCen, Z_AX) ;
@@ -639,13 +639,13 @@ PolyLine::GetAreaXY( double& dArea) const
   // verifico sia chiusa
    if ( ! IsClosed())
       return false ;
-  // calcolo l'area considerando solo XY (è la Z di Newell)
+  // calcolo l'area considerando solo XY (è la Z di Newell)
    dArea = 0 ;
    Point3d ptIni, ptFin ;
    for ( bool bFound = GetFirstLine( ptIni, ptFin) ; bFound ; bFound = GetNextLine( ptIni, ptFin)) {
       dArea += ( ptIni.x - ptFin.x) * ( ptIni.y + ptFin.y) ;   // projection on xy
    }
-  // considero anche la linea tra l'ultimo e il primo punto perchè in alcuni casi potrebbero definire area 
+  // considero anche la linea tra l'ultimo e il primo punto perchè in alcuni casi potrebbero definire area 
   // significativa anche se sono coincidenti per le nostre tolleranze
    ptIni = ptFin ;
    GetFirstPoint( ptFin) ;
@@ -746,7 +746,7 @@ DouglasPeuckerSimplification( const PNTUVECTOR& vPtU, const double dSqTol, const
       }
    }
 
-  // se la distanza massima trovata è sopra la tolleranza, allora controllo la parte di PolyLine tra
+  // se la distanza massima trovata è sopra la tolleranza, allora controllo la parte di PolyLine tra
   // (nIndStart, nMaxInd) e quella tra (nMaxInd, nIndEnd)
    if ( dMaxSqDist > dSqTol) {
      // inserisco il punto
@@ -789,7 +789,7 @@ PolyLine::RemoveAlignedPoints( double dToler)
    }
   // altrimenti chiusa
    else {
-     // cerco il punto più distante dal primo
+     // cerco il punto più distante dal primo
       double dMaxDist = 0. ;
       int nMaxInd = 0 ;
       for ( int i = 1 ; i < int( vPtU.size()) ; ++ i) {
@@ -812,6 +812,14 @@ PolyLine::RemoveAlignedPoints( double dToler)
   // ordino in senso crescente
    sort( vInd.begin(), vInd.end()) ;
 
+  // se chiusa e almeno 4 punti rimasti, controllo allineamento dell'inizio con precedente e successivo rimasti
+   if ( IsClosed()  &&  vInd.size() >= 4) {
+      if ( DistPointLine( vPtU[vInd[0]].first, vPtU[vInd[1]].first, vPtU[vInd[vInd.size()-2]].first).IsEpsilon( dToler)) {
+         vInd.erase( vInd.begin()) ;
+         vInd.back() = vInd.front() ;
+      }
+   }
+
   // rimetto in lista i soli punti rimasti
    m_lUPoints.clear() ;
    for ( auto Ind : vInd)
@@ -958,7 +966,7 @@ PolyLine::MyApproxOnSide( const Vector3d& vtN,  bool bLeftSide, double dToler)
             }
          }
       }
-     // non è stato eliminato alcunché
+     // non è stato eliminato alcunché
      // ripristino la tolleranza corrente
       dCurrToler = dToler ;
      // avanzo il terzetto di uno step
@@ -999,7 +1007,7 @@ PolyLine::MakeConvex( const Vector3d& vtN,  bool bLeftSide)
 bool
 PolyLine::MyMakeConvex( const Vector3d& vtN,  bool bLeftSide)
 {
-  // ciclo i controlli finchè non ci sono rimozioni
+  // ciclo i controlli finchè non ci sono rimozioni
    bool bRemoved = true ;
    while ( bRemoved) {
       bRemoved = false ;
@@ -1027,7 +1035,7 @@ PolyLine::MyMakeConvex( const Vector3d& vtN,  bool bLeftSide)
             bRemoved = true ;
             continue ;
          }
-        // non è stato eliminato alcunché : avanzo il terzetto di uno step
+        // non è stato eliminato alcunché : avanzo il terzetto di uno step
          precP = currP ;
          currP = nextP ;
          ++ nextP ;
@@ -1054,7 +1062,7 @@ PolyLine::Invert( bool bInvertU)
    m_lUPoints.reverse() ;
   // se richiesto, inverto anche il parametro U
    if ( bInvertU) {
-     // recupero il primo valore di U che è il vecchio finale ed è il riferimento di inversione
+     // recupero il primo valore di U che è il vecchio finale ed è il riferimento di inversione
       double dUfin = m_lUPoints.front().second ;
      // ciclo su tutti gli elementi
       for ( auto& UPoint : m_lUPoints) {
@@ -1269,7 +1277,7 @@ PolyLine::GetMinAreaRectangleXY( Point3d& ptCen, Vector3d& vtAx, double& dLen, d
 bool
 PolyLine::Trim( const Plane3d& plPlane, bool bInVsOut)
 {
-  // se vuota non faccio alcunché
+  // se vuota non faccio alcunché
    if ( m_lUPoints.size() == 0)
       return false ;
 
@@ -1371,7 +1379,7 @@ IsPointInsidePolyLine( const Point3d& ptP, const PolyLine& plPoly, double dToler
       return false ;
   // Riferimento alla lista dei punti
    PNTULIST& List = const_cast<PolyLine&>( plPoly).GetUPointList() ;
-  // Ciclo sui segmenti della polilinea per cercare il segmento più vicino al punto
+  // Ciclo sui segmenti della polilinea per cercare il segmento più vicino al punto
    double dMinSqDist = SQ_INFINITO ;
    Point3d ptMinDist ;
    auto itMinDistEnd = List.end() ;
@@ -1460,7 +1468,7 @@ GetPointParamOnPolyLine( const Point3d& ptP, const PolyLine& plPoly, double dTol
      // assegno nuovo inizio
       ptStart = ptEnd ;
    }
-  // Il punto è sulla linea se la sua distanza rispetta la tolleranza
+  // Il punto è sulla linea se la sua distanza rispetta la tolleranza
    return ( dMinSqDist < dToler * dToler) ;
 }
 
@@ -1473,7 +1481,7 @@ ChangePolyLineStart( PolyLine& plPoly, const Point3d& ptNewStart, double dToler)
       return false ;
   // Riferimento alla lista dei punti
    PNTULIST& LoopList = const_cast<PolyLine&>( plPoly).GetUPointList() ;
-  // Ciclo sui segmenti della polilinea per cercare il segmento più vicino al punto
+  // Ciclo sui segmenti della polilinea per cercare il segmento più vicino al punto
    double dMinSqDist = SQ_INFINITO ;
    auto itMinDistEnd = LoopList.end() ;
    auto itStart = LoopList.begin() ;
@@ -1517,7 +1525,7 @@ SplitPolyLineAtPoint( const PolyLine& plPoly, const Point3d& ptP, double dToler,
       return false ;
   // Riferimento alla lista dei punti
    const PNTULIST& LoopList = const_cast<PolyLine&>( plPoly).GetUPointList() ;
-  // Ciclo sui segmenti della polilinea per cercare il segmento più vicino al punto
+  // Ciclo sui segmenti della polilinea per cercare il segmento più vicino al punto
    double dMinSqDist = SQ_INFINITO ;
    auto itMinDistEnd = LoopList.end() ;
    auto itStart = LoopList.begin() ;
@@ -1581,10 +1589,6 @@ AssociatePolyLinesMinDistPoints( const PolyLine& PL1, const PolyLine& PL2, PNTIV
    int nLastJ = 0 ;
    vPnt1[0].second = 0 ;
 
-   double dFirstDist, dFirstParMinDist ;
-   DistPointPolyLine( vPnt1[0].first, PL2, dFirstDist, dFirstParMinDist) ;      
-   int nFirstMinJ = ( int)( dFirstParMinDist + 0.5) ;
-
    for ( int i = 1 ; i < nTotP1 ; ++ i) {
 
       double dDist = INFINITO ;
@@ -1601,14 +1605,10 @@ AssociatePolyLinesMinDistPoints( const PolyLine& PL1, const PolyLine& PL2, PNTIV
       }
       int nMinJ = ( int)( dMinDistPar + 0.5) ;
 
-     // eventuale correzione per i primi punti ( da forzare nel vertice 0)
-      if ( nLastJ == 0 && nFirstMinJ > 0.5 * nTotP2  && nMinJ >= nFirstMinJ) 
-         nMinJ = 0 ;
-
       if ( nMinJ < nLastJ)
          nMinJ = nLastJ ;
 
-     // verifica se è un punto interno in comune con l'altra polyline  
+     // verifica se è un punto interno in comune con l'altra polyline  
       if ( i < nTotP1 - 1  &&  dDist < EPS_SMALL  && abs( dMinDistPar - floor( dMinDistPar + 0.5)) < EPS_SMALL)
          bCommonInternalPoints = true ;
 
@@ -1619,9 +1619,6 @@ AssociatePolyLinesMinDistPoints( const PolyLine& PL1, const PolyLine& PL2, PNTIV
   // calcoli per seconda curva
    int nLastI = 0 ;
    vPnt2[0].second = 0 ;
-   
-   DistPointPolyLine( vPnt2[0].first, PL1, dFirstDist, dFirstParMinDist) ;      
-   int nFirstMinI = ( int)( dFirstParMinDist + 0.5) ;
 
    for ( int j = 1 ; j < nTotP2 ; ++ j) {
 
@@ -1640,10 +1637,6 @@ AssociatePolyLinesMinDistPoints( const PolyLine& PL1, const PolyLine& PL2, PNTIV
       }
       int nMinI = ( int)( dMinDistPar + 0.5) ;
 
-     // eventuale correzione per primi punti 
-      if ( nLastI == 0 && nFirstMinI > 0.5 * nTotP1  && nMinI >= nFirstMinI) 
-         nMinI = 0 ;
-
       if ( nMinI < nLastI)
          nMinI = nLastI ;
 

ProjectCurveSurfTm.cpp

@@ -13,8 +13,8 @@
 
 //--------------------------- Include ----------------------------------------
 #include "stdafx.h"
-#include "DistPointLine.h"
 #include "GeoConst.h"
+#include "/EgtDev/Include/EGkDistPointLine.h"
 #include "/EgtDev/Include/EGkDistPointCurve.h"
 #include "/EgtDev/Include/EGkDistPointSurfTm.h"
 #include "/EgtDev/Include/EGkIntersLineSurfTm.h"
@@ -22,6 +22,10 @@
 
 using namespace std ;
 
+//----------------------------------------------------------------------------
+// Angolo limite tra normale al triangolo e direzione di proiezione 89°
+const double COS_ANG_LIM = 0.0175 ;
+
 //----------------------------------------------------------------------------
 static bool
 PointsInTolerance( const PNT5AXVECTOR& vPt5ax, int nPrec, int nCurr, int nNext, double dSqTol)
@@ -89,9 +93,9 @@ ProjectCurveOnSurfTm( const ICurve& crCrv, const ISurfTriMesh& tmSurf, const Vec
   // controllo le tolleranze
    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 alla massima risoluzione
    PolyLine PL ;
-   if ( ! crCrv.ApproxWithLines( dLinTol, ANG_TOL_STD_DEG, ICurve::APL_STD, PL))
+   if ( ! crCrv.ApproxWithLines( EPS_SMALL, ANG_TOL_STD_DEG, ICurve::APL_STD, PL))
       return false ;
    const double MAX_SEG_LEN = min( dMaxSegmLen, 1.) ;
    if ( ! PL.AdjustForMaxSegmentLen( MAX_SEG_LEN))
@@ -112,12 +116,17 @@ ProjectCurveOnSurfTm( const ICurve& crCrv, const ISurfTriMesh& tmSurf, const Vec
    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) ;
-      if ( ! vIntRes.empty()) {
+     // 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
-         int nI = int( vIntRes.size()) - 1 ;
          Point3d ptInt ;
          if ( vIntRes[nI].nILTT == ILTT_SEGM  ||  vIntRes[nI].nILTT == ILTT_SEGM_ON_EDGE)
             ptInt = vIntRes[nI].ptI2 ;
@@ -187,9 +196,13 @@ ProjectCurveOnSurfTm( const ICurve& crCrv, const ISurfTriMesh& tmSurf, const IGe
          vtLine /= dLineLen ;
          ILSIVECTOR vIntRes ;
          if ( IntersLineSurfTm( ptP, vtLine, dLineLen, tmSurf, vIntRes, false)) {
-            if ( vIntRes.size() > 0) {
+           // 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
-               int nI = int( vIntRes.size()) - 1 ;
                Point3d ptInt ;
                if ( vIntRes[nI].nILTT == ILTT_SEGM  ||  vIntRes[nI].nILTT == ILTT_SEGM_ON_EDGE)
                   ptInt = vIntRes[nI].ptI2 ;
@@ -235,9 +248,9 @@ ProjectCurveOnSurfTm( const ICurve& crCrv, const ISurfTriMesh& tmSurf, const ICu
    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 alla massima risoluzione
    PolyLine PL ;
-   if ( ! crCrv.ApproxWithLines( dLinTol, ANG_TOL_STD_DEG, ICurve::APL_STD, PL))
+   if ( ! crCrv.ApproxWithLines( EPS_SMALL, ANG_TOL_STD_DEG, ICurve::APL_STD, PL))
       return false ;
    const double MAX_SEG_LEN = min( dMaxSegmLen, 1.) ;
    if ( ! PL.AdjustForMaxSegmentLen( MAX_SEG_LEN))
@@ -264,9 +277,13 @@ ProjectCurveOnSurfTm( const ICurve& crCrv, const ISurfTriMesh& tmSurf, const ICu
             vtLine /= dLineLen ;
             ILSIVECTOR vIntRes ;
             if ( IntersLineSurfTm( ptP, vtLine, dLineLen, tmSurf, vIntRes, false)) {
-               if ( vIntRes.size() > 0) {
+              // 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
-                  int nI = int( vIntRes.size()) - 1 ;
                   Point3d ptInt ;
                   if ( vIntRes[nI].nILTT == ILTT_SEGM  ||  vIntRes[nI].nILTT == ILTT_SEGM_ON_EDGE)
                      ptInt = vIntRes[nI].ptI2 ;
@@ -342,9 +359,13 @@ ProjectCurveOnSurfTm( const ICurve& crCrv, const ISurfTriMesh& tmSurf, const ISu
             vtLine /= dLineLen ;
             ILSIVECTOR vIntRes ;
             if ( IntersLineSurfTm( ptP, vtLine, dLineLen, tmSurf, vIntRes, false)) {
-               if ( vIntRes.size() > 0) {
+              // 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
-                  int nI = int( vIntRes.size()) - 1 ;
                   Point3d ptInt ;
                   if ( vIntRes[nI].nILTT == ILTT_SEGM  ||  vIntRes[nI].nILTT == ILTT_SEGM_ON_EDGE)
                      ptInt = vIntRes[nI].ptI2 ;

RemoveCurveDefects.cpp

@@ -15,8 +15,8 @@
 #include "stdafx.h"
 #include "GeoConst.h"
 #include "CurveComposite.h"
-#include "DistPointLine.h"
 #include "RemoveCurveDefects.h"
+#include "/EgtDev/Include/EGkDistPointLine.h"
 #include "/EgtDev/Include/EGkDistPointCurve.h"
 #include <algorithm>
 

StmFromCurves.cpp

@@ -330,9 +330,11 @@ GetSurfTriMeshSharpRectSwept( double dDimH, double dDimV, const ICurve* pGuide,
 {
   // verifico che la linea guida sia piana
    Plane3d plGuide ;
-   if ( ! pGuide->IsFlat( plGuide, false, 10 * EPS_SMALL))
+   if ( ! pGuide->IsFlat( plGuide, true, 10 * EPS_SMALL))
       return nullptr ;
-   Vector3d vtNorm = plGuide.GetVersN() ;
+   Vector3d vtNorm ; pGuide->GetExtrusion( vtNorm) ;
+   if ( vtNorm.IsSmall())
+      vtNorm = Z_AX ;
   // determino se la guida è chiusa
    bool bGuideClosed = pGuide->IsClosed() ;
   // curve di offset
@@ -355,7 +357,7 @@ GetSurfTriMeshSharpRectSwept( double dDimH, double dDimV, const ICurve* pGuide,
    PtrOwner<ISurfTriMesh> pSrfBot( pSrfTop->Clone()) ;
    if ( IsNull( pSrfBot))
       return nullptr ;
-   pSrfBot->Translate( -dDimV * vtNorm) ;
+   pSrfBot->Translate( - dDimV * vtNorm) ;
    pSrfBot->Invert() ;
    PtrOwner<ISurfTriMesh> pSrfRgt( GetSurfTriMeshByExtrusion( pCrvR, -dDimV * vtNorm, false, dLinTol)) ;
    if ( IsNull( pSrfRgt))
@@ -365,74 +367,188 @@ GetSurfTriMeshSharpRectSwept( double dDimH, double dDimV, const ICurve* pGuide,
    if ( IsNull( pSrfLft))
       return nullptr ;
   // unisco le parti
-   PtrOwner<ISurfTriMesh> pSTM( Release( pSrfTop)) ;
-   pSTM->DoSewing( *pSrfRgt) ;
-   pSTM->DoSewing( *pSrfLft) ;
-   pSTM->DoSewing( *pSrfBot) ;
-  // salvo tolleranza lineare usata e imposto angolo per smooth
-   pSTM->SetLinearTolerance( dLinTol) ;
-   pSTM->SetSmoothAngle( 20) ;
+   int nBuckets = max( 2 * ( pSrfRgt->GetVertexSize() + pSrfLft->GetVertexSize()), 1000) ;
+   StmFromTriangleSoup stmSoup ;
+   if ( ! stmSoup.Start( nBuckets))
+      return nullptr ;
+   stmSoup.AddSurfTriMesh( *pSrfTop) ;
+   stmSoup.AddSurfTriMesh( *pSrfRgt) ;
+   stmSoup.AddSurfTriMesh( *pSrfLft) ;
+   stmSoup.AddSurfTriMesh( *pSrfBot) ;
+   PtrOwner<ISurfTriMesh> pSTM ;
   // se guida aperta e tappi piatti
    if ( ! bGuideClosed  &&  nCapType == RSCAP_FLAT) {
+     // completo unione e recupero la superficie risultante
+      if ( ! stmSoup.End())
+         return nullptr ;
+      pSTM.Set( stmSoup.GetSurf()) ;
      // verifico che le due estremità siano chiuse e piatte
       POLYLINEVECTOR vPL ;
       if ( ! pSTM->GetLoops( vPL)  ||  vPL.size() != 2)
          return nullptr ;
       Plane3d plEnds ; double dArea ;
-      if ( ! vPL[0].IsClosedAndFlat( plEnds, dArea, 100 * EPS_SMALL))
+      if ( ! vPL[0].IsClosedAndFlat( plEnds, dArea, 50 * EPS_SMALL))
          return nullptr ;
-      if ( ! vPL[1].IsClosedAndFlat( plEnds, dArea, 100 * EPS_SMALL))
+      if ( ! vPL[1].IsClosedAndFlat( plEnds, dArea, 50 * EPS_SMALL))
          return nullptr ;
-     // aggiungo il cap sull'inizio
+     // calcolo il cap sull'inizio
       PtrOwner<SurfTriMesh> pSci( CreateBasicSurfTriMesh()) ;
       if ( IsNull( pSci)  ||  ! pSci->CreateByFlatContour( vPL[0]))
          return nullptr ;
       pSci->Invert() ;
-      pSTM->DoSewing( *pSci) ;
-     // aggiungo il cap sulla fine
+     // calcolo il cap sulla fine
       PtrOwner<SurfTriMesh> pSce( CreateBasicSurfTriMesh()) ;
       if ( IsNull( pSce)  ||  ! pSce->CreateByFlatContour( vPL[1]))
          return nullptr ;
       pSce->Invert() ;
-      pSTM->DoSewing( *pSce) ;
+     // cucio i tappi all'estrusione
+      if ( ! pSTM->DoSewing( *pSci)  ||  ! pSTM->DoSewing( *pSce))
+         return nullptr ;
    }
   // se altrimenti guida aperta e tappi arrotondati
    if ( ! bGuideClosed  &&  ( nCapType == RSCAP_ROUND  ||  nCapType == RSCAP_BEVEL)) {
      // step di rotazione per rispettare la tolleranza
       double dStepRotDeg = ( nCapType == RSCAP_BEVEL ? ANG_STRAIGHT / 4 : sqrt( 8 * dLinTol / dDimH) * RADTODEG) ;
-     // aggiungo il cap sull'inizio
-      Point3d ptStart ;
-      pGuide->GetStartPoint( ptStart) ;
-      Vector3d vtStart ;
-      pGuide->GetStartDir( vtStart) ;
-      vtStart.Rotate( vtNorm, 0, 1) ;
-      PolyLine PLStart ;
-      PLStart.AddUPoint( 0, ptStart) ;
-      PLStart.AddUPoint( 1, ptStart + dDimH / 2 * vtStart) ;
-      PLStart.AddUPoint( 2, ptStart + dDimH / 2 * vtStart - dDimV * vtNorm) ;
-      PLStart.AddUPoint( 3, ptStart - dDimV * vtNorm) ;
-      PtrOwner<SurfTriMesh> pSci( CreateBasicSurfTriMesh()) ;
-      if ( IsNull( pSci)  ||  ! pSci->CreateByScrewing( PLStart, ptStart, vtNorm, ANG_STRAIGHT, dStepRotDeg, 0))
+     // se l'offset interno alla guida è chiuso
+      if ( pCrvL->IsClosed()) {
+        // calcolo l'angolo di rotazione per screwing faccia Top e Bottom
+         Point3d ptRight ; pCrvR->GetEndPoint( ptRight) ;
+         Point3d ptLeft ; pCrvR->GetStartPoint( ptLeft) ;
+         Point3d ptJunction ; pCrvL->GetStartPoint( ptJunction) ;
+         Point3d ptCenter = Media( ptRight, ptLeft) ;
+         Vector3d vtRight = ptRight - ptCenter ;
+         Vector3d vtLeft = ptLeft - ptCenter ;
+         double dAng = ANG_STRAIGHT ;
+         vtRight.GetAngle( vtLeft, dAng) ;
+         vtRight.Normalize() ;
+         PolyLine plLoop ;
+        // creo il loop defininendo i punti
+         plLoop.AddUPoint( 0, ptRight) ; // primo punto
+         double dAngStep = ceil( dAng / dStepRotDeg) ; // aggiusto lo step
+         for ( int i = 1 ; i < dAngStep ; ++ i) {
+            Point3d ptRot = ptRight ;
+            ptRot.Rotate( ptCenter, vtNorm, i * ( dAng / dAngStep)) ;
+            plLoop.AddUPoint( i, ptRot) ; // punto intermedio sulla circonferenza
+         }
+         plLoop.AddUPoint( dAngStep ++, ptLeft) ; // ultimo punto
+         plLoop.AddUPoint( dAngStep ++, ptJunction) ; // punto centrale sull'offset chiuso
+         plLoop.AddUPoint( dAngStep, ptRight) ; // polyLine chiusa
+        // superificie Top
+         PtrOwner<ISurfTriMesh> pStmTop( CreateSurfTriMesh()) ;
+         if ( IsNull( pStmTop)  ||  ! pStmTop->CreateByFlatContour( plLoop))
+            return nullptr ;
+         stmSoup.AddSurfTriMesh( *pStmTop) ;
+        // superificie Bottom
+         PtrOwner<ISurfTriMesh> pStmBottom( CloneSurfTriMesh( pStmTop)) ;
+         pStmBottom->Translate( - dDimV * vtNorm) ;
+         pStmBottom->Invert() ;
+         stmSoup.AddSurfTriMesh( *pStmBottom) ;
+        // superificie perpendicolare
+        // la PolyLine che utilizzo la posso ricavare da quella calcolata sopra
+         plLoop.EraseLastUPoint() ; // apro il loop
+         plLoop.EraseLastUPoint() ; // tolgo il punto di contatto sull'offset
+         PtrOwner<ISurfTriMesh> pStmPerp( CreateSurfTriMesh()) ;
+         if ( IsNull( pStmPerp)  ||  ! pStmPerp->CreateByExtrusion( plLoop, - vtNorm * dDimV)  ||
+              ! pStmPerp->Invert())
+            return nullptr ;
+         stmSoup.AddSurfTriMesh( *pStmPerp) ;
+      }
+     // se l'offset interno della guida è aperto...
+      else {
+        // aggiungo il cap sull'inizio
+         Point3d ptStart ;
+         pGuide->GetStartPoint( ptStart) ;
+        // calcolo l'angolo di rotazione per screwing faccia Top e Bottom
+         Point3d ptSLeft ; pCrvL->GetStartPoint( ptSLeft) ;
+         Point3d ptSRight ; pCrvR->GetStartPoint( ptSRight) ;
+         Vector3d vtLeft = ptSLeft - ptStart ;
+         Vector3d vtRight = ptSRight - ptStart ;
+         double dAng = ANG_STRAIGHT ;
+         vtLeft.GetAngle( vtRight, dAng) ;
+         vtLeft.Normalize() ;
+         PolyLine plLoop ;
+        // creo il loop defininendo i punti
+         plLoop.AddUPoint( 0, ptSLeft) ;  // primo punto
+         double dAngStep = ceil( dAng / dStepRotDeg) ;
+         for ( int i = 1 ; i < dAngStep ; ++ i) {
+            Point3d ptRot = ptSLeft ;
+            ptRot.Rotate( ptStart, vtNorm, i * ( dAng / dAngStep)) ;
+            plLoop.AddUPoint( i, ptRot) ;
+         }
+         plLoop.AddUPoint( dAngStep, ptSRight) ; // ultimo punto
+         plLoop.AddUPoint( dAngStep + 1, ptSLeft) ; // polyline chiusa
+        // creo la superficie Top
+         PtrOwner<ISurfTriMesh> pStmTop_start( CreateSurfTriMesh()) ;
+         if ( IsNull( pStmTop_start)  ||  ! pStmTop_start->CreateByFlatContour( plLoop))
+            return nullptr ;
+         stmSoup.AddSurfTriMesh( *pStmTop_start) ;
+        // superificie Bottom
+         PtrOwner<ISurfTriMesh> pStmBottom_start( CloneSurfTriMesh( pStmTop_start)) ;
+         pStmBottom_start->Translate( - dDimV * vtNorm) ;
+         pStmBottom_start->Invert() ;
+         stmSoup.AddSurfTriMesh( *pStmBottom_start) ;
+        // superificie perpendicolare
+        // la PolyLine che utilizzo la posso ricavare da quella calcolata sopra
+         plLoop.EraseLastUPoint() ; // apro il loop
+         PtrOwner<ISurfTriMesh> pStmPerp_start( CreateSurfTriMesh()) ;
+         if ( IsNull( pStmPerp_start)  ||  ! pStmPerp_start->CreateByExtrusion( plLoop, - vtNorm * dDimV)  ||
+              ! pStmPerp_start->Invert())
+            return nullptr ;
+         stmSoup.AddSurfTriMesh( *pStmPerp_start) ;
+        // aggiungo il cap sulla fine
+         Point3d ptEnd ;
+         pGuide->GetEndPoint( ptEnd) ;
+        // calcolo l'angolo di rotazione per screwing faccia Top e Bottom
+         pCrvL->GetEndPoint( ptSLeft) ;
+         pCrvR->GetEndPoint( ptSRight) ;
+         vtLeft = ptSLeft - ptEnd ;
+         vtRight = ptSRight - ptEnd ;
+         dAng = ANG_STRAIGHT ;
+         vtRight.GetAngle( vtLeft, dAng) ;
+         vtRight.Normalize() ;
+         plLoop.Clear() ;
+        // creo il loop defininendo i punti
+         plLoop.AddUPoint( 0, ptSRight) ;
+         dAngStep = ceil( dAng / dStepRotDeg) ; // primo punto
+         for ( int i = 1 ; i < dAngStep ; ++ i) {
+            Point3d ptRot = ptSRight ;
+            ptRot.Rotate( ptEnd, vtNorm, i * ( dAng / dAngStep)) ;
+            plLoop.AddUPoint( i, ptRot) ;
+         }
+         plLoop.AddUPoint( dAngStep, ptSLeft) ; // ultimo punto
+         plLoop.AddUPoint( dAngStep + 1, ptSRight) ; // polyline chiusa
+        // creo la superficie Top
+         PtrOwner<ISurfTriMesh> pStmTop_end( CreateSurfTriMesh()) ;
+         if ( IsNull( pStmTop_end)  ||  ! pStmTop_end->CreateByFlatContour( plLoop))
+            return nullptr ;
+         stmSoup.AddSurfTriMesh( *pStmTop_end) ;
+        // creo la superificie Bottom
+         PtrOwner<ISurfTriMesh> pStmBottom_end( CloneSurfTriMesh( pStmTop_end)) ;
+         pStmBottom_end->Translate( - dDimV * vtNorm) ;
+         pStmBottom_end->Invert() ;
+         stmSoup.AddSurfTriMesh( *pStmBottom_end) ;
+        // creo la superificie perpendicolare alla guida
+         plLoop.EraseLastUPoint() ; // apro il loop
+         PtrOwner<ISurfTriMesh> pStmPerp_end( CreateSurfTriMesh()) ;
+         if ( IsNull( pStmPerp_end)  ||  ! pStmPerp_end->CreateByExtrusion( plLoop, - vtNorm * dDimV)  ||
+              ! pStmPerp_end->Invert())
+            return nullptr ;
+         stmSoup.AddSurfTriMesh( *pStmPerp_end) ;
+      }
+     // completo unione e recupero la superficie risultante
+      if ( ! stmSoup.End())
          return nullptr ;
-      pSci->Invert() ;
-      pSTM->DoSewing( *pSci) ;
-     // aggiungo il cap sulla fine
-      Point3d ptEnd ;
-      pGuide->GetEndPoint( ptEnd) ;
-      Vector3d vtEnd ;
-      pGuide->GetEndDir( vtEnd) ;
-      vtEnd.Rotate( vtNorm, 0, -1) ;
-      PolyLine PLEnd ;
-      PLEnd.AddUPoint( 0, ptEnd) ;
-      PLEnd.AddUPoint( 1, ptEnd + dDimH / 2 * vtEnd) ;
-      PLEnd.AddUPoint( 2, ptEnd + dDimH / 2 * vtEnd - dDimV * vtNorm) ;
-      PLEnd.AddUPoint( 3, ptEnd - dDimV * vtNorm) ;
-      PtrOwner<SurfTriMesh> pSce( CreateBasicSurfTriMesh()) ;
-      if ( IsNull( pSce)  ||  ! pSce->CreateByScrewing( PLEnd, ptEnd, vtNorm, ANG_STRAIGHT, dStepRotDeg, 0))
-         return nullptr ;
-      pSce->Invert() ;
-      pSTM->DoSewing( *pSce) ;
+      pSTM.Set( stmSoup.GetSurf()) ;
    }
+   else {
+     // completo unione e recupero la superficie risultante
+      if ( ! stmSoup.End())
+         return nullptr ;
+      pSTM.Set( stmSoup.GetSurf()) ;
+   }
+  // salvo tolleranza lineare usata e imposto angolo per smooth
+   pSTM->SetLinearTolerance( dLinTol) ;
+   pSTM->SetSmoothAngle( 20) ;
   // restituisco la superficie
    return Release( pSTM) ;
 }
@@ -445,13 +561,13 @@ GetSurfTriMeshBeveledRectSwept( double dDimH, double dDimV, double dBevelH, doub
 
   // verifico che la linea guida sia piana
    Plane3d plGuide ;
-   if ( ! pGuide->IsFlat( plGuide, false, 10 * EPS_SMALL))
+   if ( ! pGuide->IsFlat( plGuide, true, 10 * EPS_SMALL))
       return nullptr ;
-  // assegno la normale del piano
-   Vector3d vtNorm = plGuide.GetVersN() ;
+   Vector3d vtNorm ; pGuide->GetExtrusion( vtNorm) ;
+   if ( vtNorm.IsSmall())
+      vtNorm = Z_AX ;
   // determino il punto centrale della sezione
-   Point3d ptCen ;
-   pGuide->GetStartPoint( ptCen) ;
+   Point3d ptCen ; pGuide->GetStartPoint( ptCen) ;
    ptCen -= dDimV / 2 * vtNorm ;
   // determino se la guida è chiusa
    bool bGuideClosed = pGuide->IsClosed() ;
@@ -480,7 +596,6 @@ GetSurfTriMeshBeveledRectSwept( double dDimH, double dDimV, double dBevelH, doub
       for ( int i = 0 ; i < NUM_OFFS && bOk ; ++ i)
          bOk = vOffsCrv[i].Make( pGuide, vDist[i], ICurve::OFF_FILLET) ;
    }
-   
    if ( ! bOk  ||
         vOffsCrv[0].GetCurveCount() == 0  ||  vOffsCrv[1].GetCurveCount() == 0  ||
         vOffsCrv[2].GetCurveCount() == 0  ||  vOffsCrv[3].GetCurveCount() == 0)
@@ -549,10 +664,6 @@ GetSurfTriMeshBeveledRectSwept( double dDimH, double dDimV, double dBevelH, doub
       if ( ! stmSoup.End())
          return nullptr ;
       pSTM.Set( stmSoup.GetSurf()) ;
-     // preparo seconda zuppa di triangoli per inserire i tappi
-      StmFromTriangleSoup stmCapSoup ;
-      if ( ! stmCapSoup.Start( nBuckets))
-         return nullptr ;
      // verifico che le due estremità siano chiuse e piatte
       POLYLINEVECTOR vPL ;
       if ( ! pSTM->GetLoops( vPL)  ||  vPL.size() != 2)
@@ -580,42 +691,236 @@ GetSurfTriMeshBeveledRectSwept( double dDimH, double dDimV, double dBevelH, doub
    else if ( ! bGuideClosed  &&  ( nCapType == RSCAP_ROUND  ||  nCapType == RSCAP_BEVEL)) {
      // step di rotazione per rispettare il tipo o la tolleranza
       double dStepRotDeg = ( nCapType == RSCAP_BEVEL ? ANG_STRAIGHT / 4 : sqrt( 8 * dLinTol / dDimH) * RADTODEG) ;
-     // aggiungo il cap sull'inizio
-      Point3d ptStart ;
-      pGuide->GetStartPoint( ptStart) ;
-      Vector3d vtStart ;
-      pGuide->GetStartDir( vtStart) ;
-      vtStart.Rotate( vtNorm, 0, 1) ;
-      PolyLine PLStart ;
-      PLStart.AddUPoint( 0, ptStart) ;
-      PLStart.AddUPoint( 1, ptStart + ( dDimH / 2 - dBevelH) * vtStart) ;
-      PLStart.AddUPoint( 2, ptStart + dDimH / 2 * vtStart - dBevelV * vtNorm) ;
-      PLStart.AddUPoint( 3, ptStart + dDimH / 2 * vtStart - ( dDimV - dBevelV) * vtNorm) ;
-      PLStart.AddUPoint( 4, ptStart + ( dDimH / 2 - dBevelH) * vtStart - dDimV * vtNorm) ;
-      PLStart.AddUPoint( 5, ptStart - dDimV * vtNorm) ;
-      PtrOwner<SurfTriMesh> pSci( CreateBasicSurfTriMesh()) ;
-      if ( IsNull( pSci)  ||  ! pSci->CreateByScrewing( PLStart, ptStart, vtNorm, ANG_STRAIGHT, dStepRotDeg, 0))
+     // se l'offset interno della guida è chiuso...
+      if ( pCrvL->IsClosed()) {
+        // calcolo l'angolo di rotazione per screwing faccia Top e Bottom
+         Point3d ptRight ; pCrvR->GetEndPoint( ptRight) ;
+         Point3d ptLeft ; pCrvR->GetStartPoint( ptLeft) ;
+         Point3d ptJunction ; pCrvL->GetStartPoint( ptJunction) ;
+         Point3d ptCenter = Media( ptRight, ptLeft) ;
+         Vector3d vtRight = ptRight - ptCenter ;
+         Vector3d vtLeft = ptLeft - ptCenter ;
+         double dAng = ANG_STRAIGHT ;
+         vtRight.GetAngle( vtLeft, dAng) ;
+         vtRight.Normalize() ;
+         PolyLine plLoop ;
+        // creo il loop defininendo i punti
+         plLoop.AddUPoint( 0, ptRight) ; // primo punto
+         double dAngStep = ceil( dAng / dStepRotDeg) ; // aggiusto lo step
+         for ( int i = 1 ; i < dAngStep ; ++ i) {
+            Point3d ptRot = ptRight ;
+            ptRot.Rotate( ptCenter, vtNorm, i * ( dAng / dAngStep)) ;
+            plLoop.AddUPoint( i, ptRot) ; // punto intermedio sulla circonferenza
+         }
+         plLoop.AddUPoint( dAngStep ++, ptLeft) ; // ultimo punto
+         plLoop.AddUPoint( dAngStep ++, ptJunction) ; // punto centrale sull'offset chiuso
+         plLoop.AddUPoint( dAngStep, ptRight) ; // polyLine chiusa
+        // superificie Top
+         PtrOwner<ISurfTriMesh> pStmTop( CreateSurfTriMesh()) ;
+         if ( IsNull( pStmTop)  ||  ! pStmTop->CreateByFlatContour( plLoop))
+            return nullptr ;
+         stmSoup.AddSurfTriMesh( *pStmTop) ;
+        // superificie Bottom
+         PtrOwner<ISurfTriMesh> pStmBottom( CloneSurfTriMesh( pStmTop)) ;
+         if ( IsNull( pStmBottom)  ||  ! pStmBottom->Mirror( ptCen, vtNorm))
+            return nullptr ;
+         stmSoup.AddSurfTriMesh( *pStmBottom) ;
+        // calcolo l'angolo di rotazione per la faccia Top del bevel
+        // NB. Questo angolo va ricalcolato, il bevel è inclinato rispetto alla normale della guida
+         ptCenter.Translate( - dBevelV * vtNorm) ;
+         Point3d ptbRight ; pCrvRb->GetEndPoint( ptbRight) ;
+         Point3d ptbLeft ; pCrvRb->GetStartPoint( ptbLeft) ;
+         Vector3d vtbLeft = ptbLeft - ptCenter ;
+         Vector3d vtbRight = ptbRight - ptCenter ;
+         dAng = ANG_STRAIGHT ;
+         vtbRight.GetAngle( vtbLeft, dAng) ;
+         vtbRight.Normalize() ;
+        // la PolyLine che utilizzo la posso ricavare da quella calcolata sopra
+         plLoop.EraseLastUPoint() ; // apro il loop
+         plLoop.EraseLastUPoint() ; // tolgo il punto di contatto sull'offset
+        // creo il loop defininendo i punti
+         PolyLine plLoopB ;
+         plLoopB.AddUPoint( 0, ptbRight) ;
+         dAngStep = ceil( dAng / dStepRotDeg) ;
+         for ( int i = 1 ; i < dAngStep ; ++ i) {
+            Point3d ptRot = ptbRight ;
+            ptRot.Rotate( ptCenter, vtNorm, i * ( dAng / dAngStep)) ;
+            plLoopB.AddUPoint( i, ptRot) ;
+         }
+         plLoopB.AddUPoint( dAngStep, ptbLeft) ;
+        // creo la superficie Top Bevel
+         PtrOwner<ISurfTriMesh> pStmbTop_start( CreateSurfTriMesh()) ;
+         if ( IsNull( pStmbTop_start)  ||
+              ! pStmbTop_start->CreateByTwoCurves( plLoop, plLoopB, ISurfTriMesh::RLT_MINDIST)  ||
+              ! pStmbTop_start->Invert())
+            return nullptr ;
+         stmSoup.AddSurfTriMesh( *pStmbTop_start) ;
+        // creo la superificie Bottom Bevel
+         PtrOwner<ISurfTriMesh> pStmbBottom_start( CloneSurfTriMesh( pStmbTop_start)) ;
+         if ( IsNull( pStmbBottom_start)  ||  ! pStmbBottom_start->Mirror( ptCen, vtNorm))
+            return nullptr ;
+         stmSoup.AddSurfTriMesh( *pStmbBottom_start) ;
+        // creo la superficie perpendicolare alla guida
+         PolyLine plLoopB1 = plLoopB ;
+         plLoopB1.Mirror( ptCen, vtNorm) ;
+         PtrOwner<ISurfTriMesh> pStmPerp( CreateSurfTriMesh()) ;
+         if ( IsNull( pStmPerp)  ||
+              ! pStmPerp->CreateByTwoCurves( plLoopB, plLoopB1, ISurfTriMesh::RLT_MINDIST)  ||
+              ! pStmPerp->Invert())
+            return nullptr ;
+         stmSoup.AddSurfTriMesh( *pStmPerp) ;
+      }
+     // se l'offset interno della guida è aperto...
+      else {
+        // aggiungo il cap sull'inizio
+         Point3d ptStart ;
+         pGuide->GetStartPoint( ptStart) ;
+        // calcolo l'angolo di rotazione per screwing faccia Top e Bottom
+         Point3d ptSLeft ; pCrvL->GetStartPoint( ptSLeft) ;
+         Point3d ptSRight ; pCrvR->GetStartPoint( ptSRight) ;
+         Vector3d vtLeft = ptSLeft - ptStart ;
+         Vector3d vtRight = ptSRight - ptStart ;
+         double dAng = ANG_STRAIGHT ;
+         vtLeft.GetAngle( vtRight, dAng) ;
+         vtLeft.Normalize() ;
+         PolyLine plLoop ;
+        // creo il loop defininendo i punti
+         plLoop.AddUPoint( 0, ptSLeft) ;  // primo punto
+         double dAngStep = ceil( dAng / dStepRotDeg) ;
+         for ( int i = 1 ; i < dAngStep ; ++ i) {
+            Point3d ptRot = ptSLeft ;
+            ptRot.Rotate( ptStart, vtNorm, i * ( dAng / dAngStep)) ;
+            plLoop.AddUPoint( i, ptRot) ;
+         }
+         plLoop.AddUPoint( dAngStep, ptSRight) ; // ultimo punto
+         plLoop.AddUPoint( dAngStep + 1, ptSLeft) ; // polyline chiusa
+        // creo la superficie Top
+         PtrOwner<ISurfTriMesh> pStmTop_start( CreateSurfTriMesh()) ;
+         if ( IsNull( pStmTop_start)  ||  ! pStmTop_start->CreateByFlatContour( plLoop))
+            return nullptr ;
+         stmSoup.AddSurfTriMesh( *pStmTop_start) ;
+        // creo la superificie Bottom
+         PtrOwner<ISurfTriMesh> pStmBottom_start( CloneSurfTriMesh( pStmTop_start)) ;
+         if ( IsNull( pStmBottom_start)  ||  ! pStmBottom_start->Mirror( ptCen, vtNorm))
+            return nullptr ;
+         stmSoup.AddSurfTriMesh( *pStmBottom_start) ;
+        // calcolo l'angolo di rotazione per la faccia Top del bevel
+         ptStart.Translate( - dBevelV * vtNorm) ;
+         Point3d ptSbLeft ; pCrvLb->GetStartPoint( ptSbLeft) ;
+         Point3d ptSbRight ; pCrvRb->GetStartPoint( ptSbRight) ;
+         Vector3d vtbLeft = ptSbLeft - ptStart ;
+         Vector3d vtbRight = ptSbRight - ptStart ;
+         dAng = ANG_STRAIGHT ;
+         vtbLeft.GetAngle( vtbRight, dAng) ;
+         vtbLeft.Normalize() ;
+         plLoop.EraseLastUPoint() ; // apro il loop
+        // creo il loop defininendo i punti
+         PolyLine plLoopB ;
+         plLoopB.AddUPoint( 0, ptSbLeft) ;
+         dAngStep = ceil( dAng / dStepRotDeg) ; // primo punto
+         for ( int i = 1 ; i < dAngStep ; ++ i) {
+            Point3d ptRot = ptSbLeft ;
+            ptRot.Rotate( ptStart, vtNorm, i * ( dAng / dAngStep)) ;
+            plLoopB.AddUPoint( i, ptRot) ;
+         }
+         plLoopB.AddUPoint( dAngStep, ptSbRight) ; // ultimo punto
+        // creo la superficie Top Bevel
+         PtrOwner<ISurfTriMesh> pStmbTop_start( CreateSurfTriMesh()) ;
+         if ( IsNull( pStmbTop_start)  ||
+              ! pStmbTop_start->CreateByTwoCurves( plLoop, plLoopB, ISurfTriMesh::RLT_MINDIST)  ||
+              ! pStmbTop_start->Invert())
+            return nullptr ;
+         stmSoup.AddSurfTriMesh( *pStmbTop_start) ;
+        // creo la superificie Bottom Bevel
+         PtrOwner<ISurfTriMesh> pStmbBottom_start( CloneSurfTriMesh( pStmbTop_start)) ;
+         if ( IsNull( pStmbBottom_start)  ||  ! pStmbBottom_start->Mirror( ptCen, vtNorm))
+            return nullptr ;
+         stmSoup.AddSurfTriMesh( *pStmbBottom_start) ;
+        // creo la superficie perpendicolare alla guida
+         PolyLine plLoopB1 = plLoopB ;
+         plLoopB1.Mirror( ptCen, vtNorm) ;
+         PtrOwner<ISurfTriMesh> pStmPerp_start( CreateSurfTriMesh()) ;
+         if ( IsNull( pStmPerp_start)  ||
+              ! pStmPerp_start->CreateByTwoCurves( plLoopB, plLoopB1, ISurfTriMesh::RLT_MINDIST)  ||
+              ! pStmPerp_start->Invert())
+            return nullptr ;
+         stmSoup.AddSurfTriMesh( *pStmPerp_start) ;
+         // aggiungo il cap sulla fine
+         Point3d ptEnd ;
+         pGuide->GetEndPoint( ptEnd) ;
+        // calcolo l'angolo di rotazione per screwing faccia Top e Bottom
+         pCrvL->GetEndPoint( ptSLeft) ;
+         pCrvR->GetEndPoint( ptSRight) ;
+         vtLeft = ptSLeft - ptEnd ;
+         vtRight = ptSRight - ptEnd ;
+         dAng = ANG_STRAIGHT ;
+         vtRight.GetAngle( vtLeft, dAng) ;
+         vtRight.Normalize() ;
+         plLoop.Clear() ;
+        // creo il loop defininendo i punti
+         plLoop.AddUPoint( 0, ptSRight) ;
+         dAngStep = ceil( dAng / dStepRotDeg) ; // primo punto
+         for ( int i = 1 ; i < dAngStep ; ++ i) {
+            Point3d ptRot = ptSRight ;
+            ptRot.Rotate( ptEnd, vtNorm, i * ( dAng / dAngStep)) ;
+            plLoop.AddUPoint( i, ptRot) ;
+         }
+         plLoop.AddUPoint( dAngStep, ptSLeft) ; // ultimo punto
+         plLoop.AddUPoint( dAngStep + 1, ptSRight) ; // polyline chiusa
+        // creo la superficie Top
+         PtrOwner<ISurfTriMesh> pStmTop_end( CreateSurfTriMesh()) ;
+         if ( IsNull( pStmTop_end)  ||  ! pStmTop_end->CreateByFlatContour( plLoop))
+            return nullptr ;
+         stmSoup.AddSurfTriMesh( *pStmTop_end) ;
+        // creo la superificie Bottom
+         PtrOwner<ISurfTriMesh> pStmBottom_end( CreateSurfTriMesh()) ;
+        if ( IsNull( pStmBottom_end)  ||
+             ! pStmBottom_end->CopyFrom( pStmTop_end)  ||
+             ! pStmBottom_end->Mirror( ptCen, vtNorm))
+            return nullptr ;
+         stmSoup.AddSurfTriMesh( *pStmBottom_end) ;
+        // calcolo l'angolo di rotazione per la faccia Top del bevel
+         ptEnd.Translate( - dBevelV * vtNorm) ;
+         pCrvLb->GetEndPoint( ptSbLeft) ;
+         pCrvRb->GetEndPoint( ptSbRight) ;
+         vtbLeft = ptSbLeft - ptEnd ;
+         vtbRight = ptSbRight - ptEnd ;
+         dAng = ANG_STRAIGHT ;
+         vtbRight.GetAngle( vtbLeft, dAng) ;
+         vtbRight.Normalize() ;
+         plLoop.EraseLastUPoint() ; // apro il loop
+        // creo il loop defininendo i punti
+         plLoopB.Clear() ;
+         plLoopB.AddUPoint( 0, ptSbRight) ;
+         dAngStep = ceil( dAng / dStepRotDeg) ; // primo punto
+         for ( int i = 1 ; i < dAngStep ; ++ i) {
+            Point3d ptRot = ptSbRight ;
+            ptRot.Rotate( ptEnd, vtNorm, i * ( dAng / dAngStep)) ;
+            plLoopB.AddUPoint( i, ptRot) ;
+         }
+         plLoopB.AddUPoint( dAngStep, ptSbLeft) ; // ultimo punto
+        // creo la superficie Top Bevel
+         PtrOwner<ISurfTriMesh> pStmbTop_end( CreateSurfTriMesh()) ;
+         if ( IsNull( pStmbTop_end)  ||
+              ! pStmbTop_end->CreateByTwoCurves( plLoop, plLoopB, ISurfTriMesh::RLT_MINDIST)  ||
+              ! pStmbTop_end->Invert())
+            return nullptr ;
+         stmSoup.AddSurfTriMesh( *pStmbTop_end) ;
+        // creo la superificie Bottom Bevel
+         PtrOwner<ISurfTriMesh> pStmbBottom_end( CloneSurfTriMesh( pStmbTop_end)) ;
+         if ( IsNull( pStmbBottom_end)  ||  ! pStmbBottom_end->Mirror( ptCen, vtNorm))
+            return nullptr ;
+         stmSoup.AddSurfTriMesh( *pStmbBottom_end) ;
+        // creo la superficie perpendicolare alla guida
+         plLoopB1 = plLoopB ;
+         plLoopB1.Mirror( ptCen, vtNorm) ;
+         PtrOwner<ISurfTriMesh> pStmPerp_end( CreateSurfTriMesh()) ;
+         if ( IsNull( pStmPerp_end)  ||
+              ! pStmPerp_end->CreateByTwoCurves( plLoopB, plLoopB1, ISurfTriMesh::RLT_MINDIST)  ||
+              ! pStmPerp_end->Invert())
          return nullptr ;
-      pSci->Invert() ;
-      stmSoup.AddSurfTriMesh( *pSci) ;
-     // aggiungo il cap sulla fine
-      Point3d ptEnd ;
-      pGuide->GetEndPoint( ptEnd) ;
-      Vector3d vtEnd ;
-      pGuide->GetEndDir( vtEnd) ;
-      vtEnd.Rotate( vtNorm, 0, -1) ;
-      PolyLine PLEnd ;
-      PLEnd.AddUPoint( 0, ptEnd) ;
-      PLEnd.AddUPoint( 1, ptEnd + ( dDimH / 2 - dBevelH) * vtEnd) ;
-      PLEnd.AddUPoint( 2, ptEnd + dDimH / 2 * vtEnd - dBevelV * vtNorm) ;
-      PLEnd.AddUPoint( 3, ptEnd + dDimH / 2 * vtEnd - ( dDimV - dBevelV) * vtNorm) ;
-      PLEnd.AddUPoint( 4, ptEnd + ( dDimH / 2 - dBevelH) * vtEnd - dDimV * vtNorm) ;
-      PLEnd.AddUPoint( 5, ptEnd - dDimV * vtNorm) ;
-      PtrOwner<SurfTriMesh> pSce( CreateBasicSurfTriMesh()) ;
-      if ( IsNull( pSce)  ||  ! pSce->CreateByScrewing( PLEnd, ptEnd, vtNorm, ANG_STRAIGHT, dStepRotDeg, 0))
-         return nullptr ;
-      pSce->Invert() ;
-      stmSoup.AddSurfTriMesh( *pSce) ;
+         stmSoup.AddSurfTriMesh( *pStmPerp_end) ;
+      }
      // completo unione e recupero la superficie risultante
       if ( ! stmSoup.End())
          return nullptr ;
@@ -634,6 +939,7 @@ GetSurfTriMeshBeveledRectSwept( double dDimH, double dDimV, double dBevelH, doub
    return Release( pSTM) ;
 }
 
+
 //-------------------------------------------------------------------------------
 ISurfTriMesh*
 GetSurfTriMeshRectSwept( double dDimH, double dDimV, double dBevelH, double dBevelV, const ICurve* pGuide, int nCapType, double dLinTol)

SurfBezier.h

@@ -159,7 +159,7 @@ class SurfBezier : public ISurfBezier, public IGeoObjRW
    private :
       enum Status { ERR = 0, OK = 1, TO_VERIFY = 2} ;
       enum ParDir { ON_U = 1, ON_V = 2} ;
-      static const int MAXDEG = 11 ;
+      static const int MAXDEG = 21 ;
 
    private :
       bool CopyFrom( const SurfBezier& sbSrc) ;

SurfTriMesh.cpp

@@ -20,12 +20,14 @@
 #include "NgeWriter.h"
 #include "NgeReader.h"
 #include "SurfFlatRegion.h"
-#include "DistPointLine.h"
-#include "DistLineLine.h"
 #include "Triangulate.h"
 #include "GeoConst.h"
+#include "/EgtDev/Include/EGkDistPointLine.h"
+#include "/EgtDev/Include/EGkDistLineLine.h"
+#include "/EgtDev/Include/EGkDistPointSurfTm.h"
 #include "/EgtDev/Include/EGkIntersLinePlane.h"
 #include "/EgtDev/Include/EGkPointGrid3d.h"
+#include "/EgtDev/Include/EGkPolygon3d.h"
 #include "/EgtDev/Include/EGkPolyLine.h"
 #include "/EgtDev/Include/EGkStringUtils3d.h"
 #include "/EgtDev/Include/EGkUiUnits.h"
@@ -44,7 +46,7 @@ SurfTriMesh::SurfTriMesh( void)
    : m_nStatus( TO_VERIFY), m_dLinTol( STM_STD_LIN_TOL), m_dBoundaryAng( STM_STD_BOUNDARY_ANG),
      m_dSmoothAng( STM_STD_SMOOTH_ANG), m_bShowEdges( false), m_bOriented( false), m_bClosed( false),
      m_bFaceted( false), m_bFacEdged( false), m_nTimeStamp( 0), m_nTempProp{0,0}, m_dTempParam{0,0},
-     m_nMaxTFlag( 0), m_nParts( -1), m_pHGrd3d( nullptr)
+     m_nMaxTFlag( 0), m_nParts( -1), m_nShells( -1), m_pHGrd3d( nullptr)
 {
    m_dCosBndAng = cos( m_dBoundaryAng * DEGTORAD) ;
    m_dCosSmAng = cos( m_dSmoothAng * DEGTORAD) ;
@@ -84,6 +86,7 @@ SurfTriMesh::Init( int nNumVert, int nNumTria, int nNumFacet)
    m_bClosed = false ;
    m_nMaxTFlag = 0 ;
    m_nParts = -1 ;
+   m_nShells = -1 ;
 
    return true ;
 }
@@ -115,6 +118,7 @@ SurfTriMesh::Clear( void)
    m_dTempParam[1] = 0 ;
    m_nMaxTFlag = 0 ;
    m_nParts = -1 ;
+   m_nShells = -1 ;
    return true ;
 }
 
@@ -125,6 +129,7 @@ SurfTriMesh::AddVertex( const Point3d& ptVert, double dU, double dV)
    // imposto ricalcolo
    m_nStatus = TO_VERIFY ;
    m_nParts = - 1 ;
+   m_nShells = -1 ;
    m_OGrMgr.Reset() ;
    ResetHashGrids3d() ;
    // inserisco il vertice
@@ -156,6 +161,7 @@ SurfTriMesh::MoveVertex( int nInd, const Point3d& ptNewVert)
   // imposto ricalcolo
    m_nStatus = TO_VERIFY ;
    m_nParts = - 1 ;
+   m_nShells = - 1 ;
    m_bFaceted = false ;
    m_bFacEdged = false ;
    m_OGrMgr.Reset() ;
@@ -225,6 +231,7 @@ SurfTriMesh::AddTriangle( const int nIdVert[3], int nTFlag)
   // imposto ricalcolo
    m_nStatus = TO_VERIFY ;
    m_nParts = - 1 ;
+   m_nShells = -1 ;
    m_OGrMgr.Reset() ;
    ResetHashGrids3d() ;
   // inserisco il triangolo
@@ -316,6 +323,7 @@ SurfTriMesh::RemoveTriangle( int nId)
    m_bFacEdged = false ;
   // invalido calcolo connettività
    m_nParts = - 1 ;
+   m_nShells = -1 ;
    return true ;
 }
 
@@ -905,6 +913,7 @@ SurfTriMesh::CloneTriangle( int nT) const
    pSurfTM->m_dSmoothAng = m_dSmoothAng ;
    pSurfTM->m_dCosSmAng = m_dCosSmAng ;
    pSurfTM->m_nParts = 1 ;
+   pSurfTM->m_nShells = 1 ;
 
   // Copio il triangolo
    int nNewInd[3] = { pSurfTM->AddVertex( m_vVert[m_vTria[nT].nIdVert[0]].ptP),
@@ -1048,10 +1057,16 @@ bool
 SurfTriMesh::MarchAlongLoop( int nT, int nV, int nTimeStamp, PolyLine& PL) const
 {
   // mi muovo lungo il loop, un triangolo alla volta
+   int nCount = 0 ;
    bool bEnd = false ;
    while ( ! bEnd) {
+     // altro triangolo
       if ( ! MarchOneTria( nT, nV, nTimeStamp, PL, bEnd))
          return false ;
+     // per evitare loop infiniti
+      ++ nCount ;
+      if ( nCount > 3 * int( m_vTria.size()) + 10)
+         return true ;
    }
    return true ;
 }
@@ -1171,6 +1186,78 @@ SurfTriMesh::GetSilhouette( const Vector3d& vtDir, double dTol, POLYLINEVECTOR&
    return true ;
 }
 
+//----------------------------------------------------------------------------
+bool
+SurfTriMesh::GetSilhouette( const Plane3d& plPlane, double dTol, POLYLINEVECTOR& vPL, bool bAllTria) const
+{
+  // Verifico lo stato
+   if ( m_nStatus != OK)
+      return false ;
+
+  // Verifico la direzione
+   Vector3d vtVers = plPlane.GetVersN() ;
+   if ( ! vtVers.Normalize())
+      return false ;
+
+  // Controlli su tolleranza
+   dTol = max( dTol, 100 * EPS_SMALL) ;
+
+  // Determino il riferimento di proiezione
+   Frame3d frOCS ; frOCS.Set( plPlane.GetPoint(), vtVers) ;
+
+  // Ottengo la Silhouette come unione delle regioni dei triangoli proiettati (solo parti sopra il piano)
+  // calcolo la regione dei triangoli proiettati
+   PtrOwner<SurfFlatRegion> pSfr ;
+   Triangle3d Tria ;
+   int nT = GetFirstTriangle( Tria) ;
+   while ( nT != SVT_NULL) {
+     // verifico la normale
+      if ( ( bAllTria  &&  abs( Tria.GetN() * vtVers) > EPS_ZERO)  ||
+           ( ! bAllTria  &&  Tria.GetN() * vtVers > EPS_ZERO)) {
+        // ricavo il poligono equivalente al triangolo
+         Polygon3d pgTria ;
+         pgTria.FromTriangle( Tria) ;
+        // taglio il poligono con il piano
+         pgTria.Trim( plPlane, false, true, bAllTria) ;
+        // se rimasto qualcosa
+         if ( pgTria.GetSideCount() > 0) {
+           // lo proietto sul piano e creo la regione
+            pgTria.Scale( frOCS, 1, 1, 0) ;
+            PtrOwner<SurfFlatRegion> pSfrTria( GetBasicSurfFlatRegion( GetSurfFlatRegionFromPolyLine( pgTria.GetPolyLine()))) ;
+            if ( ! IsNull( pSfrTria)) {
+               if ( bAllTria  &&  Tria.GetN() * vtVers < 0)
+                  pSfrTria->Invert() ;
+               pSfrTria->Offset( dTol, ICurve::OFF_FILLET) ;
+               if ( IsNull( pSfr))
+                  pSfr.Set( pSfrTria) ;
+               else
+                  pSfr->Add( *pSfrTria) ;
+            }
+         }
+      }
+     // passo al successivo
+      nT = GetNextTriangle( nT, Tria) ;
+   }
+
+  // Se non esiste la regione
+   if ( IsNull( pSfr))
+      return false ;
+
+  // Effettuo contro-offset
+   pSfr->Offset( -dTol, ICurve::OFF_EXTEND) ;
+
+  // Recupero i contorni della regione
+   for ( int i = 0 ; i < pSfr->GetChunkCount() ; ++ i) {
+      for ( int j = 0 ; j < pSfr->GetLoopCount( i) ; ++ j) {
+         PolyLine PL ;
+         if ( pSfr->ApproxLoopWithLines( i, j, LIN_TOL_STD, ANG_TOL_STD_DEG, ICurve::APL_STD, PL))
+            vPL.emplace_back( PL) ;
+      }
+   }
+
+   return true ;
+}
+
 //----------------------------------------------------------------------------
 SurfTriMesh*
 SurfTriMesh::Clone( void) const
@@ -1225,6 +1312,7 @@ SurfTriMesh::CopyFrom( const SurfTriMesh& stmSrc)
    m_nTempProp[1] = stmSrc.m_nTempProp[1] ;
    m_nMaxTFlag = stmSrc.m_nMaxTFlag ;
    m_nParts = stmSrc.m_nParts ;
+   m_nShells = stmSrc.m_nShells ;
    m_dTempParam[0] = stmSrc.m_dTempParam[0] ;
    m_dTempParam[1] = stmSrc.m_dTempParam[1] ;
    return true ;
@@ -1262,10 +1350,12 @@ SurfTriMesh::Dump( string& sOut, bool bMM, const char* szNewLine) const
   // segnalo eventuale incongruenza di orientamento
    if ( ! m_bOriented)
       sOut += string( "Inconsistent Orientation") + szNewLine ;
-  // segnalo numero di parti se più di una
+  // segnalo numero di parti
    int nParts = GetPartCount() ;
-   if ( nParts > 1)
-      sOut += string( "Parts =") + ToString( nParts) + szNewLine ;
+   sOut += string( "Parts=") + ToString( nParts) + szNewLine ;
+  // segnalo numero di gusci
+   int nShells = GetShellCount() ;
+   sOut += string( "Shells=") + ToString( nShells) + szNewLine ;
   // numero di vertici
    sOut += "Vert : Nbr=" + ToString( GetVertexCount()) +
                 " Size=" + ToString( GetVertexSize()) + szNewLine ;
@@ -1376,6 +1466,7 @@ SurfTriMesh::Load( NgeReader& ngeIn)
    m_OGrMgr.Clear() ;
    m_nMaxTFlag = 0 ;
    m_nParts = -1 ;
+   m_nShells = -1 ;
    ResetHashGrids3d() ;
   // leggo la prossima linea ( 2 parametri : dLinTol e dSmoothAng)
   // tolleranza lineare di costruzione
@@ -1536,6 +1627,7 @@ SurfTriMesh::Validate( bool bCorrect)
 
   // invalido calcolo connessione
    m_nParts = - 1 ;
+   m_nShells = -1 ;
 
    return ( m_nStatus == OK) ;
 }
@@ -1811,6 +1903,7 @@ SurfTriMesh::AdjustTopology( void)
    m_bFacEdged = false ;
   // invalido calcolo connessione
    m_nParts = - 1 ;
+   m_nShells = - 1 ;
   // verifica indici
    if ( ! Validate( true))
       return false ;
@@ -3562,6 +3655,17 @@ SurfTriMesh::GetPartCount( void) const
    return m_nParts ;
 }
 
+//----------------------------------------------------------------------------
+int
+SurfTriMesh::GetShellCount( void) const
+{
+   if ( ! IsValid())
+      return 0 ;
+   if ( m_nShells == - 1  &&  ! VerifyConnection())
+      return 0 ;
+   return m_nShells ;
+}
+
 //----------------------------------------------------------------------------
 bool
 SurfTriMesh::VerifyConnection( void) const
@@ -3570,19 +3674,21 @@ SurfTriMesh::VerifyConnection( void) const
       return false ;
 
   // reset connessione
-   for ( auto Tria : m_vTria)
+   for ( auto& Tria : m_vTria) {
       Tria.nPart = SVT_NULL ;
+      Tria.nShell = SVT_NULL ;
+   }
 
-  // ciclo sui triangoli
-   m_nParts = 0 ;
+  // ciclo sui triangoli per determinare le shells
+   m_nShells = 0 ;
    for ( int i = 0 ; i < int( m_vTria.size()) ; ++ i) {
-     // salto triangoli cancellati o gi� assegnati
+     // salto triangoli cancellati o già assegnati
       if ( m_vTria[i].nIdVert[0] == SVT_DEL  ||
-           m_vTria[i].nPart != SVT_NULL)
+           m_vTria[i].nShell != SVT_NULL)
          continue ;
      // assegno indice di parte connessa al triangolo
-      m_vTria[i].nPart = m_nParts ;
-      ++ m_nParts ;
+      m_vTria[i].nShell = m_nShells ;
+      ++ m_nShells ;
      // set di triangoli da aggiornare
       set<int> stTria ;
       stTria.insert( i) ;
@@ -3594,13 +3700,137 @@ SurfTriMesh::VerifyConnection( void) const
         // aggiorno i triangoli adiacenti
          for ( int j = 0 ; j < 3 ; ++ j) {
             int nAdjT = m_vTria[nT].nIdAdjac[j] ;
-            if ( nAdjT != SVT_NULL  &&  m_vTria[nAdjT].nPart == SVT_NULL) {
-               m_vTria[nAdjT].nPart = m_vTria[nT].nPart ;
+            if ( nAdjT != SVT_NULL  &&  m_vTria[nAdjT].nShell == SVT_NULL) {
+               m_vTria[nAdjT].nShell = m_vTria[nT].nShell ;
                stTria.insert( nAdjT) ;
             }
          }
       }
    }
+
+  // identifico le parti
+   m_nParts = 0 ;
+   if ( m_nShells == 1) { // se ho solo una shell allora ho una sola parte
+      m_nParts = 1 ;
+      for ( auto& Tria : m_vTria)
+         Tria.nPart = Tria.nShell ;
+   }
+   else { // se ho più shells devo controllare la loro posizione
+     // salvo solo le Shells chiuse e per ognuna di esse il loro volume e l'indice della shell
+     // < stmShell, ( nShellInd, dVol) >
+     // NB. I Volumi sono sempre finiti, il segno del volume mi dice solo se sto considerando
+     // la parte interna o esterna della superficie
+      typedef pair<PtrOwner<ISurfTriMesh>, pair<int, double>> SHELLINFO ;
+      vector<SHELLINFO> vClosedShells ;
+      vClosedShells.reserve( m_nShells) ;
+      INTVECTOR vOpenShells ;    // vettore di indici di Shell aperte
+      for ( int s = 0 ; s < GetShellCount() ; ++ s) {
+        // recupero la Shell
+         PtrOwner<ISurfTriMesh> pStmShell( CloneShell( s)) ;
+         if ( IsNull( pStmShell)  ||  ! pStmShell->IsValid())
+            return false ;
+         double dVol = 0. ;
+        // se chiusa calcolo il suo volume con segno
+         if ( pStmShell->IsClosed()) {
+            pStmShell->GetVolume( dVol) ;
+            vClosedShells.emplace_back( make_pair( Release( pStmShell), make_pair( s, dVol))) ;
+         }
+        // se aperta, non rientra nella classficazione ( sarà coincidente ad una Part)
+         else
+            vOpenShells.push_back( s) ;
+      }
+
+     // ordino il vettore in senso decrescente per volume, in valore assoluto
+     // NB. tutte le shell con volume negativo sono state invertite !
+      sort( vClosedShells.begin(), vClosedShells.end(),
+            []( const SHELLINFO& a, const SHELLINFO& b) {
+               return abs( a.second.second) > abs( b.second.second) ;
+            }) ;
+
+     // se la superficie in prima posizione è negativa, allora le inverto tutte; il numero di
+     // Parts non cambia nel conto
+      if ( ! vClosedShells.empty()  &&  vClosedShells.front().second.second < - EPS_SMALL)
+         for ( auto& Shell : vClosedShells)
+            Shell.first->Invert() ;
+
+     // creo una matrice di interi ; ogni riga corrisponde ad una Parte, dove in posizione 0 c'è
+     // la shell esterna e nelle successive le rispettiva shells interne
+      INTMATRIX vnShellIndMat ;
+
+      bool bFirstShell ;         // flag per indicare se la shell corrente è esterna
+      int nIndExtShell = -1 ;    // indice del vettore di shell chiuse della shell esterna
+      INTVECTOR vIndIntShells ;  // vettore di indici di shell chiuse interne alla shell esterna corrente
+      do {
+         bFirstShell = true ;
+         for ( int i = 0 ; i < int( vClosedShells.size()) ; ++ i) {
+           // recupero l'indice della shell e verifico che sia valido
+            int j = vClosedShells[i].second.first ;
+            if ( j < 0)
+               continue ;
+           // lo inserisco come esterno...
+            if ( bFirstShell) {
+               vnShellIndMat.push_back({ j}) ;
+               vClosedShells[i].second.first = -1 ;
+               bFirstShell = false ;
+               nIndExtShell = i ;
+               vIndIntShells.clear() ;
+            }
+           // altrimento verifico se la shell è interna o no
+            else {
+              // la shell è interna se è sia interna alla shell esterna corrente ( della riga di
+              // vnShellIndMat ) e allo stesso tempo esterna a tutte le shell già inserite nella riga
+              // attuale.
+               Point3d ptCheck ;
+               vClosedShells[i].first->GetFirstVertex( ptCheck) ;
+
+              // 1) verifica rispetto alla shell esterna
+               DistPointSurfTm distCalculator( ptCheck, *vClosedShells[nIndExtShell].first) ;
+               if ( ! distCalculator.IsPointInside())
+                  continue ;
+
+              // 2) verifica rispetto a tutte la shell già trovate interne
+               bool bOk = true ;
+               for ( int k = 1 ; k < int( vIndIntShells.size()) ; ++ k) {
+                  DistPointSurfTm distCalculator( ptCheck, *vClosedShells[vIndIntShells[k]].first) ;
+                 // NB. tutte le shell hanno volume negativo !
+                  if ( ! distCalculator.IsPointInside()) {
+                     bOk = false ;
+                     break ;
+                  }
+               }
+               if ( bOk) {
+                 // inserisco nella matrice
+                  vnShellIndMat.back().push_back( j) ;
+                  vClosedShells[i].second.first = - 1 ;
+                  vIndIntShells.push_back( i) ;
+               }
+            }
+         }
+      } while ( ! bFirstShell) ;
+
+     // scorro le righe della matrice e assegno la Part ( ovvero la riga )
+      int nPart = 0 ;
+      for ( nPart = 0 ; nPart < int( vnShellIndMat.size()) ; ++ nPart) {
+         ++ m_nParts ;
+         for ( int j = 0 ; j < int( vnShellIndMat[nPart].size()) ; ++ j) {
+            for ( auto& Tria : m_vTria) {
+               if ( Tria.nShell == j)
+                  Tria.nPart = nPart ;
+            }
+         }
+      }
+
+     // per tutte le superfici aperte, assegno loro una Parte
+      for ( int i = 0 ; i < int( vOpenShells.size()) ; ++ i) {
+         ++ m_nParts ;
+         for ( auto& Tria : m_vTria) {
+            if ( Tria.nShell == i)
+               Tria.nPart = i + ( nPart ++ ) ;
+         }
+      }
+
+   }
+
    return true ;
 }
 
@@ -3636,6 +3866,38 @@ SurfTriMesh::RemovePart( int nPart)
    return true ;
 }
 
+//----------------------------------------------------------------------------
+bool
+SurfTriMesh::RemoveShell( int nShell)
+{
+   if ( ! IsValid())
+      return false ;
+
+  // Il numero delle componenti deve essere maggiore di zero o calcolabile
+   if ( m_nShells == -1  &&  ! VerifyConnection())
+      return false ;
+
+  // Se la componente non esiste, errore
+   if ( m_nShells <= 0  ||  nShell < 0  ||  nShell >= m_nShells)
+      return false ;
+
+   int nShellsOld = m_nShells ;
+  // Rimuovo i triangoli della componente nPart
+   for ( int i = 0 ; i < int( m_vTria.size()) ; ++ i) {
+      if ( m_vTria[i].nShell == nShell)
+         RemoveTriangle( i) ;
+   }
+
+  // Aggiorno il numero di componenti
+   m_nShells = nShellsOld - 1 ;
+
+  // imposto ricalcolo della grafica e di hashgrids3d
+   m_OGrMgr.Reset() ;
+   ResetHashGrids3d() ;
+
+   return true ;
+}
+
 //----------------------------------------------------------------------------
 SurfTriMesh*
 SurfTriMesh::ClonePart( int nPart) const
@@ -3663,6 +3925,7 @@ SurfTriMesh::ClonePart( int nPart) const
    pSurfTM->m_dSmoothAng = m_dSmoothAng ;
    pSurfTM->m_dCosSmAng = m_dCosSmAng ;
    pSurfTM->m_nParts = 1 ;
+   pSurfTM->m_nShells = 1 ;
 
   // Copio i triangoli
    for ( int i = 0 ; i < int( m_vTria.size()) ; ++ i) {
@@ -3682,6 +3945,53 @@ SurfTriMesh::ClonePart( int nPart) const
    return Release( pSurfTM) ;
 }
 
+//----------------------------------------------------------------------------
+SurfTriMesh* 
+SurfTriMesh::CloneShell( int nShell) const
+{
+   if ( ! IsValid())
+      return nullptr ;
+
+  // Il numero delle componenti deve essere maggiore di zero o calcolabile
+   if ( m_nShells == -1  &&  ! VerifyConnection())
+      return nullptr ;
+
+  // Se la componente non esiste, errore
+   if ( m_nShells <= 0  ||  nShell < 0  ||  nShell >= m_nShells)
+      return nullptr ;
+
+  // Creo nuovo oggetto SurfTriMesh
+   PtrOwner<SurfTriMesh> pSurfTM( new( nothrow) SurfTriMesh) ;
+   if ( IsNull( pSurfTM))
+      return nullptr ;
+
+  // Copio il valore dei membri
+   pSurfTM->m_dLinTol = m_dLinTol ;
+   pSurfTM->m_dBoundaryAng = m_dBoundaryAng ;
+   pSurfTM->m_dCosBndAng = m_dCosBndAng ;
+   pSurfTM->m_dSmoothAng = m_dSmoothAng ;
+   pSurfTM->m_dCosSmAng = m_dCosSmAng ;
+   pSurfTM->m_nParts = 1 ;
+   pSurfTM->m_nShells = 1 ;
+
+  // Copio i triangoli
+   for ( int i = 0 ; i < int( m_vTria.size()) ; ++ i) {
+      if ( m_vTria[i].nIdVert[0] != SVT_DEL  &&  m_vTria[i].nShell == nShell) {
+         int nNewInd[3] = { pSurfTM->AddVertex( m_vVert[m_vTria[i].nIdVert[0]].ptP),
+                            pSurfTM->AddVertex( m_vVert[m_vTria[i].nIdVert[1]].ptP),
+                            pSurfTM->AddVertex( m_vVert[m_vTria[i].nIdVert[2]].ptP)} ;
+         if ( pSurfTM->AddTriangle( nNewInd, m_vTria[i].nTFlag) == SVT_NULL)
+            return nullptr ;
+      }
+   }
+
+  // Aggiusto la superficie
+   pSurfTM->DoCompacting() ;
+
+  // Restituisco la nuova superficie
+   return Release( pSurfTM) ;
+}
+
 //----------------------------------------------------------------------------
 bool
 SurfTriMesh::ResetTFlags( void)
@@ -3710,6 +4020,8 @@ SurfTriMesh::SetTFlag( int nId, int nTFlag)
    if ( nId < 0  ||  nId >= GetTriangleSize()  ||  m_vTria[nId].nIdVert[0] == SVT_DEL)
       return false ;
    m_vTria[nId].nTFlag = nTFlag ;
+   m_nMaxTFlag = max( m_nMaxTFlag, abs( nTFlag)) ;
+   m_OGrMgr.Clear() ;
    return true ;
 }
 

SurfTriMesh.h

@@ -52,16 +52,16 @@ class StmTria
    public :
       StmTria( void)
          : nIdVert{ SVT_NULL, SVT_NULL, SVT_NULL}, nIdAdjac{ SVT_NULL, SVT_NULL, SVT_NULL}, nETempFlag{ 0, 0, 0},
-           vtN(), nIdFacet( SVT_NULL), nTFlag( 0), nEFlag( 0), nPart( SVT_NULL), nTemp( 0), nTempPart{ 0} {}
+         vtN(), nIdFacet( SVT_NULL), nTFlag( 0), nEFlag( 0), nShell( SVT_NULL), nPart( SVT_NULL), nTemp( 0), nTempShell{ 0}, nTempPart{ 0} {}
       StmTria( const int nIdV[3])
          : nIdVert{ nIdV[0], nIdV[1], nIdV[2]}, nIdAdjac{ SVT_NULL, SVT_NULL, SVT_NULL}, nETempFlag{ 0, 0, 0},
-           vtN(), nIdFacet( SVT_NULL), nTFlag( 0), nEFlag( 0), nPart( SVT_NULL), nTemp( 0), nTempPart{ 0} {}
+           vtN(), nIdFacet( SVT_NULL), nTFlag( 0), nEFlag( 0), nShell( SVT_NULL), nPart( SVT_NULL), nTemp( 0), nTempShell{ 0}, nTempPart{ 0} {}
       StmTria( const int nIdV[3], int nTF)
          : nIdVert{ nIdV[0], nIdV[1], nIdV[2]}, nIdAdjac{ SVT_NULL, SVT_NULL, SVT_NULL}, nETempFlag{ 0, 0, 0},
-           vtN(), nIdFacet( SVT_NULL), nTFlag( nTF), nEFlag( 0), nPart( SVT_NULL), nTemp( 0), nTempPart{ 0} {}
+           vtN(), nIdFacet( SVT_NULL), nTFlag( nTF), nEFlag( 0), nShell( SVT_NULL), nPart( SVT_NULL), nTemp( 0), nTempShell{ 0}, nTempPart{ 0} {}
       StmTria( const int nIdV[3], const int nIdA[3], const Vector3d& vtV, int nTF, int nEF)
          : nIdVert{ nIdV[0], nIdV[1], nIdV[2]}, nIdAdjac{ nIdA[0], nIdA[1], nIdA[2]}, nETempFlag{ 0, 0, 0},
-           vtN( vtV), nIdFacet( SVT_NULL), nTFlag( nTF), nEFlag( nEF), nPart( SVT_NULL), nTemp( 0), nTempPart{ 0} {}
+           vtN( vtV), nIdFacet( SVT_NULL), nTFlag( nTF), nEFlag( nEF), nShell( SVT_NULL), nPart( SVT_NULL), nTemp( 0), nTempShell{ 0}, nTempPart{ 0} {}
    public :
       int nIdVert[3] ;
       int nIdAdjac[3] ;
@@ -70,9 +70,12 @@ class StmTria
       int nIdFacet ;
       int nTFlag ;
       int nEFlag ;
+      mutable int nShell ;
       mutable int nPart ;
       mutable int nTemp ;
+      mutable int nTempShell ;
       mutable int nTempPart ;
+
 } ;
 
 //----------------------------------------------------------------------------
@@ -279,6 +282,7 @@ class SurfTriMesh : public ISurfTriMesh, public IGeoObjRW
       SurfTriMesh* CloneTriangle( int nT) const override ;
       bool GetLoops( POLYLINEVECTOR& vPL) const override ;
       bool GetSilhouette( const Vector3d& vtDir, double dTol, POLYLINEVECTOR& vPL, bool bAllTria = false) const override ;
+      bool GetSilhouette( const Plane3d& plPlane, double dTol, POLYLINEVECTOR& vPL, bool bAllTria = false) const override ;
       int  GetFacetCount( void) const override ;
       int  GetFacetSize( void) const override
               { return int( m_vFacet.size()) ; }
@@ -316,8 +320,13 @@ class SurfTriMesh : public ISurfTriMesh, public IGeoObjRW
       bool GetAllTriaOverlapBox( const BBox3d& b3Box, INTVECTOR& vT) const override ;
       const BBox3d& GetAllTriaBox( void) const override ;
       int GetPartCount( void) const override ;
+      int GetShellCount( void) const override ;
       bool RemovePart( int nPart) override ;
+      bool RemoveShell( int nShell) override ;
       SurfTriMesh* ClonePart( int nPart) const override ;
+      SurfTriMesh* CloneShell( int nShell) const override ;
+      bool SetTFlag( int nId, int nTFlag) override ;
+      bool GetTFlag( int nId, int& nFlag) const override ;
       int  GetMaxTFlag( void) const override 
               { return m_nMaxTFlag ; }
       bool ResetTFlags( void) override ;
@@ -343,10 +352,8 @@ class SurfTriMesh : public ISurfTriMesh, public IGeoObjRW
       bool ExistsTriangle( int nT) const
               { return ( nT >= 0  &&  nT < GetTriangleSize()  &&  m_vTria[nT].nIdVert[0] != SVT_DEL) ; }
       bool GetTriangleAdjacencies( int nId, int nIdAdjTriaId[3]) const ;
-      bool GetTFlag( int nId, int& nFlag) const ;
       bool GetTempInt( int nId, int& nTempInt) const ;
       bool ResetTempInts( void) const ;
-      bool SetTFlag( int nId, int nTFlag) ;
       bool SetTempInt( int nId, int nTempInt) const ;
 
    private :
@@ -404,6 +411,7 @@ class SurfTriMesh : public ISurfTriMesh, public IGeoObjRW
       bool AmbiguosTriangleManager( TRIA3DVECTORMAP& Ambiguos, SurfTriMesh& Surf) ;
       bool IntersectTriMeshTriangle( SurfTriMesh& Other) ;
       bool IdentifyParts( void) const ;
+      bool IdentifyShells( void) const ;
       bool RemoveDoubleTriangles( bool& bModified) ;
       bool RemoveTJunctions( bool& bModified) ;
       bool FlipTriangles( int nTA, int nTB) ;
@@ -434,9 +442,10 @@ class SurfTriMesh : public ISurfTriMesh, public IGeoObjRW
       int            m_nTempProp[2] ;  // vettore proprietà temporanee
       double         m_dTempParam[2] ; // vettore parametri temporanei
       int            m_nMaxTFlag ;     // massimo valore dei TFlag dei triangoli
-      mutable int    m_nParts ;        // numero di parti connesse (-1 se da calcolare)
+      mutable int    m_nParts ;        // numero di parti intese come volumi (-1 se da calcolare)
+      mutable int    m_nShells ;       // numero di gusci connessi (-1 se da calcolare)
       mutable HashGrids3d* m_pHGrd3d ; // Hash Grid 3d nel suo riferimento
-      mutable BBox3d m_b3HGrd3d ;      // Box3d collegato a Hash Grid 3d 
+      mutable BBox3d m_b3HGrd3d ;      // Box3d collegato a Hash Grid 3d
 } ;
 
 //-----------------------------------------------------------------------------

SurfTriMeshBooleans.cpp

@@ -17,13 +17,14 @@
 #include "CurveLine.h"
 #include "CurveComposite.h"
 #include "SurfFlatRegion.h"
-#include "DistPointLine.h"
 #include "Triangulate.h"
 #include "GeoConst.h"
 #include "/EgtDev/Include/EgtNumUtils.h"
 #include "/EgtDev/Include/EGkCurve.h"
+#include "/EgtDev/Include/EGkDistPointLine.h"
 #include "/EgtDev/Include/EGkDistPointCurve.h"
 #include "/EgtDev/Include/EGkDistPointTria.h"
+#include "/EgtDev/Include/EGkDistPointSurfTm.h"
 #include "/EgtDev/Include/EGkIntersLineTria.h"
 #include "/EgtDev/Include/EGkIntersLineBox.h"
 #include "/EgtDev/Include/EGkIntersPlanePlane.h"
@@ -477,7 +478,7 @@ SurfTriMesh::RetriangulationForBooleanOperation( CHAINMAP& LoopLines, TRIA3DVECT
                   Surf.m_vTria[nNewAloneTriaNum].nETempFlag[0] = 0 ;
                   Surf.m_vTria[nNewAloneTriaNum].nETempFlag[1] = 0 ;
                   Surf.m_vTria[nNewAloneTriaNum].nETempFlag[2] = 0 ;
-                  Surf.m_vTria[nNewAloneTriaNum].nTempPart = nDist[nAloneVert] ;
+                  Surf.m_vTria[nNewAloneTriaNum].nTempShell = nDist[nAloneVert] ;
                   bModif = true ;
                }
                int nNewCoupleId1[3] = { Surf.AddVertex( ptIntSt), Surf.AddVertex( trTria.GetP( ( nAloneVert + 1) % 3)), Surf.AddVertex( trTria.GetP( ( nAloneVert + 2) % 3)) } ;
@@ -486,7 +487,7 @@ SurfTriMesh::RetriangulationForBooleanOperation( CHAINMAP& LoopLines, TRIA3DVECT
                   Surf.m_vTria[nNewCoupleTriaNum1].nETempFlag[0] = 0 ;
                   Surf.m_vTria[nNewCoupleTriaNum1].nETempFlag[1] = 0 ;
                   Surf.m_vTria[nNewCoupleTriaNum1].nETempFlag[2] = 0 ;
-                  Surf.m_vTria[nNewCoupleTriaNum1].nTempPart = - nDist[nAloneVert] ;
+                  Surf.m_vTria[nNewCoupleTriaNum1].nTempShell = - nDist[nAloneVert] ;
                   bModif = true ;
                }
                int nNewCoupleId2[3] = { Surf.AddVertex( ptIntSt), Surf.AddVertex( trTria.GetP( ( nAloneVert + 2) % 3)), Surf.AddVertex( ptIntEn) } ;
@@ -495,7 +496,7 @@ SurfTriMesh::RetriangulationForBooleanOperation( CHAINMAP& LoopLines, TRIA3DVECT
                   Surf.m_vTria[nNewCoupleTriaNum2].nETempFlag[0] = 0 ;
                   Surf.m_vTria[nNewCoupleTriaNum2].nETempFlag[1] = 0 ;
                   Surf.m_vTria[nNewCoupleTriaNum2].nETempFlag[2] = 0 ;
-                  Surf.m_vTria[nNewCoupleTriaNum2].nTempPart = - nDist[nAloneVert] ;
+                  Surf.m_vTria[nNewCoupleTriaNum2].nTempShell = - nDist[nAloneVert] ;
                   bModif = true ;
                }
                continue ;
@@ -561,18 +562,18 @@ SurfTriMesh::RetriangulationForBooleanOperation( CHAINMAP& LoopLines, TRIA3DVECT
                vnDegVec.resize( 1) ;
             }
             else {
-               Surf.m_vTria[it->first].nTempPart = 0 ;
+               Surf.m_vTria[it->first].nTempShell = 0 ;
                continue ;
             }
             vbInOut.resize( 1) ;
             vbInOut[0] = true ;
             if ( ! DecomposeLoop( cvOpenChain, vnDegVec, cvBoundClosedLoopVec, vbInOut)) {
-               Surf.m_vTria[it->first].nTempPart = 0 ;
+               Surf.m_vTria[it->first].nTempShell = 0 ;
                continue ;
             }
          }
          else {
-            Surf.m_vTria[it->first].nTempPart = 0 ;
+            Surf.m_vTria[it->first].nTempShell = 0 ;
             continue ;
          }
       }
@@ -730,9 +731,9 @@ SurfTriMesh::RetriangulationForBooleanOperation( CHAINMAP& LoopLines, TRIA3DVECT
                      Surf.m_vTria[nNewTriaNum].nETempFlag[1] = 0 ;
                      Surf.m_vTria[nNewTriaNum].nETempFlag[2] = 0 ;
                      if ( vbInOut[nLoop])
-                        Surf.m_vTria[nNewTriaNum].nTempPart = 1 ;
+                        Surf.m_vTria[nNewTriaNum].nTempShell = 1 ;
                      else
-                        Surf.m_vTria[nNewTriaNum].nTempPart = - 1 ;
+                        Surf.m_vTria[nNewTriaNum].nTempShell = - 1 ;
                      bModif = true ;
                   }
                }
@@ -827,9 +828,9 @@ SurfTriMesh::RetriangulationForBooleanOperation( CHAINMAP& LoopLines, TRIA3DVECT
                      Surf.m_vTria[nNewTriaNum].nETempFlag[1] = 0 ;
                      Surf.m_vTria[nNewTriaNum].nETempFlag[2] = 0 ;
                      if ( bCodirectedNormals)
-                        Surf.m_vTria[nNewTriaNum].nTempPart = -1 ;
+                        Surf.m_vTria[nNewTriaNum].nTempShell = -1 ;
                      else
-                        Surf.m_vTria[nNewTriaNum].nTempPart = 1 ;
+                        Surf.m_vTria[nNewTriaNum].nTempShell = 1 ;
                      bModif = true ;
                   }
                }
@@ -982,9 +983,9 @@ SurfTriMesh::RetriangulationForBooleanOperation( CHAINMAP& LoopLines, TRIA3DVECT
                      Surf.m_vTria[nNewTriaNum].nETempFlag[1] = 0 ;
                      Surf.m_vTria[nNewTriaNum].nETempFlag[2] = 0 ;
                      if ( bCordirectedNormals_intLoop)
-                        Surf.m_vTria[nNewTriaNum].nTempPart = 1 ;
+                        Surf.m_vTria[nNewTriaNum].nTempShell = 1 ;
                      else
-                        Surf.m_vTria[nNewTriaNum].nTempPart = -1 ;
+                        Surf.m_vTria[nNewTriaNum].nTempShell = -1 ;
                      bModif = true ;
                   }
                }
@@ -1037,7 +1038,7 @@ SurfTriMesh::AmbiguosTriangleManager( TRIA3DVECTORMAP& Ambiguos, SurfTriMesh& Su
    for ( auto it = Ambiguos.begin() ; it != Ambiguos.end() ; ++ it) {
      // Se il triangolo ha l'indice diverso da zero vuol dire che oltre a un
      // contatto edge-edge ha avuto dei contatti che lo hanno già classificato.
-      if ( Surf.m_vTria[it->first].nTempPart != 0)
+      if ( Surf.m_vTria[it->first].nTempShell != 0)
          continue ;
      // Recupero il triangolo corrente
       Triangle3d trTria ;
@@ -1101,11 +1102,11 @@ SurfTriMesh::AmbiguosTriangleManager( TRIA3DVECTORMAP& Ambiguos, SurfTriMesh& Su
       for ( int i = 0 ; i < 3 ; ++ i) {
          if ( nTriaClassificationByEdges[i] == 0)
             continue ;
-         Surf.m_vTria[it->first].nTempPart = nTriaClassificationByEdges[i] ;
+         Surf.m_vTria[it->first].nTempShell = nTriaClassificationByEdges[i] ;
          int j ;
          for ( j = i + 1 ; j < 3 ; ++ j) {
             if ( nTriaClassificationByEdges[j] != 0  &&  nTriaClassificationByEdges[i] != nTriaClassificationByEdges[j]) {
-               Surf.m_vTria[it->first].nTempPart = 0 ;
+               Surf.m_vTria[it->first].nTempShell = 0 ;
                break ;
             }
          }
@@ -1145,13 +1146,13 @@ SurfTriMesh::IntersectTriMeshTriangle( SurfTriMesh& Other)
 
   // Setto il triangolo come né fuori né dentro
    for ( int nTA = 0 ; nTA < nTriaNumA ; ++ nTA) {
-      m_vTria[nTA].nTempPart = 0 ;
+      m_vTria[nTA].nTempShell = 0 ;
       m_vTria[nTA].nETempFlag[0] = 0 ;
       m_vTria[nTA].nETempFlag[1] = 0 ;
       m_vTria[nTA].nETempFlag[2] = 0 ;
    }
    for ( int nTB = 0 ; nTB < nTriaNumB ; ++ nTB) {
-      SurfB.m_vTria[nTB].nTempPart = 0 ;
+      SurfB.m_vTria[nTB].nTempShell = 0 ;
       SurfB.m_vTria[nTB].nETempFlag[0] = 0 ;
       SurfB.m_vTria[nTB].nETempFlag[1] = 0 ;
       SurfB.m_vTria[nTB].nETempFlag[2] = 0 ;
@@ -1259,9 +1260,9 @@ SurfTriMesh::IntersectTriMeshTriangle( SurfTriMesh& Other)
                           abs( ( trTriaA.GetP( ( nVA + 1) % 3) - trTriaB.GetP( 0)) * trTriaB.GetN()) < EPS_SMALL)
                         break ;
                   }
-                  m_vTria[nTA].nTempPart = ( ( trTriaA.GetP( nSegMaxDist) - trTriaB.GetP( 0)) * trTriaB.GetN() < - EPS_SMALL ? 1 : - 1) ;
+                  m_vTria[nTA].nTempShell = ( ( trTriaA.GetP( nSegMaxDist) - trTriaB.GetP( 0)) * trTriaB.GetN() < - EPS_SMALL ? 1 : - 1) ;
                   if ( nVA >= 0  && nVA <= 2)
-                     m_vTria[nTA].nETempFlag[nVA] = m_vTria[nTA].nTempPart ;
+                     m_vTria[nTA].nETempFlag[nVA] = m_vTria[nTA].nTempShell ;
                }
             }
            // caso di Intersezione interno-edge ( per superficie A con B)
@@ -1284,9 +1285,9 @@ SurfTriMesh::IntersectTriMeshTriangle( SurfTriMesh& Other)
                           abs( ( trTriaB.GetP( ( nVB + 1) % 3) - trTriaA.GetP( 0)) * trTriaA.GetN()) < EPS_SMALL)
                         break ;
                   }
-                  SurfB.m_vTria[vNearTria[nTB]].nTempPart = ( ( trTriaB.GetP( nSegMaxDist) - trTriaA.GetP( 0)) * trTriaA.GetN() < - EPS_SMALL ? 1 : - 1) ;
+                  SurfB.m_vTria[vNearTria[nTB]].nTempShell = ( ( trTriaB.GetP( nSegMaxDist) - trTriaA.GetP( 0)) * trTriaA.GetN() < - EPS_SMALL ? 1 : - 1) ;
                   if ( nVB >= 0  &&  nVB <= 2)
-                     SurfB.m_vTria[vNearTria[nTB]].nETempFlag[nVB] = SurfB.m_vTria[vNearTria[nTB]].nTempPart ;
+                     SurfB.m_vTria[vNearTria[nTB]].nETempFlag[nVB] = SurfB.m_vTria[vNearTria[nTB]].nTempShell ;
                }
             }
            // Intersezione edge-edge: salvo indice e vettore triangoli
@@ -1326,171 +1327,28 @@ SurfTriMesh::IntersectTriMeshTriangle( SurfTriMesh& Other)
    bool bRetriangulated = true ;
    if ( ! bModif  &&  ( int( AmbiguosA.size()) == 0  ||  int( AmbiguosB.size()) == 0)) {
       bRetriangulated = false ;
-      int nVertNum = 0 ;
-      Point3d ptFirstV ;
-      int nCurVert = GetFirstVertex( ptFirstV) ;
+     // devo assegnare a tutti i triangoli della superficie la medesima proprietà ( definita da nInOutNum)
+     // ( -1 -> esterno | 0 -> indefinito | +1 -> interno )
+     // devo farlo sia per la SurfA( *this) che per la SurfB
       int nInOutNum = 0 ;
-      while ( nInOutNum == 0  &&  nCurVert != SVT_NULL) {
-         INTVECTOR vnTriaNum ;
-         double dMinDist = DBL_MAX ;
-         for ( int nTB = 0 ; nTB < nTriaNumB ; ++ nTB) {
-           // Se il triangolo B non è valido, continuo
-            Triangle3d trTriaB ;
-            if ( ! SurfB.GetTriangle( nTB, trTriaB)  ||  ! trTriaB.Validate( true))
-               continue ;
-            double dDist ;
-            // potrei trovare più triangolo equidistanti, li salvo tutti
-            if ( DistPointTriangle( ptFirstV, trTriaB).GetDist( dDist)) {
-               if ( abs(dDist - dMinDist) < EPS_SMALL)
-                  vnTriaNum.push_back( nTB) ;
-               else if ( dDist < dMinDist){
-                  vnTriaNum.clear() ;
-                  vnTriaNum.push_back( nTB) ;
-                  dMinDist = dDist ;
-               }
-            }
-         }
-         if ( ! vnTriaNum.empty()) {
-            Triangle3d trTriaB ;
-            bool bSame = true ;
-            // controllo se rispetto a questi triangoli il punto risulta sempre fuori o sempre dentro
-            for ( int nTriaNum : vnTriaNum) {
-               SurfB.GetTriangle( nTriaNum, trTriaB) ;
-               if ( ( ptFirstV - trTriaB.GetP(0)) * trTriaB.GetN() < - EPS_SMALL) {
-                  if ( nInOutNum == 0)
-                     nInOutNum = 1 ;
-                  else if ( nInOutNum == -1) {
-                     bSame = false ;
-                     break ;
-                  }
-               }
-               else if ( ( ptFirstV - trTriaB.GetP(0)) * trTriaB.GetN() > EPS_SMALL) {
-                  if ( nInOutNum == 0)
-                     nInOutNum = -1 ;
-                  else if ( nInOutNum == 1) {
-                     bSame = false ;
-                     break ;
-                  }
-               }
-            }
-            // se le informazioni date dalle normali dei triangoli non sono concordi valuto il triangolo più vicino
-            // e ricalcolo l'informazione che mi dà la sua normale
-            if ( ! bSame) {
-               Point3d ptBar_tot ;
-               for (int nTriaNum : vnTriaNum) {
-                  SurfB.GetTriangle( nTriaNum, trTriaB) ;
-                  ptBar_tot += trTriaB.GetCentroid();
-               }
-               ptBar_tot /= int( vnTriaNum.size()) ;
-               for ( int nTriaNum : vnTriaNum) {
-                  SurfB.GetTriangle( nTriaNum, trTriaB) ;
-                  Point3d ptInters1, ptInters2 ;
-                  int nInters = IntersLineTria( ptFirstV, ptBar_tot, trTriaB, ptInters1, ptInters2, true) ;
-                  if ( nInters == ILTT_NO)
-                     continue ;
-                  else if ( nInters == ILTT_IN) {
-                     if ( ( ptFirstV - trTriaB.GetP( 0)) * trTriaB.GetN() < - EPS_SMALL)
-                        nInOutNum = 1 ;
-                     else if ( ( ptFirstV - trTriaB.GetP( 0)) * trTriaB.GetN() > EPS_SMALL)
-                        nInOutNum = -1 ;
-                     break ;
-                  }
-                  else
-                     nInOutNum = 0 ;
-               }
-            }
-         }
-         if ( nInOutNum == 0) {
-            nCurVert = GetNextVertex( nVertNum, ptFirstV) ;
-            ++ nVertNum ;
-         }
+      for ( int v = 0 ; v < int( m_vVert.size())  &&  nInOutNum == 0 ; ++ v) {
+         double dDist = 0. ;
+         DistPointSurfTm distCalculator( m_vVert[v].ptP, SurfB) ;
+         if ( distCalculator.GetDist( dDist)  &&  dDist > EPS_SMALL)
+            nInOutNum = ( distCalculator.IsPointOnLeftSide() ? 1 : -1) ;
       }
-
       for ( int nTA = 0 ; nTA < nTriaNumA ; ++ nTA)
-         m_vTria[nTA].nTempPart = nInOutNum ;
+         m_vTria[nTA].nTempShell = nInOutNum ;
 
-      nVertNum = 0 ;
-      nCurVert = SurfB.GetFirstVertex( ptFirstV) ;
       nInOutNum = 0 ;
-      while ( nInOutNum == 0  &&  nCurVert != SVT_NULL) {
-         INTVECTOR vnTriaNum ;
-         double dMinDist = DBL_MAX ;
-         for ( int nTA = 0 ; nTA < nTriaNumA ; ++ nTA) {
-            // Se il triangolo A non è valido, continuo
-            Triangle3d trTriaA ;
-            if ( ! ( GetTriangle( nTA, trTriaA)  &&  trTriaA.Validate( true)))
-               continue ;
-            DistPointTriangle DistCalculator( ptFirstV, trTriaA) ;
-            double dDist ;
-            DistCalculator.GetDist( dDist) ;
-            // potrei trovare più triangolo equidistanti, li salvo tutti
-            if ( DistPointTriangle( ptFirstV, trTriaA).GetDist( dDist)) {
-               if ( abs( dDist - dMinDist) < EPS_SMALL)
-                  vnTriaNum.push_back( nTA) ;
-               else if ( dDist < dMinDist){
-                  vnTriaNum.clear() ;
-                  vnTriaNum.push_back( nTA) ;
-                  dMinDist = dDist ;
-               }
-            }
-         }
-         if ( ! vnTriaNum.empty()) {
-            Triangle3d trTriaA ;
-            bool bSame = true ;
-            // controllo se rispetto a questi triangoli il punto risulta sempre fuori o sempre dentro
-            for ( int nTriaNum : vnTriaNum) {
-               GetTriangle( nTriaNum, trTriaA) ;
-               if ( ( ptFirstV - trTriaA.GetP(0)) * trTriaA.GetN() < - EPS_SMALL) {
-                  if ( nInOutNum == 0)
-                     nInOutNum = 1 ;
-                  else if ( nInOutNum == -1) {
-                     bSame = false ;
-                     break ;
-                  }
-               }
-               else if ( ( ptFirstV - trTriaA.GetP(0)) * trTriaA.GetN() > EPS_SMALL) {
-                  if ( nInOutNum == 0)
-                     nInOutNum = -1 ;
-                  else if ( nInOutNum == 1) {
-                     bSame = false ;
-                     break ;
-                  }
-               }
-            }
-            // se le informazioni date dalle normali dei triangoli non sono concordi valuto il triangolo più vicino
-            // e ricalcolo l'informazione che mi dà la sua normale
-            if ( ! bSame) {
-               Point3d ptBar_tot ;
-               for ( int nTriaNum : vnTriaNum) {
-                  GetTriangle( nTriaNum, trTriaA) ;
-                  ptBar_tot += trTriaA.GetCentroid();
-               }
-               ptBar_tot /= int( vnTriaNum.size()) ;
-               for ( int nTriaNum : vnTriaNum) {
-                  GetTriangle( nTriaNum, trTriaA) ;
-                  Point3d ptInters1, ptInters2 ;
-                  int nInters = IntersLineTria( ptFirstV, ptBar_tot, trTriaA, ptInters1, ptInters2, true) ;
-                  if ( nInters == ILTT_NO)
-                     continue ;
-                  else if ( nInters == ILTT_IN) {
-                     if ( ( ptFirstV - trTriaA.GetP( 0)) * trTriaA.GetN() < - EPS_SMALL)
-                        nInOutNum = 1 ;
-                     else if ( ( ptFirstV - trTriaA.GetP( 0)) * trTriaA.GetN() > EPS_SMALL)
-                        nInOutNum = -1 ;
-                     break ;
-                  }
-                  else
-                     nInOutNum = 0 ;
-               }
-            }
-         }
-         if ( nInOutNum == 0) {
-            nCurVert = SurfB.GetNextVertex( nVertNum, ptFirstV) ;
-            ++ nVertNum ;
-         }
+      for ( int v = 0 ; v < int( SurfB.m_vVert.size())  &&  nInOutNum == 0 ; ++ v) {
+         double dDist = 0. ;
+         DistPointSurfTm distCalculator( SurfB.m_vVert[v].ptP, *this) ;
+         if ( distCalculator.GetDist( dDist)  &&  dDist > EPS_SMALL)
+            nInOutNum = ( distCalculator.IsPointOnLeftSide() ? 1 : -1) ;
       }
       for ( int nTB = 0 ; nTB < nTriaNumB ; ++ nTB)
-         SurfB.m_vTria[nTB].nTempPart = nInOutNum ;
+         SurfB.m_vTria[nTB].nTempShell = nInOutNum ;
    }
 
   // Se c'è stata una ritriangolazione di almeno un triangolo, NON siamo nel caso di tutto dentro o tutto fuori.
@@ -1531,12 +1389,12 @@ SurfTriMesh::IntersectTriMeshTriangle( SurfTriMesh& Other)
             Point3d ptTempA, ptTempB ;
             int nIntTypeAB = IntersTriaTria( trTriaA, trTriaB, ptTempA, ptTempB, vTriaAB) ;
             if ( nIntTypeAB == ITTT_OVERLAPS) {
-               m_vTria[nTA].nTempPart = 2 ;
-               SurfB.m_vTria[vNearTria[nTB]].nTempPart = 2 ;
+               m_vTria[nTA].nTempShell = 2 ;
+               SurfB.m_vTria[vNearTria[nTB]].nTempShell = 2 ;
             }
             else if ( nIntTypeAB == ITTT_COUNTER_OVERLAPS) {
-               m_vTria[nTA].nTempPart = -2 ;
-               SurfB.m_vTria[vNearTria[nTB]].nTempPart = -2 ;  
+               m_vTria[nTA].nTempShell = -2 ;
+               SurfB.m_vTria[vNearTria[nTB]].nTempShell = -2 ;  
             }
          }
       }
@@ -1579,70 +1437,139 @@ SurfTriMesh::IdentifyParts( void) const
    return true ;
 }
 
+//----------------------------------------------------------------------------
+bool
+SurfTriMesh::IdentifyShells( void) const
+{
+   for ( int i = 0 ; i < int( m_vTria.size()) ; ++ i) {
+     // salto triangoli cancellati o già assegnati
+      if ( m_vTria[i].nIdVert[0] == SVT_DEL  ||
+           abs( m_vTria[i].nTempShell) != 1)
+         continue ;
+     // set di triangoli da aggiornare
+      set<int> stTria ;
+      stTria.insert( i) ;
+      while ( ! stTria.empty()) {
+        // tolgo un triangolo dal set
+         const auto iIt = stTria.begin() ;
+         int nT = *iIt ;
+         stTria.erase( iIt) ;
+        // aggiorno i triangoli adiacenti
+         for ( int j = 0 ; j < 3 ; ++ j) {
+            if ( m_vTria[nT].nETempFlag[j] != 0)
+               continue ;
+            int nAdjT = m_vTria[nT].nIdAdjac[j] ;
+            if ( nAdjT != SVT_NULL  &&  m_vTria[nAdjT].nTempShell == 0) {
+               m_vTria[nAdjT].nTempShell = m_vTria[nT].nTempShell ;
+               stTria.insert( nAdjT) ;
+            }
+         }
+      }
+   }
+
+   return true ;
+}
+
 //----------------------------------------------------------------------------
 bool
 SurfTriMesh::Add( const ISurfTriMesh& Other)
 {
-  // Le superfici devono essere valide
+
+  // le superfici devono essere valide
    if ( ! IsValid()  ||  ! Other.IsValid())
       return false ;
 
-  // Se la seconda è vuota non devo fare alcunchè
+  // se la seconda è vuota non devo fare alcunchè
    if ( Other.IsEmpty())
       return true ;
 
    m_OGrMgr.Clear() ;
 
-  // Se la prima è vuota, copio la seconda nella prima
+  // se la prima è vuota, copio la seconda nella prima
    if ( IsEmpty()) {
       CopyFrom( &Other) ;
       return true ;
    }
 
+  // clono la superficie B
    SurfTriMesh SurfB ;
    SurfB.CopyFrom( &Other) ;
 
+  // creazione del frame per scalare A e B
    Frame3d frScalingRef ;
    frScalingRef.Set( m_vVert[0].ptP, X_AX, Y_AX, Z_AX) ;
+  // scalo A e B
    Scale( frScalingRef, BOOLEAN_SCALE, BOOLEAN_SCALE, BOOLEAN_SCALE) ;
    SurfB.Scale( frScalingRef, BOOLEAN_SCALE, BOOLEAN_SCALE, BOOLEAN_SCALE) ;
 
-   IntersectTriMeshTriangle( SurfB) ;
-   IdentifyParts() ;
-   SurfB.IdentifyParts() ;
+  // tengo una copia di B ( la superficie B viene modificata durante la ritriangolazione )
+   SurfTriMesh SurfB_cl ;
+   SurfB_cl.CopyFrom( &SurfB) ;
 
+  // tengo una copia di B ( la superficie B viene modificata durante la ritriangolazione )
+   SurfTriMesh SurfA_cl ;
+   SurfA_cl.CopyFrom( this) ;
+ 
+  // ritriangolo le due superfici mediante ogni intersezione Triangolo-Triangolo
+   IntersectTriMeshTriangle( SurfB) ;
+
+  // assegno un medesimo indice ai triangoli che non interferiscono con altri
+   IdentifyShells() ;
+   SurfB.IdentifyShells() ;
+
+  // rimozione dei triangoli di A con proprietà 1 e -2 ( e gestione dei triangoli ambigui, 0)
    int nTriaNumA = GetTriangleSize() ;
    for ( int nTA = 0 ; nTA < nTriaNumA ; ++ nTA) {
-      if ( m_vTria[nTA].nTempPart == 1  ||  m_vTria[nTA].nTempPart == - 2)
+      if ( m_vTria[nTA].nTempShell == 1  ||  m_vTria[nTA].nTempShell == - 2)
          RemoveTriangle( nTA) ;
+     // se triangolo ambiguo
+      if ( m_vTria[nTA].nTempShell == 0) {
+         Triangle3d TriaA ;
+         GetTriangle( nTA, TriaA) ;
+        // rimuovo il triangolo se interno a surfB ( basta controllare un solo vertice(?) )
+         DistPointSurfTm distCalculator( TriaA.GetP( 0), SurfB_cl) ;
+         if ( distCalculator.IsPointOnLeftSide())
+            RemoveTriangle( nTA) ;
+      }
    }
+
+  // aggiunta di tutti i triangoli di B con proprietà -1 ( e gestione dei triangoli ambigui, 0)
    int nPrevMaxTFlag = m_nMaxTFlag ;
    int nTriaNumB = SurfB.GetTriangleSize() ;
    for ( int nTB = 0 ; nTB < nTriaNumB ; ++ nTB) {
-      if ( SurfB.m_vTria[nTB].nTempPart == - 1) {
-         int nNewVert[3] ;
-         for ( int nV = 0 ; nV < 3 ; ++ nV) {
-            nNewVert[nV] = AddVertex( SurfB.m_vVert[SurfB.m_vTria[nTB].nIdVert[nV]].ptP) ;
+      bool bAdd = ( SurfB.m_vTria[nTB].nTempShell == - 1) ;
+      if ( ! bAdd) {
+        // se ambiguo
+         if ( SurfB.m_vTria[nTB].nTempShell == 0) {
+           // aggiungo se il triangolo è a surfA ( basta controllare un solo vertice(?) )
+            DistPointSurfTm distCalculator( SurfB.m_vVert[SurfB.m_vTria[nTB].nIdVert[0]].ptP, SurfA_cl) ;
+            bAdd = ! distCalculator.IsPointOnLeftSide() ;
          }
+      }
+      if ( bAdd) {
+         int nNewVert[3] ;
+         for ( int nV = 0 ; nV < 3 ; ++ nV)
+            nNewVert[nV] = AddVertex( SurfB.m_vVert[SurfB.m_vTria[nTB].nIdVert[nV]].ptP) ;
          if ( nPrevMaxTFlag == m_nMaxTFlag)
             ++ m_nMaxTFlag ;
          AddTriangle( nNewVert, m_nMaxTFlag) ;
       }
    }
 
+  // sistemazioni varie
    bool bOk = ( AdjustVertices()  &&  DoCompacting()) ;
-
    bool bModified = false ;
    bOk = bOk && RemoveDoubleTriangles( bModified) ;
    if ( bModified)
       bOk = bOk && ( AdjustVertices()  &&  DoCompacting()) ;
-
    bOk = bOk && RemoveTJunctions( bModified) ;
    if ( bModified)
       bOk = bOk && ( AdjustVertices()  &&  DoCompacting()) ;
 
+  // scalo alla dimensioni originali
    Scale( frScalingRef, 1. / BOOLEAN_SCALE, 1. / BOOLEAN_SCALE, 1. / BOOLEAN_SCALE) ;
 
+  // semplifico eventuale geometria delle facce
    if ( ! SimplifyFacets())
       LOG_ERROR( GetEGkLogger(), "Error in SimplifyFacets of Stm::Add")
 
@@ -1653,6 +1580,7 @@ SurfTriMesh::Add( const ISurfTriMesh& Other)
 bool
 SurfTriMesh::Intersect( const ISurfTriMesh& Other)
 {
+
   // Le superfici devono essere valide
    if ( ! IsValid()  ||  ! Other.IsValid())
       return false ;
@@ -1681,29 +1609,51 @@ SurfTriMesh::Intersect( const ISurfTriMesh& Other)
    Scale( frScalingRef, BOOLEAN_SCALE, BOOLEAN_SCALE, BOOLEAN_SCALE) ;
    SurfB.Scale( frScalingRef, BOOLEAN_SCALE, BOOLEAN_SCALE, BOOLEAN_SCALE) ;
 
+  // tengo una copia di B ( la superficie B viene modificata durante la ritriangolazione )
+   SurfTriMesh SurfB_cl ;
+   SurfB_cl.CopyFrom( &SurfB) ;
+
+  // tengo una copia di B ( la superficie B viene modificata durante la ritriangolazione )
+   SurfTriMesh SurfA_cl ;
+   SurfA_cl.CopyFrom( this) ;
+
   // ritriangolo le due superfici mediante ogni intersezione Triangolo-Triangolo
    IntersectTriMeshTriangle( SurfB) ;
 
   // assegno un medesimo indice ai triangoli che non interferiscono con altri
-   IdentifyParts() ;
-   SurfB.IdentifyParts() ;
+   IdentifyShells() ;
+   SurfB.IdentifyShells() ;
 
-  // rimozione dei triangoli di A con proprietà -1 e -2
+  // rimozione dei triangoli di A con proprietà -1 e -2 (e gestione dei triangoli ambigui, 0)
    int nTriaNumA = GetTriangleSize() ;
    for ( int nTA = 0 ; nTA < nTriaNumA ; ++ nTA) {
-      if ( m_vTria[nTA].nTempPart == - 1  ||  m_vTria[nTA].nTempPart == - 2)
+      if ( m_vTria[nTA].nTempShell == - 1  ||  m_vTria[nTA].nTempShell == - 2)
          RemoveTriangle( nTA) ;
+     // se triangolo ambiguo
+      else if ( m_vTria[nTA].nTempShell == 0) {
+         Triangle3d TriaA ;
+         GetTriangle( nTA, TriaA) ;
+        // rimuovo il triangolo se fuori a surfB ( basta controllare un solo vertice(?) )
+         DistPointSurfTm distCalculator( TriaA.GetP( 0), SurfB_cl) ;
+         if ( ! distCalculator.IsPointOnLeftSide())
+            RemoveTriangle( nTA) ;
+      }
    }
 
-  // aggiunta dei triangoli di B con proprietà +1
+  // aggiunta dei triangoli di B con proprietà +1 ( e gestione dei triangoli ambigui, 0)
    int nPrevMaxTFlag = m_nMaxTFlag ;
    int nTriaNumB = SurfB.GetTriangleSize() ;
    for ( int nTB = 0 ; nTB < nTriaNumB ; ++ nTB) {
-      if ( SurfB.m_vTria[nTB].nTempPart == 1) {
+      bool bAdd = ( SurfB.m_vTria[nTB].nTempShell == 1) ;
+      if ( ! bAdd  &&  SurfB.m_vTria[nTB].nTempShell == 0) {
+        // aggiungo se il triangolo è interno a surfA ( basta controllare un solo vertice(?) )
+         DistPointSurfTm distCalculator( SurfB.m_vVert[SurfB.m_vTria[nTB].nIdVert[0]].ptP, SurfA_cl) ;
+         bAdd = distCalculator.IsPointOnLeftSide() ;
+      }
+      if ( bAdd) {
          int nNewVert[3] ;
          for ( int nV = 0 ; nV < 3 ; ++ nV)
             nNewVert[nV] = AddVertex( SurfB.m_vVert[SurfB.m_vTria[nTB].nIdVert[nV]].ptP) ;
-
          if ( nPrevMaxTFlag == m_nMaxTFlag)
             ++ m_nMaxTFlag ;
          AddTriangle( nNewVert, m_nMaxTFlag) ;
@@ -1720,33 +1670,6 @@ SurfTriMesh::Intersect( const ISurfTriMesh& Other)
    if ( bModified)
       bOk = bOk && ( AdjustVertices()  &&  DoCompacting()) ;
 
-  // rimuovo tutte le parti esterne all'intersezione
-   int nParts = GetPartCount() ;
-   if ( nParts > 1) {
-      for ( int i = 0 ; i < nParts ; ++ i) {
-        // recupero i triangoli della stessa Part
-         INTVECTOR vTriaPart ;
-         for ( int j = 0 ; j < int( m_vTria.size()) ; ++ j)
-            if ( m_vTria[j].nPart == i)
-               vTriaPart.push_back( j) ;
-        // controllo se il loro box interferisce con il box della superficie B
-         bool bErasePart = true ;
-         for ( int j = 0 ; j < int( vTriaPart.size())  &&  bErasePart ; ++ j) {
-           // Se il triangolo A non è valido, continuo
-            Triangle3d Tria ;
-            if ( ! GetTriangle( j, Tria)  ||  ! Tria.Validate( true))
-               continue ;
-           // Box del triangolo A
-            BBox3d b3dTriaA ; Tria.GetLocalBBox( b3dTriaA) ;
-           // Recupero i triangoli di B che interferiscono col box del triangolo di A
-            INTVECTOR vNearTria ; SurfB.GetAllTriaOverlapBox( b3dTriaA, vNearTria) ;
-            bErasePart = int( vNearTria.size() == 0) ;
-         }
-         if ( bErasePart)
-            RemovePart( i) ;
-      }
-   }
-
   // scalo alle dimensioni originali
    Scale( frScalingRef, 1. / BOOLEAN_SCALE, 1. / BOOLEAN_SCALE, 1. / BOOLEAN_SCALE) ;
 
@@ -1760,42 +1683,77 @@ SurfTriMesh::Intersect( const ISurfTriMesh& Other)
 bool
 SurfTriMesh::Subtract( const ISurfTriMesh& Other)
 {
-  // Le superfici devono essere valide
+
+  // le superfici devono essere valide
    if ( ! IsValid()  ||  ! Other.IsValid())
       return false ;
 
-  // Se una delle due è vuota non devo fare alcunchè
+  // se una delle due è vuota non devo fare alcunchè
    if ( IsEmpty()  ||  Other.IsEmpty())
       return true ;
 
+  // pulisco grafica
    m_OGrMgr.Clear() ;
 
+  // clono superficie B
    SurfTriMesh SurfB ;
    SurfB.CopyFrom( &Other) ;
 
-   Frame3d frScalingRef;
+  // creazione del frame per scalare A e B
+   Frame3d frScalingRef ;
    frScalingRef.Set( m_vVert[0].ptP, X_AX, Y_AX, Z_AX) ;
+  // scalo A e B
    Scale( frScalingRef, BOOLEAN_SCALE, BOOLEAN_SCALE, BOOLEAN_SCALE) ;
    SurfB.Scale( frScalingRef, BOOLEAN_SCALE, BOOLEAN_SCALE, BOOLEAN_SCALE) ;
 
+  // tengo una copia di B ( la superficie B viene modificata durante la ritriangolazione )
+   SurfTriMesh SurfB_cl ;
+   SurfB_cl.CopyFrom( &SurfB) ;
+
+  // tengo una copia di B ( la superficie B viene modificata durante la ritriangolazione )
+   SurfTriMesh SurfA_cl ;
+   SurfA_cl.CopyFrom( this) ;
+ 
+  // ritriangolo le due superfici mediante ogni intersezione Triangolo-Triangolo
    IntersectTriMeshTriangle( SurfB) ;
 
-   IdentifyParts() ;
-   SurfB.IdentifyParts() ;
+  // assegno un medesimo indice ai triangoli che non interferiscono con altri
+   IdentifyShells() ;
+   SurfB.IdentifyShells() ;
 
+  // rimozione dei triangoli di A con proprietà 1 e 2 ( e gestione triangoli ambigui, 0)
    int nTriaNumA = GetTriangleSize() ;
    for ( int nTA = 0 ; nTA < nTriaNumA ; ++ nTA) { 
-      if  ( m_vTria[nTA].nTempPart == 1  ||  m_vTria[nTA].nTempPart == 2)
+      if  ( m_vTria[nTA].nTempShell == 1  ||  m_vTria[nTA].nTempShell == 2)
          RemoveTriangle( nTA) ;
+     // se triangolo ambiguo...
+      if ( m_vTria[nTA].nTempShell == 0) {
+         Triangle3d TriaA ;
+         GetTriangle( nTA, TriaA) ;
+        // rimuovo il triangolo se interno a SurfB ( basta controllare un solo vertice(?) )
+         DistPointSurfTm distCalculator( TriaA.GetP( 0), SurfB_cl) ;
+         if ( distCalculator.IsPointOnLeftSide())
+            RemoveTriangle( nTA) ;
+      }
    }
+
+  // aggiunta dei triangoli di B con proprietà +1 ( e gestione triangoli ambigui, 0)
    int nPrevMaxTFlag = m_nMaxTFlag ;
    int nTriaNumB = SurfB.GetTriangleSize() ;
    for ( int nTB = 0 ; nTB < nTriaNumB ; ++ nTB) {
-      if ( SurfB.m_vTria[nTB].nTempPart == 1) {
-         int nNewVert[3] ;
-         for ( int nV = 0 ; nV < 3 ; ++ nV) {
-            nNewVert[nV] = AddVertex( SurfB.m_vVert[SurfB.m_vTria[nTB].nIdVert[nV]].ptP) ;
+      bool bAdd = SurfB.m_vTria[nTB].nTempShell == 1 ;
+      if ( ! bAdd) {
+        // se ambiguo
+         if ( SurfB.m_vTria[nTB].nTempShell == 0) {
+           // aggiungo il triangolo se interno alla SurfA ( basta controllare un solo vertice(?) )
+            DistPointSurfTm distCalculator( SurfB.m_vVert[SurfB.m_vTria[nTB].nIdVert[0]].ptP, SurfA_cl) ;
+            bAdd = distCalculator.IsPointOnLeftSide() ;
          }
+      }
+      if ( bAdd) {
+         int nNewVert[3] ;
+         for ( int nV = 0 ; nV < 3 ; ++ nV)
+            nNewVert[nV] = AddVertex( SurfB.m_vVert[SurfB.m_vTria[nTB].nIdVert[nV]].ptP) ;
          swap( nNewVert[1], nNewVert[2]) ;
          if ( nPrevMaxTFlag == m_nMaxTFlag)
             ++ m_nMaxTFlag ;
@@ -1803,22 +1761,23 @@ SurfTriMesh::Subtract( const ISurfTriMesh& Other)
       }
    }
 
+  // sistemazioni varie
    bool bOk = ( AdjustVertices()  &&  DoCompacting()) ;
-
    bool bModified = false ;
    bOk = bOk && RemoveDoubleTriangles( bModified) ;
    if ( bModified)
       bOk = bOk && ( AdjustVertices()  &&  DoCompacting()) ;
-
    bOk = bOk && RemoveTJunctions( bModified) ;
    if ( bModified)
       bOk = bOk && ( AdjustVertices()  &&  DoCompacting()) ;
 
+  // scalo alle dimensioni originali
    Scale( frScalingRef, 1. / BOOLEAN_SCALE, 1. / BOOLEAN_SCALE, 1. / BOOLEAN_SCALE) ;
 
+  // semplifico le facce
    if ( ! SimplifyFacets())
       LOG_ERROR( GetEGkLogger(), "Error in SimplifyFacets of Stm::Subtract")
-   
+
    return bOk ;
 }
 
@@ -1842,14 +1801,14 @@ SurfTriMesh::GetSurfClassification( const ISurfTriMesh& ClassifierSurf,
    SurfC.Scale( frScalingRef, BOOLEAN_SCALE, BOOLEAN_SCALE, BOOLEAN_SCALE) ;
 
    IntersectTriMeshTriangle( SurfC) ;
-   IdentifyParts() ;
+   IdentifyShells() ;
    Scale( frScalingRef, 1. / BOOLEAN_SCALE, 1. / BOOLEAN_SCALE, 1. / BOOLEAN_SCALE) ;
   
    int nTriaNum = GetTriangleSize() ;
    for ( int nT = 0 ; nT < nTriaNum ; ++ nT) {
       if ( m_vTria[nT].nIdVert[0] == SVT_DEL)
          continue ;
-      switch ( m_vTria[nT].nTempPart) {
+      switch ( m_vTria[nT].nTempShell) {
       case -2 : 
          vTriaOnM.push_back( nT) ;
          break ;
@@ -1894,13 +1853,13 @@ SurfTriMesh::CutWithOtherSurf( const ISurfTriMesh& CutterSurf, bool bInVsOut, bo
    SurfC.Scale( frScalingRef, BOOLEAN_SCALE, BOOLEAN_SCALE, BOOLEAN_SCALE) ;
 
    IntersectTriMeshTriangle( SurfC) ;
-   IdentifyParts() ;
+   IdentifyShells() ;
 
    int nPartToRemove = ( bInVsOut ? -1 : 1) ;
    int nCoplanarPartToRemove = ( bSaveOnEq ? ( nPartToRemove ? -2 : 2) : 5) ;
    int nTriaNum = GetTriangleSize() ;
    for ( int nT = 0 ; nT < nTriaNum ; ++ nT) {
-      if ( m_vTria[nT].nTempPart == nPartToRemove  ||  m_vTria[nT].nTempPart == nCoplanarPartToRemove)
+      if ( m_vTria[nT].nTempShell == nPartToRemove  ||  m_vTria[nT].nTempShell == nCoplanarPartToRemove)
          RemoveTriangle( nT) ;
    }
 

SurfTriMeshCuts.cpp

@@ -15,10 +15,10 @@
 #include "CurveComposite.h"
 #include "SurfTriMesh.h"
 #include "SurfFlatRegion.h"
-#include "DistPointLine.h"
 #include "Triangulate.h"
 #include "GeoConst.h"
 #include "/EgtDev/Include/EGkSfrCreate.h"
+#include "/EgtDev/Include/EGkDistPointLine.h"
 #include "/EgtDev/Include/EGkDistPointCurve.h"
 #include "/EgtDev/Include/EGkIntersLinePlane.h"
 #include "/EgtDev/Include/EGkIntersPlaneTria.h"
@@ -79,7 +79,7 @@ SurfTriMesh::CutByTriangles( const Plane3d& plPlane, bool bSaveOnEq, bool& bModi
       if ( m_vVert[i].nIdTria == SVT_DEL)
          continue ;
       double dDist = DistPointPlane( m_vVert[i].ptP, plPlane) ;
-      if ( abs( dDist) < EPS_SMALL)
+      if ( abs( dDist) < 1.1 * EPS_SMALL)
          m_vVert[i].nTemp = 0 ;
       else if ( dDist > 0)
          m_vVert[i].nTemp = +1 ;
@@ -781,14 +781,14 @@ SurfTriMesh::GeneralizedCut( const ICurve& cvCurve, bool bSaveOnEq)
       }
    }
 
-  // Ripristino scala originale
-   Scale( frScalingRef, 1. / CUT_SCALE, 1. / CUT_SCALE, 1. / CUT_SCALE) ;
-
   // se superficie originale a facce, cerco di semplificarle in ogni caso
    if ( nFacetOriCnt < 200  ||  double( nTriaOriCnt) / nFacetOriCnt > 4) {
-      if ( ! SimplifyFacets())
+      if ( ! SimplifyFacets( CUT_SCALE * MAX_EDGE_LEN_STD))
          LOG_ERROR( GetEGkLogger(), "Error in SimplifyFacets of Stm::GeneralizedCut")
    }
 
+  // Ripristino scala originale
+   Scale( frScalingRef, 1. / CUT_SCALE, 1. / CUT_SCALE, 1. / CUT_SCALE) ;
+
    return true ;
 }

SurfTriMeshUtilities.cpp

@@ -15,8 +15,8 @@
 #include "SurfTriMesh.h"
 #include "CurveLine.h"
 #include "Triangulate.h"
-#include "DistPointLine.h"
-#include "DistLineLine.h"
+#include "/EgtDev/Include/EGkDistPointLine.h"
+#include "/EgtDev/Include/EGkDistLineLine.h"
 #include <unordered_map>
 
 using namespace std ;

Tool.cpp

@@ -583,34 +583,78 @@ Tool::SetChiselTool( const string& sToolName, double dH, double dW, double dTh,
 
 //----------------------------------------------------------------------------
 bool
-Tool::SetAdditiveTool( const std::string& sToolName, double dH, double dR, double dRC, int nToolNum)
+Tool::SetAdditiveTool( const std::string& sToolName, double dH, double dR, double dRc, int nToolNum)
 {
-   // Impostazioni generali
+  // Impostazioni generali
    m_sName = sToolName ;
    m_nCurrentNum = nToolNum ;
    m_nType = UNDEF ;
    m_Outline.Clear() ;
    m_ArcLineApprox.Clear() ;
 
-   // verifica sulle minime dimensioni globali
-   if ( dH < EPS_SMALL  ||  dR < EPS_SMALL  ||  dRC < - EPS_SMALL)
+  // Verifica sulle dimensioni
+   if ( dH < 10 * EPS_SMALL  ||  dR < 10 * EPS_SMALL  ||  dRc < - EPS_SMALL  ||  dR < dRc - EPS_SMALL)
       return false ;
 
+  // Assegnazione dati geometrici principali
    m_dHeight = dH ;
    m_dRadius = dR ;
-   m_dRCorner = dRC ;
+   m_dRCorner = min( dRc, dR) ;
+   if ( m_dRCorner < 10 * EPS_SMALL)
+      m_dRCorner = 0 ;
    m_dTipHeight = 0 ;
    m_dTipRadius = 0 ;
    m_dRefRadius = 0 ;
    m_dCutterHeight = dH ;
 
-   bool bToolDefined = true ;
-
-   double dSquareCornerRadProj = m_dRCorner * m_dRCorner - 0.25 * m_dHeight * m_dHeight ;
-   // Utensile sfiancato
-   if ( dSquareCornerRadProj > 0) {
+  // Utensile cilindrico con eventuale raggio corner (al limite sferico)
+   if ( m_dRCorner <= m_dHeight / 2) {
+     // Se Cilindrico
+      if ( m_dRCorner < EPS_SMALL) {
+        // inizio
+         m_Outline.AddPoint( Point3d( 0, 0, 0)) ;
+        // segmento orizzontale in alto
+         m_Outline.AddLine( Point3d( m_dRadius, 0, 0)) ;
+         m_Outline.SetCurveTempProp( 0, 1, 1) ;
+        // segmento verticale
+         m_Outline.AddLine( Point3d( m_dRadius, -m_dHeight, 0)) ;
+         m_Outline.SetCurveTempProp( 1, 1, 1) ;
+        // segmento orizzontale in basso
+         m_Outline.AddLine( Point3d( 0, -m_dHeight, 0)) ;
+         m_Outline.SetCurveTempProp( 2, 1, 1) ;
+      }
+     // Generico
+      else {
+        // inizio
+         m_Outline.AddPoint( Point3d( 0, 0, 0)) ;
+         int nInd = -1 ;
+        // eventuale segmento orizzontale in alto
+         if ( m_Outline.AddLine( Point3d( m_dRadius - m_dRCorner, 0, 0)))
+            m_Outline.SetCurveTempProp( ++ nInd, 1, 1) ;
+        // raggio corner in alto
+         CurveArc cvArc ;
+         cvArc.SetC2P( Point3d( m_dRadius - m_dRCorner, -m_dRCorner, 0), Point3d( m_dRadius - m_dRCorner, 0, 0), Point3d( m_dRadius, -m_dRCorner, 0)) ;
+         m_Outline.AddCurve( cvArc, true, 2 * EPS_SMALL) ;
+         m_Outline.SetCurveTempProp( ++ nInd, 1, 1) ;
+        // eventuale segmento verticale
+         if ( m_Outline.AddLine( Point3d( m_dRadius, -m_dHeight + m_dRCorner, 0)))
+            m_Outline.SetCurveTempProp( ++ nInd, 1, 1) ;
+        // raggio corner in basso
+         cvArc.SetC2P( Point3d( m_dRadius - m_dRCorner, -m_dHeight + m_dRCorner, 0), Point3d( m_dRadius, -m_dHeight + m_dRCorner, 0), Point3d( m_dRadius - m_dRCorner, -m_dHeight, 0)) ;
+         m_Outline.AddCurve( cvArc, true, 2 * EPS_SMALL) ;
+         m_Outline.SetCurveTempProp( ++ nInd, 1, 1) ;
+        // eventuale segmento orizzontale in basso
+         if ( m_Outline.AddLine( Point3d( 0, -m_dHeight, 0)))
+            m_Outline.SetCurveTempProp( ++ nInd, 1, 1) ;
+      }
+      m_Outline.SetTempProp( 1, 1) ;
+      if ( ! SetGenTool( sToolName, &m_Outline, nToolNum))
+         return false ;
+   } 
+  // Utensile sfiancato
+   else {
       double dCenX = m_dRadius - m_dRCorner ;
-      double dCylRad = dCenX + sqrt( dSquareCornerRadProj) ;
+      double dCylRad = dCenX + sqrt( m_dRCorner * m_dRCorner - m_dHeight * m_dHeight / 4) ;
       // Utensile mal definito
       if ( dCylRad < EPS_SMALL)
          return false ;
@@ -625,41 +669,10 @@ Tool::SetAdditiveTool( const std::string& sToolName, double dH, double dR, doubl
       m_Outline.AddLine( Point3d( 0, - m_dHeight, 0)) ;
       m_Outline.SetCurveTempProp( 2, 1, 1) ;
       m_Outline.SetTempProp( 1, 1) ;
-      bToolDefined = SetGenTool( sToolName, &m_Outline, nToolNum) ;
-   }
-   // Utensile cilindrico con eventuale raggio corner 
-   else {
-      // Utensile mal definito
-      if ( m_dRadius - m_dRCorner < 0)
+      if ( ! SetGenTool( sToolName, &m_Outline, nToolNum))
          return false ;
-      // Utensile sferico
-      else if ( m_dRadius - m_dRCorner < EPS_SMALL)
-         ;
-      // Raggio corner nullo: cilindro
-      else if ( m_dRCorner < EPS_SMALL) {
-         ;
-      }
-      else {
-         // Profilo
-         m_Outline.AddPoint( Point3d( 0, 0, 0)) ;
-         m_Outline.AddLine( Point3d( m_dRadius - m_dRCorner, 0, 0)) ;
-         m_Outline.SetCurveTempProp( 0, 1, 1) ;
-         CurveArc cvArc ;
-         cvArc.SetC2P( Point3d( m_dRadius - m_dRCorner, - m_dRCorner, 0), Point3d( m_dRadius - m_dRCorner, 0, 0), Point3d( m_dRadius, - m_dRCorner, 0)) ;
-         m_Outline.AddCurve( cvArc) ;
-         m_Outline.SetCurveTempProp( 1, 1, 1) ;
-         m_Outline.AddLine( Point3d( m_dRadius, - m_dHeight + m_dRCorner, 0)) ;
-         m_Outline.SetCurveTempProp( 2, 1, 1) ;
-         cvArc.SetC2P( Point3d( m_dRadius - m_dRCorner, - m_dHeight + m_dRCorner, 0), Point3d( m_dRadius, - m_dHeight + m_dRCorner, 0), Point3d( m_dRadius - m_dRCorner, - m_dHeight, 0)) ;
-         m_Outline.AddCurve( cvArc) ;
-         m_Outline.SetCurveTempProp( 3, 1, 1) ;
-         m_Outline.AddLine( Point3d( 0, - m_dHeight, 0)) ;
-         m_Outline.SetCurveTempProp( 4, 1, 1) ;
-         m_Outline.SetTempProp( 1, 1) ;
-         bToolDefined = SetGenTool( sToolName, &m_Outline, nToolNum) ;
-      }
-   } 
+   }
 
    m_nType = ADDITIVE ;
-   return bToolDefined ;
-}
+   return true ;
+}

Tool.h

@@ -33,7 +33,7 @@ class Tool
       bool SetGenTool( const std::string& sToolName, const ICurveComposite* pToolOutline, int nToolNum) ;
       bool SetMortiserTool( const std::string& sToolName, double dH, double dW, double dTh, double dRc, int nToolNum) ;
       bool SetChiselTool( const std::string& sToolName, double dH, double dW, double dTh, int nToolNum) ;
-      bool SetAdditiveTool( const std::string& sToolName, double dH, double dR, double dRC, int nToolNum) ;
+      bool SetAdditiveTool( const std::string& sToolName, double dH, double dR, double dRc, int nToolNum) ;
       bool SetToolNum( int nToolNum)
                { m_nCurrentNum = nToolNum ; return true ; }
       int GetType() const

VolZmap.cpp

@@ -18,13 +18,14 @@
 #include "GdbGeo.h"
 #include "NgeWriter.h"
 #include "NgeReader.h"
-#include "DistPointLine.h"
 #include "GeoConst.h"
+#include "/EgtDev/Include/EGkDistPointLine.h"
 #include "/EgtDev/Include/EGkIntersLinePlane.h"
 #include "/EgtDev/Include/EGkUiUnits.h"
 #include "/EgtDev/Include/EGkIntervals.h"
 #include "/EgtDev/Include/EGkStringUtils3d.h"
 #include "/EgtDev/Include/EgtNumUtils.h"
+#include "/EgtDev/Include/EGkIntersLineSurfTm.h"
 #include <thread>
 #include <future>
 
@@ -49,7 +50,6 @@ VolZmap::VolZmap(void)
       m_nFracLin[i] = 0 ;
    }
    m_vTool.resize( 1) ;
-   m_nCurrTool = 0 ;
 }
 
 //----------------------------------------------------------------------------
@@ -1501,7 +1501,187 @@ VolZmap::GetPartMinDistFromPoint( const Point3d& ptP) const
 
 //----------------------------------------------------------------------------
 bool
-VolZmap::SetToModifyDexelBlocks( int nGrid, int nDex, int nInt) 
+VolZmap::AddSurfTm( const ISurfTriMesh* pStm)
+{
+  // controllo sulla superficie
+   double dVol ;
+   if ( pStm == nullptr  ||  ! pStm->IsValid()  ||  ! pStm->IsClosed()  ||
+        ! pStm->GetVolume( dVol)  ||  dVol < 0)
+      return false ;
+
+  // controllo se il Box3d della superficie si interseca con il Box3d dello Zmap corrente
+   BBox3d BBox_stm, BBox_curr ;
+   if ( ! pStm->GetLocalBBox( BBox_stm)  ||  ! GetLocalBBox( BBox_curr))
+      return false ;
+   BBox3d BBox_inters ;
+   if ( BBox_stm.FindIntersection( BBox_curr, BBox_inters)  &&  BBox_inters.IsEmpty())
+      return true ; // se non ci sono intersezioni, la superficie non influenza lo Zmap
+   Vector3d vtLen = BBox_curr.GetMax() - BBox_curr.GetMin() ; // dimensione massima dello spillone
+
+  // ciclo sulle griglie
+   bool bCompleted = true ;
+   for ( int g = 0 ; g < m_nMapNum ; ++ g) {
+     // definisco dei sistemi di riferimento ausiliari
+      Frame3d frMapFrame ;
+      if ( g == 0)
+         frMapFrame = m_MapFrame ;
+      else if ( g == 1)
+         frMapFrame.Set( m_MapFrame.Orig(), Y_AX, Z_AX, X_AX) ;
+      else if ( g == 2)
+         frMapFrame.Set( m_MapFrame.Orig(), Z_AX, X_AX, Y_AX) ;
+
+     // oggetto per calcolo massivo intersezioni
+      IntersParLinesSurfTm intPLSTM( frMapFrame, *pStm) ;
+
+     // numero massimo di thread
+      int nThreadMax = max( 1, int( thread::hardware_concurrency()) - 1) ;
+      vector<future<bool>> vRes ;
+      vRes.resize( nThreadMax) ;
+     // se dimensione griglia in X maggiore di dimensione Y
+      if ( m_nNx[g] > m_nNy[g]) {
+         int nDexNum = m_nNx[g] / nThreadMax ;
+         int nRemainder = m_nNx[g] % nThreadMax ;
+         int nInfI = 0 ;
+         int nSupI = 0 ;
+        // aggiungo le parti interessate alla mappa
+         for ( int nThread = 0 ; nThread < nThreadMax ; ++ nThread) {
+            nInfI = nSupI ;
+            nSupI = nInfI + ( nThread < nRemainder ? nDexNum + 1 : nDexNum) ;
+            vRes[nThread] =  async( launch::async, &VolZmap::AddMapPart, this, g,
+                                    nInfI, nSupI, 0, m_nNy[g], ref( vtLen), ref( m_MapFrame.Orig()),
+                                    ref( *pStm), ref( intPLSTM)) ;
+         }
+      }
+     // se dimensione griglia in Y maggiore di dimensione X
+      else {
+         int nDexNum = m_nNy[g] / nThreadMax ;
+         int nRemainder = m_nNy[g] % nThreadMax ;
+         int nInfJ = 0 ;
+         int nSupJ = 0 ;
+        // aggiungo le parti interessate alla mappa
+         for ( int nThread = 0 ; nThread < nThreadMax ; ++ nThread) {
+            nInfJ = nSupJ ;
+            nSupJ = nInfJ + ( nThread < nRemainder ? nDexNum + 1 : nDexNum) ;
+            vRes[nThread] =  async( launch::async, &VolZmap::AddMapPart, this, g,
+                                    0, m_nNx[g], nInfJ, nSupJ, ref( vtLen), ref( m_MapFrame.Orig()),
+                                    ref( *pStm), ref( intPLSTM)) ;
+         }
+      }
+
+     // ciclo per attendere che tutti gli async abbiano terminato.
+      int nTerminated = 0 ;
+      while ( nTerminated < nThreadMax) {
+         for ( int nL = 0 ; nL < nThreadMax ; ++ nL) {
+           // async terminato
+            if ( vRes[nL].valid()  &&  vRes[nL].wait_for( chrono::microseconds{ 1}) == future_status::ready) {
+               ++ nTerminated ;
+               bCompleted = bCompleted  &&  vRes[nL].get() ;
+            }
+         }
+      }
+
+      if ( ! bCompleted)
+         return false ;
+   }
+
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+bool
+VolZmap::MakeUniform( double dToler)
+{
+  // controllo validità dello Zmap
+   if ( ! IsValid())
+      return false ;
+  // la tolleranza deve essere minore dello step
+   dToler = min( dToler, 0.95 * m_dStep) ;
+
+  // creo lo Zmpa che andrà a sostituire il corrente
+   PtrOwner<VolZmap> pOldVolZmap( CloneBasicVolZmap( this)) ;
+   if ( IsNull( pOldVolZmap))
+      return false ;
+
+  // ciclo sulle griglie
+   for ( int nGrid = 0 ; nGrid < m_nMapNum ; ++ nGrid) {
+     // salvo lo Zmap prima di modificare gli spilloni
+      PtrOwner<VolZmap> pVolZMapCurrGrid( CloneBasicVolZmap( this)) ;
+      if ( IsNull( pVolZMapCurrGrid))
+      return false ;
+     // ciclo sul numero di dexel presenti
+      for ( int nDex = 0 ; nDex < int( m_Values[nGrid].size()) ; ++ nDex) {
+        // se il dexel corrente non ha sotto-intervalli passo al successivo
+         if ( int( m_Values[nGrid][nDex].size()) == 0)
+            continue ;
+        // indici del dexel
+         int nI = nDex % m_nNx[nGrid] ;
+         int nJ = nDex / m_nNx[nGrid] ;
+        // salvo le informazioni dei sotto-intervalli del dexel corrente
+         vector<Data> vInfo ;
+         for ( int nExtr = 0 ; nExtr < int( m_Values[nGrid][nDex].size()) ; ++ nExtr)
+            vInfo.push_back( m_Values[nGrid][nDex][nExtr]) ;
+        // per ogni sotto-intervallo, estendo a destra e a sinistra della tolleranza
+         int nSub_intervals = int( vInfo.size()) ;
+        // scorro gli intervalli
+         for ( int nInfo = 0 ; nInfo < int( pVolZMapCurrGrid->m_Values[nGrid][nDex].size()) ; ++ nInfo) {
+           // estremo inferiore
+            if ( vInfo[nInfo].dMin - dToler > m_dMinZ[nGrid]) {
+               AddIntervals( nGrid, nI, nJ,
+                             vInfo[nInfo].dMin - dToler,
+                             vInfo[nInfo].dMin + dToler,
+                             vInfo[nInfo].vtMinN, vInfo[nInfo].vtMinN, vInfo[nInfo].nToolMin,
+                             true) ;
+              // se si sono uniti degli intervalli, potrei dover aggiungere degli spilloni nelle altre due
+              // direzioni nel voxel corrispondente
+               if ( IsTriDexel()  &&  dToler > 0.5 * m_dStep - EPS_SMALL  &&
+                    nSub_intervals != int( m_Values[nGrid][nDex].size())) {
+                  nSub_intervals = int( m_Values[nGrid][nDex].size()) ; // aggiorno gli intervalli correnti
+                 // l'intervallo corrente si è unito con il precedente...
+                  AddMissingIntervalsInVoxel( pOldVolZmap, nGrid, nI, nJ, vInfo[nInfo].dMin, dToler,
+                                              vInfo[nInfo].vtMinN, vInfo[nInfo].vtMinN,
+                                              vInfo[nInfo].nToolMin) ;
+               }
+            }
+           // estremo superiore
+            if ( vInfo[nInfo].dMax + dToler < m_dMaxZ[nGrid]) {
+               AddIntervals( nGrid, nI, nJ,
+                             vInfo[nInfo].dMax - dToler,
+                             vInfo[nInfo].dMax + dToler,
+                             vInfo[nInfo].vtMaxN, vInfo[nInfo].vtMaxN, vInfo[nInfo].nToolMax,
+                             true) ;
+               if ( IsTriDexel()  &&  dToler > 0.5 * m_dStep - EPS_SMALL  &&
+                    nSub_intervals != int( m_Values[nGrid][nDex].size())) {
+                  nSub_intervals = int( m_Values[nGrid][nDex].size()) ; // aggiorno gli intervalli correnti
+                  AddMissingIntervalsInVoxel( pOldVolZmap, nGrid, nI, nJ, vInfo[nInfo].dMax, dToler,
+                                              vInfo[nInfo].vtMaxN, vInfo[nInfo].vtMaxN,
+                                              vInfo[nInfo].nToolMax) ;
+               }
+            }
+         }
+        // per ogni sotto-intervallo ricavato fino ad ora, restringo della tolleranza
+        // ( NB. avendo aggiunto intervalli, il dexel può modificare la sua struttura interna )
+         for ( int nInfo = 0 ; nInfo < int( m_Values[nGrid][nDex].size()) ; ++ nInfo) {
+           // ( NB. la rimozione di un intervallo ora va definita per intervalli a destra e a sinistra,
+           // altrimenti rimuovo parti in eccesso )
+            if ( ! pVolZMapCurrGrid->m_Values[nGrid][nDex].empty()) {
+               if ( nInfo != 0  ||
+                    pVolZMapCurrGrid->m_Values[nGrid][nDex][0].dMin - dToler > m_dMinZ[nGrid])
+                  m_Values[nGrid][nDex][nInfo].dMin += dToler ;
+               if ( nInfo != int( m_Values[nGrid][nDex].size()) - 1  ||
+                    pVolZMapCurrGrid->m_Values[nGrid][nDex].back().dMax + dToler < m_dMaxZ[nGrid])
+                  m_Values[nGrid][nDex][nInfo].dMax -= dToler ;
+            }
+         }
+      }
+   }
+
+   return true ;
+}
+
+
+//----------------------------------------------------------------------------
+bool
+VolZmap::SetToModifyDexelBlocks( int nGrid, int nDex, int nInt)
 {
   // Controllo sulla validità della griglia
    if ( nGrid < 0  ||  nGrid > 2)
@@ -2175,7 +2355,7 @@ VolZmap::SetChiselTool( const string& sToolName, double dH, double dW, double dT
 //----------------------------------------------------------------------------
 bool
 VolZmap::SetAdditiveTool( const std::string& sToolName,
-                          double dH, double dR, double dRC, int nFlag, bool bFirst)
+                          double dH, double dR, double dRc, int nFlag, bool bFirst)
 {
    if ( bFirst) {
       m_vTool.resize( 1) ;
@@ -2187,10 +2367,9 @@ VolZmap::SetAdditiveTool( const std::string& sToolName,
    if ( m_nCurrTool < 0)
       return false ;
    m_vTool[m_nCurrTool].SetTolerances( m_dToolLinTol, m_dToolAngTolDeg) ;
-   return m_vTool[m_nCurrTool].SetAdditiveTool( sToolName, dH, dR, dRC, nFlag) ;
+   return m_vTool[m_nCurrTool].SetAdditiveTool( sToolName, dH, dR, dRc, nFlag) ;
 }
 
-
 //----------------------------------------------------------------------------
 int
 VolZmap::GetToolCount( void) const

VolZmap.h

@@ -1,7 +1,7 @@
 //----------------------------------------------------------------------------
-//  EgalTech 2015-2023
+//  EgalTech 2015-2024
 //----------------------------------------------------------------------------
-// File : VolZmap.h           Data : 12.09.23          Versione : 2.5i1
+// File : VolZmap.h           Data : 22.04.24          Versione : 2.6d4
 // Contenuto : Dichiarazione della classe Volume Zmap.
 //
 //
@@ -78,7 +78,7 @@ class VolZmap : public IVolZmap, public IGeoObjRW
       bool CopyFrom( const IGeoObj* pGObjSrc) override ;
       bool Clear( void) override ;
       bool Create( const Point3d& ptO, double dDimX, double dDimY, double dDimZ, double dStep, bool bTriDex) override ;
-      bool CreateEmptyMap( const Point3d& ptO, double dLengthX, double dLengthY, double dLengthZ, double dStep, bool bTriDex) override ;
+      bool CreateEmpty( const Point3d& ptO, double dDimX, double dDimY, double dDimZ, double dStep, bool bTriDex) override ;
       bool CreateFromFlatRegion( const ISurfFlatRegion& Surf, double dDimZ, double dStep, bool bTriDex) override ;
       bool CreateFromTriMesh( const ISurfTriMesh& Surf, double dStep, bool bTriDex) override ;
       int  GetBlockCount( void) const override ;
@@ -111,7 +111,7 @@ class VolZmap : public IVolZmap, public IGeoObjRW
       bool SetChiselTool( const std::string& sToolName,
                           double dH, double dW, double dTh, int nFlag, bool bFirst) override ;
       bool SetAdditiveTool( const std::string& sToolName,
-                            double dH, double dR, double dRC, int nFlag, bool bFirst) override ;
+                            double dH, double dR, double dRc, int nFlag, bool bFirst) override ;
       int  GetToolCount( void) const override ;
       bool SetCurrTool( int nCurrTool) override ;
       bool ResetTools( void) override ;
@@ -143,6 +143,8 @@ class VolZmap : public IVolZmap, public IGeoObjRW
       VolZmap* ClonePart( int nPart) const override ;
       bool RemovePart( int nPart) override ;
       int  GetPartMinDistFromPoint( const Point3d& ptP) const override ;
+      bool AddSurfTm( const ISurfTriMesh* pStm) override ;
+      bool MakeUniform( double dToler) override ;
 
    public :                      // IGeoObjRW
       int GetNgeId( void) const override ;
@@ -243,9 +245,15 @@ class VolZmap : public IVolZmap, public IGeoObjRW
                         INTVECTOR& vAdjBlockVoxComp, INTVECTOR& vAdjBordBlockVoxComp) const ;
      // OPERAZIONI SU INTERVALLI
       bool SubtractIntervals( int nGrid, int nI, int nJ,
-                              double dMin, double dMax, const Vector3d& vtNMin, const Vector3d& vtNMax, int nToolNum) ;
+                              double dMin, double dMax, const Vector3d& vtNMin, const Vector3d& vtNMax, 
+                              int nToolNum, bool bSkipSwap = false) ;
       bool AddIntervals( int nGrid, int nI, int nJ,
-                         double dMin, double dMax, const Vector3d& vtNMin, const Vector3d& vtNMax, int nToolNum) ;
+                         double dMin, double dMax, const Vector3d& vtNMin, const Vector3d& vtNMax,
+                         int nToolNum, bool bSkipSwap = false) ;
+      bool AddMissingIntervalsInVoxel( VolZmap* VolZmapRef, int nGrid, int nI, int nJ, double dZ, double dToler,
+                                       Vector3d vtToolMin, Vector3d vtToolMax, int nToolNum) ;
+      bool AddSubIntervalInVoxel( VolZmap* VolZmapRef, int nGrid, int nI, int nJ, int nK, double& dMin, double& dMax,
+                                  Vector3d& vtMin, Vector3d& vtMax) ;
      // Spostamenti utensile
       bool MillingTranslationStep( const Point3d& ptPs, const Point3d& ptPe, const Vector3d& vtD, const Vector3d& vtA) ;
       bool MillingGeneralMotionStep( const Point3d& ptPs, const Vector3d& vtDs, const Vector3d& vtAs,
@@ -422,6 +430,8 @@ class VolZmap : public IVolZmap, public IGeoObjRW
       // Funzione per crezione solido in parallelo
       bool CreateMapPart( int nMap, int nInfI, int nSupI, int nInfJ, int nSupJ, const Vector3d& vtLen, const Point3d& ptMapOrig,
                           const ISurfTriMesh& Surf, IntersParLinesSurfTm& intPLSTM) ;
+      bool AddMapPart( int nMap, int nInfI, int nSupI, int nInfJ, int nSupJ, const Vector3d& vtLen, const Point3d& ptMapOrig,
+                       const ISurfTriMesh& Surf, IntersParLinesSurfTm& intPLSTM) ;
 
    private :
       enum Status { ERR = 0, OK = 1, TO_VERIFY = 2} ;
@@ -455,6 +465,7 @@ class VolZmap : public IVolZmap, public IGeoObjRW
       } ;
       std::vector<std::vector<Data>> m_Values[N_MAPS] ;  // dexel delle 3 griglie
 
+
       int m_nShape ;               // Forma : 0 generica, 1 box, 2 estrusione
 
       int m_nVoxNumPerBlock ;      // Numero di voxel per blocco

VolZmapCreation.cpp

@@ -1,4 +1,4 @@
-//----------------------------------------------------------------------------
+//----------------------------------------------------------------------------
 //  EgalTech 2015-2016
 //----------------------------------------------------------------------------
 // File : VolZmap.cpp         Data : 22.01.15           Versione : 1.6a4
@@ -27,10 +27,10 @@ using namespace std ;
 
 //----------------------------------------------------------------------------
 bool 
-VolZmap::Create( const Point3d& ptO, double dLengthX, double dLengthY, double dLengthZ, double dStep, bool bTriDex)
+VolZmap::Create( const Point3d& ptO, double dDimX, double dDimY, double dDimZ, double dStep, bool bTriDex)
 {
   // Controlli sull'ammissibilità delle dimensioni lineari del grezzo e del passo
-   if ( dStep < EPS_SMALL  ||  dLengthX < EPS_SMALL  ||  dLengthY < EPS_SMALL  ||  dLengthZ < EPS_SMALL)
+   if ( dStep < EPS_SMALL  ||  dDimX < EPS_SMALL  ||  dDimY < EPS_SMALL  ||  dDimZ < EPS_SMALL)
       return false ;
 
   // Il passo di discretizzazione non può essere inferiore a 100 * EPS_SMALL
@@ -47,8 +47,8 @@ VolZmap::Create( const Point3d& ptO, double dLengthX, double dLengthY, double dL
    m_MapFrame.Set( ptO, X_AX, Y_AX, Z_AX) ;
 
   // Definisco i vettori dei limiti su indici
-   m_nNx[0] = max( int( ( dLengthX + EPS_SMALL) / m_dStep + 0.5), 1) ;
-   m_nNy[0] = max( int( ( dLengthY + EPS_SMALL) / m_dStep + 0.5), 1) ;
+   m_nNx[0] = max( int( ( dDimX + EPS_SMALL) / m_dStep + 0.5), 1) ;
+   m_nNy[0] = max( int( ( dDimY + EPS_SMALL) / m_dStep + 0.5), 1) ;
 
   // Numero di componenti connesse 
    m_nConnectedCompoCount = 1 ;
@@ -56,7 +56,7 @@ VolZmap::Create( const Point3d& ptO, double dLengthX, double dLengthY, double dL
   // Se tridexel
    if ( bTriDex) {
       m_nNx[1] = m_nNy[0] ;
-      m_nNy[1] = max( int( ( dLengthZ + EPS_SMALL) / m_dStep + 0.5), 1) ;
+      m_nNy[1] = max( int( ( dDimZ + EPS_SMALL) / m_dStep + 0.5), 1) ;
       m_nNx[2] = m_nNy[1] ;
       m_nNy[2] = m_nNx[0] ;
    }
@@ -97,19 +97,19 @@ VolZmap::Create( const Point3d& ptO, double dLengthX, double dLengthY, double dL
          switch ( i) {
          case 0 :
             m_Values[i][j][0].vtMinN = - Z_AX ;
-            m_Values[i][j][0].dMax = dLengthZ ; 
+            m_Values[i][j][0].dMax = dDimZ ; 
             m_Values[i][j][0].vtMaxN = Z_AX ;
             m_Values[i][j][0].nToolMax = 0 ;
             break ;
          case 1 :
             m_Values[i][j][0].vtMinN = - X_AX ;
-            m_Values[i][j][0].dMax = dLengthX ;
+            m_Values[i][j][0].dMax = dDimX ;
             m_Values[i][j][0].vtMaxN = X_AX ;
             m_Values[i][j][0].nToolMax = 0 ;
             break ;
          case 2 :
             m_Values[i][j][0].vtMinN = - Y_AX ;
-            m_Values[i][j][0].dMax = dLengthY ;
+            m_Values[i][j][0].dMax = dDimY ;
             m_Values[i][j][0].vtMaxN = Y_AX ;
             m_Values[i][j][0].nToolMax = 0 ;
             break ; 
@@ -118,11 +118,11 @@ VolZmap::Create( const Point3d& ptO, double dLengthX, double dLengthY, double dL
 
   // Definizione delle limitazioni iniziali in Z per ogni mappa
    m_dMinZ[0] = 0 ;
-   m_dMaxZ[0] = dLengthZ ;
+   m_dMaxZ[0] = dDimZ ;
    m_dMinZ[1] = 0 ;
-   m_dMaxZ[1] = ( bTriDex ? dLengthX : 0) ;
+   m_dMaxZ[1] = ( bTriDex ? dDimX : 0) ;
    m_dMinZ[2] = 0 ;  
-   m_dMaxZ[2] = ( bTriDex ? dLengthY : 0) ;
+   m_dMaxZ[2] = ( bTriDex ? dDimY : 0) ;
 
    // Tipologia
    m_nShape = BOX ;
@@ -135,13 +135,13 @@ VolZmap::Create( const Point3d& ptO, double dLengthX, double dLengthY, double dL
 
 //----------------------------------------------------------------------------
 bool 
-VolZmap::CreateEmptyMap( const Point3d& ptO, double dLengthX, double dLengthY, double dLengthZ, double dStep, bool bTriDex)
+VolZmap::CreateEmpty( const Point3d& ptO, double dDimX, double dDimY, double dDimZ, double dStep, bool bTriDex)
 {
-  // Controlli sull'ammissibilit� delle dimensioni lineari del grezzo e del passo
-   if ( dStep < EPS_SMALL  ||  dLengthX < EPS_SMALL  ||  dLengthY < EPS_SMALL  ||  dLengthZ < EPS_SMALL)
+  // Controlli sull'ammissibilità delle dimensioni lineari del grezzo e del passo
+   if ( dStep < EPS_SMALL  ||  dDimX < EPS_SMALL  ||  dDimY < EPS_SMALL  ||  dDimZ < EPS_SMALL)
       return false ;
 
-  // Il passo di discretizzazione non pu� essere inferiore a 100 * EPS_SMALL
+  // Il passo di discretizzazione non può essere inferiore a 100 * EPS_SMALL
    m_dStep = max( dStep, 100 * EPS_SMALL) ;
 
   // Aggiorno la dimensione della mappa 1 o 3
@@ -155,8 +155,8 @@ VolZmap::CreateEmptyMap( const Point3d& ptO, double dLengthX, double dLengthY, d
    m_MapFrame.Set( ptO, X_AX, Y_AX, Z_AX) ;
 
   // Definisco i vettori dei limiti su indici
-   m_nNx[0] = max( int( ( dLengthX + EPS_SMALL) / m_dStep + 0.5), 1) ;
-   m_nNy[0] = max( int( ( dLengthY + EPS_SMALL) / m_dStep + 0.5), 1) ;
+   m_nNx[0] = max( int( ( dDimX + EPS_SMALL) / m_dStep + 0.5), 1) ;
+   m_nNy[0] = max( int( ( dDimY + EPS_SMALL) / m_dStep + 0.5), 1) ;
 
   // Numero di componenti connesse 
    m_nConnectedCompoCount = 1 ;
@@ -164,7 +164,7 @@ VolZmap::CreateEmptyMap( const Point3d& ptO, double dLengthX, double dLengthY, d
   // Se tridexel
    if ( bTriDex) {
       m_nNx[1] = m_nNy[0] ;
-      m_nNy[1] = max( int( ( dLengthZ + EPS_SMALL) / m_dStep + 0.5), 1) ;
+      m_nNy[1] = max( int( ( dDimZ + EPS_SMALL) / m_dStep + 0.5), 1) ;
       m_nNx[2] = m_nNy[1] ;
       m_nNy[2] = m_nNx[0] ;
    }
@@ -192,11 +192,11 @@ VolZmap::CreateEmptyMap( const Point3d& ptO, double dLengthX, double dLengthY, d
 
   // Definizione delle limitazioni iniziali in Z per ogni mappa
    m_dMinZ[0] = 0 ;
-   m_dMaxZ[0] = dLengthZ ;
+   m_dMaxZ[0] = dDimZ ;
    m_dMinZ[1] = 0 ;
-   m_dMaxZ[1] = ( bTriDex ? dLengthX : 0) ;
+   m_dMaxZ[1] = ( bTriDex ? dDimX : 0) ;
    m_dMinZ[2] = 0 ;  
-   m_dMaxZ[2] = ( bTriDex ? dLengthY : 0) ;
+   m_dMaxZ[2] = ( bTriDex ? dDimY : 0) ;
 
    // Tipologia
    m_nShape = GENERIC ;
@@ -597,6 +597,98 @@ VolZmap::CreateMapPart( int nMap, int nInfI, int nSupI, int nInfJ, int nSupJ, co
    return true ;
 }
 
+bool
+VolZmap::AddMapPart( int nMap, int nInfI, int nSupI, int nInfJ, int nSupJ, const Vector3d& vtLen, const Point3d& ptMapOrig,
+                     const ISurfTriMesh& Surf, IntersParLinesSurfTm& intPLSTM)
+{
+  // controllo sui parametri
+   if ( nMap < 0  ||  nMap > 2  ||
+        nInfI < 0  ||  nInfI > m_nNx[nMap]  ||
+        nSupI < 0  ||  nSupI > m_nNx[nMap]  ||
+        nInfJ < 0  ||  nInfJ > m_nNy[nMap]  ||
+        nSupJ < 0  ||  nSupJ > m_nNy[nMap])
+      return false ;
+
+  // determinazione e ridimensionamento dei dexel interni alla trimesh
+   for ( int i = nInfI ; i < nSupI ; ++ i) {
+      for ( int j = nInfJ ; j < nSupJ ; ++ j) {
+
+        // definisco la retta da intersecare con la trimesh
+         double dX = ( i + 0.5) * m_dStep ;
+         double dY = ( j + 0.5) * m_dStep ;
+         Point3d ptP0( dX, dY, 0) ;
+
+        // intersezioni della retta con la TriMesh
+         ILSIVECTOR IntersectionResults ;
+         intPLSTM.GetInters( ptP0, vtLen.v[(nMap+2)%3], IntersectionResults) ;
+
+        // rimuovo le intersezioni in eccesso
+         for ( int nI = 0 ; nI < int( IntersectionResults.size()) - 3 ; ++ nI) {
+            int nJ = nI + 1 ; // prima successiva
+            int nK = nJ + 1 ; // seconda successiva
+            int nT = nK + 1 ; // terza successiva
+           // determino i segni delle 4 intersezioni tra la linea e il trangolo della TriMesh
+            int nSgnI = IntersectionResults[nI].dCosDN > EPS_SMALL ? 1 : IntersectionResults[nI].dCosDN > -EPS_SMALL ? 0 : - 1 ;
+            int nSgnJ = IntersectionResults[nJ].dCosDN > EPS_SMALL ? 1 : IntersectionResults[nJ].dCosDN > -EPS_SMALL ? 0 : - 1 ;
+            int nSgnK = IntersectionResults[nK].dCosDN > EPS_SMALL ? 1 : IntersectionResults[nK].dCosDN > -EPS_SMALL ? 0 : - 1 ;
+            int nSgnT = IntersectionResults[nT].dCosDN > EPS_SMALL ? 1 : IntersectionResults[nT].dCosDN > -EPS_SMALL ? 0 : - 1 ;
+           // parametri dell'intersezione sulla linea
+            double dUJ = IntersectionResults[nJ].dU ;
+            double dUK = IntersectionResults[nK].dU ;
+           // controllo coerenza con segni...
+            if ( nSgnI != 0  &&  nSgnI == nSgnJ  &&
+                 nSgnK != 0  &&  nSgnK == nSgnT  &&
+                 nSgnI == - nSgnT  &&
+                 abs( dUJ - dUK) < EPS_SMALL) {
+              // ... ed elimino le intersezioni in eccesso...
+               IntersectionResults.erase( IntersectionResults.begin() + nK) ;
+               IntersectionResults.erase( IntersectionResults.begin() + nJ) ;
+            }
+         }
+
+         int nInt = int( IntersectionResults.size()) ; // numero di intersezioni valide
+         bool bInside = false ; // Flag entrata/uscita per tratto di retta
+         Point3d ptIn ; Vector3d vtInN ;
+
+        // per ogni intersezione valida trovata...
+         for ( int k = 0 ; k < nInt ; ++ k) {
+           // ricavo il tipo di intersezione
+            int nIntType = IntersectionResults[k].nILTT ;
+           // se c'è intersezione
+            if ( nIntType != ILTT_NO) {
+              // ricavo il cos tra i vettori ( normale del triangolo e tangente alla retta)
+               double dCos = IntersectionResults[k].dCosDN ;
+
+              // se entro nella superficie trimesh...
+               if ( dCos < - EPS_SMALL) {
+                  ptIn = IntersectionResults[k].ptI ; // punto di intersezione
+                  int nT = IntersectionResults[k].nT ; // triangolo di interesse
+                  int nF = Surf.GetFacetFromTria( nT) ; // faccia di interesse
+                  Surf.GetFacetNormal( nF, vtInN) ;
+                  bInside = true ; // entrata
+               }
+              // ...se esco dalla superficie trimesh ( prima sono per forza entrato)
+               else if ( dCos > EPS_SMALL  &&  bInside) {
+                  Point3d ptOut = IntersectionResults[k].ptI ; // punto di intersezione
+                  int nT = IntersectionResults[k].nT ;   // triangolo di interesse
+                  int nF = Surf.GetFacetFromTria( nT) ;  // faccia di interesse
+                  Vector3d vtOutN ; Surf.GetFacetNormal( nF, vtOutN) ; // vettore d'uscita
+
+                 // Aggiungo un tratto al dexel
+                  AddIntervals( nMap, i, j,
+                                ptIn.v[(nMap+2)%3] - ptMapOrig.v[(nMap+2)%3],
+                                ptOut.v[(nMap+2)%3] - ptMapOrig.v[(nMap+2)%3],
+                                vtInN, vtOutN, 0, true) ;
+                  bInside = false ; // uscita
+               }
+            }
+         }
+      }
+   }
+
+   return true ;
+}
+
 //----------------------------------------------------------------------------
 bool
 VolZmap::CreateFromTriMesh( const ISurfTriMesh& Surf, double dStep, bool bTriDex)

VolZmapGraphics.cpp

@@ -3327,12 +3327,12 @@ VolZmap::IsThereMat( int nI, int nJ, int nK) const
    nK *= m_nDexVoxRatio ; 
 
   // Se l'indice è alla frontiera del reticolo non vi è materiale
-   if ( nI <= - 1  ||  nI >= int( m_nNx[0])  || 
+   if ( nI <= - 1  ||  nI >= int( m_nNx[0])  ||
         nJ <= - 1  ||  nJ >= int( m_nNy[0])  ||
-        nK <= - 1  ||  nK >= int( m_nNy[1]))  
+        nK <= - 1  ||  nK >= int( m_nNy[1]))
       return false ;
 
-  // ciclo sulle griglie 
+  // ciclo sulle griglie
    int nCount = 0 ;
    int nMinPos[3] = { -1, -1, -1} ;
    int nMinIndex[3] ;
@@ -3361,18 +3361,22 @@ VolZmap::IsThereMat( int nI, int nJ, int nK) const
       }
      // verifica spillone su vertice
       double dMinDist = INFINITO ;
-      int nIndex = 0 ;   
+      int nIndex = 0 ;
       int nPos = nGrJ * m_nNx[nGrid] + nGrI ;
       nDexSize[nGrid] = int( m_Values[nGrid][nPos].size()) ;
+     // scorro i sotto-interalli dello spillone
       while ( nIndex < nDexSize[nGrid]) {
+        // distanza tra la "Z" attuale e il parametro minimo e massimo dell'intervallo nIndex-esimo
          double dDistInf = dZ[nGrid] - m_Values[nGrid][nPos][nIndex].dMin + 2 * EPS_SMALL ;
          double dDistSup = dZ[nGrid] - m_Values[nGrid][nPos][nIndex].dMax - 2 * EPS_SMALL ;
+        // se "Z" attuale a cavallo tra queste due distanze...
          if ( dDistInf > 0.  &&  dDistSup < 0.) {
-            nMinIndex[nGrid] = nIndex ;
+            nMinIndex[nGrid] = nIndex ;      // aggiorno l'indice
             ++ nCount ;
-            bInterOnNode[nGrid] = true ;
+            bInterOnNode[nGrid] = true ;    // flag T per griglia zero
             break ;
          }
+        // se "Z" attuale tutta sopra o tutta sotto
          else {
             double dDist = min( abs( dDistInf), abs( dDistSup)) ;
             if ( dDist < dMinDist) {
@@ -3380,16 +3384,22 @@ VolZmap::IsThereMat( int nI, int nJ, int nK) const
                nMinIndex[nGrid] = nIndex ;
                dMinDist = dDist ;
             }
-         } 
+         }
+        // sotto-intervallo successivo
          nIndex += 1 ;
       }
-   }
-   if ( nCount == 3)
+   } // fine ciclo sulle griglie
+
+   if ( nCount == 3) // ... se interno a tutte e 3 le griglie, allora c'è materiale ...
       return true ;
-   else if ( nCount == 2) {
+   else if ( nCount == 2) { // ... se interno solo a 2 griglie ...
+     // recupero la griglia sulla quale è esterno
       int nGrid = ( bInterOnNode[0] ? ( bInterOnNode[1] ? 2 : 1) : 0) ;
+     // se tale griglia non ha sotto-intervalli allora non c'è materiale
       if ( nDexSize[nGrid] == 0)
          return false ;
+     // se il valore in "Z" dello spillone è vicino al punto ( 1/10 dello step), aggiorno il
+     // parametro minimo e massimo e considero la presenza di materiale
       if ( dZ[nGrid] > m_Values[nGrid][nMinPos[nGrid]][nMinIndex[nGrid]].dMin - 0.1 * m_dStep  &&
            dZ[nGrid] < m_Values[nGrid][nMinPos[nGrid]][nMinIndex[nGrid]].dMax + 0.1 * m_dStep) {
          double dDistInf = abs( dZ[nGrid] - m_Values[nGrid][nMinPos[nGrid]][nMinIndex[nGrid]].dMin) ;
@@ -3401,9 +3411,9 @@ VolZmap::IsThereMat( int nI, int nJ, int nK) const
          return true ;
       }
       else
-         return false ;        
+         return false ;
    }
-   else
+   else // ... se invece interno a 1 o a nessuna delle griglie, allora non c'è materiale
       return false ;
 }
 

VolZmapVolume.cpp

@@ -30,7 +30,8 @@ using namespace std ;
 //----------------------------------------------------------------------------
 bool
 VolZmap::SubtractIntervals( int nGrid, int nI, int nJ,
-                            double dMin, double dMax, const Vector3d& vtNMin, const Vector3d& vtNMax, int nToolNum)
+                            double dMin, double dMax, const Vector3d& vtNMin, const Vector3d& vtNMax,
+                            int nToolNum, bool bSkipSwap)
 {
   // Controllo che dMin e dMax non siano quasi coincidenti
    if ( abs( dMax - dMin) < EPS_ZERO)
@@ -54,13 +55,13 @@ VolZmap::SubtractIntervals( int nGrid, int nI, int nJ,
    }
 
   // Riporto le coordinate cicliche delle normali nell'ordine di partenza (da griglia a sistema intrinseco)
-   if ( nGrid == 1) {
+   if ( !bSkipSwap  &&  nGrid == 1) {
       swap( vtNmi.x, vtNmi.z) ;
       swap( vtNmi.y, vtNmi.z) ;
       swap( vtNma.x, vtNma.z) ;
       swap( vtNma.y, vtNma.z) ;
    }
-   else if ( nGrid == 2) {
+   else if ( !bSkipSwap  &&  nGrid == 2) {
       swap( vtNmi.y, vtNmi.z) ;
       swap( vtNmi.x, vtNmi.z) ;
       swap( vtNma.y, vtNma.z) ;
@@ -250,15 +251,16 @@ VolZmap::SubtractIntervals( int nGrid, int nI, int nJ,
 //----------------------------------------------------------------------------
 bool
 VolZmap::AddIntervals( int nGrid, int nI, int nJ,
-                       double dMin, double dMax, const Vector3d& vtNMin, const Vector3d& vtNMax, int nToolNum)
+                       double dMin, double dMax, const Vector3d& vtNMin, const Vector3d& vtNMax,
+                       int nToolNum, bool bSkipSwap)
 {
   // Controllo che il numero di griglia sia entro i limiti
-   if ( nGrid < 0 || nGrid > 2)
+   if ( nGrid < 0  ||  nGrid > 2)
       return false ;
 
-   // Controllo che indici nI, nJ siano entro i limiti
-   if ( nI < 0 && nI >= m_nNx[nGrid] && 
-        nJ < 0 && nJ >= m_nNy[nGrid])
+  // Controllo che indici nI, nJ siano entro i limiti
+   if ( nI < 0  &&  nI >= m_nNx[nGrid]  && 
+        nJ < 0  &&  nJ >= m_nNy[nGrid])
       return false ;
 
   // Controllo che dMin < dMax
@@ -272,33 +274,37 @@ VolZmap::AddIntervals( int nGrid, int nI, int nJ,
   // Restringo minimo e massimo entro i limiti della mappa
    if ( dMin < m_dMinZ[nGrid]) {
       dMin = m_dMinZ[nGrid] ;
-      vtNmi = - Z_AX ;
+      if ( ! bSkipSwap)
+         vtNmi = - Z_AX ;
    }
    else if ( dMin > m_dMaxZ[nGrid]) {
       dMin = m_dMaxZ[nGrid] ;
-      vtNmi = - Z_AX ;   
+      if ( ! bSkipSwap)
+         vtNmi = - Z_AX ;
    }
    if ( dMax < m_dMinZ[nGrid]) {
       dMax = m_dMinZ[nGrid] ;
-      vtNma = Z_AX ;
+      if ( ! bSkipSwap)
+         vtNma = Z_AX ;
    }
    else if ( dMax > m_dMaxZ[nGrid]) {
       dMax = m_dMaxZ[nGrid] ;
-      vtNma = Z_AX ;   
+      if ( ! bSkipSwap)
+         vtNma = Z_AX ;
    }
 
-   // Controllo che dMin e dMax non siano quasi coincidenti
+  // Controllo che dMin e dMax non siano quasi coincidenti
    if ( abs( dMax - dMin) < EPS_SMALL)
       return true ;
 
   // Riporto le coordinate cicliche nell'ordine di partenza
-   if ( nGrid == 1) {
+   if ( !bSkipSwap  &&  nGrid == 1) {
       swap( vtNmi.x, vtNmi.z) ;
       swap( vtNmi.y, vtNmi.z) ;
       swap( vtNma.x, vtNma.z) ;
       swap( vtNma.y, vtNma.z) ;
    }
-   else if ( nGrid == 2) {
+   else if ( !bSkipSwap  &&  nGrid == 2) {
       swap( vtNmi.y, vtNmi.z) ;
       swap( vtNmi.x, vtNmi.z) ;
       swap( vtNma.y, vtNma.z) ;
@@ -306,13 +312,13 @@ VolZmap::AddIntervals( int nGrid, int nI, int nJ,
    }
 
    // Calcolo nPos
-   unsigned int nPos = nJ * m_nNx[nGrid] + nI ;
+   int nPos = nJ * m_nNx[nGrid] + nI ;
    vector<Data>& vDexel = m_Values[nGrid][nPos] ;
 
    bool bModified = false ;
 
-   // Non esistono segmenti 
-   if ( int( vDexel.size()) == 0) {
+  // Non esistono segmenti 
+   if ( vDexel.empty()) {
    
       vDexel.emplace_back() ;
       vDexel.back().dMin = dMin ;
@@ -326,10 +332,10 @@ VolZmap::AddIntervals( int nGrid, int nI, int nJ,
 
       bModified = true ;
    }
-   // Esiste almeno un segmento
+  // Esiste almeno un segmento
    else {
-      // Cerco l'ultimo intervallo a sinistra e l'ultimo intervallo a destra
-      // di quello da aggiungere, che non interferiscono con quest'ultimo.
+     // Cerco l'ultimo intervallo a sinistra e l'ultimo intervallo a destra
+     // di quello da aggiungere, che non interferiscono con quest'ultimo.
       auto itLastLeft = vDexel.end() ;
       auto itFirstRight = vDexel.end() ;
       for ( auto it = vDexel.begin() ; it != vDexel.end() ; ++ it) {
@@ -338,12 +344,12 @@ VolZmap::AddIntervals( int nGrid, int nI, int nJ,
          if ( dMax < it->dMin - EPS_SMALL  &&  itFirstRight == vDexel.end())
             itFirstRight = it ;
       }
-      // Esistono intervalli a sinistra.
+     // Esistono intervalli a sinistra.
       if ( itLastLeft != vDexel.end()) {
-         // Intervallo successivo all'ultimo a sinistra
+        // Intervallo successivo all'ultimo a sinistra
          auto itNextToLastLeft = itLastLeft ;
          ++ itNextToLastLeft ;
-         // Il successivo non esiste.
+        // Il successivo non esiste.
          if ( itNextToLastLeft == vDexel.end()) {
             // Aggiungo il nuovo semgento
             vDexel.emplace_back() ;
@@ -356,11 +362,11 @@ VolZmap::AddIntervals( int nGrid, int nI, int nJ,
             //m_Values[nGrid][nPos].back().nCompo = ;
             bModified = true ;
          }
-         // Il successivo esiste.
+        // Il successivo esiste.
          else {
-            // Il successivo è il primo a destra.
+           // Il successivo è il primo a destra.
             if ( itNextToLastLeft == itFirstRight) {
-               // Inserisco nuovo segmento-
+              // Inserisco nuovo segmento-
                Data NewSegment ;
                NewSegment.dMin = dMin ;
                NewSegment.dMax = dMax ;
@@ -373,18 +379,18 @@ VolZmap::AddIntervals( int nGrid, int nI, int nJ,
                bModified = true ;
             }
             else {
-               // Il successivo non esce a sinistra da quello da aggiungere.
+              // Il successivo non esce a sinistra da quello da aggiungere.
                if ( itNextToLastLeft->dMin > dMin + EPS_SMALL) {
                   itNextToLastLeft->dMin = dMin ;
                   itNextToLastLeft->vtMinN = vtNmi ;
                   itNextToLastLeft->nToolMin = nToolNum ;
                }
-               // Cerco l'ultimo segmento che interferisce con quello da aggiungere.
+              // Cerco l'ultimo segmento che interferisce con quello da aggiungere.
                auto itPrevToFirstRight = vDexel.end() ;
                for ( auto it = itNextToLastLeft ; it != itFirstRight ; ++ it) {
                   itPrevToFirstRight = it ;
                }
-               // L'ultimo che interferisce non esce a destra da quello da aggiungere.
+              // L'ultimo che interferisce non esce a destra da quello da aggiungere.
                if ( itPrevToFirstRight->dMax < dMax - EPS_SMALL) {
                   itNextToLastLeft->dMax = dMax ;
                   itNextToLastLeft->vtMaxN = vtNma ;
@@ -403,18 +409,18 @@ VolZmap::AddIntervals( int nGrid, int nI, int nJ,
             }
          }
       }
-      // Non esistono neanche a destra.
+     // Non esistono neanche a destra.
       else if ( itFirstRight == m_Values[nGrid][nPos].end()) {
-         // Il primo intervallo non sporge a sinistra
+        // Il primo intervallo non sporge a sinistra
          if ( vDexel.begin()->dMin > dMin + EPS_SMALL) {
             vDexel.begin()->dMin = dMin ;
             vDexel.begin()->vtMinN = vtNmi ;
             vDexel.begin()->nToolMin = nToolNum ;
             bModified = true ;
          }
-         // L'ultimo intervallo sporge a destra.
+        // L'ultimo intervallo sporge a destra.
          if ( m_Values[nGrid][nPos].back().dMax > dMax + EPS_SMALL) {
-            // Ci sono più segmenti, inglobo tutti nel primo.
+           // Ci sono più segmenti, inglobo tutti nel primo.
             if ( vDexel.back().dMax > vDexel.begin()->dMax + EPS_SMALL) {
                vDexel.begin()->dMax = vDexel.back().dMax ;
                vDexel.begin()->vtMaxN = vDexel.back().vtMaxN ;
@@ -422,7 +428,7 @@ VolZmap::AddIntervals( int nGrid, int nI, int nJ,
                bModified = true ;
             }
          }
-         // L'ultimo intervallo non sporge a destra.
+        // L'ultimo intervallo non sporge a destra.
          else {
             vDexel.begin()->dMax = dMax ;
             vDexel.begin()->vtMaxN = vtNma ;
@@ -431,11 +437,11 @@ VolZmap::AddIntervals( int nGrid, int nI, int nJ,
          }
          vDexel.erase( vDexel.begin() + 1, vDexel.end()) ;
       }
-      // A destra esistono.
+     // A destra esistono.
       else {
-         // Tutti i segmenti sono a destra di qullo da aggiungere.
+        // Tutti i segmenti sono a destra di qullo da aggiungere.
          if ( itFirstRight == vDexel.begin()) {
-            // Inserisco nuovo segmento-
+           // Inserisco nuovo segmento-
             Data NewSegment ;
             NewSegment.dMin = dMin ;
             NewSegment.dMax = dMax ;
@@ -447,19 +453,19 @@ VolZmap::AddIntervals( int nGrid, int nI, int nJ,
             bModified = true ;
          }
          else {
-            // Se il primo segmento non esce a sinistra da quello da aggiungere, cambio l'inizio.
+           // Se il primo segmento non esce a sinistra da quello da aggiungere, cambio l'inizio.
             if ( vDexel.begin()->dMin > dMin + EPS_SMALL) {
                vDexel.begin()->dMin = dMin ;
                vDexel.begin()->vtMinN = vtNmi ;
                vDexel.begin()->nToolMin = nToolNum ;
                bModified = true ;
             }
-            // Cerco l'ultimo segmento che interferisce con quello da aggiungere.
+           // Cerco l'ultimo segmento che interferisce con quello da aggiungere.
             auto itPrevToFirstRight = vDexel.begin() ;
             for ( auto it = m_Values[nGrid][nPos].begin() ; it != itFirstRight ; ++ it) {
                itPrevToFirstRight = it ;
             }
-            // L'ultimo che interferisce non esce a destra da quello da aggiungere.
+           // L'ultimo che interferisce non esce a destra da quello da aggiungere.
             if ( itPrevToFirstRight->dMax < dMax - EPS_SMALL) {
                vDexel.begin()->dMax = dMax ;
                vDexel.begin()->vtMaxN = vtNma ;
@@ -479,24 +485,24 @@ VolZmap::AddIntervals( int nGrid, int nI, int nJ,
       }  
    }
 
-   // Se nessuna modifica, esco
+  // Se nessuna modifica, esco
    if ( ! bModified)
       return true ;
 
-   // Imposto ricalcolo della grafica
+  // Imposto ricalcolo della grafica
    m_OGrMgr.Reset() ;
-   // Imposto forma generica 
+  // Imposto forma generica 
    m_nShape = GENERIC ;    
-   // Imposto ricalcolo numero di componenti connesse
+  // Imposto ricalcolo numero di componenti connesse
    m_nConnectedCompoCount = - 1 ;
 
-   // Passo da indici di dexel a indici di voxel
+  // Passo da indici di dexel a indici di voxel
    nI /= m_nDexVoxRatio ;
    nJ /= m_nDexVoxRatio ;
 
-   // Determino quali blocchi sono stati modificati
+  // Determino quali blocchi sono stati modificati
    if ( nGrid == 0) {
-      // Voxel lungo X
+     // Voxel lungo X
       int nXStop = 1 ;
       int nXBlock[2] ;
       nXBlock[0] = min( nI / m_nVoxNumPerBlock, m_nFracLin[0] - 1) ;
@@ -504,7 +510,7 @@ VolZmap::AddIntervals( int nGrid, int nI, int nJ,
          nXBlock[1] = nXBlock[0] - 1 ;
          ++ nXStop ;
       }
-      // Voxel lungo Y
+     // Voxel lungo Y
       int nYStop = 1 ;
       int nYBlock[2] ;
       nYBlock[0] = min( nJ / m_nVoxNumPerBlock, m_nFracLin[1] - 1) ;
@@ -512,13 +518,13 @@ VolZmap::AddIntervals( int nGrid, int nI, int nJ,
          nYBlock[1] = nYBlock[0] - 1 ;
          ++ nYStop ;
       }
-      // Voxel lungo Z
+     // Voxel lungo Z
       int nVoxNumZ = int( m_nNy[1] / m_nDexVoxRatio + ( m_nNy[1] % m_nDexVoxRatio == 0 ? 1 : 2)) ;  
       int nMinK = Clamp( int( floor( ( ( dMin - 0.5 * m_dStep) / ( m_nDexVoxRatio * m_dStep) - EPS_SMALL))), 0, nVoxNumZ - 2) ;
       int nMaxK = Clamp( int( floor( ( ( dMax + 0.5 * m_dStep) / ( m_nDexVoxRatio * m_dStep) + EPS_SMALL))), 0, nVoxNumZ - 2) ;
       int nMinZBlock = ( m_nMapNum == 1 ? 0 : Clamp( nMinK / int( m_nVoxNumPerBlock), 0, int( m_nFracLin[2] - 1))) ;
       int nMaxZBlock = min( int( m_nFracLin[2] - 1), nMaxK / int( m_nVoxNumPerBlock)) ;
-      // Assegno flag ai voxel
+     // Assegno flag ai voxel
       for ( int tI = 0 ; tI < nXStop ; ++ tI) {
          for ( int tJ = 0 ; tJ < nYStop ; ++ tJ) {
             for ( int k = nMinZBlock ; k <= nMaxZBlock ; ++ k) {
@@ -530,7 +536,7 @@ VolZmap::AddIntervals( int nGrid, int nI, int nJ,
    }
 
    else if ( nGrid == 1) {
-      // Voxel lungo Y
+     // Voxel lungo Y
       int nYStop = 1 ;
       int nYBlock[2] ;
       nYBlock[0] = min( nI / m_nVoxNumPerBlock, m_nFracLin[1] - 1) ;
@@ -538,7 +544,7 @@ VolZmap::AddIntervals( int nGrid, int nI, int nJ,
          nYBlock[1] = nYBlock[0] - 1 ;
          ++ nYStop ;
       }
-      // Voxel lungo Z
+     // Voxel lungo Z
       int nZStop = 1 ;
       int nZBlock[2] ;
       nZBlock[0] = min( nJ / m_nVoxNumPerBlock, m_nFracLin[2] - 1) ;
@@ -546,13 +552,13 @@ VolZmap::AddIntervals( int nGrid, int nI, int nJ,
          nZBlock[1] = nZBlock[0] - 1 ;
          ++ nZStop ;
       }
-      // Voxel lungo X
+     // Voxel lungo X
       int nVoxNumX = int( m_nNx[0] / m_nDexVoxRatio + ( m_nNx[0] % m_nDexVoxRatio == 0 ? 1 : 2)) ;
       int nMinI = Clamp( int( floor( ( ( dMin - 0.5 * m_dStep) / ( m_nDexVoxRatio * m_dStep) - EPS_SMALL))), 0, nVoxNumX - 2) ;
       int nMaxI = Clamp( int( floor( ( ( dMax + 0.5 * m_dStep) / ( m_nDexVoxRatio * m_dStep) + EPS_SMALL))), 0, nVoxNumX - 2) ;
       int nMinXBlock = Clamp( nMinI / int( m_nVoxNumPerBlock), 0, int( m_nFracLin[0] - 1)) ;
       int nMaxXBlock = min( int( m_nFracLin[0] - 1), nMaxI / int( m_nVoxNumPerBlock)) ;
-      // Assegno flag ai voxel
+     // Assegno flag ai voxel
       for ( int tI = 0 ; tI < nYStop ; ++ tI) {
          for ( int tJ = 0 ; tJ < nZStop ; ++ tJ) {
             for ( int k = nMinXBlock ; k <= nMaxXBlock ; ++ k) {
@@ -564,7 +570,7 @@ VolZmap::AddIntervals( int nGrid, int nI, int nJ,
    }
 
    else if ( nGrid == 2) {
-      // Voxel lungo X
+     // Voxel lungo X
       int nXStop = 1 ;
       int nXBlock[2] ;
       nXBlock[0] = min( nJ / m_nVoxNumPerBlock, m_nFracLin[0] - 1) ;
@@ -572,7 +578,7 @@ VolZmap::AddIntervals( int nGrid, int nI, int nJ,
          nXBlock[1] = nXBlock[0] - 1 ;
          ++ nXStop ;
       }
-      // Voxel lungo Z
+     // Voxel lungo Z
       int nZStop = 1 ;
       int nZBlock[2] ;
       nZBlock[0] = min( nI / m_nVoxNumPerBlock, m_nFracLin[2] - 1) ;
@@ -580,13 +586,13 @@ VolZmap::AddIntervals( int nGrid, int nI, int nJ,
          nZBlock[1] = nZBlock[0] - 1 ;
          ++ nZStop ;
       }
-      // Voxel lungo Y 
+     // Voxel lungo Y 
       int nVoxNumY = int( m_nNy[0] / m_nDexVoxRatio + ( m_nNy[0] % m_nDexVoxRatio == 0 ? 1 : 2)) ; 
       int nMinJ = Clamp( int( floor( ( ( dMin - 0.5 * m_dStep) / ( m_nDexVoxRatio * m_dStep) - EPS_SMALL))), 0, nVoxNumY - 2) ; 
       int nMaxJ = Clamp( int( floor( ( ( dMax + 0.5 * m_dStep) / ( m_nDexVoxRatio * m_dStep) + EPS_SMALL))), 0, nVoxNumY - 2) ; 
       int nMinYBlock = Clamp( nMinJ / int( m_nVoxNumPerBlock), 0, int( m_nFracLin[1] - 1)) ;
       int nMaxYBlock = min( int( m_nFracLin[1] - 1), nMaxJ / int( m_nVoxNumPerBlock)) ;
-      // Assegno flag ai voxel
+     // Assegno flag ai voxel
       for ( int tI = 0 ; tI < nZStop ; ++ tI) {
          for ( int tJ = 0 ; tJ < nXStop ; ++ tJ) {
             for ( int k = nMinYBlock ; k <= nMaxYBlock ; ++ k) {
@@ -600,6 +606,213 @@ VolZmap::AddIntervals( int nGrid, int nI, int nJ,
    return true ;
 }
 
+//----------------------------------------------------------------------------
+bool 
+VolZmap::AddSubIntervalInVoxel( VolZmap* VolZmapRef, int nGrid, int nI, int nJ, int nK,
+                                double& dMin, double& dMax, Vector3d& vtMin, Vector3d& vtMax)
+{
+  // se non Tridex, esco
+   if ( ! IsTriDexel())
+      return true ;
+  
+  // Controllo che il numero di griglia sia entro i limiti
+   if ( nGrid < 0  ||  nGrid > 2)
+      return false ;
+
+  // Controllo che indici nI, nJ siano entro i limiti
+   if ( nI < 0  &&  nI >= m_nNx[nGrid]  &&
+        nJ < 0  &&  nJ >= m_nNy[nGrid])
+      return false ;
+
+  // valori di default 
+   vector<double> vTdMin = { -1., -1., -1., -1.} ;
+   vector<double> vTdMax = { -1., -1., -1., -1.} ;
+   VCT3DVECTOR vtNMin = { V_INVALID, V_INVALID, V_INVALID, V_INVALID} ;
+   VCT3DVECTOR vtNMax = { V_INVALID, V_INVALID, V_INVALID, V_INVALID} ;
+
+  // se esiste un precedente lungo x-locale
+   if ( nI != 0) {
+      int nPos = nJ * m_nNx[nGrid] + ( nI - 1) ; // recupero posizione dexel
+     // cerco l'intervallo corretto sullo Zmap di riferimento
+      for ( int it = 0 ; it < int( VolZmapRef->m_Values[nGrid][nPos].size()) ; ++ it) {
+         if ( VolZmapRef->m_Values[nGrid][nPos][it].dMax + EPS_SMALL > nK * m_dStep  &&
+              VolZmapRef->m_Values[nGrid][nPos][it].dMax - EPS_SMALL < ( nK + 1) * m_dStep) {
+            vtNMax[0] = VolZmapRef->m_Values[nGrid][nPos][it].vtMaxN ;
+            vTdMax[0] = VolZmapRef->m_Values[nGrid][nPos][it].dMax ;
+         }
+         if ( VolZmapRef->m_Values[nGrid][nPos][it].dMin + EPS_SMALL > nK * m_dStep  &&
+              VolZmapRef->m_Values[nGrid][nPos][it].dMin - EPS_SMALL < ( nK + 1) * m_dStep) {
+            vtNMin[0] = VolZmapRef->m_Values[nGrid][nPos][it].vtMinN ;
+            vTdMin[0] = VolZmapRef->m_Values[nGrid][nPos][it].dMin ;
+         }
+      }
+   }
+  // se esiste il successivo lungo x-locale
+   if ( nI != m_dMaxZ[( nGrid + 2) % 3]) {
+      int nPos = nJ * m_nNx[nGrid] + ( nI + 1) ; // recupero posizione dexel
+      for ( int it = 0 ; it < int( VolZmapRef->m_Values[nGrid][nPos].size()) ; ++ it) {
+        // cerco l'intervallo corretto sullo Zmap di riferimento
+         if ( VolZmapRef->m_Values[nGrid][nPos][it].dMax + EPS_SMALL > nK * m_dStep  &&
+              VolZmapRef->m_Values[nGrid][nPos][it].dMax - EPS_SMALL < ( nK + 1) * m_dStep) {
+            vtNMax[1] = VolZmapRef->m_Values[nGrid][nPos][it].vtMaxN ;
+            vTdMax[1] = VolZmapRef->m_Values[nGrid][nPos][it].dMax ;
+         }
+         if ( VolZmapRef->m_Values[nGrid][nPos][it].dMin + EPS_SMALL > nK * m_dStep  &&
+              VolZmapRef->m_Values[nGrid][nPos][it].dMin - EPS_SMALL < ( nK + 1) * m_dStep) {
+            vtNMin[1] = VolZmapRef->m_Values[nGrid][nPos][it].vtMinN ;
+            vTdMin[1] = VolZmapRef->m_Values[nGrid][nPos][it].dMin ;
+         }
+      }
+   }
+  // se esiste il precedente lungo y-locale
+   if ( nJ != 0) {
+      int nPos = ( nJ - 1) * m_nNx[nGrid] + nI ; // recupero posizione dexel
+      for ( int it = 0 ; it < int( VolZmapRef->m_Values[nGrid][nPos].size()) ; ++ it) {
+        // cerco l'intervallo corretto sullo Zmap di riferimento
+         if ( VolZmapRef->m_Values[nGrid][nPos][it].dMax + EPS_SMALL > nK * m_dStep  &&
+              VolZmapRef->m_Values[nGrid][nPos][it].dMax - EPS_SMALL < ( nK + 1) * m_dStep) {
+            vtNMax[2] = VolZmapRef->m_Values[nGrid][nPos][it].vtMaxN ;
+            vTdMax[2] = VolZmapRef->m_Values[nGrid][nPos][it].dMax ;
+         }
+         if ( VolZmapRef->m_Values[nGrid][nPos][it].dMin + EPS_SMALL > nK * m_dStep  &&
+              VolZmapRef->m_Values[nGrid][nPos][it].dMin - EPS_SMALL < ( nK + 1) * m_dStep) {
+            vtNMin[2] = VolZmapRef->m_Values[nGrid][nPos][it].vtMinN ;
+            vTdMin[2] = VolZmapRef->m_Values[nGrid][nPos][it].dMin ;
+         }
+      }
+   }
+  // se esiste il successivo lungo y-locale
+   if ( nJ != m_dMaxZ[( nGrid + 1) % 3]) {
+      int nPos = ( nJ + 1) * m_nNx[nGrid] + nI ; // recupero posizione dexel
+     // cerco l'intervallo corretto sullo Zmap di riferimento
+      for ( int it = 0 ; it < int( VolZmapRef->m_Values[nGrid][nPos].size()) ; ++ it) {
+         if ( VolZmapRef->m_Values[nGrid][nPos][it].dMax + EPS_SMALL > nK * m_dStep  &&
+              VolZmapRef->m_Values[nGrid][nPos][it].dMax - EPS_SMALL < ( nK + 1) * m_dStep) {
+            vtNMax[3] = VolZmapRef->m_Values[nGrid][nPos][it].vtMaxN ;
+            vTdMax[3] = VolZmapRef->m_Values[nGrid][nPos][it].dMax ;
+         }
+         if ( VolZmapRef->m_Values[nGrid][nPos][it].dMin + EPS_SMALL > nK * m_dStep  &&
+              VolZmapRef->m_Values[nGrid][nPos][it].dMin - EPS_SMALL < ( nK + 1) * m_dStep) {
+            vtNMin[3] = VolZmapRef->m_Values[nGrid][nPos][it].vtMinN ;
+            vTdMin[3] = VolZmapRef->m_Values[nGrid][nPos][it].dMin ;
+         }
+      }
+   }
+
+  // scelgo le normali che si discostano il meno possibile dalla normale della faccia del voxel
+
+  // analisi dei massimi e delle normali ---------------------------
+  // 1) angolo minimo tra la normale trovata la faccia del voxel
+   double dMinAngle = ANG_FULL ;
+  // 2) vettore di riferimento per la direzione della normale
+   vtMax = ( nGrid == 0 ? m_MapFrame.VersZ() :
+           ( nGrid == 1 ? m_MapFrame.VersX() : m_MapFrame.VersY())) ;
+   Vector3d vtRef = vtMax ;
+  // 3) determino il massimo per questo intervallo
+   dMax = ( nK + 1) * m_dStep ;
+   for ( int i = 0 ; i < 4 ; ++ i) { // scorro le normali
+      if ( vtNMax[i].IsValid()) { // se normale trovata, quindi valida...
+         double dCurrAngle ; // angolo corrente tra la normale della TriMesh e quella della faccia del voxel
+         vtNMax[i].GetAngle( vtRef, dCurrAngle) ;
+         if ( abs( dCurrAngle) < dMinAngle) { // se angolo minore del minimo trovato...
+           // aggiorno i parametri
+            dMinAngle = dCurrAngle ;
+            vtMax = vtNMax[i] ;
+            dMax = vTdMax[i] ;
+         }
+      }
+   }
+
+  // analisi dei minimi e delle normali ---------------------------
+   dMinAngle = ANG_FULL ;
+   vtRef.Invert() ;
+   dMin = nK * m_dStep ;
+   vtMin = vtRef ;
+   for ( int i = 0 ; i < 4 ; ++ i) {
+      if ( vtNMin[i].IsValid()) {
+         double dCurrAngle ;
+         vtNMin[i].GetAngle( vtRef, dCurrAngle) ;
+         if ( abs( dCurrAngle) < dMinAngle) {
+            dMinAngle = dCurrAngle  ;
+            vtMin = vtNMin[i] ;
+            dMin = vTdMin[i] ;
+         }
+      }
+   }
+
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+bool
+VolZmap::AddMissingIntervalsInVoxel( VolZmap* VolZmapRef, int nGrid, int nI, int nJ, double dZ, double dToler,
+                                     Vector3d vtToolMin, Vector3d vtToolMax, int nToolNum)
+{
+
+  // se non Tridex, esco
+   if ( ! IsTriDexel())
+      return true ;
+  
+  // Controllo che il numero di griglia sia entro i limiti
+   if ( nGrid < 0  ||  nGrid > 2)
+      return false ;
+
+  // Controllo che indici nI, nJ siano entro i limiti
+   if ( nI < 0  &&  nI >= m_nNx[nGrid]  &&
+        nJ < 0  &&  nJ >= m_nNy[nGrid])
+      return false ;
+
+  // passo da indici di dexel a indici di voxel
+   nI /= m_nDexVoxRatio ;
+   nJ /= m_nDexVoxRatio ;
+
+  // numero di voxel nel dexel corrente
+   int nVoxNum = int( m_nNy[(( nGrid+1) % 3)] / m_nDexVoxRatio +
+                 ( m_nNy[(( nGrid+1) % 3)] % m_nDexVoxRatio == 0 ? 1 : 2)) ;
+
+   int nK = 0 ;
+   for ( int i = 0 ; i < nVoxNum ; ++ i) {
+     // controllo se sono nel voxel corrente
+      if ( i * m_dStep < dZ  &&  ( i + 1) * m_dStep > dZ) {
+         nK = i ;
+         break ;
+      }
+   }
+
+  // ----------- griglia successiva -----------
+   {
+      int nMyGrid = ( nGrid + 1) % 3 ;
+      int nMyI = nJ ;
+      int nMyJ = nK ;
+      int nMyK = nI ;
+      double dMyMin ;
+      double dMyMax ;
+      Vector3d vtMyMin ;
+      Vector3d vtMyMax ;
+      AddSubIntervalInVoxel( VolZmapRef, nMyGrid, nMyI, nMyJ, nMyK, dMyMin, dMyMax, vtMyMin, vtMyMax) ;
+      AddIntervals( nMyGrid, nMyI, nMyJ, dMyMin - EPS_SMALL, dMyMax + EPS_SMALL, vtMyMin, vtMyMax,
+                    nToolNum, true) ;
+   }
+
+   // ----------- griglia precedente -----------
+   {
+      int nMyGrid = ( nGrid + 2) % 3 ;
+      int nMyI = nK ;
+      int nMyJ = nI ;
+      int nMyK = nJ ;
+      double dMyMin ;
+      double dMyMax ;
+      Vector3d vtMyMin ;
+      Vector3d vtMyMax ;
+      AddSubIntervalInVoxel( VolZmapRef, nMyGrid, nMyI, nMyJ, nMyK, dMyMin, dMyMax, vtMyMin, vtMyMax) ;
+      AddIntervals( nMyGrid, nMyI, nMyJ, dMyMin - EPS_SMALL, dMyMax + EPS_SMALL, vtMyMin, vtMyMax,
+                    nToolNum, true) ;
+   }
+
+   return true ;
+
+}
+
 // ------------------------- LAVORAZIONI --------------------------------------------------------------------------------------
 
 //----------------------------------------------------------------------------       
@@ -638,9 +851,6 @@ VolZmap::MillingStep( int nCurrTool,
          return false ;
    }
 
-  // Reset numero di parti
-   m_nConnectedCompoCount = - 1 ;
-
   // Punti e vettori descriventi il moto nel sistema intrinseco dello Zmap
    Point3d ptPLs = GetToLoc( ptPs, m_MapFrame) ;
    Point3d ptPLe = GetToLoc( ptPe, m_MapFrame) ;
@@ -682,7 +892,7 @@ VolZmap::MillingGeneralMotionStep( const Point3d& ptPs, const Vector3d& vtDs, co
                                    const Point3d& ptPe, const Vector3d& vtDe, const Vector3d& vtAe)
 {
   // Divido il movimento in tratti con direzione utensile costante
-   const double ANG_STEP = 0.02 ;
+   const double ANG_STEP = 0.04 ;
    double dAngDeg ; vtDs.GetAngle( vtDe, dAngDeg) ;
    int nStepCnt = int( abs( dAngDeg) / ANG_STEP) + 1 ;
    bool bOk = true ;
@@ -837,6 +1047,8 @@ VolZmap::SelectMotion( int nGrid, const Point3d& ptLs, const Point3d& ptLe, cons
             return Conus_ZMilling( nGrid, ptLs, ptLe, vtL) ;
          case Tool::MORTISER :
             return Mrt_Milling( nGrid, ptLs, ptLe, vtL, vtAL) ;
+         case Tool::ADDITIVE :
+            return AddingMotion( nGrid, ptLs, ptLe, vtL) ;
          }
       }
    }
@@ -899,6 +1111,8 @@ VolZmap::SelectMotion( int nGrid, const Point3d& ptLs, const Point3d& ptLe, cons
             return Conus_XYMilling( nGrid, ptLs, ptLe, vtL) ;
          case Tool::MORTISER :
             return Mrt_Milling( nGrid, ptLs, ptLe, vtL, vtAL) ;
+         case Tool::ADDITIVE :
+            return AddingMotion( nGrid, ptLs, ptLe, vtL) ;
          }
       }
    }
@@ -945,6 +1159,8 @@ VolZmap::SelectMotion( int nGrid, const Point3d& ptLs, const Point3d& ptLe, cons
             if ( dSqLLong < EPS_SMALL * EPS_SMALL)
                return Chs_Milling( nGrid, ptLs, ptLe, vtL, vtAL) ;
             break ;
+         case Tool::ADDITIVE :
+            return AddingMotion( nGrid, ptLs, ptLe, vtL) ;
          }
       }
    }
@@ -4942,29 +5158,27 @@ VolZmap::CompBall_Milling( int nGrid, const Point3d& ptLs, const Point3d& ptLe,
 bool
 VolZmap::AddingMotion( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtAx)
 {
-   // Dimensioni lineari dell'utensile
+  // Dimensioni lineari dell'utensile
    double dHei = m_vTool[m_nCurrTool].GetHeigth() ;
    double dRad = m_vTool[m_nCurrTool].GetRadius() ;
    double dCornerRad = m_vTool[m_nCurrTool].GetCornRadius() ;
 
-   // Utensile sfiancato
+  // Utensile sfiancato
    if ( dCornerRad * dCornerRad - 0.25 * dHei * dHei > 0) {
-      AddingGeneral( nGrid, ptS, ptE, vtAx) ;   
+      return AddingGeneral( nGrid, ptS, ptE, vtAx) ;   
    }
-   // Utensile sferico
+  // Utensile sferico
    else if ( dRad - dCornerRad < EPS_SMALL) {
-      AddingSphere( nGrid, ptS - dRad * vtAx, ptE - dRad * vtAx, dRad) ;
+      return AddingSphere( nGrid, ptS - dRad * vtAx, ptE - dRad * vtAx, dRad) ;
    }
-   // Utensile cilindro
+  // Utensile cilindro
    else if ( dCornerRad < EPS_SMALL) {
-      AddingCylinder( nGrid, ptS, ptE, vtAx, dHei, dRad) ;
+      return AddingCylinder( nGrid, ptS, ptE, vtAx, dHei, dRad) ;
    }
-   // Utensile naso di toro
+  // Utensile naso di toro
    else {
-      AddingGeneral( nGrid, ptS, ptE, vtAx) ;
+      return AddingGeneral( nGrid, ptS, ptE, vtAx) ;
    }
-
-   return true ;
 }
 
 //----------------------------------------------------------------------------
@@ -5064,57 +5278,37 @@ VolZmap::AddingCylinder( int nGrid, const Point3d& ptS, const Point3d& ptE, cons
    double dLen1 = vtV1.Len() ;
    vtV1 /= dLen1 ;
 
-   if ( nGrid == 0) {
+   double dMyTol = 0 ;
+   Frame3d CylFrame, PolyFrame ;
+   if ( ! CylFrame.Set( ptS - dHei * vtAx, vtAx))
+      return false ;
+   if ( ! PolyFrame.Set( ptS - ( dHei + dMyTol) * vtAx, vtAx, vtV1))
+      return false ;
+   for ( int i = nStartI ; i <= nEndI ; ++ i) {
+      for ( int j = nStartJ ; j <= nEndJ ; ++ j) {
+         
+         Point3d ptC( ( i + 0.5) * m_dStep, ( j + 0.5) * m_dStep, 0) ;
+                                           
+         Point3d ptInt1, ptInt2 ;
+         Vector3d vtN1, vtN2 ;
 
-      Vector3d vtV2 = Z_AX ^ vtV1 ;
+         if ( IntersLineCylinder( ptC, Z_AX, CylFrame, dHei, dRad, true, true,
+                                  ptInt1, vtN1, ptInt2, vtN2)) {
+            AddIntervals( nGrid, i, j, ptInt1.z, ptInt2.z, - vtN1, - vtN2, CurrTool.GetToolNum()) ;
+         }
 
-      for ( int i = nStartI ; i <= nEndI ; ++ i) {
-         for ( int j = nStartJ ; j <= nEndJ ; ++ j) {
+         if ( IntersLineCylinder( ptC - dLen1 * vtV1, Z_AX, CylFrame, dHei, dRad, true, true,
+                                  ptInt1, vtN1, ptInt2, vtN2)) {
+            AddIntervals( nGrid, i, j, ptInt1.z + dLen1 * vtV1.z, ptInt2.z + dLen1 * vtV1.z, - vtN1, - vtN2, CurrTool.GetToolNum()) ;
+         }
 
-            Vector3d vtStC( ( i + 0.5) * m_dStep - ptS.x, ( j + 0.5) * m_dStep - ptS.y, 0) ;
-            Vector3d vtEnC = vtStC - dLen1 * vtV1 ;
-
-            double dX1 = vtStC * vtV1 ;
-            double dX2 = vtStC * vtV2 ;
-
-            if ( ( dX1 > 0  &&  dX1 < dLen1  &&  abs( dX2) < dRad + EPS_SMALL)  ||
-                 vtStC.SqLen() < dRad * dRad + 2 * dRad * EPS_SMALL  ||
-                 vtEnC.SqLen() < dRad * dRad + 2 * dRad * EPS_SMALL) {
-               AddIntervals( nGrid, i, j, ptS.z - dHei, ptS.z, - Z_AX, Z_AX, CurrTool.GetToolNum()) ;
-            }    
+         if ( IntersLineMyPolyhedron( ptC, Z_AX, PolyFrame, dLen1, 2 * ( dRad + dMyTol), dHei + 2 * dMyTol, 0,
+                                      ptInt1, vtN1, ptInt2, vtN2)) {
+            AddIntervals( nGrid, i, j, ptInt1.z, ptInt2.z, - vtN1, - vtN2, CurrTool.GetToolNum()) ;
          }
       }
    }
-   else {
-      double dMyTol = 0 ;
-      Frame3d CylFrame, PolyFrame ;
-      CylFrame.Set( ptS - dHei * vtAx, vtAx) ;
-      PolyFrame.Set( ptS - ( dHei + dMyTol) * vtAx, vtV1, vtAx ^ vtV1, vtAx) ;
-      for ( int i = nStartI ; i <= nEndI ; ++ i) {
-         for ( int j = nStartJ ; j <= nEndJ ; ++ j) {
-            
-            Point3d ptC( ( i + 0.5) * m_dStep, ( j + 0.5) * m_dStep, 0) ;
-                                              
-            Point3d ptInt1, ptInt2 ;
-            Vector3d vtN1, vtN2 ;
 
-            if ( IntersLineCylinder( ptC, Z_AX, CylFrame, dHei, dRad, true, true,
-                                     ptInt1, vtN1, ptInt2, vtN2)) {
-               AddIntervals( nGrid, i, j, ptInt1.z, ptInt2.z, - vtN1, - vtN2, CurrTool.GetToolNum()) ;
-            }
-
-            if ( IntersLineCylinder( ptC - dLen1 * vtV1, Z_AX, CylFrame, dHei, dRad, true, true,
-                                     ptInt1, vtN1, ptInt2, vtN2)) {
-               AddIntervals( nGrid, i, j, ptInt1.z + dLen1 * vtV1.z, ptInt2.z + dLen1 * vtV1.z, - vtN1, - vtN2, CurrTool.GetToolNum()) ;
-            }
-
-            if ( IntersLineMyPolyhedron( ptC, Z_AX, PolyFrame, dLen1, 2 * ( dRad + dMyTol), dHei + 2 * dMyTol, 0,
-                                         ptInt1, vtN1, ptInt2, vtN2)) {
-               AddIntervals( nGrid, i, j, ptInt1.z, ptInt2.z, - vtN1, - vtN2, CurrTool.GetToolNum()) ;
-            }
-         }
-      }
-   }
    return true ;
 }
 
@@ -5158,7 +5352,9 @@ VolZmap::AddingTruncatedCone( int nGrid, const Point3d& ptS, const Point3d& ptE,
    Vector3d vtV3 = vtV1 ^ vtV2 ;
 
   // Sistema di riferimento intrinseco del movimento
-   Frame3d ConusFrame ;  ConusFrame.Set( ptV, vtV2, vtV3, vtV1) ;
+   Frame3d ConusFrame ;
+   if ( ! ConusFrame.Set( ptV, vtV1, vtV2, vtV3))
+      return false ;
 
   // Dimensioni lineari movimento
    double dLongLen = 0 ;
@@ -5183,7 +5379,8 @@ VolZmap::AddingTruncatedCone( int nGrid, const Point3d& ptS, const Point3d& ptE,
   // Sistema di riferimento poliedro
    Point3d ptO = ptV + vtV1 * dl + vtV2 * ( dCos * dMinRad) ;
    Frame3d PolyFrame ;
-   PolyFrame.Set( ptO, vtV1, vtV2, vtV3) ;
+   if ( ! PolyFrame.Set( ptO, vtV1, vtV2, vtV3))
+      return false ;
 
   // Versori piani nel riferimento poliedro ( riferiti al sistema di riferimento) :
   // Sx, Dx
@@ -5423,7 +5620,7 @@ VolZmap::AddingTruncatedCone( int nGrid, const Point3d& ptS, const Point3d& ptE,
                }
             }
          }
-        // Se il poliedro � attraversato, aggiungo
+        // Se il poliedro è attraversato, aggiungo
          if ( nIntNum == 2) {
 
            // Riporto le intersezioni nel sistema griglia
@@ -5459,7 +5656,8 @@ VolZmap::AddingSphere( int nGrid, const Point3d& ptS, const Point3d& ptE, double
       vtV /= dLengthPath ;
   // Riferimento per cilindro inviluppo della sfera lungo il movimento
    Frame3d CylFrame ;
-   CylFrame.Set( ptS, vtV) ;
+   if ( ! CylFrame.Set( ptS, vtV))
+      return false ;
 
    double dSqRad = dRad * dRad ;