//----------------------------------------------------------------------------
//  EgalTech 2014-2019
//----------------------------------------------------------------------------
// File : EGkTria3d.h           Data : 17.12.19          Versione : 2.1l5
// Contenuto : Dichiarazione classe triangolo Triangle3d.
//
//
//
// Modifiche : 30.03.14 DS Creazione modulo.
//             07.02.15 DS Agg. GetArea e GetAspectRatio.
//             15.01.18 LM Agg. flag.
//             06.06.19 LM Agg. IsPointInsideTriangle.
//             17.12.19 DS Agg. controllo lunghezza di tutti e tre i lati.
//
//----------------------------------------------------------------------------

#pragma once

#include "/EgtDev/Include/EGkBBox3d.h"
#include "/EgtDev/Include/EGkPlane3d.h"
#include <vector>
#include <list>
#include <algorithm>

//-----------------------------------------------------------------------------
class Triangle3d
{
   public :
      Triangle3d( void) : m_nGrade( 0)
              { m_nAttr[0] = 0 ;  m_nAttr[1] = 0 ; m_nAttr[2] = 0 ; }
      void Set( const Point3d& ptP0, const Point3d& ptP1, const Point3d& ptP2)
              { m_ptP[0] = ptP0 ; m_ptP[1] = ptP1 ; m_ptP[2] = ptP2 ; m_vtN = V_NULL ; }
      void Set( const Point3d& ptP0, const Point3d& ptP1, const Point3d& ptP2, const Vector3d& vtV) 
              { m_ptP[0] = ptP0 ; m_ptP[1] = ptP1 ; m_ptP[2] = ptP2 ; m_vtN = vtV ; }
      bool SetP( int nInd, const Point3d& ptP)
              { if ( nInd < 0  ||  nInd >= 3)
                   return false ;
                m_ptP[nInd] = ptP ;
                m_vtN = V_NULL ;
                return true ;
              }
      bool SetGrade( int nFlag) 
              { m_nGrade = nFlag ;
                return true ;
              }
      bool SetAttrib( int nInd, int nVal) 
              { if ( nInd < 0  ||  nInd >= 3)
                   return false ;
                m_nAttr[nInd] = nVal ;
                return true ;
              }
      bool Validate( bool bOverwrite = false)
              { if ( AreSamePointApprox( m_ptP[0], m_ptP[1])  ||
                     AreSamePointApprox( m_ptP[1], m_ptP[2])  ||
                     AreSamePointApprox( m_ptP[2], m_ptP[0]))
                   return false ;
                Vector3d vtV1 = m_ptP[1] - m_ptP[0] ;
                Vector3d vtV2 = m_ptP[2] - m_ptP[1] ;
                Vector3d vtN = vtV1 ^ vtV2 ;
                double dSqN = vtN.SqLen() ;
                if ( dSqN < SQ_EPS_ZERO)
                   return false ;
                if ( dSqN < SQ_EPS_TRIA_H * std::max( { vtV1.SqLen(), vtV2.SqLen(), ( vtV1 + vtV2).SqLen()}))
                   return false ;
                vtN /= sqrt( dSqN) ;
                if ( m_vtN.IsSmall()  ||  bOverwrite) {
                   m_vtN = vtN ;
                   return true ;
                }
                return AreSameVectorApprox( vtN, m_vtN) ;
              }
      bool IsValid( void) const
              { if ( AreSamePointApprox( m_ptP[0], m_ptP[1])  ||
                     AreSamePointApprox( m_ptP[1], m_ptP[2])  ||
                     AreSamePointApprox( m_ptP[2], m_ptP[0]))
                   return false ;
                if ( m_vtN.IsSmall())
                   return false ;
                Vector3d vtV1 = m_ptP[1] - m_ptP[0] ;
                Vector3d vtV2 = m_ptP[2] - m_ptP[1] ;
                Vector3d vtN = vtV1 ^ vtV2 ;
                double dSqN = vtN.SqLen() ;
                if ( dSqN < SQ_EPS_ZERO)
                   return false ;
                if ( dSqN < SQ_EPS_TRIA_H * std::max( { vtV1.SqLen(), vtV2.SqLen(), ( vtV1 + vtV2).SqLen()}))
                   return false ;
                vtN /= sqrt( dSqN) ;
                return AreSameVectorApprox( vtN, m_vtN) ;
              }
      void Translate( const Vector3d& vtMove)
              { m_ptP[0].Translate( vtMove) ;
                m_ptP[1].Translate( vtMove) ;
                m_ptP[2].Translate( vtMove) ;
              }
      bool Rotate( const Point3d& ptAx, const Vector3d& vtAx, double dAngDeg)
              { double dAngRad = dAngDeg * DEGTORAD ;
                return Rotate( ptAx, vtAx, cos( dAngRad), sin( dAngRad)) ;
              }
      bool Rotate( const Point3d& ptAx, const Vector3d& vtAx, double dCosAng, double dSinAng)
              { return ( m_ptP[0].Rotate( ptAx, vtAx, dCosAng, dSinAng)  &&
                         m_ptP[1].Rotate( ptAx, vtAx, dCosAng, dSinAng)  &&
                         m_ptP[2].Rotate( ptAx, vtAx, dCosAng, dSinAng)  &&
                         m_vtN.Rotate( vtAx, dCosAng, dSinAng)) ;
              }
      bool Scale( const Frame3d& frRef, double dCoeffX, double dCoeffY, double dCoeffZ)
              { m_ptP[0].Scale( frRef, dCoeffX, dCoeffY, dCoeffZ) ;
                m_ptP[1].Scale( frRef, dCoeffX, dCoeffY, dCoeffZ) ;
                m_ptP[2].Scale( frRef, dCoeffX, dCoeffY, dCoeffZ) ;
                return Validate( true) ;
              }
      bool Mirror( const Point3d& ptOn, const Vector3d& vtNorm)
              { m_ptP[0].Mirror( ptOn, vtNorm) ;
                m_ptP[1].Mirror( ptOn, vtNorm) ;
                m_ptP[2].Mirror( ptOn, vtNorm) ;
                return Validate( true) ;
              }
      bool Shear( const Point3d& ptOn, const Vector3d& vtNorm, const Vector3d& vtDir, double dCoeff)
              { m_ptP[0].Shear( ptOn, vtNorm, vtDir, dCoeff) ;
                m_ptP[1].Shear( ptOn, vtNorm, vtDir, dCoeff) ;
                m_ptP[2].Shear( ptOn, vtNorm, vtDir, dCoeff) ;
                return Validate( true) ;
              }
      bool ToGlob( const Frame3d& frRef)
              { return ( m_ptP[0].ToGlob( frRef)  &&  m_ptP[1].ToGlob( frRef)  &&  m_ptP[2].ToGlob( frRef)  &&  m_vtN.ToGlob( frRef)) ; }
      bool ToLoc( const Frame3d& frRef)
              { return ( m_ptP[0].ToLoc( frRef)  &&  m_ptP[1].ToLoc( frRef)  &&  m_ptP[2].ToLoc( frRef)  &&  m_vtN.ToLoc( frRef)) ; }
      bool LocToLoc( const Frame3d& frOri, const Frame3d& frDest)
              { return ( m_ptP[0].LocToLoc( frOri, frDest)  &&
                         m_ptP[1].LocToLoc( frOri, frDest)  &&
                         m_ptP[2].LocToLoc( frOri, frDest)  &&
                         m_vtN.LocToLoc( frOri, frDest)) ;
              }
      bool GetLocalBBox( BBox3d& b3Loc) const
              { b3Loc.Reset() ;
                for ( const auto& ptP : m_ptP)
                   b3Loc.Add( ptP) ;
                return true ;
              }
      const Point3d& GetP( int nInd) const
              { if ( nInd >= 0  &&  nInd < 3)
                   return m_ptP[nInd] ;
                else if ( nInd < 0)
                   return m_ptP[0] ;
                else
                   return m_ptP[2] ;
              }
      int GetGrade( void) const 
              { return m_nGrade ; }
      int GetAttrib( int nInd) const
              { if ( nInd >= 0  &&  nInd < 3)
                   return m_nAttr[nInd] ;
                else if ( nInd < 0)
                   return m_nAttr[0] ;
                else
                   return m_nAttr[2] ;
              }
      const Vector3d& GetN( void) const
              { return m_vtN ; }
      Point3d GetCentroid( void) const
              { return ( m_ptP[0] + m_ptP[1] + m_ptP[2]) / 3 ; }
      Plane3d GetPlane( void) const
              { Plane3d plT ; plT.Set( GetCentroid(), m_vtN) ; return plT ; }
      double GetArea( void) const
              { return (( m_ptP[1] - m_ptP[0]) ^ ( m_ptP[2] - m_ptP[0])).Len() / 2  ; }
      double GetAspectRatio( void) const
              { double dSqDistA = SqDist( m_ptP[0], m_ptP[1]) ;
                double dSqDistB = SqDist( m_ptP[1], m_ptP[2]) ;
                double dSqDistC = SqDist( m_ptP[2], m_ptP[0]) ;
                double dTwoArea = (( m_ptP[1] - m_ptP[0]) ^ ( m_ptP[2] - m_ptP[0])).Len() ;
                if ( dTwoArea < SQ_EPS_SMALL)
                   return INFINITO ;
                else
                   return ( std::max( dSqDistA, std::max( dSqDistB, dSqDistC)) / dTwoArea) ;
              }
      double GetSqMinHeight( void) const
              { Vector3d vtV1 = m_ptP[1] - m_ptP[0] ;
                Vector3d vtV2 = m_ptP[2] - m_ptP[1] ;
                Vector3d vtN = vtV1 ^ vtV2 ;
                double dSqN = vtN.SqLen() ;
                if ( dSqN < SQ_EPS_ZERO)
                   return 0 ;
                return ( dSqN / std::max( { vtV1.SqLen(), vtV2.SqLen(), ( vtV1 + vtV2).SqLen()})) ;
              }

