//----------------------------------------------------------------------------
//  EgalTech 2013-2014
//----------------------------------------------------------------------------
// File : CurveByInterp.cpp         Data : 05.08.14        Versione : 1.5h1
// Contenuto : Implementazione della classe CurveByInterp, per creare
//             una curva mediante interpolazione di punti.
//
//
// Modifiche : 05.08.14 DS Creazione modulo.
//
//
//----------------------------------------------------------------------------

//--------------------------- Include ----------------------------------------
#include "stdafx.h"
#include "CalcDerivate.h"
#include "CurveBezier.h"
#include "CurveComposite.h"
#include "/EgtDev/Include/EGkCurveByInterp.h"
#include "/EgtDev/Include/EGkBiArcs.h"
#include "/EgtDev/Include/EgtPointerOwner.h"


//----------------------------------------------------------------------------
bool
CurveByInterp::Reset( void)
{
  // pulisco i diversi vettori
   m_vPnt.clear() ;
   m_vPar.clear() ;
   m_vPrevDer.clear() ;
   m_vNextDer.clear() ;
   return true ;
}

//----------------------------------------------------------------------------
bool
CurveByInterp::AddPoint( const Point3d& ptP)
{
  // se il punto coincide con il precedente, lo salto
   if ( ! m_vPnt.empty()  &&  AreSamePointApprox( ptP, m_vPnt.back()))
      return false ;
  // aggiungo il punto
   try { m_vPnt.push_back( ptP) ; }
   catch ( ...) { return false ; }
   return true ;
}

//----------------------------------------------------------------------------
ICurve*
CurveByInterp::GetCurve( int nMethod, int nType)
{
  // calcolo le tangenti
   if ( nMethod == BESSEL) {
      if ( ! CalcBesselTangents())
         return nullptr ;
   }
   else if ( nType != CUBIC_BEZIERS_LONG) {
      if ( ! CalcAkimaTangents( nMethod == AKIMA_CORNER))
         return nullptr ;
   }

  // se richiesti biarchi
   if ( nType == BIARCS) {
     // creo la curva composita
      PtrOwner<CurveComposite> pCrvCompo( CreateBasicCurveComposite()) ;
      if ( IsNull( pCrvCompo))
         return nullptr ;
     // ciclo sugli intervalli
      for ( int i = 1 ; i < int( m_vPnt.size()) ; ++ i) {
        // creo un biarco per ogni intervallo
         double dDirStartDeg ;
         m_vNextDer[i-1].ToSpherical( nullptr, nullptr, &dDirStartDeg) ;
         double dDirEndDeg ;
         m_vPrevDer[i].ToSpherical( nullptr, nullptr, &dDirEndDeg) ;
         ICurve* pCrv = GetBiArc( m_vPnt[i-1], dDirStartDeg, m_vPnt[i], dDirEndDeg, 0.5) ;
         if ( ! pCrvCompo->AddCurve( pCrv))
            return nullptr ;
      }
      return ::Release( pCrvCompo) ;
   }

  // se richieste curve di Bezier cubiche
   if ( nType == CUBIC_BEZIERS) {
     // creo la curva composita
      PtrOwner<CurveComposite> pCrvCompo( CreateBasicCurveComposite()) ;
      if ( IsNull( pCrvCompo))
         return nullptr ;
     // ciclo sugli intervalli
      for ( int i = 1 ; i < int( m_vPnt.size()) ; ++ i) {
        // creo una curva di Bezier cubica per ogni intervallo
         PtrOwner<CurveBezier> pCBez( CreateBasicCurveBezier()) ;
         if ( IsNull( pCBez))
            return nullptr ;
         pCBez->Init( 3, false) ;
         pCBez->SetControlPoint( 0, m_vPnt[i-1]) ;
         pCBez->SetControlPoint( 1, m_vPnt[i-1] + ( m_vPar[i] - m_vPar[i-1]) / 3 * m_vNextDer[i-1]) ;
         pCBez->SetControlPoint( 2, m_vPnt[i] - ( m_vPar[i] - m_vPar[i-1]) / 3 * m_vPrevDer[i]) ;
         pCBez->SetControlPoint( 3, m_vPnt[i]) ;
         if ( ! pCrvCompo->AddCurve( ::Release( pCBez)))
            return nullptr ;
      }
      return ::Release( pCrvCompo) ;
   }

   // se richieste curve di Bezier cubiche (ottenute da interpolazione con Nurbs)
   if ( nType == CUBIC_BEZIERS_LONG) {
       // creo la curva composita
       PtrOwner<ICurve> pCrv ;
       pCrv.Set( InterpolatePointSetWithBezier( m_vPnt, 50 * EPS_SMALL, 50)) ;
       if ( IsNull(pCrv)  ||  ! pCrv->IsValid())
          return nullptr ;
       return Release( pCrv) ;
   }

   return nullptr ;
}

