//----------------------------------------------------------------------------
//  EgalTech 2015-2020
//----------------------------------------------------------------------------
// File : VolZmap.cpp         Data : 05.10.20           Versione : 2.2j1
// Contenuto : Implementazione della classe Volume Zmap (tre griglie)
//
//
//
// Modifiche : 22.01.15 DS Creazione modulo.
//
//
//----------------------------------------------------------------------------

//--------------------------- Include ----------------------------------------

#include "stdafx.h"
#include "CurveLine.h"
#include "CurveArc.h"
#include "VolZmap.h"
#include "GeoConst.h"
#include "/EgtDev/Include/EGkStringUtils3d.h"
#include "/EgtDev/Include/EgtNumUtils.h"
#include "/EgtDev/Include/EgtPerfCounter.h"
#include <future>

using namespace std ;

// ------------------------- OPERAZIONI SU INTERVALLI --------------------------------------------------------------------------------------

//----------------------------------------------------------------------------
bool
VolZmap::SubtractIntervals( int nGrid, int nI, int nJ,
                            double dMin, double dMax, const Vector3d& vtNMin, const Vector3d& vtNMax,
                            int nToolNum, bool bSkipSwap)
{
  // Controllo che dMin e dMax non siano quasi coincidenti
   if ( abs( dMax - dMin) < EPS_ZERO)
      return true ;

  // Controllo che il numero di griglia sia entro i limiti
   if ( nGrid < 0  ||  nGrid > 2)
      return false ;

  // Controllo che gli indici nI, nJ siano entro i limiti
   if ( nI < 0  ||  nI >= m_nNx[nGrid]  || 
        nJ < 0  ||  nJ >= m_nNy[nGrid])
      return false ;

  // Controllo che dMin < dMax
   Vector3d vtNmi = vtNMin ;
   Vector3d vtNma = vtNMax ;
   if ( dMin > dMax) {
      swap( dMin, dMax) ;
      swap( vtNmi, vtNma) ;
   }

  // Riporto le coordinate cicliche delle normali nell'ordine di partenza (da griglia a sistema intrinseco)
   if ( !bSkipSwap  &&  nGrid == 1) {
      swap( vtNmi.x, vtNmi.z) ;
      swap( vtNmi.y, vtNmi.z) ;
      swap( vtNma.x, vtNma.z) ;
      swap( vtNma.y, vtNma.z) ;
   }
   else if ( !bSkipSwap  &&  nGrid == 2) {
      swap( vtNmi.y, vtNmi.z) ;
      swap( vtNmi.x, vtNmi.z) ;
      swap( vtNma.y, vtNma.z) ;
      swap( vtNma.x, vtNma.z) ;
   }

  // Recupero dexel da modificare
   int nPos = nJ * m_nNx[nGrid] + nI ;
   vector<Data>& vDexel = m_Values[nGrid][nPos] ;

  // Ciclo sugli intervalli del dexel (ordinati in senso crescente)
   bool bModified = false ;
   for ( int i = 0 ; i < int( vDexel.size()) ; ++i) {
     // Se interseca l'intervallo corrente
      if ( dMin < vDexel[i].dMax - EPS_ZERO  &&  dMax > vDexel[i].dMin + EPS_ZERO) {
         bModified = true ;
        // se devo limitarlo inferiormente
         if ( dMin <= vDexel[i].dMin + EPS_ZERO  &&  dMax < vDexel[i].dMax - EPS_ZERO) {
            vDexel[i].dMin = dMax ;
            vDexel[i].vtMinN = vtNma ;
            vDexel[i].nToolMin = nToolNum ;
         }
        // se devo limitarlo superiormente
         else if ( dMin > vDexel[i].dMin + EPS_ZERO  &&  dMax >= vDexel[i].dMax - EPS_ZERO) {
            vDexel[i].dMax = dMin ;
            vDexel[i].vtMaxN = vtNmi ;
            vDexel[i].nToolMax = nToolNum ;
         }
        // se devo dividerlo in due parti
         else if ( dMin > vDexel[i].dMin + EPS_ZERO  &&  dMax < vDexel[i].dMax - EPS_ZERO) {
           // inserisco nuovo intervallo (parte superiore)
            vDexel.insert( vDexel.begin() + i + 1,
                           { dMax, vtNma, nToolNum,
                             vDexel[i].dMax, vDexel[i].vtMaxN, vDexel[i].nToolMax,
                             vDexel[i].nCompo}) ;
           // aggiorno il vecchio (parte inferiore)
            vDexel[i].dMax = dMin ;
            vDexel[i].vtMaxN = vtNmi ;
            vDexel[i].nToolMax = nToolNum ;
            ++ i ;
         }
        // altrimenti devo eliminarlo
         else {
            vDexel.erase( vDexel.begin() + i) ;
            -- i ;
         }
      }
     // se è tutto minore dell'intervallo corrente, ho finito
      else if ( dMax <= vDexel[i].dMin + EPS_ZERO)
         break ;
     // altrimenti è tutto maggiore dell'intervallo corrente, passo al successivo
      else   // dMin >= vDexel[i].dMax - EPS_ZERO
         ;
   }

  // Se nessuna modifica, esco
   if ( ! bModified)
      return true ;

  // Elimino residui di intervalli inutili
   for ( int i = 0 ; i < int( vDexel.size()) ; ++ i) {
      int nMaxN = int( floor( ( vDexel[i].dMax + 2 * EPS_SMALL - 0.5 * m_dStep) / m_dStep)) ;
      int nMinN = int( floor( ( vDexel[i].dMin - 2 * EPS_SMALL - 0.5 * m_dStep) / m_dStep)) ;
      if ( nMinN == nMaxN) {
         vDexel.erase( vDexel.begin() + i) ;
         -- i ;
      }
   }

  // Imposto ricalcolo della grafica
   m_OGrMgr.Reset() ;
  // Imposto forma generica 
   m_nShape = GENERIC ;    
  // Imposto ricalcolo numero di componenti connesse
   m_nConnectedCompoCount = - 1 ; 

  // Passo da indici di dexel a indici di voxel
   nI /= m_nDexVoxRatio ;
   nJ /= m_nDexVoxRatio ;

  // Determino quali blocchi sono stati modificati
   if ( nGrid == 0) {
     // Voxel lungo X
      int nXStop = 1 ;
      int nXBlock[2] ;
      nXBlock[0] = min( nI / m_nVoxNumPerBlock, m_nFracLin[0] - 1) ;
      if ( nI % N_VOXBLOCK == 0  &&  nXBlock[0] > 0) {
         nXBlock[1] = nXBlock[0] - 1 ;
         ++ nXStop ;
      }
     // Voxel lungo Y
      int nYStop = 1 ;
      int nYBlock[2] ;
      nYBlock[0] = min( nJ / m_nVoxNumPerBlock, m_nFracLin[1] - 1) ;
      if ( nJ % N_VOXBLOCK == 0  &&  nYBlock[0] > 0) {
         nYBlock[1] = nYBlock[0] - 1 ;
         ++ nYStop ;
      }
     // Voxel lungo Z
      int nVoxNumZ = int( m_nNy[1] / m_nDexVoxRatio + ( m_nNy[1] % m_nDexVoxRatio == 0 ? 1 : 2)) ;  
      int nMinK = Clamp( int( floor( ( ( dMin - 0.5 * m_dStep) / ( m_nDexVoxRatio * m_dStep) - EPS_SMALL))), 0, nVoxNumZ - 2) ;
      int nMaxK = Clamp( int( floor( ( ( dMax + 0.5 * m_dStep) / ( m_nDexVoxRatio * m_dStep) + EPS_SMALL))), 0, nVoxNumZ - 2) ;
      int nMinZBlock = ( m_nMapNum == 1 ? 0 : Clamp( nMinK / int( m_nVoxNumPerBlock), 0, int( m_nFracLin[2] - 1))) ;
      int nMaxZBlock = min( int( m_nFracLin[2] - 1), nMaxK / int( m_nVoxNumPerBlock)) ;
     // Assegno flag ai voxel
      for ( int tI = 0 ; tI < nXStop ; ++ tI) {
         for ( int tJ = 0 ; tJ < nYStop ; ++ tJ) {
            for ( int k = nMinZBlock ; k <= nMaxZBlock ; ++ k) {
               int nBlockNum = k * m_nFracLin[0] * m_nFracLin[1] + nYBlock[tJ] * m_nFracLin[0] + nXBlock[tI] ;
               m_BlockToUpdate[nBlockNum] = true ;
            }
         }
      }
   }

   else if ( nGrid == 1) {
     // Voxel lungo Y
      int nYStop = 1 ;
      int nYBlock[2] ;
      nYBlock[0] = min( nI / m_nVoxNumPerBlock, m_nFracLin[1] - 1) ;
      if ( nI % N_VOXBLOCK == 0  &&  nYBlock[0] > 0) {
         nYBlock[1] = nYBlock[0] - 1 ;
         ++ nYStop ;
      }
     // Voxel lungo Z
      int nZStop = 1 ;
      int nZBlock[2] ;
      nZBlock[0] = min( nJ / m_nVoxNumPerBlock, m_nFracLin[2] - 1) ;
      if ( nJ % N_VOXBLOCK == 0  &&  nZBlock[0] > 0) {
         nZBlock[1] = nZBlock[0] - 1 ;
         ++ nZStop ;
      }
     // Voxel lungo X
      int nVoxNumX = int( m_nNx[0] / m_nDexVoxRatio + ( m_nNx[0] % m_nDexVoxRatio == 0 ? 1 : 2)) ;
      int nMinI = Clamp( int( floor( ( ( dMin - 0.5 * m_dStep) / ( m_nDexVoxRatio * m_dStep) - EPS_SMALL))), 0, nVoxNumX - 2) ;
      int nMaxI = Clamp( int( floor( ( ( dMax + 0.5 * m_dStep) / ( m_nDexVoxRatio * m_dStep) + EPS_SMALL))), 0, nVoxNumX - 2) ;
      int nMinXBlock = Clamp( nMinI / int( m_nVoxNumPerBlock), 0, int( m_nFracLin[0] - 1)) ;
      int nMaxXBlock = min( int( m_nFracLin[0] - 1), nMaxI / int( m_nVoxNumPerBlock)) ;
     // Assegno flag ai voxel
      for ( int tI = 0 ; tI < nYStop ; ++ tI) {
         for ( int tJ = 0 ; tJ < nZStop ; ++ tJ) {
            for ( int k = nMinXBlock ; k <= nMaxXBlock ; ++ k) {
               int nBlockNum = nZBlock[tJ] * m_nFracLin[0] * m_nFracLin[1] + nYBlock[tI] * m_nFracLin[0] + k ;
               m_BlockToUpdate[nBlockNum] = true ;
            }
         }
      } 
   }

   else if ( nGrid == 2) {
     // Voxel lungo X
      int nXStop = 1 ;
      int nXBlock[2] ;
      nXBlock[0] = min( nJ / m_nVoxNumPerBlock, m_nFracLin[0] - 1) ;
      if ( nJ % N_VOXBLOCK == 0  &&  nXBlock[0] > 0) {
         nXBlock[1] = nXBlock[0] - 1 ;
         ++ nXStop ;
      }
     // Voxel lungo Z
      int nZStop = 1 ;
      int nZBlock[2] ;
      nZBlock[0] = min( nI / m_nVoxNumPerBlock, m_nFracLin[2] - 1) ;
      if ( nI % N_VOXBLOCK == 0  &&  nZBlock[0] > 0) {
         nZBlock[1] = nZBlock[0] - 1 ;
         ++ nZStop ;
      }
     // Voxel lungo Y 
      int nVoxNumY = int( m_nNy[0] / m_nDexVoxRatio + ( m_nNy[0] % m_nDexVoxRatio == 0 ? 1 : 2)) ; 
      int nMinJ = Clamp( int( floor( ( ( dMin - 0.5 * m_dStep) / ( m_nDexVoxRatio * m_dStep) - EPS_SMALL))), 0, nVoxNumY - 2) ; 
      int nMaxJ = Clamp( int( floor( ( ( dMax + 0.5 * m_dStep) / ( m_nDexVoxRatio * m_dStep) + EPS_SMALL))), 0, nVoxNumY - 2) ; 
      int nMinYBlock = Clamp( nMinJ / int( m_nVoxNumPerBlock), 0, int( m_nFracLin[1] - 1)) ;
      int nMaxYBlock = min( int( m_nFracLin[1] - 1), nMaxJ / int( m_nVoxNumPerBlock)) ;
     // Assegno flag ai voxel
      for ( int tI = 0 ; tI < nZStop ; ++ tI) {
         for ( int tJ = 0 ; tJ < nXStop ; ++ tJ) {
            for ( int k = nMinYBlock ; k <= nMaxYBlock ; ++ k) {
               int nBlockNum = nZBlock[tI] * m_nFracLin[0] * m_nFracLin[1] + k * m_nFracLin[0] + nXBlock[tJ] ;
               m_BlockToUpdate[nBlockNum] = true ;
            }
         }
      } 
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::AddIntervals( int nGrid, int nI, int nJ,
                       double dMin, double dMax, const Vector3d& vtNMin, const Vector3d& vtNMax,
                       int nToolNum, bool bSkipSwap)
{
  // Controllo che il numero di griglia sia entro i limiti
   if ( nGrid < 0  ||  nGrid > 2)
      return false ;

  // Controllo che indici nI, nJ siano entro i limiti
   if ( nI < 0  &&  nI >= m_nNx[nGrid]  && 
        nJ < 0  &&  nJ >= m_nNy[nGrid])
      return false ;

  // Controllo che dMin < dMax
   Vector3d vtNmi = vtNMin ;
   Vector3d vtNma = vtNMax ;
   if ( dMin > dMax) {
      swap( dMin, dMax) ;
      swap( vtNmi, vtNma) ;
   }

  // Restringo minimo e massimo entro i limiti della mappa
   if ( dMin < m_dMinZ[nGrid]) {
      dMin = m_dMinZ[nGrid] ;
      if ( ! bSkipSwap)
         vtNmi = - Z_AX ;
   }
   else if ( dMin > m_dMaxZ[nGrid]) {
      dMin = m_dMaxZ[nGrid] ;
      if ( ! bSkipSwap)
         vtNmi = - Z_AX ;
   }
   if ( dMax < m_dMinZ[nGrid]) {
      dMax = m_dMinZ[nGrid] ;
      if ( ! bSkipSwap)
         vtNma = Z_AX ;
   }
   else if ( dMax > m_dMaxZ[nGrid]) {
      dMax = m_dMaxZ[nGrid] ;
      if ( ! bSkipSwap)
         vtNma = Z_AX ;
   }

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

  // Riporto le coordinate cicliche nell'ordine di partenza
   if ( !bSkipSwap  &&  nGrid == 1) {
      swap( vtNmi.x, vtNmi.z) ;
      swap( vtNmi.y, vtNmi.z) ;
      swap( vtNma.x, vtNma.z) ;
      swap( vtNma.y, vtNma.z) ;
   }
   else if ( !bSkipSwap  &&  nGrid == 2) {
      swap( vtNmi.y, vtNmi.z) ;
      swap( vtNmi.x, vtNmi.z) ;
      swap( vtNma.y, vtNma.z) ;
      swap( vtNma.x, vtNma.z) ;
   }

   // Calcolo nPos
   int nPos = nJ * m_nNx[nGrid] + nI ;
   vector<Data>& vDexel = m_Values[nGrid][nPos] ;

   bool bModified = false ;

  // Non esistono segmenti 
   if ( vDexel.empty()) {
   
      vDexel.emplace_back() ;
      vDexel.back().dMin = dMin ;
      vDexel.back().vtMinN = vtNmi ;
      vDexel.back().nToolMin = nToolNum ;
      vDexel.back().dMax = dMax ;
      vDexel.back().vtMaxN = vtNma ;
      vDexel.back().nToolMax = nToolNum ;
      
      m_OGrMgr.Reset() ;

      bModified = true ;
   }
  // Esiste almeno un segmento
   else {
     // Cerco l'ultimo intervallo a sinistra e l'ultimo intervallo a destra
     // di quello da aggiungere, che non interferiscono con quest'ultimo.
      auto itLastLeft = vDexel.end() ;
      auto itFirstRight = vDexel.end() ;
      for ( auto it = vDexel.begin() ; it != vDexel.end() ; ++ it) {
         if ( dMin > it->dMax + EPS_SMALL)
            itLastLeft = it ;
         if ( dMax < it->dMin - EPS_SMALL  &&  itFirstRight == vDexel.end())
            itFirstRight = it ;
      }
     // Esistono intervalli a sinistra.
      if ( itLastLeft != vDexel.end()) {
        // Intervallo successivo all'ultimo a sinistra
         auto itNextToLastLeft = itLastLeft ;
         ++ itNextToLastLeft ;
        // Il successivo non esiste.
         if ( itNextToLastLeft == vDexel.end()) {
            // Aggiungo il nuovo semgento
            vDexel.emplace_back() ;
            vDexel.back().dMin = dMin ;
            vDexel.back().dMax = dMax ;
            vDexel.back().vtMinN = vtNmi ;
            vDexel.back().vtMaxN = vtNma;
            vDexel.back().nToolMin = nToolNum ;
            vDexel.back().nToolMax = nToolNum ;
            //m_Values[nGrid][nPos].back().nCompo = ;
            bModified = true ;
         }
        // Il successivo esiste.
         else {
           // Il successivo è il primo a destra.
            if ( itNextToLastLeft == itFirstRight) {
              // Inserisco nuovo segmento-
               Data NewSegment ;
               NewSegment.dMin = dMin ;
               NewSegment.dMax = dMax ;
               NewSegment.vtMinN = vtNmi ;
               NewSegment.vtMaxN = vtNma ;
               NewSegment.nToolMin = nToolNum ;
               NewSegment.nToolMax = nToolNum ;
               //NewSegment.nCompo = ;
               vDexel.insert( itFirstRight, NewSegment) ;
               bModified = true ;
            }
            else {
              // Il successivo non esce a sinistra da quello da aggiungere.
               if ( itNextToLastLeft->dMin > dMin + EPS_SMALL) {
                  itNextToLastLeft->dMin = dMin ;
                  itNextToLastLeft->vtMinN = vtNmi ;
                  itNextToLastLeft->nToolMin = nToolNum ;
               }
              // Cerco l'ultimo segmento che interferisce con quello da aggiungere.
               auto itPrevToFirstRight = vDexel.end() ;
               for ( auto it = itNextToLastLeft ; it != itFirstRight ; ++ it) {
                  itPrevToFirstRight = it ;
               }
              // L'ultimo che interferisce non esce a destra da quello da aggiungere.
               if ( itPrevToFirstRight->dMax < dMax - EPS_SMALL) {
                  itNextToLastLeft->dMax = dMax ;
                  itNextToLastLeft->vtMaxN = vtNma ;
                  itNextToLastLeft->nToolMax = nToolNum ;
                  bModified = true ;
               }
               else {
                  itNextToLastLeft->dMax = itPrevToFirstRight->dMax ;
                  itNextToLastLeft->vtMaxN = itPrevToFirstRight->vtMaxN ;
                  itNextToLastLeft->nToolMax = nToolNum ;
                  bModified = true ;
               }
               auto itFirstToCancel = itNextToLastLeft ;
               ++ itFirstToCancel ;
               vDexel.erase( itFirstToCancel, itFirstRight) ;
            }
         }
      }
     // Non esistono neanche a destra.
      else if ( itFirstRight == m_Values[nGrid][nPos].end()) {
        // Il primo intervallo non sporge a sinistra
         if ( vDexel.begin()->dMin > dMin + EPS_SMALL) {
            vDexel.begin()->dMin = dMin ;
            vDexel.begin()->vtMinN = vtNmi ;
            vDexel.begin()->nToolMin = nToolNum ;
            bModified = true ;
         }
        // L'ultimo intervallo sporge a destra.
         if ( m_Values[nGrid][nPos].back().dMax > dMax + EPS_SMALL) {
           // Ci sono più segmenti, inglobo tutti nel primo.
            if ( vDexel.back().dMax > vDexel.begin()->dMax + EPS_SMALL) {
               vDexel.begin()->dMax = vDexel.back().dMax ;
               vDexel.begin()->vtMaxN = vDexel.back().vtMaxN ;
               vDexel.begin()->nToolMax = nToolNum ;
               bModified = true ;
            }
         }
        // L'ultimo intervallo non sporge a destra.
         else {
            vDexel.begin()->dMax = dMax ;
            vDexel.begin()->vtMaxN = vtNma ;
            vDexel.begin()->nToolMax = nToolNum ;
            bModified = true ;
         }
         vDexel.erase( vDexel.begin() + 1, vDexel.end()) ;
      }
     // A destra esistono.
      else {
        // Tutti i segmenti sono a destra di qullo da aggiungere.
         if ( itFirstRight == vDexel.begin()) {
           // Inserisco nuovo segmento-
            Data NewSegment ;
            NewSegment.dMin = dMin ;
            NewSegment.dMax = dMax ;
            NewSegment.vtMinN = vtNmi ;
            NewSegment.vtMaxN = vtNma ;
            NewSegment.nToolMin = nToolNum ;
            NewSegment.nToolMax = nToolNum ;
            vDexel.insert( vDexel.begin(), NewSegment) ;
            bModified = true ;
         }
         else {
           // Se il primo segmento non esce a sinistra da quello da aggiungere, cambio l'inizio.
            if ( vDexel.begin()->dMin > dMin + EPS_SMALL) {
               vDexel.begin()->dMin = dMin ;
               vDexel.begin()->vtMinN = vtNmi ;
               vDexel.begin()->nToolMin = nToolNum ;
               bModified = true ;
            }
           // Cerco l'ultimo segmento che interferisce con quello da aggiungere.
            auto itPrevToFirstRight = vDexel.begin() ;
            for ( auto it = m_Values[nGrid][nPos].begin() ; it != itFirstRight ; ++ it) {
               itPrevToFirstRight = it ;
            }
           // L'ultimo che interferisce non esce a destra da quello da aggiungere.
            if ( itPrevToFirstRight->dMax < dMax - EPS_SMALL) {
               vDexel.begin()->dMax = dMax ;
               vDexel.begin()->vtMaxN = vtNma ;
               vDexel.begin()->nToolMax = nToolNum ;
               bModified = true ;
            }
            else {
               vDexel.begin()->dMax = itPrevToFirstRight->dMax ;
               vDexel.begin()->vtMaxN = itPrevToFirstRight->vtMaxN ;
               vDexel.begin()->nToolMax = nToolNum ;
               bModified = true ;
            }
            auto itFirstToCancel = vDexel.begin() ;
            ++ itFirstToCancel ;
            vDexel.erase( itFirstToCancel, itFirstRight) ;
         }
      }  
   }

  // Se nessuna modifica, esco
   if ( ! bModified)
      return true ;

  // Imposto ricalcolo della grafica
   m_OGrMgr.Reset() ;
  // Imposto forma generica 
   m_nShape = GENERIC ;    
  // Imposto ricalcolo numero di componenti connesse
   m_nConnectedCompoCount = - 1 ;

  // Passo da indici di dexel a indici di voxel
   nI /= m_nDexVoxRatio ;
   nJ /= m_nDexVoxRatio ;

  // Determino quali blocchi sono stati modificati
   if ( nGrid == 0) {
     // Voxel lungo X
      int nXStop = 1 ;
      int nXBlock[2] ;
      nXBlock[0] = min( nI / m_nVoxNumPerBlock, m_nFracLin[0] - 1) ;
      if ( nI % N_VOXBLOCK == 0  &&  nXBlock[0] > 0) {
         nXBlock[1] = nXBlock[0] - 1 ;
         ++ nXStop ;
      }
     // Voxel lungo Y
      int nYStop = 1 ;
      int nYBlock[2] ;
      nYBlock[0] = min( nJ / m_nVoxNumPerBlock, m_nFracLin[1] - 1) ;
      if ( nJ % N_VOXBLOCK == 0  &&  nYBlock[0] > 0) {
         nYBlock[1] = nYBlock[0] - 1 ;
         ++ nYStop ;
      }
     // Voxel lungo Z
      int nVoxNumZ = int( m_nNy[1] / m_nDexVoxRatio + ( m_nNy[1] % m_nDexVoxRatio == 0 ? 1 : 2)) ;  
      int nMinK = Clamp( int( floor( ( ( dMin - 0.5 * m_dStep) / ( m_nDexVoxRatio * m_dStep) - EPS_SMALL))), 0, nVoxNumZ - 2) ;
      int nMaxK = Clamp( int( floor( ( ( dMax + 0.5 * m_dStep) / ( m_nDexVoxRatio * m_dStep) + EPS_SMALL))), 0, nVoxNumZ - 2) ;
      int nMinZBlock = ( m_nMapNum == 1 ? 0 : Clamp( nMinK / int( m_nVoxNumPerBlock), 0, int( m_nFracLin[2] - 1))) ;
      int nMaxZBlock = min( int( m_nFracLin[2] - 1), nMaxK / int( m_nVoxNumPerBlock)) ;
     // Assegno flag ai voxel
      for ( int tI = 0 ; tI < nXStop ; ++ tI) {
         for ( int tJ = 0 ; tJ < nYStop ; ++ tJ) {
            for ( int k = nMinZBlock ; k <= nMaxZBlock ; ++ k) {
               int nBlockNum = k * m_nFracLin[0] * m_nFracLin[1] + nYBlock[tJ] * m_nFracLin[0] + nXBlock[tI] ;
               m_BlockToUpdate[nBlockNum] = true ;
            }
         }
      }
   }

   else if ( nGrid == 1) {
     // Voxel lungo Y
      int nYStop = 1 ;
      int nYBlock[2] ;
      nYBlock[0] = min( nI / m_nVoxNumPerBlock, m_nFracLin[1] - 1) ;
      if ( nI % N_VOXBLOCK == 0  &&  nYBlock[0] > 0) {
         nYBlock[1] = nYBlock[0] - 1 ;
         ++ nYStop ;
      }
     // Voxel lungo Z
      int nZStop = 1 ;
      int nZBlock[2] ;
      nZBlock[0] = min( nJ / m_nVoxNumPerBlock, m_nFracLin[2] - 1) ;
      if ( nJ % N_VOXBLOCK == 0  &&  nZBlock[0] > 0) {
         nZBlock[1] = nZBlock[0] - 1 ;
         ++ nZStop ;
      }
     // Voxel lungo X
      int nVoxNumX = int( m_nNx[0] / m_nDexVoxRatio + ( m_nNx[0] % m_nDexVoxRatio == 0 ? 1 : 2)) ;
      int nMinI = Clamp( int( floor( ( ( dMin - 0.5 * m_dStep) / ( m_nDexVoxRatio * m_dStep) - EPS_SMALL))), 0, nVoxNumX - 2) ;
      int nMaxI = Clamp( int( floor( ( ( dMax + 0.5 * m_dStep) / ( m_nDexVoxRatio * m_dStep) + EPS_SMALL))), 0, nVoxNumX - 2) ;
      int nMinXBlock = Clamp( nMinI / int( m_nVoxNumPerBlock), 0, int( m_nFracLin[0] - 1)) ;
      int nMaxXBlock = min( int( m_nFracLin[0] - 1), nMaxI / int( m_nVoxNumPerBlock)) ;
     // Assegno flag ai voxel
      for ( int tI = 0 ; tI < nYStop ; ++ tI) {
         for ( int tJ = 0 ; tJ < nZStop ; ++ tJ) {
            for ( int k = nMinXBlock ; k <= nMaxXBlock ; ++ k) {
               int nBlockNum = nZBlock[tJ] * m_nFracLin[0] * m_nFracLin[1] + nYBlock[tI] * m_nFracLin[0] + k ;
               m_BlockToUpdate[nBlockNum] = true ;
            }
         }
      } 
   }

   else if ( nGrid == 2) {
     // Voxel lungo X
      int nXStop = 1 ;
      int nXBlock[2] ;
      nXBlock[0] = min( nJ / m_nVoxNumPerBlock, m_nFracLin[0] - 1) ;
      if ( nJ % N_VOXBLOCK == 0  &&  nXBlock[0] > 0) {
         nXBlock[1] = nXBlock[0] - 1 ;
         ++ nXStop ;
      }
     // Voxel lungo Z
      int nZStop = 1 ;
      int nZBlock[2] ;
      nZBlock[0] = min( nI / m_nVoxNumPerBlock, m_nFracLin[2] - 1) ;
      if ( nI % N_VOXBLOCK == 0  &&  nZBlock[0] > 0) {
         nZBlock[1] = nZBlock[0] - 1 ;
         ++ nZStop ;
      }
     // Voxel lungo Y 
      int nVoxNumY = int( m_nNy[0] / m_nDexVoxRatio + ( m_nNy[0] % m_nDexVoxRatio == 0 ? 1 : 2)) ; 
      int nMinJ = Clamp( int( floor( ( ( dMin - 0.5 * m_dStep) / ( m_nDexVoxRatio * m_dStep) - EPS_SMALL))), 0, nVoxNumY - 2) ; 
      int nMaxJ = Clamp( int( floor( ( ( dMax + 0.5 * m_dStep) / ( m_nDexVoxRatio * m_dStep) + EPS_SMALL))), 0, nVoxNumY - 2) ; 
      int nMinYBlock = Clamp( nMinJ / int( m_nVoxNumPerBlock), 0, int( m_nFracLin[1] - 1)) ;
      int nMaxYBlock = min( int( m_nFracLin[1] - 1), nMaxJ / int( m_nVoxNumPerBlock)) ;
     // Assegno flag ai voxel
      for ( int tI = 0 ; tI < nZStop ; ++ tI) {
         for ( int tJ = 0 ; tJ < nXStop ; ++ tJ) {
            for ( int k = nMinYBlock ; k <= nMaxYBlock ; ++ k) {
               int nBlockNum = nZBlock[tI] * m_nFracLin[0] * m_nFracLin[1] + k * m_nFracLin[0] + nXBlock[tJ] ;
               m_BlockToUpdate[nBlockNum] = true ;
            }
         }
      } 
   }

   return true ;
}

//----------------------------------------------------------------------------
bool 
VolZmap::AddSubIntervalInVoxel( VolZmap* VolZmapRef, int nGrid, int nI, int nJ, int nK,
                                double& dMin, double& dMax, Vector3d& vtMin, Vector3d& vtMax)
{
  // se non Tridex, esco
   if ( ! IsTriDexel())
      return true ;
  
  // Controllo che il numero di griglia sia entro i limiti
   if ( nGrid < 0  ||  nGrid > 2)
      return false ;

  // Controllo che indici nI, nJ siano entro i limiti
   if ( nI < 0  &&  nI >= m_nNx[nGrid]  &&
        nJ < 0  &&  nJ >= m_nNy[nGrid])
      return false ;

  // valori di default 
   vector<double> vTdMin = { -1., -1., -1., -1.} ;
   vector<double> vTdMax = { -1., -1., -1., -1.} ;
   VCT3DVECTOR vtNMin = { V_INVALID, V_INVALID, V_INVALID, V_INVALID} ;
   VCT3DVECTOR vtNMax = { V_INVALID, V_INVALID, V_INVALID, V_INVALID} ;

  // se esiste un precedente lungo x-locale
   if ( nI != 0) {
      int nPos = nJ * m_nNx[nGrid] + ( nI - 1) ; // recupero posizione dexel
     // cerco l'intervallo corretto sullo Zmap di riferimento
      for ( int it = 0 ; it < int( VolZmapRef->m_Values[nGrid][nPos].size()) ; ++ it) {
         if ( VolZmapRef->m_Values[nGrid][nPos][it].dMax + EPS_SMALL > nK * m_dStep  &&
              VolZmapRef->m_Values[nGrid][nPos][it].dMax - EPS_SMALL < ( nK + 1) * m_dStep) {
            vtNMax[0] = VolZmapRef->m_Values[nGrid][nPos][it].vtMaxN ;
            vTdMax[0] = VolZmapRef->m_Values[nGrid][nPos][it].dMax ;
         }
         if ( VolZmapRef->m_Values[nGrid][nPos][it].dMin + EPS_SMALL > nK * m_dStep  &&
              VolZmapRef->m_Values[nGrid][nPos][it].dMin - EPS_SMALL < ( nK + 1) * m_dStep) {
            vtNMin[0] = VolZmapRef->m_Values[nGrid][nPos][it].vtMinN ;
            vTdMin[0] = VolZmapRef->m_Values[nGrid][nPos][it].dMin ;
         }
      }
   }
  // se esiste il successivo lungo x-locale
   if ( nI != m_dMaxZ[( nGrid + 2) % 3]) {
      int nPos = nJ * m_nNx[nGrid] + ( nI + 1) ; // recupero posizione dexel
      for ( int it = 0 ; it < int( VolZmapRef->m_Values[nGrid][nPos].size()) ; ++ it) {
        // cerco l'intervallo corretto sullo Zmap di riferimento
         if ( VolZmapRef->m_Values[nGrid][nPos][it].dMax + EPS_SMALL > nK * m_dStep  &&
              VolZmapRef->m_Values[nGrid][nPos][it].dMax - EPS_SMALL < ( nK + 1) * m_dStep) {
            vtNMax[1] = VolZmapRef->m_Values[nGrid][nPos][it].vtMaxN ;
            vTdMax[1] = VolZmapRef->m_Values[nGrid][nPos][it].dMax ;
         }
         if ( VolZmapRef->m_Values[nGrid][nPos][it].dMin + EPS_SMALL > nK * m_dStep  &&
              VolZmapRef->m_Values[nGrid][nPos][it].dMin - EPS_SMALL < ( nK + 1) * m_dStep) {
            vtNMin[1] = VolZmapRef->m_Values[nGrid][nPos][it].vtMinN ;
            vTdMin[1] = VolZmapRef->m_Values[nGrid][nPos][it].dMin ;
         }
      }
   }
  // se esiste il precedente lungo y-locale
   if ( nJ != 0) {
      int nPos = ( nJ - 1) * m_nNx[nGrid] + nI ; // recupero posizione dexel
      for ( int it = 0 ; it < int( VolZmapRef->m_Values[nGrid][nPos].size()) ; ++ it) {
        // cerco l'intervallo corretto sullo Zmap di riferimento
         if ( VolZmapRef->m_Values[nGrid][nPos][it].dMax + EPS_SMALL > nK * m_dStep  &&
              VolZmapRef->m_Values[nGrid][nPos][it].dMax - EPS_SMALL < ( nK + 1) * m_dStep) {
            vtNMax[2] = VolZmapRef->m_Values[nGrid][nPos][it].vtMaxN ;
            vTdMax[2] = VolZmapRef->m_Values[nGrid][nPos][it].dMax ;
         }
         if ( VolZmapRef->m_Values[nGrid][nPos][it].dMin + EPS_SMALL > nK * m_dStep  &&
              VolZmapRef->m_Values[nGrid][nPos][it].dMin - EPS_SMALL < ( nK + 1) * m_dStep) {
            vtNMin[2] = VolZmapRef->m_Values[nGrid][nPos][it].vtMinN ;
            vTdMin[2] = VolZmapRef->m_Values[nGrid][nPos][it].dMin ;
         }
      }
   }
  // se esiste il successivo lungo y-locale
   if ( nJ != m_dMaxZ[( nGrid + 1) % 3]) {
      int nPos = ( nJ + 1) * m_nNx[nGrid] + nI ; // recupero posizione dexel
     // cerco l'intervallo corretto sullo Zmap di riferimento
      for ( int it = 0 ; it < int( VolZmapRef->m_Values[nGrid][nPos].size()) ; ++ it) {
         if ( VolZmapRef->m_Values[nGrid][nPos][it].dMax + EPS_SMALL > nK * m_dStep  &&
              VolZmapRef->m_Values[nGrid][nPos][it].dMax - EPS_SMALL < ( nK + 1) * m_dStep) {
            vtNMax[3] = VolZmapRef->m_Values[nGrid][nPos][it].vtMaxN ;
            vTdMax[3] = VolZmapRef->m_Values[nGrid][nPos][it].dMax ;
         }
         if ( VolZmapRef->m_Values[nGrid][nPos][it].dMin + EPS_SMALL > nK * m_dStep  &&
              VolZmapRef->m_Values[nGrid][nPos][it].dMin - EPS_SMALL < ( nK + 1) * m_dStep) {
            vtNMin[3] = VolZmapRef->m_Values[nGrid][nPos][it].vtMinN ;
            vTdMin[3] = VolZmapRef->m_Values[nGrid][nPos][it].dMin ;
         }
      }
   }

  // scelgo le normali che si discostano il meno possibile dalla normale della faccia del voxel

  // analisi dei massimi e delle normali ---------------------------
  // 1) angolo minimo tra la normale trovata la faccia del voxel
   double dMinAngle = ANG_FULL ;
  // 2) vettore di riferimento per la direzione della normale
   vtMax = ( nGrid == 0 ? m_MapFrame.VersZ() :
           ( nGrid == 1 ? m_MapFrame.VersX() : m_MapFrame.VersY())) ;
   Vector3d vtRef = vtMax ;
  // 3) determino il massimo per questo intervallo
   dMax = ( nK + 1) * m_dStep ;
   for ( int i = 0 ; i < 4 ; ++ i) { // scorro le normali
      if ( vtNMax[i].IsValid()) { // se normale trovata, quindi valida...
         double dCurrAngle ; // angolo corrente tra la normale della TriMesh e quella della faccia del voxel
         vtNMax[i].GetAngle( vtRef, dCurrAngle) ;
         if ( abs( dCurrAngle) < dMinAngle) { // se angolo minore del minimo trovato...
           // aggiorno i parametri
            dMinAngle = dCurrAngle ;
            vtMax = vtNMax[i] ;
            dMax = vTdMax[i] ;
         }
      }
   }

  // analisi dei minimi e delle normali ---------------------------
   dMinAngle = ANG_FULL ;
   vtRef.Invert() ;
   dMin = nK * m_dStep ;
   vtMin = vtRef ;
   for ( int i = 0 ; i < 4 ; ++ i) {
      if ( vtNMin[i].IsValid()) {
         double dCurrAngle ;
         vtNMin[i].GetAngle( vtRef, dCurrAngle) ;
         if ( abs( dCurrAngle) < dMinAngle) {
            dMinAngle = dCurrAngle  ;
            vtMin = vtNMin[i] ;
            dMin = vTdMin[i] ;
         }
      }
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::AddMissingIntervalsInVoxel( VolZmap* VolZmapRef, int nGrid, int nI, int nJ, double dZ, double dToler,
                                     Vector3d vtToolMin, Vector3d vtToolMax, int nToolNum)
{

  // se non Tridex, esco
   if ( ! IsTriDexel())
      return true ;
  
  // Controllo che il numero di griglia sia entro i limiti
   if ( nGrid < 0  ||  nGrid > 2)
      return false ;

  // Controllo che indici nI, nJ siano entro i limiti
   if ( nI < 0  &&  nI >= m_nNx[nGrid]  &&
        nJ < 0  &&  nJ >= m_nNy[nGrid])
      return false ;

  // passo da indici di dexel a indici di voxel
   nI /= m_nDexVoxRatio ;
   nJ /= m_nDexVoxRatio ;

  // numero di voxel nel dexel corrente
   int nVoxNum = int( m_nNy[(( nGrid+1) % 3)] / m_nDexVoxRatio +
                 ( m_nNy[(( nGrid+1) % 3)] % m_nDexVoxRatio == 0 ? 1 : 2)) ;

   int nK = 0 ;
   for ( int i = 0 ; i < nVoxNum ; ++ i) {
     // controllo se sono nel voxel corrente
      if ( i * m_dStep < dZ  &&  ( i + 1) * m_dStep > dZ) {
         nK = i ;
         break ;
      }
   }

  // ----------- griglia successiva -----------
   {
      int nMyGrid = ( nGrid + 1) % 3 ;
      int nMyI = nJ ;
      int nMyJ = nK ;
      int nMyK = nI ;
      double dMyMin ;
      double dMyMax ;
      Vector3d vtMyMin ;
      Vector3d vtMyMax ;
      AddSubIntervalInVoxel( VolZmapRef, nMyGrid, nMyI, nMyJ, nMyK, dMyMin, dMyMax, vtMyMin, vtMyMax) ;
      AddIntervals( nMyGrid, nMyI, nMyJ, dMyMin - EPS_SMALL, dMyMax + EPS_SMALL, vtMyMin, vtMyMax,
                    nToolNum, true) ;
   }

   // ----------- griglia precedente -----------
   {
      int nMyGrid = ( nGrid + 2) % 3 ;
      int nMyI = nK ;
      int nMyJ = nI ;
      int nMyK = nJ ;
      double dMyMin ;
      double dMyMax ;
      Vector3d vtMyMin ;
      Vector3d vtMyMax ;
      AddSubIntervalInVoxel( VolZmapRef, nMyGrid, nMyI, nMyJ, nMyK, dMyMin, dMyMax, vtMyMin, vtMyMax) ;
      AddIntervals( nMyGrid, nMyI, nMyJ, dMyMin - EPS_SMALL, dMyMax + EPS_SMALL, vtMyMin, vtMyMax,
                    nToolNum, true) ;
   }

   return true ;

}

// ------------------------- LAVORAZIONI --------------------------------------------------------------------------------------

//----------------------------------------------------------------------------       
bool
VolZmap::MillingStep( int nCurrTool,
                      const Point3d& ptPs, const Vector3d& vtDs,
                      const Point3d& ptPe, const Vector3d& vtDe) 
{
   return MillingStep( nCurrTool, ptPs, vtDs, Vector3d(), ptPe, vtDe, Vector3d()) ;
}

//----------------------------------------------------------------------------       
bool
VolZmap::MillingStep( int nCurrTool,
                      const Point3d& ptPs, const Vector3d& vtDs, const Vector3d& vtAs,
                      const Point3d& ptPe, const Vector3d& vtDe, const Vector3d& vtAe)
{
  // Controllo utensile
   if ( nCurrTool < 0  ||  nCurrTool >= int( m_vTool.size()))
      return false ;
   m_nCurrTool = nCurrTool ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

  // Se non è definito l'utensile, non posso fare alcunchè
   if ( CurrTool.GetType() == Tool::UNDEF)
      return false ;
  // Controllo definizione vettori direzione
   if ( vtDs.IsSmall()  ||  vtDe.IsSmall())
      return false ;
  // Controllo sull'effettiva esistenza del movimento
   if ( AreSamePointApprox( ptPs, ptPe)  &&  AreSameVectorApprox( vtDs, vtDe))
      return true ;
  // Se utensile chisel o mortise, controllo definizione vettori ausiliari
   if ( CurrTool.GetType() == Tool::MORTISER  ||  CurrTool.GetType() == Tool::CHISEL) {
      if ( vtAs.IsSmall()  ||  vtAe.IsSmall())
         return false ;
   }

  // Punti e vettori descriventi il moto nel sistema intrinseco dello Zmap
   Point3d ptPLs = GetToLoc( ptPs, m_MapFrame) ;
   Point3d ptPLe = GetToLoc( ptPe, m_MapFrame) ;
   Vector3d vtDLs = GetToLoc( vtDs, m_MapFrame) ;
   vtDLs.Normalize() ;
   Vector3d vtDLe = GetToLoc( vtDe, m_MapFrame) ;
   vtDLe.Normalize() ;
   Vector3d vtALs = GetToLoc( vtAs, m_MapFrame) ;
   vtALs.Normalize() ;
   Vector3d vtALe = GetToLoc( vtAe, m_MapFrame) ;
   vtALe.Normalize() ;

   //static PerformanceCounter Counter ;
   //{
   //string sOut = "Draw=" + ToString( Counter.Stop(), 3) ;
   //LOG_INFO( GetEGkLogger(), sOut.c_str())
   //Counter.Start() ;
   //}

  // Se pura traslazione
   bool bOk ;
   if ( AreSameVectorApprox( vtDLs, vtDLe))
      bOk = MillingTranslationStep( ptPLs, ptPLe, vtDLs, vtALs) ;
   else
      bOk = MillingGeneralMotionStep( ptPLs, vtDLs, vtALs, ptPLe, vtDLe, vtALe) ;

   //{
   //string sOut = "Calc=" + ToString( Counter.Stop(), 3) ;
   //LOG_INFO( GetEGkLogger(), sOut.c_str())
   //Counter.Start() ;
   //}

   return bOk ;
}

//----------------------------------------------------------------------------
bool
VolZmap::MillingGeneralMotionStep( const Point3d& ptPs, const Vector3d& vtDs, const Vector3d& vtAs,
                                   const Point3d& ptPe, const Vector3d& vtDe, const Vector3d& vtAe)
{
  // Divido il movimento in tratti con direzione utensile costante
   const double ANG_STEP = 0.04 ;
   double dAngDeg ; vtDs.GetAngle( vtDe, dAngDeg) ;
   int nStepCnt = int( abs( dAngDeg) / ANG_STEP) + 1 ;
   bool bOk = true ;
   Point3d ptSt = ptPs ;
   for ( int i = 0 ; i <= nStepCnt && bOk ; ++ i) {
      double dPosCoeff, dDirCoeff ;
      if ( i < nStepCnt) {
         dPosCoeff = ( i + 0.5) / nStepCnt ;
         dDirCoeff = double( i) / nStepCnt ;
      }
      else {
         dPosCoeff = 1 ;
         dDirCoeff = 1 ;
      }
      Point3d ptEn = Media( ptPs, ptPe, dPosCoeff) ;
      Vector3d vtD = Media( vtDs, vtDe, dDirCoeff) ; vtD.Normalize() ;
      Vector3d vtA = Media( vtAs, vtAe, dDirCoeff) ; vtA.Normalize() ;
      bOk = bOk && MillingTranslationStep( ptSt, ptEn, vtD, vtA) ;
     // aggiorno prossimo inizio
      ptSt = ptEn ;
   }
   return bOk ;
}

//----------------------------------------------------------------------------
bool
VolZmap::MillingTranslationStep( const Point3d& ptPs, const Point3d& ptPe, const Vector3d& vtD, const Vector3d& vtA)
{
   Point3d ptLs[N_MAPS] ;
   Point3d ptLe[N_MAPS] ;
   Vector3d vtLs[N_MAPS] ;
   Vector3d vtALs[N_MAPS] ;
   InitializePointsAndVectors( ptPs, ptPe, vtD, vtA, ptLs, ptLe, vtLs, vtALs) ;
  // Ciclo sulle mappe (scommentare solo per DEBUG)
   //{
   //   bool bOk = true ;
   //   for ( int i = 0 ; i < m_nMapNum ; ++ i) {
   //      bOk = SelectMotion( i, ptLs[i], ptLe[i], vtLs[i], vtALs[i]) && bOk ; 
   //   }
   //   return true ;
   //}
  // Ciclo sulle mappe
   vector< future<bool>> vRes ;
   vRes.resize( m_nMapNum) ;
   for ( int i = 0 ; i < m_nMapNum ; ++ i) {
      vRes[i] = async( launch::async, &VolZmap::SelectMotion, this, i, cref( ptLs[i]), cref( ptLe[i]), cref( vtLs[i]), cref( vtALs[i])) ;
   }
   bool bOk = true ;
   int nTerminated = 0 ;
   while ( nTerminated < m_nMapNum) {
      for ( int i = 0 ; i < m_nMapNum ; ++ i) {
         if ( vRes[i].valid()  &&  vRes[i].wait_for( chrono::nanoseconds{ 1}) == future_status::ready) {
            bOk = vRes[i].get() && bOk ;
            ++ nTerminated ;
         }
      }
   }
   return bOk ;
}

//----------------------------------------------------------------------------
bool
VolZmap::InitializePointsAndVectors( const Point3d& ptPs, const Point3d& ptPe, const Vector3d& vtDs, const Vector3d& vtAs,
                                     Point3d ptLs[3], Point3d ptLe[3], Vector3d vtLs[3], Vector3d vtALs[3])
{
  // Punti posizione
   ptLs[0] = ptPs ;
   ptLe[0] = ptPe ;
   if ( m_nMapNum > 1) {
      ptLs[1].x = ptLs[0].y ;  ptLs[1].y = ptLs[0].z ;  ptLs[1].z = ptLs[0].x ;
      ptLs[2].x = ptLs[0].z ;  ptLs[2].y = ptLs[0].x ;  ptLs[2].z = ptLs[0].y ;
      ptLe[1].x = ptLe[0].y ;  ptLe[1].y = ptLe[0].z ;  ptLe[1].z = ptLe[0].x ;
      ptLe[2].x = ptLe[0].z ;  ptLe[2].y = ptLe[0].x ;  ptLe[2].z = ptLe[0].y ;
   }

  // Vettori asse utensile
   vtLs[0] = vtDs ;
   if ( m_nMapNum > 1) {
      vtLs[1].x = vtLs[0].y ;  vtLs[1].y = vtLs[0].z ;  vtLs[1].z = vtLs[0].x ;
      vtLs[2].x = vtLs[0].z ;  vtLs[2].y = vtLs[0].x ;  vtLs[2].z = vtLs[0].y ;
   }

  // Vettori ausiliari 
   vtALs[0] = vtAs ;
   if ( m_nMapNum > 1) {
      vtALs[1].x = vtALs[0].y ;  vtALs[1].y = vtALs[0].z ;  vtALs[1].z = vtALs[0].x ;
      vtALs[2].x = vtALs[0].z ;  vtALs[2].y = vtALs[0].x ;  vtALs[2].z = vtALs[0].y ;
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::SelectMotion( int nGrid, const Point3d& ptLs, const Point3d& ptLe, const Vector3d& vtL, const Vector3d& vtAL)
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

  // Versore utensile parallelo all'asse Z (coincide con spillone)
   if ( vtL.SqLenXY() < EPS_ZERO * EPS_ZERO) {

      Vector3d vtMove = ptLe - ptLs ;

     // Foratura
      if ( vtMove.SqLenXY() < EPS_SMALL * EPS_SMALL) {
         switch ( CurrTool.GetType()) {
         case Tool::GEN :
            return GenTool_ZDrilling( nGrid, ptLs, ptLe, vtL) ;
         case Tool::CYLMILL :
         case Tool::BALLMILL :
            return CylBall_ZDrilling( nGrid, ptLs, ptLe, vtL) ;
         case Tool::CONEMILL :
            return Conus_ZDrilling( nGrid, ptLs, ptLe, vtL) ;
         case Tool::MORTISER :
            return Mrt_ZDrilling( nGrid, ptLs, ptLe, vtL, vtAL) ;
         case Tool::CHISEL :
            return Chs_ZDrilling( nGrid, ptLs, ptLe, vtL, vtAL) ;
         }
      }

     // Fresatura con vettore movimento perpendicolare all'utensile
      else if ( abs( vtMove.z) < EPS_SMALL) {
         switch ( CurrTool.GetType()) {
         case Tool::GEN :
            return GenTool_ZMilling( nGrid, ptLs, ptLe, vtL) ;
         case Tool::CYLMILL :
         case Tool::BALLMILL :
            return CylBall_ZPerp( nGrid, ptLs, ptLe, vtL) ;
         case Tool::CONEMILL :
            return Conus_ZPerp( nGrid, ptLs, ptLe, vtL) ;
         case Tool::MORTISER :
            return Mrt_Milling( nGrid, ptLs, ptLe, vtL, vtAL) ;
         case Tool::CHISEL :
            return Chs_Milling( nGrid, ptLs, ptLe, vtL, vtAL) ;
         case Tool::ADDITIVE :
            return AddingMotion( nGrid, ptLs, ptLe, vtL) ;
         }
      }

     // Fresatura con vettore movimento generico rispetto all'utensile
      else {
         switch ( CurrTool.GetType()) {
         case Tool::GEN :
            return GenTool_ZMilling( nGrid, ptLs, ptLe, vtL) ;
         case Tool::CYLMILL :
         case Tool::BALLMILL :
            return CylBall_ZMilling( nGrid, ptLs, ptLe, vtL) ;
         case Tool::CONEMILL :
            return Conus_ZMilling( nGrid, ptLs, ptLe, vtL) ;
         case Tool::MORTISER :
            return Mrt_Milling( nGrid, ptLs, ptLe, vtL, vtAL) ;
         case Tool::ADDITIVE :
            return AddingMotion( nGrid, ptLs, ptLe, vtL) ;
         }
      }
   }

  // Versore utensile nel piano XY (perpendicolare allo spillone)
   else if ( abs( vtL.z) < EPS_ZERO) {

      Vector3d vtMove = ptLe - ptLs ;
      Vector3d vtMLong = ( vtMove * vtL) * vtL ;
      Vector3d vtMOrt = vtMove - vtMLong ;

      double dSqLLong = vtMLong.SqLen() ;
      double dSqLOrt = vtMOrt.SqLen() ;

     // Foratura
      if ( dSqLOrt < EPS_SMALL * EPS_SMALL) {
         switch ( CurrTool.GetType()) {
         case Tool::GEN :
            return GenTool_Drilling( nGrid, ptLs, ptLe, vtL) ;
         case Tool::CYLMILL :
         case Tool::BALLMILL :
            return CylBall_XYDrilling( nGrid, ptLs, ptLe, vtL) ;
         case Tool::CONEMILL :
            return Conus_XYDrilling( nGrid, ptLs, ptLe, vtL) ;
         case Tool::MORTISER :
            return Mrt_Drilling( nGrid, ptLs, ptLe, vtL, vtAL) ;
         case Tool::CHISEL :
            return Chs_Drilling( nGrid, ptLs, ptLe, vtL, vtAL) ;
         }
      }
     // Fresatura con vettore movimento perpendicolare all'utensile
      else if ( dSqLLong < EPS_SMALL * EPS_SMALL) {
         switch ( CurrTool.GetType()) {
         case Tool::GEN :
            return GenTool_Milling( nGrid, ptLs, ptLe, vtL) ;
         case Tool::CYLMILL :
         case Tool::BALLMILL :
            return CylBall_XYPerp( nGrid, ptLs, ptLe, vtL) ;
         case Tool::CONEMILL :
           // Usiamo la generica per via dell'intsabilità di Conus_XYPerp
           //return  Conus_XYPerp( i, ptLs[i], ptLe[i], vtLs[i]) ;
            return Conus_Milling( nGrid, ptLs, ptLe, vtL) ;
         case Tool::MORTISER :
            return Mrt_Milling( nGrid, ptLs, ptLe, vtL, vtAL) ;
         case Tool::CHISEL :
            return Chs_Milling( nGrid, ptLs, ptLe, vtL, vtAL) ;
         case Tool::ADDITIVE :
            return AddingMotion( nGrid, ptLs, ptLe, vtL) ;
         }
      }
     // Fresatura con vettore movimento generico rispetto all'utensile
      else {
         switch ( CurrTool.GetType()) {
         case Tool::GEN :
            return GenTool_Milling( nGrid, ptLs, ptLe, vtL);
         case Tool::CYLMILL :
         case Tool::BALLMILL :
            return CylBall_XYMilling( nGrid, ptLs, ptLe, vtL) ;
         case Tool::CONEMILL :
            return Conus_XYMilling( nGrid, ptLs, ptLe, vtL) ;
         case Tool::MORTISER :
            return Mrt_Milling( nGrid, ptLs, ptLe, vtL, vtAL) ;
         case Tool::ADDITIVE :
            return AddingMotion( nGrid, ptLs, ptLe, vtL) ;
         }
      }
   }

  // Versore utensile con direzione generica
   else {
      Vector3d vtMove = ptLe - ptLs ;
      Vector3d vtMLong = ( vtMove * vtL) * vtL ;
      Vector3d vtMOrt = vtMove - vtMLong ;

      double dSqLLong = vtMLong.SqLen() ;
      double dSqLOrt = vtMOrt.SqLen() ;

      // Foratura
      if ( dSqLOrt < EPS_SMALL * EPS_SMALL) {
         switch ( CurrTool.GetType()) {
         case Tool::GEN :
            return GenTool_Drilling( nGrid, ptLs, ptLe, vtL) ;
         case Tool::CYLMILL :
         case Tool::BALLMILL :
            return CylBall_Drilling( nGrid, ptLs, ptLe, vtL) ;
         case Tool::CONEMILL :
            return Conus_Drilling( nGrid, ptLs, ptLe, vtL) ;
         case Tool::MORTISER :
            return Mrt_Drilling( nGrid, ptLs, ptLe, vtL, vtAL) ;
         case Tool::CHISEL :
            return Chs_Drilling( nGrid, ptLs, ptLe, vtL, vtAL) ;
         }
      }
     // Fresatura con vettore movimento generico rispetto all'utensile
      else {
         switch ( CurrTool.GetType()) {
         case Tool::GEN :
            return GenTool_Milling( nGrid, ptLs, ptLe, vtL) ;
         case Tool::CYLMILL :
         case Tool::BALLMILL :
            return CylBall_Milling( nGrid, ptLs, ptLe, vtL) ;
         case Tool::CONEMILL :
            return Conus_Milling( nGrid, ptLs, ptLe, vtL) ;
         case Tool::MORTISER :
            return Mrt_Milling( nGrid, ptLs, ptLe, vtL, vtAL) ;
         case Tool::CHISEL :
           // ammesso solo movimento perpendicolare all'asse utensile
            if ( dSqLLong < EPS_SMALL * EPS_SMALL)
               return Chs_Milling( nGrid, ptLs, ptLe, vtL, vtAL) ;
            break ;
         case Tool::ADDITIVE :
            return AddingMotion( nGrid, ptLs, ptLe, vtL) ;
         }
      }
   }
   return false ;
}


// ---------- VERSORE UTENSILE DERETTO COME Z --------------------------------


// ---------- Cilindro e sfera -----------------------------------------------
//----------------------------------------------------------------------------     
bool
VolZmap::CylBall_ZDrilling( int nGrid, const Point3d & ptS, const Point3d & ptE, const Vector3d & vtToolDir)
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

  // Verifica sull'interferenza utensile Zmap
   int nStartI, nStartJ, nEndI, nEndJ ;
   if ( ! TestToolBBox( nGrid, ptS, ptE, vtToolDir, nStartI, nStartJ, nEndI, nEndJ))
      return true ;

  // Proiezione dei punti sul piano
   Point3d ptSxy( ptS.x, ptS.y, 0) ;

  // Parametri geometrici dell'utensile
   double dStemHeigth = CurrTool.GetHeigth() - CurrTool.GetTipHeigth() ;
   double dSqRad = CurrTool.GetRadius() * CurrTool.GetRadius() ;

  // Punte del gambo
   Point3d ptTStemS = ptS - vtToolDir * dStemHeigth ;
   Point3d ptTStemE = ptE - vtToolDir * dStemHeigth ;

  // Quote estreme del gambo
   double dMinStemZ = min( min( ptS.z, ptTStemS.z), min( ptE.z, ptTStemE.z)) ;
   double dMaxStemZ = max( max( ptS.z, ptTStemS.z), max( ptE.z, ptTStemE.z)) ;

   for ( int i = nStartI ; i <= nEndI ; ++ i) {
      for( 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 - ptSxy ;
         double dSqLen = vtC.SqLen() ;

        // Se il punto si trova dentro il cerchio taglio
         if ( dSqLen < dSqRad - 2 * CurrTool.GetRadius() * EPS_SMALL) {
           // utensile cilindrico
            if ( CurrTool.GetType() == Tool::CYLMILL)
               SubtractIntervals( nGrid, i, j, dMinStemZ, dMaxStemZ, Z_AX, -Z_AX, CurrTool.GetToolNum()) ;
           // utensile sferico
            else if ( CurrTool.GetType() == Tool::BALLMILL) {
               double dH = sqrt( dSqRad - dSqLen) ;
               if ( vtToolDir.z > 0) {
                   Vector3d vtNorm = Point3d( ptSxy.x, ptSxy.y, dMinStemZ) - Point3d( dX, dY, dMinStemZ - dH) ;
                   vtNorm.Normalize() ;
                   SubtractIntervals( nGrid, i, j, dMinStemZ - dH, dMaxStemZ, vtNorm, -vtToolDir, CurrTool.GetToolNum()) ;
               }
               else {
                   Vector3d vtNorm = Point3d( ptSxy.x, ptSxy.y, dMinStemZ) - Point3d( dX, dY, dMinStemZ + dH) ;
                   vtNorm.Normalize() ;
                   SubtractIntervals( nGrid, i, j, dMinStemZ, dMaxStemZ + dH, -vtToolDir, vtNorm, CurrTool.GetToolNum()) ;
               }
            }
         }
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::CylBall_ZPerp( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir)
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

  // Verifica sull'interferenza utensile Zmap
   int nStartI, nStartJ, nEndI, nEndJ ;
   if ( ! TestToolBBox( nGrid, ptS, ptE, vtToolDir, nStartI, nStartJ, nEndI, nEndJ))
      return true ;

  // Parametri geometrici dell'utensile
   double dStemHeigth = CurrTool.GetHeigth() - CurrTool.GetTipHeigth() ;
   double dSafeRad = CurrTool.GetRadius() - EPS_SMALL ;
   double dSafeSqRad = dSafeRad * dSafeRad ;

  // Punte del gambo
   Point3d ptTStemS = ptS - vtToolDir * dStemHeigth ;
   Point3d ptTStemE = ptE - vtToolDir * dStemHeigth ;

  // Quote estreme del gambo
   double dMinStemZ = min( ptS.z, ptTStemS.z) ;
   double dMaxStemZ = max( ptS.z, ptTStemS.z) ;

  // Vettore movimento e sua lunghezza
   Vector3d vtMove = ptE - ptS ;
   double dLen = vtMove.LenXY() ;

  // Definizione di un sistema di riferimento ad hoc
   Point3d ptSxy( ptS.x, ptS.y, 0) ;
   Vector3d vtV1 = vtMove ;
   vtV1.Normalize() ;   // se |vtMove| < EPS è un buco con dz = 0
   Vector3d vtV2 = vtV1 ;
   vtV2.Rotate( Z_AX, 0, 1) ;

   for ( int i = nStartI ; i <= nEndI ; ++ i) {
      for( 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 - ptSxy ;

         double dP1 = vtC * vtV1 ;           double dP2 = vtC * vtV2 ;

        // Utensile cilindrico
         if ( CurrTool.GetType() == Tool::CYLMILL) { 
           // Se il punto cade nella zona di interesse taglio 
            if ( ( dP1 * dP1 + dP2 * dP2 < dSafeSqRad) || 
                 ( ( dP1 - dLen) * ( dP1 - dLen) + dP2 * dP2) < dSafeSqRad ||
                 ( dP1 > 0 && dP1 < dLen && abs( dP2) < dSafeRad))
               
                 SubtractIntervals( nGrid, i, j, dMinStemZ, dMaxStemZ, Z_AX, -Z_AX, CurrTool.GetToolNum()) ;
         }

        // Utensile sferico
         else if ( CurrTool.GetType() == Tool::BALLMILL) {

            if ( abs( dP2) < dSafeRad) {

               if ( dP1 < 0) {
                 
                  if ( dP1 * dP1 + dP2 * dP2 < dSafeSqRad) {

                     double dH = sqrt( dSafeSqRad - dP1 * dP1 - dP2 * dP2) ;

                     if ( vtToolDir.z > 0) {
                        Point3d ptInt( dX, dY, ptTStemS.z - dH) ;
                        Vector3d vtN = ptTStemS - ptInt ;
                        vtN.Normalize() ;
                        SubtractIntervals( nGrid, i, j, ptTStemS.z - dH, ptS.z, vtN, -Z_AX, CurrTool.GetToolNum()) ;
                     }

                     else {
                        Point3d ptInt( dX, dY, ptTStemS.z + dH) ;
                        Vector3d vtN = ptTStemS - ptInt ;
                        vtN.Normalize() ;
                        SubtractIntervals( nGrid, i, j, ptS.z, ptTStemS.z + dH, Z_AX, vtN, CurrTool.GetToolNum()) ;
                     }
                  }
               }

               else if ( dP1 < dLen) {

                  double dH = sqrt( dSafeSqRad - dP2 * dP2) ;
                         
                  if ( vtToolDir.z > 0) {
                     Point3d ptInt( dX, dY, ptTStemS.z - dH) ;
                     Vector3d vtN = - ( ptInt - ptTStemS) + ( ptInt - ptTStemS) * vtV1 * vtV1 ;
                     vtN.Normalize() ;
                     SubtractIntervals( nGrid, i, j, ptTStemS.z - dH, ptS.z, vtN, -Z_AX, CurrTool.GetToolNum()) ;
                  }

                  else {
                     Point3d ptInt( dX, dY, ptTStemS.z + dH) ;
                     Vector3d vtN = - ( ptInt - ptTStemS) + ( ptInt - ptTStemS) * vtV1 * vtV1 ;
                     vtN.Normalize() ;
                     SubtractIntervals( nGrid, i, j, ptS.z, ptTStemS.z + dH, Z_AX, vtN, CurrTool.GetToolNum()) ;
                  }
               }

               else {
                  if ( ( dP1 - dLen) * ( dP1 - dLen) + dP2 * dP2 < dSafeSqRad) {
                  
                     double dH = sqrt( dSafeSqRad - ( dP1 - dLen) * ( dP1 - dLen) - dP2 * dP2) ;

                     if ( vtToolDir.z > 0) {
                        Point3d ptInt( dX, dY, ptTStemE.z - dH) ;
                        Vector3d vtN = ptTStemE - ptInt ;
                        vtN.Normalize() ;
                        SubtractIntervals( nGrid, i, j, ptTStemE.z - dH, ptS.z, vtN, -Z_AX, CurrTool.GetToolNum()) ;
                     }

                     else  {
                        Point3d ptInt( dX, dY, ptTStemE.z + dH) ;
                        Vector3d vtN = ptTStemE - ptInt ;
                        vtN.Normalize() ;
                        SubtractIntervals( nGrid, i, j, ptS.z, ptTStemE.z + dH, Z_AX, vtN, CurrTool.GetToolNum()) ;
                     }
                  }
               }
            }
         }
      }
   }
   return true ;
}

//----------------------------------------------------------------------------     
bool
VolZmap::CylBall_ZMilling( int nGrid, const Point3d & ptS, const Point3d & ptE, const Vector3d & vtToolDir)
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

  // Verifica sull'interferenza utensile Zmap
   int nStartI, nStartJ, nEndI, nEndJ ;
   if ( ! TestToolBBox( nGrid, ptS, ptE, vtToolDir, nStartI, nStartJ, nEndI, nEndJ))
      return true ;

  // Parametri geometrici dell'utensile
   double dStemHeigth = CurrTool.GetHeigth() - CurrTool.GetTipHeigth() ;
   double dSafeRadius = CurrTool.GetRadius() - EPS_SMALL ;
   double dSqRad = CurrTool.GetRadius() * CurrTool.GetRadius() ;
   double dSafeSqRad = dSqRad - 2 * CurrTool.GetRadius() * EPS_SMALL ;

  // Studio delle simmetrie
   Point3d ptI = ( ptS.z < ptE.z ? ptS : ptE) ;
   Point3d ptF = ( ptS.z < ptE.z ? ptE : ptS) ;
   Point3d ptIT = ptI - vtToolDir * dStemHeigth ;
   Point3d ptFT = ptF - vtToolDir * dStemHeigth ;
   Point3d ptIxy( ptI.x, ptI.y, 0) ;

   Point3d ptIUp( ptI.x, ptI.y, max( ptI.z, ptIT.z)) ;
   Point3d ptIDw( ptI.x, ptI.y, min( ptI.z, ptIT.z)) ;
  
  // Quote iniziali e finali massime e
  // minime del gambo dell'utensile e DeltaZ
   double dZMaxI = max( ptI.z, ptIT.z) ;
   double dZMaxF = max( ptF.z, ptFT.z) ;
   double dZMinI = dZMaxI - dStemHeigth ;
   double dDeltaZ = dZMaxF - dZMaxI ;

  // Vettori caratterizzanti il moto 
   Vector3d vtMove = ptF - ptI ;
   Vector3d vtMoveXY( vtMove.x, vtMove.y, 0) ;
   double dLenXY = vtMoveXY.LenXY() ;
   vtMove.Normalize() ;

  // Definizione di un sistema di riferimento ad hoc
   Vector3d vtV1, vtV2 ;

  // Se la lunghezza è troppo piccola lo allungo
   if ( dLenXY < EPS_SMALL) 
      vtV1 = ( 1 / dLenXY) * vtMoveXY ;
   else
      vtV1 = vtMoveXY ;

  // Normalizzo vtV1
   vtV1.Normalize() ;
  // Definisco vtV2
   vtV2 = vtV1 ;
   vtV2.Rotate( Z_AX, 0, 1) ;

   double dMin, dMax ;

   for ( int i = nStartI ; i <= nEndI ; ++ i) {
      for ( 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 - ptIxy ;

         double dX1 = vtC * vtV1 ;
         double dX2 = vtC * vtV2 ;

        // Se il punto appartiene alla proiezione del volume spazzato valuto massimo e minimo
         if ( ( dX1 > 0 && dX1 < dLenXY && abs( dX2) < dSafeRadius)    ||
              ( dX1 - dLenXY) * ( dX1 - dLenXY) + dX2 * dX2 < dSafeSqRad || 
                dX1 * dX1 + dX2 * dX2 < dSafeSqRad) {

            double dX1_0 = sqrt( dSqRad - dX2 * dX2) ; 
            Vector3d vtNmin, vtNmax ;

           // Massimo
            if ( ( dX1 - dLenXY) * ( dX1 - dLenXY) + dX2 * dX2 < dSqRad) {
               vtNmax = - Z_AX ;
               dMax = dZMaxF ;
            }

            else { 
               //vtNmax = - ( dX2 / dX1_0) * vtMove ^ vtV1 - vtMove ^ vtV2 ; 
               dMax = dZMaxI + dDeltaZ * ( dX1 + dX1_0) / dLenXY ;

               Vector3d vtCirc = dX1_0 * vtV1 - dX2 * vtV2 ; // Punta verso il centro
               Vector3d vtTan( - vtCirc.y, vtCirc.x, 0) ;    // Tangente alla circonferenza
               Vector3d vtCross = vtTan ^ vtMove ; 
               
               vtNmax = ( vtCross * vtCirc > - EPS_ZERO ? vtCross : - vtCross) ;
               vtNmax.Normalize() ;
            }
           // Minimo
            if ( dX1 * dX1 + dX2 * dX2 < dSqRad) {
               vtNmin = Z_AX ;
               dMin = dZMinI ;
            }

            else {
               //vtNmin = - ( dX2 / dX1_0) * vtMove ^ vtV1 + vtMove ^ vtV2 ;
               dMin = dZMinI + dDeltaZ * ( dX1 - dX1_0) / dLenXY ;

               Vector3d vtCirc = - dX1_0 * vtV1 - dX2 * vtV2 ; // Punta verso il centro
               Vector3d vtTan( - vtCirc.y, vtCirc.x, 0) ;    // Tangente alla circonferenza
               Vector3d vtCross = vtTan ^ vtMove ; 

               vtNmin = ( vtCross * vtCirc > - EPS_ZERO ? vtCross : - vtCross) ;
               vtNmin.Normalize() ;
            }

            SubtractIntervals( nGrid, i, j, dMin, dMax, vtNmin, vtNmax, CurrTool.GetToolNum()) ;
         }
      }
   }

   if ( CurrTool.GetType() == Tool::BALLMILL)

      CompBall_Milling( nGrid, ptIT, ptFT, CurrTool.GetRadius(), CurrTool.GetToolNum()) ;

   return true ;
}


// ---------- Coni -----------------------------------------------------------
//----------------------------------------------------------------------------
bool
VolZmap::Conus_ZDrilling( int nGrid, const Point3d & ptS, const Point3d & ptE, const Vector3d & vtToolDir)
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

  // Verifica sull'interferenza utensile Zmap
   int nStartI, nStartJ, nEndI, nEndJ ;
   if ( ! TestToolBBox( nGrid, ptS, ptE, vtToolDir, nStartI, nStartJ, nEndI, nEndJ))
      return true ;

  // Parametri geometrici dell'utensile
   double dStemHeigth = CurrTool.GetHeigth() - CurrTool.GetTipHeigth() ;
   double dMinRad = min( CurrTool.GetRadius(), CurrTool.GetTipRadius()) ;
   double dMaxRad = max( CurrTool.GetRadius(), CurrTool.GetTipRadius()) ;
   double dDeltaRad = dMaxRad - dMinRad ;
   double dSqMinRad = dMinRad * dMinRad ;
   double dSqMaxRad = dMaxRad * dMaxRad ;

  // Proiezione delle posizioni sul piano
   Point3d ptO( ptS.x, ptS.y, 0) ;

  // Quote massime e minime dell'utensile durante il moto
   double dZMax = max( max( ptS.z, ptS.z - vtToolDir.z * CurrTool.GetHeigth()),
                       max( ptE.z, ptE.z - vtToolDir.z * CurrTool.GetHeigth())) ;
   double dZMin = min( min( ptS.z, ptS.z - vtToolDir.z * CurrTool.GetHeigth()), 
                       min( ptE.z, ptE.z - vtToolDir.z * CurrTool.GetHeigth())) ;

  // Parametri geometrici per 
  // determinare il vettore normale
   double dL = CurrTool.GetTipHeigth() * dMaxRad / dDeltaRad ;

  // Trapano
   if ( CurrTool.GetTipRadius() < CurrTool.GetRadius()) {

     // Ciclo sui punti
      for ( int i = nStartI ; i <= nEndI ; ++ i) {
         for ( 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.SqLenXY() ;

            if ( dSqDist < dSqMinRad)

               SubtractIntervals( nGrid, i, j, dZMin, dZMax, Z_AX, -Z_AX, CurrTool.GetToolNum()) ;

            else if ( dSqDist < dSqMaxRad) {

               double dr = sqrt( dSqDist) ;
               double dl = dr * dMaxRad / dL ;

               if ( vtToolDir.z > 0) {

                  double dMin = dZMin + CurrTool.GetTipHeigth() * ( dr - dMinRad) / dDeltaRad ;
                  double dMax = dZMax ;
                  
                  Point3d ptInt( dX, dY, dMin) ;
                  Point3d ptPn( ptO.x, ptO.y, dMin + dl) ;

                  Vector3d vtNmin = ptPn - ptInt ;
                  vtNmin.Normalize() ;

                  SubtractIntervals( nGrid, i, j, dMin, dMax, vtNmin, -Z_AX, CurrTool.GetToolNum()) ;
               }

               else {
                  double dMin = dZMin ;
                  double dMax = dZMax -  CurrTool.GetTipHeigth() * ( dr - dMinRad) / dDeltaRad ;

                  Point3d ptInt( dX, dY, dMax) ;
                  Point3d ptPn( ptO.x, ptO.y, dMax - dl) ;

                  Vector3d vtNmax = ptPn - ptInt ;
                  vtNmax.Normalize() ;

                  SubtractIntervals( nGrid, i, j, dMin, dMax, Z_AX, vtNmax, CurrTool.GetToolNum()) ;
               }
            }
         }
      }
   }

  // Coda di rondine
   else {

     // Ciclo sui punti
      for ( int i = nStartI ; i <= nEndI ; ++ i) {
         for ( 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.SqLenXY() ;
            
            if ( dSqDist < dSqMinRad) {
               SubtractIntervals( nGrid, i, j, dZMin, dZMax, Z_AX, -Z_AX, CurrTool.GetToolNum()) ;
            }

            else if ( dSqDist < dSqMaxRad) {

               double dr = sqrt( dSqDist) ;
               double dl = dr * dMaxRad / dL ;

               if ( vtToolDir.z > 0) {  
                  double dMin = dZMin ;
                  double dMax = dZMax - dStemHeigth - CurrTool.GetTipHeigth() * ( dr - dMinRad) / dDeltaRad ;

                  Point3d ptInt( dX, dY, dMax) ;
                  Point3d ptPn( ptO.x, ptO.y, dMax - dl) ;

                  Vector3d vtNmax = ptPn - ptInt ;
                  vtNmax.Normalize() ;
                 
                  SubtractIntervals( nGrid, i, j, dMin, dMax, Z_AX, vtNmax, CurrTool.GetToolNum()) ;  
               }

               else {
                  double dMin = dZMin + dStemHeigth + CurrTool.GetTipHeigth() * ( dr - dMinRad) / dDeltaRad ;
                  double dMax = dZMax ;

                  Point3d ptInt( dX, dY, dMin) ;
                  Point3d ptPn( ptO.x, ptO.y, dMin + dl) ;

                  Vector3d vtNmin = ptPn - ptInt ;
                  vtNmin.Normalize() ;

                  SubtractIntervals( nGrid, i, j, dMin, dMax, vtNmin, -Z_AX, CurrTool.GetToolNum()) ;  
               }
            }
         }
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::Conus_ZPerp( int nGrid, const Point3d & ptS, const Point3d & ptE, const Vector3d & vtToolDir)
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

  // Verifica sull'interferenza utensile Zmap
   int nStartI, nStartJ, nEndI, nEndJ ;
   if ( ! TestToolBBox( nGrid, ptS, ptE, vtToolDir, nStartI, nStartJ, nEndI, nEndJ))
      return true ;

  // Parametri geometrici dell'utensile
   double dStemHeigth = CurrTool.GetHeigth() - CurrTool.GetTipHeigth() ;
   double dMinRad = min( CurrTool.GetRadius(), CurrTool.GetTipRadius()) ;
   double dMaxRad = max( CurrTool.GetRadius(), CurrTool.GetTipRadius()) ;
   double dDeltaRad = dMaxRad - dMinRad ;
   double dSqMinRad = dMinRad * dMinRad ;
   double dSqMaxRad = dMaxRad * dMaxRad ;

   double dSafeMinRad = dMinRad - EPS_SMALL ;
   double dSafeMaxRad = dMaxRad - EPS_SMALL ;
   double dSafeSqMaxRad = dSafeMaxRad * dSafeMaxRad ;
   double dSafeSqMinRad = dSafeMinRad * dSafeMinRad ;
   
   Point3d ptSxy( ptS.x, ptS.y, 0) ;
   Point3d ptExy( ptE.x, ptE.y, 0) ;

   Vector3d vtV1( ptE.x - ptS.x, ptE.y - ptS.y, 0) ;
   Vector3d vtV2( - vtV1.y, vtV1.x, 0) ;

   double dLen = vtV1.LenXY() ;

   vtV1.Normalize() ;
   vtV2.Normalize() ;

   double dBaseZ = ptS.z ;
   double dStemZ = ptS.z - vtToolDir.z * dStemHeigth ;
   double dTipZ = ptS.z -  vtToolDir.z * CurrTool.GetHeigth() ;

  // Lunghezza cono
   double dL = CurrTool.GetTipHeigth() * dMaxRad / dDeltaRad ;

  // Punta a trapano
   if ( CurrTool.GetTipRadius() < CurrTool.GetRadius()) {

      Vector3d vtV = vtToolDir ;

      Point3d ptVS = ptS - vtToolDir * ( dStemHeigth + dL) ;
      Point3d ptVE = ptE - vtToolDir * ( dStemHeigth + dL) ;

      for ( int i = nStartI ; i <= nEndI ; ++ i) {
         for( 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 vtCS = ptC - ptSxy ;
            Vector3d vtCE = ptC - ptExy ;

            double dP1 = vtCS * vtV1 ;
            double dP2 = vtCS * vtV2 ;

            double dSqLenS = vtCS.SqLenXY() ;
            double dSqLenE = vtCE.SqLenXY() ;

            if ( dSqLenS < dSafeSqMaxRad  ||  dSqLenE < dSafeSqMaxRad  ||
                 ( abs( dP2) < dSafeMaxRad  &&  dP1 >= 0  &&  dP1 <= dLen)) {

               if ( dSqLenS < dSqMinRad  ||  dSqLenE < dSqMinRad  || 
                    ( abs( dP2) < dMinRad  &&  dP1 >= 0  &&  dP1 <= dLen)) {

                  double dMin = min( dBaseZ, dTipZ) ;
                  double dMax = max( dBaseZ, dTipZ) ;

                  SubtractIntervals( nGrid, i, j, dMin, dMax, Z_AX, -Z_AX, CurrTool.GetToolNum()) ;
               }

               else {
                  Vector3d vtNmin, vtNmax ;

                  if ( dP1 < 0) {

                     double dr = sqrt( dSqLenS) ;

                     double dMin = min( dBaseZ, dTipZ + ( dStemZ - dTipZ) * ( dr - dMinRad) / dDeltaRad) ;
                     double dMax = max( dBaseZ, dTipZ + ( dStemZ - dTipZ) * ( dr - dMinRad) / dDeltaRad) ;

                     if ( vtToolDir.z > 0) {
                        Point3d ptInt( dX, dY, dMin) ;
                        Vector3d vtU = ( ptInt - ptVS) - ( ptInt - ptVS) * vtV * vtV ;
                        vtU.Normalize() ;
                        vtNmin = dDeltaRad * vtV - CurrTool.GetTipHeigth() * vtU ;
                        vtNmin.Normalize() ;
                        vtNmax = - Z_AX ;
                     }

                     else {
                        Point3d ptInt( dX, dY, dMax) ;
                        Vector3d vtU = ( ptInt - ptVS) - ( ptInt - ptVS) * vtV * vtV ;
                        vtU.Normalize() ;
                        vtNmin = Z_AX ;
                        vtNmax = dDeltaRad * vtV - CurrTool.GetTipHeigth() * vtU ;
                        vtNmax.Normalize() ;
                     }
                     SubtractIntervals( nGrid, i, j, dMin, dMax, vtNmin, vtNmax, CurrTool.GetToolNum()) ;
                  }

                  else if ( dP1 <= dLen) {
                     double dr = abs( dP2) ;
                     double dMin = min( dBaseZ, dTipZ + ( dStemZ - dTipZ) * ( dr - dMinRad) / dDeltaRad) ;
                     double dMax = max( dBaseZ, dTipZ + ( dStemZ - dTipZ) * ( dr - dMinRad) / dDeltaRad) ;
                     Point3d ptVtemp = ptVS + vtV1 * dP1 ;

                     if ( vtToolDir.z > 0) {
                        Point3d ptInt( dX, dY, dMin) ;
                        Vector3d vtU = ( ptInt - ptVtemp) - ( ptInt - ptVtemp) * vtV * vtV ;
                        vtU.Normalize() ;
                        vtNmin = dDeltaRad * vtV - CurrTool.GetTipHeigth() * vtU ;
                        vtNmin.Normalize() ;
                        vtNmax = - Z_AX ;
                     }

                     else {
                        Point3d ptInt( dX, dY, dMax) ;
                        Vector3d vtU = ( ptInt - ptVtemp) - ( ptInt - ptVtemp) * vtV * vtV ;
                        vtU.Normalize() ;
                        vtNmin = Z_AX ;
                        vtNmax = dDeltaRad * vtV - CurrTool.GetTipHeigth() * vtU ;
                        vtNmax.Normalize() ;
                     }

                     SubtractIntervals( nGrid, i, j, dMin, dMax, vtNmin, vtNmax, CurrTool.GetToolNum()) ;
                  }

                  else {

                     double dr = sqrt( dSqLenE) ;

                     double dMin = min( dBaseZ, dTipZ + ( dStemZ - dTipZ) * ( dr - dMinRad) / dDeltaRad) ;
                     double dMax = max( dBaseZ, dTipZ + ( dStemZ - dTipZ) * ( dr - dMinRad) / dDeltaRad) ;

                     if ( vtToolDir.z > 0) {
                        Point3d ptInt( dX, dY, dMin) ;
                        Vector3d vtU = ( ptInt - ptVE) - ( ptInt - ptVE) * vtV * vtV ; 
                        vtU.Normalize() ;
                        vtNmin = dDeltaRad * vtV - CurrTool.GetTipHeigth() * vtU ;
                        vtNmin.Normalize() ;
                        vtNmax = -Z_AX ;
                     }

                     else {
                        Point3d ptInt( dX, dY, dMax) ;
                        Vector3d vtU = ( ptInt - ptVE) - ( ptInt - ptVE) * vtV * vtV ; 
                        vtU.Normalize() ;
                        vtNmin =  Z_AX ;
                        vtNmax = dDeltaRad * vtV - CurrTool.GetTipHeigth() * vtU ;
                        vtNmax.Normalize() ;
                     }

                     SubtractIntervals( nGrid, i, j, dMin, dMax, vtNmin, vtNmax, CurrTool.GetToolNum()) ;
                  }
               }
            }
         }
      }
   }

  // Coda di rondine
   else {

      Vector3d vtV = - vtToolDir ;

      Point3d ptVS = ptS - vtToolDir * ( CurrTool.GetHeigth() - dL) ;
      Point3d ptVE = ptE - vtToolDir * ( CurrTool.GetHeigth() - dL) ;

      for ( int i = nStartI ; i <= nEndI ; ++ i) {
         for( 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 vtCS = ptC - ptSxy ;
            Vector3d vtCE = ptC - ptExy ;

            double dP1 = vtCS * vtV1 ;
            double dP2 = vtCS * vtV2 ;

            double dSqLenS = vtCS.SqLenXY() ;
            double dSqLenE = vtCE.SqLenXY() ;

            if ( dSqLenS < dSqMaxRad  ||  dSqLenE < dSqMaxRad  || 
                 ( abs( dP2) < dMaxRad  &&  dP1 >= 0  &&  dP1 <= dLen)) {

               if ( dSqLenS < dSafeSqMinRad  ||  dSqLenE < dSafeSqMinRad  || 
                    ( abs( dP2) < dSafeMinRad  &&  dP1 >= 0  &&  dP1 <= dLen)) {
                  double dMin = min( dBaseZ, dTipZ) ;
                  double dMax = max( dBaseZ, dTipZ) ;
                  SubtractIntervals( nGrid, i, j, dMin, dMax, Z_AX, -Z_AX, CurrTool.GetToolNum()) ;
               }

               else {
                   Vector3d vtNmin, vtNmax ;

                  if ( dP1 < 0) {
                     double dr = sqrt( dSqLenS) ;
                     double dMin = min( dTipZ, dStemZ + ( dTipZ - dStemZ) * ( dr - dMinRad) / dDeltaRad) ;
                     double dMax = max( dTipZ, dStemZ + ( dTipZ - dStemZ) * ( dr - dMinRad) / dDeltaRad) ;

                     if ( vtToolDir.z > 0) {
                        Point3d ptInt( dX, dY, dMax) ;
                        Vector3d vtU = ( ptInt - ptVS) - ( ptInt - ptVS) * vtV * vtV ;
                        vtU.Normalize() ;
                        vtNmin = Z_AX ;
                        vtNmax = dDeltaRad * vtV - CurrTool.GetTipHeigth() * vtU ;
                        vtNmax.Normalize() ;
                     }

                     else {
                        Point3d ptInt( dX, dY, dMin) ;
                        Vector3d vtU = ( ptInt - ptVS) - ( ptInt - ptVS) * vtV * vtV ;
                        vtU.Normalize() ;
                        vtNmax = -Z_AX ;
                        vtNmin = dDeltaRad * vtV - CurrTool.GetTipHeigth() * vtU ;
                        vtNmin.Normalize() ;
                     }
                     SubtractIntervals( nGrid, i, j, dMin, dMax, vtNmin, vtNmax, CurrTool.GetToolNum()) ;
                  }

                  else if ( dP1 <= dLen) {
                      double dr = abs( dP2) ;
                      Point3d ptVtemp = ptVS + vtV1 * dP1 ;
                      double dMin = min( dTipZ, dStemZ + ( dTipZ - dStemZ) * ( dr - dMinRad) / dDeltaRad) ;
                      double dMax = max( dTipZ, dStemZ + ( dTipZ - dStemZ) * ( dr - dMinRad) / dDeltaRad) ;

                     if ( vtToolDir.z > 0) {
                        Point3d ptInt( dX, dY, dMax) ;
                        Vector3d vtU = ( ptInt - ptVtemp) - ( ptInt - ptVtemp) * vtV * vtV ;
                        vtU.Normalize() ;
                        vtNmin = Z_AX ;
                        vtNmax = dDeltaRad * vtV - CurrTool.GetTipHeigth() * vtU ;
                        vtNmax.Normalize() ;
                     }

                     else {
                        Point3d ptInt( dX, dY, dMin) ;
                        Vector3d vtU = ( ptInt - ptVtemp) - ( ptInt - ptVtemp) * vtV * vtV ;
                        vtU.Normalize() ;
                        vtNmax = -Z_AX ;
                        vtNmin = dDeltaRad * vtV - CurrTool.GetTipHeigth() * vtU ;
                        vtNmin.Normalize() ;
                     }

                     SubtractIntervals( nGrid, i, j, dMin, dMax, vtNmin, vtNmax, CurrTool.GetToolNum()) ;
                  }

                  else {
                     double dr = sqrt( dSqLenE) ;
                     double dMin = min( dTipZ, dStemZ + ( dTipZ - dStemZ) * ( dr - dMinRad) / dDeltaRad) ;
                     double dMax = max( dTipZ, dStemZ + ( dTipZ - dStemZ) * ( dr - dMinRad) / dDeltaRad) ;

                     if ( vtToolDir.z > 0) {
                        Point3d ptInt( dX, dY, dMax) ;
                        Vector3d vtU = ( ptInt - ptVE) - ( ptInt - ptVE) * vtV * vtV ;
                        vtU.Normalize() ;
                        vtNmin = Z_AX ;
                        vtNmax = dDeltaRad * vtV - CurrTool.GetTipHeigth() * vtU ;
                        vtNmax.Normalize() ;
                     }

                     else {
                        Point3d ptInt( dX, dY, dMin) ;
                        Vector3d vtU = ( ptInt - ptVE) - ( ptInt - ptVE) * vtV * vtV ;
                        vtU.Normalize() ;
                        vtNmax = -Z_AX ;
                        vtNmin = dDeltaRad * vtV - CurrTool.GetTipHeigth() * vtU ;
                        vtNmin.Normalize() ;
                     }

                     SubtractIntervals( nGrid, i, j, dMin, dMax, vtNmin, vtNmax, CurrTool.GetToolNum()) ;
                  }
               }
            }
         }
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
bool 
VolZmap::Conus_ZMilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir)
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

   double dStemH = CurrTool.GetHeigth() - CurrTool.GetTipHeigth() ;

   CompCyl_ZMilling( nGrid, ptS, ptE, vtToolDir, dStemH, CurrTool.GetRadius(), CurrTool.GetToolNum()) ;

   if ( CurrTool.GetTipRadius() < CurrTool.GetRadius()) {
      Point3d ptSC = ptS - vtToolDir * dStemH ;
      Point3d ptEC = ptE - vtToolDir * dStemH ;
      CompConus_ZMilling( nGrid, ptSC, ptEC, vtToolDir,
                          CurrTool.GetTipHeigth(), CurrTool.GetRadius(), CurrTool.GetTipRadius(), V_NULL, V_NULL, CurrTool.GetToolNum()) ;
   }

   else {
      Point3d ptSC = ptS - vtToolDir * CurrTool.GetHeigth() ;
      Point3d ptEC = ptE - vtToolDir * CurrTool.GetHeigth() ;
      CompConus_ZMilling( nGrid, ptSC, ptEC, - vtToolDir,
                          CurrTool.GetTipHeigth(), CurrTool.GetTipRadius(), CurrTool.GetRadius(), V_NULL, V_NULL, CurrTool.GetToolNum()) ;
   }

   return true ;
}

// --------- Mortasatrice -----------------------------------------------------

// ----------------------------------------------------------------------------
bool 
VolZmap::Mrt_ZDrilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, const Vector3d& vtAux)
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

  // Proiezione della traiettoria sul piano dei movimenti possibili
   Vector3d vtMoveOnP = ( ptE - ptS) * vtToolDir * vtToolDir ;
   Point3d ptEOnP = ptS + vtMoveOnP ;

  // Scompongo la mortasatrice in solidi semplici

  // Parallelepipedo di base
   double dLenX = CurrTool.GetMrtChsWidth() ;
   double dLenY =  CurrTool.GetMrtChsThickness() ;
   double dLenZ = CurrTool.GetHeigth() - CurrTool.GetCornRadius() ;

   CompPar_ZDrilling( nGrid, dLenX, dLenY, dLenZ, ptS, ptEOnP, vtToolDir, vtAux, CurrTool.GetToolNum()) ;
   
  // Se la punta è di tipo bull-nose
   if ( abs( CurrTool.GetMrtChsWidth() - 2 * CurrTool.GetCornRadius()) > EPS_SMALL) {

     // Parallelepipedo di punta 
      Point3d ptTipS = ptS - dLenZ * vtToolDir ;
      Point3d ptTipE = ptEOnP - dLenZ * vtToolDir ;

      dLenX = abs( CurrTool.GetMrtChsWidth() - 2 * CurrTool.GetCornRadius()) ;
      dLenZ = CurrTool.GetCornRadius() ;

      CompPar_ZDrilling( nGrid, dLenX, dLenY, dLenZ, ptTipS, ptTipE, vtToolDir, vtAux, CurrTool.GetToolNum()) ;

      Vector3d vtVOnP = vtToolDir ^ vtAux ;

     // Cilindri 
      Point3d ptSminus = ptTipS - ( 0.5 * dLenX) * vtVOnP + 0.5 * dLenY * vtAux ;
      Point3d ptEminus = ptTipE - ( 0.5 * dLenX) * vtVOnP + 0.5 * dLenY * vtAux ;  
      Point3d ptSplus = ptTipS + ( 0.5 * dLenX) * vtVOnP + 0.5 * dLenY * vtAux ;
      Point3d ptEplus = ptTipE + ( 0.5 * dLenX) * vtVOnP + 0.5 * dLenY * vtAux ;  

      CompCyl_Milling( nGrid, ptSminus, ptEminus, vtAux, dLenY, CurrTool.GetCornRadius(), false, false, CurrTool.GetToolNum()) ;
      CompCyl_Milling( nGrid, ptSplus, ptEplus, vtAux, dLenY, CurrTool.GetCornRadius(), false, false, CurrTool.GetToolNum()) ;
   }

  // se la punta è di tipo sfera
   else {
     // Cilindro
      Point3d ptCylS = ptS - dLenZ * vtToolDir + 0.5 * dLenY * vtAux ;
      Point3d ptCylE = ptEOnP - dLenZ * vtToolDir + 0.5 * dLenY * vtAux ;

      CompCyl_Milling( nGrid, ptCylS, ptCylE, vtAux, dLenY, CurrTool.GetCornRadius(), false, false, CurrTool.GetToolNum()) ;
   }

   return true ;
}


// --------- Chisel ----------------------------------------------------------

//----------------------------------------------------------------------------
bool 
VolZmap::Chs_ZDrilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, const Vector3d& vtAux) 
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

   CompPar_ZDrilling( nGrid, CurrTool.GetMrtChsWidth(), CurrTool.GetMrtChsThickness(), CurrTool.GetHeigth(),
                      ptS, ptE, vtToolDir, vtAux, CurrTool.GetToolNum()) ;
   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::Chs_ZMilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, const Vector3d& vtAux) 
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

  // Proiezione della traiettoria sul piano dei moti possibili
   Point3d ptSp = ptS ;
   Vector3d vtMoveOnP = ( ptE - ptS) - ( ptE - ptS) * vtToolDir * vtToolDir ;
   Point3d ptEp = ptS + vtMoveOnP ;

   CompPar_ZMilling( nGrid, CurrTool.GetMrtChsWidth(), CurrTool.GetMrtChsThickness(), CurrTool.GetHeigth(),
                     ptSp, ptEp, vtToolDir, vtAux, CurrTool.GetToolNum()) ;
   return true ;
}

// --------- Utensile generico ------------------------------------------------ 

//----------------------------------------------------------------------------
bool
VolZmap::GenTool_ZDrilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir)
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

  // Descrizione geometrica del moto  
   Point3d ptI = ptS ;
   Point3d ptF = ptE ;
   Vector3d vtMove = ptE - ptS ;
  // Vettore delle normali agli archi
   const VCT3DVECTOR& vArcNorm = CurrTool.GetArcNormalVec() ;
  // Poinché l'asse utensile è parallelo all'asse Z, definisco un sistema di
  // riferimento ad hoc in cui le normali agli archi giacciano nel piano XZ.
   Frame3d frNormFrame ;
   frNormFrame.Set( ORIG, X_AX, - Z_AX, Y_AX) ;
  // Ciclo sulle curve
   const CurveComposite& ToolProfile = CurrTool.GetApproxOutline() ;
   int i = - 1 ;
   const ICurve* pCurve = ToolProfile.GetCurve( ++ i) ;
   while ( pCurve != nullptr) {

      double dHeight = 0 ;

     // Se segmento
      if ( pCurve->GetType() == CRV_LINE) {
        // Recupero gli estremi
         const ICurveLine* pLine = GetCurveLine( pCurve) ;
         Point3d ptStart = pLine->GetStart() ;
         Point3d ptEnd = pLine->GetEnd() ;
         int nNormNum = pLine->GetTempProp() ;
         Vector3d vtNormSt, vtNormEn ;
         if ( nNormNum != 0) {
            vtNormSt = vArcNorm[nNormNum - 1] ;
            vtNormEn = vArcNorm[nNormNum] ;
            vtNormSt.ToLoc( frNormFrame) ;
            vtNormEn.ToLoc( frNormFrame) ; 
         }
        // Ne determino l'altezza
         dHeight = abs( ptStart.y - ptEnd.y) ;
         if ( dHeight > EPS_SMALL) {
           // Se X costante, è un cilindro
            if ( abs( ptStart.x - ptEnd.x) < EPS_SMALL) {
               double dRadius = ptStart.x ;
               if ( dRadius > 10 * EPS_SMALL)
                  CompCyl_ZDrilling( nGrid, ptI, ptF, vtToolDir, dHeight, dRadius, CurrTool.GetToolNum()) ;
            }
           // Se X crescente, è un cono con vettore equiverso a quello dell'utensile
            else if ( ptStart.x > ptEnd.x) {
               double dMaxRad = ptStart.x ;
               double dMinRad = ptEnd.x ;
               CompConus_ZDrilling( nGrid, ptI, ptF, vtToolDir, dHeight, dMaxRad, dMinRad, vtNormSt, vtNormEn, CurrTool.GetToolNum()) ;
            }
           // Se X decrescente, è un cono con vettore opposto a quello dell'utensile
            else if ( ptStart.x < ptEnd.x) {
               double dMaxRad = ptEnd.x ; 
               double dMinRad = ptStart.x ;
               Point3d ptIn = ptI - vtToolDir * dHeight ;
               Point3d ptFn = ptIn + vtMove ;
               vtNormEn.z *= -1 ;
               vtNormSt.z *= -1 ;
               CompConus_ZDrilling( nGrid, ptIn, ptFn, - vtToolDir, dHeight, dMaxRad, dMinRad, vtNormEn, vtNormSt, CurrTool.GetToolNum()) ;
            }
         }
      }

     // se altrimenti arco
      else if ( pCurve->GetType() == CRV_ARC) {
        // Recupero estremi, centro e raggio
         const ICurveArc* pArc = GetCurveArc( pCurve) ;
         Point3d ptStart ; pArc->GetStartPoint( ptStart) ;
         Point3d ptEnd ; pArc->GetEndPoint( ptEnd) ;
         Point3d ptCen = pArc->GetCenter() ;
         double dRadius = pArc->GetRadius() ;
        // Determino le posizioni iniziale e finale del centro della sfera
         Point3d ptCenS = ptI - vtToolDir * ( ptStart.y - ptCen.y) ;
         Point3d ptCenE = ptCenS + vtMove ;
        // Eseguo l'asportazione del materiale
         CompBall_Milling( nGrid, ptCenS, ptCenE, dRadius, CurrTool.GetToolNum()) ;
        // aggiorno l'altezza
         dHeight = abs( ptStart.y - ptEnd.y) ;
      }

     // Determino le posizioni iniziale e finale del componente successivo
      ptI = ptI - vtToolDir * dHeight ;
      ptF = ptI + vtMove ;

     // Passo alla curva successiva del profilo
      pCurve = ToolProfile.GetCurve( ++ i) ;
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::GenTool_ZMilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir)
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

  // Descrizione geometrica del moto  
   Point3d ptI = ptS ;
   Point3d ptF = ptE ;
   Vector3d vtMove = ptE - ptS ;
  // Vettore delle normali agli archi
   const VCT3DVECTOR& vArcNorm = CurrTool.GetArcNormalVec() ;
  // Poinché l'asse utensile è parallelo all'asse Z, definisco un sistema di
  // riferimento ad hoc in cui le normali agli archi giacciano nel piano XZ.
   Frame3d frNormFrame;
   frNormFrame.Set( ORIG, X_AX, -Z_AX, Y_AX) ;
  // Ciclo sulle curve del profilo utensile
   const CurveComposite& ToolProfile = CurrTool.GetApproxOutline() ;
   int i = - 1 ;
   const ICurve* pCurve = ToolProfile.GetCurve( ++ i) ;
   while ( pCurve != nullptr) {

      double dHeight = 0 ;

     // Se segmento
      if ( pCurve->GetType() == CRV_LINE) {
        // Recupero gli estremi
         const ICurveLine* pLine = GetCurveLine( pCurve) ;
         Point3d ptStart = pLine->GetStart() ;
         Point3d ptEnd = pLine->GetEnd() ;
         int nNormNum = pLine->GetTempProp() ;
         Vector3d vtNormSt, vtNormEn ;
         if ( nNormNum != 0) {
            vtNormSt = vArcNorm[nNormNum - 1] ;
            vtNormEn = vArcNorm[nNormNum] ;
            vtNormSt.ToLoc( frNormFrame) ;
            vtNormEn.ToLoc( frNormFrame) ;
         }
        // Ne determino l'altezza
         dHeight = abs( ptStart.y - ptEnd.y) ;
         if ( dHeight > EPS_SMALL) {
           // Se X costante, è un cilindro
            if ( abs( ptStart.x - ptEnd.x) < EPS_SMALL) {
               double dRadius = ptStart.x ;
               if ( dRadius > 10 * EPS_SMALL)
                  CompCyl_ZMilling( nGrid, ptI, ptF, vtToolDir, dHeight, dRadius, CurrTool.GetToolNum()) ;
            }
           // Se X crescente, è un cono con vettore equiverso a quello dell'utensile
            else if ( ptStart.x > ptEnd.x) {
               double dMaxRad = ptStart.x ;
               double dMinRad = ptEnd.x ;
               CompConus_ZMilling( nGrid, ptI, ptF, vtToolDir, dHeight, dMaxRad, dMinRad, vtNormSt, vtNormEn, CurrTool.GetToolNum()) ;
            }
           // Se X decrescente, è un cono con vettore opposto a quello dell'utensile
            else if ( ptStart.x < ptEnd.x) {
               double dMaxRad = ptEnd.x ; 
               double dMinRad = ptStart.x ;
               Point3d ptIn = ptI - vtToolDir * dHeight ;
               Point3d ptFn = ptIn + vtMove ;
               vtNormEn.z *= -1 ;
               vtNormSt.z *= -1 ;
               CompConus_ZMilling( nGrid, ptIn, ptFn, - vtToolDir, dHeight, dMaxRad, dMinRad, vtNormEn, vtNormSt, CurrTool.GetToolNum()) ;
            }
         }
      }

     // se altrimenti arco
      else if ( pCurve->GetType() == CRV_ARC) { 
        // Recupero estremi, centro e raggio
         const ICurveArc* pArc = GetCurveArc( pCurve) ;
         Point3d ptStart ; pArc->GetStartPoint( ptStart) ;
         Point3d ptEnd ; pArc->GetEndPoint( ptEnd) ;
         Point3d ptCen = pArc->GetCenter() ;
         double dRadius = pArc->GetRadius() ;
        // Determino le posizioni iniziale e finale del centro della sfera
         Point3d ptCenS = ptI - vtToolDir * ( ptStart.y - ptCen.y) ;
         Point3d ptCenE = ptCenS + vtMove ;
        // Eseguo l'asportazione del materiale
         CompBall_Milling( nGrid, ptCenS, ptCenE, dRadius, CurrTool.GetToolNum()) ;
        // aggiorno l'altezza
         dHeight = abs( ptStart.y - ptEnd.y) ;
      } 

     // Determino le posizioni iniziale e finale del componente successivo
      ptI = ptI - vtToolDir * dHeight ;
      ptF = ptI + vtMove ;

     // Passo alla curva successiva del profilo
      pCurve = ToolProfile.GetCurve( ++ i) ;
   }

   return true ;
}

