#ifndef STIM_CYLINDER_H #define STIM_CYLINDER_H #include #include #include namespace stim { template class cylinder { private: stim::circle s; std::vector< stim::vec > pos; std::vector< stim::vec > mags; std::vector< T > L; void init() { } void init(std::vector > inP, std::vector > inM) { pos = inP; mags = inM; L.resize(pos.size()-1); for(int i = 0; i < L.size(); i++) { L[i] += (pos[i] - pos[i+1]).len(); } } stim::vec d(int idx) { return (pos[idx] - pos[idx+1]).norm(); } T getl(int j) { for(int i = 0; i < j-1; ++i) { L += (pos[i] -pos[i+1]).len(); } } int findIdx(T l) { int i = pos.size()/2; while(1) { if(L[i] < l) { i = i/2; } else if(L[i] < l && L[i+1] > l) { break; } else { i = i+i/2; } } return i; } public: cylinder() { } ///constructor to create a cylinder from a set of points, radii, and the number of sides for the cylinder. ///The higher the number of sides, the more rectangeles compose the surface of 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 > inP, std::vector > inM) { init(inP, inM); } ///Returns a position vector at the given p-value (p value ranges from 0 to 1). stim::vec p(T pvalue) { T l = pvalue*L[L.size()-1]; int idx = findIdx(l); return (pos[idx] + (pos[idx+1]-pos[idx])*((l-L[idx])/(L[idx+1]- L[idx]))); } stim::vec p(T l, int idx) { 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). T r(T pvalue) { 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]))); } T r(T l, int idx) { return (mags[idx] + (mags[idx+1]-mags[idx])*((l-L[idx])/(L[idx+1]- L[idx]))); } ///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 stim::vec surf(T pvalue, T theta) { T l = pvalue*L[L.size()-1]; int idx = findIdx(l); stim::vec ps = p(l, idx); T m = r(l, idx); stim::vec dr = d(idx); s = stim::circle(ps, m, dr); return(s.p(theta)); } std::vector > > getPoints(int sides) { if(pos.size() < 2) { return; } else { std::vector > > points; points.resize(pos.size()); stim::vec d = (pos[0] - pos[1]).norm(); s = stim::circle(pos[0], mags[0][0], d); points[0] = s.getPoints(sides); for(int i = 1; i < pos.size(); i++) { d = (pos[i] - pos[i-1]).norm(); s.center(pos[i]); s.normal(d); s.scale(mags[i][0]/mags[i-1][0], mags[i][0]/mags[i-1][0]); points[i] = s.getPoints(sides); } return points; } } }; } #endif