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

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

#include "stdafx.h"
#include "CurveLine.h"
#include "VolZmap.h"
#include "GeoConst.h"
#include "IntersLineSurfTm.h"
#include "\EgtDev\Include\EgtNumUtils.h"


using namespace std ;

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

   dU1 = - INFINITO ;
   dU2 = INFINITO ;

  // confronto con piani YZ (perpendicolari ad asse X)
   if ( vtV.x > EPS_ZERO) {
      dU1 = max( dU1, ( ptMin.x - ptP.x) / vtV.x) ;
      dU2 = min( dU2, ( ptMax.x - ptP.x) / vtV.x) ;
   }
   else if ( vtV.x < - EPS_ZERO) {
      dU1 = max( dU1, ( ptMax.x - ptP.x) / vtV.x) ;
      dU2 = min( dU2, ( ptMin.x - ptP.x) / vtV.x) ;
   }
   else if ( ptP.x < ptMin.x - EPS_SMALL  ||  ptP.x > ptMax.x + EPS_SMALL)
      return false ;

  // confronto con piani ZX (perpendicolari ad asse Y)
   if ( vtV.y > EPS_ZERO) {
      dU1 = max( dU1, ( ptMin.y - ptP.y) / vtV.y) ;
      dU2 = min( dU2, ( ptMax.y - ptP.y) / vtV.y) ;
   }
   else if ( vtV.y < - EPS_ZERO) {
      dU1 = max( dU1, ( ptMax.y - ptP.y) / vtV.y) ;
      dU2 = min( dU2, ( ptMin.y - ptP.y) / vtV.y) ;
   }
   else if ( ptP.y < ptMin.y - EPS_SMALL  ||  ptP.y > ptMax.y + EPS_SMALL)
      return false ;

  // confronto con piani XZ (perpendicolari ad asse Z)
   if ( vtV.z > EPS_ZERO) {
      dU1 = max( dU1, ( ptMin.z - ptP.z) / vtV.z) ;
      dU2 = min( dU2, ( ptMax.z - ptP.z) / vtV.z) ;
   }
   else if ( vtV.z < - EPS_ZERO) {
      dU1 = max( dU1, ( ptMax.z - ptP.z) / vtV.z) ;
      dU2 = min( dU2, ( ptMin.z - ptP.z) / vtV.z) ;
   }
   else if ( ptP.z < ptMin.z - EPS_SMALL  ||  ptP.z > ptMax.z + EPS_SMALL)
      return false ;

   return ( dU2 >= dU1) ;
}

