TNT header files were missing. They are templates and don't have to be in the makefile
This commit is contained in:
Christian Schmitt
parent
ce1f89f832
commit
d10991cbb0
6 changed files with 1359 additions and 6 deletions
|
|
@ -18,10 +18,5 @@ add_library(Polygon STATIC
|
|||
tg_unique_vec3f.hxx
|
||||
trinodes.cxx
|
||||
trinodes.hxx
|
||||
TNT/jama_qr.h
|
||||
TNT/tnt_array1d.h
|
||||
TNT/tnt_array2d.h
|
||||
TNT/tnt_i_refvec.h
|
||||
TNT/tnt_math_utils.h
|
||||
point3d.hxx
|
||||
)
|
||||
)
|
||||
|
|
|
|||
326
src/Lib/Polygon/TNT/jama_qr.h
Normal file
326
src/Lib/Polygon/TNT/jama_qr.h
Normal file
|
|
@ -0,0 +1,326 @@
|
|||
#ifndef JAMA_QR_H
|
||||
#define JAMA_QR_H
|
||||
|
||||
#include "tnt_array1d.h"
|
||||
#include "tnt_array2d.h"
|
||||
#include "tnt_math_utils.h"
|
||||
|
||||
namespace JAMA
|
||||
{
|
||||
|
||||
/**
|
||||
<p>
|
||||
Classical QR Decompisition:
|
||||
for an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n
|
||||
orthogonal matrix Q and an n-by-n upper triangular matrix R so that
|
||||
A = Q*R.
|
||||
<P>
|
||||
The QR decompostion always exists, even if the matrix does not have
|
||||
full rank, so the constructor will never fail. The primary use of the
|
||||
QR decomposition is in the least squares solution of nonsquare systems
|
||||
of simultaneous linear equations. This will fail if isFullRank()
|
||||
returns 0 (false).
|
||||
|
||||
<p>
|
||||
The Q and R factors can be retrived via the getQ() and getR()
|
||||
methods. Furthermore, a solve() method is provided to find the
|
||||
least squares solution of Ax=b using the QR factors.
|
||||
|
||||
<p>
|
||||
(Adapted from JAMA, a Java Matrix Library, developed by jointly
|
||||
by the Mathworks and NIST; see http://math.nist.gov/javanumerics/jama).
|
||||
*/
|
||||
|
||||
template <class Real>
|
||||
class QR {
|
||||
|
||||
|
||||
/** Array for internal storage of decomposition.
|
||||
@serial internal array storage.
|
||||
*/
|
||||
|
||||
TNT::Array2D<Real> QR_;
|
||||
|
||||
/** Row and column dimensions.
|
||||
@serial column dimension.
|
||||
@serial row dimension.
|
||||
*/
|
||||
int m, n;
|
||||
|
||||
/** Array for internal storage of diagonal of R.
|
||||
@serial diagonal of R.
|
||||
*/
|
||||
TNT::Array1D<Real> Rdiag;
|
||||
|
||||
|
||||
public:
|
||||
|
||||
/**
|
||||
Create a QR factorization object for A.
|
||||
|
||||
@param A rectangular (m>=n) matrix.
|
||||
*/
|
||||
QR(const TNT::Array2D<Real> &A) /* constructor */
|
||||
{
|
||||
QR_ = A.copy();
|
||||
m = A.dim1();
|
||||
n = A.dim2();
|
||||
Rdiag = TNT::Array1D<Real>(n);
|
||||
int i=0, j=0, k=0;
|
||||
|
||||
// Main loop.
|
||||
for (k = 0; k < n; k++) {
|
||||
// Compute 2-norm of k-th column without under/overflow.
|
||||
Real nrm = 0;
|
||||
for (i = k; i < m; i++) {
|
||||
nrm = TNT::hypot(nrm,QR_[i][k]);
|
||||
}
|
||||
|
||||
if (nrm != 0.0) {
|
||||
// Form k-th Householder vector.
|
||||
if (QR_[k][k] < 0) {
|
||||
nrm = -nrm;
|
||||
}
|
||||
for (i = k; i < m; i++) {
|
||||
QR_[i][k] /= nrm;
|
||||
}
|
||||
QR_[k][k] += 1.0;
|
||||
|
||||
// Apply transformation to remaining columns.
|
||||
for (j = k+1; j < n; j++) {
|
||||
Real s = 0.0;
|
||||
for (i = k; i < m; i++) {
|
||||
s += QR_[i][k]*QR_[i][j];
|
||||
}
|
||||
s = -s/QR_[k][k];
|
||||
for (i = k; i < m; i++) {
|
||||
QR_[i][j] += s*QR_[i][k];
|
||||
}
|
||||
}
|
||||
}
|
||||
Rdiag[k] = -nrm;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
Flag to denote the matrix is of full rank.
|
||||
|
||||
@return 1 if matrix is full rank, 0 otherwise.
|
||||
*/
|
||||
int isFullRank() const
|
||||
{
|
||||
for (int j = 0; j < n; j++)
|
||||
{
|
||||
if (Rdiag[j] == 0)
|
||||
return 0;
|
||||
}
|
||||
return 1;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
/**
|
||||
|
||||
Retreive the Householder vectors from QR factorization
|
||||
@returns lower trapezoidal matrix whose columns define the reflections
|
||||
*/
|
||||
|
||||
TNT::Array2D<Real> getHouseholder (void) const
|
||||
{
|
||||
TNT::Array2D<Real> H(m,n);
|
||||
|
||||
/* note: H is completely filled in by algorithm, so
|
||||
initializaiton of H is not necessary.
|
||||
*/
|
||||
for (int i = 0; i < m; i++)
|
||||
{
|
||||
for (int j = 0; j < n; j++)
|
||||
{
|
||||
if (i >= j) {
|
||||
H[i][j] = QR_[i][j];
|
||||
} else {
|
||||
H[i][j] = 0.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
return H;
|
||||
}
|
||||
|
||||
|
||||
|
||||
/** Return the upper triangular factor, R, of the QR factorization
|
||||
@return R
|
||||
*/
|
||||
|
||||
TNT::Array2D<Real> getR() const
|
||||
{
|
||||
TNT::Array2D<Real> R(n,n);
|
||||
for (int i = 0; i < n; i++) {
|
||||
for (int j = 0; j < n; j++) {
|
||||
if (i < j) {
|
||||
R[i][j] = QR_[i][j];
|
||||
} else if (i == j) {
|
||||
R[i][j] = Rdiag[i];
|
||||
} else {
|
||||
R[i][j] = 0.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
return R;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
/**
|
||||
Generate and return the (economy-sized) orthogonal factor
|
||||
@param Q the (ecnomy-sized) orthogonal factor (Q*R=A).
|
||||
*/
|
||||
|
||||
TNT::Array2D<Real> getQ() const
|
||||
{
|
||||
int i=0, j=0, k=0;
|
||||
|
||||
TNT::Array2D<Real> Q(m,n);
|
||||
for (k = n-1; k >= 0; k--) {
|
||||
for (i = 0; i < m; i++) {
|
||||
Q[i][k] = 0.0;
|
||||
}
|
||||
Q[k][k] = 1.0;
|
||||
for (j = k; j < n; j++) {
|
||||
if (QR_[k][k] != 0) {
|
||||
Real s = 0.0;
|
||||
for (i = k; i < m; i++) {
|
||||
s += QR_[i][k]*Q[i][j];
|
||||
}
|
||||
s = -s/QR_[k][k];
|
||||
for (i = k; i < m; i++) {
|
||||
Q[i][j] += s*QR_[i][k];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
return Q;
|
||||
}
|
||||
|
||||
|
||||
/** Least squares solution of A*x = b
|
||||
@param B m-length array (vector).
|
||||
@return x n-length array (vector) that minimizes the two norm of Q*R*X-B.
|
||||
If B is non-conformant, or if QR.isFullRank() is false,
|
||||
the routine returns a null (0-length) vector.
|
||||
*/
|
||||
|
||||
TNT::Array1D<Real> solve(const TNT::Array1D<Real> &b) const
|
||||
{
|
||||
if (b.dim1() != m) /* arrays must be conformant */
|
||||
return TNT::Array1D<Real>();
|
||||
|
||||
if ( !isFullRank() ) /* matrix is rank deficient */
|
||||
{
|
||||
return TNT::Array1D<Real>();
|
||||
}
|
||||
|
||||
TNT::Array1D<Real> x = b.copy();
|
||||
|
||||
// Compute Y = transpose(Q)*b
|
||||
for (int k = 0; k < n; k++)
|
||||
{
|
||||
Real s = 0.0;
|
||||
for (int i = k; i < m; i++)
|
||||
{
|
||||
s += QR_[i][k]*x[i];
|
||||
}
|
||||
s = -s/QR_[k][k];
|
||||
for (int i = k; i < m; i++)
|
||||
{
|
||||
x[i] += s*QR_[i][k];
|
||||
}
|
||||
}
|
||||
// Solve R*X = Y;
|
||||
for (int k = n-1; k >= 0; k--)
|
||||
{
|
||||
x[k] /= Rdiag[k];
|
||||
for (int i = 0; i < k; i++) {
|
||||
x[i] -= x[k]*QR_[i][k];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* return n x nx portion of X */
|
||||
TNT::Array1D<Real> x_(n);
|
||||
for (int i=0; i<n; i++)
|
||||
x_[i] = x[i];
|
||||
|
||||
return x_;
|
||||
}
|
||||
|
||||
/** Least squares solution of A*X = B
|
||||
@param B m x k Array (must conform).
|
||||
@return X n x k Array that minimizes the two norm of Q*R*X-B. If
|
||||
B is non-conformant, or if QR.isFullRank() is false,
|
||||
the routine returns a null (0x0) array.
|
||||
*/
|
||||
|
||||
TNT::Array2D<Real> solve(const TNT::Array2D<Real> &B) const
|
||||
{
|
||||
if (B.dim1() != m) /* arrays must be conformant */
|
||||
return TNT::Array2D<Real>(0,0);
|
||||
|
||||
if ( !isFullRank() ) /* matrix is rank deficient */
|
||||
{
|
||||
return TNT::Array2D<Real>(0,0);
|
||||
}
|
||||
|
||||
int nx = B.dim2();
|
||||
TNT::Array2D<Real> X = B.copy();
|
||||
int i=0, j=0, k=0;
|
||||
|
||||
// Compute Y = transpose(Q)*B
|
||||
for (k = 0; k < n; k++) {
|
||||
for (j = 0; j < nx; j++) {
|
||||
Real s = 0.0;
|
||||
for (i = k; i < m; i++) {
|
||||
s += QR_[i][k]*X[i][j];
|
||||
}
|
||||
s = -s/QR_[k][k];
|
||||
for (i = k; i < m; i++) {
|
||||
X[i][j] += s*QR_[i][k];
|
||||
}
|
||||
}
|
||||
}
|
||||
// Solve R*X = Y;
|
||||
for (k = n-1; k >= 0; k--) {
|
||||
for (j = 0; j < nx; j++) {
|
||||
X[k][j] /= Rdiag[k];
|
||||
}
|
||||
for (i = 0; i < k; i++) {
|
||||
for (j = 0; j < nx; j++) {
|
||||
X[i][j] -= X[k][j]*QR_[i][k];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* return n x nx portion of X */
|
||||
TNT::Array2D<Real> X_(n,nx);
|
||||
for (i=0; i<n; i++)
|
||||
for (j=0; j<nx; j++)
|
||||
X_[i][j] = X[i][j];
|
||||
|
||||
return X_;
|
||||
}
|
||||
|
||||
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
// namespace JAMA
|
||||
|
||||
#endif
|
||||
// JAMA_QR__H
|
||||
|
||||
301
src/Lib/Polygon/TNT/tnt_array1d.h
Normal file
301
src/Lib/Polygon/TNT/tnt_array1d.h
Normal file
|
|
@ -0,0 +1,301 @@
|
|||
/*
|
||||
*
|
||||
* Template Numerical Toolkit (TNT)
|
||||
*
|
||||
* Mathematical and Computational Sciences Division
|
||||
* National Institute of Technology,
|
||||
* Gaithersburg, MD USA
|
||||
*
|
||||
*
|
||||
* This software was developed at the National Institute of Standards and
|
||||
* Technology (NIST) by employees of the Federal Government in the course
|
||||
* of their official duties. Pursuant to title 17 Section 105 of the
|
||||
* United States Code, this software is not subject to copyright protection
|
||||
* and is in the public domain. NIST assumes no responsibility whatsoever for
|
||||
* its use by other parties, and makes no guarantees, expressed or implied,
|
||||
* about its quality, reliability, or any other characteristic.
|
||||
*
|
||||
*/
|
||||
|
||||
|
||||
|
||||
#ifndef TNT_ARRAY1D_H
|
||||
#define TNT_ARRAY1D_H
|
||||
|
||||
//#include <cstdlib>
|
||||
#include <iostream>
|
||||
|
||||
#ifdef TNT_BOUNDS_CHECK
|
||||
#include <assert.h>
|
||||
#endif
|
||||
|
||||
|
||||
#include "tnt_i_refvec.h"
|
||||
|
||||
namespace TNT
|
||||
{
|
||||
|
||||
template <class T>
|
||||
class Array1D
|
||||
{
|
||||
|
||||
private:
|
||||
|
||||
/* ... */
|
||||
i_refvec<T> v_;
|
||||
int n_;
|
||||
T* data_; /* this normally points to v_.begin(), but
|
||||
* could also point to a portion (subvector)
|
||||
* of v_.
|
||||
*/
|
||||
|
||||
void copy_(T* p, const T* q, int len) const;
|
||||
void set_(T* begin, T* end, const T& val);
|
||||
|
||||
|
||||
public:
|
||||
|
||||
typedef T value_type;
|
||||
|
||||
|
||||
Array1D();
|
||||
explicit Array1D(int n);
|
||||
Array1D(int n, const T &a);
|
||||
Array1D(int n, T *a);
|
||||
inline Array1D(const Array1D &A);
|
||||
inline operator T*();
|
||||
inline operator const T*();
|
||||
inline Array1D & operator=(const T &a);
|
||||
inline Array1D & operator=(const Array1D &A);
|
||||
inline Array1D & ref(const Array1D &A);
|
||||
Array1D copy() const;
|
||||
Array1D & inject(const Array1D & A);
|
||||
inline T& operator[](int i);
|
||||
inline T& operator()(int i);
|
||||
inline const T& operator[](int i) const;
|
||||
inline const T& operator()(int i) const;
|
||||
inline int dim1() const;
|
||||
inline int dim() const;
|
||||
~Array1D();
|
||||
|
||||
|
||||
/* ... extended interface ... */
|
||||
|
||||
inline int ref_count() const;
|
||||
inline Array1D<T> subarray(int i0, int i1);
|
||||
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
template <class T>
|
||||
Array1D<T>::Array1D() : v_(), n_(0), data_(0) {}
|
||||
|
||||
template <class T>
|
||||
Array1D<T>::Array1D(const Array1D<T> &A) : v_(A.v_), n_(A.n_),
|
||||
data_(A.data_)
|
||||
{
|
||||
#ifdef TNT_DEBUG
|
||||
std::cout << "Created Array1D(const Array1D<T> &A) \n";
|
||||
#endif
|
||||
|
||||
}
|
||||
|
||||
template <class T>
|
||||
Array1D<T>::Array1D(int n) : v_(n), n_(n), data_(v_.begin())
|
||||
{
|
||||
#ifdef TNT_DEBUG
|
||||
std::cout << "Created Array1D(int n) \n";
|
||||
#endif
|
||||
}
|
||||
|
||||
template <class T>
|
||||
Array1D<T>::Array1D(int n, const T &val) : v_(n), n_(n), data_(v_.begin())
|
||||
{
|
||||
#ifdef TNT_DEBUG
|
||||
std::cout << "Created Array1D(int n, const T& val) \n";
|
||||
#endif
|
||||
set_(data_, data_+ n, val);
|
||||
|
||||
}
|
||||
|
||||
template <class T>
|
||||
Array1D<T>::Array1D(int n, T *a) : v_(a), n_(n) , data_(v_.begin())
|
||||
{
|
||||
#ifdef TNT_DEBUG
|
||||
std::cout << "Created Array1D(int n, T* a) \n";
|
||||
#endif
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline Array1D<T>::operator T*()
|
||||
{
|
||||
return &(v_[0]);
|
||||
}
|
||||
|
||||
|
||||
template <class T>
|
||||
inline Array1D<T>::operator const T*()
|
||||
{
|
||||
return &(v_[0]);
|
||||
}
|
||||
|
||||
|
||||
|
||||
template <class T>
|
||||
inline T& Array1D<T>::operator[](int i)
|
||||
{
|
||||
#ifdef TNT_BOUNDS_CHECK
|
||||
assert(i>= 0);
|
||||
assert(i < n_);
|
||||
#endif
|
||||
return data_[i];
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline const T& Array1D<T>::operator[](int i) const
|
||||
{
|
||||
#ifdef TNT_BOUNDS_CHECK
|
||||
assert(i>= 0);
|
||||
assert(i < n_);
|
||||
#endif
|
||||
return data_[i];
|
||||
}
|
||||
|
||||
|
||||
|
||||
template <class T>
|
||||
inline T& Array1D<T>::operator()(int i)
|
||||
{
|
||||
#ifdef TNT_BOUNDS_CHECK
|
||||
assert(i>= 0);
|
||||
assert(i < n_);
|
||||
#endif
|
||||
return data_[i-1];
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline const T& Array1D<T>::operator()(int i) const
|
||||
{
|
||||
#ifdef TNT_BOUNDS_CHECK
|
||||
assert(i>= 0);
|
||||
assert(i < n_);
|
||||
#endif
|
||||
return data_[i-1];
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
template <class T>
|
||||
Array1D<T> & Array1D<T>::operator=(const T &a)
|
||||
{
|
||||
set_(data_, data_+n_, a);
|
||||
return *this;
|
||||
}
|
||||
|
||||
template <class T>
|
||||
Array1D<T> Array1D<T>::copy() const
|
||||
{
|
||||
Array1D A( n_);
|
||||
copy_(A.data_, data_, n_);
|
||||
|
||||
return A;
|
||||
}
|
||||
|
||||
|
||||
template <class T>
|
||||
Array1D<T> & Array1D<T>::inject(const Array1D &A)
|
||||
{
|
||||
if (A.n_ == n_)
|
||||
copy_(data_, A.data_, n_);
|
||||
|
||||
return *this;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template <class T>
|
||||
Array1D<T> & Array1D<T>::ref(const Array1D<T> &A)
|
||||
{
|
||||
if (this != &A)
|
||||
{
|
||||
v_ = A.v_; /* operator= handles the reference counting. */
|
||||
n_ = A.n_;
|
||||
data_ = A.data_;
|
||||
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
template <class T>
|
||||
Array1D<T> & Array1D<T>::operator=(const Array1D<T> &A)
|
||||
{
|
||||
return ref(A);
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline int Array1D<T>::dim1() const { return n_; }
|
||||
|
||||
template <class T>
|
||||
inline int Array1D<T>::dim() const { return n_; }
|
||||
|
||||
template <class T>
|
||||
Array1D<T>::~Array1D() {}
|
||||
|
||||
|
||||
/* ............................ exented interface ......................*/
|
||||
|
||||
template <class T>
|
||||
inline int Array1D<T>::ref_count() const
|
||||
{
|
||||
return v_.ref_count();
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline Array1D<T> Array1D<T>::subarray(int i0, int i1)
|
||||
{
|
||||
if (((i0 >= 0) && (i1 < n_)) || (i0 <= i1))
|
||||
{
|
||||
Array1D<T> X(*this); /* create a new instance of this array. */
|
||||
X.n_ = i1-i0+1;
|
||||
X.data_ += i0;
|
||||
|
||||
return X;
|
||||
}
|
||||
else
|
||||
{
|
||||
return Array1D<T>();
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* private internal functions */
|
||||
|
||||
|
||||
template <class T>
|
||||
void Array1D<T>::set_(T* begin, T* end, const T& a)
|
||||
{
|
||||
for (T* p=begin; p<end; p++)
|
||||
*p = a;
|
||||
|
||||
}
|
||||
|
||||
template <class T>
|
||||
void Array1D<T>::copy_(T* p, const T* q, int len) const
|
||||
{
|
||||
T *end = p + len;
|
||||
while (p<end )
|
||||
*p++ = *q++;
|
||||
|
||||
}
|
||||
|
||||
|
||||
} /* namespace TNT */
|
||||
|
||||
#endif
|
||||
/* TNT_ARRAY1D_H */
|
||||
|
||||
451
src/Lib/Polygon/TNT/tnt_array2d.h
Normal file
451
src/Lib/Polygon/TNT/tnt_array2d.h
Normal file
|
|
@ -0,0 +1,451 @@
|
|||
/*
|
||||
*
|
||||
* Template Numerical Toolkit (TNT)
|
||||
*
|
||||
* Mathematical and Computational Sciences Division
|
||||
* National Institute of Technology,
|
||||
* Gaithersburg, MD USA
|
||||
*
|
||||
*
|
||||
* This software was developed at the National Institute of Standards and
|
||||
* Technology (NIST) by employees of the Federal Government in the course
|
||||
* of their official duties. Pursuant to title 17 Section 105 of the
|
||||
* United States Code, this software is not subject to copyright protection
|
||||
* and is in the public domain. NIST assumes no responsibility whatsoever for
|
||||
* its use by other parties, and makes no guarantees, expressed or implied,
|
||||
* about its quality, reliability, or any other characteristic.
|
||||
*
|
||||
*/
|
||||
|
||||
|
||||
|
||||
#ifndef TNT_ARRAY2D_H
|
||||
#define TNT_ARRAY2D_H
|
||||
|
||||
#include <cstdlib>
|
||||
#include <iostream>
|
||||
#ifdef TNT_BOUNDS_CHECK
|
||||
#include <assert.h>
|
||||
#endif
|
||||
|
||||
#include "tnt_array1d.h"
|
||||
|
||||
namespace TNT
|
||||
{
|
||||
|
||||
|
||||
/**
|
||||
Ttwo-dimensional numerical array which
|
||||
looks like a conventional C multiarray.
|
||||
Storage corresponds to C (row-major) ordering.
|
||||
Elements are accessed via A[i][j] notation for 0-based indexing,
|
||||
and A(i,j) for 1-based indexing..
|
||||
|
||||
<p>
|
||||
Array assignment is by reference (i.e. shallow assignment).
|
||||
That is, B=A implies that the A and B point to the
|
||||
same array, so modifications to the elements of A
|
||||
will be reflected in B. If an independent copy
|
||||
is required, then B = A.copy() can be used. Note
|
||||
that this facilitates returning arrays from functions
|
||||
without relying on compiler optimizations to eliminate
|
||||
extensive data copying.
|
||||
|
||||
<p>
|
||||
The indexing and layout of this array object makes
|
||||
it compatible with C and C++ algorithms that utilize
|
||||
the familiar C[i][j] notation. This includes numerous
|
||||
|
||||
textbooks, such as Numercial Recipes, and various
|
||||
public domain codes.
|
||||
|
||||
*/
|
||||
template <class T>
|
||||
class Array2D
|
||||
{
|
||||
|
||||
|
||||
private:
|
||||
|
||||
|
||||
|
||||
Array1D<T> data_;
|
||||
Array1D<T*> v_;
|
||||
int m_;
|
||||
int n_;
|
||||
|
||||
public:
|
||||
|
||||
|
||||
|
||||
/**
|
||||
Used to determined the data type of array entries. This is most
|
||||
commonly used when requiring scalar temporaries in templated algorithms
|
||||
that have TNT arrays as input. For example,
|
||||
<pre>
|
||||
template < class ArrayTwoD >
|
||||
void foo(ArrayTwoD &A)
|
||||
{
|
||||
A::value_type first_entry = A[0][0];
|
||||
...
|
||||
}
|
||||
</pre>
|
||||
*/
|
||||
typedef T value_type;
|
||||
|
||||
/**
|
||||
Create a null array. This is <b>not</b> the same
|
||||
as Array2D(0,0), which consumes some memory overhead.
|
||||
*/
|
||||
|
||||
Array2D();
|
||||
|
||||
|
||||
/**
|
||||
Create a new (m x n) array, without initalizing elements. (This
|
||||
encurs an O(1) operation cost, rather than a O(m*n) cost.)
|
||||
|
||||
@param m the first (row) dimension of the new matrix.
|
||||
@param n the second (column) dimension of the new matrix.
|
||||
*/
|
||||
Array2D(int m, int n);
|
||||
|
||||
|
||||
/**
|
||||
Create a new (m x n) array, as a view of an existing one-dimensional
|
||||
array stored in row-major order, i.e. right-most dimension varying fastest.
|
||||
Note that the storage for this pre-existing array will
|
||||
never be destroyed by TNT.
|
||||
|
||||
@param m the first (row) dimension of the new matrix.
|
||||
@param n the second (column) dimension of the new matrix.
|
||||
@param a the one dimensional C array to use as data storage for
|
||||
the array.
|
||||
*/
|
||||
Array2D(int m, int n, T *a);
|
||||
|
||||
|
||||
|
||||
/**
|
||||
Create a new (m x n) array, initializing array elements to
|
||||
constant specified by argument. Most often used to
|
||||
create an array of zeros, as in A(m, n, 0.0).
|
||||
|
||||
@param m the first (row) dimension of the new matrix.
|
||||
@param n the second (column) dimension of the new matrix.
|
||||
@param val the constant value to set all elements of the new array to.
|
||||
*/
|
||||
|
||||
Array2D(int m, int n, const T &val);
|
||||
|
||||
|
||||
/**
|
||||
Copy constructor. Array data is NOT copied, but shared.
|
||||
Thus, in Array2D B(A), subsequent changes to A will
|
||||
be reflected in B. For an indepent copy of A, use
|
||||
Array2D B(A.copy()), or B = A.copy(), instead.
|
||||
*/
|
||||
inline Array2D(const Array2D &A);
|
||||
|
||||
|
||||
/**
|
||||
Convert 2D array into a regular multidimensional C pointer. Most often
|
||||
called automatically when calling C interfaces that expect things like
|
||||
double** rather than Array2D<dobule>.
|
||||
|
||||
*/
|
||||
inline operator T**();
|
||||
|
||||
/**
|
||||
Convert a const 2D array into a const multidimensional C pointer.
|
||||
Most often called automatically when calling C interfaces that expect
|
||||
things like "const double**" rather than "const Array2D<dobule>&".
|
||||
|
||||
*/
|
||||
|
||||
inline operator const T**() const;
|
||||
|
||||
|
||||
/**
|
||||
Assign all elements of array the same value.
|
||||
|
||||
@param val the value to assign each element.
|
||||
*/
|
||||
inline Array2D & operator=(const T &val);
|
||||
|
||||
|
||||
/**
|
||||
Assign one Array2D to another. (This is a shallow-assignement operation,
|
||||
and it is the identical semantics to ref(A).
|
||||
|
||||
@param A the array to assign this one to.
|
||||
*/
|
||||
inline Array2D & operator=(const Array2D &A);
|
||||
|
||||
|
||||
inline Array2D & ref(const Array2D &A);
|
||||
Array2D copy() const;
|
||||
Array2D & inject(const Array2D & A);
|
||||
inline T* operator[](int i);
|
||||
inline const T* operator[](int i) const;
|
||||
inline int dim1() const;
|
||||
inline int dim2() const;
|
||||
~Array2D();
|
||||
|
||||
/* extended interface (not part of the standard) */
|
||||
|
||||
|
||||
inline int ref_count();
|
||||
inline int ref_count_data();
|
||||
inline int ref_count_dim1();
|
||||
Array2D subarray(int i0, int i1, int j0, int j1);
|
||||
|
||||
};
|
||||
|
||||
|
||||
/**
|
||||
Create a new (m x n) array, WIHOUT initializing array elements.
|
||||
To create an initialized array of constants, see Array2D(m,n,value).
|
||||
|
||||
<p>
|
||||
This version avoids the O(m*n) initialization overhead and
|
||||
is used just before manual assignment.
|
||||
|
||||
@param m the first (row) dimension of the new matrix.
|
||||
@param n the second (column) dimension of the new matrix.
|
||||
*/
|
||||
template <class T>
|
||||
Array2D<T>::Array2D() : data_(), v_(), m_(0), n_(0) {}
|
||||
|
||||
|
||||
template <class T>
|
||||
Array2D<T>::Array2D(const Array2D<T> &A) : data_(A.data_), v_(A.v_),
|
||||
m_(A.m_), n_(A.n_) {}
|
||||
|
||||
|
||||
|
||||
|
||||
template <class T>
|
||||
Array2D<T>::Array2D(int m, int n) : data_(m*n), v_(m), m_(m), n_(n)
|
||||
{
|
||||
if (m>0 && n>0)
|
||||
{
|
||||
T* p = &(data_[0]);
|
||||
for (int i=0; i<m; i++)
|
||||
{
|
||||
v_[i] = p;
|
||||
p += n;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
template <class T>
|
||||
Array2D<T>::Array2D(int m, int n, const T &val) : data_(m*n), v_(m),
|
||||
m_(m), n_(n)
|
||||
{
|
||||
if (m>0 && n>0)
|
||||
{
|
||||
data_ = val;
|
||||
T* p = &(data_[0]);
|
||||
for (int i=0; i<m; i++)
|
||||
{
|
||||
v_[i] = p;
|
||||
p += n;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template <class T>
|
||||
Array2D<T>::Array2D(int m, int n, T *a) : data_(m*n, a), v_(m), m_(m), n_(n)
|
||||
{
|
||||
if (m>0 && n>0)
|
||||
{
|
||||
T* p = &(data_[0]);
|
||||
|
||||
for (int i=0; i<m; i++)
|
||||
{
|
||||
v_[i] = p;
|
||||
p += n;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template <class T>
|
||||
inline T* Array2D<T>::operator[](int i)
|
||||
{
|
||||
#ifdef TNT_BOUNDS_CHECK
|
||||
assert(i >= 0);
|
||||
assert(i < m_);
|
||||
#endif
|
||||
|
||||
return v_[i];
|
||||
|
||||
}
|
||||
|
||||
|
||||
template <class T>
|
||||
inline const T* Array2D<T>::operator[](int i) const
|
||||
{
|
||||
#ifdef TNT_BOUNDS_CHECK
|
||||
assert(i >= 0);
|
||||
assert(i < m_);
|
||||
#endif
|
||||
|
||||
return v_[i];
|
||||
|
||||
}
|
||||
|
||||
template <class T>
|
||||
Array2D<T> & Array2D<T>::operator=(const T &a)
|
||||
{
|
||||
/* non-optimzied, but will work with subarrays in future verions */
|
||||
|
||||
for (int i=0; i<m_; i++)
|
||||
for (int j=0; j<n_; j++)
|
||||
v_[i][j] = a;
|
||||
return *this;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
template <class T>
|
||||
Array2D<T> Array2D<T>::copy() const
|
||||
{
|
||||
Array2D A(m_, n_);
|
||||
|
||||
for (int i=0; i<m_; i++)
|
||||
for (int j=0; j<n_; j++)
|
||||
A[i][j] = v_[i][j];
|
||||
|
||||
|
||||
return A;
|
||||
}
|
||||
|
||||
|
||||
template <class T>
|
||||
Array2D<T> & Array2D<T>::inject(const Array2D &A)
|
||||
{
|
||||
if (A.m_ == m_ && A.n_ == n_)
|
||||
{
|
||||
for (int i=0; i<m_; i++)
|
||||
for (int j=0; j<n_; j++)
|
||||
v_[i][j] = A[i][j];
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
template <class T>
|
||||
Array2D<T> & Array2D<T>::ref(const Array2D<T> &A)
|
||||
{
|
||||
if (this != &A)
|
||||
{
|
||||
v_ = A.v_;
|
||||
data_ = A.data_;
|
||||
m_ = A.m_;
|
||||
n_ = A.n_;
|
||||
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
|
||||
|
||||
template <class T>
|
||||
Array2D<T> & Array2D<T>::operator=(const Array2D<T> &A)
|
||||
{
|
||||
return ref(A);
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline int Array2D<T>::dim1() const { return m_; }
|
||||
|
||||
template <class T>
|
||||
inline int Array2D<T>::dim2() const { return n_; }
|
||||
|
||||
|
||||
template <class T>
|
||||
Array2D<T>::~Array2D() {}
|
||||
|
||||
|
||||
|
||||
|
||||
template <class T>
|
||||
inline Array2D<T>::operator T**()
|
||||
{
|
||||
return &(v_[0]);
|
||||
}
|
||||
template <class T>
|
||||
inline Array2D<T>::operator const T**() const
|
||||
{
|
||||
return static_cast<const T**>(&(v_[0]));
|
||||
}
|
||||
|
||||
/* ............... extended interface ............... */
|
||||
/**
|
||||
Create a new view to a subarray defined by the boundaries
|
||||
[i0][i0] and [i1][j1]. The size of the subarray is
|
||||
(i1-i0) by (j1-j0). If either of these lengths are zero
|
||||
or negative, the subarray view is null.
|
||||
|
||||
*/
|
||||
template <class T>
|
||||
Array2D<T> Array2D<T>::subarray(int i0, int i1, int j0, int j1)
|
||||
{
|
||||
Array2D<T> A;
|
||||
int m = i1-i0+1;
|
||||
int n = j1-j0+1;
|
||||
|
||||
/* if either length is zero or negative, this is an invalide
|
||||
subarray. return a null view.
|
||||
*/
|
||||
if (m<1 || n<1)
|
||||
return A;
|
||||
|
||||
A.data_ = data_;
|
||||
A.m_ = m;
|
||||
A.n_ = n;
|
||||
A.v_ = Array1D<T*>(m);
|
||||
T* p = &(data_[0]) + i0 * n_ + j0;
|
||||
for (int i=0; i<m; i++)
|
||||
{
|
||||
A.v_[i] = p + i*n_;
|
||||
|
||||
}
|
||||
return A;
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline int Array2D<T>::ref_count()
|
||||
{
|
||||
return ref_count_data();
|
||||
}
|
||||
|
||||
|
||||
|
||||
template <class T>
|
||||
inline int Array2D<T>::ref_count_data()
|
||||
{
|
||||
return data_.ref_count();
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline int Array2D<T>::ref_count_dim1()
|
||||
{
|
||||
return v_.ref_count();
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
} /* namespace TNT */
|
||||
|
||||
#endif
|
||||
/* TNT_ARRAY2D_H */
|
||||
|
||||
243
src/Lib/Polygon/TNT/tnt_i_refvec.h
Normal file
243
src/Lib/Polygon/TNT/tnt_i_refvec.h
Normal file
|
|
@ -0,0 +1,243 @@
|
|||
/*
|
||||
*
|
||||
* Template Numerical Toolkit (TNT)
|
||||
*
|
||||
* Mathematical and Computational Sciences Division
|
||||
* National Institute of Technology,
|
||||
* Gaithersburg, MD USA
|
||||
*
|
||||
*
|
||||
* This software was developed at the National Institute of Standards and
|
||||
* Technology (NIST) by employees of the Federal Government in the course
|
||||
* of their official duties. Pursuant to title 17 Section 105 of the
|
||||
* United States Code, this software is not subject to copyright protection
|
||||
* and is in the public domain. NIST assumes no responsibility whatsoever for
|
||||
* its use by other parties, and makes no guarantees, expressed or implied,
|
||||
* about its quality, reliability, or any other characteristic.
|
||||
*
|
||||
*/
|
||||
|
||||
|
||||
|
||||
#ifndef TNT_I_REFVEC_H
|
||||
#define TNT_I_REFVEC_H
|
||||
|
||||
#include <cstdlib>
|
||||
#include <iostream>
|
||||
|
||||
#ifdef TNT_BOUNDS_CHECK
|
||||
#include <assert.h>
|
||||
#endif
|
||||
|
||||
#ifndef NULL
|
||||
#define NULL 0
|
||||
#endif
|
||||
|
||||
namespace TNT
|
||||
{
|
||||
/*
|
||||
Internal representation of ref-counted array. The TNT
|
||||
arrays all use this building block.
|
||||
|
||||
<p>
|
||||
If an array block is created by TNT, then every time
|
||||
an assignment is made, the left-hand-side reference
|
||||
is decreased by one, and the right-hand-side refernce
|
||||
count is increased by one. If the array block was
|
||||
external to TNT, the refernce count is a NULL pointer
|
||||
regardless of how many references are made, since the
|
||||
memory is not freed by TNT.
|
||||
|
||||
|
||||
|
||||
*/
|
||||
template <class T>
|
||||
class i_refvec
|
||||
{
|
||||
|
||||
|
||||
private:
|
||||
T* data_;
|
||||
int *ref_count_;
|
||||
|
||||
|
||||
public:
|
||||
|
||||
i_refvec();
|
||||
explicit i_refvec(int n);
|
||||
inline i_refvec(T* data);
|
||||
inline i_refvec(const i_refvec &v);
|
||||
inline T* begin();
|
||||
inline const T* begin() const;
|
||||
inline T& operator[](int i);
|
||||
inline const T& operator[](int i) const;
|
||||
inline i_refvec<T> & operator=(const i_refvec<T> &V);
|
||||
void copy_(T* p, const T* q, const T* e);
|
||||
void set_(T* p, const T* b, const T* e);
|
||||
inline int ref_count() const;
|
||||
inline int is_null() const;
|
||||
inline void destroy();
|
||||
~i_refvec();
|
||||
|
||||
};
|
||||
|
||||
template <class T>
|
||||
void i_refvec<T>::copy_(T* p, const T* q, const T* e)
|
||||
{
|
||||
for (T* t=p; q<e; t++, q++)
|
||||
*t= *q;
|
||||
}
|
||||
|
||||
template <class T>
|
||||
i_refvec<T>::i_refvec() : data_(NULL), ref_count_(NULL) {}
|
||||
|
||||
/**
|
||||
In case n is 0 or negative, it does NOT call new.
|
||||
*/
|
||||
template <class T>
|
||||
i_refvec<T>::i_refvec(int n) : data_(NULL), ref_count_(NULL)
|
||||
{
|
||||
if (n >= 1)
|
||||
{
|
||||
#ifdef TNT_DEBUG
|
||||
std::cout << "new data storage.\n";
|
||||
#endif
|
||||
data_ = new T[n];
|
||||
ref_count_ = new int;
|
||||
*ref_count_ = 1;
|
||||
}
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline i_refvec<T>::i_refvec(const i_refvec<T> &V): data_(V.data_),
|
||||
ref_count_(V.ref_count_)
|
||||
{
|
||||
if (V.ref_count_ != NULL)
|
||||
(*(V.ref_count_))++;
|
||||
}
|
||||
|
||||
|
||||
template <class T>
|
||||
i_refvec<T>::i_refvec(T* data) : data_(data), ref_count_(NULL) {}
|
||||
|
||||
template <class T>
|
||||
inline T* i_refvec<T>::begin()
|
||||
{
|
||||
return data_;
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline const T& i_refvec<T>::operator[](int i) const
|
||||
{
|
||||
return data_[i];
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline T& i_refvec<T>::operator[](int i)
|
||||
{
|
||||
return data_[i];
|
||||
}
|
||||
|
||||
|
||||
template <class T>
|
||||
inline const T* i_refvec<T>::begin() const
|
||||
{
|
||||
return data_;
|
||||
}
|
||||
|
||||
|
||||
|
||||
template <class T>
|
||||
i_refvec<T> & i_refvec<T>::operator=(const i_refvec<T> &V)
|
||||
{
|
||||
if (this == &V)
|
||||
return *this;
|
||||
|
||||
|
||||
if (ref_count_ != NULL)
|
||||
{
|
||||
(*ref_count_) --;
|
||||
if ((*ref_count_) == 0)
|
||||
destroy();
|
||||
}
|
||||
|
||||
data_ = V.data_;
|
||||
ref_count_ = V.ref_count_;
|
||||
|
||||
if (V.ref_count_ != NULL)
|
||||
(*(V.ref_count_))++;
|
||||
|
||||
return *this;
|
||||
}
|
||||
|
||||
template <class T>
|
||||
void i_refvec<T>::destroy()
|
||||
{
|
||||
if (ref_count_ != NULL)
|
||||
{
|
||||
#ifdef TNT_DEBUG
|
||||
std::cout << "destorying data... \n";
|
||||
#endif
|
||||
delete ref_count_;
|
||||
|
||||
#ifdef TNT_DEBUG
|
||||
std::cout << "deleted ref_count_ ...\n";
|
||||
#endif
|
||||
if (data_ != NULL)
|
||||
delete []data_;
|
||||
#ifdef TNT_DEBUG
|
||||
std::cout << "deleted data_[] ...\n";
|
||||
#endif
|
||||
data_ = NULL;
|
||||
}
|
||||
}
|
||||
|
||||
/*
|
||||
* return 1 is vector is empty, 0 otherwise
|
||||
*
|
||||
* if is_null() is false and ref_count() is 0, then
|
||||
*
|
||||
*/
|
||||
template<class T>
|
||||
int i_refvec<T>::is_null() const
|
||||
{
|
||||
return (data_ == NULL ? 1 : 0);
|
||||
}
|
||||
|
||||
/*
|
||||
* returns -1 if data is external,
|
||||
* returns 0 if a is NULL array,
|
||||
* otherwise returns the positive number of vectors sharing
|
||||
* this data space.
|
||||
*/
|
||||
template <class T>
|
||||
int i_refvec<T>::ref_count() const
|
||||
{
|
||||
if (data_ == NULL)
|
||||
return 0;
|
||||
else
|
||||
return (ref_count_ != NULL ? *ref_count_ : -1) ;
|
||||
}
|
||||
|
||||
template <class T>
|
||||
i_refvec<T>::~i_refvec()
|
||||
{
|
||||
if (ref_count_ != NULL)
|
||||
{
|
||||
(*ref_count_)--;
|
||||
|
||||
if (*ref_count_ == 0)
|
||||
destroy();
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
} /* namespace TNT */
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
#endif
|
||||
/* TNT_I_REFVEC_H */
|
||||
|
||||
37
src/Lib/Polygon/TNT/tnt_math_utils.h
Normal file
37
src/Lib/Polygon/TNT/tnt_math_utils.h
Normal file
|
|
@ -0,0 +1,37 @@
|
|||
#ifndef MATH_UTILS_H
|
||||
#define MATH_UTILS_H
|
||||
|
||||
|
||||
/* needed for abs(), sqrt() below */
|
||||
#include <cmath>
|
||||
|
||||
|
||||
|
||||
namespace TNT
|
||||
{
|
||||
|
||||
using namespace std;
|
||||
/**
|
||||
@returns hypotenuse of real (non-complex) scalars a and b by
|
||||
avoiding underflow/overflow
|
||||
using (a * sqrt( 1 + (b/a) * (b/a))), rather than
|
||||
sqrt(a*a + b*b).
|
||||
*/
|
||||
template <class Real>
|
||||
Real hypot(const Real &a, const Real &b)
|
||||
{
|
||||
|
||||
if (a== 0)
|
||||
return abs(b);
|
||||
else
|
||||
{
|
||||
Real c = b/a;
|
||||
return abs(a) * sqrt(1 + c*c);
|
||||
}
|
||||
}
|
||||
} /* TNT namespace */
|
||||
|
||||
|
||||
|
||||
#endif
|
||||
/* MATH_UTILS_H */
|
||||
Loading…
Add table
Reference in a new issue