Interpolator.hh

/*
This file is part of the BIAS library (Basic ImageAlgorithmS).

Copyright (C) 2003-2009    (see file CONTACT for details)
Multimediale Systeme der Informationsverarbeitung
Institut fuer Informatik
Christian-Albrechts-Universitaet Kiel


BIAS is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.

BIAS is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with BIAS; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/


#ifndef _BIAS_INTERPOLATOR_
#define _BIAS_INTERPOLATOR_

#include <bias_config.h>

#include <math.h>

#include <Base/Common/BIASpragmaStart.hh>

#include <vector>
#include "Base/Debug/Debug.hh"
#include "Base/Math/Matrix.hh"
#include "Base/Math/Vector.hh"
#include "Base/Math/Vector2.hh"
#include "Base/Math/Vector3.hh"
#include "Base/Image/Image.hh"


namespace BIAS {
  /** @class Interpolator
      @ingroup g_mathalgo
      @brief this class interpolates a function y=f(t)
      between given control points (the y-values)

      this class interpolates a function y=f(t)
      with different algorithms between given control points (the y-values),
      which can be one, two or three dimensional.
      Interpolation means, the curve reaches all control points.

      usage:
      -# SetControlPoints()

      -# optional: SetKnotPoints(), if not set the t values range form 0 to 1
      such that the controlpoints are uniformly spaced

      -# call a Interpolation function, like Spline() or Bezier3()

      see also 'Computer Graphics & Geometric Modeling' from David Salomon
      @author Daniel Grest, Juli 2002
  */
