codebase / stimlib

#ifndef STIM_CYLINDER_H
#define STIM_CYLINDER_H
#include <iostream>
#include <stim/math/circle.h>
#include <stim/math/vector.h>


namespace stim
{
template<typename T>
class cylinder
{
	private:
		stim::circle<T> s;			//an arbitrary circle
		std::vector< stim::vec<T> > pos;	//positions of the cylinder.
		std::vector< stim::vec<T> > mags;	//radii at each position
		std::vector< T > L;			//length of the cylinder at each position.
	
		///default init	
		void
		init()
		{

		}

		///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)
		{
			pos = inP;
			mags = inM;

			//calculate each L.
			L.resize(pos.size());
			T temp = (T)0;
			L[0] = 0;
			for(int i = 1; i < L.size(); i++)
			{
				temp += (pos[i-1] - pos[i]).len();
				L[i] = temp;
			}
		}
		
		///returns the direction vector at point idx.
		stim::vec<T>
		d(int idx)
		{
			if(idx == 0)
			{
				return (pos[idx+1] - pos[idx]).norm();
			}
			else if(idx == pos.size()-1)
			{
				return (pos[idx] - pos[idx-1]).norm();
			}
			else
			{
				stim::vec<float> v1 = (pos[idx]-pos[idx-1]).norm();
				stim::vec<float> v2 = (pos[idx+1]-pos[idx]).norm();
				return (v1+v2).norm();			
			}
			
		}

		///returns the total length of the line at index j.
		T
		getl(int j)
		{
			T temp = (T) 0;
			for(int i = 0; i < j; ++i)
			{
				temp += (pos[i] - pos[i+1]).len();
				L[i] = temp;
			}
		}

		///finds the index of the point closest to the length l on the lower bound.
		///binary search.
		int
		findIdx(T l)
		{
			unsigned int i = L.size()/2;
			unsigned int max = L.size()-1;
			unsigned int min = 0;
//			std::cerr << "Index initially: " << i << std::endl;
//			std::cerr << "l initially: " << l << std::endl;
			while(i > 0 && i < L.size()-1)
			{
//				std::cerr << "L[i] =  " << L[i] << " i= "  << i << " maxval= " << L.size()-1 << std::endl;
				if(l < L[i])
				{
					max = i;
//				std::cerr << l << " < " << L[i] << "so max is set to " << i << " and the new i is " << min+(max-min)/2 << std::endl;
					i = min+(max-min)/2;
				}
				else if(L[i] <= l && L[i+1] >= l)
				{
//				std::cerr << L[i] << " < " << l << " < " << L[i+1] << std::endl;
					break;
				}
				else
				{
					min = i;
//				std::cerr << l << " > " << L[i] << "so min is set to " << i << " and the new i is " << min+(max-min)/2 << std::endl;
					i = min+(max-min)/2;
				}
			}
			return i;
		}

	public:
		///default constructor
		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)
		{
			init(inP, inM);
		}


		///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>
		p(T pvalue)
		{
			if(pvalue < 0.0 || pvalue > 1.0)
			{
				return stim::vec<float>(-1,-1,-1);
			}
			T l = pvalue*L[L.size()-1];
			int idx = findIdx(l);
				T rat = (l-L[idx])/(L[idx+1]-L[idx]);
			return(	pos[idx] + (pos[idx+1]-pos[idx])*rat);
		}

		///Returns a position vector at the given length into the fiber (based on the pvalue).
		///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>
		p(T l, int idx)
		{
				T rat = (l-L[idx])/(L[idx+1]-L[idx]);
			return(	pos[idx] + (pos[idx+1]-pos[idx])*rat);
//			return(
//			return (pos[idx] + (pos[idx+1]-pos[idx])*((l-L[idx])/(L[idx+1]- L[idx])));
		}

		///Returns a radius at the given p-value (p value ranges from 0 to 1).
		///interpolates the radius along the line.
		///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1).
		T
		r(T pvalue)
		{
			if(pvalue < 0.0 || pvalue > 1.0)
				return;
			T l = pvalue*L[L.size()-1];
			int idx = findIdx(l);
			return (mags[idx] + (mags[idx+1]-mags[idx])*((l-L[idx])/(L[idx+1]- L[idx])))[0];
		}

		///Returns a radius at the given length into the fiber (based on the pvalue).
		///Interpolates the position 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.
		T
		r(T l, int idx)
		{
				T rat = (l-L[idx])/(L[idx+1]-L[idx]);
			return(	mags[idx] + (mags[idx+1]-mags[idx])*rat)[0];
//			return (mags[idx] + (mags[idx+1]-mags[idx])*((l-L[idx])/(L[idx+1]- L[idx])))[0];
		}


		///returns the position of the point with a given pvalue and theta on the surface
		///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>
		surf(T pvalue, T theta)
		{
			if(pvalue < 0.0 || pvalue > 1.0)
			{
				return stim::vec<float>(-1,-1,-1);
			} else {
//			std::cout << "AAAAAAAAAAAAAAAAAAAAAAA             " << pvalue*L[L.size()-1] << " " << L[L.size()-1] << std::endl;
			T l = pvalue*L[L.size()-1];
//			std::cerr << "The l value is: " << l << std::endl;
			int idx = findIdx(l);
//			std::cerr << "Index: " << idx << std::endl;
			stim::vec<T> ps = p(l, idx); 
			T m = r(l, idx);
			stim::vec<T> dr = d(idx);
			s = stim::circle<T>(ps, m, dr, stim::vec<T>(1,0,0));
			return(s.p(theta));
			}
		}

		///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> > >
		getPoints(int sides)
		{
			if(pos.size() < 2)
			{
				return;
			} else {
				std::vector<std::vector <vec<T> > > points;
				points.resize(pos.size());
				stim::vec<T> dr = (pos[1] - pos[0]).norm();
				s = stim::circle<T>(pos[0], mags[0][0], dr, stim::vec<T>(1,0,0));
				points[0] = s.getPoints(sides);
				for(int i = 1; i < pos.size(); i++)
				{
					dr = d(i);
//					dr = (pos[i-1] - pos[i]).norm();
					
//					s = stim::circle<T>(pos[i], mags[i][0], dr, stim::vec<T>(1,0,0));
					s.center(pos[i]);
					s.normal(dr);
//					s.scale(mags[i][0]/mags[i-1][0], mags[i][0]/mags[i-1][0]);
					s.scale(mags[i][0]/mags[i-1][0]);
					points[i] = s.getPoints(sides);
				}
				return points;
			}
		}

		void
		print(int idx)
		{
			std::cout << d(idx) << std::endl;
		}
		
};

}
#endif