changer to the base stim/math/*, stim/grid/* and /stim/gl/* classes in order to …

…allow the new vector to function. Added general convenience push methods to the mathvec class, as well as changed a convenience constructor to a size constructor, many small fixes

changer to the base stim/math/, stim/grid/ and /stim/gl/* classes in order to …
…allow the new vector to function. Added general convenience push methods to the mathvec class, as well as changed a convenience constructor to a size constructor, many small fixes
Pavel Govyadinov
2 parents 84583942 487a9b49
Showing 14 changed files with 839 additions and 690 deletions Show diff stats
stim/cuda/cost.h
stim/gl/gl_spider.h
stim/grids/grid.h
stim/grids/image_stack.h
stim/math/mathvec.h
stim/math/matrix.h
stim/math/quaternion.h
stim/math/rect.h
stim/math/spharmonics.h
stim/math/triangle.h
stim/math/vector.h
stim/visualization/camera.h
stim/visualization/sph_harmonics.h → stim/visualization/gl_spharmonics.h
stim/visualization/obj.h
@@ -67,8 +67,8 @@ void get_diff (float *result)
 	float valIn		= tex2D(texIn, x, y)/255.0;
 	float valTemp		= Template(x);
-	//result[idx]		= abs(valIn-valTemp);
-	result[idx]		= abs(valIn);
+	result[idx]		= abs(valIn-valTemp);
+	//result[idx]		= abs(valIn);
 }
@@ -10,7 +10,7 @@
 #include "stim/gl/gl_texture.h"
 #include "stim/visualization/camera.h"
 #include "stim/gl/error.h"
-#include "stim/math/vector.h"
+#include "stim/math/mathvec.h"
 #include "stim/math/rect.h"
 #include "stim/math/matrix.h"
 #include "stim/cuda/cost.h"
@@ -64,19 +64,22 @@ class gl_spider
 		findOptimalDirection()
 		{
 			//genTemplate(dV, 0);
+			cout << "Direction Before: " << d << endl;
 			setMatrix();
 			glCallList(dList);
 			int best = getCost();
-			stim::vec<float, 4> next;
- 			next[0] = dV[best][0]*S[0]*R[0];
-			next[1] = dV[best][1]*S[1]*R[1];
-			next[2] = dV[best][2]*S[2]*R[2];
-			next[3] = 1;
+			stim::vec<float> next(
+ 			dV[best][0]*S[0]*R[0],
+			dV[best][1]*S[1]*R[1],
+			dV[best][2]*S[2]*R[2],
+			1);
+			//next = (cT*next).norm();
 			next = (cT*next).norm();
 			setPosition(	p[0]+next[0]*m[0]/stepsize,
 					p[1]+next[1]*m[0]/stepsize,
 					p[2]+next[2]*m[0]/stepsize);
 			setDirection(next[0], next[1], next[2]);
+			cout << "Direction After: " << d << endl;
 		}
 		/// Method for finding the best d for the spider.
@@ -88,16 +91,16 @@ class gl_spider
 			setMatrix();
 			glCallList(dList+1);
 			int best = getCost();
-			stim::vec<float, 4> next;
- 			next[0] = pV[best][0];
-			next[1] = pV[best][1];
-			next[2] = pV[best][2];
-			next[3] = 1;
+			stim::vec<float> next(
+ 			pV[best][0],
+			pV[best][1],
+			pV[best][2],
+			1);
 			next = cT*next;
-			std::cout << "Optimal p:"<< next << std::endl;
 			setPosition(	next[0]*S[0]*R[0],
 					next[1]*S[1]*R[1],
 					next[2]*S[2]*R[2]);
+			std::cout << "Optimal p:"<< p << std::endl;
 		}
 		/// Method for finding the best scale for the spider.
@@ -109,11 +112,11 @@ class gl_spider
 			setMatrix();
 			glCallList(dList+2);
 			int best = getCost();
-			stim::vec<float, 4> next;
- 			next[0] = mV[best][0]*S[0]*R[0];
-			next[1] = mV[best][1]*S[1]*R[1];
-			next[2] = mV[best][2]*S[2]*R[2];
-			next[3] = 1;
+			stim::vec<float> next(
+ 			mV[best][0]*S[0]*R[0],
+			mV[best][1]*S[1]*R[1],
+			mV[best][2]*S[2]*R[2],
+			1);
 			next = cT*next;
 			std::cout << "Optimal Scale:"<< next << std::endl;
 			setMagnitude(next[0]);
@@ -123,7 +126,7 @@ class gl_spider
 		branchDetection()
 		{
 			Bind();
-			positionTemplate();
+			setMatrix();
 			glCallList(dList+3);
 			int best = getCost();
@@ -151,16 +154,17 @@ class gl_spider
 		void
 		genDirectionVectors(float solidAngle = M_PI)
 		{
+			//ofstream file;
+			//file.open("dvectors.txt");
 			//Set up the vectors necessary for Rectangle creation.
 			vec<float> Y(1.0,0.0,0.0);
 			vec<float> pos(0.0,0.0,0.0);
 			vec<float> mag(1.0, 1.0, 1.0);
 			vec<float> dir(0.0, 0.0, 1.0);
-			vec<float> d_s = direction.cart2sph().norm();
 			//Set up the variable necessary for vector creation.
 			vec<float> d_s = d.cart2sph().norm();
-			vec<float> temp;
+			vec<float> temp(0,0,0);
 			int dim = (sqrt(numSamples)-1)/2;
 			float p0  	= -M_PI;
 			float dt  	= solidAngle/(2.0 * ((float)dim + 1.0));
@@ -171,7 +175,7 @@ class gl_spider
 			int idx = 0;
 			for(int i = -dim; i <= dim; i++){
 				for(int j = -dim; j <= dim; j++){
-				
+					
 					//Create linear index
 					idx = (j+dim)+(i+dim)*((dim*2)+1);	
 					temp[0] = d_s[0]; 			//rotate vector
@@ -179,6 +183,7 @@ class gl_spider
 					temp[2] = d_s[2]+dt*(float) j;
 					temp = (temp.sph2cart()).norm();	//back to cart
+			//		file << temp[0] << "," << temp[1] << "," << temp[2] << endl;
 					dV.push_back(temp);
 				 	if(cos(Y.dot(temp))< 0.087){                                                                             Y[0] = 0.0; Y[1] = 1.0;}
 					else{Y[0] = 1.0; Y[1] = 0.0;}
@@ -202,6 +207,8 @@ class gl_spider
 		void
 		genPositionVectors(float delta = 0.2)
 		{
+			ofstream file;
+			file.open("pvectors.txt");
 			//Set up the vectors necessary for Rectangle creation.
 			vec<float> Y(1.0,0.0,0.0);
 			vec<float> pos(0.0,0.0,0.0);
@@ -209,10 +216,11 @@ class gl_spider
 			vec<float> dir(0.0, 0.0, 1.0);
 			//Set up the variable necessary for vector creation.
-			vec<float> temp;
+			vec<float> temp(0,0,0);
 			int dim = (sqrt(numSamples)-1)/2;
 			stim::rect<float> samplingPlane =
-				 stim::rect<float>(m[0]*delta, p, d);
+				 stim::rect<float>(p, d);
+			samplingPlane.scale(m[0]*delta, m[0]*delta);
 			float step = 1.0/(dim);
 			//Loop over the samples, keeping the original p sample
@@ -228,6 +236,7 @@ class gl_spider
 							0.5+step*i,
 								 0.5+step*j
 										);
+					file << temp[0] << "," << temp[1] << "," << temp[2] << endl;
 					pV.push_back(temp);
                 			hor = stim::rect<float>(mag,
 						 temp, dir,
@@ -262,7 +271,7 @@ class gl_spider
 			float max 	= 1.0+delta;
 			float step	= (max-min)/(numSamples-1);
 			float factor;
-			vec<float> temp;
+			vec<float> temp(0.0,0.0,0.0);
 			glNewList(dList+2, GL_COMPILE);
 			for(int i = 0; i < numSamples; i++){
@@ -284,7 +293,6 @@ class gl_spider
 		}
 		///@param v_x x-coordinate in buffer-space,
 		///@param v_y y-coordinate in buffer-space.
-					 ilVertex2f(0.0, j*10.0);
 		///Samples the texturespace and places a sample in the provided coordinates
 		///of bufferspace.
 		void
@@ -434,7 +442,7 @@ class gl_spider
 		///Based on the p of the spider in real space (arbitrary).
 		void setMatrix()
 		{
-			stim::vec<float, 4> rot = getRotation(d);
+			stim::vec<float> rot = getRotation(d);
 			glMatrixMode(GL_TEXTURE);
 			glLoadIdentity();
@@ -514,20 +522,30 @@ class gl_spider
 		gl_spider
 		(int samples = 1089)
 		{
-			setPosition(0.0,0.0,0.0);
-			setDirection(0.0,0.0,1.0);
-			setMagnitude(1.0);
+			p.push(0.0,0.0,0.0);
+			d.push(0.0,0.0,1.0);
+			m.push(1.0, 1.0);
+			S.push(1.0,1.0,1.0);
+			R.push(1.0,1.0,1.0);
+			//setPosition(0.0,0.0,0.0);
+			//setDirection(0.0,0.0,1.0);
+			//setMagnitude(1.0);
 			numSamples = samples;
 		}
 		///temporary constructor for convenience, will be removed in further updates.	
 		gl_spider
 		(float pos_x, float pos_y, float pos_z, float dir_x, float dir_y, float dir_z,
-			float mag_x)
-		{
-			setPosition(pos_x, pos_y, pos_z);
-			setDirection(dir_x, dir_y, dir_z);
-			setMagnitude(mag_x);
+			float mag_x, int numSamples = 1089)
+		{
+			p.push(pos_x, pos_y, pos_z);
+			d.push(dir_x, dir_y, dir_z);
+			m.push(mag_x, mag_x);
+			S.push(1.0,1.0,1.0);
+			R.push(1.0,1.0,1.0);
+			//setPosition(pos_x, pos_y, pos_z);
+			//setDirection(dir_x, dir_y, dir_z);
+			//setMagnitude(mag_x);
 		}
@@ -551,7 +569,7 @@ class gl_spider
 			GenerateFBO(20, numSamples*10);
 			setDims(0.6, 0.6, 1.0);
 			setSize(512.0, 512.0, 426.0);
-
+			setMatrix();
 			dList = glGenLists(3);
 			glListBase(dList);
 			Bind();
@@ -664,10 +682,10 @@ class gl_spider
 		///@param dir, the vector to which we are rotating
 		///given a vector to align to, finds the required
 		///axis and angle for glRotatef
-		stim::vec<float, 4>
+		stim::vec<float>
 		getRotation(stim::vec<float> dir)
 		{
-			stim::vec<float, 4> out;
+			stim::vec<float> out(0,0,0,0);;
 			stim::vec<float> from(0,0,1);
 			out[0] = acos(from.dot(dir))*M_PI/180;
 			if(out[0] < 0.0001){
@@ -718,7 +736,7 @@ class gl_spider
 			findOptimalDirection();
 			findOptimalPosition();
 			findOptimalScale();
-			branchDetection();
+//			branchDetection();
 			Unbind();
 		}
@@ -7,7 +7,7 @@
 #include <fstream>
 #include <cstdarg>
-#include "../math/vector.h"
+#include "../math/mathvec.h"
 namespace stim{
@@ -20,8 +20,8 @@ class grid{
 protected:
-	stim::vec<unsigned long, D> R;		//elements in each dimension
-	stim::vec<float, D> S;
+	stim::vec<unsigned long> R;		//elements in each dimension
+	stim::vec<float> S;
 	T* ptr;									//pointer to the data (on the GPU or CPU)
 	///Return the total number of values in the binary file
@@ -56,7 +56,7 @@ public:
 	///Constructor used to specify the grid size as a vector
 	/// @param _R is a vector describing the grid resolution
-	grid( stim::vec<unsigned long, D> _R){
+	grid( stim::vec<unsigned long> _R){
 		//set the grid resolution
 		R = _R;
@@ -65,7 +65,7 @@ public:
 	}
 	void
-	setDim(stim::vec<float, D> s)
+	setDim(stim::vec<float> s)
 	{
 		S = s;
 	}
@@ -106,7 +106,7 @@ public:
 	/// @param filename is the name of the file containing the binary data
 	/// @param S is the size of the binary file along each dimension
 	/// @param header is the size of the header in bytes
-	void read(std::string filename, stim::vec<unsigned long, D> S, unsigned long header = 0){
+	void read(std::string filename, stim::vec<unsigned long> S, unsigned long header = 0){
 		R = S;	//set the sample resolution
@@ -177,7 +177,7 @@ public:
 		result<<"]"<<std::endl;
 		//calculate the number of values to output
-		unsigned long nV = min(R[0], (unsigned long)10);
+		unsigned long nV = min((unsigned long long)R[0], (unsigned long long)10);
 		for(unsigned long v = 0; v<nV; v++){
 			result<<ptr[v];
@@ -51,10 +51,10 @@ public:
 		stim::image<T> I(file_list[0].str());
 		//set the image resolution and number of channels
-		R[0] = I.channels();
-		R[1] = I.width();
-		R[2] = I.height();
-		R[3] = file_list.size();
+		R.push(I.channels());
+		R.push(I.width());
+		R.push(I.height());
+		R.push(file_list.size());
 		//allocate storage space
 		ptr = (T*)malloc(sizeof(T) * samples());
+#ifndef RTS_VECTOR_H
+#define RTS_VECTOR_H
+
+#include <iostream>
+#include <cmath>
+#include <sstream>
+#include <vector>
+#include "../cuda/callable.h"
+
+namespace stim
+{
+
+
+
+template <class T>
+struct vec : public std::vector<T>
+{
+	using std::vector<T>::size;
+	using std::vector<T>::at;
+	using std::vector<T>::resize;
+	using std::vector<T>::push_back;
+	vec(){
+
+	}
+
+//	//efficiency constructors, makes construction easier for 1D-4D vectors
+	vec(T rhs)
+	{
+		resize(1, rhs);
+//		cerr << " Created a vector " << endl;
+//		//for(int i=0; i<N; i++)
+//		//	v[i] = rhs;
+	}
+
+//	vec(int s)
+//	{
+//		resize(s,0);
+//	}	
+
+	vec(T x, T y)
+	{
+		push_back(x);
+		push_back(y);
+	}
+	vec(T x, T y, T z)
+	{
+		push_back(x);
+		push_back(y);
+		push_back(z);
+	}
+	vec(T x, T y, T z, T w)
+	{
+		push_back(x);
+		push_back(y);
+		push_back(z);
+		push_back(w);
+	}
+
+	
+	
+	//copy constructor
+	vec( const vec<T>& other){
+		unsigned int N = other.size();
+		for(int i=0; i<N; i++)
+			push_back(other[i]);
+	}
+
+	vec<T> push(T x)
+	{
+		push_back(x);
+		return *this;
+	}
+	
+	vec<T> push(T x, T y)
+	{
+		push_back(x);
+		push_back(y);
+		return *this;
+	}
+	vec<T> push(T x, T y, T z)
+	{
+		push_back(x);
+		push_back(y);
+		push_back(z);
+		return *this;
+	}
+	vec<T> push(T x, T y, T z, T w)
+	{
+		push_back(x);
+		push_back(y);
+		push_back(z);
+		push_back(w);
+		return *this;
+	}
+	template< typename U >
+	operator vec<U>(){
+		unsigned int N = size();
+		vec<U> result;
+		for(int i=0; i<N; i++)
+			result.push_back(at(i));
+
+		return result;
+	}
+
+	//template<class U>
+	//friend vec<U, N>::operator vec<T, N>();
+
+	T len() const
+	{
+		unsigned int N = size();
+
+        //compute and return the vector length
+        T sum_sq = (T)0;
+        for(int i=0; i<N; i++)
+        {
+            sum_sq += pow( at(i), 2 );
+        }
+        return sqrt(sum_sq);
+
+	}
+
+	vec<T> cart2sph() const
+	{
+		//convert the vector from cartesian to spherical coordinates
+		//x, y, z -> r, theta, phi (where theta = 0 to 2*pi)
+
+		vec<T> sph;
+		sph.push_back(std::sqrt(at(0)*at(0) + at(1)*at(1) + at(2)*at(2)));
+		sph.push_back(std::atan2(at(1), at(0)));
+
+		if(sph[0] == 0)
+			sph.push_back(0);
+		else
+			sph.push_back(std::acos(at(2) / sph[0]));
+
+		return sph;
+	}
+
+	vec<T> sph2cart() const
+	{
+		//convert the vector from cartesian to spherical coordinates
+		//r, theta, phi -> x, y, z (where theta = 0 to 2*pi)
+
+		vec<T> cart;
+		cart.push_back(at(0) * std::cos(at(1)) * std::sin(at(2)));
+		cart.push_back(at(0) * std::sin(at(1)) * std::sin(at(2)));
+		cart.push_back(at(0) * std::cos(at(2)));
+
+		return cart;
+	}
+
+	vec<T> norm() const
+	{
+		unsigned int N = size();
+
+        //compute and return the unit vector
+        vec<T> result;
+
+        //compute the vector length
+        T l = len();
+
+        //normalize
+        for(int i=0; i<N; i++)
+        {
+            result.push_back(at(i) / l);
+        }
+
+        return result;
+	}
+
+	vec<T> cross(const vec<T> rhs) const
+	{
+		vec<T> result;
+
+		//compute the cross product (only valid for 3D vectors)
+		result.push_back(at(1) * rhs[2] - at(2) * rhs[1]);
+		result.push_back(at(2) * rhs[0] - at(0) * rhs[2]);
+		result.push_back(at(0) * rhs[1] - at(1) * rhs[0]);
+
+		return result;
+	}
+
+    T dot(vec<T> rhs) const
+    {
+        T result = (T)0;
+	unsigned int N = size();
+        for(int i=0; i<N; i++)
+            result += at(i) * rhs[i];
+
+        return result;
+
+    }
+
+	//arithmetic
+	vec<T> operator+(vec<T> rhs) const
+	{
+		vec<T> result;
+
+		unsigned int N = size();
+		for(int i=0; i<N; i++)
+		    result.push_back(at(i) + rhs[i]);
+
+		return result;
+	}
+	vec<T> operator+(T rhs) const
+	{
+		vec<T> result;
+		unsigned int N = size();
+		for(int i=0; i<N; i++)
+		    result.push_back(at(i) + rhs);
+
+		return result;
+	}
+	vec<T> operator-(vec<T> rhs) const
+	{
+        vec<T> result;
+
+		unsigned int N = size();
+        for(int i=0; i<N; i++)
+            result.push_back(at(i) - rhs[i]);
+
+        return result;
+	}
+	vec<T> operator*(T rhs) const
+	{
+        vec<T> result;
+
+		unsigned int N = size();
+        for(int i=0; i<N; i++)
+            result.push_back(at(i) * rhs);
+
+        return result;
+	}
+	vec<T> operator/(T rhs) const
+	{
+        vec<T> result;
+
+		unsigned int N = size();
+        for(int i=0; i<N; i++)
+            result.push_back(at(i) / rhs);
+
+        return result;
+	}
+	vec<T> operator*=(T rhs){
+
+		unsigned int N = size();
+		for(int i=0; i<N; i++)
+			at(i) = at(i) * rhs;
+		return *this;
+	}
+	vec<T> operator+=(vec<T> rhs){
+		unsigned int N = size();
+		for(int i=0; i<N; i++)
+			at(i) += rhs[i];
+		return *this;
+	}
+	vec<T> & operator=(T rhs){
+
+		unsigned int N = size();
+		for(int i=0; i<N; i++)
+			at(i) = rhs;
+		return *this;
+	}
+
+	template<typename Y>
+	vec<T> & operator=(vec<Y> rhs){
+		unsigned int N = rhs.size();
+		resize(N);
+
+		for(int i=0; i<N; i++)
+			at(i) = rhs[i];
+		return *this;
+	}
+	//unary minus
+	vec<T> operator-() const{
+		vec<T> r;
+
+		//negate the vector
+		unsigned int N = size();
+		for(int i=0; i<N; i++)
+		    r.push_back(-at(i));
+
+		return r;
+	}
+
+
+	std::string str() const
+	{
+		std::stringstream ss;
+
+		unsigned int N = size();
+
+		ss<<"[";
+		for(int i=0; i<N; i++)
+		{
+			ss<<at(i);
+			if(i != N-1)
+				ss<<", ";
+		}
+		ss<<"]";
+
+		return ss.str();
+	}
+
+};
+
+
+}	//end namespace rts
+
+template <typename T>
+std::ostream& operator<<(std::ostream& os, stim::vec<T> v)
+{
+    os<<v.str();
+    return os;
+}
+
+
+
+template <typename T>
+stim::vec<T> operator*(T lhs, stim::vec<T> rhs)
+{
+	stim::vec<T> r;
+
+    return rhs * lhs;
+}
+
+//#if __GNUC__ > 3 && __GNUC_MINOR__ > 7
+//template<class T, int N> using rtsVector = rts::vector<T, N>;
+//#endif
+
+#endif
@@ -4,7 +4,7 @@
 //#include "rts/vector.h"
 #include <string.h>
 #include <iostream>
-#include "vector.h"
+#include "mathvec.h"
 #include "../cuda/callable.h"
 namespace stim{
@@ -52,9 +52,12 @@ struct matrix
 	template<typename Y>
-	CUDA_CALLABLE vec<Y, N> operator*(vec<Y, N> rhs)
+	CUDA_CALLABLE vec<Y> operator*(vec<Y> rhs)
 	{
-		vec<Y, N> result;
+		unsigned int N = rhs.size();
+
+		vec<Y> result;
+		result.resize(N);
 		for(int r=0; r<N; r++)
 			for(int c=0; c<N; c++)
@@ -26,13 +26,13 @@ public:
 	CUDA_CALLABLE void CreateRotation(T theta, T ux, T uy, T uz){
-		vec<T, 3> u(ux, uy, uz);
+		vec<T> u(ux, uy, uz);
 		CreateRotation(theta, u);		
 	}
-	CUDA_CALLABLE void CreateRotation(T theta, vec<T, 3> u){
+	CUDA_CALLABLE void CreateRotation(T theta, vec<T> u){
-		vec<T, 3> u_hat = u.norm();
+		vec<T> u_hat = u.norm();
 		//assign the given Euler rotation to this quaternion
 		w = (T)cos(theta/2);
@@ -41,7 +41,7 @@ public:
 		z = u_hat[2]*(T)sin(theta/2);
 	}
-	CUDA_CALLABLE void CreateRotation(vec<T, 3> from, vec<T, 3> to){
+	CUDA_CALLABLE void CreateRotation(vec<T> from, vec<T> to){
 		vec<T> r = from.cross(to);			//compute the rotation vector
 		T theta = asin(r.len());				//compute the angle of the rotation about r
@@ -3,7 +3,7 @@
 //enable CUDA_CALLABLE macro
 #include "../cuda/callable.h"
-#include "../math/vector.h"
+#include "../math/mathvec.h"
 #include "../math/triangle.h"
 #include "../math/quaternion.h"
 #include <iostream>
@@ -13,7 +13,7 @@
 namespace stim{
 //template for a rectangle class in ND space
-template <class T, int N = 3>
+template <class T>
 struct rect
 {
 	/*
@@ -28,9 +28,9 @@ struct rect
 private:
-	stim::vec<T, N> C;
-	stim::vec<T, N> X;
-	stim::vec<T, N> Y;
+	stim::vec<T> C;
+	stim::vec<T> X;
+	stim::vec<T> Y;
 	CUDA_CALLABLE void scale(T factor){
 		X *= factor;
@@ -38,10 +38,10 @@ private:
 	}
-	CUDA_CALLABLE void normal(vec<T, N> n){		//orient the rectangle along the specified normal
+	CUDA_CALLABLE void normal(vec<T> n){		//orient the rectangle along the specified normal
 		n = n.norm();								//normalize, just in case
-		vec<T, N> n_current = X.cross(Y).norm();	//compute the current normal
+		vec<T> n_current = X.cross(Y).norm();	//compute the current normal
 		quaternion<T> q;							//create a quaternion
 		q.CreateRotation(n_current, n);				//initialize a rotation from n_current to n
@@ -51,9 +51,9 @@ private:
 	}
 	CUDA_CALLABLE void init(){
-		C = vec<T, N>(0, 0, 0);
-		X = vec<T, N>(1, 0, 0);
-		Y = vec<T, N>(0, 1, 0);
+		C = vec<T>(0, 0, 0);
+		X = vec<T>(1, 0, 0);
+		Y = vec<T>(0, 1, 0);
 	}
 public:
@@ -69,30 +69,22 @@ public:
 		C[2] = z_pos;
 	}
-	//create a rectangle from a center point, normal, and size
-	CUDA_CALLABLE rect(T size, vec<T, N> c, vec<T, N> n = vec<T, N>(0, 0, 1)){
-		init();			//start with the default setting
-		C = c;
-		scale(size);	//scale the rectangle
-		normal(n);		//orient
-	}
 	//create a rectangle from a center point, normal, and size
-	CUDA_CALLABLE rect(vec<T,2> size, vec<T, N> c, vec<T, N> n = vec<T, N>(0, 0, 1)){
+	CUDA_CALLABLE rect(vec<T> c, vec<T> n = vec<T>(0, 0, 1)){
 		init();			//start with the default setting
 		C = c;
-		scale(size);	//scale the rectangle
 		normal(n);		//orient
 	}
-	CUDA_CALLABLE rect(vec<T, N> center, vec<T, N> directionX, vec<T, N> directionY )
+	CUDA_CALLABLE rect(vec<T> center, vec<T> directionX, vec<T> directionY )
 	{
 		C = center;
 		X = directionX;
 		Y = directionY;
 	}
-	CUDA_CALLABLE rect(T size, vec<T, N> center, vec<T, N> directionX, vec<T, N> directionY )
+	CUDA_CALLABLE rect(T size, vec<T> center, vec<T> directionX, vec<T> directionY )
 	{	
 		C = center;
 		X = directionX;
@@ -100,7 +92,7 @@ public:
 		scale(size);
 	}
-	CUDA_CALLABLE rect(vec<T,N> size, vec<T, N> center, vec<T, N> directionX, vec<T, N> directionY )
+	CUDA_CALLABLE rect(vec<T> size, vec<T> center, vec<T> directionX, vec<T> directionY )
 	{	
 		C = center;
 		X = directionX;
@@ -114,7 +106,7 @@ public:
 	}
 	//boolean comparison
-	bool operator==(const rect<T, N> & rhs)
+	bool operator==(const rect<T> & rhs)
 	{
 		if(C == rhs.C && X == rhs.X && Y == rhs.Y)
 			return true;
@@ -125,24 +117,24 @@ public:
 	/*******************************************
 	Return the normal for the rect
 	*******************************************/
