EgtGeomKernel 1.8j4 :

- aggiunta classe Polygon3d (da EgtExchange) - razionalizzata classe Plane3d - corretta funzione IntersLineTria.
2017-10-16 07:56:04 +00:00
parent 31658bb165
commit 5bcd4bb67d
26 changed files with 969 additions and 105 deletions
@@ -31,12 +31,12 @@ CDBoxPlane( const BBox3d& Box, const Plane3d& Plane)
      return false ;

  // Compute the projection interval radius of b onto L(t) = ptCen + t * Plane.vtN
-   double dR = vtExt.x * abs( Plane.vtN.x) +
-               vtExt.y * abs( Plane.vtN.y) +
-               vtExt.z * abs( Plane.vtN.z) ;
+   double dR = vtExt.x * abs( Plane.GetVersN().x) +
+               vtExt.y * abs( Plane.GetVersN().y) +
+               vtExt.z * abs( Plane.GetVersN().z) ;

  // Compute distance of box center from plane
-   double dS = ( ptCen - ORIG) * Plane.vtN - Plane.dDist ;
+   double dS = ( ptCen - ORIG) * Plane.GetVersN() - Plane.GetDist() ;

  // Intersection occurs when distance dS falls within [-dR,+dR] interval
   return ( abs( dS) < dR + EPS_SMALL) ;
@@ -148,7 +148,6 @@ CDBoxTria( const BBox3d& b3Box, const Triangle3d& trTria)

  // Test separating axis corresponding to triangle face normal (category 2)
   Plane3d Plane ;
-   Plane.vtN = trTria.GetN() ;
-   Plane.dDist = Plane.vtN * ( trTria.GetP( 0) - ORIG) ;
+   Plane.Set( trTria.GetP( 0), trTria.GetN()) ;
   return CDBoxPlane( b3Box, Plane) ;
 }
@@ -789,14 +789,15 @@ CurveArc::IsFlat( Plane3d& plPlane, double dToler) const

  // assegno dati possibile piano
   bool bFlat = IsFlat( dToler) ;
+   Vector3d vtN ;
   if ( m_VtExtr.IsSmall()) {
-      plPlane.vtN = m_VtN ;
+      vtN = m_VtN ;
   }
   else {
-      plPlane.vtN = m_VtExtr ;
+      vtN = m_VtExtr ;
      bFlat = bFlat  &&  AreSameOrOppositeVectorApprox( m_VtExtr, m_VtN) ; 
   }
-   plPlane.dDist = ( m_PtCen - ORIG) * plPlane.vtN ;
+   plPlane.Set( m_PtCen, vtN) ;
  // ritorno conferma
   return bFlat ;
 } 