// ---------- VERSORE UTENSILE NEL PIANO XY ----------------------------------

// --------- Cilindro e sfera ------------------------------------------------ 
//----------------------------------------------------------------------------
bool
VolZmap::CylBall_XYDrilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir)
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

  // Verifica sull'interferenza utensile Zmap
   int nStartI, nStartJ, nEndI, nEndJ ;
   if ( ! TestToolBBox( nGrid, ptS, ptE, vtToolDir, nStartI, nStartJ, nEndI, nEndJ))
      return true ;

  // Parametri geometrici dell'utensile e quota Z del movimento
   double dStemHeigth = CurrTool.GetHeigth() - CurrTool.GetTipHeigth() ;
   double dSqRad = CurrTool.GetRadius() * CurrTool.GetRadius() ;
   double dZ = ptS.z ;
  
  // Vettore movimento e sua lunghezza
   Vector3d vtMove = ptE - ptS ;  double dLen = vtMove.LenXY() ;

  // Definizione di un sistema di riferimento ad hoc
   Point3d ptI = ( vtMove * vtToolDir > 0 ? ptE : ptS) ;
   Point3d ptF = ( vtMove * vtToolDir > 0 ? ptS - dStemHeigth * vtToolDir : ptE - dStemHeigth * vtToolDir) ;
   Point3d ptIxy( ptI.x, ptI.y, 0) ;
   Point3d ptFxy( ptF.x, ptF.y, 0) ;

   Vector3d vtV1( - vtToolDir.x, - vtToolDir.y, 0) ;
   vtV1.Normalize() ;
   Vector3d vtV2 = vtV1 ;
   vtV2.Rotate( Z_AX, 0, 1) ;

   for ( int i = nStartI ; i <= nEndI ; ++ i) {
      for ( 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 - ptIxy ;

         double dP1 = vtC * vtV1 ;           double dP2 = vtC * vtV2 ;

         Vector3d vtBall = ptC - ptFxy ;     double dSqLen = vtBall.SqLenXY() ;

        // Zona lavorata dalla parte cilindrica
         if ( dP1 > EPS_SMALL && dP1 < dStemHeigth + dLen - EPS_SMALL && 
              abs( dP2) < CurrTool.GetRadius()) {

            double dH = sqrt( dSqRad - dP2 * dP2) ;
            double dMin = dZ - dH ;
            double dMax = dZ + dH ;

            Point3d ptIntMin( dX, dY, dMin) ;
            Point3d ptIntMax( dX, dY, dMax) ;

            Vector3d vtMin = ( ptI - ptIntMin) -
                             ( ptI - ptIntMin) * vtToolDir * vtToolDir ;
            Vector3d vtMax = ( ptI - ptIntMax) - 
                             ( ptI - ptIntMax) * vtToolDir * vtToolDir ;

            vtMin.Normalize() ;
            vtMax.Normalize() ;

            SubtractIntervals( nGrid, i, j, dMin, dMax, vtMin, vtMax, CurrTool.GetToolNum()) ;
         }

        // Se l'utensile è sferico sottraggo anche la punta
         if ( CurrTool.GetType() == Tool::BALLMILL)
            if ( dSqLen < dSqRad) { // LA SOLUZIONE MOMENTANEA è CREARE UTENSILE GENERICO SE LO STELO è PIù CORTO DEL RAGGIO

               double dH = sqrt( dSqRad - dSqLen) ;
               double dMin = dZ - dH ;
               double dMax = dZ + dH ;

               Vector3d vtMin = ptF - Point3d( dX, dY, dMin) ;
               Vector3d vtMax = ptF - Point3d( dX, dY, dMax) ;

               vtMin.Normalize() ;
               vtMax.Normalize() ;

               SubtractIntervals( nGrid, i, j, dZ - dH, dZ + dH, vtMin, vtMax, CurrTool.GetToolNum()) ;
            }
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::CylBall_XYPerp( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d & vtToolDir) 
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

  // Verifica sull'interferenza utensile Zmap
   int nStartI, nStartJ, nEndI, nEndJ ;
   if ( ! TestToolBBox( nGrid, ptS, ptE, vtToolDir, nStartI, nStartJ, nEndI, nEndJ))
      return true ;

  // Parametri geometrici dell'utensile
   double dStemHeigth = CurrTool.GetHeigth() - CurrTool.GetTipHeigth() ;
   double dSqRad = CurrTool.GetRadius() * CurrTool.GetRadius() ;
   double dSafeSqRad = dSqRad - 2 * CurrTool.GetRadius() * EPS_SMALL ;

  // Studio simmetrie del problema
   Point3d ptI = ( ptS.z <= ptE.z ? ptS : ptE) ;
   Point3d ptF = ( ptS.z <= ptE.z ? ptE : ptS) ;
   // elimino eventuale piccolo errore di perpendicolarità del movimento rispetto all'utensile
   Vector3d vtErr = 0.5 * ( ( ptF - ptI) * vtToolDir) * vtToolDir ;
   ptI += vtErr ;
   ptF -= vtErr ;
   Point3d ptIxy( ptI.x, ptI.y, 0) ;

   Point3d ptIStemT = ptI - vtToolDir * dStemHeigth ;
   Point3d ptFStemT = ptF - vtToolDir * dStemHeigth ;

  // Quote punti iniziale e finale
   double dZI = ptI.z ;
   double dZF = ptF.z ;
   double dDeltaZ = ptF.z - ptI.z ;

  // Vettori caratterizzanti il moto 
   Vector3d vtMove = ptF - ptI ;
   double dLenXY = vtMove.LenXY() ;
   vtMove.Normalize() ;
  
  // Sistema di riferimento ad hoc
   Vector3d vtV1( - vtToolDir.x, - vtToolDir.y, 0) ;
   vtV1.Normalize() ;
   Vector3d vtV2 = vtV1 ;
   vtV2.Rotate( Z_AX, 0, 1) ;
   if ( vtV2 * vtMove < 0)
      vtV2 = - vtV2 ;

  // Vettori e punti determinanti i piani
   Vector3d vtP = vtMove ;  // Se dLen < EPS_SMALL non si usa
   vtP.Rotate( vtToolDir, 0, 1) ;
   Point3d ptUp = ptI + CurrTool.GetRadius() * ( vtP.z > 0 ? vtP : - vtP) ;
   Point3d ptDw = ptI + CurrTool.GetRadius() * ( vtP.z > 0 ? - vtP : vtP) ;

   Vector3d vtPXY( vtP.x, vtP.y, 0) ;

   Vector3d vtUp = ptUp - ORIG ;  double dDotUp = vtUp * vtP ;
   Vector3d vtDw = ptDw - ORIG ;  double dDotDw = vtDw * vtP ;

   double dSmall = CurrTool.GetRadius() * vtPXY.LenXY() ;

  // Parte sferica
   double dCos = vtMove.z ; // vtMove.z > 0 : ptF.z >= ptI.z
   double dSin = ( dCos < 1 ? sqrt( 1 - dCos * dCos) : 0) ;

   double dInfZ, dSupZ ;

   if ( dLenXY < EPS_SMALL) {

      for ( int i = nStartI ; i <= nEndI ; ++ i) {
         for ( 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 - ptIxy ;

            double dP1 = vtC * vtV1 ;
            double dP2 = vtC * vtV2 ;

            Vector3d vtNmin, vtNmax ;

            if ( dP1 > EPS_SMALL  &&  dP1 < dStemHeigth - EPS_SMALL  && 
                 dP2 > - CurrTool.GetRadius() + EPS_SMALL  &&
                 dP2 < CurrTool.GetRadius() - EPS_SMALL) {

               dInfZ = ptI.z - sqrt( dSqRad - dP2 * dP2) ;
               dSupZ = ptF.z + sqrt( dSqRad - dP2 * dP2) ;

               Point3d ptIntInf( dX, dY, dInfZ) ;
               Point3d ptIntSup( dX, dY, dSupZ) ;

               vtNmin = - ( ptIntInf - ptI) + ( ptIntInf - ptI) * vtV1 * vtV1 ;
               vtNmax = - ( ptIntSup - ptF) + ( ptIntSup - ptF) * vtV1 * vtV1 ;

               vtNmin.Normalize() ;
               vtNmax.Normalize() ;

               SubtractIntervals( nGrid, i, j, dInfZ, dSupZ, vtNmin, vtNmax, CurrTool.GetToolNum()) ;
            }

            if ( CurrTool.GetType() == Tool::BALLMILL) {

               if ( dP1 > dStemHeigth - 2 * EPS_SMALL && ( dP1 - dStemHeigth) * ( dP1 - dStemHeigth) + dP2 * dP2 < dSafeSqRad) {

                  dInfZ = ptI.z - sqrt( dSqRad - ( dP1 - dStemHeigth) * ( dP1 - dStemHeigth) - dP2 * dP2) ;
                  dSupZ = ptF.z + sqrt( dSqRad - ( dP1 - dStemHeigth) * ( dP1 - dStemHeigth) - dP2 * dP2) ;

                  Point3d ptIntInf( dX, dY, dInfZ) ;
                  Point3d ptIntSup( dX, dY, dSupZ) ;

                  vtNmin = ptIStemT - ptIntInf ;
                  vtNmax = ptFStemT - ptIntSup ;

                  vtNmin.Normalize() ;
                  vtNmax.Normalize() ;

                  SubtractIntervals( nGrid, i, j, dInfZ, dSupZ, vtNmin, vtNmax, CurrTool.GetToolNum()) ;
               }
            }
         }
      }
   }

   else {
      for ( int i = nStartI ; i <= nEndI ; ++ i) {
         for ( 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 - ptIxy ;

            double dP1 = vtC * vtV1 ;
            double dP2 = vtC * vtV2 ;

           // Parte cilindrica 
            if ( dP1 > EPS_SMALL  &&  dP1 < dStemHeigth - EPS_SMALL) {

               if ( dP2 > - CurrTool.GetRadius()  &&  dP2 < dLenXY + CurrTool.GetRadius()) {

                  Vector3d vtNmin, vtNmax ;

                 // Massimi 
                  if ( dP2 < - dSmall + EPS_SMALL) {
                     double dHsq = dSqRad - dP2 * dP2 ;
                     double dH = ( dHsq > 0 ? sqrt( dHsq) : 0) ;
                     dSupZ = dZI + dH ; 
                     Point3d ptInt( dX, dY, dSupZ) ;
                     vtNmax = - ( ptInt - ptI - ( ptInt - ptI) * vtToolDir * vtToolDir) ;
                     vtNmax.Normalize() ;
                  }

                  else if ( dP2 < dLenXY - dSmall - EPS_SMALL) {
                     dSupZ = ( dDotUp - dX * vtP.x - dY * vtP.y) / vtP.z ;
                     Vector3d vtT = vtToolDir ^ vtMove ; 
                     vtNmax = ( vtT.z < 0 ? vtT : - vtT) ;
                     vtNmax.Normalize() ;
                  }

                  else {
                     double dH = sqrt( dSqRad - ( dP2 - dLenXY) * ( dP2 - dLenXY)) ;
                     dSupZ = dZF + dH ;
                     Point3d ptInt( dX, dY, dSupZ) ;
                     vtNmax = - ( ptInt - ptF - ( ptInt - ptF) * vtToolDir * vtToolDir) ;
                     vtNmax.Normalize() ;
                  }

                 // Minimi
                  if ( dP2 < dSmall + EPS_SMALL) {
                     double dHsq = dSqRad - dP2 * dP2 ;
                     double dH = ( dHsq > 0 ? sqrt( dHsq) : 0) ;
                     dInfZ = dZI - dH ;    
                     Point3d ptInt( dX, dY, dInfZ) ;
                     vtNmin = - ( ( ptInt - ptI) - ( ptInt - ptI) * vtToolDir * vtToolDir) ;
                     vtNmin.Normalize() ;
                  }

                  else if ( dP2 < dLenXY + dSmall - EPS_SMALL) {
                     dInfZ = ( dDotDw - dX * vtP.x - dY * vtP.y) / vtP.z ;
                     Vector3d vtT = vtToolDir ^ vtMove ; 
                     vtNmin = ( vtT.z < 0 ? - vtT : vtT) ;
                     vtNmin.Normalize() ;
                  }

                  else {
                     double dH = sqrt( dSqRad - ( dP2 - dLenXY) * ( dP2 - dLenXY)) ;
                     dInfZ = dZF - dH ;
                     Point3d ptInt( dX, dY, dInfZ) ;
                     vtNmin = - ( ( ptInt - ptF) -   ( ptInt - ptF) * vtToolDir * vtToolDir) ;
                     vtNmin.Normalize() ;
                  }

                  SubtractIntervals( nGrid, i, j, dInfZ, dSupZ, vtNmin, vtNmax, CurrTool.GetToolNum()) ;
               }
            }

           // Se l'utensile è ball-end sottraggo la punta
            if ( CurrTool.GetType() == Tool::BALLMILL) {

               if ( ( ( dP1 - dStemHeigth) * ( dP1 - dStemHeigth) + dP2 * dP2 < dSqRad  ||
                      ( dP1 - dStemHeigth) * ( dP1 - dStemHeigth) + ( dP2 - dLenXY) * ( dP2 - dLenXY) < dSqRad  ||
                      ( dP2 > 0 && dP2 < dLenXY  &&  dP1 < CurrTool.GetHeigth()))  &&  ( dP1 > dStemHeigth - 2 * EPS_SMALL)) {

                  double dSqRoot = sqrt( dSqRad - ( dP1 - dStemHeigth) * ( dP1 - dStemHeigth)) ;
                  double dP2_0 = dCos * dSqRoot ;
                  double dH0 = dSin * dSqRoot ;

                  double dMin, dMax ;

                  Vector3d vtNmin, vtNmax ;

                 // Massimo
                  if ( dP2 < - dP2_0) {
                     dMax = dZI + sqrt( dSqRad - ( dP1 - dStemHeigth) * ( dP1 - dStemHeigth) - dP2 * dP2) ;
                     Point3d ptInt( dX, dY, dMax) ;
                     vtNmax = ptIStemT - ptInt ;
                     vtNmax.Normalize() ;
                  }

                  else if ( dP2 < dLenXY - dP2_0) {
                     dMax = dZI + dH0 + dDeltaZ * ( dP2 + dP2_0) / dLenXY ;
                     Point3d ptInt( dX, dY, dMax) ;
                     vtNmax = - ( ptInt - ptIStemT) + ( ptInt - ptIStemT) * vtMove * vtMove ;
                     vtNmax.Normalize() ;
                  }

                  else { 
                     dMax = dZF + sqrt( dSqRad - ( dP1 - dStemHeigth) * ( dP1 - dStemHeigth) - ( dP2 - dLenXY) * ( dP2 - dLenXY)) ;
                     Point3d ptInt( dX, dY, dMax) ;
                     vtNmax = ptFStemT - ptInt;
                     vtNmax.Normalize() ;
                  }

                 // Minimo
                  if ( dP2 < dP2_0) {
                     dMin = dZI - sqrt( dSqRad - ( dP1 - dStemHeigth) * ( dP1 - dStemHeigth) - dP2 * dP2) ;
                     Point3d ptInt( dX, dY, dMin) ;
                     vtNmin = ptIStemT - ptInt ;
                     vtNmin.Normalize() ;
                  }

                  else if ( dP2 < dLenXY + dP2_0) {
                     dMin =  dZI - dH0 + dDeltaZ * ( dP2 - dP2_0) / dLenXY ;
                     Point3d ptInt( dX, dY, dMin) ;
                     vtNmin = - ( ptInt - ptIStemT) + ( ptInt - ptIStemT) * vtMove * vtMove ;
                     vtNmin.Normalize() ;
                  }

                  else {
                     dMin = dZF - sqrt( dSqRad - ( dP1 - dStemHeigth) * ( dP1 - dStemHeigth) - ( dP2 - dLenXY) * ( dP2 - dLenXY)) ;
                     Point3d ptInt( dX, dY, dMin) ;
                     vtNmin = ptFStemT - ptInt ;
                     vtNmin.Normalize() ;
                  }

                  SubtractIntervals( nGrid, i, j, dMin, dMax, vtNmin, vtNmax, CurrTool.GetToolNum()) ;
               }
            }
         }
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::CylBall_XYMilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir)
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

   if ( CurrTool.GetType() == Tool::CYLMILL)
      return CompCyl_Milling( nGrid, ptS, ptE, vtToolDir, CurrTool.GetHeigth(), CurrTool.GetRadius(), false, false, CurrTool.GetToolNum()) ;

   else if ( CurrTool.GetType() == Tool::BALLMILL) {
      double dHei = CurrTool.GetHeigth() - CurrTool.GetTipHeigth() ; 
      CompCyl_Milling( nGrid, ptS, ptE, vtToolDir, dHei, CurrTool.GetRadius(), false, false, CurrTool.GetToolNum()) ;
      CompBall_Milling( nGrid, ptS - dHei * vtToolDir, ptE - dHei * vtToolDir, CurrTool.GetRadius(), CurrTool.GetToolNum()) ;
      return true ;
   }

   else
      return false ;
}


// --------- Coni ------------------------------------------------------------
//----------------------------------------------------------------------------
bool
VolZmap::Conus_XYDrilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir)
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

  // Verifica sull'interferenza utensile Zmap
   int nStartI, nStartJ, nEndI, nEndJ ;
   if ( ! TestToolBBox( nGrid, ptS, ptE, vtToolDir, nStartI, nStartJ, nEndI, nEndJ))
      return true ;

  // Parametri geometrici dell'utensile
   double dStemHeigth = CurrTool.GetHeigth() - CurrTool.GetTipHeigth() ;
   double dMinRad = min( CurrTool.GetRadius(), CurrTool.GetTipRadius()) ;
   double dMaxRad = max( CurrTool.GetRadius(), CurrTool.GetTipRadius()) ;
   double dDeltaRad = dMaxRad - dMinRad ;

  // Geometria del moto
   double dLenXY = ( ptE - ptS).LenXY() ;
   Point3d ptI = ( vtToolDir * ( ptE - ptS) < 0 ? ptS : ptE) ;
   double dMatStemLen = ( CurrTool.GetRadius() > CurrTool.GetTipRadius() ? dStemHeigth + dLenXY : dStemHeigth) ;
   double dSqTipRad = CurrTool.GetTipRadius() * CurrTool.GetTipRadius() ;
   double dSqRad = CurrTool.GetRadius() * CurrTool.GetRadius() ;

  // Determinazione del vertice del cono
   double dL = CurrTool.GetTipHeigth() * dMaxRad / dDeltaRad ;

   Point3d ptV ;  // Vertice
   Vector3d vtV ; // Vettore riferimento con origine nel vertice

   if ( CurrTool.GetRadius() > CurrTool.GetTipRadius()) {
      vtV = vtToolDir ;
      ptV = ( vtToolDir * ( ptE - ptS) < 0 ? ptE : ptS) - vtToolDir * ( dStemHeigth + dL) ;
   }

   else {
      vtV = - vtToolDir ;
      ptV = ( vtToolDir * ( ptE - ptS) < 0 ? ptS : ptE) - vtToolDir * ( CurrTool.GetHeigth() - dL) ;
   }

  // Sistema di riferimento sul piano
   Vector3d vtV1 = - vtToolDir ;
   Vector3d vtV2 = vtV1 ;
   vtV2.Rotate( Z_AX, 0, 1) ;

  // Proiezione di ptI sul piano
   Point3d ptIxy( ptI.x, ptI.y, 0) ;

  // Ciclo sui punti
   for ( int i = nStartI ; i <= nEndI ; ++ i) {
      for ( 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 - ptIxy ;

         double dX1 = vtC * vtV1 ;
         double dX2 = vtC * vtV2 ;

         double dr = CurrTool.GetRadius() + ( dX1 - dMatStemLen) * ( CurrTool.GetTipRadius() - CurrTool.GetRadius()) / CurrTool.GetTipHeigth() ; 

         if ( dX1 > EPS_SMALL  &&  dX1 < dMatStemLen  &&  abs( dX2) < CurrTool.GetRadius()) {

            double dH = sqrt( dSqRad - dX2 * dX2) ;
            double dMin = ptI.z - dH ;
            double dMax = ptI.z + dH ;

            Point3d ptIntMin( dX, dY, dMin) ;
            Vector3d vtNmin = ( ptI - ptIntMin) - ( ptI - ptIntMin) * vtV1 * vtV1 ;
            vtNmin.Normalize() ;

            Point3d ptIntMax( dX, dY, dMax) ;
            Vector3d vtNmax = ( ptI - ptIntMax) - ( ptI - ptIntMax) * vtV1 * vtV1 ;
            vtNmax.Normalize() ;

            SubtractIntervals( nGrid, i, j, ptI.z - dH, ptI.z + dH, vtNmin, vtNmax, CurrTool.GetToolNum()) ;
         }

         else if ( dX1 >= dMatStemLen  && 
                   dX1 < dMatStemLen + CurrTool.GetTipHeigth() - EPS_SMALL  &&
                   abs( dX2) < dr) {

            double dH = sqrt( dr * dr - dX2 * dX2) ;
            double dMin = ptI.z - dH ;
            double dMax = ptI.z + dH ;

            Point3d ptIntMin( dX, dY, dMin) ;
            Vector3d vtUmin = ( ptIntMin - ptV) - ( ptIntMin - ptV) * vtV * vtV ;
            vtUmin.Normalize() ;
            Vector3d vtNmin = dDeltaRad * vtV - CurrTool.GetTipHeigth() * vtUmin ;
            vtNmin.Normalize() ;

            Point3d ptIntMax( dX, dY, dMax) ;
            Vector3d vtUmax = ( ptIntMax - ptV) - ( ptIntMax - ptV) * vtV * vtV ;
            vtUmax.Normalize() ;
            Vector3d vtNmax = dDeltaRad * vtV - CurrTool.GetTipHeigth() * vtUmax ;
            vtNmax.Normalize() ;

            SubtractIntervals( nGrid, i, j, dMin, dMax, vtNmin, vtNmax, CurrTool.GetToolNum()) ;
         }

         if ( CurrTool.GetTipRadius() >= CurrTool.GetRadius()) {

            if ( dX1 > dMatStemLen + CurrTool.GetTipHeigth() - 2 * EPS_SMALL  && 
                 dX1 < dMatStemLen + CurrTool.GetTipHeigth() + dLenXY - EPS_SMALL  &&
                 abs( dX2) < dSqTipRad) {

               double dH = sqrt( dSqTipRad - dX2 * dX2) ;
               double dMin = ptI.z - dH ;
               double dMax = ptI.z + dH ;

               Point3d ptIntMin( dX, dY, dMin) ;
               Vector3d vtNmin = ( ptI - ptIntMin) - ( ptI - ptIntMin) * vtV1 * vtV1 ;
               vtNmin.Normalize() ;

               Point3d ptIntMax( dX, dY, dMax) ;
               Vector3d vtNmax = ( ptI - ptIntMax) - ( ptI - ptIntMax) * vtV1 * vtV1 ;
               vtNmax.Normalize() ;

               SubtractIntervals( nGrid, i, j, ptI.z - dH, ptI.z + dH, vtNmin, vtNmax, CurrTool.GetToolNum()) ;
            }
         }
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::Conus_XYPerp( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir) 
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

  // Verifica sull'interferenza utensile Zmap
   int nStartI, nStartJ, nEndI, nEndJ ;
   if ( ! TestToolBBox( nGrid, ptS, ptE, vtToolDir, nStartI, nStartJ, nEndI, nEndJ))
      return true ;

  // Parametri geometrici dell'utensile
   double dStemHeigth = CurrTool.GetHeigth() - CurrTool.GetTipHeigth()  ;
   double dMinRad = min( CurrTool.GetRadius(), CurrTool.GetTipRadius()) ;
   double dMaxRad = max( CurrTool.GetRadius(), CurrTool.GetTipRadius()) ;
   double dSqRad = CurrTool.GetRadius() * CurrTool.GetRadius() ;
   double dDeltaRad = dMaxRad - dMinRad ;

  // Studio delle simmetrie del moto
   Point3d ptI = ( ptS.z < ptE.z ? ptS : ptE) ;
   Point3d ptF = ( ptS.z < ptE.z ? ptE : ptS) ;
   // elimino eventuale piccolo errore di perpendicolarità del movimento rispetto all'utensile
   Vector3d vtErr = 0.5 * ( ( ptF - ptI) * vtToolDir) * vtToolDir ;
   ptI += vtErr ;
   ptF -= vtErr ;
   Point3d ptIxy( ptI.x, ptI.y, 0) ;
   Point3d ptFxy( ptF.x, ptF.y, 0) ;

  // Cinematica del moto
   Vector3d vtMove = ptF - ptI ;
   double dLengthPathXY = vtMove.LenXY() ;
   double dDeltaZ = ptF.z - ptI.z ;
   double dZI = ptI.z ;
   double dZF = ptF.z ;
   vtMove.Normalize() ;

  // Riferimento coni
   double dL = CurrTool.GetTipHeigth() * dMaxRad / dDeltaRad ;
   Vector3d vtV = ( CurrTool.GetRadius() > CurrTool.GetTipRadius() ? vtToolDir : - vtToolDir) ; 
   Point3d ptVI = ptI - ( CurrTool.GetRadius() > CurrTool.GetTipRadius() ? dStemHeigth + dL : CurrTool.GetHeigth() - dL) * vtToolDir ;
   Point3d ptVF = ptF - ( CurrTool.GetRadius() > CurrTool.GetTipRadius() ? dStemHeigth + dL : CurrTool.GetHeigth() - dL) * vtToolDir ;

  // Movimento verticale 
   if ( dLengthPathXY < EPS_SMALL) {
   
     // Riferimento sul piano
      Vector3d vtV1 = - vtToolDir ;
      Vector3d vtV2( - vtV1.y, vtV1.x, 0) ;

      for ( int i = nStartI ; i <= nEndI ; ++ i) {
         for ( 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 - ptIxy ;

            double dP1 = vtC * vtV1 ;
            double dP2 = vtC * vtV2 ;
          
           // Parte cilindrica 
            if ( dP1 > EPS_SMALL  &&  dP1 <= dStemHeigth  && 
                 dP2 > - CurrTool.GetRadius() + EPS_SMALL  &&  dP2 < CurrTool.GetRadius() - EPS_SMALL) {

               double dMin = dZI - sqrt( dSqRad - dP2 * dP2) ;
               double dMax = dZF + sqrt( dSqRad - dP2 * dP2) ;

               Point3d ptIntMin( dX, dY, dMin) ;
               Vector3d vtNmin = - ( ptIntMin - ptI) + ( ptIntMin - ptI) * vtV1 * vtV1 ;
               vtNmin.Normalize() ;

               Point3d ptIntMax( dX, dY, dMax) ;
               Vector3d vtNmax = - ( ptIntMax - ptF) + ( ptIntMax - ptF) * vtV1 * vtV1 ;
               vtNmax.Normalize() ;

               SubtractIntervals( nGrid, i, j, dMin, dMax, vtNmin, vtNmax, CurrTool.GetToolNum()) ;
            }

           // Parte conica
            else if ( dP1 > dStemHeigth && dP1 < CurrTool.GetHeigth() - EPS_SMALL && abs( dP2) < CurrTool.GetRadius() + 
                    ( CurrTool.GetTipRadius() - CurrTool.GetRadius()) * ( dP1 - dStemHeigth) / CurrTool.GetTipHeigth() - EPS_SMALL) {

               double dr = dP2 ;
               double dMr = CurrTool.GetRadius() + ( CurrTool.GetTipRadius() - CurrTool.GetRadius()) * ( dP1 - dStemHeigth) / CurrTool.GetTipHeigth() ;

               double dMin = dZI - sqrt( dMr * dMr - dr * dr) ;
               double dMax = dZF + sqrt( dMr * dMr - dr * dr) ;

               Point3d ptIntMin( dX, dY, dMin) ;
               Vector3d vtUmin = ( ptIntMin - ptVI) - ( ptIntMin - ptVI) * vtV * vtV ;
               vtUmin.Normalize() ;
               Vector3d vtNmin = dDeltaRad * vtV - CurrTool.GetTipHeigth() * vtUmin ;
               vtNmin.Normalize() ;

               Point3d ptIntMax( dX, dY, dMax) ;
               Vector3d vtUmax = ( ptIntMax - ptVF) - ( ptIntMax - ptVF) * vtV * vtV ;
               vtUmax.Normalize() ;
               Vector3d vtNmax = dDeltaRad * vtV - CurrTool.GetTipHeigth() * vtUmax ;
               vtNmax.Normalize() ;

               SubtractIntervals( nGrid, i, j, dMin, dMax, vtNmin, vtNmax, CurrTool.GetToolNum()) ;
            }
         }
      }
   }

   else {
     // Riferimento sul piano
      Vector3d vtV1 = - vtToolDir ;
      Vector3d vtV2( vtMove.x, vtMove.y, 0) ;
      vtV2.Normalize() ;
     // Vettore per individuare i punti di tangenza 
     // dei piani con il cilindro
      Vector3d vtCross = vtToolDir ^ vtMove ;
      if ( vtCross.z < 0)  vtCross = - vtCross ;
     // Punti di tangenza
      Point3d ptUp = ptI - vtToolDir * dStemHeigth + CurrTool.GetRadius() * vtCross ;
      Point3d ptDw = ptI - vtToolDir * dStemHeigth - CurrTool.GetRadius() * vtCross ;
      double dSmallLength = CurrTool.GetRadius() * vtCross.LenXY() ;
     // Punti di tangenza in corrispondenza della punta 
      Point3d ptTipUp = ptI - vtToolDir * CurrTool.GetHeigth() + CurrTool.GetTipRadius() * vtCross ;
      Point3d ptTipDw = ptI - vtToolDir * CurrTool.GetHeigth() - CurrTool.GetTipRadius() * vtCross ;

      Vector3d vtUpTan = ptTipUp - ptUp ;
      Vector3d vtDwTan = ptTipDw - ptDw ;

      Vector3d vtUpTanXY( vtUpTan.x, vtUpTan.y, 0) ;
      double dDeltaSmallAbs =  abs( vtUpTanXY * vtV2) ;
      double dDeltaSmall = ( CurrTool.GetRadius() > CurrTool.GetTipRadius() ? dDeltaSmallAbs : - dDeltaSmallAbs) ;

      vtUpTan.Normalize() ;
      vtDwTan.Normalize() ;

      Vector3d vtUpCross = vtMove ^ vtUpTan ;
      Vector3d vtDwCross = - vtMove ^ vtDwTan ;

      if ( vtUpCross.z > 0) 
         vtUpCross = - vtUpCross ;

      if ( vtDwCross.z < 0) 
         vtDwCross = - vtDwCross ;

     // Descrizione piani tangenti al cono
      Vector3d vtR0Up = ptUp - ORIG ;
      Vector3d vtR0Dw = ptDw - ORIG ;

      double dDotUp = vtR0Up * vtUpCross ;
      double dDotDw = vtR0Dw * vtDwCross ;

      double dMin, dMax ;
      Vector3d vtNmin, vtNmax ;
      Point3d ptInt ;

      for ( int i = nStartI ; i <= nEndI ; ++ i) {
         for ( 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 - ptIxy ;

            double dP1 = vtC * vtV1 ;
            double dP2 = vtC * vtV2 ;

           // Parte cilindrica
            if ( dP1 > EPS_SMALL && dP1 <= dStemHeigth &&  // vedere se conviene dP1 < dStemHeigth + eps oppure  dP1 <= dStemHeigth
                 dP2 > - CurrTool.GetRadius() && dP2 < dLengthPathXY + CurrTool.GetRadius()) {

              // Massimi
               if ( dP2 < - dSmallLength) {
                  dMax = dZI + sqrt( dSqRad - dP2 * dP2) ;
                  ptInt.Set( dX, dY, dMax) ;
                  vtNmax = - ( ptInt - ptI) + ( ptInt - ptI) * vtV1 * vtV1 ; 
                  vtNmax.Normalize() ;
               }  

               else if ( dP2 < dLengthPathXY - dSmallLength) {
                  dMax = ptUp.z + dDeltaZ * ( dP2 + dSmallLength) / dLengthPathXY ;
                  vtNmax = - vtCross ;
               }

               else {
                  dMax = dZF + sqrt( dSqRad - ( dP2 - dLengthPathXY) * ( dP2 - dLengthPathXY)) ;
                  ptInt.Set( dX, dY, dMax) ;
                  vtNmax = - ( ptInt - ptF) + ( ptInt - ptF) * vtV1 * vtV1 ;
                  vtNmax.Normalize() ;
               }

              // Minimi 
               if ( dP2 < dSmallLength) {
                  dMin = dZI - sqrt( dSqRad - dP2 * dP2) ;
                  ptInt.Set( dX, dY, dMin) ;
                  vtNmin = - ( ptInt - ptI) + ( ptInt - ptI) * vtV1 * vtV1 ;
                  vtNmin.Normalize() ;
               }

               else if ( dP2 < dLengthPathXY + dSmallLength) {
                  dMin = ptDw.z + dDeltaZ * ( dP2 - dSmallLength) / dLengthPathXY ;
                  vtNmin = vtCross ;
               }

               else {
                  dMin = dZF - sqrt( dSqRad - ( dP2 - dLengthPathXY) * ( dP2 - dLengthPathXY)) ;
                  ptInt.Set( dX, dY, dMin) ;
                  vtNmin = - ( ptInt - ptF) + ( ptInt - ptF) * vtV1 * vtV1 ;
                  vtNmin.Normalize() ;
               }

               SubtractIntervals( nGrid, i, j, dMin, dMax, vtNmin, vtNmax, CurrTool.GetToolNum()) ;
            }

           // Parte conica
            else if ( dP1 > dStemHeigth && dP1 < CurrTool.GetHeigth() - EPS_SMALL && 
                      dP2 > - CurrTool.GetRadius() - ( CurrTool.GetTipRadius() - CurrTool.GetRadius()) * ( dP1 - dStemHeigth) / CurrTool.GetTipHeigth() &&
                      dP2 <   CurrTool.GetRadius() + dLengthPathXY + 
                    ( CurrTool.GetTipRadius() - CurrTool.GetRadius()) * ( dP1 - dStemHeigth) / CurrTool.GetTipHeigth()) {

              // Massimi
               if ( dP2 < - dSmallLength + dDeltaSmall * ( dP1 - dStemHeigth) / CurrTool.GetTipHeigth()) {

                  double dr = dP2 ;
                  double dMr = CurrTool.GetRadius() + ( CurrTool.GetTipRadius() - CurrTool.GetRadius()) * ( dP1 - dStemHeigth) / CurrTool.GetTipHeigth() ;

                  dMax = dZI + sqrt( dMr * dMr - dr * dr) ;

                  ptInt.Set( dX, dY, dMax) ;

                  Vector3d vtU = ( ptInt - ptVI) - ( ptInt - ptVI) * vtV * vtV ;
                  vtU.Normalize() ;
                  vtNmax = dDeltaRad * vtV - CurrTool.GetTipHeigth() * vtU ;
                  vtNmax.Normalize() ;
               }

               else if ( dP2 < dLengthPathXY - dSmallLength + 
                         dDeltaSmall * ( dP1 - dStemHeigth) / CurrTool.GetTipHeigth()) {
                  dMax = ( dDotUp - dX * vtUpCross.x - dY * vtUpCross.y) / vtUpCross.z ;
                  vtNmax = vtUpCross ;
               }

               else {
                  double dr = dP2 - dLengthPathXY ; 
                  double dMr = CurrTool.GetRadius() + ( CurrTool.GetTipRadius() - CurrTool.GetRadius()) * ( dP1 - dStemHeigth) / CurrTool.GetTipHeigth() ;

                  dMax = dZF + sqrt( dMr * dMr - dr * dr) ;

                  ptInt.Set( dX, dY, dMax) ;

                  Vector3d vtU = ( ptInt - ptVF) - ( ptInt - ptVF) * vtV * vtV ;
                  vtU.Normalize() ;
                  vtNmax = dDeltaRad * vtV - CurrTool.GetTipHeigth() * vtU ;
                  vtNmax.Normalize() ;
               }

              // Minimi
               if ( dP2 < dSmallLength - dDeltaSmall * ( dP1 - dStemHeigth) / CurrTool.GetTipHeigth()) {

                  double dr = dP2 ;
                  double dMr = CurrTool.GetRadius() + ( CurrTool.GetTipRadius() - CurrTool.GetRadius()) * ( dP1 - dStemHeigth) / CurrTool.GetTipHeigth() ;

                  dMin = dZI - sqrt( dMr * dMr - dr * dr) ;

                  ptInt.Set( dX, dY, dMin) ;

                  Vector3d vtU = ( ptInt - ptVI) - ( ptInt - ptVI) * vtV * vtV ;
                  vtU.Normalize() ;
                  vtNmin = dDeltaRad * vtV - CurrTool.GetTipHeigth() * vtU ;
                  vtNmin.Normalize() ;
               }

               else if ( dP2 < dLengthPathXY + dSmallLength - dDeltaSmall * ( dP1 - dStemHeigth) / CurrTool.GetTipHeigth()) {
                  dMin = ( dDotDw - dX * vtDwCross.x - dY * vtDwCross.y) / vtDwCross.z ;
                  vtNmin = vtDwCross ;
               }

               else {

                  double dr = dP2 - dLengthPathXY ; 
                  double dMr = CurrTool.GetRadius() + ( CurrTool.GetTipRadius() - CurrTool.GetRadius()) * ( dP1 - dStemHeigth) / CurrTool.GetTipHeigth() ;

                  dMin = dZF - sqrt( dMr * dMr - dr * dr) ;

                  ptInt.Set( dX, dY, dMin) ;

                  Vector3d vtU = ( ptInt - ptVF) - ( ptInt - ptVF) * vtV * vtV ;
                  vtU.Normalize() ;
                  vtNmin = dDeltaRad * vtV - CurrTool.GetTipHeigth() * vtU ;
                  vtNmin.Normalize() ;
               }

               SubtractIntervals( nGrid, i, j, dMin, dMax, vtNmin, vtNmax, CurrTool.GetToolNum()) ;
            }
         }
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::Conus_XYMilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir)
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

   double dStemHeigth = CurrTool.GetHeigth() - CurrTool.GetTipHeigth() ;

   CompCyl_Milling( nGrid, ptS, ptE, vtToolDir, dStemHeigth, CurrTool.GetRadius(), false, true, CurrTool.GetToolNum()) ;

   if ( CurrTool.GetTipRadius() < CurrTool.GetRadius()) {
      Point3d ptSTip = ptS - dStemHeigth * vtToolDir ;
      Point3d ptETip = ptE - dStemHeigth * vtToolDir ;
      CompConus_Milling( nGrid, ptSTip, ptETip, vtToolDir,
                         CurrTool.GetTipHeigth(), CurrTool.GetRadius(), CurrTool.GetTipRadius(),
                         true, false, V_NULL, V_NULL, CurrTool.GetToolNum()) ;
   }

   else {
      Point3d ptSTip = ptS - CurrTool.GetHeigth() * vtToolDir ;
      Point3d ptETip = ptE - CurrTool.GetHeigth() * vtToolDir ;
      CompConus_Milling( nGrid, ptSTip, ptETip, - vtToolDir,
                         CurrTool.GetTipHeigth(), CurrTool.GetTipRadius(), CurrTool.GetRadius(),
                         false, true, V_NULL, V_NULL, CurrTool.GetToolNum()) ;
   }

   return true ;
}

// ---------- VERSORE UTENSILE CON ORIENTAZIONE GENERICA ---------------------

// ---------- Cilindro e sfera -----------------------------------------------
//----------------------------------------------------------------------------
bool
VolZmap::CylBall_Drilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir)
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

  // Altezza cilindro
   double dStemHeigth = CurrTool.GetHeigth() - CurrTool.GetTipHeigth() ;
  // Sottraggo cilindro
   CompCyl_Drilling( nGrid, ptS, ptE, vtToolDir, dStemHeigth, CurrTool.GetRadius(), false, false, CurrTool.GetToolNum()) ;
  // Se è sfera la sottraggo 
   if ( CurrTool.GetType() == Tool::BALLMILL) {
      Point3d ptSBall = ptS - dStemHeigth * vtToolDir ;
      Point3d ptEBall = ptE - dStemHeigth * vtToolDir ;
      CompBall_Milling( nGrid, ptSBall, ptEBall, CurrTool.GetRadius(), CurrTool.GetToolNum()) ;
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::CylBall_Milling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir)
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

  // Altezza cilindro
   double dStemHeigth = CurrTool.GetHeigth() - CurrTool.GetTipHeigth() ;
  // Sottraggo cilindro
   CompCyl_Milling( nGrid, ptS, ptE, vtToolDir, dStemHeigth, CurrTool.GetRadius(), false, false, CurrTool.GetToolNum()) ;
  // Se è sfera la sottraggo
   if ( CurrTool.GetType() == Tool::BALLMILL) {
      Point3d ptSBall = ptS - dStemHeigth * vtToolDir ;
      Point3d ptEBall = ptE - dStemHeigth * vtToolDir ;
      CompBall_Milling( nGrid, ptSBall, ptEBall, CurrTool.GetRadius(), CurrTool.GetToolNum()) ;
   }
   return true ;
}

// ---------- Coni -----------------------------------------------------------
//----------------------------------------------------------------------------
bool
VolZmap::Conus_Drilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir)
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

   double dStemHeigth = CurrTool.GetHeigth() - CurrTool.GetTipHeigth() ;

   CompCyl_Drilling( nGrid, ptS, ptE, vtToolDir, dStemHeigth, CurrTool.GetRadius(), false, true, CurrTool.GetToolNum()) ;

  // Trapano
   if ( CurrTool.GetTipRadius() < CurrTool.GetRadius()) {
      Point3d ptSCone = ptS - dStemHeigth * vtToolDir ;
      Point3d ptECone = ptE - dStemHeigth * vtToolDir ;
      CompConus_Drilling( nGrid, ptSCone, ptECone, vtToolDir,
                          CurrTool.GetTipHeigth(), CurrTool.GetRadius(), CurrTool.GetTipRadius(),
                          true, false, V_NULL, V_NULL, CurrTool.GetToolNum()) ;
   }

   else {
      Point3d ptSCone = ptS - CurrTool.GetHeigth() * vtToolDir ;
      Point3d ptECone = ptE - CurrTool.GetHeigth() * vtToolDir ;
      CompConus_Drilling( nGrid, ptSCone, ptECone, - vtToolDir,
                          CurrTool.GetTipHeigth(), CurrTool.GetTipRadius(), CurrTool.GetRadius(),
                          false, true, V_NULL, V_NULL, CurrTool.GetToolNum()) ;
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::Conus_Milling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir)
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

   double dStemHeigth = CurrTool.GetHeigth() - CurrTool.GetTipHeigth() ;

   CompCyl_Milling( nGrid, ptS, ptE, vtToolDir, dStemHeigth, CurrTool.GetRadius(), false, true, CurrTool.GetToolNum()) ;

  // Trapano
   if ( CurrTool.GetTipRadius() < CurrTool.GetRadius()) {
      Point3d ptSBall = ptS - dStemHeigth * vtToolDir ;
      Point3d ptEBall = ptE - dStemHeigth * vtToolDir ;
      CompConus_Milling( nGrid, ptSBall, ptEBall, vtToolDir, CurrTool.GetTipHeigth(),
                         CurrTool.GetRadius(), CurrTool.GetTipRadius(), true, false, V_NULL, V_NULL, CurrTool.GetToolNum()) ;  
   }

   else {
      Point3d ptSBall = ptS - CurrTool.GetHeigth() * vtToolDir ;
      Point3d ptEBall = ptE - CurrTool.GetHeigth() * vtToolDir ;
      CompConus_Milling( nGrid, ptSBall, ptEBall, - vtToolDir, CurrTool.GetTipHeigth(),
                         CurrTool.GetTipRadius(), CurrTool.GetRadius(), false, true, V_NULL, V_NULL, CurrTool.GetToolNum()) ;
   }

   return true ;
}

