numerics/bezier.h

// -*- c++ -*-

// +--------------------------------------------------------------------+
// | This file is part of MathTL - the Mathematical Template Library    |
// |                                                                    |
// | Copyright (c) 2002-2009                                            |
// | Thorsten Raasch, Manuel Werner                                     |
// +--------------------------------------------------------------------+

#ifndef _MATHTL_BEZIER_H
#define _MATHTL_BEZIER_H

#include <cmath>
#include <utils/function.h>
#include <utils/tiny_tools.h>

namespace MathTL
{
  template <int d>
  double EvaluateBernsteinPolynomial(const int k, const double x)
  {
    return ((k < 0 || k > d)
            ? 0
            : (1-x) * EvaluateBernsteinPolynomial<d-1>(k,  x)
            + x   * EvaluateBernsteinPolynomial<d-1>(k-1,x));
  }
  
  template <>
  double EvaluateBernsteinPolynomial<0>(const int k, const double x)
  {
    return (k == 0 ? 1.0 : 0.0);
  }

  template <int d>
  double EvaluateBernsteinPolynomial_x(const int k, const double x)
  {
    return d * (EvaluateBernsteinPolynomial<d-1>(k-1, x)
                - EvaluateBernsteinPolynomial<d-1>(k, x));
  }

  template <int d>
  double EvaluateBernsteinPolynomial_xx(const int k, const double x)
  {
    return d*(d-1) * (EvaluateBernsteinPolynomial<d-2>(k-2, x)
                      - 2*EvaluateBernsteinPolynomial<d-2>(k-1, x)
                      + EvaluateBernsteinPolynomial<d-2>(k, x));
  }

  inline
  double EvaluateBernsteinRepresentation(const double b0, const double b1, const double b2,
                                         const double x)
  {
    const double b10 = (1.-x)*b0+x*b1;
    const double b11 = (1.-x)*b1+x*b2;
    
    return (1.-x)*b10+x*b11;
  }
  
  inline
  double EvaluateBernsteinRepresentation(const double b0, const double b1,
                                         const double b2, const double b3,
                                         const double x)
  {
    const double b10 (b0+x*(b1-b0)); // = (1.-x)*b0+x*b1;
    const double b11 (b1+x*(b2-b1)); // = (1.-x)*b1+x*b2;
    const double b12 (b2+x*(b3-b2)); // = (1.-x)*b2+x*b3;
    
    const double b20 (b10+x*(b11-b10)); // = (1.-x)*b10+x*b11;
    const double b21 (b11+x*(b12-b11)); // = (1.-x)*b11+x*b12;

    return b20+x*(b21-b20); // = (1.-x)*b20+x*b21;
  }

  inline
  double EvaluateBernsteinRepresentation_x(const double b0, const double b1,
                                           const double b2, const double b3,
                                           const double x)
  {
    return 3*EvaluateBernsteinRepresentation(b1-b0, b2-b1, b3-b2, x);
  }

  inline
  double EvaluateBernsteinRepresentation_xx(const double b0, const double b1,
                                            const double b2, const double b3,
                                            const double x)
  {
    return 6*((1-x)*(b2-2*b1+b0) + x*(b3-2*b2+b1));
  }
  
  double EvaluateHermiteSpline(const int i, const double x) {
    if (i == 0) {
      // phi_0, Bezier coefficients are {0,0,1,1,1,0,0}
      if (x <= -1. || x >= 1.) {
        return 0.;
      } else {
        return (x <= 0.
                ? EvaluateBernsteinRepresentation(0., 0., 1., 1., x+1.)
//              (x+1)*(x+1)*(1-2*x)
                : EvaluateBernsteinRepresentation(1., 1., 0., 0., x)
//              (1-x)*(1-x)*(2*x+1)
                );
      }
    } else {
      // phi_1, Bezier coefficients are {0,0,-1/3,0,1/3,0,0}
      if (x <= -1. || x >= 1.) {
        return 0.;
      } else {
        return (x <= 0.
                ? EvaluateBernsteinRepresentation(0., 0., -1./3., 0., x+1.)
                : EvaluateBernsteinRepresentation(0., 1./3., 0., 0., x));
      }
    }
  }

  inline
  double EvaluateHermiteSpline_td(const int i, const int j, const int k, const double x) {
#if 0
    const double factor(1<<j);
    return sqrt(factor) * EvaluateHermiteSpline(i, factor * x - k);
#else
    return twotothejhalf(j) * EvaluateHermiteSpline(i, (1<<j) * x - k);
#endif
  }
  
