cylinder.h 4.73 KB
``````#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()-1);
T temp = (T)0;
for(int i = 0; i < L.size()-1; i++)
{
temp += (pos[i] - pos[i+1]).len();
L[i] = temp;
}
}

///returns the direction vector at point idx.
stim::vec<T>
d(int idx)
{
return (pos[idx] - pos[idx+1]).norm();

}

///returns the total length of the line at index j.
T
getl(int j)
{
for(int i = 0; i < j-1; ++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)
{
int i = pos.size()/2;
while(i > 0 && i < pos.size())
{
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:
///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;
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])));
}

///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)
{
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])));
}

///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)
{
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.
///@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;
T l = pvalue*L[L.size()-1];
int idx = findIdx(l);
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);
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> d = (pos[0] - pos[1]).norm();
s = stim::circle<T>(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
``````