-	CUDA_CALLABLE stim::vec<T, N> n()
+	CUDA_CALLABLE stim::vec<T> n()
 	{
         return (X.cross(Y)).norm();
 	}
 	//get the world space value given the planar coordinates a, b in [0, 1]
-	CUDA_CALLABLE stim::vec<T, N> p(T a, T b)
+	CUDA_CALLABLE stim::vec<T> p(T a, T b)
 	{
-		stim::vec<T, N> result;
+		stim::vec<T> result;
 		//given the two parameters a, b = [0 1], returns the position in world space
-		vec<T, N> A = C - X * (T)0.5 - Y * (T)0.5;
+		vec<T> A = C - X * (T)0.5 - Y * (T)0.5;
 		result = A + X * a + Y * b;
 		return result;
 	}
 	//parenthesis operator returns the world space given rectangular coordinates a and b in [0 1]
-	CUDA_CALLABLE stim::vec<T, N> operator()(T a, T b)
+	CUDA_CALLABLE stim::vec<T> operator()(T a, T b)
 	{
 		return p(a, b);
 	}
@@ -150,7 +142,7 @@ public:
 	std::string str()
 	{
 		std::stringstream ss;
-		vec<T, N> A = C - X * (T)0.5 - Y * (T)0.5;
+		vec<T> A = C - X * (T)0.5 - Y * (T)0.5;
 		ss<<std::left<<"B="<<std::setfill('-')<<std::setw(20)<<A + Y<<">"<<"C="<<A + Y + X<<std::endl;
 		ss<<std::setfill(' ')<<std::setw(23)<<"|"<<"|"<<std::endl<<std::setw(23)<<"|"<<"|"<<std::endl;
 		ss<<std::left<<"A="<<std::setfill('-')<<std::setw(20)<<A<<">"<<"D="<<A + X;
@@ -160,12 +152,12 @@ public:
 	}
 	//scales the rectangle by a value rhs