  inline
  double EvaluateHermiteSpline_x(const int i, const double x) {
    if (i == 0) {
      // phi_0, Bezier coefficients are {0,0,1,1,1,0,0}
      if (x <= -1 || x >= 1) {
        return 0;
      } else {
        return (x <= 0
                ? EvaluateBernsteinRepresentation_x(0, 0, 1, 1, x+1.0)
                : EvaluateBernsteinRepresentation_x(1, 1, 0, 0, x));
      }
    } else {
      // phi_1, Bezier coefficients are {0,0,-1/3,0,1/3,0,0}
      if (x <= -1 || x >= 1) {
        return 0;
      } else {
        return (x <= 0
                ? EvaluateBernsteinRepresentation_x(0, 0, -1./3., 0, x+1.0)
                : EvaluateBernsteinRepresentation_x(0, 1./3., 0, 0, x));
      }
    }
  }

  inline
  double EvaluateHermiteSpline_xx(const int i, const double x) {
    if (i == 0) {
      // phi_0, Bezier coefficients are {0,0,1,1,1,0,0}
      if (x <= -1 || x >= 1) {
        return 0;
      } else {
        return (x <= 0
                ? EvaluateBernsteinRepresentation_xx(0, 0, 1, 1, x+1.0)
                : EvaluateBernsteinRepresentation_xx(1, 1, 0, 0, x));
      }
    } else {
      // phi_1, Bezier coefficients are {0,0,-1/3,0,1/3,0,0}
      if (x <= -1 || x >= 1) {
        return 0;
      } else {
        return (x <= 0
                ? EvaluateBernsteinRepresentation_xx(0, 0, -1./3., 0, x+1.0)
                : EvaluateBernsteinRepresentation_xx(0, 1./3., 0, 0, x));
      }
    }
  }

  inline
  double EvaluateHermiteSpline_td_x(const int i, const int j, const int k, const double x) {
#if 0
    return twotothejhalf(3*j) * EvaluateHermiteSpline_x(i, (1<<j) * x - k);
#else
    return (double)(1<<j) * (double)twotothejhalf(j) * EvaluateHermiteSpline_x(i, (1<<j) * x - k);
#endif
  }
  
  inline
  double EvaluateHermiteSpline_td_xx(const int i, const int j, const int k, const double x) {
#if 0
    return twotothejhalf(5*j) * EvaluateHermiteSpline_xx(i, (1<<j) * x - k);
#else
    const double factor = (double)(1<<j);
    return factor * factor * twotothejhalf(j) * EvaluateHermiteSpline_xx(i, (1<<j) * x - k);
#endif
  }

  class CubicHermiteInterpolant_td : public Function<1>
  {
  public:
    CubicHermiteInterpolant_td(const int j, const int c, const int k)
      : j_(j), c_(c), k_(k) {}
    
    virtual ~CubicHermiteInterpolant_td() {}
    
    inline double value(const Point<1>& p,
                        const unsigned int component = 0) const
    {
      return EvaluateHermiteSpline_td(c_, j_, k_, p[0]);
    }
    
    void vector_value(const Point<1> &p,
                      Vector<double>& values) const
    {
      values.resize(1, false);
      values[0] = value(p);
    }
    
    void set_j(const int j) { j_ = j; }
    
    void set_c(const int c) { c_ = c; }

    void set_k(const int k) { k_ = k; }
    
  protected:
    int j_, c_, k_;
  };
  
  class CubicHermiteInterpolant2D_td : public Function<2>
  {
  public:
    CubicHermiteInterpolant2D_td(const int j, const int c0, const int c1, const int k0, const int k1)
      : j_(j), c0_(c0), c1_(c1), k0_(k0), k1_(k1) {}
    
    virtual ~CubicHermiteInterpolant2D_td() {}
    
    inline double value(const Point<2>& p,
                        const unsigned int component = 0) const
    {
      return EvaluateHermiteSpline_td(c0_, j_, k0_, p[0])
        * EvaluateHermiteSpline_td(c1_, j_, k1_, p[1]);
    }
    
    void vector_value(const Point<2> &p,
                      Vector<double>& values) const
    {
      values.resize(1, false);
      values[0] = value(p);
    }
    
    void set_j(const int j) { j_ = j; }
    
    void set_c0(const int c0) { c0_ = c0; }

    void set_c1(const int c1) { c1_ = c1; }

    void set_k0(const int k0) { k0_ = k0; }
    
    void set_k1(const int k1) { k1_ = k1; }
    
  protected:
    int j_, c0_, c1_, k0_, k1_;
  };
  
}

#include <numerics/bezier.cpp>

#endif