nrutil_nr.h

Go to the documentation of this file.

#ifndef _NR_UTIL_H_
#define _NR_UTIL_H_

#include <string>
#include <cmath>
#include <complex>
#include <iostream>
#include <stdlib.h>
using namespace std;

typedef double DP;

template<class T>
inline const T SQR(const T a) {return a*a;}

template<class T>
inline const T SIGN(const T &a, const T &b) {return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);}

template<class T>
inline void SWAP(T &a, T &b) {T dum=a; a=b; b=dum;}

namespace NR {
        inline void nrerror(const string error_text)
        // Numerical Recipes standard error handler
        {
                cerr << "Numerical Recipes run-time error..." << endl;
                cerr << error_text << endl;
                cerr << "...now exiting to system..." << endl;
                exit(1);
        }
}

template <class T>
class NRVec {
private:
        int nn; // size of array. upper index is nn-1
        T *v;
public:
        NRVec();
        explicit NRVec(int n);          // Zero-based array
        NRVec(const T &a, int n);       //initialize to constant value
        NRVec(const T *a, int n);       // Initialize to array
        NRVec(const NRVec &rhs);        // Copy constructor
        NRVec & operator=(const NRVec &rhs);    //assignment
        NRVec & operator=(const T &a);  //assign a to every element
        inline T & operator[](const int i);     //i'th element
        inline const T & operator[](const int i) const;
        inline int size() const;
        ~NRVec();
};

template <class T>
NRVec<T>::NRVec() : nn(0), v(0) {}

template <class T>
NRVec<T>::NRVec(int n) : nn(n), v(new T[n]) {}

template <class T>
NRVec<T>::NRVec(const T& a, int n) : nn(n), v(new T[n])
{
        for(int i=0; i<n; i++)
                v[i] = a;
}

template <class T>
NRVec<T>::NRVec(const T *a, int n) : nn(n), v(new T[n])
{
        for(int i=0; i<n; i++)
                v[i] = *a++;
}

template <class T>
NRVec<T>::NRVec(const NRVec<T> &rhs) : nn(rhs.nn), v(new T[nn])
{
        for(int i=0; i<nn; i++)
                v[i] = rhs[i];
}

template <class T>
NRVec<T> & NRVec<T>::operator=(const NRVec<T> &rhs)
// postcondition: normal assignment via copying has been performed;
//              if vector and rhs were different sizes, vector
//              has been resized to match the size of rhs
{
        if (this != &rhs)
        {
                if (nn != rhs.nn) {
                        if (v != 0) delete [] (v);
                        nn=rhs.nn;
                        v= new T[nn];
                }
                for (int i=0; i<nn; i++)
                        v[i]=rhs[i];
        }
        return *this;
}

template <class T>
NRVec<T> & NRVec<T>::operator=(const T &a)      //assign a to every element
{
        for (int i=0; i<nn; i++)
                v[i]=a;
        return *this;
}

template <class T>
inline T & NRVec<T>::operator[](const int i)    //subscripting
{
        return v[i];
}

template <class T>
inline const T & NRVec<T>::operator[](const int i) const        //subscripting
{
        return v[i];
}

template <class T>
inline int NRVec<T>::size() const
{
        return nn;
}

template <class T>
NRVec<T>::~NRVec()
{
        if (v != 0)
                delete[] (v);
}

template <class T>
class NRMat {
private:
        int nn;
        int mm;
        T **v;
public:
        NRMat();
        NRMat(int n, int m);                    // Zero-based array
        NRMat(const T &a, int n, int m);        //Initialize to constant
        NRMat(const T *a, int n, int m);        // Initialize to array
        NRMat(const NRMat &rhs);                // Copy constructor
        NRMat & operator=(const NRMat &rhs);    //assignment
        NRMat & operator=(const T &a);          //assign a to every element
        inline T* operator[](const int i);      //subscripting: pointer to row i
        inline const T* operator[](const int i) const;
        inline int nrows() const;
        inline int ncols() const;
        ~NRMat();
};

template <class T>
NRMat<T>::NRMat() : nn(0), mm(0), v(0) {}

template <class T>
NRMat<T>::NRMat(int n, int m) : nn(n), mm(m), v(new T*[n])
{
        v[0] = new T[m*n];
        for (int i=1; i< n; i++)
                v[i] = v[i-1] + m;
}