// ---------- Mortasatrice ---------------------------------------------------
//----------------------------------------------------------------------------
bool
VolZmap::Mrt_Drilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, const Vector3d& vtAux)
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

  // Proiezione della traiettoria sul piano dei movimenti possibili
   Point3d ptEOnP = ptS + ( ptE - ptS) * vtToolDir * vtToolDir ;

  // Scompongo la mortasatrice in solidi semplici

  // Parallelepipedo di base
   double dLenX = CurrTool.GetMrtChsWidth() ;
   double dLenY = CurrTool.GetMrtChsThickness() ;
   double dLenZ = CurrTool.GetHeigth() - CurrTool.GetCornRadius() ;
   CompPar_Drilling( nGrid, dLenX, dLenY, dLenZ, ptS, ptEOnP, vtToolDir, vtAux, CurrTool.GetToolNum()) ;
  
  // Se la punta è di tipo bull-nose
   if ( abs( CurrTool.GetMrtChsWidth() - 2 * CurrTool.GetCornRadius()) > EPS_SMALL) {

     // Parallelepipedo di punta 
      Point3d ptTipS = ptS - dLenZ * vtToolDir ; 
      Point3d ptTipE = ptEOnP - dLenZ * vtToolDir ;

      dLenX = abs( CurrTool.GetMrtChsWidth() - 2 * CurrTool.GetCornRadius()) ;
      dLenZ = CurrTool.GetCornRadius() ;

      CompPar_Drilling( nGrid, dLenX, dLenY, dLenZ, ptTipS, ptTipE, vtToolDir, vtAux, CurrTool.GetToolNum()) ;

      Vector3d vtVOnP = vtToolDir ^ vtAux ;

     // Cilindri 
      Point3d ptSminus = ptTipS - ( 0.5 * dLenX) * vtVOnP + 0.5 * dLenY * vtAux ;
      Point3d ptEminus = ptTipE - ( 0.5 * dLenX) * vtVOnP + 0.5 * dLenY * vtAux ;  
      Point3d ptSplus = ptTipS + ( 0.5 * dLenX) * vtVOnP + 0.5 * dLenY * vtAux ;
      Point3d ptEplus = ptTipE + ( 0.5 * dLenX) * vtVOnP + 0.5 * dLenY * vtAux ;  

      CompCyl_Milling( nGrid, ptSminus, ptEminus, vtAux, dLenY, CurrTool.GetCornRadius(), false, false, CurrTool.GetToolNum()) ;
      CompCyl_Milling( nGrid, ptSplus, ptEplus, vtAux, dLenY, CurrTool.GetCornRadius(), false, false, CurrTool.GetToolNum()) ;
   }

  // se la punta è di tipo sfera
   else {
     // Cilindro
      Point3d ptCylS = ptS - dLenZ * vtToolDir + 0.5 * dLenY * vtAux ;
      Point3d ptCylE = ptEOnP - dLenZ * vtToolDir + 0.5 * dLenY * vtAux ;
      CompCyl_Milling( nGrid, ptCylS, ptCylE, vtAux, dLenY, CurrTool.GetCornRadius(), false, false, CurrTool.GetToolNum()) ;
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::Mrt_Milling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, const Vector3d& vtAux)
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

  // Scompongo la mortasatrice in solidi semplici

  // Parallelepipedo di base
   double dLenX = CurrTool.GetMrtChsWidth() ;
   double dLenY = CurrTool.GetMrtChsThickness() ;
   double dLenZ = CurrTool.GetHeigth() - CurrTool.GetCornRadius() ;

   Point3d ptBasePS = ptS - 0.5 * dLenZ * vtToolDir + 0.5 * dLenY * vtAux ;
   Point3d ptBasePE = ptE - 0.5 * dLenZ * vtToolDir + 0.5 * dLenY * vtAux ;

   CompPar_Milling( nGrid, dLenX, dLenZ, dLenY, ptBasePS, ptBasePE, vtAux, vtToolDir, CurrTool.GetToolNum()) ;
  // Se la punta è di tipo bull-nose
   if ( abs( CurrTool.GetMrtChsWidth() - 2 * CurrTool.GetCornRadius()) > EPS_SMALL) {
   
     // Parallelepipedo di punta 
      Point3d ptTipPS = ptBasePS - 0.5 * dLenZ * vtToolDir ;
      Point3d ptTipPE = ptBasePE - 0.5 * dLenZ * vtToolDir ;

      dLenX = abs( CurrTool.GetMrtChsWidth() - 2 * CurrTool.GetCornRadius()) ; 
      dLenZ = CurrTool.GetCornRadius() ;

      ptTipPS -= 0.5 * dLenZ * vtToolDir ;
      ptTipPE -= 0.5 * dLenZ * vtToolDir ;

      CompPar_Milling( nGrid, dLenX, dLenZ, dLenY, ptTipPS, ptTipPE, vtAux, vtToolDir, CurrTool.GetToolNum()) ;
     // Cilindri 
     // Vettore nel piano
      Vector3d vtVOnP = vtToolDir ^ vtAux ;
      if ( vtVOnP * ( ptE - ptS) < 0)
         vtVOnP *= -1 ;
      Point3d ptSminus = ptTipPS + 0.5 * dLenZ * vtToolDir - ( 0.5 * dLenX) * vtVOnP ;
      Point3d ptEminus = ptTipPE + 0.5 * dLenZ * vtToolDir - ( 0.5 * dLenX) * vtVOnP ;
      CompCyl_Milling( nGrid, ptSminus, ptEminus, vtAux, dLenY, CurrTool.GetCornRadius(), false, false, CurrTool.GetToolNum()) ;
      Point3d ptSplus = ptTipPS + 0.5 * dLenZ * vtToolDir + ( 0.5 * dLenX) * vtVOnP ;
      Point3d ptEplus = ptTipPE + 0.5 * dLenZ * vtToolDir + ( 0.5 * dLenX) * vtVOnP ;
      CompCyl_Milling( nGrid, ptSplus, ptEplus, vtAux, dLenY, CurrTool.GetCornRadius(), false, false, CurrTool.GetToolNum()) ;
   }

  // se la punta è di tipo sfera
   else {
     // Cilindro
      Point3d ptCylS = ptBasePS - 0.5 * dLenZ * vtToolDir ;
      Point3d ptCylE = ptBasePE - 0.5 * dLenZ * vtToolDir ;
      CompCyl_Milling( nGrid, ptCylS, ptCylE, vtAux, dLenY, CurrTool.GetCornRadius(), false, false, CurrTool.GetToolNum()) ;
   }

   return true ;
}

