//----------------------------------------------------------------------------
//  EgalTech 2015-2016
//----------------------------------------------------------------------------
// File : VolZmap.cpp         Data : 22.01.15           Versione : 1.6a4
// Contenuto : Implementazione della classe Volume Zmap (tre griglie)
//
//
//
// Modifiche : 22.01.15 DS Creazione modulo.
//
//
//----------------------------------------------------------------------------

//--------------------------- Include ----------------------------------------
#include "stdafx.h"
#include "CurveLine.h"
#include "VolZmap.h"
#include "GeoConst.h"
#include "IntersLineBox.h"
#include "IntersLineSurfStd.h"
#include "/EgtDev/Include/EGkIntersLineTria.h"
#include "/EgtDev/Include/EGkIntersLinePlane.h"
#include "/EgtDev/Include/EGkIntersLineSphere.h"
#include "/EgtDev/Include/EGkIntersPlaneBox.h"
#include "/EgtDev/Include/EGkIntersPlaneTria.h"
#include "/EgtDev/Include/EGkChainCurves.h"
#include "/EgtDev/Include/ENkPolynomialRoots.h"
#include "/EgtDev/Include/EGnStringUtils.h"
#include "/EgtDev/Include/EgtNumUtils.h"
#include <future>
#include <atomic>

using namespace std ;

//----------------------------------------------------------------------------
// La retta è rappresentata da un punto e dal versore espressi nel riferimento dello Zmap.
// Se non vi è intersezione, viene restituito falso
bool 
VolZmap::IntersLineZMapLattice( const Point3d& ptP, const Vector3d& vtV, double& dU1, double& dU2) const 
{
  // Punti estremi della diagonale del box del reticolo
   Point3d ptMin( - 0.5 * m_dStep, - 0.5 * m_dStep, - 0.5 * m_dStep) ;
   Point3d ptMax( ( m_nNx[0] + 0.5) * m_dStep, ( m_nNy[0] + 0.5) * m_dStep, ( m_nNy[1] + 0.5) * m_dStep) ; 

   return IntersLineBox( ptP, vtV, ptMin, ptMax, dU1, dU2) ;
}

//----------------------------------------------------------------------------
// Trova le intersezioni di una retta, definita tramite un suo punto e il vettore della direzione,
// e il bounding box dello Zmap. Punto e vettore devono essere espressi nel sistema di riferimento
// intrinseco dello Zmap. Se vi è intersezione viene restitutito true, false altrimenti.
// Se vi è intersezione in dU1 e dU2 vengono restituiti i parametri a cui avvengono le intersezioni.
bool
VolZmap::IntersLineZMapBBox( const Point3d& ptP, const Vector3d& vtV, double& dU1, double& dU2) const 
{
  // Punti estremi del box dello Zmap
   Point3d ptMin = ORIG ;
   Point3d ptMax ;
   if ( m_nMapNum == 3)
      ptMax = ptMin + Vector3d( m_dMaxZ[1], m_dMaxZ[2], m_dMaxZ[0]) ;
   else
      ptMax = ptMin + Vector3d( ( m_nNx[0] + 0.5) * m_dStep, ( m_nNy[0] + 0.5) * m_dStep, m_dMaxZ[0]) ;

   return IntersLineBox( ptP, vtV, ptMin, ptMax, dU1, dU2) ;
}

//----------------------------------------------------------------------------
// Trova le intersezioni fra una retta e un dexel parallelepipedale (gli intervalli sono 
// parallelepipdi). Punto e vettore devono essere espressi nel sistema grilgia che coincide 
// con quello intrinseco dello Zmap solo nel caso della prima griglia. Per le altre griglie 
// è necessario permutare ciclicamente le coordinate.
bool 
VolZmap::IntersLineDexel( const Point3d& ptP, const Vector3d& vtV, int nGrid, int nI, int nJ,
                          double& dU1, double& dU2) const 
{ 
  // Determino l'indice del dexel e il numero di suoi intervalli
   int nDexelPos = nJ * m_nNx[nGrid] + nI ;
   int nDexelSize = int( m_Values[nGrid][nDexelPos].size()) ;

 // Se non c'è materiale non devo fare alcunché
  if ( nDexelSize == 0) 
     return false ;

  // Determino estremi nel piano XY intrinseco del dexel
   double dXmin = nI * m_dStep ;
   double dYmin = nJ * m_dStep ;
   double dXmax = ( nI + 1) * m_dStep ;
   double dYmax = ( nJ + 1) * m_dStep ;

  // ciclo sugli intervalli
   dU1 = INFINITO ;
   dU2 = - INFINITO ;
   bool bInters = false ;
   for ( int nIndex = 0 ; nIndex < nDexelSize ; nIndex += 1) {
     // estremi del box del singolo intervallo 
      Point3d ptE1( dXmin, dYmin, m_Values[nGrid][nDexelPos][nIndex].dMin) ;
      Point3d ptE2( dXmax, dYmax, m_Values[nGrid][nDexelPos][nIndex].dMax) ;
      double dP1, dP2 ;
      if ( IntersLineBox( ptP, vtV, ptE1, ptE2, dP1, dP2)) {
         bInters = true ;
         dU1 = min( dU1, dP1) ;
         dU2 = max( dU2, dP2) ;
      }
   }

   return bInters ;
}

//----------------------------------------------------------------------------
// Trova le intersezioni fra una semi-retta e un dexel parallelepipedale (gli intervalli sono 
// parallelepipdi). Punto e vettore devono essere espressi nel sistema grilgia che coincide 
// con quello intrinseco dello Zmap solo nel caso della prima griglia. Per le altre griglie 
// è necessario permutare ciclicamente le coordinate.
bool 
VolZmap::IntersRayDexel( const Point3d& ptP, const Vector3d& vtV, int nGrid, int nI, int nJ,
                         double& dU1, double& dU2) const 
{ 
  // Determino l'indice del dexel e il numero di suoi intervalli
   int nDexelPos = nJ * m_nNx[nGrid] + nI ;
   int nDexelSize = int( m_Values[nGrid][nDexelPos].size()) ;

 // Se non c'è materiale non devo fare alcunché
  if ( nDexelSize == 0) 
     return false ;

  // Determino estremi nel piano XY intrinseco del dexel
   double dXmin = nI * m_dStep ;
   double dYmin = nJ * m_dStep ;
   double dXmax = ( nI + 1) * m_dStep ;
   double dYmax = ( nJ + 1) * m_dStep ;

  // ciclo sugli intervalli
   dU1 = INFINITO ;
   dU2 = - INFINITO ;
   bool bInters = false ;
   for ( int nIndex = 0 ; nIndex < nDexelSize ; nIndex += 1) {
     // estremi del box del singolo intervallo 
      Point3d ptE1( dXmin, dYmin, m_Values[nGrid][nDexelPos][nIndex].dMin) ;
      Point3d ptE2( dXmax, dYmax, m_Values[nGrid][nDexelPos][nIndex].dMax) ;
      double dP1, dP2 ;
      if ( IntersLineBox( ptP, vtV, ptE1, ptE2, dP1, dP2)) {
         if ( dP2 >= 0) {
            bInters = true ;
            if ( dP1 >= 0)
               dU1 = min( dU1, dP1) ;
            else
               dU1 = - 1 ;
            dU2 = max( dU2, dP2) ;
         }
      }
   }

   return bInters ;
}

//----------------------------------------------------------------------------
// Calcola la profondità del materiale lungo un raggio, identificato da un punto e un versore.
// Punto e versore devono essere espressi nel sistema di riferimento locale, in cui è immerso quello dello 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::GetDepth( const Point3d& ptPLoc, const Vector3d& vtDLoc, double& dInLength, double& dOutLength, bool bExact) const
{
   if ( bExact  &&  m_nMapNum == 3)
      return GetDepthWithVoxel( ptPLoc, vtDLoc, dInLength, dOutLength) ;
   else {
      return GetDepthWithDexel( ptPLoc, vtDLoc, dInLength, dOutLength) ;
   }
}

//----------------------------------------------------------------------------
// Calcola la profondità del materiale usando i parallelepipedi corrispondenti ai segmenti di dexel di una griglia.
// Punto e versore devono essere espressi nel sistema di riferimento Zmap. Il vettore deve essere normalizzato.
// 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::GetDepthWithDexel( const Point3d& ptPLoc, const Vector3d& vtDLoc, double& dInLength, double& dOutLength) const
{ 
  // Punto e direzione nel riferimento intrinseco
   Point3d ptP = ptPLoc ;
   Vector3d vtD = vtDLoc ;
   ptP.ToLoc( m_MapFrame) ;
   vtD.ToLoc( m_MapFrame) ;
   vtD.Normalize() ;

  // Intersezione fra semiretta e BBox dello Zmap
   double dU1, dU2 ;
   bool bTest = IntersLineZMapBBox( ptP, vtD, dU1, dU2)  &&  ( dU1 > 0  ||  dU2 > 0) ;

  // Semiretta esterna al box dello Zmap quindi esterna anche allo Zmap
   if ( ! bTest) {
      dInLength = -2 ; 
      dOutLength = -2 ;
      return true ;
   }

   double dInLen[N_MAPS] ;
   double dOutLen[N_MAPS] ;

  // Ciclo sulle griglie
   for ( int nGrid = 0 ; nGrid < m_nMapNum ; ++ nGrid) {

      Point3d ptP0 = ptP ;
      Vector3d vtV0 = vtD ;

      Point3d ptI, ptF ;
     // Una sola intersezione valida ( punto interno, intersezione valida 2)
      if ( dU1 < 0  &&  dU2 > 0) {
         ptI = ptP ;
         ptF = ptP + dU2 * vtD ;
      }
     // due soluzioni valide ( punto esterno)
      else {
         ptI = ptP + dU1 * vtD ;
         ptF = ptP + dU2 * vtD ;
      }
     // Passo dal sistema intrinseco alla griglia
      if ( nGrid == 1) {
         swap( ptP0.y, ptP0.z) ;
         swap( ptP0.x, ptP0.z) ;
         swap( vtV0.y, vtV0.z) ;
         swap( vtV0.x, vtV0.z) ;
         swap( ptI.y, ptI.z) ;
         swap( ptI.x, ptI.z) ;
         swap( ptF.y, ptF.z) ;
         swap( ptF.x, ptF.z) ;
      }
      else if ( nGrid == 2) {
         swap( ptP0.x, ptP0.z) ;
         swap( ptP0.y, ptP0.z) ;
         swap( vtV0.x, vtV0.z) ;
         swap( vtV0.y, vtV0.z) ;
         swap( ptI.x, ptI.z) ;
         swap( ptI.y, ptI.z) ;
         swap( ptF.x, ptF.z) ;
         swap( ptF.y, ptF.z) ;
      }

     // Inizializzo distanze
      dInLen[nGrid] = INFINITO ;
      dOutLen[nGrid] = - INFINITO ;

     // Determino asse di spostamento maggiore
      bool bOnX = ( abs( ptF.y - ptI.y) <= abs( ptF.x - ptI.x)) ;

     // Movimento crescente su asse di movimento principale
      if ( ( bOnX  &&  ptF.x < ptI.x)  ||  ( ! bOnX  &&  ptF.y < ptI.y) )
         swap( ptI, ptF) ;

     // Pendenza
      double dDeltaX = ( ptF.x - ptI.x) ;
      double dDeltaY = ( ptF.y - ptI.y) ;
      double dM = 0 ;
      if ( bOnX  &&  abs( dDeltaX) > EPS_SMALL)
         dM = dDeltaY / dDeltaX ;
      else if ( ! bOnX  &&  abs( dDeltaY) > EPS_SMALL)
         dM = dDeltaX / dDeltaY ;

     // Determinazione degli indici i j dei punti ptI e ptF
      int nIi = Clamp( int( floor( ptI.x / m_dStep)), 0, m_nNx[nGrid] - 1) ;
      int nIj = Clamp( int( floor( ptI.y / m_dStep)), 0, m_nNy[nGrid] - 1) ;
      int nFi = Clamp( int( floor( ptF.x / m_dStep)), 0, m_nNx[nGrid] - 1) ;
      int nFj = Clamp( int( floor( ptF.y / m_dStep)), 0, m_nNy[nGrid] - 1) ;

     // mi muovo
      int i = nIi ; int j = nIj ;
      while ( ( bOnX  &&  i <= nFi)  ||  ( ! bOnX  &&  j <= nFj)) {

        // Eseguo controllo
         double dV1, dV2 ;
         if ( IntersRayDexel( ptP0, vtV0, nGrid, i, j, dV1, dV2)) {
            dInLen[nGrid] = min( dInLen[nGrid], dV1) ;
            dOutLen[nGrid] = max( dOutLen[nGrid], dV2) ;
         }

        // Calcolo spostamento (a destra o sopra)
         if ( bOnX) {
            double dMoveX = ( ( i + 1) * m_dStep - ptI.x) ;
            double dMoveY = dMoveX * dM ;
            double dY = ptI.y + dMoveY ;
            int nNewJ = Clamp( int( floor( dY / m_dStep)), 0, m_nNy[nGrid] - 1) ;
            if ( nNewJ != j)
               j = nNewJ ;
            else
               ++ i ;
         }
         else {
            double dMoveY = ( ( j + 1) * m_dStep - ptI.y) ;
            double dMoveX = dMoveY * dM ;
            double dX = ptI.x + dMoveX ;
            int nNewI = Clamp( int( floor( dX / m_dStep)), 0, m_nNx[nGrid] - 1) ;
            if ( nNewI != i)
               i = nNewI ;
            else
               ++ j ;
         }
      }

     // Se non abbiamo incontrato materiale
      if ( dInLen[nGrid] > dOutLen[nGrid] - EPS_SMALL) {
         dInLen[nGrid] = -2 ;
         dOutLen[nGrid] = -2 ;
      }
     // Se parto dall'interno
      else if ( dInLen[nGrid] < - EPS_SMALL)
         dInLen[nGrid] = -1 ;
   }

   if ( m_nMapNum == 1) {
      dInLength = dInLen[0] ;
      dOutLength = dOutLen[0] ;
      return true ;
   }

   if ( abs( dInLen[0] + 2) < EPS_SMALL  &&
        abs( dInLen[1] + 2) < EPS_SMALL  &&
        abs( dInLen[2] + 2) < EPS_SMALL) {
      dInLength = -2 ;
      dOutLength = -2 ;
      return true ;
   }
   else {
     dInLength = INFINITO ;
     dOutLength = - INFINITO ;
      for ( int nGr = 0 ; nGr < m_nMapNum ; ++ nGr) {
         if ( dInLen[nGr] > -2  &&  dInLen[nGr] < dInLength)
            dInLength = dInLen[nGr] ;
         if ( dOutLen[nGr] > -2  &&  dOutLen[nGr] > dOutLength) 
            dOutLength = dOutLen[nGr] ;
      }
   }

   return true ;
}

