EgtGeomKernel/IntersLineCyl.cpp

//----------------------------------------------------------------------------
//  EgalTech 2023-2023
//----------------------------------------------------------------------------
// File : IntersLineCyl.cpp        Data : 16.05.23          Versione : 2.5e3
// Contenuto : Implementazione della intersezione linea/cilindro.
//
//
//
// Modifiche : 16.05.23 DS Creazione modulo.
//
//
//----------------------------------------------------------------------------

//--------------------------- Include ----------------------------------------
#include "stdafx.h"
#include "IntersLineCyl.h"
#include "/EgtDev/Include/EGkFrame3d.h"
#include "/EgtDev/Include/ENkPolynomialRoots.h"

using namespace std ;

//----------------------------------------------------------------------------
// Linea e cilindro sono nel medesimo riferimento.
// Il cilindro è centrato sull'asse Z e appoggiato sul piano XY.
// In caso di intersezione viene restituito true e i parametri in dU1 e dU2.
//----------------------------------------------------------------------------
bool
IntersLineCyl( const Point3d& ptL, const Vector3d& vtL,
               double dRad, double dHeight,
               double& dU1, double& dU2)
{
  // Verifico il versore
   if ( vtL.IsSmall())
      return false ;

  // Verifico il cilindro
   if ( dRad < EPS_SMALL  ||  dHeight < EPS_SMALL)
      return false ;

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

  // Se la linea interseca entrambe le basi, si sono trovate le due intersezioni
   if ( nBasInt == 3) {
      if ( dU1 > dU2)
         swap( dU1, dU2) ;
     // Trovate intersezioni
      return true ;
   }
   
  // Determino le intersezioni con la superficie laterale del cilindro
   DBLVECTOR vdCoeff{ ptL.x * ptL.x + ptL.y * ptL.y - dRad * dRad,
                      2 * ( ptL.x * vtL.x + ptL.y * vtL.y),
                      vtL.x * vtL.x + vtL.y * vtL.y} ;
   DBLVECTOR vdRoots ;
   int nRoot = PolynomialRoots( 2, vdCoeff, vdRoots) ;

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

  // Due soluzioni: la retta interseca due volte la superficie laterale
   if ( nRoot == 2) {
      dU1 = vdRoots[0] ;
      dU2 = vdRoots[1] ;
      if ( dU1 > dU2)
         swap( dU1, dU2) ;
     // Trovate intersezioni
      return true ;
   }
   
  // Una soluzione : la retta interseca la superficie laterale e un piano
   else if ( nRoot == 1) {
     // Se piano superiore
      if ( nBasInt == 2) {
         dU1 = vdRoots[0] ;
      }
     // altrimenti piano inferiore
      else if ( nBasInt == 1) {
         dU2 = vdRoots[0] ;
      }
     // altrimenti niente
      else
         return false ;
      if ( dU1 > dU2)
         swap( dU1, dU2) ;
     // Trovate intersezioni
      return true ;
   }

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

//----------------------------------------------------------------------------
// Linea e cilindro sono nel medesimo riferimento.
// Il cilindro è definito con centro della base, asse, raggio e altezza.
// In caso di intersezione viene restituito true e i parametri in dU1 e dU2.
//----------------------------------------------------------------------------
bool
IntersLineCyl( const Point3d& ptL, const Vector3d& vtL,
               const Point3d& ptCyl, const Vector3d& vtCyl, double dRad, double dHeight,
               double& dU1, double& dU2) 
{
  // Riferimento intrinseco del cilindro
   Frame3d frCyl ;
   if ( ! frCyl.Set( ptCyl, vtCyl))
      return false ;
  // Ora eseguo i conti nel riferimento intrinseco
   return IntersLineCyl( GetToLoc( ptL, frCyl), GetToLoc( vtL, frCyl), dRad, dHeight, dU1, dU2) ;
}

//----------------------------------------------------------------------------
// Linea e cilindro sono nel medesimo riferimento.
// Il cilindro è definito con centri delle due basi, e raggio.
// In caso di intersezione viene restituito true e i parametri in dU1 e dU2.
//----------------------------------------------------------------------------
bool
IntersLineCyl( const Point3d& ptL, const Vector3d& vtL,
               const Point3d& ptCyl1, const Point3d& ptCyl2, double dRad,
               double& dU1, double& dU2) 
{
  // Determino asse ed altezza del cilindro
   Vector3d vtCyl = ptCyl2 - ptCyl1 ;
   double dHeight = vtCyl.Len() ;
   if ( dHeight < EPS_SMALL)
      return false ;
   vtCyl /= dHeight ;
  // Riferimento intrinseco del cilindro
   Frame3d frCyl ;
   if ( ! frCyl.Set( ptCyl1, vtCyl))
      return false ;
  // Ora eseguo i conti nel riferimento intrinseco
   return IntersLineCyl( GetToLoc( ptL, frCyl), GetToLoc( vtL, frCyl), dRad, dHeight, dU1, dU2) ;
}