diff --git a/stim/biomodels/centerline.h b/stim/biomodels/centerline.h index 40d56eb..63d8f38 100644 --- a/stim/biomodels/centerline.h +++ b/stim/biomodels/centerline.h @@ -3,7 +3,6 @@ #include #include -//#include namespace stim{ @@ -12,195 +11,123 @@ namespace stim{ * class to describe an interconnected (often biological) network. */ template -class centerline{ +class centerline : public std::vector< stim::vec3 >{ protected: - unsigned int N; //number of points in the fiber - double **c; //centerline (array of double pointers) -// ANNkd_tree* kdt; //kd-tree stores all points in the fiber for fast searching - /// Initialize an empty fiber - void init() - { - N=0; - c=NULL; -// kdt = NULL; - } - - /// Initialize a fiber with N centerline points (all located at [0, 0, 0] with radius 0) - void init(unsigned int n) - { + std::vector L; //stores the integrated length along the fiber (used for parameterization) - N = n; //set the number of points -// kdt = NULL; - c = (double**) malloc(sizeof(double*) * N); //allocate the array pointer + ///Return the normalized direction vector at point i (average of the incoming and outgoing directions) + vec3 d(size_t i) { + if (size() <= 1) return vec3(0, 0, 0); //if there is insufficient information to calculate the direction, return a null vector + if (size() == 2) return (at(1) - at(0)).norm(); //if there are only two points, the direction vector at both is the direction of the line segment + if (i == 0) return (at(1) - at(0)).norm(); //the first direction vector is oriented towards the first line segment + if (i == size() - 1) return at(size() - 1) - at(size() - 2); //the last direction vector is oriented towards the last line segment - for(unsigned int i = 0; i < N; i++) //allocate space for each point - c[i] = (double*) malloc(sizeof(double) * 3); + //all other direction vectors are the average direction of the two joined line segments + return ((at(i) - at(i - 1)).norm() + (at(i + 1) - at(i)).norm()).norm(); } - - /// Copies an existing fiber to the current fiber - - /// @param cpy stores the new copy of the fiber - void copy( const stim::centerline& cpy, bool kd = 0){ - - ///allocate space for the new fiber - init(cpy.N); - - ///copy the points - for(unsigned int i = 0; i < N; i++){ - for(unsigned int d = 0; d < 3; d++) //for each dimension - c[i][d] = cpy.c[i][d]; //copy the coordinate + //initializes the integrated length vector to make parameterization easier, starting with index idx (all previous indices are assumed to be correct) + void update_L(size_t start = 0) { + L.resize(size()); //allocate space for the L array + if (start == 0) { + L[0] = 0; //initialize the length value for the first point to zero (0) + start++; } -// if(kd) -// gen_kdtree(); //generate the kd tree for the new fiber - } - /// generate a KD tree for points on fiber -// void gen_kdtree() -// { -// int n_data = N; //create an array of data points -// ANNpointArray pts = (ANNpointArray)c; //cast the centerline list to an ANNpointArray -// kdt = new ANNkd_tree(pts, n_data, 3); //build a KD tree -// } - - /// find distance between two points - double dist(double* p0, double* p1){ - - double sum = 0; // initialize variables - float v; - for(unsigned int d = 0; d < 3; d++) - { - v = p1[d] - p0[d]; - sum +=v * v; + stim::vec3 d; + for (size_t i = start; i < size(); i++) { //for each line segment in the centerline + d = at(i) - at(i - 1); + L[i] = L[i - 1] + d.len(); //calculate the running length total } - return sqrt(sum); } - /// This function retreives the index for the fiber point closest to q - - /// @param q is a reference point used to find the closest point on the fiber center line -// unsigned int ann( stim::vec q ){ -// ANNidxArray idx = new ANNidx[1]; //variable used to hold the nearest point -// ANNdistArray sq_dist = new ANNdist[1]; //variable used to hold the squared distance to the nearest point -// kdt->annkSearch(q.data(), 1, idx, sq_dist); //search the KD tree for the nearest neighbor -// return *idx; -// } + void init() { + if (size() == 0) return; //return if there aren't any points + update_L(); + } /// Returns a stim::vec representing the point at index i /// @param i is an index of the desired centerline point stim::vec get_vec(unsigned i){ - stim::vec3 r; - r.resize(3); - r[0] = c[i][0]; - r[1] = c[i][1]; - r[2] = c[i][2]; + return std::vector< stim::vec3 >::at(i); + } + + ///finds the index of the point closest to the length l on the lower bound. + ///binary search. + size_t findIdx(T l) { + size_t i = L.size() / 2; + size_t max = L.size() - 1; + size_t min = 0; + while (i > 0 && i < L.size() - 1){ + 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; + } - return r; + ///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 = at(idx); + stim::vec3 v2 = at(idx + 1); + return(v1 + (v2 - v1)*rat); } public: - centerline(){ - init(); - } + using std::vector< stim::vec3 >::at; + using std::vector< stim::vec3 >::size; - /// Copy constructor - centerline(const stim::centerline &obj){ - copy(obj); + centerline() : std::vector< stim::vec3 >() { + init(); } - - //temp constructor for graph visualization - centerline(int n) - { - init(n); + centerline(size_t n) : std::vector< stim::vec3 >(n){ + init(); } - /// Constructor takes a list of stim::vec points, the radius at each point is set to zero - centerline(std::vector< stim::vec > p, bool kd = 0){ - init(p.size()); //initialize the fiber - - //for each point, set the centerline position and radius - for(unsigned int i = 0; i < N; i++){ - - //set the centerline position - for(unsigned int d = 0; d < 3; d++) - c[i][d] = (double) p[i][d]; - - //set the radius - } - //generate a kd tree -// if(kd) -// gen_kdtree(); + //overload the push_back function to update the length vector + void push_back(stim::vec3 p) { + std::vector< stim::vec3 >::push_back(p); + update_L(size() - 1); } - /// constructor takes a list of points - centerline(std::vector< stim::vec3< T > > pos, bool kd = 0){ - init(pos.size()); //initialize the fiber + ///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) return at(0); //return the first element + if (pvalue >= 1.0) return back(); //return the last element - //for each point, set the centerline position and radius - for(unsigned int i = 0; i < N; i++){ - //set the centerline position - for(unsigned int d = 0; d < 3; d++) - c[i][d] = (double) pos[i][d]; - //set the radius - } - - //generate a kd tree - //if(kd) - // gen_kdtree(); + T l = pvalue*L[L.size() - 1]; + int idx = findIdx(l); + return p(l, idx); } - /// Assignment operation - centerline& operator=(const centerline &rhs){ - if(this == &rhs) return *this; //test for and handle self-assignment - copy(rhs); - return *this; - } - - - /// Return the point on the fiber closest to q - /// @param q is the query point used to locate the nearest point on the fiber centerline -// stim::vec nearest(stim::vec q){ -// -// stim::vec temp( (double) q[0], (double) q[1], (double) q[2]); -// -// unsigned int idx = ann(temp); //determine the index of the nearest neighbor -// -// return stim::vec((T) c[idx][0], (T) c[idx][1], (T) c[idx][2]); //return the nearest centerline point -// } - - /// Return the point index on the fiber closest to q - /// @param q is the query point used to locate the nearest point on the fiber centerline -// unsigned int nearest_idx(stim::vec q){ -// -// stim::vec temp((double) q[0], (double) q[1], (double) q[2]); -// -// unsigned int idx = ann(temp); //determine the index of the nearest neighbor -// -// return idx; //return the nearest centerline point index -// } - - /// Returns the fiber centerline as an array of stim::vec points - std::vector< stim::vec > get_centerline(){ - - //create an array of stim vectors - std::vector< stim::vec3 > pts(N); - - //cast each point to a stim::vec, keeping only the position information - for(unsigned int i = 0; i < N; i++) - pts[i] = stim::vec3((T) c[i][0], (T) c[i][1], (T) c[i][2]); - - //return the centerline array - return pts; + ///Return the length of the entire centerline + T length() { + return L.back(); } /// Split the fiber at the specified index. If the index is an end point, only one fiber is returned std::vector< stim::centerline > split(unsigned int idx){ - std::vector< stim::centerline > fl; //create an array to store up to two fibers + std::vector< stim::centerline > fl; //create an array to store up to two fibers + size_t N = size(); //if the index is an end point, only the existing fiber is returned if(idx == 0 || idx == N-1){ @@ -216,123 +143,84 @@ public: fl.resize(2); //set the array size to 2 - fl[0].init(N1); //set the size of each fiber - fl[1].init(N2); + fl[0] = stim::centerline(N1); //set the size of each fiber + fl[1] = stim::centerline(N2); //copy both halves of the fiber unsigned int i, d; //first half - for(i = 0; i < N1; i++){ //for each centerline point - for(d = 0; d < 3; d++) - fl[0].c[i][d] = c[i][d]; //copy each coordinate - } + for(i = 0; i < N1; i++) //for each centerline point + fl[0][i] = std::vector< stim::vec3 >::at(i); + fl[0].init(); //initialize the length vector //second half - for(i = 0; i < N2; i++){ - for(d = 0; d < 3; d++) - fl[1].c[i][d] = c[idx + i][d]; - - } + for(i = 0; i < N2; i++) + fl[1][i] = std::vector< stim::vec3 >::at(idx+i); + fl[1].init(); //initialize the length vector } return fl; //return the array } - /// Calculates the set of fibers resulting from a connection between the current fiber and a fiber f - - /// @param f is the fiber that will be connected to the current fiber -/* std::vector< stim::centerline > connect( stim::centerline &f, double dist){ - - double min_dist; - unsigned int idx0, idx1; - - //go through each point in the query fiber, looking for the indices for the closest points - for(unsigned int i = 0; i < f.n_pts(); i++){ - //Run through all points and find the index with the closest point, then partition the fiber and return two fibers. - - } - - - - } -*/ /// Outputs the fiber as a string std::string str(){ std::stringstream ss; - - //create an iterator for the point list - //typename std::list< point >::iterator i; - for(unsigned int i = 0; i < N; i++){ - ss<<" [ "; - for(unsigned int d = 0; d < 3; d++){ - ss< >::size(); + ss << "---------[" << N << "]---------" << std::endl; + for (size_t i = 0; i < N; i++) + ss << std::vector< stim::vec3 >::at(i) << std::endl; + ss << "--------------------" << std::endl; return ss.str(); } - /// Returns the number of centerline points in the fiber - unsigned int size(){ - return N; - } - - - /// Bracket operator returns the element at index i - - /// @param i is the index of the element to be returned as a stim::vec - stim::vec operator[](unsigned i){ - return get_vec(i); - } /// Back method returns the last point in the fiber - stim::vec back(){ - return get_vec(N-1); + stim::vec3 back(){ + return std::vector< stim::vec3 >::back(); } + ////resample a fiber in the network stim::centerline resample(T spacing) { std::cout<<"fiber::resample()"< v(3); //v-direction vector of the segment - stim::vec p(3); //- intermediate point to be added - stim::vec p1(3); // p1 - starting point of an segment on the fiber, - stim::vec p2(3); // p2 - ending point, + stim::vec3 v; //v-direction vector of the segment + stim::vec3 p; //- intermediate point to be added + stim::vec3 p1; // p1 - starting point of an segment on the fiber, + stim::vec3 p2; // p2 - ending point, double sum=0; //distance summation - std::vector > fiberPositions = centerline(); - std::vector > newPointList; // initialize list of new resampled points on the fiber + + size_t N = size(); + + centerline new_c; // initialize list of new resampled points on the fiber // for each point on the centerline (skip if it is the last point on centerline) - //unsigned int N = fiberPositions.size(); // number of points on the fiber for(unsigned int f=0; f< N-1; f++) - { + { + p1 = at(f); + p2 = at(f+1); + v = p2 - p1; - p1 = fiberPositions[f]; p2 = fiberPositions[f + 1]; v = p2 - p1; - for(unsigned int d = 0; d < 3; d++){ - sum +=v[d] * v[d];} //length of segment-distance between starting and ending point - - T lengthSegment = sqrt(sum); //find Length of the segment as distance between the starting and ending points of the segment - - if(lengthSegment >= spacing) // if length of the segment is greater than standard deviation resample - { - // repeat resampling until accumulated stepsize is equsl to length of the segment - for(T step=0.0; step= spacing){ // if length of the segment is greater than standard deviation resample + + // repeat resampling until accumulated stepsize is equsl to length of the segment + for(T step=0.0; step +#include + +namespace stim{ + +/** This class stores information about a single fiber represented as a set of geometric points + * between two branch or end points. This class is used as a fundamental component of the stim::network + * class to describe an interconnected (often biological) network. + */ +template +class centerline{ + +protected: + unsigned int N; //number of points in the fiber + double **c; //centerline (array of double pointers) + + /// Initialize an empty fiber + void init(){ + N=0; + c=NULL; + } + + /// Initialize a fiber with N centerline points (all located at [0, 0, 0] with radius 0) + void init(unsigned int n){ + + N = n; //set the number of points + c = (double**) malloc(sizeof(double*) * N); //allocate the array pointer + + for(unsigned int i = 0; i < N; i++) //allocate space for each point + c[i] = (double*) malloc(sizeof(double) * 3); + } + + /// Copies an existing fiber to the current fiber + + /// @param cpy stores the new copy of the fiber + void copy( const stim::centerline& cpy, bool kd = 0){ + + ///allocate space for the new fiber + init(cpy.N); + + ///copy the points + for(unsigned int i = 0; i < N; i++){ + for(unsigned int d = 0; d < 3; d++) //for each dimension + c[i][d] = cpy.c[i][d]; //copy the coordinate + } + } + + /// find distance between two points + double dist(double* p0, double* p1){ + + double sum = 0; // initialize variables + float v; + for(unsigned int d = 0; d < 3; d++) + { + v = p1[d] - p0[d]; + sum +=v * v; + } + return sqrt(sum); + } + + /// Returns a stim::vec representing the point at index i + + /// @param i is an index of the desired centerline point + stim::vec get_vec(unsigned i){ + stim::vec3 r; + r.resize(3); + r[0] = c[i][0]; + r[1] = c[i][1]; + r[2] = c[i][2]; + + return r; + } + + +public: + + centerline(){ + init(); + } + + /// Copy constructor + centerline(const stim::centerline &obj){ + copy(obj); + } + + //initialize a centerline with n points + centerline(int n){ + init(n); + } + + /// Constructor takes a list of stim::vec points, the radius at each point is set to zero + centerline(std::vector< stim::vec > p, bool kd = 0){ + init(p.size()); //initialize the fiber + + //for each point, set the centerline position and radius + for(unsigned int i = 0; i < N; i++){ + + //set the centerline position + for(unsigned int d = 0; d < 3; d++) + c[i][d] = (double) p[i][d]; + } + } + + /// constructor takes a list of points + centerline(std::vector< stim::vec3< T > > pos, bool kd = 0){ + init(pos.size()); //initialize the fiber + + //for each point, set the centerline position and radius + for(unsigned int i = 0; i < N; i++){ + //set the centerline position + for(unsigned int d = 0; d < 3; d++) + c[i][d] = (double) pos[i][d]; + } + } + + /// Assignment operation + centerline& operator=(const centerline &rhs){ + if(this == &rhs) return *this; //test for and handle self-assignment + copy(rhs); + return *this; + } + + + /// Returns the fiber centerline as an array of stim::vec points + std::vector< stim::vec > get_centerline(){ + + //create an array of stim vectors + std::vector< stim::vec3 > pts(N); + + //cast each point to a stim::vec, keeping only the position information + for(unsigned int i = 0; i < N; i++) + pts[i] = stim::vec3((T) c[i][0], (T) c[i][1], (T) c[i][2]); + + //return the centerline array + return pts; + } + + /// Split the fiber at the specified index. If the index is an end point, only one fiber is returned + std::vector< stim::centerline > split(unsigned int idx){ + + std::vector< stim::centerline > fl; //create an array to store up to two fibers + + //if the index is an end point, only the existing fiber is returned + if(idx == 0 || idx == N-1){ + fl.resize(1); //set the size of the fiber to 1 + fl[0] = *this; //copy the current fiber + } + + //if the index is not an end point + else{ + + unsigned int N1 = idx + 1; //calculate the size of both fibers + unsigned int N2 = N - idx; + + fl.resize(2); //set the array size to 2 + + fl[0].init(N1); //set the size of each fiber + fl[1].init(N2); + + //copy both halves of the fiber + unsigned int i, d; + + //first half + for(i = 0; i < N1; i++){ //for each centerline point + for(d = 0; d < 3; d++) + fl[0].c[i][d] = c[i][d]; //copy each coordinate + } + + //second half + for(i = 0; i < N2; i++){ + for(d = 0; d < 3; d++) + fl[1].c[i][d] = c[idx + i][d]; + + } + } + + return fl; //return the array + + } + + /// Outputs the fiber as a string + std::string str(){ + std::stringstream ss; + + //create an iterator for the point list + //typename std::list< point >::iterator i; + for(unsigned int i = 0; i < N; i++){ + ss<<" [ "; + for(unsigned int d = 0; d < 3; d++){ + ss< operator[](unsigned i){ + return get_vec(i); + } + + /// Back method returns the last point in the fiber + stim::vec back(){ + return get_vec(N-1); + } + ////resample a fiber in the network + stim::centerline resample(T spacing) + { + std::cout<<"fiber::resample()"< v(3); //v-direction vector of the segment + stim::vec p(3); //- intermediate point to be added + stim::vec p1(3); // p1 - starting point of an segment on the fiber, + stim::vec p2(3); // p2 - ending point, + double sum=0; //distance summation + std::vector > fiberPositions = centerline(); + std::vector > newPointList; // initialize list of new resampled points on the fiber + // for each point on the centerline (skip if it is the last point on centerline) + for(unsigned int f=0; f< N-1; f++) + { + + p1 = fiberPositions[f]; p2 = fiberPositions[f + 1]; v = p2 - p1; + for(unsigned int d = 0; d < 3; d++){ + sum +=v[d] * v[d];} //length of segment-distance between starting and ending point + + T lengthSegment = sqrt(sum); //find Length of the segment as distance between the starting and ending points of the segment + + if(lengthSegment >= spacing) // if length of the segment is greater than standard deviation resample + { + // repeat resampling until accumulated stepsize is equsl to length of the segment + for(T step=0.0; step private: - stim::vec3 Y; + //stim::vec3 Y; + T R; //radius of the circle - CUDA_CALLABLE void + /*CUDA_CALLABLE void init() { Y = U.cross(N).norm(); - } + }*/ public: using stim::plane::n; @@ -42,26 +43,28 @@ public: init(); } - ///create a rectangle given a size and position in Z space. + ///create a circle given a size and position in Z space. ///@param size: size of the rectangle in ND space. ///@param z_pos z coordinate of the rectangle. CUDA_CALLABLE - circle(T size, T z_pos = (T)0) : plane() + circle(T radius, T z_pos = (T)0) : plane(z_pos) { - center(stim::vec3(0,0,z_pos)); - scale(size); - init(); + //center(stim::vec3(0, 0, z_pos)); + //scale(size); + //init(); + R = radius; } ///create a rectangle from a center point, normal ///@param c: x,y,z location of the center. ///@param n: x,y,z direction of the normal. CUDA_CALLABLE - circle(vec3 c, vec3 n = vec3(0,0,1)) : plane() + circle(vec3 c, vec3 n = vec3(0,0,1)) : plane(n, c) { - center(c); - normal(n); - init(); + //center(c); + //normal(n); + //init(); + R = (T)1; } /*///create a rectangle from a center point, normal, and size @@ -84,14 +87,18 @@ public: ///@param n: x,y,z direction of the normal. ///@param u: x,y,z direction for the zero vector (from where the rotation starts) CUDA_CALLABLE - circle(vec3 c, T s, vec3 n = vec3(0,0,1), vec3 u = vec3(1, 0, 0)) : plane() + circle(vec3 c, T r, vec3 n = vec3(0,0,1), vec3 u = vec3(1, 0, 0)) : plane() { - init(); + P = c; + N = n; setU(u); + R = r; + //init(); + //setU(u); // U = u; - center(c); - normal(n); - scale(s); + //center(c); + //normal(n); + //scale(s); } ///scales the circle by a certain factor @@ -99,23 +106,23 @@ public: CUDA_CALLABLE void scale(T factor) { - U *= factor; - Y *= factor; + //U *= factor; + //Y *= factor; + R *= factor; } ///sets the normal for the cirlce ///@param n: x,y,z direction of the normal. CUDA_CALLABLE void - normal(vec3 n) - { - rotate(n, Y); + normal(vec3 n){ + rotate(n); } ///sets the center of the circle. ///@param n: x,y,z location of the center. CUDA_CALLABLE void center(vec3 p){ - this->P = p; + P = p; } ///boolean comparison @@ -128,57 +135,63 @@ public: return false; } + //returns the point in world space corresponding to the polar coordinate (r, theta) + CUDA_CALLABLE stim::vec3 + p(T r, T theta) { + T u = r * cos(theta); //calculate the coordinates in the planar space defined by the circle + T v = r * sin(theta); + + vec3 V = U.cross(N); //calculate the orthogonal vector V + return P + U * u + V * v; //calculate the cartesian coordinate of the specified point + } + + //returns the point in world space corresponding to the value theta at radius R + CUDA_CALLABLE stim::vec3 + p(T theta) { + return p(R, theta); + } + + //get the world space value given the polar coordinates (r, theta) + ///get the world space value given the planar coordinates a, b in [0, 1] - CUDA_CALLABLE stim::vec3 p(T a, T b) + /*CUDA_CALLABLE stim::vec3 p(T a, T b) { stim::vec3 result; vec3 A = this->P - this->U * (T)0.5 - Y * (T)0.5; result = A + this->U * a + Y * b; return result; - } + }*/ ///parenthesis operator returns the world space given rectangular coordinates a and b in [0 1] - CUDA_CALLABLE stim::vec3 operator()(T a, T b) + CUDA_CALLABLE stim::vec3 operator()(T r, T theta) { - return p(a,b); + return p(r, theta); + } + + //parenthesis operator returns the world space coordinate at the edge of the circle given theta + CUDA_CALLABLE stim::vec3 operator()(T theta) { + return p(theta); } ///returns a vector with the points on the initialized circle. ///connecting the points results in a circle. ///@param n: integer for the number of points representing the circle. - std::vector > - getPoints(int n) - { - std::vector > result; - stim::vec3 point; - T x,y; - float step = 360.0/(float) n; - for(float j = 0; j <= 360.0; j += step) - { - y = 0.5*cos(j*stim::TAU/360.0)+0.5; - x = 0.5*sin(j*stim::TAU/360.0)+0.5; - result.push_back(p(x,y)); - } + std::vector< stim::vec3 > points(unsigned n) { + std::vector< stim::vec3 > result(n); //initialize a vector of n points + + float dt = stim::TAU / n; + for (unsigned i = 0; i < n; i++) + result[i] = p(i * dt); //calculate a point on the edge of the circle return result; } - ///returns a vector with the points on the initialized circle. - ///connecting the points results in a circle. - ///@param n: integer for the number of points representing the circle. - CUDA_CALLABLE stim::vec3 - p(T theta) - { - T x,y; - y = 0.5*cos(theta*STIM_TAU/360.0)+0.5; - x = 0.5*sin(theta*STIM_TAU/360.0)+0.5; - return p(x,y); - } + std::string str() const { std::stringstream ss; - ss << "(P=" << P.str() << ", N=" << N.str() << ", U=" << U.str() << ", Y=" << Y.str(); + ss << "r = "< +#include +#include +#include +#include +#include +#include +#include + +namespace stim{ + +template +class circle : plane +{ + +private: + + stim::vec3 Y; + + CUDA_CALLABLE void + init() + { + Y = U.cross(N).norm(); + } + +public: + using stim::plane::n; + using stim::plane::P; + using stim::plane::N; + using stim::plane::U; + using stim::plane::rotate; + using stim::plane::setU; + + ///base constructor + ///@param th value of the angle of the starting point from 0 to 360. + CUDA_CALLABLE + circle() : plane() + { + init(); + } + + ///create a rectangle given a size and position in Z space. + ///@param size: size of the rectangle in ND space. + ///@param z_pos z coordinate of the rectangle. + CUDA_CALLABLE + circle(T size, T z_pos = (T)0) : plane() + { + center(stim::vec3(0,0,z_pos)); + scale(size); + init(); + } + + ///create a rectangle from a center point, normal + ///@param c: x,y,z location of the center. + ///@param n: x,y,z direction of the normal. + CUDA_CALLABLE + circle(vec3 c, vec3 n = vec3(0,0,1)) : plane() + { + center(c); + normal(n); + init(); + } + + /*///create a rectangle from a center point, normal, and size + ///@param c: x,y,z location of the center. + ///@param s: size of the rectangle. + ///@param n: x,y,z direction of the normal. + CUDA_CALLABLE + circle(vec3 c, T s, vec3 n = vec3(0,0,1)) : plane() + { + init(); + center(c); + rotate(n, U, Y); + scale(s); + } + */ + + ///create a rectangle from a center point, normal, and size + ///@param c: x,y,z location of the center. + ///@param s: size of the rectangle. + ///@param n: x,y,z direction of the normal. + ///@param u: x,y,z direction for the zero vector (from where the rotation starts) + CUDA_CALLABLE + circle(vec3 c, T s, vec3 n = vec3(0,0,1), vec3 u = vec3(1, 0, 0)) : plane() + { + init(); + setU(u); +// U = u; + center(c); + normal(n); + scale(s); + } + + ///scales the circle by a certain factor + ///@param factor: the factor by which the dimensions of the shape are scaled. + CUDA_CALLABLE + void scale(T factor) + { + U *= factor; + Y *= factor; + } + + ///sets the normal for the cirlce + ///@param n: x,y,z direction of the normal. + CUDA_CALLABLE void + normal(vec3 n) + { + rotate(n, Y); + } + + ///sets the center of the circle. + ///@param n: x,y,z location of the center. + CUDA_CALLABLE void + center(vec3 p){ + this->P = p; + } + + ///boolean comparison + bool + operator==(const circle & rhs) + { + if(P == rhs.P && U == rhs.U && Y == rhs.Y) + return true; + else + return false; + } + + ///get the world space value given the planar coordinates a, b in [0, 1] + CUDA_CALLABLE stim::vec3 p(T a, T b) + { + stim::vec3 result; + + vec3 A = this->P - this->U * (T)0.5 - Y * (T)0.5; + result = A + this->U * a + Y * b; + return result; + } + + ///parenthesis operator returns the world space given rectangular coordinates a and b in [0 1] + CUDA_CALLABLE stim::vec3 operator()(T a, T b) + { + return p(a,b); + } + + ///returns a vector with the points on the initialized circle. + ///connecting the points results in a circle. + ///@param n: integer for the number of points representing the circle. + std::vector > + getPoints(int n) + { + std::vector > result; + stim::vec3 point; + T x,y; + float step = 360.0/(float) n; + for(float j = 0; j <= 360.0; j += step) + { + y = 0.5*cos(j*stim::TAU/360.0)+0.5; + x = 0.5*sin(j*stim::TAU/360.0)+0.5; + result.push_back(p(x,y)); + } + return result; + } + + ///returns a vector with the points on the initialized circle. + ///connecting the points results in a circle. + ///@param n: integer for the number of points representing the circle. + CUDA_CALLABLE stim::vec3 + p(T theta) + { + T x,y; + y = 0.5*cos(theta*STIM_TAU/360.0)+0.5; + x = 0.5*sin(theta*STIM_TAU/360.0)+0.5; + return p(x,y); + } + + std::string str() const + { + std::stringstream ss; + ss << "(P=" << P.str() << ", N=" << N.str() << ", U=" << U.str() << ", Y=" << Y.str(); + return ss.str(); + } + +}; +} +#endif diff --git a/stim/visualization/cylinder.h b/stim/visualization/cylinder.h index ff3fb82..1f1fb10 100644 --- a/stim/visualization/cylinder.h +++ b/stim/visualization/cylinder.h @@ -3,56 +3,163 @@ #include #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; +class cylinder : public centerline { +protected: + //stim::circle s; //an arbitrary circle + //std::vector > e; //an array of circles that store the centerline + + //std::vector > norms; + std::vector< stim::vec3 > U; //stores the array of U vectors defining the Frenet frame + std::vector< std::vector > M; //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; + using stim::centerline::size; + using stim::centerline::at; + using stim::centerline::findIdx; - ///default init - void - init() - { - + //calculates the U values for each point to initialize the frenet frame + // this function assumes that the centerline has already been set + void init() { + U.resize(size()); //allocate space for the frenet frame vectors + + stim::circle c; //create a circle + stim::vec3 d0, d1; + for (size_t i = 0; i < size() - 1; i++) { //for each line segment in the centerline + c.rotate(d(i)); //rotate the circle to match that normal + U[i] = c.U; //save the U vector from the circle } - - ///Calculate the physical length of the fiber at each point in the fiber. - void - calculateL() - { -/* L.resize(inP.size()); - T temp = (T)0; - L[0] = 0; - for(size_t i = 1; i < L.size(); i++) - { - temp += (inP[i-1] - inP[i]).len(); - L[i] = temp; + U[size() - 1] = c.U; //for the last point, duplicate the final frenet frame vector + } + + ///adds a magnitude to each point in the cylinder + void add_mag(T val = 0) { + if (M.size() == 0) M.resize(size()); //if the magnitude vector isn't initialized, resize it to match the centerline + for (size_t i = 0; i < size(); i++) //for each point + M[i].push_back(val); //add this value to the magnitude vector at each point + } + + //adds a magnitude based on a list of magnitudes for each point + void add_mag(std::vector val) { + if (M.size() == 0) M.resize(size()); //if the magnitude vector isn't initialized, resize it to match the centerline + for (size_t i = 0; i < size(); i++) //for each point + M[i].push_back(val[i]); //add this value to the magnitude vector at each point + } + +public: + std::string str() { + std::stringstream ss; + size_t N = std::vector< stim::vec3 >::size(); + ss << "---------[" << N << "]---------" << std::endl; + for (size_t i = 0; i < N; i++) + ss << std::vector< stim::vec3 >::at(i) << " r = "<< M[i][0]<<" u = "< r) : centerline(c){ + init(); + add_mag(r); + } + + ///Returns magnitude i 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 m(T l, int idx, size_t i = 0) { + T a = (l - L[idx]) / (L[idx + 1] - L[idx]); + T v2 = (M[idx][i] + (M[idx + 1][i] - M[idx][i])*a); + + return v2; + } + + ///Returns the ith magnitude 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 m(T pvalue, size_t i = 0) { + if (pvalue <= 0.0) return M[0][i]; + if (pvalue >= 1.0) return M[size() - 1][i]; + + T l = pvalue*L[L.size() - 1]; + int idx = findIdx(l); + return m(l, idx, i); + } + + ///Returns a circle representing the cylinder cross section at point i + stim::circle circ(size_t i, size_t m = 0) { + return stim::circle(at(i), M[i][m], d(i), U[i]); + } + + ///Returns an OBJ object representing the cylinder with a radial tesselation value of N using magnitude m + stim::obj OBJ(size_t N, size_t m = 0) { + stim::obj out; //create an OBJ object + stim::circle c0, c1; + std::vector< stim::vec3 > p0, p1; //allocate space for the point sets representing the circles bounding each cylinder segment + T u0, u1, v0, v1; //allocate variables to store running texture coordinates + for (size_t i = 1; i < size(); i++) { //for each line segment in the cylinder + c0 = circ(i - 1); //get the two circles bounding the line segment + c1 = circ(i); + + p0 = c0.points(N); //get t points for each of the end caps + p1 = c1.points(N); + + u0 = L[i - 1] / length(); //calculate the texture coordinate (u, v) where u runs along the cylinder length + u1 = L[i] / length(); + + for (size_t n = 1; n < N; n++) { //for each point in the circle + v0 = (double)(n-1) / (double)(N - 1); //v texture coordinate runs around the cylinder + v1 = (double)(n) / (double)(N - 1); + out.Begin(OBJ_FACE); //start a face (quad) + out.TexCoord(u0, v0); + out.Vertex(p0[n - 1][0], p0[n - 1][1], p0[n - 1][2]); //output the points composing a strip of quads wrapping the cylinder segment + out.TexCoord(u1, v0); + out.Vertex(p1[n - 1][0], p1[n - 1][1], p1[n - 1][2]); + + out.TexCoord(u0, v1); + out.Vertex(p1[n][0], p1[n][1], p1[n][2]); + out.TexCoord(u1, v1); + out.Vertex(p0[n][0], p0[n][1], p0[n][2]); + out.End(); } -*/ } + v0 = (double)(N - 2) / (double)(N - 1); //v texture coordinate runs around the cylinder + v1 = 1.0; + out.Begin(OBJ_FACE); + out.TexCoord(u0, v0); + out.Vertex(p0[N - 1][0], p0[N - 1][1], p0[N - 1][2]); //output the points composing a strip of quads wrapping the cylinder segment + out.TexCoord(u1, v0); + out.Vertex(p1[N - 1][0], p1[N - 1][1], p1[N - 1][2]); + + out.TexCoord(u0, v1); + out.Vertex(p1[0][0], p1[0][1], p1[0][2]); + out.TexCoord(u1, v1); + out.Vertex(p0[0][0], p0[0][1], p0[0][2]); + out.End(); + } + return out; + } - ///inits the cylinder from a list of points (std::vector of stim::vec3 --inP) and radii (inM) + + /* + ///inits the cylinder from a list of points (std::vector of stim::vec3 --inP) and magnitudes (inM) void - init(std::vector > inP, std::vector > inM) - { - mags = inM; + init(centerline inP, std::vector< std::vector > inM) { + M = inM; //the magnitude vector can be copied directly + (*this) = inP; //the centerline can be copied to this class directly stim::vec3 v1; stim::vec3 v2; e.resize(inP.size()); @@ -246,13 +353,18 @@ class cylinder init(inP, inM); } + //assignment operator creates a cylinder from a centerline (default radius is zero) + cylinder& operator=(stim::centerline c) { + init(c); + } + ///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). @@ -266,22 +378,7 @@ class cylinder 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. @@ -321,10 +418,7 @@ class cylinder 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. @@ -506,7 +600,7 @@ class cylinder return cylinder(result); - } + }*/ }; diff --git a/stim/visualization/cylinder_dep.h b/stim/visualization/cylinder_dep.h new file mode 100644 index 0000000..26d27bd --- /dev/null +++ b/stim/visualization/cylinder_dep.h @@ -0,0 +1,498 @@ +#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 -- libgit2 0.21.4