@@ -706,14 +706,14 @@ ProjectCurveOnPlane( const ICurve& crCrv, const Plane3d& plPlane)
   Plane3d plCrv ;
   if ( crCrv.IsFlat( plCrv)) {
     // se il piano della curva è parallelo a quello di proiezione
-      if ( AreSameOrOppositeVectorApprox( plCrv.vtN, plPlane.vtN)) {
+      if ( AreSameOrOppositeVectorApprox( plCrv.GetVersN(), plPlane.GetVersN())) {
        // copio la curva
         PtrOwner<ICurve> pCrv( crCrv.Clone()) ;
         if ( IsNull( pCrv))
            return nullptr ;
        // se non coincidenti, basta eseguire una traslazione
-         Point3d ptOC = ORIG + plCrv.dDist * plCrv.vtN ;
-         Point3d ptOP = ORIG + plPlane.dDist * plPlane.vtN ;
+         Point3d ptOC = ORIG + plCrv.GetDist() * plCrv.GetVersN() ;
+         Point3d ptOP = ORIG + plPlane.GetDist() * plPlane.GetVersN() ;
         if ( ! AreSamePointApprox( ptOC, ptOP))
            pCrv->Translate( ptOP - ptOC) ;
        // restituisco la nuova curva
@@ -727,7 +727,7 @@ ProjectCurveOnPlane( const ICurve& crCrv, const Plane3d& plPlane)
      return nullptr ;
  // calcolo un riferimento con piano XY coincidente con il piano di proiezione
   Frame3d frRef ;
-   if ( ! frRef.Set( ORIG + plPlane.dDist * plPlane.vtN, plPlane.vtN))
+   if ( ! frRef.Set( ORIG + plPlane.GetDist() * plPlane.GetVersN(), plPlane.GetVersN()))
      return nullptr ;
  // eseguo scalatura con fattori X e Y unitari e fattore Z nullo
   if ( ! pCrv->Scale( frRef, 1, 1, 0))
@@ -585,34 +585,26 @@ CurveBezier::IsFlat( Plane3d& plPlane, double dToler) const
   // se punto
   switch ( nRank) {
   case 0 :      // punto
-      if ( bFlat) {
-         plPlane.vtN = ( m_VtExtr.IsSmall() ? Z_AX : m_VtExtr) ;
-         plPlane.dDist = ( ptCen - ORIG) * plPlane.vtN ;
-      }
-      else {
-         plPlane.vtN = ( m_VtExtr.IsSmall() ? V_NULL : m_VtExtr) ;
-         plPlane.dDist = 0 ;
-      }
+      if ( bFlat)
+         plPlane.Set( ptCen, ( m_VtExtr.IsSmall() ? Z_AX : m_VtExtr)) ;
+      else
+         plPlane.Set( 0, ( m_VtExtr.IsSmall() ? V_NULL : m_VtExtr)) ;
      break ;
   case 1 :      // linea
-      if ( m_VtExtr.IsSmall()) {
-         plPlane.vtN = FromUprightOrtho( vtDir) ;
-      }
+      if ( m_VtExtr.IsSmall())
+         plPlane.Set( ptCen, FromUprightOrtho( vtDir)) ;
      else {
-         plPlane.vtN = m_VtExtr ;
+         plPlane.Set( ptCen, m_VtExtr) ;
         bFlat = bFlat && ( AreOrthoApprox( m_VtExtr, vtDir)) ;
      }
-      plPlane.dDist = ( ptCen - ORIG) * plPlane.vtN ;
      break ;
   default :     // piana o 3d
-      if ( m_VtExtr.IsSmall()) {
-         plPlane.vtN = vtDir ;
-      }
+      if ( m_VtExtr.IsSmall())
+         plPlane.Set( ptCen, vtDir) ;
      else {
-         plPlane.vtN = m_VtExtr ;
+         plPlane.Set( ptCen, m_VtExtr) ;
         bFlat = bFlat && ( AreSameOrOppositeVectorApprox( m_VtExtr, vtDir)) ;
      }
-      plPlane.dDist = ( ptCen - ORIG) * plPlane.vtN ;
      break ;
   }

@@ -860,34 +860,26 @@ CurveComposite::IsFlat( Plane3d& plPlane, double dToler) const
   // se punto
   switch ( nRank) {
   case 0 :      // punto
-      if ( bFlat) {
-         plPlane.vtN = ( m_VtExtr.IsSmall() ? Z_AX : m_VtExtr) ;
-         plPlane.dDist = ( ptCen - ORIG) * plPlane.vtN ;
-      }
-      else {
-         plPlane.vtN = ( m_VtExtr.IsSmall() ? V_NULL : m_VtExtr) ;
-         plPlane.dDist = 0 ;
-      }
+      if ( bFlat)
+         plPlane.Set( ptCen, ( m_VtExtr.IsSmall() ? Z_AX : m_VtExtr)) ;
+      else
+         plPlane.Set( 0, ( m_VtExtr.IsSmall() ? V_NULL : m_VtExtr)) ;
      break ;
   case 1 :      // linea
-      if ( m_VtExtr.IsSmall()) {
-         plPlane.vtN = FromUprightOrtho( vtDir) ;
-      }
+      if ( m_VtExtr.IsSmall())
+         plPlane.Set( ptCen, FromUprightOrtho( vtDir)) ;
      else {
-         plPlane.vtN = m_VtExtr ;
+         plPlane.Set( ptCen, m_VtExtr) ;
         bFlat = bFlat && ( AreOrthoApprox( m_VtExtr, vtDir)) ;
      }
-      plPlane.dDist = ( ptCen - ORIG) * plPlane.vtN ;
      break ;
   default :     // piana o 3d
-      if ( m_VtExtr.IsSmall()) {
-         plPlane.vtN = vtDir ;
-      }
+      if ( m_VtExtr.IsSmall())
+         plPlane.Set( ptCen, vtDir) ;
      else {
-         plPlane.vtN = m_VtExtr ;
+         plPlane.Set( ptCen, m_VtExtr) ;
         bFlat = bFlat && ( AreSameOrOppositeVectorApprox( m_VtExtr, vtDir)) ;
      }
-      plPlane.dDist = ( ptCen - ORIG) * plPlane.vtN ;
      break ;
   }

@@ -271,15 +271,16 @@ CurveLine::IsFlat( Plane3d& plPlane, double dToler) const

  // assegno dati piano
   bool bFlat ;
+   Vector3d vtN ;
   if ( m_VtExtr.IsSmall()) {
-      plPlane.vtN = FromUprightOrtho( m_PtEnd - m_PtStart) ;
+      vtN = FromUprightOrtho( m_PtEnd - m_PtStart) ;
      bFlat = true ;
   }
   else {
-      plPlane.vtN = m_VtExtr ;
+      vtN = m_VtExtr ;
      bFlat = ( AreOrthoApprox( m_VtExtr, ( m_PtEnd - m_PtStart))) ;
   }
-   plPlane.dDist = ( 0.5 * ( m_PtStart + m_PtEnd) - ORIG) * plPlane.vtN ;
+   plPlane.Set( 0.5 * ( m_PtStart + m_PtEnd), vtN) ;

  // ritorno conferma
   return bFlat ;
@@ -267,6 +267,7 @@ copy $(TargetPath) \EgtProg\Dll64</Command>
    <ClCompile Include="CurveByApprox.cpp" />
    <ClCompile Include="CurveByInterp.cpp" />
    <ClCompile Include="CurveCompositeOffset.cpp" />
+    <ClCompile Include="IntersPlanePlane.cpp" />
    <ClCompile Include="OffsetCurve.cpp" />
    <ClCompile Include="DistPointArc.cpp" />
    <ClCompile Include="DistPointCrvAux.cpp" />
@@ -324,6 +325,7 @@ copy $(TargetPath) \EgtProg\Dll64</Command>
    <ClCompile Include="NgeReader.cpp" />
    <ClCompile Include="FontOs.cpp" />
    <ClCompile Include="OffsetCurveOnX.cpp" />
+    <ClCompile Include="Polygon3d.cpp" />
    <ClCompile Include="RemoveCurveOverlaps.cpp" />
    <ClCompile Include="RemoveCurveSpikes.cpp" />
    <ClCompile Include="SelfIntersCurve.cpp" />
@@ -405,6 +407,7 @@ copy $(TargetPath) \EgtProg\Dll64</Command>
    <ClInclude Include="..\Include\EGkIntersLinePlane.h" />
    <ClInclude Include="..\Include\EGkIntersLineSurfTm.h" />
    <ClInclude Include="..\Include\EGkIntersLineTria.h" />
+    <ClInclude Include="..\Include\EGkIntersPlanePlane.h" />
    <ClInclude Include="..\Include\EGkLinePerpTwoCurves.h" />
    <ClInclude Include="..\Include\EGkLinePntMinDistCurve.h" />
    <ClInclude Include="..\Include\EGkLinePntPerpCurve.h" />
@@ -357,6 +357,12 @@
    <ClCompile Include="RemoveCurveSpikes.cpp">
      <Filter>File di origine\GeoInters</Filter>
    </ClCompile>
+    <ClCompile Include="Polygon3d.cpp">
+      <Filter>File di origine\Base</Filter>
+    </ClCompile>
+    <ClCompile Include="IntersPlanePlane.cpp">
+      <Filter>File di origine\GeoInters</Filter>
+    </ClCompile>
  </ItemGroup>
  <ItemGroup>
    <ClInclude Include="stdafx.h">
@@ -782,6 +788,9 @@
    <ClInclude Include="RemoveCurveSpikes.h">
      <Filter>File di intestazione</Filter>
    </ClInclude>
+    <ClInclude Include="..\Include\EGkIntersPlanePlane.h">
+      <Filter>File di intestazione\Include</Filter>
+    </ClInclude>
  </ItemGroup>
  <ItemGroup>
    <ResourceCompile Include="EgtGeomKernel.rc">
@@ -17,7 +17,7 @@

 //----------------------------------------------------------------------------
 int
-IntersLinePlane( const Point3d& ptL1, const Point3d& ptL2, const Plane3d& plPlane, Point3d& ptInt)
+IntersLinePlane( const Point3d& ptL1, const Point3d& ptL2, const Plane3d& plPlane, Point3d& ptInt, bool bFinite)
 {
   Vector3d vtL = ptL2 - ptL1 ;
   double dLen = vtL.Len() ;
@@ -25,20 +25,20 @@ IntersLinePlane( const Point3d& ptL1, const Point3d& ptL2, const Plane3d& plPlan
      vtL /= dLen ;
   else
      vtL = V_NULL ;
-   return IntersLinePlane( ptL1, vtL, dLen, plPlane, ptInt) ;
+   return IntersLinePlane( ptL1, vtL, dLen, plPlane, ptInt, bFinite) ;
 }

 //----------------------------------------------------------------------------
 int
-IntersLinePlane( const Point3d& ptL, const Vector3d& vtL, double dLen, const Plane3d& plPlane, Point3d& ptInt)
+IntersLinePlane( const Point3d& ptL, const Vector3d& vtL, double dLen, const Plane3d& plPlane, Point3d& ptInt, bool bFinite)
 {
   Vector3d vtDir = vtL ;
   if ( dLen > EPS_SMALL  &&  ! vtDir.Normalize( EPS_ZERO))
      return ILPT_NO ;
  // coseno dell'angolo tra direzione linea e normale al piano
-   double dCosAng = vtDir * plPlane.vtN ;
+   double dCosAng = vtDir * plPlane.GetVersN() ;
  // posizione del punto rispetto al piano
-   double dPosL = ( ptL - ORIG) * plPlane.vtN - plPlane.dDist ;
+   double dPosL = ( ptL - ORIG) * plPlane.GetVersN() - plPlane.GetDist() ;
  // verifico se linea parallela al piano (ovvero ortogonale alla normale)
   if ( fabs( dCosAng) < COS_ORTO_ANG_ZERO) {
     // se il punto giace nel piano, vi giace anche la linea
@@ -51,6 +51,9 @@ IntersLinePlane( const Point3d& ptL, const Vector3d& vtL, double dLen, const Pla
  // determino il punto di intersezione
   double dDistI = - dPosL / dCosAng ;
   ptInt = ptL + vtDir * dDistI ;
+  // se linea infinita, ho trovato l'intersezione
+   if ( ! bFinite)
+      return ILPT_YES ;
  // intersezione sul primo estremo
   if ( fabs( dDistI) < EPS_SMALL)
      return ILPT_START ;
@@ -46,8 +46,7 @@ IntersLineTria( const Point3d& ptL, const Vector3d& vtL, double dLen, const Tria
 {
  // determino il piano del triangolo
   Plane3d plPlane ;
-   plPlane.vtN = trTria.GetN() ;
-   plPlane.dDist = ( trTria.GetP(0) - ORIG) * plPlane.vtN ;
+   plPlane.Set( trTria.GetP( 0), trTria.GetN()) ;
  // calcolo l'intersezione tra il segmento di linea e il piano del triangolo
   int nRes = IntersLinePlane( ptL, vtL, dLen, plPlane, ptInt) ;
  // se non c'è intersezione
@@ -83,8 +82,17 @@ IntersCoplanarLineTria( const Point3d& ptL, const Vector3d& vtL, double dLen, co
   ptInt2 = ptL + vtL * dLen ;
   int nPTT = PointInTria( ptInt, trTria) ;
   int nPTT2 = PointInTria( ptInt2, trTria) ;
-   if ( nPTT != PTT_OUT  &&  nPTT2 != PTT_OUT)
+   if ( ( nPTT == PTT_IN  &&  nPTT2 != PTT_OUT)  ||  ( nPTT != PTT_OUT  &&  nPTT2 == PTT_IN))
      return ILTT_SEGM ;
+   if ( nPTT == PTT_VERT  &&  nPTT2 == PTT_VERT)
+      return ILTT_SEGM_ON_EDGE ;
+   if ( ( nPTT == PTT_VERT  ||  nPTT == PTT_EDGE)  &&  ( nPTT2 == PTT_VERT  ||  nPTT2 == PTT_EDGE)) {
+      Point3d ptMed = Media( ptInt, ptInt2, 0.5) ;
+      if ( PointInTria( ptMed, trTria) == PTT_IN)
+         return ILTT_SEGM ;
+      else
+         return ILTT_SEGM_ON_EDGE ;
+   }
  // se segmento nullo, basta verificare la posizione del punto iniziale
   if ( dLen < EPS_SMALL  ||  vtL.IsSmall()) {
      switch ( nPTT) {
@@ -140,13 +148,18 @@ IntersCoplanarLineTria( const Point3d& ptL, const Vector3d& vtL, double dLen, co
      else {
        // se l'inizio è nel triangolo
         if ( nPTT == PTT_IN) {
-            ptInt2 = ptInt ;
            ptInt = ptL ; 
+            ptInt2 = ptInt ;
+            return ILTT_SEGM ;
         }
-        // altrimenti è la fine nel triangolo
-         else
+        // se altrimenti è la fine nel triangolo
+         else if ( nPTT2 == PTT_IN) {
            ptInt2 = ptL + vtL * dLen ;
-         return ILTT_SEGM ;
+            return ILTT_SEGM ;
+         }
+        // altrimenti entrambi esterni, è un vertice
+         else
+            return ILTT_VERT ;
      }
   }
  // se 2 intersezioni
@@ -0,0 +1,46 @@
+//----------------------------------------------------------------------------
+//  EgalTech 2017-2017
+//----------------------------------------------------------------------------
+// File : IntersPlanePlane.cpp       Data : 15.10.17         Versione : 1.8j3
+// Contenuto : Implementazione della intersezione piano/piano.
+//
+//
+//
+// Modifiche : 15.12.17 DS Creazione modulo.
+//
+//
+//----------------------------------------------------------------------------
+
+//--------------------------- Include ----------------------------------------
+#include "stdafx.h"
+#include "/EgtDev/Include/EGkIntersLinePlane.h"
+#include "/EgtDev/Include/EGkIntersPlanePlane.h"
+
+//----------------------------------------------------------------------------
+int
+IntersPlanePlane( const Plane3d& plPlane1, const Plane3d& plPlane2, Point3d& ptInt, Vector3d& vtDir)
+{
+  // Direzione linea di intersezione
+   vtDir = plPlane1.GetVersN() ^ plPlane2.GetVersN() ;
+  // Verifico se piani praticamente paralleli
+   double dDenom = vtDir * vtDir ;
+   if ( dDenom < EPS_ZERO * EPS_ZERO)
+      return ( abs( plPlane1.GetDist() - plPlane2.GetDist()) < EPS_SMALL ? IPPT_OVERLAPS : IPPT_NO) ;
+  // Calcolo un punto sulla retta di intersezione
+   ptInt = ORIG + ( ( plPlane1.GetDist() * plPlane2.GetVersN() - plPlane2.GetDist() * plPlane1.GetVersN()) ^ vtDir) / dDenom ;
+  // Normalizzo la direzione
+   vtDir.Normalize() ;
+   return IPPT_YES ;
+}
+
+//----------------------------------------------------------------------------
+int
+Inters3Planes( const Plane3d& plPlane1, const Plane3d& plPlane2, const Plane3d& plPlane3, Point3d& ptInt)
+{
+  // linea di intersezione tra i primi due piani
+   Point3d ptL ; Vector3d vtL ;
+   if ( IntersPlanePlane( plPlane1, plPlane2, ptL, vtL) != IPPT_YES)
+      return IPPT_NO ;
+  // intersezione della linea con il terzo piano
+   return ( IntersLinePlane( ptL, vtL, 1, plPlane3, ptInt) == ILPT_YES ? IPPT_YES : IPPT_NO) ;
+}
@@ -63,11 +63,11 @@ OffsetCurve::Make( const ICurve* pCrv, double dDist, int nType)
  // se esiste estrusione, verifico sia perpendicolare al piano della curva
   Vector3d vtExtr ;
   if ( pCrv->GetExtrusion( vtExtr)  &&  ! vtExtr.IsSmall()) {
-      if ( ! AreSameOrOppositeVectorApprox( plPlane.vtN, vtExtr))
+      if ( ! AreSameOrOppositeVectorApprox( plPlane.GetVersN(), vtExtr))
         return false ;
   }
   else
-      vtExtr = plPlane.vtN ;
+      vtExtr = plPlane.GetVersN() ;
  // recupero spessore
   double dThick ;
   pCrv->GetThickness( dThick) ;
@@ -56,11 +56,11 @@ OffsetCurveOnX::Make( const ICurve* pCrv, double dDist)
  // se esiste estrusione, verifico sia perpendicolare al piano della curva
   Vector3d vtExtr ;
   if ( pCrv->GetExtrusion( vtExtr)  &&  ! vtExtr.IsSmall()) {
-      if ( ! AreSameOrOppositeVectorApprox( plPlane.vtN, vtExtr))
+      if ( ! AreSameOrOppositeVectorApprox( plPlane.GetVersN(), vtExtr))
         return false ;
   }
   else
-      vtExtr = plPlane.vtN ;
+      vtExtr = plPlane.GetVersN() ;

  // determino se necessario cambiare riferimento ( dal vettore estrusione)
   bool bNeedRef = ( ! vtExtr.IsZplus()) ;
@@ -532,11 +532,10 @@ PolyLine::IsFlat( int& nRank, Point3d& ptCen, Vector3d& vtDir, double dToler) co
      vtDir.Invert() ;
   // Plane calculation
   Plane3d plPlane ;
-   plPlane.vtN = vtDir ;
-   plPlane.dDist = ( ptCen - ORIG) * plPlane.vtN ;
+   plPlane.Set( ptCen, vtDir) ;
   // Test each vertex to see if it is farther from plane than allowed max distance
   for ( bool bFound = GetFirstPoint( ptP) ; bFound ; bFound = GetNextPoint( ptP)) {
-      double dDist = ( ( ptP - ORIG) * plPlane.vtN) - plPlane.dDist ;
+      double dDist = ( ( ptP - ORIG) * plPlane.GetVersN()) - plPlane.GetDist() ;
      if ( fabs( dDist) > dToler)
         return false ;
   }
@@ -559,7 +558,7 @@ PolyLine::IsClosedAndFlat( Plane3d& plPlane, double& dArea, double dToler) const
      return false ;
  // Test each vertex to see if it is farther from plane than allowed max distance
   for ( bool bFound = GetFirstPoint( ptP) ; bFound ; bFound = GetNextPoint( ptP)) {
-      double dDist = ( ( ptP - ORIG) * plPlane.vtN) - plPlane.dDist ;
+      double dDist = ( ( ptP - ORIG) * plPlane.GetVersN()) - plPlane.GetDist() ;
      if ( fabs( dDist) > dToler)
         return false ;
   }
@@ -0,0 +1,306 @@
+//----------------------------------------------------------------------------
+//  EgalTech 2015-2017
+//----------------------------------------------------------------------------
+// File : Polygon3d.cpp         Data : 15.10.17         Versione : 1.8j3
+// Contenuto : Implementazione della classe Polygon3d.
+//
+//
+//
+// Modifiche : 30.08.15 DS Creazione modulo.
+//
+//
+//----------------------------------------------------------------------------
+
+//--------------------------- Include ----------------------------------------
+#include "stdafx.h"
+#include "/EgtDev/Include/EgkPolygon3d.h"
+
+//----------------------------------------------------------------------------
+bool
+Polygon3d::FromRectangle( double dDimX, double dDimY)
+{
+  // verifico i dati
+   if ( dDimX < EPS_SMALL  ||  dDimY < EPS_SMALL)
+      return false ;
+  // assegno il piano
+   m_Plane.Set( 0, Z_AX) ;
+  // assegno i vertici del poligono
+   m_vVert.reserve( 4) ;
+   m_vVert.emplace_back( + 0.5 * dDimX, + 0.5 * dDimY, 0) ;
+   m_vVert.emplace_back( - 0.5 * dDimX, + 0.5 * dDimY, 0) ;
+   m_vVert.emplace_back( - 0.5 * dDimX, - 0.5 * dDimY, 0) ;
+   m_vVert.emplace_back( + 0.5 * dDimX, - 0.5 * dDimY, 0) ;
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+bool
+Polygon3d::FromPolyLine( const PolyLine& PL)
+{
+  // annullo il piano del poligono
+   m_Plane.Reset() ;
+  // verifico sia chiusa e piana
+   Plane3d plPlane ;
+   double dArea ;
+   if ( ! PL.IsClosedAndFlat( plPlane, dArea))
+      return false ;
+  // assegno il piano
+   m_Plane = plPlane ;
+  // assegno i punti
+   Point3d ptP ;
+   bool bPoint = PL.GetFirstPoint( ptP, true) ;
+   // inserisco i punti
+   while ( bPoint) {
+      m_vVert.emplace_back( ptP) ;
+      bPoint = PL.GetNextPoint( ptP, true) ;
+   }
+
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+bool
+Polygon3d::FromPlaneTrimmedWithBox( const Plane3d& plPlane, const Point3d& ptMin, const Point3d& ptMax)
+{
+  // il piano e il box devono essere validi
+   if ( plPlane.GetVersN().IsSmall()  ||  AreSamePointApprox( ptMin, ptMax))
+      return false ;
+  // assegno il piano
+   m_Plane = plPlane ;
+  // centro del box proiettato nel piano
+   Point3d ptCen = Media( ptMin, ptMax, 0.5) ;
+   ptCen -= (( ptCen - ORIG) * plPlane.GetVersN() - plPlane.GetDist()) * plPlane.GetVersN() ;
+  // raggio del box
+   double dDiam = Dist( ptMin, ptMax) ;
+  // poligono nel suo piano
+   FromRectangle( dDiam, dDiam) ;
+  // riferimento del piano
+   Frame3d frRef ;
+   frRef.Set( ptCen, plPlane.GetVersN(), FromUprightOrtho( plPlane.GetVersN())) ;
+   ToGlob( frRef) ;
+  // lo trimmo con le facce del box
+   Plane3d plFace ;
+   plFace.Set( ptMin.y, -Y_AX) ;       // faccia 1 (Norm = Y-) 
+   Trim( plFace, true, true) ;
+   plFace.Set( ptMin.z, -Z_AX) ;       // faccia 2 (Norm = Z-) 
+   Trim( plFace, true, true) ;
+   plFace.Set( ptMax.y, Y_AX) ;        // faccia 3 (Norm = Y+) 
+   Trim( plFace, true, true) ;
+   plFace.Set( ptMax.z, Z_AX) ;        // faccia 4 (Norm = Z+) 
+   Trim( plFace, true, true) ;
+   plFace.Set( ptMin.x, -X_AX) ;       // faccia 5 (Norm = X-) 
+   Trim( plFace, true, true) ;
+   plFace.Set( ptMax.x, X_AX) ;        // faccia 6 (Norm = X+) 
+   Trim( plFace, true, true) ;
+
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+bool
+Polygon3d::FromPlaneTrimmedWithBox( const Point3d& ptOn, const Vector3d& vtN, const Point3d& ptMin, const Point3d& ptMax)
+{
+   Plane3d plPlane ;
+   if ( plPlane.Set( ptOn, vtN))
+      return false ;
+   return FromPlaneTrimmedWithBox( plPlane, ptMin, ptMax) ;
+}
+
+//----------------------------------------------------------------------------
+bool
+Polygon3d::Trim( const Plane3d& plPlane, bool bInVsOut, bool bOnEq)
+{
+  // calcolo le distanze dei punti del poligono dal piano
+   DBLVECTOR vDist ;
+   vDist.reserve( m_vVert.capacity()) ;
+   for ( auto& ptP : m_vVert)
+      vDist.push_back( ( ptP - ORIG) * plPlane.GetVersN() - plPlane.GetDist()) ;
+
+  // verifico se il poligono giace nel piano
+   bool bOnPlane = true ;
+   for ( auto& dDist : vDist) {
+      if ( fabs( dDist) > EPS_SMALL)
+         bOnPlane = false ;
+   }
+   if ( bOnPlane) {
+     // se non ammesso equiverso oppure controverso, lo annullo
+      if ( ! bOnEq  ||  m_Plane.GetVersN() * plPlane.GetVersN() < 0)
+         m_vVert.clear() ;
+      return true ;
+   }
+
+  // determino le intersezioni dei lati con il piano
+   for ( size_t i = 0 ; i < m_vVert.size() ; ++ i) {
+     // indice del punto finale del lato
+      size_t j = ( ( i + 1 < m_vVert.size()) ? i + 1 : 0) ;
+     // se il lato attraversa il piano, inserisco il punto di intersezione nel poligono
+      if ( ( vDist[i] > EPS_SMALL  &&  vDist[j] < - EPS_SMALL)  ||
+           ( vDist[i] < - EPS_SMALL  &&  vDist[j] > EPS_SMALL)) {
+         double dCoeff = fabs( vDist[i]) / ( fabs( vDist[i]) + fabs( vDist[j])) ;
+         Point3d ptInt = Media( m_vVert[i], m_vVert[j], dCoeff) ;
+         m_vVert.insert( m_vVert.begin() + i + 1, ptInt) ;
+         vDist.insert( vDist.begin() + i + 1, 0.0) ;
+         ++ i ;
+      }
+   }
+
+  // elimino i punti che non rispettano la posizione rispetto al piano
+   for ( size_t i = 0 ; i < m_vVert.size() ; ++ i) {
+      if ( ( bInVsOut  &&  vDist[i] > EPS_SMALL)  ||
+           ( ! bInVsOut  &&  vDist[i] < - EPS_SMALL)) {
+         m_vVert.erase( m_vVert.begin() + i) ;
+         vDist.erase( vDist.begin() + i) ;
+         -- i ;
+      }
+   }
+
+  // se i punti rimasti giacciono tutti sul piano (quindi su una linea), annullo il poligono
+   bool bIsLine = true ;
+   for ( auto& dDist : vDist) {
+      if ( fabs( dDist) > EPS_SMALL)
+         bIsLine = false ;
+   }
+   if ( bIsLine)
+      m_vVert.clear() ;
+
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+bool
+Polygon3d::Trim( const Polygon3d& plyOther, bool bInVsOut, bool bOnEq)
+{
+   if ( GetSideCount() == 0  ||  plyOther.GetSideCount() == 0)
+      return true ;
+   return Trim( plyOther.m_Plane, bInVsOut, bOnEq) ;
+}
+
+//----------------------------------------------------------------------------
+PolyLine
+Polygon3d::GetPolyLine( void) const
+{
+   PolyLine PL ;
+  // se il piano non è definito, errore
+   if ( m_Plane.GetVersN().IsSmall())
+      return PL ;
+  // converto il poligono in PolyLine
+   for ( size_t i = 0 ; i < m_vVert.size() ; ++ i)
+      PL.AddUPoint( int( i), m_vVert[i]) ;
+   PL.Close() ;
+   return PL ;
+}
+
+//----------------------------------------------------------------------------
+void
+Polygon3d::Translate( const Vector3d& vtMove)
+{
+  // traslo il piano
+   m_Plane.Translate( vtMove) ;
+  // traslo i punti
+   for ( auto& ptP : m_vVert)
+      ptP.Translate( vtMove) ;
+}
+
+//----------------------------------------------------------------------------
+bool
+Polygon3d::Rotate( const Point3d& ptAx, const Vector3d& vtAx, double dCosAng, double dSinAng)
+{
+  // ruoto il piano
+   if ( ! m_Plane.Rotate( ptAx, vtAx, dCosAng, dSinAng))
+      return false ;
+  // ruoto i punti
+   for ( auto& ptP : m_vVert)
+      ptP.Rotate( ptAx, vtAx, dCosAng, dSinAng) ;
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+bool
+Polygon3d::Scale( const Frame3d& frRef, double dCoeffX, double dCoeffY, double dCoeffZ)
+{
+  // verifico validità della scalatura
+   if ( fabs( dCoeffX) < EPS_ZERO  &&  fabs( dCoeffY) < EPS_ZERO  &&  fabs( dCoeffZ) < EPS_ZERO)
+      return false ;
+  // scalo il piano
+   if ( ! m_Plane.Scale( frRef, dCoeffX, dCoeffY, dCoeffZ))
+      return false ;
+  // scalo i punti
+   for ( auto& ptP : m_vVert)
+      ptP.Scale( frRef, dCoeffX, dCoeffY, dCoeffZ) ;
+  // determino se contiene anche un mirror (numero dispari di coefficienti negativi)
+   bool bMirror = ( dCoeffX < 0) ;
+   bMirror = ( bMirror ? ( dCoeffY > 0) : ( dCoeffY < 0)) ; 
+   bMirror = ( bMirror ? ( dCoeffZ > 0) : ( dCoeffZ < 0)) ; 
+   if ( bMirror)
+      reverse( m_vVert.begin(), m_vVert.end()) ;
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+bool
+Polygon3d::Mirror( const Point3d& ptOn, const Vector3d& vtNorm)
+{
+  // specchio il piano
+   if ( ! m_Plane.Mirror( ptOn, vtNorm))
+      return false ;
+  // specchio i punti e ne inverto l'ordine
+   for ( auto& ptP : m_vVert)
+      ptP.Mirror( ptOn, vtNorm) ;
+   reverse( m_vVert.begin(), m_vVert.end()) ;
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+bool
+Polygon3d::Shear( const Point3d& ptOn, const Vector3d& vtNorm, const Vector3d& vtDir, double dCoeff)
+{
+  // stiro il piano
+   if ( ! m_Plane.Shear( ptOn, vtNorm, vtDir, dCoeff))
+      return false ;
+  // stiro i punti
+   for ( auto& ptP : m_vVert)
+      ptP.Shear( ptOn, vtNorm, vtDir, dCoeff) ;
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+bool
+Polygon3d::ToGlob( const Frame3d& frRef)
+{
+  // trasformo il piano
+   if ( ! m_Plane.ToGlob( frRef))
+      return false ;
+  // trasformo i punti
+   for ( auto& ptP : m_vVert)
+      ptP.ToGlob( frRef) ;
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+bool
+Polygon3d::ToLoc( const Frame3d& frRef)
+{
+  // trasformo il piano
+   if ( ! m_Plane.ToLoc( frRef))
+      return false ;
+  // trasformo i punti
+   for ( auto& ptP : m_vVert)
+      ptP.ToLoc( frRef) ;
+   return true ;
+}
+
+//----------------------------------------------------------------------------
+bool
+Polygon3d::LocToLoc( const Frame3d& frOri, const Frame3d& frDest)
+{
+  // se i due riferimenti coincidono, non devo fare alcunché
+   if ( AreSameFrame( frOri, frDest))
+      return true ;
+  // trasformo il piano
+   if ( ! m_Plane.LocToLoc( frOri, frDest))
+      return false ;
+  // trasformo i punti
+   for ( auto& ptP : m_vVert)
+      ptP.LocToLoc( frOri, frDest) ;
+   return true ;
+}
@@ -133,9 +133,7 @@ PolygonPlane::GetPlane( Plane3d& plPlane)
   if ( ! Finalize())
      return false ;
  // assegno i dati al piano
-   plPlane.vtN = m_vtN ;
-   plPlane.dDist = ( ( m_ptMid - ORIG) * plPlane.vtN) ;
-   return true ;
+   return plPlane.Set( m_ptMid, m_vtN) ;
 }

 //----------------------------------------------------------------------------
@@ -284,12 +284,12 @@ SurfFlatRegionByContours::Prepare( void)
         return false ;
     // imposto la normale del primo contorno come riferimento
      if ( i == 0)
-         vtN0 = plPlane.vtN ;
+         vtN0 = plPlane.GetVersN() ;
     // verifico che le normali siano molto vicine
-      if ( ! AreSameOrOppositeVectorApprox( plPlane.vtN, vtN0))
+      if ( ! AreSameOrOppositeVectorApprox( plPlane.GetVersN(), vtN0))
         return false ;
     // assegno il segno all'area secondo il verso della normale
-      if ( ( plPlane.vtN * vtN0) > 0)
+      if ( ( plPlane.GetVersN() * vtN0) > 0)
         m_vArea.emplace_back( i, dArea) ;
      else
         m_vArea.emplace_back( i, - dArea) ;
@@ -87,7 +87,7 @@ GetSurfTriMeshByExtrusion( const ICurve* pCurve, const Vector3d& vtExtr,
      double dArea ;
      if ( PL.IsClosedAndFlat( plPlane, dArea, 50 * EPS_SMALL)) {
        // componente dell'estrusione perpendicolare al piano della curva
-         double dOrthoExtr = plPlane.vtN * vtExtr ;
+         double dOrthoExtr = plPlane.GetVersN() * vtExtr ;
         if ( ( fabs( dOrthoExtr) > EPS_SMALL)) {
            bDoCapEnds = true ;
           // se negativa, inverto il senso del contorno
@@ -354,12 +354,12 @@ CalcRegionPolyLines( const CICURVEPVECTOR& vpCurve, double dLinTol,
         return false ;
     // imposto la normale del primo contorno come riferimento
      if ( i == 0)
-         vtN0 = plPlane.vtN ;
+         vtN0 = plPlane.GetVersN() ;
     // verifico che le normali siano molto vicine
-      if ( ! AreSameOrOppositeVectorApprox( plPlane.vtN, vtN0))
+      if ( ! AreSameOrOppositeVectorApprox( plPlane.GetVersN(), vtN0))
         return false ;
     // assegno il segno all'area secondo il verso della normale
-      if ( ( plPlane.vtN * vtN0) > 0)
+      if ( ( plPlane.GetVersN() * vtN0) > 0)
         vArea.emplace_back( i, dArea) ;
      else
         vArea.emplace_back( i, - dArea) ;
@@ -89,7 +89,7 @@ SurfFlatRegion::AddExtLoop( ICurve* pCrv)
   Plane3d plPlane ;
   if ( ! pCrv->GetArea( plPlane, dArea))
      return false ;
-   pCrv->SetExtrusion( plPlane.vtN) ;
+   pCrv->SetExtrusion( plPlane.GetVersN()) ;
   pCrv->SetThickness( 0) ;
  // rimuovo eventuali sovrapposizioni (calcolate nel suo piano)
   ICURVEPLIST CrvLst ;
@@ -132,7 +132,7 @@ SurfFlatRegion::AddSimpleExtLoop( ICurve* pCrv)
  // se sto costruendo il primo chunk
   if ( m_vExtInd.size() == 0) {
     // assegno il riferimento intrinseco
-      if ( ! m_frF.Set( ORIG + plPlane.vtN * plPlane.dDist, plPlane.vtN))
+      if ( ! m_frF.Set( ORIG + plPlane.GetVersN() * plPlane.GetDist(), plPlane.GetVersN()))
         return false ;
     // sistemo il senso di rotazione (deve essere CCW -> area > 0)
      if ( dArea < 0)
@@ -142,11 +142,11 @@ SurfFlatRegion::AddSimpleExtLoop( ICurve* pCrv)
   else {
     // verifico che il piano della curva coincida con quello XY intrinseco
      plPlane.ToLoc( m_frF) ;
-      if ( ! ( plPlane.vtN.IsZplus()  ||  plPlane.vtN.IsZminus())  ||  fabs( plPlane.dDist) > EPS_SMALL)
+      if ( ! ( plPlane.GetVersN().IsZplus()  ||  plPlane.GetVersN().IsZminus())  ||  fabs( plPlane.GetDist()) > EPS_SMALL)
         return false ;
     // sistemo il senso di rotazione (deve essere CCW -> area > 0)
-      if ( ( plPlane.vtN.IsZplus()  &&  dArea < 0)  ||
-           ( plPlane.vtN.IsZminus()  &&  dArea > 0))
+      if ( ( plPlane.GetVersN().IsZplus()  &&  dArea < 0)  ||
+           ( plPlane.GetVersN().IsZminus()  &&  dArea > 0))
         pMyCrv->Invert() ;
   }
  // porto la curva nel riferimento intrinseco
@@ -235,7 +235,7 @@ SurfFlatRegion::AddIntLoop( ICurve* pCrv)
   Plane3d plPlane ;
   if ( ! pCrv->GetArea( plPlane, dArea))
      return false ;
-   pCrv->SetExtrusion( plPlane.vtN) ;
+   pCrv->SetExtrusion( plPlane.GetVersN()) ;
   pCrv->SetThickness( 0) ;
  // rimuovo eventuali sovrapposizioni (calcolate nel suo piano)
   ICURVEPLIST CrvLst ;
@@ -284,10 +284,10 @@ SurfFlatRegion::AddSimpleIntLoop( ICurve* pCrv)
   Plane3d plPlane ;
   if ( ! pMyCrv->GetArea( plPlane, dArea))
      return false ;
-   if ( ! ( plPlane.vtN.IsZplus()  ||  plPlane.vtN.IsZminus())  ||  fabs( plPlane.dDist) > EPS_SMALL)
+   if ( ! ( plPlane.GetVersN().IsZplus()  ||  plPlane.GetVersN().IsZminus())  ||  fabs( plPlane.GetDist()) > EPS_SMALL)
      return false ;
  // sistemo il senso di rotazione (deve essere CW -> se N==Z+ area < 0, se N==Z- area>0)
-   if ( ( plPlane.vtN.IsZplus()  &&  dArea > 0)  ||  ( plPlane.vtN.IsZminus()  &&  dArea < 0))
+   if ( ( plPlane.GetVersN().IsZplus()  &&  dArea > 0)  ||  ( plPlane.GetVersN().IsZminus()  &&  dArea < 0))
      pMyCrv->Invert() ;
  // verifico non abbia auto-intersezioni
   SelfIntersCurve sInt( *pMyCrv) ;
@@ -2867,7 +2867,7 @@ SurfTriMesh::Cut( const Plane3d& plPlane, bool bSaveOnEq)
      double dDist = DistPointPlane( m_vVert[i].ptP, plPlane) ;
      if ( abs( dDist) < EPS_SMALL) {
         m_vVert[i].nTemp = 0 ;
-         m_vVert[i].ptP -= plPlane.vtN * dDist ;
+         m_vVert[i].ptP -= plPlane.GetVersN() * dDist ;
      }
      else if ( dDist > 0)
         m_vVert[i].nTemp = +1 ;
@@ -2891,7 +2891,7 @@ SurfTriMesh::Cut( const Plane3d& plPlane, bool bSaveOnEq)
           m_vVert[ m_vTria[i].nIdVert[0]].nTemp == 0  &&
           m_vVert[ m_vTria[i].nIdVert[1]].nTemp == 0  &&
           m_vVert[ m_vTria[i].nIdVert[2]].nTemp == 0  &&
-           AreSameVectorApprox( m_vTria[i].vtN, plPlane.vtN))
+           AreSameVectorApprox( m_vTria[i].vtN, plPlane.GetVersN()))
         ;
     // se giace sul piano o dalla parte esterna del piano, lo cancello
      else if ( m_vVert[ m_vTria[i].nIdVert[0]].nTemp != -1  &&
@@ -60,7 +60,7 @@ SurfTriMesh::UpdateFaceting( void)
         ++ nFacet ;
        // piano del triangolo
         Plane3d plPlane ;
-         if ( ! SetPlane( m_vVert[m_vTria[i].nIdVert[0]].ptP, m_vTria[i].vtN, plPlane))
+         if ( ! plPlane.Set( m_vVert[m_vTria[i].nIdVert[0]].ptP, m_vTria[i].vtN))
            return false ;
        // aggiorno i triangoli adiacenti
         for ( int j = 0 ; j < 3 ; ++ j) {
@@ -88,7 +88,7 @@ bool
 SurfTriMesh::UpdateTriaFaceting( int nRefT, int nFacet, const Plane3d& plPlane, int nT)
 {
  // verifica scostamento della normale da quella del piano
-   if ( plPlane.vtN * m_vTria[nT].vtN < m_dCosBndAng)
+   if ( plPlane.GetVersN() * m_vTria[nT].vtN < m_dCosBndAng)
      return true ;
  // verifica scostamento del vertice opposto al lato in comune dal piano
   int nV = SVT_NULL ;
@@ -411,7 +411,7 @@ SurfTriMesh::GetFacetLoops( int nF, POLYLINEVECTOR& vPL) const
      if ( ! vPL[i].IsClosedAndFlat( plPlane, dArea))
         return false ;
     // se loop esterno
-      if ( vtN * plPlane.vtN > 0) {
+      if ( vtN * plPlane.GetVersN() > 0) {
        // se non c'è ancora loop esterno in prima posizione
         if ( ! bOutFirst) {
           // lo sposto in prima posizione
@@ -44,7 +44,7 @@ Triangulate::Make( const PolyLine& PL, PNTVECTOR& vPt, INTVECTOR& vTr)
      return false ;
  // determino il piano ottimale di proiezione e il relativo senso di rotazione
   bool bCCW ;
-   if ( ! CalcProjPlane( plPlane.vtN, m_nPlane, bCCW))
+   if ( ! CalcProjPlane( plPlane.GetVersN(), m_nPlane, bCCW))
      return false ;
  // riempio il vettore con i vertici del poligono da triangolare
   vPt.clear() ;
@@ -103,7 +103,7 @@ Triangulate::Make( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr)
      return false ;
  // determino il piano ottimale di proiezione e il relativo senso di rotazione
   bool bCCW ;
-   if ( ! CalcProjPlane( plExtPlane.vtN, m_nPlane, bCCW))
+   if ( ! CalcProjPlane( plExtPlane.GetVersN(), m_nPlane, bCCW))
      return false ;
  // verifico le altre polilinee
   for ( int i = 1 ; i < int( vPL.size()) ; ++i) {
@@ -111,7 +111,7 @@ Triangulate::Make( const POLYLINEVECTOR& vPL, PNTVECTOR& vPt, INTVECTOR& vTr)
      double dArea ;
      Plane3d plPlane ;
      if ( ! vPL[i].IsClosedAndFlat( plPlane, dArea, 50 * EPS_SMALL)  ||
-           ! AreOppositeVectorApprox( plExtPlane.vtN, plPlane.vtN))
+           ! AreOppositeVectorApprox( plExtPlane.GetVersN(), plPlane.GetVersN()))
         return false ;
   }
  // se non CCW inverto tutte le polilinee
@@ -18,16 +18,25 @@
 #include "GeoConst.h"
 #include "IntersLineSurfTm.h"
 #include "/EgtDev/Include/EgtNumUtils.h"
+#include "/EgtDev/Include/EGkIntersLineTria.h"

 using namespace std ;

+//----------------------------------------------------------------------------
+bool
+VolZmap::IntersLineBox( const Point3d& ptP, const Vector3d& vtV, const Point3d& ptMin, const Point3d& ptMax) const 
+{
+  // Il box è allineato agli assi
+   double dU1, dU2 ;
+   return IntersLineBox( ptP, vtV, ptMin, ptMax, dU1, dU2) ;
+}
+
 //----------------------------------------------------------------------------
 bool
 VolZmap::IntersLineBox( const Point3d& ptP, const Vector3d& vtV,
                        const Point3d& ptMin, const Point3d& ptMax, double& dU1, double& dU2) const
 {
  // Il box è allineato agli assi
-
   dU1 = - INFINITO ;
   dU2 = INFINITO ;

@@ -168,6 +177,11 @@ VolZmap::IntersRayDexel( const Point3d& ptP, const Vector3d& vtV, unsigned int n
 bool
 VolZmap::GetDepth( const Point3d& ptPGlob, const Vector3d& vtDir, double& dInLength, double& dOutLength) const
 {
+#if 1
+   GetDepthWithVoxel( ptPGlob, vtDir, dInLength, dOutLength) ;
+   return true ;
+#else
+
   int nGrid = 0 ;

  // Porto il raggio nel riferimento intrinseco
@@ -271,6 +285,190 @@ VolZmap::GetDepth( const Point3d& ptPGlob, const Vector3d& vtDir, double& dInLen
   if ( dInLength < - EPS_SMALL)
      dInLength = - 1 ;
     
+   return true ;
+#endif
+}
+
+//----------------------------------------------------------------------------
+// Punti e vettori devono essere espressi nel sistema locale (quello in cui è immerso lo Zmap)
+// InLength = distanza di ingresso (se -1 il punto è interno, se -2 il punto è esterno e il raggio non interseca lo Zmap)
+// OutLength = distanza di uscita
+bool
+VolZmap::GetDepthWithVoxel( const Point3d& ptPLoc, const Vector3d& vtDir, double& dInLength, double& dOutLength) const 
+{
+  // Porto punto e vettore della retta nel sistema Zmap
+   Point3d ptP = ptPLoc ;
+   ptP.ToLoc( m_MapFrame[0]) ;
+   Vector3d vtD = vtDir ;
+   vtD.ToLoc( m_MapFrame[0]) ;
+   vtD.Normalize() ;
+
+  // Parametri di intersezione retta BBox dello Zmap
+   double dU1, dU2 ;
+  
+  // Intersezione fra semiretta e BBox dello Zmap    
+   bool bLineBBoxInters = IntersLineZMapBBox( ptP, vtD, 0, dU1, dU2)  &&  ( dU1 > 0  ||  dU2 > 0) ;
+  // Semiretta esterna al box dello Zmap quindi esterna anche allo Zmap
+   if ( ! bLineBBoxInters) {
+      dInLength = - 2 ; 
+      dOutLength = - 2 ;
+      return true ;
+   }
+
+  // Indici del voxel corrente
+   int nVoxI, nVoxJ, nVoxK ;
+  // Determino il voxel di partenza
+   if ( ! GetPointVoxel( ptP, nVoxI, nVoxJ, nVoxK)) { 
+      if ( ! GetPointVoxel( ptP + dU1 * vtD, nVoxI, nVoxJ, nVoxK)) 
+         return false ;
+   }
+  
+  // Indici dell'ultimo voxel
+   int nVxEndI, nVxEndJ, nVxEndK ;
+  // Determino il voxel finale
+   if ( ! GetPointVoxel( ptP + dU2 * vtD, nVxEndI, nVxEndJ, nVxEndK))
+       return false ;
+
+  // Struttura studio dell'intersezione
+   struct LineTriaInt {
+      double dPar ;
+      double dDot ; 
+   } ;
+   std::vector<std::vector<LineTriaInt>> IntMatrix ;
+
+   int nStI = min( nVoxI, nVxEndI) ;
+   int nStJ = min( nVoxJ, nVxEndJ) ;
+   int nStK = min( nVoxK, nVxEndK) ; 
+   int nEnI = max( nVoxI, nVxEndI) ;
+   int nEnJ = max( nVoxJ, nVxEndJ) ;
+   int nEnK = max( nVoxK, nVxEndK) ; 
+ 
+   for ( int nI = nStI ; nI <= nEnI ; ++ nI) {
+      for ( int nJ = nStJ ; nJ <= nEnJ ; ++ nJ) {
+         for ( int nK = nStK ; nK <= nEnK ; ++ nK) {
+
+            Point3d ptMin( ( nI + 0.5) * m_dStep, 
+                           ( nJ + 0.5) * m_dStep, 
+                           ( nK + 0.5) * m_dStep) ;
+            Point3d ptMax( ( nI + 1.5) * m_dStep, 
+                           ( nJ + 1.5) * m_dStep, 
+                           ( nK + 1.5) * m_dStep) ;
+
+            if ( IsValidVoxel( nI, nJ, nK)  &&  
+                 IntersLineBox( ptP, vtD, ptMin, ptMax)) {
+            
+              // Analisi del voxel
+               int nCbType ;
+               TRIA3DLIST lstTria ;
+               ProcessCube( nI, nJ, nK, nCbType, lstTria) ;
+
+              // Se il voxel contiene triangoli
+               if ( nCbType == VOX_ON_BOUNDARY) {
+      
+                 // Ciclo sui triangoli del voxel
+                  for ( auto it = lstTria.begin() ; it != lstTria.end() ; ++it ) {
+        
+                     Triangle3d CurrTria = *it ;
+                     Point3d ptLineTria1, ptLineTria2 ;
+                    // Studio dell'intersezione della retta con il triangolo corrente
+                     int nIntType = IntersLineTria( ptP, vtD, 1.5 * dU2, CurrTria, ptLineTria1, ptLineTria2) ;
+                    // Se non ci sono intersezioni passo al prossimo triangolo
+                     if ( nIntType == ILTT_NO)     
+                        continue ;
+                    // se altrimenti c'è una sola intersezione
+                     else if ( nIntType == ILTT_VERT  ||   
+                               nIntType == ILTT_EDGE  ||
+                               nIntType == ILTT_IN) {           
+                        LineTriaInt NewInt ;
+                        NewInt.dPar = ( ptLineTria1 - ptP) * vtD ;
+                        NewInt.dDot = vtD * CurrTria.GetN() ;
+                        std::vector<LineTriaInt> vSing ;
+                        vSing.emplace_back( NewInt) ;
+              
+                        IntMatrix.emplace_back( vSing) ;
+                     }
+                    // altrimenti ci sono due intersezioni
+                     else {
+                        LineTriaInt NewInt1, NewInt2 ;
+                        std::vector<LineTriaInt> vCouple ;
+                        double dP1 = ( ptLineTria1 - ptP) * vtD ;
+                        double dP2 = ( ptLineTria2 - ptP) * vtD ;
+                        NewInt1.dPar = ( dP1 < dP2 ? dP1 : dP2) ;
+                        NewInt2.dPar = ( dP1 < dP2 ? dP2 : dP1) ;        
+                        NewInt1.dDot = vtD * CurrTria.GetN() ;
+                        NewInt2.dDot = NewInt1.dDot ;
+                                   
+                        vCouple.emplace_back( NewInt1) ;
+                        vCouple.emplace_back( NewInt2) ;
+                        IntMatrix.emplace_back( vCouple) ;
+                     }
+                  }
+               }
+            }
+         }
+      }
+   }
+  
+  // Ordino le intersezioni in base al parametro distanza con segno da ptP
+   for ( int nN = 0 ; nN < int( IntMatrix.size()) - 1 ; ++ nN) {
+
+      double dFP = ( IntMatrix[nN].size() == 2 ? 0.5 * ( IntMatrix[nN][0].dPar + IntMatrix[nN][1].dPar) : 
+                     IntMatrix[nN][0].dPar) ;
+      double dLP = ( IntMatrix[nN+1].size() == 2 ? 0.5 * ( IntMatrix[nN+1][0].dPar + IntMatrix[nN+1][1].dPar) : 
+                     IntMatrix[nN+1][0].dPar) ;
+      
+      if ( dFP > dLP) 
+         swap( IntMatrix[nN], IntMatrix[nN+1]) ;
+   }
+
+   std::vector<LineTriaInt> vInt ;
+   
+   for ( int nN = 0 ; nN < int( IntMatrix.size()) ; ++ nN) {
+      vInt.emplace_back( IntMatrix[nN][0]) ;
+      if ( IntMatrix[nN].size() == 2) 
+         vInt.emplace_back( IntMatrix[nN][1]) ;
+   }
+
+  // Inizializzo le distanze di ingresso e uscita:
+  // dInLength diminuisce, dOutLength aumenta.
+   dInLength = INFINITO ; 
+   dOutLength = - INFINITO ;
+
+   int nFirstPosN ;
+   int nN = 0 ;
+   for ( ; nN < int( vInt.size()) ; ++ nN) {  
+      if ( vInt[nN].dPar >= 0) {
+         nFirstPosN = nN ;
+         break ;   
+      }
+   }
+
+   if ( nN == int( vInt.size())) {
+      dInLength = - 2 ;
+      dOutLength = - 2 ;
+      return true ;
+   }
+
+   if ( nFirstPosN > 0) {   
+      if ( vInt[nFirstPosN - 1].dDot < EPS_ZERO)
+         dInLength = -1 ;    
+   }
+   else if ( nFirstPosN == 0) {
+      if ( vInt[nFirstPosN].dDot > - EPS_ZERO)
+         dInLength = -1 ;    
+   }
+
+   for ( int nN = nFirstPosN ; nN < int( vInt.size()) ; ++ nN) {
+   
+      if ( vInt[nN].dDot > - EPS_ZERO  && 
+           dOutLength < vInt[nN].dPar)
+         dOutLength = vInt[nN].dPar ;
+     
+      if ( vInt[nN].dDot < EPS_ZERO  && 
+           dInLength > vInt[nN].dPar)
+         dInLength = vInt[nN].dPar ; 
+   }
+
   return true ;
 }

@@ -1161,3 +1359,5 @@ VolZmap::GetPartVolume( int nPart, double& dVol) const

   return true ;
 }
+
+
@@ -614,6 +614,285 @@ VolZmap::AddDexelSideFace( int nPos, int nPosAdj, const Point3d& ptP, const Poin
   return true ;
 }

+//----------------------------------------------------------------------------
+// I triangoli sono espressi nel sistema Zmap
+bool
+VolZmap::ProcessCube( int nVoxI, int nVoxJ, int nVoxK, int& nCubeIndex, TRIA3DLIST& lstTria) const 
+{
+  // Se il voxel non esiste, vi è un errore.
+   if ( nVoxI + 1 < 0  ||  nVoxI + 1 > int( m_nNx[0])  ||
+        nVoxJ + 1 < 0  ||  nVoxJ + 1 > int( m_nNy[0])  ||
+        nVoxK + 1 < 0  ||  nVoxK + 1 > int( m_nNy[1]))
+      return false ;
+
+  // Classificazione dei vertici: interni o esterni al materiale o di frontiera 
+   int nIndex = CalcIndex( nVoxI, nVoxJ, nVoxK) ;
+   if ( nIndex == 0)
+      nCubeIndex = VOX_EXTERN ;
+   else if ( nIndex == 255)
+      nCubeIndex = VOX_INNER ;
+   else 
+      nCubeIndex = VOX_ON_BOUNDARY ;
+
+   if ( nCubeIndex == VOX_ON_BOUNDARY) {
+
+     // Indici i,j,k dei vertici
+      int IndexCorner[8][3] = {
+         { nVoxI, nVoxJ, nVoxK}, 
+         { nVoxI + 1, nVoxJ, nVoxK}, 
+         { nVoxI + 1, nVoxJ + 1, nVoxK}, 
+         { nVoxI, nVoxJ + 1, nVoxK}, 
+         { nVoxI, nVoxJ, nVoxK + 1}, 
+         { nVoxI + 1, nVoxJ, nVoxK + 1}, 
+         { nVoxI + 1, nVoxJ + 1, nVoxK + 1}, 
+         { nVoxI, nVoxJ + 1, nVoxK + 1} 
+      } ;
+
+      static int intersections[12][2] = {
+         { 0, 1 }, { 1, 2 }, { 3, 2 }, { 0, 3 }, { 4, 5 }, { 5, 6 },
+         { 7, 6 }, { 4, 7 }, { 0, 4 }, { 1, 5 }, { 2, 6 }, { 3, 7 }
+      } ;
+
+     // Array di strutture punto di intersezione e normale alla superficie in esso.
+      VectorField VecField[12] ;
+
+     // Flag di regolarità dei campi scalare e vettoriale  
+      bool bReg = true ;
+
+     // Ciclo sui segmenti
+      for ( int EdgeIndex = 0 ; EdgeIndex < 12 ; ++ EdgeIndex) {
+        // Se il segmento non attraversa la superficie passo al successivo
+         if ( ( EdgeTable[nIndex] & ( 1 << EdgeIndex)) == 0)
+            continue ;
+        // Indici per linee di griglia sui vertici
+         int n1 = intersections[EdgeIndex][0] ;
+         int n2 = intersections[EdgeIndex][1] ;
+        // Flag posizione corner
+         bool bN1 = ( ( nIndex & ( 1 << n1)) != 0) ;
+        // Determino con precisione il punto di intersezione sullo spigolo,
+        // se i campi scalare e vettoriale non sono regolari bReg diviene falso.
+         if ( ! IntersPos( IndexCorner[n1], IndexCorner[n2], bN1,
+                           VecField[EdgeIndex].ptInt, VecField[EdgeIndex].vtNorm))
+            bReg = false ;
+      }
+
+     // Determino il numero di componenti connesse nel voxel in caso di configurazione standard
+      int nComponents = TriangleTableEn[nIndex][1][0] ;
+
+     // Matrici di campi vettoriali:
+     // CompoVert[i] ha i vertici della base del triangle fan della (i+1)-esima componente connessa;
+     // CompoTriVert[i] ha i vertici di tutti i triangoli, nel caso di assenza di sharp feature, 
+     // della (i+1)-esima componente connessa.
+      VectorField CompoVert[6][12] ;
+      VectorField CompoTriVert[6][17] ;
+
+     // Arrey numero di vertici della base del fan per componente
+     // connessa: nVertComp[i] contiene il numero di vertici 
+     // della base del fan della (i+1)-esima componente connessa. 
+      int nVertComp[6] ;
+
+     
+      int nExtTabOff = nComponents ;
+      int nStdTabOff = 0 ;
+
+     // Carico le matrici CompoVert e CompoTriVert  
+      for ( int nCompCount = 1 ; nCompCount <= nComponents ; ++ nCompCount) {
+
+        // Numero vertici per componenti
+         nVertComp[nCompCount - 1] = TriangleTableEn[nIndex][1][nCompCount] ;
+  
+        // Riempio il nCompCount-esimo vettore di vertici dei triangoli in assenza di  
+        // sharp feature: in una mesh di triangoli con n vertici vi sono n - 2 triangoli.
+         for ( int nVertCount = 0 ; nVertCount < 3 * ( nVertComp[nCompCount - 1] - 2) ; nVertCount += 3) { 
+            CompoTriVert[nCompCount - 1][nVertCount] =   VecField[TriangleTableEn[nIndex][0][nStdTabOff + nVertCount+2]] ;
+            CompoTriVert[nCompCount - 1][nVertCount+1] = VecField[TriangleTableEn[nIndex][0][nStdTabOff + nVertCount+1]] ;
+            CompoTriVert[nCompCount - 1][nVertCount+2] =   VecField[TriangleTableEn[nIndex][0][nStdTabOff + nVertCount]] ;
+         }
+
+        // Aggiorno gli offsets per raggiungere i
+        // vertici della componente successiva.
+         nExtTabOff += nVertComp[nCompCount - 1] ;
+         nStdTabOff += 3 * ( nVertComp[nCompCount - 1] - 2) ;
+      }
+
+     // Test sulla topologia: dal momento che il nostro test
+     // si fonda sugli angoli compresi fra le normali, esso ha 
+     // senso solo se il campo è regolare.    
+      if ( bReg  &&  nAllConfig[nIndex] == 3) { 
+
+         Vector3d vtCmpAvg0, vtCmpAvg1 ; 
+
+        // Verifico se i versori delle componenti sono tutti 
+        // più o meno concordi (per esserlo non devono esserci
+        // due vettori di una medesima componente con prodotto 
+        // scalare inferiore a 0.7).
+         bool bTest0 = DotTest( CompoVert[0], 3, vtCmpAvg0, 0.7) ;
+         bool bTest1 = DotTest( CompoVert[1], 3, vtCmpAvg1, 0.7) ;
+
+        // Se i versori di entrambe le componenti sono concordi
+        // ha senso parlare di vettori medi, altrimenti non ha 
+        // senso. Se non ha senso parlare di vettori medi non 
+        // ha senso parlare di prodotti scalari fra loro, 
+        // quindi pongo il loro prodotto a un valore assurdo -2
+        // (il prodotto scalare fra versori ha modulo non superiore
+        // a uno).
+         double dScProd = - 2 ;
+
+         if ( bTest0  &&  bTest1) 
+            dScProd = vtCmpAvg0 * vtCmpAvg1 ;
+
+         double dThreshold = 0.7 ; 
+
+         if ( ! ( bTest0  &&  bTest1)  ||  ( bTest0  &&  bTest1  &&  dScProd > dThreshold)) {
+
+            int nt = 0 ;
+
+            while ( nIndexVsIndex3[nt][0] != nIndex)
+               ++ nt ;
+
+            int nRotCase = nIndexVsIndex3[nt][1] ;
+
+            nComponents = Cases3Plus[nRotCase][1][0] ;
+
+           // Riaggiorno gli offsets
+            nExtTabOff = nComponents ;
+            nStdTabOff = 0 ;
+
+           // Modifico le matrici
+            for ( int nC = 1 ; nC <= nComponents ; ++ nC) {
+
+              // Numero vertici per componenti
+               nVertComp[nC - 1] = Cases3Plus[nRotCase][1][nC] ;            
+
+              // Matrici dei vertici dei triangoli in assenza di sharp feature 
+               for ( int nTriVert = 0 ; nTriVert < 3 * ( nVertComp[nC - 1] - 2) ; nTriVert += 3) { 
+
+                  CompoTriVert[nC - 1][nTriVert] =   VecField[Cases3Plus[nRotCase][0][nStdTabOff + nTriVert+2]] ;
+                  CompoTriVert[nC - 1][nTriVert+1] = VecField[Cases3Plus[nRotCase][0][nStdTabOff + nTriVert+1]] ;
+                  CompoTriVert[nC - 1][nTriVert+2] =   VecField[Cases3Plus[nRotCase][0][nStdTabOff + nTriVert]] ;
+               }
+
+              // Aggiorno gli offsets per raggiungere i 
+              // vertici della componente successiva.
+               nExtTabOff += nVertComp[nC - 1] ;
+               nStdTabOff += 3 * ( nVertComp[nC - 1] - 2) ;
+            }
+         }
+      }
+
+      else if ( bReg  &&  nAllConfig[nIndex] == 6) {
+
+        // Procedura analoga a quella della configurazione 3
+         Vector3d vtCmpAvg0, vtCmpAvg1 ;
+
+         bool bTest0 = DotTest( CompoVert[0], 4, vtCmpAvg0, 0.7) ;
+         bool bTest1 = DotTest( CompoVert[1], 3, vtCmpAvg1, 0.7) ;
+
+         double dScProd = - 2 ;
+
+         if ( bTest0  &&  bTest1) 
+            dScProd = vtCmpAvg0 * vtCmpAvg1 ;
+
+         double dThreshold = 0.7 ; 
+
+         if ( ! ( bTest0  &&  bTest1)  ||  ( bTest0  &&  bTest1  &&  dScProd > dThreshold)) {
+
+            int nt = 0 ;
+
+            while ( nIndexVsIndex6[nt][0] != nIndex)
+               ++ nt ;
+
+            int nRotCase = nIndexVsIndex6[nt][1] ;
+
+            nComponents = 1 ;
+
+           // Riaggiorno gli offsets
+            nExtTabOff = nComponents ;
+            nStdTabOff = 0 ;
+
+           // Modifico le matrici
+            for ( int nC = 1 ; nC <= nComponents ; ++ nC) {
+
+              // Numero vertici per componenti
+               nVertComp[nC - 1] = 7 ;            
+
+              // Matrici dei vertici dei triangoli in assenza di sharp feature 
+               for ( int nTriVert = 0 ; nTriVert < 3 * ( nVertComp[nC - 1] - 2) ; nTriVert += 3) { 
+
+                  CompoTriVert[nC - 1][nTriVert] =   VecField[Cases3Plus[nRotCase][0][nStdTabOff + nTriVert+2]] ;
+                  CompoTriVert[nC - 1][nTriVert+1] = VecField[Cases3Plus[nRotCase][0][nStdTabOff + nTriVert+1]] ;
+                  CompoTriVert[nC - 1][nTriVert+2] =   VecField[Cases3Plus[nRotCase][0][nStdTabOff + nTriVert]] ;
+               }
+
+              // Aggiorno gli offsets per raggiungere i 
+              // vertici della componente successiva.
+               nExtTabOff += nVertComp[nC - 1] ;
+               nStdTabOff += 3 * ( nVertComp[nC - 1] - 2) ;
+            }
+         }
+      }
+
+      else if ( bReg  &&  nAllConfig[nIndex] == 10) {
+
+         Vector3d vtCmpAvg0, vtCmpAvg1 ; 
+
+        // Verifico se i versori delle componenti sono tutti 
+        // più o meno concordi (per esserlo non devono esserci
+        // due vettori di una medesima componente con prodotto 
+        // scalare inferiore a 0).  decidere se 0.0 o 0.7
+         bool bTest0 = DotTest( CompoVert[0], 4, vtCmpAvg0) ;
+         bool bTest1 = DotTest( CompoVert[1], 4, vtCmpAvg1) ;
+
+         if ( ! ( bTest0 && bTest1)) {
+
+            int nt = 0 ;
+            while ( nIndexVsIndex10[nt][0] != nIndex)
+               ++ nt ;
+
+            int nRotCase = nIndexVsIndex10[nt][1] ;
+
+           // Riaggiorno gli offsets
+            nExtTabOff = 2 ;
+            nStdTabOff = 0 ;
+
+           // Modifico le matrici
+            for ( int nC = 1 ; nC <= 2 ; ++ nC) {
+              // Numero vertici per componenti
+               nVertComp[nC - 1] = Cases10Plus[nRotCase][1][nC] ;
+
+               
+              // Matrici dei vertici dei triangoli in assenza di sharp feature 
+               for ( int nTriVert = 0 ; nTriVert < 6 ; nTriVert += 3) { 
+                  CompoTriVert[nC - 1][nTriVert] =   VecField[Cases10Plus[nRotCase][0][nStdTabOff + nTriVert+2]] ;
+                  CompoTriVert[nC - 1][nTriVert+1] = VecField[Cases10Plus[nRotCase][0][nStdTabOff + nTriVert+1]] ;
+                  CompoTriVert[nC - 1][nTriVert+2] = VecField[Cases10Plus[nRotCase][0][nStdTabOff + nTriVert]] ;
+               }
+
+              // Aggiorno gli offsets per raggiungere i vertici della componente successiva.
+               nExtTabOff += nVertComp[nC - 1] ;
+               nStdTabOff += 3 * ( nVertComp[nC - 1] - 2) ;
+            }
+         }
+      }  
+     // Ciclo sulle componenti
+      for ( int nCompCount = 1 ; nCompCount <= nComponents ; ++ nCompCount) {
+        // Costruzione dei triangoli
+         for ( int TriIndex =  0; TriIndex < ( nVertComp[nCompCount - 1] - 2) * 3 ; TriIndex += 3) {
+           // Il triangolo è pronto
+            Triangle3d CurrentTriangle ;
+            CurrentTriangle.Set( CompoTriVert[nCompCount - 1][TriIndex].ptInt,
+                                 CompoTriVert[nCompCount - 1][TriIndex+1].ptInt, 
+                                 CompoTriVert[nCompCount - 1][TriIndex+2].ptInt) ;
+            bool bV = CurrentTriangle.Validate( true) ;
+           // Aggiungo alla lista
+            lstTria.emplace_back( CurrentTriangle) ;
+         }          
+      }
+   }
+
+   return true ;
+}
+
 //----------------------------------------------------------------------------
 bool
 VolZmap::MarchingCubes( int nBlock, TRIA3DLIST& lstTria) const
@@ -1289,7 +1568,7 @@ VolZmap::ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold)
               }

              // Flag ExtMC
-               bool bExtMC = ( nFeatureType != NO_FEATURE) && bGridControl ;
+               bool bExtMC = ( nFeatureType != NO_FEATURE  &&  bGridControl) ;

              // Extended MC
               if ( bExtMC) {
@@ -2510,6 +2789,22 @@ VolZmap::IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, Point3d& ptInt,
   return bFound ;
 }

+//----------------------------------------------------------------------------
+bool
+VolZmap::IsValidVoxel( int nN) const
+{
+   return ( nN >= 0  &&  nN < int( ( m_nNx[0] + 1) * ( m_nNy[0] + 1) * ( m_nNy[1] + 1))) ;
+}
+
+//----------------------------------------------------------------------------
+bool 
+VolZmap::IsValidVoxel( int nI, int nJ, int nK) const 
+{
+   return ( nI >= - 1  &&  nI < int( m_nNx[0])  &&
+            nJ >= - 1  &&  nJ < int( m_nNy[0])  &&
+            nK >= - 1  &&  nK < int( m_nNy[1])  ) ;
+}
+
 //----------------------------------------------------------------------------
 bool
 VolZmap::GetVoxIJKFromN( int nN, int& nI, int& nJ, int& nK) const
@@ -24,7 +24,7 @@
 #include "/EgtDev/Include/EgkIntersLinesurfTm.h"
 #include <unordered_map>
 #include <stack>
-      
+
 //----------------------------------------------------------------------------
 struct TriaStruct {
   int i, j, k ;
@@ -37,7 +37,7 @@ typedef std::vector<TriaStruct> TriHolder ;

 // Vettori di TriHolder con sharp feature di frontiera:
 // il primo indice individua il blocco, il secondo il voxel
-typedef std::vector<TriHolder> TriaMatrix ;     
+typedef std::vector<TriHolder> TriaMatrix ;

 //----------------------------------------------------------------------------
 class VolZmap : public IVolZmap, public IGeoObjRW
@@ -97,6 +97,7 @@ class VolZmap : public IVolZmap, public IGeoObjRW
      bool MillingStep( const Point3d& ptPs, const Vector3d& vtDs, const Point3d& ptPe, const Vector3d& vtDe) override ;
      bool MillingStep( const Point3d& ptPs, const Vector3d& vtDs, const Vector3d& vtAs, const Point3d& ptPe, const Vector3d& vtDe, const Vector3d& vtAe) override ;
      bool GetDepth( const Point3d& ptP, const Vector3d& vtDir, double& dInLength, double& dOutLength) const override ;
+      bool GetDepthWithVoxel( const Point3d& ptPLoc, const Vector3d& vtDir, double& dInLength, double& dOutLength) const ;
      bool AvoidBox( const Frame3d& frBox, const Vector3d& vtDiag) const override ;
      VolZmap* ClonePart( int nPart) const override ;
      bool RemovePart( int nPart) override ;
@@ -124,6 +125,9 @@ class VolZmap : public IVolZmap, public IGeoObjRW
                      CONUSMILL = 4,               // con parte terminale conica
                      MORTISER = 5,                // mortasatrice
                      CHISEL = 6} ;                // scalpello
+      enum CubeType { VOX_EXTERN = 1,
+                      VOX_ON_BOUNDARY = 0,
+                      VOX_INNER = -1} ;
      typedef std::unordered_map <int, double> VoxelContainer ;
      
   private :
@@ -134,6 +138,7 @@ class VolZmap : public IVolZmap, public IGeoObjRW
      bool CalcDexelPrisms( int nPos1, int nPos2, TRIA3DLIST& lstTria) const ;
      bool AddDexelSideFace( int nPos, int nPosAdj, const Point3d& ptP, const Point3d& ptQ,
                             const Vector3d& vtZ, const Vector3d& vtNorm, TRIA3DLIST& lstTria) const ;
+      bool ProcessCube( int nVoxI, int nVoxJ, int nVoxK, int& nCubeIndex, TRIA3DLIST& lstTria) const ;
      bool MarchingCubes( int nBlock, TRIA3DLIST& lstTria) const ;
      bool ExtMarchingCubes( int nBlock, TRIA3DLIST& lstTria, TriHolder& triHold) const ; 
      bool FlipEdgesII( TriHolder& TriHold) const ;
@@ -144,7 +149,6 @@ class VolZmap : public IVolZmap, public IGeoObjRW
      bool IntersPos( int nVec1[], int nVec2[], bool bFirstCorner, Point3d& ptInt, Vector3d& vtNormal) const ;
      bool IsPointInsideVoxelApprox( int nI, int nJ, int nK, const Point3d& ptP, double dPrec = EPS_SMALL) const ;
      bool GetPointVoxel( const Point3d& ptP, int& nVoxI, int& nVoxJ, int& nVoxK) const ;
-           
     // OPERAZIONI SU INTERVALLI
      bool SubtractIntervals( unsigned int nGrid, unsigned int nI, unsigned int nJ,
                              double dMin, double dMax, const Vector3d& vtNMin, const Vector3d& vtNMax) ;
@@ -232,6 +236,7 @@ class VolZmap : public IVolZmap, public IGeoObjRW
                               const Point3d& ptS, const Point3d& ptE, const Vector3d& vtD, const Vector3d& vtA,
                               unsigned int& nStI, unsigned int& nStJ, unsigned int& nEnI, unsigned int& nEnJ) ;
    // Intersezioni
+      bool IntersLineBox( const Point3d& ptP, const Vector3d& vtV, const Point3d& ptMin, const Point3d& ptMax) const ;
      bool IntersLineBox( const Point3d& ptP, const Vector3d& vtV, const Point3d& ptMin, const Point3d& ptMax,
                          double& dU1, double& dU2) const ;
      bool IntersLineZMapBBox( const Point3d& ptP, const Vector3d& vtV, unsigned int nGrid, double& dU1, double& dU2) const ;
@@ -260,7 +265,9 @@ class VolZmap : public IVolZmap, public IGeoObjRW
      bool IntersLineMyPolyhedron( const Point3d& ptLineSt, const Vector3d& vtLineDir,
                                   const Frame3d& PolyFrame, double dLenX, double dLenY, double dLenZ, double dDeltaX,
                                   Point3d& ptInt1, Point3d& ptInt2, Vector3d& vtN1, Vector3d& vtN2) ;
-     // Passaggio da N a ijk per i voxel
+     // Voxel: esistenza e passaggio da N a ijk per i voxel
+      bool IsValidVoxel( int nN) const ;
+      bool IsValidVoxel( int nI, int nJ, int nK) const ;     
      bool GetVoxIJKFromN( int nN, int& nI, int& nJ, int& nK) const ;
      bool GetVoxNFromIJK( int nI, int nJ, int nK, int& nN) const ;
     // Funzioni di gestione dei blocchi
@@ -271,7 +278,7 @@ class VolZmap : public IVolZmap, public IGeoObjRW
      bool GetAdjBlockToBlock( int nBlockN, int nDeltaI, int nDeltaJ, int nDeltaK, int& nAdjBlockN) const ;
      bool IsAVoxelOnBoundary( const int nLimits[], const int nIJK[], bool bType) const ;
      bool IsATriangleOnBorder( const Triangle3d& trTria, const Point3d& ptVert, 
-                                const int nBlockLimits[], const int nVoxIJK[]) const ;
+                                const int nBlockLimits[], const int nVoxIJK[]) const ;    
     // Funzioni per facce canoniche con grandi triangoli
      bool ProcessVoxContXY( VoxelContainer& VoxContXY, bool bPlus, TRIA3DLIST& lstTria) const ;
      bool ProcessVoxContYZ( VoxelContainer& VoxContYZ, bool bPlus, TRIA3DLIST& lstTria) const ;