   protected :
      Point3d  m_ptP[3] ;
      Vector3d m_vtN ;
      int      m_nGrade ;
      int      m_nAttr[3] ;
} ;

//----------------------------------------------------------------------------
// Raccolte di Triangle3d
typedef std::vector<Triangle3d>  TRIA3DVECTOR ;      // vettore di Triangle3d
typedef std::list<Triangle3d>    TRIA3DLIST ;        // lista di Triangle3d
typedef std::vector<TRIA3DLIST>  TRIA3DLISTVECTOR ;  // vettore di liste di Triangle3d

//-----------------------------------------------------------------------------
// Flags indicanti se i lati sono parte del contorno di un poligono di più triangoli
class TriFlags3d
{
   public :
      TriFlags3d( void)
         { bFlag[0] = true ; bFlag[1] = true ; bFlag[2] = true ; }
      TriFlags3d( bool bF)
         { bFlag[0] = bF ; bFlag[1] = bF ; bFlag[2] = bF ; }
      TriFlags3d( bool bF0, bool bF1, bool bF2)
         { bFlag[0] = bF0 ; bFlag[1] = bF1 ; bFlag[2] = bF2 ; }

   public :
      bool bFlag[3] ;
} ;

//-----------------------------------------------------------------------------
// Normali sui vertici, ottenute mediando opportunamente coi triangoli vicini
class TriNormals3d
{
   public :
      TriNormals3d( void)
         {}
      TriNormals3d( const Vector3d& vtNrm)
         { vtN[0] = vtNrm ; vtN[1] = vtNrm ; vtN[2] = vtNrm ; }
      TriNormals3d( const Vector3d& vtN0, const Vector3d& vtN1, const Vector3d& vtN2)
         { vtN[0] = vtN0 ; vtN[1] = vtN1 ; vtN[2] = vtN2 ; }