// ---------- Chisel ---------------------------------------------------------
//----------------------------------------------------------------------------
bool
VolZmap::Chs_Drilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, const Vector3d& vtAux)
{  
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

  // Proiezione della traiettoria sulla varietà del movimento
   Point3d ptProjE = ptS + ( ptE - ptS) * vtToolDir * vtToolDir ;

   CompPar_Drilling( nGrid, CurrTool.GetMrtChsWidth(), CurrTool.GetMrtChsThickness(), CurrTool.GetHeigth(),
                     ptS, ptProjE, vtToolDir, vtAux, CurrTool.GetToolNum()) ;

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::Chs_Milling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, const Vector3d& vtAux)
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

  // Proiezione della traiettoria sul piano dei movimenti possibili
   Vector3d vtMoveOnP = ( ptE - ptS) - ( ptE - ptS) * vtToolDir * vtToolDir ;
   Point3d ptProjE = ptS + vtMoveOnP ;

   CompPar_Milling( nGrid, CurrTool.GetMrtChsWidth(), CurrTool.GetMrtChsThickness(), CurrTool.GetHeigth(),
                    ptS, ptProjE, vtToolDir, vtAux, CurrTool.GetToolNum()) ;

   return true ;
}

// ---------- Utensile generico ----------------------------------------------
//----------------------------------------------------------------------------
bool
VolZmap::GenTool_Drilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir)
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

  // Descrizione geometrica del moto  
   Point3d ptI = ptS ;
   Point3d ptF = ptE ;
   Vector3d vtMove = ptE - ptS ;
  // Vettore delle normali agli archi
   const VCT3DVECTOR& vArcNorm = CurrTool.GetArcNormalVec() ;
  // Poinché l'asse utensile è parallelo all'asse Z, definisco un sistema di
  // riferimento ad hoc in cui le normali agli archi giacciano nel piano XZ.
   Frame3d frNormFrame ;
   frNormFrame.Set( ORIG, X_AX, -Z_AX, Y_AX) ;
  // Ciclo sulle curve del profilo
   const CurveComposite& ToolProfile = CurrTool.GetApproxOutline() ;
   int i = - 1 ;
   const ICurve* pPrevCurve = nullptr ;
   const ICurve* pCurve = ToolProfile.GetCurve( ++ i) ;
   while ( pCurve != nullptr) {

      double dHeight = 0 ;

     // Se segmento
      if ( pCurve->GetType() == CRV_LINE) {
        // Recupero gli estremi
         const ICurveLine* pLine = GetCurveLine( pCurve) ;
         Point3d ptStart = pLine->GetStart() ;
         Point3d ptEnd = pLine->GetEnd() ;
         int nNormNum = pLine->GetTempProp();
         Vector3d vtNormSt, vtNormEn;
         if ( nNormNum != 0) {
            vtNormSt = vArcNorm[nNormNum - 1] ;
            vtNormEn = vArcNorm[nNormNum] ;
            vtNormSt.ToLoc(frNormFrame);
            vtNormEn.ToLoc(frNormFrame);
         }
        // Ne determino l'altezza
         dHeight = abs( ptStart.y - ptEnd.y) ;
         if ( dHeight > EPS_SMALL) {
           // Verifiche curva precedente per eventuale tappo sopra
            bool bTapT = false ;
            if ( pPrevCurve != nullptr  &&  pPrevCurve->GetType() == CRV_LINE) {
               const ICurveLine* pOthLine = GetCurveLine( pPrevCurve) ;
               Point3d ptOthStart = pOthLine->GetStart() ;
               Point3d ptOthEnd = pOthLine->GetEnd() ;  
               if ( abs( ptOthStart.y - ptOthEnd.y) < EPS_SMALL  &&  ptOthStart.x < ptOthEnd.x)
                  bTapT = true ;
            }
           // Verifiche curva successiva per eventuale tappo sotto
            bool bTapB = false ;
            const ICurve* pNextCurve = ToolProfile.GetCurve( ++ i) ;
            if ( pNextCurve != nullptr  &&  pNextCurve->GetType() == CRV_LINE) {
               const ICurveLine* pOthLine = GetCurveLine( pNextCurve) ;
               Point3d ptOthStart = pOthLine->GetStart() ;
               Point3d ptOthEnd = pOthLine->GetEnd() ;  
               if ( abs( ptOthStart.y - ptOthEnd.y) < EPS_SMALL  &&  ptOthStart.x > ptOthEnd.x)
                  bTapB = true ;
            }
           // Se X costante, è un cilindro
            if ( abs( ptStart.x - ptEnd.x) < EPS_SMALL) {
               double dRadius = ptStart.x ;
               if ( dRadius > 10 * EPS_SMALL)
                  CompCyl_Drilling( nGrid, ptI, ptF, vtToolDir, dHeight, dRadius, bTapB, bTapT, CurrTool.GetToolNum()) ;
            }
           // Se X crescente, è un cono con vettore equiverso a quello dell'utensile
            else if ( ptStart.x > ptEnd.x) {
               double dMaxRad = ptStart.x ;
               double dMinRad = ptEnd.x ;
               CompConus_Drilling( nGrid, ptI, ptF, vtToolDir, dHeight, dMaxRad, dMinRad, bTapB, bTapT, vtNormSt, vtNormEn, CurrTool.GetToolNum()) ;
            }
           // Se X decrescente, è un cono con vettore opposto a quello dell'utensile
            else if ( ptStart.x < ptEnd.x) {
               double dMaxRad = ptEnd.x ; 
               double dMinRad = ptStart.x ;
               Point3d ptIn = ptI - vtToolDir * dHeight ;
               Point3d ptFn = ptIn + vtMove ;
               vtNormEn.z *= -1 ;
               vtNormSt.z *= -1 ;
               CompConus_Drilling( nGrid, ptIn, ptFn, - vtToolDir, dHeight, dMaxRad, dMinRad, bTapT, bTapB, vtNormEn, vtNormSt, CurrTool.GetToolNum()) ;
            }
           // Passo alla curva successiva
            pPrevCurve = pCurve ;
            pCurve = pNextCurve ;
         }
         else {
           // Passo alla curva successiva
            pPrevCurve = pCurve ;
            pCurve = ToolProfile.GetCurve( ++ i) ;
         }
      }

     // Se arco
      else if ( pCurve->GetType() == CRV_ARC) { 
        // Recupero estremi, centro e raggio
         const ICurveArc* pArc = GetCurveArc( pCurve) ;
         Point3d ptStart ; pArc->GetStartPoint( ptStart) ;
         Point3d ptEnd ; pArc->GetEndPoint( ptEnd) ;
         Point3d ptCen = pArc->GetCenter() ;
         double dRadius = pArc->GetRadius() ;
        // Determino le posizioni iniziale e finale del centro della sfera
         Point3d ptCenS = ptI - vtToolDir * ( ptStart.y - ptCen.y) ;
         Point3d ptCenE = ptCenS + vtMove ;
        // Eseguo l'asportazione del materiale
         CompBall_Milling( nGrid, ptCenS, ptCenE, dRadius, CurrTool.GetToolNum()) ;
        // aggiorno l'altezza
         dHeight = abs( ptStart.y - ptEnd.y) ;
        // Passo alla curva successiva
         pPrevCurve = pCurve ;
         pCurve = ToolProfile.GetCurve( ++ i) ;
      } 

     // Determino le posizioni iniziale e finale del componente successivo
      ptI = ptI - vtToolDir * dHeight ;
      ptF = ptI + vtMove ;
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::GenTool_Milling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir)
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

  // Descrizione geometrica del moto  
   Point3d ptI = ptS ;
   Point3d ptF = ptE ;
   Vector3d vtMove = ptE - ptS ;
  // Vettore delle normali agli archi
   const VCT3DVECTOR& vArcNorm = CurrTool.GetArcNormalVec() ;
  // Poinché l'asse utensile è parallelo all'asse Z, definisco un sistema di
  // riferimento ad hoc in cui le normali agli archi giacciano nel piano XZ.
   Frame3d frNormFrame ;
   frNormFrame.Set( ORIG, X_AX, -Z_AX, Y_AX) ;
  // Ciclo sulle curve del profilo
   const CurveComposite& ToolProfile = CurrTool.GetApproxOutline() ;
   int i = - 1 ;
   const ICurve* pPrevCurve = nullptr ;
   const ICurve* pCurve = ToolProfile.GetCurve( ++ i) ;
   while ( pCurve != nullptr) {
   
      double dHeight = 0 ;

     // Se segmento
      if ( pCurve->GetType() == CRV_LINE) {
        // Recupero gli estremi
         const ICurveLine* pLine = GetCurveLine( pCurve) ;
         Point3d ptStart = pLine->GetStart() ;
         Point3d ptEnd = pLine->GetEnd() ;
         int nNormNum = pLine->GetTempProp() ;
         Vector3d vtNormSt, vtNormEn ;
         if ( nNormNum != 0) {
            vtNormSt = vArcNorm[nNormNum - 1] ;
            vtNormEn = vArcNorm[nNormNum] ;
            vtNormSt.ToLoc( frNormFrame) ;
            vtNormEn.ToLoc( frNormFrame) ;
         }
        // Ne determino l'altezza
         dHeight = abs( ptStart.y - ptEnd.y) ;
         if ( dHeight > EPS_SMALL) {
           // verifiche curva precedente per eventuale tappo sopra
            bool bTapT = true ;
            if ( pPrevCurve != nullptr  &&  pPrevCurve->GetType() == CRV_LINE) {
               const ICurveLine* pOthLine = GetCurveLine( pPrevCurve) ;
               Point3d ptOthStart = pOthLine->GetStart() ;
               Point3d ptOthEnd = pOthLine->GetEnd() ;  
               if ( abs( ptOthStart.y - ptOthEnd.y) < EPS_SMALL  &&  ptOthStart.x < ptOthEnd.x)
                  bTapT = false ;
            }
           // verifiche curva successiva per eventuale tappo sotto
            bool bTapB = true ;
            const ICurve* pNextCurve = ToolProfile.GetCurve( ++ i) ;
            if ( pNextCurve != nullptr  &&  pNextCurve->GetType() == CRV_LINE) {
               const ICurveLine* pOthLine = GetCurveLine( pNextCurve) ;
               Point3d ptOthStart = pOthLine->GetStart() ;
               Point3d ptOthEnd = pOthLine->GetEnd() ;  
               if ( abs( ptOthStart.y - ptOthEnd.y) < EPS_SMALL  &&  ptOthStart.x > ptOthEnd.x)
                  bTapB = false ;
            }
           // Se X costante, è un cilindro
            if ( abs( ptStart.x - ptEnd.x) < EPS_SMALL) {
               double dRadius = ptStart.x ;
               if ( dRadius > 10 * EPS_SMALL)
                  CompCyl_Milling( nGrid, ptI, ptF, vtToolDir, dHeight, dRadius, bTapB, bTapT, CurrTool.GetToolNum()) ;
            }
           // se altrimenti X decrescente, è un cono con vettore equiverso a quello dell'utensile
            else if ( ptStart.x > ptEnd.x) {
               double dMaxRad = ptStart.x ;
               double dMinRad = ptEnd.x ;
               CompConus_Milling( nGrid, ptI, ptF, vtToolDir, dHeight, dMaxRad, dMinRad,
                                  bTapB, bTapT, vtNormSt, vtNormEn, CurrTool.GetToolNum()) ;
            }
           // altrimenti X crescente, è un cono con vettore opposto a quello dell'utensile
            else {
               double dMaxRad = ptEnd.x ; 
               double dMinRad = ptStart.x ;
               Point3d ptIn = ptI - vtToolDir * dHeight ;
               Point3d ptFn = ptIn + vtMove ;
               vtNormEn.z = -vtNormEn.z ;
               vtNormSt.z = -vtNormSt.z ;
               CompConus_Milling( nGrid, ptIn, ptFn, - vtToolDir, dHeight, dMaxRad, dMinRad,
                                  bTapT, bTapB, vtNormEn, vtNormSt, CurrTool.GetToolNum()) ;
            }  
           // Passo alla curva successiva
            pPrevCurve = pCurve ;
            pCurve = pNextCurve ;
         }
         else {
           // Passo alla curva successiva
            pPrevCurve = pCurve ;
            pCurve = ToolProfile.GetCurve( ++ i) ;
         }
      }

     // Se arco
      else if ( pCurve->GetType() == CRV_ARC) { 
        // Recupero estremi, centro e raggio
         const ICurveArc* pArc = GetCurveArc( pCurve) ;
         Point3d ptStart ; pArc->GetStartPoint( ptStart) ;
         Point3d ptEnd ; pArc->GetEndPoint( ptEnd) ;
         Point3d ptCen = pArc->GetCenter() ;
         double dRadius = pArc->GetRadius() ;
        // Determino le posizioni iniziale e finale del centro della sfera
         Point3d ptCenS = ptI - vtToolDir * ( ptStart.y - ptCen.y) ;
         Point3d ptCenE = ptCenS + vtMove ;
        // Eseguo l'asportazione del materiale
         CompBall_Milling( nGrid, ptCenS, ptCenE, dRadius, CurrTool.GetToolNum()) ;
        // aggiorno l'altezza
         dHeight = abs( ptStart.y - ptEnd.y) ;
        // Passo alla curva successiva
         pPrevCurve = pCurve ;
         pCurve = ToolProfile.GetCurve( ++ i) ;
      } 

     // Determino le posizioni iniziale e finale del componente successivo
      ptI = ptI - vtToolDir * dHeight ;
      ptF = ptI + vtMove ;
   }

   return true ;
}

