base_splines.hpp

/* base_splines.hpp */
/* $Id: base_splines.hpp,v 1.3 2007/09/10 08:42:07 guser1 Exp $ */
#ifndef __base_splines_hpp
#define __base_splines_hpp

#include <vector>
#include <math.h>

#include "parspline.hpp"

using namespace std;

#define U3D_EPSILON 1.e-6
#define EQUAL_REALS(a,b) ( (fabs((a)-(b)) < U3D_EPSILON ? true : false) ) 

//------------------------------------------------------------------------------------------
#define U3D_BASE_MSG_LEN 2048

class Spl_U3dMessage
{
public:
  static void info(const char* format, ...);
  
  static void warning(const char* format, ...);
  
  static void error(const char* format, ...);
  
  static void debug(const char* format, ...);
};


//------------------------------------------------------------------------------------------
class Spl_U3dPoint
{
public:
  
  Spl_U3dPoint ();
  
  virtual ~Spl_U3dPoint(){};
  
  Spl_U3dPoint (const Spl_U3dPoint& p);
  
  Spl_U3dPoint (double x, double y, double z);
  
  /* Create point with given point coordinates. */
  //Spl_U3dPoint (Spl_U3dVertex vertex);
  
  Spl_U3dPoint& operator= (const Spl_U3dPoint& p);
  
  int operator==(const Spl_U3dPoint&) const;
  
  int operator!=(const Spl_U3dPoint&) const;
  
  // Point and Vector Operations (always valid) 
  Spl_U3dPoint  operator+(Spl_U3dPoint);      // +translate
  Spl_U3dPoint  operator-(Spl_U3dPoint);      // -translate
  //Spl_U3dPoint& operator+=(Spl_U3dVector);     // inc translate
  //Spl_U3dPoint& operator-=(Spl_U3dVector);     // dec translate
  
  // Scalar Multiplication
  friend Spl_U3dPoint operator*(int, Spl_U3dPoint);
  friend Spl_U3dPoint operator*(Spl_U3dPoint, int);
  friend Spl_U3dPoint operator*(double, Spl_U3dPoint);
  friend Spl_U3dPoint operator*(Spl_U3dPoint, double); 
  
  // Scalar Division
  friend Spl_U3dPoint operator/(Spl_U3dPoint, int);
  friend Spl_U3dPoint operator/(Spl_U3dPoint, double);
  
  double    getX () const;
  
  double    getY () const;
  
  double    getZ () const;
 
  /* Check the point is inside given box. */
  //virtual bool isInside (const Spl_U3dBox& mbox) const;
  
  void      setXYZ (double x, double y, double z);
  
  void      setX (double x);
  
  void      setY (double y);
  
  void      setZ (double z);
  
  //Spl_U3dVertex getVertex();
  
protected:
  
  double    fX, fY, fZ;
  
};


/* INLINE implementation */

inline Spl_U3dPoint::Spl_U3dPoint ()
  : fX(0.0), fY(0.0), fZ(0.0)
{
  /* dummy */
}

inline Spl_U3dPoint::Spl_U3dPoint (const Spl_U3dPoint& p)
  : fX(p.getX()), fY(p.getY()), fZ(p.getZ())
{
  /* dummy */
}

inline Spl_U3dPoint::Spl_U3dPoint (double x, double y, double z)
  : fX(x), fY(y), fZ(z)
{
  /* dummy */
}
/*
  inline Spl_U3dPoint::Spl_U3dPoint (Spl_U3dVertex vertex)
  : fX(vertex.x), fY(vertex.y), fZ(vertex.z)
  {
 
  }*/

inline Spl_U3dPoint&
Spl_U3dPoint::operator= (const Spl_U3dPoint& p)
{
  fX = p.getX();
  fY = p.getY();
  fZ = p.getZ();
  return *this;
}

inline double
Spl_U3dPoint::getX () const
{
  return fX;
}   

inline double
Spl_U3dPoint::getY () const
{
  return fY;
}

