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


namespace stim
{
template<typename T>
class cylinder
 : public centerline<T>
{
	private:
		stim::circle<T> s;							//an arbitrary circle
		std::vector<stim::circle<T> > e;			//an array of circles that store the centerline

		std::vector<stim::vec3<T> > norms;
		std::vector<stim::vec<T> > Us;
		std::vector<stim::vec<T> > mags;			//stores a list of magnitudes for each point in the centerline (assuming mags[0] is the radius)
		std::vector< T > L;							//length of the cylinder at each position (pre-integration)


		using stim::centerline<T>::c;
		using stim::centerline<T>::N;
	
		///default init	
		void
		init()
		{

		}

		///inits the cylinder from a list of points (std::vector of stim::vec3 --inP) and radii (inM)
		void
		init(std::vector<stim::vec3<T> > inP, std::vector<stim::vec<T> > inM)
		{
			mags = inM;
			stim::vec3<float> v1;
			stim::vec3<float> v2;
			e.resize(inP.size());

			norms.resize(inP.size());
			Us.resize(inP.size());

			if(inP.size() < 2)
				return;

			//calculate each L.
			L.resize(inP.size());						//the number of precomputed lengths will equal the number of points
			T temp = (T)0;								//length up to that point
			L[0] = temp;
			for(size_t i = 1; i < L.size(); i++)
			{
				temp += (inP[i-1] - inP[i]).len();
				L[i] = temp;
			}

			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));
			norms[0] = s.N;
			e[0] = s;
			Us[0] = e[0].U;
			for(size_t i = 1; i < inP.size()-1; i++)
			{
				s.center(inP[i]);
				v1 = (inP[i] - inP[i-1]).norm();
				v2 = (inP[i+1] - inP[i]).norm();
				dr = (v1+v2).norm();
				s.normal(dr);

				norms[i] = s.N;

				s.scale(inM[i][0]/inM[i-1][0]);
				e[i] = s;
				Us[i] = e[i].U;
			}
			
			int j = inP.size()-1;
			s.center(inP[j]);
			dr = (inP[j] - inP[j-1]).norm();
			s.normal(dr);

			norms[j] = s.N;

