Biomechanical Joint Model  Author: Anderson Maciel

---

vec.h

Go to the documentation of this file.

#ifndef AMAURY_VEC_H
#define AMAURY_VEC_H

#include <iostream>
#include <sstream>
#include <assert.h>
#include <math.h>


namespace LinAlg{

#define REAL double

 const double epsilon=1E-8;
 const double infinity=1E16;
 const int dimspace = 3;   // most vectors have dimension 3

class Vector 
{


public:
    typedef         double   value_type;
    typedef         double*  pointer;
    typedef         double*  iterator;
    typedef         double&  reference;
    typedef const   double*  const_iterator;
    typedef const   double&  const_reference;
  
    // constructors
    Vector();
    Vector(const Vector &A);
    Vector(unsigned N, const double& value = double(0));
    Vector(unsigned N, const double* v);
    Vector(const float *v, unsigned N);
    Vector(unsigned N, const char *s);

    // destructor
    virtual ~Vector();

    // assignments
    virtual Vector& operator=(const Vector &A);
    Vector& operator=(const double& scalar);
    Vector& operator+=(const Vector &A);
    Vector& operator-=(const Vector &A);
    Vector& operator*=(const double& scalar);
    Vector& operator/=(const double& scalar);

    // methods
    unsigned dim() const;
    double norm2() const;
    double norm() const;
    double normalize();

    bool operator >=(const double& scalar) const;
    bool operator <=(const double& scalar) const;
    bool operator >(const double& scalar) const;
    bool operator <(const double& scalar) const;
    
    // to build & manipulate vectors of another dimension
    Vector& insert(const Vector &A, int pos=0);
    Vector& add(const Vector &A, int pos=0);
    Vector extract(int pos=0, unsigned N=dimspace) const;
   
    // compatibility with other libraries that use floats only
    Vector& setVecf(const float *A);
    float *getVecf();
    operator float *();

    // access
    iterator begin() { return v_;}
    iterator end()   { return v_ + n_; }
    const_iterator begin() const { return v_;}
    const_iterator end() const  { return v_ + n_; }
    reference operator[](int i);
    const_reference operator[](int i) const;
  
    // binary operators
    friend Vector operator+(const Vector &A, const Vector &B);
    friend Vector operator-(const Vector &A, const Vector &B);
    friend Vector operator*(const Vector &A, const Vector &B);
    friend Vector operator*(const Vector &A, double s);
    friend Vector operator*(double s, const Vector &A);
    friend Vector operator/(const Vector &A, double s);
    friend bool operator ==(const Vector &A, const Vector &B);
    friend bool operator !=(const Vector &A, const Vector &B);

    friend bool operator >=(const Vector &A, const Vector &B);
    friend bool operator <=(const Vector &A, const Vector &B);
    friend bool operator >(const Vector &A, const Vector &B);
    friend bool operator <(const Vector &A, const Vector &B);

    friend double dot_prod(const Vector &A, const Vector &B);
    friend Vector cross_prod(const Vector &A, const Vector &B);

    friend double angle_between(const Vector &A, const Vector &B);

    // unary operators
    friend Vector& Sqrt(Vector &A);
    friend Vector& abs(Vector &A);
    friend Vector operator-(const Vector &A);

protected:
    double* v_;
    float *f_;
    unsigned n_, nf_;
    // internal helper function to create the array
    // of row pointers
    virtual void initialize(unsigned N);
    void copy(const double*  v);
    void set(const double& val);
    virtual void destroy(); 
    
    Vector& div(const Vector &A);
    Vector& div(const double& scalar);
};
 
  typedef Vector VEC;

  std::ostream& operator<<(std::ostream &s, const Vector &A);

  // friend function prototypes
  Vector operator+(const Vector &A, const Vector &B);
  Vector operator-(const Vector &A, const Vector &B);
  Vector operator*(const Vector &A, const Vector &B);
  Vector operator*(const Vector &A, double s);
  Vector operator*(double s, const Vector &A);
  Vector operator/(const Vector &A, double s);
  bool operator ==(const Vector &A, const Vector &B);
  bool operator !=(const Vector &A, const Vector &B);
  