   public :
      Vector3d vtN[3] ;
} ;

//-----------------------------------------------------------------------------
// Classe derivata da Triangle3d che include anche i flag dei lati e le normali dei vertici
class Triangle3dEx : public Triangle3d
{
   public :
      Triangle3dEx( void)
              {}
      Triangle3dEx( const Triangle3d& trTria)
              { *(static_cast<Triangle3d*>(this)) = trTria ;
                for ( int i = 0 ; i < 3 ; ++ i)
                   m_triN.vtN[i] = m_vtN ;
              }
      bool SetEdgeFlag( int nInd, bool bVal)
              { if ( nInd < 0  ||  nInd >= 3)
                   return false ;
                m_triF.bFlag[nInd] = bVal ;
                return true ;
              }
      bool GetEdgeFlag( int nInd) const
              { if ( nInd >= 0  &&  nInd < 3)
                   return m_triF.bFlag[nInd] ;
                else if ( nInd < 0)
                   return m_triF.bFlag[0] ;
                else
                   return m_triF.bFlag[2] ;
              }
      const TriFlags3d& GetTriFlags( void) const
              { return m_triF ; }
      void ResetEdgeFlags( void)
              { m_triF.bFlag[0] = m_triF.bFlag[1] = m_triF.bFlag[2] = true ; }
      bool SetVertexNorm( int nInd, const Vector3d& vtN)
              { if ( nInd < 0  ||  nInd >= 3)
                   return false ;
                m_triN.vtN[nInd] = vtN ;
                return true ;
              }
      const Vector3d& GetVertexNorm( int nInd) const
              { if ( nInd >= 0  &&  nInd < 3)
                   return m_triN.vtN[nInd] ;
                else if ( nInd < 0)
                   return m_triN.vtN[0] ;
                else
                   return m_triN.vtN[2] ;
              }
      const TriNormals3d& GetTriNormals( void) const
              { return m_triN ; }
      bool Validate( bool bOverwrite = false)
              { if ( ! Triangle3d::Validate( bOverwrite))
                   return false ;
                for ( int i = 0 ; i < 3 ; ++ i) {
                   if ( m_triN.vtN[i].IsZero())
                      m_triN.vtN[i] = m_vtN ;
                }
                return true ;
              }
      bool Rotate( const Point3d& ptAx, const Vector3d& vtAx, double dAngDeg)
              { double dAngRad = dAngDeg * DEGTORAD ;
                return Rotate( ptAx, vtAx, cos( dAngRad), sin( dAngRad)) ; }
      bool Rotate( const Point3d& ptAx, const Vector3d& vtAx, double dCosAng, double dSinAng)
              { return ( m_ptP[0].Rotate( ptAx, vtAx, dCosAng, dSinAng)  &&
                         m_ptP[1].Rotate( ptAx, vtAx, dCosAng, dSinAng)  &&
                         m_ptP[2].Rotate( ptAx, vtAx, dCosAng, dSinAng)  &&
                         m_vtN.Rotate( vtAx, dCosAng, dSinAng)  &&
                         m_triN.vtN[0].Rotate( vtAx, dCosAng, dSinAng)  &&
                         m_triN.vtN[1].Rotate( vtAx, dCosAng, dSinAng)  &&
                         m_triN.vtN[2].Rotate( vtAx, dCosAng, dSinAng)) ;
              }
      bool Scale( const Frame3d& frRef, double dCoeffX, double dCoeffY, double dCoeffZ)
              { m_ptP[0].Scale( frRef, dCoeffX, dCoeffY, dCoeffZ) ;
                m_ptP[1].Scale( frRef, dCoeffX, dCoeffY, dCoeffZ) ;
                m_ptP[2].Scale( frRef, dCoeffX, dCoeffY, dCoeffZ) ;
                for ( int i = 0 ; i < 3 ; ++ i)
                   m_triN.vtN[i] = V_NULL ;
                return Validate( true) ;
              }
      bool Mirror( const Point3d& ptOn, const Vector3d& vtNorm)
              { m_ptP[0].Mirror( ptOn, vtNorm) ;
                m_ptP[1].Mirror( ptOn, vtNorm) ;
                m_ptP[2].Mirror( ptOn, vtNorm) ;
                for ( int i = 0 ; i < 3 ; ++ i)
                   m_triN.vtN[i].Mirror( vtNorm) ;
                return Validate( true) ;
              }
      bool Shear( const Point3d& ptOn, const Vector3d& vtNorm, const Vector3d& vtDir, double dCoeff)
              { m_ptP[0].Shear( ptOn, vtNorm, vtDir, dCoeff) ;
                m_ptP[1].Shear( ptOn, vtNorm, vtDir, dCoeff) ;
                m_ptP[2].Shear( ptOn, vtNorm, vtDir, dCoeff) ;
                for ( int i = 0 ; i < 3 ; ++ i)
                   m_triN.vtN[i] = V_NULL ;
                return Validate( true) ;
              }
      bool ToGlob( const Frame3d& frRef)
              { return ( m_ptP[0].ToGlob( frRef)  &&  m_ptP[1].ToGlob( frRef)  &&  m_ptP[2].ToGlob( frRef)  &&  m_vtN.ToGlob( frRef)  &&
                         m_triN.vtN[0].ToGlob( frRef)  &&  m_triN.vtN[1].ToGlob( frRef)  &&  m_triN.vtN[2].ToGlob( frRef)) ; }
      bool ToLoc( const Frame3d& frRef)
              { return ( m_ptP[0].ToLoc( frRef)  &&  m_ptP[1].ToLoc( frRef)  &&  m_ptP[2].ToLoc( frRef)  &&  m_vtN.ToLoc( frRef)  &&
                         m_triN.vtN[0].ToLoc( frRef)  &&  m_triN.vtN[1].ToLoc( frRef)  &&  m_triN.vtN[2].ToLoc( frRef)) ; }
      bool LocToLoc( const Frame3d& frOri, const Frame3d& frDest)
              { return ( m_ptP[0].LocToLoc( frOri, frDest)  &&
                         m_ptP[1].LocToLoc( frOri, frDest)  &&
                         m_ptP[2].LocToLoc( frOri, frDest)  &&
                         m_vtN.LocToLoc( frOri, frDest)  &&
                         m_triN.vtN[0].LocToLoc( frOri, frDest)  &&
                         m_triN.vtN[1].LocToLoc( frOri, frDest)  &&
                         m_triN.vtN[2].LocToLoc( frOri, frDest)) ; }

