fixed class compatibility with stim::vec3

David Mayerich
1 parent 814eb271
Showing 9 changed files with 142 additions and 127 deletions Show diff stats
stim/biomodels/network.h
stim/math/circle.h
stim/math/matrix.h
stim/math/plane.h
stim/math/vec3.h
stim/math/vector.h
stim/visualization/aaboundingbox.h
stim/visualization/camera.h
stim/visualization/cylinder.h
@@ -8,7 +8,7 @@
 #include <algorithm>
 #include <string.h>
 #include <math.h>
-#include <stim/math/vector.h>
+#include <stim/math/vec3.h>
 #include <stim/visualization/obj.h>
 #include <stim/visualization/cylinder.h>
 #include <ANN/ANN.h>
@@ -37,7 +37,7 @@ class network{
 		/// Constructor - creates an edge from a list of points by calling the stim::fiber constructor
 		///@param p is an array of positions in space
-		edge(std::vector< stim::vec<T> > p) : cylinder<T>(p){}
+		edge(std::vector< stim::vec3<T> > p) : cylinder<T>(p){}
 		/// Copy constructor creates an edge from a fiber
 		edge(stim::cylinder<T> f) : cylinder<T>(f) {}
@@ -61,20 +61,20 @@ class network{
 	};
 	///Node class that stores the physical position of the node as well as the edges it is connected to (edges that connect to it), As well as any additional data necessary.
-	class vertex : public stim::vec<T>
+	class vertex : public stim::vec3<T>
 	{
 		public:
 			//std::vector<unsigned int> edges;		//indices of edges connected to this node.
 			std::vector<unsigned int> e[2];			//indices of edges going out (e[0]) and coming in (e[1])
-			//stim::vec<T> p;						//position of this node in physical space.
+			//stim::vec3<T> p;						//position of this node in physical space.
 			//constructor takes a stim::vec
-			vertex(stim::vec<T> p) : stim::vec<T>(p){}
+			vertex(stim::vec3<T> p) : stim::vec3<T>(p){}
 			/// Output the vertex information as a string
 			std::string	str(){
 				std::stringstream ss;
-				ss<<"\t(x, y, z) = "<<stim::vec<T>::str();
+				ss<<"\t(x, y, z) = "<<stim::vec3<T>::str();
 				if(e[0].size() > 0){
 					ss<<"\t> ";
@@ -129,7 +129,11 @@ public:
 			std::vector< stim::vec<T> > c;				//allocate an array of points for the vessel centerline
 			O.getLine(l, c);							//get the fiber centerline
-			edge new_edge = c;							//create an edge from the given centerline
+			std::vector< stim::vec3<T> > c3(c.size());
+			for(size_t j = 0; j < c.size(); j++)
+				c3[j] = c[j];
+
+			edge new_edge = c3;							//create an edge from the given centerline
 			unsigned int I = new_edge.size();			//calculate the number of points on the centerline
 			//get the first and last vertex IDs for the line
@@ -222,7 +226,7 @@ public:
 	float gaussianFunction(float x, float std=25){ return exp(-x/(2*std*std));} // by default std = 25
     // stim 3d vector to annpoint of 3 dimensions
-	void stim2ann(ANNpoint &a, stim::vec<T> b){
+	void stim2ann(ANNpoint &a, stim::vec3<T> b){
 		a[0] = b[0];
 		a[1] = b[1];
 		a[2] = b[2];
@@ -278,10 +282,9 @@ public:
 		ANNdistArray dists = new ANNdist[1];     // near neighbor distances
 		ANNidxArray nnIdx = new ANNidx[1];				// near neighbor indices // allocate near neigh indices
-		stim::vec<T> p0, p1;
-		float m0, m1;
+		stim::vec3<T> p0, p1;
+		float m1;
 		float M = 0;											//stores the total metric value
-		float l;												//stores the segment length
 		float L = 0;											//stores the total network length
 		ANNpoint queryPt = annAllocPt(3);
 		for(unsigned e = 0; e < R.E.size(); e++){					//for each edge in A
@@ -292,7 +295,7 @@ public:
 				p1 = R.E[e][p];									//get the next point in the edge
 				stim2ann(queryPt, p1);
 				kdt->annkSearch( queryPt, 1, nnIdx, dists, eps);	//find the distance between A and the current network
-				m1 = 1.0f - gaussianFunction(dists[0], sigma);		//calculate the metric value based on the distance
+				m1 = 1.0f - gaussianFunction((float)dists[0], sigma);		//calculate the metric value based on the distance
 				R.E[e].set_mag(m1, p, 1);						//set the error for the second point in the segment
 			}
@@ -17,7 +17,7 @@ class circle : plane&lt;T&gt;
 private:
-	stim::vec<T> Y;
+	stim::vec3<T> Y;
 	CUDA_CALLABLE void
 	init()
@@ -48,7 +48,7 @@ public:
 	circle(T size, T z_pos = (T)0) : plane<T>()
 	{
 		init();
-		center(stim::vec<T>(0,0,z_pos));
+		center(stim::vec3<T>(0,0,z_pos));
 		scale(size);
 	}
@@ -56,7 +56,7 @@ public:
 	///@param c: x,y,z location of the center.
 	///@param n: x,y,z direction of the normal.	
 	CUDA_CALLABLE
-	circle(vec<T> c, vec<T> n = vec<T>(0,0,1)) : plane<T>()
+	circle(vec3<T> c, vec3<T> n = vec3<T>(0,0,1)) : plane<T>()
 	{
 		center(c);
 		normal(n);
@@ -68,7 +68,7 @@ public:
 	///@param s: size of the rectangle.
 	///@param n: x,y,z direction of the normal.
 	CUDA_CALLABLE 
-	circle(vec<T> c, T s, vec<T> n = vec<T>(0,0,1)) : plane<T>()
+	circle(vec3<T> c, T s, vec3<T> n = vec3<T>(0,0,1)) : plane<T>()
 	{
 		init();
 		center(c);
@@ -82,7 +82,7 @@ public:
 	///@param n: x,y,z direction of the normal.
 	///@param u: x,y,z direction for the zero vector (from where the rotation starts)
 	CUDA_CALLABLE
-	circle(vec<T> c, T s, vec<T> n = vec<T>(0,0,1), vec<T> u = vec<T>(1, 0, 0)) : plane<T>()
+	circle(vec3<T> c, T s, vec3<T> n = vec3<T>(0,0,1), vec3<T> u = vec3<T>(1, 0, 0)) : plane<T>()
 	{
 		init();
 		setU(u);
@@ -103,16 +103,15 @@ public:
 	///sets the normal for the cirlce
 	///@param n: x,y,z direction of the normal.
 	CUDA_CALLABLE void
-	normal(vec<T> n)
+	normal(vec3<T> n)
 	{
 		rotate(n, Y);
 	}
 	///sets the center of the circle.
 	///@param n: x,y,z location of the center.
-	CUDA_CALLABLE T
-	center(vec<T> p)
-	{
+	CUDA_CALLABLE void
+	center(vec3<T> p){
 		this->P = p;
 	}
@@ -127,17 +126,17 @@ public:
 	}
 	///get the world space value given the planar coordinates a, b in [0, 1]
-	CUDA_CALLABLE stim::vec<T> p(T a, T b)
+	CUDA_CALLABLE stim::vec3<T> p(T a, T b)
 	{
-		stim::vec<T> result;
+		stim::vec3<T> result;
-		vec<T> A = this->P - this->U * (T)0.5 - Y * (T)0.5;
+		vec3<T> A = this->P - this->U * (T)0.5 - Y * (T)0.5;
 		result = A + this->U * 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> operator()(T a, T b)
+	CUDA_CALLABLE stim::vec3<T> operator()(T a, T b)
 	{
 		return p(a,b);
 	}
@@ -145,11 +144,11 @@ public:
 	///returns a vector with the points on the initialized circle.
 	///connecting the points results in a circle.
 	///@param n: integer for the number of points representing the circle.
-	std::vector<stim::vec<T> >
+	std::vector<stim::vec3<T> >
 	getPoints(int n)
 	{
-		std::vector<stim::vec<T> > result;
-		stim::vec<T> point;
+		std::vector<stim::vec3<T> > result;
+		stim::vec3<T> point;
 		T x,y;
 		float step = 360.0/(float) n;
 		for(float j = 0; j <= 360.0; j += step)
@@ -164,7 +163,7 @@ public:
 	///returns a vector with the points on the initialized circle.
 	///connecting the points results in a circle.
 	///@param n: integer for the number of points representing the circle.
-	stim::vec<T>
+	stim::vec3<T>
 	p(T theta)
 	{
 		T x,y;
@@ -5,6 +5,7 @@
 #include <string.h>
 #include <iostream>
 #include <stim/math/vector.h>
+#include <stim/math/vec3.h>
 #include <stim/cuda/cudatools/callable.h>
 namespace stim{
@@ -54,7 +55,7 @@ struct matrix
 	vec<Y> operator*(vec<Y> rhs){
 		unsigned int N = rhs.size();
-		vec<Y> result;
+		vec3Y> result;
 		result.resize(N);
 		for(int r=0; r<N; r++)
@@ -2,7 +2,7 @@
 #define STIM_PLANE_H
 #include <iostream>
-#include <stim/math/vector.h>
+#include <stim/math/vec3.h>
 #include <stim/cuda/cudatools/callable.h>
 #include <stim/math/quaternion.h>
@@ -22,17 +22,17 @@ template &lt;typename T&gt;
 class plane
 {
 	protected:
-		stim::vec<T> P;
-		stim::vec<T> N;
-		stim::vec<T> U;
+		stim::vec3<T> P;
+		stim::vec3<T> N;
+		stim::vec3<T> U;
 		///Initializes the plane with standard coordinates.
 		///
 		CUDA_CALLABLE void init()
 		{
-			P = stim::vec<T>(0, 0, 0);
-			N = stim::vec<T>(0, 0, 1);
-			U = stim::vec<T>(1, 0, 0);
+			P = stim::vec3<T>(0, 0, 0);
+			N = stim::vec3<T>(0, 0, 1);
+			U = stim::vec3<T>(1, 0, 0);
 		}
 	public:
@@ -42,7 +42,7 @@ class plane
 			init();
 		}
-		CUDA_CALLABLE plane(vec<T> n, vec<T> p = vec<T>(0, 0, 0))
+		CUDA_CALLABLE plane(vec3<T> n, vec3<T> p = vec3<T>(0, 0, 0))
 		{
 			init();
 			P = p;
@@ -56,11 +56,11 @@ class plane
 		}
 		//create a plane from three points (a triangle)
-		CUDA_CALLABLE plane(vec<T> a, vec<T> b, vec<T> c)
+		CUDA_CALLABLE plane(vec3<T> a, vec3<T> b, vec3<T> c)
 		{
 			init();
 			P = c;
-			stim::vec<T> n = (c - a).cross(b - a);
+			stim::vec3<T> n = (c - a).cross(b - a);
 			try
 			{
 				if(n.len() != 0)
@@ -84,17 +84,17 @@ class plane
 		}
-		CUDA_CALLABLE vec<T> n()
+		CUDA_CALLABLE vec3<T> n()
 		{
 			return N;
 		}
-		CUDA_CALLABLE vec<T> p()
+		CUDA_CALLABLE vec3<T> p()
 		{
 			return P;
 		}
-		CUDA_CALLABLE vec<T> u()
+		CUDA_CALLABLE vec3<T> u()
 		{
 			return U;
 		}
@@ -107,7 +107,7 @@ class plane
 		}
 		//determines how a vector v intersects the plane (1 = intersects front, 0 = within plane,     -1 = intersects back)
-		CUDA_CALLABLE int face(vec<T> v){
+		CUDA_CALLABLE int face(vec3<T> v){
 			T dprod = v.dot(N);             //get the dot product between v and N
@@ -121,46 +121,46 @@ class plane
 		}
 		//determine on which side of the plane a point lies (1 = front, 0 = on the plane, -1 = bac    k)
-		CUDA_CALLABLE int side(vec<T> p){
+		CUDA_CALLABLE int side(vec3<T> p){
-			vec<T> v = p - P;    //get the vector from P to the query point p
+			vec3<T> v = p - P;    //get the vector from P to the query point p
 			return face(v);
 		}
 		//compute the component of v that is perpendicular to the plane
-		CUDA_CALLABLE vec<T> perpendicular(vec<T> v){
+		CUDA_CALLABLE vec3<T> perpendicular(vec3<T> v){
 			return N * v.dot(N);
 		}
 		//compute the projection of v in the plane
-		CUDA_CALLABLE vec<T> parallel(vec<T> v){
+		CUDA_CALLABLE vec3<T> parallel(vec3<T> v){
 			return v - perpendicular(v);
 		}
-		CUDA_CALLABLE void setU(vec<T> v)
+		CUDA_CALLABLE void setU(vec3<T> v)
 		{
 			U = (parallel(v.norm())).norm();		
 		}
-		CUDA_CALLABLE void decompose(vec<T> v, vec<T>& para, vec<T>& perp){
+		CUDA_CALLABLE void decompose(vec3<T> v, vec3<T>& para, vec3<T>& perp){
 			perp = N * v.dot(N);
 			para = v - perp;
 		}
 		//get both the parallel and perpendicular components of a vector v w.r.t. the plane
-		CUDA_CALLABLE void project(vec<T> v, vec<T> &v_par, vec<T> &v_perp){
+		CUDA_CALLABLE void project(vec3<T> v, vec3<T> &v_par, vec3<T> &v_perp){
 			v_perp = v.dot(N);
 			v_par = v - v_perp;
 		}
 		//compute the reflection of v off of the plane
-		CUDA_CALLABLE vec<T> reflect(vec<T> v){
+		CUDA_CALLABLE vec3<T> reflect(vec3<T> v){
 			//compute the reflection using N_prime as the plane normal
-			vec<T> par = parallel(v);
-			vec<T> r = (-v) + par * 2;
+			vec3<T> par = parallel(v);
+			vec3<T> r = (-v) + par * 2;
 			return r;
 		}
@@ -184,7 +184,7 @@ class plane
 		}
-		CUDA_CALLABLE void rotate(vec<T> n)
+		CUDA_CALLABLE void rotate(vec3<T> n)
 		{
 			quaternion<T> q;
 			q.CreateRotation(N, n);
@@ -194,7 +194,7 @@ class plane
 		}
-		CUDA_CALLABLE void rotate(vec<T> n, vec<T> &Y)
+		CUDA_CALLABLE void rotate(vec3<T> n, vec3<T> &Y)
 		{
 			quaternion<T> q;
 			q.CreateRotation(N, n);
@@ -205,7 +205,7 @@ class plane
 		}
-		CUDA_CALLABLE void rotate(vec<T> n, vec<T> &X, vec<T> &Y)
+		CUDA_CALLABLE void rotate(vec3<T> n, vec3<T> &X, vec3<T> &Y)
 		{
 			quaternion<T> q;
 			q.CreateRotation(N, n);
@@ -217,12 +217,12 @@ public:
 	std::string str() const{
 		std::stringstream ss;
-		size_t N = size();
+		const size_t N = 3;
 		ss<<"[";
 		for(size_t i=0; i<N; i++)
 		{
-			ss<<at(i);
+			ss<<ptr[i];
 			if(i != N-1)
 				ss<<", ";
 		}
@@ -230,7 +230,10 @@ public:
 		return ss.str();
 	}
-	};						//end class triple
+
+	size_t size(){ return 3; }
+
+	};						//end class vec3
 }							//end namespace stim
 /// Multiply a vector by a constant when the vector is on the right hand side
@@ -317,6 +317,15 @@ struct vec : public std::vector&lt;T&gt;
 		return *this;
 	}
+	/// Cast to a vec3
+	operator vec3<T>(){
+		vec3<T> r;
+		size_t N = std::min<size_t>(size(), 3);
+		for(size_t i = 0; i < N; i++)
+			r[i] = at(i);
+		return r;
+	}
+
 	/// Casting and assignment
 	template<typename Y>
 	vec<T> & operator=(vec<Y> rhs){
@@ -10,8 +10,8 @@ class aaboundingbox{
 public:
 	bool set;				//has the bounding box been set to include any points?
-	stim::vec<T> A;			//minimum point in the bounding box
-	stim::vec<T> B;			//maximum point in the bounding box
+	stim::vec3<T> A;			//minimum point in the bounding box
+	stim::vec3<T> B;			//maximum point in the bounding box
 	aaboundingbox(){					//constructor generates an empty bounding box
 		set = false;
@@ -21,7 +21,7 @@ public:
 	/// Test if a point is inside of the bounding box and returns true if it is.
 	/// @param p is the point to be tested
-	bool test(stim::vec<T> p){
+	bool test(stim::vec3<T> p){
 		for(unsigned d = 0; d < p.size(); p++){		//for each dimension
 			if(p[d] < A[d]) return false;			//if the point is less than the minimum bound, return false
@@ -33,7 +33,7 @@ public:
 	/// Expand the bounding box to include the specified point.
 	/// @param p is the point to be included
-	void expand(stim::vec<T> p){
+	void expand(stim::vec3<T> p){
 		if(!set){							//if the bounding box is empty, fill it with the current point
 			A = B = p;
@@ -47,12 +47,12 @@ public:
 	}
 	/// Return the center point of the bounding box as a stim::vec
-	stim::vec<T> center(){
+	stim::vec3<T> center(){
 		return (B + A) * 0.5;
 	}
 	/// Return the size of the bounding box as a stim::vec
-	stim::vec<T> size(){
+	stim::vec3<T> size(){
 		return (B - A);
 	}
@@ -11,32 +11,32 @@ namespace stim{
 class camera
 {
-	vec<float> d;	//direction that the camera is pointing
-	vec<float> p;	//position of the camera
-	vec<float> up;	//"up" direction
+	vec3<float> d;	//direction that the camera is pointing
+	vec3<float> p;	//position of the camera
+	vec3<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> side = up.cross(d);
+		vec3<float> side = up.cross(d);
 		up = d.cross(side);
 		up = up.norm();
 		d = d.norm();
 	}
 public:
-	void setPosition(vec<float> pos)
+	void setPosition(vec3<float> pos)
 	{
 		p = pos;
 	}
-	void setPosition(float x, float y, float z){setPosition(vec<float>(x, y, z));}
+	void setPosition(float x, float y, float z){setPosition(vec3<float>(x, y, z));}
 	void setFocalDistance(float distance){focus = distance;}
 	void setFOV(float field_of_view){fov = field_of_view;}
-	void LookAt(vec<float> pos)
+	void LookAt(vec3<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>(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 LookAt(float px, float py, float pz){LookAt(vec3<float>(px, py, pz));}
+	void LookAt(vec3<float> pos, vec3<float> new_up){up = new_up; LookAt(pos);}
+	void LookAt(float px, float py, float pz, float ux, float uy, float uz){LookAt(vec3<float>(px, py, pz), vec3<float>(ux, uy, uz));}
 	void LookAtDolly(float lx, float ly, float lz)
 	{
 		//find the current focus point
-		vec<float> f = p + focus*d;
-		vec<float> T = vec<float>(lx, ly, lz) - f;
+		vec3<float> f = p + focus*d;
+		vec3<float> T = vec3<float>(lx, ly, lz) - f;
 		p = p + T;
 	}
-	void Dolly(vec<float> direction)
+	void Dolly(vec3<float> direction)
 	{
 		p = p+direction;
 	}
-	void Dolly(float x, float y, float z){Dolly(vec<float>(x, y, z));}
+	void Dolly(float x, float y, float z){Dolly(vec3<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> side = up.cross(d);
+		vec3<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> focal_point = p + focus*d;
+		vec3<float> focal_point = p + focus*d;
 		//center the coordinate system on the focal point
-		vec<float> centered = p - (focal_point - vec<float>(0, 0, 0));
+		vec3<float> centered = p - (focal_point - vec3<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>(0, 0, 0) + qx.toMatrix3()*(centered - vec<float>(0, 0, 0));
+		centered = vec3<float>(0, 0, 0) + qx.toMatrix3()*(centered - vec3<float>(0, 0, 0));
 		//get a side vector for theta_y rotation
-		vec<float> side = up.cross((vec<float>(0, 0, 0) - centered).norm());
+		vec3<float> side = up.cross((vec3<float>(0, 0, 0) - centered).norm());
 		quaternion<float> qy;
 		qy.CreateRotation(theta_y, side[0], side[1], side[2]);
-		centered = vec<float>(0, 0, 0) + qy.toMatrix3()*(centered - vec<float>(0, 0, 0));
+		centered = vec3<float>(0, 0, 0) + qy.toMatrix3()*(centered - vec3<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>(0, 0, 0));
+		p = centered + (focal_point - vec3<float>(0, 0, 0));
 		//make sure we are looking at the focal point
 		LookAt(focal_point);
@@ -151,17 +151,17 @@ public:
 	void Slide(float u, float v)
 	{
-		vec<float> V = up.norm();
-		vec<float> U = up.cross(d).norm();
+		vec3<float> V = up.norm();
+		vec3<float> U = up.cross(d).norm();
 		p = p + (V * v) + (U * u);
 	}
 	//accessor methods
-	vec<float> getPosition(){return p;}
-	vec<float> getUp(){return up;}
-	vec<float> getDirection(){return d;}
-	vec<float> getLookAt(){return p + focus*d;}
+	vec3<float> getPosition(){return p;}
+	vec3<float> getUp(){return up;}
+	vec3<float> getDirection(){return d;}
+	vec3<float> getLookAt(){return p + focus*d;}
 	float getFOV(){return fov;}
 	//output the camera settings
@@ -182,9 +182,9 @@ public:
 	//constructor
 	camera()
 	{
-		p = vec<float>(0, 0, 0);
-		d = vec<float>(0, 0, 1);
-		up = vec<float>(0, 1, 0);
+		p = vec3<float>(0, 0, 0);
+		d = vec3<float>(0, 0, 1);
+		up = vec3<float>(0, 1, 0);
 		focus = 1;
 	}
@@ -2,7 +2,7 @@
 #define STIM_CYLINDER_H
 #include <iostream>
 #include <stim/math/circle.h>
-#include <stim/math/vector.h>
+#include <stim/math/vec3.h>
 namespace stim
@@ -25,11 +25,11 @@ class cylinder
 		///inits the cylinder from a list of points (inP) and radii (inM)
 		void
-		init(std::vector<stim::vec<T> > inP, std::vector<stim::vec<T> > inM)
+		init(std::vector<stim::vec3<T> > inP, std::vector<stim::vec<T> > inM)
 		{
 			mags = inM;
-			stim::vec<float> v1;
-			stim::vec<float> v2;
+			stim::vec3<float> v1;
+			stim::vec3<float> v2;
 			e.resize(inP.size());
 			if(inP.size() < 2)
 				return;
@@ -38,16 +38,16 @@ class cylinder
 			L.resize(inP.size());
 			T temp = (T)0;
 			L[0] = 0;
-			for(int i = 1; i < L.size(); i++)
+			for(size_t i = 1; i < L.size(); i++)
 			{
 				temp += (inP[i-1] - inP[i]).len();
 				L[i] = temp;
 			}
-			stim::vec<T> dr = (inP[1] - inP[0]).norm();
-			s = stim::circle<T>(inP[0], inM[0][0], dr, stim::vec<T>(1,0,0));
+			stim::vec3<T> dr = (inP[1] - inP[0]).norm();
+			s = stim::circle<T>(inP[0], inM[0][0], dr, stim::vec3<T>(1,0,0));
 			e[0] = s;
-			for(int i = 1; i < inP.size()-1; i++)
+			for(size_t i = 1; i < inP.size()-1; i++)
 			{
 				s.center(inP[i]);
 				v1 = (inP[i] - inP[i-1]).norm();
@@ -67,7 +67,7 @@ class cylinder
 		}
 		///returns the direction vector at point idx.
-		stim::vec<T>
+		stim::vec3<T>
 		d(int idx)
 		{
 			if(idx == 0)
@@ -81,15 +81,15 @@ class cylinder
 			else
 			{
 //				return (e[idx+1].P - e[idx].P).norm();
-				stim::vec<float> v1 = (e[idx].P-e[idx-1].P).norm();
-				stim::vec<float> v2 = (e[idx+1].P-e[idx].P).norm();
+				stim::vec3<float> v1 = (e[idx].P-e[idx-1].P).norm();
+				stim::vec3<float> v2 = (e[idx+1].P-e[idx].P).norm();
 				return (v1+v2).norm();			
 			} 
 	//		return e[idx].N;	
 		}
-		stim::vec<T>
+		stim::vec3<T>
 		d(T l, int idx)
 		{
 			if(idx == 0 || idx == e.size()-1)
@@ -144,13 +144,13 @@ class cylinder
 		///constructor to create a cylinder from a set of points, radii, and the number of sides for the cylinder.
 		///@param inP:  Vector of stim vecs composing the points of the centerline.
 		///@param inM:  Vector of stim vecs composing the radii of the centerline.
-		cylinder(std::vector<stim::vec<T> > inP, std::vector<stim::vec<T> > inM){
+		cylinder(std::vector<stim::vec3<T> > inP, std::vector<stim::vec3<T> > inM){
 			init(inP, inM);
 		}
 		///Constructor defines a cylinder with centerline inP and magnitudes of zero
 		///@param inP: Vector of stim vecs composing the points of the centerline
-		cylinder(std::vector< stim::vec<T> > inP){
+		cylinder(std::vector< stim::vec3<T> > inP){
 			std::vector< stim::vec<T> > inM;						//create an array of arbitrary magnitudes
 			stim::vec<T> zero;
@@ -171,12 +171,12 @@ class cylinder
 		///Returns a position vector at the given p-value (p value ranges from 0 to 1).
 		///interpolates the position along the line.
 		///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1).
-		stim::vec<T>
+		stim::vec3<T>
 		p(T pvalue)
 		{
 			if(pvalue < 0.0 || pvalue > 1.0)
 			{
-				return stim::vec<float>(-1,-1,-1);
+				return stim::vec3<float>(-1,-1,-1);
 			}
 			T l = pvalue*L[L.size()-1];
 			int idx = findIdx(l);
@@ -188,7 +188,7 @@ class cylinder
 		///Interpolates the radius along the line.
 		///@param l: the location of the in the cylinder.
 		///@param idx: integer location of the point closest to l but prior to it.
-		stim::vec<T>
+		stim::vec3<T>
 		p(T l, int idx)
 		{
 				T rat = (l-L[idx])/(L[idx+1]-L[idx]);
@@ -252,16 +252,16 @@ class cylinder
 		///in x, y, z coordinates. Theta is in degrees from 0 to 360.
 		///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1).
 		///@param theta: the angle to the point of a circle.
-		stim::vec<T>
+		stim::vec3<T>
 		surf(T pvalue, T theta)
 		{
 			if(pvalue < 0.0 || pvalue > 1.0)
 			{
-				return stim::vec<float>(-1,-1,-1);
+				return stim::vec3<float>(-1,-1,-1);
 			} else {
 			T l = pvalue*L[L.size()-1];
 			int idx = findIdx(l);
-			stim::vec<T> ps = p(l, idx); 
+			stim::vec3<T> ps = p(l, idx); 
 			T m = r(l, idx);
 			s = e[idx];
 			s.center(ps);
@@ -273,10 +273,10 @@ class cylinder
 		///returns a vector of points necessary to create a circle at every position in the fiber.
 		///@param sides: the number of sides of each circle.	
-		std::vector<std::vector<vec<T> > >
+		std::vector<std::vector<vec3<T> > >
 		getPoints(int sides)
 		{
-			std::vector<std::vector <vec<T> > > points;
+			std::vector<std::vector <vec3<T> > > points;
 			points.resize(e.size());
 			for(int i = 0; i < e.size(); i++)
 			{
@@ -293,7 +293,7 @@ class cylinder
 		}
 		/// Allows a point on the centerline to be accessed using bracket notation
-		vec<T> operator[](unsigned int i){
+		vec3<T> operator[](unsigned int i){
 			return e[i].P;
 		}
@@ -309,7 +309,7 @@ class cylinder
 			T M = 0;						//initialize the integral to zero
 			T m0, m1;						//allocate space for both magnitudes in a single segment
-			//vec<T> p0, p1;					//allocate space for both points in a single segment
+			//vec3<T> p0, p1;					//allocate space for both points in a single segment
 			m0 = mags[0][m];				//initialize the first point and magnitude to the first point in the cylinder
 			//p0 = pos[0];
@@ -325,7 +325,7 @@ class cylinder
 				if(p > 1) len = (L[p-1] - L[p-2]);		//calculate the segment length using the L array
 				//add the average magnitude, weighted by the segment length
-				M += (m0 + m1)/2.0 * len;
+				M += (m0 + m1)/(T)2.0 * len;
 				m0 = m1;								//move to the next segment by shifting points
 			}
@@ -345,21 +345,21 @@ class cylinder
 		/// @param spacing is the maximum spacing allowed between sample points
 		cylinder<T> resample(T spacing){
-			std::vector< vec<T> > result;
+			std::vector< vec3<T> > result;
-			vec<T> p0 = e[0].P;								//initialize p0 to the first point on the centerline
-			vec<T> p1;
+			vec3<T> p0 = e[0].P;								//initialize p0 to the first point on the centerline
+			vec3<T> p1;
 			unsigned N = size();							//number of points in the current centerline
 			//for each line segment on the centerline
 			for(unsigned int i = 1; i < N; i++){
 				p1 = e[i].P;								//get the second point in the line segment
-				vec<T> v = p1 - p0;							//calculate the vector between these two points
+				vec3<T> v = p1 - p0;							//calculate the vector between these two points
 				T d = v.len();								//calculate the distance between these two points (length of the line segment)
-				unsigned nsteps = d / spacing+1;		//calculate the number of steps to take along the segment to meet the spacing criteria
-				T stepsize = 1.0 / nsteps;			//calculate the parametric step size between new centerline points
+				size_t nsteps = (size_t)std::ceil(d / spacing);		//calculate the number of steps to take along the segment to meet the spacing criteria
+				T stepsize = (T)1.0 / nsteps;			//calculate the parametric step size between new centerline points
 				//for each step along the line segment
 				for(unsigned s = 0; s < nsteps; s++){