#ifndef STIM_CYLINDER_H #define STIM_CYLINDER_H #include #include #include namespace stim { template class cylinder : public centerline { private: stim::circle s; //an arbitrary circle std::vector > e; //an array of circles that store the centerline std::vector > norms; std::vector > Us; std::vector > 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::c; using stim::centerline::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 > inP, std::vector > inM) { mags = inM; stim::vec3 v1; stim::vec3 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 dr = (inP[1] - inP[0]).norm(); s = stim::circle(inP[0], inM[0][0], dr, stim::vec3(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 d(int idx) { if(idx == 0) { stim::vec3 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 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 v1 = (e[idx].P-e[idx-1].P).norm(); stim::vec3 v1( c[idx][0]-c[idx-1][0], c[idx][1]-c[idx-1][1], c[idx][2]-c[idx-1][2] ); // stim::vec3 v2 = (e[idx+1].P-e[idx].P).norm(); stim::vec3 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 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() { } ///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 > inP, std::vector > inM) : centerline(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 > inP, std::vector< T > inM) : centerline(inP) { std::vector > temp; stim::vec 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 > inP) : centerline(inP) { std::vector< stim::vec > inM; //create an array of arbitrary magnitudes stim::vec 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 p(T pvalue) { if(pvalue < 0.0 || pvalue > 1.0) { return stim::vec3(-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 v1( c[idx][0], c[idx][1], c[idx][2] ); stim::vec3 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 p(T l, int idx) { T rat = (l-L[idx])/(L[idx+1]-L[idx]); stim::vec3 v1( c[idx][0], c[idx][1], c[idx][2] ); stim::vec3 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 "<= 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 surf(T pvalue, T theta) { if(pvalue < 0.0 || pvalue > 1.0) { return stim::vec3(-1,-1,-1); } else { T l = pvalue*L[L.size()-1]; int idx = findIdx(l); stim::vec3 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 > > getPoints(int sides) { std::vector > > 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 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 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 resample(T spacing){ std::vector< vec3 > result; vec3 p0 = e[0].P; //initialize p0 to the first point on the centerline vec3 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 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(result); } }; } #endif