template <class T>
NRMat<T>::NRMat(const T &a, int n, int m) : nn(n), mm(m), v(new T*[n])
{
        int i,j;
        v[0] = new T[m*n];
        for (i=1; i< n; i++)
                v[i] = v[i-1] + m;
        for (i=0; i< n; i++)
                for (j=0; j<m; j++)
                        v[i][j] = a;
}

template <class T>
NRMat<T>::NRMat(const T *a, int n, int m) : nn(n), mm(m), v(new T*[n])
{
        int i,j;
        v[0] = new T[m*n];
        for (i=1; i< n; i++)
                v[i] = v[i-1] + m;
        for (i=0; i< n; i++)
                for (j=0; j<m; j++)
                        v[i][j] = *a++;
}

template <class T>
NRMat<T>::NRMat(const NRMat &rhs) : nn(rhs.nn), mm(rhs.mm), v(new T*[nn])
{
        int i,j;
        v[0] = new T[mm*nn];
        for (i=1; i< nn; i++)
                v[i] = v[i-1] + mm;
        for (i=0; i< nn; i++)
                for (j=0; j<mm; j++)
                        v[i][j] = rhs[i][j];
}

template <class T>
NRMat<T> & NRMat<T>::operator=(const NRMat<T> &rhs)
// postcondition: normal assignment via copying has been performed;
//              if matrix and rhs were different sizes, matrix
//              has been resized to match the size of rhs
{
        if (this != &rhs) {
                int i,j;
                if (nn != rhs.nn || mm != rhs.mm) {
                        if (v != 0) {
                                delete[] (v[0]);
                                delete[] (v);
                        }
                        nn=rhs.nn;
                        mm=rhs.mm;
                        v = new T*[nn];
                        v[0] = new T[mm*nn];
                }
                for (i=1; i< nn; i++)
                        v[i] = v[i-1] + mm;
                for (i=0; i< nn; i++)
                        for (j=0; j<mm; j++)
                                v[i][j] = rhs[i][j];
        }
        return *this;
}

template <class T>
NRMat<T> & NRMat<T>::operator=(const T &a)      //assign a to every element
{
        for (int i=0; i< nn; i++)
                for (int j=0; j<mm; j++)
                        v[i][j] = a;
        return *this;
}

template <class T>
inline T* NRMat<T>::operator[](const int i)     //subscripting: pointer to row i
{
        return v[i];
}

template <class T>
inline const T* NRMat<T>::operator[](const int i) const
{
        return v[i];
}

template <class T>
inline int NRMat<T>::nrows() const
{
        return nn;
}

template <class T>
inline int NRMat<T>::ncols() const
{
        return mm;
}

template <class T>
NRMat<T>::~NRMat()
{
        if (v != 0) {
                delete[] (v[0]);
                delete[] (v);
        }
}

template <class T>
class NRMat3d {
private:
        int nn;
        int mm;
        int kk;
        T ***v;
public:
        NRMat3d();
        NRMat3d(int n, int m, int k);
        inline T** operator[](const int i);     //subscripting: pointer to row i
        inline const T* const * operator[](const int i) const;
        inline int dim1() const;
        inline int dim2() const;
        inline int dim3() const;
        ~NRMat3d();
};

template <class T>
NRMat3d<T>::NRMat3d(): nn(0), mm(0), kk(0), v(0) {}

template <class T>
NRMat3d<T>::NRMat3d(int n, int m, int k) : nn(n), mm(m), kk(k), v(new T**[n])
{
        int i,j;
        v[0] = new T*[n*m];
        v[0][0] = new T[n*m*k];
        for(j=1; j<m; j++)
                v[0][j] = v[0][j-1] + k;
        for(i=1; i<n; i++) {
                v[i] = v[i-1] + m;
                v[i][0] = v[i-1][0] + m*k;
                for(j=1; j<m; j++)
                        v[i][j] = v[i][j-1] + k;
        }
}

template <class T>
inline T** NRMat3d<T>::operator[](const int i) //subscripting: pointer to row i
{
        return v[i];
}

template <class T>
inline const T* const * NRMat3d<T>::operator[](const int i) const
{
        return v[i];
}

template <class T>
inline int NRMat3d<T>::dim1() const
{
        return nn;
}

template <class T>
inline int NRMat3d<T>::dim2() const
{
        return mm;
}