-	CUDA_CALLABLE rect<T, N> operator*(T rhs)
+	CUDA_CALLABLE rect<T> operator*(T rhs)
 	{
 		//scales the plane by a scalar value
 		//create the new rectangle
-		rect<T, N> result = *this;
+		rect<T> result = *this;
 		result.scale(rhs);
 		return result;
@@ -173,15 +165,15 @@ public:
 	}
 	//computes the distance between the specified point and this rectangle
-	CUDA_CALLABLE T dist(vec<T, N> p)
+	CUDA_CALLABLE T dist(vec<T> p)
 	{
         //compute the distance between a point and this rect
-		vec<T, N> A = C - X * (T)0.5 - Y * (T)0.5;
+		vec<T> A = C - X * (T)0.5 - Y * (T)0.5;
         //first break the rect up into two triangles
-        triangle<T, N> T0(A, A+X, A+Y);
-        triangle<T, N> T1(A+X+Y, A+X, A+Y);
+        triangle<T> T0(A, A+X, A+Y);
+        triangle<T> T1(A+X+Y, A+X, A+Y);
         T d0 = T0.dist(p);
@@ -193,9 +185,9 @@ public:
             return d1;
 	}
-	CUDA_CALLABLE T dist_max(vec<T, N> p)
+	CUDA_CALLABLE T dist_max(vec<T> p)
 	{
-		vec<T, N> A = C - X * (T)0.5 - Y * (T)0.5;
+		vec<T> A = C - X * (T)0.5 - Y * (T)0.5;
         T da = (A - p).len();
         T db = (A+X - p).len();
         T dc = (A+Y - p).len();
@@ -208,7 +200,7 @@ public:
 }	//end namespace rts
 template <typename T, int N>
