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

//--------------------------- Include ----------------------------------------
#include "stdafx.h"
#include "VolZmap.h"
#include "GeoObjFactory.h"
#include "NgeWriter.h"
#include "NgeReader.h"
#include "\EgtDev\Include\EGkIntervals.h"

using namespace std ;

//----------------------------------------------------------------------------
GEOOBJ_REGISTER( VOL_ZMAP, NGE_V_ZMP, VolZmap) ;

//----------------------------------------------------------------------------
VolZmap::VolZmap(void)
   : m_nStatus( TO_VERIFY), m_nTempProp()
{
   m_dStep = 0 ;
   m_nNx = 0 ;
   m_nNy = 0 ;
   m_nDim = 0 ;
}

//----------------------------------------------------------------------------
VolZmap::~VolZmap( void)
{
}

//----------------------------------------------------------------------------
VolZmap*
VolZmap::Clone( void) const
{
  // alloco oggetto
   VolZmap* pVzm = new(nothrow) VolZmap ;
   if ( pVzm != nullptr) {
      if ( ! pVzm->CopyFrom( *this)) {
         delete pVzm ;
         return nullptr ;
      }
   }

   return pVzm ;
}

//----------------------------------------------------------------------------
bool
VolZmap::CopyFrom( const IGeoObj* pGObjSrc)
{
   const VolZmap* pVzm = dynamic_cast<const VolZmap*>( pGObjSrc) ;
   if ( pVzm == nullptr)
      return false ;
   return CopyFrom( *pVzm) ;
}

//----------------------------------------------------------------------------
bool
VolZmap::CopyFrom( const VolZmap& vzmSrc)
{
   if ( &vzmSrc == this)
      return true ;
   m_OGrMgr.Reset() ;
   m_LocalFrame = vzmSrc.m_LocalFrame ;
   m_dStep = vzmSrc.m_dStep ;
   m_nDim = vzmSrc.m_nDim ;
   m_nNx = vzmSrc.m_nNx ;
   m_nNy = vzmSrc.m_nNy ;
   m_ZValues = vzmSrc.m_ZValues ;
   m_nStatus = vzmSrc.m_nStatus ;
   m_nTempProp = vzmSrc.m_nTempProp ;
   return true ;
}

//----------------------------------------------------------------------------
GeoObjType
VolZmap::GetType( void) const
{
   return static_cast<GeoObjType>( GEOOBJ_GETTYPE( VolZmap)) ;
}

//----------------------------------------------------------------------------
const string&
VolZmap::GetTitle( void) const
{
   static const string sTitle = "Zmap" ;
   return sTitle ;
}

//----------------------------------------------------------------------------
bool
VolZmap::Dump( string& sOut, bool bMM, const char* szNewLine) const
{
   return false ;
}

//----------------------------------------------------------------------------
int
VolZmap::GetNgeId( void) const
{
   return GEOOBJ_GETNGEID( VolZmap) ;
}

//----------------------------------------------------------------------------
bool
VolZmap::Save( NgeWriter& ngeOut) const
{
  // parametri di scrittura: sistema di riferimento, minimo incremento, numero di passi
  //             in direzione x e y, e per ogni casella, numero di valori e valori
   if ( ! ngeOut.WriteFrame( m_LocalFrame, ";", true))
      return false ;
   if ( ! ngeOut.WriteDouble( m_dStep, ",", false))
      return false ;
   if ( ! ngeOut.WriteInt( m_nNx, ",", false))
      return false ;
   if ( ! ngeOut.WriteInt( m_nNy, ";", true))
      return false ;
  // ciclo sui dexel
   for ( unsigned int i = 0 ; i < m_nDim ; ++ i) {
     // numero di estremi
      int nDim = int( m_ZValues[i].size()) ;
      if ( ! ngeOut.WriteInt( nDim, ",", false))
         return false ;
     // se dexel nullo
      if ( nDim == 0) {
        // scrivo un valore dummy
         if ( ! ngeOut.WriteDouble( 0, ";", true))
            return false ;
      }
     // altrimenti 
      else {
         for ( unsigned int k = 0 ; k < m_ZValues[i].size() ; ++ k) {
            bool bEndL = ( k == m_ZValues[i].size() - 1) ;
            if ( ! ngeOut.WriteDouble( m_ZValues[i][k], ( bEndL ? ";" : ","), bEndL))
               return false ;
         }
      }
   }
         
   return true ;
}
    