  bool operator >=(const Vector &A, const Vector &B);
  bool operator <=(const Vector &A, const Vector &B);
  bool operator >(const Vector &A, const Vector &B);
  bool operator <(const Vector &A, const Vector &B);
  
  double dot_prod(const Vector &A, const Vector &B);
  Vector cross_prod(const Vector &A, const Vector &B);
  double angle_between(const Vector &A, const Vector &B);
  
  // unary operators
  Vector& Sqrt(Vector &A);
  Vector& abs(Vector &A);
  Vector operator-(const Vector &A);


 // private members


inline void Vector::copy(const double*  v)
{
  memcpy(v_, v, n_*sizeof(REAL));
}


inline void Vector::set(const double& val)
{
  unsigned N = n_;
  unsigned i;

  for (i=0; i< N; i++)
    v_[i] = val;        
}


inline Vector& Vector::div(const double& scalar)
{
  assert(scalar!=0);
  for (int i=0; i<n_; i++)
    v_[i] /= scalar;
  return *this;
}
 
inline Vector::~Vector()
{ 
   destroy(); 
}

// private members

inline void Vector::destroy()
{       
   delete [] v_;
   delete [] f_;
   f_ = NULL;
   v_ = NULL;
   n_ = nf_ = 0;
}


inline void Vector::initialize(unsigned N)
{
   v_ = new double[N];
   n_ = N;
}

// assignments

inline Vector& Vector::operator=(const Vector &A)
{
   if (v_ == A.v_)
      return *this;

   if (n_ == A.n_)         // no need to re-alloc
      copy(A.v_);
   else
   {
      destroy();
      initialize(A.n_);
      copy(A.v_);
   }

   return *this;
}

inline Vector operator*(const Vector &A, const Vector &B)
{
    unsigned N = A.n_;

    assert(N==B.n_);

    Vector tmp(N);
    unsigned i;

    for (i=0; i<N; i++)
            tmp.v_[i] = A.v_[i] * B.v_[i];

    return tmp;
}


/**************[ implementation of class Vector ]********************/

// constructors


inline Vector::Vector() : v_(0), n_(0), nf_(0), f_(0)
{}


inline Vector::Vector(const Vector &A) : v_(0), n_(0), nf_(0), f_(0)
{
  initialize(A.n_);
  copy(A.v_);
}


inline Vector::Vector(unsigned N, const double& value) :  v_(0),  n_(0), nf_(0), f_(0)
{
  initialize(N);
  set(value);
}


inline Vector::Vector(unsigned N, const double* v) :  v_(0), n_(0), nf_(0), f_(0)
{
  initialize(N);
  copy(v);
}
 

inline Vector::Vector(const  float *v, unsigned N) :  v_(0), n_(0), nf_(0), f_(0)
{
  initialize(N);
  for (int i=0; i< N; i++)
    v_[i] = v[i];
}
 

inline Vector::Vector(unsigned N, const char *s) :  v_(0),  n_(0), nf_(0), f_(0)
{
  initialize(N);
  std::istringstream ins(s);

  unsigned i;

  for (i=0; i<N; i++)
    ins >> v_[i];
}


// methods

// assignments
//

        

inline Vector& Vector::operator=(const double& scalar)
{ 
  set(scalar);  
  return *this;
}


inline Vector& Vector::operator+=(const Vector &A)
{
  assert(n_==A.n_);
  for (int i=0; i<n_; i++)
    v_[i] += A.v_[i];
  return *this;
}


inline Vector& Vector::operator-=(const Vector &A)
{
  assert(n_==A.n_);
  for (int i=0; i<n_; i++)
    v_[i] -= A.v_[i];
  return *this;
}


inline Vector& Vector::operator*=(const double& scalar)
{
  for (int i=0; i<n_; i++)
    v_[i] *= scalar;
  return *this;
}


inline Vector& Vector::operator/=(const double& scalar)
{
  return div(scalar);
}


inline unsigned Vector::dim() const 
{
  return  n_; 
}


inline double Vector::norm2() const 
{
  double sum = 0;

  for (int i=0; i<n_; i++)
    sum += v_[i] * v_[i];
  return sum; 
}


inline double Vector::norm() const 
{
  double sum = norm2();
  return sqrt(sum);
}


inline double Vector::normalize()
{
  double length = norm();
  if (length>epsilon) div(length);
  return length;
}


inline bool Vector::operator >= (const double& scalar) const
{
    for (int i=0; i<n_; i++)
      if (v_[i] < scalar) return false;
    return true;
}


inline bool Vector::operator <= (const double& scalar) const
{
    for (int i=0; i<n_; i++)
      if (v_[i] > scalar)  return false;
    return true;
}


inline bool Vector::operator > (const double& scalar) const     
{
    for (int i=0; i<n_; i++)
      if (v_[i] <= scalar) return false;
    return true;
}


inline bool Vector::operator < (const double& scalar) const      
{
    for (int i=0; i<n_; i++)
      if (v_[i] >= scalar)  return false;
    return true;
}
   

inline Vector& Vector::setVecf(const float *A)
{
  for (int i=0; i<n_; i++)
    v_[i] = A[i];
  return *this;
}


inline Vector::reference Vector::operator[](int i)
{ 
#ifdef BOUNDS_CHECK
  assert(0<=i);
  assert(i < n_);
#endif
  return v_[i]; 
}


inline Vector::const_reference Vector::operator[](int i) const
{
#ifdef BOUNDS_CHECK
  assert(0<=i);
  assert(i < n_);
#endif
  return v_[i]; 
}


inline Vector& Vector::insert(const Vector &A, int pos)
{
    unsigned N = A.n_;
#ifdef BOUNDS_CHECK
    assert(0<=pos);
    assert(N+pos <= n_);
#endif
    for (int i=0; i<N; i++)
      v_[i+pos] = A.v_[i];

    return *this;
}


inline Vector& Vector::add(const Vector &A, int pos)
{
    unsigned N = A.n_;
#ifdef BOUNDS_CHECK
    assert(0<=pos);
    assert(N+pos <= n_);
#endif
    for (int i=0; i<N; i++)
      v_[i+pos] += A.v_[i];

    return *this;
}

inline Vector Vector::extract(int pos, unsigned N) const
{
#ifdef BOUNDS_CHECK
    assert(0<=pos);
    assert( (N+pos) <= n_);
#endif
    Vector tmp(N);

    for (int i=0; i<N; i++)
      tmp.v_[i] = v_[i+pos];

    return tmp;
}

// to be compatible with other libraries

inline float* Vector::getVecf()
{   
  if (nf_ != n_)
    {
      delete [] f_;
      f_ = new float[n_];
    }
  
  for (int i=0; i<n_; i++) 
    f_[i] = v_[i];
        return f_;    
}


inline Vector::operator float *()
{ return getVecf(); }

/****************************  I/O  ********************************/


inline std::ostream& operator<<(std::ostream &s, const Vector &A)
{
    unsigned N=A.dim();

    s <<  "[ ";

    for (unsigned i=0; i<N; i++)
        s   << A[i] << " ";
    s << "]";

    return s;
}


inline Vector operator+(const Vector &A, const Vector &B)
{
  unsigned N = A.dim();  
  assert(N==B.dim());
  
  Vector tmp(N);
  unsigned i;
  
  for (i=0; i<N; i++)
      tmp.v_[i] = A.v_[i] + B.v_[i];
  
  return tmp;
}


inline Vector operator-(const Vector &A, const Vector &B)
{
    unsigned N = A.n_;

    assert(N==B.n_);

    Vector tmp(N);
    unsigned i;

    for (i=0; i<N; i++)
      tmp.v_[i] = A.v_[i] - B.v_[i];

    return tmp;
}



inline Vector operator-(const Vector &A)
{
    unsigned N = A.n_;

    Vector tmp(N);
    unsigned i;

    for (i=0; i<N; i++)
            tmp.v_[i] = -A.v_[i];

    return tmp;
}


inline Vector operator*(const Vector &A, double s)
{
    unsigned N = A.n_;

    Vector tmp(N);
    unsigned i;

    for (i=0; i<N; i++)
            tmp.v_[i] = A.v_[i] * s;

    return tmp;
}


inline Vector operator*(double s, const Vector &A)
{
    return A*s;
}


inline Vector operator/(const Vector &A, double s)
{
    unsigned N = A.n_;    
    Vector tmp(N);
    unsigned i;
    assert(s!=0);

    for (i=0; i<N; i++)
            tmp.v_[i] = A.v_[i] / s;

    return tmp;
}


inline double dot_prod(const Vector &A, const Vector &B)
{
    unsigned N = A.n_;
    assert(N == B.n_);

    unsigned i;
    double sum = 0;

    for (i=0; i<N; i++)
        sum += A.v_[i] * B.v_[i];

    return sum;
}


inline double angle_between(const Vector &A, const Vector &B)
{
    VEC U(A), V(B);
    U.normalize();
    V.normalize();

    return acosf( dot_prod(U,V) );
}

// defined for dim-3 vectors only

inline Vector cross_prod(const Vector &A, const Vector &B)
{
    assert(A.n_==3 && B.n_==3);
    Vector tmp(3);

    tmp[0] =   A[1]*B[2] - A[2]*B[1] ;
    tmp[1] = - A[0]*B[2] + A[2]*B[0] ;
    tmp[2] =   A[0]*B[1] - A[1]*B[0] ;

    return tmp;
}


inline bool operator ==(const Vector &A, const Vector &B)
{
    unsigned N = A.n_;
    if (N == B.n_)
    {
      unsigned i;
      for ( i=0; i<N; i++)
        if (A.v_[i] != B.v_[i]) break;
      return (i==N);
    }
    return false;
}


inline bool operator!=(const Vector &A, const  Vector &B)      
{
  return !(A==B);
}


inline bool operator >=(const Vector &A, const Vector &B)
{
    unsigned N = A.n_;
    if (N != B.n_) return false;

    for (int i=0; i<N; i++)
      if (A.v_[i] < B.v_[i]) return false;
    return true;
}


inline bool operator <=(const Vector &A, const Vector &B)
{
    unsigned N = A.n_;
    if (N != B.n_) return false;

    for (int i=0; i<N; i++)
      if (A.v_[i] > B.v_[i]) return false;
    return true;
}


inline bool operator >(const Vector &A, const Vector &B)      
{
    unsigned N = A.n_;
    if (N != B.n_) return false;

    for (int i=0; i<N; i++)
    {
       REAL delta = A.v_[i] - B.v_[i];
       if (delta <= 0) return false;
       if (fabs(delta)<1E-7) return false;
    }
    return true;
}


inline bool operator <(const Vector &A, const Vector &B)      
{
    unsigned N = A.n_;
    if (N != B.n_) return false;

    for (int i=0; i<N; i++)
    {
       REAL delta = A.v_[i] - B.v_[i];
       if (delta >= 0) return false;
       if (fabs(delta)<1E-7) return false;
    }
    return true;
}


inline Vector& Sqrt(Vector &A)
{  
  assert(A>=0);

  for (int i=0; i<A.n_; i++)
    A.v_[i] = sqrt(A.v_[i]);
  return A;
}

inline Vector& abs(Vector &A)
{  
  for (int i=0; i<A.n_; i++)
    A.v_[i] = fabs(A.v_[i]);
  return A;
}
   

}


#endif
// VECdoubleOR_H

---