//----------------------------------------------------------------------------
bool
CurveByInterp::CalcAkimaTangents( bool bDetectCorner)
{
  // pulisco i vettori dei parametri e delle tangenti
   m_vPar.clear() ;
   m_vPrevDer.clear() ;
   m_vNextDer.clear() ;

  // numero di punti
   int nSize = int( m_vPnt.size()) ;

  // sono necessari almeno due punti
   if ( nSize < 2)
      return false ;

  // calcolo le distanze tra i punti per derivarne i parametri
   m_vPar.reserve( nSize) ;
   double dPar = 0 ;
   m_vPar.push_back( dPar) ;
   for ( int i = 1 ; i < nSize ; ++ i) {
      double dDist = Dist( m_vPnt[i-1], m_vPnt[i]) ;
      dPar += dDist ;
      m_vPar.push_back( dPar) ;
   }

  // calcolo le derivate
   m_vPrevDer.reserve( nSize) ;
   m_vNextDer.reserve( nSize) ;
  // se ci sono solo 2 punti, le tangenti devono essere dirette lungo la linea che li unisce
   if ( nSize == 2) {
     // non esiste derivata prima del primo punto
      m_vPrevDer.emplace_back( 0, 0, 0) ;
      m_vNextDer.push_back( ( m_vPnt[1] - m_vPnt[0]) / ( m_vPar[1] - m_vPar[0])) ;
      m_vPrevDer.push_back( m_vNextDer[0]) ;
     // non esiste derivata dopo il secondo e ultimo punto
      m_vNextDer.emplace_back( 0, 0, 0) ;
      return true ;
   }
  // verifico se curva chiusa (primo e ultimo punto coincidono)
   bool bClosed = AreSamePointApprox( m_vPnt.front(), m_vPnt.back()) ;
  // calcolo le derivate
   for ( int i = 0 ; i < nSize ; ++ i) {
      Vector3d vtPrevDer ;
      Vector3d vtNextDer ;
     // primo punto
      if ( i == 0) {
        // se curva chiusa, come precedente uso il penultimo punto
         if ( bClosed) {
           // se non ci sono almeno 5 punti
            if ( nSize < 5) {
               if ( ! CalcCircleMidDer( m_vPar[nSize-2] - m_vPar[nSize-1], m_vPnt[nSize-2], m_vPar[i], m_vPnt[i],
                                        m_vPar[i+1], m_vPnt[i+1], vtNextDer))
                  return false ;
               vtPrevDer = vtNextDer ;
            }
           // altrimenti
            else {
               if ( ! CalcAkimaMidDer( m_vPar[nSize-3] - m_vPar[nSize-1], m_vPnt[nSize-3], m_vPar[nSize-2] - m_vPar[nSize-1], m_vPnt[nSize-2],
                                       m_vPar[i], m_vPnt[i], m_vPar[i+1], m_vPnt[i+1],
                                       m_vPar[i+2], m_vPnt[i+2], bDetectCorner,
                                       vtPrevDer, vtNextDer))
                  return false ;
            }
         }
        // altrimenti, uso arco sui primi tre punti
         else {
            if ( ! CalcCircleStartDer( m_vPar[i], m_vPnt[i], m_vPar[i+1], m_vPnt[i+1],
                                       m_vPar[i+2], m_vPnt[i+2], vtNextDer))
               return false ;
            vtPrevDer = Vector3d( 0, 0, 0) ;
         }
      }
     // ultimo punto
      else if ( i == nSize - 1) {
        // se curva chiusa, le tg devono coincidere con quelle del primo
         if ( bClosed) {
            vtPrevDer = m_vPrevDer[0] ;
            vtNextDer = m_vNextDer[0] ;
         }
        // altrimenti, uso arco sugli ultimi tre punti
         else {
            if ( ! CalcCircleEndDer( m_vPar[i-2], m_vPnt[i-2], m_vPar[i-1], m_vPnt[i-1],
                                       m_vPar[i], m_vPnt[i], vtPrevDer))
               return false ;
            vtNextDer = Vector3d( 0, 0, 0) ;
         }
      }
     // punti intermedi
      else {
        // se secondo punto
         if ( i == 1) {
           // se curva aperta o non ci sono almeno 5 punti
            if ( ! bClosed  ||  nSize < 5) {
               if ( ! CalcCircleMidDer( m_vPar[i-1], m_vPnt[i-1], m_vPar[i], m_vPnt[i],
                                        m_vPar[i+1], m_vPnt[i+1], vtPrevDer))
                  return false ;
               vtNextDer = vtPrevDer ;
            }
           // altrimenti
            else {
               if ( ! CalcAkimaMidDer( m_vPar[nSize-2] - m_vPar[nSize-1], m_vPnt[nSize-2], m_vPar[i-1], m_vPnt[i-1],
                                       m_vPar[i], m_vPnt[i], m_vPar[i+1], m_vPnt[i+1],
                                       m_vPar[i+2], m_vPnt[i+2], bDetectCorner,
                                       vtPrevDer, vtNextDer))
                  return false ;
            }
         }
        // se penultimo punto
         else if ( i == nSize - 2) {
           // se curva aperta o non ci sono almeno 5 punti
            if ( ! bClosed  ||  nSize < 5) {
               if ( ! CalcCircleMidDer( m_vPar[i-1], m_vPnt[i-1], m_vPar[i], m_vPnt[i],
                                        m_vPar[i+1], m_vPnt[i+1], vtPrevDer))
                  return false ;
               vtNextDer = vtPrevDer ;
            }
           // altrimenti
            else {
               if ( ! CalcAkimaMidDer( m_vPar[i-2], m_vPnt[i-2], m_vPar[i-1], m_vPnt[i-1],
                                       m_vPar[i], m_vPnt[i], m_vPar[i+1], m_vPnt[i+1],
                                       m_vPar[1] + m_vPar[i+1], m_vPnt[1], bDetectCorner,
                                       vtPrevDer, vtNextDer))
                  return false ;
            }
         }
        // altrimenti
         else {
            if ( ! CalcAkimaMidDer( m_vPar[i-2], m_vPnt[i-2], m_vPar[i-1], m_vPnt[i-1],
                                    m_vPar[i], m_vPnt[i], m_vPar[i+1], m_vPnt[i+1],
                                    m_vPar[i+2], m_vPnt[i+2], bDetectCorner,
                                    vtPrevDer, vtNextDer))
               return false ;
         }
      }
     // salvo la derivata
      m_vPrevDer.push_back( vtPrevDer) ;
      m_vNextDer.push_back( vtNextDer) ;
   }

   return true ;
}