//----------------------------------------------------------------------------       
bool
VolZmap::Load( NgeReader& ngeIn)
{
   m_nStatus = TO_VERIFY ;
  // parametri di lettura: sistema di riferimento, minimo incremento, numero di passi
  //             in direzione x e y, e per ogni casella, numero di valori e valori
   if ( ! ngeIn.ReadFrame( m_LocalFrame, ";", true))
      return false ;
   if ( ! ngeIn.ReadDouble( m_dStep, ",", false))
      return false ;
   int nTemp ;
   if ( ! ngeIn.ReadInt( nTemp, ",", false))
      return false ;
   m_nNx = nTemp ;
   if ( ! ngeIn.ReadInt( nTemp, ";", true))
      return false ;
   m_nNy = nTemp ;

  // dimensione del vettore di dexel 
   m_nDim = m_nNx * m_nNy ;
   m_ZValues.resize(m_nDim) ;

  // ciclo sui dexel
   for ( unsigned int i = 0 ; i < m_nDim ; ++ i)  {
     // leggo il numero di estremi nel dexel
      if ( ! ngeIn.ReadInt( nTemp, ",", false))
         return false ;
     // devono essere pari
      if ( ( nTemp % 2) != 0)
         return false ;
     // se dexel nullo
      if ( nTemp == 0) {
        // leggo un valore dummy
         double dDummy ;
         if ( ! ngeIn.ReadDouble( dDummy, ",", true))
            return false ;
      }
     // altrimenti 
      else {
        // dimensiono l'array
         m_ZValues[i].resize(nTemp) ;
        // leggo i valori
         for ( unsigned int k = 0 ; k < m_ZValues[i].size() ; ++ k) {
            bool bEndL = ( k == m_ZValues[i].size() - 1) ;
            if ( ! ngeIn.ReadDouble( m_ZValues[i][k], ( bEndL ? ";" : ","), bEndL))
               return false ;
         }
      }
   }
   
   m_nStatus = OK ;
   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::GetLocalBBox( BBox3d& b3Loc, int nFlag) const
{
  // verifico lo stato
   if ( m_nStatus != OK)
      return false ;
  // reset box
   b3Loc.Reset() ;
  // ciclo sui dexel
   double dY = 0.5 * m_dStep ;
   for ( size_t i = 0 ; i < m_nNy ; ++ i) {
      double dX = 0.5 * m_dStep ;
      for ( size_t j = 0 ; j < m_nNx ; ++ j) {
         size_t nPos = j + i * m_nNx ;
         if ( m_ZValues[nPos].size() > 0) {
            Point3d ptP = m_LocalFrame.Orig() + dX * m_LocalFrame.VersX() + dY * m_LocalFrame.VersY() ;
            b3Loc.Add( ptP + m_ZValues[nPos][0] * m_LocalFrame.VersZ()) ;
            b3Loc.Add( ptP + m_ZValues[nPos][m_ZValues[nPos].size()-1] * m_LocalFrame.VersZ()) ;
         }
         dX += m_dStep ;
      }
      dY += m_dStep ;
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::GetBBox( const Frame3d& frRef, BBox3d& b3Ref, int nFlag) const
{
  // verifico lo stato
   if ( m_nStatus != OK)
      return false ;
  // reset box
   b3Ref.Reset() ;
  // trasformo il riferimento locale tramite quello passato
   Frame3d frUse = m_LocalFrame ;
   frUse.ToGlob( frRef) ;
  // ciclo sui dexel
   double dY = 0.5 * m_dStep ;
   for ( size_t i = 0 ; i < m_nNy ; ++ i) {
      double dX = 0.5 * m_dStep ;
      for ( size_t j = 0 ; j < m_nNx ; ++ j) {
         size_t nPos = j + i * m_nNx ;
         if ( m_ZValues[nPos].size() > 0) {
            Point3d ptP = frUse.Orig() + dX * frUse.VersX() + dY * frUse.VersY() ;
            b3Ref.Add( ptP + m_ZValues[nPos][0] * frUse.VersZ()) ;
            b3Ref.Add( ptP + m_ZValues[nPos][m_ZValues[nPos].size()-1] * frUse.VersZ()) ;
         }
         dX += m_dStep ;
      }
      dY += m_dStep ;
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::Translate( const Vector3d& vtMove)
{
  // verifico lo stato
   if ( m_nStatus != OK)
      return false ;
  // imposto ricalcolo della grafica
   m_OGrMgr.Reset() ;
  // traslo il riferimento
   m_LocalFrame.Translate( vtMove) ;
   return true ;
}

//----------------------------------------------------------------------------
bool 
VolZmap::Rotate( const Point3d& ptAx, const Vector3d& vtAx, double dCosAng, double dSinAng)
{
  // verifico lo stato
   if ( m_nStatus != OK)
      return false ;
  // imposto ricalcolo della grafica
   m_OGrMgr.Reset() ;
  // ruoto il riferimento
   return m_LocalFrame.Rotate( ptAx, vtAx, dCosAng, dSinAng) ;
}

//----------------------------------------------------------------------------
bool 
VolZmap::Scale( const Frame3d& frRef, double dCoeffX, double dCoeffY, double dCoeffZ)
{
  // verifico lo stato
   if ( m_nStatus != OK)
      return false ;
   return false ;
}

//----------------------------------------------------------------------------
bool 
VolZmap::Mirror( const Point3d& ptOn, const Vector3d& vtNorm)
{
  // verifico lo stato
   if ( m_nStatus != OK)
      return false ;
   return false ;
}

//----------------------------------------------------------------------------
bool 
VolZmap::Shear( const Point3d& ptOn, const Vector3d& vtNorm, const Vector3d& vtDir, double dCoeff)
{
  // verifico lo stato
   if ( m_nStatus != OK)
      return false ;
   return false ;
}

//----------------------------------------------------------------------------
bool 
VolZmap::ToGlob( const Frame3d& frRef)
{
  // verifico lo stato
   if ( m_nStatus != OK)
      return false ;
   return false ;
}

//----------------------------------------------------------------------------
bool 
VolZmap::ToLoc( const Frame3d& frRef)
{
  // verifico lo stato
   if ( m_nStatus != OK)
      return false ;
   return false ;
}

//----------------------------------------------------------------------------
bool 
VolZmap::LocToLoc( const Frame3d& frOri, const Frame3d& frDest)
{
  // verifico lo stato
   if ( m_nStatus != OK)
      return false ;
   return false ;
}

//----------------------------------------------------------------------------
bool
VolZmap::CreateMap( const Point3d& ptO, double dPrec, double dLengthX, double dLengthY, double dLengthZ)
{
  // Controlli sui parametri
   if ( dPrec < EPS_SMALL || dLengthX < EPS_SMALL || dLengthY < EPS_SMALL || dLengthZ < EPS_SMALL)
      return false ; 

   m_LocalFrame.Set( ptO, X_AX, Y_AX, Z_AX) ;

   m_dStep = dPrec ;

   m_nNx = static_cast <unsigned int> ( ceil( dLengthX / m_dStep)) ;
   m_nNy = static_cast <unsigned int> ( ceil( dLengthY / m_dStep)) ;

   m_nDim = m_nNx * m_nNy ;

   m_ZValues.resize( m_nDim) ;

   for ( int i = 0 ; i < int( m_nDim) ; i++) {
      m_ZValues[i].resize(2) ;

      m_ZValues[i][0] = 0 ;
      m_ZValues[i][1] = dLengthZ ;
   }

   m_dMinZ = 0 ;
   m_dMaxZ = dLengthZ ;

   m_nStatus = OK ;
   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::SubtractIntervals( unsigned int nI, unsigned int nJ, double dMin, double dMax)    
{
   unsigned int nPos ;
   //unsigned int nKD ;
   //unsigned int nKU ;
   //unsigned int nDelta ; // Variazione di dimensione del vettore 

   // Controllo che dMin e dMax non siano quasi coincidenti
   if ( abs( dMax - dMin) < EPS_SMALL)
   
      return true ;
   // Controllo che dMin < dMax

   if ( dMax < dMin ) {
   
         double dTemp = dMax ;
   
         dMax = dMin ;
         dMin = dTemp ; // cambiato il 3/2
   }
   
   // Calcolo nPos
   nPos = nJ*m_nNx + nI ;  // O nJ*m_Nx + nI + 1 ?

   unsigned int i = 0 ;

   while ( i < m_ZValues[nPos].size() - 1) {

      if ( m_ZValues[nPos].size() == 0)
         return true ;
      // Casi:
      // Intervallo da sottrarre è tutto a sinistra di quello corrente, non vi è intersezione 
      if ( m_ZValues[nPos][i] > dMax - EPS_SMALL) {
      
      }
      // Intersezione
      else if ( m_ZValues[nPos][i + 1] > dMax + EPS_SMALL) {
         // L'intervallo corrente corrente viene limitato a sinistra
         if ( m_ZValues[nPos][i] > dMin - EPS_SMALL) {
         
            m_ZValues[nPos][i] = dMax ;
         }
         // L'intervallo si divide in due intervalli
         else {

            m_ZValues[nPos].resize( m_ZValues[nPos].size() + 2) ;

            for ( size_t j = m_ZValues[nPos].size() - 1 ; j >= i + 3 ; -- j)

               m_ZValues[nPos][j] = m_ZValues[nPos][j - 2] ;

            m_ZValues[nPos][i + 1] = dMin ;
            m_ZValues[nPos][i + 2] = dMax ;

            i = i + 2 ;
         }
      }
      else {
         // L'intervallo corrente viene eliminato
         if ( m_ZValues[nPos][i] > dMin - EPS_SMALL) {
         
            for ( unsigned int j = i ; j < m_ZValues[nPos].size() - 2 ; ++ j)

               m_ZValues[nPos][j] = m_ZValues[nPos][j + 2] ;

            m_ZValues[nPos].resize( m_ZValues[nPos].size() - 2) ;

            i = i - 2 ;
         }
         // L'intervallo corrente viene limitato a destra
         else if ( m_ZValues[nPos][i + 1] > dMin + EPS_SMALL) {
         
            m_ZValues[nPos][i + 1] = dMin ;
         }
         // L'intervallo da sottrarre è tutto a destra di quello corrente, non vi è intersezione
         else {
         
         }
      }
      
      i = i + 2 ;
   }

    return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::SubtractIntervals( const Point3d& ptP, double dMin, double dMax) // ptP è espresso nel sistema locale e viene convertito in quello intrinseco (localFrame) 
{
   Point3d ptPL = ptP ;
   ptPL.ToLoc( m_LocalFrame) ;  // ptPL è ora espresso rispetto al sistema di riferimento intrinseco.

   double dX, dY, dZ ;    // Coordinate di ptPL nel sistema intrinseco
   double dhMin, dhMax ;  // Altezze dMin e dMax RIESPRESSE nel sistema intrinseco (dMin e dMax sono altezze rispetto a ptP)

   dX = ptPL.x ; dY = ptPL.y ; dZ = ptPL.z ;
   dhMin = dZ + dMin ; dhMax = dZ + dMax ;

   // Cerco il punto della griglia più vicino

   double integerPartX = floor( dX / m_dStep) ;
   double integerPartY = floor( dY / m_dStep) ;
   
   unsigned int i = static_cast <unsigned int> (integerPartX) ; // Indici del punto di griglia più vicino.
   unsigned int j = static_cast <unsigned int> (integerPartY) ; // i = 0, 1, ..., m_Nx - 1 ; j = 0, 1, ..., m_Ny - 1

   // Controllo se le coordinate x e y del punto dato siano all'interno della griglia:
   // se sono dentro la griglia chiamo l'altra subtract
   if ( dX < m_dStep*m_nNx && dY < m_dStep*m_nNy
        && dX >= 0 && dY >= 0) { // Mettendo > - qlc può sempre capitare un punto compreso fra - qlc e 0 e si esce dai limiti dell vector e vi è errore runtime
      
      return SubtractIntervals( i, j, dhMin, dhMax) ;
   }
   // altrimenti non succede niente
   else {
   
      return true ;
   }

}

//----------------------------------------------------------------------------
bool
VolZmap::AddIntervals( unsigned int nI, unsigned int nJ, double dMin, double dMax)     
{
   unsigned int nPos ;
   //unsigned int nKD ;
   //unsigned int nKU ;
   //unsigned int nDelta ; // Variazione di dimensione del vettore 

   // Controllo che dMin e dMax non siano quasi coincidenti
   if ( abs( dMax - dMin) < EPS_SMALL)
   
      return true ;
   
   // Controllo che dMin < dMax
   if ( dMax < dMin ) {
   
      double dTemp = dMax ;
   
      dMax = dMin ;
      dMin = dTemp ; // cambiato il 3/2
   }
   
   // Calcolo nPos
   nPos = nJ*m_nNx + nI ;  // O nJ*m_Nx + nI + 1 ?

   if ( m_ZValues[nPos].size() == 0) {
   
      m_ZValues[nPos].resize( 2) ;
   
      m_ZValues[nPos][0] = dMin ;
      m_ZValues[nPos][1] = dMax ;

      return true ;
   }

   unsigned int i = 0 ;

   while ( i < m_ZValues[nPos].size() - 1) {

      if ( i > 0) {
         // Caso in cui un intervallo precedente sconfini in quello corrente
         if ( m_ZValues[nPos][i] < m_ZValues[nPos][i - 1] + EPS_SMALL) {
            // Se l'intervallo corrente non è contenuto totalmente si esegue l'istruzione successiva
            if ( m_ZValues[nPos][i - 1] < m_ZValues[nPos][i + 1] + EPS_SMALL)
         
               m_ZValues[nPos][i - 1] = m_ZValues[nPos][i + 1] ;
         
            for ( unsigned int j = i ; j < m_ZValues[nPos].size() - 2 ; ++ j) 
            
               m_ZValues[nPos][j] = m_ZValues[nPos][j + 2] ;
         
            m_ZValues[nPos].resize( m_ZValues[nPos].size() - 2) ; 
            
            i = i - 2 ; 
         }
      }

      // Caso in cui devo aggiungere un intervallo a sinistra dell'intervallo corrente
      if ( m_ZValues[nPos][i] > dMax + EPS_SMALL) {
         
            m_ZValues[nPos].resize( m_ZValues[nPos].size() + 2) ;
         
            for ( size_t j = m_ZValues[nPos].size() - 1 ; j >= i + 2 ; -- j) 
               
               m_ZValues[nPos][j] = m_ZValues[nPos][j - 2] ;
         
            m_ZValues[nPos][i] = dMin ;
            m_ZValues[nPos][i + 1] = dMax ;

            i = i + 2 ;
      } 
      // Casi d'intersezione:
      else if ( m_ZValues[nPos][i + 1] > dMax - EPS_SMALL) {
         // Se l'intervallo da aggiungere sconfina a sinistra modifico il minimo dell'intervalo corrente
         if (m_ZValues[nPos][i] > dMin - EPS_SMALL) { 

            m_ZValues[nPos][i] = dMin ;
         }
      }
      else {
         // Se l'intervallo corrente è tutto contenuto nell'intervallo da aggungere modifico gli estremi    
         if ( m_ZValues[nPos][i] > dMin + EPS_SMALL) {
      
               m_ZValues[nPos][i] = dMin ;
               m_ZValues[nPos][i + 1] = dMax ;
         }
         // Se l'intervallo da aggiungere sconfina a destra modifico il massimo dell'intervallo corrente 
         else if ( m_ZValues[nPos][i + 1] > dMin - EPS_SMALL) {
      
            m_ZValues[nPos][i + 1] = dMax ;
         }
         else {
            // Aggiungo intervallo a destra dell'ultimo intervallo
            if ( i == m_ZValues[nPos].size() - 2) {
            
               m_ZValues[nPos].resize( m_ZValues[nPos].size() + 2) ;
               
               m_ZValues[nPos][i + 2] = dMin ;
               m_ZValues[nPos][i + 3] = dMax ;

               i = i + 2 ;
            }
         }
      }

      i = i + 2 ;
   }

   return true ;    
}

//----------------------------------------------------------------------------
bool
VolZmap::AddIntervals( const Point3d& ptP, double dMin, double dMax) // ptP è espresso nel sistema locale e viene convertito in quello intrinseco (localFrame) 
{ 
   Point3d ptPL = ptP ;
   ptPL.ToLoc( m_LocalFrame) ;  // ptPL è ora espresso rispetto al sistema di riferimento intrinseco.

   double dX, dY, dZ ;    // Coordinate di ptPL nel sistema intrinseco
   double dhMin, dhMax ;  // Altezze dMin e dMax RIESPRESSE nel sistema intrinseco (dMin e dMax sono altezze rispetto a ptP)

   dX = ptPL.x ; dY = ptPL.y ; dZ = ptPL.z ;
   dhMin = dZ + dMin ; dhMax = dZ + dMax ;

   // Cerco il punto della griglia più vicino

   double integerPartX = floor( dX / m_dStep) ;
   double integerPartY = floor( dY / m_dStep) ;
   
   unsigned int i = static_cast <unsigned int> (integerPartX) ; // Indici del punto di griglia più vicino.
   unsigned int j = static_cast <unsigned int> (integerPartY) ; // i = 0, 1, ..., m_Nx - 1 ; j = 0, 1, ..., m_Ny - 1

   // Controllo se le coordinate x e y del punto dato siano all'interno della griglia:
   // se sono dentro la griglia chiamo l'altra subtract
   if ( dX < m_dStep*m_nNx && dY < m_dStep*m_nNy
        && dX >= 0 && dY >= 0) {  // Mettendo > - qlc può sempre capitare un punto compreso fra - qlc e 0 e si esce dai limiti dell vector e vi è errore runtime
      
      return AddIntervals( i, j, dhMin, dhMax) ;
   }
   // altrimenti non succede niente
   else {
   
      return false ;
   }     
}

//----------------------------------------------------------------------------
bool
VolZmap::GetDexelLines( int nDir, int nPos1, int nPos2, POLYLINELIST& lstPL) const
{
  // per ora solo perpendicolari a XY (1)
   if ( nDir != 1)
      return false ;
  // verifiche sugli indici
   if ( nPos1 < 0  ||  nPos1 >= int( m_nNx)  ||  nPos2 < 0  ||  nPos2 >= int( m_nNy))
      return false ;
   int nPos = nPos1 + nPos2 * m_nNx ;
   if ( nPos < 0  ||  nPos >= int( m_ZValues.size()))
      return false ;
  // calcolo coordinate punto
   double dX = m_dStep * ( 0.5 + nPos1) ;
   double dY = m_dStep * ( 0.5 + nPos2) ;
   Point3d ptP = m_LocalFrame.Orig() + dX * m_LocalFrame.VersX() + dY * m_LocalFrame.VersY() ;
  // creo le polilinee
   for ( int i = 1 ; i < int( m_ZValues[nPos].size()) ; i += 2) {
     // aggiungo polilinea a lista
      lstPL.emplace_back() ;
     // inserisco punti estremi
      lstPL.back().AddUPoint( 0, ptP + m_ZValues[nPos][i-1] * m_LocalFrame.VersZ()) ;
      lstPL.back().AddUPoint( 1, ptP + m_ZValues[nPos][i] * m_LocalFrame.VersZ()) ;
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::GetAllTriangles( TRIA3DLIST& lstTria) const
{
   const int DIM_CHUNK = 4 ;
   for ( int i = 0 ; i < int( m_nNx) ; i += DIM_CHUNK) {
      int nDimChunkX = min( DIM_CHUNK, int( m_nNx) - i) ;
      for ( int j = 0 ; j < int( m_nNy) ; j += DIM_CHUNK) {
         int nDimChunkY = min( DIM_CHUNK, int( m_nNy) - j) ;
         GetChunkPrisms( i, j, nDimChunkX, nDimChunkY, lstTria) ;
      } 
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::GetChunkPrisms( int nPos1, int nPos2, int nDim1, int nDim2, TRIA3DLIST& lstTria) const
{
  // determino se è un semplice parallelepipedo
   bool bIsSimple = true ;
   double dBotZ ;
   double dTopZ ;
   for ( int i = 0 ; i < nDim1 && bIsSimple ; ++ i) {
      for ( int j = 0 ; j < nDim2 && bIsSimple ; ++ j) {
         int nPos = ( nPos1 + i) + ( nPos2 + j) * m_nNx ;
         if ( nPos > int( m_nDim)  ||
              int( m_ZValues[nPos].size()) != 2)
            bIsSimple = false ;  
         else if ( i == 0  &&  j == 0) {
            dBotZ = m_ZValues[nPos][0] ;
            dTopZ = m_ZValues[nPos][1] ;
         }
         else if ( abs( m_ZValues[nPos][0] - dBotZ) > EPS_SMALL  ||
                   abs( m_ZValues[nPos][1] - dTopZ) > EPS_SMALL)
            bIsSimple = false ;  
      }
   }

  // se semplice parallelepipedo
   if ( bIsSimple) {
      CalcChunkPrisms( nPos1, nPos2, nDim1, nDim2, lstTria) ;
   }
  // altrimenti
   else {
     // elaboro ogni singolo dexel
      for ( int i = 0 ; i < nDim1 ; ++ i) {
         for ( int j = 0 ; j < nDim2 ; ++ j) {
            CalcDexelPrisms( nPos1 + i, nPos2 + j, lstTria) ;
         } 
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::CalcChunkPrisms( int nPos1, int nPos2, int nDim1, int nDim2, TRIA3DLIST& lstTria) const
{
  // verifiche sugli indici
   if ( nPos1 < 0  ||  nPos1 + nDim1 > int( m_nNx)  ||  nPos2 < 0  ||  nPos2 + nDim2 > int( m_nNy))
      return false ;
   int nPos = nPos1 + nPos2 * m_nNx ;
   if ( nPos < 0  ||  nPos >= int( m_nDim))
      return false ;

  // calcolo coordinate punti
   double dX = m_dStep * nPos1 ;
   double dY = m_dStep * nPos2 ;
   Point3d ptP1 = m_LocalFrame.Orig() + dX * m_LocalFrame.VersX() + dY * m_LocalFrame.VersY() ;
   Point3d ptP2 = ptP1 + nDim1 * m_dStep * m_LocalFrame.VersX() ;
   Point3d ptP3 = ptP2 + nDim2 * m_dStep * m_LocalFrame.VersY() ;
   Point3d ptP4 = ptP1 + nDim2 * m_dStep * m_LocalFrame.VersY() ;

  // creo le facce sopra e sotto
   Vector3d vtDZt = m_ZValues[nPos][1] * m_LocalFrame.VersZ() ;
   Vector3d vtDZb = m_ZValues[nPos][0] * m_LocalFrame.VersZ() ;
   // faccia superiore P1t->P2t->P3t->P4t : sempre visibile
   lstTria.emplace_back() ;
   lstTria.back().Set( ptP1 + vtDZt, ptP2 + vtDZt, ptP3 + vtDZt, m_LocalFrame.VersZ()) ;
   lstTria.emplace_back() ;
   lstTria.back().Set( ptP3 + vtDZt, ptP4 + vtDZt, ptP1 + vtDZt, m_LocalFrame.VersZ()) ;
   // faccia inferiore P1b->P4b->P3b->P2b : sempre visibile
   lstTria.emplace_back() ;
   lstTria.back().Set( ptP1 + vtDZb, ptP4 + vtDZb, ptP3 + vtDZb, - m_LocalFrame.VersZ()) ;
   lstTria.emplace_back() ;
   lstTria.back().Set( ptP3 + vtDZb, ptP2 + vtDZb, ptP1 + vtDZb, - m_LocalFrame.VersZ()) ;

  // creo le facce laterali
   for ( int j = 0 ; j < nDim2 ; ++ j) {
      int nPosD = nPos + nDim1 - 1 + j * m_nNx ;
      int nPosEst = ( nPos1 + nDim1 - 1 < int( m_nNx - 1) ? nPosD + 1 : - 1) ;
      Point3d ptP2D = ptP2 + j * m_dStep * m_LocalFrame.VersY() ;
      Point3d ptP3D = ptP2D + m_dStep * m_LocalFrame.VersY() ;
      AddDexelSideFace( nPosD, nPosEst, ptP2D, ptP3D, m_LocalFrame.VersZ(), m_LocalFrame.VersX(), lstTria) ;
   }
   for ( int i = 0 ; i < nDim1 ; ++ i) {
      int nPosD = nPos + ( nDim2 - 1) * m_nNx + i ;
      int nPosNord = ( nPos2 + nDim2 - 1 < int( m_nNy - 1) ? nPosD + m_nNx : - 1) ;
      Point3d ptP4D = ptP4 + i * m_dStep * m_LocalFrame.VersX() ;
      Point3d ptP3D = ptP4D + m_dStep * m_LocalFrame.VersX() ;
      AddDexelSideFace( nPosD, nPosNord, ptP3D, ptP4D, m_LocalFrame.VersZ(), m_LocalFrame.VersY(), lstTria) ;
   }
   for ( int j = 0 ; j < nDim2 ; ++ j) {
      int nPosD = nPos + j * m_nNx ;
      int nPosWest = ( nPos1 > 0 ? nPosD - 1 : - 1) ;
      Point3d ptP1D = ptP1 + j * m_dStep * m_LocalFrame.VersY() ;
      Point3d ptP4D = ptP1D + m_dStep * m_LocalFrame.VersY() ;
      AddDexelSideFace( nPosD, nPosWest, ptP4D, ptP1D, m_LocalFrame.VersZ(), - m_LocalFrame.VersX(), lstTria) ;
   }
   for ( int i = 0 ; i < nDim1 ; ++ i) {
      int nPosD = nPos + i ;
      int nPosSud = ( nPos2 > 0 ? nPosD - m_nNx : - 1) ;
      Point3d ptP1D = ptP1 + i * m_dStep * m_LocalFrame.VersX() ;
      Point3d ptP2D = ptP1D + m_dStep * m_LocalFrame.VersX() ;
      AddDexelSideFace( nPosD, nPosSud, ptP1D, ptP2D, m_LocalFrame.VersZ(), - m_LocalFrame.VersY(), lstTria) ;
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::CalcDexelPrisms( int nPos1, int nPos2, TRIA3DLIST& lstTria) const
{
  // verifiche sugli indici
   if ( nPos1 < 0  ||  nPos1 >= int( m_nNx)  ||  nPos2 < 0  ||  nPos2 >= int( m_nNy))
      return false ;
   int nPos = nPos1 + nPos2 * m_nNx ;
   if ( nPos < 0  ||  nPos >= int( m_nDim))
      return false ;

  // calcolo coordinate punto
   double dX = m_dStep * nPos1 ;
   double dY = m_dStep * nPos2 ;
   Point3d ptP1 = m_LocalFrame.Orig() + dX * m_LocalFrame.VersX() + dY * m_LocalFrame.VersY() ;
   Point3d ptP2 = ptP1 + m_dStep * m_LocalFrame.VersX() ;
   Point3d ptP3 = ptP2 + m_dStep * m_LocalFrame.VersY() ;
   Point3d ptP4 = ptP1 + m_dStep * m_LocalFrame.VersY() ;

  // creo le facce sopra e sotto di ogni intervallo (sempre visibili)
   for ( int i = 1 ; i < int( m_ZValues[nPos].size()) ; i += 2) {
      Vector3d vtDZt = m_ZValues[nPos][i] * m_LocalFrame.VersZ() ;
      Vector3d vtDZb = m_ZValues[nPos][i-1] * m_LocalFrame.VersZ() ;
     // faccia superiore P1t->P2t->P3t->P4t : sempre visibile
      lstTria.emplace_back() ;
      lstTria.back().Set( ptP1 + vtDZt, ptP2 + vtDZt, ptP3 + vtDZt, m_LocalFrame.VersZ()) ;
      lstTria.emplace_back() ;
      lstTria.back().Set( ptP3 + vtDZt, ptP4 + vtDZt, ptP1 + vtDZt, m_LocalFrame.VersZ()) ;
     // faccia inferiore P1b->P4b->P3b->P2b : sempre visibile
      lstTria.emplace_back() ;
      lstTria.back().Set( ptP1 + vtDZb, ptP4 + vtDZb, ptP3 + vtDZb, - m_LocalFrame.VersZ()) ;
      lstTria.emplace_back() ;
      lstTria.back().Set( ptP3 + vtDZb, ptP2 + vtDZb, ptP1 + vtDZb, - m_LocalFrame.VersZ()) ;
   }

  // creo le facce laterali
   int nPosEst = ( nPos1 < int( m_nNx - 1) ? nPos + 1 : - 1) ;
   AddDexelSideFace( nPos, nPosEst, ptP2, ptP3, m_LocalFrame.VersZ(), m_LocalFrame.VersX(), lstTria) ;
   int nPosNord = ( nPos2 < int( m_nNy - 1) ? nPos + m_nNx : - 1) ;
   AddDexelSideFace( nPos, nPosNord, ptP3, ptP4, m_LocalFrame.VersZ(), m_LocalFrame.VersY(), lstTria) ;
   int nPosWest = ( nPos1 > 0 ? nPos - 1 : - 1) ;
   AddDexelSideFace( nPos, nPosWest, ptP4, ptP1, m_LocalFrame.VersZ(), - m_LocalFrame.VersX(), lstTria) ;
   int nPosSud = ( nPos2 > 0 ? nPos - m_nNx : - 1) ;
   AddDexelSideFace( nPos, nPosSud, ptP1, ptP2, m_LocalFrame.VersZ(), - m_LocalFrame.VersY(), lstTria) ;

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::AddDexelSideFace( int nPos, int nPosAdj, const Point3d& ptP, const Point3d& ptQ,
                           const Vector3d& vtZ, const Vector3d& vtNorm, TRIA3DLIST& lstTria) const
{
   Intervals intFace ;
   for ( int i = 1 ; i < int( m_ZValues[nPos].size()) ; i += 2)
       intFace.Add( m_ZValues[nPos][i-1], m_ZValues[nPos][i]) ;
   if ( nPosAdj > 0) {
      for ( int i = 1 ; i < int( m_ZValues[nPosAdj].size()) ; i += 2)
          intFace.Subtract( m_ZValues[nPosAdj][i-1], m_ZValues[nPosAdj][i]) ;
   }
   double dMin, dMax ;
   bool bFound = intFace.GetFirst( dMin, dMax) ;
   while ( bFound) {
      Vector3d vtDZt = dMax * vtZ ;
      Vector3d vtDZb = dMin * vtZ ;
      lstTria.emplace_back() ;
      lstTria.back().Set( ptP + vtDZb, ptQ + vtDZb, ptQ + vtDZt, vtNorm) ;
      lstTria.emplace_back() ;
      lstTria.back().Set( ptQ + vtDZt, ptP + vtDZt, ptP + vtDZb, vtNorm) ;
      bFound = intFace.GetNext( dMin, dMax) ;
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::SetTool( const string& pToolName, const CurveComposite* pToolOutline)
{
   m_sToolName = pToolName ;

   m_nToolType = GenericTool ;

   m_ToolOutline.CopyFrom( pToolOutline) ;

   // Dimensioni dell'utensile
   BBox3d Bounding ;   
   m_ToolOutline.GetLocalBBox( Bounding) ;     
   double m_dHeight = Bounding.GetMax().y - Bounding.GetMin().y ;   
   double m_dRadius = Bounding.GetMax().x - Bounding.GetMin().x ; 

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::SetStdTool( const string& pToolName, unsigned int nToolType, double dH, double dTH, double dR, double dTR, double dRc)
{
   // Controllo sull'ammissibilità del tipo di utensile
   if ( nToolType == GenericTool || nToolType > 4)

      return false ;
   else {
      // Controllo sull'ammissibilità dei parametri geometrici dell'utensile
      if ( dH < EPS_SMALL || dR < EPS_SMALL)
         
         return false ;

      m_sToolName = pToolName ;
      m_nToolType = nToolType ;

      m_dRadius = dR ;
      m_dHeight = dH ;

      if ( nToolType == CylindricalMill) {
      
         m_dTipHeight = 0 ;
         m_dTipRadius = m_dRadius ;
      }
      else if ( nToolType == BallEndMill) {
      
         m_dTipHeight = m_dRadius ;
         m_dTipRadius = 0 ;
      }
      else if ( nToolType == BullNoseMill) { 

         if ( dR < dRc)

            return false ; 
      
         m_dRCorner = dRc ;
         m_dTipHeight = dRc ;
         m_dTipRadius = dR - dRc ;

         if ( dRc < EPS_SMALL) {
            
            m_dRCorner = 0 ;
            m_nToolType = CylindricalMill ; 
         }

         if ( dRc > dR - EPS_SMALL) {
         
            m_dRCorner = m_dRadius ;
            m_nToolType = BallEndMill ;
         }
      }
      else if ( nToolType == ConusMill) {

         if ( dTH < 0 || dTR < 0)
            
            return false ;

         if ( abs( dR - dTR) < EPS_SMALL) {
         
            m_nToolType = CylindricalMill ;

            m_dTipHeight = 0 ;
            m_dTipRadius = m_dRadius ; 
         }
         else {
         
            m_dTipHeight = dTH ;
            m_dTipRadius = dTR ;
         }
      }
      return true ;
   }
}

//----------------------------------------------------------------------------
bool
VolZmap::MillingStep( const Point3d& ptPs, const Point3d& ptPe, Vector3d& vtDs, Vector3d& vtDe) 
{                                                                                               
                                                                                                
  // Porto i dati del movimento nel riferimento intrinseco                                      
   Point3d ptLs = ptPs ;            
   ptLs.ToLoc( m_LocalFrame) ;      
   Point3d ptLe = ptPe ; 
   ptLe.ToLoc( m_LocalFrame) ;
   Vector3d vtLs = vtDs ;
   vtLs.ToLoc( m_LocalFrame) ;
   Vector3d vtLe = vtDe ;
   vtLe.ToLoc( m_LocalFrame) ;

   if ( AreSamePointApprox( ptLs, ptLe) && AreSameVectorApprox( vtLs, vtLe)) 
      return true ;

  // Normalizzo i vettori 
   vtLs.Normalize() ;
   vtLe.Normalize() ;

  // Direzione utensile costante 
   if ( AreSameVectorApprox( vtLs, vtLe)) {  
     // Versori della direzione utensile paralleli/ antiparalleli a Z
      if ( vtLs.IsZplus() || vtLs.IsZminus()) {
        // Movimento diretto come direzione utensile
         if ( AreSamePointXYApprox( ptLs, ptLe)) 
            return MillingDrillZ( ptLs, ptLe, vtLs) ;
        // Movimento perpendicolare a direzione utensile
         else if ( abs( ptLe.z - ptLs.z) < EPS_SMALL)
            return MillingPerpZ( ptLs, ptLe, vtLs) ; 
        // Movimento generico
         else
            return MillingZ( ptLs, ptLe, vtLs) ;
      }
     // Versori della direzione utensile perpendicolari all'asse Z
      else if ( abs(vtLs * Z_AX) < EPS_SMALL) {
        // Movimento con Z costante (con vettore movimento parallelo od ortogonale al versore dell'utensile)
         if ( abs( ptLe.z - ptLs.z) < EPS_SMALL) {
           // Movimento parallelo alla direzione dell'utensile (foratura) 
            Vector3d vtTest = ptLe - ptLs ;  vtTest.Normalize() ;
            if ( abs( ( vtTest * vtLs) - 1) < EPS_SMALL)   
                return MillingDrillXY( ptLs, ptLe, vtLs) ;
           // Movimento perpendicolare alla direzione dell'utensile 
            else if ( abs( (ptLe - ptLs) * vtLs) < EPS_SMALL) 
               return MillingPerpXY( ptLs, ptLe, vtLs) ;
           // Movimento nel piano generico
            else
               return MillingXYPlaneGen( ptLs, ptLe, vtLs) ;
         }
        // Movimento con Z non costante
         else {
           // Grandezze geometriche per selezione
            Vector3d vtMove = ptLe - ptLs ;
            Vector3d vtTest = vtMove - ( vtMove * vtLs) * vtLs ; 

            if ( abs( vtTest * X_AX) < EPS_SMALL && abs( vtTest * Y_AX) < EPS_SMALL)
               return MillingXYLongVert( ptLs, ptLe, vtLs) ;
            else {
               if ( abs( vtMove * vtLs) < EPS_SMALL) { 
                  if ( m_nToolType == 1)
                     return MillingXYCyl( ptLs, ptLe, vtLs) ;
                  else if ( m_nToolType == 2)
                     return MillingXYBall( ptLs, ptLe, vtLs) ;
               }
               else { 
                  if ( m_nToolType ==1)
                     return MillingXYPlusCyl( ptLs, ptLe, vtLs) ;
                  else if ( m_nToolType == 2)
                     return MillingXYPlusBall( ptLs, ptLe, vtLs) ;
               }
            }
         }
      }
   }
  // Altri casi, non gestiti
   return false ;
}

// Frese: cylindrical, ball-end e bull-nose //////////////////////////////////
// Versore utensile parallelo all'asse Z 
//----------------------------------------------------------------------------
bool
VolZmap::MillingDrillZ( const Point3d& ptLs, const Point3d& ptLe, const Vector3d& vtToolDir)
{
   if ( m_nToolType == 4)
      return ConusDrillingZ( ptLs, ptLe, vtToolDir) ;

  // Bounding box
   double dMinX = min( ptLs.x, ptLe.x) - m_dRadius ;
   double dMinY = min( ptLs.y, ptLe.y) - m_dRadius ;
   double dMaxX = max( ptLs.x, ptLe.x) + m_dRadius ;
   double dMaxY = max( ptLs.y, ptLe.y) + m_dRadius ;

  // Verifico interferisca con lo Zmap
   if ( dMaxX < EPS_SMALL  ||  dMinX > m_nNx * m_dStep - EPS_SMALL)
      return true ;
   if ( dMaxY < EPS_SMALL  ||  dMinY > m_nNy * m_dStep - EPS_SMALL)
      return true ;

  // Determino i limiti sugli indici
   unsigned int nStartI = ( dMinX < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinX / m_dStep)) ;
   unsigned int nEndI = ( dMaxX > m_nNx * m_dStep - EPS_SMALL ? m_nNx - 1 : static_cast <unsigned int> ( dMaxX / m_dStep)) ;
   unsigned int nStartJ = ( dMinY < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinY / m_dStep)) ;
   unsigned int nEndJ = ( dMaxY > m_nNy * m_dStep - EPS_SMALL ? m_nNy - 1 : static_cast <unsigned int> ( dMaxY / m_dStep)) ;

  // Determino quote estreme del tagliente            
   double dMax = ptLs.z - m_dHeight ;  // Volendo si può effettuare un controllo su quota punto finale e iniziale,
   double dMin = ptLe.z - m_dHeight ;  // ma non è necessaria dal momento che Subtract intervals scambia min con max se max < min 
   double dZCutBase = ptLs.z ; // Quota della base del tagliente nella posizione iniziale
   double dDeltaZ = ptLe.z - ptLs.z ; // Differenza delle quote fra le posizioni finale e iniziale della base del tagliente 

  // Limite sul quadrato del raggio
   double dSqRad = ( m_dRadius + EPS_SMALL) * ( m_dRadius + EPS_SMALL) ;

  // Ciclo sui punti nei limiti
   for ( unsigned int i = nStartI ; i <= nEndI ; ++ i) {
      double dX = ( i + 0.5) * m_dStep ;
      for ( unsigned int j = nStartJ ; j <= nEndJ ; ++ j) {
         double dY = ( j + 0.5) * m_dStep ;
        // punto
         Point3d ptQ( dX, dY, 0) ;
        // determino il quadrato della distanza
         double dSqDist = SqDistXY( ptQ, ptLe) ;
        // se distanza nei limiti, taglio
         if ( dSqDist < dSqRad)
            GetMinMaxZ( i, j, dZCutBase, dDeltaZ, dSqDist, vtToolDir) ;
      }
   }

   return true ; 
}

//----------------------------------------------------------------------------
bool
VolZmap::MillingPerpZ( const Point3d& ptLs, const Point3d& ptLe, const Vector3d& vtToolDir)  
{
   if ( m_nToolType == 4)
      return ConusPerpZ( ptLs, ptLe, vtToolDir) ;

  // Bounding box
   double dMinX = min( ptLs.x, ptLe.x) - m_dRadius ;
   double dMinY = min( ptLs.y, ptLe.y) - m_dRadius ;
   double dMaxX = max( ptLs.x, ptLe.x) + m_dRadius ;
   double dMaxY = max( ptLs.y, ptLe.y) + m_dRadius ;

  // Verifico interferisca con lo Zmap
   if ( dMaxX < EPS_SMALL  ||  dMinX > m_nNx * m_dStep - EPS_SMALL)
      return true ;
   if ( dMaxY < EPS_SMALL  ||  dMinY > m_nNy * m_dStep - EPS_SMALL)
      return true ;

  // Determino i limiti sugli indici
   unsigned int nStartI = ( dMinX < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinX / m_dStep)) ;
   unsigned int nEndI = ( dMaxX > m_nNx * m_dStep - EPS_SMALL ? m_nNx - 1 : static_cast <unsigned int> ( dMaxX / m_dStep)) ;
   unsigned int nStartJ = ( dMinY < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinY / m_dStep)) ;
   unsigned int nEndJ = ( dMaxY > m_nNy * m_dStep - EPS_SMALL ? m_nNy - 1 : static_cast <unsigned int> ( dMaxY / m_dStep)) ;

  // Determino quote estreme del tagliente            
   double dZCutBase = ptLs.z ; // Quota della base del tagliente nella posizione iniziale
   double dDeltaZ = ptLe.z - ptLs.z ; // Differenza delle quote fra le posizioni finale e iniziale della base del tagliente 

  // Limite sul quadrato del raggio
   double dSqRad = ( m_dRadius + EPS_SMALL) * ( m_dRadius + EPS_SMALL) ;

  // Segmento di movimento (nel piano griglia)
   Point3d ptStart( ptLs.x, ptLs.y, 0) ;
   Point3d ptEnd( ptLe.x, ptLe.y, 0) ;
   double dLen ;
   Vector3d vtDir ;
   DirDist( ptStart, ptEnd, vtDir, dLen) ;

  // Ciclo sui punti nei limiti
   for ( unsigned int i = nStartI ; i <= nEndI ; ++ i) {
      double dX = ( i + 0.5) * m_dStep ;
      for ( unsigned int j = nStartJ ; j <= nEndJ ; ++ j) {
         double dY = ( j + 0.5) * m_dStep ;
        // punto
         Point3d ptQ( dX, dY, 0) ;
        // determino il quadrato della distanza del punto dal segmento
         double dProiez = vtDir * ( ptQ - ptStart) ;
         if ( dProiez < 0)
            dProiez = 0 ;
         else if ( dProiez > dLen)
            dProiez = dLen ;
         Point3d ptMinDist = ptStart + vtDir * dProiez ;
         double dSqDist = SqDistXY( ptQ, ptMinDist) ;
        // se distanza nei limiti, taglio
         if ( dSqDist < dSqRad)
            GetMinMaxZ( i, j, dZCutBase, dDeltaZ, dSqDist, vtToolDir) ;
      }
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::MillingZ( const Point3d& ptLs, const Point3d& ptLe, const Vector3d& vtToolDir) {

   if ( m_nToolType == 4)
      if ( m_dRadius > m_dTipRadius)
         return ConusMillingZDr( ptLs, ptLe, vtToolDir) ;
      else
         return ConusMillingZSw( ptLs, ptLe, vtToolDir) ;

  // Setto il fattore per l'orientazione in z
   double dFactor = ( vtToolDir.z < 0 ?  1 : - 1) ;

  // Bounding box
   double dMinX = min( ptLs.x, ptLe.x) - m_dRadius ;
   double dMaxX = max( ptLs.x, ptLe.x) + m_dRadius ;
   double dMinY = min( ptLs.y, ptLe.y) - m_dRadius ;
   double dMaxY = max( ptLs.y, ptLe.y) + m_dRadius ;
   double dMinZ = min( min( ptLs.z, ptLs.z + dFactor * m_dHeight), min( ptLe.z, ptLe.z + dFactor * m_dHeight)) ;
   double dMaxZ = max( max( ptLs.z, ptLs.z + dFactor * m_dHeight), max( ptLe.z, ptLe.z + dFactor * m_dHeight)) ;

  // Verifico interferisca con lo Zmap
   if ( dMaxX < EPS_SMALL  ||  dMinX > m_nNx * m_dStep - EPS_SMALL)
      return true ;
   if ( dMaxY < EPS_SMALL  ||  dMinY > m_nNy * m_dStep - EPS_SMALL)
      return true ;
   if ( dMaxZ < m_dMinZ + EPS_SMALL || dMinZ > m_dMaxZ - EPS_SMALL)
      return true ;

  // Determino i limiti sugli indici
   unsigned int nStartI = ( dMinX < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinX / m_dStep)) ;
   unsigned int nEndI = ( dMaxX > m_nNx * m_dStep - EPS_SMALL ? m_nNx - 1 : static_cast <unsigned int> ( dMaxX / m_dStep)) ;
   unsigned int nStartJ = ( dMinY < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinY / m_dStep)) ;
   unsigned int nEndJ = ( dMaxY > m_nNy * m_dStep - EPS_SMALL ? m_nNy - 1 : static_cast <unsigned int> ( dMaxY / m_dStep)) ;
  //
   double dZCutBase = ptLs.z ; // Quota della base del tagliente nella posizione iniziale
   double dDeltaZ = ptLe.z - ptLs.z ; // Differenza delle quote fra le posizioni finale e iniziale della base del tagliente 

  // Limite sul quadrato del raggio
   double dSqRad = ( m_dRadius + EPS_SMALL) * ( m_dRadius + EPS_SMALL) ;

  // Segmento di movimento (nel piano griglia)
   Point3d ptStart( ptLs.x, ptLs.y, 0) ;
   Point3d ptEnd( ptLe.x, ptLe.y, 0) ;
   double dLen ;
   Vector3d vtDir ;
   DirDist( ptStart, ptEnd, vtDir, dLen) ;

  // Ciclo sui punti nei limiti
   for ( unsigned int i = nStartI ; i <= nEndI ; ++ i) {
      double dX = ( i + 0.5) * m_dStep ;
      for ( unsigned int j = nStartJ ; j <= nEndJ ; ++ j) {
         double dY = ( j + 0.5) * m_dStep ;
        // punto
         Point3d ptQ( dX, dY, 0) ;
        // determino il quadrato della distanza del punto dal segmento
         double dProj = vtDir * ( ptQ - ptStart) ;
         double dProiez ;

         if ( dProj < 0)
            dProiez = 0 ;
         else if ( dProj < dLen)
            dProiez = dProj ;
         else 
            dProiez = dLen ;

         Point3d ptMinDist = ptStart + vtDir * dProiez ;
         double dSqDist = SqDistXY( ptQ, ptMinDist) ;

        // se distanza nei limiti, taglio
         if ( dSqDist < dSqRad)
            GetMinMaxZGen( i, j, dProj, dSqDist, dLen, dZCutBase, dDeltaZ, vtToolDir) ;
      }
   }
   return true ;   
}

//----------------------------------------------------------------------------
inline bool
VolZmap::GetMinMaxZ( unsigned int nI, unsigned int nJ, double dZCutBase, double dDeltaZ, double dSqDist, const Vector3d& vtToolDir) {

  // Definisco variabili 
   double dSqRad = m_dRadius * m_dRadius ;
   double dFactor = ( vtToolDir.z < 0 ? - 1 : 1) ;
   double dZtip = ( vtToolDir.z < 0 ? dZCutBase + m_dHeight : dZCutBase - m_dHeight) ;

   double dMin ; double dMax ;
   
   // Caso utensile generico al momento non gestito
   if ( m_nToolType == 0)

      return false;
   // Caso Cylindrical Mill
   else if ( m_nToolType == 1) { 
    
      if ( abs(vtToolDir.z * dDeltaZ) < EPS_SMALL) {            // Non è meglio fare if (vt * delta < - Eps_small ) else if ( vt * delta < Eps) else ?
         dMin = min(dZCutBase, dZtip) ;
         dMax = max(dZCutBase, dZtip) ;
      }
      else if ( vtToolDir.z * dDeltaZ < 0) {
         dMin = min(dZtip, dZtip + dDeltaZ) ;
         dMax = max(dZtip, dZtip + dDeltaZ) ;
      }
      else {
         dMin = min(dZCutBase, dZCutBase + dDeltaZ) ;
         dMax = max(dZCutBase, dZCutBase + dDeltaZ) ;
      }
      return SubtractIntervals( nI, nJ, dMin, dMax) ; 
   }
   // Caso Ball-End Mill
   else if ( m_nToolType == 2) { 

      if ( abs(vtToolDir.z * dDeltaZ) < EPS_SMALL) {
         dMin = min(dZCutBase, dZtip + dFactor*(m_dRadius - sqrt(dSqRad - dSqDist))) ;
         dMax = max(dZCutBase, dZtip + dFactor*(m_dRadius - sqrt(dSqRad - dSqDist))) ;
      }
      else if ( vtToolDir.z * dDeltaZ < 0) {
         dMin = min(dZtip + dFactor*(m_dRadius - sqrt(dSqRad - dSqDist)), dZtip + dFactor*(m_dRadius - sqrt(dSqRad - dSqDist)) + dDeltaZ) ;// ocio
         dMax = max(dZtip + dFactor*(m_dRadius - sqrt(dSqRad - dSqDist)), dZtip + dFactor*(m_dRadius - sqrt(dSqRad - dSqDist)) + dDeltaZ) ;// ocio
      }
      else {
         dMin = min(dZCutBase, dZCutBase + dDeltaZ) ;
         dMax = max(dZCutBase, dZCutBase + dDeltaZ) ;
      }
      return SubtractIntervals( nI, nJ, dMin, dMax) ; 
   }
   // Caso Bull-Nose Mill
   else if ( m_nToolType == 3) {

      double dDeltaR = m_dRadius - m_dRCorner ;

      if ( dSqDist < dDeltaR*dDeltaR) {

         if ( abs(vtToolDir.z * dDeltaZ) < EPS_SMALL) {
            dMin = min(dZCutBase, dZtip) ;
            dMax = max(dZCutBase, dZtip) ;
         }
         else if ( vtToolDir.z * dDeltaZ < 0) {
            dMin = min(dZtip, dZtip + dDeltaZ) ;
            dMax = max(dZtip, dZtip + dDeltaZ) ;
         }
         else {
            dMin = min(dZCutBase, dZCutBase + dDeltaZ) ;
            dMax = max(dZCutBase, dZCutBase + dDeltaZ) ;
         }
         return SubtractIntervals( nI, nJ, dMin, dMax) ; 
      }
      else { 

         double dSqRadC = m_dRCorner*m_dRCorner ;
         double dSqd = dSqDist + dDeltaR*dDeltaR - 2*sqrt(dSqDist)*dDeltaR ; 

         if ( abs(vtToolDir.z * dDeltaZ) < EPS_SMALL) {
            dMin = min(dZCutBase, dZtip + dFactor*(m_dRCorner - sqrt(dSqRadC - dSqd))) ;
            dMax = max(dZCutBase, dZtip + dFactor*(m_dRCorner - sqrt(dSqRadC - dSqd))) ;
         }
         else if ( vtToolDir.z * dDeltaZ < 0) {
            dMin = min(dZtip + dFactor*(m_dRCorner - sqrt(dSqRadC - dSqd)), dZtip + dFactor*(m_dRCorner - sqrt(dSqRadC - dSqd)) + dDeltaZ) ;
            dMax = max(dZtip + dFactor*(m_dRCorner - sqrt(dSqRadC - dSqd)), dZtip + dFactor*(m_dRCorner - sqrt(dSqRadC - dSqd)) + dDeltaZ) ;
         }
         else {
            dMin = min(dZCutBase, dZCutBase + dDeltaZ) ;
            dMax = max(dZCutBase, dZCutBase + dDeltaZ) ;
         }
         return SubtractIntervals( nI, nJ, dMin, dMax) ; 
      }
   }
   // Caso di utensile inesistente ( m_nToolType fuori dai valori concessi)
   else
      return false ;
} 

//----------------------------------------------------------------------------
inline bool
VolZmap::GetMinMaxZGen( unsigned int nI, unsigned int nJ, double dProj, double dSqd, double dLenPath, double dZheight, double dDelta, const Vector3d& vtToolDir) {

  // Controllo sul tipo di utensile: se 0 utensile generico (al momento non gestito) se maggiore di 3 utensile fuori dai tipi consentiti
   if ( m_nToolType == 0 || m_nToolType > 3)
      return false ;

  // Definisco variabili quota della punta, minimo e massimo dell'intervallo da sottrarre
   double dZTip = ( vtToolDir.z < 0 ? dZheight + m_dHeight : dZheight - m_dHeight) ;
   double dMin, dMax ; double dStart ; // dStart è un parametro che esprime l'ascissa in cui la retta congiungente le due posizioni iniziale e finale 
   double dRSqDist1, dRSqDist2 ;       // di un punto del tagliente (nel sistema dell'asse dell'utensile nella posizione iniziale) assume la quota dZheight

  // Nei conti è comodo che dSqd assuma sempre il significato di distanza del punto dall'asse del movimento al quadrato
   if ( dProj < 0)
      dSqd = dSqd - dProj * dProj ;
   else if ( dProj > dLenPath)
      dSqd = dSqd - ( dProj - dLenPath) * ( dProj - dLenPath) ; 
  
  // Caso di cylindrical mill
   if ( m_nToolType == 1) {

      double dZ1, dZ2 ; 
     // Se lavora la punta
      if ( vtToolDir.z * dDelta < 0) {
         dZ1 = dZTip ;
         dZ2 = dZheight ;
      }
     // Se lavora il fondo
      else {
         dZ1 = dZheight ;
         dZ2 = dZTip ;
      }

      dStart = sqrt( m_dRadius * m_dRadius - dSqd) ;
      dRSqDist1 = dProj * dProj + dSqd ;
      dRSqDist2 = ( dLenPath - dProj) * ( dLenPath - dProj) + dSqd ;

      if ( dRSqDist1 < (m_dRadius + EPS_SMALL) * (m_dRadius + EPS_SMALL)) {

         dMin = min( max( dZ1 + ( dDelta * ( dProj + dStart)) / dLenPath, dZ1 + dDelta), dZ2) ;
         dMax = max( min( dZ1 + ( dDelta * ( dProj + dStart)) / dLenPath, dZ1 + dDelta), dZ2) ;
      }
      else if ( dRSqDist2 < (m_dRadius + EPS_SMALL) * (m_dRadius + EPS_SMALL)) {
            
         dMin = min( min(dZ2 + ( dDelta * ( dProj - dStart)) / dLenPath, dZ2 + dDelta), dZ1 + dDelta) ;
         dMax = max( max(dZ2 + ( dDelta * ( dProj - dStart)) / dLenPath, dZ2 + dDelta), dZ1 + dDelta) ;
      }
      else {

         dMin = min( dZ1 + ( dDelta * ( dProj + dStart)) / dLenPath, dZ2 + ( dDelta * ( dProj - dStart)) / dLenPath) ;
         dMax = max( dZ1 + ( dDelta * ( dProj + dStart)) / dLenPath, dZ2 + ( dDelta * ( dProj - dStart)) / dLenPath) ;
      }
   }
  // Caso di ball-end mill
   else if ( m_nToolType == 2) {         

      if ( dDelta < 0) {
            
         dDelta = - dDelta ;
         dProj = dLenPath - dProj ;
         dZheight = dZheight - dDelta ;
         dZTip = dZTip - dDelta ;
      }

      dStart = sqrt( m_dRadius * m_dRadius - dSqd) ;
      dRSqDist1 = dProj * dProj + dSqd ;                             
      dRSqDist2 = ( dLenPath - dProj) * ( dLenPath - dProj) + dSqd ;

      if ( vtToolDir.z > 0) {

        // Semi-asse ellisse 
         double dSemiAxMin = m_dRadius * sqrt( 1 - dLenPath * dLenPath / ( dLenPath * dLenPath + dDelta * dDelta)) ;
         double dSqrSemiAxMin = m_dRadius * m_dRadius * ( 1 - dLenPath * dLenPath / ( dLenPath * dLenPath + dDelta * dDelta)) ;
         double dXD2 = ( dProj - dLenPath) * ( dProj - dLenPath) ;

         if ( dRSqDist2 < m_dRadius * m_dRadius) {

            dMax = dZheight + dDelta ;

            double dTest = ( 1 - dSqd / ( m_dRadius * m_dRadius) > 0 ? sqrt( 1 - dSqd / ( m_dRadius * m_dRadius)) : 0) ;
          
            if ( dProj - dLenPath > dSemiAxMin * dTest) {

               double dSqrRad = ( dProj - dLenPath) * ( dProj - dLenPath) + dSqd ; 
               double dH = ( m_dRadius * m_dRadius - dSqrRad > 0 ? sqrt( m_dRadius * m_dRadius - dSqrRad) : 0) ;
               dMin = dZTip + m_dRadius + dDelta - dH ;
            }
            else if ( dProj > dSemiAxMin * dTest) {
              // Determino l'altezza del punto sull'ellisse da cui passa la retta
               double dPr0 = ( ( 1 - dSqd / ( m_dRadius * m_dRadius)) > 0 ? dLenPath + dSemiAxMin * sqrt( ( 1 - dSqd / (m_dRadius * m_dRadius))) : dLenPath) ;
               double dPar = ( m_dRadius * m_dRadius - dSqrSemiAxMin > 0 ? sqrt( m_dRadius * m_dRadius - dSqrSemiAxMin) : 0) ;
               double dZ0 = ( ( dPr0 - dLenPath) * dPar) / dSemiAxMin ; 
              
               dMin = dZTip + dDelta + m_dRadius - dZ0 + ( dDelta / dLenPath) * ( dProj - dPr0) ;  
            }
            else {

               double dSqrRad = dProj * dProj + dSqd ; 
               double dH = ( m_dRadius * m_dRadius - dSqrRad > 0 ? sqrt( m_dRadius * m_dRadius - dSqrRad) : 0) ;               
               dMin = dZTip + m_dRadius - dH ;
            }
         }
         else { 

            dMax = dZheight + ( dDelta * ( dProj + dStart)) / dLenPath ;

            double dTest = ( 1 - dSqd / ( m_dRadius * m_dRadius) > 0 ? sqrt( 1 - dSqd / ( m_dRadius * m_dRadius)) : 0) ;

            if ( dProj  > dSemiAxMin * dTest) {            
              // Determino l'altezza del punto sull'ellisse da cui passa la retta
               double dPr0 = ( ( 1 - dSqd / ( m_dRadius * m_dRadius)) > 0 ? dLenPath + dSemiAxMin * sqrt( ( 1 - dSqd / (m_dRadius * m_dRadius))) : dLenPath) ;
               double dPar = ( m_dRadius * m_dRadius - dSqrSemiAxMin > 0 ? sqrt( m_dRadius * m_dRadius - dSqrSemiAxMin) : 0) ;
               double dZ0 = ( ( dPr0 - dLenPath) * dPar) / dSemiAxMin ; 
              
               dMin = dZTip + dDelta + m_dRadius - dZ0 + ( dDelta / dLenPath) * ( dProj - dPr0) ;    
            }
            else {

               double dSqrRad = dProj * dProj + dSqd ; 
               double dH = ( m_dRadius * m_dRadius - dSqrRad > 0 ? sqrt( m_dRadius * m_dRadius - dSqrRad) : 0) ;               
               
               dMin = dZTip + m_dRadius - dH ;
            }
         }         
      }     
      else {

        // Semi-asse ellisse 
         double dSemiAxMin = m_dRadius * sqrt( 1 - dLenPath * dLenPath / ( dLenPath * dLenPath + dDelta * dDelta)) ;
         double dSqrSemiAxMin = m_dRadius * m_dRadius * ( 1 - dLenPath * dLenPath / ( dLenPath * dLenPath + dDelta * dDelta)) ;
         double dXD2 = ( dProj - dLenPath) * ( dProj - dLenPath) ;

         if ( dRSqDist1 < m_dRadius * m_dRadius) {

            dMin = dZheight ;

            double dTest = ( 1 - dSqd / ( m_dRadius * m_dRadius) > 0 ? sqrt( 1 - dSqd / ( m_dRadius * m_dRadius)) : 0) ;
          
            if ( dProj < - dSemiAxMin * dTest) {

               double dSqrRad = dProj * dProj + dSqd ; 
               double dH = ( m_dRadius * m_dRadius - dSqrRad > 0 ? sqrt( m_dRadius * m_dRadius - dSqrRad) : 0) ;
               
               dMax = dZTip - m_dRadius + dH ;
            }
            else if ( dProj - dLenPath < - dSemiAxMin * dTest) {
              // Determino l'altezza del punto sull'ellisse da cui passa la retta
               double dPr0 = ( ( 1 - dSqd / ( m_dRadius * m_dRadius)) > 0 ? dSemiAxMin * sqrt( ( 1 - dSqd / (m_dRadius * m_dRadius))) : 0) ;
               double dPar = ( m_dRadius * m_dRadius - dSqrSemiAxMin > 0 ? sqrt( m_dRadius * m_dRadius - dSqrSemiAxMin) : 0) ;
               double dZ0 = ( dPr0 * dPar) / dSemiAxMin ; 
              
               dMax = dZTip - m_dRadius + dZ0 + ( dDelta / dLenPath) * ( dProj + dPr0) ;  
            }
            else {

               double dSqrRad = dProj * dProj + dSqd ; 
               double dH = ( m_dRadius * m_dRadius - dSqrRad > 0 ? sqrt( m_dRadius * m_dRadius - dSqrRad) : 0) ;               
               dMax = dZTip - m_dRadius + dDelta + dH ;
            }
         }
         else { 

            dMin = dZheight + ( dDelta * ( dProj - dStart)) / dLenPath ;

            double dTest = ( 1 - dSqd / ( m_dRadius * m_dRadius) > 0 ? sqrt( 1 - dSqd / ( m_dRadius * m_dRadius)) : 0) ;

            if ( dProj - dLenPath < - dSemiAxMin * dTest) {            
              // Determino l'altezza del punto sull'ellisse da cui passa la retta
               double dPr0 = ( ( 1 - dSqd / ( m_dRadius * m_dRadius)) > 0 ? dSemiAxMin * sqrt( ( 1 - dSqd / (m_dRadius * m_dRadius))) : 0) ;
               double dPar = ( m_dRadius * m_dRadius - dSqrSemiAxMin > 0 ? sqrt( m_dRadius * m_dRadius - dSqrSemiAxMin) : 0) ;
               double dZ0 = ( dPr0 * dPar) / dSemiAxMin ; 
              
               dMax = dZTip - m_dRadius + dZ0 + ( dDelta / dLenPath) * ( dProj + dPr0) ;    
            }
            else {

               double dSqrRad = ( dProj - dLenPath) * ( dProj - dLenPath) + dSqd ; 
               double dH = ( m_dRadius * m_dRadius - dSqrRad > 0 ? sqrt( m_dRadius * m_dRadius - dSqrRad) : 0) ;               
               dMax = dZTip + dDelta - m_dRadius + dH ;
            }
         }        
      }
   }
  // Caso di bull-nose mill
   else
      return true ;

   return SubtractIntervals( nI, nJ, dMin, dMax) ;
}

// Versore utensile nel piano XY 
//----------------------------------------------------------------------------
bool
VolZmap::MillingDrillXY( const Point3d& ptLs, const Point3d& ptLe, const Vector3d& vtToolDir) {
      
   // Punti per la determinazione del materiale asportato
    Vector3d vtMove = ptLe - ptLs ;
   // Punti di riferimento dell'asportazione
    Point3d ptLNs ;
    Point3d ptLNe ;
   // Parametro relativo all'utensile altezza della parte non cilindrica e fattore determinante
   // la lunghezza della parte lavorata a seconda che lavori la punta o il fondo 
    double dCylH ; double dFactor ;                                                     

   // Caso utensile generico (al momento non gestito)
    if ( m_nToolType == 0)
       return false ;
   // Caso Cylindrical Mill
    else if ( m_nToolType == 1)
       dCylH = m_dHeight ;
   // Caso Ball-end Mill
    else if ( m_nToolType == 2)
       dCylH = m_dHeight - m_dRadius ;
   // Caso Bull-nose Mill
    else if ( m_nToolType == 3)
       dCylH = m_dHeight - m_dRCorner ;

   // Normalizzo tale vettore e ne determino la lunghezza: 
    double dLenPath = vtMove.Len() ; vtMove.Normalize() ;

   // Prodotto scalare fra versore direzione utensile e direzione movimento
    double dScProd = vtMove * vtToolDir ;
   // Se lavora la punta
    if ( dScProd < 0) {
      // Trovo i punti di riferimento per la lavorazione
       ptLNs = ptLs + ( vtMove * dCylH) ;
       ptLNe = ptLe + ( vtMove * dCylH) ;
       dFactor = 1 ;
    } 
   // Se lavora il fondo
    else {
      ptLNs = ptLs ;
      ptLNe = ptLe ;
      dFactor = 0 ;
    }

   // Quota z dei punti iniziale e finale
    double dHz = ptLNs.z ;

   // Bounding box
    double dMinX = min( ptLNs.x, ptLNe.x) - m_dRadius ;
    double dMaxX = max( ptLNs.x, ptLNe.x) + m_dRadius ;
    double dMinY = min( ptLNs.y, ptLNe.y) - m_dRadius ;
    double dMaxY = max( ptLNs.y, ptLNe.y) + m_dRadius ;
    double dMinZ = dHz - m_dRadius ;
    double dMaxZ = dHz + m_dRadius ;

   // Verifico se il movimento intersca lo Zmap
    if ( dMaxX < EPS_SMALL || dMinX > m_nNx * m_dStep - EPS_SMALL)
       return true ;
    if ( dMaxY < EPS_SMALL || dMinY > m_nNy * m_dStep - EPS_SMALL)
       return true ;
    if ( dMaxZ < m_dMinZ + EPS_SMALL || dMinZ > m_dMaxZ - EPS_SMALL)
       return true ;

   // Determino i limiti sugli indici
    unsigned int nStartI = ( dMinX < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinX / m_dStep)) ;
    unsigned int nEndI = ( dMaxX > m_nNx * m_dStep - EPS_SMALL ? m_nNx - 1 : static_cast <unsigned int> ( dMaxX / m_dStep)) ; 
    unsigned int nStartJ = ( dMinY < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinY / m_dStep)) ;
    unsigned int nEndJ = ( dMaxY > m_nNy * m_dStep - EPS_SMALL ? m_nNy - 1 : static_cast <unsigned int> ( dMaxY / m_dStep)) ;

   //Proietto ptLNs sul piano XY:
    Point3d ptStart( ptLNs.x, ptLNs.y, 0) ; 

   //Ciclo sui punti
     for ( unsigned int i = nStartI ; i <= nEndI ; ++ i)
        for ( unsigned int j = nStartJ ; j <= nEndJ ; ++ j) {
          // Punto da valutare 
           Point3d ptC( (i + 0.5)*m_dStep, (j + 0.5)*m_dStep, 0) ;
          // Vettore spostamento da ptLe a ptC 
           Vector3d vtC = ( ptC - ptStart) ;
          // Componenti parallela e perpendicolare a vtMove
          // Vettore ortogonale a vtMove
           Vector3d vtOrt = vtMove ;
          // Ruoto vtOrt affinché sia ortogonale
           vtOrt.Rotate( Z_AX, -90) ;

           double dProj = vtC * vtOrt ;
           Vector3d vtPara = vtOrt * dProj ; // Parallelo alla perpendicolare al movimento
           Vector3d vtPerp = vtC - vtPara ;   // Perpendicolare alla perpendicolare al movimento ( serve per unire le getminmax di drill e perp)
          // Distanza di ptC dall'asse dell'untensile
           double dSqDist = vtPerp.SqLen() ;
           double dLimitMill = dLenPath + dFactor * ( m_dHeight - dCylH) ;
          // Se dTestProj è positivo è vtC è dalla parte giusta
           double dTestProj = vtC * vtMove ;

           if ( dTestProj > 0 && dSqDist < dLimitMill * dLimitMill)
              if ( dProj > - m_dRadius && dProj < m_dRadius)
                 GetMinMaxXY( i, j, dProj, dHz, dSqDist, 0, dLenPath, dScProd) ;        
        }
    return true ;
}   

//----------------------------------------------------------------------------
bool
VolZmap::MillingPerpXY( const Point3d& ptLs, const Point3d& ptLe, const Vector3d& vtToolDir) {
   
  // Determinazione delle posizioni iniziali e finali della punta dell'utensile
   Point3d ptTLs = ptLs - vtToolDir * m_dHeight ;
   Point3d ptTLe = ptLe - vtToolDir * m_dHeight ;

  // Quota Z
   double dZH = ptLs.z ;
   
  // Estremi Bounding box
   double dMinX = min( min( ptLs.x, ptLe.x), min( ptTLs.x, ptTLe.x)) - m_dRadius ;
   double dMaxX = max( max( ptLs.x, ptLe.x), max( ptTLs.x, ptTLe.x)) + m_dRadius ;
   double dMinY = min( min( ptLs.y, ptLe.y), min( ptTLs.y, ptTLe.y)) - m_dRadius ;
   double dMaxY = max( max( ptLs.y, ptLe.y), max( ptTLs.y, ptTLe.y)) + m_dRadius ;
   double dMinZ = dZH - m_dRadius ;
   double dMaxZ = dZH + m_dRadius ;
  
  // Verifica dell'interferenza dell'utensile con lo Zmap
   if ( dMaxX < EPS_SMALL || dMinX > m_nNx * m_dStep - EPS_SMALL)
      return true ;
   if ( dMaxY < EPS_SMALL || dMinY > m_nNy * m_dStep - EPS_SMALL)
      return true ;
   if ( dMaxZ < m_dMinZ + EPS_SMALL || dMinZ > m_dMaxZ - EPS_SMALL)
      return true ;

  // Determinazione limiti sugli indici
   unsigned int nStartI = ( dMinX < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinX / m_dStep)) ;
   unsigned int nEndI = ( dMaxX > m_nNx * m_dStep - EPS_SMALL ? m_nNx - 1 : static_cast <unsigned int> ( dMaxX / m_dStep)) ;
   unsigned int nStartJ = ( dMinY < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinY / m_dStep)) ;
   unsigned int nEndJ = ( dMaxY > m_nNy * m_dStep - EPS_SMALL ? m_nNy - 1 : static_cast <unsigned int> ( dMaxY / m_dStep)) ;

  // Proiezione sul piano XY delle grandezze di interesse
   Point3d ptStart( ptLs.x, ptLs.y, 0) ;
   Point3d ptEnd( ptLe.x, ptLe.y, 0) ;
   Vector3d vtMove = ptEnd - ptStart ;
   double dLenPath = vtMove.Len() ;
  // Normalizzo il vettore vtMove congiungente le posizioni iniziale e finale della base dell'utensile
   vtMove.Normalize() ;

   for ( unsigned int i = nStartI ; i <= nEndI ; ++ i)
      for ( unsigned int j = nStartJ ; j <= nEndJ ; ++ j) {
        // Punto su cui ciclare e vettore congiungente posizione iniziale della base a tale punto
         Point3d ptC( (i + 0.5) * m_dStep, (j + 0.5) * m_dStep, 0) ;
         Vector3d vtC = ptC - ptStart ;
        // Proiezione di vtC sulla direzione del movimento
         double dProj = vtC * vtMove ;
        // Componente di vtC ortogonale al movimento
         Vector3d vtPerp = vtC - vtMove * dProj ;
        // Lunghezza quadrata del precedente vettore
         double dSqDist = vtPerp.SqLen() ;

         if ( dProj > - m_dRadius && dProj < dLenPath + m_dRadius && vtPerp * vtToolDir < 0 && dSqDist < m_dHeight * m_dHeight)
            GetMinMaxXY( i, j, dProj, dZH, dSqDist, dLenPath, 0, 0) ;   
      }   
   return true ;
}

//----------------------------------------------------------------------------
bool 
VolZmap::MillingXYPlaneGen( const Point3d& ptLs, const Point3d& ptLe, const Vector3d& vtToolDir) {

   double dMinZ = min( ptLs.z, ptLe.z) - m_dRadius ;
   double dMaxZ = max( ptLs.z, ptLe.z) + m_dRadius ;

  // Prima verifica sull'interferenza dell'utensile con lo Zmap
   if ( dMaxZ < m_dMinZ + EPS_SMALL || dMinZ > m_dMaxZ - EPS_SMALL)
      return true ;

   Point3d ptI = ptLs ;  Point3d ptF = ptLe ;

  // Quote dei punti ptI e ptF
   double dZ = ptI.z ;

   Point3d ptIT = ptI - vtToolDir * m_dHeight ;
   Point3d ptFT = ptF - vtToolDir * m_dHeight ;

  // Bounding box 
   double dMinX = min( min( ptI.x, ptIT.x), min( ptF.x, ptFT.x)) - m_dRadius ;
   double dMaxX = max( max( ptI.x, ptIT.x), max( ptF.x, ptFT.x)) + m_dRadius ;
   double dMinY = min( min( ptI.y, ptIT.y), min( ptF.y, ptFT.y)) - m_dRadius ;
   double dMaxY = max( max( ptI.y, ptIT.y), max( ptF.y, ptFT.y)) + m_dRadius ;

  // Seconda verifica dell'interferenza dell'utensile con lo Zmap
   if ( dMaxX < EPS_SMALL || dMinX > m_nNx * m_dStep - EPS_SMALL)
      return true ;
   if ( dMaxY < EPS_SMALL || dMinY > m_nNy * m_dStep - EPS_SMALL)
      return true ;

   Vector3d vtMove = ptF - ptI ; 

   if ( vtMove * vtToolDir > 0) {
      Point3d ptTemp = ptI ;
      ptI = ptF ;
      ptF = ptTemp ;
      vtMove = - vtMove ;
   }

   Vector3d vtMoveOrt = vtMove - ( vtMove * vtToolDir) * vtToolDir ;     double dLen2 = vtMoveOrt.Len() ;
   Vector3d vtMoveLong = ( vtMove * vtToolDir) * vtToolDir ;             double dLen1 = vtMoveLong.Len() ;

  // Determino sistema di riferimento del movimento: costruisco un sistrema di vettori ortonormali e destrorsi spiccati dal punto iniziale del movimento
   Vector3d vtV1 = vtToolDir ; 
   Vector3d vtV2 = vtMoveOrt ; vtV2.Normalize() ;
  
  // Punti iniziale e finale proiettati sul piano
   Point3d ptIxy( ptI.x, ptI.y, 0) ;  Point3d ptFxy( ptF.x, ptF.y, 0) ;

  // Determinazione limiti sugli indici
   unsigned int nStartI = ( dMinX < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinX / m_dStep)) ;
   unsigned int nEndI = ( dMaxX > m_nNx * m_dStep - EPS_SMALL ? m_nNx - 1 : static_cast <unsigned int> ( dMaxX / m_dStep)) ;
   unsigned int nStartJ = ( dMinY < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinY / m_dStep)) ;
   unsigned int nEndJ = ( dMaxY > m_nNy * m_dStep - EPS_SMALL ? m_nNy - 1 : static_cast <unsigned int> ( dMaxY / m_dStep)) ;

   // Ciclo
   for ( unsigned int i = nStartI ; i <= nEndI ; ++ i) 
    for ( unsigned int j = nStartJ ; j <= nEndJ ; ++ j) {

       if ( m_nToolType == 1) {

          Point3d ptC( ( i + 0.5) * m_dStep, ( j + 0.5) * m_dStep, 0) ;

          Vector3d vtC = ptC - ptIxy ;

          double dProj1 = vtC * vtV1 ;
          double dProj2 = vtC * vtV2 ;

          GetMMPlaneGenCyl( i, j, dZ, dLen1, dLen2, dProj1, dProj2) ;
       }
       else if ( m_nToolType == 2) 

          GetMMPlaneGenBall( i, j, dZ, dLen1, dLen2, ptIxy, vtMove, vtV1, vtV2) ; 
    }
    return true ;
}

//----------------------------------------------------------------------------
inline bool
VolZmap::GetMMPlaneGenCyl( unsigned int i, unsigned int j, double dZ, double dLen1, double dLen2, double dProj1, double dProj2) {

   double dMin, dMax ;

   if ( dProj2 > - m_dRadius && dProj2 < 0) {

      if ( dProj1 < 0 && dProj1 > - m_dHeight) {

         double dH = sqrt( m_dRadius * m_dRadius - dProj2 * dProj2) ;

         dMin = dZ - dH ; dMax = dZ + dH ;

         SubtractIntervals( i, j, dMin, dMax) ;    
      }
      else if ( dProj1 <= - m_dHeight) {
          
         if ( ( dProj2 < dLen2 - m_dRadius && dProj1 > - m_dHeight - ( dLen1 / dLen2) * ( dProj2 + m_dRadius))
              || ( dProj2 >= dLen2 - m_dRadius && dProj1 > - m_dHeight - dLen1)) {   // In questo costrutto if-else non c'è bisogno di specificare nient'altro perché già siamo nella regione - m_dRadius < dProj2 < 0

            double dPar = m_dHeight + ( dLen1 / dLen2) * ( dProj2 + m_dRadius) + dProj1 ;
            double dL = m_dRadius - ( dLen2 / dLen1) * dPar ;
            double dH = ( m_dRadius * m_dRadius - dL * dL > 0 ? sqrt( m_dRadius * m_dRadius - dL * dL) : 0) ;

            dMin = dZ - dH ; dMax = dZ + dH ;

            SubtractIntervals( i, j, dMin, dMax) ;
         }
      }
   }
   else if ( dProj2 >= 0 && dProj2 < m_dRadius) {  

      if ( dProj1 < 0 && ( ( dProj2 < dLen2 && dProj1 > - ( dLen1 / dLen2) * dProj2) || ( dProj2 >= dLen2 && dProj1 > - dLen1))) {

         double dPar = dProj2 + ( dLen2 / dLen1) * dProj1 ;
         /* Oppure
         double dPar1 = ( dLen1 / dLen2) * dProj2 + dProj1 ;
         double dPar2 = ( dLen2 / dLen1) * dPar1 ; */ 
         double dH = ( m_dRadius * m_dRadius - dPar * dPar > 0 ? sqrt( m_dRadius * m_dRadius - dPar * dPar) : 0) ;

         dMin = dZ - dH ; dMax = dZ + dH ;
             
         SubtractIntervals( i, j, dMin, dMax) ;
      }
      else if ( dProj2 < dLen2 && dProj1 <= - ( dLen1 / dLen2) * dProj2) {

         if ( dProj1 > - ( dLen1 / dLen2) * dProj2 - m_dHeight) {
             
            dMin = dZ - m_dRadius ; dMax = dZ + m_dRadius ;

            SubtractIntervals( i, j, dMin, dMax) ;
         }
         else if ( dProj1 <= - ( dLen1 / dLen2) * dProj2 - m_dHeight) {  

            if ( ( dProj2 < dLen2 - m_dRadius && dProj1 > - ( dLen1 / dLen2) * ( dProj2 + m_dRadius) - m_dHeight)
                 || ( dProj2 >= dLen2 - m_dRadius && dProj1 > - m_dHeight - dLen1)) {     // modificato qui

               double dPar = m_dHeight + ( dLen1 / dLen2) * ( dProj2 + m_dRadius) + dProj1 ;
               double dL = m_dRadius - ( dLen2 / dLen1) * dPar ;
               double dH = ( m_dRadius * m_dRadius - dL * dL > 0 ? sqrt( m_dRadius * m_dRadius - dL * dL) : 0) ;

               dMin = dZ - dH ; dMax = dZ + dH ;

               SubtractIntervals( i, j, dMin, dMax) ;
            }
         }
      }
      else if ( dProj2 >= dLen2 && dProj1 < - dLen1) {

         if ( dProj1 > - m_dHeight - dLen1) {
             
            double dH = sqrt( m_dRadius * m_dRadius - ( dProj2 - dLen2) * ( dProj2 - dLen2)) ;

            dMin = dZ - dH ; dMax = dZ + dH ;

            SubtractIntervals( i, j, dMin, dMax) ;
         }
      }
   }
   else if ( dProj2 >= m_dRadius && dProj2 < dLen2) {
       
      if ( dProj1 < - ( dLen1 / dLen2) * ( dProj2 - m_dRadius) && dProj1 > - ( dLen1 / dLen2) * dProj2) {
          
         double dL = ( dLen2 / dLen1) * ( ( dLen1 / dLen2) * dProj2 + dProj1) ; 
         double dH = ( m_dRadius * m_dRadius - dL * dL > 0 ? sqrt( m_dRadius * m_dRadius - dL * dL) : 0) ;
         
         dMin = dZ - dH ; dMax = dZ + dH ;

         SubtractIntervals( i, j, dMin, dMax) ;
      }
      else if ( dProj1 <= - ( dLen1 / dLen2) * dProj2 && dProj1 > - ( dLen1 / dLen2) * dProj2 - m_dHeight) {
          
         dMin = dZ - m_dRadius ;   dMax = dZ + m_dRadius ;

         SubtractIntervals( i, j, dMin , dMax) ;
      }
      else if ( dProj1 <= - ( dLen1 / dLen2) * dProj2 - m_dHeight) {
          
         if ( ( dProj2 < dLen2 - m_dRadius && dProj1 > - ( dLen1 / dLen2) * ( dProj2 + m_dRadius) - m_dHeight) 
              || ( dProj2 >= dLen2 - m_dRadius && dProj1 > - m_dHeight - dLen1)) {

            double dPar = m_dHeight + ( dLen1 / dLen2) * ( dProj2 + m_dRadius) + dProj1 ;
            double dL = m_dRadius - ( dLen2 / dLen1) * dPar ;
            double dH = ( m_dRadius * m_dRadius - dL * dL > 0 ? sqrt( m_dRadius * m_dRadius - dL * dL) : 0) ;

            dMin = dZ - dH ; dMax = dZ + dH ;

            SubtractIntervals( i, j, dMin, dMax) ;
         }
      }
   }  
   else if ( dProj2 >=  dLen2 && dProj2 < dLen2 + m_dRadius) {
       
      if ( dProj1 < -  ( dLen1 / dLen2) * ( dProj2 - m_dRadius) && dProj1 >= - dLen1) {

         double dL = ( dLen2 / dLen1) * ( ( dLen1 / dLen2) * dProj2 + dProj1) ; 
         double dH = ( m_dRadius * m_dRadius - dL * dL > 0 ? sqrt( m_dRadius * m_dRadius - dL * dL) : 0) ;

         dMin = dZ - dH ; dMax = dZ + dH ;

         SubtractIntervals( i, j, dMin, dMax) ;
      }
      else if ( dProj1 < - dLen1 && dProj1 > - dLen1 - m_dHeight) {
          
         double dL = dProj2 - dLen2 ;
         double dH = sqrt( m_dRadius * m_dRadius - dL * dL) ;

         dMin = dZ - dH ; dMax = dZ + dH ;

         SubtractIntervals( i, j, dMin, dMax) ;
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
inline bool
VolZmap::GetMMPlaneGenBall( unsigned int i, unsigned int j, double dZ, double dLen1, double dLen2, Point3d ptIxy, Vector3d vtMove, Vector3d vtV1, Vector3d vtV2) {

   double dMin, dMax ;

   Point3d ptC( ( i + 0.5) * m_dStep, ( j + 0.5) * m_dStep, 0) ;

   Vector3d vtC = ptC - ptIxy ;

   double dProj1 = vtC * vtV1 ;     // vtV1, vtV2 sono paralleli al piano
   double dProj2 = vtC * vtV2 ;

   double dCH = m_dHeight - m_dRadius ;
   double dLMove = vtMove.LenXY() ;

   Point3d ptCS = ptIxy - dCH * vtV1 ;   Point3d ptCE = ptCS + vtMove ;   // vtMove è orizzontale

   Vector3d vtCS = ptC - ptCS ;   Vector3d vtCE = ptC - ptCE ;

   double dProjMove = ( vtCS * vtMove) / dLMove ;

   Vector3d vtCSP = vtCS - ( ( vtCS * vtMove) / ( dLMove * dLMove)) * vtMove ;

   double dSQDist = vtCSP.SqLenXY() ;
   double dSQDistS = vtCS.SqLenXY() ; 
   double dSQDistE = vtCE.SqLenXY() ;

  // parte cilindrica

   if ( dProj2 > - m_dRadius && dProj2 < 0) {
   
      if ( dProj1 > - dCH && dProj1 < 0) {
      
         double dH = sqrt( m_dRadius * m_dRadius - dProj2 * dProj2) ;

         dMin = dZ - dH ; dMax = dZ + dH ;

         SubtractIntervals( i, j, dMin, dMax) ;
      }
   }
   else if ( dProj2 >= 0 && dProj2 < m_dRadius) {
   
      if ( dProj1 < 0 && ( ( dProj2 < dLen2 && dProj1 > - ( dLen1 / dLen2) * dProj2) || ( dProj2 >= dLen2 && dProj1 > - dLen1))) {

         double dPar = dProj2 + ( dLen2 / dLen1) * dProj1 ;
         /* Oppure
         double dPar1 = ( dLen1 / dLen2) * dProj2 + dProj1 ;
         double dPar2 = ( dLen2 / dLen1) * dPar1 ; */ 
         double dH = ( m_dRadius * m_dRadius - dPar * dPar > 0 ? sqrt( m_dRadius * m_dRadius - dPar * dPar) : 0) ;

         dMin = dZ - dH ; dMax = dZ + dH ;
             
         SubtractIntervals( i, j, dMin, dMax) ;
      }
      else if ( dProj2 < dLen2 && dProj1 <= - ( dLen1 / dLen2) * dProj2 && dProj1 > - ( dLen1 / dLen2) * dProj2 - dCH) {

         dMin = dZ - m_dRadius ;  dMax = dZ + m_dRadius ;

         SubtractIntervals( i, j, dMin, dMax) ;
      }
      else if ( dProj2 >= dLen2 && dProj1 <= - dLen1 && dProj1 > - dLen1 - dCH) {
      
         double dL = dProj2 - dLen2 ;
         double dH = sqrt( m_dRadius * m_dRadius - dL * dL) ;

         dMin = dZ - dH ;  dMax = dZ + dH ;

         SubtractIntervals( i, j, dMin, dMax) ;
      }
   }
   else if ( dProj2 >= m_dRadius && dProj1 < - ( dLen1 / dLen2) * ( dProj2 - m_dRadius)) {

      if ( dProj2 < dLen2 && dProj1 > - ( dLen1 / dLen2) * dProj2) {

         double dPar = dProj2 + ( dLen2 / dLen1) * dProj1 ; 
         double dH = ( m_dRadius * m_dRadius - dPar * dPar > 0 ? sqrt( m_dRadius * m_dRadius - dPar * dPar) : 0) ;

         dMin = dZ - dH ; dMax = dZ + dH ;
             
         SubtractIntervals( i, j, dMin, dMax) ;
      }
      else if ( dProj2 < dLen2 && dProj1 <= - ( dLen1 / dLen2) * dProj2 && dProj1 > - ( dLen1 / dLen2) * dProj2 - dCH) {
      
         dMin = dZ - m_dRadius ;  dMax = dZ + m_dRadius ;

         SubtractIntervals( i, j, dMin, dMax) ;
      }
      else if ( dProj2 >= dLen2 && dProj2 < dLen2 + m_dRadius) {
      
         if ( dProj1 > - dLen1) {
         
            double dPar = dProj2 + ( dLen2 / dLen1) * dProj1 ; 
            double dH = ( m_dRadius * m_dRadius - dPar * dPar > 0 ? sqrt( m_dRadius * m_dRadius - dPar * dPar) : 0) ;

            dMin = dZ - dH ; dMax = dZ + dH ;
             
            SubtractIntervals( i, j, dMin, dMax) ;
         }
         else if ( dProj1 <= - dLen1 && dProj1 > - dLen1 - dCH) {
         
            double dL = dProj2 - dLen2 ;
            double dH = sqrt( m_dRadius * m_dRadius - dL * dL) ;

            dMin = dZ - dH ; dMax = dZ + dH ;

            SubtractIntervals( i, j, dMin, dMax) ;
         }
      }
   }
  // parte non cilindrica
  
   
   if ( dProjMove > - m_dRadius && dProjMove < 0) {

      if ( dSQDistS < m_dRadius * m_dRadius) {
      
         double dH = sqrt( m_dRadius * m_dRadius - dSQDistS) ;

         dMin = dZ - dH ; dMax = dZ + dH ;

         SubtractIntervals( i, j, dMin, dMax) ;
      }
   }
   else if ( dProjMove >= 0 && dProjMove < dLMove) {

      if ( dSQDist < m_dRadius * m_dRadius) {
      
         double dH = sqrt( m_dRadius * m_dRadius - dSQDist) ;

         dMin = dZ - dH ; dMax = dZ + dH ;

         SubtractIntervals( i, j, dMin, dMax) ;
      }
   }
   else {

      if ( dSQDistE < m_dRadius * m_dRadius) {
         
         double dH = sqrt( m_dRadius * m_dRadius - dSQDistE) ;

         dMin = dZ - dH ; dMax = dZ + dH ;

         SubtractIntervals( i, j, dMin, dMax) ;
      }
   }
   
   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::MillingXYLongVert( const Point3d& ptLs, const Point3d& ptLe, const Vector3d& vtToolDir) {

  // Determino posizioni punte
   Point3d ptTLs = ptLs - vtToolDir * m_dHeight ;
   Point3d ptTLe = ptLe - vtToolDir * m_dHeight ;

  // Variabili utili
   double dDelta = ptLe.z - ptLs.z ;  double dMin, dMax ;

  // Bounding box
   double dMinX = min( min( ptLs.x, ptLe.x), min( ptTLs.x, ptTLe.x)) - m_dRadius ;
   double dMaxX = max( max( ptLs.x, ptLe.x), max( ptTLs.x, ptTLe.x)) + m_dRadius ;
   double dMinY = min( min( ptLs.y, ptLe.y), min( ptTLs.y, ptTLe.y)) - m_dRadius ;
   double dMaxY = max( max( ptLs.y, ptLe.y), max( ptTLs.y, ptTLe.y)) + m_dRadius ;
   double dMinZ = min( ptLs.z, ptLs.z + dDelta) - m_dRadius ;
   double dMaxZ = max( ptLs.z, ptLs.z + dDelta) + m_dRadius ;

  // Controllo su interferenza utensile-Zmap
   if ( dMaxX < EPS_SMALL || dMinX > m_nNx * m_dStep - EPS_SMALL)
      return true ;
   if ( dMaxY < EPS_SMALL || dMinY > m_nNy * m_dStep - EPS_SMALL)
      return true ;
   if ( dMaxZ < m_dMinZ + EPS_SMALL || dMinZ > m_dMaxZ - EPS_SMALL)
      return true ;

  // Determinazione limiti sugli indici
   unsigned int nStartI = ( dMinX < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinX / m_dStep)) ;
   unsigned int nEndI = ( dMaxX > m_nNx * m_dStep - EPS_SMALL ? m_nNx - 1 : static_cast <unsigned int> ( dMaxX / m_dStep)) ;
   unsigned int nStartJ = ( dMinY < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinY / m_dStep)) ;
   unsigned int nEndJ = ( dMaxY > m_nNy * m_dStep - EPS_SMALL ? m_nNy - 1 : static_cast <unsigned int> ( dMaxY / m_dStep)) ;

  // Definizione delle grandezze geometriche essenziali
   Vector3d vtMove = ptLe - ptLs ;
   Vector3d vtPlaneMove = ( vtMove * vtToolDir) * vtToolDir ; //( vtMove * vtToolDir < 0 ? ( - vtMove * vtToolDir) * vtToolDir : ( vtMove * vtToolDir) * vtToolDir) ;
   Vector3d vtZMove = vtMove - vtPlaneMove ;

   double dLenPlanePath = sqrt( vtPlaneMove * vtPlaneMove) ;

   if ( m_nToolType == 1) {

      if ( abs ( ( vtMove * vtToolDir) * dDelta) < EPS_SMALL) {

         for ( unsigned int i = nStartI ; i <= nEndI ; ++ i) 
            for ( unsigned int j = nStartJ ; j <= nEndJ ; ++ j) {

               Point3d ptC( ( i + 0.5) * m_dStep, ( j + 0.5) * m_dStep, 0) ;
               Point3d ptStart( ptLs.x, ptLs.y, 0) ;
               Vector3d vtC = ptC - ptStart ;
               double dProj = vtC * vtToolDir ;
               Vector3d vtP = vtC - dProj * vtToolDir ;
               double dSqD = vtP * vtP ;
               double dZH = ptLs.z ;

               if ( dSqD < m_dRadius * m_dRadius && dProj < 0 && dProj > - m_dHeight) {

                  double dFactor = ( dDelta > 0 ? 1 : - 1) ;
                  dMin = min( dZH, dZH + dDelta + dFactor * sqrt( m_dRadius * m_dRadius - dSqD)) ;
                  dMax = max( dZH, dZH + dDelta + dFactor * sqrt( m_dRadius * m_dRadius - dSqD)) ; 

                  SubtractIntervals( i, j, dMin, dMax) ;
               }
            }
      }
      else {

         Point3d ptSt ;

         if ( vtMove * vtToolDir < 0 && dDelta > 0)
            ptSt = ptLs ;
         else if ( vtMove * vtToolDir < 0 && dDelta < 0) {
            ptSt = ptTLe ; vtPlaneMove = - vtPlaneMove ; dDelta = - dDelta ;
         }
         else if ( vtMove * vtToolDir > 0 && dDelta < 0) {
            ptSt = ptLe ; vtPlaneMove = - vtPlaneMove ; dDelta = - dDelta ;
         }
         else {
            ptSt = ptTLs ;// vtPlaneMove = - vtPlaneMove ;
         }

         double dZH = ptSt.z ;
         Point3d ptStart ( ptSt.x, ptSt.y, 0) ;

         for ( unsigned int i = nStartI ; i <= nEndI ; ++ i)
            for ( unsigned int j = nStartJ ; j <= nEndJ ; ++ j) {

               vtPlaneMove.Normalize() ;

               Point3d ptC( ( i + 0.5) * m_dStep, ( j + 0.5) * m_dStep, 0) ;
               Vector3d vtC = ptC - ptStart ;

               double dProj = vtPlaneMove * vtC ;

               Vector3d vtP = vtC - dProj * vtPlaneMove ;

               double dSqD = vtP * vtP ;

               if ( dSqD < m_dRadius * m_dRadius) {

                  if ( dProj > 0 && dProj < m_dHeight) {
               
                     dMin = dZH - sqrt( m_dRadius * m_dRadius - dSqD) ;

                     if ( dProj < dLenPlanePath)
                        dMax = dZH + sqrt( m_dRadius * m_dRadius - dSqD) + ( dDelta/dLenPlanePath) * dProj ;
                     else
                        dMax = dZH + dDelta + sqrt( m_dRadius * m_dRadius - dSqD) ;  

                     SubtractIntervals( i, j, dMin, dMax) ;
                  }
                  else if ( dProj >= m_dHeight && dProj < m_dHeight + dLenPlanePath) { 
               
                     dMin = dZH - sqrt( m_dRadius * m_dRadius - dSqD) + ( dDelta/dLenPlanePath) * ( dProj - m_dHeight) ;

                     if ( dProj < dLenPlanePath)
                        dMax = dZH + sqrt( m_dRadius * m_dRadius - dSqD) + ( dDelta/dLenPlanePath) * dProj ;
                     else 
                        dMax = dZH + dDelta + sqrt( m_dRadius * m_dRadius - dSqD) ;

                     SubtractIntervals( i, j, dMin, dMax) ;
                  }
               }            
            }
      }
      return true ;
   }
   else if ( m_nToolType == 2) {

      if ( abs ( ( vtMove * vtToolDir) * dDelta) < EPS_SMALL) {

         for ( unsigned int i = nStartI ; i <= nEndI ; ++ i) 
            for ( unsigned int j = nStartJ ; j <= nEndJ ; ++ j) {

               Point3d ptC( ( i + 0.5) * m_dStep, ( j + 0.5) * m_dStep, 0) ;
               Point3d ptStart( ptLs.x, ptLs.y, 0) ;
               Vector3d vtC = ptC - ptStart ;
               double dProj = vtC * vtToolDir ;
               Vector3d vtP = vtC - dProj * vtToolDir ;
               double dSqD = vtP * vtP ;
               double dZH = ptLs.z ;


               if ( dProj < 0 && dProj > - m_dHeight + m_dRadius) {
               
                  if ( dSqD < m_dRadius * m_dRadius) {

                     dMin = min( dZH - sqrt( m_dRadius * m_dRadius - dSqD), dZH + dDelta - sqrt( m_dRadius * m_dRadius - dSqD)) ;
                     dMax = max( dZH + sqrt( m_dRadius * m_dRadius - dSqD), dZH + dDelta + sqrt( m_dRadius * m_dRadius - dSqD)) ; 

                     SubtractIntervals( i, j, dMin, dMax) ;
                  }
               }        // Chiedere per =
               else if ( dProj <= - m_dHeight + m_dRadius  && dProj > - m_dHeight) {     
                     
                  double dCylH = m_dHeight - m_dRadius ;                            
                  double dL = abs( - dProj - dCylH) ;
                  double dR = ( m_dRadius * m_dRadius - dL * dL > 0 ? sqrt( m_dRadius * m_dRadius - dL * dL) : 0) ;

                  if ( dSqD < dR * dR) {

                     double dH = ( dR * dR - dSqD > 0 ? sqrt ( dR * dR - dSqD) : 0) ;

                     dMin = min( dZH - dH, dZH - dH + dDelta) ;
                     dMax = max( dZH + dH, dZH + dH + dDelta) ;

                     SubtractIntervals( i, j, dMin, dMax) ;
                  }
               }   
            }
      }
      else {  // Movimento obliquo

         Point3d ptStart = ptLs ;

         double dA = m_dRadius / sqrt( 1 + ( dDelta * dDelta) / (dLenPlanePath * dLenPlanePath)) ;
         double dSemiAxMin = ( abs( dDelta) * dA) / dLenPlanePath ;

         if ( ( vtToolDir * vtMove) * dDelta < 0) {

            if ( dDelta < 0) {
               vtPlaneMove = - vtPlaneMove ;
               dDelta = - dDelta ;
               ptStart = ptLe ;
            }

            for ( unsigned int i = nStartI ; i <= nEndI ; ++ i) 
               for ( unsigned int j = nStartJ ; j <= nEndJ ; ++ j) {
               
                  Point3d ptS( ptStart.x, ptStart.y, 0) ;
                  Point3d ptC( ( i + 0.5) * m_dStep, ( j + 0.5) * m_dStep, 0) ;
                  Vector3d vtC = ptC - ptS ;
                  double dProj = vtC * vtToolDir ;
                  Vector3d vtP = vtC - ( vtC * vtToolDir) * vtToolDir ;
                  double dSqD = vtP * vtP ;
                  double dCylH = m_dHeight - m_dRadius ;
                  double dZH = ptStart.z ;
                   
                  if ( dSqD < m_dRadius * m_dRadius) {

                     if ( dProj < 0 && dProj > - m_dHeight - dLenPlanePath) {
                  
                        if ( dProj > -dCylH) {
                          // parte cilindrica orizzontale inferiore
                           double dH = ( m_dRadius * m_dRadius - dSqD > 0 ? sqrt( m_dRadius * m_dRadius - dSqD) : 0) ;
                           dMin = dZH - dH ; 
                        }               
                        else if ( dProj > - dCylH - dSemiAxMin * sqrt( 1 - dSqD / ( m_dRadius * m_dRadius))) {
                          // parte sferica inferiore
                           double dL = - dProj - dCylH ;
                           double dR = ( m_dRadius * m_dRadius - dL * dL > 0 ? sqrt( m_dRadius * m_dRadius - dL * dL) : 0) ;
                           double dH = ( dR * dR - dSqD > 0 ? sqrt( dR * dR - dSqD) : 0) ;

                           dMin = dZH - dH ;
                        }
                        else if ( dProj > - dCylH - dLenPlanePath - dSemiAxMin * sqrt( 1 - dSqD / (m_dRadius * m_dRadius))) {
                          // parte cilindrica obliqua inferiore
                           double dl = dSemiAxMin * sqrt( 1 - dSqD / (m_dRadius * m_dRadius)) ; 

                           dMin = dZH - dA * sqrt( 1 - dSqD / (m_dRadius * m_dRadius)) - ( dDelta / dLenPlanePath) * ( dProj + dCylH + dl) ;
                        }
                        else if ( dProj > - m_dHeight - dLenPlanePath) {
                          // seconda parte sferica inferiore 
                           double dL = - dProj - dCylH - dLenPlanePath;
                           double dR = ( m_dRadius * m_dRadius - dL * dL > 0 ? sqrt( m_dRadius * m_dRadius - dL * dL) : 0) ;
                           double dH = ( dR * dR - dSqD > 0 ? sqrt( dR * dR - dSqD) : 0) ;

                           dMin = dZH + dDelta - dH ;
                        }

                        if ( dProj > - dLenPlanePath) {
                          // parte cilindrica obliqua superiore
                           double dH = sqrt( m_dRadius * m_dRadius - dSqD) ;

                           dMax = dZH + dH - ( dDelta / dLenPlanePath) * dProj ;
                        }
                        else if ( dProj > - dLenPlanePath - dCylH) {
                          // parte cilindrica orizzontale superiore
                           double dH = sqrt( m_dRadius * m_dRadius - dSqD) ;

                           dMax = dZH + dDelta + dH ;
                        }
                        else if ( dProj > - dLenPlanePath - m_dHeight) {
                          // parte sferica superiore
                           double dL = - dProj - dCylH - dLenPlanePath;
                           double dR = ( m_dRadius * m_dRadius - dL * dL > 0 ? sqrt( m_dRadius * m_dRadius - dL * dL) : 0) ;
                           double dH = ( dR * dR - dSqD > 0 ? sqrt( dR * dR - dSqD) : 0) ;

                           dMax = dZH + dDelta + dH ;
                        }
                     
                        SubtractIntervals( i, j, dMin, dMax) ;
                     }
                  }
               }
               return true ;
         }
         else {

            if ( dDelta > 0) {
               vtPlaneMove = - vtPlaneMove ;
               dDelta = - dDelta ;
               ptStart = ptLe ;
            }

            for ( unsigned int i = nStartI ; i <= nEndI ; ++ i) 
               for ( unsigned int j = nStartJ ; j <= nEndJ ; ++ j) {
               
                  Point3d ptS( ptStart.x, ptStart.y, 0) ;
                  Point3d ptC( ( i + 0.5) * m_dStep, ( j + 0.5) * m_dStep, 0) ;
                  Vector3d vtC = ptC - ptS ;
                  double dProj = vtC * vtToolDir ;
                  Vector3d vtP = vtC - ( vtC * vtToolDir) * vtToolDir ;
                  double dSqD = vtP * vtP ;
                  double dCylH = m_dHeight - m_dRadius ;
                  double dZH = ptStart.z ;
                   
                  if ( dSqD < m_dRadius * m_dRadius) {

                     if ( dProj < 0 && dProj > - m_dHeight - dLenPlanePath) {

                        if ( dProj > - dLenPlanePath) {
                     
                           double dH = sqrt( m_dRadius * m_dRadius - dSqD) ;

                           dMin = dZH - dH - ( dDelta / dLenPlanePath) * dProj ;
                        }
                        else if ( dProj > - dLenPlanePath - dCylH) {

                           double dH = sqrt( m_dRadius * m_dRadius - dSqD) ;

                           dMin = dZH + dDelta - dH ;
                        }
                        else if ( dProj > - dLenPlanePath - m_dHeight) {
                     
                           double dL = - dProj - dCylH - dLenPlanePath;
                           double dR = ( m_dRadius * m_dRadius - dL * dL > 0 ? sqrt( m_dRadius * m_dRadius - dL * dL) : 0) ;
                           double dH = ( dR * dR - dSqD > 0 ? sqrt( dR * dR - dSqD) : 0) ;

                           dMin = dZH + dDelta - dH ;
                        }
                  
                        if ( dProj > - dCylH) {
                     
                           double dH = ( m_dRadius * m_dRadius - dSqD > 0 ? sqrt( m_dRadius * m_dRadius - dSqD) : 0) ;
                           dMax = dZH + dH ; 
                        }               
                        else if ( dProj > - dCylH - dSemiAxMin * sqrt( 1 - dSqD / ( m_dRadius * m_dRadius))) {
                        
                           double dL = - dProj - dCylH ;
                           double dR = ( m_dRadius * m_dRadius - dL * dL > 0 ? sqrt( m_dRadius * m_dRadius - dL * dL) : 0) ;
                           double dH = ( dR * dR - dSqD > 0 ? sqrt( dR * dR - dSqD) : 0) ;

                           dMax = dZH + dH ;
                        }
                        else if ( dProj > - dCylH - dLenPlanePath - dSemiAxMin * sqrt( 1 - dSqD / (m_dRadius * m_dRadius))) {

                           double dl = dSemiAxMin * sqrt( 1 - dSqD / (m_dRadius * m_dRadius)) ; 

                           dMax = dZH + dA * sqrt( 1 - dSqD / (m_dRadius * m_dRadius)) - ( dDelta / dLenPlanePath) * ( dProj + dCylH + dl) ;
                        }
                        else if ( dProj > - m_dHeight - dLenPlanePath) {
                           
                           double dL = - dProj - dCylH - dLenPlanePath;
                           double dR = ( m_dRadius * m_dRadius - dL * dL > 0 ? sqrt( m_dRadius * m_dRadius - dL * dL) : 0) ;
                           double dH = ( dR * dR - dSqD > 0 ? sqrt( dR * dR - dSqD) : 0) ;

                           dMax = dZH + dDelta + dH ;
                        }
                     
                        SubtractIntervals( i, j, dMin, dMax) ;
                     }
                  }
               }
            return true ;
         }
      }

      return true ;
   }
   else if ( m_nToolType == 3)
      return true ;
   else 
      return true ; 
} 

//----------------------------------------------------------------------------
bool
VolZmap::MillingXYCyl( const Point3d& ptLs, const Point3d& ptLe, const Vector3d& vtToolDir) {


   double dMinZ = min( ptLs.z, ptLe.z) - m_dRadius ;
   double dMaxZ = max( ptLs.z, ptLe.z) + m_dRadius ;

  // Prima verifica sull'interferenza dell'utensile con lo Zmap
   if ( dMaxZ < m_dMinZ + EPS_SMALL || dMinZ > m_dMaxZ - EPS_SMALL)
      return true ;

   Point3d ptI, ptF ;

   if ( ptLs.z < ptLe.z) {
      ptI = ptLs ;
      ptF = ptLe ;
   }
   else {
      ptI = ptLe ;
      ptF = ptLs ;
   }

  // Quote dei punti ptI e ptF
   double dZI = ptI.z ; double dZF = ptF.z ;

   Point3d ptIT = ptI - vtToolDir * m_dHeight ;
   Point3d ptFT = ptF - vtToolDir * m_dHeight ;

  // Bounding box 
   double dMinX = min( min( ptI.x, ptIT.x), min( ptF.x, ptFT.x)) - m_dRadius ;
   double dMaxX = max( max( ptI.x, ptIT.x), max( ptF.x, ptFT.x)) + m_dRadius ;
   double dMinY = min( min( ptI.y, ptIT.y), min( ptF.y, ptFT.y)) - m_dRadius ;
   double dMaxY = max( max( ptI.y, ptIT.y), max( ptF.y, ptFT.y)) + m_dRadius ;

  // Seconda verifica dell'interferenza dell'utensile con lo Zmap
   if ( dMaxX < EPS_SMALL || dMinX > m_nNx * m_dStep - EPS_SMALL)
      return true ;
   if ( dMaxY < EPS_SMALL || dMinY > m_nNy * m_dStep - EPS_SMALL)
      return true ;
   
   Vector3d vtMove = ptF - ptI ; //double dLenPath = vtMove.Len() ;

  // Determino sistema di riferimento del movimento: costruisco un sistrema di vettori ortonormali e destrorsi spiccati dal punto iniziale del movimento
   Vector3d vtV1 = vtToolDir ; 
   Vector3d vtV2 = vtMove ; vtV2.Normalize() ;
   Vector3d vtV3  = vtV1 ^ vtV2 ;

  // Determinazione punti notevoli del volume spazzato dall'utensile
   Point3d ptPlaneSup, ptPlaneInf ;

   if ( vtV3.z > 0) {
      ptPlaneInf = ptI - m_dRadius * vtV3 ;
      ptPlaneSup = ptI + m_dRadius * vtV3 ;
   }
   else {
      ptPlaneInf = ptI + m_dRadius * vtV3 ;
      ptPlaneSup = ptI - m_dRadius * vtV3 ;
   }

  // Prodotti scalari per costruire i piani passanti per i punti notevoli 
   Vector3d vtR0Inf = ptPlaneInf - ORIG ;
   Vector3d vtR0Sup = ptPlaneSup - ORIG ;

   double dPInf = vtR0Inf * vtV3 ; double dPSup = vtR0Sup * vtV3 ;

  // Determinazione delle proiezioni sul piano delle entità geometriche fondamentali
   Vector3d vtPlaneMove( vtMove.x, vtMove.y, 0) ; double dPlanePath = vtPlaneMove.Len() ; vtPlaneMove.Normalize() ;
   Vector3d vtPlaneLim( m_dRadius * vtV3.x, m_dRadius * vtV3.y, 0) ; double dPlaneLim = vtPlaneLim.Len() ;
   Point3d ptIxy( ptI.x, ptI.y, 0) ;
    
  // Determinazione limiti sugli indici
   unsigned int nStartI = ( dMinX < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinX / m_dStep)) ;
   unsigned int nEndI = ( dMaxX > m_nNx * m_dStep - EPS_SMALL ? m_nNx - 1 : static_cast <unsigned int> ( dMaxX / m_dStep)) ;
   unsigned int nStartJ = ( dMinY < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinY / m_dStep)) ;
   unsigned int nEndJ = ( dMaxY > m_nNy * m_dStep - EPS_SMALL ? m_nNy - 1 : static_cast <unsigned int> ( dMaxY / m_dStep)) ;
 
  // Ciclo
   for ( unsigned int i = nStartI ; i <= nEndI ; ++ i) 
      for ( unsigned int j = nStartJ ; j <= nEndJ ; ++ j) {

         double dMin, dMax ;

         double dX = ( i + 0.5) * m_dStep ; double dY = ( j + 0.5) * m_dStep ; 
         double dZInf = ( dPInf - vtV3.x * dX - vtV3.y * dY) / vtV3.z ;  // Quota Z Piano inferiore come funzione di x e y
         double dZSup = ( dPSup - vtV3.x * dX - vtV3.y * dY) / vtV3.z ;  // Quota Z Piano inferiore come funzione di x e y

         Point3d ptC( dX, dY, 0) ;
         Vector3d vtC = ptC - ptIxy ;

         double dProj1 = vtC * vtV1 ;         // vtV1 è vtToolDir che per il momento giace nel piano
         double dProj2 = vtC * vtPlaneMove ;  // vtPlaneMove è stato normalizzato dopo averne calcolato la lunghezza

         if ( dProj1 < 0 && dProj1 > - m_dHeight && dProj2 > - m_dRadius && dProj2 < dPlanePath + m_dRadius) {
            
            
           // Minimi 
            if ( dProj2 > - m_dRadius && dProj2 < dPlaneLim) { 
               
               double dH = sqrt( m_dRadius * m_dRadius - dProj2 * dProj2) ;

               dMin = dZI - dH ;
            }
            else if ( dProj2 >= dPlaneLim && dProj2 < dPlanePath + dPlaneLim) { 
               
               dMin = dZInf ;
            }
            else if ( dProj2 >= dPlanePath + dPlaneLim && dProj2 < dPlanePath + m_dRadius) { 
               
               double dH = sqrt( m_dRadius * m_dRadius - ( dProj2 - dPlanePath) * ( dProj2 - dPlanePath)) ;
               
               dMin = dZF - dH ;
            }
           // Massimi
            if ( dProj2 > - m_dRadius && dProj2 < - dPlaneLim) { 
               
               double dH = sqrt( m_dRadius * m_dRadius - dProj2 * dProj2) ;

               dMax = dZI + dH ; 
            }
            else if ( dProj2 >= - dPlaneLim && dProj2 < dPlanePath - dPlaneLim) { 
               
               dMax = dZSup ;
            }
            else if ( dProj2 >= dPlanePath - dPlaneLim && dProj2 < dPlanePath + m_dRadius) { 
               
               double dH = sqrt( m_dRadius * m_dRadius - ( dProj2 - dPlanePath) * ( dProj2 - dPlanePath)) ;

               dMax = dZF + dH ; 
            }

            SubtractIntervals( i, j, dMin, dMax) ; 
         }
      } // Fine ciclo
      return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::MillingXYBall( const Point3d& ptLs, const Point3d& ptLe, const Vector3d& vtToolDir) {

   double dMinZ = min( ptLs.z, ptLe.z) - m_dRadius ;
   double dMaxZ = max( ptLs.z, ptLe.z) + m_dRadius ;

  // Prima verifica sull'interferenza dell'utensile con lo Zmap
   if ( dMaxZ < m_dMinZ + EPS_SMALL || dMinZ > m_dMaxZ - EPS_SMALL)
      return true ;

   Point3d ptI, ptF ;

   if ( ptLs.z < ptLe.z) {
      ptI = ptLs ;
      ptF = ptLe ;
   }
   else {
      ptI = ptLe ;
      ptF = ptLs ;
   }

  // Quote dei punti ptI e ptF e DeltaZ
   double dZI = ptI.z ; double dZF = ptF.z ; double dDeltaZ = dZF - dZI ;

   Point3d ptIT = ptI - vtToolDir * m_dHeight ;
   Point3d ptFT = ptF - vtToolDir * m_dHeight ;

  // Bounding box 
   double dMinX = min( min( ptI.x, ptIT.x), min( ptF.x, ptFT.x)) - m_dRadius ;
   double dMaxX = max( max( ptI.x, ptIT.x), max( ptF.x, ptFT.x)) + m_dRadius ;
   double dMinY = min( min( ptI.y, ptIT.y), min( ptF.y, ptFT.y)) - m_dRadius ;
   double dMaxY = max( max( ptI.y, ptIT.y), max( ptF.y, ptFT.y)) + m_dRadius ;

  // Seconda verifica dell'interferenza dell'utensile con lo Zmap
   if ( dMaxX < EPS_SMALL || dMinX > m_nNx * m_dStep - EPS_SMALL)
      return true ;
   if ( dMaxY < EPS_SMALL || dMinY > m_nNy * m_dStep - EPS_SMALL)
      return true ;

   Vector3d vtMove = ptF - ptI ; 

  // Determino sistema di riferimento del movimento: costruisco un sistrema di vettori ortonormali e destrorsi spiccati dal punto iniziale del movimento
   Vector3d vtV1 = vtToolDir ; 
   Vector3d vtV2 = vtMove ; vtV2.Normalize() ;
   Vector3d vtV3  = vtV1 ^ vtV2 ;

  // Determinazione punti notevoli del volume spazzato dall'utensile
   Point3d ptPlaneSup, ptPlaneInf ;

   if ( vtV3.z > 0) {
      ptPlaneInf = ptI - m_dRadius * vtV3 ;
      ptPlaneSup = ptI + m_dRadius * vtV3 ;
   }
   else {
      ptPlaneInf = ptI + m_dRadius * vtV3 ;
      ptPlaneSup = ptI - m_dRadius * vtV3 ;
   }

  // Prodotti scalari per costruire i piani passanti per i punti notevoli 
   Vector3d vtR0Inf = ptPlaneInf - ORIG ;
   Vector3d vtR0Sup = ptPlaneSup - ORIG ;

   double dPInf = vtR0Inf * vtV3 ; double dPSup = vtR0Sup * vtV3 ;

  // Determinazione delle proiezioni sul piano delle entità geometriche fondamentali
   Vector3d vtPlaneMove( vtMove.x, vtMove.y, 0) ; double dPlanePath = vtPlaneMove.Len() ; vtPlaneMove.Normalize() ;
   Vector3d vtPlaneLim( m_dRadius * vtV3.x, m_dRadius * vtV3.y, 0) ; double dPlaneLim = vtPlaneLim.Len() ;
   Point3d ptIxy( ptI.x, ptI.y, 0) ;
    
  // Determinazione limiti sugli indici
   unsigned int nStartI = ( dMinX < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinX / m_dStep)) ;
   unsigned int nEndI = ( dMaxX > m_nNx * m_dStep - EPS_SMALL ? m_nNx - 1 : static_cast <unsigned int> ( dMaxX / m_dStep)) ;
   unsigned int nStartJ = ( dMinY < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinY / m_dStep)) ;
   unsigned int nEndJ = ( dMaxY > m_nNy * m_dStep - EPS_SMALL ? m_nNy - 1 : static_cast <unsigned int> ( dMaxY / m_dStep)) ;

    // Ciclo
   for ( unsigned int i = nStartI ; i <= nEndI ; ++ i) 
      for ( unsigned int j = nStartJ ; j <= nEndJ ; ++ j) {

         double dMin, dMax ;
         double dCylH = m_dHeight - m_dRadius ;

         double dX = ( i + 0.5) * m_dStep ; double dY = ( j + 0.5) * m_dStep ; 
         double dZInf = ( dPInf - vtV3.x * dX - vtV3.y * dY) / vtV3.z ;  // Quota Z Piano inferiore come funzione di x e y
         double dZSup = ( dPSup - vtV3.x * dX - vtV3.y * dY) / vtV3.z ;  // Quota Z Piano inferiore come funzione di x e y

         Point3d ptC( dX, dY, 0) ;
         Vector3d vtC = ptC - ptIxy ;

         double dProj1 = vtC * vtV1 ;         // vtV1 è vtToolDir che per il momento giace nel piano
         double dProj2 = vtC * vtPlaneMove ;  // vtPlaneMove è stato normalizzato dopo averne calcolato la lunghezza

        // Parte cilindrica
         if ( dProj1 < 0 && dProj1 > - dCylH && dProj2 > - m_dRadius && dProj2 < dPlanePath + m_dRadius) {

             // Minimi 
            if ( dProj2 > - m_dRadius && dProj2 < dPlaneLim) { 
               
               double dH = sqrt( m_dRadius * m_dRadius - dProj2 * dProj2) ;

               dMin = dZI - dH ;
            }
            else if ( dProj2 >= dPlaneLim && dProj2 < dPlanePath + dPlaneLim) { 
               
               dMin = dZInf ;
            }
            else if ( dProj2 >= dPlanePath + dPlaneLim && dProj2 < dPlanePath + m_dRadius) { 
               
               double dH = sqrt( m_dRadius * m_dRadius - ( dProj2 - dPlanePath) * ( dProj2 - dPlanePath)) ;
               
               dMin = dZF - dH ;
            }
           // Massimi
            if ( dProj2 > - m_dRadius && dProj2 < - dPlaneLim) { 
               
               double dH = sqrt( m_dRadius * m_dRadius - dProj2 * dProj2) ;

               dMax = dZI + dH ; 
            }
            else if ( dProj2 >= - dPlaneLim && dProj2 < dPlanePath - dPlaneLim) { 
               
               dMax = dZSup ;
            }
            else if ( dProj2 >= dPlanePath - dPlaneLim && dProj2 < dPlanePath + m_dRadius) { 
               
               double dH = sqrt( m_dRadius * m_dRadius - ( dProj2 - dPlanePath) * ( dProj2 - dPlanePath)) ;

               dMax = dZF + dH ; 
            }

            SubtractIntervals( i, j, dMin, dMax) ; 
         }
        // Parte sferica
         else if ( ( dProj2 > - m_dRadius && dProj2 < 0 && dProj2 * dProj2 + (- dProj1 - dCylH) * (- dProj1 - dCylH) < m_dRadius * m_dRadius) ||
                   ( dProj2 >= 0 && dProj2 < dPlanePath && dProj1 > - m_dHeight && dProj1 <= - dCylH) ||
                   ( dProj2 >= dPlanePath && dProj2 < dPlanePath + m_dRadius && ( dProj2 - dPlanePath) * ( dProj2 - dPlanePath) + (- dProj1 - dCylH) * (- dProj1 - dCylH) < m_dRadius * m_dRadius)) {


            double dSemiMin = dPlaneLim ; double dSemiMax = m_dRadius ;   // Semi-assi dell'ellisse
            double dl = - dProj1 - dCylH ; 

           // Massimi
            if ( dProj2 < - dSemiMin * sqrt( 1 - ( dl * dl) / ( dSemiMax * dSemiMax))) {

               double dr = sqrt( dProj2 * dProj2 + dl * dl) ;
               double dH = sqrt( m_dRadius * m_dRadius - dr * dr) ;

               dMax = dZI + dH ;
            }
            else if ( dProj2 < dPlanePath - dSemiMin * sqrt( 1 - ( dl * dl) / ( dSemiMax * dSemiMax))) {
            
               double dCos = abs( vtV3 * Z_AX) ;
               
               dMax = dZI + sqrt( m_dRadius * m_dRadius - dl * dl) * dCos + ( dDeltaZ / dPlanePath) * ( dProj2 + dSemiMin * sqrt( 1 - ( dl * dl) / ( dSemiMax * dSemiMax))) ;
            }
            else {

               double dL = dProj2 - dPlanePath ;
               double dr = sqrt( dL * dL + dl * dl) ;

               double dH = sqrt( m_dRadius * m_dRadius - dr * dr) ;

               dMax = dZF + dH ;
            }

           // Minimi
            if ( dProj2 < dSemiMin * sqrt( 1 - ( dl * dl) / ( dSemiMax * dSemiMax))) {

               double dr = sqrt( dProj2 * dProj2 + dl * dl) ;
               double dH = sqrt( m_dRadius * m_dRadius - dr * dr) ;

               dMin = dZI - dH ;
            }
            else if ( dProj2 < dPlanePath + dSemiMin * sqrt( 1 - ( dl * dl) / ( dSemiMax * dSemiMax))) {

               double dCos = abs( vtV3 * Z_AX) ;
               
               dMin = dZI - sqrt( m_dRadius * m_dRadius - dl * dl) * dCos + ( dDeltaZ / dPlanePath) * ( dProj2 - dSemiMin * sqrt( 1 - ( dl * dl) / ( dSemiMax * dSemiMax))) ;
            }
            else {

               double dr = sqrt( ( dProj2 - dPlanePath) * ( dProj2 - dPlanePath) + dl * dl) ;
               double dH = sqrt( m_dRadius * m_dRadius - dr * dr) ;

               dMin = dZF - dH ;
            }
         
            SubtractIntervals( i, j, dMin, dMax) ;
         }
      } // Fine ciclo
      return true ;
}
/*
//  DA TENERE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
//----------------------------------------------------------------------------
bool
VolZmap::MillingXYPlusCyl( const Point3d& ptLs, const Point3d& ptLe, const Vector3d& vtToolDir) {

   double dMinZ = min( ptLs.z, ptLe.z) - m_dRadius ;
   double dMaxZ = max( ptLs.z, ptLe.z) + m_dRadius ;

  // Prima verifica sull'interferenza dell'utensile con lo Zmap
   if ( dMaxZ < m_dMinZ + EPS_SMALL || dMinZ > m_dMaxZ - EPS_SMALL)
      return true ;

   Point3d ptI, ptF ;

   if ( ptLs.z <= ptLe.z) {
      ptI = ptLs ;
      ptF = ptLe ;
   }
   else {
      ptI = ptLe ;
      ptF = ptLs ;
   }

  // Quote dei punti ptI e ptF e DeltaZ
   double dZI = ptI.z ; double dZF = ptF.z ; double dDeltaZ = dZF - dZI ;

   Point3d ptIT = ptI - vtToolDir * m_dHeight ;
   Point3d ptFT = ptF - vtToolDir * m_dHeight ;

  // Bounding box 
   double dMinX = min( min( ptI.x, ptIT.x), min( ptF.x, ptFT.x)) - m_dRadius ;
   double dMaxX = max( max( ptI.x, ptIT.x), max( ptF.x, ptFT.x)) + m_dRadius ;
   double dMinY = min( min( ptI.y, ptIT.y), min( ptF.y, ptFT.y)) - m_dRadius ;
   double dMaxY = max( max( ptI.y, ptIT.y), max( ptF.y, ptFT.y)) + m_dRadius ;

  // Seconda verifica dell'interferenza dell'utensile con lo Zmap
   if ( dMaxX < EPS_SMALL || dMinX > m_nNx * m_dStep - EPS_SMALL)
      return true ;
   if ( dMaxY < EPS_SMALL || dMinY > m_nNy * m_dStep - EPS_SMALL)
      return true ;
  
  // Vettori di riferimento nello spazio e controllo per simmetria
   Vector3d vtMove = ptF - ptI ;
   Vector3d vtTool = vtToolDir ;

   if ( vtToolDir * vtMove > 0) {
   
      Point3d ptTemp = ptI ;
      ptI = ptIT ; ptIT = ptTemp ;
      ptTemp = ptF ;
      ptF = ptFT ;
      ptFT = ptTemp ;
      vtTool = - vtTool ;
   }
          
   Vector3d vtMoveOrt = vtMove - ( vtMove * vtTool) * vtTool ;     double dLen2 = vtMoveOrt.Len() ;
   Vector3d vtMoveLong = ( vtMove * vtTool) * vtTool ;             double dLen1 = vtMoveLong.Len() ;
 
  // Vettori di riferimento nel piano
   Vector3d vtPlaneMove( vtMove.x, vtMove.y, 0) ;                  double dPLen = vtPlaneMove.LenXY() ;
   Vector3d vtPlaneMoveLong( vtMoveLong.x, vtMoveLong.y, 0) ;      double dPLen1 = vtPlaneMoveLong.LenXY() ; 
   Vector3d vtPlaneMoveOrt( vtMoveOrt.x, vtMoveOrt.y, 0) ;         double dPLen2 = vtPlaneMoveOrt.LenXY() ;   vtPlaneMoveOrt.Normalize() ;
   
  // Punti iniziale e finale proiettati sul piano
   Point3d ptIxy( ptI.x, ptI.y, 0) ;  Point3d ptFxy( ptF.x, ptF.y, 0) ; 
   
  // Determino i sistemi di riferimento del movimento: costruisco un sistrema di vettori ortonormali e destrorsi spiccati dal punto iniziale del movimento
   Vector3d vtV1 = vtTool ; 
   Vector3d vtV2 = vtMoveOrt ;  vtV2.Normalize() ;
   Vector3d vtV3  = vtV1 ^ vtV2 ;

  // Determinazione punti notevoli del volume spazzato dall'utensile
   Point3d ptPlaneSup, ptPlaneInf ;

   if ( vtV3.z > 0) {
      ptPlaneInf = ptI - m_dRadius * vtV3 ;
      ptPlaneSup = ptI + m_dRadius * vtV3 ;
   }
   else {
      ptPlaneInf = ptI + m_dRadius * vtV3 ;
      ptPlaneSup = ptI - m_dRadius * vtV3 ;
   }

   Point3d ptPlaneSupxy( ptPlaneSup.x, ptPlaneSup.y, 0) ;
   Point3d ptPlaneInfxy( ptPlaneInf.x, ptPlaneInf.y, 0) ;

   double dr = sqrt( ( ptPlaneSupxy - ptIxy) * ( ptPlaneSupxy - ptIxy)) ; 
  
  // Determinazione degli analoghi punti sulla punta dell'utensile e delle loro proiezioni sul piano XY 
   Point3d ptPlTInf = ptPlaneInf - m_dHeight * vtTool ;  Point3d ptPlTInfxy( ptPlTInf.x, ptPlTInf.y, 0) ;
   Point3d ptPlTSup = ptPlaneSup - m_dHeight * vtTool ;  Point3d ptPlTSupxy( ptPlTSup.x, ptPlTSup.y, 0) ;

  // Prodotti scalari per costruire i piani passanti per i punti notevoli  
   Vector3d vtR0Inf = ptPlaneInf - ORIG ;
   Vector3d vtR0Sup = ptPlaneSup - ORIG ;
   Vector3d vtR0 = ptI - ORIG ;
   double dPInf = vtR0Inf * vtV3 ; double dPSup = vtR0Sup * vtV3 ; double dP = vtR0 * vtV3 ;
  
  // Determinazione limiti sugli indici
   unsigned int nStartI = ( dMinX < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinX / m_dStep)) ;
   unsigned int nEndI = ( dMaxX > m_nNx * m_dStep - EPS_SMALL ? m_nNx - 1 : static_cast <unsigned int> ( dMaxX / m_dStep)) ;
   unsigned int nStartJ = ( dMinY < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinY / m_dStep)) ;
   unsigned int nEndJ = ( dMaxY > m_nNy * m_dStep - EPS_SMALL ? m_nNy - 1 : static_cast <unsigned int> ( dMaxY / m_dStep)) ;

  // Ciclo
   for ( unsigned int i = nStartI ; i <= nEndI ; ++ i) 
      for ( unsigned int j = nStartJ ; j <= nEndJ ; ++ j) {

         double dMin, dMax ;

         double dX = ( i + 0.5) * m_dStep ; double dY = ( j + 0.5) * m_dStep ; 
         double dZPInf = ( dPInf - vtV3.x * dX - vtV3.y * dY) / vtV3.z ;  // Quota Z Piano inferiore come funzione di x e y
         double dZPSup = ( dPSup - vtV3.x * dX - vtV3.y * dY) / vtV3.z ;  // Quota Z Piano superiore come funzione di x e y
       
        // Punto e vettori del ciclo
         Point3d ptC( dX, dY, 0) ;  Vector3d vtCi = ptC - ptIxy ;    
        // Proiezione fondamentale
         double dPro1 = vtCi * vtV1 ;  double dPro2 = vtCi * vtPlaneMoveOrt ;


        // Se il punto cade nella proiezione sul piano XY del volume spazzato si taglia 
         if ( ((dPro2 > - m_dRadius && dPro2 < m_dRadius) && (((dPro2 < dPLen2 - m_dRadius && (dPro1 > - m_dHeight - (dPLen1 / dPLen2) * (dPro2 + m_dRadius) && dPro1 < 0))) || (dPro2 >= dPLen2 - m_dRadius && (dPro1 > - m_dHeight - dPLen1 && dPro1 < 0))))
             || ((dPro2 >= m_dRadius && dPro2 < dPLen2 + m_dRadius) && ((dPro2 < dPLen2 - m_dRadius && (dPro1 > - (dPLen1 / dPLen2) * (dPro2 + m_dRadius) - m_dHeight && dPro1 < - (dPLen1 / dPLen2) * (dPro2 - m_dRadius))) || (dPro2 >= dPLen2 - m_dRadius && (dPro1 > - m_dHeight - dPLen1 && dPro1 < - (dPLen1 / dPLen2) * (dPro2 - m_dRadius)))))) {
         
           // Massimi //////////////////////////////////////////////////////////////////////////
           // Prima zona cilindrica superiore 
            if ( ( dPro2 > - m_dRadius && dPro2 < - dr) && ( dPro1 > - m_dHeight && dPro1 < 0)) {
            
               double dH = sqrt( m_dRadius * m_dRadius - dPro2 * dPro2) ;

               dMax = dZI + dH ;
            } 

           // Vettore per seconda zona cilindrica superiore
            Vector3d vtCf = ptC - ptFxy ;
           // Proiezione per seconda zona cilindrica sueriore
            double dPr1 = vtCf * vtV1 ;   double dPr2 = vtCf * vtPlaneMoveOrt ;
         
           // Seconda zona cilindrica superiore 
            if ( ( dPr2 >= - dr && dPr2 < m_dRadius) && ( dPr1 > - m_dHeight && dPr1 <= 0)) {
            
               double dH = sqrt( m_dRadius * m_dRadius - dPr2 * dPr2) ;

               dMax = dZF + dH ;
            }

           // Vettore per Piano superiore e zona di fondo superiore
            Vector3d vtCS = ptC - ptPlaneSupxy ;
           // Proiezioni 
            double dPrS1 = vtCS * vtV1 ; double dPrS2 = vtCS * vtPlaneMoveOrt ;

           // Piano superiore 
            if ( dPrS2 >= 0 && dPrS2 < dPLen2 && dPrS1 > - m_dHeight - ( dPLen1/dPLen2) * dPrS2 && dPrS1 <= - ( dPLen1/dPLen2) * dPrS2) 
            
               dMax = dZPSup ;
            
           // Vettore per zona di punta superiore
            Vector3d vtCTS = ptC - ptPlTSupxy ;
           // Proiezioni 
            double dPrTS1 = vtCTS * vtV1 ;  double dPrTS2 = vtCTS * vtPlaneMoveOrt ;

           // Zona di punta superiore
            if ( dPrTS1 <= 0 && dPrTS1 > - dPLen1)
               if ( dPrTS2 <= - ( dPLen2 / dPLen1) * dPrTS1  && dPrTS2 > - ( m_dRadius - dr) - ( dPLen2 / dPLen1) * dPrTS1) {

                 double dDist = - ( dPLen2 / dPLen1) * dPrTS1 - dPrTS2 ;
                 double dL = dDist + dr ;
                 double dH = ( m_dRadius * m_dRadius - dL * dL > 0 ? sqrt( m_dRadius * m_dRadius - dL * dL) : 0) ;
                 double dl = sqrt( dPrTS1 * dPrTS1 + ( ( dPLen2 / dPLen1) * dPrTS1) * ( ( dPLen2 / dPLen1) * dPrTS1)) ;

                 dMax = dZI + ( dDeltaZ / dPLen) * dl + dH ;
               }

           // Zona di fondo superiore
            if ( dPrS1 < 0 && dPrS1 > - dPLen1) 
               if ( dPrS2 > - ( dPLen2 / dPLen1) * dPrS1 && dPrS2 < dr - ( dPLen2 / dPLen1) * dPrS1) {
               
                 double dL = dr - ( dPLen2 / dPLen1) * dPrS1 - dPrS2 ;
                 double dH = ( m_dRadius * m_dRadius - dL * dL > 0 ? sqrt( m_dRadius * m_dRadius - dL * dL) : 0) ;
                 double dl = sqrt( dPrS1 * dPrS1 + ( ( dPLen2 / dPLen1) * dPrS1) * ( ( dPLen2 / dPLen1) * dPrS1)) ;

                 dMax = dZI + ( dDeltaZ / dPLen) * dl + dH ;  
               }
               else if ( dPrS2 >= dr - ( dPLen2 / dPLen1) * dPrS1 && dPrS2 < dr + m_dRadius - ( dPLen2 / dPLen1) * dPrS1) {
               
                  double dL = dPrS2 - ( dr - ( dPLen2 / dPLen1) * dPrS1) ;
                  double dH = ( m_dRadius * m_dRadius - dL * dL > 0 ? sqrt( m_dRadius * m_dRadius - dL * dL) : 0) ;
                  double dl = sqrt( dPrS1 * dPrS1 + ( ( dPLen2 / dPLen1) * dPrS1) * ( ( dPLen2 / dPLen1) * dPrS1)) ;

                  dMax = dZI + ( dDeltaZ / dPLen) * dl + dH ;      
               }   

           // Minimi //////////////////////////////////////////////////////////////////////////
           // Prima zona cilindrica inferiore
            if ( ( dPro2 > - m_dRadius && dPro2 < dr) && ( dPro1 > - m_dHeight && dPro1 < 0)) {
            
               double dH = sqrt( m_dRadius * m_dRadius - dPro2 * dPro2) ;

               dMin = dZI - dH ;
            }

           // Seconda zona cilindrica inferiore
            if ( ( dPr2 >= dr && dPr2 < m_dRadius) && ( dPr1 > - m_dHeight && dPr1 <= 0)) {
            
               double dH = sqrt( m_dRadius * m_dRadius - dPr2 * dPr2) ;

               dMin = dZF - dH ;
            }

           // Vettore per piano inferiore e zona di fondo inferiore
            Vector3d vtCI = ptC - ptPlaneInfxy ;
           // Proiezioni
            double dPrI1 = vtCI * vtV1 ;  double dPrI2 = vtCI * vtPlaneMoveOrt ;

           // Piano inferiore
            if ( dPrI2 >= 0 && dPrI2 < dPLen2 && dPrI1 > - m_dHeight - ( dPLen1/dPLen2) * dPrI2 && dPrI1 <= - ( dPLen1/dPLen2) * dPrI2) 
            
               dMin = dZPInf ;
            
           // Zona di fondo inferiore
            if ( dPrI1 <= 0 && dPrI1 > - dPLen1)
               if ( dPrI2 > - ( dPLen2 / dPLen1) * dPrI1 && dPrI2 < ( m_dRadius - dr) - ( dPLen2 / dPLen1) * dPrI1) {
               
                 double dDist = dPrI2 + ( dPLen2 / dPLen1) * dPrI1 ;
                 double dL = dDist + dr ;
                 double dH = ( m_dRadius * m_dRadius - dL * dL > 0 ? sqrt( m_dRadius * m_dRadius - dL * dL) : 0) ;
                 double dl = sqrt( dPrI1 * dPrI1 + ( ( dPLen2 / dPLen1) * dPrI1) * ( ( dPLen2 / dPLen1) * dPrI1)) ;

                 dMin = dZI + ( dDeltaZ / dPLen) * dl - dH ;
               }

           // Vettore per zona di punta inferiore
            Vector3d vtCTI = ptC - ptPlTInfxy ;
           // Proiezioni
            double dPrTI1 = vtCTI * vtV1 ;  double dPrTI2 = vtCTI * vtPlaneMoveOrt ;

           // zona di punta inferiore
            if ( dPrTI1 <= 0 && dPrTI1 > - dPLen1) 
               if ( dPrTI2 > - m_dRadius - dr - ( dPLen2 / dPLen1) * dPrTI1 && dPrTI2 <= - dr - ( dPLen2 / dPLen1) * dPrTI1) {
               
                  double dL = - dr - ( dPLen2 / dPLen1) * dPrTI1 - dPrTI2 ;
                  double dH = ( m_dRadius * m_dRadius - dL * dL > 0 ? sqrt( m_dRadius * m_dRadius - dL * dL) : 0) ;
                  double dl = sqrt( dPrTI1 * dPrTI1 + ( ( dPLen2 / dPLen1) * dPrTI1) * ( ( dPLen2 / dPLen1) * dPrTI1)) ;

                  dMin = dZI + ( dDeltaZ / dPLen) * dl - dH ;  
               }
               else if ( dPrTI2 > - dr - ( dPLen2 / dPLen1) * dPrTI1 && dPrTI2 < - ( dPLen2 / dPLen1) * dPrTI1) {

                  double dL = dPrTI2 - ( - dr - ( dPLen2 / dPLen1) * dPrTI1 ) ;
                  double dH = ( m_dRadius * m_dRadius - dL * dL > 0 ? sqrt( m_dRadius * m_dRadius - dL * dL) : 0) ;
                  double dl = sqrt( dPrTI1 * dPrTI1 + ( ( dPLen2 / dPLen1) * dPrTI1) * ( ( dPLen2 / dPLen1) * dPrTI1)) ;

                  dMin = dZI + ( dDeltaZ / dPLen) * dl - dH ;    
               } 

     
            SubtractIntervals( i, j, dMin, dMax) ;
         }
      } // Fine ciclo
      return true ;
}*/

//----------------------------------------------------------------------------
bool
VolZmap::MillingXYPlusCyl( const Point3d& ptLs, const Point3d& ptLe, const Vector3d& vtToolDir) {

   double dMinZ = min( ptLs.z, ptLe.z) - m_dRadius ;
   double dMaxZ = max( ptLs.z, ptLe.z) + m_dRadius ;

  // Prima verifica sull'interferenza dell'utensile con lo Zmap
   if ( dMaxZ < m_dMinZ + EPS_SMALL || dMinZ > m_dMaxZ - EPS_SMALL)
      return true ;

   Point3d ptI, ptF ;

   if ( ptLs.z <= ptLe.z) {
      ptI = ptLs ;
      ptF = ptLe ;
   }
   else {
      ptI = ptLe ;
      ptF = ptLs ;
   }

  // Quote dei punti ptI e ptF e DeltaZ
   double dZI = ptI.z ; double dZF = ptF.z ; double dDeltaZ = dZF - dZI ;

   Point3d ptIT = ptI - vtToolDir * m_dHeight ;
   Point3d ptFT = ptF - vtToolDir * m_dHeight ;

  // Bounding box 
   double dMinX = min( min( ptI.x, ptIT.x), min( ptF.x, ptFT.x)) - m_dRadius ;
   double dMaxX = max( max( ptI.x, ptIT.x), max( ptF.x, ptFT.x)) + m_dRadius ;
   double dMinY = min( min( ptI.y, ptIT.y), min( ptF.y, ptFT.y)) - m_dRadius ;
   double dMaxY = max( max( ptI.y, ptIT.y), max( ptF.y, ptFT.y)) + m_dRadius ;

  // Seconda verifica dell'interferenza dell'utensile con lo Zmap
   if ( dMaxX < EPS_SMALL || dMinX > m_nNx * m_dStep - EPS_SMALL)
      return true ;
   if ( dMaxY < EPS_SMALL || dMinY > m_nNy * m_dStep - EPS_SMALL)
      return true ;
  
  // Vettori di riferimento nello spazio e controllo per simmetria
   Vector3d vtMove = ptF - ptI ;
   Vector3d vtTool = vtToolDir ;

   if ( vtToolDir * vtMove > 0) {
   
      Point3d ptTemp = ptI ;
      ptI = ptIT ; ptIT = ptTemp ;
      ptTemp = ptF ;
      ptF = ptFT ;
      ptFT = ptTemp ;
      vtTool = - vtTool ;
   }
   
   Vector3d vtMoveLong = ( vtMove * vtTool) * vtTool ;        double dLen1 = vtMoveLong.Len() ;
   Vector3d vtMoveOrt = vtMove - vtMoveLong ;                 double dLen2 = vtMoveOrt.Len() ;
 
  // Vettori di riferimento nel piano
   Vector3d vtPlaneMove( vtMove.x, vtMove.y, 0) ;                  double dPLen = vtPlaneMove.LenXY() ;
   Vector3d vtPlaneMoveLong( vtMoveLong.x, vtMoveLong.y, 0) ;      double dPLen1 = vtPlaneMoveLong.LenXY() ; 
   Vector3d vtPlaneMoveOrt( vtMoveOrt.x, vtMoveOrt.y, 0) ;         double dPLen2 = vtPlaneMoveOrt.LenXY() ;   vtPlaneMoveOrt.Normalize() ;
   
  // Punti iniziale e finale proiettati sul piano
   Point3d ptIxy( ptI.x, ptI.y, 0) ;  Point3d ptFxy( ptF.x, ptF.y, 0) ; 
   
  // Determino i sistemi di riferimento del movimento: costruisco un sistrema di vettori ortonormali e destrorsi spiccati dal punto iniziale del movimento
   Vector3d vtV1 = vtTool ; 
   Vector3d vtV2 = vtMoveOrt ;  vtV2.Normalize() ;
   Vector3d vtV3  = vtV1 ^ vtV2 ;

  // Determinazione punti notevoli del volume spazzato dall'utensile
   Point3d ptPlaneSup, ptPlaneInf ;

   if ( vtV3.z > 0) {
      ptPlaneInf = ptI - m_dRadius * vtV3 ;
      ptPlaneSup = ptI + m_dRadius * vtV3 ;
   }
   else {
      ptPlaneInf = ptI + m_dRadius * vtV3 ;
      ptPlaneSup = ptI - m_dRadius * vtV3 ;
   }

   Point3d ptPlaneSupxy( ptPlaneSup.x, ptPlaneSup.y, 0) ;
   Point3d ptPlaneInfxy( ptPlaneInf.x, ptPlaneInf.y, 0) ;

   double dr = sqrt( ( ptPlaneSupxy - ptIxy) * ( ptPlaneSupxy - ptIxy)) ; 
  
  // Determinazione degli analoghi punti sulla punta dell'utensile e delle loro proiezioni sul piano XY 
   Point3d ptPlTInf = ptPlaneInf - m_dHeight * vtTool ;  Point3d ptPlTInfxy( ptPlTInf.x, ptPlTInf.y, 0) ;
   Point3d ptPlTSup = ptPlaneSup - m_dHeight * vtTool ;  Point3d ptPlTSupxy( ptPlTSup.x, ptPlTSup.y, 0) ;

  // Prodotti scalari per costruire i piani passanti per i punti notevoli  
   Vector3d vtR0Inf = ptPlaneInf - ORIG ;
   Vector3d vtR0Sup = ptPlaneSup - ORIG ;
   Vector3d vtR0 = ptI - ORIG ;
   double dPInf = vtR0Inf * vtV3 ; double dPSup = vtR0Sup * vtV3 ; double dP = vtR0 * vtV3 ;
  
  // Determinazione limiti sugli indici
   unsigned int nStartI = ( dMinX < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinX / m_dStep)) ;
   unsigned int nEndI = ( dMaxX > m_nNx * m_dStep - EPS_SMALL ? m_nNx - 1 : static_cast <unsigned int> ( dMaxX / m_dStep)) ;
   unsigned int nStartJ = ( dMinY < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinY / m_dStep)) ;
   unsigned int nEndJ = ( dMaxY > m_nNy * m_dStep - EPS_SMALL ? m_nNy - 1 : static_cast <unsigned int> ( dMaxY / m_dStep)) ;

  // Ciclo
   for ( unsigned int i = nStartI ; i <= nEndI ; ++ i) 
      for ( unsigned int j = nStartJ ; j <= nEndJ ; ++ j) {

         double dMin, dMax ;

         double dX = ( i + 0.5) * m_dStep ; double dY = ( j + 0.5) * m_dStep ; 
         double dZPInf = ( dPInf - vtV3.x * dX - vtV3.y * dY) / vtV3.z ;  // Quota Z Piano inferiore come funzione di x e y
         double dZPSup = ( dPSup - vtV3.x * dX - vtV3.y * dY) / vtV3.z ;  // Quota Z Piano superiore come funzione di x e y
       
        // Punto e vettori del ciclo
         Point3d ptC( dX, dY, 0) ;  Vector3d vtCi = ptC - ptIxy ;    
        // Proiezione fondamentale
         double dPro1 = vtCi * vtV1 ;  double dPro2 = vtCi * vtPlaneMoveOrt ;


        // Se il punto cade nella proiezione sul piano XY del volume spazzato si taglia 
         if ( ((dPro2 > - m_dRadius && dPro2 < m_dRadius) && (((dPro2 < dPLen2 - m_dRadius && (dPro1 > - m_dHeight - (dPLen1 / dPLen2) * (dPro2 + m_dRadius) && dPro1 < 0))) || (dPro2 >= dPLen2 - m_dRadius && (dPro1 > - m_dHeight - dPLen1 && dPro1 < 0))))
             || ((dPro2 >= m_dRadius && dPro2 < dPLen2 + m_dRadius) && ((dPro2 < dPLen2 - m_dRadius && (dPro1 > - (dPLen1 / dPLen2) * (dPro2 + m_dRadius) - m_dHeight && dPro1 < - (dPLen1 / dPLen2) * (dPro2 - m_dRadius))) || (dPro2 >= dPLen2 - m_dRadius && (dPro1 > - m_dHeight - dPLen1 && dPro1 < - (dPLen1 / dPLen2) * (dPro2 - m_dRadius)))))) {
         
           // Massimi //////////////////////////////////////////////////////////////////////////
           // Prima zona cilindrica superiore 
            if ( ( dPro2 > - m_dRadius && dPro2 < - dr) && ( dPro1 > - m_dHeight && dPro1 < 0)) {
            
               double dH = sqrt( m_dRadius * m_dRadius - dPro2 * dPro2) ;

               dMax = dZI + dH ;
            } 

           // Vettore per seconda zona cilindrica superiore
            Vector3d vtCf = ptC - ptFxy ;
           // Proiezione per seconda zona cilindrica sueriore
            double dPr1 = vtCf * vtV1 ;   double dPr2 = vtCf * vtPlaneMoveOrt ;
         
           // Seconda zona cilindrica superiore 
            if ( ( dPr2 >= - dr && dPr2 < m_dRadius) && ( dPr1 > - m_dHeight && dPr1 <= 0)) {
            
               double dH = sqrt( m_dRadius * m_dRadius - dPr2 * dPr2) ;

               dMax = dZF + dH ;
            }

           // Vettore per Piano superiore e zona di fondo superiore
            Vector3d vtCS = ptC - ptPlaneSupxy ;
           // Proiezioni 
            double dPrS1 = vtCS * vtV1 ; double dPrS2 = vtCS * vtPlaneMoveOrt ;

           // Piano superiore 
            if ( dPrS2 >= 0 && dPrS2 < dPLen2 && dPrS1 > - m_dHeight - ( dPLen1/dPLen2) * dPrS2 && dPrS1 <= - ( dPLen1/dPLen2) * dPrS2) 
            
               dMax = dZPSup ;
            
           // Vettore per zona di punta superiore
            Vector3d vtCTS = ptC - ptPlTSupxy ;
           // Proiezioni 
            double dPrTS1 = vtCTS * vtV1 ;  double dPrTS2 = vtCTS * vtPlaneMoveOrt ;

           // Zona di punta superiore
            if ( dPrTS1 <= 0 && dPrTS1 > - dPLen1)
               if ( dPrTS2 <= - ( dPLen2 / dPLen1) * dPrTS1  && dPrTS2 > - ( m_dRadius - dr) - ( dPLen2 / dPLen1) * dPrTS1) {

                 double dDist = - ( dPLen2 / dPLen1) * dPrTS1 - dPrTS2 ;
                 double dL = dDist + dr ;
                 double dH = ( m_dRadius * m_dRadius - dL * dL > 0 ? sqrt( m_dRadius * m_dRadius - dL * dL) : 0) ;
                 double dl = sqrt( dPrTS1 * dPrTS1 + ( ( dPLen2 / dPLen1) * dPrTS1) * ( ( dPLen2 / dPLen1) * dPrTS1)) ;

                 dMax = dZI + ( dDeltaZ / dPLen) * dl + dH ;
               }

           // Zona di fondo superiore
            if ( dPrS1 < 0 && dPrS1 > - dPLen1) 
               if ( dPrS2 > - ( dPLen2 / dPLen1) * dPrS1 && dPrS2 < m_dRadius + dr - ( dPLen2 / dPLen1) * dPrS1) {
               
                 double dL = abs( dr - ( dPLen2 / dPLen1) * dPrS1 - dPrS2) ;
                 double dH = ( m_dRadius * m_dRadius - dL * dL > 0 ? sqrt( m_dRadius * m_dRadius - dL * dL) : 0) ;
                 double dl = sqrt( dPrS1 * dPrS1 + ( ( dPLen2 / dPLen1) * dPrS1) * ( ( dPLen2 / dPLen1) * dPrS1)) ;

                 dMax = dZI + ( dDeltaZ / dPLen) * dl + dH ;      
               }   

           // Minimi //////////////////////////////////////////////////////////////////////////
           // Prima zona cilindrica inferiore
            if ( ( dPro2 > - m_dRadius && dPro2 < dr) && ( dPro1 > - m_dHeight && dPro1 < 0)) {
            
               double dH = sqrt( m_dRadius * m_dRadius - dPro2 * dPro2) ;

               dMin = dZI - dH ;
            }

           // Seconda zona cilindrica inferiore
            if ( ( dPr2 >= dr && dPr2 < m_dRadius) && ( dPr1 > - m_dHeight && dPr1 <= 0)) {
            
               double dH = sqrt( m_dRadius * m_dRadius - dPr2 * dPr2) ;

               dMin = dZF - dH ;
            }

           // Vettore per piano inferiore e zona di fondo inferiore
            Vector3d vtCI = ptC - ptPlaneInfxy ;
           // Proiezioni
            double dPrI1 = vtCI * vtV1 ;  double dPrI2 = vtCI * vtPlaneMoveOrt ;

           // Piano inferiore
            if ( dPrI2 >= 0 && dPrI2 < dPLen2 && dPrI1 > - m_dHeight - ( dPLen1/dPLen2) * dPrI2 && dPrI1 <= - ( dPLen1/dPLen2) * dPrI2) 
            
               dMin = dZPInf ;
            
           // Zona di fondo inferiore
            if ( dPrI1 <= 0 && dPrI1 > - dPLen1)
               if ( dPrI2 > - ( dPLen2 / dPLen1) * dPrI1 && dPrI2 < ( m_dRadius - dr) - ( dPLen2 / dPLen1) * dPrI1) {
               
                 double dDist = dPrI2 + ( dPLen2 / dPLen1) * dPrI1 ;
                 double dL = dDist + dr ;
                 double dH = ( m_dRadius * m_dRadius - dL * dL > 0 ? sqrt( m_dRadius * m_dRadius - dL * dL) : 0) ;
                 double dl = sqrt( dPrI1 * dPrI1 + ( ( dPLen2 / dPLen1) * dPrI1) * ( ( dPLen2 / dPLen1) * dPrI1)) ;

                 dMin = dZI + ( dDeltaZ / dPLen) * dl - dH ;
               }

           // Vettore per zona di punta inferiore
            Vector3d vtCTI = ptC - ptPlTInfxy ;
           // Proiezioni
            double dPrTI1 = vtCTI * vtV1 ;  double dPrTI2 = vtCTI * vtPlaneMoveOrt ;

           // zona di punta inferiore
            if ( dPrTI1 <= 0 && dPrTI1 > - dPLen1) 
               if ( dPrTI2 > - m_dRadius - dr - ( dPLen2 / dPLen1) * dPrTI1 && dPrTI2 < - ( dPLen2 / dPLen1) * dPrTI1) {
               
                  double dL = abs( - dr - ( dPLen2 / dPLen1) * dPrTI1 - dPrTI2) ;
                  double dH = ( m_dRadius * m_dRadius - dL * dL > 0 ? sqrt( m_dRadius * m_dRadius - dL * dL) : 0) ;
                  double dl = sqrt( dPrTI1 * dPrTI1 + ( ( dPLen2 / dPLen1) * dPrTI1) * ( ( dPLen2 / dPLen1) * dPrTI1)) ;

                  dMin = dZI + ( dDeltaZ / dPLen) * dl - dH ;    
               } 

            SubtractIntervals( i, j, dMin, dMax) ;
         }
      } // Fine ciclo
      return true ;
}
 
/*
//----------------------------------------------------------------------------
bool
VolZmap::MillingXYPlusBall( const Point3d& ptLs, const Point3d& ptLe, const Vector3d& vtToolDir) {

  // Parte cilindrica

   m_dHeight = m_dHeight - m_dRadius ;

   MillingXYPlusCyl( ptLs, ptLe, vtToolDir) ;

   m_dHeight = m_dHeight + m_dRadius ;
  ////////////////////////////////////////
   Point3d ptI, ptF ;

   if ( ptLs.z < ptLe.z) {
   
      ptI = ptLs ;
      ptF = ptLe ;
   }
   else {
   
      ptI = ptLe ;
      ptF = ptLs ;
   }

   double dZI = ptI.z ;  double dZF = ptF.z ;  double dDeltaZ = dZF - dZI ;

   double dCylH = m_dHeight - m_dRadius ;

   Point3d ptCI = ptI - dCylH * vtToolDir ;
   Point3d ptCF = ptF - dCylH * vtToolDir ;

   Point3d ptCIxy( ptCI.x, ptCI.y, 0) ;

  // Bounding box 
   double dMinX = min( ptCI.x, ptCF.x) - m_dRadius ;
   double dMaxX = max( ptCI.x, ptCF.x) + m_dRadius ;
   double dMinY = min( ptCI.y, ptCF.y) - m_dRadius ;
   double dMaxY = max( ptCI.y, ptCF.y) + m_dRadius ;

   Vector3d vtMove = ptF - ptI ;

  // Sistema di riferimento
   Vector3d vtV2 = vtMove ;    vtV2.Normalize() ; 
   Vector3d vtV1 ;

   if ( abs( vtV2.y) < EPS_SMALL) 
   
      vtV1 = Y_AX ;
   else 
      vtV1 = X_AX - ( vtV2.x / vtV2.y) * Y_AX ;

   if ( vtV1 * vtToolDir < 0)

      vtV1 = - vtV1 ;

   vtV1.Normalize() ;

   Vector3d vtV3 = vtV1 ^ vtV2 ;

  // Determino il semi-asse minore 
   Vector3d vtSemiMin( vtV3.x, vtV3.y, 0) ;
   double dSemiAxMin = m_dRadius * vtSemiMin.LenXY() ;

   Vector3d vtMovexy( vtMove.x, vtMove.y, 0) ;   double dPlaneLen = vtMovexy.LenXY() ;
   Vector3d vtU1 = vtV1 ;
   Vector3d vtU2 ;

   if ( abs( vtV2.y) < EPS_SMALL) {

      if ( vtV2.x > 0) 
         
         vtU2 = X_AX ;
      else 
         vtU2 = - X_AX ;
   }
   else {

      vtU2 = vtMovexy ;  vtU2.Normalize() ;
   }

   

  // Determinazione limiti sugli indici
   unsigned int nStartI = ( dMinX < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinX / m_dStep)) ;
   unsigned int nEndI = ( dMaxX > m_nNx * m_dStep - EPS_SMALL ? m_nNx - 1 : static_cast <unsigned int> ( dMaxX / m_dStep)) ;
   unsigned int nStartJ = ( dMinY < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinY / m_dStep)) ;
   unsigned int nEndJ = ( dMaxY > m_nNy * m_dStep - EPS_SMALL ? m_nNy - 1 : static_cast <unsigned int> ( dMaxY / m_dStep)) ;

  // Ciclo
   for ( unsigned int i = nStartI ; i <= nEndI ; ++ i)
      for ( unsigned int j = nStartJ ; j <= nEndJ ; ++ j) {

         double dMin, dMax ;

         double dX = ( i + 0.5) * m_dStep ;  double dY = ( j + 0.5) * m_dStep ;
      
         Point3d ptC( dX, dY, 0) ;

         Vector3d vtC = ptC - ptCIxy ;

         double dProj1 = vtC * vtU1 ;  double dProj2 = vtC * vtU2 ;

         double dSqRadDistI = dProj1 * dProj1 + dProj2 * dProj2 ;
         double dSqRadDistF = dProj1 * dProj1 + ( dProj2 - dPlaneLen) * ( dProj2 - dPlaneLen) ;
         double dSqAxDist = dProj1 * dProj1 ; 

         if ( ( dProj2 < 0 && dSqRadDistI < m_dRadius * m_dRadius) || 
              ( dProj2 >= 0 && dProj2 < dPlaneLen && dSqAxDist < m_dRadius * m_dRadius) || 
              ( dProj2 >= dPlaneLen && dSqRadDistF < m_dRadius * m_dRadius)) {

           // Massimi
            if ( dProj2 < - dSemiAxMin * sqrt( 1 - ( dProj1 * dProj1) / ( m_dRadius * m_dRadius))) {
            
               double dH = sqrt( m_dRadius * m_dRadius - dSqRadDistI) ;

               dMax = dZI + dH ;
            }
            else if ( dProj2 >=  - dSemiAxMin * sqrt( 1 - ( dProj1 * dProj1) / ( m_dRadius * m_dRadius)) &&
                      dProj2 < dPlaneLen - dSemiAxMin * sqrt( 1 - ( dProj1 * dProj1) / ( m_dRadius * m_dRadius))) {

               double dProj0 = - dSemiAxMin * sqrt( 1 - ( dProj1 * dProj1) / ( m_dRadius * m_dRadius)) ;
               double dZ0 = dZI + ( dPlaneLen / dDeltaZ) * ( - dProj0) ;

               dMax = dZ0  + ( dDeltaZ / dPlaneLen) * ( dProj2 - dProj0) ;
            }
            else {

               double dH = sqrt( m_dRadius * m_dRadius - dSqRadDistF) ;

               dMax = dZF + dH ;
            }

           // Minimi
            if ( dProj2 < dSemiAxMin * sqrt( 1 - ( dProj1 * dProj1) / ( m_dRadius * m_dRadius))) {
            
               double dH = sqrt( m_dRadius * m_dRadius - dSqRadDistI) ;

               dMin = dZI - dH ;
            }
            else if ( dProj2 >= dSemiAxMin * sqrt( 1 - ( dProj1 * dProj1) / ( m_dRadius * m_dRadius)) &&
                      dProj2 < dPlaneLen + dSemiAxMin * sqrt( 1 - ( dProj1 * dProj1) / ( m_dRadius * m_dRadius))) {

               double dProj0 = dSemiAxMin * sqrt( 1 - ( dProj1 * dProj1) / ( m_dRadius * m_dRadius)) ;
               double dZ0 = dZI + ( dPlaneLen / dDeltaZ) * ( - dProj0) ;

               dMin = dZ0  + ( dDeltaZ / dPlaneLen) * ( dProj2 - dProj0) ;
            }
            else {

               double dH = sqrt( m_dRadius * m_dRadius - dSqRadDistF) ;

               dMin = dZF - dH ;
            }

            SubtractIntervals( i, j, dMin, dMax) ;
         }
      } // Fine ciclo

      return true ;
}
*/

//----------------------------------------------------------------------------
bool
VolZmap::MillingXYPlusBall( const Point3d& ptLs, const Point3d& ptLe, const Vector3d& vtToolDir) {

  // Parte cilindrica

   double dCylH = m_dHeight - m_dRadius ;

   m_dHeight = dCylH ;

   MillingXYPlusCyl( ptLs, ptLe, vtToolDir) ;

   m_dHeight = m_dHeight + m_dRadius ;
  ////////////////////////////////////////
   Point3d ptI, ptF ;

   if ( ptLs.z < ptLe.z) {
   
      ptI = ptLs ;
      ptF = ptLe ;
   }
   else {
   
      ptI = ptLe ;
      ptF = ptLs ;
   }

   double dZI = ptI.z ;  double dZF = ptF.z ;  double dDeltaZ = dZF - dZI ;

   Point3d ptCI = ptI - dCylH * vtToolDir ;
   Point3d ptCF = ptF - dCylH * vtToolDir ;

   Point3d ptCIxy( ptCI.x, ptCI.y, 0) ;

  // Bounding box 
   double dMinX = min( ptCI.x, ptCF.x) - m_dRadius ;
   double dMaxX = max( ptCI.x, ptCF.x) + m_dRadius ;
   double dMinY = min( ptCI.y, ptCF.y) - m_dRadius ;
   double dMaxY = max( ptCI.y, ptCF.y) + m_dRadius ;

   Vector3d vtMove = ptF - ptI ;
   Vector3d vtPlaneMove( vtMove.x, vtMove.y, 0) ;   double dPLen = vtPlaneMove.LenXY() ;  vtMove.Normalize() ;

  // Sistema di riferimento nel piano
   Vector3d vtV2 = vtPlaneMove ;  vtV2.Normalize() ;  double dComp1 = vtV2 * X_AX ;  double dComp2 = vtV2 * Y_AX ;
   Vector3d vtV1 = dComp2 * X_AX - dComp1 * Y_AX ;

  // Determino il semi-asse minore 
   double dOriz = vtMove * Z_AX ;  double dVert = vtMove * vtV2 ;

   double dSemiAxMin = m_dRadius * dOriz ;

  // Determinazione limiti sugli indici
   unsigned int nStartI = ( dMinX < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinX / m_dStep)) ;
   unsigned int nEndI = ( dMaxX > m_nNx * m_dStep - EPS_SMALL ? m_nNx - 1 : static_cast <unsigned int> ( dMaxX / m_dStep)) ;
   unsigned int nStartJ = ( dMinY < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinY / m_dStep)) ;
   unsigned int nEndJ = ( dMaxY > m_nNy * m_dStep - EPS_SMALL ? m_nNy - 1 : static_cast <unsigned int> ( dMaxY / m_dStep)) ;

  // Ciclo
   for ( unsigned int i = nStartI ; i <= nEndI ; ++ i)
      for ( unsigned int j = nStartJ ; j <= nEndJ ; ++ j) {

         double dMin, dMax ;

         double dX = ( i + 0.5) * m_dStep ;  double dY = ( j + 0.5) * m_dStep ;
      
         Point3d ptC( dX, dY, 0) ;

         Vector3d vtC = ptC - ptCIxy ;

         double dProj1 = vtC * vtV1 ;  double dProj2 = vtC * vtV2 ;

         double dSqRadDistI = dProj1 * dProj1 + dProj2 * dProj2 ;
         double dSqRadDistF = dProj1 * dProj1 + ( dProj2 - dPLen) * ( dProj2 - dPLen) ;
         double dSqAxDist = dProj1 * dProj1 ; 

         if ( ( dProj2 < 0 && dSqRadDistI < m_dRadius * m_dRadius) || 
              ( dProj2 >= 0 && dProj2 < dPLen && dSqAxDist < m_dRadius * m_dRadius) || 
              ( dProj2 >= dPLen && dSqRadDistF < m_dRadius * m_dRadius)) {

           // Massimi
            if ( dProj2 < - dSemiAxMin * sqrt( 1 - ( dSqAxDist) / ( m_dRadius * m_dRadius))) {
            
               double dH = sqrt( m_dRadius * m_dRadius - dSqRadDistI) ;

               dMax = dZI + dH ;
            }
            else if ( dProj2 >=  - dSemiAxMin * sqrt( 1 - ( dSqAxDist) / ( m_dRadius * m_dRadius)) &&
                      dProj2 < dPLen - dSemiAxMin * sqrt( 1 - ( dSqAxDist) / ( m_dRadius * m_dRadius))) {

               double dProj0 = - dSemiAxMin * sqrt( 1 - ( dSqAxDist) / ( m_dRadius * m_dRadius)) ;
               double dZ0 = dZI + dVert * sqrt( m_dRadius * m_dRadius - dSqAxDist) ;

               dMax = dZ0  + ( dDeltaZ / dPLen) * ( dProj2 - dProj0) ;
            }
            else {

               double dH = sqrt( m_dRadius * m_dRadius - dSqRadDistF) ;

               dMax = dZF + dH ;
            }

           // Minimi
            if ( dProj2 < dSemiAxMin * sqrt( 1 - ( dSqAxDist) / ( m_dRadius * m_dRadius))) {
            
               double dH = sqrt( m_dRadius * m_dRadius - dSqRadDistI) ;

               dMin = dZI - dH ;
            }
            else if ( dProj2 >= dSemiAxMin * sqrt( 1 - ( dSqAxDist) / ( m_dRadius * m_dRadius)) &&
                      dProj2 < dPLen + dSemiAxMin * sqrt( 1 - ( dSqAxDist) / ( m_dRadius * m_dRadius))) {

               double dProj0 = dSemiAxMin * sqrt( 1 - ( dSqAxDist) / ( m_dRadius * m_dRadius)) ;
               double dZ0 = dZI - dVert * sqrt( m_dRadius * m_dRadius - dSqAxDist) ;

               dMin = dZ0  + ( dDeltaZ / dPLen) * ( dProj2 - dProj0) ;
            }
            else {

               double dH = sqrt( m_dRadius * m_dRadius - dSqRadDistF) ;

               dMin = dZF - dH ;
            }

            SubtractIntervals( i, j, dMin, dMax) ;
         }
      } // Fine ciclo

      return true ;
}

//----------------------------------------------------------------------------
inline bool
VolZmap::GetMinMaxXY( unsigned int nI, unsigned int nJ, double dProj, double dZheight, double dSqD, double dPathPerp, double dPathPar, double dScProd) {

  // Definisco la variabile altezza della parte cilindrica dell'utensile
   double dCylH ;
  // Definisco variabili per determinazione intervallo da sottrarre
   double dL, dR, dH ;

  // Caso di utensile generico ( per ora non gestito) 
   if ( m_nToolType == 0)
      return false ;
  // Se utensile standard setto altezza della parte cilindrica
   else if ( m_nToolType == 1)
      dCylH = m_dHeight ;
   else if ( m_nToolType == 2)
      dCylH = m_dHeight - m_dRadius ;
   else if ( m_nToolType == 3)
      dCylH = m_dHeight - m_dRCorner ;
  // Caso di utensile non definito 
   else  
      return false ;

  // Parametri per isolare la lavorazione della parte cilindrica comune a tutti gli utensili standard  
  // Limite nella direziona parallela all'asse dell'utensile e parametro perpendicolare
   double dLimPar, dParPerp ;

  // Tagli perpendicolari all'asse dell'utensile
   if ( abs( dScProd) < EPS_SMALL) {
      dLimPar = dCylH ;
      dParPerp = dPathPerp ;
   }
  // Tagli paralleli all'asse dell'utensile
   else {
      dLimPar = dPathPar ;
      dParPerp = 0 ;
   }
   
  // Parte cilindrica della lavorazione
   if ( dSqD < dLimPar * dLimPar) {
     // Qui il raggio del semicerchio è m_dRadius
      if ( dProj < 0) 
        // dH = sqrt( m_dRadius^2 - dProj^2)
         dH = sqrt( m_dRadius * m_dRadius - ( dProj * dProj)) ;
     // Qui l'altezza è costante
      else if ( dProj < dParPerp) 
        // dH = m_dRadius
         dH = m_dRadius ;
     // Qui il raggio del semicerchio è m_dRadius 
      else if ( dProj < dParPerp + m_dRadius) 
        // dH = sqrt( m_dRadius^2 - ( dProj - dParPerp)^2)
         dH = sqrt( m_dRadius * m_dRadius - ( dProj - dParPerp) * ( dProj - dParPerp)) ;
     // Eseguo il taglio
      return SubtractIntervals( nI, nJ, dZheight - dH, dZheight + dH) ;      
   }
  // Parte non cilindrica, questa parte esiste solo nei tagli (dScProd = 0) e nel foro con la punta quindi (dScProd < 0) 
  // e solo con frese non cilindriche
   else {
      if ( dScProd < EPS_SMALL) {
        // Definisco variabili 
        // Caso fresa ball-end
         if ( m_nToolType == 2) {
           // Il raggio è sqrt( m_dRadius^2 - ( sqrt( dSqD) - dLimPar)^2)
            if ( dProj < 0) {
               dL = sqrt( dSqD) - dLimPar ;
               dR = sqrt( m_dRadius * m_dRadius - dL * dL) ;
               dH = sqrt( dR * dR - (dProj * dProj)) ;
            }
           // Qui dH non dipende da dProj
            else if ( dProj < dParPerp) {
               dL = sqrt( dSqD) - dLimPar ; 
               dH = sqrt( m_dRadius * m_dRadius - dL * dL) ;
            }
           // E' analogo al primo caso con la sostituzione di dProj con dProj - dParPerp
            else if ( dProj < dParPerp + m_dRadius) { // In questo caso è equivalente a else
               dL = sqrt( dSqD) - dLimPar ;
               dR = sqrt( m_dRadius * m_dRadius - dL * dL) ;
               dH = sqrt( dR * dR - ((dProj - dParPerp) * (dProj - dParPerp))) ; 
            }
         }
        // Caso di fresa bull-nose
         else if ( m_nToolType == 3) {
           // Raggio semicerchio m_dRadius - m_dRCorner + sqrt( m_dRCorner^2 - ( sqrt( dSqD) - dLimPar)^2)   
            if ( dProj < 0) {
               dL = sqrt( dSqD) - dLimPar ; 
               dR = m_dRadius - m_dRCorner + sqrt( m_dRCorner * m_dRCorner - dL * dL) ;
               dH = sqrt( dR * dR - ( dProj * dProj)) ;
            }
           // Qui dH non dipende da dProj
            else if ( dProj < dParPerp) {
               dL = sqrt( dSqD) - dLimPar ;
               dH = m_dRadius - m_dRCorner + sqrt( m_dRCorner * m_dRCorner - dL * dL) ;
            }
           // E' analogo al primo caso con la sostituzione di di dProj con dProj - dParPerp 
            else if ( dProj < dParPerp + m_dRadius) { // In questo caso equivalente a else
               dL = sqrt( dSqD) - dLimPar ; 
               dR = m_dRadius - m_dRCorner + sqrt( m_dRCorner * m_dRCorner - dL * dL) ;
               dH = sqrt( dR * dR - ( ( dProj - dParPerp) * ( dProj - dParPerp))) ;
            }         
         }
      return SubtractIntervals( nI, nJ, dZheight - dH, dZheight + dH) ;
      }
   } 
   return true ;
}
 
// Frese conus
//----------------------------------------------------------------------------
bool
VolZmap::ConusDrillingZ( const Point3d ptLs, const Point3d ptLe, const Vector3d vtToolDir) {

   double dMinZ = min( min( ptLs.z, ptLs.z - vtToolDir.z * m_dHeight), min( ptLe.z, ptLe.z - vtToolDir.z * m_dHeight)) ;
   double dMaxZ = max( max( ptLs.z, ptLs.z - vtToolDir.z * m_dHeight), max( ptLe.z, ptLe.z - vtToolDir.z * m_dHeight)) ;

  // Prima verifica sull'interferenza dell'utensile con lo Zmap
   if ( dMaxZ < m_dMinZ + EPS_SMALL || dMinZ > m_dMaxZ - EPS_SMALL)
      return true ;

   double dMaxRad = max( m_dRadius, m_dTipRadius) ;
   double dMinRad = min( m_dRadius, m_dTipRadius) ;

  // Bounding box
   double dMinX = ptLs.x - dMaxRad ;
   double dMaxX = ptLs.x + dMaxRad ;
   double dMinY = ptLs.y - dMaxRad ;
   double dMaxY = ptLs.y + dMaxRad ;

  // Seconda verifica dell'interferenza dell'utensile con lo Zmap
   if ( dMaxX < EPS_SMALL || dMinX > m_nNx * m_dStep - EPS_SMALL)
      return true ;
   if ( dMaxY < EPS_SMALL || dMinY > m_nNy * m_dStep - EPS_SMALL)
      return true ;

   double dDeltaZ = abs( ptLe.z - ptLs.z) ;

   Point3d ptO( ptLs.x, ptLs.y, 0) ;

  // Limiti su indici
   unsigned int nStartI = ( dMinX < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinX / m_dStep)) ;
   unsigned int nEndI = ( dMaxX > m_nNx * m_dStep - EPS_SMALL ? m_nNx - 1 : static_cast <unsigned int> ( dMaxX / m_dStep)) ;
   unsigned int nStartJ = ( dMinY < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinY / m_dStep)) ;
   unsigned int nEndJ = ( dMaxY > m_nNy * m_dStep - EPS_SMALL ? m_nNy - 1 : static_cast <unsigned int> ( dMaxY / m_dStep)) ;


   if ( m_dRadius < m_dTipRadius) 

      return GetMinMaxZSw( ptO, nStartI, nEndI, nStartJ, nEndJ, dMinZ, dMaxZ, dMinRad, dMaxRad, vtToolDir.z, dDeltaZ) ;

   else 
      
      return GetMinMaxZDr( ptO, nStartI, nEndI, nStartJ, nEndJ, dMinZ, dMaxZ, dMinRad, dMaxRad, vtToolDir.z, dDeltaZ) ;
}

//----------------------------------------------------------------------------
inline bool
VolZmap::GetMinMaxZSw( const Point3d ptO, unsigned int nStartI, unsigned int nEndI, unsigned int nStartJ, unsigned int nEndJ, double dMinZ, double dMaxZ, double dMinRad, double dMaxRad, double dDir, double dDeltaZ) {

   double dMin, dMax ;

   for ( unsigned int i = nStartI ; i <= nEndI ; ++ i)
         for ( unsigned int j = nStartJ ; j <= nEndJ ; ++ j) {

            double dX = ( i + 0.5) * m_dStep ;  double dY = ( j + 0.5) * m_dStep ; 

            Point3d ptC( dX, dY, 0) ;    Vector3d vtC = ptC - ptO ;  

            double dSqDist = vtC * vtC ;

            if ( dSqDist < dMaxRad * dMaxRad) {

               if ( dSqDist < dMinRad * dMinRad) {

                  dMin = dMinZ ;
                  dMax = dMaxZ ;               
               }
               else {

                  if ( dDir > 0) {

                     dMin = dMinZ ;
                     dMax = dMinZ + dDeltaZ + ( m_dTipHeight * ( dMaxRad - sqrt( dSqDist))) / ( dMaxRad - dMinRad) ;
                  }
                  else {

                     dMin = dMaxZ - dDeltaZ + ( m_dTipHeight * ( sqrt( dSqDist) - dMaxRad)) / ( dMaxRad - dMinRad) ;
                     dMax = dMaxZ ;        
                  }
               }

               SubtractIntervals( i, j, dMin, dMax) ;
            }       
      } // Fine ciclo
      
      return true ;
}

//----------------------------------------------------------------------------
inline bool
VolZmap::GetMinMaxZDr( const Point3d ptO, unsigned int nStartI, unsigned int nEndI, unsigned int nStartJ, unsigned int nEndJ, double dMinZ, double dMaxZ, double dMinRad, double dMaxRad, double dDir, double dDeltaZ)  {

   double dMin, dMax ;

   for ( unsigned int i = nStartI ; i <= nEndI ; ++ i)
         for ( unsigned int j = nStartJ ; j <= nEndJ ; ++ j) {

            double dX = ( i + 0.5) * m_dStep ;  double dY = ( j + 0.5) * m_dStep ; 

            Point3d ptC( dX, dY, 0) ;    Vector3d vtC = ptC - ptO ;  

            double dSqDist = vtC * vtC ;

            if ( dSqDist < dMaxRad * dMaxRad) {

               if ( dSqDist < dMinRad * dMinRad) {

                  dMin = dMinZ ;
                  dMax = dMaxZ ;
               }
               else {

                  if ( dDir > 0) {

                     dMin = dMinZ + ( m_dTipHeight * ( sqrt( dSqDist) - dMinRad)) / ( dMaxRad - dMinRad) ;           
                     dMax = dMaxZ ;
                  }
                  else {

                     dMin = dMinZ ;
                     dMax = dMaxZ + ( m_dTipHeight * ( sqrt( dSqDist) - dMinRad)) / ( dMinRad - dMaxRad) ;              
                  }
               }

               SubtractIntervals( i, j, dMin, dMax) ;
            }    
      } // Fine ciclo
      return true ;
}

//----------------------------------------------------------------------------
bool 
VolZmap::ConusPerpZ( const Point3d ptLs, const Point3d ptLe, const Vector3d vtToolDir) {

   double dMinZ = min( min( ptLs.z, ptLs.z - vtToolDir.z * m_dHeight), min( ptLe.z, ptLe.z - vtToolDir.z * m_dHeight)) ;
   double dMaxZ = max( max( ptLs.z, ptLs.z - vtToolDir.z * m_dHeight), max( ptLe.z, ptLe.z - vtToolDir.z * m_dHeight)) ;

  // Prima verifica sull'interferenza dell'utensile con lo Zmap
   if ( dMaxZ < m_dMinZ + EPS_SMALL || dMinZ > m_dMaxZ - EPS_SMALL)
      return true ;

   double dMaxRad = max( m_dRadius, m_dTipRadius) ;
   double dMinRad = min( m_dRadius, m_dTipRadius) ;

  // Bounding box
   double dMinX = min( ptLs.x, ptLe.x) - dMaxRad ;
   double dMaxX = max( ptLs.x, ptLe.x) + dMaxRad ;
   double dMinY = min( ptLs.y, ptLe.y) - dMaxRad ;
   double dMaxY = max( ptLs.y, ptLe.y) + dMaxRad ;   // Ricordati del caso balordo da mettere nella documentazione

  // Seconda verifica dell'interferenza dell'utensile con lo Zmap
   if ( dMaxX < EPS_SMALL || dMinX > m_nNx * m_dStep - EPS_SMALL)
      return true ;
   if ( dMaxY < EPS_SMALL || dMinY > m_nNy * m_dStep - EPS_SMALL)
      return true ;

   Point3d ptO( ptLs.x, ptLs.y, 0) ;
   Vector3d vtMove = ptLe - ptLs ;  double dLen = vtMove.LenXY() ;  
   Vector3d vtV1 = vtMove ;  vtV1.Normalize() ; 

   double dZBase = ptLs.z ; 
   double dZStem = ptLs.z - ( m_dHeight - m_dTipHeight) * vtToolDir.z ; 
   double dZTip = ptLs.z - m_dHeight * vtToolDir.z ; 

  // Limiti su indici
   unsigned int nStartI = ( dMinX < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinX / m_dStep)) ;
   unsigned int nEndI = ( dMaxX > m_nNx * m_dStep - EPS_SMALL ? m_nNx - 1 : static_cast <unsigned int> ( dMaxX / m_dStep)) ;
   unsigned int nStartJ = ( dMinY < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinY / m_dStep)) ;
   unsigned int nEndJ = ( dMaxY > m_nNy * m_dStep - EPS_SMALL ? m_nNy - 1 : static_cast <unsigned int> ( dMaxY / m_dStep)) ;

   for( unsigned int i = nStartI ; i <= nEndI ; ++ i)
      for ( unsigned int j = nStartJ ; j <= nEndJ ; ++ j) {

         double dMin, dMax ;
         double dX = ( i + 0.5) * m_dStep ;  double dY = ( j + 0.5) * m_dStep ;
      
         Point3d ptC( dX, dY, 0) ;  Vector3d vtC = ptC - ptO ;

         double dProj1 = vtC * vtV1 ;
         Vector3d vtOrt = vtC - dProj1 * vtV1 ;

         if ( dProj1 < 0) {

            double dSqDist = vtC * vtC ; 

            if ( dSqDist < dMinRad * dMinRad) {
            
               dMin = min( dZBase, dZTip) ;
               dMax = max( dZBase, dZTip) ;

               SubtractIntervals( i, j, dMin, dMax) ;
            }
            else if ( dSqDist >= dMinRad * dMinRad && dSqDist < dMaxRad * dMaxRad) {

               if ( m_dRadius < m_dTipRadius) {

                  dMin = min( dZTip, dZStem + ( sqrt( dSqDist) - dMinRad) * ( dZTip - dZStem) / ( dMaxRad - dMinRad)) ;
                  dMax = max( dZTip, dZStem+ ( sqrt( dSqDist) - dMinRad) * ( dZTip - dZStem) / ( dMaxRad - dMinRad)) ;
               }
               else {

                  dMin = min( dZBase, dZTip + ( sqrt( dSqDist) - dMinRad) * ( dZStem - dZTip) / ( dMaxRad - dMinRad)) ;
                  dMax = max( dZBase, dZTip + ( sqrt( dSqDist) - dMinRad) * ( dZStem - dZTip) / ( dMaxRad - dMinRad)) ;
               }
               SubtractIntervals( i, j, dMin, dMax) ;
            }
         }
         else if ( dProj1 >= 0 && dProj1 < dLen) {

            double dSqDist = vtOrt * vtOrt ;

            if ( dSqDist < dMinRad * dMinRad) {

               dMin = min( dZBase, dZTip) ;
               dMax = max( dZBase, dZTip) ; 

               SubtractIntervals( i, j, dMin, dMax) ;
            }
            else if ( dSqDist >= dMinRad * dMinRad && dSqDist < dMaxRad * dMaxRad) {

               if ( m_dRadius < m_dTipRadius) {

                   dMin = min( dZTip, dZStem + ( sqrt( dSqDist) - dMinRad) * ( dZTip - dZStem) / ( dMaxRad - dMinRad)) ;
                   dMax = max( dZTip, dZStem + ( sqrt( dSqDist) - dMinRad) * ( dZTip - dZStem) / ( dMaxRad - dMinRad)) ;
               }
               else {
               
                   dMin = min( dZBase, dZTip + ( sqrt( dSqDist) - dMinRad) * ( dZStem - dZTip) / ( dMaxRad - dMinRad)) ;
                   dMax = max( dZBase, dZTip + ( sqrt( dSqDist) - dMinRad) * ( dZStem - dZTip) / ( dMaxRad - dMinRad)) ;
               }
               SubtractIntervals( i, j, dMin, dMax) ;
            }
         }
         else if ( dProj1 >= dLen) {

            Vector3d vtCn = vtC - vtMove ;

            double dSqDist = vtCn * vtCn ;

            if ( dSqDist < dMinRad * dMinRad) {

               dMin = min( dZBase, dZTip) ;
               dMax = max( dZBase, dZTip) ;

               SubtractIntervals( i, j, dMin, dMax) ;
            }
            else if ( dSqDist >= dMinRad * dMinRad && dSqDist < dMaxRad * dMaxRad) {

               if ( m_dRadius < m_dTipRadius) {

                   dMin = min( dZTip, dZStem + ( sqrt( dSqDist) - dMinRad) * ( dZTip - dZStem) / ( dMaxRad - dMinRad)) ;
                   dMax = max( dZTip, dZStem + ( sqrt( dSqDist) - dMinRad) * ( dZTip - dZStem) / ( dMaxRad - dMinRad)) ;
               }
               else {

                   dMin = min( dZBase, dZTip + ( sqrt( dSqDist) - dMinRad) * ( dZStem - dZTip) / ( dMaxRad - dMinRad)) ;
                   dMax = max( dZBase, dZTip + ( sqrt( dSqDist) - dMinRad) * ( dZStem - dZTip) / ( dMaxRad - dMinRad)) ;
               }
               SubtractIntervals( i, j, dMin, dMax) ;
            }
         }
      } // Fine ciclo
      return true ;
}

//----------------------------------------------------------------------------
bool 
VolZmap::ConusMillingZDr( const Point3d ptLs, const Point3d ptLe, const Vector3d vtToolDir) {

   double dMinZ = min( min( ptLs.z, ptLs.z - vtToolDir.z * m_dHeight), min( ptLe.z, ptLe.z - vtToolDir.z * m_dHeight)) ;
   double dMaxZ = max( max( ptLs.z, ptLs.z - vtToolDir.z * m_dHeight), max( ptLe.z, ptLe.z - vtToolDir.z * m_dHeight)) ;

  // Prima verifica sull'interferenza dell'utensile con lo Zmap
   if ( dMaxZ < m_dMinZ + EPS_SMALL || dMinZ > m_dMaxZ - EPS_SMALL)
      return true ;

   double dMaxRad = max( m_dRadius, m_dTipRadius) ;
   double dMinRad = min( m_dRadius, m_dTipRadius) ;

  // Bounding box
   double dMinX = min( ptLs.x, ptLe.x) - dMaxRad ;
   double dMaxX = max( ptLs.x, ptLe.x) + dMaxRad ;
   double dMinY = min( ptLs.y, ptLe.y) - dMaxRad ;
   double dMaxY = max( ptLs.y, ptLe.y) + dMaxRad ;  

  // Seconda verifica dell'interferenza dell'utensile con lo Zmap
   if ( dMaxX < EPS_SMALL || dMinX > m_nNx * m_dStep - EPS_SMALL)
      return true ;
   if ( dMaxY < EPS_SMALL || dMinY > m_nNy * m_dStep - EPS_SMALL)
      return true ;

  // Punti iniziale e finale e proiezione sul piano del punto iniziale
   Point3d ptI, ptF, ptO ;

   if( vtToolDir.z > 0) {

      ptI = ( ptLs.z < ptLe.z ? ptLs : ptLe) ;  ptO.x = ptI.x ;  ptO.y = ptI.y ;  ptO.z = 0 ;
      ptF = ( ptLs.z < ptLe.z ? ptLe : ptLs) ; 
   }
   else {

      ptI = ( ptLs.z < ptLe.z ? ptLe : ptLs) ;  ptO.x = ptI.x ;  ptO.y = ptI.y ;  ptO.z = 0 ;
      ptF = ( ptLs.z < ptLe.z ? ptLs : ptLe) ;
   }
   

  // Quote iniziale e finale della base dell'utensile e DeltaZ
   double dZI = ptI.z ;  double dZF = ptF.z ;  double dDeltaZ = dZF - dZI ;
  
  // Vettori di movimento
   Vector3d vtMove = ptF - ptI ;
   Vector3d vtMoveXY( vtMove.x, vtMove.y, 0) ;   double dLen = vtMoveXY.LenXY() ;

  // Sistema di riferimento sul cono e vertice del cono
   Vector3d vtV1 = vtToolDir ;
   Vector3d vtV2 = vtMoveXY ; vtV2.Normalize() ;
   Vector3d vtV3 = vtV1 ^ vtV2 ;
   Point3d ptV = ptI - ( m_dHeight + m_dTipHeight * dMinRad / ( dMaxRad - dMinRad)) * vtV1 ;

  // Apertura del cono e parametri per determinare i piani
   double dTanAlpha = ( dMaxRad - dMinRad) / m_dTipHeight ;
   double dRatio = ( vtMove * vtV1) / ( vtMove * vtV2) ;
   double dCos = dTanAlpha * dRatio ;
  
  // Versori normali e prodotti scalari per per determinare i piani
   Vector3d vtNs = - ( dTanAlpha / sqrt( 1 + dTanAlpha * dTanAlpha)) * vtV1 + ( dCos / sqrt( 1 + dTanAlpha * dTanAlpha)) * vtV2 + ( sqrt( 1 - dCos * dCos) / sqrt( 1 + dTanAlpha * dTanAlpha)) * vtV3 ;
   Vector3d vtNd = - ( dTanAlpha / sqrt( 1 + dTanAlpha * dTanAlpha)) * vtV1 + ( dCos / sqrt( 1 + dTanAlpha * dTanAlpha)) * vtV2 - ( sqrt( 1 - dCos * dCos) / sqrt( 1 + dTanAlpha * dTanAlpha)) * vtV3 ;
   Vector3d vtR0 = ptV - ORIG ;
   double dDots = vtR0 * vtNs ;
   double dDotd = vtR0 * vtNd ;

  // Limiti su indici
   unsigned int nStartI = ( dMinX < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinX / m_dStep)) ;
   unsigned int nEndI = ( dMaxX > m_nNx * m_dStep - EPS_SMALL ? m_nNx - 1 : static_cast <unsigned int> ( dMaxX / m_dStep)) ;
   unsigned int nStartJ = ( dMinY < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinY / m_dStep)) ;
   unsigned int nEndJ = ( dMaxY > m_nNy * m_dStep - EPS_SMALL ? m_nNy - 1 : static_cast <unsigned int> ( dMaxY / m_dStep)) ;

   for ( unsigned int i = nStartI ; i <= nEndI ; ++ i)
      for ( unsigned int j = nStartJ ; j <= nEndJ ; ++ j) {

         double dMin, dMax ;
         double dX = ( i + 0.5) * m_dStep ;  double dY = ( j + 0.5) * m_dStep ; 

         Point3d ptC( dX, dY, 0) ;   Vector3d vtC = ptC - ptO ;  Vector3d vtCf = vtC - vtMoveXY ;
         
         Vector3d vtUC = vtC ;  Vector3d vtUCf = vtCf ;   vtUC.Normalize() ;  vtUCf.Normalize() ; 

         double dCCos = vtUC * vtV2 ;  double dCCosf = vtUCf * vtV2 ;

         double dProj = vtC * vtV2 ;
         Vector3d vtOrt = vtC - dProj * vtV2 ;  

         double dSqDistI = vtC * vtC ;  double dSqDistM = vtOrt * vtOrt ;  double dSqDistF = vtCf * vtCf ;

         if ( ( dProj < 0 && dSqDistI < dMaxRad * dMaxRad) 
            ||( dProj >= 0 && dProj < dLen && dSqDistM < dMaxRad * dMaxRad) 
            ||( dProj >= dLen && dSqDistF < dMaxRad * dMaxRad)) {

            if ( vtToolDir.z > 0) {
              // Massimi
               double dPMaxI = ( dMaxRad * dMaxRad - dSqDistM > 0 ? sqrt( dMaxRad * dMaxRad - dSqDistM) : 0) ;

               if ( dProj < dLen - dPMaxI) 

                  dMax = dZI + ( dDeltaZ / dLen) * ( dProj + dPMaxI) ;                    
               else            
                  dMax = dZF ;

              // Minimi
               if ( dSqDistI < dMinRad * dMinRad) 

                  dMin = dZI - m_dHeight ;
                  
               else {

                  if ( ( vtMove * vtV1) / ( vtMove * vtV2) <= 1 / dTanAlpha) {
                  
                     if ( dSqDistI < dMaxRad * dMaxRad && dCCos < dCos) 

                        dMin = dZI - m_dHeight + ( ( sqrt( dSqDistI) - dMinRad) * m_dTipHeight) / ( dMaxRad - dMinRad) ;
                        
                     else if ( dCCos >= dCos && dCCosf < dCos && dSqDistM < dMaxRad * dMaxRad * ( 1 - dCos * dCos)) {

                        double dMinSql = dMinRad * dMinRad * ( 1 - dCos * dCos) ;
                        double dMaxSql = dMaxRad * dMaxRad * ( 1 - dCos * dCos) ;
                        double dPMinI = ( dMinRad * dMinRad - dSqDistM > 0 ? sqrt( dMinRad * dMinRad - dSqDistM) : 0) ;

                        if ( dSqDistM <= dMinSql)

                           dMin = dZI - m_dHeight + ( dProj - dPMinI) * dDeltaZ / dLen ; 

                        else if ( dSqDistM < dMaxSql) {

                           if ( vtC * vtV3 > 0)

                              dMin = ( dDots - dX * vtNs.x - dY * vtNs.y) / vtNs.z ;
                           else 
                              dMin = ( dDotd - dX * vtNd.x - dY * vtNd.y) / vtNd.z ;
                        }    
                     }
                     else if ( dCCosf >= dCos) {

                        double dPMinI = ( dMinRad * dMinRad - dSqDistM > 0 ? sqrt( dMinRad * dMinRad - dSqDistM) : 0) ;

                        if ( dSqDistF < dMinRad * dMinRad)

                           dMin = dZI - m_dHeight + ( dProj - dPMinI) * dDeltaZ / dLen ;
                        else 
                           dMin = dZF - m_dHeight + ( ( sqrt( dSqDistF) - dMinRad) * m_dTipHeight) / ( dMaxRad - dMinRad) ;
                     }
                     else 

                        dMin = dZI - m_dHeight + m_dTipHeight + ( dProj - dPMaxI) * dDeltaZ / dLen ; 
                  }
                  else {

                     if ( dSqDistI < dMaxRad * dMaxRad) 

                        dMin = dZI - m_dHeight + ( sqrt( dSqDistI) - dMinRad) * m_dTipHeight / ( dMaxRad - dMinRad) ;
                     
                     else if ( dProj >= dPMaxI && dProj < dPMaxI + dLen) 

                        dMin = dZI - m_dHeight + m_dTipHeight + ( dProj - dPMaxI) * dDeltaZ / dLen ;
                  }
               }
            }
            else {
              // Massimi
               double dPMaxI = ( dMaxRad * dMaxRad - dSqDistM > 0 ? sqrt( dMaxRad * dMaxRad - dSqDistM) : 0) ;

               if ( dSqDistI < dMinRad * dMinRad)

                  dMax = dZI + m_dHeight ;
                  
               else {

                  if ( - ( vtMove * vtV1) / ( vtMove * vtV2) <= 1 / dTanAlpha) {

                     if ( dSqDistI < dMaxRad * dMaxRad && dCCos < dCos) 
                     
                        dMax = dZI + m_dHeight - ( sqrt( dSqDistI) - dMinRad) * m_dTipHeight / ( dMaxRad - dMinRad) ;
                     
                     else if ( dCCos >= dCos && dCCosf < dCos && dSqDistM < dMaxRad * dMaxRad * ( 1 - dCos * dCos)) {

                        double dMinSql = dMinRad * dMinRad * ( 1 - dCos * dCos) ;
                        double dMaxSql = dMaxRad * dMaxRad * ( 1 - dCos * dCos) ;
                        double dPMinI = ( dMinRad * dMinRad - dSqDistM > 0 ? sqrt( dMinRad * dMinRad - dSqDistM) : 0) ;

                        if ( dSqDistM <= dMinSql)

                           dMax = dZI + m_dHeight + ( dProj - dPMinI) * dDeltaZ / dLen ; 

                        else if ( dSqDistM < dMaxSql) {

                           if ( vtC * vtV3 > 0)

                              dMax = ( dDots - dX * vtNs.x - dY * vtNs.y) / vtNs.z ;
                           else 
                              dMax = ( dDotd - dX * vtNd.x - dY * vtNd.y) / vtNd.z ;
                        }    
                     }
                     else if ( dCCosf >= dCos) {

                        double dPMinI = ( dMinRad * dMinRad - dSqDistM > 0 ? sqrt( dMinRad * dMinRad - dSqDistM) : 0) ;

                        if ( dSqDistF < dMinRad * dMinRad)

                           dMax = dZI + m_dHeight + ( dProj - dPMinI) * dDeltaZ / dLen ;
                        else 
                           dMax = dZF + m_dHeight - ( ( sqrt( dSqDistF) - dMinRad) * m_dTipHeight) / ( dMaxRad - dMinRad) ;
                     }
                     else 

                        dMax = dZI + m_dHeight - m_dTipHeight + ( dProj - dPMaxI) * dDeltaZ / dLen ; 

                  }
                  else {

                     if ( dSqDistI < dMaxRad * dMaxRad) 

                        dMax = dZI + m_dHeight - ( ( sqrt( dSqDistI) - dMinRad) * m_dTipHeight) / ( dMaxRad - dMinRad) ;
                     
                     else if ( dProj >= dPMaxI && dProj < dPMaxI + dLen) 
                     
                        dMax = dZI + m_dHeight - m_dTipHeight + ( dProj - dPMaxI) * dDeltaZ / dLen ;
                     
                  }   
               }
              // Minimi
               double dPMaxCirc = sqrt( dMaxRad * dMaxRad - dSqDistM) ;

               if ( dProj > dLen - dPMaxCirc) 

                  dMin = dZF ;       
               else
                  dMin = dZI + ( dDeltaZ / dLen) * ( dProj + dPMaxCirc) ; 
            }
            SubtractIntervals( i, j, dMin, dMax) ;
         }
      } 
   return true ;
}

