util.h

Go to the documentation of this file.

#ifndef UTIL_H
#define UTIL_H

#include <iostream>
#include <list>
#include <vector>
#include <cassert>
#include <cmath>
#include <string>
#include <string.h>

#include "persistance-noop.h"

#define bzero(b,len) (memset((b), '\0', (len)), (void) 0)

/* These are values and weights for gausian quadrature */
extern const double ALPHA[7];
extern const double BETA[7];
extern const double WEIGHT[7];

typedef double Radians;
extern const double PI;
extern const double TWOPI;
extern const double HALFPI;

namespace Tumble {
  inline bool is_zero(double a) {
    return fabs(a) < 1e-15;
  }
  inline bool are_equal(double a, double b) {
    return is_zero(a - b);
  }
};


class Point2D {
public:
    double coords[ 2 ]; 
    /* constructors */
    Point2D()
    {
        coords[ 0 ] = 0.0;
        coords[ 1 ] = 0.0;
    }

    Point2D( double x, double y )
    {
        coords[ 0 ] = x;
        coords[ 1 ] = y;
    }

    Point2D( const Point2D &o)
    {
        coords[ 0 ] = o.coords[ 0 ];
        coords[ 1 ] = o.coords[ 1 ];
    }

    /* accessors and setters */
    Point2D& assign( double x, double y )
    {
        coords[ 0 ] = x;
        coords[ 1 ] = y;
        return *this;
    }

    double operator[] ( int i ) const { return coords[ i ]; }
    double x() const { return coords[ 0 ]; }
    double y() const { return coords[ 1 ]; }

    bool operator<( const Point2D &o ) const
    {
        return ( magsq() < o.magsq() );
    }

    bool operator<=( const Point2D &o ) const
    {
        return ( magsq() <= o.magsq() );
    }

    bool operator>( const Point2D &o ) const
    {
        return ( magsq() > o.magsq() );
    }

    bool operator>=( const Point2D &o ) const
    {
        return ( magsq() >= o.magsq() );
    }

    Point2D& operator= ( const Point2D &o )
    {
        coords[ 0 ] = o.coords[ 0 ];
        coords[ 1 ] = o.coords[ 1 ];
        return *this;
    }

    /* in place operations (perfer these)*/
    Point2D& operator+=( const Point2D &o )
    {
        coords[ 0 ] += o.coords[ 0 ];
        coords[ 1 ] += o.coords[ 1 ];
        return *this;
    }

    Point2D& operator-=( const Point2D &o )
    {
        coords[ 0 ] -= o.coords[ 0 ];
        coords[ 1 ] -= o.coords[ 1 ];
        return *this;
    }

    Point2D& operator*=( double s )
    {
        coords[ 0 ] *= s;
        coords[ 1 ] *= s;
        return *this;
    }

    Point2D& operator/=( double s )
    {
        coords[ 0 ] /= s;
        coords[ 1 ] /= s;
        return *this;
    }

    /* temporary generating operations (avoid these) */
    Point2D operator+ ( const Point2D &o ) const
    {
        return Point2D(coords[ 0 ] + o.coords[ 0 ], coords[ 1 ] + o.coords[ 1 ]);
    }

    Point2D operator- ( const Point2D &o ) const
    {
        return Point2D(coords[ 0 ] - o.coords[ 0 ], coords[ 1 ] - o.coords[ 1 ]);
    }

    Point2D operator* ( double s ) const
    {
        return Point2D(coords[ 0 ] *s, coords[ 1 ] * s);
    }

    Point2D operator/ ( double s ) const
    {
        return Point2D(coords[ 0 ] / s, coords[ 1 ] / s);
    }

    /* vector arithmetic operations */
    double dot( const Point2D &o ) const
    {
        return coords[ 0 ] * o.coords[ 0 ] + coords[ 1 ] * o.coords[ 1 ];
    }

    double magsq() const
    {
        return coords[ 0 ] * coords[ 0 ] + coords[ 1 ] * coords[ 1 ];
    }

    double mag() const
    {
        return sqrt( coords[ 0 ] * coords[ 0 ] + coords[ 1 ] * coords[ 1 ] );
    }

