Splines

/*--------------------------------------------------------------------------*\
 |                                                                          |
 |  Copyright (C) 2016                                                      |
 |                                                                          |
 |         , __                 , __                                        |
 |        /|/  \               /|/  \                                       |
 |         | __/ _   ,_         | __/ _   ,_                                |
 |         |   \|/  /  |  |   | |   \|/  /  |  |   |                        |
 |         |(__/|__/   |_/ \_/|/|(__/|__/   |_/ \_/|/                       |
 |                           /|                   /|                        |
 |                           \|                   \|                        |
 |                                                                          |
 |      Enrico Bertolazzi                                                   |
 |      Dipartimento di Ingegneria Industriale                              |
 |      Universita` degli Studi di Trento                                   |
 |      email: enrico.bertolazzi@unitn.it                                   |
 |                                                                          |
\*--------------------------------------------------------------------------*/
 
/*\
 |   ____        _ _          __     __
 |  / ___| _ __ | (_)_ __   __\ \   / /__  ___
 |  \___ \| '_ \| | | '_ \ / _ \ \ / / _ \/ __|
 |   ___) | |_) | | | | | |  __/\ V /  __/ (__
 |  |____/| .__/|_|_|_| |_|\___| \_/ \___|\___|
 |        |_|
\*/
 
namespace Splines {
 
  class SplineVec {
 
  protected:
 
    #ifndef DOXYGEN_SHOULD_SKIP_THIS
 
    string const m_name;
 
    Utils::Malloc<real_type>  m_mem;
    Utils::Malloc<real_type*> m_mem_p;
 
    integer m_dim{0};
    integer m_npts{0};
    bool    m_curve_is_closed{false};
    bool    m_curve_can_extend{true};
 
    real_type *  m_X{nullptr};
    real_type ** m_Y{nullptr};
    real_type ** m_Yp{nullptr};
 
    #ifdef SPLINES_USE_THREADS
    mutable std::mutex                                         m_last_interval_mutex;
    mutable std::map<std::thread::id,std::shared_ptr<integer>> m_last_interval;
    #else
    mutable integer m_last_interval;
    #endif
 
    void init_last_interval();
    void allocate( integer dim, integer npts );
    void compute_chords();
 
    #endif
 
  public:
 
    SplineVec( SplineVec const & ) = delete;
    SplineVec const & operator = ( SplineVec const & ) = delete;
 
 
    SplineVec( string const & name = "SplineVec" );
 
    virtual
    ~SplineVec();
 
 
    void search( std::pair<integer,real_type> & res ) const;
 
    string const & name() const { return m_name; }
 
 
    bool is_closed() const { return m_curve_is_closed; }
 
    void make_closed() { m_curve_is_closed = true; }
 
    void make_open() { m_curve_is_closed = false; }
 
    bool can_extend() const { return m_curve_can_extend; }
 
    void make_unbounded()   { m_curve_can_extend = true; }
 
    void make_buonded()     { m_curve_can_extend = false; }
 
 
    integer num_points() const { return m_npts; }
 
    integer dimension() const { return m_dim; }
 
    real_type const * x_nodes() const { return m_X; }
 
    real_type x_node( integer npt ) const { return m_X[size_t(npt)]; }
 
    real_type const * y_nodes( integer j ) const { return m_Y[size_t(j)]; }
 
    real_type
    y_node( integer npt, integer j ) const
    { return m_Y[size_t(j)][size_t(npt)]; }
 
    real_type x_min() const { return m_X[0]; }
 
    real_type x_max() const { return m_X[size_t(m_npts-1)]; }
 
 
 
    real_type
    eval( real_type x, integer i ) const;
 
    real_type
    operator () ( real_type x, integer i ) const
    { return this->eval( x, i ); }
 
    real_type
    D( real_type x, integer i ) const;
 
    real_type
    eval_D( real_type x, integer i ) const
    { return this->D(x,i); }
 
    real_type
    DD( real_type x, integer i ) const;
 
    real_type
    eval_DD( real_type x, integer i ) const
    { return this->DD(x,i); }
 
    real_type
    DDD( real_type x, integer i ) const;
 
    real_type
    eval_DDD( real_type x, integer i ) const
    { return this->DDD(x,i); }
 
    real_type
    DDDD( real_type x, integer i ) const;
 
    real_type
    eval_DDDD( real_type x, integer i ) const
    { return this->DDDD(x,i); }
 
    real_type
    DDDDD( real_type x, integer i ) const;
 