   private :
     TriFlags3d   m_triF ;
     TriNormals3d m_triN ;
} ;

//----------------------------------------------------------------------------
// Raccolte di Triangle3dEx
typedef std::vector<Triangle3dEx>  TRIA3DEXVECTOR ;      // vettore di Triangle3dEx
typedef std::list<Triangle3dEx>    TRIA3DEXLIST ;        // lista di Triangle3dEx
typedef std::vector<TRIA3DEXLIST>  TRIA3DEXLISTVECTOR ;  // vettore di liste di Triangle3dEx

//-----------------------------------------------------------------------------
enum PlaneType { PL_NULL = 0,
                 PL_XY = 1,
                 PL_YZ = 2,
                 PL_ZX = 3} ;

//----------------------------------------------------------------------------
// Piano canonico di miglior proiezione (perpendicolare alla componente maggiore della normale)
inline bool
CalcProjPlane( const Vector3d& vtN, int& nPlane, bool& bCCW)
{
  // verifico che la normale non sia nulla
   if ( vtN.IsZero())
      return false ;
  // proiezione sul piano XY (Nz con valore maggiore)
   if ( abs( vtN.z) > abs( vtN.x)  &&
        abs( vtN.z) > abs( vtN.y)) {
      nPlane = PL_XY ;
      bCCW = ( vtN.z > 0) ;
   }
  // proiezione sul piano YZ (Nx con valore maggiore)
   else if ( abs( vtN.x) > abs( vtN.y)) {
      nPlane = PL_YZ ;
      bCCW = ( vtN.x > 0) ;
   }
  // proiezione sul piano ZX (Ny con valore maggiore)
   else {
      nPlane = PL_ZX ;
      bCCW = ( vtN.y > 0) ;
   }
   return true ;
}

