removed dependencies to theICTlab repo; separated GUI components; inc…

…luded bare minimum files for compiling and running TC

Algebra.h

@@ -0,0 +1,259 @@
+#ifndef __ALGEBRA_H__
+#define __ALGEBRA_H__
+
+#include <iostream>
+#include <limits>
+#include <Eigen/Eigen>
+using namespace Eigen;
+
+#define XX 0
+#define YY 1
+#define ZZ 2
+
+#define EMAX 0
+#define EMID 1
+#define EMIN 2
+static short eig_sys3d(float a[][3], float d[])	{
+	#define SIGN(a,b) ((b)<0 ? -fabs(a) : fabs(a))
+/* Adopted From Numerical Recipes in C (2nd Edition) pp. 474-475 :
+   Householder reduction of a real, symmetric matrix a[0..2][0..2]. On output,
+   a is replaced by the orthogonal matrix Q effecting the transformation.
+   d[0..2] returns the diagonal elements of the tridiagonal matrix,
+   and e[0..2] the off-diagonal elements, with e[0] = 0. Several
+   statements, as noted in comments, can be omitted if only eigenvalues are
+   to be found, in which case a contains no useful information on output.
+   Otherwise they are to be included. */
+
+    static float scale,hh,h,g,f, e[3];
+
+    // i = 2; l = 1;
+    h = 0.0;
+    scale = (float) (fabs(a[2][0])+fabs(a[2][1]));
+    if (scale == 0.0) e[2] = a[2][1];
+    else {
+      a[2][0] /= scale;
+      a[2][1] /= scale;
+      h = a[2][0]*a[2][0]+a[2][1]*a[2][1];
+      f = a[2][1];
+      g = f>=0 ? ((float) -sqrt(h)) : ((float) sqrt(h));
+      e[2] = scale*g;
+      h -= f*g;
+      a[2][1] = f-g;
+
+      a[0][2] = a[2][0]/h;
+      e[0] = (a[0][0]*a[2][0]+a[1][0]*a[2][1])/h;
+      a[1][2] = a[2][1]/h;
+      e[1] = (a[1][0]*a[2][0]+a[1][1]*a[2][1])/h;
+      f = e[0]*a[2][0]+e[1]*a[2][1];
+
+      hh = f/(h+h);
+      f = a[2][0];
+      e[0] = g = e[0]-hh*f;
+      a[0][0] -= f*e[0]+g*a[2][0];
+      f = a[2][1];
+      e[1] = g = e[1]-hh*f;
+      a[1][0] -= f*e[0]+g*a[2][0];
+      a[1][1] -= f*e[1]+g*a[2][1];
+    }
+    d[2] = h;
+
+    // i = 1; l = 0;
+    e[1] = a[1][0];
+    d[1] = 0.0;
+
+    // Next statement can be omitted if eigenvectors not wanted
+    // i=0; l=-1;
+    d[0] = a[0][0];
+    a[0][0] = 1.0;
+
+    // i=1; l=0;
+    if (d[1]) a[0][0] -= (a[1][0]*a[0][0])*a[0][1];
+    d[1] = a[1][1];
+    a[1][1] = 1.0;
+    a[0][1]=a[1][0]=0.0;
+
+    // i=2; l=1;
+    if (d[2]) {
+      g = a[2][0]*a[0][0]+a[2][1]*a[1][0];
+      a[0][0] -= g*a[0][2];
+      a[1][0] -= g*a[1][2];
+      g = a[2][0]*a[0][1]+a[2][1]*a[1][1];
+      a[0][1] -= g*a[0][2];
+      a[1][1] -= g*a[1][2];
+    }
+    d[2] = a[2][2];
+    a[2][2] = 1.0;
+    a[0][2] = a[2][0] = 0.0;
+    a[1][2] = a[2][1] = 0.0;
+
+/* Adopted From Numerical Recipes in C (2nd Edition) pp. 480-481 :
+   QL algorithm with implicit shifts, to determine the eigenvalues and
+   eigenvectors of a real, symmetric, tridiagonal matrix, or of a real,
+   symmetric matrix previously reduced by tred2.  On input, d[1..n] contains
+   the diagonal elements of the tridiagonal matrix. On output, it returns the
+   eigenvalues.  The vector e[1..n] inputs the subdiagonal elements of the
+   tridiagonal matrix, with e[1] arbitrary.  On output e is destroyed.  When
+   finding only the eigenvalues, several lines may be omitted, as noted in the
+   comments.  If the eigenvectors of a tridiagonal matrix are desired, the
+   matrix a[1..n][1..n] is input as the identity matrix.  If the eigenvectors
+   of a matrix that has been reduced by tred2 are required, then z is input as
+   the matrix output by tred2. In either case, the kth column of z returns the
+   normalized eigenvector corresponding to d[k]. */
+
+    static int m,l,iter,i,k;
+    static float s,r,p,dd,c,b;
+
+    e[0] = e[1];
+    e[1] = e[2];
+    e[2]=0.0;
+
+    for (l=0;l<=2;l++) {
+        iter=0;
+        do {
+            for (m=l;m<=1;m++) {
+                dd= (float) (fabs(d[m])+fabs(d[m+1]));
+                if (fabs(e[m])+dd == dd) break;
+            }
+            if (m != l) {
+                if (iter++ == 50) {
+                    std::cout << "Too many iterations in TQLI" << std::endl;
+                    return(0);
+                }
+                g=(d[l+1]-d[l])/(2.0f*e[l]);
+                r= (float) sqrt((g*g)+1.0);
+                g=d[m]-d[l]+e[l]/(g+(float) SIGN(r,g));
+                s=c=1.0;
+                p=0.0;
+                for (i=m-1;i>=l;i--) {
+                    f=s*e[i];
+                    b=c*e[i];
+                    if (fabs(f) >= fabs(g)) {
+                        c=g/f;
+                        r= (float) sqrt((c*c)+1.0);
+                        e[i+1]=f*r;
+                        c *= (s=1.0f/r);
+                    } else {
+                        s=f/g;
+                        r= (float) sqrt((s*s)+1.0);
+                        e[i+1]=g*r;
+                        s *= (c=1.0f/r);
+                    }
+                    g=d[i+1]-p;
+                    r=(d[i]-g)*s+2.0f*c*b;
+                    p=s*r;
+                    d[i+1]=g+p;
+                    g=c*r-b;
+                    // Next loop can be omitted if eigenvectors not wanted
+                    for (k=0;k<=2;k++) {
+                        f=a[k][i+1];
+                        a[k][i+1]=s*a[k][i]+c*f;
+                        a[k][i]=c*a[k][i]-s*f;
+                    }
+                }
+                d[l]=d[l]-p;
+                e[l]=g;
+                e[m]=0.0;
+            }
+        } while (m != l);
+    }
+
+	// sort the eigenvalues and eigenvectors
+	static float max_val;
+	static int max_index;
+	for(i=0;i<=1;i++) {
+	  max_val = d[i];
+	  max_index = i;
+	  for(k=i+1;k<=2;k++)
+		if (max_val < d[k]) { max_val = d[k]; max_index = k; }
+	  if (max_index != i) {
+		e[0] = d[i]; d[i] = d[max_index]; d[max_index] = e[0];
+		e[0] = a[0][i]; a[0][i] = a[0][max_index]; a[0][max_index] = e[0];
+		e[0] = a[1][i]; a[1][i] = a[1][max_index]; a[1][max_index] = e[0];
+		e[0] = a[2][i]; a[2][i] = a[2][max_index]; a[2][max_index] = e[0];
+	  }
+	}
+
+    return(1);
+}
+
+static void rotationMatrixAboutArbitraryAxis(Matrix3f *rotation_matrix, Vector3f const &axis,float angle, bool left_handed = true)	{
+	float c = (float) cos(angle);
+	float s = (float) sin(angle);
+	float x = axis(0);
+	float y = axis(1);
+	float z = axis(2);
+	float t = 1.0f - c;
+
+	if (!left_handed)	{
+		///Create the matrix right-handed
+		(*rotation_matrix)(0,0) = c + t*x*x;
+		(*rotation_matrix)(0,1) = t*y*x + s*z;
+		(*rotation_matrix)(0,2) = t*z*x - s*y;
+
+		(*rotation_matrix)(1,0) = t*x*y - s*z;
+		(*rotation_matrix)(1,1) = c + t*y*y;
+		(*rotation_matrix)(1,2) = t*z*y + s*x;
+
+		(*rotation_matrix)(2,0) = t*x*z + s*y;
+		(*rotation_matrix)(2,1) = t*y*z - s*x;
+		(*rotation_matrix)(2,2) = c + t*z*z;
+	}
+	else	{
+		///Create the matrix left-handed
+		(*rotation_matrix)(0,0) = c + t*x*x;
+		(*rotation_matrix)(0,1) = t*y*x - s*z;
+		(*rotation_matrix)(0,2) = t*z*x + s*y;
+
+		(*rotation_matrix)(1,0) = t*x*y + s*z;
+		(*rotation_matrix)(1,1) = c + t*y*y;
+		(*rotation_matrix)(1,2) = t*z*y - s*x;
+
+		(*rotation_matrix)(2,0) = t*x*z - s*y;
+		(*rotation_matrix)(2,1) = t*y*z + s*x;
+		(*rotation_matrix)(2,2) = c + t*z*z;
+	}
+
+	return;
+}
+
+static void eulerAnglesFromRotationMatrix(float *eulerAngles,Matrix3f const &matrix)	{
+
+	//Given the rotation matrix this method retrieves the euler angles
+	//angle a is around X axis
+	//angle b is around Y axis
+	//angle c is around Z axis
+
+	float cosb;
+
+	eulerAngles[1] = (float) -asin(matrix(2,0));	// also PI + asin(matrix(3,1))
+
+	cosb = (float) cos(eulerAngles[1]);
+
+	if (cosb > 0.0f)	{
+		eulerAngles[0] = (float) atan2(matrix(2,1),matrix(2,2));
+		eulerAngles[2] = (float) atan2(matrix(1,0),matrix(0,0));
+	}
+	else	{
+		if (cosb < 0.0f)	{
+			eulerAngles[0] = (float) atan2(-matrix(2,1),-matrix(2,2));
+			eulerAngles[2] = (float) atan2(-matrix(1,0),-matrix(0,0));
+		}
+		else	{ //cosb == 0.0f
+				  //You can write a as a function of c in this case,
+				  //there are infinitely many solutions (gimbal lock)
+				  //We just choose c = 0.0f
+				 eulerAngles[2] = 0.0f;
+				 if (matrix(2,0) == 1.0f)	{
+					 eulerAngles[1] = (float) M_PI/2.0f;
+					 eulerAngles[0] = (float) atan2(matrix(0,1),matrix(0,2));
+				 }
+				 else	{ //% (M(3,1) == -1)
+					 eulerAngles[1] = (float) -M_PI/2.0f;
+					 eulerAngles[0] = (float) atan2(-matrix(0,1),-matrix(0,2));
+				 }
+		}
+	}
+
+}
+
+#endif //__ALGEBRA_H__