//----------------------------------------------------------------------------
bool 
VolZmap::IntersLineVoxel( const Point3d& ptP, const Vector3d& vtV, int nIndI, int nIndJ, int nIndK,
                          int& nFaceF, int& nFaceL, double& dUF, double& dUL) const 
{
  // Controllo sull'ammissibilità del voxel
   if ( nIndI < - 1 || nIndI >= int( m_nNx[0]) || 
        nIndJ < - 1 || nIndJ >= int( m_nNy[0]) ||
        nIndK < - 1 || nIndK >= int( m_nNy[1]))
        return false ;

   Point3d ptInt, ptIntF, ptIntL ;
   double dSqEps = EPS_SMALL * EPS_SMALL ;
   int nIntNum  = 0 ;
   
   double dU1 = ( ( nIndJ + 0.5) * m_dStep - ptP.y) / vtV.y ;
   double dU2 = ( ( nIndI + 0.5) * m_dStep - ptP.x) / vtV.x ;
   double dU3 = ( ( nIndJ + 1.5) * m_dStep - ptP.y) / vtV.y ;
   double dU4 = ( ( nIndI + 1.5) * m_dStep - ptP.x) / vtV.x ;
   double dU5 = ( ( nIndK + 0.5) * m_dStep - ptP.z) / vtV.z ;
   double dU6 = ( ( nIndK + 1.5) * m_dStep - ptP.z) / vtV.z ;
   
  // Intersezione con le facce 1 e 3
   if ( abs( vtV.y) > EPS_ZERO) {
   
     // Intersezione con la prima faccia
      ptInt = ptP + dU1 * vtV ;

      if ( ptInt.x >= ( nIndI + 0.5) * m_dStep && 
           ptInt.x <= ( nIndI + 1.5) * m_dStep &&
           ptInt.z >= ( nIndK + 0.5) * m_dStep && 
           ptInt.z <= ( nIndK + 1.5) * m_dStep) {
         
         dUF = dU1 ;
         nFaceF = 1 ;
         ptIntF = ptInt ;    
         ++ nIntNum ;
      }
    
     // Intersezione con la terza faccia
      ptInt = ptP + dU3 * vtV ;

      if ( ptInt.x >= ( nIndI + 0.5) * m_dStep && 
           ptInt.x <= ( nIndI + 1.5) * m_dStep &&
           ptInt.z >= ( nIndK + 0.5) * m_dStep && 
           ptInt.z <= ( nIndK + 1.5) * m_dStep) {
           
         if ( nIntNum == 0) {
         
            dUF = dU3 ;
            nFaceF = 3 ;
            ptIntF = ptInt ;
            ++ nIntNum ;
         }
         else if ( nIntNum == 1 && 
                 ( ptIntF - ptInt).SqLen() > dSqEps) {
         
            dUL = dU3 ;
            nFaceL = 3 ;
            ptIntL = ptInt ;
            ++ nIntNum ;
         }  
      }
   }

  // Intersezione con le facce 2 e 4
   if ( abs( vtV.x) > EPS_ZERO) {
  
     // Intersezione con la seconda faccia
      ptInt = ptP + dU2 * vtV ;

      if ( ptInt.y >= ( nIndJ + 0.5) * m_dStep && 
           ptInt.y <= ( nIndJ + 1.5) * m_dStep &&
           ptInt.z >= ( nIndK + 0.5) * m_dStep && 
           ptInt.z <= ( nIndK + 1.5) * m_dStep) {
           
         if ( nIntNum == 0) {
            
            dUF = dU2 ;
            nFaceF = 2 ;
            ptIntF = ptInt ;
            ++ nIntNum ;
         }
         else if ( nIntNum == 1 &&
                 ( ptIntF - ptInt).SqLen() > dSqEps) {
         
            dUL = dU2 ;
            nFaceL = 2 ;
            ptIntL = ptInt ;
            ++ nIntNum ;
         }
      }

     // Intersezione con la quarta faccia
      ptInt = ptP + dU4 * vtV ;

      if ( ptInt.y >= ( nIndJ + 0.5) * m_dStep && 
           ptInt.y <= ( nIndJ + 1.5) * m_dStep &&
           ptInt.z >= ( nIndK + 0.5) * m_dStep && 
           ptInt.z <= ( nIndK + 1.5) * m_dStep) { 
           
         if ( nIntNum == 0) {
            
            dUF = dU4 ;
            nFaceF = 4 ;
            ptIntF = ptInt ;
            ++ nIntNum ;
         }
         else if ( nIntNum == 1 &&
                 ( ptIntF - ptInt).SqLen() > dSqEps) {
         
            dUL = dU4 ;
            nFaceL = 4 ;
            ptIntL = ptInt ;
            ++ nIntNum ;
         }    
      }
   }

  // Intersezione con le facce 5 e 6
   if ( abs( vtV.z) > EPS_ZERO) {
   
     // Intersezione con la quinta faccia
      ptInt = ptP + dU5 * vtV ;

      if ( ptInt.x >= ( nIndI + 0.5) * m_dStep && 
           ptInt.x <= ( nIndI + 1.5) * m_dStep &&
           ptInt.y >= ( nIndJ + 0.5) * m_dStep && 
           ptInt.y <= ( nIndJ + 1.5) * m_dStep) {
      
         if ( nIntNum == 0) {
            
            dUF = dU5 ;
            nFaceF = 5 ;
            ptIntF = ptInt ;
            ++ nIntNum ;
         }
         else if ( nIntNum == 1 &&
                 ( ptIntF - ptInt).SqLen() > dSqEps) {
         
            dUL = dU5 ;
            nFaceL = 5 ;
            ptIntL = ptInt ;
            ++ nIntNum ;
         }    
      } 
           
     // Intersezione con la sesta faccia
      ptInt = ptP + dU6 * vtV ;

      if ( ptInt.x >= ( nIndI + 0.5) * m_dStep && 
           ptInt.x <= ( nIndI + 1.5) * m_dStep &&
           ptInt.y >= ( nIndJ + 0.5) * m_dStep && 
           ptInt.y <= ( nIndJ + 1.5) * m_dStep) {
      
      
          if ( nIntNum == 0) {
            
            dUF = dU6 ;
            nFaceF = 6 ;
            ptIntF = ptInt ;
            ++ nIntNum ;
         }
         else if ( nIntNum == 1 &&
                 ( ptIntF - ptInt).SqLen() > dSqEps) {
         
            dUL = dU6 ;
            nFaceL = 6 ;
            ptIntL = ptInt ;
            ++ nIntNum ;
         }  
      }            
   }

   if ( dUF > dUL) {
   
      swap( dUF, dUL) ;
      swap( nFaceF, nFaceL) ;
   }
   return ( nIntNum == 2) ;
} 

//----------------------------------------------------------------------------
bool
VolZmap::IntersLineZMapBBox( unsigned int nGrid, const Point3d& ptP, const Vector3d& vtV, double& dU1, double& dU2) 
{
  // Punti estremi del box dello Zmap
   Point3d ptMin = ORIG ;
   Point3d ptMax = ptMin + Vector3d( m_nNx[nGrid] * m_dStep, m_nNy[nGrid] * m_dStep, m_dMaxZ[nGrid]) ;

   return ( IntersLineBox( ptP, vtV, ptMin, ptMax, dU1, dU2)  &&  ( dU1 > 0  ||  dU2 > 0)) ;
}

//----------------------------------------------------------------------------
bool 
VolZmap::IntersLineDexel( unsigned int nGrid, const Point3d& ptP, const Vector3d& vtV, unsigned int nI,
                          unsigned int nJ, double& dU1, double& dU2) 
{ 
  // Determino l'indice del dexel e il doppio del numero di suo intervalli
   unsigned int nDexelPos = nJ * m_nNx[nGrid] + nI ;
   unsigned int nDexelSize = unsigned int( m_Values[nGrid][nDexelPos].size()) ;

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

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

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

   return bInters ;
}