//----------------------------------------------------------------------------
// Prodotto vettoriale nel piano di proiezione ( positivo se i 3 punti in ordine CCW)
inline double
TwoAreaInPlane( int nPlane, const Point3d& ptA, const Point3d& ptB, const Point3d& ptC)
{
   switch ( nPlane) {
   default : // PL_XY
      return ( ptB.x - ptA.x) * ( ptC.y - ptB.y) - ( ptB.y - ptA.y) * ( ptC.x - ptB.x) ;
   case PL_YZ :
      return ( ptB.y - ptA.y) * ( ptC.z - ptB.z) - ( ptB.z - ptA.z) * ( ptC.y - ptB.y) ;
   case PL_ZX :
      return ( ptB.z - ptA.z) * ( ptC.x - ptB.x) - ( ptB.x - ptA.x) * ( ptC.z - ptB.z) ;
   }
}

//----------------------------------------------------------------------------
// Prodotto scalare nel piano di proiezione
inline double
ProScaInPlane( int nPlane, const Point3d& ptP, const Point3d& ptA, const Point3d& ptB)
{
   switch ( nPlane) {
   default : // PL_XY
      return ( ptP.x - ptA.x) * ( ptB.x - ptA.x) + ( ptP.y - ptA.y) * ( ptB.y - ptA.y) ;
   case PL_YZ :
      return ( ptP.y - ptA.y) * ( ptB.y - ptA.y) + ( ptP.z - ptA.z) * ( ptB.z - ptA.z) ;
   case PL_ZX :
      return ( ptP.z - ptA.z) * ( ptB.z - ptA.z) + ( ptP.x - ptA.x) * ( ptB.x - ptA.x) ;
   }
}

