From d10991cbb002d5cd2ebe546c2b7bcb2976bfa70d Mon Sep 17 00:00:00 2001
From: Christian Schmitt <chris@ilovelinux.de>
Date: Sat, 15 Dec 2012 11:25:32 +0100
Subject: [PATCH] TNT header files were missing. They are templates and don't
 have to be in the makefile

---
 src/Lib/Polygon/CMakeLists.txt       |   7 +-
 src/Lib/Polygon/TNT/jama_qr.h        | 326 +++++++++++++++++++
 src/Lib/Polygon/TNT/tnt_array1d.h    | 301 ++++++++++++++++++
 src/Lib/Polygon/TNT/tnt_array2d.h    | 451 +++++++++++++++++++++++++++
 src/Lib/Polygon/TNT/tnt_i_refvec.h   | 243 +++++++++++++++
 src/Lib/Polygon/TNT/tnt_math_utils.h |  37 +++
 6 files changed, 1359 insertions(+), 6 deletions(-)
 create mode 100644 src/Lib/Polygon/TNT/jama_qr.h
 create mode 100644 src/Lib/Polygon/TNT/tnt_array1d.h
 create mode 100644 src/Lib/Polygon/TNT/tnt_array2d.h
 create mode 100644 src/Lib/Polygon/TNT/tnt_i_refvec.h
 create mode 100644 src/Lib/Polygon/TNT/tnt_math_utils.h

diff --git a/src/Lib/Polygon/CMakeLists.txt b/src/Lib/Polygon/CMakeLists.txt
index 3ad344f2..5ce94618 100644
--- a/src/Lib/Polygon/CMakeLists.txt
+++ b/src/Lib/Polygon/CMakeLists.txt
@@ -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
-)
\ No newline at end of file
+)
diff --git a/src/Lib/Polygon/TNT/jama_qr.h b/src/Lib/Polygon/TNT/jama_qr.h
new file mode 100644
index 00000000..98d37f53
--- /dev/null
+++ b/src/Lib/Polygon/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
+{
+
+/** 
+<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
+
diff --git a/src/Lib/Polygon/TNT/tnt_array1d.h b/src/Lib/Polygon/TNT/tnt_array1d.h
new file mode 100644
index 00000000..c4c155fe
--- /dev/null
+++ b/src/Lib/Polygon/TNT/tnt_array1d.h
@@ -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 */
+
diff --git a/src/Lib/Polygon/TNT/tnt_array2d.h b/src/Lib/Polygon/TNT/tnt_array2d.h
new file mode 100644
index 00000000..71eeee50
--- /dev/null
+++ b/src/Lib/Polygon/TNT/tnt_array2d.h
@@ -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 &lt class ArrayTwoD &gt
+  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 */
+
diff --git a/src/Lib/Polygon/TNT/tnt_i_refvec.h b/src/Lib/Polygon/TNT/tnt_i_refvec.h
new file mode 100644
index 00000000..5a67eb57
--- /dev/null
+++ b/src/Lib/Polygon/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 <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 */
+
diff --git a/src/Lib/Polygon/TNT/tnt_math_utils.h b/src/Lib/Polygon/TNT/tnt_math_utils.h
new file mode 100644
index 00000000..4d81950f
--- /dev/null
+++ b/src/Lib/Polygon/TNT/tnt_math_utils.h
@@ -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 */