@@ -43,12 +33,12 @@
// lookup node elevations for each point in the point_list. Returns
// average of all points. Doesn't modify the original list.
-double tgAverageElevation( const string &root, const string_list elev_src,
+double tgAverageElevation( const std::string &root, const string_list elev_src,
const point_list points_source );
// lookup node elevations for each point in the specified nurbs++
// matrix.
-void tgCalcElevations( const string &root, const string_list elev_src, SimpleMatrix &Pts, double average );
+void tgCalcElevations( const std::string &root, const string_list elev_src, SimpleMatrix &Pts, double average );
// clamp all elevations to the specified range
void tgClampElevations( SimpleMatrix &Pts, double center_m, double max_clamp_m );
diff --git a/src/Lib/Geometry/TNT/jama_qr.h b/src/Lib/Geometry/TNT/jama_qr.h
new file mode 100644
index 00000000..98d37f53
--- /dev/null
+++ b/src/Lib/Geometry/TNT/jama_qr.h
@@ -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
+{
+
+/**
+
+ 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.
+
+ 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).
+
+
+ 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.
+
+
+ (Adapted from JAMA, a Java Matrix Library, developed by jointly
+ by the Mathworks and NIST; see http://math.nist.gov/javanumerics/jama).
+*/
+
+template
+class QR {
+
+
+ /** Array for internal storage of decomposition.
+ @serial internal array storage.
+ */
+
+ TNT::Array2D 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 Rdiag;
+
+
+public:
+
+/**
+ Create a QR factorization object for A.
+
+ @param A rectangular (m>=n) matrix.
+*/
+ QR(const TNT::Array2D &A) /* constructor */
+ {
+ QR_ = A.copy();
+ m = A.dim1();
+ n = A.dim2();
+ Rdiag = TNT::Array1D(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 getHouseholder (void) const
+ {
+ TNT::Array2D 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 getR() const
+ {
+ TNT::Array2D 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 getQ() const
+ {
+ int i=0, j=0, k=0;
+
+ TNT::Array2D 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 solve(const TNT::Array1D &b) const
+ {
+ if (b.dim1() != m) /* arrays must be conformant */
+ return TNT::Array1D();
+
+ if ( !isFullRank() ) /* matrix is rank deficient */
+ {
+ return TNT::Array1D();
+ }
+
+ TNT::Array1D 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 x_(n);
+ for (int i=0; i solve(const TNT::Array2D &B) const
+ {
+ if (B.dim1() != m) /* arrays must be conformant */
+ return TNT::Array2D(0,0);
+
+ if ( !isFullRank() ) /* matrix is rank deficient */
+ {
+ return TNT::Array2D(0,0);
+ }
+
+ int nx = B.dim2();
+ TNT::Array2D 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 X_(n,nx);
+ for (i=0; i
+#include
+
+#ifdef TNT_BOUNDS_CHECK
+#include
+#endif
+
+
+#include "tnt_i_refvec.h"
+
+namespace TNT
+{
+
+template
+class Array1D
+{
+
+ private:
+
+ /* ... */
+ i_refvec 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 subarray(int i0, int i1);
+
+};
+
+
+
+
+template
+Array1D::Array1D() : v_(), n_(0), data_(0) {}
+
+template
+Array1D::Array1D(const Array1D &A) : v_(A.v_), n_(A.n_),
+ data_(A.data_)
+{
+#ifdef TNT_DEBUG
+ std::cout << "Created Array1D(const Array1D &A) \n";
+#endif
+
+}
+
+template
+Array1D::Array1D(int n) : v_(n), n_(n), data_(v_.begin())
+{
+#ifdef TNT_DEBUG
+ std::cout << "Created Array1D(int n) \n";
+#endif
+}
+
+template
+Array1D::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
+Array1D::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
+inline Array1D::operator T*()
+{
+ return &(v_[0]);
+}
+
+
+template
+inline Array1D::operator const T*()
+{
+ return &(v_[0]);
+}
+
+
+
+template
+inline T& Array1D::operator[](int i)
+{
+#ifdef TNT_BOUNDS_CHECK
+ assert(i>= 0);
+ assert(i < n_);
+#endif
+ return data_[i];
+}
+
+template
+inline const T& Array1D::operator[](int i) const
+{
+#ifdef TNT_BOUNDS_CHECK
+ assert(i>= 0);
+ assert(i < n_);
+#endif
+ return data_[i];
+}
+
+
+
+template
+inline T& Array1D::operator()(int i)
+{
+#ifdef TNT_BOUNDS_CHECK
+ assert(i>= 0);
+ assert(i < n_);
+#endif
+ return data_[i-1];
+}
+
+template
+inline const T& Array1D::operator()(int i) const
+{
+#ifdef TNT_BOUNDS_CHECK
+ assert(i>= 0);
+ assert(i < n_);
+#endif
+ return data_[i-1];
+}
+
+
+
+
+template
+Array1D & Array1D::operator=(const T &a)
+{
+ set_(data_, data_+n_, a);
+ return *this;
+}
+
+template
+Array1D Array1D::copy() const
+{
+ Array1D A( n_);
+ copy_(A.data_, data_, n_);
+
+ return A;
+}
+
+
+template
+Array1D & Array1D::inject(const Array1D &A)
+{
+ if (A.n_ == n_)
+ copy_(data_, A.data_, n_);
+
+ return *this;
+}
+
+
+
+
+
+template
+Array1D & Array1D::ref(const Array1D &A)
+{
+ if (this != &A)
+ {
+ v_ = A.v_; /* operator= handles the reference counting. */
+ n_ = A.n_;
+ data_ = A.data_;
+
+ }
+ return *this;
+}
+
+template
+Array1D & Array1D::operator=(const Array1D &A)
+{
+ return ref(A);
+}
+
+template
+inline int Array1D::dim1() const { return n_; }
+
+template
+inline int Array1D::dim() const { return n_; }
+
+template
+Array1D::~Array1D() {}
+
+
+/* ............................ exented interface ......................*/
+
+template
+inline int Array1D::ref_count() const
+{
+ return v_.ref_count();
+}
+
+template
+inline Array1D Array1D::subarray(int i0, int i1)
+{
+ if ((i0 >= 0) && (i1 < n_) || (i0 <= i1))
+ {
+ Array1D X(*this); /* create a new instance of this array. */
+ X.n_ = i1-i0+1;
+ X.data_ += i0;
+
+ return X;
+ }
+ else
+ {
+ return Array1D();
+ }
+}
+
+
+/* private internal functions */
+
+
+template
+void Array1D::set_(T* begin, T* end, const T& a)
+{
+ for (T* p=begin; p
+void Array1D::copy_(T* p, const T* q, int len) const
+{
+ T *end = p + len;
+ while (p
+#include
+#ifdef TNT_BOUNDS_CHECK
+#include
+#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..
+
+
+ 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.
+
+
+ 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 Array2D
+{
+
+
+ private:
+
+
+
+ Array1D data_;
+ Array1D 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,
+
+ template < class ArrayTwoD >
+ void foo(ArrayTwoD &A)
+ {
+ A::value_type first_entry = A[0][0];
+ ...
+ }
+
+*/
+ typedef T value_type;
+
+/**
+ Create a null array. This is not 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.
+
+*/
+ 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&".
+
+*/
+
+ 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).
+
+
+ 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
+Array2D::Array2D() : data_(), v_(), m_(0), n_(0) {}
+
+
+template
+Array2D::Array2D(const Array2D &A) : data_(A.data_), v_(A.v_),
+ m_(A.m_), n_(A.n_) {}
+
+
+
+
+template
+Array2D::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
+Array2D::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
+Array2D::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
+inline T* Array2D::operator[](int i)
+{
+#ifdef TNT_BOUNDS_CHECK
+ assert(i >= 0);
+ assert(i < m_);
+#endif
+
+return v_[i];
+
+}
+
+
+template
+inline const T* Array2D::operator[](int i) const
+{
+#ifdef TNT_BOUNDS_CHECK
+ assert(i >= 0);
+ assert(i < m_);
+#endif
+
+return v_[i];
+
+}
+
+template
+Array2D & Array2D::operator=(const T &a)
+{
+ /* non-optimzied, but will work with subarrays in future verions */
+
+ for (int i=0; i
+Array2D Array2D::copy() const
+{
+ Array2D A(m_, n_);
+
+ for (int i=0; i
+Array2D & Array2D::inject(const Array2D &A)
+{
+ if (A.m_ == m_ && A.n_ == n_)
+ {
+ for (int i=0; i
+Array2D & Array2D::ref(const Array2D &A)
+{
+ if (this != &A)
+ {
+ v_ = A.v_;
+ data_ = A.data_;
+ m_ = A.m_;
+ n_ = A.n_;
+
+ }
+ return *this;
+}
+
+
+
+template
+Array2D & Array2D::operator=(const Array2D &A)
+{
+ return ref(A);
+}
+
+template
+inline int Array2D::dim1() const { return m_; }
+
+template
+inline int Array2D::dim2() const { return n_; }
+
+
+template
+Array2D::~Array2D() {}
+
+
+
+
+template
+inline Array2D::operator T**()
+{
+ return &(v_[0]);
+}
+template
+inline Array2D::operator const T**() const
+{
+ return static_cast(&(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
+Array2D Array2D::subarray(int i0, int i1, int j0, int j1)
+{
+ Array2D 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(m);
+ T* p = &(data_[0]) + i0 * n_ + j0;
+ for (int i=0; i
+inline int Array2D::ref_count()
+{
+ return ref_count_data();
+}
+
+
+
+template
+inline int Array2D::ref_count_data()
+{
+ return data_.ref_count();
+}
+
+template
+inline int Array2D::ref_count_dim1()
+{
+ return v_.ref_count();
+}
+
+
+
+
+} /* namespace TNT */
+
+#endif
+/* TNT_ARRAY2D_H */
+
diff --git a/src/Lib/Geometry/TNT/tnt_i_refvec.h b/src/Lib/Geometry/TNT/tnt_i_refvec.h
new file mode 100644
index 00000000..5a67eb57
--- /dev/null
+++ b/src/Lib/Geometry/TNT/tnt_i_refvec.h
@@ -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
+#include
+
+#ifdef TNT_BOUNDS_CHECK
+#include
+#endif
+
+#ifndef NULL
+#define NULL 0
+#endif
+
+namespace TNT
+{
+/*
+ Internal representation of ref-counted array. The TNT
+ arrays all use this building block.
+
+
+ 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 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 & operator=(const i_refvec &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
+void i_refvec::copy_(T* p, const T* q, const T* e)
+{
+ for (T* t=p; q
+i_refvec::i_refvec() : data_(NULL), ref_count_(NULL) {}
+
+/**
+ In case n is 0 or negative, it does NOT call new.
+*/
+template
+i_refvec::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
+inline i_refvec::i_refvec(const i_refvec &V): data_(V.data_),
+ ref_count_(V.ref_count_)
+{
+ if (V.ref_count_ != NULL)
+ (*(V.ref_count_))++;
+}
+
+
+template
+i_refvec::i_refvec(T* data) : data_(data), ref_count_(NULL) {}
+
+template
+inline T* i_refvec::begin()
+{
+ return data_;
+}
+
+template
+inline const T& i_refvec::operator[](int i) const
+{
+ return data_[i];
+}
+
+template
+inline T& i_refvec::operator[](int i)
+{
+ return data_[i];
+}
+
+
+template
+inline const T* i_refvec::begin() const
+{
+ return data_;
+}
+
+
+
+template
+i_refvec & i_refvec::operator=(const i_refvec &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
+void i_refvec::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
+int i_refvec::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
+int i_refvec::ref_count() const
+{
+ if (data_ == NULL)
+ return 0;
+ else
+ return (ref_count_ != NULL ? *ref_count_ : -1) ;
+}
+
+template
+i_refvec::~i_refvec()
+{
+ if (ref_count_ != NULL)
+ {
+ (*ref_count_)--;
+
+ if (*ref_count_ == 0)
+ destroy();
+ }
+}
+
+
+} /* namespace TNT */
+
+
+
+
+
+#endif
+/* TNT_I_REFVEC_H */
+
diff --git a/src/Lib/Geometry/TNT/tnt_math_utils.h b/src/Lib/Geometry/TNT/tnt_math_utils.h
new file mode 100644
index 00000000..4d81950f
--- /dev/null
+++ b/src/Lib/Geometry/TNT/tnt_math_utils.h
@@ -0,0 +1,37 @@
+#ifndef MATH_UTILS_H
+#define MATH_UTILS_H
+
+
+/* needed for abs(), sqrt() below */
+#include
+
+
+
+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
+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 */