			s.scale(inM[j][0]/inM[j-1][0]);
			e[j] = s;
			Us[j] = e[j].U;
		}
		
		///returns the direction vector at point idx.
		stim::vec3<T>
		d(int idx)
		{
			if(idx == 0)
			{
				stim::vec3<T> temp(
						c[idx+1][0]-c[idx][0],	
						c[idx+1][1]-c[idx][1],	
						c[idx+1][2]-c[idx][2]
						);	
//				return (e[idx+1].P - e[idx].P).norm();
				return (temp.norm());
			}
			else if(idx == N-1)
			{
				stim::vec3<T> temp(
						c[idx][0]-c[idx+1][0],	
						c[idx][1]-c[idx+1][1],	
						c[idx][2]-c[idx+1][2]
						);	
			//	return (e[idx].P - e[idx-1].P).norm();
				return (temp.norm());
			}
			else
			{
//				return (e[idx+1].P - e[idx].P).norm();
//				stim::vec3<float> v1 = (e[idx].P-e[idx-1].P).norm();
				stim::vec3<T> v1(
						c[idx][0]-c[idx-1][0],	
						c[idx][1]-c[idx-1][1],	
						c[idx][2]-c[idx-1][2]
						);	
						
//				stim::vec3<float> v2 = (e[idx+1].P-e[idx].P).norm();
				stim::vec3<T> v2(
						c[idx+1][0]-c[idx][0],	
						c[idx+1][1]-c[idx][1],	
						c[idx+1][2]-c[idx][2]
						);
					
				return (v1.norm()+v2.norm()).norm();			
			} 
	//		return e[idx].N;	

		}

		stim::vec3<T>
		d(T l, int idx)
		{
			if(idx == 0 || idx == N-1)
			{
				return norms[idx];
//				return e[idx].N;
			}
			else
			{
				
				T rat = (l-L[idx])/(L[idx+1]-L[idx]);
				return(	norms[idx] + (norms[idx+1] - norms[idx])*rat);
//				return(	e[idx].N + (e[idx+1].N - e[idx].N)*rat);
			} 
		}


		///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;
			while(i > 0 && i < L.size()-1)
			{
//				std::cerr << "Trying " << i << std::endl;
//				std::cerr << "l is " << l << ", L[" << i << "]" << L[i] << std::endl;
				if(l < L[i])
				{
					max = i;
					i = min+(max-min)/2;
				}
				else if(L[i] <= l && L[i+1] >= l)
				{
					break;
				}
				else
				{
					min = i;
					i = min+(max-min)/2;
				}
			}
			return i;
		}

	public:
		///default constructor
		cylinder()
		// : centerline<T>()
		{

		}

		///constructor to create a cylinder from a set of points, radii, and the number of sides for the cylinder.
		///@param inP:  Vector of stim vec3 composing the points of the centerline.
		///@param inM:  Vector of stim vecs composing the radii of the centerline.
		cylinder(std::vector<stim::vec3<T> > inP, std::vector<stim::vec<T> > inM)
			: centerline<T>(inP)
		{
			init(inP, inM);
		}

		///constructor to create a cylinder from a set of points, radii, and the number of sides for the cylinder.
		///@param inP:  Vector of stim vec3 composing the points of the centerline.
		///@param inM:  Vector of stim vecs composing the radii of the centerline.
		cylinder(std::vector<stim::vec3<T> > inP, std::vector< T > inM)
			: centerline<T>(inP)
		{
			std::vector<stim::vec<T> > temp;
			stim::vec<T> zero(0.0,0.0);
			temp.resize(inM.size(), zero);
			for(int i = 0; i < inM.size(); i++)
				temp[i][0] = inM[i];
			init(inP, temp);
		}


		///Constructor defines a cylinder with centerline inP and magnitudes of zero
		///@param inP: Vector of stim vec3 composing the points of the centerline
		cylinder(std::vector< stim::vec3<T> > inP)
			: centerline<T>(inP)
		{
			std::vector< stim::vec<T> > inM;						//create an array of arbitrary magnitudes

			stim::vec<T> zero;
			zero.push_back(0);

			inM.resize(inP.size(), zero);								//initialize the magnitude values to zero
			init(inP, inM);
		}

		///Returns the number of points on the cylinder centerline

		unsigned int size(){
			return N;
		}


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