//----------------------------------------------------------------------------
bool
CurveByInterp::CalcBesselTangents( void)
{
  // pulisco i vettori dei parametri e delle tangenti
   m_vPar.clear() ;
   m_vPrevDer.clear() ;
   m_vNextDer.clear() ;

  // numero di punti
   int nSize = int( m_vPnt.size()) ;

  // sono necessari almeno due punti
   if ( nSize < 2)
      return false ;

  // calcolo le distanze tra i punti per derivarne i parametri
   m_vPar.reserve( nSize) ;
   double dPar = 0 ;
   m_vPar.push_back( dPar) ;
   for ( int i = 1 ; i < nSize ; ++ i) {
      double dDist = Dist( m_vPnt[i-1], m_vPnt[i]) ;
      dPar += dDist ;
      m_vPar.push_back( dPar) ;
   }

  // calcolo le derivate
   m_vPrevDer.reserve( nSize) ;
   m_vNextDer.reserve( nSize) ;
  // se ci sono solo 2 punti, le tangenti devono essere dirette lungo la linea che li unisce
   if ( nSize == 2) {
     // non esiste derivata prima del primo punto
      m_vPrevDer.emplace_back( 0, 0, 0) ;
      m_vNextDer.push_back( ( m_vPnt[1] - m_vPnt[0]) / ( m_vPar[1] - m_vPar[0])) ;
      m_vPrevDer.push_back( m_vNextDer[0]) ;
     // non esiste derivata dopo il secondo e ultimo punto
      m_vNextDer.emplace_back( 0, 0, 0) ;
      return true ;
   }
  // verifico se curva chiusa (primo e ultimo punto coincidono)
   bool bClosed = AreSamePointApprox( m_vPnt.front(), m_vPnt.back()) ;
  // calcolo le derivate
   for ( int i = 0 ; i < nSize ; ++ i) {
      Vector3d vtPrevDer ;
      Vector3d vtNextDer ;
     // primo punto
      if ( i == 0) {
        // se curva chiusa, come precedente uso il penultimo punto
         if ( bClosed) {
            if ( ! CalcBesselMidDer( m_vPar[nSize-2] - m_vPar[nSize-1], m_vPnt[nSize-2], m_vPar[i], m_vPnt[i],
                                     m_vPar[i+1], m_vPnt[i+1], vtNextDer))
               return false ;
            vtPrevDer = vtNextDer ;
         }
        // altrimenti, uso i primi tre punti
         else {
            if ( ! CalcBesselStartDer( m_vPar[i], m_vPnt[i], m_vPar[i+1], m_vPnt[i+1],
                                       m_vPar[i+2], m_vPnt[i+2], vtNextDer))
               return false ;
            vtPrevDer = Vector3d( 0, 0, 0) ;
         }
      }
     // ultimo punto
      else if ( i == nSize - 1) {
        // se curva chiusa, le tg devono coincidere con quelle del primo
         if ( bClosed) {
            vtPrevDer = m_vPrevDer[0] ;
            vtNextDer = m_vNextDer[0] ;
         }
        // altrimenti, uso gli ultimi tre punti
         else {
            if ( ! CalcBesselEndDer( m_vPar[i-2], m_vPnt[i-2], m_vPar[i-1], m_vPnt[i-1],
                                     m_vPar[i], m_vPnt[i], vtPrevDer))
               return false ;
            vtNextDer = Vector3d( 0, 0, 0) ;
         }
      }
     // punti intermedi
      else {
         if ( ! CalcBesselMidDer( m_vPar[i-1], m_vPnt[i-1], m_vPar[i], m_vPnt[i],
                                  m_vPar[i+1], m_vPnt[i+1], vtPrevDer))
            return false ;
         vtNextDer = vtPrevDer ;
      }
     // salvo la derivata
      m_vPrevDer.push_back( vtPrevDer) ;
      m_vNextDer.push_back( vtNextDer) ;
   }

   return true ;
}