//----------------------------------------------------------------------------
bool
VolZmap::GetDepth( const Point3d& ptPGlob, const Vector3d& vtDir, double& dInLength, double& dOutLength)
{
  // Porto il raggio nel riferimento intrinseco
   Point3d ptP = ptPGlob ;
   ptP.ToLoc( m_MapFrame[0]) ;
   Vector3d vtV = vtDir ;
   vtV.ToLoc( m_MapFrame[0]) ;
   vtV.Normalize() ;
 
  // Studio dell'intersezione fra semiretta e BBox dello Zmap
   double dU1, dU2 ;
   bool bTest = IntersLineZMapBBox( 0, ptP, vtV, dU1, dU2) ;

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

   Point3d ptI, ptF ;
  // Una sola intersezione valida ( punto interno, intersezione valida 2)
   if ( dU1 < 0  &&  dU2 > 0) {
      ptI = ptP ;
      ptF = ptP + dU2 * vtV ;
   }
  // due soluzioni valide ( punto esterno)
   else {
      ptI = ptP + dU1 * vtV ;
      ptF = ptP + dU2 * vtV ;
   }

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

  // Variazioni
   double dDeltaX = ptF.x - ptI.x ;  
   double dDeltaY = ptF.y - ptI.y ;

  // se inclinazione da asse X minore di 45 gradi (in assoluto)
   if ( abs( dDeltaY) <= abs( dDeltaX)) {
     // mi muovo lungo X (i)
      int nIncrI = ( nFi >= nIi ? 1 : - 1) ;
      for ( int i = nIi, j = nIj ;
            i != nFi + nIncrI ;
            i += nIncrI) {

        // Controllo con nuovo i e j corrente (considero il bordo sinistro del dexel)
         double dU1, dU2 ;
         if ( IntersLineDexel( 0, ptP, vtV, i, j, dU1, dU2)) {
            dInLength = min( dInLength, dU1) ;
            dOutLength = max( dOutLength, dU2) ;
         }
           
        // Mi sposto sul bordo destro del dexel
         double dMoveX = ( ( i + max( nIncrI, 0)) * m_dStep - ptI.x) ;
         double dMoveY = dMoveX * dDeltaY / dDeltaX ;
         double dY = ptI.y + dMoveY ;
         int OldJ = j ;
         j = Clamp( int( floor( dY / m_dStep)), 0, m_nNy[0] - 1) ;

        // Analisi del dexel
         if ( j != OldJ) {
            double dU1, dU2 ;
            if ( IntersLineDexel( 0, ptP, vtV, i, j, dU1, dU2)) {
               dInLength = min( dInLength, dU1) ;
               dOutLength = max( dOutLength, dU2) ;
            }
         }
      }
   }

  // altrimenti
   else {
     // mi muovo lungo Y (j)
      int nIncrJ = ( nFj >= nIj ? 1 : - 1) ;
      for ( int i = nIi, j = nIj ;
            j != nFj + nIncrJ ;
            j += nIncrJ) {

        // Controllo con nuovo j e i corrente (considero il bordo sotto del dexel)
         double dU1, dU2 ;
         if ( IntersLineDexel( 0, ptP, vtV, i, j, dU1, dU2)) {
            dInLength = min( dInLength, dU1) ;
            dOutLength = max( dOutLength, dU2) ;
         }
           
        // Mi sposto sul bordo sopra del dexel
         double dMoveY = ( ( j + max( nIncrJ, 0)) * m_dStep - ptI.y) ;
         double dMoveX = dMoveY * dDeltaX / dDeltaY ;
         double dX = ptI.x + dMoveX ;
         int OldI = i ;
         i = Clamp( int( floor( dX / m_dStep)), 0, m_nNx[0] - 1) ;

       // Analisi del dexel
         if ( i != OldI) {
            double dU1, dU2 ;
            if ( IntersLineDexel( 0, ptP, vtV, i, j, dU1, dU2)) {
               dInLength = min( dInLength, dU1) ;
               dOutLength = max( dOutLength, dU2) ;
            }
         }
      }   
   }

  // Se non abbiamo incontrato materiale
   if ( dInLength > dOutLength - EPS_SMALL) {
      dInLength = - 2 ;
      dOutLength = - 2 ;
      return true ;
   }

  // Se parto dall'interno
   if ( dInLength < - EPS_SMALL)
      dInLength = - 1 ;
     
   return true ;
}

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

  // lo porto nel riferimento intrinseco dello Zmap
   b3Box.LocToLoc( frBox, m_MapFrame[0]) ;

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

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

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

  // Vettore direzione dei dexel nel riferimento del Box
   Vector3d vtK = Z_AX ;  vtK.LocToLoc( m_MapFrame[0], frBox) ;

  // Riferimento intrinseco dei dexel nel riferimento del box
   Point3d ptO = ORIG ;   ptO.LocToLoc( m_MapFrame[0], frBox) ;
   Vector3d vtX = X_AX ;  vtX.LocToLoc( m_MapFrame[0], frBox) ;
   Vector3d vtY = Y_AX ;  vtY.LocToLoc( m_MapFrame[0], frBox) ;

  // Ciclo di intersezione dei dexel con il BBox
   for ( int i = nStI ; i <= nEnI ; ++ i) {

      for ( int j = nStJ ; j <= nEnJ ; ++ j) {

         int nPos = j * m_nNx[0] + i ;
         int nSize = int( m_Values[0][nPos].size()) ;
         if ( nSize == 0)
            continue ;

         Point3d ptC = ptO + ( i + 0.5) * m_dStep * vtX + ( j + 0.5) * m_dStep * vtY ;

         double dZmin, dZmax ;
         if ( IntersLineBox( ptC, vtK, ORIG, ORIG + vtDiag, dZmin, dZmax)) {

            for ( int nIndex = 0 ; nIndex < nSize ; nIndex += 2) {
               if ( ! ( dZmax < m_Values[0][nPos][nIndex].dZVal - EPS_SMALL  ||
                        dZmin > m_Values[0][nPos][nIndex + 1].dZVal + EPS_SMALL))
                  return false ;
            }
         }
      }
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
VolZmap::IntersLineCylinder( const Point3d& ptLineSt, const Vector3d& vtLineDir,
                             const Frame3d& CylFrame, double dL, double dR,
                             Point3d& ptInt1, Point3d& ptInt2, Vector3d& vtN1, Vector3d& vtN2,
                             bool bTapO, bool bTapL) 
{  
  // NB: L'origine del sistema di riferimento deve essere 
  // nel centro della circonferenza di base e l'asse di simmetria 
  // deve coincidere con l'asse x. 
  // La funzione restituisce true in caso di intersezione,
  // false altrimenti.

   Point3d ptP = ptLineSt ;
   Vector3d vtV = vtLineDir ; 

  // Trasformazione delle coordinate:
  // l'asse del cilindro corrisponde con
  // l'asse x del sistema di riferimento
   ptP.ToLoc( CylFrame) ;
   vtV.ToLoc( CylFrame) ;

   DBLVECTOR vdCoef(3) ;
   DBLVECTOR vdRoots ;

   double dSqRad = dR * dR - 2 * dR * EPS_SMALL ;

   vdCoef[0] = ptP.y * ptP.y + ptP.z * ptP.z - dSqRad ;
   vdCoef[1] = 2 * ( ptP.y * vtV.y + ptP.z * vtV.z) ;
   vdCoef[2] = vtV.y * vtV.y + vtV.z * vtV.z ;

  // Computo radici
   int nRoot = PolynomialRoots( 2, vdCoef, vdRoots) ;

  // Nessuna soluzione 
   if ( nRoot == 0 || nRoot == 1) {

      if ( abs( vtV.x) > EPS_ZERO) {
      
         ptInt1 = ptP - ( ptP.x / vtV.x) * vtV ;
         ptInt2 = ptP + ( ( dL - ptP.x) / vtV.x) * vtV ;

         vtN1 = X_AX ;
         vtN2 = - X_AX ;

         if ( ptInt1.y * ptInt1.y + ptInt1.z * ptInt1.z <= dSqRad  &&
              ptInt2.y * ptInt2.y + ptInt2.z * ptInt2.z <= dSqRad) {

            ptInt1.ToGlob( CylFrame) ;
            ptInt2.ToGlob( CylFrame) ;

            vtN1.ToGlob( CylFrame) ;
            vtN2.ToGlob( CylFrame) ;
           
            return true ;
         }
        // Nessuna intersezione
         else 
            return false ;
      }
     // Nessuna intersezione
      else 
         return false ;
   }

  // L'equazione ammette o due soluzioni (eventualmente
  // coincidenti) oppure nessuna o infinite se la la retta
  // appartiene alla superficie

   if ( nRoot == 2) {

      double dEpsO = ( bTapO ? - EPS_SMALL : EPS_SMALL) ;
      double dEpsL = ( bTapL ? EPS_SMALL : - EPS_SMALL) ;

      ptInt1 = ptP + vdRoots[0] * vtV ;
      ptInt2 = ptP + vdRoots[1] * vtV ;

      if ( ptInt1.x > ptInt2.x)
         swap( ptInt1, ptInt2) ; 

      vtN1.Set( 0, ( ORIG - ptInt1).y, ( ORIG - ptInt1).z) ;
      vtN2.Set( 0, ( ORIG - ptInt2).y, ( ORIG - ptInt2).z) ;
      

      if ( ptInt1.x < dL + dEpsL) {
      
         if (  ptInt1.x > dEpsO) {

            if ( ptInt2.x > dL + dEpsL) {

               ptInt2 = ptP + ( ( dL - ptP.x) / vtV.x) * vtV ;  
               vtN2.Set( -1, 0, 0) ;
            }
         }
         else {

            if ( ptInt2.x > dL + dEpsL) {
            
               ptInt1 = ptP - ( ptP.x / vtV.x) * vtV ;
               ptInt2 = ptP + ( ( dL - ptP.x) / vtV.x) * vtV ;
               vtN1.Set( 1, 0, 0) ;
               vtN2.Set( -1, 0, 0) ;
            }
            else if ( ptInt2.x > dEpsO) {
            
               ptInt1 = ptP - ( ptP.x / vtV.x) * vtV ;                 
               vtN1.Set( 1, 0, 0) ;
            }
            else 
               return false ;
         }
      }
      else 
         return false ;

     // Riporto le coordinate nel sistema di riferimento griglia
      ptInt1.ToGlob( CylFrame) ;
      ptInt2.ToGlob( CylFrame) ;
      vtN1.ToGlob( CylFrame) ;
      vtN2.ToGlob( CylFrame) ;

      vtN1.Normalize() ;
      vtN2.Normalize() ;
   }                  

   return true ;
}

//----------------------------------------------------------------------------
bool 
VolZmap::IntersZLineCylinder( const Point3d& ptLine,
                              const Point3d& ptBase, const Point3d& ptTop, double dCylR,
                              double& dInfZ, double& dSupZ) 
{
  // NB: Le coordinate sono espresse nel sistema griglia
  // La funzione restituisce true in caso di intersezione,
  // false altrimenti.

   double dSqRad = dCylR * dCylR ;

  // Cilindro verticale 
   if ( AreSamePointXYApprox( ptBase, ptTop)) {

     // Intersezione
      if ( SqDistXY( ptLine, ptBase) < dSqRad) {
         dInfZ = min( ptBase.z, ptTop.z) ;
         dSupZ = max( ptBase.z, ptTop.z) ;
         return true ;
      }

     // Non vi è intersezione
      else 
         return false ;
   }

  // Cilindro non verticale
   else {

     // Studio delle simmetrie
      Point3d ptS = ( ptBase.z < ptTop.z ? ptBase : ptTop) ;
      Point3d ptE = ( ptBase.z < ptTop.z ? ptTop : ptBase) ;

      Vector3d vtAx = ptE - ptS ;  

      Vector3d vtV1( vtAx.x, vtAx.y, 0) ;

      double dLenXY = vtV1.LenXY() ;
      double dSZ = ptS.z ;
      double dEZ = ptE.z ;
      double dDeltaZ = dEZ - dSZ ;

      Vector3d vtL( ptLine.x - ptS.x, ptLine.y - ptS.y, 0) ;

     // vtV1 e vtV2 formano un sistema ortonormale 
     // sul piano e insieme a ptSxy formano un sistema
     // di riferimento bidimensionale
      vtV1.Normalize() ;
      Vector3d vtV2 = vtV1 ;
      vtV2.Rotate( Z_AX, 90) ;

      double dLen = vtAx.Len() ;

     // Sono seno e coseno dell'angolo complementare
     // rispetto a quello formato dal vettore movimento
     // con il piano, per questo motivo si ha dCos con 
     // dDeltaZ e dSin con dLenXY
      double dCos = dDeltaZ / dLen ;
      double dSin = dLenXY / dLen ;

     // Nuove coordinate piane del punto
      double dLocX1 = vtL * vtV1 ;
      double dLocX2 = vtL * vtV2 ;

      double dSqRoot = sqrt( dSqRad - dLocX2 * dLocX2) ;
      double dX1_0 = dCos * dSqRoot ;

      if ( dLocX1 >= - dX1_0  &&  dLocX1 <= dLenXY + dX1_0  &&
           abs( dLocX2) < dCylR) {

        // Minimi
         if ( dLocX1 < dX1_0) {

            double dDotS = vtAx * ( ptS - ORIG) ;
           // Qui usiamo ptLine perché servono coordinate griglia
            dInfZ = ( dDotS - vtAx.x * ptLine.x - vtAx.y * ptLine.y) / vtAx.z ;
         }
         else {

            double dZ0 = - dSin * dSqRoot ;
                          
            dInfZ = dSZ + dZ0 + ( dLocX1 - dX1_0) * dDeltaZ / dLenXY ;
         } 
            
        // Massimi
         if ( dLocX1 <= dLenXY - dX1_0)  {

            double dZ0 = dSin * dSqRoot ;

            dSupZ = dSZ + dZ0 + ( dLocX1 + dX1_0) * dDeltaZ / dLenXY ;       
         }
         else {

            double dDotE = vtAx * ( ptE - ORIG) ;
           // Qui usiamo ptLine perché servono coordinate griglia
            dSupZ = ( dDotE - vtAx.x * ptLine.x - vtAx.y * ptLine.y) / vtAx.z ;
         }

         return true ;
      } 

      return false ;
   }
}

//----------------------------------------------------------------------------
bool
VolZmap::IntersLineConus( const Point3d& ptLineSt, const Vector3d& vtLineDir,
                          const Frame3d& ConusFrame, double dTan, double dl, double dL, 
                          Point3d& ptInt1, Point3d& ptInt2, Vector3d& vtN1, Vector3d& vtN2,
                          bool bTapLow, bool bTapUp) 
{
  // NB: L'origine del sistema di riferimento deve essere 
  // nel vertice del cono e l'asse di simmetria deve coincidere 
  // con l'asse x.
  // La funzione restituisce true in caso di intersezione,
  // false altrimenti.

   Point3d ptP = ptLineSt ;
   Vector3d vtV = vtLineDir ; 

  // Trasformazione delle coordinate
   ptP.ToLoc( ConusFrame) ;
   vtV.ToLoc( ConusFrame) ;

   DBLVECTOR vdCoef(3) ;
   DBLVECTOR vdRoots ;

   double dSqTan = dTan * dTan ;
   double dMinRad = dTan * dl ;
   double dMaxRad = dTan * dL ;
   double dDeltaR = dMaxRad - dMinRad ;
   double dHei = dL - dl ;
          
   vdCoef[0] = dSqTan * ptP.x * ptP.x - ptP.y * ptP.y - ptP.z * ptP.z ;
   vdCoef[1] = 2 * ( dSqTan * ptP.x * vtV.x - ptP.y * vtV.y - ptP.z * vtV.z) ;
   vdCoef[2] = dSqTan * vtV.x * vtV.x - vtV.y * vtV.y - vtV.z * vtV.z ;

  // Computo radici
   int nRoot = PolynomialRoots( 2, vdCoef, vdRoots) ;
  
  // Nessuna soluzione 
   if ( nRoot == 0)
      return false ;

   double dEpsLow = ( bTapLow ? - EPS_SMALL : EPS_SMALL) ;
   double dEpsUp = ( bTapUp ? EPS_SMALL : - EPS_SMALL) ;
   
  // Una soluzione: la retta iterseca superficie 
  // laterale e un piano
   if ( nRoot == 1) {

      ptInt1 = ptP + vdRoots[0] * vtV ;
      
      Vector3d vtU = ( ptInt1 - ORIG) - ( ptInt1 - ORIG).x * X_AX ;

      vtU.Normalize() ;

      vtN1 = dDeltaR * X_AX - dHei * vtU ;

      vtN1.Normalize() ;

      if ( ptInt1.x < dL + dEpsUp) {

         if ( ptInt1.x > dl + dEpsLow) {
      
            ptInt2 = ptP + ( ( dL - ptP.x) / vtV.x) * vtV ;

            vtN2 = - X_AX ;

         }
         else if ( ptInt1.x > - EPS_SMALL) {

            ptInt1 = ptP + ( ( dl - ptP.x) / vtV.x) * vtV ;
            ptInt2 = ptP + ( ( dL - ptP.x) / vtV.x) * vtV ;  

            vtN1 = X_AX ;
            vtN2 = - X_AX ;

            if ( ptInt2.y * ptInt2.y + ptInt2.z * ptInt2.z > dMaxRad * dMaxRad)

                 return false ;
         }
         else 
            return false ;

         ptInt1.ToGlob( ConusFrame) ;
         ptInt2.ToGlob( ConusFrame) ;

         vtN1.ToGlob( ConusFrame) ;
         vtN2.ToGlob( ConusFrame) ;

         return true ;
      }
      else 
         return false ;
      
   }
  // Due soluzioni: la retta interseca due volte la 
  // superficie laterale
   else if ( nRoot == 2) {

      ptInt1 = ptP + vdRoots[0] * vtV ;
      ptInt2 = ptP + vdRoots[1] * vtV ;

      if ( ptInt1.x > ptInt2.x) 
         swap( ptInt1, ptInt2) ;

      Vector3d vtU1 = ( ptInt1 - ORIG) - ( ptInt1 - ORIG).x * X_AX ;
      Vector3d vtU2 = ( ptInt2 - ORIG) - ( ptInt2 - ORIG).x * X_AX ;

      vtU1.Normalize() ;
      vtU2.Normalize() ;

      vtN1 = dDeltaR * X_AX - dHei * vtU1 ;
      vtN2 = dDeltaR * X_AX - dHei * vtU2 ;

      vtN1.Normalize() ;
      vtN2.Normalize() ;

      if ( abs( vtV.x) < EPS_ZERO) {

         if ( ptInt1.x > dl + dEpsLow && ptInt1.x < dL + dEpsUp) {
         
            ptInt1.ToGlob( ConusFrame) ;
            ptInt2.ToGlob( ConusFrame) ;

            vtN1.ToGlob( ConusFrame) ;
            vtN2.ToGlob( ConusFrame) ;

            vtN1.Normalize() ;
            vtN2.Normalize() ;

            return true ;
         }
         else 
            return false ;     
      }
     

      if ( ptInt1.x < dL + dEpsUp) {
      
         if ( ptInt1.x > dl + dEpsLow) {

            if ( ptInt2.x > dL + dEpsUp) { 

               ptInt2 = ptP + ( ( dL - ptP.x) / vtV.x) * vtV ; 

               vtN2 = - X_AX ;
            }
         }
         else if ( ptInt1.x > - EPS_SMALL) {

            if ( ptInt2.x > dL + dEpsUp) {

               ptInt1 = ptP + ( ( dl - ptP.x) / vtV.x) * vtV ;
               ptInt2 = ptP + ( ( dL - ptP.x) / vtV.x) * vtV ;  

               vtN1 = X_AX ;
               vtN2 = - X_AX ;
            }
            else if ( ptInt2.x > dl + dEpsLow) { 

               ptInt1 = ptP + ( ( dl - ptP.x) / vtV.x) * vtV ;
               vtN1 = X_AX ;
            }
            else 
               return false ;
         }
         else {

            if ( ptInt2.x < 0) 

               return false ;

            else if ( ptInt2.x < dl + dEpsLow) {

               ptInt1 = ptP + ( ( dl - ptP.x) / vtV.x) * vtV ;
               ptInt2 = ptP + ( ( dL - ptP.x) / vtV.x) * vtV ; 

               vtN1 = X_AX ;
               vtN2 = - X_AX ;
            }          
            else if ( ptInt2.x < dL + dEpsUp) { 

               ptInt1 = ptP + ( ( dL - ptP.x) / vtV.x) * vtV ;
               vtN1 = - X_AX ;
            }
            else 
               return false ;
         }

         ptInt1.ToGlob( ConusFrame) ;
         ptInt2.ToGlob( ConusFrame) ;

         vtN1.ToGlob( ConusFrame) ;
         vtN2.ToGlob( ConusFrame) ;

         vtN1.Normalize() ;
         vtN2.Normalize() ;

         return true ;
      }
      else
         return false ;
   }
   return false ;
}

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

   double dObCoef = dOrtMvLen / dLongMvLen ;
   double dSqCoef = dObCoef * dObCoef ;

   Point3d ptP = ptLineSt ;
   Vector3d vtV = vtLineDir ; 

  // Asse cilindro ellittico
   Vector3d vtAx( dLongMvLen, dOrtMvLen, 0) ;
   vtAx.Normalize() ;

  // Trasformazione delle coordinate
   ptP.ToLoc( CircFrame) ;
   vtV.ToLoc( CircFrame) ;

   std::vector <double> vdCoef(3) ;
   std::vector <double> vdRoots ;

   vdCoef[0] = dSqCoef * ptP.x * ptP.x + ptP.y * ptP.y + ptP.z * ptP.z - 2 * dObCoef * ptP.x * ptP.y - dSqRad ;
   vdCoef[1] = 2 * ( dSqCoef * vtV.x * ptP.x + vtV.y * ptP.y + vtV.z * ptP.z - dObCoef * ( vtV.x * ptP.y + vtV.y * ptP.x)) ;
   vdCoef[2] = dSqCoef * vtV.x * vtV.x + vtV.y * vtV.y + vtV.z * vtV.z - 2 * dObCoef * vtV.x * vtV.y ;

   int nRoot = PolynomialRoots( 2, vdCoef, vdRoots) ;


   // Nessuna soluzione 
   if ( nRoot == 0 || nRoot == 1) {

      if ( abs( vtV.x) > EPS_ZERO) {
      
         ptInt1 = ptP - ( ptP.x / vtV.x) * vtV ;
         ptInt2 = ptP + ( ( dLongMvLen - ptP.x) / vtV.x) * vtV ;

         if ( ptInt1.y * ptInt1.y + ptInt1.z * ptInt1.z < dSqRad &&
            ( ptInt2.y - dOrtMvLen) * ( ptInt2.y - dOrtMvLen) + ptInt2.z * ptInt2.z < dSqRad) {

            ptInt1.ToGlob( CircFrame) ;
            ptInt2.ToGlob( CircFrame) ;

            vtN1 = X_AX ;
            vtN2 = - X_AX ;

            vtN1.ToGlob( CircFrame) ;
            vtN2.ToGlob( CircFrame) ;

            return true ;
         }
        // Nessuna intersezione
         else 
            return false ;
      }
     // Nessuna intersezione
      else 
         return false ;
   }

   double dEpsLow = ( bTapLow ? - EPS_SMALL : EPS_SMALL) ;
   double dEpsUp = ( bTapUp ?  EPS_SMALL : - EPS_SMALL) ;

  // L'equazione ammette o due soluzioni (eventualmente
  // coincidenti) oppure nessuna o infinite se la la retta
  // appartiene alla superficie

   Vector3d vtMv( dLongMvLen, dOrtMvLen, 0) ;

   if ( nRoot == 2) {

      ptInt1 = ptP + vdRoots[0] * vtV ;
      ptInt2 = ptP + vdRoots[1] * vtV ;

   
      if ( ptInt1.x > ptInt2.x) 
         
         swap( ptInt1, ptInt2) ;

      Vector3d vtTest1 = ( ptInt1 - ORIG) - ( ptInt1 - ORIG) * vtAx * vtAx ;
      Vector3d vtTest2 = ( ptInt2 - ORIG) - ( ptInt2 - ORIG) * vtAx * vtAx ;

      double dY0_1, dY0_2 ;
      
      if ( vtTest1.y > 0) {

         dY0_1 = ( dSqRad - ptInt1.z * ptInt1.z > 0 ? sqrt( dSqRad - ptInt1.z * ptInt1.z) : 0) ;  
      }
      else {
      
         dY0_1 = ( dSqRad - ptInt1.z * ptInt1.z > 0 ? - sqrt( dSqRad - ptInt1.z * ptInt1.z) : 0) ;  
      }

      Vector3d vtCirc1( 0, - dY0_1, - ptInt1.z) ;
      Vector3d vtTan1( 0, - vtCirc1.z, vtCirc1.y) ;
      Vector3d vtCross1 = vtTan1 ^ vtMv ;

      vtN1 = ( vtCross1 * vtCirc1 > - EPS_ZERO ? vtCross1 : - vtCross1) ;

      if ( vtTest2.y > 0) {

         dY0_2 = ( dSqRad - ptInt2.z * ptInt2.z > 0 ? sqrt( dSqRad - ptInt2.z * ptInt2.z) : 0) ;  
      }
      else {
      
         dY0_2 = ( dSqRad - ptInt2.z * ptInt2.z > 0 ? - sqrt( dSqRad - ptInt2.z * ptInt2.z) : 0) ;  
      }

      Vector3d vtCirc2( 0, - dY0_2, - ptInt2.z) ;
      Vector3d vtTan2( 0, - vtCirc2.z, vtCirc2.y) ;
      Vector3d vtCross2 = vtTan2 ^ vtMv ;

      vtN2 = ( vtCross2 * vtCirc2 > - EPS_ZERO ? vtCross2 : - vtCross2) ;



      if ( ptInt1.x < dLongMvLen + dEpsUp) {
      
         if (  ptInt1.x > + dEpsLow) {

            if ( ptInt2.x > dLongMvLen + dEpsUp) {

               ptInt2 = ptP + ( ( dLongMvLen - ptP.x) / vtV.x) * vtV ;  
               vtN2 = - X_AX ;
            }
         }
         else {

            if ( ptInt2.x > dLongMvLen + dEpsUp) {
            
               ptInt1 = ptP - ( ptP.x / vtV.x) * vtV ;
               ptInt2 = ptP + ( ( dLongMvLen - ptP.x) / vtV.x) * vtV ;

               vtN1.Set( 1, 0, 0) ;
               vtN2.Set( -1, 0, 0) ;

               if ( ptInt1.y * ptInt1.y + ptInt1.z * ptInt1.z > dSqRad &&
                    ptInt2.y * ptInt2.y + ptInt2.z * ptInt2.z > dSqRad)

                    return false ;
            }
            else if ( ptInt2.x > dEpsLow) { 
            
               ptInt1 = ptP - ( ptP.x / vtV.x) * vtV ;
               vtN1.Set( 1, 0, 0) ;
            }
            else 
               return false ;
         }
      }
      else 
         return false ;

     // Riporto le coordinate nel sistema di riferimento
     // griglia
      ptInt1.ToGlob( CircFrame) ;
      ptInt2.ToGlob( CircFrame) ;

      vtN1.ToGlob( CircFrame) ;
      vtN2.ToGlob( CircFrame) ;

      vtN1.Normalize() ;
      vtN2.Normalize() ;
   }   

   return true ;
}