    double cross( const Point2D &a ) const
    {
        return coords[ 0 ] * a.coords[ 1 ] - coords[ 1 ] * a.coords[ 0 ];
    }

    /* geometry tests */

    bool machine_equal(const Point2D& other) const {
      return x() == other.x() && y() == other.y();
    }

    bool is_left_of( const Point2D &a, const Point2D &b) const { return line_side_test(a,b)>0.0; }

    bool is_right_of( const Point2D &a, const Point2D &b) const { return line_side_test(a,b)<0.0; }

    double line_side_test( Point2D a, Point2D b) const;

    bool in_circle_test( Point2D a, Point2D b, Point2D c) const;

    double dist_from_line(const Point2D &a, const Point2D &b ) const;

    double angle(const Point2D &a, const Point2D &c) const
    {
        Point2D ab = a - *this;
        Point2D cb = c - *this;
    return (ab.dot(cb) / (ab.mag() * cb.mag()));
    }

    /* information retrieval */
    void print(){ std::cout<<*this;}
    friend std::ostream &operator<<( std::ostream&, const Point2D &p );
};

class LinearData {
public:
    double *fields; 
    int length; 
    /* constructors */
    LinearData(); 
    LinearData(int len); 
    LinearData( double *data, int len); 
    LinearData( const std::vector<double> &data); 
    LinearData( const LinearData &ld); 
    ~LinearData(); 
    /* data modification */
    int update( const std::vector<double> &data ); 
    int update( double *data, int len ); 
    void resize( int len); 
    void zero();  
    void set(double d, int i) {
        assert(i>=0 && i < length);
        fields[i] = d;
    }

    double operator[]( int i ) const {
        assert(i>=0 && i<length);
        return fields[ i ];
    }

    double get(int i) const {
        assert(i>=0 && i < length);
        return fields[i];
    }

    LinearData& operator= ( const LinearData& d);

    /* in place operations (stongly prefer these) */
    LinearData& operator+= ( const LinearData& d);
    LinearData& operator-= ( const LinearData& d);
    LinearData& operator*= ( double s);
    LinearData& operator/= ( double s);

    /* temporary creation operations (stongly avoid these) */
    LinearData operator+ ( const LinearData& d) const;
    LinearData operator- ( const LinearData& d) const;
    LinearData operator* ( double s) const;
    LinearData operator/ ( double s) const;

    /* information retrieval */
    void print(){ std::cout<<*this;}
    friend std::ostream &operator<<( std::ostream&, const LinearData& );
};

typedef PersistantList<Point2D>::iterator ControlPoint;
std::ostream& operator << (std::ostream& os, const ControlPoint& cp);

typedef PersistantList<LinearData>::iterator DataPoint;
std::ostream& operator << (std::ostream& os, const DataPoint& cp);

class Vec3D {
    public:
        double coords[ 3 ]; 
        /* constructors */
        Vec3D()
        {
      coords[0] = 0;
      coords[1] = 0;
      coords[2] = 0;
        }

        Vec3D( double x, double y, double z )
        {
            coords[ 0 ] = x;
            coords[ 1 ] = y;
            coords[ 2 ] = z;
        }

        Vec3D( const Point2D &p, double z=0.0)
        {
            coords[ 0 ] = p.x();
            coords[ 1 ] = p.y();
            coords[ 2 ] = z;
        }

        Vec3D( const Vec3D &o)
        {
      coords[0] = o.coords[0];
      coords[1] = o.coords[1];
      coords[2] = o.coords[2];
        }

        /* accessors and setters */
        Vec3D& assign( double x, double y, double z )
        {
            coords[ 0 ] = x;
            coords[ 1 ] = y;
            coords[ 2 ] = z;
            return *this;
        }

        Vec3D& assign( const Point2D &p, double z )
        {
            coords[ 0 ] = p.x();
            coords[ 1 ] = p.y();
            coords[ 2 ] = z;
            return *this;
        }

        double operator[] ( int i ) const { return coords[ i ]; }
        double x() const { return coords[ 0 ]; }
        double y() const { return coords[ 1 ]; }
        double z() const { return coords[ 2 ]; }

        bool operator<( const Vec3D &o ) const
        {
            return ( magsq() < o.magsq() );
        }

        bool operator<=( const Vec3D &o ) const
        {
            return ( magsq() <= o.magsq() );
        }