    real_type
    eval_DDDDD( real_type x, integer i ) const
    { return this->DDDDD(x,i); }
 
 
 
    void
    eval(
      real_type x,
      real_type vals[],
      integer   inc
    ) const;
 
    void
    eval_D(
      real_type x,
      real_type vals[],
      integer   inc
    ) const;
 
    void
    eval_DD(
      real_type x,
      real_type vals[],
      integer   inc
    ) const;
 
    void
    eval_DDD(
      real_type x,
      real_type vals[],
      integer   inc
    ) const;
 
    void
    eval_DDDD(
      real_type x,
      real_type vals[],
      integer   inc
    ) const;
 
    void
    eval_DDDDD(
      real_type x,
      real_type vals[],
      integer   inc
    ) const;
 
 
    void eval( real_type x, vector<real_type> & vals ) const;
 
    void eval_D( real_type x, vector<real_type> & vals ) const;
 
    void eval_DD( real_type x, vector<real_type> & vals ) const;
 
    void eval_DDD( real_type x, vector<real_type> & vals ) const;
 
    void eval_DDDD( real_type x, vector<real_type> & vals ) const;
 
    void eval_DDDDD( real_type x, vector<real_type> & vals ) const;
 
 
    void
    eval( real_type x, GenericContainer & vals ) const
    { eval( x, vals.set_vec_real(m_dim) ); }
 
    void
    eval_D( real_type x, GenericContainer & vals ) const
    { eval_D( x, vals.set_vec_real(m_dim) ); }
 
    void
    eval_DD( real_type x, GenericContainer & vals ) const
    { eval_DD( x, vals.set_vec_real(m_dim) ); }
 
    void
    eval_DDD( real_type x, GenericContainer & vals ) const
    { eval_DDD( x, vals.set_vec_real(m_dim) ); }
 
    void
    eval_DDDD( real_type x, GenericContainer & vals ) const
    { eval_DDDD( x, vals.set_vec_real(m_dim) ); }
 
    void
    eval_DDDDD( real_type x, GenericContainer & vals ) const
    { eval_DDDDD( x, vals.set_vec_real(m_dim) ); }
 
 
    void
    eval( vec_real_type const & x, GenericContainer & vals ) const;
 
    void
    eval_D( vec_real_type const & x, GenericContainer & vals ) const;
 
    void
    eval_DD( vec_real_type const & x, GenericContainer & vals ) const;
 
    void
    eval_DDD( vec_real_type const & x, GenericContainer & vals ) const;
 
    void
    eval_DDDD( vec_real_type const & x, GenericContainer & vals ) const;
 
    void
    eval_DDDDD( vec_real_type const & x, GenericContainer & vals ) const;
 
 
    void
    setup(
      integer           dim,
      integer           npts,
      real_type const * Y[]
    );
 
    void
    setup(
      integer         dim,
      integer         npts,
      real_type const Y[],
      integer         ldY
    );
 
    void set_knots( real_type const X[] );
 
    void set_knots_chord_length();
 
    void set_knots_centripetal();
 
    void set_knots_foley();
 
    void catmull_rom();
 
    void
    build( GenericContainer const & gc )
    { setup(gc); }
 
    void setup( GenericContainer const & gc );
 
 
    real_type curvature( real_type x ) const;
 
    real_type curvature_D( real_type x ) const;
 
    SplineType1D type() const { return SplineType1D::SPLINE_VEC; }
 
    string info() const;
 
    void
    info( ostream_type & stream ) const
    { stream << this->info() << '\n'; }
 
    void
    dump_table( ostream_type & s, integer num_points ) const;
 
    #ifdef SPLINES_BACK_COMPATIBILITY
    integer numPoints() const { return m_npts; }
    real_type const * xNodes() const { return m_X; }
    real_type xNode( integer npt ) const { return m_X[size_t(npt)]; }
    real_type const * yNodes( integer j ) const { return m_Y[size_t(j)]; }
    real_type yNode( integer npt, integer j ) const { return y_node(npt,j); }
    real_type xMin() const { return m_X[0]; }
    real_type xMax() const { return m_X[size_t(m_npts-1)]; }
    void setKnots( real_type const X[] ) { this->set_knots( X ); }
    void setKnotsChordLength() { this->set_knots_chord_length(); }
    void setKnotsCentripetal() { this->set_knots_centripetal() }
    void setKnotsFoley() { this->set_knots_foley(); }
    void CatmullRom() { this->catmull_rom(); }
    #endif
 
  };
 
}
 
// EOF: SplineVec.hxx