#define D_INTERPOL_SPLINE3 1<<0
#define D_INTERPOL_SPLINE2 1<<1
#define D_INTERPOL_COMMON  1<<2
#define D_INTERPOL_BEZIER  1<<3

  class BIASMathAlgo_EXPORT Interpolator : public Debug
  {

  public:
    Interpolator();
    Interpolator(const Interpolator &ip)
    { operator=(ip); }
    ~Interpolator();

    /** resets everything.
        clears all lists and calls constructor
    */
    void Clear();

    /** set the control points, which control the interpolating curve
     */
    void SetControlPoints(const std::vector<double> &cPnt1){
      listPntDim1_=cPnt1;};
    void SetControlPoints(const std::vector<BIAS::Vector2<double> > &cPnt);
    void SetControlPoints(const std::vector<BIAS::Vector3<double> > &cPnt);
    void GetControlPoints(std::vector<double> &cPnt) const {
      cPnt=listPntDim1_;};
    void GetControlPoints(std::vector<BIAS::Vector2<double> > &cPnt) const;
    void GetControlPoints(std::vector<BIAS::Vector3<double> > &cPnt) const;

    /** set the additional knot points,
        if you want nonuniform interpolation.
        Otherwise they will be uniformly spaced.
        The values in kPnt must be increasing,
        such that kPnt[i]<kPnt[i+1]
    */
    void SetKnotPoints(const std::vector<double> &kPnt);
    void GetKnotPoints(std::vector<double> &kPnt) const{
      kPnt=KPts_;};


    inline void SetStartTangent(double startTangent){
      startTangentX_=startTangent;
    };
    inline void SetStartTangent(BIAS::Vector2<double> startTangent){
      startTangentX_=startTangent[0];startTangentY_=startTangent[1];
    };
    inline void SetStartTangent(BIAS::Vector3<double> startTangent){
      startTangentX_=startTangent[0]; startTangentY_=startTangent[1];
      startTangentZ_=startTangent[2];
    };
    inline void ClearStartTangent(){
      startTangentX_=startTangentY_=startTangentZ_=-10000;
    };

    inline void SetEndTangent(double endTangent){endTangentX_=endTangent;};
    inline void SetEndTangent(BIAS::Vector2<double> endTangent){
      endTangentX_=endTangent[0];endTangentY_=endTangent[1];};
    inline void SetEndTangent(BIAS::Vector3<double> endTangent){
      endTangentX_=endTangent[0]; endTangentY_=endTangent[1];
      endTangentZ_=endTangent[2];
    };
    inline void ClearEndTangent(){
      endTangentX_=endTangentY_=endTangentZ_=-10000;
    };

    /** for hermite splines, a tangent can be given for each control point.
     */
    void SetTangentsHermite(const std::vector<double> &cTan1){
      listTangDim1_=cTan1;};

    /** call this for restart at t= first knot point
        initiates recalculation of all polynom coefficients
        at first call of Spline()
    */
    void InitSpline();


    /** these functions do the Spline interpolation
        which reaches each control point.
        k is the degree of the interpolation polynom.
        implemented for k=1,2,3,4 (linear, quadratic splines,
        cubic splines, akima cubic splines)
        If no start or end tangents are set, the end condition is
        relaxed for cubic and akima splines.
        for k=5 and a 1D function hermite splines are implemented (todo for vector2/3).
    */
    int Spline(double &res, double t, unsigned int k=3);
    int Spline(BIAS::Vector2<double> &res, double t, unsigned int k=3);
    int Spline(BIAS::Vector3<double> &res, double t, unsigned int k=3);


    /** here the calculation of all coefficients for the bezier interpolation
        is done and the nmatrix is set
        is called by Bezier3() if not called manually
        @author Ingo Schiller
        @param dimension of Control Points
    */
    void InitBezier(int dim_of_CP);


    /** these functions do the bezier interpolation as described
        in David Salomon 4.14.19 which reaches each control point
        /////////////////////////////////////////////
        IMPORTANT NOTICE:
        if the knotpoints are not set, the value of t must be between 0 and 1!
        if t > 1 the value of the last controlpoint is returned!
        /////////////////////////////////////////////
        calculates P(t)= tvector(t)*nmatrix*cpoints
        @author Ingo Schiller
        @param reference to the result(return value)
        @param time t of which value is to be calculated
    */

    int Bezier3(double &res, double t);
    int Bezier3(BIAS::Vector2<double> &res, double t);
    int Bezier3(BIAS::Vector3<double> &res, double t);

    void DrawBezierCoords2(Image<unsigned char> &img);

    /**
       @brief Prepare Look-up-Table in range [min,max] with N samples

       Prepare a Look-up-Table for spline interpolation in the range
       [min,max] with N samples. Use SplineFromLuT() afterwards;
       @author evers
       @date 2005-04-15
    */
    int PrepareLuTSpline(double min, double max, unsigned int N,
                         unsigned int k=3);

    /**
       @brief Get spline interpolation using Look-up-Table.
       return code is negativ if argument is out of bound os LuT
       @author evers
       @date 2005-04-15
    */
    int SplineFromLuT(double &res, double t) const;

    /**
     * Compare LuT data structures if matching perfectly.
     * Used for ProjectionParametersSpherical compare.
     * @author bartczak
     * @date 06/2006
     */
    const bool DoLuTSplinesDiffer(const Interpolator& I) const;

  protected:
    void InitSpline(const unsigned int k, const unsigned int dim);
    void InitLinear(std::vector<BIAS::Vector<double> > &listP,
                    const std::vector<double> &listPnt);
    // precalculates all coefficients for the segments
    // for polynoms of degree 3
    void InitSpline3(std::vector<BIAS::Vector<double> > &listP,
                     const std::vector<double>& listPnt,
                     const double startTangent=0, const double endTangent=0);

    // precalculates all coefficients for the segments
    // for polynoms of degree 2
    void InitSpline2(std::vector<BIAS::Vector<double> > &listP,
                     const std::vector<double>& listPnt,
                     double startTangent=0);

    // precalculates all coefficients for an akima spline
    // of degree 3
    void InitSplineAkima(std::vector<BIAS::Vector<double> > &listP,
                         const std::vector<double>& listPnt,
                         const double startTangent,const double endTangent);

    void InitSplineHermite(std::vector<double> &tangents,
                           const std::vector<double>& listPnt,
                           const double startTangent, const double endTangent);

    // validates important conditions
    bool validate(double &t, bool debug=false);

    //subroutine for akima spline
    void fakima(BIAS::Vector<double> &a,
                BIAS::Vector<double> Px, BIAS::Vector<double> Py,
                int from, int to, int intervall);

    std::vector<double> KPts_; // knot values for control points
    std::vector<double> T_; // knot value intervals for control points
    //double tOffset_;  // T_[0] has to be 0, for coefficient calculation

    std::vector<double> listPntDim1_; // control points dimension 1
    std::vector<double> listPntDim2_; // control points dimension 2
    std::vector<double> listPntDim3_; // control points dimension 3

    std::vector<double> listTangDim1_; // hermite tangents dimension 1;

    // list of polynom coefficients for dimension 1,2,3
    std::vector<BIAS::Vector<double> > listPolynoms1_;
    std::vector<BIAS::Vector<double> > listPolynoms2_;
    std::vector<BIAS::Vector<double> > listPolynoms3_;

    // lists of bezier coefficients for all dimensions
    std::vector<double> bXK_1_;
    std::vector<double> bXK_2_;
    std::vector<double> bXK_3_;

    std::vector<double> bYK_1_;
    std::vector<double> bYK_2_;
    std::vector<double> bYK_3_;
    //nmatrix for the calcutaion of the bezier interpolation
    BIAS::Matrix<double> nmatrix_;

    double startTangentX_;
    double startTangentY_;
    double startTangentZ_;
    double endTangentX_;
    double endTangentY_;
    double endTangentZ_;

    int actIndex_; // the index of the actual intervall for t
    // the degree of the polynom for which  was initialized
    unsigned int initializedPoly_;
    //indicates if InitBezier has been called
    bool bezierInitialized_;

    std::vector<double> SplineLuT1_;
    double LuTMin1_,LuTMax1_, LuTSampleDist1_;

  };

}

#include <Base/Common/BIASpragmaEnd.hh>

#endif