//----------------------------------------------------------------------------
bool 
VolZmap::IntersLineMyPolyhedron( const Point3d& ptLineSt, const Vector3d& vtLineDir,
                                      const Frame3d& PolyFrame, double dLenX, double dLenY, double dLenZ, double dDeltaX,
                                      Point3d& ptInt1, Point3d& ptInt2, Vector3d& vtN1, Vector3d& vtN2) 
{

   double SqIndet = EPS_SMALL * EPS_SMALL ;

   Point3d ptP = ptLineSt ;
   Vector3d vtV = vtLineDir ;

  // Trasformazione delle coordinate
   ptP.ToLoc( PolyFrame) ;
   vtV.ToLoc( PolyFrame) ;

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

  // Servono per descrivere i piani obliqui
   Vector3d vtFacet5 = ptFacet135 - ptP ;
   Vector3d vtFacet6 = ptFacet246 - ptP ;

   Vector3d vtOb( dLenY, - dDeltaX, 0) ;

   vtOb.Normalize() ;

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

   if ( abs( vtV.z) < EPS_ZERO && 
        abs( ptP.z) > dLenZ / 2 - EPS_SMALL) 
   
        return false ;

   int nIntNum = 0 ;

  // Intersezione con la prima faccia
   if ( ptI1.y >= 0 && ptI1.y <= dLenY && 
        ptI1.x * dLenY >= dDeltaX * ptI1.y && ( ptI1.x - dLenX) * dLenY <= dDeltaX * ptI1.y) {

      ptInt1 = ptI1 ;
      vtN1 = - Z_AX ;
      ++ nIntNum ;
   }
  
  // Intersezione con la seconda faccia
   if ( ptI2.y >= 0 && ptI2.y <= dLenY && 
        ptI2.x * dLenY >= dDeltaX * ptI2.y && ( ptI2.x - dLenX) * dLenY <= dDeltaX * ptI2.y) {

      if ( nIntNum == 0) {
      
         ptInt1 = ptI2 ;
         vtN1 = Z_AX ;
         ++ nIntNum ;
      }
      else if ( ( ptInt1 - ptI2).SqLen() > SqIndet) {
      
         ptInt2 = ptI2 ;
         vtN2 = Z_AX ;
         ++ nIntNum ;
      }
   }

  // Intersezione con la terza faccia
   if ( nIntNum < 2 &&
        ptI3.x >= 0 && ptI3.x <= dLenX && 
        ptI3.z >= - ptFacet135.z && ptI3.z <= ptFacet135.z) {

      if ( nIntNum == 0) {
      
         ptInt1 = ptI3 ;
         vtN1 = Y_AX ;
         ++ nIntNum ;
      }
      else if ( ( ptInt1 - ptI3).SqLen() > SqIndet) {
      
         ptInt2 = ptI3 ;
         vtN2 = Y_AX ;
         ++ nIntNum ;
      }
   }

  // Intersezione con la quarta faccia
   if ( nIntNum < 2 &&
        ptI4.x >= dDeltaX && ptI4.x <= dLenX + dDeltaX && 
        ptI4.z >= - ptFacet135.z && ptI4.z <= ptFacet135.z) {

      if ( nIntNum == 0) {
      
         ptInt1 = ptI4 ;
         vtN1 = - Y_AX ;
         ++ nIntNum ;
      }
      else if ( ( ptInt1 - ptI4).SqLen() > SqIndet) {
      
         ptInt2 = ptI4 ;
         vtN2 = - Y_AX ;
         ++ nIntNum ;
      }
   }

  // Intersezione con la quinta faccia
   if ( nIntNum < 2 &&
        ptI5.y >= 0 && ptI5.y <= dLenY && 
        ptI5.z >= - ptFacet135.z && ptI5.z <= ptFacet135.z) {

      if ( nIntNum == 0) {
      
         ptInt1 = ptI5 ;
         vtN1 = vtOb ;
         ++ nIntNum ;
      }
      else if ( ( ptInt1 - ptI5).SqLen() > SqIndet) {
      
         ptInt2 = ptI5 ;
         vtN2 = vtOb ;
         ++ nIntNum ;
      }  
   }

  // Intersezione con la sesta faccia
   if ( nIntNum < 2 &&
        ptI6.y >= 0 && ptI6.y <= dLenY && 
        ptI6.z >= - ptFacet135.z && ptI6.z <= ptFacet135.z) {

      if ( nIntNum == 0) {
      
         ptInt1 = ptI6;
         vtN1 = - vtOb ;
         ++ nIntNum ;
      }
      else if ( ( ptInt1 - ptI6).SqLen() > SqIndet) {
      
         ptInt2 = ptI6;
         vtN2 = - vtOb ;
         ++ nIntNum ;
      } 
   }

   
   if ( nIntNum == 2) {
   
      ptInt1.ToGlob( PolyFrame) ;
      ptInt2.ToGlob( PolyFrame) ;

      vtN1.ToGlob( PolyFrame) ;
      vtN2.ToGlob( PolyFrame) ;

      return true ;
   }
   else 
      return false ;
}