EgtGeomKernel/ProjectCurveSurf.cpp

//----------------------------------------------------------------------------
//  EgalTech 2023-2023
//----------------------------------------------------------------------------
// File : ProjectCurveSurfTm.cpp    Data : 16.11.23        Versione : 2.5kh3
// Contenuto : Implementazione funzioni proiezione curve su superficie Trimesh.
//
//
//
// Modifiche : 31.08.23 DS Creazione modulo.
//
//
//----------------------------------------------------------------------------

//--------------------------- Include ----------------------------------------
#include "stdafx.h"
#include "SurfTriMesh.h"
#include "SurfBezier.h"
#include "GeoConst.h"
#include "/EgtDev/Include/EGkDistPointLine.h"
#include "/EgtDev/Include/EGkDistPointCurve.h"
#include "/EgtDev/Include/EGkDistPointSurfTm.h"
#include "/EgtDev/Include/EGkIntersLineSurfTm.h"
#include "/EgtDev/Include/EGkProjectCurveSurf.h"

using namespace std ;

//----------------------------------------------------------------------------
// Angolo limite tra normale al triangolo e direzione di proiezione 89°
const double COS_ANG_LIM = 0.0175 ;

//----------------------------------------------------------------------------
bool
ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const Vector3d& vtDir,
                    double dLinTol, double dMaxSegmLen, PNT5AXVECTOR& vPt5ax)
{
  // sistemazioni per tipo di superficie
   const SurfTriMesh* pSurfTm = nullptr ;
   switch ( sfSurf.GetType()) {
   case SRF_TRIMESH :
      pSurfTm = GetBasicSurfTriMesh( &sfSurf) ;
      break ;
   case SRF_BEZIER :
      pSurfTm = GetBasicSurfBezier( &sfSurf)->GetAuxSurf() ;
      break ;
   case SRF_FLATRGN :
      pSurfTm = GetBasicSurfFlatRegion( &sfSurf)->GetAuxSurf() ;
      break ;
   }
   if ( pSurfTm == nullptr)
      return false ;

  // controllo le tolleranze
   dLinTol = max( dLinTol, LIN_TOL_MIN) ;
   dMaxSegmLen = max( dMaxSegmLen, 10 * EPS_SMALL) ;
  // approssimo la curva con una polilinea alla massima risoluzione
   PolyLine PL ;
   if ( ! crCrv.ApproxWithLines( EPS_SMALL, ANG_TOL_STD_DEG, ICurve::APL_STD, PL))
      return false ;
   const double MAX_SEG_LEN = min( dMaxSegmLen, 0.977) ;
   if ( ! PL.AdjustForMaxSegmentLen( MAX_SEG_LEN))
      return false ;

  // Oggetto per calcolo massivo intersezioni tra linee di proiezione e superficie
   Frame3d frRefLine ;
   if ( ! frRefLine.Set( ORIG, vtDir))
      return false ;
   IntersParLinesSurfTm intPLSTM( frRefLine, *pSurfTm) ;

  // Vettore locale dei punti risultanti
   PNT5AXVECTOR vMyPt5ax ;
   vMyPt5ax.reserve( PL.GetPointNbr()) ;

  // proietto i punti della polilinea sulla superficie
   double dPar ;
   Point3d ptP ;
   bool bFound = PL.GetFirstUPoint( &dPar, &ptP) ;
   while ( bFound) {
     // intersezione retta di proiezione con superficie
      Point3d ptL = GetToLoc( ptP, frRefLine) ;
      ILSIVECTOR vIntRes ;
      intPLSTM.GetInters( ptL, 1, vIntRes, false) ;
     // cerco la prima intersezione valida a partire dall'ultima (è la più alta)
      int nI = int( vIntRes.size()) - 1 ;
      while ( nI >= 0  &&  abs( vIntRes[nI].dCosDN) < COS_ANG_LIM)
         --nI ;
     // se trovata
      if ( nI >= 0) {
        // calcolo il punto
         Point3d ptInt ;
         if ( vIntRes[nI].nILTT == ILTT_SEGM  ||  vIntRes[nI].nILTT == ILTT_SEGM_ON_EDGE)
            ptInt = vIntRes[nI].ptI2 ;
         else
            ptInt = vIntRes[nI].ptI ;
        // calcolo la normale (si calcola smooth, in caso di errore si prende quella del triangolo)
         Triangle3dEx trTria ;
         if ( ! pSurfTm->GetTriangle( vIntRes[nI].nT, trTria))
            return false ;
         Vector3d vtN ;
         if ( ! CalcNormal( ptInt, trTria, vtN))
            vtN = trTria.GetN() ;
        // aggiungo al vettore dei proiettati
         vMyPt5ax.emplace_back( ptInt, vtN, vtDir, dPar, 1) ;
      }
      bFound = PL.GetNextUPoint( &dPar, &ptP) ;
   }

   // eventuale rimozione punti in eccesso rispetto alle tolleranze lasciata alla funzione chiamante

  // copio i punti nel vettore di ritorno
   vPt5ax.clear() ;
   for ( const auto& Pt5ax : vMyPt5ax) {
      if ( Pt5ax.nFlag != -1)
         vPt5ax.emplace_back( Pt5ax) ;
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const IGeoPoint3d& gpRef,
                    double dLinTol, double dMaxSegmLen, PNT5AXVECTOR& vPt5ax)
{
  // sistemazioni per tipo di superficie
   const SurfTriMesh* pSurfTm = nullptr ;
   switch ( sfSurf.GetType()) {
   case SRF_TRIMESH :
      pSurfTm = GetBasicSurfTriMesh( &sfSurf) ;
      break ;
   case SRF_BEZIER :
      pSurfTm = GetBasicSurfBezier( &sfSurf)->GetAuxSurf() ;
      break ;
   case SRF_FLATRGN :
      pSurfTm = GetBasicSurfFlatRegion( &sfSurf)->GetAuxSurf() ;
      break ;
   }
   if ( pSurfTm == nullptr)
      return false ;

  // controllo le tolleranze
   dLinTol = max( dLinTol, LIN_TOL_MIN) ;
   dMaxSegmLen = max( dMaxSegmLen, 10 * EPS_SMALL) ;

  // approssimo la curva con una polilinea entro la metà della tolleranza
   PolyLine PL ;
   if ( ! crCrv.ApproxWithLines( dLinTol, ANG_TOL_STD_DEG, ICurve::APL_STD, PL))
      return false ;
   const double MAX_SEG_LEN = min( dMaxSegmLen, 0.977) ;
   if ( ! PL.AdjustForMaxSegmentLen( MAX_SEG_LEN))
      return false ;

  // Vettore locale dei punti risultanti
   PNT5AXVECTOR vMyPt5ax ;
   vMyPt5ax.reserve( PL.GetPointNbr()) ;

  // proietto i punti della polilinea sulla superficie con direzione data dal punto di riferimento
   double dPar ;
   Point3d ptP ;
   bool bFound = PL.GetFirstUPoint( &dPar, &ptP) ;
   while ( bFound) {
     // punto di riferimento
      Point3d ptMin = gpRef.GetPoint() ;
     // intersezione della retta di minima distanza con la superficie
      Vector3d vtLine = ptP - ptMin ;
      double dLineLen = vtLine.Len() ;
      if ( dLineLen > EPS_SMALL) {        
         vtLine /= dLineLen ;
         ILSIVECTOR vIntRes ;
         if ( IntersLineSurfTm( ptP, vtLine, dLineLen, *pSurfTm, vIntRes, false)) {
           // cerco la prima intersezione valida a partire dall'ultima (è la più alta)
            int nI = int( vIntRes.size()) - 1 ;
            while ( nI >= 0  &&  abs( vIntRes[nI].dCosDN) < COS_ANG_LIM)
               --nI ;
           // se trovata
            if ( nI >= 0) {
              // calcolo il punto
               Point3d ptInt ;
               if ( vIntRes[nI].nILTT == ILTT_SEGM  ||  vIntRes[nI].nILTT == ILTT_SEGM_ON_EDGE)
                  ptInt = vIntRes[nI].ptI2 ;
               else
                  ptInt = vIntRes[nI].ptI ;
              // calcolo la normale (si calcola smooth, in caso di errore si prende quella del triangolo)
               Triangle3dEx trTria ;
               if ( ! pSurfTm->GetTriangle( vIntRes[nI].nT, trTria))
                  return false ;
               Vector3d vtN ;
               if ( ! CalcNormal( ptInt, trTria, vtN))
                  vtN = trTria.GetN() ;
              // aggiungo al vettore dei proiettati
               vMyPt5ax.emplace_back( ptInt, vtN, vtLine, dPar, 1) ;
            }
         }
      }
      bFound = PL.GetNextUPoint( &dPar, &ptP) ;
   }

  // eventuale rimozione punti in eccesso rispetto alle tolleranze lasciata alla funzione chiamante

  // copio i punti nel vettore di ritorno
   vPt5ax.clear() ;
   for ( const auto& Pt5ax : vMyPt5ax) {
      if ( Pt5ax.nFlag != -1)
         vPt5ax.emplace_back( Pt5ax) ;
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const ICurve& crRef,
                    double dLinTol, double dMaxSegmLen, PNT5AXVECTOR& vPt5ax)
{
  // sistemazioni per tipo di superficie
   const SurfTriMesh* pSurfTm = nullptr ;
   switch ( sfSurf.GetType()) {
   case SRF_TRIMESH :
      pSurfTm = GetBasicSurfTriMesh( &sfSurf) ;
      break ;
   case SRF_BEZIER :
      pSurfTm = GetBasicSurfBezier( &sfSurf)->GetAuxSurf() ;
      break ;
   case SRF_FLATRGN :
      pSurfTm = GetBasicSurfFlatRegion( &sfSurf)->GetAuxSurf() ;
      break ;
   }
   if ( pSurfTm == nullptr)
      return false ;

  // controllo le tolleranze
   dLinTol = max( dLinTol, LIN_TOL_MIN) ;
   dMaxSegmLen = max( dMaxSegmLen, 10 * EPS_SMALL) ;

  // approssimo la curva con una polilinea alla massima risoluzione
   PolyLine PL ;
   if ( ! crCrv.ApproxWithLines( EPS_SMALL, ANG_TOL_STD_DEG, ICurve::APL_STD, PL))
      return false ;
   const double MAX_SEG_LEN = min( dMaxSegmLen, 0.977) ;
   if ( ! PL.AdjustForMaxSegmentLen( MAX_SEG_LEN))
      return false ;

  // Vettore locale dei punti risultanti
   PNT5AXVECTOR vMyPt5ax ;
   vMyPt5ax.reserve( PL.GetPointNbr()) ;

  // proietto i punti della polilinea sulla superficie con direzione normale alla curva di riferimento
   double dPar ;
   Point3d ptP ;
   bool bFound = PL.GetFirstUPoint( &dPar, &ptP) ;
   while ( bFound) {
     // punto sulla curva a minima distanza
      DistPointCurve dPC( ptP, crRef) ;
      Point3d ptMin ;
      int nFlag ;
      if ( dPC.GetMinDistPoint( 0, ptMin, nFlag)) {
        // intersezione della retta di minima distanza con la superficie
         Vector3d vtLine = ptP - ptMin ;
         double dLineLen = vtLine.Len() ;
         if ( dLineLen > EPS_SMALL) {        
            vtLine /= dLineLen ;
            ILSIVECTOR vIntRes ;
            if ( IntersLineSurfTm( ptP, vtLine, dLineLen, *pSurfTm, vIntRes, false)) {
              // cerco la prima intersezione valida a partire dall'ultima (è la più alta)
               int nI = int( vIntRes.size()) - 1 ;
               while ( nI >= 0  &&  abs( vIntRes[nI].dCosDN) < COS_ANG_LIM)
                  --nI ;
              // se trovata
               if ( nI >= 0) {
                 // calcolo il punto
                  Point3d ptInt ;
                  if ( vIntRes[nI].nILTT == ILTT_SEGM  ||  vIntRes[nI].nILTT == ILTT_SEGM_ON_EDGE)
                     ptInt = vIntRes[nI].ptI2 ;
                  else
                     ptInt = vIntRes[nI].ptI ;
                 // calcolo la normale (si calcola smooth, in caso di errore si prende quella del triangolo)
                  Triangle3dEx trTria ;
                  if ( ! pSurfTm->GetTriangle( vIntRes[nI].nT, trTria))
                     return false ;
                  Vector3d vtN ;
                  if ( ! CalcNormal( ptInt, trTria, vtN))
                     vtN = trTria.GetN() ;
                 // aggiungo al vettore dei proiettati
                  vMyPt5ax.emplace_back( ptInt, vtN, vtLine, dPar, 1) ;
               }
            }
         }
      }
      bFound = PL.GetNextUPoint( &dPar, &ptP) ;
   }

  // eventuale rimozione punti in eccesso rispetto alle tolleranze lasciata alla funzione chiamante

  // copio i punti nel vettore di ritorno
   vPt5ax.clear() ;
   for ( const auto& Pt5ax : vMyPt5ax) {
      if ( Pt5ax.nFlag != -1)
         vPt5ax.emplace_back( Pt5ax) ;
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
ProjectCurveOnSurf( const ICurve& crCrv, const ISurf& sfSurf, const ISurf& sfRef,
                    double dLinTol, double dMaxSegmLen, PNT5AXVECTOR& vPt5ax)
{
  // sistemazioni per tipo di superficie
   const SurfTriMesh* pSurfTm = nullptr ;
   switch ( sfSurf.GetType()) {
   case SRF_TRIMESH :
      pSurfTm = GetBasicSurfTriMesh( &sfSurf) ;
      break ;
   case SRF_BEZIER :
      pSurfTm = GetBasicSurfBezier( &sfSurf)->GetAuxSurf() ;
      break ;
   case SRF_FLATRGN :
      pSurfTm = GetBasicSurfFlatRegion( &sfSurf)->GetAuxSurf() ;
      break ;
   }
   if ( pSurfTm == nullptr)
      return false ;

  // sistemazioni per tipo di superficie di riferimento
   const SurfTriMesh* pRefTm = nullptr ;
   switch ( sfRef.GetType()) {
   case SRF_TRIMESH :
      pRefTm = GetBasicSurfTriMesh( &sfRef) ;
      break ;
   case SRF_BEZIER :
      pRefTm = GetBasicSurfBezier( &sfRef)->GetAuxSurf() ;
      break ;
   case SRF_FLATRGN :
      pRefTm = GetBasicSurfFlatRegion( &sfRef)->GetAuxSurf() ;
      break ;
   }
   if ( pRefTm == nullptr)
      return false ;

  // controllo le tolleranze
   dLinTol = max( dLinTol, LIN_TOL_MIN) ;
   dMaxSegmLen = max( dMaxSegmLen, 10 * EPS_SMALL) ;

  // approssimo la curva con una polilinea entro la metà della tolleranza
   PolyLine PL ;
   if ( ! crCrv.ApproxWithLines( dLinTol, ANG_TOL_STD_DEG, ICurve::APL_STD, PL))
      return false ;
   const double MAX_SEG_LEN = min( dMaxSegmLen, 0.977) ;
   if ( ! PL.AdjustForMaxSegmentLen( MAX_SEG_LEN))
      return false ;

  // Vettore locale dei punti risultanti
   PNT5AXVECTOR vMyPt5ax ;
   vMyPt5ax.reserve( PL.GetPointNbr()) ;

  // proietto i punti della polilinea sulla superficie con direzione normale alla curva di riferimento
   double dPar ;
   Point3d ptP ;
   bool bFound = PL.GetFirstUPoint( &dPar, &ptP) ;
   while ( bFound) {
     // punto sulla superficie guida a minima distanza
      DistPointSurfTm dPS( ptP, *pRefTm) ;
      Point3d ptMin ;
      int nTriaMin ;
      if ( dPS.GetMinDistPoint( ptMin)  &&  dPS.GetMinDistTriaIndex ( nTriaMin)) {
        // recupero direzione della retta di minima distanza, altrimenti normale alla superficie
         Vector3d vtLine = ptP - ptMin ;
         double dLineLen = vtLine.Len() ;
         if ( dLineLen > EPS_SMALL)       
            vtLine /= dLineLen ;
         else {
           // calcolo la normale della superficie guida
            Triangle3dEx trGuide ;
            if ( ! pRefTm->GetTriangle( nTriaMin, trGuide))
               return false ;
            if ( ! CalcNormal( ptMin, trGuide, vtLine))
               vtLine = trGuide.GetN() ;
            dLineLen = 100 ;
         }
        // intersezione della retta con la superficie
         ILSIVECTOR vIntRes ;
         if ( IntersLineSurfTm( ptP, vtLine, dLineLen, *pSurfTm, vIntRes, false)) {
           // cerco la prima intersezione valida a partire dall'ultima (è la più alta)
            int nI = int( vIntRes.size()) - 1 ;
            while ( nI >= 0  &&  abs( vIntRes[nI].dCosDN) < COS_ANG_LIM)
               --nI ;
           // se trovata
            if ( nI >= 0) {
              // calcolo il punto
               Point3d ptInt ;
               if ( vIntRes[nI].nILTT == ILTT_SEGM  ||  vIntRes[nI].nILTT == ILTT_SEGM_ON_EDGE)
                  ptInt = vIntRes[nI].ptI2 ;
               else
                  ptInt = vIntRes[nI].ptI ;
              // calcolo la normale (si calcola smooth, in caso di errore si prende quella del triangolo)
               Triangle3dEx trTria ;
               if ( ! pSurfTm->GetTriangle( vIntRes[nI].nT, trTria))
                  return false ;
               Vector3d vtN ;
               if ( ! CalcNormal( ptMin, trTria, vtN))
                  vtN = trTria.GetN() ;
              // calcolo la normale della superficie guida
               Triangle3dEx trGuide ;
               if ( ! pRefTm->GetTriangle( nTriaMin, trGuide))
                  return false ;
               Vector3d vtN2 ;
               if ( ! CalcNormal( ptMin, trGuide, vtN2))
                  vtN2 = trGuide.GetN() ;
              // aggiungo al vettore dei proiettati
               vMyPt5ax.emplace_back( ptInt, vtN, vtN2, dPar, 1) ;
            }
         }
      }
      bFound = PL.GetNextUPoint( &dPar, &ptP) ;
   }

  // eventuale rimozione punti in eccesso rispetto alle tolleranze lasciata alla funzione chiamante

  // copio i punti nel vettore di ritorno
   vPt5ax.clear() ;
   for ( const auto& Pt5ax : vMyPt5ax) {
      if ( Pt5ax.nFlag != -1)
         vPt5ax.emplace_back( Pt5ax) ;
   }

   return true ;
}

//----------------------------------------------------------------------------
static bool
PointsInTolerance( const PNT5AXVECTOR& vPt5ax, int nPrec, int nCurr, int nNext, double dSqTol)
{
   for ( int i = nPrec + 1 ; i < nCurr ; ++ i) {
      double dSqDist ;
      if ( ! DistPointLine( vPt5ax[i].ptP, vPt5ax[nPrec].ptP, vPt5ax[nNext].ptP).GetSqDist( dSqDist)  ||  dSqDist > dSqTol)
         return false ;
   }
   return true ;
}

//----------------------------------------------------------------------------
bool
RemovePointsInExcess( PNT5AXVECTOR& vMyPt5ax, double dLinTol, double dMaxSegmLen, bool bTestDir)
{
  // rimuovo i punti allineati entro la tolleranza e non più lontani tra loro del massimo
   double dSqMaxLen = dMaxSegmLen * dMaxSegmLen ;
   double dSqTol = dLinTol * dLinTol ;
   const double LENREF = 100 ;
   double dCosAngLim = 1 - dSqTol / ( 2 * LENREF * LENREF) ;
   int nPrec = 0 ;
   int nCurr = 1 ;
   int nNext = 2 ;
   while ( nNext < int( vMyPt5ax.size())) {
      bool bRemove = false ;
     // lunghezza del segmento che unisce gli adiacenti
      double dSqLen = SqDist( vMyPt5ax[nPrec].ptP, vMyPt5ax[nNext].ptP) ;
     // se lunghezza inferiore al massimo, passo agli altri controlli
      if ( dSqLen <= dSqMaxLen) {
        // distanza del punto corrente dal segmento che unisce gli adiacenti
         DistPointLine dPL( vMyPt5ax[nCurr].ptP, vMyPt5ax[nPrec].ptP, vMyPt5ax[nNext].ptP) ;
         double dSqDist ;
        // se distanza inferiore a tolleranza lineare
         if ( dPL.GetSqDist( dSqDist)  &&  dSqDist < dSqTol  &&  PointsInTolerance( vMyPt5ax, nPrec, nCurr, nNext, dSqTol)) {
           // verifico se errore angolare inferiore a limite
            double dPar ; dPL.GetParamAtMinDistPoint( dPar) ;
            if ( bTestDir) {
               Vector3d vtNew = Media( vMyPt5ax[nPrec].vtDir, vMyPt5ax[nNext].vtDir, dPar) ;
               if ( vtNew.Normalize()  &&  vtNew * vMyPt5ax[nCurr].vtDir > dCosAngLim)
                  bRemove = true ;
            }
            else {
               Vector3d vtNew = Media( vMyPt5ax[nPrec].vtDir2, vMyPt5ax[nNext].vtDir2, dPar) ;
               if ( vtNew.Normalize()  &&  vtNew * vMyPt5ax[nCurr].vtDir2 > dCosAngLim)
                  bRemove = true ;
            }
         }
      }
     // se da eliminare
      if ( bRemove) {   
         // dichiaro da eliminare il punto
         vMyPt5ax[nCurr].nFlag = -1 ;
        // avanzo con corrente e successivo
         nCurr = nNext ;
         ++ nNext ;
      }
     // altrimenti da tenere
      else {
        // avanzo il terzetto di uno step
         nPrec = nCurr ;
         nCurr = nNext ;
         ++ nNext ;
      }
   }
   return true ;
}