template <class T>
inline int NRMat3d<T>::dim3() const
{
        return kk;
}

template <class T>
NRMat3d<T>::~NRMat3d()
{
        if (v != 0) {
                delete[] (v[0][0]);
                delete[] (v[0]);
                delete[] (v);
        }
}

//The next 3 classes are used in artihmetic coding, Huffman coding, and
//wavelet transforms respectively. This is as good a place as any to put them!

class arithcode {
private:
        NRVec<unsigned long> *ilob_p,*iupb_p,*ncumfq_p;
public:
        NRVec<unsigned long> &ilob,&iupb,&ncumfq;
        unsigned long jdif,nc,minint,nch,ncum,nrad;
        arithcode(unsigned long n1, unsigned long n2, unsigned long n3)
                : ilob_p(new NRVec<unsigned long>(n1)),
                iupb_p(new NRVec<unsigned long>(n2)),
                ncumfq_p(new NRVec<unsigned long>(n3)),
                ilob(*ilob_p),iupb(*iupb_p),ncumfq(*ncumfq_p) {}
        ~arithcode() {
                if (ilob_p != 0) delete ilob_p;
                if (iupb_p != 0) delete iupb_p;
                if (ncumfq_p != 0) delete ncumfq_p;
        }
};

class huffcode {
private:
        NRVec<unsigned long> *icod_p,*ncod_p,*left_p,*right_p;
public:
        NRVec<unsigned long> &icod,&ncod,&left,&right;
        int nch,nodemax;
        huffcode(unsigned long n1, unsigned long n2, unsigned long n3,
                unsigned long n4) :
                icod_p(new NRVec<unsigned long>(n1)),
                ncod_p(new NRVec<unsigned long>(n2)),
                left_p(new NRVec<unsigned long>(n3)),
                right_p(new NRVec<unsigned long>(n4)),
                icod(*icod_p),ncod(*ncod_p),left(*left_p),right(*right_p) {}
        ~huffcode() {
                if (icod_p != 0) delete icod_p;
                if (ncod_p != 0) delete ncod_p;
                if (left_p != 0) delete left_p;
                if (right_p != 0) delete right_p;
        }
};

class wavefilt {
private:
        NRVec<DP> *cc_p,*cr_p;
public:
        int ncof,ioff,joff;
        NRVec<DP> &cc,&cr;
        wavefilt() : cc(*cc_p),cr(*cr_p) {}
        wavefilt(const DP *a, const int n) :  //initialize to array
                cc_p(new NRVec<DP>(n)),cr_p(new NRVec<DP>(n)),
                ncof(n),ioff(-(n >> 1)),joff(-(n >> 1)),cc(*cc_p),cr(*cr_p) {
                        int i;
                        for (i=0; i<n; i++)
                                cc[i] = *a++;
                        DP sig = -1.0;
                        for (i=0; i<n; i++) {
                                cr[n-1-i]=sig*cc[i];
                                sig = -sig;
                        }
        }
        ~wavefilt() {
                if (cc_p != 0) delete cc_p;
                if (cr_p != 0) delete cr_p;
        }
};

//Overloaded complex operations to handle mixed float and double
//This takes care of e.g. 1.0/z, z complex<float>

inline const complex<float> operator+(const double &a,
        const complex<float> &b) { return float(a)+b; }

inline const complex<float> operator+(const complex<float> &a,
        const double &b) { return a+float(b); }

inline const complex<float> operator-(const double &a,
        const complex<float> &b) { return float(a)-b; }

inline const complex<float> operator-(const complex<float> &a,
        const double &b) { return a-float(b); }

inline const complex<float> operator*(const double &a,
        const complex<float> &b) { return float(a)*b; }

inline const complex<float> operator*(const complex<float> &a,
        const double &b) { return a*float(b); }

inline const complex<float> operator/(const double &a,
        const complex<float> &b) { return float(a)/b; }

inline const complex<float> operator/(const complex<float> &a,
        const double &b) { return a/float(b); }

//some compilers choke on pow(float,double) in single precision. also atan2

inline float pow (float x, double y) {return pow(double(x),y);}
inline float pow (double x, float y) {return pow(x,double(y));}
inline float atan2 (float x, double y) {return atan2(double(x),y);}
inline float atan2 (double x, float y) {return atan2(x,double(y));}
#endif /* _NR_UTIL_H_ */