- Main Page
- Namespaces
- Classes
- Files

- File List
- File Members

# C:/AngioLocalizeNoAuto/MevisLabNoAutoReg/STUMevisApp/Angiolocalize/Sources/ML/STUChamferMatching/nr3.h

    00001 #ifndef _NR3_H_
    00002 #define _NR3_H_
    00003
    00004 //#define _CHECKBOUNDS_ 1
    00005 //#define _USESTDVECTOR_ 1
    00006 //#define _USENRERRORCLASS_ 1
    00007 //#define _TURNONFPES_ 1
    00008
    00009 // all the system #include's we'll ever need
    00010 #include <fstream>
    00011 #include <cmath>
    00012 #include <complex>
    00013 #include <iostream>
    00014 #include <iomanip>
    00015 #include <vector>
    00016 #include <limits>
    00017 #include <stdlib.h>
    00018 #include <stdio.h>
    00019 #include <time.h>
    00020 #include <fcntl.h>
    00021 #include <string.h>
    00022 #include <ctype.h>
    00023 #include <limits>
    00024
    00025 using namespace std;
    00026
    00027 // macro-like inline functions
    00028
    00029 template<class T>
    00030 inline T SQR(const T a) {return a*a;}
    00031
    00032 template<class T>
    00033 inline const T &MAX(const T &a, const T &b)
    00034         {return b > a ? (b) : (a);}
    00035
    00036 inline float MAX(const double &a, const float &b)
    00037         {return b > a ? (b) : float(a);}
    00038
    00039 inline float MAX(const float &a, const double &b)
    00040         {return b > a ? float(b) : (a);}
    00041
    00042 template<class T>
    00043 inline const T &MIN(const T &a, const T &b)
    00044         {return b < a ? (b) : (a);}
    00045
    00046 inline float MIN(const double &a, const float &b)
    00047         {return b < a ? (b) : float(a);}
    00048
    00049 inline float MIN(const float &a, const double &b)
    00050         {return b < a ? float(b) : (a);}
    00051
    00052 template<class T>
    00053 inline T SIGN(const T &a, const T &b)
    00054         {return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);}
    00055
    00056 inline float SIGN(const float &a, const double &b)
    00057         {return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);}
    00058
    00059 inline float SIGN(const double &a, const float &b)
    00060         {return (float)(b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a));}
    00061
    00062 template<class T>
    00063 inline void SWAP(T &a, T &b)
    00064         {T dum=a; a=b; b=dum;}
    00065
    00066 // exception handling
    00067
    00068 #ifndef _USENRERRORCLASS_
    00069 #define throw(message) \
    00070 {printf("ERROR: %s\n     in file %s at line %d\n", message,__FILE__,__LINE__); throw(1);}
    00071 #else
    00072 struct NRerror {
    00073         char *message;
    00074         char *file;
    00075         int line;
    00076         NRerror(char *m, char *f, int l) : message(m), file(f), line(l) {}
    00077 };
    00078 #define throw(message) throw(NRerror(message,__FILE__,__LINE__));
    00079 void NRcatch(NRerror err) {
    00080         printf("ERROR: %s\n     in file %s at line %d\n",
    00081                 err.message, err.file, err.line);
    00082         exit(1);
    00083 }
    00084 #endif
    00085
    00086 // usage example:
    00087 //
    00088 //      try {
    00089 //              somebadroutine();
    00090 //      }
    00091 //      catch(NRerror s) {NRcatch(s);}
    00092 //
    00093 // (You can of course substitute any other catch body for NRcatch(s).)
    00094
    00095
    00096 // Vector and Matrix Classes
    00097
    00098 #ifdef _USESTDVECTOR_
    00099 #define NRvector vector
    00100 #else
    00101
    00102 template <class T>
    00103 class NRvector {
    00104 private:
    00105         int nn; // size of array. upper index is nn-1
    00106         T *v;
    00107 public:
    00108         NRvector();
    00109         explicit NRvector(int n);               // Zero-based array
    00110         NRvector(int n, const T &a);    //initialize to constant value
    00111         NRvector(int n, const T *a);    // Initialize to array
    00112         NRvector(const NRvector &rhs);  // Copy constructor
    00113         NRvector & operator=(const NRvector &rhs);      //assignment
    00114         typedef T value_type; // make T available externally
    00115         inline T & operator[](const int i);     //i'th element
    00116         inline const T & operator[](const int i) const;
    00117         inline int size() const;
    00118         void resize(int newn); // resize (contents not preserved)
    00119         void assign(int newn, const T &a); // resize and assign a constant value
    00120         ~NRvector();
    00121 };
    00122
    00123 // NRvector definitions
    00124
    00125 template <class T>
    00126 NRvector<T>::NRvector() : nn(0), v(NULL) {}
    00127
    00128 template <class T>
    00129 NRvector<T>::NRvector(int n) : nn(n), v(n>0 ? new T[n] : NULL) {}
    00130
    00131 template <class T>
    00132 NRvector<T>::NRvector(int n, const T& a) : nn(n), v(n>0 ? new T[n] : NULL)
    00133 {
    00134         for(int i=0; i<n; i++) v[i] = a;
    00135 }
    00136
    00137 template <class T>
    00138 NRvector<T>::NRvector(int n, const T *a) : nn(n), v(n>0 ? new T[n] : NULL)
    00139 {
    00140         for(int i=0; i<n; i++) v[i] = *a++;
    00141 }
    00142
    00143 template <class T>
    00144 NRvector<T>::NRvector(const NRvector<T> &rhs) : nn(rhs.nn), v(nn>0 ? new T[nn] : NULL)
    00145 {
    00146         for(int i=0; i<nn; i++) v[i] = rhs[i];
    00147 }
    00148
    00149 template <class T>
    00150 NRvector<T> & NRvector<T>::operator=(const NRvector<T> &rhs)
    00151 // postcondition: normal assignment via copying has been performed;
    00152 //              if vector and rhs were different sizes, vector
    00153 //              has been resized to match the size of rhs
    00154 {
    00155         if (this != &rhs)
    00156         {
    00157                 if (nn != rhs.nn) {
    00158                         if (v != NULL) delete [] (v);
    00159                         nn=rhs.nn;
    00160                         v= nn>0 ? new T[nn] : NULL;
    00161                 }
    00162                 for (int i=0; i<nn; i++)
    00163                         v[i]=rhs[i];
    00164         }
    00165         return *this;
    00166 }
    00167
    00168 template <class T>
    00169 inline T & NRvector<T>::operator[](const int i) //subscripting
    00170 {
    00171 #ifdef _CHECKBOUNDS_
    00172 if (i<0 || i>=nn) {
    00173         throw("NRvector subscript out of bounds");
    00174 }
    00175 #endif
    00176         return v[i];
    00177 }
    00178
    00179 template <class T>
    00180 inline const T & NRvector<T>::operator[](const int i) const     //subscripting
    00181 {
    00182 #ifdef _CHECKBOUNDS_
    00183 if (i<0 || i>=nn) {
    00184         throw("NRvector subscript out of bounds");
    00185 }
    00186 #endif
    00187         return v[i];
    00188 }
    00189
    00190 template <class T>
    00191 inline int NRvector<T>::size() const
    00192 {
    00193         return nn;
    00194 }
    00195
    00196 template <class T>
    00197 void NRvector<T>::resize(int newn)
    00198 {
    00199         if (newn != nn) {
    00200                 if (v != NULL) delete[] (v);
    00201                 nn = newn;
    00202                 v = nn > 0 ? new T[nn] : NULL;
    00203         }
    00204 }
    00205
    00206 template <class T>
    00207 void NRvector<T>::assign(int newn, const T& a)
    00208 {
    00209         if (newn != nn) {
    00210                 if (v != NULL) delete[] (v);
    00211                 nn = newn;
    00212                 v = nn > 0 ? new T[nn] : NULL;
    00213         }
    00214         for (int i=0;i<nn;i++) v[i] = a;
    00215 }
    00216
    00217 template <class T>
    00218 NRvector<T>::~NRvector()
    00219 {
    00220         if (v != NULL) delete[] (v);
    00221 }
    00222
    00223 // end of NRvector definitions
    00224
    00225 #endif //ifdef _USESTDVECTOR_
    00226
    00227 template <class T>
    00228 class NRmatrix {
    00229 private:
    00230         int nn;
    00231         int mm;
    00232         T **v;
    00233 public:
    00234         NRmatrix();
    00235         NRmatrix(int n, int m);                 // Zero-based array
    00236         NRmatrix(int n, int m, const T &a);     //Initialize to constant
    00237         NRmatrix(int n, int m, const T *a);     // Initialize to array
    00238         NRmatrix(const NRmatrix &rhs);          // Copy constructor
    00239         NRmatrix & operator=(const NRmatrix &rhs);      //assignment
    00240         typedef T value_type; // make T available externally
    00241         inline T* operator[](const int i);      //subscripting: pointer to row i
    00242         inline const T* operator[](const int i) const;
    00243         inline int nrows() const;
    00244         inline int ncols() const;
    00245         void resize(int newn, int newm); // resize (contents not preserved)
    00246         void assign(int newn, int newm, const T &a); // resize and assign a constant value
    00247         ~NRmatrix();
    00248 };
    00249
    00250 template <class T>
    00251 NRmatrix<T>::NRmatrix() : nn(0), mm(0), v(NULL) {}
    00252
    00253 template <class T>
    00254 NRmatrix<T>::NRmatrix(int n, int m) : nn(n), mm(m), v(n>0 ? new T*[n] : NULL)
    00255 {
    00256         int i,nel=m*n;
    00257         if (v) v[0] = nel>0 ? new T[nel] : NULL;
    00258         for (i=1;i<n;i++) v[i] = v[i-1] + m;
    00259 }
    00260
    00261 template <class T>
    00262 NRmatrix<T>::NRmatrix(int n, int m, const T &a) : nn(n), mm(m), v(n>0 ? new T*[n] : NULL)
    00263 {
    00264         int i,j,nel=m*n;
    00265         if (v) v[0] = nel>0 ? new T[nel] : NULL;
    00266         for (i=1; i< n; i++) v[i] = v[i-1] + m;
    00267         for (i=0; i< n; i++) for (j=0; j<m; j++) v[i][j] = a;
    00268 }
    00269
    00270 template <class T>
    00271 NRmatrix<T>::NRmatrix(int n, int m, const T *a) : nn(n), mm(m), v(n>0 ? new T*[n] : NULL)
    00272 {
    00273         int i,j,nel=m*n;
    00274         if (v) v[0] = nel>0 ? new T[nel] : NULL;
    00275         for (i=1; i< n; i++) v[i] = v[i-1] + m;
    00276         for (i=0; i< n; i++) for (j=0; j<m; j++) v[i][j] = *a++;
    00277 }
    00278
    00279 template <class T>
    00280 NRmatrix<T>::NRmatrix(const NRmatrix &rhs) : nn(rhs.nn), mm(rhs.mm), v(nn>0 ? new T*[nn] : NULL)
    00281 {
    00282         int i,j,nel=mm*nn;
    00283         if (v) v[0] = nel>0 ? new T[nel] : NULL;
    00284         for (i=1; i< nn; i++) v[i] = v[i-1] + mm;
    00285         for (i=0; i< nn; i++) for (j=0; j<mm; j++) v[i][j] = rhs[i][j];
    00286 }
    00287
    00288 template <class T>
    00289 NRmatrix<T> & NRmatrix<T>::operator=(const NRmatrix<T> &rhs)
    00290 // postcondition: normal assignment via copying has been performed;
    00291 //              if matrix and rhs were different sizes, matrix
    00292 //              has been resized to match the size of rhs
    00293 {
    00294         if (this != &rhs) {
    00295                 int i,j,nel;
    00296                 if (nn != rhs.nn || mm != rhs.mm) {
    00297                         if (v != NULL) {
    00298                                 delete[] (v[0]);
    00299                                 delete[] (v);
    00300                         }
    00301                         nn=rhs.nn;
    00302                         mm=rhs.mm;
    00303                         v = nn>0 ? new T*[nn] : NULL;
    00304                         nel = mm*nn;
    00305                         if (v) v[0] = nel>0 ? new T[nel] : NULL;
    00306                         for (i=1; i< nn; i++) v[i] = v[i-1] + mm;
    00307                 }
    00308                 for (i=0; i< nn; i++) for (j=0; j<mm; j++) v[i][j] = rhs[i][j];
    00309         }
    00310         return *this;
    00311 }
    00312
    00313 template <class T>
    00314 inline T* NRmatrix<T>::operator[](const int i)  //subscripting: pointer to row i
    00315 {
    00316 #ifdef _CHECKBOUNDS_
    00317 if (i<0 || i>=nn) {
    00318         throw("NRmatrix subscript out of bounds");
    00319 }
    00320 #endif
    00321         return v[i];
    00322 }
    00323
    00324 template <class T>
    00325 inline const T* NRmatrix<T>::operator[](const int i) const
    00326 {
    00327 #ifdef _CHECKBOUNDS_
    00328 if (i<0 || i>=nn) {
    00329         throw("NRmatrix subscript out of bounds");
    00330 }
    00331 #endif
    00332         return v[i];
    00333 }
    00334
    00335 template <class T>
    00336 inline int NRmatrix<T>::nrows() const
    00337 {
    00338         return nn;
    00339 }
    00340
    00341 template <class T>
    00342 inline int NRmatrix<T>::ncols() const
    00343 {
    00344         return mm;
    00345 }
    00346
    00347 template <class T>
    00348 void NRmatrix<T>::resize(int newn, int newm)
    00349 {
    00350         int i,nel;
    00351         if (newn != nn || newm != mm) {
    00352                 if (v != NULL) {
    00353                         delete[] (v[0]);
    00354                         delete[] (v);
    00355                 }
    00356                 nn = newn;
    00357                 mm = newm;
    00358                 v = nn>0 ? new T*[nn] : NULL;
    00359                 nel = mm*nn;
    00360                 if (v) v[0] = nel>0 ? new T[nel] : NULL;
    00361                 for (i=1; i< nn; i++) v[i] = v[i-1] + mm;
    00362         }
    00363 }
    00364
    00365 template <class T>
    00366 void NRmatrix<T>::assign(int newn, int newm, const T& a)
    00367 {
    00368         int i,j,nel;
    00369         if (newn != nn || newm != mm) {
    00370                 if (v != NULL) {
    00371                         delete[] (v[0]);
    00372                         delete[] (v);
    00373                 }
    00374                 nn = newn;
    00375                 mm = newm;
    00376                 v = nn>0 ? new T*[nn] : NULL;
    00377                 nel = mm*nn;
    00378                 if (v) v[0] = nel>0 ? new T[nel] : NULL;
    00379                 for (i=1; i< nn; i++) v[i] = v[i-1] + mm;
    00380         }
    00381         for (i=0; i< nn; i++) for (j=0; j<mm; j++) v[i][j] = a;
    00382 }
    00383
    00384 template <class T>
    00385 NRmatrix<T>::~NRmatrix()
    00386 {
    00387         if (v != NULL) {
    00388                 delete[] (v[0]);
    00389                 delete[] (v);
    00390         }
    00391 }
    00392
    00393 template <class T>
    00394 class NRMat3d {
    00395 private:
    00396         int nn;
    00397         int mm;
    00398         int kk;
    00399         T ***v;
    00400 public:
    00401         NRMat3d();
    00402         NRMat3d(int n, int m, int k);
    00403         inline T** operator[](const int i);     //subscripting: pointer to row i
    00404         inline const T* const * operator[](const int i) const;
    00405         inline int dim1() const;
    00406         inline int dim2() const;
    00407         inline int dim3() const;
    00408         ~NRMat3d();
    00409 };
    00410
    00411 template <class T>
    00412 NRMat3d<T>::NRMat3d(): nn(0), mm(0), kk(0), v(NULL) {}
    00413
    00414 template <class T>
    00415 NRMat3d<T>::NRMat3d(int n, int m, int k) : nn(n), mm(m), kk(k), v(new T**[n])
    00416 {
    00417         int i,j;
    00418         v[0] = new T*[n*m];
    00419         v[0][0] = new T[n*m*k];
    00420         for(j=1; j<m; j++) v[0][j] = v[0][j-1] + k;
    00421         for(i=1; i<n; i++) {
    00422                 v[i] = v[i-1] + m;
    00423                 v[i][0] = v[i-1][0] + m*k;
    00424                 for(j=1; j<m; j++) v[i][j] = v[i][j-1] + k;
    00425         }
    00426 }
    00427
    00428 template <class T>
    00429 inline T** NRMat3d<T>::operator[](const int i) //subscripting: pointer to row i
    00430 {
    00431         return v[i];
    00432 }
    00433
    00434 template <class T>
    00435 inline const T* const * NRMat3d<T>::operator[](const int i) const
    00436 {
    00437         return v[i];
    00438 }
    00439
    00440 template <class T>
    00441 inline int NRMat3d<T>::dim1() const
    00442 {
    00443         return nn;
    00444 }
    00445
    00446 template <class T>
    00447 inline int NRMat3d<T>::dim2() const
    00448 {
    00449         return mm;
    00450 }
    00451
    00452 template <class T>
    00453 inline int NRMat3d<T>::dim3() const
    00454 {
    00455         return kk;
    00456 }
    00457
    00458 template <class T>
    00459 NRMat3d<T>::~NRMat3d()
    00460 {
    00461         if (v != NULL) {
    00462                 delete[] (v[0][0]);
    00463                 delete[] (v[0]);
    00464                 delete[] (v);
    00465         }
    00466 }
    00467
    00468
    00469 // basic type names (redefine if your bit lengths don't match)
    00470
    00471 typedef int Int; // 32 bit integer
    00472 typedef unsigned int Uint;
    00473
    00474 #ifdef _MSC_VER
    00475 typedef __int64 Llong; // 64 bit integer
    00476 typedef unsigned __int64 Ullong;
    00477 #else
    00478 typedef long long int Llong; // 64 bit integer
    00479 typedef unsigned long long int Ullong;
    00480 #endif
    00481
    00482 typedef char Char; // 8 bit integer
    00483 typedef unsigned char Uchar;
    00484
    00485 typedef double Doub; // default floating type
    00486 typedef long double Ldoub;
    00487
    00488 typedef complex<double> Complex; // default complex type
    00489
    00490 typedef bool Bool;
    00491
    00492 // NaN: uncomment one of the following 3 methods of defining a global NaN
    00493 // you can test by verifying that (NaN != NaN) is true
    00494
    00495 static const Doub NaN = numeric_limits<Doub>::quiet_NaN();
    00496
    00497 //Uint proto_nan[2]={0xffffffff, 0x7fffffff};
    00498 //double NaN = *( double* )proto_nan;
    00499
    00500 //Doub NaN = sqrt(-1.);
    00501
    00502 // vector types
    00503
    00504 typedef const NRvector<Int> VecInt_I;
    00505 typedef NRvector<Int> VecInt, VecInt_O, VecInt_IO;
    00506
    00507 typedef const NRvector<Uint> VecUint_I;
    00508 typedef NRvector<Uint> VecUint, VecUint_O, VecUint_IO;
    00509
    00510 typedef const NRvector<Llong> VecLlong_I;
    00511 typedef NRvector<Llong> VecLlong, VecLlong_O, VecLlong_IO;
    00512
    00513 typedef const NRvector<Ullong> VecUllong_I;
    00514 typedef NRvector<Ullong> VecUllong, VecUllong_O, VecUllong_IO;
    00515
    00516 typedef const NRvector<Char> VecChar_I;
    00517 typedef NRvector<Char> VecChar, VecChar_O, VecChar_IO;
    00518
    00519 typedef const NRvector<Char*> VecCharp_I;
    00520 typedef NRvector<Char*> VecCharp, VecCharp_O, VecCharp_IO;
    00521
    00522 typedef const NRvector<Uchar> VecUchar_I;
    00523 typedef NRvector<Uchar> VecUchar, VecUchar_O, VecUchar_IO;
    00524
    00525 typedef const NRvector<Doub> VecDoub_I;
    00526 typedef NRvector<Doub> VecDoub, VecDoub_O, VecDoub_IO;
    00527
    00528 typedef const NRvector<Doub*> VecDoubp_I;
    00529 typedef NRvector<Doub*> VecDoubp, VecDoubp_O, VecDoubp_IO;
    00530
    00531 typedef const NRvector<Complex> VecComplex_I;
    00532 typedef NRvector<Complex> VecComplex, VecComplex_O, VecComplex_IO;
    00533
    00534 typedef const NRvector<Bool> VecBool_I;
    00535 typedef NRvector<Bool> VecBool, VecBool_O, VecBool_IO;
    00536
    00537 // matrix types
    00538
    00539 typedef const NRmatrix<Int> MatInt_I;
    00540 typedef NRmatrix<Int> MatInt, MatInt_O, MatInt_IO;
    00541
    00542 typedef const NRmatrix<Uint> MatUint_I;
    00543 typedef NRmatrix<Uint> MatUint, MatUint_O, MatUint_IO;
    00544
    00545 typedef const NRmatrix<Llong> MatLlong_I;
    00546 typedef NRmatrix<Llong> MatLlong, MatLlong_O, MatLlong_IO;
    00547
    00548 typedef const NRmatrix<Ullong> MatUllong_I;
    00549 typedef NRmatrix<Ullong> MatUllong, MatUllong_O, MatUllong_IO;
    00550
    00551 typedef const NRmatrix<Char> MatChar_I;
    00552 typedef NRmatrix<Char> MatChar, MatChar_O, MatChar_IO;
    00553
    00554 typedef const NRmatrix<Uchar> MatUchar_I;
    00555 typedef NRmatrix<Uchar> MatUchar, MatUchar_O, MatUchar_IO;
    00556
    00557 typedef const NRmatrix<Doub> MatDoub_I;
    00558 typedef NRmatrix<Doub> MatDoub, MatDoub_O, MatDoub_IO;
    00559
    00560 typedef const NRmatrix<Bool> MatBool_I;
    00561 typedef NRmatrix<Bool> MatBool, MatBool_O, MatBool_IO;
    00562
    00563 // 3D matrix types
    00564
    00565 typedef const NRMat3d<Doub> Mat3DDoub_I;
    00566 typedef NRMat3d<Doub> Mat3DDoub, Mat3DDoub_O, Mat3DDoub_IO;
    00567
    00568 // Floating Point Exceptions for Microsoft compilers
    00569
    00570 #ifdef _TURNONFPES_
    00571 #ifdef _MSC_VER
    00572 struct turn_on_floating_exceptions {
    00573         turn_on_floating_exceptions() {
    00574                 int cw = _controlfp( 0, 0 );
    00575                 cw &=~(EM_INVALID | EM_OVERFLOW | EM_ZERODIVIDE );
    00576                 _controlfp( cw, MCW_EM );
    00577         }
    00578 };
    00579 turn_on_floating_exceptions yes_turn_on_floating_exceptions;
    00580 #endif /* _MSC_VER */
    00581 #endif /* _TURNONFPES */
    00582
    00583 #endif /* _NR3_H_ */
    00584

All Classes Namespaces Files Functions Variables Enumerator Friends