        bool operator>( const Vec3D &o ) const
        {
            return ( magsq() > o.magsq() );
        }

        bool operator>=( const Vec3D &o ) const
        {
            return ( magsq() >= o.magsq() );
        }

        /* assignment */
        Vec3D& operator= ( const Vec3D &o )
        {
      coords [0] = o.coords[0];
      coords [1] = o.coords[1];
      coords [2] = o.coords[2];
      return *this;
        }

        /* in place operations (perfer these)*/
        Vec3D& operator+=( const Vec3D &o )
        {
            coords[ 0 ] += o.coords[ 0 ];
            coords[ 1 ] += o.coords[ 1 ];
            coords[ 2 ] += o.coords[ 2 ];
            return *this;
        }

        Vec3D& operator-=( const Vec3D &o )
        {
            coords[ 0 ] -= o.coords[ 0 ];
            coords[ 1 ] -= o.coords[ 1 ];
            coords[ 2 ] -= o.coords[ 2 ];
            return *this;
        }

        Vec3D& operator*=( double s )
        {
            coords[ 0 ] *= s;
            coords[ 1 ] *= s;
            coords[ 2 ] *= s;
            return *this;
        }

        Vec3D& operator/=( double s )
        {
            coords[ 0 ] /= s;
            coords[ 1 ] /= s;
            coords[ 2 ] /= s;
            return *this;
        }

        /* temporary generating operations (avoid these) */
        Vec3D operator+ ( const Vec3D &o ) const
        {
            return Vec3D(coords[ 0 ] + o.coords[ 0 ], coords[ 1 ] + o.coords[ 1 ], coords[ 2 ] + o.coords[ 2 ]);
        }

        Vec3D operator- ( const Vec3D &o ) const
        {
            return Vec3D(coords[ 0 ] - o.coords[ 0 ], coords[ 1 ] - o.coords[ 1 ], coords[ 2 ] - o.coords[ 2 ]);
        }

        Vec3D operator* ( double s ) const
        {
            return Vec3D(coords[ 0 ] *s, coords[ 1 ] * s, coords[ 2 ] * s);
        }

        Vec3D operator/ ( double s ) const
        {
            return Vec3D(coords[ 0 ] / s, coords[ 1 ] / s, coords[ 2 ] / s);
        }

        /* vector arithmetic operations */
        double dot( const Vec3D &o ) const
        {
            return coords[ 0 ] * o.coords[ 0 ] + coords[ 1 ] * o.coords[ 1 ] + coords[ 2 ] * o.coords[ 2 ];
        }

        /* vector arithmetic operations */
        double operator*( const Vec3D &o ) const
        {
            return coords[ 0 ] * o.coords[ 0 ] + coords[ 1 ] * o.coords[ 1 ] + coords[ 2 ] * o.coords[ 2 ];
        }

        double magsq() const
        {
            return coords[ 0 ] * coords[ 0 ] + coords[ 1 ] * coords[ 1 ] + coords[ 2 ] * coords[ 2 ];
        }

        double mag() const
        {
            return sqrt( coords[ 0 ] * coords[ 0 ] + coords[ 1 ] * coords[ 1 ] + coords[ 2 ] * coords[ 2 ]);
        }

        double norm() const { return mag(); }

        Vec3D cross( const Vec3D &o ) const
        {
            return Vec3D( coords[1]*o.coords[2]-coords[2]*o.coords[1],
                          coords[0]*o.coords[2]-coords[2]*o.coords[0],
                          coords[0]*o.coords[1]-coords[1]*o.coords[0]);
        }

        double angle(const Vec3D &o)
        {
           double n = norm();
           double on = o.norm();
           if ((0.0==n) || (0.0==on))  return PI;
           double costheta = dot(o) / (n*on);
           if (costheta <= -1.0) return PI;
           if (costheta >= 1.0 ) return 0;
           return acos(costheta);
        }

        /* information retrieval */
        void print(){ std::cout<<*this;}
        friend std::ostream &operator<<( std::ostream&, const Vec3D &p );
};

class Matrix {
private:
    /* This allocates a new Matrix which you must delete */
    Matrix* submat( int i, int j );