// ------------------------- SOLIDI ELEMENTARI -----------------------------------------------------------------------------
// Asse di simmetria diretto come l'asse Z: FORATURA 

//----------------------------------------------------------------------------     
bool
VolZmap::CompCyl_ZDrilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, double dHei, double dRad, int nToolNum) 
{
  // Verifica sull'interferenza con lo Zmap
   int nStartI, nStartJ, nEndI, nEndJ ;
   if ( ! TestCompoBBox( nGrid, ptS, ptE, vtToolDir, dRad, dRad, dHei, nStartI, nStartJ, nEndI, nEndJ))
      return true ;

  // Proiezione dei punti sul piano
   Point3d ptSxy( ptS.x, ptS.y, 0) ;

  // Parametri geometrici dell'utensile
   double dSafeSqRad = dRad * dRad - 2 * dRad * EPS_SMALL ;

  // Punte del gambo
   Point3d ptTStemS = ptS - vtToolDir * dHei ;
   Point3d ptTStemE = ptE - vtToolDir * dHei ;

  // Quote estreme del gambo
   double dMinStemZ = min( min( ptS.z, ptTStemS.z), min( ptE.z, ptTStemE.z)) ;
   double dMaxStemZ = max( max( ptS.z, ptTStemS.z), max( ptE.z, ptTStemE.z)) ;

  // Ciclo sui punti
   for ( int i = nStartI ; i <= nEndI ; ++ i) {
      for ( 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 - ptSxy ;

         double dSqLen = vtC.SqLen() ;

        // Se il punto si trova dentro il cerchio taglio
         if ( dSqLen < dSafeSqRad)
            SubtractIntervals( nGrid, i, j, dMinStemZ, dMaxStemZ, Z_AX, -Z_AX, nToolNum) ;
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::CompConus_ZDrilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir,
                              double dHei, double dMaxRad, double dMinRad,
                              const Vector3d& vtArcNormMaxR, const Vector3d& vtArcNormMinR, int nToolNum)
{
  // Verifica sull'interferenza con lo Zmap
   int nStartI, nStartJ, nEndI, nEndJ ;
   if ( ! TestCompoBBox( nGrid, ptS, ptE, vtToolDir, dMaxRad, dMinRad, dHei, nStartI, nStartJ, nEndI, nEndJ))
      return true ;

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

   double dZMin, dZMax ;
   double dAngC = dHei / ( dMaxRad - dMinRad) ; 
   double dSqMinRad = dMinRad * dMinRad ;
   double dSqMaxRad = dMaxRad * dMaxRad ;
   double dSafeSqMaxRad = dSqMaxRad - 2 * dMaxRad * EPS_SMALL ; // Questa variabile è sperimentale: serve per evitare il taglio di un dexel dalla parte cilindrica del volume spazzato dalla traslazione del cono.
   double dDeltaR = dMaxRad - dMinRad ;                         // Per tornare alla versione precedente basta sostituire dSafeSqMaxRad con dSqMaxRad. Per risolvere il problema in modo forse più sicuro, ma  
                                                                // computazionalmente più pesante è sottrarre prima il cilindro con dSafeSqMaxRad e dopo il cono con dSqMaxRad.
  // Studio delle simmetrie
   if ( vtToolDir.z > 0) {
      dZMin = ( ptS.z < ptE.z ? ptS.z - dHei : ptE.z - dHei) ;
      dZMax = ( ptS.z < ptE.z ? ptE.z : ptS.z) ;
   }

   else {
      dZMin = ( ptS.z < ptE.z ? ptS.z : ptE.z) ;
      dZMax = ( ptS.z < ptE.z ? ptE.z + dHei : ptS.z + dHei) ;
   }

   double dL = dMaxRad * dAngC ;

   Point3d ptV = ( vtToolDir * ( ptE - ptS) < 0 ? ptE : ptS) - vtToolDir * dL ;

   double dMin, dMax ;
   Vector3d vtMin, vtMax ;

  // Ciclo sui punti
   for ( int i = nStartI ; i <= nEndI ; ++ i) {
      for ( 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 < dSqMinRad)  {
            SubtractIntervals( nGrid, i, j, dZMin, dZMax, Z_AX, -Z_AX, nToolNum) ;
         }

         else if ( dSqDist < dSafeSqMaxRad) {             // dSafeSqMaxRad è sperimentale 

            double dr = sqrt( dSqDist) ;

            if ( vtToolDir.z > 0) { 
               dMin = dZMin + dAngC * ( dr - dMinRad) ;
               dMax = dZMax ;

               vtMax = - Z_AX ;

               Point3d ptInt( dX, dY, dMin) ;

               Vector3d vtU = ( ptInt - ptV) - ( ptInt - ptV) * vtToolDir * vtToolDir ;
               vtU.Normalize() ;

               if ( vtArcNormMaxR.IsSmall()  ||  vtArcNormMinR.IsSmall()) {
                  vtMin = dDeltaR * vtToolDir - dHei * vtU ;
                  vtMin.Normalize() ;
               }
               else {
                  Vector3d vtOriginalN = ( ( dDeltaR - dr + dMinRad) / dDeltaR) * vtArcNormMinR + ( ( dr - dMinRad) / dDeltaR) * vtArcNormMaxR ;
                  vtOriginalN.Normalize() ;
                  vtMin = - vtOriginalN.z * vtToolDir - vtOriginalN.x * vtU ;
                  vtMin.Normalize() ;
               }
               
               SubtractIntervals( nGrid, i, j, dMin, dMax, vtMin, vtMax, nToolNum) ;
            }

            else {
               dMin = dZMin ;
               dMax = dZMax - dAngC * ( dr - dMinRad) ;

               vtMin = Z_AX ;

               Point3d ptInt( dX, dY, dMax) ;

               Vector3d vtU = ( ptInt - ptV) - ( ptInt - ptV) * vtToolDir * vtToolDir ;
               vtU.Normalize()  ;

               if ( vtArcNormMinR.IsSmall()  ||  vtArcNormMaxR.IsSmall()) {
                  vtMax = dDeltaR * vtToolDir - dHei * vtU ;
                  vtMax.Normalize() ;
               }
               else {
                  Vector3d vtOriginalN = ( ( dDeltaR - dr + dMinRad) / dDeltaR) * vtArcNormMinR + ( ( dr - dMinRad) / dDeltaR) * vtArcNormMaxR ;
                  vtOriginalN.Normalize() ;
                  vtMax = - vtOriginalN.z * vtToolDir - vtOriginalN.x * vtU ;
                  vtMax.Normalize() ;
               }

               SubtractIntervals( nGrid, i, j, dMin, dMax, vtMin, vtMax, nToolNum) ;
            }
         }
      }
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::CompPar_ZDrilling( int nGrid, double dLenX, double dLenY, double dLenZ,
                            const Point3d& ptS, const Point3d& ptE,
                            const Vector3d& vtToolDir, const Vector3d& vtAux, int nToolNum)
{
   Point3d ptMyS = ptS - EPS_SMALL * vtToolDir ;
   Point3d ptMyE = ptE - EPS_SMALL * vtToolDir ;

   dLenX -= ( 2 * EPS_SMALL) ;
   dLenY -= ( 2 * EPS_SMALL) ;
   dLenZ -= ( 2 * EPS_SMALL) ;

  // Controllo sull'interferenza utensile-solido
   int nStartI, nStartJ, nEndI, nEndJ ;
   if ( ! TestParaBBox( nGrid, ptMyS, ptMyE, vtToolDir, vtAux, dLenX, dLenY, dLenZ, nStartI, nStartJ, nEndI, nEndJ))
      return true ;

  // Costruisco un sistema di riferimento nel piano
   Point3d ptOXY( ptMyS.x, ptMyS.y, 0) ;
   Vector3d vtV1 = vtToolDir ^ vtAux ;
   Vector3d vtV2 = vtAux ;

  // Quote estreme del volume asportato
   double dMinZ = min( min( ptMyS.z, ptMyS.z - vtToolDir.z * dLenZ), min( ptMyE.z, ptMyE.z - vtToolDir.z * dLenZ)) ;
   double dMaxZ = max( max( ptMyS.z, ptMyS.z - vtToolDir.z * dLenZ), max( ptMyE.z, ptMyE.z - vtToolDir.z * dLenZ)) ;

   for ( int i = nStartI ; i <= nEndI ; ++ i) {
      for ( 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 vtR = ptC - ptOXY ;

         double dR1 = vtR * vtV1 ;
         double dR2 = vtR * vtV2 ;
         
         if ( abs( dR1) < 0.5 * dLenX  &&  abs( dR2) < 0.5 * dLenY) 
            SubtractIntervals( nGrid, i, j, dMinZ, dMaxZ, Z_AX, - Z_AX, nToolNum) ;
      }
   }

   return true ;
}

// Asse di simmetria diretto come l'asse Z: FRESATURA

//----------------------------------------------------------------------------
bool
VolZmap::CompCyl_ZMilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir,
                           double dHei, double dRad, int nToolNum) 
{
  // Verifica sull'interferenza con lo Zmap
   int nStartI, nStartJ, nEndI, nEndJ ;
   if ( ! TestCompoBBox( nGrid, ptS, ptE, vtToolDir, dRad, dRad, dHei, nStartI, nStartJ, nEndI, nEndJ))
      return true ;

  // Parametri geometrici 
   double dSqRad = dRad * dRad ;

  // Studio delle simmetrie
   Point3d ptI = ( ptS.z < ptE.z ? ptS : ptE) ;
   Point3d ptF = ( ptS.z < ptE.z ? ptE : ptS) ;
   Point3d ptIT = ptI - vtToolDir * dHei ;
   Point3d ptFT = ptF - vtToolDir * dHei ;
   Point3d ptIxy( ptI.x, ptI.y, 0) ;

   Point3d ptIUp( ptI.x, ptI.y, max( ptI.z, ptIT.z)) ;
   Point3d ptIDw( ptI.x, ptI.y, min( ptI.z, ptIT.z)) ;
  
  // Quote iniziali e finali massime e 
  // minime del gambo dell'utensile e DeltaZ 
   double dZMaxI = max( ptI.z, ptIT.z) ;
   double dZMaxF = max( ptF.z, ptFT.z) ;
   double dZMinI = dZMaxI - dHei ;
   double dDeltaZ = dZMaxF - dZMaxI ;

  // Vettori caratterizzanti il moto 
   Vector3d vtMove = ptF - ptI ;
   Vector3d vtMoveXY( vtMove.x, vtMove.y, 0) ;
   double dLenXY = vtMoveXY.LenXY() ;
   vtMove.Normalize() ;

   double dSafeRad = dRad - EPS_SMALL ;

  // Definizione di un sistema di riferimento ad hoc
   Vector3d vtV1, vtV2 ;

  // Se la lunghezza è troppo piccola lo allungo
   if ( dLenXY < EPS_SMALL) 
      vtV1 = ( 1 / dLenXY) * vtMoveXY ;
   else 
      vtV1 = vtMoveXY ;

  // Normalizzo vtV1
   vtV1.Normalize() ;
  // Definisco vtV2
   vtV2 = vtV1 ;  
   vtV2.Rotate( Z_AX, 0, 1) ;

   double dMin, dMax ;
   Vector3d vtMin, vtMax ;

   for ( int i = nStartI ; i <= nEndI ; ++ i) {
      for ( 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 - ptIxy ;

         double dX1 = vtC * vtV1 ;
         double dX2 = vtC * vtV2 ;

         Point3d ptInt ;

        // Se il punto appartiene alla proiezione del volume spazzato valuto massimo e minimo
         if ( ( dX1 > 0 && dX1 < dLenXY && abs( dX2) < dSafeRad)    ||
              ( dX1 - dLenXY) * ( dX1 - dLenXY) + dX2 * dX2 < dSqRad || 
                dX1 * dX1 + dX2 * dX2 < dSqRad) {

            double dX1_0 = sqrt( dSqRad - dX2 * dX2) ;
           // Massimo
            if ( ( dX1 - dLenXY) * ( dX1 - dLenXY) + dX2 * dX2 < dSqRad) {
               dMax = dZMaxF ;
               vtMax = - Z_AX ;
            }

            else {
               dMax = dZMaxI + dDeltaZ * ( dX1 + dX1_0) / dLenXY ;
               if ( abs( vtMove * vtToolDir) < EPS_SMALL) {
                  vtMax = - Z_AX ;
               }
               else {
                  Vector3d vtCirc = dX1_0 * vtV1 - dX2 * vtV2 ;
                  Vector3d vtTan( - vtCirc.y, vtCirc.x, 0) ;
                  Vector3d vtCross = vtTan ^ vtMove ;
                  vtMax = ( vtCross * vtCirc > 0 ? vtCross : - vtCross) ;
                  vtMax.Normalize() ;
               }
            }
           // Minimo
            if ( dX1 * dX1 + dX2 * dX2 < dSqRad) {
               dMin = dZMinI ;
               vtMin = Z_AX ;
            }
            else {
               dMin = dZMinI + dDeltaZ * ( dX1 - dX1_0) / dLenXY ;
               if ( abs( vtMove * vtToolDir) < EPS_SMALL) {
                  vtMin = Z_AX ;
               }
               else {
                  Vector3d vtCirc = - dX1_0 * vtV1 - dX2 * vtV2 ;
                  Vector3d vtTan( - vtCirc.y, vtCirc.x, 0) ;
                  Vector3d vtCross = vtTan ^ vtMove ;
                  vtMin = ( vtCross * vtCirc > 0 ? vtCross : - vtCross) ;
                  vtMin.Normalize() ;
               }
            }
            SubtractIntervals( nGrid, i, j, dMin, dMax, vtMin, vtMax, nToolNum) ;
         }    
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::CompConus_ZMilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir,
                             double dHei, double dMaxRad, double dMinRad,
                             const Vector3d& vtArcNormMaxR, const Vector3d& vtArcNormMinR, int nToolNum)
{
   int nStartI, nStartJ, nEndI, nEndJ ;
   if ( ! TestCompoBBox( nGrid, ptS, ptE, vtToolDir, dMaxRad, dMinRad, dHei, nStartI, nStartJ, nEndI, nEndJ))
      return true ;

   Point3d ptI = ( vtToolDir * ( ptE - ptS) > 0 ? ptS : ptE) ;
   Point3d ptF = ( vtToolDir * ( ptE - ptS) > 0 ? ptE : ptS) ;

   Point3d ptIxy( ptI.x, ptI.y, 0) ;
   Point3d ptFxy( ptF.x, ptF.y, 0) ;

   Vector3d vtMove = ptF - ptI ;
   Vector3d vtMLong = ( vtMove * vtToolDir) * vtToolDir ;
   Vector3d vtMOrt = vtMove - vtMLong ;  double dLOrt = vtMOrt.Len() ;

   Vector3d vtV1 = vtToolDir ;  
   Vector3d vtV2 = vtMOrt ;  vtV2.Normalize() ;
   Vector3d vtV3 = vtV1 ^ vtV2 ;

   double dZI = ptI.z ;
   double dZTI = ptI.z - vtV1.z * dHei ;
   double dDeltaZ = ptF.z - ptI.z ;
   double dDeltaR = dMaxRad - dMinRad ;

   double dTan = dDeltaR / dHei ;
   double dRatio = ( vtMove * vtV1) / ( vtMove * vtV2) ;

   bool bCase = ( dRatio * dTan > 1) ;

   double dCos = dTan * dRatio ;
   double dSin = ( 1 - dCos * dCos > 0 ? sqrt( 1 - dCos * dCos) : 0) ;

   double dDen = sqrt( 1 + dTan * dTan) ;
    
   Point3d ptV = ptI - vtV1 * ( dHei * dMaxRad / dDeltaR) ;
   Point3d ptVF = ptV + vtMove ;

   Vector3d vtNs = - ( dTan / dDen) * vtV1 + ( dCos / dDen) * vtV2 + ( dSin / dDen) * vtV3 ;
   Vector3d vtNd = - ( dTan / dDen) * vtV1 + ( dCos / dDen) * vtV2 - ( dSin / dDen) * vtV3 ;
   Vector3d vtR0 = ptV - ORIG ;

   vtNs.Normalize() ;
   vtNd.Normalize() ;

   double dDots = vtR0 * vtNs ;
   double dDotd = vtR0 * vtNd ;

   double dMin = 0, dMax = 0, dPLim = 0, dMLim = 0 ;
   Vector3d vtMin, vtMax, vtP, vtM ;

   Vector3d vtUmv = vtMove ;
   vtUmv.Normalize() ;

   Point3d ptInt ;

   for ( int i = nStartI ; i <= nEndI ; ++ i) {

      for ( 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 vtCI = ptC - ptIxy ;  double dSqDI = vtCI.SqLenXY() ;
         Vector3d vtCF = ptC - ptFxy ;  double dSqDF = vtCF.SqLenXY() ;
         
         double dIDO = vtCI * vtV3 ;
         double dIDL = vtCI * vtV2 ;
         double dIVarCos = dIDL / sqrt( dSqDI) ;

         double dFDL = vtCF * vtV2 ;
         double dFVarCos = dFDL / sqrt( dSqDF) ;

         vtCI.Normalize() ;
         vtCF.Normalize() ;

         if ( dSqDI < dMaxRad * dMaxRad  ||  dSqDF < dMaxRad * dMaxRad  || 
              ( abs( dIDO) < dMaxRad  &&  dIDL > 0  &&  dIDL < dLOrt)) {

           // Caso dTan * dRatio > 1
            if ( bCase) {
              // Limiti nella direzione positiva di vtV1 
               if ( dSqDF < dMaxRad * dMaxRad) {
                  dPLim = dZI + dDeltaZ ;
                  vtP = - vtV1 ;
               }

               else { 
                  double dIDL_0 = - sqrt( dMaxRad * dMaxRad - dIDO * dIDO) ;
                  dPLim = dZI + ( dIDL - dIDL_0) * dDeltaZ / dLOrt ;
                  if ( abs( vtV1 * vtUmv) < EPS_SMALL) {
                     vtP = - vtV1 ;
                  }
                  else {
                     Vector3d vtCirc = - dIDL_0 * vtV2 - dIDO * vtV3 ;
                     Vector3d vtTan( - vtCirc.y, vtCirc.x, 0) ;
                     Vector3d vtCross = vtTan ^ vtUmv ;
                     vtP = ( vtCross * vtCirc > 0 ? vtCross : - vtCross) ;
                     vtP.Normalize() ;
                  }
               }

              // Limiti nella direzione negativa di vtV1
               if ( dSqDI < dMinRad * dMinRad) {
                  dMLim = dZTI ;
                  vtM = vtToolDir ;
               }

               else if ( dSqDI < dMaxRad * dMaxRad) { 
                  dMLim = dZTI + ( sqrt( dSqDI) - dMinRad) * ( dZI - dZTI) / dDeltaR ;
                  ptInt.Set( dX, dY, dMLim) ;
                  Vector3d vtU = ( ptInt - ptV) - ( ptInt - ptV) * vtToolDir * vtToolDir ;
                  vtU.Normalize() ;
                  if ( vtArcNormMaxR.IsSmall()  ||  vtArcNormMinR.IsSmall()) {
                     vtM = dDeltaR * vtToolDir - dHei * vtU ;
                     vtM.Normalize() ;
                  }
                  else {
                     Vector3d vtOrigMaxR = - vtArcNormMaxR.x * vtCI - vtArcNormMaxR.z * vtV1 ;
                     Vector3d vtOrigMinR = - vtArcNormMinR.x * vtCI - vtArcNormMinR.z * vtV1 ;
                     vtOrigMaxR.Normalize() ;
                     vtOrigMinR.Normalize() ;
                     vtM = ( ( dDeltaR - sqrt( dSqDI) + dMinRad) / dDeltaR) * vtOrigMinR + ( ( sqrt(dSqDI) - dMinRad) / dDeltaR) * vtOrigMaxR ;
                     vtM.Normalize() ;
                  }          
               }

               else {
                  double dIDL_0 = sqrt( dMaxRad * dMaxRad - dIDO * dIDO) ;
                  dMLim = dZI + ( dIDL - dIDL_0) * dDeltaZ / dLOrt ;
                  if ( abs( vtUmv * vtV1) < EPS_SMALL) {
                     vtM = vtV1 ;
                  }
                  else {
                     Vector3d vtCirc = - dIDL_0 * vtV2 - dIDO * vtV3 ;
                     Vector3d vtTan( - vtCirc.y, vtCirc.x, 0) ;
                     Vector3d vtCross = vtTan ^ vtUmv ;
                     vtM = ( vtCross * vtCirc > 0 ? vtCross : - vtCross) ;
                     vtM.Normalize() ;
                  }
               }
            }

            else {
              // Limiti nella direzione positiva di vtV1 
               if ( dSqDF < dMaxRad * dMaxRad) { 
                  dPLim = dZI + dDeltaZ ;
                  vtP = - vtV1 ;
               }

               else {
                  double dIDL_0 = - sqrt( dMaxRad * dMaxRad - dIDO * dIDO) ;
                  dPLim = dZI + ( dIDL - dIDL_0) * dDeltaZ / dLOrt ; 
                  if ( abs( vtUmv * vtToolDir) < EPS_SMALL) {
                     vtP = - vtV1 ;
                  }
                  else {
                     Vector3d vtCirc = - dIDL_0 * vtV2 - dIDO * vtV3 ;
                     Vector3d vtTan( - vtCirc.y, vtCirc.x, 0) ;
                     Vector3d vtCross = vtTan ^ vtUmv ;
                     vtP = ( vtCross * vtCirc > 0 ? vtCross : - vtCross) ; // vtCross * vtCirc o vtCross * vtMove?
                     vtP.Normalize() ;
                  }
               }

              // Limiti nella direzione negativa di vtV1
               if ( dSqDI < dMinRad * dMinRad) {
                  dMLim = dZTI ;
                  vtM = vtV1 ;
               }

               else if ( dSqDI >= dMinRad * dMinRad && dSqDI < dMaxRad * dMaxRad && dIVarCos < dCos) { 
                  dMLim = dZTI + ( sqrt( dSqDI) - dMinRad) * ( dZI - dZTI) / dDeltaR ;
                  ptInt.Set( dX, dY, dMLim) ;
                  Vector3d vtU = ( ptInt - ptV) - ( ptInt - ptV) * vtToolDir * vtToolDir ;
                  vtU.Normalize() ;
                  if ( vtArcNormMaxR.IsSmall()  ||  vtArcNormMinR.IsSmall()) {
                     vtM = dDeltaR * vtToolDir - dHei * vtU ;
                     vtM.Normalize() ;
                  }
                  else {
                     Vector3d vtOrigMaxR = - vtArcNormMaxR.x * vtCI - vtArcNormMaxR.z * vtV1 ;
                     Vector3d vtOrigMinR = - vtArcNormMinR.x * vtCI - vtArcNormMinR.z * vtV1 ;
                     vtOrigMaxR.Normalize() ;
                     vtOrigMinR.Normalize() ;
                     vtM = ( ( dDeltaR - sqrt( dSqDI) + dMinRad) / dDeltaR) * vtOrigMinR + ( ( sqrt( dSqDI) - dMinRad) / dDeltaR) * vtOrigMaxR ;
                     vtM.Normalize() ;
                  }
               }

               else if ( dSqDI >= dMinRad * dMinRad && dIVarCos >= dCos && dFVarCos < dCos && abs( dIDO) < dMaxRad * dSin) {

                  if ( dIDO > - dMaxRad * dSin && dIDO <= - dMinRad * dSin) { 
                     dMLim = ( dDotd - dX * vtNd.x - dY * vtNd.y) / vtNd.z ; 
                     if ( vtArcNormMaxR.IsSmall()  ||  vtArcNormMinR.IsSmall()) {
                        vtM = - vtNd ;
                     }
                     else {
                        double dLatMin = dMinRad * dSin ;
                        double dDeltaLat = dDeltaR * dSin ;
                        Vector3d vtRadial = dCos * vtV2 - dSin * vtV3 ;
                        vtRadial.Normalize() ;
                        Vector3d vtOrigMaxR = - vtArcNormMaxR.x * vtRadial - vtArcNormMaxR.z * vtV1 ;
                        Vector3d vtOrigMinR = - vtArcNormMinR.x * vtRadial - vtArcNormMinR.z * vtV1 ;
                        vtOrigMaxR.Normalize() ;
                        vtOrigMinR.Normalize() ;                    
                        vtM = ( ( dDeltaLat - abs( dIDO) + dLatMin) / dDeltaLat) * vtOrigMinR + ( ( abs( dIDO) - dLatMin) / dDeltaLat) * vtOrigMaxR ;
                        vtM.Normalize() ;
                     }
                  }

                  else if ( dIDO > - dMinRad * dSin && dIDO < dMinRad * dSin) {
                     double dIDL_0 = sqrt( dMinRad * dMinRad - dIDO * dIDO) ;
                     dMLim = dZTI + ( dIDL - dIDL_0) * dDeltaZ / dLOrt ;
                     if ( abs( vtUmv * vtV1) < EPS_SMALL) {
                        vtM = vtV1 ;
                     }
                     else {
                        Vector3d vtCirc = - dIDL_0 * vtV2 - dIDO * vtV3 ;
                        Vector3d vtTan( - vtCirc.y, vtCirc.x, 0) ;
                        Vector3d vtCross = vtTan ^ vtUmv ;
                        double dDotCrossCirc = vtCross * vtCirc ;
                        vtM = ( dDotCrossCirc > 0 ? vtCross : - vtCross) ;
                        vtM.Normalize() ;
                     }
                  }

                  else if ( dIDO >= dMinRad * dSin && dIDO < dMaxRad * dSin) {
                     dMLim = ( dDots - dX * vtNs.x - dY * vtNs.y) / vtNs.z ;
                     if ( vtArcNormMaxR.IsSmall()  ||  vtArcNormMinR.IsSmall()) {
                        vtM = - vtNs ;
                     }
                     else {
                        double dLatMin = dMinRad * dSin ;
                        double dDeltaLat = dDeltaR * dSin ;
                        Vector3d vtRadial = dCos * vtV2 + dSin * vtV3 ;
                        vtRadial.Normalize() ;
                        Vector3d vtOrigMaxR = - vtArcNormMaxR.x * vtRadial - vtArcNormMaxR.z * vtV1 ;
                        Vector3d vtOrigMinR = - vtArcNormMinR.x * vtRadial - vtArcNormMinR.z * vtV1 ;
                        vtOrigMaxR.Normalize() ;
                        vtOrigMinR.Normalize() ;
                        vtM = ( ( dDeltaLat -  dIDO + dLatMin) / dDeltaLat) * vtOrigMinR + ( ( dIDO - dLatMin) / dDeltaLat) * vtOrigMaxR ;                       
                        vtM.Normalize() ;
                     }
                  }
               }

               else if ( dFVarCos >= dCos) {
            
                  if ( dSqDF < dMinRad * dMinRad) { 

                     double dIDL_0 = sqrt( dMinRad * dMinRad - dIDO * dIDO) ;

                     dMLim = dZTI + ( dIDL - dIDL_0) * dDeltaZ / dLOrt ;

                     if ( abs( vtUmv * vtV1) < EPS_SMALL) {
                        vtM = vtV1 ;
                     }
                     else {
                        Vector3d vtCirc = - dIDL_0 * vtV2 - dIDO * vtV3 ;
                        Vector3d vtTan( - vtCirc.y, vtCirc.x, 0) ;
                        Vector3d vtCross = vtTan ^ vtUmv ;
                        vtM = ( vtCross * vtCirc > 0 ? vtCross : - vtCross) ;
                        vtM.Normalize() ;
                     }
                  }

                  else {
                     dMLim = dZTI + dDeltaZ + ( sqrt( dSqDF) - dMinRad) * ( dZI - dZTI) / dDeltaR ;
                     ptInt.Set( dX, dY, dMLim) ;
                     Vector3d vtU = ( ptInt - ptVF) - ( ptInt - ptVF) * vtToolDir * vtToolDir ;
                     vtU.Normalize() ;
                     if ( vtArcNormMaxR.IsSmall()  ||  vtArcNormMinR.IsSmall()) {
                        vtM = dDeltaR * vtToolDir - dHei * vtU ;
                        vtM.Normalize() ;
                     }
                     else {
                        Vector3d vtOrigMaxR = - vtArcNormMaxR.x * vtCF - vtArcNormMaxR.z * vtV1 ;
                        Vector3d vtOrigMinR = - vtArcNormMinR.x * vtCF - vtArcNormMinR.z * vtV1 ;
                        vtOrigMaxR.Normalize() ;
                        vtOrigMinR.Normalize() ;
                        vtM = ( ( dDeltaR - sqrt( dSqDF) + dMinRad) / dDeltaR) * vtOrigMinR + ( ( sqrt( dSqDF) - dMinRad) / dDeltaR) * vtOrigMaxR ;
                        vtM.Normalize() ;
                     }
                  }
               }

               else { 
                  double dIDL_0 = sqrt( dMaxRad * dMaxRad - dIDO * dIDO) ;
                  dMLim =  dZI + ( dIDL - dIDL_0) * dDeltaZ / dLOrt ;
                  if ( abs( vtUmv * vtV1) < EPS_SMALL) {
                     vtM = vtV1 ;
                  }
                  else {
                     Vector3d vtCirc = - dIDL * vtV2 - dIDO * vtV3 ;
                     Vector3d vtTan( - vtCirc.y, vtCirc.x, 0) ;
                     Vector3d vtCross = vtTan ^ vtMove ;
                     vtM = ( vtCross * vtMove > 0 ? vtCross : - vtCross) ;
                     vtM.Normalize() ;
                  }
               }
            }

            if ( dMLim < dPLim) {
               dMin = dMLim ;
               dMax = dPLim ;
               vtMin = vtM ;
               vtMax = vtP ;
            }

            else {
               dMin = dPLim ;
               dMax = dMLim ;
               vtMin = vtP ;
               vtMax = vtM ;
            }
            
            SubtractIntervals( nGrid, i, j, dMin, dMax, vtMin, vtMax, nToolNum) ;
         }
      }
   }
   return true ;
} 

//----------------------------------------------------------------------------
bool //  E' in realtà  MillingPerp
VolZmap::CompPar_ZMilling( int nGrid, double dLenX, double dLenY, double dLenZ,
                           const Point3d& ptS, const Point3d& ptE,
                           const Vector3d& vtToolDir, const Vector3d& vtAux, int nToolNum)
{  
   Point3d ptMyS = ptS - EPS_SMALL * vtToolDir ;
   Point3d ptMyE = ptE - EPS_SMALL * vtToolDir ;

   dLenX -= ( 2 * EPS_SMALL) ;
   dLenY -= ( 2 * EPS_SMALL) ;
   dLenZ -= ( 2 * EPS_SMALL) ;

  // Controllo sull'interferenza utensile-solido
   int nStartI, nStartJ, nEndI, nEndJ ;
   if ( ! TestParaBBox( nGrid, ptMyS, ptMyE, vtToolDir, vtAux, dLenX, dLenY, dLenZ, nStartI, nStartJ, nEndI, nEndJ))
      return true ;

  // Quote estreme Z 
   double dMinZ = min( min( ptMyS.z, ptMyS.z - vtToolDir.z * dLenZ), min( ptMyE.z, ptMyE.z - vtToolDir.z * dLenZ)) ;
   double dMaxZ = max( max( ptMyS.z, ptMyS.z - vtToolDir.z * dLenZ), max( ptMyE.z, ptMyE.z - vtToolDir.z * dLenZ)) ;

  // Vettore di movimento
   Vector3d vtMove = ptMyE - ptMyS ;

  // Sistemi di riferimento del solido nella posizione iniziale e finale
   Vector3d vtV2 = vtAux ; 
   Vector3d vtV1 = vtV2 ^ vtToolDir ;
   vtV1.Normalize() ;
   Point3d ptSXY( ptMyS.x, ptMyS.y, 0) ;
   Point3d ptEXY( ptMyE.x, ptMyE.y, 0) ;
  
  // Studio del volume asportato durante il moto
   Frame3d MotionFrame ;
   Vector3d vtDiagA = dLenX * vtV1 + dLenY * vtV2 ;
   Vector3d vtDiagB = dLenX * vtV1 - dLenY * vtV2 ;
   Vector3d vtW1 = ( abs( vtDiagA * vtMove) <= abs( vtDiagB * vtMove) ? vtDiagA : vtDiagB) ;
   Point3d ptOrigXY = ptSXY - vtW1 / 2 ;
   double dLenW1 = vtW1.Len() ;
   vtW1 /= dLenW1 ;
   double dDeltaW1 = vtMove * vtW1 ;
   Vector3d vtW2 = vtMove - dDeltaW1 * vtW1 ;
   double dLenW2 = vtW2.Len() ;
   vtW2 /= dLenW2 ;

   for ( int i = nStartI ; i <= nEndI ; ++ i) {
      for ( 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 vtRS = ptC - ptSXY ;
         Vector3d vtRE = ptC - ptEXY ;

         double dRS1 = vtRS * vtV1 ;
         double dRS2 = vtRS * vtV2 ;

         double dRE1 = vtRE * vtV1 ;
         double dRE2 = vtRE * vtV2 ;
       
        // Asportazione materiale nella posizione iniziale
         if ( abs( dRS1) < 0.5 * dLenX  &&  abs( dRS2) < 0.5 * dLenY)
            SubtractIntervals( nGrid, i, j, dMinZ, dMaxZ, Z_AX, - Z_AX, nToolNum) ;

        // Asportazione materiale nella posizione finale
         if ( abs( dRE1) < 0.5 * dLenX  &&  abs( dRE2) < 0.5 * dLenY)
            SubtractIntervals( nGrid, i, j, dMinZ, dMaxZ, Z_AX, - Z_AX, nToolNum) ;

        // Asportazione materiale nel moto
         Vector3d vtR = ptC - ptOrigXY ;
         double dR1 = vtR * vtW1 ;
         double dR2 = vtR * vtW2 ;
         if ( dR2 > 0  &&  dR2 < dLenW2  && 
              dR1 * dLenW2 > dDeltaW1 * dR2  && 
              dR1 * dLenW2 < dLenW1 * dLenW2 + dDeltaW1 * dR2)
            SubtractIntervals( nGrid, i, j, dMinZ, dMaxZ, Z_AX, - Z_AX, nToolNum) ;
      }
   }
   return true ;
}

// Asse di simmetria con orientazione generica: FORATURA
//----------------------------------------------------------------------------
bool 
VolZmap::CompCyl_Drilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir,
                           double dHei, double dRad, bool bTapB, bool bTapT, int nToolNum)
{
  // Verifico che il cilindro con il suo movimento intersechi la griglia
   int nStartI, nEndI, nStartJ, nEndJ ;
   if ( ! TestCompoBBox( nGrid, ptS, ptE, vtToolDir, dRad, dRad, dHei, nStartI, nStartJ, nEndI, nEndJ))
      return true ;

  // Studio delle simmetrie
   Point3d ptI = ( ( ptE - ptS) * vtToolDir > 0 ? ptE : ptS) ;
   Point3d ptF = ( ( ptE - ptS) * vtToolDir > 0 ? ptS - vtToolDir * dHei : ptE - vtToolDir * dHei) ;

  // Altezza cilindro totale altezza + moto
   double dH = ( ptF - ptI).Len() ;

  // Sistema di riferimento del cilindro
   Frame3d CylFrame ; CylFrame.Set( ptF, vtToolDir) ;

   for ( int i = nStartI ; i <= nEndI ; ++ i) {
      for ( int j = nStartJ ; j <= nEndJ ; ++ j) {
         Point3d ptC( ( i + 0.5) * m_dStep, ( j + 0.5) * m_dStep, 0) ;
         Point3d ptInt1, ptInt2 ;
         Vector3d vtN1, vtN2 ;
         if ( IntersLineCylinder( ptC, Z_AX, CylFrame, dH, dRad, bTapB, bTapT, ptInt1, vtN1, ptInt2, vtN2)) {
            SubtractIntervals( nGrid, i, j, ptInt1.z, ptInt2.z, vtN1, vtN2, nToolNum) ;
         }
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
bool 
VolZmap::CompConus_Drilling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, 
                             double dHei, double dMaxRad, double dMinRad, bool bTapB, bool bTapT,
                             const Vector3d& vtArcNormMaxR, const Vector3d& vtArcNormMinR, int nToolNum)
{
   int nStartI, nStartJ, nEndI, nEndJ ;
   if ( ! TestCompoBBox( nGrid, ptS, ptE, vtToolDir, dMaxRad, dMinRad, dHei, nStartI, nStartJ, nEndI, nEndJ))
      return true ;

  // Apertura del cono
   double dDeltaR = dMaxRad - dMinRad ;
   double dTan = dDeltaR / dHei ;
   double dConeMaxH = ( ( dMaxRad * dHei) / dDeltaR) ;
   double dConeMinH = dConeMaxH - dHei ;

   Point3d ptVertex = ( vtToolDir * ( ptE - ptS) > 0 ? ptS : ptE) - vtToolDir * dConeMaxH ;

  // Sistemi di riferimento del cono e del cilindro
   Frame3d ConusFrame ; ConusFrame.Set( ptVertex, vtToolDir) ;
   Frame3d CylFrame = ConusFrame ; CylFrame.Translate( vtToolDir * dConeMaxH) ;

  // L'altezza del cilindro è il movimento
   double dH = ( ptE - ptS).Len() ;

   for ( int i = nStartI ; i <= nEndI ; ++ i) { 
      for ( 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) ;

         Point3d ptInt1, ptInt2 ;
         Vector3d vtN1, vtN2 ;

        // Cilindro
         if ( IntersLineCylinder( ptC, Z_AX, CylFrame, dH, dMaxRad, true, bTapT,
                                  ptInt1, vtN1, ptInt2, vtN2)) {
            SubtractIntervals( nGrid, i, j, ptInt1.z, ptInt2.z, vtN1, vtN2, nToolNum) ;
         }

        // Cono
         if ( IntersLineConus( ptC, Z_AX, ConusFrame, dTan, dConeMinH, dConeMaxH, bTapB, true,
                               ptInt1, vtN1, ptInt2, vtN2)) {
            if ( ! ( vtArcNormMaxR.IsSmall()  ||  vtArcNormMinR.IsSmall())) {
               if ( ! AreSameOrOppositeVectorEpsilon(vtN1, vtToolDir, 0.1 * EPS_SMALL)) {
                  Vector3d vtL1 = ptInt1 - ptVertex ;
                  vtL1 -= ( vtL1 * vtToolDir) * vtToolDir ;
                  double dL1 = vtL1.Len() ;
                  vtL1 /= dL1 ;
                  Vector3d vtOriginalN1 = ( ( dDeltaR - dL1 + dMinRad) / dDeltaR) * vtArcNormMinR + ( ( dL1 - dMinRad) / dDeltaR) * vtArcNormMaxR ;
                  vtOriginalN1.Normalize() ;
                  vtN1 = - vtOriginalN1.z * vtToolDir - vtOriginalN1.x * vtL1 ;
               }
               if ( ! AreSameOrOppositeVectorEpsilon( vtN2, vtToolDir, 0.1 * EPS_SMALL)) {
                  Vector3d vtL2 = ptInt2 - ptVertex ;
                  vtL2 -= ( vtL2 * vtToolDir) * vtToolDir ;
                  double dL2 = vtL2.Len() ;
                  vtL2 /= dL2 ;
                  Vector3d vtOriginalN2 = ( ( dDeltaR - dL2 + dMinRad) / dDeltaR) * vtArcNormMinR + ( ( dL2 - dMinRad) / dDeltaR) * vtArcNormMaxR ;
                  vtOriginalN2.Normalize() ;
                  vtN2 = - vtOriginalN2.z * vtToolDir - vtOriginalN2.x * vtL2 ;
               }
            }
            SubtractIntervals( nGrid, i, j, ptInt1.z, ptInt2.z, vtN1, vtN2, nToolNum) ;
         }
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::CompPar_Drilling( int nGrid, double dLenX, double dLenY, double dLenZ,
                           const Point3d& ptS, const Point3d& ptE,
                           const Vector3d& vtToolDir, const Vector3d& vtAux, int nToolNum)
{
   Point3d ptMyS = ptS - EPS_SMALL * vtToolDir ;
   Point3d ptMyE = ptE - EPS_SMALL * vtToolDir ;

   dLenX -= ( 2 * EPS_SMALL) ; 
   dLenY -= ( 2 * EPS_SMALL) ; 
   dLenZ -= ( 2 * EPS_SMALL) ;

  // Controllo sull'interferenza utensile-solido
   int nStartI, nStartJ, nEndI, nEndJ ;
   if ( ! TestParaBBox( nGrid, ptMyS, ptMyE, vtToolDir, vtAux, dLenX, dLenY, dLenZ, nStartI, nStartJ, nEndI, nEndJ))
      return true ;

  // Costruzione di un sistema di riferimento per il volume
  // asportato dal solido nelle posizioni iniziale e finale.
   Vector3d vtV1 = vtToolDir ;
   Vector3d vtV2 = vtAux ^ vtV1 ;
   Vector3d vtV3 = vtAux ;

   double dLenMove = ( ptMyE - ptMyS).Len() ;
   Point3d ptO = ptMyS - ( ( ptMyE - ptMyS) * vtV1 > 0 ? dLenZ : dLenZ + dLenMove) * vtV1 - ( 0.5 * dLenX) * vtV2 ;
   Frame3d ParaFrame ; ParaFrame.Set( ptO, vtV2, vtV3, vtV1) ;

   for ( int i = nStartI ; i <= nEndI ; ++ i) {
      for ( 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) ;

         Point3d ptInt1, ptInt2 ;
         Vector3d vtN1, vtN2 ;

         if ( IntersLineMyPolyhedron(  ptC, Z_AX, ParaFrame, dLenX, dLenY, dLenZ + dLenMove, 0, ptInt1, vtN1, ptInt2, vtN2)) {
            SubtractIntervals( nGrid, i, j, ptInt1.z, ptInt2.z, vtN1, vtN2, nToolNum) ;
         }
      }
   }

   return true ;
}

// Asse di simmetria con orientazione generica: FRESATURA
//----------------------------------------------------------------------------
bool
VolZmap::CompCyl_Milling( int nGrid, const Point3d& ptS, const Point3d& ptE,
                          const Vector3d& vtToolDir, double dHei, double dRad, bool bTapB, bool bTapT, int nToolNum) 
{
  // Verifica sull'interferenza utensile Zmap
   int nStartI, nStartJ, nEndI, nEndJ ;
   if ( ! TestCompoBBox( nGrid, ptS, ptE, vtToolDir, dRad, dRad, dHei, nStartI, nStartJ, nEndI, nEndJ))
      return true ;

  // Studio delle simmetrie
   Point3d ptI = ( vtToolDir * ( ptE - ptS) > 0 ? ptS : ptE) ;
   Point3d ptF = ( vtToolDir * ( ptE - ptS) > 0 ? ptE : ptS) ;
   Point3d ptITip = ptI - vtToolDir * dHei ;

  // Definizione terne vettoriali e sistemi di riferimento intrinseci al movimento
   Vector3d vtMove = ptF - ptI ;
   Vector3d vtMoveLong = ( vtMove * vtToolDir) * vtToolDir ;
   Vector3d vtMoveOrt = vtMove - vtMoveLong ;

   Vector3d vtV1 = vtToolDir ;
   Vector3d vtV2 = vtMoveOrt ;  vtV2.Normalize() ;
   Vector3d vtV3 = vtV1 ^ vtV2 ;
   Frame3d CylFrame ;  CylFrame.Set( ptITip, vtV2, vtV3, vtV1) ;
 
  // Parametri geometrici di moto e cilindro
   double dLongLen = vtMoveLong.Len() ;
   double dOrtLen = vtMoveOrt.Len() ;

   for ( int i = nStartI ; i <= nEndI ; ++ i) {
      for ( int j = nStartJ ; j <= nEndJ ; ++ j) {
      
         Point3d ptC( ( i + 0.5) * m_dStep, ( j + 0.5) * m_dStep, 0) ;
            
         Point3d ptInt1, ptInt2 ;
         Vector3d vtN1, vtN2 ;

        // Cilindro iniziale
         CylFrame.ChangeOrig( ptITip) ;
         if ( IntersLineCylinder( ptC, Z_AX, CylFrame, dHei, dRad, bTapB, bTapT, ptInt1, vtN1, ptInt2, vtN2)) {
            SubtractIntervals( nGrid, i, j, ptInt1.z, ptInt2.z, vtN1, vtN2, nToolNum) ;
         }

        // Cilindro finale:L'unica differenza rispetto a prima è l'origine 
        // del sistema di riferimento, quindi usiamo lo stesso sistema sommando a ptC
        // il vettore che congiunge le due origini.
         CylFrame.ChangeOrig( ptITip + vtMove) ;
         if ( IntersLineCylinder( ptC, Z_AX, CylFrame, dHei, dRad, bTapB, bTapT, ptInt1, vtN1, ptInt2, vtN2)) {
            SubtractIntervals( nGrid, i, j, ptInt1.z, ptInt2.z, vtN1, vtN2, nToolNum) ;
         }

        // Poliedro interno
         CylFrame.ChangeOrig( ptITip) ;
         if ( IntersLineMyPolyhedron( ptC, Z_AX, CylFrame, dOrtLen, 2 * dRad, dHei,
                                      dLongLen, ptInt1, vtN1, ptInt2, vtN2)) {
            SubtractIntervals( nGrid, i, j, ptInt1.z, ptInt2.z, vtN1, vtN2, nToolNum) ;
         }
        
        // Se movimento non ortogonale all'asse cilindro ellittico di punta
         if ( dLongLen > EPS_SMALL) {
            CylFrame.ChangeOrig( ptITip) ;
            if ( IntersLineEllipticalCylinder( ptC, Z_AX, CylFrame, dRad, dLongLen, dOrtLen, 
                                               bTapB, true, ptInt1, vtN1, ptInt2, vtN2)) {
               SubtractIntervals( nGrid, i, j, ptInt1.z, ptInt2.z, vtN1, vtN2, nToolNum) ;
            }
           
           // Cilindro ellittico di base: L'unica differenza rispetto a prima è l'origine
           // del sistema di riferimento, quindi usiamo lo stesso sistema sommando a ptC
           // il vettore che congiunge le due origini.
            CylFrame.ChangeOrig( ptITip + dHei * vtV1) ;
            if ( IntersLineEllipticalCylinder( ptC, Z_AX, CylFrame, dRad, dLongLen, dOrtLen,
                                               true, bTapT, ptInt1, vtN1, ptInt2, vtN2)) {
               SubtractIntervals( nGrid, i, j, ptInt1.z, ptInt2.z, vtN1, vtN2, nToolNum) ;
            }
         }
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
static Vector3d
AdjustConeNormal( const Point3d& ptInt, const Vector3d& vtN, const Point3d& ptV, const Vector3d& vtToolDir,
                  double dMinRad, double dDeltaR, const Vector3d& vtArcNormMinR, const Vector3d& vtArcNormMaxR)
{
   if ( AreSameOrOppositeVectorEpsilon( vtN, vtToolDir, 0.1 * EPS_SMALL))
      return vtN ;
   Vector3d vtL = ( ptInt - ptV) - (( ptInt - ptV) * vtToolDir) * vtToolDir ;
   double dL = vtL.Len() ;
   vtL /= dL ;
   Vector3d vtOriginalN = ( ( dDeltaR - dL + dMinRad) / dDeltaR) * vtArcNormMinR + ((dL - dMinRad) / dDeltaR) * vtArcNormMaxR ;
   Vector3d vtNewN = - vtOriginalN.z * vtToolDir - vtOriginalN.x * vtL ;
   vtNewN.Normalize() ;
   return vtNewN ;
}

//----------------------------------------------------------------------------
bool
VolZmap::CompConus_Milling( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtToolDir, 
                            double dHei, double dMaxRad, double dMinRad, bool bTapB, bool bTapT,
                            const Vector3d& vtArcNormMaxR, const Vector3d& vtArcNormMinR, int nToolNum)
{
  // Verifico interferenza
   int nStartI, nStartJ, nEndI, nEndJ ;
   if ( ! TestCompoBBox( nGrid, ptS, ptE, vtToolDir, dMaxRad, dMinRad, dHei, nStartI, nStartJ, nEndI, nEndJ))
      return true ;

  // Geometria del cono
   double dDeltaR = dMaxRad - dMinRad ;
  
  // Studio simmetrie
   Point3d ptI = ptS ;
   Point3d ptF = ptE ;
   if ( vtToolDir * ( ptE - ptS) <= 0)
      swap( ptI, ptF) ;

   double dL = ( dMaxRad * dHei) / dDeltaR ;
   double dl = dL - dHei ;

   Point3d ptV = ptI - vtToolDir * dL ;

  // Vettori caratteristici del movimento
   Vector3d vtMove = ptF - ptI ;
   Vector3d vtMvLong = ( vtMove * vtToolDir) * vtToolDir ;
   Vector3d vtMvOrt = vtMove - vtMvLong ;

  // Terna destrorsa e unitaria
   Vector3d vtV1 = vtToolDir ;
   Vector3d vtV2 = vtMvOrt ;  vtV2.Normalize() ;
   Vector3d vtV3 = vtV1 ^ vtV2 ;

  // Sistema di riferimento intrinseco del movimento
   Frame3d ConusFrame ;  ConusFrame.Set( ptV, vtV2, vtV3, vtV1) ;

  // Dimensioni lineari movimento
   double dLongLen = vtMvLong.Len() ;
   double dOrtLen = vtMvOrt.Len() ;

  // Apertura del cono 
   double dTan = dDeltaR / dHei ;
   double dRatio = dLongLen / dOrtLen ;

  // Per costruire piani laterali poliedro interno
   double dCos = dTan * dRatio ;
   double dSin = ( 1 - dCos * dCos > 0 ? sqrt( 1 - dCos * dCos) : 0) ;

  // Dimensioni lineari descriventi il poliedro interno
   double dLenX = dLongLen ;
   double dLenY = dOrtLen ;
   double dLenZ = dSin * dMinRad ;
   double dDeltaX = dHei ;
   double dDeltaY = dCos * dDeltaR ;
   double dDeltaZ = dSin * dDeltaR ;

  // Sistema di riferimento poliedro
   Point3d ptO = ptV + vtV1 * dl + vtV2 * ( dCos * dMinRad) ;
   Frame3d PolyFrame ;
   PolyFrame.Set( ptO, vtV1, vtV2, vtV3) ;

  // Versori piani nel riferimento poliedro ( riferiti al sistema di riferimento) :
  // Sx, Dx
   Vector3d vtNs( - dTan, dCos, dSin) ;
   vtNs.Normalize() ;
   Vector3d vtNd( - dTan, dCos, - dSin) ;
   vtNd.Normalize() ;
  // Iniziale e finale
   Vector3d vtIF( - dDeltaY, dDeltaX, 0) ;
   vtIF.Normalize() ;
  // Up e Down
   Vector3d vtUD( - dLenY, dLenX, 0) ;
   vtUD.Normalize() ;

  // Punti dei piani (sempre espressi nel sistema PolyFrame)
   Point3d ptFacet135( 0, 0, dLenZ) ;
   Point3d ptFacet246( dLenX + dDeltaX, dLenY + dDeltaY, - dLenZ - dDeltaZ) ;

  // Necessità ricalcolo normali (perchè variabili per approx curve)
   bool bRecalNorm = ( ! vtArcNormMaxR.IsSmall()  &&  ! vtArcNormMinR.IsSmall()) ;

   if ( dRatio * dTan <= 1) {
    
      for ( int i = nStartI ; i <= nEndI ; ++ i) {
         for ( int j = nStartJ ; j <= nEndJ ; ++ j) {

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

            Point3d ptInt1, ptInt2 ;
            Vector3d vtN1, vtN2 ;

           // Cono iniziale
            ConusFrame.ChangeOrig( ptV) ;
            if ( IntersLineConus( ptC, Z_AX, ConusFrame, dTan, dl, dL, bTapB, bTapT, ptInt1, vtN1, ptInt2, vtN2)) {
               if ( bRecalNorm) {
                  vtN1 = AdjustConeNormal( ptInt1, vtN1, ptV, vtToolDir, dMinRad, dDeltaR, vtArcNormMinR, vtArcNormMaxR) ;
                  vtN2 = AdjustConeNormal( ptInt2, vtN2, ptV, vtToolDir, dMinRad, dDeltaR, vtArcNormMinR, vtArcNormMaxR) ;
               }
               SubtractIntervals( nGrid, i, j, ptInt1.z, ptInt2.z, vtN1, vtN2, nToolNum) ;
            }

           // Cono finale
            ConusFrame.ChangeOrig( ptV + vtMove) ;
            if ( IntersLineConus( ptC, Z_AX, ConusFrame, dTan, dl, dL, bTapB, bTapT, ptInt1, vtN1, ptInt2, vtN2)) {
               if ( bRecalNorm) {
                  vtN1 = AdjustConeNormal( ptInt1, vtN1, ptV + vtMove, vtToolDir, dMinRad, dDeltaR, vtArcNormMinR, vtArcNormMaxR) ;
                  vtN2 = AdjustConeNormal( ptInt2, vtN2, ptV + vtMove, vtToolDir, dMinRad, dDeltaR, vtArcNormMinR, vtArcNormMaxR) ;
               }
               SubtractIntervals( nGrid, i, j, ptInt1.z, ptInt2.z, vtN1, vtN2, nToolNum) ;
            }
            
           // Solido interno
            Point3d ptPoly = GetToLoc( ptC, PolyFrame) ;
            Vector3d vtPoly = GetToLoc( Z_AX, PolyFrame) ;

           // Intervallo di intersezione (infinito) e normali (nulle)
            bool bValid = true ;
            double dPar1 = -INFINITO ;
            double dPar2 = +INFINITO ;
            vtN1 = V_NULL ;
            vtN2 = V_NULL ;

           // Verifica con facce iniziale e finale
            if ( bValid) {
              // Distanza con segno del punto di riferimento del dexel dal piano delle facce iniziale e finale
               double dDistI = ( ptPoly - ptFacet135) * -vtIF ;
               double dDistF = ( ptPoly - ptFacet246) * vtIF ;
              // Se dexel non parallelo alle facce
               if ( abs( vtPoly * vtIF) > COS_ORTO_ANG_ZERO) {
                 // posizione parametrica delle intersezioni
                  double dParI = -dDistI / ( vtPoly * -vtIF) ;
                  double dParF = -dDistF / ( vtPoly * vtIF) ;
                 // se intervallo tra intersezioni praticamente nullo
                  if ( abs( dParI - dParF) < EPS_ZERO)
                     bValid = false ;
                 // altrimenti
                  else {
                     if ( dParI < dParF) {
                        dPar1 = dParI ;
                        vtN1 = vtIF ;
                        dPar2 = dParF ;
                        vtN2 = -vtIF ;
                     }
                     else {
                        dPar1 = dParF ;
                        vtN1 = -vtIF ;
                        dPar2 = dParI ;
                        vtN2 = vtIF ;
                     }
                  }
               }
              // altrimenti praticamente parallelo
               else {
                 // se esterno ad almeno uno invalida tutto
                  if ( dDistI > 0  ||  dDistF > 0)
                     bValid = false ;
                 // altrimenti non cambia niente
               }
            }

           // Verifica con facce sopra e sotto
            if ( bValid) {
              // Distanza con segno del punto di riferimento del dexel dal piano delle facce sopra e sotto
               double dDistU = ( ptPoly - ptFacet246) * -vtUD ;
               double dDistD = ( ptPoly - ptFacet135) * vtUD ;
              // Se dexel non parallelo alle facce
               if ( abs( vtPoly * vtUD) > COS_ORTO_ANG_ZERO) {
                 // posizione parametrica delle intersezioni
                  double dParU = -dDistU / ( vtPoly * -vtUD) ;
                  double dParD = -dDistD / ( vtPoly * vtUD) ;
                 // se intervallo tra intersezioni praticamente nullo
                  if ( abs( dParU - dParD) < EPS_ZERO)
                     bValid = false ;
                 // altrimenti
                  else {
                     if ( dParU < dParD) {
                        if ( dParD < dPar1 + EPS_ZERO  ||  dParU > dPar2 - EPS_ZERO)
                           bValid = false ;
                        else {
                           if ( dParU > dPar1) {
                              dPar1 = dParU ;
                              vtN1 = vtUD ;
                           }
                           if ( dParD < dPar2) {
                              dPar2 = dParD ;
                              vtN2 = -vtUD ;
                           }
                        }
                     }
                     else {
                        if ( dParU < dPar1 + EPS_ZERO  ||  dParD > dPar2 - EPS_ZERO)
                           bValid = false ;
                        else {
                           if ( dParD > dPar1) {
                              dPar1 = dParD ;
                              vtN1 = -vtUD ;
                           }
                           if ( dParU < dPar2) {
                              dPar2 = dParU ;
                              vtN2 = vtUD ;
                           }
                        }
                     }
                  }
               }
              // altrimenti praticamente parallelo
               else {
                 // se esterno ad almeno uno invalida tutto
                  if ( dDistU > 0  ||  dDistD > 0)
                     bValid = false ;
                 // altrimenti non cambia niente
               }
            }

           // Taglio con la faccia sinistra
            if ( bValid) {
              // Distanza con segno del punto di riferimento del dexel dal piano della faccia sinistra
               double dDistS = ( ptPoly - ptFacet135) * vtNs ;
              // Se dexel non parallelo alla faccia
               if ( abs( vtPoly * vtNs) > COS_ORTO_ANG_ZERO) {
                 // posizione parametrica della intersezione
                  double dParS = -dDistS / ( vtPoly * vtNs) ;
                 // verifico limitazioni su inizio e fine dell'intervallo
                  int nLimit = 0 ;
                 // se limita inizio
                  if ( vtPoly * vtNs < 0) {
                    // se oltre la fine, invalida tutto
                     if ( dParS > dPar2 - EPS_ZERO)
                        bValid = false ;
                    // se altrimenti solo oltre inizio, riduce
                     else if ( dParS >= dPar1) {
                        dPar1 = dParS ;
                        nLimit = 1 ;
                     }
                  }
                 // altrimenti limita fine
                  else {
                    // se prima dell'inizio, invalida tutto
                     if ( dParS < dPar1 + EPS_ZERO)
                        bValid = false ;
                    // se altrimenti solo prima della fine, riduce
                     else if ( dParS <= dPar2) {
                        dPar2 = dParS ;
                        nLimit = 2 ; 
                     }        
                  }
                 // se limita, devo aggiornare la normale
                  if ( nLimit != 0) {
                     Vector3d vtNewN = -vtNs ;
                     if ( bRecalNorm) {
                        Vector3d vtRadial( 0, dMinRad * dCos, dMinRad * dSin) ;
                        vtRadial.Normalize() ;
                        Vector3d vtOrigMaxR = - vtArcNormMaxR.x * vtRadial - vtArcNormMaxR.z * X_AX ;
                        Vector3d vtOrigMinR = - vtArcNormMinR.x * vtRadial - vtArcNormMinR.z * X_AX ;
                        vtOrigMaxR.Normalize() ;
                        vtOrigMinR.Normalize() ;
                        Point3d ptInt = ptPoly + dParS * vtPoly ;
                        vtNewN = ( ( dDeltaZ - ptInt.z + dLenZ) / dDeltaZ) * vtOrigMinR + ( ( ptInt.z - dLenZ) / dDeltaZ) * vtOrigMaxR ;
                        vtNewN.Normalize() ;
                     }
                     if ( nLimit == 1)
                        vtN1 = vtNewN ;
                     else
                        vtN2 = vtNewN ;
                  }
               }
              // altrimenti praticamente parallelo
               else {
                 // se esterno invalida tutto
                  if ( dDistS > 0)
                     bValid = false ;
                 // altrimenti non cambia niente
               }
            }

           // Taglio con la faccia destra
            if ( bValid) {
              // Distanza con segno del punto di riferimento del dexel dal piano della faccia destra
               double dDistD = ( ptPoly - ptFacet246) * vtNd ;
              // Se dexel non parallelo alla faccia
               if ( abs( vtPoly * vtNd) > COS_ORTO_ANG_ZERO) {
                 // posizione parametrica della intersezione
                  double dParD = -dDistD / ( vtPoly * vtNd) ;
                 // verifico limitazioni su inizio e fine dell'intervallo
                  int nLimit = 0 ;
                 // se limita inizio
                  if ( vtPoly * vtNd < 0) {
                    // se oltre la fine, invalida tutto
                     if ( dParD > dPar2 - EPS_ZERO)
                        bValid = false ;
                    // se altrimenti solo oltre inizio, riduce
                     else if ( dParD >= dPar1) {
                        dPar1 = dParD ;
                        nLimit = 1 ;
                     }
                  }
                 // altrimenti limita fine
                  else {
                    // se prima dell'inizio, invalida tutto
                     if ( dParD < dPar1 + EPS_ZERO)
                        bValid = false ;
                    // se altrimenti solo prima della fine, riduce
                     else if ( dParD <= dPar2) {
                        dPar2 = dParD ;
                        nLimit = 2 ; 
                     }        
                  }
                 // se limita, devo aggiornare la normale
                  if ( nLimit != 0) {
                     Vector3d vtNewN = -vtNd ;
                     if ( bRecalNorm) {
                        Vector3d vtRadial( 0, dMinRad * dCos, -dMinRad * dSin) ;
                        vtRadial.Normalize() ;
                        Vector3d vtOrigMaxR = - vtArcNormMaxR.x * vtRadial - vtArcNormMaxR.z * X_AX ;
                        Vector3d vtOrigMinR = - vtArcNormMinR.x * vtRadial - vtArcNormMinR.z * X_AX ;
                        vtOrigMaxR.Normalize() ;
                        vtOrigMinR.Normalize() ;
                        Point3d ptInt = ptPoly + dParD * vtPoly ;
                        vtNewN = ( ( dDeltaZ - abs( ptInt.z) + dLenZ) / dDeltaZ) * vtOrigMinR + ( ( abs( ptInt.z) - dLenZ) / dDeltaZ) * vtOrigMaxR ;
                        vtNewN.Normalize() ;
                     }
                     if ( nLimit == 1)
                        vtN1 = vtNewN ;
                     else
                        vtN2 = vtNewN ;
                  }
               }
              // altrimenti praticamente parallelo
               else {
                 // se esterno invalida tutto
                  if ( dDistD > 0)
                     bValid = false ;
                 // altrimenti non cambia niente
               }
            }

           // Se rimasto qualcosa lo sottraggo al dexel
            if ( bValid) {
              // Punti di intersezione e normali nel sistema griglia
               ptInt1 = ptC + dPar1 * Z_AX ;
               vtN1.ToGlob( PolyFrame) ;
               ptInt2 = ptC + dPar2 * Z_AX ; 
               vtN2.ToGlob( PolyFrame) ;
              // Eseguo sottrazione
               SubtractIntervals( nGrid, i, j, ptInt1.z, ptInt2.z, vtN1, vtN2, nToolNum) ;
            }

           // Se movimento non ortogonale all'asse sottraggo i cilindri ellittici
            if ( dLongLen > EPS_SMALL) {
            
              // Traslazione ellisse di punta
               ConusFrame.ChangeOrig( ptV + vtV1 * dl) ;
               if ( IntersLineEllipticalCylinder( ptC, Z_AX, ConusFrame, dMinRad, dLongLen, dOrtLen,
                                                  bTapB, true, ptInt1, vtN1, ptInt2, vtN2)) {
                  SubtractIntervals( nGrid, i, j, ptInt1.z, ptInt2.z, vtN1, vtN2, nToolNum) ;
               }
            
              // Traslazione ellisse di base
               ConusFrame.ChangeOrig( ptV + vtV1 * dL) ;
               if ( IntersLineEllipticalCylinder( ptC, Z_AX, ConusFrame, dMaxRad, dLongLen, dOrtLen,
                                                  true, bTapT, ptInt1, vtN1, ptInt2, vtN2)) {
                  SubtractIntervals( nGrid, i, j, ptInt1.z, ptInt2.z, vtN1, vtN2, nToolNum) ;
               }
            }
         }
      }
      return true ;
   }

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

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

            Point3d ptInt1, ptInt2 ;
            Vector3d vtN1, vtN2 ;

           // Cono 
            ConusFrame.ChangeOrig( ptV) ;
            if ( IntersLineConus( ptC, Z_AX, ConusFrame, dTan, dl, dL, bTapB, true, ptInt1, vtN1, ptInt2, vtN2)) {
               if ( bRecalNorm) {
                  vtN1 = AdjustConeNormal( ptInt1, vtN1, ptV, vtToolDir, dMinRad, dDeltaR, vtArcNormMinR, vtArcNormMaxR) ;
                  vtN2 = AdjustConeNormal( ptInt2, vtN2, ptV, vtToolDir, dMinRad, dDeltaR, vtArcNormMinR, vtArcNormMaxR) ;
               }
               SubtractIntervals( nGrid, i, j, ptInt1.z, ptInt2.z, vtN1, vtN2, nToolNum) ;
            }

           // Traslazione ellisse
            ConusFrame.ChangeOrig( ptV + vtV1 * dL) ;
            if ( IntersLineEllipticalCylinder( ptC, Z_AX, ConusFrame, dMaxRad, dLongLen, dOrtLen,
                                               true, bTapT, ptInt1, vtN1, ptInt2, vtN2)) {
               SubtractIntervals( nGrid, i, j, ptInt1.z, ptInt2.z, vtN1, vtN2, nToolNum) ;
            }
         }
      }
      return true ;
   }
}

//----------------------------------------------------------------------------
bool 
VolZmap::CompPar_Milling( int nGrid, double dLenX, double dLenY, double dLenZ,
                          const Point3d& ptS, const Point3d& ptE,
                          const Vector3d& vtToolDir, const Vector3d& vtAux, int nToolNum)
{
   Point3d ptMyS = ptS - EPS_SMALL * vtToolDir ;
   Point3d ptMyE = ptE - EPS_SMALL * vtToolDir ;

   dLenX -= ( 2 * EPS_SMALL) ;
   dLenY -= ( 2 * EPS_SMALL) ;
   dLenZ -= ( 2 * EPS_SMALL) ;

  // Controllo sull'interferenza utensile-solido
   int nStartI, nStartJ, nEndI, nEndJ ;
   if ( ! TestParaBBox( nGrid, ptMyS, ptMyE, vtToolDir, vtAux, dLenX, dLenY, dLenZ, nStartI, nStartJ, nEndI, nEndJ))
      return true ;

  // Vettore di movimento
   Vector3d vtMove = ptMyE - ptMyS ;

  // Sistema di riferimento del solido nella posizione iniziale (con vtMove derivo quello finale)
   Frame3d ParaFrame ;
   Vector3d vtV3 = vtToolDir ;
   Vector3d vtV1 = vtAux ^ vtV3 ;
   Vector3d vtV2 = vtV3 ^ vtV1 ;
   ParaFrame.Set( ORIG, vtV1, vtV2, vtV3) ;
   Point3d ptV = ptMyS - dLenZ * ParaFrame.VersZ() - ( 0.5 * dLenX) * ParaFrame.VersX() ;
   ParaFrame.ChangeOrig( ptV) ;

  // Studio del volume asportato durante il moto
   Frame3d MotionFrame ;
   Vector3d vtDiagA = dLenX * ParaFrame.VersX() + dLenY * ParaFrame.VersY() ;
   Vector3d vtDiagB = dLenX * ParaFrame.VersX() - dLenY * ParaFrame.VersY() ;
   Vector3d vtW3 = ( abs( vtDiagA * vtMove) <= abs( vtDiagB * vtMove) ? vtDiagA : vtDiagB) ;
   Point3d ptOrig = ptMyS - vtW3 / 2 - vtToolDir * ( dLenZ / 2) ;
   double dLenW3 = vtW3.Len() ;
   vtW3 /= dLenW3 ;
   double dDeltaW3 = vtMove * vtW3 ;
   Vector3d vtW1 = vtMove - dDeltaW3 * vtW3 ;
   double dLenW1 = vtW1.Len() ;
   vtW1 /= dLenW1 ;
   Vector3d vtW2 = vtW3 ^ vtW1 ;
   MotionFrame.Set( ptOrig, vtW1, vtW2, vtW3) ;

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

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

         Point3d ptInt1, ptInt2 ;
         Vector3d vtN1, vtN2 ;

        // Solido nella posizione iniziale
         ParaFrame.ChangeOrig( ptV) ;
         if ( IntersLineMyPolyhedron( ptC, Z_AX, ParaFrame, dLenX, dLenY, dLenZ, 0, ptInt1, vtN1, ptInt2, vtN2)) {
            SubtractIntervals( nGrid, i, j, ptInt1.z, ptInt2.z, vtN1, vtN2, nToolNum) ; 
         }

        // Solido nella posizione finale
         ParaFrame.ChangeOrig( ptV + vtMove) ;
         if ( IntersLineMyPolyhedron( ptC, Z_AX, ParaFrame, dLenX, dLenY, dLenZ, 0, ptInt1, vtN1, ptInt2, vtN2)) {
            SubtractIntervals( nGrid, i, j, ptInt1.z, ptInt2.z, vtN1, vtN2, nToolNum) ;
         }

        // Volume spazzato nel moto
         if ( IntersLineMyPolyhedron( ptC, Z_AX, MotionFrame, dLenW1, dLenZ, dLenW3, dDeltaW3, ptInt1, vtN1, ptInt2, vtN2)) {
            SubtractIntervals( nGrid, i, j, ptInt1.z, ptInt2.z, vtN1, vtN2, nToolNum) ; 
         }
      }
   }

   return true ;
}

// ---------- SFERA ----------------------------------------------------------

//---------------------------------------------------------------------------- pizza
bool
VolZmap::CompBall_Milling( int nGrid, const Point3d& ptLs, const Point3d& ptLe, double dRad, int nToolNum)
{
  // Verifico interferisca
   int nStartI, nStartJ, nEndI, nEndJ ;
   if ( ! TestCompoBBox( nGrid, ptLs, ptLe, V_NULL, dRad, 0, 0, nStartI, nStartJ, nEndI, nEndJ))
      return true ;
  // Vettore modivemnto
   Vector3d vtV = ptLe - ptLs ;
   double dLengthPath = vtV.Len() ;
   if ( dLengthPath > EPS_ZERO)
      vtV /= dLengthPath ;
  // Riferimento per cilindro inviluppo della sfera lungo il movimento
   Frame3d CylFrame ;
   CylFrame.Set( ptLs, vtV) ;

   double dSqRad = dRad * dRad ;

   for ( int i = nStartI ; i <= nEndI ; ++ i) {
      for ( 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) ;

        // Sfera in posizione start
         double dStSqDXY = SqDistXY( ptC, ptLs) ;
         if ( dStSqDXY < dSqRad) {
            double dMin = ptLs.z - sqrt( dSqRad - dStSqDXY) ;
            Vector3d vtNmin = ptLs - Point3d( dX, dY, dMin) ;
            vtNmin.Normalize() ;
            double dMax = ptLs.z + sqrt( dSqRad - dStSqDXY) ;
            Vector3d vtNmax = ptLs - Point3d( dX, dY, dMax) ;
            vtNmax.Normalize() ;
            SubtractIntervals( nGrid, i, j, dMin, dMax, vtNmin, vtNmax, nToolNum) ;
         }

        // Sfera in posizione end
         double dEnSqDXY = SqDistXY( ptC, ptLe) ;
         if ( dEnSqDXY < dSqRad) {
            double dMin = ptLe.z - sqrt( dSqRad - dEnSqDXY) ;
            Vector3d vtNmin = ptLe - Point3d( dX, dY, dMin) ;
            vtNmin.Normalize() ;
            double dMax = ptLe.z + sqrt( dSqRad - dEnSqDXY) ;
            Vector3d vtNmax = ptLe - Point3d( dX, dY, dMax) ;
            vtNmax.Normalize() ;
            SubtractIntervals( nGrid, i, j, dMin, dMax, vtNmin, vtNmax, nToolNum) ;
         }

        // Cilindro inviluppo della sfera
         Point3d ptInt1, ptInt2 ;
         Vector3d vtN1, vtN2 ;
         if ( IntersLineCylinder( ptC, Z_AX, CylFrame, dLengthPath, dRad, false, false, ptInt1, vtN1, ptInt2, vtN2)) {
            SubtractIntervals( nGrid, i, j, ptInt1.z, ptInt2.z, vtN1, vtN2, nToolNum) ;
         }
      }
   }

   return true ;
}

// ------------------------- Utensili additivi ---------------------------------------------------------------------------------

//----------------------------------------------------------------------------
bool
VolZmap::AddingMotion( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtAx)
{
  // Dimensioni lineari dell'utensile
   double dHei = m_vTool[m_nCurrTool].GetHeigth() ;
   double dRad = m_vTool[m_nCurrTool].GetRadius() ;
   double dCornerRad = m_vTool[m_nCurrTool].GetCornRadius() ;

  // Utensile sfiancato
   if ( dCornerRad * dCornerRad - 0.25 * dHei * dHei > 0) {
      return AddingGeneral( nGrid, ptS, ptE, vtAx) ;   
   }
  // Utensile sferico
   else if ( dRad - dCornerRad < EPS_SMALL) {
      return AddingSphere( nGrid, ptS - dRad * vtAx, ptE - dRad * vtAx, dRad) ;
   }
  // Utensile cilindro
   else if ( dCornerRad < EPS_SMALL) {
      return AddingCylinder( nGrid, ptS, ptE, vtAx, dHei, dRad) ;
   }
  // Utensile naso di toro
   else {
      return AddingGeneral( nGrid, ptS, ptE, vtAx) ;
   }
}

//----------------------------------------------------------------------------
bool
VolZmap::AddingGeneral( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtAx)
{
  // Descrizione geometrica del moto  
   Point3d ptI = ptS ;
   Point3d ptF = ptE ;
   Vector3d vtMove = ptE - ptS ;
  // Vettore delle normali agli archi
   const VCT3DVECTOR& vArcNorm = m_vTool[m_nCurrTool].GetArcNormalVec() ;
  // Poinché l'asse utensile è parallelo all'asse Z, definisco un sistema di
  // riferimento ad hoc in cui le normali agli archi giacciano nel piano XZ.
   Frame3d frNormFrame ;
   frNormFrame.Set( ORIG, X_AX, -Z_AX, Y_AX) ;
  // Ciclo sulle curve del profilo
   const CurveComposite& ToolProfile = m_vTool[m_nCurrTool].GetApproxOutline() ;
   int i = - 1 ;
   const ICurve* pCurve = ToolProfile.GetCurve( ++ i) ;
   while ( pCurve != nullptr) {
   
      double dHeight = 0 ;

     // Se segmento
      if ( pCurve->GetType() == CRV_LINE) {
        // Recupero gli estremi
         const ICurveLine* pLine = GetCurveLine( pCurve) ;
         Point3d ptStart = pLine->GetStart() ;
         Point3d ptEnd = pLine->GetEnd() ;
         int nNormNum = pLine->GetTempProp();
         Vector3d vtNormSt, vtNormEn;
         if ( nNormNum != 0) {
            vtNormSt = vArcNorm[nNormNum - 1] ;
            vtNormEn = vArcNorm[nNormNum] ;
            vtNormSt.ToLoc(frNormFrame);
            vtNormEn.ToLoc(frNormFrame);
         }
        // Ne determino l'altezza
         dHeight = abs( ptStart.y - ptEnd.y) ;
         if ( dHeight > EPS_SMALL) {
           // Se X costante, è un cilindro
            if ( abs( ptStart.x - ptEnd.x) < EPS_SMALL) {
               double dRadius = ptStart.x ;
               if (dRadius > 10 * EPS_SMALL)
                  AddingCylinder( nGrid, ptI, ptF, vtAx, dHeight, dRadius) ;
            }
           // Se X crescente, è un cono con vettore equiverso a quello dell'utensile
            else if ( ptStart.x > ptEnd.x) {
               double dMaxRad = ptStart.x ;
               double dMinRad = ptEnd.x ;
               AddingTruncatedCone( nGrid, ptI, ptF, vtAx, dMaxRad, dMinRad, dHeight, vtNormSt, vtNormEn) ;
            }
           // Se X decrescente, è un cono con vettore opposto a quello dell'utensile
            else if ( ptStart.x < ptEnd.x) {
               double dMaxRad = ptEnd.x ; 
               double dMinRad = ptStart.x ;
               Point3d ptIn = ptI - vtAx * dHeight ;
               Point3d ptFn = ptIn + vtMove ;
               vtNormEn.z *= -1 ;
               vtNormSt.z *= -1 ;
               AddingTruncatedCone( nGrid, ptIn, ptFn, - vtAx, dMaxRad, dMinRad, dHeight, vtNormEn, vtNormSt) ;
               
            }  
           // Passo alla curva successiva
            pCurve = ToolProfile.GetCurve( ++ i) ;
         }
         else {
           // Passo alla curva successiva
            pCurve = ToolProfile.GetCurve( ++ i) ;
         }
      }

     // Determino le posizioni iniziale e finale del componente successivo
      ptI = ptI - vtAx * dHeight ;
      ptF = ptI + vtMove ;
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::AddingCylinder( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtAx, double dHei, double dRad)
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

   // Verifica sull'interferenza utensile Zmap
   int nStartI, nStartJ, nEndI, nEndJ ;
   if ( ! TestCompoBBox( nGrid, ptS, ptE, vtAx, dRad, dRad, dHei, nStartI, nStartJ, nEndI, nEndJ))
      return true ;

   Vector3d vtV1 = ptE - ptS ;
   double dLen1 = vtV1.Len() ;
   vtV1 /= dLen1 ;

   double dMyTol = 0 ;
   Frame3d CylFrame, PolyFrame ;
   if ( ! CylFrame.Set( ptS - dHei * vtAx, vtAx))
      return false ;
   if ( ! PolyFrame.Set( ptS - ( dHei + dMyTol) * vtAx, vtAx, vtV1))
      return false ;
   for ( int i = nStartI ; i <= nEndI ; ++ i) {
      for ( int j = nStartJ ; j <= nEndJ ; ++ j) {
         
         Point3d ptC( ( i + 0.5) * m_dStep, ( j + 0.5) * m_dStep, 0) ;
                                           
         Point3d ptInt1, ptInt2 ;
         Vector3d vtN1, vtN2 ;

         if ( IntersLineCylinder( ptC, Z_AX, CylFrame, dHei, dRad, true, true,
                                  ptInt1, vtN1, ptInt2, vtN2)) {
            AddIntervals( nGrid, i, j, ptInt1.z, ptInt2.z, - vtN1, - vtN2, CurrTool.GetToolNum()) ;
         }

         if ( IntersLineCylinder( ptC - dLen1 * vtV1, Z_AX, CylFrame, dHei, dRad, true, true,
                                  ptInt1, vtN1, ptInt2, vtN2)) {
            AddIntervals( nGrid, i, j, ptInt1.z + dLen1 * vtV1.z, ptInt2.z + dLen1 * vtV1.z, - vtN1, - vtN2, CurrTool.GetToolNum()) ;
         }

         if ( IntersLineMyPolyhedron( ptC, Z_AX, PolyFrame, dLen1, 2 * ( dRad + dMyTol), dHei + 2 * dMyTol, 0,
                                      ptInt1, vtN1, ptInt2, vtN2)) {
            AddIntervals( nGrid, i, j, ptInt1.z, ptInt2.z, - vtN1, - vtN2, CurrTool.GetToolNum()) ;
         }
      }
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::AddingTruncatedCone( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtAx,
                              double dMaxRad, double dMinRad, double dHei, 
                              const Vector3d& vtArcNormMaxR, const Vector3d& vtArcNormMinR)
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

   // Verifico interferenza
   int nStartI, nStartJ, nEndI, nEndJ ;
   if ( ! TestCompoBBox( nGrid, ptS, ptE, vtAx, dMaxRad, dMinRad, dHei, nStartI, nStartJ, nEndI, nEndJ))
      return true ;

  // Geometria del cono
   double dDeltaR = dMaxRad - dMinRad ;
  
  // Studio simmetrie
   Point3d ptI = ( vtAx * ( ptE - ptS) > 0 ? ptS : ptE) ;
   Point3d ptF = ( vtAx * ( ptE - ptS) > 0 ? ptE : ptS) ;

   double dL = ( dMaxRad * dHei) / dDeltaR ;
   double dl = dL - dHei ;

   Point3d ptV = ptI - vtAx * dL ;

  // Vettori caratteristici del movimento
  // Elimino eventuali componenti del moto lungo l'asse.
   Vector3d vtMove = ptF - ptI ;
   Vector3d vtMvLong = ( vtMove * vtAx) * vtAx ;
   Vector3d vtMvOrt = vtMove - vtMvLong ;

  // Terna destrorsa e unitaria
   Vector3d vtV1 = vtAx ;
   Vector3d vtV2 = vtMvOrt ;  vtV2.Normalize() ;
   Vector3d vtV3 = vtV1 ^ vtV2 ;

  // Sistema di riferimento intrinseco del movimento
   Frame3d ConusFrame ;
   if ( ! ConusFrame.Set( ptV, vtV1, vtV2, vtV3))
      return false ;

  // Dimensioni lineari movimento
   double dLongLen = 0 ;
   double dOrtLen = vtMvOrt.Len() ;

  // Apertura del cono 
   double dTan = dDeltaR / dHei ;
   double dRatio = dLongLen / dOrtLen ;

  // Per costruire piani laterali poliedro interno
   double dCos = dTan * dRatio ;
   double dSin = ( 1 - dCos * dCos > 0 ? sqrt( 1 - dCos * dCos) : 0) ;

  // Dimensioni lineari descriventi il poliedro interno
   double dLenX = dLongLen ;
   double dLenY = dOrtLen ;
   double dLenZ = dSin * dMinRad ;
   double dDeltaX = dHei ;
   double dDeltaY = dCos * dDeltaR ;
   double dDeltaZ = dSin * dDeltaR ;

  // Sistema di riferimento poliedro
   Point3d ptO = ptV + vtV1 * dl + vtV2 * ( dCos * dMinRad) ;
   Frame3d PolyFrame ;
   if ( ! PolyFrame.Set( ptO, vtV1, vtV2, vtV3))
      return false ;

  // Versori piani nel riferimento poliedro ( riferiti al sistema di riferimento) :
  // Sx, Dx
   Vector3d vtNs( - dTan, dCos, dSin) ;
   vtNs.Normalize() ;
   Vector3d vtNd( - dTan, dCos, - dSin) ;
   vtNd.Normalize() ;
  // Iniziale e finale
   Vector3d vtIF( - dDeltaY, dDeltaX, 0) ;
   vtIF.Normalize() ;
  // Up e Down
   Vector3d vtUD( - dLenY, dLenX, 0) ;
   vtUD.Normalize() ;

  // Punti dei piani (sempre espressi nel sistema PolyFrame)
   Point3d ptFacet135( 0, 0, dLenZ) ;
   Point3d ptFacet246( dLenX + dDeltaX, dLenY + dDeltaY, - dLenZ - dDeltaZ) ;

   Vector3d vtUmv = vtMove ;  vtUmv.Normalize() ;

   for ( int i = nStartI ; i <= nEndI ; ++ i) {
      for ( int j = nStartJ ; j <= nEndJ ; ++ j) {
         
         Point3d ptC( ( i + 0.5) * m_dStep, ( j + 0.5) * m_dStep, 0) ;

         Point3d ptInt1, ptInt2 ;
         Vector3d vtN1, vtN2 ;

        // Cono iniziale
         ConusFrame.ChangeOrig( ptV) ;
         if ( IntersLineConus( ptC, Z_AX, ConusFrame, dTan, dl, dL, true, true, ptInt1, vtN1, ptInt2, vtN2)) {
            vtN1 *= - 1 ;
            vtN2 *= - 1 ;
            if ( ! ( vtArcNormMaxR.IsSmall()  ||  vtArcNormMinR.IsSmall())) {
               if ( ! AreSameOrOppositeVectorEpsilon( vtN1, vtAx, 0.1 * EPS_SMALL)) {
                  Vector3d vtL1 = ptInt1 - ptV ;
                  vtL1 -= ( vtL1 * vtAx) * vtAx ;
                  double dL1 = vtL1.Len() ;
                  vtL1 /= dL1 ;
                  Vector3d vtOriginalN1 = ( ( dDeltaR - dL1 + dMinRad) / dDeltaR) * vtArcNormMinR + ((dL1 - dMinRad) / dDeltaR) * vtArcNormMaxR;
                  vtOriginalN1.Normalize() ;
                  vtN1 = vtOriginalN1.z * vtAx + vtOriginalN1.x * vtL1 ;
                  vtN1.Normalize() ;
               }
               if ( ! AreSameOrOppositeVectorEpsilon( vtN2, vtAx, 0.1 * EPS_SMALL)) {
                  Vector3d vtL2 = ptInt2 - ptV ;
                  vtL2 -= ( vtL2 * vtAx) * vtAx ;
                  double dL2 = vtL2.Len() ;
                  vtL2 /= dL2 ;
                  Vector3d vtOriginalN2 = ( ( dDeltaR - dL2 + dMinRad) / dDeltaR) * vtArcNormMinR + ( ( dL2 - dMinRad) / dDeltaR) * vtArcNormMaxR ;
                  vtOriginalN2.Normalize() ;
                  vtN2 = vtOriginalN2.z * vtAx + vtOriginalN2.x * vtL2 ;
                  vtN2.Normalize() ;
               }
            }
            AddIntervals( nGrid, i, j, ptInt1.z, ptInt2.z, vtN1, vtN2, CurrTool.GetToolNum()) ;
         }
      
        // Cono finale
         ConusFrame.ChangeOrig( ptV + vtMove) ;
         if ( IntersLineConus( ptC, Z_AX, ConusFrame, dTan, dl, dL, true, true, ptInt1, vtN1, ptInt2, vtN2)) {
            vtN1 *= - 1 ;
            vtN2 *= - 1 ;
            if ( ! ( vtArcNormMaxR.IsSmall()  ||  vtArcNormMinR.IsSmall())) {
               if ( ! AreSameOrOppositeVectorEpsilon( vtN1, vtAx, 0.1 * EPS_SMALL)) {
                  Vector3d vtL1 = ptInt1 - ptV - vtMove ;
                  vtL1 -= ( vtL1 * vtAx) * vtAx ;
                  double dL1 = vtL1.Len() ;
                  vtL1 /= dL1 ;
                  Vector3d vtOriginalN1 = ( ( dDeltaR - dL1 + dMinRad) / dDeltaR) * vtArcNormMinR + ( ( dL1 - dMinRad) / dDeltaR) * vtArcNormMaxR ;
                  vtOriginalN1.Normalize() ;
                  vtN1 = vtOriginalN1.z * vtAx + vtOriginalN1.x * vtL1 ;
                  vtN1.Normalize() ;
               }
               if ( ! AreSameOrOppositeVectorEpsilon(vtN2, vtAx, 0.1 * EPS_SMALL)) {
                  Vector3d vtL2 = ptInt2 - ptV - vtMove ;
                  vtL2 -= (vtL2 * vtAx) * vtAx;
                  double dL2 = vtL2.Len() ;
                  vtL2 /= dL2 ;
                  Vector3d vtOriginalN2 = ( ( dDeltaR - dL2 + dMinRad) / dDeltaR) * vtArcNormMinR + ( ( dL2 - dMinRad) / dDeltaR) * vtArcNormMaxR ;
                  vtOriginalN2.Normalize() ;
                  vtN2 = vtOriginalN2.z * vtAx + vtOriginalN2.x * vtL2 ;
                  vtN2.Normalize() ;
               }
            }
            AddIntervals( nGrid, i, j, ptInt1.z, ptInt2.z, vtN1, vtN2, CurrTool.GetToolNum()) ;
         }
         
        // Solido interno
         Point3d ptPoly = ptC ;
         Vector3d vtPoly = Z_AX ;

         ptPoly.ToLoc( PolyFrame) ;
         vtPoly.ToLoc( PolyFrame) ;

         Point3d ptPoly1 = ptPoly + ( ( ( ptFacet135 - ptPoly) * vtNs) / ( vtPoly * vtNs)) * vtPoly ;
         Point3d ptPoly2 = ptPoly + ( ( ( ptFacet246 - ptPoly) * vtNd) / ( vtPoly * vtNd)) * vtPoly ;
         Point3d ptPoly3 = ptPoly + ( ( ( ptFacet135 - ptPoly) * vtIF) / ( vtPoly * vtIF)) * vtPoly ;
         Point3d ptPoly4 = ptPoly + ( ( ( ptFacet246 - ptPoly) * vtIF) / ( vtPoly * vtIF)) * vtPoly ;
         Point3d ptPoly5 = ptPoly + ( ( ( ptFacet135 - ptPoly) * vtUD) / ( vtPoly * vtUD)) * vtPoly ;
         Point3d ptPoly6 = ptPoly + ( ( ( ptFacet246 - ptPoly) * vtUD) / ( vtPoly * vtUD)) * vtPoly ;

         int nIntNum = 0 ;

        // Intersezione con la prima faccia
         if ( abs( vtPoly * vtNs) > COS_ORTO_ANG_ZERO) {
            if ( dLenY * ( ptPoly1.x + EPS_SMALL) > dLenX * ptPoly1.y  && 
                 dLenY * ( ptPoly1.x - dDeltaX - EPS_SMALL) < dLenX * ( ptPoly1.y - dDeltaY)  && 
                 dDeltaX * ( ptPoly1.y + EPS_SMALL) > dDeltaY * ptPoly1.x  && 
                 dDeltaX * ( ptPoly1.y - dLenY - EPS_SMALL) < dDeltaY * ( ptPoly1.x - dLenX)) {
               ptInt1 = ptPoly1 ;
               vtN1 = - vtNs ;
               if ( ! ( vtArcNormMaxR.IsSmall()  ||  vtArcNormMinR.IsSmall())) {
                  Vector3d vtRadial( 0, dMinRad * dCos, dMinRad * dSin) ;
                  vtRadial.Normalize() ;
                  Vector3d vtOrigMaxR = - vtArcNormMaxR.x * vtRadial - vtArcNormMaxR.z * X_AX ;
                  Vector3d vtOrigMinR = - vtArcNormMinR.x * vtRadial - vtArcNormMinR.z * X_AX ;
                  vtOrigMaxR.Normalize() ;
                  vtOrigMinR.Normalize() ;
                  vtN1 = - ( ( dDeltaZ - ptInt1.z + dLenZ) / dDeltaZ) * vtOrigMinR - ( ( ptInt1.z - dLenZ) / dDeltaZ) * vtOrigMaxR ;
                  vtN1.Normalize() ;
               }
               ++ nIntNum ;
            }
         }
        // Intersezione con la seconda faccia
         if ( abs( vtPoly * vtNd) > COS_ORTO_ANG_ZERO) {
            if ( dLenY * ( ptPoly2.x + EPS_SMALL) > dLenX * ptPoly2.y  && 
                 dLenY * ( ptPoly2.x - dDeltaX - EPS_SMALL) < dLenX * ( ptPoly2.y - dDeltaY)  && 
                 dDeltaX * ( ptPoly2.y + EPS_SMALL) > dDeltaY * ptPoly2.x  && 
                 dDeltaX * ( ptPoly2.y - dLenY - EPS_SMALL) < dDeltaY * ( ptPoly2.x - dLenX)) {

               if ( nIntNum == 0) {
                  ptInt1 = ptPoly2 ;
                  vtN1 = - vtNd ;
                  if ( ! ( vtArcNormMaxR.IsSmall()  ||  vtArcNormMinR.IsSmall())) {
                     Vector3d vtRadial( 0, dMinRad * dCos, - dMinRad * dSin) ;
                     vtRadial.Normalize() ;
                     Vector3d vtOrigMaxR = - vtArcNormMaxR.x * vtRadial - vtArcNormMaxR.z * X_AX ;
                     Vector3d vtOrigMinR = - vtArcNormMinR.x * vtRadial - vtArcNormMinR.z * X_AX ;
                     vtOrigMaxR.Normalize() ;
                     vtOrigMinR.Normalize() ;
                     vtN1 = - ( ( dDeltaZ - abs( ptInt1.z) + dLenZ) / dDeltaZ) * vtOrigMinR - ( ( abs( ptInt1.z) - dLenZ) / dDeltaZ) * vtOrigMaxR ;
                     vtN1.Normalize() ;
                  }
                  ++ nIntNum ;
               }         
               else if ( ( ptInt1 - ptPoly2).SqLen() > SQ_EPS_SMALL) {
                  ptInt2 = ptPoly2 ;
                  vtN2 = - vtNd ;
                  if ( ! ( vtArcNormMaxR.IsSmall()  ||  vtArcNormMinR.IsSmall())) {
                     Vector3d vtRadial( 0, dMinRad * dCos, -dMinRad * dSin) ;
                     vtRadial.Normalize() ;
                     Vector3d vtOrigMaxR = -vtArcNormMaxR.x * vtRadial - vtArcNormMaxR.z * X_AX ;
                     Vector3d vtOrigMinR = -vtArcNormMinR.x * vtRadial - vtArcNormMinR.z * X_AX ;
                     vtOrigMaxR.Normalize() ;
                     vtOrigMinR.Normalize() ;
                     vtN2 = - ( ( dDeltaZ - abs( ptInt2.z) + dLenZ) / dDeltaZ) * vtOrigMinR - ( ( abs( ptInt2.z) - dLenZ) / dDeltaZ) * vtOrigMaxR ;
                     vtN2.Normalize() ;
                  }
                  ++ nIntNum ;
               }
            }
         }
        // Intersezione con la terza faccia
         if ( abs( vtPoly * vtIF) > COS_ORTO_ANG_ZERO) {
            if ( nIntNum < 2  &&
                 ptPoly3.x > - EPS_SMALL  &&  ptPoly3.x < dDeltaX + EPS_SMALL  && 
                 dDeltaX * abs( ptPoly3.z) < dDeltaX * dLenZ + dDeltaZ * ptPoly3.x + dDeltaX * EPS_SMALL) {

               if ( nIntNum == 0) {
                  ptInt1 = ptPoly3 ;
                  vtN1 = - vtIF ; 
                  ++ nIntNum ;
               }
               else if ( ( ptInt1 - ptPoly3).SqLen() > SQ_EPS_SMALL) {
                  ptInt2 = ptPoly3 ;
                  vtN2 = - vtIF ;
                  ++ nIntNum ;
               }
            }
         }

        // Intersezione con la quarta faccia
         if ( abs( vtPoly * vtIF) > COS_ORTO_ANG_ZERO) {
            if ( nIntNum < 2  &&
                 ptPoly4.x > dLenX - EPS_SMALL  &&  ptPoly4.x < dLenX + dDeltaX + EPS_SMALL && 
                 dDeltaX * abs( ptPoly4.z) < dDeltaX * dLenZ + dDeltaZ * ( ptPoly4.x - dLenX) + dDeltaX * EPS_SMALL) {

               if ( nIntNum == 0) {
                  ptInt1 = ptPoly4 ;
                  vtN1 = vtIF ; 
                  ++ nIntNum ;
               }
               else if ( ( ptInt1 - ptPoly4).SqLen() > SQ_EPS_SMALL) {
                  ptInt2 = ptPoly4 ;
                  vtN2 = vtIF ;
                  ++ nIntNum ;
               }
            }
         }

        // Intersezione con la quinta faccia
         if ( abs( vtPoly * vtUD) > COS_ORTO_ANG_ZERO) {
            if ( nIntNum < 2  &&
                 ptPoly5.y >= 0  &&  ptPoly5.y <= dLenY  &&
                 abs( ptPoly5.z) <= dLenZ) {

               if ( nIntNum == 0) {
                  ptInt1 = ptPoly5 ;
                  vtN1 = vtUD ;
                  ++ nIntNum ;
               }
               else if ( ( ptInt1 - ptPoly5).SqLen() > SQ_EPS_SMALL) {
                  ptInt2 = ptPoly5 ;
                  vtN2 = vtUD ;
                  ++ nIntNum ;
               }
            }
         }

        // Intersezione con la sesta faccia
         if ( abs( vtPoly * vtUD) > COS_ORTO_ANG_ZERO) {
            if ( nIntNum < 2  &&
                 ptPoly6.y >= dDeltaY  &&  ptPoly6.y <= dLenY + dDeltaY  &&
                 abs( ptPoly6.z) <= dLenZ + dDeltaZ) {

               if ( nIntNum == 0) {
                  ptInt1 = ptPoly6;
                  vtN1 = - vtUD ;
                  ++ nIntNum ;
               }
               else if ( ( ptInt1 - ptPoly6).SqLen() > SQ_EPS_SMALL) {
                  ptInt2 = ptPoly6;
                  vtN2 = - vtUD ;
                  ++ nIntNum ;
               }
            }
         }
        // Se il poliedro è attraversato, aggiungo
         if ( nIntNum == 2) {

           // Riporto le intersezioni nel sistema griglia
            ptInt1.ToGlob( PolyFrame) ;
            vtN1.ToGlob( PolyFrame) ;
            ptInt2.ToGlob( PolyFrame) ;
            vtN2.ToGlob( PolyFrame) ;

            AddIntervals( nGrid, i, j, ptInt1.z, ptInt2.z, vtN1, vtN2, CurrTool.GetToolNum()) ;
         }
      }
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::AddingSphere( int nGrid, const Point3d& ptS, const Point3d& ptE, double dRad)
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

   // Verifico interferisca
   int nStartI, nStartJ, nEndI, nEndJ ;
   if ( ! TestCompoBBox( nGrid, ptS, ptE, V_NULL, dRad, 0, 0, nStartI, nStartJ, nEndI, nEndJ))
      return true ;
  // Vettore movimento
   Vector3d vtV = ptE - ptS ;
   double dLengthPath = vtV.Len() ;
   if ( dLengthPath > EPS_ZERO)
      vtV /= dLengthPath ;
  // Riferimento per cilindro inviluppo della sfera lungo il movimento
   Frame3d CylFrame ;
   if ( ! CylFrame.Set( ptS, vtV))
      return false ;

   double dSqRad = dRad * dRad ;

   for ( int i = nStartI ; i <= nEndI ; ++ i) {
      for ( 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) ;

        // Sfera in posizione start
         double dStSqDXY = SqDistXY( ptC, ptS) ;
         if ( dStSqDXY < dSqRad) {
            double dMin = ptS.z - sqrt( dSqRad - dStSqDXY) ;
            Vector3d vtNmin = Point3d( dX, dY, dMin) - ptS ;
            vtNmin.Normalize() ;
            double dMax = ptS.z + sqrt( dSqRad - dStSqDXY) ;
            Vector3d vtNmax = Point3d( dX, dY, dMax) - ptS ;
            vtNmax.Normalize() ;
            AddIntervals( nGrid, i, j, dMin, dMax, vtNmin, vtNmax, CurrTool.GetToolNum()) ;
         }

        // Sfera in posizione end
         double dEnSqDXY = SqDistXY( ptC, ptE) ;
         if ( dEnSqDXY < dSqRad) {
            double dMin = ptE.z - sqrt( dSqRad - dEnSqDXY) ;
            Vector3d vtNmin = Point3d( dX, dY, dMin) - ptE ;
            vtNmin.Normalize() ;
            double dMax = ptE.z + sqrt( dSqRad - dEnSqDXY) ;
            Vector3d vtNmax = Point3d( dX, dY, dMax) - ptE ;
            vtNmax.Normalize() ;
            AddIntervals( nGrid, i, j, dMin, dMax, vtNmin, vtNmax, CurrTool.GetToolNum()) ;
         }

        // Cilindro inviluppo della sfera
         Point3d ptInt1, ptInt2 ;
         Vector3d vtN1, vtN2 ;
         if ( IntersLineCylinder( ptC, Z_AX, CylFrame, dLengthPath, dRad, false, false, ptInt1, vtN1, ptInt2, vtN2)) {
            AddIntervals( nGrid, i, j, ptInt1.z, ptInt2.z, - vtN1, - vtN2, CurrTool.GetToolNum()) ;
         }
      }
   }

   return true ;
}


// ------------------------- BOUNDING BOX --------------------------------------------------------------------------------------

//----------------------------------------------------------------------------
inline BBox3d
GetCylMoveBBox( const Point3d& ptP1, const Point3d& ptP2, const Vector3d& vtV, double dRad, double dH)
{
  // Determinazione delle posizioni della punta del componente nelle posizioni iniziale e finale
   Point3d ptP1T = ptP1 - dH * vtV ;
   Point3d ptP2T = ptP2 - dH * vtV ;
  // Calcolo del box del movimento dell'asse
   BBox3d b3Box ;
   b3Box.Add( ptP1) ;
   b3Box.Add( ptP1T) ;
   b3Box.Add( ptP2) ;
   b3Box.Add( ptP2T) ;
  // Aggiungo ingombro raggio
   if ( AreSameOrOppositeVectorApprox( vtV, X_AX))
      b3Box.Expand( 0, dRad, dRad) ;
   else if ( AreSameOrOppositeVectorApprox( vtV, Y_AX))
      b3Box.Expand( dRad, 0, dRad) ;
   else if ( AreSameOrOppositeVectorApprox( vtV, Z_AX))
      b3Box.Expand( dRad, dRad, 0) ;
   else {
      double dExpandX = dRad * sqrt( 1 - vtV.x * vtV.x) ;
      double dExpandY = dRad * sqrt( 1 - vtV.y * vtV.y) ;
      double dExpandZ = dRad * sqrt( 1 - vtV.z * vtV.z) ;
      b3Box.Expand( dExpandX, dExpandY, dExpandZ) ;
   }
  // Restituisco il box calcolato
   return b3Box ;
}

//----------------------------------------------------------------------------
inline BBox3d
GetSphereMoveBBox( const Point3d& ptP1, const Point3d& ptP2, double dRad)
{
  // Calcolo del box del movimento del centro
   BBox3d b3Box ;
   b3Box.Add( ptP1) ;
   b3Box.Add( ptP2) ;
  // Aggiungo ingombro raggio
   b3Box.Expand( dRad) ;
  // Restituisco il box calcolato
   return b3Box ;
}

//----------------------------------------------------------------------------
inline bool
VolZmap::TestToolBBox( int nGrid, const Point3d& ptP1, const Point3d& ptP2, const Vector3d& vtV,
                       int& nStI, int& nStJ, int& nEnI, int& nEnJ)
{
  // Controllo utensile
   if ( m_nCurrTool < 0  ||  m_nCurrTool >= int( m_vTool.size()))
      return false ;
   Tool& CurrTool = m_vTool[m_nCurrTool] ;

   return TestCompoBBox( nGrid, ptP1, ptP2, vtV, CurrTool.GetRadius(), CurrTool.GetTipRadius(), CurrTool.GetHeigth(),
                         nStI, nStJ, nEnI, nEnJ) ;
}

//----------------------------------------------------------------------------
inline bool
VolZmap::TestCompoBBox( int nGrid, const Point3d& ptP1, const Point3d& ptP2, const Vector3d& vtV,
                        double dRad, double dTipRad, double dHei,
                        int& nStI, int& nStJ, int& nEnI, int& nEnJ)
{
  // I punti e i vettori devono essere nel sistema di riferimento opportuno

  // Controllo sull'ammissibilità del numero di griglia
   if ( nGrid < 0  ||  nGrid > 2)
      return false ;

  // BBox dello Zmap
   BBox3d b3Zmap( 0, 0, m_dMinZ[nGrid], m_nNx[nGrid] * m_dStep, m_nNy[nGrid] * m_dStep, m_dMaxZ[nGrid]) ;

  // BBox dell'utensile nel suo movimento
   double dMaxRad = max( dRad, dTipRad) ;
   BBox3d b3Box = ( vtV.IsSmall() ? GetSphereMoveBBox( ptP1, ptP2, dRad) : GetCylMoveBBox( ptP1, ptP2, vtV, dMaxRad, dHei)) ;

  // Verifica dell'interferenza dell'utensile con lo Zmap
   if ( ! b3Zmap.FindIntersection( b3Box, b3Box))
      return false ;

  // Limiti su indici
   nStI = max( 0, int( b3Box.GetMin().x / m_dStep)) ;
   nEnI = min( m_nNx[nGrid] - 1, int( b3Box.GetMax().x / m_dStep)) ;
   nStJ = max( 0, int( b3Box.GetMin().y / m_dStep)) ;
   nEnJ = min( m_nNy[nGrid] - 1, int( b3Box.GetMax().y / m_dStep)) ;
   
   return true ;
}

//----------------------------------------------------------------------------
inline bool 
VolZmap::TestParaBBox( int nGrid, const Point3d& ptS, const Point3d& ptE, const Vector3d& vtD, const Vector3d& vtA,
                       double dLenX, double dLenY, double dLenZ,
                       int& nStI, int& nStJ, int& nEnI, int& nEnJ) 
{
  // I punti e i vettori devono essere nel sistema di riferimento opportuno

  // Determino le posizioni iniziale e  
  // finale della punta dell'utensile.
   Point3d ptSTip = ptS - vtD * dLenZ ;
   Point3d ptETip = ptE - vtD * dLenZ ; 
 
   double dSemiDiag = sqrt( dLenX * dLenX + dLenY * dLenY) / 2 ;

  // Determinazione dei limiti del più piccolo parallelepipedo contenente il movimento
   double dMinX = min( min( ptS.x, ptSTip.x), min( ptE.x, ptETip.x)) - dSemiDiag ;
   double dMinY = min( min( ptS.y, ptSTip.y), min( ptE.y, ptETip.y)) - dSemiDiag ;
   double dMinZ = min( min( ptS.z, ptSTip.z), min( ptE.z, ptETip.z)) - dSemiDiag ;
   double dMaxX = max( max( ptS.x, ptSTip.x), max( ptE.x, ptETip.x)) + dSemiDiag ;
   double dMaxY = max( max( ptS.y, ptSTip.y), max( ptE.y, ptETip.y)) + dSemiDiag ;
   double dMaxZ = max( max( ptS.z, ptSTip.z), max( ptE.z, ptETip.z)) + dSemiDiag ;

   int nMaxNx = m_nNx[nGrid] ;
   int nMaxNy = m_nNy[nGrid] ;

   double dMaxXValue = nMaxNx * m_dStep ;
   double dMaxYValue = nMaxNy * m_dStep ;
   double dMinZValue = m_dMinZ[nGrid] ;
   double dMaxZValue = m_dMaxZ[nGrid] ;

   // Verifica dell'interferenza dell'utensile con lo Zmap
   if ( dMaxX < EPS_SMALL  ||  dMinX > dMaxXValue - EPS_SMALL)
      return false ;
   if ( dMaxY < EPS_SMALL  ||  dMinY > dMaxYValue - EPS_SMALL)
      return false ;
   if ( dMaxZ < dMinZValue + EPS_SMALL  ||  dMinZ > dMaxZValue - EPS_SMALL)
      return false ;

  // Limiti su indici
   nStI = ( dMinX < EPS_SMALL ? 0 : int( dMinX / m_dStep)) ;
   nEnI = ( dMaxX > dMaxXValue - EPS_SMALL ? nMaxNx - 1 : int( dMaxX / m_dStep)) ;
   nStJ = ( dMinY < EPS_SMALL ? 0 : int( dMinY / m_dStep)) ;
   nEnJ = ( dMaxY > dMaxYValue - EPS_SMALL ? nMaxNy - 1 : int ( dMaxY / m_dStep)) ;

   return true ;
}