-std::ostream& operator<<(std::ostream& os, stim::rect<T, N> R)
+std::ostream& operator<<(std::ostream& os, stim::rect<T> R)
 {
     os<<R.str();
     return os;
+#ifndef STIM_SPH_HARMONICS
+#define STIM_SPH_HARMONICS
+
+
+#include <boost/math/special_functions/spherical_harmonic.hpp>
+#include <vector>
+
+#define PI 3.14159
+#define WIRE_SCALE 1.001
+namespace stim{
+
+template<class T>
+class spharmonics{
+
+protected:
+
+	std::vector<T> C;	//list of SH coefficients
+
+	unsigned int mcN;	//number of Monte-Carlo samples
+
+	//calculate the value of the SH basis function (l, m) at (theta, phi)
+		//here, theta = [0, PI], phi = [0, 2*PI]
+	double SH(int l, int m, double theta, double phi){
+		return boost::math::spherical_harmonic_r(l, m, phi, theta);
+	}
+
+	unsigned int coeff_1d(unsigned int l, int m){
+		return pow(l + 1, 2) - (l - m) - 1;
+	}
+
+	
+
+
+public:
+
+	void push(double c){
+		C.push_back(c);
+	}
+
+	void resize(unsigned int n){
+		C.resize(n);
+	}
+
+	void setc(unsigned int l, int m, T value){
+		unsigned int c = coeff_1d(l, m);
+		C[c] = value;
+	}
+
+	void setc(unsigned int c, T value){
+		C[c] = value;
+	}
+
+	/// Initialize Monte-Carlo sampling of a function using N spherical harmonics coefficients
+
+	/// @param N is the number of spherical harmonics coefficients used to represent the user function
+	void mcBegin(unsigned int coefficients){
+		C.resize(coefficients, 0);
+		mcN = 0;
+	}
+
+	void mcBegin(unsigned int l, int m){
+		unsigned int c = pow(l + 1, 2) - (l - m);
+		mcBegin(c);
+	}
+
+	void mcSample(double theta, double phi, double val){
+
+		int l, m;
+		double sh;
+
+		l = m = 0;
+		for(unsigned int i = 0; i < C.size(); i++){
+
+			sh = SH(l, m, theta, phi);
+			C[i] += sh * val;
+
+			m++;			//increment m
+
+			//if we're in a new tier, increment l and set m = -l
+			if(m > l){		
+				l++;
+				m = -l;
+			}
+		}	//end for all coefficients
+
+		//increment the number of samples
+		mcN++;
+
+	}	//end mcSample()
+
+	void mcEnd(){
+
+		//divide all coefficients by the number of samples
+		for(unsigned int i = 0; i < C.size(); i++)
+			C[i] /= mcN;
+	}
+
+	/// Generates a PDF describing the probability distribution of points on a spherical surface
+
+	/// @param sph_pts is a list of points in spherical coordinates (theta, phi) where theta = [0, 2pi] and phi = [0, pi]
+	/// @param l is the maximum degree of the spherical harmonic function
+	/// @param m is the maximum order
+	void pdf(std::vector<stim::vec<double>> sph_pts, unsigned int l, int m){
+			
+		mcBegin( l, m );		//begin spherical harmonic sampling
+
+		unsigned int nP = sph_pts.size();
+
+		for(unsigned int p = 0; p < nP; p++){
+			mcSample(sph_pts[p][1], sph_pts[p][2], 1.0);
+		}
+
+		mcEnd();
+	}
+
+	std::string str(){
+
+		std::stringstream ss;
+
+		int l, m;
+		l = m = 0;
+		for(unsigned int i = 0; i < C.size(); i++){
+				
+			ss<<C[i]<<'\t';
+
+			m++;			//increment m
+
+			//if we're in a new tier, increment l and set m = -l
+			if(m > l){
+				l++;
+				m = -l;
+
+				ss<<std::endl;
+					
+			}
+		}
+
+		return ss.str();
+
+
+	}
+
+	/// Returns the value of the function at the coordinate (theta, phi)
+
+	/// @param theta = [0, 2pi]
+	/// @param phi = [0, pi]
+	double operator()(double theta, double phi){
+
+		double fx = 0;
+
+		int l = 0;
+		int m = 0;
+		for(unsigned int i = 0; i < C.size(); i++){
+			fx += C[i] * SH(l, m, theta, phi);
+			m++;
+			if(m > l){
+				l++;
+				m = -l;					
+			}
+
+		}
+
+		return fx;
+	}
+
+};		//end class sph_harmonics
+
+
+
+
+}
+
+
+#endif
@@ -3,12 +3,12 @@
 //enable CUDA_CALLABLE macro
 #include "../cuda/callable.h"
-#include "../math/vector.h"
+#include "../math/mathvec.h"
 #include <iostream>
 namespace stim{
-template <class T, int N=3>
+template <class T>
 struct triangle
 {
     /*
@@ -23,15 +23,15 @@ struct triangle
     */
     private:
-    vec<T, N> A;
-    vec<T, N> B;
-    vec<T, N> C;
+    vec<T> A;
+    vec<T> B;
+    vec<T> C;
-    CUDA_CALLABLE vec<T, N> _p(T s, T t)
+    CUDA_CALLABLE vec<T> _p(T s, T t)
     {
         //This function returns the point specified by p = A + s(B-A) + t(C-A)
-        vec<T, N> E0 = B-A;
-        vec<T, N> E1 = C-A;
+        vec<T> E0 = B-A;
+        vec<T> E1 = C-A;
         return A + s*E0 + t*E1;
     }
@@ -46,26 +46,26 @@ struct triangle
 	}
-	CUDA_CALLABLE triangle(vec<T, N> a, vec<T, N> b, vec<T, N> c)
+	CUDA_CALLABLE triangle(vec<T> a, vec<T> b, vec<T> c)
 	{
 		A = a;
 		B = b;
 		C = c;
 	}
-	CUDA_CALLABLE stim::vec<T, N> operator()(T s, T t)
+	CUDA_CALLABLE stim::vec<T> operator()(T s, T t)
 	{
         return _p(s, t);
 	}
-	CUDA_CALLABLE vec<T, N> nearest(vec<T, N> p)
+	CUDA_CALLABLE vec<T> nearest(vec<T> p)
 	{
         //comptue the distance between a point and this triangle
         //  This code is adapted from: http://www.geometrictools.com/Documentation/DistancePoint3Triangle3.pdf
-        vec<T, N> E0 = B-A;
-        vec<T, N> E1 = C-A;
-        vec<T, N> D = A - p;
+        vec<T> E0 = B-A;
+        vec<T> E1 = C-A;
+        vec<T> D = A - p;
         T a = E0.dot(E0);
         T b = E0.dot(E1);
@@ -175,9 +175,9 @@ struct triangle
 	}
-	CUDA_CALLABLE T dist(vec<T, N> p)
+	CUDA_CALLABLE T dist(vec<T> p)
 	{
-        vec<T, N> n = nearest(p);
+        vec<T> n = nearest(p);
         return (p - n).len();
 	}
-#ifndef RTS_VECTOR_H
-#define RTS_VECTOR_H
-
-#include <iostream>
-#include <cmath>
-#include <sstream>
-//#include "rts/math/point.h"
-#include "../cuda/callable.h"
-
-
-namespace stim
-{
-
-
-
-template <class T, int N=3>
-struct vec
-{
-	T v[N];
-
-	CUDA_CALLABLE vec()
-	{
-		//memset(v, 0, sizeof(T) * N);
-		for(int i=0; i<N; i++)
-			v[i] = 0;
-	}
-
-	//efficiency constructor, makes construction easier for 1D-4D vectors
-	CUDA_CALLABLE vec(T rhs)
-	{
-		for(int i=0; i<N; i++)
-			v[i] = rhs;
-	}
-	CUDA_CALLABLE vec(T x, T y)
-	{
-		v[0] = x;
-		v[1] = y;
-	}
-	CUDA_CALLABLE vec(T x, T y, T z)
-	{
-		v[0] = x;
-		v[1] = y;
-		v[2] = z;
-	}
-	CUDA_CALLABLE vec(T x, T y, T z, T w)
-	{
-		v[0] = x;
-		v[1] = y;
-		v[2] = z;
-		v[3] = w;
-	}
-
-	//copy constructor
-	CUDA_CALLABLE vec( const vec<T, N>& other){
-		for(int i=0; i<N; i++)
-			v[i] = other.v[i];
-	}
-
-	
-	template< typename U >
-	CUDA_CALLABLE operator vec<U, N>(){
-		vec<U, N> result;
-		for(int i=0; i<N; i++)
-			result.v[i] = v[i];
-
-		return result;
-	}
-
-	//template<class U>
-	//friend vec<U, N>::operator vec<T, N>();
-
-	CUDA_CALLABLE T len() const
-	{
-        //compute and return the vector length
-        T sum_sq = (T)0;
-        for(int i=0; i<N; i++)
-        {
-            sum_sq += v[i] * v[i];
-        }
-        return sqrt(sum_sq);
-
-	}
-
-	CUDA_CALLABLE vec<T, N> cart2sph() const
-	{
-		//convert the vector from cartesian to spherical coordinates
-		//x, y, z -> r, theta, phi (where theta = 0 to 2*pi)
-
-		vec<T, N> sph;
-		sph[0] = std::sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
-		sph[1] = std::atan2(v[1], v[0]);
-		sph[2] = std::acos(v[2] / sph[0]);
-
-		return sph;
-	}
-
-	CUDA_CALLABLE vec<T, N> sph2cart() const
-	{
-		//convert the vector from cartesian to spherical coordinates
-		//r, theta, phi -> x, y, z (where theta = 0 to 2*pi)
-
-		vec<T, N> cart;
-		cart[0] = v[0] * std::cos(v[1]) * std::sin(v[2]);
-		cart[1] = v[0] * std::sin(v[1]) * std::sin(v[2]);
-		cart[2] = v[0] * std::cos(v[2]);
-
-		return cart;
-	}
-
-	CUDA_CALLABLE vec<T, N> norm() const
-	{
-        //compute and return the vector norm
-        vec<T, N> result;
-
-        //compute the vector length
-        T l = len();
-
-        //normalize
-        for(int i=0; i<N; i++)
-        {
-            result.v[i] = v[i] / l;
-        }
-
-        return result;
-	}
-
-	CUDA_CALLABLE vec<T, 3> cross(const vec<T, 3> rhs) const
-	{
-		vec<T, 3> result;
-
-		//compute the cross product (only valid for 3D vectors)
-		result[0] = v[1] * rhs.v[2] - v[2] * rhs.v[1];
-		result[1] = v[2] * rhs.v[0] - v[0] * rhs.v[2];
-		result[2] = v[0] * rhs.v[1] - v[1] * rhs.v[0];
-
-		return result;
-	}
-
-    CUDA_CALLABLE T dot(vec<T, N> rhs) const
-    {
-        T result = (T)0;
-
-        for(int i=0; i<N; i++)
-            result += v[i] * rhs.v[i];
-
-        return result;
-
-    }
-
-	//arithmetic
-	CUDA_CALLABLE vec<T, N> operator+(vec<T, N> rhs) const
-	{
-		vec<T, N> result;
-
-		for(int i=0; i<N; i++)
-		    result.v[i] = v[i] + rhs.v[i];
-
-		return result;
-	}
-	CUDA_CALLABLE vec<T, N> operator+(T rhs) const
-	{
-		vec<T, N> result;
-
-		for(int i=0; i<N; i++)
-		    result.v[i] = v[i] + rhs;
-
-		return result;
-	}
-	CUDA_CALLABLE vec<T, N> operator-(vec<T, N> rhs) const
-	{
-        vec<T, N> result;
-
-        for(int i=0; i<N; i++)
-            result.v[i] = v[i] - rhs.v[i];
-
-        return result;
-	}
-	CUDA_CALLABLE vec<T, N> operator*(T rhs) const
-	{
-        vec<T, N> result;
-
-        for(int i=0; i<N; i++)
-            result.v[i] = v[i] * rhs;
-
-        return result;
-	}
-	CUDA_CALLABLE vec<T, N> operator/(T rhs) const
-	{
-        vec<T, N> result;
-
-        for(int i=0; i<N; i++)
-            result.v[i] = v[i] / rhs;
-
-        return result;
-	}
-	CUDA_CALLABLE vec<T, N> operator*=(T rhs){
-		for(int i=0; i<N; i++)
-			v[i] = v[i] * rhs;
-		return *this;
-	}
-	CUDA_CALLABLE vec<T, N> & operator=(T rhs){
-		for(int i=0; i<N; i++)
-			v[i] = rhs;
-		return *this;
-	}
-	template<typename Y>
-	CUDA_CALLABLE vec<T, N> & operator=(vec<Y, N> rhs){
-		for(int i=0; i<N; i++)
-			v[i] = rhs.v[i];
-		return *this;
-	}
-	//unary minus
-	CUDA_CALLABLE vec<T, N> operator-() const{
-		vec<T, N> r;
-
-		//negate the vector
-		for(int i=0; i<N; i++)
-		    r.v[i] = -v[i];
-
-		return r;
-	}
-
-	CUDA_CALLABLE bool operator==(vec<T, N> rhs) const
-	{
-        if ( (rhs.v[0] == v[0]) && (rhs.v[1] == v[1]) && (rhs.v[2] == v[2]) )
-            return true;
-
-        return false;
-	}
-
-	std::string str() const
-	{
-		std::stringstream ss;
-
-		ss<<"[";
-		for(int i=0; i<N; i++)
-		{
-			ss<<v[i];
-			if(i != N-1)
-				ss<<", ";
-		}
-		ss<<"]";
-
-		return ss.str();
-	}
-
-	//bracket operator - allows assignment to the vector
-	CUDA_CALLABLE T& operator[](const unsigned int i)
-	{
-        return v[i];
-    }
-
-};
-
-
-}	//end namespace rts
-
-template <typename T, int N>
-std::ostream& operator<<(std::ostream& os, stim::vec<T, N> v)
-{
-    os<<v.str();
-    return os;
-}
-
-
-
-template <typename T, int N>
-CUDA_CALLABLE stim::vec<T, N> operator*(T lhs, stim::vec<T, N> rhs)
-{
-	stim::vec<T, N> r;
-
-    return rhs * lhs;
-}
-
-//#if __GNUC__ > 3 && __GNUC_MINOR__ > 7
-//template<class T, int N> using rtsVector = rts::vector<T, N>;
-//#endif
-
-#endif
-#include "../math/vector.h"
+#include "../math/mathvec.h"
 #include "../math/quaternion.h"
 #include "../math/matrix.h"
@@ -11,32 +11,32 @@ namespace stim{
 class camera
 {
-	vec<float, 3> d;	//direction that the camera is pointing
-	vec<float, 3> p;	//position of the camera
-	vec<float, 3> up;	//"up" direction
+	vec<float> d;	//direction that the camera is pointing
+	vec<float> p;	//position of the camera
+	vec<float> up;	//"up" direction
 	float focus;		//focal length of the camera
 	float fov;
 	//private function makes sure that the up vector is orthogonal to the direction vector and both are normalized
 	void stabalize()
 	{
-		vec<float, 3> side = up.cross(d);
+		vec<float> side = up.cross(d);
 		up = d.cross(side);
 		up = up.norm();
 		d = d.norm();
 	}
 public:
-	void setPosition(vec<float, 3> pos)
+	void setPosition(vec<float> pos)
 	{
 		p = pos;
 	}
-	void setPosition(float x, float y, float z){setPosition(vec<float, 3>(x, y, z));}
+	void setPosition(float x, float y, float z){setPosition(vec<float>(x, y, z));}
 	void setFocalDistance(float distance){focus = distance;}
 	void setFOV(float field_of_view){fov = field_of_view;}
-	void LookAt(vec<float, 3> pos)
+	void LookAt(vec<float> pos)
 	{
 		//find the new direction
 		d = pos - p;
@@ -47,22 +47,22 @@ public:
 		//stabalize the camera
 		stabalize();
 	}
-	void LookAt(float px, float py, float pz){LookAt(vec<float, 3>(px, py, pz));}
-	void LookAt(vec<float, 3> pos, vec<float, 3> new_up){up = new_up; LookAt(pos);}
-	void LookAt(float px, float py, float pz, float ux, float uy, float uz){LookAt(vec<float, 3>(px, py, pz), vec<float, 3>(ux, uy, uz));}
+	void LookAt(float px, float py, float pz){LookAt(vec<float>(px, py, pz));}
+	void LookAt(vec<float> pos, vec<float> new_up){up = new_up; LookAt(pos);}
+	void LookAt(float px, float py, float pz, float ux, float uy, float uz){LookAt(vec<float>(px, py, pz), vec<float>(ux, uy, uz));}
 	void LookAtDolly(float lx, float ly, float lz)
 	{
 		//find the current focus point
-		vec<float, 3> f = p + focus*d;
-		vec<float, 3> T = vec<float, 3>(lx, ly, lz) - f;
+		vec<float> f = p + focus*d;
+		vec<float> T = vec<float>(lx, ly, lz) - f;
 		p = p + T;
 	}
-	void Dolly(vec<float, 3> direction)
+	void Dolly(vec<float> direction)
 	{
 		p = p+direction;
 	}
-	void Dolly(float x, float y, float z){Dolly(vec<float, 3>(x, y, z));}
+	void Dolly(float x, float y, float z){Dolly(vec<float>(x, y, z));}
 	void Push(float delta)
 	{
 		if(delta > focus)
@@ -80,7 +80,7 @@ public:
 		qx.CreateRotation(theta_x, up[0], up[1], up[2]);
 		//y rotation is around the side axis
-		vec<float, 3> side = up.cross(d);
+		vec<float> side = up.cross(d);
 		quaternion<float> qy;
 		qy.CreateRotation(theta_y, side[0], side[1], side[2]);
@@ -118,28 +118,28 @@ public:
 	void OrbitFocus(float theta_x, float theta_y)
 	{
 		//find the focal point
-		vec<float, 3> focal_point = p + focus*d;
+		vec<float> focal_point = p + focus*d;
 		//center the coordinate system on the focal point
-		vec<float, 3> centered = p - (focal_point - vec<float, 3>(0, 0, 0));
+		vec<float> centered = p - (focal_point - vec<float>(0, 0, 0));
 		//create the x rotation (around the up vector)
 		quaternion<float> qx;
 		qx.CreateRotation(theta_x, up[0], up[1], up[2]);
-		centered = vec<float, 3>(0, 0, 0) + qx.toMatrix3()*(centered - vec<float, 3>(0, 0, 0));
+		centered = vec<float>(0, 0, 0) + qx.toMatrix3()*(centered - vec<float>(0, 0, 0));
 		//get a side vector for theta_y rotation
-		vec<float, 3> side = up.cross((vec<float, 3>(0, 0, 0) - centered).norm());
+		vec<float> side = up.cross((vec<float>(0, 0, 0) - centered).norm());
 		quaternion<float> qy;
 		qy.CreateRotation(theta_y, side[0], side[1], side[2]);
-		centered = vec<float, 3>(0, 0, 0) + qy.toMatrix3()*(centered - vec<float, 3>(0, 0, 0));
+		centered = vec<float>(0, 0, 0) + qy.toMatrix3()*(centered - vec<float>(0, 0, 0));
 		//perform the rotation on the centered camera position
 		//centered = final.toMatrix()*centered;
 		//re-position the camera
-		p = centered + (focal_point - vec<float, 3>(0, 0, 0));
+		p = centered + (focal_point - vec<float>(0, 0, 0));
 		//make sure we are looking at the focal point
 		LookAt(focal_point);
@@ -151,21 +151,21 @@ public:
 	void Slide(float u, float v)
 	{
-		vec<float, 3> V = up.norm();
-		vec<float, 3> U = up.cross(d).norm();
+		vec<float> V = up.norm();
+		vec<float> U = up.cross(d).norm();
 		p = p + (V * v) + (U * u);
 	}
 	//accessor methods
-	vec<float, 3> getPosition(){return p;}
-	vec<float, 3> getUp(){return up;}
-	vec<float, 3> getDirection(){return d;}
-	vec<float, 3> getLookAt(){return p + focus*d;}
+	vec<float> getPosition(){return p;}
+	vec<float> getUp(){return up;}
+	vec<float> getDirection(){return d;}
+	vec<float> getLookAt(){return p + focus*d;}
 	float getFOV(){return fov;}
 	//output the camera settings
-	const void print(std::ostream& output)
+	void print(std::ostream& output)
 	{
 		output<<"Position: "<<p.str()<<std::endl;
@@ -182,9 +182,9 @@ public:
 	//constructor
 	camera()
 	{
-		p = vec<float, 3>(0, 0, 0);
-		d = vec<float, 3>(0, 0, 1);
-		up = vec<float, 3>(0, 1, 0);
+		p = vec<float>(0, 0, 0);
+		d = vec<float>(0, 0, 1);
+		up = vec<float>(0, 1, 0);
 		focus = 1;
 	}
-#ifndef STIM_SPH_HARMONICS
-#define STIM_SPH_HARMONICS
+#ifndef STIM_GL_SPHARMONICS_H
+#define STIM_GL_SPHARMONICS_H
 #include <GL/gl.h>
 #include <stim/gl/error.h>
 #include <stim/visualization/colormap.h>
-#include <vector>
+#include <stim/math/spharmonics.h>
-#define PI 3.14159
-#define WIRE_SCALE 1.001
 namespace stim{
-	class sph_harmonics{
+	
+template <typename T>
+class gl_spharmonics : public spharmonics<T>{
-	private:
+protected:
+	
+	T* func;		//stores the raw function data (samples at each point)
-		double* func;		//stores the raw function data (samples at each point)
+	GLuint color_tex;	//texture map that acts as a colormap for the spherical function
+		
+	unsigned int N;		//resolution of the spherical grid
-		GLuint color_tex;	//texture map that acts as a colormap for the spherical function
-		
-		unsigned int N;		//resolution of the spherical grid
+	void gen_function(){
-		std::vector<double> C;	//list of SH coefficients
+		//initialize the function to zero
+		memset(func, 0, sizeof(double) * N * N);
+		double theta, phi;
+		double result;
+		int l, m;
-		//evaluates an associated Legendre polynomial (-l <= m <= l)
-		double P(int l,int m,double x)
-		{
-			// evaluate an Associated Legendre Polynomial P(l,m,x) at x
-			double pmm = 1.0;
-			if(m>0) {
-				double somx2 = sqrt((1.0-x)*(1.0+x));
-				double fact = 1.0;
-				for(int i=1; i<=m; i++) {
-					pmm *= (-fact) * somx2;
-					fact += 2.0;
-				}
-			}
-			if(l==m) return pmm;
-			double pmmp1 = x * (2.0*m+1.0) * pmm;
-			if(l==m+1) return pmmp1;
-			double pll = 0.0;
-			for(int ll=m+2; ll<=l; ++ll) {
-				pll = ( (2.0*ll-1.0)*x*pmmp1-(ll+m-1.0)*pmm ) / (ll-m);
-				pmm = pmmp1;
-				pmmp1 = pll;
-			}
-			return pll;
-		}
-		//recursively calculate a factorial given a positive integer n
-		unsigned int factorial(unsigned int n) {
-			if (n == 0)
-			   return 1;
-			return n * factorial(n - 1);
-		}
-
-		//calculate the SH scaling constant
-		double K(int l, int m){
+		l = m = 0;
+		//for each coefficient
+		for(unsigned int c = 0; c < C.size(); c++){				
-			// renormalisation constant for SH function
-			double temp = ((2.0*l+1.0)*factorial(l-m)) / (4.0*PI*factorial(l+m));
-			return sqrt(temp);
-		}
+			//iterate through the entire 2D grid representing the function
+			for(unsigned int xi = 0; xi < N; xi++){
+				for(unsigned int yi = 0; yi < N; yi++){
-		//calculate the value of the SH basis function (l, m) at (theta, phi)
-			//here, theta = [0, PI], phi = [0, 2*PI]
-		double SH(int l, int m, double theta, double phi){
-			// return a point sample of a Spherical Harmonic basis function
-			// l is the band, range [0..N]
-			// m in the range [-l..l]
-			// theta in the range [0..Pi]
-			// phi in the range [0..2*Pi]
-			const double sqrt2 = sqrt(2.0);
-			if(m==0) return K(l,0)*P(l,m,cos(theta));
-			else if(m>0) return sqrt2*K(l,m)*cos(m*phi)*P(l,m,cos(theta));
-			else return sqrt2*K(l,-m)*sin(-m*phi)*P(l,-m,cos(theta));
-		}
+					//get the spherical coordinates for each grid point
+					theta = (2 * PI) * ((double)xi / (N-1));
+					phi = PI * ((double)yi / (N-1));
-		void gen_function(){
+					//sum the contribution of the current spherical harmonic based on the coefficient
+					result = C[c] * SH(l, m, theta, phi);
-			//initialize the function to zero
-			memset(func, 0, sizeof(double) * N * N);
-
-			double theta, phi;
-			double result;
-			int l, m;
-
-			l = m = 0;
-			for(unsigned int c = 0; c < C.size(); c++){
-				
-
-				for(unsigned int xi = 0; xi < N; xi++)
-					for(unsigned int yi = 0; yi < N; yi++){
-
-						theta = (2 * PI) * ((double)xi / (N-1));
-						phi = PI * ((double)yi / (N-1));
-						result = C[c] * SH(l, m, phi, theta);		//phi and theta are reversed here (damn physicists)
-						func[yi * N + xi] += result;
-					}
-
-				m++;			//increment m
-
-				//if we're in a new tier, increment l and set m = -l
-				if(m > l){		
-					l++;
-					m = -l;
+					//store the result in a 2D array (which will later be used as a texture map)
+					func[yi * N + xi] += result;
 				}
 			}
-		}
-
-		void gl_prep_draw(){
-			//enable depth testing
-				//this has to be used instead of culling because the sphere can have negative values
-			glEnable(GL_DEPTH_TEST);
-			glDepthMask(GL_TRUE);
-			glEnable(GL_TEXTURE_2D);	//enable 2D texture mapping
+			//keep track of m and l here
+			m++;			//increment m
+			//if we're in a new tier, increment l and set m = -l
+			if(m > l){		
+				l++;
+				m = -l;
+			}
 		}
+	}
-		//draw a texture mapped sphere representing the function surface
-		void gl_draw_sphere() {
-
-			//PI is used to convert from spherical to cartesian coordinates
-			//const double PI = 3.14159;
-
-			//bind the 2D texture representing the color map
-			glBindTexture(GL_TEXTURE_2D, color_tex);
-
-			//Draw the Sphere
-			int i, j;
-
-			for(i = 1; i <= N-1; i++) {
-				double phi0 = PI * ((double) (i - 1) / (N-1));
-				double phi1 = PI * ((double) i / (N-1));
-
-				glBegin(GL_QUAD_STRIP);
-				for(j = 0; j <= N; j++) {
-
-					//calculate the indices into the function array
-					int phi0_i = i-1;
-					int phi1_i = i;
-					int theta_i = j;
-					if(theta_i == N)
-						theta_i = 0;
-
-					double v0 = func[phi0_i * N + theta_i];
-					double v1 = func[phi1_i * N + theta_i];
-
-					v0 = fabs(v0);
-					v1 = fabs(v1);
-
-
-					double theta = 2 * PI * (double) (j - 1) / N;
-					double x0 = v0 * cos(theta) * sin(phi0);
-					double y0 = v0 * sin(theta) * sin(phi0);
-					double z0 = v0 * cos(phi0);
-
-					double x1 = v1 * cos(theta) * sin(phi1);
-					double y1 = v1 * sin(theta) * sin(phi1);
-					double z1 = v1 * cos(phi1);
-
-					glTexCoord2f(theta / (2 * PI), phi0 / PI);
-					glVertex3f(x0, y0, z0);
-
-					glTexCoord2f(theta / (2 * PI), phi1 / PI);
-					glVertex3f(x1, y1, z1);
-				}
-				glEnd();
-			}
-		}
+	void gl_prep_draw(){
-		//draw a wire frame sphere representing the function surface
-		void gl_draw_wireframe() {
+		//enable depth testing
+			//this has to be used instead of culling because the sphere can have negative values
+		glEnable(GL_DEPTH_TEST);
+		glDepthMask(GL_TRUE);
+		glEnable(GL_TEXTURE_2D);	//enable 2D texture mapping
+	}
-					//PI is used to convert from spherical to cartesian coordinates
-					//const double PI = 3.14159;
+	//draw a texture mapped sphere representing the function surface
+	void gl_draw_sphere() {
-					//bind the 2D texture representing the color map
-					glDisable(GL_TEXTURE_2D);
-					glColor3f(0.0f, 0.0f, 0.0f);
+		//bind the 2D texture representing the color map
+		glBindTexture(GL_TEXTURE_2D, color_tex);
-					//Draw the Sphere
-					int i, j;
+		//Draw the Sphere
+		int i, j;
-					for(i = 1; i <= N-1; i++) {
-						double phi0 = PI * ((double) (i - 1) / (N-1));
-						double phi1 = PI * ((double) i / (N-1));
+		for(i = 1; i <= N-1; i++) {
+			double phi0 = PI * ((double) (i - 1) / (N-1));
+			double phi1 = PI * ((double) i / (N-1));
-						glBegin(GL_LINE_STRIP);
-						for(j = 0; j <= N; j++) {
+			glBegin(GL_QUAD_STRIP);
+			for(j = 0; j <= N; j++) {
-							//calculate the indices into the function array
-							int phi0_i = i-1;
-							int phi1_i = i;
-							int theta_i = j;
-							if(theta_i == N)
-								theta_i = 0;
+				//calculate the indices into the function array
+				int phi0_i = i-1;
+				int phi1_i = i;
+				int theta_i = j;
+				if(theta_i == N)
+					theta_i = 0;
-							double v0 = func[phi0_i * N + theta_i];
-							double v1 = func[phi1_i * N + theta_i];
+				double v0 = func[phi0_i * N + theta_i];
+				double v1 = func[phi1_i * N + theta_i];
-							v0 = fabs(v0);
-							v1 = fabs(v1);
+				v0 = fabs(v0);
+				v1 = fabs(v1);
-							double theta = 2 * PI * (double) (j - 1) / N;
-							double x0 = WIRE_SCALE * v0 * cos(theta) * sin(phi0);
-							double y0 = WIRE_SCALE * v0 * sin(theta) * sin(phi0);
-							double z0 = WIRE_SCALE * v0 * cos(phi0);
+				double theta = 2 * PI * (double) (j - 1) / N;
+				double x0 = v0 * cos(theta) * sin(phi0);
+				double y0 = v0 * sin(theta) * sin(phi0);
+				double z0 = v0 * cos(phi0);
-							double x1 = WIRE_SCALE * v1 * cos(theta) * sin(phi1);
-							double y1 = WIRE_SCALE * v1 * sin(theta) * sin(phi1);
-							double z1 = WIRE_SCALE * v1 * cos(phi1);
+				double x1 = v1 * cos(theta) * sin(phi1);
+				double y1 = v1 * sin(theta) * sin(phi1);
+				double z1 = v1 * cos(phi1);
-							glTexCoord2f(theta / (2 * PI), phi0 / PI);
-							glVertex3f(x0, y0, z0);
+				glTexCoord2f(theta / (2 * PI), phi0 / PI);
+				glVertex3f(x0, y0, z0);
-							glTexCoord2f(theta / (2 * PI), phi1 / PI);
-							glVertex3f(x1, y1, z1);
-						}
-						glEnd();
-					}
-				}
+				glTexCoord2f(theta / (2 * PI), phi1 / PI);
+				glVertex3f(x1, y1, z1);
+			}
+			glEnd();
+		}
+	}
-		void init(unsigned int n){
+	void gl_init(unsigned int n){
-			//set the sphere resolution
-			N = n;
+		//set the sphere resolution
+		N = n;
-			//allocate space for the color map
-			unsigned int bytes = N * N * sizeof(unsigned char) * 3;
-			unsigned char* color_image;
-			color_image = (unsigned char*) malloc(bytes);
+		//allocate space for the color map
+		unsigned int bytes = N * N * sizeof(unsigned char) * 3;
+		unsigned char* color_image;
+		color_image = (unsigned char*) malloc(bytes);
-			//allocate space for the function
-			func = (double*) malloc(N * N * sizeof(double));
+		//allocate space for the function
+		func = (double*) malloc(N * N * sizeof(double));
-			//generate a function (temporary)
-			gen_function();
+		//generate a functional representation that will be used for the texture map and vertices
+		gen_function();
-			//generate a colormap from the function
-			stim::cpu2cpu<double>(func, color_image, N*N, stim::cmBrewer);
+		//generate a colormap from the function
+		stim::cpu2cpu<double>(func, color_image, N*N, stim::cmBrewer);
-			//prep everything for drawing
-			gl_prep_draw();			
+		//prep everything for drawing
+		gl_prep_draw();			
-			//generate an OpenGL texture map in the current context
-			glGenTextures(1, &color_tex);
-			//bind the texture
-			glBindTexture(GL_TEXTURE_2D, color_tex);
+		//generate an OpenGL texture map in the current context
+		glGenTextures(1, &color_tex);
+		//bind the texture
+		glBindTexture(GL_TEXTURE_2D, color_tex);
-			//copy the color data from the buffer to the GPU
-			glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, N, N, 0, GL_RGB, GL_UNSIGNED_BYTE, color_image);
+		//copy the color data from the buffer to the GPU
+		glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, N, N, 0, GL_RGB, GL_UNSIGNED_BYTE, color_image);
-			//initialize all of the texture parameters
-			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
-			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);
-			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
-			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
-			glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);
+		//initialize all of the texture parameters
+		glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
+		glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);
+		glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
+		glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
+		glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);
-			//free the buffer
-			free(color_image);
-		}
-
-
-	public:
+		//free the buffer
+		free(color_image);
+	}
-		void glRender(){
-			//set all OpenGL parameters required for drawing
-			gl_prep_draw();
-
-			//draw the sphere
-			gl_draw_sphere();
-			//gl_draw_wireframe();
-
-		}
-
-		void glInit(unsigned int n){
-			init(n);
-		}
-
-		void push(double c){
-			C.push_back(c);
-		}
+public:
+	void glRender(){
+		//set all OpenGL parameters required for drawing
+		gl_prep_draw();
+		//draw the sphere
+		gl_draw_sphere();
+	}
+	void glInit(unsigned int n = 256){
+		gl_init(n);
+		gen_function();
+	}
-	};		//end class sph_harmonics
+};		//end gl_spharmonics
+};		//end namespace stim
-}
-#endif
+#endif
 \ No newline at end of file
@@ -6,6 +6,7 @@
 #include <fstream>
 #include <stdlib.h>
 #include <stim/parser/parser.h>
+#include <stim/math/mathvec.h>
 namespace stim{
@@ -34,37 +35,14 @@ protected:
 	enum token_type { OBJ_INVALID, OBJ_V, OBJ_VT, OBJ_VN, OBJ_P, OBJ_L, OBJ_F };
-	struct vertex : public std::vector<T>{
+	struct vertex : public stim::vec<T>{
-		using std::vector<T>::push_back;
-		using std::vector<T>::size;
-		using std::vector<T>::at;
+		using vec<T>::push_back;
+		using vec<T>::size;
+		using vec<T>::at;
 		vertex(){}
-		//constructors for quickly building vertices with arbitrary numbers of components
-		vertex(T x){
-			push_back(x);
-		}
-
-		vertex(T x, T y){
-			push_back(x);
-			push_back(y);
-		}
-
-		vertex(T x, T y, T z){
-			push_back(x);
-			push_back(y);
-			push_back(z);
-		}
-
-		vertex(T x, T y, T z, T w){
-			push_back(x);
-			push_back(y);
-			push_back(z);
-			push_back(w);
-		}
-
 		//constructor creates a vertex from a line string
 		vertex(std::string line){
@@ -305,6 +283,11 @@ protected:
 	}
 public:
+	/// Constructor loads a Wavefront OBJ file
+	obj(std::string filename){
+		load(filename);
+	}
+
 	//functions for setting the texture coordinate for the next vertex
 	void TexCoord(T x){ update_vt(vertex(x));}
 	void TexCoord(T x, T y){ update_vt(vertex(x, y));}
@@ -504,6 +487,47 @@ public:
 	}
+	/// Return the number of position vertices in the OBJ model
+	unsigned int numV(){
+		return V.size();
+	}
+
+	/// Retrieve the vertex stored in index i
+
+	/// @param i is the desired vertex index
+	stim::vec<T> getV(unsigned int i){
+
+		stim::vec<T> v = V[i];
+		return v;
+	}
+
+	stim::vec<T> centroid(){
+
+		//get the number of coordinates
+		unsigned int N = V[0].size();
+
+		//allocate space for the minimum and maximum coordinate points (bounding box corners)
+		stim::vec<float> vmin, vmax;
+		vmin.resize(N);
+		vmax.resize(N);
+
+		//find the minimum and maximum value for each coordinate
+		unsigned int NV = V.size();
+		for(unsigned int v = 0; v < NV; v++)
+			for(unsigned int i = 0; i < N; i++){
+
+				if(V[v][i] < vmin[i])
+					vmin[i] = V[v][i];
+				if(V[v][i] > vmax[i])
+					vmax[i] = V[v][i];
+			}
+
+		//find the centroid using the min and max points
+		stim::vec<T> c = vmin * 0.5 + vmax * 0.5;
+
+		return c;
+	}
+
 };