//----------------------------------------------------------------------------
bool 
VolZmap::ConusMillingZSw( const Point3d ptLs, const Point3d ptLe, const Vector3d vtToolDir) {

   double dMinZ = min( min( ptLs.z, ptLs.z - vtToolDir.z * m_dHeight), min( ptLe.z, ptLe.z - vtToolDir.z * m_dHeight)) ;
   double dMaxZ = max( max( ptLs.z, ptLs.z - vtToolDir.z * m_dHeight), max( ptLe.z, ptLe.z - vtToolDir.z * m_dHeight)) ;

  // Prima verifica sull'interferenza dell'utensile con lo Zmap
   if ( dMaxZ < m_dMinZ + EPS_SMALL || dMinZ > m_dMaxZ - EPS_SMALL)
      return true ;

   double dMaxRad = max( m_dRadius, m_dTipRadius) ;
   double dMinRad = min( m_dRadius, m_dTipRadius) ;

  // Bounding box
   double dMinX = min( ptLs.x, ptLe.x) - dMaxRad ;
   double dMaxX = max( ptLs.x, ptLe.x) + dMaxRad ;
   double dMinY = min( ptLs.y, ptLe.y) - dMaxRad ;
   double dMaxY = max( ptLs.y, ptLe.y) + dMaxRad ;  

  // Seconda verifica dell'interferenza dell'utensile con lo Zmap
   if ( dMaxX < EPS_SMALL || dMinX > m_nNx * m_dStep - EPS_SMALL)
      return true ;
   if ( dMaxY < EPS_SMALL || dMinY > m_nNy * m_dStep - EPS_SMALL)
      return true ;

  // Punti iniziale e finale e proiezione sul piano del punto iniziale
   Point3d ptI, ptF, ptO ;

   if( vtToolDir.z > 0) {

      ptI = ( ptLs.z < ptLe.z ? ptLe : ptLs) ;  ptO.x = ptI.x ;  ptO.y = ptI.y ;  ptO.z = 0 ;
      ptF = ( ptLs.z < ptLe.z ? ptLs : ptLe) ; 
   }
   else {

      ptI = ( ptLs.z < ptLe.z ? ptLs : ptLe) ;  ptO.x = ptI.x ;  ptO.y = ptI.y ;  ptO.z = 0 ;
      ptF = ( ptLs.z < ptLe.z ? ptLe : ptLs) ;
   }

   Point3d ptICyl = ( ptLs.z < ptLe.z ? ptLs : ptLe) ;

   double dDeltaH = m_dHeight - m_dTipHeight ;
   double dMinZCyl = min( ptICyl.z, ptICyl.z - vtToolDir.z * dDeltaH) ;
  

  // Punti contatto cono cilindro
   Point3d ptIS = ptI - ( m_dHeight - m_dTipHeight) * vtToolDir ;  
   Point3d ptFS = ptF - ( m_dHeight - m_dTipHeight) * vtToolDir ;    
 
  // Quote iniziali e finali e DeltaZ
   double dZI = ptI.z ;  double dZF = ptF.z ;   double dDeltaZ = dZF - dZI ; double dADeltaZ = abs( dDeltaZ) ;
   double dZIS = ptIS.z ;  double dZFS = ptFS.z ;  
  
  // Vettori di movimento
   Vector3d vtMove = ptF - ptI ;
   Vector3d vtMoveXY( vtMove.x, vtMove.y, 0) ;   double dLen = vtMoveXY.LenXY() ;

  // Sistema di riferimento sul cono e vertice del cono
   Vector3d vtV1 = - vtToolDir ;
   Vector3d vtV2 = vtMoveXY ; vtV2.Normalize() ;
   Vector3d vtV3 = vtV1 ^ vtV2 ;
   Point3d ptV = ptI + ( m_dHeight - m_dTipHeight * ( 1 + dMinRad / ( dMaxRad - dMinRad))) * vtV1 ;

  // Apertura del cono e parametri per determinare i piani
   double dTanAlpha = ( dMaxRad - dMinRad) / m_dTipHeight ;
   double dRatio = ( vtMove * vtV1) / ( vtMove * vtV2) ;
   double dCos = dTanAlpha * dRatio ;
  
  // Versori normali e prodotti scalari per per determinare i piani
   Vector3d vtNp = - ( dTanAlpha / sqrt( 1 + dTanAlpha * dTanAlpha)) * vtV1 + ( dCos / sqrt( 1 + dTanAlpha * dTanAlpha)) * vtV2 + ( sqrt( 1 - dCos * dCos) / sqrt( 1 + dTanAlpha * dTanAlpha)) * vtV3 ;
   Vector3d vtNm = - ( dTanAlpha / sqrt( 1 + dTanAlpha * dTanAlpha)) * vtV1 + ( dCos / sqrt( 1 + dTanAlpha * dTanAlpha)) * vtV2 - ( sqrt( 1 - dCos * dCos) / sqrt( 1 + dTanAlpha * dTanAlpha)) * vtV3 ;
   Vector3d vtR0 = ptV - ORIG ;
   double dDotp = vtR0 * vtNp ;
   double dDotm = vtR0 * vtNm ;

  // Limiti su indici
   unsigned int nStartI = ( dMinX < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinX / m_dStep)) ;
   unsigned int nEndI = ( dMaxX > m_nNx * m_dStep - EPS_SMALL ? m_nNx - 1 : static_cast <unsigned int> ( dMaxX / m_dStep)) ;
   unsigned int nStartJ = ( dMinY < EPS_SMALL ? 0 : static_cast <unsigned int> ( dMinY / m_dStep)) ;
   unsigned int nEndJ = ( dMaxY > m_nNy * m_dStep - EPS_SMALL ? m_nNy - 1 : static_cast <unsigned int> ( dMaxY / m_dStep)) ;

   for ( unsigned int i = nStartI ; i <= nEndI ; ++i)
      for ( unsigned int j = nStartJ ; j <= nEndJ ; ++j) {

         double dMin, dMax ;
         double dX = ( i + 0.5) * m_dStep ;  double dY = ( j + 0.5) * m_dStep ; 

         Point3d ptC( dX, dY, 0) ;   Vector3d vtC = ptC - ptO ;  Vector3d vtCf = vtC - vtMoveXY ;
         
         Vector3d vtUC = vtC ;  Vector3d vtUCf = vtCf ;   vtUC.Normalize() ;  vtUCf.Normalize() ; 

         double dCCos = vtUC * vtV2 ;  double dCCosf = vtUCf * vtV2 ;

         double dProj = vtC * vtV2 ;
         Vector3d vtOrt = vtC - dProj * vtV2 ;  

         double dSqDistI = vtC * vtC ;  double dSqDistM = vtOrt * vtOrt ;  double dSqDistF = vtCf * vtCf ;

         if ( ( dProj < 0 && dSqDistI < dMaxRad * dMaxRad) 
            ||( dProj >= 0 && dProj < dLen && dSqDistM < dMaxRad * dMaxRad) 
            ||( dProj >= dLen && dSqDistF < dMaxRad * dMaxRad)) {

            if ( vtV1.z < 0) {

               double dPMaxI = ( dMaxRad * dMaxRad - dSqDistM > 0 ? sqrt( dMaxRad * dMaxRad - dSqDistM) : 0) ; 
              
              // Massimi
               if ( dRatio <= 1 / dTanAlpha) {

                  if ( dSqDistI < dMinRad * dMinRad) 

                     dMax = dZIS ;

                  else {
               
                     if ( dSqDistI < dMaxRad * dMaxRad && dCCos < dCos) 
                  
                        dMax = dZIS - ( sqrt( dSqDistI) - dMinRad) * m_dTipHeight / ( dMaxRad - dMinRad) ;
                  
                     else if ( dCCos >= dCos && dCCosf < dCos && dSqDistM < dMaxRad * dMaxRad * ( 1 - dCos * dCos)) {

                        double dMinSql = dMinRad * dMinRad * ( 1 - dCos * dCos) ;
                        double dMaxSql = dMaxRad * dMaxRad * ( 1 - dCos * dCos) ;
                        double dPMinI = ( dMinRad * dMinRad - dSqDistM > 0 ? sqrt( dMinRad * dMinRad - dSqDistM) : 0) ;

                        if ( dSqDistM <= dMinSql)

                           dMax = dZIS + ( dProj - dPMinI) * dDeltaZ / dLen ; 

                        else if ( dSqDistM < dMaxSql) {

                           if ( vtC * vtV3 > 0)

                              dMax = ( dDotp - dX * vtNp.x - dY * vtNp.y) / vtNp.z ;
                           else 
                              dMax = ( dDotm - dX * vtNm.x - dY * vtNm.y) / vtNm.z ;
                        }    
                     }
                     else if ( dCCosf >= dCos) {

                        double dPMinI = ( dMinRad * dMinRad - dSqDistM > 0 ? sqrt( dMinRad * dMinRad - dSqDistM) : 0) ;

                        if ( dSqDistF < dMinRad * dMinRad)

                           dMax = dZIS + ( dProj - dPMinI) * dDeltaZ / dLen ;
                        else 
                           dMax = dZFS - ( ( sqrt( dSqDistF) - dMinRad) * m_dTipHeight) / ( dMaxRad - dMinRad) ;
                     }
                     else 
                        dMax = dZIS - m_dTipHeight + ( dProj - dPMaxI) * dDeltaZ / dLen ;
                  }
               }
               else {
               
                  if ( dSqDistI < dMinRad * dMinRad)

                     dMax = dZIS ;
                  
                  else if ( dSqDistI < dMaxRad * dMaxRad)

                     dMax = dZIS - ( sqrt( dSqDistI) - dMinRad) * m_dTipHeight / ( dMaxRad - dMinRad) ;

                  else {
                  
                     if ( dProj >= dPMaxI)

                        dMax = dZIS - m_dTipHeight + ( dProj - dPMaxI) * dDeltaZ / dLen ;
                  }
               }
              // Minimi
               
               if ( dSqDistF < dMaxRad * dMaxRad)

                  dMin = dZFS - m_dTipHeight ;

               else if ( dProj <= dLen - dPMaxI)

                  dMin = dZIS - m_dTipHeight + ( dProj + dPMaxI) * dDeltaZ / dLen ;
               
            } // Fine vtV1 < 0
            else {

               double dPMaxI = ( dMaxRad * dMaxRad - dSqDistM > 0 ? sqrt( dMaxRad * dMaxRad - dSqDistM) : 0) ; 

              // Massimi 
               if ( dSqDistF < dMaxRad * dMaxRad) 

                  dMax = dZFS + m_dTipHeight ;
               
               else if ( dProj <= dLen - dPMaxI)

                  dMax = dZIS + m_dTipHeight + ( dProj + dPMaxI) * dDeltaZ / dLen ;

              // Minimi
               if ( dRatio <= 1 / dTanAlpha) {
                  
                  if ( dSqDistI < dMinRad * dMinRad) 

                     dMin = dZIS ;

                  else {

                     if ( dSqDistI < dMaxRad * dMaxRad && dCCos < dCos)

                        dMin = dZIS + ( sqrt( dSqDistI) - dMinRad) * m_dTipHeight / ( dMaxRad - dMinRad) ;

                     else if ( dCCos >= dCos && dCCosf < dCos && dSqDistM < dMaxRad * dMaxRad * ( 1 - dCos * dCos)) {

                        double dMinSql = dMinRad * dMinRad * ( 1 - dCos * dCos) ;
                        double dMaxSql = dMaxRad * dMaxRad * ( 1 - dCos * dCos) ;
                        double dPMinI = ( dMinRad * dMinRad - dSqDistM > 0 ? sqrt( dMinRad * dMinRad - dSqDistM) : 0) ;

                        if ( dSqDistM <= dMinSql)

                           dMin = dZIS + ( dProj - dPMinI) * dDeltaZ / dLen ; 

                        else if ( dSqDistM < dMaxSql) {

                           if ( vtC * vtV3 > 0)

                              dMin = ( dDotp - dX * vtNp.x - dY * vtNp.y) / vtNp.z ;
                           else 
                              dMin = ( dDotm - dX * vtNm.x - dY * vtNm.y) / vtNm.z ;
                        }    
                     }
                     else if ( dCCosf >= dCos) {

                        double dPMinI = ( dMinRad * dMinRad - dSqDistM > 0 ? sqrt( dMinRad * dMinRad - dSqDistM) : 0) ;

                        if ( dSqDistF < dMinRad * dMinRad)

                           dMin = dZIS + ( dProj - dPMinI) * dDeltaZ / dLen ;
                        else 
                           dMin = dZFS + ( ( sqrt( dSqDistF) - dMinRad) * m_dTipHeight) / ( dMaxRad - dMinRad) ;
                     }
                     else 
                        dMin = dZIS + m_dTipHeight + ( dProj - dPMaxI) * dDeltaZ / dLen ;
                  }
               }
               else {

                  if ( dSqDistI < dMinRad * dMinRad) 

                     dMin = dZIS ;

                  else if ( dSqDistI < dMaxRad * dMaxRad)

                     dMin = dZIS + ( sqrt( dSqDistI) - dMinRad) * m_dTipHeight / ( dMaxRad - dMinRad) ;

                  else 

                     dMin = dZIS + m_dTipHeight + ( dProj - dPMaxI) * dDeltaZ / dLen ;
               
               }
            }
            SubtractIntervals( i, j, dMin, dMax) ;
         }

        // Parte cilindrica
         Vector3d vtCyl = ( vtV1.z < 0 ? vtCf : vtC) ;
         Vector3d vtCylf = ( vtV1.z < 0 ? vtC : vtCf) ;
         Vector3d vtMot = ( vtV1.z < 0 ? - vtV2 : vtV2) ;

         double dCylProj = vtCyl * vtMot ;
         double dCylSqDistI = vtCyl * vtCyl ;
         double dCylSqDistF = vtCylf * vtCylf ;
         double dCylSqDistM = ( vtCyl -  dCylProj * vtMot) * ( vtCyl -  dCylProj * vtMot) ;
         

         if (  dCylSqDistI < dMinRad * dMinRad 
            || dCylSqDistF < dMinRad * dMinRad  
            ||   ( dCylProj > 0 && dCylProj < dLen && dCylSqDistM < dMinRad * dMinRad)) {

            double dSt = sqrt( dMinRad * dMinRad - dCylSqDistM) ;
           
           // Minimi 
            if ( dCylSqDistI < dMinRad * dMinRad )

               dMin = dMinZCyl ;

            else if ( dCylProj >= dSt)

               dMin = dMinZCyl + ( dCylProj - dSt) *  dADeltaZ / dLen ;
           
           // Massimi
            if ( dCylSqDistF < dMinRad * dMinRad)

               dMax = dMinZCyl + dDeltaH + dADeltaZ ;

            else if ( dCylProj <= dLen - dSt)

               dMax = dMinZCyl + dDeltaH + ( dCylProj + dSt) * dADeltaZ / dLen ;

            SubtractIntervals( i, j, dMin, dMax) ;
         
         }
      }
      return true ;
}