inline double
Spl_U3dPoint::getZ () const
{
  return fZ;
}

inline void
Spl_U3dPoint::setXYZ (double x, double y, double z)
{
  fX = x;
  fY = y;
  fZ = z;
}

inline void
Spl_U3dPoint::setX (double x)
{
  fX = x;
}

inline void
Spl_U3dPoint::setY (double y)
{
  fY = y;
}

inline void
Spl_U3dPoint::setZ (double z)
{
  fZ = z;
}

typedef vector<Spl_U3dPoint*> Spl_U3dPointsAr; 
typedef vector<Spl_U3dPointsAr> Spl_U3dPointsMatrix; 


/*=========================== GLOBAL FUNCTIONS ========================*/

inline double pDistance(Spl_U3dPoint p0, Spl_U3dPoint p1)
{
  return sqrt ( (p0.getX()-p1.getX())*(p0.getX()-p1.getX()) + 
        (p0.getY()-p1.getY())*(p0.getY()-p1.getY()) + 
        (p0.getZ()-p1.getZ())*(p0.getZ()-p1.getZ()) );
};

//------------------------------------------------------------------------------------------
typedef std::vector<double> Spl_U3dMathDArray;

typedef std::vector<Spl_U3dMathDArray> Spl_U3dMathDMatrix;

class Spl_U3dMath
{
public:
  static double compositionMul(const Spl_U3dPoint &p1, const Spl_U3dPoint &p2, const Spl_U3dPoint &p3);
  
  static double scalarMul(const Spl_U3dPoint &p1, const Spl_U3dPoint &p2);
  
  static void matrix4x4Mul(const double (&m1)[4][4], const double (&m2)[4][4], double (&m3)[4][4], 
               bool tr1 = false, bool tr2 = false);
  
  static void matrix2vec4x4Mul(const double (&m)[4][4], const double (&u)[4], double (&v)[4], 
                   bool tr = false);
  
  static void vec2matrix4x4Mul(const double (&u)[4], const double (&m)[4][4], double (&v)[4], 
                   bool tr = false);

  static void getLineCoeffs(double x1, double y1, double x2, double y2, 
                double &k, double &b) {
    k = (y1 - y2) / (x1 - x2);
    b = y1 - k*x1;
  }
  static int solveSquareEq(double a, double b, double c, double &x1, double &x2); 

  static double m_d_prec;
};

//------------------------------------------------------------------------------------------
using namespace std;

#define U3D_MATH_BSPL_PARAM_PREC 1.e-10

//------------------------------------------------------------------------------------------
class Spl_U3dBspl3AbstractMatrix
{
public:
  
  Spl_U3dBspl3AbstractMatrix(){m_is_constructed = false;};
  
  Spl_U3dBspl3AbstractMatrix(const Spl_U3dBspl3AbstractMatrix &obj);            
  
  virtual ~Spl_U3dBspl3AbstractMatrix(){};
  
  virtual double getB(int i, double t) const; 
  
  virtual void info() const;
  
  bool isConstructed() const
  {return m_is_constructed;}
  
  const double (&get() const)[4][4];
 
protected:

  bool m_is_constructed;   

  double m_bm[4][4];
};

//------------------------------------------------------------------------------------------
class Spl_U3dBspl3BetaMatrix : public Spl_U3dBspl3AbstractMatrix
{
public:
  
  Spl_U3dBspl3BetaMatrix(double beta1, double beta2);
  
  Spl_U3dBspl3BetaMatrix(const Spl_U3dBspl3BetaMatrix &obj);
  
  virtual ~Spl_U3dBspl3BetaMatrix(){};
  
protected:
  
  void reconstruct(double beta1, double beta2); 
};

//------------------------------------------------------------------------------------------
class Spl_U3dBspl3CRomMatrix : public Spl_U3dBspl3AbstractMatrix
{
public:
  
  Spl_U3dBspl3CRomMatrix(double t);
  
