// Gnurbs - A curve and surface library
// Copyright (C) 2008-2026 Eric Bechet
//
// See the LICENSE file for contributions and license information.
// Please report all bugs and problems to <bechet@cadxfem.org>.
//
// Contributor: Leblanc Christophe.
// cleblancad@gmail.com

#ifndef __REPARAMETRIZATION_H
#define __REPARAMETRIZATION_H

#include "nbspline.h"

#include <stdexcept>
#include <set>
#include <limits>
#include <cmath>

/// \class Reparametrization
/// \brief Reparametrize a b-spline/nurbs by another b-spline.
class Reparametrization
{
private:
  nbspline m_curve; //!< b-spline/nurbs curve to be reparametrized.
  nbspline m_reparametrization; //!< Reparametrization b-spline: contains values at coord 0 of its control points. Should be a bijection.
  bool m_done; //!< True if curve has been reparametrized, false elsewhere.
  bool m_use_inverse; //!< Reparametrization with the inverse function.

public:
  /// \brief Constructor.
  /// \paraml[in] curve curve to be reparametrized.
  /// \param[in] param reparametrization function as a b-spline.
  Reparametrization(const nbspline &curve, const nbspline &param):
    m_curve(curve), m_reparametrization(param), m_done(false), m_use_inverse(false)
  { 
    if(!is_bspline(param)) 
      throw std::runtime_error("Reparmatrization only implemented for b-spline reparametrizing curves."); 
  }

  /// \brief Computes the reparametrization.
  void apply();

  /// \brief Returns the reparametrized curve.
  const nbspline& get()
  { if(!m_done) apply(); return m_curve; }

  /// \brief Sets if the inverse reparametrization function should be used.
  void use_inverse(const bool inv)
  { m_use_inverse = inv; }

  /// \brief Tells if the inverse reparametrization function is ised.
  bool is_inverse_used() const
  { return m_use_inverse; }

private:
  /// \brief Test if the current curve is a b-spline.
  /// \param[in] curve curve to test.
  /// \return True if the current curve is a b-spline, false otherwise (nurbs).
  bool is_bspline(const nbspline &curve) const;

  /// \brief Sets the knots of the reparametrized curve.
  /// \param[in,out] curve reparametrized curve with its nodal sequence updated.
  ///  Warning: control points are not set. see function 'set_control_points'.
  /// \param[out] saturated_curve input curve in form of Bézier segments
  ///  and with the knots ui = f(si) inserted.
  void set_nodal_sequence(nbspline &curve, nbspline &saturated_curve) const;

  /// \brief Returns the distincts internal knots of a curve.
  /// \param[in] curve curve to consider.
  /// \param[out] knots distincts internal knots.
  void get_distinct_internal_knots(const nbspline &curve, 
    std::set<double> &knots) const;

  /// \brief Returns the distincts knots of a curve.
  /// \param[in] curve curve to consider.
  /// \param[out] knots distincts knots.
  void get_distinct_knots(const nbspline &curve, 
    std::set<double> &knots) const;

  /// \brief Compute the 'inverted' knots si = g(ui) = f^{(-1)}(si).
  /// \param[in] ui knots to 'invert'.
  /// \param[out] si 'inverted' knots.
  void get_inverse_knots(const std::set<double> &ui,
    std::vector<double> &si) const;

  /// \brief Sets the control points of the reparametrized curve.
  /// \param[in,out] reparametrized_curve In: curve with reparametrized nodal
  ///   sequence (as done by 'set_nodal_sequence'). Out: control points of 
  ///   the reparametrized curve are computed.
  /// \param[in] saturated_curve saturated curve as returned by 'set_nodal_sequence'.
  void set_control_points(nbspline &reparametrized_curve,
    const nbspline &saturated_curve) const;

  /// \brief Compute equation 6.40 in Nurbs Book, p248.
  /// \param[in] C derivatives of the curve at saturated knots.
  /// \param[in] f derivatives of the reparametrization at saturated (and merged) knots.
  /// \param[in] knot_index knot index value at which C and f derivatives have to be evaluated.
  /// \param[in] n derivative order.
  /// \warning: slow for high degrees (typically > 10).
  /// \return Value of formula 6.40.
  npoint compute_equation640(const std::vector< std::vector<npoint> > &C,
    const std::vector< std::vector<double> > &f, const size_t knot_index,
    const size_t n) const;

  /// \brief Helper function for compute_equation640.
  /// \param[in] C derivatives of the curve at saturated knots.
  /// \param[in] f derivatives of the reparametrization at saturated (and merged) knots.
  /// \param[in] knot_index knot index value at which C and f derivatives have to be evaluated.
  /// \param[in] j current sum index in formula 6.40.
  /// \param[in] k vector of k's value satisfiying the sum constraints in formula 6.40.
  /// \return Value of one term of the sum over the k's of formula 6.40 in Nurbs Book, p248.
  double compute_equation640_helper(const std::vector< std::vector<double> > &f, 
    const int knot_index, const int j, const std::vector<size_t> &k) const;

  /// \brief Checks if a set of integers k1, ..., kn satisfies
  ///   the constraints of equation 6.40 in Nurbs Book, p248.
  /// \param[in] k vector of integers.
  /// \param[in] j constraint bound.
  /// \param[out] smaller_than_n_constraint_satisfied true if the
  ///   The relaxed constraint k_1+2*k_2+...+n*k_n <= n
  ///   (from constraint k_1+2*k_2+...+n*k_n = n) is satisfied.
  ///   False otherwise.
  /// \return True if constraints satisfied, false otherwise.
  inline bool constraints640_satisfied(const std::vector<size_t> &k,
    const size_t j, bool &smaller_than_n_constraint_satisfied) const;

  /// \brief Factorials table for the long double type (from 0! to 170!).
  /// (From the boost library).
  /// \param[in] n input number (between 0 and 170 included).
  /// \return n!
  inline long double factorial(const int n) const;
};

#endif