BoundingBox.cpp

@@ -0,0 +1,57 @@
+////////////////////////////////////////////////////////////////////////////////////
+// Copyright © Charalambos "Charis" Poullis, [email protected]    	  //
+// This work can only be used under an exclusive license of the author.           //
+//////////////////////////////////////////////////////////////////////////////////// 
+
+#ifndef __BOUNDING_BOX_CPP__
+#define __BOUNDING_BOX_CPP__
+
+#include "BoundingBox.h"
+
+BoundingBox::BoundingBox()	{
+	min_pt = Vector3f(FLT_MAX, FLT_MAX, FLT_MAX);
+	max_pt = Vector3f(-FLT_MAX, -FLT_MAX, -FLT_MAX);
+	onb = 0x00;
+}
+
+BoundingBox::BoundingBox(Vector3f const &_min_pt, Vector3f const &_max_pt, ONB *_onb)	{
+	min_pt = _min_pt;
+	max_pt = _max_pt;
+	onb = _onb;
+}
+
+BoundingBox::~BoundingBox()	{
+	if (onb)	delete onb;
+}
+
+
+Vector3f BoundingBox::getMinPt() const	{
+	return min_pt;
+}
+
+Vector3f BoundingBox::getMaxPt() const	{
+	return max_pt;
+}
+
+ONB *BoundingBox::getONB() const	{
+	return onb;
+}
+
+void BoundingBox::setMinPt(Vector3f const _min_pt)	{
+	min_pt = _min_pt;
+	return;
+}
+
+void BoundingBox::setMaxPt(Vector3f const _max_pt)	{
+	max_pt = _max_pt;
+	return;
+}
+
+void BoundingBox::setONB(ONB *_onb)	{
+	if (onb)	delete onb;
+	onb = _onb;
+	return;
+}
+
+
+#endif