//----------------------------------------------------------------------------
// Calcola la profondità del materiale usando i triangoli della rappresentazione grafica.
// Punto e versore devono essere espressi nel sistema di riferimento 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& ptP, const Vector3d& vtD, double& dInLength, double& dOutLength) const 
{
   ILZIVECTOR vIntersInfo ;
   if ( ! GetLineIntersection( ptP, vtD, vIntersInfo))
      return false ;

   if ( vIntersInfo.empty()) {
      dInLength = -2. ;
      dOutLength = -2. ;
      return true ;
   }
   
   // Inizializzo le distanze di ingresso e uscita:
   // dInLength diminuisce, dOutLength aumenta.
   dInLength = INFINITO ;
   dOutLength = -INFINITO ;

   int nFirstPosN = 0 ;
   for ( int nN = 0; nN < int( vIntersInfo.size()) ; ++ nN) {
      if ( vIntersInfo[nN].dU > - EPS_SMALL) {
         nFirstPosN = nN ;
         break ;
      }
   }

   if ( nFirstPosN == int( vIntersInfo.size())) {
      dInLength = -2 ;
      dOutLength = -2 ;
      return true ;
   }

   if ( nFirstPosN > 0) {
      if ( vIntersInfo[nFirstPosN - 1].trTria.GetN() * vtD < EPS_ZERO)
         dInLength = -1 ;
   }
   else if ( nFirstPosN == 0) {
      if ( vIntersInfo[nFirstPosN].trTria.GetN() * vtD > EPS_ZERO)
         dInLength = -1 ;
      else {
         Point3d ptPi = ptP ;
         ptPi.ToLoc( m_MapFrame) ;
         int nVoxIJK[3] ;
         GetBlockIJKFromN( vIntersInfo[0].nVox, nVoxIJK) ;
         if ( GetPointVoxel( ptPi, nVoxIJK[0], nVoxIJK[1], nVoxIJK[2])) {
            int nCubeType = CalcIndex( nVoxIJK[0], nVoxIJK[1], nVoxIJK[2]) ;
            if ( nCubeType == 255)
               dInLength = -1 ;
         }
      }
   }

   for ( int nN = nFirstPosN ; nN < int( vIntersInfo.size()) ; ++ nN) {
      if ( vIntersInfo[nN].trTria.GetN() * vtD > - EPS_ZERO) {
         if ( ( vIntersInfo[nN].nILTT == ILTT_SEGM  ||  vIntersInfo[nN].nILTT == ILTT_SEGM_ON_EDGE)  && 
              dOutLength < vIntersInfo[nN].dU2)
            dOutLength = vIntersInfo[nN].dU2 ;
         else if ( dOutLength < vIntersInfo[nN].dU)
            dOutLength = vIntersInfo[nN].dU ;
      }
      if ( vIntersInfo[nN].trTria.GetN() * vtD < EPS_ZERO  &&
           dInLength > vIntersInfo[nN].dU)
         dInLength = vIntersInfo[nN].dU ;
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::AvoidSimpleBox( const Frame3d& frBox, const Vector3d& vtDiag, bool bPrecise) const
{
  // BBox
   BBox3d b3BoxL( ORIG, ORIG + vtDiag) ;

  // Porto il box nel riferimento intrinseco dello Zmap
   Frame3d frBoxInt = GetToLoc( frBox, m_MapFrame) ;
   BBox3d b3Box = GetToGlob( b3BoxL, frBoxInt) ;

  // BBox dello Zmap nel suo riferimento intrinseco
   BBox3d b3Zmap( ORIG, Point3d( m_nNx[0] * m_dStep, m_nNy[0] * m_dStep, m_dMaxZ[0])) ;

  // Se non interferiscono, posso uscire
   if ( ! b3Zmap.Overlaps( b3Box)  ||  ! b3Zmap.Overlaps( frBoxInt, b3BoxL))
      return true ;
   BBox3d b3Int ;
   if ( ! b3Zmap.FindIntersection( b3Box, b3Int))
      return true ;

  // Se verifico solo prima mappa
   if ( ! bPrecise) {

     // Limiti su indici
      int nStI = Clamp( int( b3Int.GetMin().x / m_dStep), 0, m_nNx[0] - 1) ;
      int nEnI = Clamp( int( b3Int.GetMax().x / m_dStep), 0, m_nNx[0] - 1) ;
      int nStJ = Clamp( int( b3Int.GetMin().y / m_dStep), 0, m_nNy[0] - 1) ;
      int nEnJ = Clamp( int( b3Int.GetMax().y / m_dStep), 0, m_nNy[0] - 1) ;

     // Vettore direzione dei dexel nel riferimento del Box
      Vector3d vtK = Z_AX ;  vtK.ToLoc( frBoxInt) ;

     // Riferimento intrinseco dei dexel nel riferimento del box
      Point3d ptO = ORIG ;   ptO.ToLoc( frBoxInt) ;
      Vector3d vtX = X_AX ;  vtX.ToLoc( frBoxInt) ;
      Vector3d vtY = Y_AX ;  vtY.ToLoc( frBoxInt) ;

     // Ciclo di intersezione dei dexel con il BBox
      for ( int i = nStI ; i <= nEnI ; ++ i) {
         for ( int j = nStJ ; j <= nEnJ ; ++ j) {
            int nPos = j * m_nNx[0] + i ;
            int nSize = int( m_Values[0][nPos].size()) ;
            if ( nSize == 0)
               continue ;
            for ( int k = 0 ; k < 5 ; ++ k) {
               Point3d ptT = ptO + ( i + 0.5) * m_dStep * vtX + ( j + 0.5) * m_dStep * vtY ;
               if ( k == 0)
                  ;
               else if ( k == 1)
                  ptT += - 0.4 * m_dStep * vtX - 0.4 * m_dStep * vtY ;
               else if ( k == 2)
                  ptT += + 0.4 * m_dStep * vtX - 0.4 * m_dStep * vtY ;
               else if ( k == 3)
                  ptT += + 0.4 * m_dStep * vtX + 0.4 * m_dStep * vtY ;
               else if ( k == 4)
                  ptT += - 0.4 * m_dStep * vtX + 0.4 * m_dStep * vtY ;
               double dZmin, dZmax ;
               if ( IntersLineBox( ptT, vtK, ORIG, ORIG + vtDiag, dZmin, dZmax)) {
                  for ( int nIndex = 0 ; nIndex < nSize ; nIndex += 1) {
                     if ( dZmax > m_Values[0][nPos][nIndex].dMin - EPS_SMALL  &&
                          dZmin < m_Values[0][nPos][nIndex].dMax + EPS_SMALL)
                        return false ;
                  }
               }
            }
         }
      }
   }

  // altrimenti verifico con tutte e tre le mappe
   else {
     // Ciclo di intersezione dei dexel con il BBox
      for ( int nMap = 0 ; nMap < m_nMapNum ; ++ nMap) {
         Point3d ptO = ORIG ;
         Vector3d vtX = X_AX ;
         Vector3d vtY = Y_AX ;
         Vector3d vtK = Z_AX ;
        // Estremi del box
         Point3d ptBoxInf = b3Int.GetMin() ;
         Point3d ptBoxSup = b3Int.GetMax() ;
         if ( nMap == 1) {
            swap( ptBoxInf.x, ptBoxInf.z) ;
            swap( ptBoxInf.x, ptBoxInf.y) ;
            swap( ptBoxSup.x, ptBoxSup.z) ;
            swap( ptBoxSup.x, ptBoxSup.y) ;
            vtX = Y_AX ;
            vtY = Z_AX ;
            vtK = X_AX ;
         }
         else if ( nMap == 2) {
            swap( ptBoxInf.y, ptBoxInf.z) ;
            swap( ptBoxInf.x, ptBoxInf.y) ;
            swap( ptBoxSup.y, ptBoxSup.z) ;
            swap( ptBoxSup.x, ptBoxSup.y) ;
            vtX = Z_AX ;
            vtY = X_AX ;
            vtK = Y_AX ;
         }
        // Passo da sistema griglia a sistema BBox
         ptO.ToLoc( frBoxInt) ;
         vtX.ToLoc( frBoxInt) ;
         vtY.ToLoc( frBoxInt) ;
         vtK.ToLoc( frBoxInt) ;
         // Limiti su indici
         int nStI = Clamp( int( ptBoxInf.x / m_dStep), 0, m_nNx[nMap] - 1) ;
         int nEnI = Clamp( int( ptBoxSup.x / m_dStep), 0, m_nNx[nMap] - 1) ;
         int nStJ = Clamp( int( ptBoxInf.y / m_dStep), 0, m_nNy[nMap] - 1) ;
         int nEnJ = Clamp( int( ptBoxSup.y / m_dStep), 0, m_nNy[nMap] - 1) ;
         // Ciclo sui dexel.
         for ( int i = nStI ; i <= nEnI ; ++ i) {
            for ( int j = nStJ ; j <= nEnJ ; ++ j) {
               int nPos = j * m_nNx[nMap] + i;
               int nSize = int( m_Values[nMap][nPos].size()) ;
               if ( nSize == 0)
                  continue ;
               Point3d ptT = ptO + ( i + 0.5) * m_dStep * vtX + ( j + 0.5) * m_dStep * vtY ;
               double dMinU, dMaxU ;
               // La retta associata al dexel interseca il box.
               if ( IntersLineBox( ptT, vtK, ORIG, ORIG + vtDiag, dMinU, dMaxU)) {
                  // Ciclo sui segmenti del dexel
                  for ( int nIndex = 0 ; nIndex < nSize ; nIndex += 1) {
                     // Se il segmento è interno all'intervallo d'intersezione, ho finito.
                     if ( dMaxU > m_Values[nMap][nPos][nIndex].dMin - EPS_SMALL  &&
                          dMinU < m_Values[nMap][nPos][nIndex].dMax + EPS_SMALL)
                        return false ;
                  }
               }
            }
         }
      }
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::AvoidBox( const Frame3d& frBox, const Vector3d& vtDiag, double dSafeDist, bool bPrecise) const
{
  // Se distanza di sicurezza nulla
   if ( dSafeDist < EPS_SMALL)
      return AvoidSimpleBox( frBox, vtDiag, bPrecise) ;

  // Verifica preliminare con box esteso
   Frame3d frEst = frBox ; frEst.Translate( -dSafeDist * ( frBox.VersX() + frBox.VersY() + frBox.VersZ())) ;
   if ( AvoidSimpleBox( frEst, vtDiag + 2 * Vector3d( dSafeDist, dSafeDist, dSafeDist), bPrecise))
      return true ;

  // Tre box aumentati con distanza di sicurezza in un sola dimensione
   Frame3d frTmp = frBox ; frTmp.Translate( -dSafeDist * frBox.VersX()) ;
   if ( ! AvoidSimpleBox( frTmp, vtDiag + 2 * dSafeDist * X_AX, bPrecise))
      return false ;
   frTmp = frBox ; frTmp.Translate( -dSafeDist * frBox.VersY()) ;
   if ( ! AvoidSimpleBox( frTmp, vtDiag + 2 * dSafeDist * Y_AX, bPrecise))
      return false ;
   frTmp = frBox ; frTmp.Translate( -dSafeDist * frBox.VersZ()) ;
   if ( ! AvoidSimpleBox( frTmp, vtDiag + 2 * dSafeDist * Z_AX, bPrecise))
      return false ;

  // Vertici del box
   Point3d ptVert0 = frBox.Orig() ;
   Point3d ptVert1 = ptVert0 + vtDiag.x * frBox.VersX() ;
   Point3d ptVert2 = ptVert0 + vtDiag.y * frBox.VersY() ;
   Point3d ptVert3 = ptVert1 + vtDiag.y * frBox.VersY() ;
   Point3d ptVert4 = ptVert0 + vtDiag.z * frBox.VersZ() ;
   Point3d ptVert5 = ptVert1 + vtDiag.z * frBox.VersZ() ;
   Point3d ptVert6 = ptVert2 + vtDiag.z * frBox.VersZ() ;
   Point3d ptVert7 = ptVert3 + vtDiag.z * frBox.VersZ() ;

  // Sfere centrate negli otto vertici
   if ( ! AvoidSimpleSphere( ptVert0, dSafeDist, bPrecise))
      return false ;
   if ( ! AvoidSimpleSphere( ptVert1, dSafeDist, bPrecise))
      return false ;
   if ( ! AvoidSimpleSphere( ptVert2, dSafeDist, bPrecise))
      return false ;
   if ( ! AvoidSimpleSphere( ptVert3, dSafeDist, bPrecise))
      return false ;
   if ( ! AvoidSimpleSphere( ptVert4, dSafeDist, bPrecise))
      return false ;
   if ( ! AvoidSimpleSphere( ptVert5, dSafeDist, bPrecise))
      return false ;
   if ( ! AvoidSimpleSphere( ptVert6, dSafeDist, bPrecise))
      return false ;
   if ( ! AvoidSimpleSphere( ptVert7, dSafeDist, bPrecise))
      return false ;

  // Cilindri centrati sui dodici spigoli
   frTmp.Set( ptVert0, frBox.VersX()) ;
   if ( ! AvoidSimpleCylinder( frTmp, dSafeDist, vtDiag.x, bPrecise))
      return false ;
   frTmp.Set( ptVert2, frBox.VersX()) ;
   if ( ! AvoidSimpleCylinder( frTmp, dSafeDist, vtDiag.x, bPrecise))
      return false ;
   frTmp.Set( ptVert0, frBox.VersY()) ;
   if ( ! AvoidSimpleCylinder( frTmp, dSafeDist, vtDiag.y, bPrecise))
      return false ;
   frTmp.Set( ptVert1, frBox.VersY()) ;
   if ( ! AvoidSimpleCylinder( frTmp, dSafeDist, vtDiag.y, bPrecise))
      return false ;
   frTmp.Set( ptVert4, frBox.VersX()) ;
   if ( ! AvoidSimpleCylinder( frTmp, dSafeDist, vtDiag.x, bPrecise))
      return false ;
   frTmp.Set( ptVert6, frBox.VersX()) ;
   if ( ! AvoidSimpleCylinder( frTmp, dSafeDist, vtDiag.x, bPrecise))
      return false ;
   frTmp.Set( ptVert4, frBox.VersY()) ;
   if ( ! AvoidSimpleCylinder( frTmp, dSafeDist, vtDiag.y, bPrecise))
      return false ;
   frTmp.Set( ptVert5, frBox.VersY()) ;
   if ( ! AvoidSimpleCylinder( frTmp, dSafeDist, vtDiag.y, bPrecise))
      return false ;
   frTmp.Set( ptVert0, frBox.VersZ()) ;
   if ( ! AvoidSimpleCylinder( frTmp, dSafeDist, vtDiag.z, bPrecise))
      return false ;
   frTmp.Set( ptVert1, frBox.VersZ()) ;
   if ( ! AvoidSimpleCylinder( frTmp, dSafeDist, vtDiag.z, bPrecise))
      return false ;
   frTmp.Set( ptVert2, frBox.VersZ()) ;
   if ( ! AvoidSimpleCylinder( frTmp, dSafeDist, vtDiag.z, bPrecise))
      return false ;
   frTmp.Set( ptVert3, frBox.VersZ()) ;
   if ( ! AvoidSimpleCylinder( frTmp, dSafeDist, vtDiag.z, bPrecise))
      return false ;

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::AvoidSimpleSphere( const Point3d& ptCenter, double dRad, bool bPrecise) const
{
  // Porto la sfera nel riferimento intrinseco dello Zmap
   Point3d ptC = ptCenter ;
   ptC.ToLoc( m_MapFrame) ;

  // BBox della sfera
   BBox3d b3Box( ptC) ;
   b3Box.Expand( dRad) ;

  // BBox dello Zmap nel suo riferimento intrinseco
   BBox3d b3Zmap( ORIG, Point3d( m_nNx[0] * m_dStep, m_nNy[0] * m_dStep, m_dMaxZ[0])) ;

  // Se non interferiscono, posso uscire
   BBox3d b3Int ;
   if ( ! b3Zmap.FindIntersection( b3Box, b3Int))
      return true ;

  // Se verifico solo prima mappa
   if ( ! bPrecise) {

     // Limiti su indici
      int nStI = Clamp( int( b3Int.GetMin().x / m_dStep), 0, m_nNx[0] - 1) ;
      int nEnI = Clamp( int( b3Int.GetMax().x / m_dStep), 0, m_nNx[0] - 1) ;
      int nStJ = Clamp( int( b3Int.GetMin().y / m_dStep), 0, m_nNy[0] - 1) ;
      int nEnJ = Clamp( int( b3Int.GetMax().y / m_dStep), 0, m_nNy[0] - 1) ;

     // Ciclo di intersezione dei dexel con la sfera (nel riferimento intrinseco)
      for ( int i = nStI ; i <= nEnI ; ++ i) {
         for ( int j = nStJ ; j <= nEnJ ; ++ j) {
            int nPos = j * m_nNx[0] + i ;
            int nSize = int( m_Values[0][nPos].size()) ;
            if ( nSize == 0)
               continue ;
            if ( m_Values[0][nPos][nSize-1].dMax < b3Int.GetMin().z  ||  m_Values[0][nPos][0].dMin > b3Int.GetMax().z)
               continue ;
            for ( int k = 0 ; k < 5 ; ++ k) {
               Point3d ptT = ORIG + ( i + 0.5) * m_dStep * X_AX + ( j + 0.5) * m_dStep * Y_AX ;
               switch ( k) {
               case 0 : break ;
               case 1 : ptT += -0.4 * m_dStep * X_AX - 0.4 * m_dStep * Y_AX ; break ;
               case 2 : ptT += +0.4 * m_dStep * X_AX - 0.4 * m_dStep * Y_AX ; break ;
               case 3 : ptT += +0.4 * m_dStep * X_AX + 0.4 * m_dStep * Y_AX ; break ;
               case 4 : ptT += -0.4 * m_dStep * X_AX + 0.4 * m_dStep * Y_AX ; break ;
               }
               Point3d ptI1, ptI2 ;
               if ( ::IntersLineSphere( ptT, Z_AX, ptC, dRad, ptI1, ptI2) != ILST_NO) {
                  double dZmin = min( ptI1.z, ptI2.z) ;
                  double dZmax = max( ptI1.z, ptI2.z) ;
                  for ( int nIndex = 0 ; nIndex < nSize ; nIndex += 1) {
                     if ( dZmax > m_Values[0][nPos][nIndex].dMin - EPS_SMALL  &&
                          dZmin < m_Values[0][nPos][nIndex].dMax + EPS_SMALL)
                        return false ;
                  }
               }
            }
         }
      }
   }

  // altrimenti verifico con tutte e tre le mappe
   else {
     // Ciclo di intersezione dei dexel con la sfera (nel riferimento intrinseco)
      for ( int nMap = 0 ; nMap < m_nMapNum ; ++ nMap) {
         Vector3d vtX = X_AX ;
         Vector3d vtY = Y_AX ;
         Vector3d vtK = Z_AX ;
        // Estremi del box
         Point3d ptBoxInf = b3Int.GetMin() ;
         Point3d ptBoxSup = b3Int.GetMax() ;
         if ( nMap == 1) {
            swap( ptBoxInf.x, ptBoxInf.z) ;
            swap( ptBoxInf.x, ptBoxInf.y) ;
            swap( ptBoxSup.x, ptBoxSup.z) ;
            swap( ptBoxSup.x, ptBoxSup.y) ;
            vtX = Y_AX ;
            vtY = Z_AX ;
            vtK = X_AX ;
         }
         else if ( nMap == 2) {
            swap( ptBoxInf.y, ptBoxInf.z) ;
            swap( ptBoxInf.x, ptBoxInf.y) ;
            swap( ptBoxSup.y, ptBoxSup.z) ;
            swap( ptBoxSup.x, ptBoxSup.y) ;
            vtX = Z_AX ;
            vtY = X_AX ;
            vtK = Y_AX ;
         }
        // Limiti su indici
         int nStI = Clamp( int( ptBoxInf.x / m_dStep), 0, m_nNx[nMap] - 1) ;
         int nEnI = Clamp( int( ptBoxSup.x / m_dStep), 0, m_nNx[nMap] - 1) ;
         int nStJ = Clamp( int( ptBoxInf.y / m_dStep), 0, m_nNy[nMap] - 1) ;
         int nEnJ = Clamp( int( ptBoxSup.y / m_dStep), 0, m_nNy[nMap] - 1) ;
        // Ciclo sui dexel.
         for ( int i = nStI ; i <= nEnI ; ++ i) {
            for ( int j = nStJ ; j <= nEnJ ; ++ j) {
               int nPos = j * m_nNx[nMap] + i ;
               int nSize = int( m_Values[nMap][nPos].size()) ;
               if ( nSize == 0)
                  continue ;
               if ( m_Values[nMap][nPos][nSize-1].dMax < b3Int.GetMin().z  ||  m_Values[nMap][nPos][0].dMin > b3Int.GetMax().z)
                  continue ;
               Point3d ptT = ORIG + ( i + 0.5) * m_dStep * vtX + ( j + 0.5) * m_dStep * vtY ;
               Point3d ptI1, ptI2 ;
              // La linea del dexel interseca la sfera.
               if ( ::IntersLineSphere( ptT, vtK, ptC, dRad, ptI1, ptI2) != ILST_NO) {
                  double dMinU, dMaxU ;
                  if ( nMap == 0) {
                     dMinU = min( ptI1.z, ptI2.z) ;
                     dMaxU = max( ptI1.z, ptI2.z) ;
                  }
                  else if ( nMap == 1) {
                     dMinU = min( ptI1.x, ptI2.x) ;
                     dMaxU = max( ptI1.x, ptI2.x) ;
                  }
                  else {
                     dMinU = min( ptI1.y, ptI2.y) ;
                     dMaxU = max( ptI1.y, ptI2.y) ;
                  }
                // Ciclo sui segmenti del dexel.
                  for ( int nIndex = 0 ; nIndex < nSize ; nIndex += 1) {
                    // Se il segmento è interno all'intervallo d'intersezione, ho finito.
                     if ( dMaxU > m_Values[nMap][nPos][nIndex].dMin - EPS_SMALL  &&
                          dMinU < m_Values[nMap][nPos][nIndex].dMax + EPS_SMALL)
                        return false ;
                  }
               }
            }
         }
      }
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::AvoidSphere( const Point3d& ptCenter, double dRad, double dSafeDist, bool bPrecise) const
{
   return AvoidSimpleSphere( ptCenter, dRad + max( 0., dSafeDist), bPrecise) ;
}

//----------------------------------------------------------------------------
bool
VolZmap::AvoidSimpleCylinder( const Frame3d& frCyl, double dR, double dH, bool bPrecise) const
{
  // BBox del cilindro in locale
   BBox3d b3CylL( Point3d( -dR, -dR, 0), Point3d( dR, dR, dH)) ;

  // BBox del cilindro nel riferimento intrinseco dello Zmap
   Frame3d frCylInt = GetToLoc( frCyl, m_MapFrame) ;
   BBox3d b3CylI = GetToGlob( b3CylL, frCylInt) ;

  // BBox dello Zmap nel suo riferimento intrinseco
   BBox3d b3Zmap( ORIG, Point3d( m_nNx[0] * m_dStep, m_nNy[0] * m_dStep, m_dMaxZ[0])) ;

  // Se non interferiscono, posso uscire
   if ( ! b3Zmap.Overlaps( b3CylI)  ||  ! b3Zmap.Overlaps( frCylInt, b3CylL))
      return true ;

  // BBox del cilindro ottimizzato nel riferimento intrinseco dello Zmap
   Point3d ptMyCen = frCylInt.Orig() ;
   Vector3d vtMyAx = frCylInt.VersZ() ;
   BBox3d b3Cyl( ptMyCen) ;
   b3Cyl.Add( ptMyCen + vtMyAx * dH) ;
   if ( vtMyAx.IsX())
      b3Cyl.Expand( 0, dR, dR) ;
   else if ( vtMyAx.IsY())
      b3Cyl.Expand( dR, 0, dR) ;
   else if ( vtMyAx.IsZ())
      b3Cyl.Expand( dR, dR, 0) ;
   else {
      double dExpandX = dR * sqrt( 1 - vtMyAx.x * vtMyAx.x) ;
      double dExpandY = dR * sqrt( 1 - vtMyAx.y * vtMyAx.y) ;
      double dExpandZ = dR * sqrt( 1 - vtMyAx.z * vtMyAx.z) ;
      b3Cyl.Expand( dExpandX, dExpandY, dExpandZ) ;
   }

  // Se non interferiscono, posso uscire
   BBox3d b3Int ;
   if ( ! b3Zmap.FindIntersection( b3Cyl, b3Int))
      return true ;

  // Se verifico solo prima mappa
   if ( ! bPrecise) {

     // Limiti su indici
      int nStI = Clamp( int( b3Int.GetMin().x / m_dStep), 0, m_nNx[0] - 1) ;
      int nEnI = Clamp( int( b3Int.GetMax().x / m_dStep), 0, m_nNx[0] - 1) ;
      int nStJ = Clamp( int( b3Int.GetMin().y / m_dStep), 0, m_nNy[0] - 1) ;
      int nEnJ = Clamp( int( b3Int.GetMax().y / m_dStep), 0, m_nNy[0] - 1) ;

     // Ciclo di intersezione dei dexel con il cilindro (nel riferimento intrinseco)
      for ( int i = nStI ; i <= nEnI ; ++ i) {
         for ( int j = nStJ ; j <= nEnJ ; ++ j) {
            int nPos = j * m_nNx[0] + i ;
            int nSize = int( m_Values[0][nPos].size()) ;
            if ( nSize == 0)
               continue ;
            if ( m_Values[0][nPos][nSize-1].dMax < b3Int.GetMin().z  ||  m_Values[0][nPos][0].dMin > b3Int.GetMax().z)
               continue ;
            for ( int k = 0 ; k < 5 ; ++ k) {
               Point3d ptT = ORIG + ( i + 0.5) * m_dStep * X_AX + ( j + 0.5) * m_dStep * Y_AX ;
               switch ( k) {
               case 0 : break ;
               case 1 : ptT += -0.4 * m_dStep * X_AX - 0.4 * m_dStep * Y_AX ; break ;
               case 2 : ptT += +0.4 * m_dStep * X_AX - 0.4 * m_dStep * Y_AX ; break ;
               case 3 : ptT += +0.4 * m_dStep * X_AX + 0.4 * m_dStep * Y_AX ; break ;
               case 4 : ptT += -0.4 * m_dStep * X_AX + 0.4 * m_dStep * Y_AX ; break ;
               }
               Point3d ptI1, ptI2 ;
               Vector3d vtN1, vtN2 ;
               if ( IntersLineCylinder( ptT, Z_AX, frCylInt, dH, dR, true, true, ptI1, vtN1, ptI2, vtN2)) {
                  double dZmin = min( ptI1.z, ptI2.z) ;
                  double dZmax = max( ptI1.z, ptI2.z) ;
                  for ( int nIndex = 0 ; nIndex < nSize ; nIndex += 1) {
                     if ( dZmax > m_Values[0][nPos][nIndex].dMin - EPS_SMALL  &&
                          dZmin < m_Values[0][nPos][nIndex].dMax + EPS_SMALL)
                        return false ;
                  }
               }
            }
         }
      }
   }

  // altrimenti verifico con tutte e tre le mappe
   else {
     // Ciclo di intersezione dei dexel con il cilindro (nel riferimento intrinseco)
      for ( int nMap = 0 ; nMap < m_nMapNum ; ++ nMap) {
         Vector3d vtX = X_AX ;
         Vector3d vtY = Y_AX ;
         Vector3d vtK = Z_AX ;
        // Estremi del box
         Point3d ptBoxInf = b3Int.GetMin() ;
         Point3d ptBoxSup = b3Int.GetMax() ;
         if ( nMap == 1) {
            swap( ptBoxInf.x, ptBoxInf.z) ;
            swap( ptBoxInf.x, ptBoxInf.y) ;
            swap( ptBoxSup.x, ptBoxSup.z) ;
            swap( ptBoxSup.x, ptBoxSup.y) ;
            vtX = Y_AX ;
            vtY = Z_AX ;
            vtK = X_AX ;
         }
         else if ( nMap == 2) {
            swap( ptBoxInf.y, ptBoxInf.z) ;
            swap( ptBoxInf.x, ptBoxInf.y) ;
            swap( ptBoxSup.y, ptBoxSup.z) ;
            swap( ptBoxSup.x, ptBoxSup.y) ;
            vtX = Z_AX ;
            vtY = X_AX ;
            vtK = Y_AX ;
         }
        // Limiti su indici
         int nStI = Clamp( int( ptBoxInf.x / m_dStep), 0, m_nNx[nMap] - 1) ;
         int nEnI = Clamp( int( ptBoxSup.x / m_dStep), 0, m_nNx[nMap] - 1) ;
         int nStJ = Clamp( int( ptBoxInf.y / m_dStep), 0, m_nNy[nMap] - 1) ;
         int nEnJ = Clamp( int( ptBoxSup.y / m_dStep), 0, m_nNy[nMap] - 1) ;
        // Ciclo sui dexel.
         for ( int i = nStI ; i <= nEnI ; ++ i) {
            for ( int j = nStJ ; j <= nEnJ ; ++ j) {
               int nPos = j * m_nNx[nMap] + i ;
               int nSize = int( m_Values[nMap][nPos].size()) ;
               if ( nSize == 0)
                  continue ;
               if ( m_Values[nMap][nPos][nSize-1].dMax < b3Int.GetMin().z  ||  m_Values[nMap][nPos][0].dMin > b3Int.GetMax().z)
                  continue ;
               Point3d ptT = ORIG + ( i + 0.5) * m_dStep * vtX + ( j + 0.5) * m_dStep * vtY ;
               Point3d ptI1, ptI2 ;
               Vector3d vtN1, vtN2 ;
              // La linea del dexel interseca il cilindro.
               if ( IntersLineCylinder( ptT, vtK, frCylInt, dH, dR, true, true, ptI1, vtN1, ptI2, vtN2)) {
                  double dMinU, dMaxU ;
                  if ( nMap == 0) {
                     dMinU = min( ptI1.z, ptI2.z) ;
                     dMaxU = max( ptI1.z, ptI2.z) ;
                  }
                  else if ( nMap == 1) {
                     dMinU = min( ptI1.x, ptI2.x) ;
                     dMaxU = max( ptI1.x, ptI2.x) ;
                  }
                  else {
                     dMinU = min( ptI1.y, ptI2.y) ;
                     dMaxU = max( ptI1.y, ptI2.y) ;
                  }
                 // Ciclo sui segmenti del dexel. 
                  for ( int nIndex = 0 ; nIndex < nSize ; nIndex += 1) {
                    // Se il segmento è interno all'intervallo d'intersezione, ho finito.
                     if ( dMaxU > m_Values[nMap][nPos][nIndex].dMin - EPS_SMALL  &&
                          dMinU < m_Values[nMap][nPos][nIndex].dMax + EPS_SMALL)
                        return false ;
                  }
               }
            }
         }
      }
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::AvoidCylinder( const Frame3d& frCyl, double dR, double dH, double dSafeDist, bool bPrecise) const
{
  // Se altezza negativa, sposto riferimento da faccia sopra a quella sotto
   Frame3d frMyCyl = frCyl ;
   if ( dH < 0) {
      frMyCyl.Translate( dH * frMyCyl.VersZ()) ;
      dH = - dH ;
   }

  // Se distanza di sicurezza nulla
   if ( dSafeDist < EPS_SMALL)
      return AvoidSimpleCylinder( frMyCyl, dR, dH, bPrecise) ;

  // Verifica preliminare con cilindro esteso
   Frame3d frEst = frMyCyl ; frEst.Translate( -dSafeDist * frMyCyl.VersZ()) ;
   if ( AvoidSimpleCylinder( frEst, dR + dSafeDist, dH + 2 * dSafeDist, bPrecise))
      return true ;

  // Cilindro allargato
   if ( ! AvoidSimpleCylinder( frMyCyl, dR + dSafeDist, dH, bPrecise))
      return false ;
  // Cilindro allungato
   Frame3d frTmp = frMyCyl ; frTmp.Translate( - dSafeDist * frMyCyl.VersZ()) ;
   if ( ! AvoidSimpleCylinder( frTmp, dR, dH + 2 * dSafeDist, bPrecise))
      return false ;
  // Toro inferiore
   if ( ! AvoidSimpleTorus( frMyCyl, dR, dSafeDist, bPrecise))
      return false ;
  // Toro superiore
   frTmp = frMyCyl ; frTmp.Translate( dH * frMyCyl.VersZ()) ;
   if ( ! AvoidSimpleTorus( frTmp, dR, dSafeDist, bPrecise))
      return false ;

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::SingleMapDexelConeCollision( int nStI, int nEnI, int nStJ, int nEnJ, const Point3d& ptRefPoint, const Vector3d& vtRefAx,
                                      double dMinRad, double dMaxRad, double dHeight, double dMinBoxH, double dMaxBoxH) const
{
  // Ciclo di intersezione dei dexel con il cono (nel riferimento intrinseco)
   for ( int i = nStI ; i <= nEnI ; ++ i) {
      for ( int j = nStJ ; j <= nEnJ ; ++ j) {
         int nPos = j * m_nNx[0] + i ;
         int nSize = int( m_Values[0][nPos].size()) ;
         if ( nSize == 0)
            continue ;
         double dParMin = m_Values[0][nPos][0].dMin ;
         double dParMax = m_Values[0][nPos][nSize-1].dMax ;
         if ( dParMax < dMinBoxH  ||  dParMin > dMaxBoxH)
            continue ;
         for ( int k = 0 ; k < 5 ; ++ k) {
           // se richiesta interruzione, esco
            if ( m_bBreak)
               return false ;
           // Calcolo spillone
            Point3d ptLineSt = ORIG + ( i + 0.5) * m_dStep * X_AX + ( j + 0.5) * m_dStep * Y_AX ;
            switch ( k) {
            case 0 : break ;
            case 1 : ptLineSt += -0.4 * m_dStep * X_AX - 0.4 * m_dStep * Y_AX ; break ;
            case 2 : ptLineSt += +0.4 * m_dStep * X_AX - 0.4 * m_dStep * Y_AX ; break ;
            case 3 : ptLineSt += +0.4 * m_dStep * X_AX + 0.4 * m_dStep * Y_AX ; break ;
            case 4 : ptLineSt += -0.4 * m_dStep * X_AX + 0.4 * m_dStep * Y_AX ; break ;
            }
           // Cono proprio
            if ( dMinRad < EPS_SMALL) {
               Point3d ptSegSt = ptLineSt + dParMin * Z_AX ;
               double dU1, dU2 ;
               int nIntType = SegmentCone( ptSegSt, Z_AX, dParMax - dParMin, ptRefPoint, vtRefAx,
                                           dMaxRad, dHeight, dU1, dU2) ;
               if ( nIntType == LinCompCCIntersType::CC_ERROR_INT)
                  return true ;
               else if ( nIntType != LinCompCCIntersType::CC_NO_INTERS) {
                  for ( int nIndex = 0 ; nIndex < nSize ; nIndex += 1) {
                     if ( m_Values[0][nPos][nIndex].dMax >= dU1  &&  m_Values[0][nPos][nIndex].dMin <= dU2)
                        return true ;
                  }
               }
               nIntType = SegmentDisc( ptSegSt, Z_AX, dParMax - dParMin, ptRefPoint + dHeight * vtRefAx,
                                       vtRefAx, dMaxRad, dU1, dU2) ;
               if ( nIntType == LinCompDiscIntersType::D_ERROR_INT)
                  return true ;
               else if ( nIntType != LinCompDiscIntersType::D_NO_INTERS) {
                  for ( int nIndex = 0 ; nIndex < nSize ; nIndex += 1) {
                     if ( m_Values[0][nPos][nIndex].dMax >= dU1  &&  m_Values[0][nPos][nIndex].dMin <= dU2)
                        return true ;
                  }
               }
            }
           // Tronco di cono
            else {
               Point3d ptSegSt = ptLineSt + dParMin * Z_AX ;
               double dU1, dU2 ;
               int nIntType = SegmentConeFrustum( ptSegSt, Z_AX, dParMax - dParMin, ptRefPoint, vtRefAx,
                                                  dMinRad, dMaxRad, dHeight, dU1, dU2) ;
               if ( nIntType == LinCompCCIntersType::CC_ERROR_INT)
                  return true ;
               else if ( nIntType != LinCompCCIntersType::CC_NO_INTERS) {
                  for ( int nIndex = 0 ; nIndex < nSize ; nIndex += 1) {
                     if ( m_Values[0][nPos][nIndex].dMax >= dU1  &&  m_Values[0][nPos][nIndex].dMin <= dU2)
                        return true ;
                  }
               }
               nIntType = SegmentDisc( ptSegSt, Z_AX, dParMax - dParMin, ptRefPoint + dHeight * vtRefAx,
                                       vtRefAx, dMaxRad, dU1, dU2) ;
               if ( nIntType == LinCompDiscIntersType::D_ERROR_INT)
                  return true ;
               else if ( nIntType != LinCompDiscIntersType::D_NO_INTERS) {
                  for ( int nIndex = 0 ; nIndex < nSize ; nIndex += 1) {
                     if ( m_Values[0][nPos][nIndex].dMax >= dU1  &&  m_Values[0][nPos][nIndex].dMin <= dU2)
                        return true ;
                  }
               }
               nIntType = SegmentDisc( ptSegSt, Z_AX, dParMax - dParMin, ptRefPoint, vtRefAx, dMinRad, dU1, dU2) ;
               if ( nIntType == LinCompDiscIntersType::D_ERROR_INT)
                  return true ;
               else if ( nIntType != LinCompDiscIntersType::D_NO_INTERS) {
                  for ( int nIndex = 0 ; nIndex < nSize ; nIndex += 1) {
                     if ( m_Values[0][nPos][nIndex].dMax >= dU1  &&  m_Values[0][nPos][nIndex].dMin <= dU2)
                        return true ;
                  }
               }
            }
         }
      }
   }

   return false ;
}

//----------------------------------------------------------------------------
bool
VolZmap::AvoidSimpleConeFrustum( const Frame3d& frCone, double dMinRad, double dMaxRad, double dHeight, bool bPrecise) const
{
  // BBox del tronco di cono in locale
   BBox3d b3ConeL( Point3d( -dMaxRad, -dMaxRad, 0), Point3d( dMaxRad, dMaxRad, dHeight)) ;

  // BBox del tronco di cono nel riferimento intrinseco dello Zmap
   Frame3d frConeInt = GetToLoc( frCone, m_MapFrame) ;
   BBox3d b3ConeI = GetToGlob( b3ConeL, frConeInt) ;

  // BBox dello Zmap nel suo riferimento intrinseco
   BBox3d b3Zmap( ORIG, Point3d( m_nNx[0] * m_dStep, m_nNy[0] * m_dStep, m_dMaxZ[0])) ;

  // Se non interferiscono, posso uscire
   if ( ! b3Zmap.Overlaps( b3ConeI)  ||  ! b3Zmap.Overlaps( frConeInt, b3ConeL))
      return true ;

  // BBox del tronco di cono ottimizzato nel riferimento intrinseco dello Zmap 
   Point3d ptRefPoint = frConeInt.Orig() ;
   Vector3d vtRefAx = frConeInt.VersZ() ;
   BBox3d b3Cone( ptRefPoint) ;
   b3Cone.Add( ptRefPoint + vtRefAx * dHeight) ;
   if ( vtRefAx.IsX())
      b3Cone.Expand( 0, dMaxRad, dMaxRad) ;
   else if ( vtRefAx.IsY())
      b3Cone.Expand( dMaxRad, 0, dMaxRad) ;
   else if ( vtRefAx.IsZ())
      b3Cone.Expand( dMaxRad, dMaxRad, 0) ;
   else {
      double dExpandX = dMaxRad * sqrt( 1 - vtRefAx.x * vtRefAx.x) ;
      double dExpandY = dMaxRad * sqrt( 1 - vtRefAx.y * vtRefAx.y) ;
      double dExpandZ = dMaxRad * sqrt( 1 - vtRefAx.z * vtRefAx.z) ;
      b3Cone.Expand( dExpandX, dExpandY, dExpandZ) ;
   }

  // Se non interferiscono, posso uscire
   BBox3d b3Int ;
   if ( ! b3Zmap.FindIntersection( b3Cone, b3Int))
      return true ;

  // Uso solo la prima mappa
   if ( ! bPrecise) {
     // Limiti su indici
      int nStI = Clamp( int( b3Int.GetMin().x / m_dStep), 0, m_nNx[0] - 1) ;
      int nEnI = Clamp( int( b3Int.GetMax().x / m_dStep), 0, m_nNx[0] - 1) ;
      int nStJ = Clamp( int( b3Int.GetMin().y / m_dStep), 0, m_nNy[0] - 1) ;
      int nEnJ = Clamp( int( b3Int.GetMax().y / m_dStep), 0, m_nNy[0] - 1) ;
     // Limiti su Z
      double dZmin = b3Int.GetMin().z ;
      double dZmax = b3Int.GetMax().z ;
     // Numero massimo di thread
      int nThreadMax = max( 1, int( thread::hardware_concurrency()) - 1) ;
     // se un solo thread
      if ( nThreadMax == 1) {
         m_bBreak = false ;
         bool bCollision = SingleMapDexelConeCollision( nStI, nEnI, nStJ, nEnJ, ptRefPoint, vtRefAx,
                                                        dMinRad, dMaxRad, dHeight, dZmin, dZmax) ;
         return ( ! bCollision) ;
      }
     // altrimenti esecuzione in parallelo dei calcoli
      else {
         //string sOut = "I=" + ToString( nStI) + "," + ToString( nEnI) + " J=" + ToString( nStJ) + "," + ToString( nEnJ) ;
         //LOG_INFO( GetEGkLogger(), sOut.c_str())
         m_bBreak = false ;
        // lancio dei thread
         int nSpanI = ( nEnI - nStI + 1) / nThreadMax + 1 ;
         int nSpanJ = ( nEnJ - nStJ + 1) / nThreadMax + 1 ;
         bool bOnI = ( nSpanI >= nSpanJ) ;
         int nThreadTot = 0 ;
         vector< future<bool>> vRes( nThreadMax) ;
         for ( int nT = 0 ; nT < nThreadMax ; ++ nT) {
            int nMyStI = ( bOnI ? nStI + nT * nSpanI : nStI) ;
            int nMyEnI = ( bOnI ? min( nMyStI + nSpanI - 1, nEnI) : nEnI) ;
            int nMyStJ = ( bOnI ? nStJ : nStJ + nT * nSpanJ) ;
            int nMyEnJ = ( bOnI ? nEnJ : min( nMyStJ + nSpanJ - 1, nEnJ)) ;
            if ( nMyStI > nEnI  ||  nMyStJ > nEnJ)
               break ;
            //string sOut = "MyI=" + ToString( nMyStI) + "," + ToString( nMyEnI) + " MyJ=" + ToString( nMyStJ) + "," + ToString( nMyEnJ) ;
            //LOG_INFO( GetEGkLogger(), sOut.c_str())
            vRes[nT] = async( launch::async, &VolZmap::SingleMapDexelConeCollision, this,
                              nMyStI, nMyEnI, nMyStJ, nMyEnJ, cref( ptRefPoint), cref( vtRefAx),
                              dMinRad, dMaxRad, dHeight, dZmin, dZmax) ;
            ++ nThreadTot ;
         }
        // recupero i risultati dei thread alla loro terminazione
         bool bCollision = false ;
         int nTerminated = 0 ;
         while ( nTerminated < nThreadTot) {
            for ( int nT = 0 ; nT < nThreadTot ; ++ nT) {
              // Async terminato
               if ( vRes[nT].valid()  &&  vRes[nT].wait_for( chrono::nanoseconds{ 1}) == future_status::ready) {
                  ++ nTerminated ;
                 // Se c'è collisione ...
                  if ( vRes[nT].get()) {
                     bCollision = true ;
                     m_bBreak = true ;
                  }
               }
            }
         }
         return ( ! bCollision) ;
      }
   }

  // Uso tutte le mappe
   else {
     // Ciclo sulle mappe.
      for ( int nMap = 0 ; nMap < m_nMapNum ; ++ nMap) {
         Point3d ptInfIntBox = b3Int.GetMin() ;
         Point3d ptSupIntBox = b3Int.GetMax() ;
        // Dal sistema intrinseco al sistema griglia (per la prima griglia coincidono).
         if ( nMap == 1) {
            swap( ptInfIntBox.x, ptInfIntBox.z) ;
            swap( ptInfIntBox.x, ptInfIntBox.y) ;
            swap( ptSupIntBox.x, ptSupIntBox.z) ;
            swap( ptSupIntBox.x, ptSupIntBox.y) ;
         }
         else if ( nMap == 2) {
            swap( ptInfIntBox.y, ptInfIntBox.z) ;
            swap( ptInfIntBox.x, ptInfIntBox.y) ;
            swap( ptSupIntBox.y, ptSupIntBox.z) ;
            swap( ptSupIntBox.x, ptSupIntBox.y) ;
         }
        // Limiti su indici
         int nStI = Clamp( int( ptInfIntBox.x / m_dStep), 0, m_nNx[nMap] - 1) ;
         int nEnI = Clamp( int( ptSupIntBox.x / m_dStep), 0, m_nNx[nMap] - 1) ;
         int nStJ = Clamp( int( ptInfIntBox.y / m_dStep), 0, m_nNy[nMap] - 1) ;
         int nEnJ = Clamp( int( ptSupIntBox.y / m_dStep), 0, m_nNy[nMap] - 1) ;
         // Ciclo sui dexel.
         for ( int nDex = 0 ; nDex < int( m_Values[nMap].size()) ; ++ nDex) {
            int nDexSize = (int)m_Values[nMap][nDex].size() ;
            if ( nDexSize == 0  ||
                 m_Values[nMap][nDex][ nDexSize- 1].dMax < ptInfIntBox.z  ||
                 m_Values[nMap][nDex][0].dMin > ptSupIntBox.z)
               continue ;
           // Indici del dexel.
            int nI = nDex % m_nNx[nMap] ;
            int nJ = nDex / m_nNx[nMap] ;
           // Se fuori dalla regione ammissibile salto l'iterazione
            if ( nI < nStI  ||  nI > nEnI  ||  nJ < nStJ ||  nJ > nEnJ)
               continue ;
           // Posizione del dexel.
            double dX = ( nI + 0.5) * m_dStep ;
            double dY = ( nJ + 0.5) * m_dStep ;
            Point3d ptLineSt( dX, dY, 0.) ;
            Vector3d vtLineDir( 0., 0., 1.) ;
           // Dal sistema griglia al sistema intrinseco (per la prima griglia coincidono).
            if ( nMap == 1) {
               swap( ptLineSt.x, ptLineSt.y) ;
               swap( ptLineSt.x, ptLineSt.z) ;
               swap( vtLineDir.x, vtLineDir.y) ;
               swap( vtLineDir.x, vtLineDir.z) ;
            }
            else if ( nMap == 2) {
               swap( ptLineSt.x, ptLineSt.y) ;
               swap( ptLineSt.y, ptLineSt.z) ;
               swap( vtLineDir.x, vtLineDir.y) ;
               swap( vtLineDir.y, vtLineDir.z) ;
            }
           // Cono proprio
            if ( dMinRad < EPS_SMALL) {
               double dMinPar = m_Values[nMap][nDex][0].dMin ;
               double dMaxPar = m_Values[nMap][nDex][nDexSize - 1].dMax ;
               Point3d ptSegSt = ptLineSt + dMinPar * vtLineDir ;
               double dU1, dU2 ;
               int nIntType = SegmentCone( ptSegSt, vtLineDir, dMaxPar - dMinPar, ptRefPoint, vtRefAx,
                                           dMaxRad, dHeight, dU1, dU2) ;
               if ( nIntType == LinCompCCIntersType::CC_ERROR_INT)
                  return false ;
               else if ( nIntType != LinCompCCIntersType::CC_NO_INTERS) {
                  for ( int nIndex = 0 ; nIndex < nDexSize ; nIndex += 1) {
                     if ( m_Values[nMap][nDex][nIndex].dMax >= dU1  &&  m_Values[nMap][nDex][nIndex].dMin <= dU2)
                        return false ;
                  }
               }
               nIntType = SegmentDisc( ptSegSt, vtLineDir, dMaxPar - dMinPar, ptRefPoint + dHeight * vtRefAx,
                                       vtRefAx, dMaxRad, dU1, dU2) ;
               if ( nIntType == LinCompDiscIntersType::D_ERROR_INT)
                  return false ;
               else if ( nIntType != LinCompDiscIntersType::D_NO_INTERS) {
                  for ( int nIndex = 0 ; nIndex < nDexSize ; nIndex += 1) {
                     if ( m_Values[nMap][nDex][nIndex].dMax >= dU1  &&  m_Values[nMap][nDex][nIndex].dMin <= dU2)
                        return false ;
                  }
               }
            }
           // Tronco di cono
            else {
               double dMinPar = m_Values[nMap][nDex][0].dMin ;
               double dMaxPar = m_Values[nMap][nDex][nDexSize - 1].dMax ;
               Point3d ptSegSt = ptLineSt + dMinPar * vtLineDir ;
               double dU1, dU2 ;
               int nIntType = SegmentConeFrustum( ptSegSt, vtLineDir, dMaxPar - dMinPar, ptRefPoint, vtRefAx,
                                                  dMinRad, dMaxRad, dHeight, dU1, dU2) ;
               if ( nIntType == LinCompCCIntersType::CC_ERROR_INT)
                  return false ;
               else if ( nIntType != LinCompCCIntersType::CC_NO_INTERS) {
                  for ( int nIndex = 0 ; nIndex < nDexSize ; nIndex += 1) {
                     if ( m_Values[nMap][nDex][nIndex].dMax >= dU1  &&  m_Values[nMap][nDex][nIndex].dMin <= dU2)
                        return false ;
                  }
               }
               nIntType = SegmentDisc( ptSegSt, vtLineDir, dMaxPar - dMinPar, ptRefPoint + dHeight * vtRefAx,
                                       vtRefAx, dMaxRad, dU1, dU2) ;
               if ( nIntType == LinCompDiscIntersType::D_ERROR_INT)
                  return false ;
               else if ( nIntType != LinCompDiscIntersType::D_NO_INTERS) {
                  for ( int nIndex = 0 ; nIndex < nDexSize ; nIndex += 1) {
                     if ( m_Values[nMap][nDex][nIndex].dMax >= dU1  &&  m_Values[nMap][nDex][nIndex].dMin <= dU2)
                        return false ;
                  }
               }
               nIntType = SegmentDisc( ptSegSt, vtLineDir, dMaxPar - dMinPar, ptRefPoint, vtRefAx, dMinRad, dU1, dU2) ;
               if ( nIntType == LinCompDiscIntersType::D_ERROR_INT)
                  return false ;
               else if ( nIntType != LinCompDiscIntersType::D_NO_INTERS) {
                  for ( int nIndex = 0 ; nIndex < nDexSize ; nIndex += 1) {
                     if ( m_Values[nMap][nDex][nIndex].dMax >= dU1  &&  m_Values[nMap][nDex][nIndex].dMin <= dU2)
                        return false ;
                  }
               }
            }
         }
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
// Ha come asse Z l'asse del cono/tronco di cono. Se è un cono proprio l'origine 
// è nel vertice del cono e l'asse Z è diretto verso l'apertura.
// Se è un tronco di cono l'origine è nel centro della base Bot e l'azze Z è diretto verso
// la base Top, a prescindere da quale base abbia raggio maggiore.
bool
VolZmap::AvoidConeFrustum( const Frame3d& frCone, double dRadBot, double dRadTop, double dHeight,
                           double dSafeDist, bool bPrecise) const
{
  // Se cilindro
   if ( abs( dRadBot - dRadTop) < EPS_SMALL)
      return AvoidCylinder( frCone, max( dRadBot, dRadTop), dHeight, dSafeDist, bPrecise) ;

  // Verifico che la base minore sia in basso nel suo riferimento
   Frame3d frMyCone = frCone ;
   if ( dRadBot > dRadTop) {
      frMyCone.Translate( dHeight * frMyCone.VersZ()) ;
      frMyCone.Rotate( frMyCone.Orig(), frMyCone.VersX(), -1, 0) ;
   }
   double dMinRad = min( dRadBot, dRadTop) ;
   double dMaxRad = max( dRadBot, dRadTop) ;

  // Se distanza di sicurezza nulla
   if ( dSafeDist < EPS_SMALL)
      return AvoidSimpleConeFrustum( frMyCone, dMinRad, dMaxRad, dHeight, bPrecise) ;

  // Se vero e proprio cono
   if ( dMinRad < EPS_SMALL)  {
     // Se non c'è collisione con l'offset esteso, ho finito.
      double dDeltaVert = dSafeDist * sqrt( 1 + ( dHeight * dHeight / (dMaxRad * dMaxRad))) ;
      double dHeightExt = dHeight + dSafeDist + dDeltaVert ;
      double dRadExt = dMaxRad * dHeightExt / dHeight ;
      Frame3d frFrame = frMyCone ;
      frFrame.Translate( -dDeltaVert * frFrame.VersZ()) ;
      if ( AvoidSimpleConeFrustum( frFrame, 0., dRadExt, dHeightExt, bPrecise))
         return true ;
     // Sfera nel vertice in basso
      frFrame = frMyCone ;
      if ( ! AvoidSimpleSphere( frMyCone.Orig(), dSafeDist, bPrecise))
         return false ;
     // Tronco di cono intermedio
      double dHypo = sqrt( dMaxRad * dMaxRad + dHeight * dHeight ) ;
      double dDeltaH = dSafeDist * dMaxRad / dHypo ;
      double dDeltaR = dSafeDist * dHeight / dHypo ;
      frFrame = frMyCone ; frFrame.Translate( -dDeltaH * frFrame.VersZ()) ;
      if ( ! AvoidSimpleConeFrustum( frFrame, dDeltaR, dMaxRad + dDeltaR, dHeight, bPrecise))
         return false ;
     // Cilindro nel toro in alto
      frFrame = frMyCone ; frFrame.Translate( ( dHeight - dSafeDist) * frFrame.VersZ()) ;
      if ( ! AvoidSimpleCylinder( frFrame, dMaxRad, 2 * dSafeDist, bPrecise))
         return false ;
     // Toro in alto
      frFrame = frMyCone ; frFrame.Translate( dHeight * frFrame.VersZ()) ;
      if ( ! AvoidSimpleTorus( frFrame, dMaxRad, dSafeDist, bPrecise))
         return false ;
      return true ;
   }

  // Tronco di cono
   // Definisco un sistema di riferimento che ha per origine il centro della base
   // minore e asse Z l'asse di simmetria rivolto verso la direzione di apertura.
  // Se non c'è collisione con l'offset esteso, ho finito.
   double dDiffRad = dMaxRad - dMinRad ;
   double dDeltaVert = dSafeDist * sqrt( 1 + ( dHeight * dHeight / ( dDiffRad * dDiffRad))) ;
   double dDeltaMaxRad = ( dDeltaVert + dSafeDist) * dDiffRad / dHeight ;
   double dDeltaMinRad = ( dDeltaVert - dSafeDist) * dDiffRad / dHeight ;
   Frame3d frFrame = frMyCone ; frFrame.Translate( -dSafeDist * frFrame.VersZ()) ;
   if ( AvoidSimpleConeFrustum( frFrame, dMinRad + dDeltaMinRad, dMaxRad + dDeltaMaxRad, dHeight + 2 * dSafeDist, bPrecise))
      return true ;
  // Tronco di cono intermedio
   double dHypo = sqrt( dDiffRad * dDiffRad + dHeight * dHeight ) ;
   double dDeltaH = dSafeDist * dDiffRad / dHypo ;
   double dDeltaR = dSafeDist * dHeight / dHypo ;
   frFrame = frMyCone ; frFrame.Translate( -dDeltaH * frFrame.VersZ()) ;
   if ( ! AvoidSimpleConeFrustum( frFrame, dMinRad + dDeltaR, dMaxRad + dDeltaR, dHeight, bPrecise))
      return false ;
  // Cilindro sotto
   frFrame = frMyCone ; frFrame.Translate( - dSafeDist * frFrame.VersZ()) ;
   if ( ! AvoidSimpleCylinder( frFrame, dMinRad, 2 * dSafeDist, bPrecise))
      return false ;
  // Cilindro sopra
   frFrame.Translate( dHeight * frFrame.VersZ()) ;
   if ( ! AvoidSimpleCylinder( frFrame, dMaxRad, 2 * dSafeDist, bPrecise))
      return false ;
  // Toro sotto
   frFrame = frMyCone ;
   if ( ! AvoidSimpleTorus( frFrame, dMinRad, dSafeDist, bPrecise))
      return false ;
  // Toro sopra
   frFrame.Translate( dHeight * frFrame.VersZ()) ;
   if ( ! AvoidSimpleTorus( frFrame, dMaxRad, dSafeDist, bPrecise))
      return false ;

   return true ;
}

//----------------------------------------------------------------------------
static bool
RectPrismoidSegmentCollision( const Frame3d& frPrismoid, double dLenghtBaseX, double dLenghtBaseY,
                              double dLenghtTopX, double dLenghtTopY, double dHeight,
                              const Point3d& ptSt, const Point3d& ptEn)
{
  // Se il solido non è ben definito, non ha senso continuare.
   if ( max( dLenghtBaseX, dLenghtTopX) < EPS_SMALL  ||
        max( dLenghtBaseY, dLenghtTopY) < EPS_SMALL  ||
        dHeight < EPS_SMALL)
      return false ;

  // Porto il segmento nel sistema del tronco.
   Point3d ptMySt = ptSt ;
   ptMySt.ToLoc( frPrismoid) ;
   Point3d ptMyEn = ptEn ;
   ptMyEn.ToLoc( frPrismoid) ;

  // Se il segmento non è ben definito, non ha senso continuare.
   Vector3d vtMySeg = ptMyEn - ptMySt ;
   if ( ! vtMySeg.Normalize())
      return false ;

   double dHalfBaseX = 0.5 * dLenghtBaseX ; 
   double dHalfBaseY = 0.5 * dLenghtBaseY ;
   double dHalfTopX = 0.5 * dLenghtTopX ;
   double dHalfTopY = 0.5 * dLenghtTopY ;

  // Se almeno un vertice collide ho finito
   if ( (   ptMySt.z > - EPS_SMALL  &&  ptMySt.z < dHeight + EPS_SMALL  &&
          ( ptMySt.x + dHalfBaseX + EPS_SMALL) * dHeight > ( dHalfBaseX - dHalfTopX) * ptMySt.z  &&
          ( ptMySt.x - dHalfBaseX - EPS_SMALL) * dHeight < ( dHalfTopX - dHalfBaseX) * ptMySt.z  &&
          ( ptMySt.y + dHalfBaseY + EPS_SMALL) * dHeight > ( dHalfBaseY - dHalfTopY) * ptMySt.z  &&
          ( ptMySt.y - dHalfBaseY - EPS_SMALL) * dHeight < ( dHalfTopY - dHalfBaseY) * ptMySt.z)  ||
        (   ptMyEn.z > - EPS_SMALL  &&  ptMySt.z < dHeight + EPS_SMALL  &&
          ( ptMyEn.x + dHalfBaseX + EPS_SMALL) * dHeight > ( dHalfBaseX - dHalfTopX) * ptMyEn.z  &&
          ( ptMyEn.x - dHalfBaseX - EPS_SMALL) * dHeight < ( dHalfTopX - dHalfBaseX) * ptMyEn.z  &&
          ( ptMyEn.y + dHalfBaseY + EPS_SMALL) * dHeight > ( dHalfBaseY - dHalfTopY) * ptMyEn.z  &&
          ( ptMyEn.y - dHalfBaseY - EPS_SMALL) * dHeight < ( dHalfTopY - dHalfBaseY) * ptMyEn.z))
      return true ;
   
  // Se c'è collisione con almeno un triangolo delle facce ho finito
   Triangle3d trFaceTria1, trFaceTria2 ;
   Point3d ptInt, ptInt2 ;
  // Faccia base
   trFaceTria1.Set( Point3d( - dHalfBaseX, - dHalfBaseY, 0.),
                    Point3d( - dHalfBaseX,   dHalfBaseY, 0.),
                    Point3d(   dHalfBaseX,   dHalfBaseY, 0.)) ;
   if ( trFaceTria1.Validate()  &&  IntersLineTria( ptMySt, ptMyEn, trFaceTria1, ptInt, ptInt2) != ILTT_NO)
      return true ;
   trFaceTria2.Set( Point3d( - dHalfBaseX, - dHalfBaseY, 0.),
                    Point3d(   dHalfBaseX,   dHalfBaseY, 0.),
                    Point3d(   dHalfBaseX, - dHalfBaseY, 0.)) ;
   if ( trFaceTria2.Validate()  &&  IntersLineTria( ptMySt, ptMyEn, trFaceTria2, ptInt, ptInt2) != ILTT_NO)
      return true ;
  // Faccia top
   trFaceTria1.Set( Point3d( - dHalfTopX, - dHalfTopY, dHeight),
                    Point3d(   dHalfTopX,   dHalfTopY, dHeight),
                    Point3d( - dHalfTopX,   dHalfTopY, dHeight)) ;
   if ( trFaceTria1.Validate()  &&  IntersLineTria( ptMySt, ptMyEn, trFaceTria1, ptInt, ptInt2) != ILTT_NO)
      return true ;
   trFaceTria2.Set( Point3d( - dHalfTopX, - dHalfTopY, dHeight),
                    Point3d(   dHalfTopX, - dHalfTopY, dHeight),
                    Point3d(   dHalfTopX,   dHalfTopY, dHeight)) ;
   if ( trFaceTria2.Validate()  &&  IntersLineTria( ptMySt, ptMyEn, trFaceTria2, ptInt, ptInt2) != ILTT_NO)
      return true ;
  // Faccia laterale 1
   trFaceTria1.Set( Point3d( - dHalfBaseX, - dHalfBaseY,      0.),
                    Point3d(   dHalfTopX , - dHalfTopY , dHeight),
                    Point3d( - dHalfTopX , - dHalfTopY , dHeight)) ;
   if ( trFaceTria1.Validate()  &&  IntersLineTria( ptMySt, ptMyEn, trFaceTria1, ptInt, ptInt2) != ILTT_NO)
      return true ;
   trFaceTria2.Set( Point3d( - dHalfBaseX, - dHalfBaseY,      0.),
                    Point3d(   dHalfBaseX, - dHalfBaseY,      0.),
                    Point3d(   dHalfTopX,  - dHalfTopY , dHeight)) ;
   if ( trFaceTria2.Validate()  &&  IntersLineTria( ptMySt, ptMyEn, trFaceTria2, ptInt, ptInt2) != ILTT_NO)
      return true ;
  // Faccia laterale 2
   trFaceTria1.Set( Point3d( dHalfBaseX, - dHalfBaseY,      0.),
                    Point3d( dHalfTopX ,   dHalfTopY , dHeight),
                    Point3d( dHalfTopX , - dHalfTopY , dHeight)) ;
   if ( trFaceTria1.Validate()  &&  IntersLineTria( ptMySt, ptMyEn, trFaceTria1, ptInt, ptInt2) != ILTT_NO)
      return true ;
   trFaceTria2.Set( Point3d( dHalfBaseX, - dHalfBaseY,      0.),
                    Point3d( dHalfBaseX,   dHalfBaseY,      0.),
                    Point3d( dHalfTopX ,   dHalfTopY , dHeight)) ;
   if ( trFaceTria2.Validate()  &&  IntersLineTria( ptMySt, ptMyEn, trFaceTria2, ptInt, ptInt2) != ILTT_NO)
      return true ;
  // Faccia laterale 3
   trFaceTria1.Set( Point3d(   dHalfBaseX, dHalfBaseY,      0.),
                    Point3d( - dHalfTopX , dHalfTopY , dHeight),
                    Point3d(   dHalfTopX , dHalfTopY , dHeight)) ;
   if ( trFaceTria1.Validate()  &&  IntersLineTria( ptMySt, ptMyEn, trFaceTria1, ptInt, ptInt2) != ILTT_NO)
      return true ;
   trFaceTria2.Set( Point3d(   dHalfBaseX, dHalfBaseY,      0.),
                    Point3d( - dHalfBaseX, dHalfBaseY,      0.),
                    Point3d( - dHalfTopX , dHalfTopY , dHeight)) ;
   if ( trFaceTria2.Validate()  &&  IntersLineTria( ptMySt, ptMyEn, trFaceTria2, ptInt, ptInt2) != ILTT_NO)
      return true ;
  // Faccia laterale 4
   trFaceTria1.Set( Point3d( - dHalfBaseX,   dHalfBaseY,      0.),
                    Point3d( - dHalfTopX , - dHalfTopY , dHeight),
                    Point3d( - dHalfTopX ,   dHalfTopY , dHeight)) ;
   if ( trFaceTria1.Validate()  &&  IntersLineTria( ptMySt, ptMyEn, trFaceTria1, ptInt, ptInt2) != ILTT_NO)
      return true ;
   trFaceTria2.Set( Point3d( - dHalfBaseX,   dHalfBaseY,      0.),
                    Point3d( - dHalfBaseX, - dHalfBaseY,      0.),
                    Point3d( - dHalfTopX , - dHalfTopY , dHeight)) ;
   if ( trFaceTria2.Validate()  &&  IntersLineTria( ptMySt, ptMyEn, trFaceTria2, ptInt, ptInt2) != ILTT_NO)
      return true ;

   return false ;
}

//----------------------------------------------------------------------------
static bool
IntersSegmentTrianglePlus( const Point3d& ptTr1, const Point3d& ptTr2, const Point3d& ptTr3,
                           const Point3d& ptLnSt, const Vector3d& vtLnDir, double dLnLen,
                           double& dStU, double& dEnU)
{
   Triangle3d trFaceTria ;
   trFaceTria.Set( ptTr1, ptTr2, ptTr3) ;
   if ( trFaceTria.Validate()) {
      Point3d ptMyIntSt, ptMyIntEn ;
      int nIntType = IntersLineTria( ptLnSt, vtLnDir, dLnLen, trFaceTria, ptMyIntSt, ptMyIntEn, true) ;
      if ( nIntType != ILTT_NO) {
         if ( nIntType == ILTT_VERT  || nIntType == ILTT_EDGE  ||  nIntType == ILTT_IN) {
            double dCurU = ( ptMyIntSt - ptLnSt) * vtLnDir ;
            if ( dCurU < dStU)
               dStU = dCurU ;
            if ( dCurU > dEnU)
               dEnU = dCurU ;
         }
         else {
            double dCurStU = ( ptMyIntSt - ptLnSt) * vtLnDir ;
            double dCurEnU = ( ptMyIntEn - ptLnSt) * vtLnDir ;
            if ( dCurStU <  dStU)
               dStU = dCurStU ;
            if ( dCurEnU > dEnU)
               dEnU = dCurEnU ;
         }
      }
      return true ;
   }

   return false ;
}

//----------------------------------------------------------------------------
static bool
RectPrismoidSegmentCollisionPlus( const Frame3d& frPrismoid, double dLenghtBaseX, double dLenghtBaseY,
                                  double dLenghtTopX, double dLenghtTopY, double dHeight,
                                  const Point3d& ptSt, const Point3d& ptEn,
                                  double& dStU, double& dEnU)
{
  // Se il solido non è ben definito, non ha senso continuare.
   if ( max( dLenghtBaseX, dLenghtTopX) < EPS_SMALL  ||
        max( dLenghtBaseY, dLenghtTopY) < EPS_SMALL  ||
        dHeight < EPS_SMALL)
      return false ;

  // Porto il segmento nel sistema del prismoide a base rettangolare
   Point3d ptMySt = ptSt ;
   ptMySt.ToLoc( frPrismoid) ;
   Point3d ptMyEn = ptEn ;
   ptMyEn.ToLoc( frPrismoid) ;

  // Se il segmento non è ben definito, non ha senso continuare.
   Vector3d vtMySeg = ptMyEn - ptMySt ;
   double dSegLen = vtMySeg.Len() ;
   if ( dSegLen < EPS_SMALL)
      return false ;
   vtMySeg /= dSegLen ;

  // Semidimensioni delle basi
   double dHalfBaseX = 0.5 * dLenghtBaseX ; 
   double dHalfBaseY = 0.5 * dLenghtBaseY ;
   double dHalfTopX = 0.5 * dLenghtTopX ;
   double dHalfTopY = 0.5 * dLenghtTopY ;

  // Inizializzo gli estremi della parte di retta che interseca il prismoide
   dStU = INFINITO ;
   dEnU = -INFINITO ;
  // Interseco la retta con le facce e salvo i punti d'intersezione
  // Faccia base
   IntersSegmentTrianglePlus( Point3d( -dHalfBaseX, -dHalfBaseY, 0.),
                              Point3d( -dHalfBaseX,  dHalfBaseY, 0.),
                              Point3d(  dHalfBaseX,  dHalfBaseY, 0.),
                              ptMySt, vtMySeg, dSegLen,
                              dStU, dEnU) ;
   IntersSegmentTrianglePlus( Point3d( -dHalfBaseX, -dHalfBaseY, 0.),
                              Point3d(  dHalfBaseX,  dHalfBaseY, 0.),
                              Point3d(  dHalfBaseX, -dHalfBaseY, 0.),
                              ptMySt, vtMySeg, dSegLen,
                              dStU, dEnU) ;
  // Faccia top
   IntersSegmentTrianglePlus( Point3d( -dHalfTopX, -dHalfTopY, dHeight),
                              Point3d(  dHalfTopX,  dHalfTopY, dHeight),
                              Point3d( -dHalfTopX,  dHalfTopY, dHeight),
                              ptMySt, vtMySeg, dSegLen,
                              dStU, dEnU) ;
   IntersSegmentTrianglePlus( Point3d( -dHalfTopX, -dHalfTopY, dHeight),
                              Point3d(  dHalfTopX, -dHalfTopY, dHeight),
                              Point3d(  dHalfTopX,  dHalfTopY, dHeight),
                              ptMySt, vtMySeg, dSegLen,
                              dStU, dEnU) ;
  // Faccia laterale 1
   IntersSegmentTrianglePlus( Point3d( -dHalfBaseX, -dHalfBaseY,      0.),
                              Point3d(  dHalfTopX , -dHalfTopY , dHeight),
                              Point3d( -dHalfTopX , -dHalfTopY , dHeight),
                              ptMySt, vtMySeg, dSegLen,
                              dStU, dEnU) ;
   IntersSegmentTrianglePlus( Point3d( -dHalfBaseX, -dHalfBaseY,      0.),
                              Point3d(  dHalfBaseX, -dHalfBaseY,      0.),
                              Point3d(  dHalfTopX,  -dHalfTopY , dHeight),
                              ptMySt, vtMySeg, dSegLen,
                              dStU, dEnU) ;
  // Faccia laterale 2
   IntersSegmentTrianglePlus( Point3d( dHalfBaseX, -dHalfBaseY,      0.),
                              Point3d( dHalfTopX ,  dHalfTopY , dHeight),
                              Point3d( dHalfTopX , -dHalfTopY , dHeight),
                              ptMySt, vtMySeg, dSegLen,
                              dStU, dEnU) ;
   IntersSegmentTrianglePlus( Point3d( dHalfBaseX, -dHalfBaseY,      0.),
                              Point3d( dHalfBaseX,  dHalfBaseY,      0.),
                              Point3d( dHalfTopX ,  dHalfTopY , dHeight),
                              ptMySt, vtMySeg, dSegLen,
                              dStU, dEnU) ;
  // Faccia laterale 3
   IntersSegmentTrianglePlus( Point3d(  dHalfBaseX, dHalfBaseY,      0.),
                              Point3d( -dHalfTopX , dHalfTopY , dHeight),
                              Point3d(  dHalfTopX , dHalfTopY , dHeight),
                              ptMySt, vtMySeg, dSegLen,
                              dStU, dEnU) ;
   IntersSegmentTrianglePlus( Point3d(  dHalfBaseX, dHalfBaseY,      0.),
                              Point3d( -dHalfBaseX, dHalfBaseY,      0.),
                              Point3d( -dHalfTopX , dHalfTopY , dHeight),
                              ptMySt, vtMySeg, dSegLen,
                              dStU, dEnU) ;
  // Faccia laterale 4
   IntersSegmentTrianglePlus( Point3d( -dHalfBaseX,  dHalfBaseY,      0.),
                              Point3d( -dHalfTopX , -dHalfTopY , dHeight),
                              Point3d( -dHalfTopX ,  dHalfTopY , dHeight),
                              ptMySt, vtMySeg, dSegLen,
                              dStU, dEnU) ;
   IntersSegmentTrianglePlus( Point3d( -dHalfBaseX,  dHalfBaseY,      0.),
                              Point3d( -dHalfBaseX, -dHalfBaseY,      0.),
                              Point3d( -dHalfTopX , -dHalfTopY , dHeight),
                              ptMySt, vtMySeg, dSegLen,
                              dStU, dEnU) ;

   return ( dEnU > dStU - EPS_ZERO  &&  dStU < dSegLen + EPS_SMALL  &&  dEnU > -EPS_SMALL) ;
}

//----------------------------------------------------------------------------
bool
VolZmap::AvoidSimpleRectPrismoid( const Frame3d& frPrismoid, double dLenghtBaseX, double dLenghtBaseY,
                                  double dLenghtTopX, double dLenghtTopY, double dHeight, bool bPrecise) const
{
  // Box del tronco di prismoide nel suo sistema locale
   double dMaxLenX = max( dLenghtBaseX, dLenghtTopX) ;
   double dMaxLenY = max( dLenghtBaseY, dLenghtTopY) ;
   BBox3d b3PrismL( Point3d( -dMaxLenX / 2, -dMaxLenY / 2, 0.),
                    Point3d(  dMaxLenX / 2,  dMaxLenY / 2, dHeight)) ;

  // BBox del tronco di prismoide nel riferimento intrinseco dello Zmap
   Frame3d frPrismInt = GetToLoc( frPrismoid, m_MapFrame) ;
   BBox3d b3PrismI = GetToGlob( b3PrismL, frPrismInt) ;

  // BBox dello Zmap nel suo riferimento intrinseco
   BBox3d b3Zmap( ORIG, Point3d( m_nNx[0] * m_dStep, m_nNy[0] * m_dStep, m_dMaxZ[0])) ;

   // Se i box non interferiscono, posso uscire
   if ( ! b3Zmap.Overlaps( b3PrismI)  ||  ! b3Zmap.Overlaps( frPrismInt, b3PrismL))
      return true ;
   BBox3d b3Int ;
   if ( ! b3Zmap.FindIntersection( b3PrismI, b3Int))
      return true ;

   if ( ! bPrecise) {
     // Limiti su indici
      int nStI = Clamp( int( b3Int.GetMin().x / m_dStep), 0, m_nNx[0] - 1) ;
      int nEnI = Clamp( int( b3Int.GetMax().x / m_dStep), 0, m_nNx[0] - 1) ;
      int nStJ = Clamp( int( b3Int.GetMin().y / m_dStep), 0, m_nNy[0] - 1) ;
      int nEnJ = Clamp( int( b3Int.GetMax().y / m_dStep), 0, m_nNy[0] - 1) ;
     // Ciclo di intersezione dei dexel con il cilindro (nel riferimento intrinseco)
      for ( int i = nStI ; i <= nEnI ; ++ i) {
         for ( int j = nStJ ; j <= nEnJ ; ++ j) {
            int nPos = j * m_nNx[0] + i ;
            int nSize = int( m_Values[0][nPos].size()) ;
            if ( nSize == 0)
               continue ;
            double dParMin = m_Values[0][nPos][0].dMin ;
            double dParMax = m_Values[0][nPos][nSize-1].dMax ;
            if ( dParMax < b3Int.GetMin().z  ||  dParMin > b3Int.GetMax().z)
               continue ;
            for  ( int k = 0 ; k < 5 ; ++ k) {
               Point3d ptLineSt = ORIG + ( i + 0.5) * m_dStep * X_AX + ( j + 0.5) * m_dStep * Y_AX ;
               switch ( k) {
               case 0 : break ;
               case 1 : ptLineSt += -0.4 * m_dStep * X_AX - 0.4 * m_dStep * Y_AX ; break ;
               case 2 : ptLineSt += +0.4 * m_dStep * X_AX - 0.4 * m_dStep * Y_AX ; break ;
               case 3 : ptLineSt += +0.4 * m_dStep * X_AX + 0.4 * m_dStep * Y_AX ; break ;
               case 4 : ptLineSt += -0.4 * m_dStep * X_AX + 0.4 * m_dStep * Y_AX ; break ;
               }
               double dStU, dEnU ;
               Point3d ptSegSt = ptLineSt + m_Values[0][nPos][0].dMin * Z_AX ;
               Point3d ptSegEn = ptLineSt + m_Values[0][nPos][nSize-1].dMax * Z_AX ;
               if ( RectPrismoidSegmentCollisionPlus( frPrismInt, dLenghtBaseX, dLenghtBaseY, dLenghtTopX, dLenghtTopY,
                                                      dHeight, ptSegSt, ptSegEn, dStU, dEnU)) {       
                  for ( int nIndex = 0 ; nIndex < nSize ; nIndex += 1) {
                     if ( m_Values[0][nPos][nIndex].dMax >= dStU  &&  m_Values[0][nPos][nIndex].dMin <= dEnU)
                        return false ;
                  }
               }
            }
         }
      }
   }
   else { 
     // Ciclo sulle mappe
      for ( int nMap = 0 ; nMap < m_nMapNum ; ++ nMap) {
         Point3d ptInfIntBox = b3Int.GetMin();
         Point3d ptSupIntBox = b3Int.GetMax();
        // Dal sistema intrinseco al sistema griglia (per la prima griglia coincidono).
         if ( nMap == 1) {
            swap( ptInfIntBox.x, ptInfIntBox.z) ;
            swap( ptInfIntBox.x, ptInfIntBox.y) ;
            swap( ptSupIntBox.x, ptSupIntBox.z) ;
            swap( ptSupIntBox.x, ptSupIntBox.y) ;
         }
         else if ( nMap == 2) {
            swap( ptInfIntBox.y, ptInfIntBox.z) ;
            swap( ptInfIntBox.x, ptInfIntBox.y) ;
            swap( ptSupIntBox.y, ptSupIntBox.z) ;
            swap( ptSupIntBox.x, ptSupIntBox.y) ;
         }
        // Limiti su indici
         int nStI = Clamp( int( ptInfIntBox.x / m_dStep), 0, m_nNx[nMap] - 1) ;
         int nEnI = Clamp( int( ptSupIntBox.x / m_dStep), 0, m_nNx[nMap] - 1) ;
         int nStJ = Clamp( int( ptInfIntBox.y / m_dStep), 0, m_nNy[nMap] - 1) ;
         int nEnJ = Clamp( int( ptSupIntBox.y / m_dStep), 0, m_nNy[nMap] - 1) ;
         for ( int nDex = 0 ; nDex < int( m_Values[nMap].size()) ; ++ nDex) {
            int nSize = int( m_Values[nMap][nDex].size()) ;
            if ( nSize == 0)
               continue ;
           // Indici del dexel
            int nI = nDex % m_nNx[nMap] ;
            int nJ = nDex / m_nNx[nMap] ;
           // Se fuori dalla regione ammissibile salto l'iterazione
            if ( nI < nStI  ||  nI > nEnI  ||  nJ < nStJ  ||  nJ > nEnJ)
               continue ;
           // Posizione del dexel
            double dX = ( nI + 0.5) * m_dStep ;
            double dY = ( nJ + 0.5) * m_dStep ;
            Point3d ptLineSt( dX, dY, 0.) ;
            Vector3d vtLineDir( 0., 0., 1.) ;
           // Dal sistema griglia al sistema intrinseco (per la prima griglia coincidono).
            if ( nMap == 1) {
               swap( ptLineSt.x, ptLineSt.y) ;
               swap( ptLineSt.x, ptLineSt.z) ;
               swap( vtLineDir.x, vtLineDir.y) ;
               swap( vtLineDir.x, vtLineDir.z) ;
            }
            else if ( nMap == 2) {
               swap( ptLineSt.x, ptLineSt.y) ;
               swap( ptLineSt.y, ptLineSt.z) ;
               swap( vtLineDir.x, vtLineDir.y) ;
               swap( vtLineDir.y, vtLineDir.z) ;
            }
            double dStU, dEnU ;
            Point3d ptSegSt = ptLineSt + m_Values[nMap][nDex][0].dMin * vtLineDir ;
            Point3d ptSegEn = ptLineSt + m_Values[nMap][nDex][nSize-1].dMax * vtLineDir ;
            if ( RectPrismoidSegmentCollisionPlus( frPrismInt, dLenghtBaseX, dLenghtBaseY, dLenghtTopX, dLenghtTopY,
                                                   dHeight, ptSegSt, ptSegEn, dStU, dEnU)) {
               for ( int nIndex = 0 ; nIndex < nSize ; nIndex += 1) {
                  if ( m_Values[nMap][nDex][nIndex].dMax >= dStU  &&  m_Values[nMap][nDex][nIndex].dMin <= dEnU)
                     return false ;
               }
            }
         }
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::AvoidRectPrismoid( const Frame3d& frPrismoid, double dLenghtBaseX, double dLenghtBaseY,
                            double dLenghtTopX, double dLenghtTopY, double dHeight,
                            double dSafeDist, bool bPrecise) const
{
  // Se il tronco di piramide non è ben definito non procedo
   if ( max( dLenghtBaseX, dLenghtTopX) < EPS_SMALL  ||
        max( dLenghtBaseY, dLenghtTopY) < EPS_SMALL  ||
        dHeight < EPS_SMALL)
      return false ;

  // Se distanza di sicurezza nulla
   if ( dSafeDist < EPS_SMALL)
      return AvoidSimpleRectPrismoid( frPrismoid, dLenghtBaseX, dLenghtBaseY, dLenghtTopX, dLenghtTopY, dHeight, bPrecise) ;

  // Verifica preliminare con offset esteso
   double dHDiffX = ( dLenghtBaseX - dLenghtTopX) / 2 ;
   double dTgAx = dHDiffX / dHeight ;
   double dSecAx = sqrt( 1 + dTgAx * dTgAx) ;
   double dOffsBaseX = dSafeDist * ( dSecAx + dTgAx) ;
   double dOffsTopX = dSafeDist * ( dSecAx - dTgAx) ;
   double dHDiffY = ( dLenghtBaseY - dLenghtTopY) / 2 ;
   double dTgAy = dHDiffY / dHeight ;
   double dSecAy = sqrt( 1 + dTgAy * dTgAy) ;
   double dOffsBaseY = dSafeDist * ( dSecAy + dTgAy) ;
   double dOffsTopY = dSafeDist * ( dSecAy - dTgAy) ;
   Frame3d frFrame = frPrismoid ; frFrame.Translate( -dSafeDist * frFrame.VersZ()) ;
   if ( AvoidSimpleRectPrismoid( frFrame, dLenghtBaseX + 2 * dOffsBaseX, dLenghtBaseY + 2 * dOffsBaseY,
                                 dLenghtTopX + 2 * dOffsTopX, dLenghtTopY + 2 * dOffsTopY, dHeight + 2 * dSafeDist, bPrecise))
      return true ;

  // Offset fine
  // Sfere centrate nei vertici
   double dHalfBaseX = dLenghtBaseX / 2 ;
   double dHalfBaseY = dLenghtBaseY / 2 ;
   double dHalfTopX = dLenghtTopX / 2 ;
   double dHalfTopY = dLenghtTopY / 2 ;
   PNTVECTOR vVert = { Point3d( -dHalfBaseX, -dHalfBaseY, 0.),
                       Point3d( dHalfBaseX, -dHalfBaseY, 0.),
                       Point3d( dHalfBaseX, dHalfBaseY, 0.),
                       Point3d( -dHalfBaseX, dHalfBaseY, 0.),
                       Point3d( -dHalfTopX, -dHalfTopY, dHeight),
                       Point3d( dHalfTopX, -dHalfTopY, dHeight),
                       Point3d( dHalfTopX, dHalfTopY, dHeight),
                       Point3d( -dHalfTopX, dHalfTopY, dHeight)} ;
   for ( auto& ptV : vVert) {
      ptV.ToGlob( frPrismoid) ;
      if ( ! AvoidSimpleSphere( ptV, dSafeDist, bPrecise))
         return false ;
   }
   // Cilindri con i segmenti come asse
   frFrame.Set( vVert[0], frPrismoid.VersX()) ;
   if ( ! AvoidSimpleCylinder( frFrame, dSafeDist, dLenghtBaseX, bPrecise))
      return false ;
   frFrame.Set( vVert[1], frPrismoid.VersY()) ;
   if ( ! AvoidSimpleCylinder( frFrame, dSafeDist, dLenghtBaseY, bPrecise))
      return false ;
   frFrame.Set( vVert[2], -frPrismoid.VersX()) ;
   if ( ! AvoidSimpleCylinder( frFrame, dSafeDist, dLenghtBaseX, bPrecise))
      return false ;
   frFrame.Set( vVert[3], -frPrismoid.VersY()) ;
   if ( ! AvoidSimpleCylinder( frFrame, dSafeDist, dLenghtBaseY, bPrecise))
      return false ;
   frFrame.Set( vVert[4], frPrismoid.VersX()) ;
   if ( ! AvoidSimpleCylinder( frFrame, dSafeDist, dLenghtTopX, bPrecise))
      return false ;
   frFrame.Set( vVert[5], frPrismoid.VersY()) ;
   if ( ! AvoidSimpleCylinder( frFrame, dSafeDist, dLenghtTopY, bPrecise))
      return false ;
   frFrame.Set( vVert[6], -frPrismoid.VersX()) ;
   if ( ! AvoidSimpleCylinder( frFrame, dSafeDist, dLenghtTopX, bPrecise))
      return false ;
   frFrame.Set( vVert[7], -frPrismoid.VersY()) ;
   if ( ! AvoidSimpleCylinder( frFrame, dSafeDist, dLenghtTopY, bPrecise))
      return false ;
   Vector3d vtSeg04 = vVert[4] - vVert[0] ;
   double dLenSeg04 = vtSeg04.Len() ;
   frFrame.Set( vVert[0], vtSeg04) ;
   if ( ! AvoidSimpleCylinder( frFrame, dSafeDist, dLenSeg04, bPrecise))
      return false ;
   Vector3d vtSeg15 = vVert[5] - vVert[1] ;
   double dLenSeg15 = vtSeg15.Len() ;
   frFrame.Set( vVert[1], vtSeg15) ;
   if ( ! AvoidSimpleCylinder( frFrame, dSafeDist, dLenSeg15, bPrecise))
      return false ;
   Vector3d vtSeg26 = vVert[6] - vVert[2] ;
   double dLenSeg26 = vtSeg26.Len() ;
   frFrame.Set( vVert[2], vtSeg26) ;
   if ( ! AvoidSimpleCylinder( frFrame, dSafeDist, dLenSeg26, bPrecise))
      return false ;
   Vector3d vtSeg37 = vVert[7] - vVert[3] ;
   double dLenSeg37 = vtSeg37.Len();
   frFrame.Set( vVert[3], vtSeg37) ;
   if ( ! AvoidSimpleCylinder( frFrame, dSafeDist, dLenSeg37, bPrecise))
      return false ;
  // Box sotto
   frFrame = frPrismoid ; frFrame.Translate( -dSafeDist * frFrame.VersZ()) ;
   if ( ! AvoidSimpleBox( frFrame, Vector3d( dLenghtBaseX, dLenghtBaseY, dSafeDist), bPrecise))
      return false ;
  // Box sopra
   frFrame = frPrismoid ; frFrame.Translate( dHeight * frFrame.VersZ()) ;
   if ( ! AvoidSimpleBox( frFrame, Vector3d( dLenghtBaseX, dLenghtBaseY, dSafeDist), bPrecise))
      return false ;
  // Prismoide allungato in X
   double dHypoX = sqrt( dHDiffX * dHDiffX + dHeight * dHeight) ;
   double dOffsX = dSafeDist * dHeight / dHypoX ;
   double dMoveXZ = dSafeDist * dHDiffX / dHypoX ;
   frFrame = frPrismoid ; frFrame.Translate( dMoveXZ * frFrame.VersZ()) ;
   if ( ! AvoidSimpleRectPrismoid( frFrame, dLenghtBaseX + 2 * dOffsX, dLenghtBaseY,
                                   dLenghtTopX + 2 * dOffsX, dLenghtTopY, dHeight, bPrecise))
      return false ;
  // Prismoide allungato in Y
   double dHypoY = sqrt( dHDiffY * dHDiffY + dHeight * dHeight) ;
   double dOffsY = dSafeDist * dHeight / dHypoY ;
   double dMoveYZ = dSafeDist * dHDiffY / dHypoY ;
   frFrame = frPrismoid ; frFrame.Translate( dMoveYZ * frFrame.VersZ()) ;
   if ( ! AvoidSimpleRectPrismoid( frFrame, dLenghtBaseX, dLenghtBaseY + 2 * dOffsY,
                                   dLenghtTopX, dLenghtTopY + 2 * dOffsY, dHeight, bPrecise))
      return false ;

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::AvoidSimpleTorus( const Frame3d& frTorus, double dMaxRad, double dMinRad, bool bPrecise) const
{
  // BBox del toro in locale
   BBox3d b3TorusL( Point3d( -dMaxRad - dMinRad, -dMaxRad - dMinRad, -dMinRad),
                    Point3d(  dMaxRad + dMinRad,  dMaxRad + dMinRad,  dMinRad)) ;

  // BBox del toro nel riferimento intrinseco dello Zmap
   Frame3d frTorusInt = GetToLoc( frTorus, m_MapFrame) ;
   BBox3d b3TorusI = GetToGlob( b3TorusL, frTorusInt) ;

  // BBox dello Zmap nel suo riferimento intrinseco
   BBox3d b3Zmap( ORIG, Point3d( m_nNx[0] * m_dStep, m_nNy[0] * m_dStep, m_dMaxZ[0])) ;

  // Se non interferiscono, posso uscire
   if ( ! b3Zmap.Overlaps( b3TorusI)  ||  ! b3Zmap.Overlaps( frTorusInt, b3TorusL))
      return true ;

  // BBox del toro ottimizzato nel riferimento intrinseco dello Zmap
   Point3d ptMyCen = frTorusInt.Orig() ;
   Vector3d vtMyAx = frTorusInt.VersZ() ;
   BBox3d b3Torus( ptMyCen) ;
   b3Torus.Add( ptMyCen + vtMyAx * dMinRad) ;
   b3Torus.Add( ptMyCen - vtMyAx * dMinRad) ;
   double dTotRad = dMaxRad + dMinRad ;
   if ( vtMyAx.IsX())
      b3Torus.Expand( 0, dTotRad, dTotRad) ;
   else if ( vtMyAx.IsY())
      b3Torus.Expand( dTotRad, 0, dTotRad) ;
   else if ( vtMyAx.IsZ())
      b3Torus.Expand( dTotRad, dTotRad, 0) ;
   else {
      double dExpandX = dTotRad * sqrt( 1 - vtMyAx.x * vtMyAx.x) ;
      double dExpandY = dTotRad * sqrt( 1 - vtMyAx.y * vtMyAx.y) ;
      double dExpandZ = dTotRad * sqrt( 1 - vtMyAx.z * vtMyAx.z) ;
      b3Torus.Expand( dExpandX, dExpandY, dExpandZ) ;
   }

  // Se non interferiscono, posso uscire
   BBox3d b3Int ;
   if ( ! b3Zmap.FindIntersection( b3Torus, b3Int))
      return true ;

  // Se verifico solo prima mappa
   if ( ! bPrecise) {
     // Limiti su indici
      int nStI = Clamp( int( b3Int.GetMin().x / m_dStep), 0, m_nNx[0] - 1) ;
      int nEnI = Clamp( int( b3Int.GetMax().x / m_dStep), 0, m_nNx[0] - 1) ;
      int nStJ = Clamp( int( b3Int.GetMin().y / m_dStep), 0, m_nNy[0] - 1) ;
      int nEnJ = Clamp( int( b3Int.GetMax().y / m_dStep), 0, m_nNy[0] - 1) ;
     // Ciclo di intersezione dei dexel con il toro (nel riferimento intrinseco)
      for ( int i = nStI ; i <= nEnI ; ++ i) {
         for ( int j = nStJ ; j <= nEnJ ; ++ j) {
            int nPos = j * m_nNx[0] + i ;
            int nSize = int( m_Values[0][nPos].size()) ;
            if ( nSize == 0)
               continue ;
            double dParMin = m_Values[0][nPos][0].dMin ;
            double dParMax = m_Values[0][nPos][nSize-1].dMax ;
            if ( dParMax < b3Int.GetMin().z  ||  dParMin > b3Int.GetMax().z)
               continue ;
            for ( int k = 0 ; k < 5 ; ++ k) {
               Point3d ptLineSt = ORIG + ( i + 0.5) * m_dStep * X_AX + ( j + 0.5) * m_dStep * Y_AX ;
               switch ( k) {
               case 0 : break ;
               case 1 : ptLineSt += -0.4 * m_dStep * X_AX - 0.4 * m_dStep * Y_AX ; break ;
               case 2 : ptLineSt += +0.4 * m_dStep * X_AX - 0.4 * m_dStep * Y_AX ; break ;
               case 3 : ptLineSt += +0.4 * m_dStep * X_AX + 0.4 * m_dStep * Y_AX ; break ;
               case 4 : ptLineSt += -0.4 * m_dStep * X_AX + 0.4 * m_dStep * Y_AX ; break ;
               }
              // Intersezione segmento toro
               Point3d ptSegSt = ptLineSt + dParMin * Z_AX ;
               BOOLVECTOR vbType ;
               DBLVECTOR  vdPar ;
               int nIntType = SegmentTorus( ptSegSt, Z_AX, dParMax - dParMin, ptMyCen, vtMyAx, dMinRad, dMaxRad,
                                            vbType, vdPar) ;
               if ( nIntType == LinCompTorusIntersType::T_ERROR)
                  return false ;
               else if ( nIntType != LinCompTorusIntersType::T_NO_INT) {
                  double dUmin = vdPar.front() ;
                  double dUmax = vdPar.back() ;
                  for ( int nIndex = 0 ; nIndex < nSize ; nIndex += 1) {
                     if ( m_Values[0][nPos][nIndex].dMax >= dUmin  &&  m_Values[0][nPos][nIndex].dMin <= dUmax)
                        return false ;
                  }
               }
            }
         }
      }
   }
   else {
     // Ciclo sulle mappe
      for ( int nMap = 0 ; nMap < m_nMapNum ; ++ nMap) {
         Point3d ptInfIntBox = b3Int.GetMin() ;
         Point3d ptSupIntBox = b3Int.GetMax() ;
        // Dal sistema intrinseco al sistema griglia (per la prima griglia coincidono).
         if ( nMap == 1) {
            swap( ptInfIntBox.x, ptInfIntBox.z) ;
            swap( ptInfIntBox.x, ptInfIntBox.y) ;
            swap( ptSupIntBox.x, ptSupIntBox.z) ;
            swap( ptSupIntBox.x, ptSupIntBox.y) ;
         }
         else if (nMap == 2) {
            swap( ptInfIntBox.y, ptInfIntBox.z) ;
            swap( ptInfIntBox.x, ptInfIntBox.y) ;
            swap( ptSupIntBox.y, ptSupIntBox.z) ;
            swap( ptSupIntBox.x, ptSupIntBox.y) ;
         }
        // Limiti su indici
         int nStI = Clamp( int( ptInfIntBox.x / m_dStep), 0, m_nNx[nMap] - 1) ;
         int nEnI = Clamp( int( ptSupIntBox.x / m_dStep), 0, m_nNx[nMap] - 1) ;
         int nStJ = Clamp( int( ptInfIntBox.y / m_dStep), 0, m_nNy[nMap] - 1) ;
         int nEnJ = Clamp( int( ptSupIntBox.y / m_dStep), 0, m_nNy[nMap] - 1) ;
         for ( int nDex = 0 ; nDex < int( m_Values[nMap].size()) ; ++ nDex) {
            int nSize = int( m_Values[nMap][nDex].size()) ;
            if ( nSize == 0)
               continue ;
            double dParMin = m_Values[nMap][nDex][0].dMin ;
            double dParMax = m_Values[nMap][nDex][nSize-1].dMax ;
            if ( dParMax < ptInfIntBox.z  ||  dParMin > ptSupIntBox.z)
               continue ;
           // Indici del dexel
            int nI = nDex % m_nNx[nMap] ;
            int nJ = nDex / m_nNx[nMap] ;
           // Se fuori dalla regione ammissibile salto l'iterazione
            if ( nI < nStI  ||  nI > nEnI  ||  nJ < nStJ  ||  nJ > nEnJ)
               continue ;
           // Posizione del dexel
            double dX = ( nI + 0.5) * m_dStep ;
            double dY = ( nJ + 0.5) * m_dStep ;
            Point3d ptLineSt( dX, dY, 0.) ;
            Vector3d vtLineDir( 0., 0., 1.) ;
           // Dal sistema griglia al sistema intrinseco (per la prima griglia coincidono).
            if ( nMap == 1) {
               swap( ptLineSt.x, ptLineSt.y) ;
               swap( ptLineSt.x, ptLineSt.z) ;
               swap( vtLineDir.x, vtLineDir.y) ;
               swap( vtLineDir.x, vtLineDir.z) ;
            }
            else if ( nMap == 2) {
               swap( ptLineSt.x, ptLineSt.y) ;
               swap( ptLineSt.y, ptLineSt.z) ;
               swap( vtLineDir.x, vtLineDir.y) ;
               swap( vtLineDir.y, vtLineDir.z) ;
            }
            BOOLVECTOR vbType ;
            DBLVECTOR  vdPar ;
            Point3d ptSegSt = ptLineSt + dParMin * vtLineDir ;
            int nIntType = SegmentTorus( ptSegSt, vtLineDir, dParMax - dParMin, ptMyCen, vtMyAx, dMinRad, dMaxRad,
                                         vbType, vdPar) ;
            if ( nIntType == LinCompTorusIntersType::T_ERROR)
               return false ;
            else if ( nIntType != LinCompTorusIntersType::T_NO_INT) {
               double dUmin = vdPar.front() ;
               double dUmax = vdPar.back() ;
               for ( int nIndex = 0 ; nIndex < nSize ; nIndex += 1) {
                  if ( m_Values[nMap][nDex][nIndex].dMax >= dUmin  &&  m_Values[nMap][nDex][nIndex].dMin <= dUmax)
                     return false ;
               }
            }
         }
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::AvoidTorus( const Frame3d& frTorus, double dMaxRad, double dMinRad,
                     double dSafeDist, bool bPrecise) const
{
  // I raggi devono essere non nulli
   if ( dMaxRad < EPS_SMALL  ||  dMinRad < EPS_SMALL)
      return true ;

   return AvoidSimpleTorus( frTorus, dMaxRad, dMinRad + max( 0., dSafeDist), bPrecise) ;
}

//----------------------------------------------------------------------------
// La superficie deve essere valida e chiusa.
// La superficie deve essere espressa nel sistema locale, in cui è definito quello intrinseco.
// I conti sui box per scremare i triangoli e i dexel vengono eseguiti nel sistema intrinsco,
// quelli di intersezione segmento triangolo vengono eseguiti nel sistema locale. Questo perché
// è più veloce trasformare le coordinate degli estremi del segmento piuttosto che quelle del
// triangolo.
bool
VolZmap::AvoidSurfTm( const ISurfTriMesh& tmSurf, double dSafeDist, bool bPrecise) const
{
   // Controllo sulla validità della superficie ed eventualmente sulla sua chiusura
   if ( ! ( tmSurf.IsValid()  &&  tmSurf.IsClosed()))
      return false ;
   // Bounding box della superficie espresso nel sistema locale
   BBox3d b3SurfBox ;
   tmSurf.GetLocalBBox( b3SurfBox) ;
   // Se la distanza di sicurezza è significativa, espando il box
   if ( dSafeDist > EPS_SMALL )
      b3SurfBox.Expand( dSafeDist) ;
   // Box dello Zmap nel suo sistema locale
   BBox3d b3ZmapBox( ORIG, Point3d( m_nNx[0] * m_dStep, m_nNy[0] * m_dStep, m_dMaxZ[0])) ;
   b3ZmapBox.ToGlob( m_MapFrame) ;
   // Box intersezione: se non c'è intersezione ho finito.
   BBox3d b3IntBox ;
   if ( ! b3ZmapBox.FindIntersection( b3SurfBox, b3IntBox))
      return true ;
   // Recupero i triangoli della superficie che cadono nel box intersezione.
   INTVECTOR vTriaIndex ;
   tmSurf.GetAllTriaOverlapBox( b3IntBox, vTriaIndex) ;
   // Non è richeista precisione
   if ( ! bPrecise) {
      // Ciclo sui triangoli che cadono nel box
      for ( int nT : vTriaIndex) {
         Triangle3d trTria ;
         if ( ! ( tmSurf.GetTriangle( nT, trTria)  &&  trTria.Validate()))
            continue ;
         BBox3d b3TriaBox ;
         trTria.GetLocalBBox( b3TriaBox) ;
         // Se è necessario, espando il box di una costante additiva pari alla distanza di sicurezza.
         if ( dSafeDist > EPS_SMALL)
            b3TriaBox.Expand( dSafeDist) ;
         // Porto il bounding-box del triangolo nel sistema intrinseco.
         b3TriaBox.ToLoc( m_MapFrame) ;
         // Limiti su indici
         int nStI = Clamp( int( b3TriaBox.GetMin().x / m_dStep), 0, m_nNx[0] - 1) ;
         int nEnI = Clamp( int( b3TriaBox.GetMax().x / m_dStep), 0, m_nNx[0] - 1) ;
         int nStJ = Clamp( int( b3TriaBox.GetMin().y / m_dStep), 0, m_nNy[0] - 1) ;
         int nEnJ = Clamp( int( b3TriaBox.GetMax().y / m_dStep), 0, m_nNy[0] - 1) ;
         // Ciclo di intersezione dei dexel con il toro (nel riferimento intrinseco)
         for ( int i = nStI ; i <= nEnI ; ++ i) {
            for ( int j = nStJ ; j <= nEnJ ; ++ j) {
               int nPos = j * m_nNx[0] + i ;
               int nSize = int( m_Values[0][nPos].size()) ;
               if ( nSize == 0)
                  continue ;
               for ( int k = 0 ; k < 5 ; ++ k) {
                  Point3d ptLineSt = ORIG + ( i + 0.5) * m_dStep * X_AX + ( j + 0.5) * m_dStep * Y_AX ;
                  if ( k == 0)
                     ;
                  else if ( k == 1)
                     ptLineSt += - 0.4 * m_dStep * X_AX - 0.4 * m_dStep * Y_AX ;
                  else if ( k == 2)
                     ptLineSt += + 0.4 * m_dStep * X_AX - 0.4 * m_dStep * Y_AX ;
                  else if ( k == 3)
                     ptLineSt += + 0.4 * m_dStep * X_AX + 0.4 * m_dStep * Y_AX ;
                  else if ( k == 4)
                     ptLineSt += - 0.4 * m_dStep * X_AX + 0.4 * m_dStep * Y_AX ;
                  Vector3d vtLineDir = Z_AX ;
                  // Porto punto iniziale e vettore nel sistema locale
                  ptLineSt.ToGlob( m_MapFrame) ;
                  vtLineDir.ToGlob( m_MapFrame) ;
                  for ( int nIndex = 0 ; nIndex < nSize ; nIndex += 1) {
                     // Segmento ormai nel sistema locale
                     Point3d ptSegSt = ptLineSt + m_Values[0][nPos][nIndex].dMin * vtLineDir ;
                     double dSegLen = m_Values[0][nPos][nIndex].dMax - m_Values[0][nPos][nIndex].dMin ;
                     // Copia del triangolo che può essere evenutalmente traslata
                     Triangle3d trNewTria = trTria ;
                     // Se la distanza di sicurezza è significativa
                     if ( dSafeDist > EPS_SMALL) {
                        // Valuto sfere nei vertici e cilindri lungo gli edge.
                        for ( int nV = 0 ; nV < 3 ; ++ nV) {
                           Point3d ptVertP = trTria.GetP( nV) ;
                           Vector3d vtEdgeV = trTria.GetP( ( nV + 1)) - ptVertP ;
                           double dEdgeLen = vtEdgeV.Len() ;
                           vtEdgeV /= dEdgeLen ;
                           double dU1, dU2 ;
                           int nIntersType = SegmentSphere( ptSegSt, vtLineDir, dSegLen, ptVertP, dSafeDist, dU1, dU2) ;
                           if ( nIntersType != LinCompSphereIntersType::S_NO_INTERS)
                              return false ;
                           nIntersType = IntersSegmentCylinder( ptSegSt, vtLineDir, dSegLen, ptVertP, vtEdgeV,
                                                                 dSafeDist, dEdgeLen, dU1, dU2) ;
                           if ( nIntersType != LinCompCCIntersType::CC_NO_INTERS )
                              return false ;
                        }
                        // Traslo il triangolo.
                        trNewTria.Translate( dSafeDist * trTria.GetN()) ;
                     }
                     // Intersezione segento triangolo
                     Point3d ptInt, ptInt2 ;
                     int nIntersType = IntersLineTria( ptSegSt, vtLineDir, dSegLen, trNewTria, ptInt, ptInt2) ;
                     // Collisione
                     if ( nIntersType != IntLineTriaType::ILTT_NO)
                        return false ;
                  }
               }
            }
         }
      }
   }
   // È richiesta precisione
   else {
      // Ciclo sui triangoli che cadono nel box
      for ( int nT : vTriaIndex) {
         Triangle3d trTria ;
         if ( ! ( tmSurf.GetTriangle( nT, trTria)  &&  trTria.Validate()))
            continue ;
         BBox3d b3TriaBox ;
         trTria.GetLocalBBox( b3TriaBox) ;
         // Se è necessario, espando il box di una costante additiva pari alla distanza di sicurezza.
         if ( dSafeDist > EPS_SMALL)
            b3TriaBox.Expand( dSafeDist) ;
         // Porto il bounding-box del triangolo nel sistema intrinseco.
         b3TriaBox.ToLoc( m_MapFrame) ;
         // Ciclo sulle mappe
         for ( int nMap = 0 ; nMap < m_nMapNum ; ++ nMap) {
            Point3d ptInfIntBox = b3TriaBox.GetMin() ;
            Point3d ptSupIntBox = b3TriaBox.GetMax() ;
            // Dal sistema intrinseco al sistema griglia (per la prima griglia coincidono).
            if ( nMap == 1) {
               swap( ptInfIntBox.x, ptInfIntBox.z) ;
               swap( ptInfIntBox.x, ptInfIntBox.y) ;
               swap( ptSupIntBox.x, ptSupIntBox.z) ;
               swap( ptSupIntBox.x, ptSupIntBox.y) ;
            }
            else if (nMap == 2) {
               swap( ptInfIntBox.y, ptInfIntBox.z) ;
               swap( ptInfIntBox.x, ptInfIntBox.y) ;
               swap( ptSupIntBox.y, ptSupIntBox.z) ;
               swap( ptSupIntBox.x, ptSupIntBox.y) ;
            }
            // Limiti su indici
            int nStI = Clamp( int( ptInfIntBox.x / m_dStep), 0, m_nNx[nMap] - 1) ;
            int nEnI = Clamp( int( ptSupIntBox.x / m_dStep), 0, m_nNx[nMap] - 1) ;
            int nStJ = Clamp( int( ptInfIntBox.y / m_dStep), 0, m_nNy[nMap] - 1) ;
            int nEnJ = Clamp( int( ptSupIntBox.y / m_dStep), 0, m_nNy[nMap] - 1) ;
            for ( int nDex = 0 ; nDex < int( m_Values[nMap].size()) ; ++ nDex) {
               // Indici del dexel
               int nI = nDex % m_nNx[nMap] ;
               int nJ = nDex / m_nNx[nMap] ;
               // Se fuori dalla regione ammissibile salto l'iterazione
               if ( nI < nStI  ||  nI > nEnI  ||  nJ < nStJ  ||  nJ > nEnJ)
                  continue ;
               // Posizione del dexel
               double dX = ( nI + 0.5) * m_dStep ;
               double dY = ( nJ + 0.5) * m_dStep ;
               Point3d ptLineSt( dX, dY, 0.) ;
               Vector3d vtLineDir( 0., 0., 1.) ;
               // Dal sistema griglia al sistema intrinseco (per la prima griglia coincidono).
               if ( nMap == 1) {
                  swap( ptLineSt.x, ptLineSt.y) ;
                  swap( ptLineSt.x, ptLineSt.z) ;
                  swap( vtLineDir.x, vtLineDir.y) ;
                  swap( vtLineDir.x, vtLineDir.z) ;
               }
               else if ( nMap == 2) {
                  swap( ptLineSt.x, ptLineSt.y) ;
                  swap( ptLineSt.y, ptLineSt.z) ;
                  swap( vtLineDir.x, vtLineDir.y) ;
                  swap( vtLineDir.y, vtLineDir.z) ;
               }
               // Porto punto iniziale e vettore nel sistema locale
               ptLineSt.ToGlob( m_MapFrame) ;
               vtLineDir.ToGlob( m_MapFrame) ;
               // Ciclo sui segmenti del dexel.
               for ( int nInt = 0 ; nInt < int( m_Values[nMap][nDex].size()) ; ++ nInt) {
                  // Segmento ormai nel sistema locale
                  Point3d ptSegSt = ptLineSt + m_Values[nMap][nDex][nInt].dMin * vtLineDir ;
                  double dSegLen = m_Values[nMap][nDex][nInt].dMax - m_Values[nMap][nDex][nInt].dMin ;
                  // Copia del triangolo che può essere evenutalmente traslata
                  Triangle3d trNewTria = trTria ;
                  // Se la distanza di sicurezza è significativa
                  if ( dSafeDist > EPS_SMALL) {
                     // Valuto sfere nei vertici e cilindri lungo gli edge.
                     for ( int nV = 0 ; nV < 3 ; ++ nV) {
                        Point3d ptVertP = trTria.GetP( nV) ;
                        Vector3d vtEdgeV = trTria.GetP( ( nV + 1)) - ptVertP ;
                        double dEdgeLen = vtEdgeV.Len() ;
                        vtEdgeV /= dEdgeLen ;
                        double dU1, dU2 ;
                        int nIntersType = SegmentSphere( ptSegSt, vtLineDir, dSegLen, ptVertP, dSafeDist, dU1, dU2) ;
                        if ( nIntersType != LinCompSphereIntersType::S_NO_INTERS)
                           return false ;
                        nIntersType = IntersSegmentCylinder( ptSegSt, vtLineDir, dSegLen, ptVertP, vtEdgeV,
                           dSafeDist, dEdgeLen, dU1, dU2) ;
                        if ( nIntersType != LinCompCCIntersType::CC_NO_INTERS )
                           return false ;
                     }
                     // Traslo il triangolo.
                     trNewTria.Translate( dSafeDist * trTria.GetN()) ;
                  }
                  // Intersezione segento triangolo
                  Point3d ptInt, ptInt2 ;
                  int nIntersType = IntersLineTria( ptSegSt, vtLineDir, dSegLen, trNewTria, ptInt, ptInt2) ;
                  // Collisione
                  if ( nIntersType != IntLineTriaType::ILTT_NO)
                     return false ;
               }
            }
         }
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
// Riferimento con origine nel centro della base e asse di simmetria coincidente con l'asse Z.
// La funzione restituisce true in caso di intersezione, false altrimenti.
//----------------------------------------------------------------------------
bool
VolZmap::IntersLineCylinder( const Point3d& ptLineSt, const Vector3d& vtLineDir,
                             const Frame3d& CylFrame, double dH, double dRad, bool bTapLow, bool bTapUp,
                             Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2) const
{
  // Porto la linea nel riferimento del cilindro
   Point3d ptP = GetToLoc( ptLineSt, CylFrame) ;
   Vector3d vtV = GetToLoc( vtLineDir, CylFrame) ;

  // Determino le eventuali intersezioni con le due basi a quota minima e massima (solo se linea non parallela ad esse)
   int nBasInt = 0 ;
   if ( abs( vtV.z) > EPS_ZERO) {
      // le linee tangenti al cilindro non sono considerate intersecanti
      double EpsRad = ( vtV.IsZeroXY() ? - EPS_SMALL : EPS_SMALL) ;
      ptInt1 = ptP + ( ( 0 - ptP.z) / vtV.z) * vtV ;
      if ( ptInt1.x * ptInt1.x + ptInt1.y * ptInt1.y < dRad * dRad + 2 * dRad * EpsRad) {
         nBasInt += 1 ;
         vtN1 = Z_AX ;
      }
      ptInt2 = ptP + ( ( dH - ptP.z) / vtV.z) * vtV ;  
      if ( ptInt2.x * ptInt2.x + ptInt2.y * ptInt2.y < dRad * dRad + 2 * dRad * EpsRad) {
         nBasInt += 2 ;
         vtN2 = - Z_AX ;
      }
   }

  // Se la linea interseca entrambe le basi, si sono trovate le due intersezioni
   if ( nBasInt == 3) {
     // Porto i punti e i versori nel riferimento globale
      ptInt1.ToGlob( CylFrame) ;
      vtN1.ToGlob( CylFrame) ;
      ptInt2.ToGlob( CylFrame) ;
      vtN2.ToGlob( CylFrame) ;
     // Trovate intersezioni
      return true ;
   }
   
  // Determino le intersezioni con la superficie laterale del cilindro
   DBLVECTOR vdCoeff{ ptP.x * ptP.x + ptP.y * ptP.y - dRad * dRad,
                      2 * ( ptP.x * vtV.x + ptP.y * vtV.y),
                      vtV.x * vtV.x + vtV.y * vtV.y} ;
   DBLVECTOR vdRoots ;
   int nRoot = PolynomialRoots( 2, vdCoeff, vdRoots) ;

  // Epsilon per piani di tappo
   double dEpsLow = ( bTapLow ? - EPS_SMALL : EPS_SMALL) ;
   double dEpsUp = ( bTapUp ? EPS_SMALL : - EPS_SMALL) ;

  // Elimino le soluzioni cha danno intersezioni fuori dai limiti in Z del cilindro
   if ( nRoot == 2) {
      double dIntZ2 = ptP.z + vdRoots[1] * vtV.z ;
      if ( dIntZ2 < 0 + dEpsLow  ||  dIntZ2 > dH + dEpsUp)
         -- nRoot ;
   }
   if ( nRoot >= 1) {
      double dIntZ1 = ptP.z + vdRoots[0] * vtV.z ;
      if ( dIntZ1 < 0 + dEpsLow  ||  dIntZ1 > dH + dEpsUp) {
         if ( nRoot == 2)
            vdRoots[0] = vdRoots[1] ;
         -- nRoot ;
      }
   }

  // Due soluzioni: la retta interseca due volte la superficie laterale
   if ( nRoot == 2) {
     // Punti di intersezione con la superficie del cilindro
      ptInt1 = ptP + vdRoots[0] * vtV ;
      ptInt2 = ptP + vdRoots[1] * vtV ;
     // Determino le normali
      vtN1.Set( -ptInt1.x, -ptInt1.y, 0) ;
      vtN1.Normalize() ;
      vtN2.Set( -ptInt2.x, -ptInt2.y, 0) ;
      vtN2.Normalize() ;
     // Porto i punti e i versori nel riferimento globale
      ptInt1.ToGlob( CylFrame) ;
      vtN1.ToGlob( CylFrame) ;
      ptInt2.ToGlob( CylFrame) ;
      vtN2.ToGlob( CylFrame) ;
     // Trovate intersezioni
      return true ;
   }
   
  // Una soluzione : la retta interseca la superficie laterale e un piano
   else if ( nRoot == 1) {
     // Se piano superiore
      if ( nBasInt == 2) {
        // Punto di intersezione
         ptInt1 = ptP + vdRoots[0] * vtV ;
        // Normale alla superficie del cilindro verso l'interno
         vtN1.Set( -ptInt1.x, -ptInt1.y, 0) ;
         vtN1.Normalize() ;
      }
     // altrimenti piano inferiore
      else if ( nBasInt == 1) {
        // Punto di intersezione
         ptInt2 = ptP + vdRoots[0] * vtV ;
        // Normale alla superficie del cilindro verso l'interno
         vtN2.Set( -ptInt2.x, -ptInt2.y, 0) ;
         vtN2.Normalize() ;
      }
     // altrimenti niente
      else
         return false ;
     // Porto i punti e i versori nel riferimento globale
      ptInt1.ToGlob( CylFrame) ;
      vtN1.ToGlob( CylFrame) ;
      ptInt2.ToGlob( CylFrame) ;
      vtN2.ToGlob( CylFrame) ;
     // Trovate intersezioni
      return true ;
   }

  // Nessuna soluzione : nessuna intersezione 
   else
      return false ;
}

//----------------------------------------------------------------------------
// Riferimento con origine nel vertice del cono e asse di simmetria coincidente con l'asse Z.
// La funzione restituisce true in caso di intersezione, false altrimenti.
//----------------------------------------------------------------------------
bool
VolZmap::IntersLineConus( const Point3d& ptLineSt, const Vector3d& vtLineDir,
                          const Frame3d& ConusFrame, double dTan, double dMinH, double dMaxH, bool bTapLow, bool bTapUp,
                          Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2) const
{
  // Porto la linea nel riferimento del cono
   Point3d ptP = GetToLoc( ptLineSt, ConusFrame) ;
   Vector3d vtV = GetToLoc( vtLineDir, ConusFrame) ;

  // Raggi delle due basi
   double dMinRad = dTan * dMinH ;
   double dMaxRad = dTan * dMaxH ;

  // Epsilon per piani di tappo
   double dEpsLow = ( bTapLow ? - EPS_SMALL : EPS_SMALL) ;
   double dEpsUp = ( bTapUp ? EPS_SMALL : - EPS_SMALL) ;

  // Determino le eventuali intersezioni con le due basi a quota minima e massima (solo se linea non parallela ad esse)
   int nBasInt = 0 ;
   if ( abs( vtV.z) > EPS_ZERO) {
      ptInt1 = ptP + ( ( dMinH - ptP.z) / vtV.z) * vtV ;
      if ( ptInt1.x * ptInt1.x + ptInt1.y * ptInt1.y < dMinRad * dMinRad + 2 * dMinRad * dTan * dEpsLow) {
         nBasInt += 1 ;
         vtN1 = Z_AX ;
      }
      ptInt2 = ptP + ( ( dMaxH - ptP.z) / vtV.z) * vtV ;  
      if ( ptInt2.x * ptInt2.x + ptInt2.y * ptInt2.y < dMaxRad * dMaxRad + 2 * dMaxRad * dTan * dEpsUp) {
         nBasInt += 2 ;
         vtN2 = - Z_AX ;
      }
   }

  // Se la linea interseca entrambe le basi, si sono trovate le due intersezioni
   if ( nBasInt == 3) {
     // Porto i punti e i versori nel riferimento globale
      ptInt1.ToGlob( ConusFrame) ;
      vtN1.ToGlob( ConusFrame) ;
      ptInt2.ToGlob( ConusFrame) ;
      vtN2.ToGlob( ConusFrame) ;
     // Trovate intersezioni
      return true ;
   }

  // Determino le intersezioni con la superficie laterale del cono
   double dSqTan = dTan * dTan ;
   DBLVECTOR vdCoeff{ dSqTan * ptP.z * ptP.z - ptP.x * ptP.x - ptP.y * ptP.y,
                      2 * ( dSqTan * ptP.z * vtV.z - ptP.x * vtV.x - ptP.y * vtV.y),
                      dSqTan * vtV.z * vtV.z - vtV.x * vtV.x - vtV.y * vtV.y} ;
   DBLVECTOR vdRoots ;
   int nRoot = PolynomialRoots( 2, vdCoeff, vdRoots) ;

  // Elimino le soluzioni cha danno intersezioni fuori dai limiti in Z del tronco
   if ( nRoot == 2) {
      double dIntZ2 = ptP.z + vdRoots[1] * vtV.z ;
      if ( dIntZ2 < dMinH + dEpsLow  ||  dIntZ2 > dMaxH + dEpsUp)
         -- nRoot ;
   }
   if ( nRoot >= 1) {
      double dIntZ1 = ptP.z + vdRoots[0] * vtV.z ;
      if ( dIntZ1 < dMinH + dEpsLow  ||  dIntZ1 > dMaxH + dEpsUp) {
         if ( nRoot == 2)
            vdRoots[0] = vdRoots[1] ;
         -- nRoot ;
      }
   }

  // Due soluzioni: la retta interseca due volte la superficie laterale
   if ( nRoot == 2) {
     // Punti di intersezione con la superficie del cono
      ptInt1 = ptP + vdRoots[0] * vtV ;
      ptInt2 = ptP + vdRoots[1] * vtV ;
      if ( ptInt1.z > ptInt2.z) 
         swap( ptInt1, ptInt2) ;
     // Determino le normali
      vtN1.Set( -ptInt1.x, -ptInt1.y, ( ptInt1.x * ptInt1.x + ptInt1.y * ptInt1.y) / ptInt1.z) ;
      vtN1.Normalize() ;
      vtN2.Set( -ptInt2.x, -ptInt2.y, ( ptInt2.x * ptInt2.x + ptInt2.y * ptInt2.y) / ptInt2.z) ;
      vtN2.Normalize() ;
     // Porto i punti e i versori nel riferimento globale
      ptInt1.ToGlob( ConusFrame) ;
      vtN1.ToGlob( ConusFrame) ;
      ptInt2.ToGlob( ConusFrame) ;
      vtN2.ToGlob( ConusFrame) ;
     // Trovate intersezioni
      return true ;
   }
   
  // Una soluzione : la retta interseca la superficie laterale e un piano
   else if ( nRoot == 1) {
     // Se piano superiore
      if ( nBasInt == 2) {
        // Punto di intersezione
         ptInt1 = ptP + vdRoots[0] * vtV ;
        // Normale alla superficie del cono verso l'interno
         vtN1.Set( -ptInt1.x, -ptInt1.y, ( ptInt1.x * ptInt1.x + ptInt1.y * ptInt1.y) / ptInt1.z) ;
         vtN1.Normalize() ;
      }
     // altrimenti piano inferiore
      else if( nBasInt == 1)  {
        // Punto di intersezione
         ptInt2 = ptP + vdRoots[0] * vtV ;
        // Normale alla superficie del cono verso l'interno
         vtN2.Set( -ptInt2.x, -ptInt2.y, ( ptInt2.x * ptInt2.x + ptInt2.y * ptInt2.y) / ptInt2.z) ;
         vtN2.Normalize() ;
      }
     // altrimenti niente
      else
         return false ;
     // Porto i punti e i versori nel riferimento globale
      ptInt1.ToGlob( ConusFrame) ;
      vtN1.ToGlob( ConusFrame) ;
      ptInt2.ToGlob( ConusFrame) ;
      vtN2.ToGlob( ConusFrame) ;
     // Trovate intersezioni
      return true ;
   }

  // Nessuna soluzione : nessuna intersezione 
   else
      return false ;
}

//----------------------------------------------------------------------------
// L'origine del sistema di riferimento deve essere 
// nel centro della circonferenza di base, la cui traslazione obliqua
// genera il cilindro ellittico, e l'asse z deve essere l'asse 
// di simmetria di tale circonferenza.
// La funzione restituisce true in caso di intersezione, false altrimenti.
// NB: dSqRad è il quadrato del raggio della circonferenza la cui
// traslazione obliqua genera il cilindro ellittico, dLongMvLen e 
// dOrtMvLen sono rispettivamente le lunghezze delle proiezioni del 
// movimento su z e x del sistema di riferimento CircFrame. 
bool 
VolZmap::IntersLineEllipticalCylinder( const Point3d& ptLineSt, const Vector3d& vtLineDir,
                                       const Frame3d& CircFrame, double dRad, double dLongMvLen, double dOrtMvLen,
                                       bool bTapLow, bool bTapUp,
                                       Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2) const
{
  // Se il cilindrico ellittico degenera in un piano, non bisogna tagliare
   if ( abs( dLongMvLen) < EPS_SMALL)
      return false ;

  // Porto la linea nel riferimento del cilindro
   Point3d ptP = GetToLoc( ptLineSt, CircFrame) ;
   Vector3d vtV = GetToLoc( vtLineDir, CircFrame) ; 

  // Quadrato del raggio
   double dSqRad = dRad * dRad ;
  // Asse cilindro ellittico
   Vector3d vtAx( dOrtMvLen, 0, dLongMvLen) ;
   vtAx.Normalize() ;

  // Retta parallela all'asse del cilindro
   if ( AreSameOrOppositeVectorApprox( vtV, vtAx)) {
     // Interseco la retta con i piani delle circonferenze    
      Point3d ptOLsCirc( dOrtMvLen, 0, dLongMvLen) ;
      ptInt1 = ptP - ( ptP.z / vtV.z) * vtV ;
      ptInt2 = ptP - ( ( ptP.z - dLongMvLen) / vtV.z) * vtV ;
      double dSafeSqRad = dSqRad - 2 * dRad * EPS_SMALL ;   
      if ( ( ptInt1 - ORIG).SqLenXY() < dSafeSqRad  &&
           ( ptInt2 - ptOLsCirc).SqLenXY() < dSafeSqRad) {
         vtN1 = Z_AX ;
         vtN2 = - Z_AX ;
         ptInt1.ToGlob( CircFrame) ;
         ptInt2.ToGlob( CircFrame) ;
         vtN1.ToGlob( CircFrame) ;
         vtN2.ToGlob( CircFrame) ;
         return true ;
      }
      else 
         return false ; 
   }

  // Coefficiente angolare della retta di movimento nel 
  // piano ZX del sistema di riferimento del movimento e suo quadrato
   double dObCoef = dOrtMvLen / dLongMvLen ;
   double dSqCoef = dObCoef * dObCoef ;
  // Setto i coeficienti dell'equazione
   DBLVECTOR vdCoef(3) ;
   vdCoef[0] = dSqCoef * ptP.z * ptP.z + ptP.x * ptP.x + ptP.y * ptP.y 
               - 2 * dObCoef * ptP.z * ptP.x - dSqRad ;
   vdCoef[1] = 2 * ( dSqCoef * vtV.z * ptP.z + vtV.x * ptP.x + vtV.y * ptP.y 
               - dObCoef * ( vtV.z * ptP.x + vtV.x * ptP.z)) ;
   vdCoef[2] = dSqCoef * vtV.z * vtV.z + vtV.x * vtV.x + vtV.y * vtV.y
               - 2 * dObCoef * vtV.z * vtV.x ;
  // Numero di soluzioni
   DBLVECTOR vdRoots ;
   int nRoot = PolynomialRoots( 2, vdCoef, vdRoots) ;
  
  // L'equazione ammette o due soluzioni (eventualmente
  // coincidenti) oppure nessuna o infinite se la la retta
  // appartiene alla superficie         

  // Se numero di soluzioni diverso da due le eventuali intersezioni sono già state trovate
   if ( nRoot != 2)
      return false ;

  // Due soluzioni trovate
  // Vettore di movimento
   Vector3d vtMv( dOrtMvLen, 0, dLongMvLen) ;
  // Punti di intersezione
   ptInt1 = ptP + vdRoots[0] * vtV ;
   ptInt2 = ptP + vdRoots[1] * vtV ;
  // Simmetria del problema
   if ( ptInt1.z > ptInt2.z)        
      swap( ptInt1, ptInt2) ;
  // Determino le normali alla superficie nei punti d'intersezione
   // Vettore 1 
   Vector3d vtTest1 = ( ptInt1 - ORIG) - ( ptInt1 - ORIG) * vtAx * vtAx ;
   double dX0_1 = ( vtTest1.x > 0 ? 1 : -1) * sqrt( max( dSqRad - ptInt1.y * ptInt1.y, 0.)) ;  
   Vector3d vtCirc1( - dX0_1, - ptInt1.y, 0) ;
   Vector3d vtTan1( vtCirc1.y, - vtCirc1.x, 0) ;
   Vector3d vtCross1 = vtTan1 ^ vtMv ;
   vtN1 = ( vtCross1 * vtCirc1 > - EPS_ZERO ? vtCross1 : - vtCross1) ;
   vtN1.Normalize() ;
   // Vettore 2
   Vector3d vtTest2 = ( ptInt2 - ORIG) - ( ptInt2 - ORIG) * vtAx * vtAx ;
   double dX0_2 = ( vtTest2.x > 0 ? 1 : -1) * sqrt( max( dSqRad - ptInt2.y * ptInt2.y, 0.)) ;  
   Vector3d vtCirc2( - dX0_2, - ptInt2.y, 0) ;
   Vector3d vtTan2( vtCirc2.y, - vtCirc2.x, 0) ;
   Vector3d vtCross2 = vtTan2 ^ vtMv ;
   vtN2 = ( vtCross2 * vtCirc2 > - EPS_ZERO ? vtCross2 : - vtCross2) ;
   vtN2.Normalize() ;

  // Flag per i tappi
   double dEpsLow = ( bTapLow ? - EPS_SMALL : EPS_SMALL) ;
   double dEpsUp = ( bTapUp ?  EPS_SMALL : - EPS_SMALL) ;

  // Studio le soluzioni: se una è fuori dalla regione 
  // ammissibile, vuol dire che la retta esce da un tappo.
   if ( ptInt1.z < dLongMvLen + dEpsUp) {     
      if (  ptInt1.z > + dEpsLow) {
        // ptInt2 è sul tappo
         if ( ptInt2.z > dLongMvLen + dEpsUp) {
            ptInt2 = ptP + ( ( dLongMvLen - ptP.z) / vtV.z) * vtV ;  
            vtN2 = - Z_AX ;
         }
      }
      else {
        // Entrambe le soluzioni sono su un tappo
         if ( ptInt2.z > dLongMvLen + dEpsUp) {        
            ptInt1 = ptP - ( ptP.z / vtV.z) * vtV ;
            ptInt2 = ptP + ( ( dLongMvLen - ptP.z) / vtV.z) * vtV ;
            vtN1.Set( 0, 0, 1)  ;
            vtN2.Set( 0, 0, -1) ;
            if ( ptInt1.x * ptInt1.x + ptInt1.y * ptInt1.y > dSqRad &&
                 ptInt2.x * ptInt2.x + ptInt2.y * ptInt2.y > dSqRad)
                 return false ;
         }
        // La prima soluzione è sul tappo
         else if ( ptInt2.z > dEpsLow) {            
            ptInt1 = ptP - ( ptP.z / vtV.z) * vtV ;
            vtN1.Set( 0, 0, 1) ;
         }
         else 
            return false ;
      }
   }
   else 
      return false ;

  // Se il secondo punto è sceso in Z sotto il primo o il primo è salito sopra il secondo vuol
  // dire che l'intersezione è al limite del cilindro ellittico; non c'è reale intersezione. 
   if ( ptInt1.z > ptInt2.z)
      return false ;

  // Porto i punti e i versori nel riferimento globale
   ptInt1.ToGlob( CircFrame) ;
   vtN1.ToGlob( CircFrame) ;
   ptInt2.ToGlob( CircFrame) ;
   vtN2.ToGlob( CircFrame) ; 
       
   return true ;
}

//----------------------------------------------------------------------------
bool 
VolZmap::IntersLineMyPolyhedron( const Point3d& ptLineSt, const Vector3d& vtLineDir,
                                 const Frame3d& PolyFrame, double dLenX, double dLenY, double dLenZ, double dDeltaZ,
                                 Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2) const
{
  // Controllo sulle dimensioni lineari affinché sia valido il poliedro
   if ( dLenX <= EPS_SMALL  ||  dLenY <= EPS_SMALL  ||  dLenZ <= EPS_SMALL)
      return false ;
  
  // Porto la linea nel riferimento del poliedro
   Point3d ptP = GetToLoc( ptLineSt, PolyFrame) ;
   Vector3d vtV = GetToLoc( vtLineDir, PolyFrame) ;

  // Facce 1 e 2 parallele a XY
  // Facce 3 e 4 parallele a XZ
  // Facce 5 e 6 oblique
   Point3d ptFacet135( 0, dLenY / 2, 0) ;
   Point3d ptFacet246( dLenX, - dLenY / 2, dLenZ + dDeltaZ) ;

  // Vettori per descrizione piani obliqui
   Vector3d vtFacet5 = ptFacet135 - ptP ;
   Vector3d vtFacet6 = ptFacet246 - ptP ;
   Vector3d vtOb( - dDeltaZ, 0, dLenX) ;
   vtOb.Normalize() ;

  // Controlli affinché non vengano tagliati dexel a filo
  // con il passaggio dell'utensile:
  // Controllo sulle facce 1 e 2
   if ( abs( vtV.y) < EPS_ZERO  &&  abs( ptP.y) > ptFacet135.y)   
      return false ;
  // Controllo sulle facce 3 e 4
   if ( abs( vtV.x) < EPS_ZERO  &&  ( ptP.x < ptFacet135.x  ||  ptP.x > ptFacet246.x))
      return false ;
  // Controllo sulle facce 5 e 6 
   double dP1 = abs ( ( ptFacet135 - ptP) * vtOb) ;
   double dP2 = abs ( ( ptFacet246 - ptP) * vtOb) ;
   if ( abs( vtV * vtOb) < EPS_ZERO  &&  ( dP1 < EPS_SMALL  ||  dP2  < EPS_SMALL))
      return false ;

  // Punti notevoli
   Point3d ptI1 = ptP + ( ( ptFacet135.y - ptP.y) / vtV.y) * vtV ;
   Point3d ptI2 = ptP + ( ( ptFacet246.y - ptP.y) / vtV.y) * vtV ;
   Point3d ptI3 = ptP + ( ( ptFacet135.x - ptP.x) / vtV.x) * vtV ;
   Point3d ptI4 = ptP + ( ( ptFacet246.x - ptP.x) / vtV.x) * vtV ;
   Point3d ptI5 = ptP + ( ( vtFacet5 * vtOb) / ( vtV * vtOb)) * vtV ;
   Point3d ptI6 = ptP + ( ( vtFacet6 * vtOb) / ( vtV * vtOb)) * vtV ;

  // Ricerca intersezioni con le facce 
   int nIntNum = 0 ;

  // Intersezione con la prima faccia
   if ( ptI1.x >= 0  &&  ptI1.x <= dLenX  && 
        ptI1.z * dLenX >= dDeltaZ * ptI1.x  &&  ( ptI1.z - dLenZ) * dLenX <= dDeltaZ * ptI1.x) {
      ptInt1 = ptI1 ;
      vtN1 = - Y_AX ;
      ++ nIntNum ;
   }
  
  // Intersezione con la seconda faccia
   if ( ptI2.x >= 0  &&  ptI2.x <= dLenX  && 
        ptI2.z * dLenX >= dDeltaZ * ptI2.x  &&  ( ptI2.z - dLenZ) * dLenX <= dDeltaZ * ptI2.x) {
      if ( nIntNum == 0) {
         ptInt1 = ptI2 ;
         vtN1 = Y_AX ;
         ++ nIntNum ;
      }
      else if ( ( ptInt1 - ptI2).SqLen() > SQ_EPS_SMALL) {
         ptInt2 = ptI2 ;
         vtN2 = Y_AX ;
         ++ nIntNum ;
      }
   }

  // Intersezione con la terza faccia
   if ( nIntNum < 2  &&
        ptI3.z >= 0  &&  ptI3.z <= dLenZ  && 
        ptI3.y >= - ptFacet135.y  &&  ptI3.y <= ptFacet135.y) {
      if ( nIntNum == 0) {
         ptInt1 = ptI3 ;
         vtN1 = X_AX ;
         ++ nIntNum ;
      }
      else if ( ( ptInt1 - ptI3).SqLen() > SQ_EPS_SMALL) {
         ptInt2 = ptI3 ;
         vtN2 = X_AX ;
         ++ nIntNum ;
      }
   }

  // Intersezione con la quarta faccia
   if ( nIntNum < 2  &&
        ptI4.z >= dDeltaZ  &&  ptI4.z <= dLenZ + dDeltaZ  && 
        ptI4.y >= - ptFacet135.y  &&  ptI4.y <= ptFacet135.y) {
      if ( nIntNum == 0) {
         ptInt1 = ptI4 ;
         vtN1 = - X_AX ;
         ++ nIntNum ;
      }
      else if ( ( ptInt1 - ptI4).SqLen() > SQ_EPS_SMALL) {
         ptInt2 = ptI4 ;
         vtN2 = - X_AX ;
         ++ nIntNum ;
      }
   }

  // Intersezione con la quinta faccia
   if ( nIntNum < 2  &&
        ptI5.x >= 0  &&  ptI5.x <= dLenX  && 
        ptI5.y >= - ptFacet135.y  &&  ptI5.y <= ptFacet135.y) {
      if ( nIntNum == 0) {
         ptInt1 = ptI5 ;
         vtN1 = vtOb ;
         ++ nIntNum ;
      }
      else if ( ( ptInt1 - ptI5).SqLen() > SQ_EPS_SMALL) {
         ptInt2 = ptI5 ;
         vtN2 = vtOb ;
         ++ nIntNum ;
      }  
   }

  // Intersezione con la sesta faccia
   if ( nIntNum < 2  &&
        ptI6.x >= 0  &&  ptI6.x <= dLenX  && 
        ptI6.y >= - ptFacet135.y  &&  ptI6.y <= ptFacet135.y) {
      if ( nIntNum == 0) {
         ptInt1 = ptI6;
         vtN1 = - vtOb ;
         ++ nIntNum ;
      }
      else if ( ( ptInt1 - ptI6).SqLen() > SQ_EPS_SMALL) {
         ptInt2 = ptI6;
         vtN2 = - vtOb ;
         ++ nIntNum ;
      } 
   }
   
   if ( nIntNum != 2)
      return false ;

  // Porto i punti e i versori nel riferimento globale
   ptInt1.ToGlob( PolyFrame) ;
   ptInt2.ToGlob( PolyFrame) ;
   vtN1.ToGlob( PolyFrame) ;
   vtN2.ToGlob( PolyFrame) ;
   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::IntersLineTruncatedPyramid( const Point3d& ptLineSt, const Vector3d& vtLineDir,
                                     const Frame3d& frTruncPyramFrame, double dSegMin, double dSegMax, double dHeight,
                                     Point3d& ptInt1, Vector3d& vtN1, Point3d& ptInt2, Vector3d& vtN2) const
{
  // Controllo sulle dimensioni lineari affinché sia valido il solido
   if ( dSegMin < EPS_SMALL  ||  dSegMax < EPS_SMALL  ||  dHeight < EPS_SMALL)
      return false ;

  // Porto la linea nel riferimento del solido
   Point3d ptP = GetToLoc( ptLineSt, frTruncPyramFrame) ;
   Vector3d vtV = GetToLoc( vtLineDir, frTruncPyramFrame) ;
 
  // Se la retta sta sopra o sotto il solido non vi può essere intersezione
   if ( abs( vtV.z) < EPS_ZERO  &&  ( ptP.z < EPS_SMALL  ||  ptP.z > dHeight + EPS_SMALL))
      return false ;

   double dHalfMax = 0.5 * dSegMax ;
   double dHalfMin = 0.5 * dSegMin ;

  // Cerco le intersezioni con i piani delle facce
   int nIntNum = 0 ;
  // Base maggiore
   Point3d ptIntPlaneMax = ptP + ( - ptP.z / vtV.z) * vtV ;
   if ( abs( ptIntPlaneMax.x) <= dHalfMax &&  abs( ptIntPlaneMax.y) <= dHalfMax) {
      ptInt1 = ptIntPlaneMax ;
      vtN1 = - Z_AX ;
      ++ nIntNum ;
   }
  // Base minore
   Point3d ptIntPlaneMin = ptP + ( ( dHeight - ptP.z) / vtV.z) * vtV ;
   if ( abs( ptIntPlaneMin.x) <= dHalfMin &&  abs( ptIntPlaneMin.y) <= dHalfMin) {
      if ( nIntNum == 0) {
         ptInt1 = ptIntPlaneMin ;
         vtN1 =  Z_AX ;
      }
      else {
         ptInt2 = ptIntPlaneMin ;
         vtN2 = Z_AX ;
      }
      ++ nIntNum ;
   }

  // Se la retta intersca i piani delle basi al loro interno, ho finito
   if ( nIntNum == 2) {
      if ( ( ptInt1 - ptP) * vtV > ( ptInt2 - ptP) * vtV) {
         swap( ptInt1, ptInt2) ;
         swap( vtN1, vtN2) ;
      }
      ptInt1.ToGlob( frTruncPyramFrame) ;
      ptInt2.ToGlob( frTruncPyramFrame) ;
      vtN1.ToGlob( frTruncPyramFrame) ;
      vtN2.ToGlob( frTruncPyramFrame) ;
      return true ;
   }

  // Altrimenti se le intersezioni sono fuori dalle basi e dallo stesso lato rispetto a queste ultime non c'è intersezione
   if ( nIntNum == 0  &&
        ( ( ptIntPlaneMax.x < - dHalfMax - EPS_SMALL  &&  ptIntPlaneMin.x < - dHalfMin - EPS_SMALL)  ||
          ( ptIntPlaneMax.x >   dHalfMax + EPS_SMALL  &&  ptIntPlaneMin.x >   dHalfMin + EPS_SMALL)  ||
          ( ptIntPlaneMax.y < - dHalfMax - EPS_SMALL  &&  ptIntPlaneMin.y < - dHalfMin - EPS_SMALL)  ||
          ( ptIntPlaneMax.y >   dHalfMax + EPS_SMALL  &&  ptIntPlaneMin.y >   dHalfMin + EPS_SMALL)))
      return false ;

  // Cerco le eventuali intersezioni con i piani delle facce laterali
   Point3d ptPlaneP( - dHalfMax, - dHalfMax, 0.) ;
   Vector3d vtPlaneN( 0., - dHeight, ( dHalfMax - dHalfMin)) ;
   vtPlaneN.Normalize() ;
   if ( nIntNum < 2  &&  abs( vtV * vtPlaneN) > EPS_ZERO) {
      if ( nIntNum == 0) {
         ptInt1 = ptP + ( ( ( ptPlaneP - ptP) * vtPlaneN) / ( vtV * vtPlaneN)) * vtV ;
         vtN1 = vtPlaneN ;
         if (   ptInt1.z >= 0  &&  ptInt1.z <= dHeight  &&
              ( ptInt1.x + dHalfMax) * dHeight >= ( dHalfMax - dHalfMin) * ptInt1.z  &&
              ( ptInt1.x - dHalfMax) * dHeight <= ( dHalfMin - dHalfMax) * ptInt1.z)
            ++ nIntNum ;
      }
      else {
         ptInt2 = ptP + ( ( ( ptPlaneP - ptP) * vtPlaneN) / ( vtV * vtPlaneN)) * vtV ;
         vtN2 = vtPlaneN ;
         if ( ( ptInt2 - ptInt2).SqLen() > SQ_EPS_SMALL  &&
                ptInt2.z >= 0  &&  ptInt1.z <= dHeight  &&
              ( ptInt2.x + dHalfMax) * dHeight >= ( dHalfMax - dHalfMin) * ptInt2.z  &&
              ( ptInt2.x - dHalfMax) * dHeight <= ( dHalfMin - dHalfMax) * ptInt2.z)
            ++ nIntNum ;
      }
   }
   vtPlaneN = Vector3d( - dHeight, 0., ( dHalfMax - dHalfMin)) ;
   vtPlaneN.Normalize() ;
   if ( nIntNum < 2  &&  abs( vtV * vtPlaneN) > EPS_ZERO) {
      if ( nIntNum == 0) {
         ptInt1 = ptP + ( ( ( ptPlaneP - ptP) * vtPlaneN) / ( vtV * vtPlaneN)) * vtV ;
         vtN1 = vtPlaneN ;
         if (   ptInt1.z >= 0  &&  ptInt1.z <= dHeight  &&
              ( ptInt1.y + dHalfMax) * dHeight >= ( dHalfMax - dHalfMin) * ptInt1.z  &&
              ( ptInt1.y - dHalfMax) * dHeight <= ( dHalfMin - dHalfMax) * ptInt1.z)
            ++ nIntNum ;
      }
      else {
         ptInt2 = ptP + ( ( ( ptPlaneP - ptP) * vtPlaneN) / ( vtV * vtPlaneN)) * vtV ;
         vtN2 = vtPlaneN ;
         if ( ( ptInt1 - ptInt2).SqLen() > SQ_EPS_SMALL  &&
                ptInt2.z >= 0  &&  ptInt2.z <= dHeight  &&
              ( ptInt2.y + dHalfMax) * dHeight >= ( dHalfMax - dHalfMin) * ptInt2.z  &&
              ( ptInt2.y - dHalfMax) * dHeight <= ( dHalfMin - dHalfMax) * ptInt2.z)
            ++ nIntNum ;
      }
   }
   ptPlaneP = Point3d( dHalfMax, dHalfMax, 0.) ;
   vtPlaneN = Vector3d( dHeight, 0., ( dHalfMax - dHalfMin)) ;
   vtPlaneN.Normalize() ;
   if ( nIntNum < 2  &&  abs( vtV * vtPlaneN) > EPS_ZERO) {
      if ( nIntNum == 0) {
         ptInt1 = ptP + ( ( ( ptPlaneP - ptP) * vtPlaneN) / ( vtV * vtPlaneN)) * vtV ;
         vtN1 = vtPlaneN ;
         if (   ptInt1.z >= 0  &&  ptInt1.z <= dHeight  &&
              ( ptInt1.y + dHalfMax) * dHeight >= ( dHalfMax - dHalfMin) * ptInt1.z  &&
              ( ptInt1.y - dHalfMax) * dHeight <= ( dHalfMin - dHalfMax) * ptInt1.z)
            ++ nIntNum ;
      }
      else {
         ptInt2 = ptP + ( ( ( ptPlaneP - ptP) * vtPlaneN) / ( vtV * vtPlaneN)) * vtV ;
         vtN2 = vtPlaneN ;
         if ( ( ptInt1 - ptInt2).SqLen() > SQ_EPS_SMALL  &&
                ptInt2.z >= 0  &&  ptInt2.z <= dHeight  &&
              ( ptInt2.y + dHalfMax) * dHeight >= ( dHalfMax - dHalfMin) * ptInt2.z  &&
              ( ptInt2.y - dHalfMax) * dHeight <= ( dHalfMin - dHalfMax) * ptInt2.z)
            ++ nIntNum ;
      }
   }
   vtPlaneN = Vector3d( 0., dHeight, ( dHalfMax - dHalfMin)) ;
   vtPlaneN.Normalize() ;
   if ( nIntNum < 2  &&  abs( vtV * vtPlaneN) > EPS_ZERO) {
      if ( nIntNum == 0) {
         ptInt1 = ptP + ( ( ( ptPlaneP - ptP) * vtPlaneN) / ( vtV * vtPlaneN)) * vtV ;
         vtN1 = vtPlaneN ;
         if (   ptInt1.z >= 0 && ptInt1.z <= dHeight  &&
              ( ptInt1.x + dHalfMax) * dHeight >= ( dHalfMax - dHalfMin) * ptInt1.z &&
              ( ptInt1.x - dHalfMax) * dHeight <= ( dHalfMin - dHalfMax) * ptInt1.z)
            ++ nIntNum ;
      }
      else {
         ptInt2 = ptP + ( ( ( ptPlaneP - ptP) * vtPlaneN) / ( vtV * vtPlaneN)) * vtV ;
         vtN2 = vtPlaneN ;
         if ( ( ptInt1 - ptInt2).SqLen() > SQ_EPS_SMALL  &&
                ptInt2.z >= 0 && ptInt1.z <= dHeight  &&
              ( ptInt2.x + dHalfMax) * dHeight >= ( dHalfMax - dHalfMin) * ptInt2.z  &&
              ( ptInt2.x - dHalfMax) * dHeight <= ( dHalfMin - dHalfMax) * ptInt2.z)
            ++ nIntNum ;
      }
   }

  // La retta interseca il tronco di piramide
   if ( nIntNum == 2) {
      if ( ( ptInt1 - ptP) * vtV > ( ptInt2 - ptP) * vtV) {
         swap( ptInt1, ptInt2) ;
         swap( vtN1, vtN2) ;
      }
      ptInt1.ToGlob( frTruncPyramFrame) ;
      ptInt2.ToGlob( frTruncPyramFrame) ;
      vtN1.ToGlob( frTruncPyramFrame) ;
      vtN2.ToGlob( frTruncPyramFrame) ;
      return true ;
   }

   return false ;
}

//----------------------------------------------------------------------------
bool
VolZmap::GetVolume( double& dVol) const 
{
  // verifico lo stato
   if ( m_nStatus != OK)
      return false ;
  // Eseguo il calcolo della lunghezza totale degli spilloni
   double dLen = 0 ;
   for ( int nMap = 0 ; nMap < int( m_nMapNum) ; ++ nMap) {
      for ( int nDex = 0 ; nDex < int( m_Values[nMap].size()) ; ++ nDex) {
         for ( int nInt = 0 ; nInt < int( m_Values[nMap][nDex].size()) ; ++ nInt) {
            dLen += m_Values[nMap][nDex][nInt].dMax - m_Values[nMap][nDex][nInt].dMin ;
         }
      }
   }
  // Il volume si trova moltiplicando per l'area di ogni spillone e dividendo per il numero di mappe
   dVol = dLen * ( m_dStep * m_dStep) / m_nMapNum ;

   return true ;  
}

//----------------------------------------------------------------------------
bool
VolZmap::GetPartVolume( int nPart, double& dVol) const
{ 
  // verifico lo stato
   if ( m_nStatus != OK)
      return false ;
  // Se una sola mappa o il numero di componenti è indefinito, errore
   if (  m_nMapNum == 1  ||  m_nConnectedCompoCount == -1) 
      return false ;
  // Se la componente richiesta non esiste, errore
   if ( nPart < 0  ||  nPart > m_nConnectedCompoCount - 1)
      return false ;

  // Eseguo il calcolo della lunghezza totale degli spilloni
   double dLen = 0 ;
   for ( int nMap = 0 ; nMap < m_nMapNum ; ++ nMap) {
      for ( int nDex = 0 ; nDex < m_nDim[nMap] ; ++ nDex) {
         for ( int nInt = 0 ; nInt < int( m_Values[nMap][nDex].size()) ; ++ nInt) {
           // se il pezzo di spillone appartiene alla parte, ne aggiungo la lunghezza
            if ( m_Values[nMap][nDex][nInt].nCompo == nPart + 1) {
               dLen += m_Values[nMap][nDex][nInt].dMax - m_Values[nMap][nDex][nInt].dMin ;
            }
         }
      }
   }
  // Il volume si trova moltiplicando per l'area di ogni spillone e dividendo per il numero di mappe
   dVol = dLen * ( m_dStep * m_dStep) / 3 ;

   return true ;
}

//----------------------------------------------------------------------------
// La retta deve essere espressa nel sistema Zmap.
// Per riferimento viene restituito un vettore di intersezioni della retta con i triangoli della superficie del solido.
// Si restituisce false in caso di errore, true altrimenti.
bool
VolZmap::GetLineIntersection( const Point3d& ptP, const Vector3d& vtD, ILZIVECTOR& vIntersInfo) const
{
  // Direzione normalizzata
   Vector3d vtDir = vtD ;
   vtDir.Normalize() ;

  // Punto e vettore espressi nel riferimento intrinseco dello Zmap
   Point3d ptLocP = ptP ;
   Vector3d vtDirL = vtDir ;
   ptLocP.ToLoc( m_MapFrame) ;
   vtDirL.ToLoc( m_MapFrame) ;

  // Intersezione fra semiretta e BBox dello Zmap
   double dU1, dU2 ;
   bool bLineBBoxInters = IntersLineZMapBBox( ptLocP, vtDirL, dU1, dU2) && ( dU1 > 0  ||  dU2 > 0) ;

  // Semiretta esterna al box dello Zmap quindi esterna anche allo Zmap
   if ( ! bLineBBoxInters)
      return true ;

  // Lancio eventuale aggiornamento pendente della grafica
   if ( m_nMapNum == 1)
      UpdateSingleMapGraphics() ;
   else
      UpdateTripleMapGraphics() ;

  // Ciclo sui blocchi
   for ( int nB = 0 ; nB < m_nNumBlock ; ++ nB) {
     // Determino indici IJK del blocco; se non trovo il blocco, errore
      int nBlockIJK[3] ;
      if ( ! GetBlockIJKFromN( nB, nBlockIJK))
         return false ;
     // Costruisco il bounding-box del blocco; se non è possibile, errore
      BBox3d b3BlockBox ;
      if ( ! GetBlockBox( nBlockIJK, b3BlockBox))
         return false ;
     // Se c'è intersezione valuto tutti i voxel interni
      if ( TestIntersLineBox( ptLocP, vtDirL, b3BlockBox.GetMin(), b3BlockBox.GetMax())) {
         // Ciclo sui voxel del blocco.
         // Triangoli smooth
         for ( int nV = 0 ; nV < int( m_BlockSmoothTria[nB].size()) ; ++ nV) {
           // Box del voxel
            BBox3d b3Vox ;
            int nCurVoxIJK[3] = { m_BlockSmoothTria[nB][nV].i,
                                  m_BlockSmoothTria[nB][nV].j,
                                  m_BlockSmoothTria[nB][nV].k } ;
            GetVoxelBox( nCurVoxIJK[0], nCurVoxIJK[1], nCurVoxIJK[2], b3Vox) ;
           // Se non c'è intersezione col voxel, passo al successivo.
            if ( ! TestIntersLineBox( ptLocP, vtDirL, b3Vox.GetMin(), b3Vox.GetMax()))
               continue ;
            for ( int nT = 0 ; nT < int( m_BlockSmoothTria[nB][nV].vTria.size()) ; ++ nT) {
               Triangle3d trTria = m_BlockSmoothTria[nB][nV].vTria[nT] ;
               if ( ! trTria.Validate( true))
                  continue ;
               Point3d ptLineTria1, ptLineTria2 ;
              // Studio dell'intersezione della retta con il triangolo corrente
               int nIntType = IntersLineTria( ptP, vtDir, 1.5 * dU2, trTria, 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) {
                  int nNumVox ;  
                  GetVoxNFromIJK( nCurVoxIJK[0], nCurVoxIJK[1], nCurVoxIJK[2], nNumVox) ;
                  vIntersInfo.emplace_back( nIntType, ( ptLineTria1 - ptP) * vtDir,
                                            nNumVox, nB, ptLineTria1, trTria) ;
               }
              // altrimenti ci sono due intersezioni
               else {
                  int nNumVox ;
                  GetVoxNFromIJK( nCurVoxIJK[0], nCurVoxIJK[1], nCurVoxIJK[2], nNumVox) ;
                  double dP1 = ( ptLineTria1 - ptP) * vtDir ;
                  double dP2 = ( ptLineTria2 - ptP) * vtDir ;
                  vIntersInfo.emplace_back( nIntType, ( dP1 < dP2 ? dP1 : dP2), ( dP1 < dP2 ? dP2 : dP1),
                                            nNumVox, nB, ptLineTria1, ptLineTria2, trTria) ;
               }
            }
         }
         // Triangoli sharp interni al blocco
         for ( int nV = 0 ; nV < int( m_BlockSharpTria[nB].size()) ; ++ nV) {
            int nCurVoxIJK[3] = { m_BlockSharpTria[nB][nV].i,
                                  m_BlockSharpTria[nB][nV].j,
                                  m_BlockSharpTria[nB][nV].k } ;

            // Ciclo sulle componenti connesse
            for ( int nC = 0 ; nC < int( m_BlockSharpTria[nB][nV].vCompoTria.size()) ; ++ nC) {
               for ( int nT = 0 ; nT < int( m_BlockSharpTria[nB][nV].vCompoTria[nC].size()) ; ++ nT) {
                  Triangle3d trTria = m_BlockSharpTria[nB][nV].vCompoTria[nC][nT] ;
                  if ( ! trTria.Validate( true))
                     continue ;
                  Point3d ptLineTria1, ptLineTria2 ;
                 // Studio dell'intersezione della retta con il triangolo corrente
                  int nIntType = IntersLineTria( ptP, vtDir, 1.5 * dU2, trTria, 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) {
                     int nNumVox ;
                     GetVoxNFromIJK( nCurVoxIJK[0], nCurVoxIJK[1], nCurVoxIJK[2], nNumVox) ;
                     vIntersInfo.emplace_back( nIntType, ( ptLineTria1 - ptP) * vtDir,
                                               nNumVox, nB, ptLineTria1, trTria) ;
                  }
                  // altrimenti ci sono due intersezioni
                  else {
                     int nNumVox ;
                     GetVoxNFromIJK( nCurVoxIJK[0], nCurVoxIJK[1], nCurVoxIJK[2], nNumVox) ;
                     double dP1 = ( ptLineTria1 - ptP) * vtDir ;
                     double dP2 = ( ptLineTria2 - ptP) * vtDir ;
                     vIntersInfo.emplace_back( nIntType, ( dP1 < dP2 ? dP1 : dP2), ( dP1 < dP2 ? dP2 : dP1),
                                               nNumVox, nB, ptLineTria1, ptLineTria2, trTria) ;
                  }
               }
            }
         }
         // Triangoli grandi del blocco
         for ( int nT = 0 ; nT < int( m_BlockBigTria[nB].size()) ; ++ nT) {
            Triangle3d trTria = m_BlockBigTria[nB][nT] ;
            Point3d ptLineTria1, ptLineTria2 ;
           // Studio dell'intersezione della retta con il triangolo corrente
            int nIntType = IntersLineTria( ptP, vtDir, 1.5 * dU2, trTria, 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) {
               vIntersInfo.emplace_back( nIntType, ( ptLineTria1 - ptP) * vtDir,
                                         -1, nB, ptLineTria1, trTria) ;
            }
           // altrimenti ci sono due intersezioni
            else {
               double dP1 = ( ptLineTria1 - ptP) * vtDir ;
               double dP2 = ( ptLineTria2 - ptP) * vtDir ;
               vIntersInfo.emplace_back( nIntType, ( dP1 < dP2 ? dP1 : dP2), ( dP1 < dP2 ? dP2 : dP1),
                                         -1, nB, ptLineTria1, ptLineTria2, trTria) ;
            }
         }
      }
      // In ogni caso valuto i triangoli sharp fra blocchi
      for ( int nV = 0 ; nV < int( m_InterBlockSharpTria[nB].size()) ; ++ nV) {
         int nCurVoxIJK[3] = { m_InterBlockSharpTria[nB][nV].i,
                               m_InterBlockSharpTria[nB][nV].j,
                               m_InterBlockSharpTria[nB][nV].k } ;
         // Ciclo sulle componenti connesse
         for ( int nC = 0 ; nC < int( m_InterBlockSharpTria[nB][nV].vCompoTria.size()) ; ++ nC) {
            for ( int nT = 0 ; nT < int( m_InterBlockSharpTria[nB][nV].vCompoTria[nC].size()) ; ++ nT) {
               Triangle3d trTria = m_InterBlockSharpTria[nB][nV].vCompoTria[nC][nT] ;
               if ( ! trTria.Validate( true))
                  continue ;
               Point3d ptLineTria1, ptLineTria2 ;
              // Studio dell'intersezione della retta con il triangolo corrente
               int nIntType = IntersLineTria( ptP, vtDir, 1.5 * dU2, trTria, 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) {
                  int nNumVox ;
                  GetVoxNFromIJK( nCurVoxIJK[0], nCurVoxIJK[1], nCurVoxIJK[2], nNumVox) ;
                  vIntersInfo.emplace_back( nIntType, ( ptLineTria1 - ptP) * vtDir,
                                            nNumVox, nB, ptLineTria1, trTria) ;
               }
               // altrimenti ci sono due intersezioni
               else {
                  int nNumVox ;
                  GetVoxNFromIJK( nCurVoxIJK[0], nCurVoxIJK[1], nCurVoxIJK[2], nNumVox) ;
                  double dP1 = ( ptLineTria1 - ptP) * vtDir ;
                  double dP2 = ( ptLineTria2 - ptP) * vtDir ;
                  vIntersInfo.emplace_back( nIntType, ( dP1 < dP2 ? dP1 : dP2), ( dP1 < dP2 ? dP2 : dP1),
                                            nNumVox, nB, ptLineTria1, ptLineTria2, trTria) ;
               }
            }
         }
      }
   }

   // Ordino le intersezioni in base al parametro distanza con segno da ptP
   sort( vIntersInfo.begin(), vIntersInfo.end(),
         []( const IntLineZmapInfo& a, const IntLineZmapInfo& b)
            {  double dFP = (( a.nILTT == ILTT_SEGM || a.nILTT == ILTT_SEGM_ON_EDGE) ? 0.5 * ( a.dU + a.dU2) : a.dU) ;
               double dLP = (( b.nILTT == ILTT_SEGM || b.nILTT == ILTT_SEGM_ON_EDGE) ? 0.5 * ( b.dU + b.dU2) : b.dU) ;
               return ( dLP > dFP) ; }) ;

   return true ;
}

//----------------------------------------------------------------------------
// Il piano è espresso nel sistema locale, la funzione lo esprime nel sistema Zmap.
// I loop della flat region ottenuta dall'intersezione vengono salvati nello 
// standard vector vLoop passato per riferimento.
// Se il processo va a buon fine viene restituito true, false altrimenti.
bool
VolZmap::GetPlaneIntersection( const Plane3d& plPlane, ICURVEPOVECTOR& vpLoop) const 
{
  // Verifico validità parametri
   if ( ! plPlane.IsValid())
     return false ;

  // Porto il piano nel sistema Zmap
   Plane3d plPlaneLoc = plPlane ;
   plPlaneLoc.ToLoc( m_MapFrame) ;

  // Se non c'è intersezione fra piano e bounding-box del solido, ho finito.
   if ( ! TestIntersPlaneZmapBBox( plPlaneLoc))
      return true ;
  // Lancio eventuale aggiornamento pendente della grafica
   if ( m_nMapNum == 1)
      UpdateSingleMapGraphics() ;
   else
      UpdateTripleMapGraphics() ;
  // Vettore di segmenti 
   vector<CurveLine> vLine ;
  // Ciclo sui blocchi
   for ( int nB = 0 ; nB < m_nNumBlock ; ++ nB) {
     // Determino indici IJK del blocco; Se non trovo il blocco, errore 
      int nBlockIJK[3] ;
      if ( ! GetBlockIJKFromN( nB, nBlockIJK))
         return false ;
     // Costruisco il bounding-bx del blocco; se non è possiblie, errore
      BBox3d b3BlockBox ;
      if ( ! GetBlockBox( nBlockIJK, b3BlockBox))
         return false ;
     // Se c'è intersezione valuto tutti i voxel interni
      if ( TestIntersPlaneBox( plPlaneLoc, b3BlockBox)) {    
        // Ciclo sui voxel del blocco.    
        // Triangoli smooth
         for ( int nV = 0 ; nV < int( m_BlockSmoothTria[nB].size()) ; ++ nV) {
           // Box del voxel
            BBox3d b3Vox ;
            GetVoxelBox( m_BlockSmoothTria[nB][nV].i, m_BlockSmoothTria[nB][nV].j, m_BlockSmoothTria[nB][nV].k, b3Vox) ;
            // Se non c'è intersezione col voxel, passo al successivo.
            if ( ! TestIntersPlaneBox( plPlaneLoc, b3Vox))
               continue ; 
            for ( int nT = 0 ; nT < int( m_BlockSmoothTria[nB][nV].vTria.size()) ; ++ nT) {
               Triangle3d trTria = m_BlockSmoothTria[nB][nV].vTria[nT] ;
               Point3d ptSt, ptEn ;
               int nIntersType = IntersPlaneTria( plPlane, trTria, ptSt, ptEn) ;
               if ( nIntersType == IPTT_EDGE  ||  nIntersType == IPTT_YES) {
                 // Costruisco il tratto di curva
                  CurveLine cvLine ;
                  if ( cvLine.Set( ptSt, ptEn))
                     vLine.emplace_back( cvLine) ;
               }
            }
         }
        // Triangoli sharp interni al blocco    
         for ( int nV = 0 ; nV < int( m_BlockSharpTria[nB].size()) ; ++ nV) {
           // Box del voxel     
            BBox3d b3Vox ;
            GetVoxelBox( m_BlockSharpTria[nB][nV].i, m_BlockSharpTria[nB][nV].j, m_BlockSharpTria[nB][nV].k, b3Vox) ;
           // Se non c'è intersezione col voxel, passo al successivo.
           // Ciclo sulle componenti connesse
            for ( int nC = 0 ; nC < int( m_BlockSharpTria[nB][nV].vCompoTria.size()) ; ++ nC) {
               for ( int nT = 0 ; nT < int( m_BlockSharpTria[nB][nV].vCompoTria[nC].size()) ; ++ nT) {
                  Triangle3d trTria = m_BlockSharpTria[nB][nV].vCompoTria[nC][nT] ;
                  Point3d ptSt, ptEn ;
                  int nIntersType = IntersPlaneTria(plPlane, trTria, ptSt, ptEn) ;
                  if (nIntersType == IPTT_EDGE || nIntersType == IPTT_YES) {
                    // Costruisco il tratto di curva
                     CurveLine cvLine ;
                     if ( cvLine.Set(ptSt, ptEn))
                        vLine.emplace_back( cvLine) ;
                  }
               }
            }
         }
        // Triangoli grandi del blocco
         for ( int nT = 0 ; nT < int( m_BlockBigTria[nB].size()) ; ++ nT) {
            Triangle3d trTria = m_BlockBigTria[nB][nT] ;
            Point3d ptSt, ptEn ;
            int nIntersType = IntersPlaneTria(plPlane, trTria, ptSt, ptEn) ;
            if ( nIntersType == IPTT_EDGE  ||  nIntersType == IPTT_YES) {
               // Costruisco il tratto di curva
               CurveLine cvLine ;
               if ( cvLine.Set(ptSt, ptEn))
                  vLine.emplace_back( cvLine) ;
            }
         }
      }
     // In ogni caso valuto i triangoli sharp fra blocchi
      for ( int nV = 0 ; nV < int( m_InterBlockSharpTria[nB].size()) ; ++ nV) {
        // Ciclo sulle componenti connesse
         for ( int nC = 0 ; nC < int( m_InterBlockSharpTria[nB][nV].vCompoTria.size()) ; ++ nC) {
            for ( int nT = 0 ; nT < int( m_InterBlockSharpTria[nB][nV].vCompoTria[nC].size()) ; ++ nT) {
               Triangle3d trTria = m_InterBlockSharpTria[nB][nV].vCompoTria[nC][nT] ;
               Point3d ptSt, ptEn ;
               int nIntersType = IntersPlaneTria( plPlane, trTria, ptSt, ptEn) ;
               if ( nIntersType == IPTT_EDGE  ||  nIntersType == IPTT_YES) {
                  // Costruisco il tratto di curva
                  CurveLine cvLine ;
                  if ( cvLine.Set(ptSt, ptEn))
                     vLine.emplace_back( cvLine) ;
               }
            }
         }
      }   
   }
   
  // Creo i loop
   ChainCurves LoopCreator ;
   LoopCreator.Init( false, EPS_SMALL, int( vLine.size())) ;
  // Carico le curve per concatenarle
   for ( int nCv = 0 ; nCv < int( vLine.size()) ; ++ nCv) {
      Point3d ptSt = vLine[nCv].GetStart() ;
      Point3d ptEn = vLine[nCv].GetEnd() ;
      Vector3d vtDir; vLine[nCv].GetStartDir(vtDir) ;
      LoopCreator.AddCurve( nCv + 1, ptSt, vtDir, ptEn, vtDir) ;
   }
   // Recupero i concatenamenti
   INTVECTOR vIds ;
   Point3d ptNearStart ;
   while ( LoopCreator.GetChainFromNear( ptNearStart, false, vIds)) {
      PtrOwner<CurveComposite> pLoop( CreateBasicCurveComposite()) ;
      if ( IsNull( pLoop))
         return false ;
      for ( auto i : vIds) {
        // Aggiungo la linea alla curva composta.
         if ( ! pLoop->AddCurve( vLine[i - 1], true, 10 * EPS_SMALL))
            return false ;
      }
      pLoop->SetExtrusion( plPlane.GetVersN()) ;
      pLoop->MergeCurves( 10 * EPS_SMALL, ANG_TOL_STD_DEG) ;
      // Inserisco la curva composita nella raccolta da ritornare
      vpLoop.emplace_back( Release( pLoop)) ;
   }
   
   return true ;
}

//---------------------------------------------------------------------------- 
// Il piano deve essere espresso nel sistema di riferimento intrinseco dello Zmap.
bool
VolZmap::TestIntersPlaneZmapBBox( const Plane3d& plPlane) const 
{
   BBox3d b3Zmap( ORIG, Point3d( m_nNx[0] * m_dStep, m_nNy[0] * m_dStep, m_dMaxZ[0])) ;
   return TestIntersPlaneBox( plPlane, b3Zmap) ;
}