  Spl_U3dBspl3CRomMatrix(const Spl_U3dBspl3CRomMatrix &obj);
  
  virtual ~Spl_U3dBspl3CRomMatrix(){};
  
protected:
  
  void reconstruct(double t); 
};

//------------------------------------------------------------------------------------------
class Spl_U3dBspl3BSMatrix : public Spl_U3dBspl3AbstractMatrix
{
public:
  
  Spl_U3dBspl3BSMatrix();
  
  Spl_U3dBspl3BSMatrix(const Spl_U3dBspl3BSMatrix &obj);
  
  virtual ~Spl_U3dBspl3BSMatrix(){};
  
protected:
  
  void reconstruct(); 
};

//------------------------------------------------------------------------------------------
class Spl_U3dBspl3BezierMatrix : public Spl_U3dBspl3AbstractMatrix
{
public:
  
  Spl_U3dBspl3BezierMatrix();
  
  Spl_U3dBspl3BezierMatrix(const Spl_U3dBspl3BezierMatrix &obj);
  
  virtual ~Spl_U3dBspl3BezierMatrix(){};
  
protected:
  
  void reconstruct(); 
};


//------------------------------------------------------------------------------------------
class Spl_U3dSimpleBspl3Curve
{
public:
  
  Spl_U3dSimpleBspl3Curve(const Spl_U3dBspl3AbstractMatrix *bm, const double (&cp)[4]);
  
  virtual ~Spl_U3dSimpleBspl3Curve(){};
  
  double get(double p) const;
  
  double getSp(double p) const;
  
  double getSpp(double p) const;
  
  void getCoeffs(double (&c)[4]) const {
    c[0] = m_rv[0];
    c[1] = m_rv[1];
    c[2] = m_rv[2];
    c[3] = m_rv[3];
  }
  
  double getMaxLineDeviation(double k, double b, double p1, double p2);
  
  double getLineDeviation(double k, double b, double p);
  
protected:
  
  const Spl_U3dBspl3AbstractMatrix *m_bm;
  
  double m_rv[4];
  
private:

  void cp2rv(const Spl_U3dBspl3AbstractMatrix *bm, const double (&cp)[4]);
};

typedef vector<Spl_U3dSimpleBspl3Curve*> Spl_U3dSimpleBspl3CurveAr;

//------------------------------------------------------------------------------------------
class Spl_U3dSimpleBspl3Surface
{
public:
  
  Spl_U3dSimpleBspl3Surface(const Spl_U3dBspl3AbstractMatrix *bm, const double (&cp)[4][4]);
  
  virtual ~Spl_U3dSimpleBspl3Surface(){};
  
  double get(double p, double q) const;
  
  double getSp(double p, double q) const;
  
  double getSq(double p, double q) const;
  
  double getSpp(double p, double q) const;
  
  double getSqq(double p, double q) const;
  
  double getSpq(double p, double q) const;
  
protected:
  
  const Spl_U3dBspl3AbstractMatrix *m_bm;
  
  double m_rm[4][4];
  
  double calc(double *pp, double *qq) const;
  
private:
  
  void cp2rm(const Spl_U3dBspl3AbstractMatrix *bm, const double (&cp)[4][4]);
};


//------------------------------------------------------------------------------------------
class Spl_U3dSimpleBspl3Surf3d
{
public:
  
  Spl_U3dSimpleBspl3Surf3d(const Spl_U3dBspl3AbstractMatrix *bm, const double (&cp)[4][4][4]);
  
  virtual ~Spl_U3dSimpleBspl3Surf3d(){};
  
  double get(double p, double q, double r) const;
  
protected:
  
  double m_rt[4][4][4];
  
  const Spl_U3dBspl3AbstractMatrix *m_bm;
  
  void make_tensor(const Spl_U3dBspl3AbstractMatrix *bm, const double (&cp)[4][4][4]);

  double calc(double *pp, double *qq, double *rr) const;
};