BoundingBox.h

@@ -0,0 +1,57 @@
+////////////////////////////////////////////////////////////////////////////////////
+// Copyright © Charalambos "Charis" Poullis, [email protected]    	  //
+// This work can only be used under an exclusive license of the author.           //
+//////////////////////////////////////////////////////////////////////////////////// 
+
+
+#ifndef __BOUNDING_BOX_H__
+#define __BOUNDING_BOX_H__
+
+/**		Bounding Box
+* This class defines a bounding box and related functions
+* 
+*/
+
+
+#include <float.h>
+#include "ONB.h"
+
+
+class BoundingBox	{
+	public: 
+		///Constructor
+		BoundingBox();
+		///Constructor
+		BoundingBox(Vector3f const &_min_pt, Vector3f const &_max_pt, ONB *_onb = 0x00);
+		///Destructor
+		~BoundingBox();
+
+		///Get min point
+		Vector3f getMinPt() const;
+
+		///Get max point
+		Vector3f getMaxPt() const;
+
+		///Get ONB
+		ONB *getONB()	const;
+
+		///Set min point
+		void setMinPt(Vector3f const _min_pt);
+
+		///Set max point
+		void setMaxPt(Vector3f const _max_pt);
+
+		///Set the orientation
+		void setONB(ONB *_onb);
+
+	protected:
+		///The minimum point
+		Vector3f min_pt;
+		///The maximum point
+		Vector3f max_pt;
+		///The orthonormal basis of the bbox
+		ONB *onb;
+};
+
+
+#endif