//----------------------------------------------------------------------------
// Coordinate baricentriche di un punto giacente nel piano del triangolo
inline bool
BarycentricCoord( const Point3d& ptP, const Triangle3d& Tria,
                  double& dU, double& dV, double& dW)
{
  // calcolo del piano ottimale di proiezione
   int nPlane ;
   bool bCCW ;
   if ( ! CalcProjPlane( Tria.GetN(), nPlane, bCCW))
      return false ;
  // verifico che l'area (doppia) non sia nulla
   double d2Area = TwoAreaInPlane( nPlane, Tria.GetP( 0), Tria.GetP( 1), Tria.GetP( 2)) ;
   if ( abs( d2Area) < SQ_EPS_SMALL)
      return false ;
  // calcolo delle coordinate baricentriche
   double dInv2Area = 1 / d2Area ;
   dU = TwoAreaInPlane( nPlane, ptP, Tria.GetP( 1), Tria.GetP( 2)) * dInv2Area ;
   dV = TwoAreaInPlane( nPlane, Tria.GetP( 0), ptP, Tria.GetP( 2)) * dInv2Area ;
   dW = TwoAreaInPlane( nPlane, Tria.GetP( 0), Tria.GetP( 1), ptP) * dInv2Area ;
  // devono dare somma unitaria
   double dSumm = dU + dV + dW ;
   if ( abs( dSumm) < EPS_ZERO)
      return false ;
   if ( abs( dSumm - 1) > EPS_ZERO) {
      double dInvSumm = 1 / dSumm ; 
      dU *= dInvSumm ;
      dV *= dInvSumm ;
      dW *= dInvSumm ;
   }
   return true ;
}

//----------------------------------------------------------------------------
// Normale in un punto del triangolo, come media baricentrica delle normali nei vertici
inline bool
CalcNormal( const Point3d& ptP, const Triangle3d& Tria, const TriNormals3d& Tnorms, Vector3d& vtNorm)
{
  // calcolo le coordinate baricentriche del punto
   double dU, dV, dW ;
   if ( ! BarycentricCoord( ptP, Tria, dU, dV, dW))
      return false ;
  // calcolo la media
   vtNorm = dU * Tnorms.vtN[0] + dV * Tnorms.vtN[1] + dW * Tnorms.vtN[2] ;
  // la normalizzo
   return vtNorm.Normalize( EPS_ZERO) ;
}

//----------------------------------------------------------------------------
enum TriangleType { OPEN = -1, EXACT = 0, CLOSED = 1 };