// Utensile generico
//----------------------------------------------------------------------------
inline bool
VolZmap::GetMinMaxZGenTool( unsigned int nI, unsigned int nJ, double dZCutBase, double dDeltaZ, double dSqDist, const Vector3d& vtToolDir) {
  
  // Estremi dell'intervallo da asportare
   double dMin, dMax ;

   unsigned int CurveNum = m_ToolOutline.GetCurveCount() ;  // Non dovrebbe servire
  // Puntatore a un oggetto costante di tipo ICurve
   const ICurve* pCurve ;
  // L'inizio del ciclo dipende dal fatto che lavori la punta o il fondo dell'utensile
   if ( vtToolDir.z * dDeltaZ < 0)
      pCurve = m_ToolOutline.GetLastCurve() ;
   else
      pCurve = m_ToolOutline.GetFirstCurve() ;
 
  // Ciclo
   while ( pCurve != nullptr) {
   
      int nCrvType = pCurve -> GetType() ;

      const ICurve* pTest = pCurve ;
         
      if ( vtToolDir.z * dDeltaZ < 0)
         pTest = m_ToolOutline.GetNextCurve() ;
      else 
         pTest = m_ToolOutline.GetPrevCurve() ;

     // Caso di segmento
      if ( nCrvType == CRV_LINE) {
         
         Point3d ptStart, ptEnd ;
         pCurve->GetStartPoint( ptStart) ;
         pCurve->GetEndPoint( ptEnd) ;
         double dZStart, dMDPrev, angC, angCPrev ;

         if ( abs( ptEnd.x - ptStart.x) > EPS_SMALL) {

            if ( pTest == nullptr) {

               if ( vtToolDir.z * dDeltaZ < 0) {
                  dZStart = ptEnd.y ;
                  angC = ( ptStart.y - ptEnd.y) / ( ptStart.x - ptEnd.x) ;
                  angCPrev = angC ;
               }
               else {
                  dZStart = ptStart.y ;
                  angC = ( ptEnd.y - ptStart.y) / ( ptEnd.x - ptStart.x) ;
                  angCPrev = angC ;
               }  
               
               dMin = min( dZStart + angC * dSqDist, dZStart + angC * dSqDist + dDeltaZ) ;
               dMax = max( dZStart + angC * dSqDist, dZStart + angC * dSqDist + dDeltaZ) ;
               dMDPrev = max( ptStart.x, ptEnd.x) ;
              // Eseguo l'asportazione
               SubtractIntervals( nI, nJ, dMin, dMax) ;
            }
            else {

               double dMD ;

               if ( vtToolDir.z * dDeltaZ < 0) {

                  dZStart = ptEnd.y ;
                  angC = ( ptStart.y - ptEnd.y) / ( ptStart.x - ptEnd.x) ;
                  dMD = max( ptStart.x, ptEnd.x) ;

                  if ( angCPrev < 0 && dMDPrev > dMD) {
                      dMin = min( dZStart + angC * dSqDist, dZStart + angC * dSqDist + dDeltaZ) ;
                      dMax = max( dZStart + angC * dSqDist, dZStart + angC * dSqDist + dDeltaZ) ;
                      dMDPrev = max( ptStart.x, ptEnd.x) ;
                  }

                  angCPrev = angC ;
                  dMDPrev = dMD ;
               }
               else {

                  dZStart = ptStart.y ;
                  angC = ( ptEnd.y - ptStart.y) / ( ptEnd.x - ptStart.x) ;
                  dMD = max( ptStart.x, ptEnd.x) ;

                   if ( angCPrev < 0 && dMDPrev > max( ptStart.x, ptEnd.x)) {
                      dMin = min( dZStart + angC * dSqDist, dZStart + angC * dSqDist + dDeltaZ) ;
                      dMax = max( dZStart + angC * dSqDist, dZStart + angC * dSqDist + dDeltaZ) ;
                      dMDPrev = max( ptStart.x, ptEnd.x) ;
                  }

                  angCPrev = angC ;
                  dMDPrev = dMD ;
               }     
            }
         }
      }
     // Caso di arco
      else if ( nCrvType == CRV_ARC) {
          return true ;
      }
      else 
     // Caso di Bézier et al..
         return true ;
      
      }
     // Decremento/incremento il puntatore
      if ( vtToolDir.z * dDeltaZ < 0)
         pCurve = m_ToolOutline.GetPrevCurve() ;
      else
         pCurve = m_ToolOutline.GetNextCurve() ;

      // Ricorda di guardare che tutti i casi abbiano il relatrivo enunciato di return
}