NurbsBase.hh

/*
 * NurbsBase.hh
 *
 *  Created on: Jul 1, 2010
 *      Author: jordt
 */

#ifndef NURBSBASE_HH_
#define NURBSBASE_HH_

#include <vector>

//The following switches are meant to be for debuging and experimenting:

// Evaluate the nurbs recursivly according to their mathematical
// definition (very slow)
//#define NURBS_BASE_RECURSIVE_EVALUATION

// This switch makes the evaluation a little slower but thread safe.
// To avoid this draw back call the BaseFunc with an additional Cache
// pointer (This is fast and thread safe)
#define NURBS_BASE_THREAD_SAFE

namespace BIAS{

  /** Base class for NURBS classes. The NURBS Base function is implemented here
    * @author jordt
    */
  class NurbsBase {
  public:

    /** Standard Constructor
      * @param knotsPerDim The number of knots in the knot vector (without the dublications)
      * @param degree The degree of the b-splines
      * @param dimension The dimension of the parameter space! (not the control point space)
      */
    NurbsBase(int knotsPerDim = 4, int degree = 2, int dimension = 1);

    /// Empty Destructor (Virtual to make this whole thing polymorphic)
    ~NurbsBase(){}

    /// Returns the number of knots of the NURBS
    int GetNumKnots() const;

    /// Returns the degree of the NURBS
    int GetDegree() const;

    /// Returns the number of control points per dimension
    int GetNumCtrlPts() const;

    /// Returns the position of a knot with a given index
    /// of a given dimension.
    double GetKnot(int dim, int index) const;

    /// Set the position of a knot with a given index and
    /// of a given dimension
    void SetKnot(int dim, int index, double value);

  protected:

    /// The vector of knot vectors (One for each dimension)
    std::vector<std::vector<double> > knotVector_;

    /// Number of knots in each dimension
    int knotsPerDim_;

    /// Degree of the b-splines
    int degree_;

    /// The number of control points per dimension
    int numCtrlPoints_;

    #ifndef NURBS_BASE_THREAD_SAFE
      /// The function value cache for the serial evaluation
      std::vector<double> evalCache_;
    #endif //NURBS_BASE_THREAD_SAFE

    int baseOps_;

  #ifdef NURBS_BASE_RECURSIVE_EVALUATION
    /// This is standard implementation of the recursive nurbs evaluation.
    inline double BaseFunc(const int &dim,const  int &knot,const int &degree,const double &u) const{
      if (degree == 0){
        if (u < knotVector_[dim][knot]) return 0.0;
        if (u > knotVector_[dim][knot+1]) return 0.0;
        if (u == knotVector_[dim][knot+1] && u != 1.0) return 0.0;
        return 1.0;
      }

      double recursiveA = 0.0;
      if (knotVector_[dim][knot+degree] != knotVector_[dim][knot])
       recursiveA = (u - knotVector_[dim][knot]) / (knotVector_[dim][knot+degree]-knotVector_[dim][knot]) * BaseFunc(dim,knot,degree - 1,u);

      double recursiveB  = 0.0;
      if (knotVector_[dim][knot+degree+1] != knotVector_[dim][knot+1])
         recursiveB = (knotVector_[dim][knot+degree+1]-u) / (knotVector_[dim][knot+degree+1]-knotVector_[dim][knot+1]) * BaseFunc(dim,knot+1,degree - 1,u);

      return recursiveA + recursiveB;
    }
  #else //NURBS_BASE_RECURSIVE_EVALUATION
     /** This is a serial implementation of the recursive nurbs evaluation.
       * Recursive function values are stored in the evaluation cache (evalCache_)
       * @param dim The dimension (determines which knot vector is used)
       * @param knot The index of the knot in the knot vector
       * @param degree The current degree the function is evaluated in
       * @param u The parmeter for which the base function is evaluated */

    inline double BaseFunc(const int &dim, const int &knot, const int &degree, const double &u) const{

      #ifdef NURBS_BASE_THREAD_SAFE
        double* evalCache_ = new double[degree+1];
      #endif

      // baseOps: The number of function calls on this level of recursion.
      int baseOps = degree+1;
      for(int k = 0; k < baseOps; k++)
        evalCache_[k] = (u < knotVector_[dim][knot+k] || u > knotVector_[dim][knot+k+1] || (u == knotVector_[dim][knot+k+1] && u!=1.0) ? 0.0 : 1.0);

      // After computing the deepest recursion level, the levels based on the cache values are rendered:
      for(int level = 1; level <= degree; level++){
        baseOps = baseOps -1;
        for(int k = 0; k < baseOps; k++){
          evalCache_[k] =
              ((knotVector_[dim][knot+level+k] == knotVector_[dim][knot+k]) ? 0.0 :(u - knotVector_[dim][knot+k]) / (knotVector_[dim][knot+level+k]-knotVector_[dim][knot+k]) * evalCache_[k])
         +
              ((knotVector_[dim][knot+level+1+k] == knotVector_[dim][knot+1+k]) ? 0.0 :(knotVector_[dim][knot+level+1+k]-u) / (knotVector_[dim][knot+level+1+k]-knotVector_[dim][knot+1+k]) * evalCache_[k+1]);
        }
      }

      // The result then lies in the first cache array element:
       #ifdef NURBS_BASE_THREAD_SAFE
        double returnValue = evalCache_[0];
        delete[] evalCache_;
        return returnValue;
      #else
        return evalCache_[0];
      #endif

    }
  #endif //NURBS_BASE_RECURSIVE_EVALUATION


    /** This is a serial implementation of the recursive nurbs evaluation.
      * Recursive function values are stored in the evaluation cache (evalCache_)
      * In this function, the cache pointer has to be provided in order to save the
      * memory allocation time.
      * @param dim The dimension (determines which knot vector is used)
      * @param knot The index of the knot in the knot vector
      * @param degree The current degree the function is evaluated in
      * @param u The parmeter for which the base function is evaluated
      * @param cachePointer The pointer to a cache, this function is allowed to use (already allocated). The size should be degree+1*/
    inline double BaseFunc(const int &dim, const int &knot, const int &degree, const double &u, double* &cachePointer) const{

      // baseOps: The number of function calls on this level of recursion.
      int baseOps = degree+1;
      for(int k = 0; k < baseOps; k++)
        cachePointer[k] = (u < knotVector_[dim][knot+k] || u > knotVector_[dim][knot+k+1] || (u == knotVector_[dim][knot+k+1] && u!=1.0)? 0.0 : 1.0);

      // After computing the deepest recursion level, the levels based on the cache values are rendered:
      for(int level = 1; level <= degree; level++){
        baseOps = baseOps -1;
        for(int k = 0; k < baseOps; k++){
          cachePointer[k] =
              ((knotVector_[dim][knot+level+k] == knotVector_[dim][knot+k]) ? 0.0 :(u - knotVector_[dim][knot+k]) / (knotVector_[dim][knot+level+k]-knotVector_[dim][knot+k]) * cachePointer[k])
         +
              ((knotVector_[dim][knot+level+1+k] == knotVector_[dim][knot+1+k]) ? 0.0 :(knotVector_[dim][knot+level+1+k]-u) / (knotVector_[dim][knot+level+1+k]-knotVector_[dim][knot+1+k]) * cachePointer[k+1]);
        }
      }

      // The result then lies in the first cache array element:
      return cachePointer[0];
    }

  };

}

#endif /* NURBSBASE_HH_ */