//----------------------------------------------------------------------------
// Stabilisce se un punto appartiene al triangolo aperto, esatto o chiuso.
inline bool
IsPointInsideTriangle( const Point3d& ptP, const Triangle3d& trTria, int nTriType)
{
  // Se il punto non è nel box del triangolo, è sicuramente esterno
   BBox3d b3Tria( trTria.GetP( 0)) ;
   b3Tria.Add( trTria.GetP( 1)) ;
   b3Tria.Add( trTria.GetP( 2)) ;
   if ( ! b3Tria.Encloses( ptP))
      return false ;

  // Se il punto non sta sul piano del triangolo, è sicuramente esterno
   if ( abs( ( ptP - trTria.GetP(0)) * trTria.GetN()) > EPS_SMALL)
      return false ;

  // Tolleranza di appartenenza dipendente dal tipo di triangolo
   double dEps ;
   switch ( nTriType) {
   case OPEN :
      dEps = EPS_SMALL ;
      break ;
   case EXACT :
      dEps = 0 ;
      break ;
   case CLOSED :
      dEps = -EPS_SMALL ;
      break ;
   default :
      return false ;
   }

  // Verifico se il punto è a destra del primo lato
   Vector3d vtV0 = trTria.GetP( 1) - trTria.GetP( 0) ;
   vtV0.Normalize() ;
   double dProd0 = ( vtV0 ^ (ptP - trTria.GetP( 0))) * trTria.GetN() ;
   if ( dProd0 < dEps)
      return false ;
  // Verifico se il punto è a destra del secondo lato
   Vector3d vtV1 = trTria.GetP( 2) - trTria.GetP( 1) ;
   vtV1.Normalize() ;
   double dProd1 = ( vtV1 ^ (ptP - trTria.GetP( 1))) * trTria.GetN() ;
   if ( dProd1 < dEps)
      return false ;
  // Verifico se il punto è a destra del terzo lato
   Vector3d vtV2 = trTria.GetP( 0) - trTria.GetP( 2) ;
   vtV2.Normalize() ;
   double dProd2 = ( vtV2 ^ (ptP - trTria.GetP(2))) * trTria.GetN() ;
   if ( dProd2 < dEps)
      return false ;
  // Punto a sinistra di tutti i lati, quindi interno
   return true ;
}

//----------------------------------------------------------------------------
// Stabilisce se un punto appartiene a un vertice del triangolo
inline bool
IsPointOnTriangleVertex( const Point3d& ptP, const Triangle3d& trTria)
{
   return ( AreSamePointApprox( trTria.GetP( 0), ptP)  ||
            AreSamePointApprox( trTria.GetP( 1), ptP)  ||
            AreSamePointApprox( trTria.GetP( 2), ptP)) ;
}

//----------------------------------------------------------------------------
inline bool
GetTriaVertexNearestToLine( const Triangle3d& trTria, const Point3d& ptL, const Vector3d& vtL,
                            int& nVert, double& dSqDist)
{
  // Imposto valori di default
   nVert = - 1 ;
   dSqDist = SQ_INFINITO ;
  // Verifico validità triangolo
   if ( ! trTria.IsValid())
      return false ;
  // Verifico validità linea
   Vector3d vtDir = vtL ;
   if ( ! vtDir.Normalize())
      return false ;
  // Cerco il vertice più vicino
   for ( int i = 0 ; i < 3 ; ++ i) {
      Point3d ptP = trTria.GetP( i) ;
      double dCurrSqDist = (( ptP - ptL) - vtL * ( vtL * ( ptP - ptL))).SqLen() ;
      if ( dCurrSqDist < dSqDist) {
         dSqDist = dCurrSqDist ;
         nVert = i ;
      }
   }
   return ( nVert != -1) ;
}

//----------------------------------------------------------------------------
inline bool
GetTriaMidEdgeNearestToLine( const Triangle3d& trTria, const Point3d& ptL, const Vector3d& vtL,
                             Point3d& ptMid, double& dSqDist)
{
  // Imposto valori di default
   dSqDist = SQ_INFINITO ;
  // Verifico validità triangolo
   if ( ! trTria.IsValid())
      return false ;
  // Verifico validità linea
   Vector3d vtDir = vtL ;
   if ( ! vtDir.Normalize())
      return false ;
  // Cerco il punto medio dei lati del triangolo più vicino
   bool bFound = false ;
   for ( int i = 0 ; i < 3 ; ++ i) {
      Point3d ptP = Media( trTria.GetP( i), trTria.GetP( (i + 1) % 3)) ;
      double dCurrSqDist = (( ptP - ptL) - vtL * ( vtL * ( ptP - ptL))).SqLen() ;
      if ( dCurrSqDist < dSqDist) {
         dSqDist = dCurrSqDist ;
         ptMid = ptP ;
         bFound = true ;
      }
   }
   return bFound ;
}