//------------------------------------------------------------------------------------------
class Spl_U3dBspl3Curve
{
public:
  Spl_U3dBspl3Curve (const Spl_U3dBspl3AbstractMatrix *bm, const std::vector<double> &ps);
  
  virtual ~Spl_U3dBspl3Curve();
  
  double get(double p) const;
  
protected:
  
  Spl_U3dSimpleBspl3CurveAr m_spline;
  
  bool m_is_interpolated;
  
private:
  
  bool ips2spline(const Spl_U3dBspl3AbstractMatrix *bm, const std::vector<double> &ps);
};

//------------------------------------------------------------------------------------------

class Spl_U3dBspl3Curve3d
{
public:
  Spl_U3dBspl3Curve3d(const Spl_U3dBspl3AbstractMatrix &bm, const Spl_U3dPointsAr &ps);
  
  ~Spl_U3dBspl3Curve3d();
  
  
  bool getPoint(double p, Spl_U3dPoint &point) const;
  
  bool getXYZp(double p, Spl_U3dPoint &point) const;
  
  bool getXYZpp(double p, Spl_U3dPoint &point) const;
  
  bool getCurvature(double p, double &c);
  
  int getControlPointsNum() const;
  
  bool getControlPoint(int ind, Spl_U3dPoint &point) const;
  
  
  bool add(Spl_U3dPoint &point, int ind = -1);
  
  bool remove(int ind);
  
  bool move(int ind);
  
  
  bool isReady() const
  {return m_is_ready;}
  
  
  bool getLineSegments(double ds, std::vector<double> &segs);
  
  bool getLineDeviations(double p1, double p2, double p, Spl_U3dPoint &dev);
  
  bool getMaxLineDeviation(double p1, double p2, double &md);
  
  
  void info();
  
protected:
  
  struct Bspl3dCurve 
  {
    Spl_U3dSimpleBspl3Curve *x, *y, *z;
  }; 
      
  std::vector<Bspl3dCurve> m_bspls;
  
  Spl_U3dBspl3AbstractMatrix *m_bm;
  
  Spl_U3dPointsAr m_ps;
  
  bool m_is_ready;
  
  void reconstruct();
  
  void clear();
  
  bool check_add_bps(const Spl_U3dPointsAr &ps);
  
  bool find_element(double p, int &nn, double &pp) const;
  
  bool get_line_distance(double p, Spl_U3dPoint kP, Spl_U3dPoint bP, Spl_U3dPoint &dist);
  
  bool get_line_distance(double pp, int nn, Spl_U3dPoint kP, Spl_U3dPoint bP, Spl_U3dPoint &dist);
  
  bool construct_isect_line(double p1, double p2, Spl_U3dPoint &kP, Spl_U3dPoint &bP);
  
  void devs_update(double p, int n,  Spl_U3dPoint &kP, Spl_U3dPoint &bP, 
           Spl_U3dPoint &dev, double &smx, double &smy, double &smz) {
    get_line_distance(p, n, kP, bP, dev);
    if (dev.getX() > smx)
      smx = dev.getX();
    if (dev.getY() > smy)
      smy = dev.getY();
    if (dev.getZ() > smz)
      smz = dev.getZ();
  }
  
  void next_segment (double p1, double p2, std::vector<double> &segs);
};

//------------------------------------------------------------------------------------------

class BaseSpline2D : public ParSpline
{
public:
  
  BaseSpline2D();
  
  BaseSpline2D(double *xx, double *zz, double *ddx, double *ddz, double *ccurv,bool *flag, int nn);
  
  virtual ~BaseSpline2D();
  
protected:
  
  virtual int IntervalSmoothParSpline(int N1,double *x1,double *z1,bool *flag,int N2, double *xf,double *zf,
                      double rc, double sp, double_ar *xx1, double_ar *zz1, double_ar *dx, double_ar *dz,
                      double_ar *curvity,int_ar *NumNodeBreak);
};

#endif /* base_splines.hpp */