    double matminor( int i, int j );
    double cofactor( int i, int j );
    int lu_decompose(int *P );
public:
    double** m;
    int rows, cols;

    /* Constructors */
    Matrix( int, int );
    Matrix( const Matrix& );
    ~Matrix();

    /* Accessors */
    double operator() (const int i,const int j) const {
        return m[i][j];
    }
    void set(int i, int j, double val){
        if(i>=0 && j>=0 && i<rows && j<cols) m[i][j]=val;
    }

    /* detrminate */
    double det();

    /* ***********************************************
     * WARNING: these functions allocate a matrix or vector
     * and return it.  You must be careful to delete it before
     * throwing away the reference
     */
    Matrix* transpose_times_self();
    Matrix* inverse();
    double* transpose_times_vector( double* );
    double* times_vector( double* );
    /************************************************/

    void print(){ std::cout<<*this<<std::endl; }
    friend std::ostream& operator<<( std::ostream&, Matrix );
};


class Matrix5 {
private:
    void factor();
public:
    int size;
    double *d1, *d2, *d3, *d4, *d5;
    bool factored;
    Matrix5( const Matrix& );
    Matrix5( const Matrix5& );
    Matrix5( int );
    ~Matrix5();
    double* LUsolve( double* );
   // friend std::ostream& operator<<( std::ostream&, Matrix5 );
};




template <class T>
class PQueueEntry {
public:
    T i;
    double p;
    PQueueEntry *prev, *next;
    PQueueEntry( const T &_i, double _p ) : i( _i ), p( _p ), prev( NULL ), next( NULL ) {}
};


template <class T>
class PQueue {
public:
    int sz;
    PQueueEntry<T> *front, *back;

    PQueue() {
        sz=0;
        front = NULL;
        back = NULL;
    }
    int size() {
        return sz;
    }
    bool pop( T &val );
    void push( const T &val, double p );
    bool remove( const T &val );
    void print(){
        std::cout<<"PQueue: [";
        PQueueEntry<T> *temp=front;
        while(temp && temp->next!=NULL){
            std::cout<<"{"<<temp->i<<", "<<temp->p<<"}, ";
            temp=temp->next;
        }
        if(temp!=NULL) std::cout<<"{"<<temp->i<<", "<<temp->p<<"}";
        std::cout<<"]"<<std::endl;
    }
};

template<class T>
bool PQueue<T>::pop( T &val ) {
    PQueueEntry<T> * e = front;
    if ( e == NULL ) return false;
    val = e->i;
    front = e->next;
    if ( front == NULL ) back = NULL;
    else front->prev = NULL;
    sz--;
    delete e;
    return true;
}

template<class T>
void PQueue<T>::push( const T &val, double p ) {
    PQueueEntry<T> * n = new PQueueEntry<T>( val, p );
    PQueueEntry<T> * temp;
    sz++;
    if ( front == NULL ) {
        front = n;
        back = n;
    } else if ( p <= front->p ) {
        n->next = front;
        front->prev = n;
        front = n;
    } else {
        temp = front;
        while ( temp && p > temp->p ) {temp = temp->next;}
        n->next = temp;
        if ( temp ){
            assert(temp != front);
            n->prev = temp->prev;
            temp->prev = n;
        }
        else{
            assert(back != NULL);
            n->prev = back;
            back = n;
        }
        assert(n->prev != NULL);
        n->prev->next = n;
    }
}

template<class T>
bool PQueue<T>::remove( const T &val ) {
    PQueueEntry<T> * temp=front;
    while ( temp && temp->i != val ) temp = temp->next;
    if ( temp == NULL ) return false;
    sz--;
    if ( front == temp ) {
        front = temp->next;
        if ( front != NULL ) front->prev = NULL;
    } else {
        temp->prev->next = temp->next;
    }

    if ( back == temp ) {
        back = temp->prev;
        if ( back != NULL ) back->next = NULL;
    } else {
        temp->next->prev = temp->prev;
    }

    delete temp;
    return true;
}

int add_filename_extension(char *filename, const char *ext, char *buf, int buf_len);

bool poly_convex(Point2D poly[], int size);

bool point_in_poly_kernel(Point2D poly[], int size, const Point2D &point);

#endif /* UTIL_H */

---