diff --git a/stim/biomodels/fiber.h b/stim/biomodels/fiber.h deleted file mode 100644 index c901dd4..0000000 --- a/stim/biomodels/fiber.h +++ /dev/null @@ -1,483 +0,0 @@ -#ifndef STIM_FIBER_H -#define STIM_FIBER_H - -#include -#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 fiber{ - -protected: - unsigned int N; //number of points in the fiber - double **c; //centerline (array of double pointers) - - T* r; // array of fiber radii - ANNkd_tree* kdt; //kd-tree stores all points in the fiber for fast searching - - /// Initialize an empty fiber - void init() - { - kdt = NULL; - c=NULL; - r=NULL; - N=0; - } - - /// 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 - kdt = NULL; - 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); - - r = (T*) malloc(sizeof(T) * N); //allocate space for the radii - } - - /// Copies an existing fiber to the current fiber - - /// @param cpy stores the new copy of the fiber - void copy( const stim::fiber& cpy ){ - - ///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 - - r[i] = cpy.r[i]; //copy the radius - } - - 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; - - } - 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; - } - - /// 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::vec r; - r.resize(3); - r[0] = c[i][0]; - r[1] = c[i][1]; - r[2] = c[i][2]; - - return r; - } - - -public: - - fiber(){ - init(); - } - - /// Copy constructor - fiber(const stim::fiber &obj){ - - copy(obj); - - } - - //temp constructor for graph visualization - fiber(int n) - { - init(n); - } - - /// Constructor takes a list of stim::vec points, the radius at each point is set to zero - fiber(std::vector< stim::vec > p){ - 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 - r[i] = 0; - } - - //generate a kd tree - gen_kdtree(); - } - - /// constructor takes a list of points and radii - fiber(std::vector< stim::vec< T > > pos, std::vector< T > radii){ - 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]; - - //set the radius - r[i] = radii[i]; - } - - //generate a kd tree - gen_kdtree(); - } - - /// constructor takes an array of points and radii - // this function is used when the radii are represented as a stim::vec, - // since this may be easier when importing OBJs - fiber(std::vector< stim::vec > pos, std::vector< stim::vec > radii){ - - init(pos.size()); - - //for each point, set the position and radius - for(unsigned int i = 0; i < N; i++){ - //at(i) = (double*)malloc(sizeof(double) * 3); - for(unsigned int d = 0; d < 3; d++) - c[i][d] = (double) pos[i][d]; - - r[i] = radii[i][(unsigned int)0]; - } - - gen_kdtree(); - } - - /// Assignment operation - fiber& operator=(const fiber &rhs){ - - if(this == &rhs) return *this; //test for and handle self-assignment - - copy(rhs); - } - - /// Calculate the length of the fiber and return it. - double length(){ - - double* p0; - double *p1; - double l = 0; //initialize the length to zero - - //for each point - //typename std::list< point >::iterator i; //create a point iterator - for(unsigned int i = 0; i < N; i++){ //for each point in the fiber - - if(i == 0) //if this is the first point, just store it - p1 = c[0]; - else{ //if this is any other point - p0 = p1; //shift p1->p0 - p1 = c[i]; //set p1 to the new point - l += dist(p0, p1); //add the length of p1 - p0 to the running sum - } - } - - return l; //return the length - } - - /// Calculates the length and average radius of the fiber - - /// @param length is filled with the fiber length - T radius(T& length){ - - double* p0; //temporary variables to store point positions - double* p1; - T r0, r1; //temporary variables to store radii at points - double l; - T r_mean; //temporary variable to store the length and average radius of a fiber segment - double length_sum = 0; //initialize the length to zero - T radius_sum = 0; //initialize the radius sum to zero - - //for each point - //typename std::list< point >::iterator i; //create a point iterator - for(unsigned int i = 0; i < N; i++){ //for each point in the fiber - - if(i == 0){ //if this is the first point, just store it - p1 = c[0]; - r1 = r[0]; - } - else{ //if this is any other point - p0 = p1; //shift p1->p0 and r1->r0 - r0 = r1; - p1 = c[i]; //set p1 to the new point - r1 = r[i]; - - l = dist(p0, p1); //calculate the length of the p0-p1 segment - r_mean = (r0 + r1) / 2; //calculate the average radius of the segment - - radius_sum += r_mean * (T) l; //add the radius scaled by the length to a running sum - length_sum += l; //add the length of p1 - p0 to the running sum - } - } - - length = length_sum; //store the total length - - //if the total length is zero, store a radius of zero - if(length == 0) - return 0; - else - return (T)(radius_sum / length); //return the average radius of the fiber - } - T average_radius() - { - T r_sum = 0.; - for(unsigned int i = 0; i < N; i++) - { - r_sum = r_sum + r[i]; - } - return r_sum/((T) N); - } - - /// Calculates the average radius of the fiber - T radius(){ - T length; - return radius(length); - } - - /// Returns the radius at index idx. - T radius(int idx){ - return r[idx]; - } - - /// 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 > centerline(){ - - //create an array of stim vectors - std::vector< stim::vec > 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::vec((T) c[i][0], (T) c[i][1], (T) c[i][2]); - - //return the centerline array - return pts; - } - - /// Returns the fiber centerline magnitudes as an array of stim::vec points - std::vector< stim::vec > centerlinemag(){ - - //create an array of stim vectors - std::vector< stim::vec > 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::vec(r[i], r[i]);; - - //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::fiber > split(unsigned int idx){ - - std::vector< stim::fiber > 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 - fl[0].r[i] = r[i]; //copy the corresponding radius - } - - //second half - for(i = 0; i < N2; i++){ - for(d = 0; d < 3; d++) - fl[1].c[i][d] = c[idx + i][d]; - fl[1].r[i] = r[idx + i]; - } - } - - 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::fiber > connect( stim::fiber &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< 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) - //unsigned int N = fiberPositions.size(); // number of points on the fiber - 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