//			return(	e[idx].P + (e[idx+1].P-e[idx].P)*rat);
			return( v1 + (v2 - v1)*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::vec3<T>
		p(T l, int idx)
		{
				T rat = (l-L[idx])/(L[idx+1]-L[idx]);
				stim::vec3<T> v1(
						c[idx][0],	
						c[idx][1],	
						c[idx][2]
						);	
						
				stim::vec3<T> v2(
						c[idx+1][0],	
						c[idx+1][1],	
						c[idx+1][2]
						);
//			return(	e[idx].P + (e[idx+1].P-e[idx].P)*rat);
			return(	v1 + (v2-v1)*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){
				std::cerr<<"Error, value "<<pvalue<<" is outside of [0 1]."<<std::endl;
				exit(1);
			}
			T l = pvalue*L[L.size()-1];
			int idx = findIdx(l);
			return (r(l,idx));
/*			T rat = (l-L[idx])/(L[idx+1]-L[idx]);
//			return (mags[idx][0] + (mags[idx+1][0]-mags[idx][0])*rat);
			return (e[idx].U.len() + (e[idx+1].U.len() - e[idx].U.len())*rat);
*/		}

		///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]);
			T v1 = (e[idx].U.len() + (e[idx+1].U.len() - e[idx].U.len())*rat);
			T v3 = (Us[idx].len() + (Us[idx+1].len() - Us[idx].len())*rat);
			T v2 = (mags[idx][0] + (mags[idx+1][0]-mags[idx][0])*rat);
//			std::cout << (float)v1 = (float) v2 << std::endl;
			if(abs(v3 - v1) >= 10e-6)
			{
				std::cout << "-------------------------" << std::endl;
				std::cout << e[idx].str() << std::endl << std::endl;
				std::cout << Us[idx].str() << std::endl;
				std::cout << (float)v1 - (float) v2 << std::endl;
				std::cout << "failed" << std::endl;
			}
//			std::cout << e[idx].U.len() << " " << mags[idx][0] << std::endl;
//			std::cout << v2 << std::endl;
			return(v2);
//			return (mags[idx][0] + (mags[idx+1][0]-mags[idx][0])*rat);
	//	(
		}

		///	Returns the magnitude at the given index
		///	@param i is the index of the desired point
		/// @param m is the index of the magnitude value
		T ri(unsigned i, unsigned m = 0){
			return mags[i][m];
		}

		/// Adds a new magnitude value to all points
		/// @param m is the starting value for the new magnitude
		void add_mag(T m = 0){
			for(unsigned int p = 0; p < N; p++)
				mags[p].push_back(m);
		}

		/// Sets a magnitude value
		/// @param val is the new value for the magnitude
		/// @param p is the point index for the magnitude to be set
		/// @param m is the index for the magnitude
		void set_mag(T val, unsigned p, unsigned m = 0){
			mags[p][m] = val;
		}

		/// Returns the number of magnitude values at each point
		unsigned nmags(){
			return mags[0].size();
		}

		///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::vec3<T>
		surf(T pvalue, T theta)
		{
			if(pvalue < 0.0 || pvalue > 1.0)
			{
				return stim::vec3<float>(-1,-1,-1);
			} else {
			T l = pvalue*L[L.size()-1];
			int idx = findIdx(l);
			stim::vec3<T> ps = p(l, idx); 
			T m = r(l, idx);
			s = e[idx];
			s.center(ps);
			s.normal(d(l, idx));
			s.scale(m/e[idx].U.len());
			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<vec3<T> > >
		getPoints(int sides)
		{
			std::vector<std::vector <vec3<T> > > points;
			points.resize(N);
			for(int i = 0; i < N; i++)
			{
				points[i] = e[i].getPoints(sides);
			}
			return points;
		}

		///returns the total length of the line at index j.
		T
		getl(int j)
		{
			return (L[j]);
		}
		/// Allows a point on the centerline to be accessed using bracket notation

		vec3<T> operator[](unsigned int i){
			return e[i].P;
		}

		/// Returns the total length of the cylinder centerline
		T length(){
			return L.back();
		}

		/// Integrates a magnitude value along the cylinder.
		/// @param m is the magnitude value to be integrated (this is usually the radius)
		T integrate(unsigned m = 0){

			T M = 0;						//initialize the integral to zero
			T m0, m1;						//allocate space for both magnitudes 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];

			T len = L[0];						//allocate space for the segment length

			//for every consecutive point in the cylinder
			for(unsigned p = 1; p < N; p++){

				//p1 = pos[p];							//get the position and magnitude for the next point
				m1 = mags[p][m];

				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)/(T)2.0 * len;

				m0 = m1;								//move to the next segment by shifting points
			}
			return M;			//return the integral
		}

		/// Averages a magnitude value across the cylinder
		/// @param m is the magnitude value to be averaged (this is usually the radius)
		T average(unsigned m = 0){			

			//return the average magnitude
			return integrate(m) / L.back();
		}

		/// Resamples the cylinder to provide a maximum distance of "spacing" between centerline points. All current
		///		centerline points are guaranteed to exist in the new cylinder
		/// @param spacing is the maximum spacing allowed between sample points
		cylinder<T> resample(T spacing){

			std::vector< vec3<T> > result;

			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

				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)

				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++){
					T alpha = stepsize * s;							//calculate the fraction of the distance along the line segment covered
					result.push_back(p0 + alpha * v);				//push the point at alpha position along the line segment
				}

				p0 = p1;											//shift the points to move to the next line segment
			}

			result.push_back(e[size() - 1].P);						//push the last point in the centerline

			return cylinder<T>(result);

		}

		
};

}
#endif