diff --git a/stim/biomodels/fiber.h b/stim/biomodels/fiber.h new file mode 100644 index 0000000..2ee5a96 --- /dev/null +++ b/stim/biomodels/fiber.h @@ -0,0 +1,502 @@ +#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 +// // fibers this fiber intersects. + 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 + } + + void gen_kdtree(){ + + //create an array of data points + int n_data = N; + + ANNpointArray pts = (ANNpointArray)c; //cast the centerline list to an ANNpointArray + + kdt = new ANNkd_tree(pts, n_data, 3); //build a KD tree + } + + + double dist(double* p0, double* p1){ + + double sum = 0; + 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(); + } + + ///// This constructor takes a list of points and radii + //fiber(std::list< stim::vec< T > > pos, std::list< T > radii){ + + // init(pos.size()); //initialize the array size + + // //create an iterator for each list + // typename std::list< stim::vec< T > >::iterator pi = pos.begin(); + // typename std::list< T >::iterator ri = radii.begin(); + + // //create a counter for indexing into the fiber array + // unsigned int i = 0; + + // //for each point, set the position and radius + // for(pi = pos.begin(), ri = radii.begin(); pi != pos.end(); pi++, ri++){ + + // //set the centerline position + // for(unsigned int d = 0; d < 3; d++) + // c[i][d] = (double) (*pi)[d]; + // + // r[i] = *ri; //set the radius + // i++; //increment the array index + // } + + // gen_kdtree(); //generate a kd tree + + //} + + /// 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]; + } + /// get index of a node on a fiber + // by matching the node on fiber to already set vertices (both strings) + // used in obj file conversion + int + getIndexes(std::string* input, std::string searched, int sizeV) + { + int result = 0; + for (int i = 0; i < sizeV; i++) + { + if (input[i] == searched) + { + result = i + 1; + } + } + return result; + } + // strObj returns a string of fiber indices corresponding to vectors of positions in the fiber including intermediate nodes + std::string + strObj(std::string* strArray, int sizeV) + { + std::stringstream ss; + std::stringstream oss; + for(unsigned int i = 0; i < N; i++) + { + ss.str(std::string()); + for(unsigned int d = 0; d < 3; d++) + { + ss< 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<::N<<")\tl = "< - class t_node{ - - public: - - unsigned int id; - - //lists of edge indices for capillaries - //the "in" and "out" just indicate how the geometry is defined: - // edges in the "in" list are geometrically oriented such that the terminal node is last - // edges in the "out" list are geometrically oriented such that the terminal node is first - std::list< fiber_i > in; //edge indices for incoming capillaries - std::list< fiber_i > out; //edge indices for outgoing capillaries - std::string str(){ - - std::stringstream ss; - - ss<::iterator f; - - for(f = in.begin(); f != in.end(); f++){ - - if(f != in.begin()) - ss<<", "; - ss<<(*f)->n[0]->id; - } - - //if there are nodes in both lists, separate them by a comma - if(out.size() > 0 && in.size() > 0) - ss<<", "; - - for(f = out.begin(); f != out.end(); f++){ + }; + + ///Node class that stores the physical position of the node as well as the edges it is connected to (edges that connect to it), As well as any additional data necessary. + class vertex : public stim::vec + { + public: + //std::vector edges; //indices of edges connected to this node. + std::vector e[2]; //indices of edges going out (e[0]) and coming in (e[1]) + //stim::vec p; //position of this node in physical space. + + //constructor takes a stim::vec + vertex(stim::vec p) : stim::vec(p){} + + /// Output the vertex information as a string + std::string str(){ + std::stringstream ss; + ss<<"\t(x, y, z) = "<::str(); + + if(e[0].size() > 0){ + ss<<"\t> "; + for(unsigned int o = 0; o < e[0].size(); o++) + ss< 0){ + ss<<"\t< "; + for(unsigned int i = 0; i < e[1].size(); i++) + ss<n[1]->id; + return ss.str(); } - - - return ss.str(); - - - - } + }; + private: -//---------------NETWORK CLASS----------------------------- - -protected: - - //list of terminal nodes - std::list N; - - //list of fibers - std::list F; - - /// Sets a unique ID for each terminal node and fiber - void set_names(){ - - unsigned int i; - - i = 0; - for(t_node_i ti = N.begin(); ti != N.end(); ti++) - ti->id = i++; + std::vector E; //list of edges + std::vector V; //list of vertices. + + public: - i = 0; - for(fiber_i fi = F.begin(); fi != F.end(); fi++) - fi->id = i++; + ///Returns the number of edges in the network. + unsigned int edges(){ + return E.size(); } -public: - - std::string str(){ - - - //assign names to elements of the network - set_names(); - - //create a stringstream for output - std::stringstream ss; - - //output the nodes - ss<<"Nodes-----------------------------------------------"<str()<radius(length); - - //output the IDs of the terminal nodes - ss<n[0]->id<<" -- "<n[1]->id<<": length = "<& n0, std::vector& n1, std::vector& length, std::vector& radius){ - - //determine the number of fibers in the network - unsigned int nF = F.size(); - - //allocate the necessary space to store the fiber information - n0.resize(nF); - n1.resize(nF); - length.resize(nF); - radius.resize(nF); - - //assign names (identifiers) to the network components - set_names(); - - //fill the arrays - unsigned int i = 0; - T l, r; - for(fiber_i f = F.begin(); f != F.end(); f++){ - n0[i] = f->n[0]->id; //get the identifiers for the first and second nodes for the current fiber - n1[i] = f->n[1]->id; - - r = f->radius(l); //get the length and radius of the capillary (calculated at the same time) - - radius[i] = r; //store the radius in the output array - length[i] = l; //store the length in the output array - - i++; //increment the array index + ss<<"Edges ("< +#include +#include +#include + +namespace stim{ + +/** This class provides an interface for dealing with biological networks. + * It takes the following aspects into account: + * 1) Network geometry and centerlines + * 2) Network connectivity (a graph structure can be extracted) + * 3) Network surface structure (the surface is represented as a triangular mesh and referenced to the centerline) + */ + +template +class network{ + + //-------------------HELPER CLASSES----------------------- + /// Stores information about a geometric point on the network centerline (including point position and radius) + //template + class point : public stim::vec{ + + public: + T r; + + point() : stim::vec(){} + + //casting constructor + point(stim::vec rhs) : stim::vec(rhs){} + }; + + //template + class t_node; + class fiber; + + //create typedefs for the iterators to simplify the network code + typedef typename std::list< fiber >::iterator fiber_i; + typedef typename std::list< t_node >::iterator t_node_i; + + /// Stores information about a single capillary (a length of vessel between two branch or end points) + //template + class fiber : public std::list< point >{ + + using std::list< point >::begin; + using std::list< point >::end; + using std::list< point >::size; + + + public: + //std::list< point > P; //geometric point positions + + typename std::list< t_node >::iterator n[2]; //indices to terminal nodes + unsigned int id; + + public: + + /// Calculate the length of the fiber and return it. + T length(){ + + point p0, p1; + T l = 0; //initialize the length to zero + + //for each point + typename std::list< point >::iterator i; //create a point iterator + for(i = begin(); i != end(); i++){ //for each point in the fiber + + if(i == begin()) //if this is the first point, just store it + p1 = *i; + else{ //if this is any other point + p0 = p1; //shift p1->p0 + p1 = *i; //set p1 to the new point + l += (p1 - p0).len(); //add the length of p1 - p0 to the running sum + } + } + + return l; //return the length + } + + T radius(T& length){ + + point p0, p1; //temporary variables to store point positions + T r0, r1; //temporary variables to store radii at points + T l, r; //temporary variable to store the length and average radius of a fiber segment + T 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(i = begin(); i != end(); i++){ //for each point in the fiber + + if(i == begin()){ //if this is the first point, just store it + p1 = *i; + r1 = i->r; + } + else{ //if this is any other point + p0 = p1; //shift p1->p0 and r1->r0 + r0 = r1; + p1 = *i; //set p1 to the new point + r1 = i->r; //and r1 + + l = (p1 - p0).len(); //calculate the length of the p0-p1 segment + r = (r0 + r1) / 2; //calculate the average radius of the segment + + radius_sum += r * 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 radius_sum / length; //return the average radius of the fiber + } + + std::vector< stim::vec > geometry(){ + + std::vector< stim::vec > result; //create an array to store the fiber geometry + result.resize( size() ); //pre-allocate the array + + typename std::list< point >::iterator p; //create a list iterator + unsigned int pi = 0; //create an index into the result array + + //for each geometric point on the fiber centerline + for(p = begin(); p != end(); p++){ + result[pi] = *p; + pi++; + } + + return result; //return the geometry array + + } + + std::string str(){ + std::stringstream ss; + + //create an iterator for the point list + typename std::list::iterator i; + for(i = begin(); i != end(); i++){ + ss<str()<<" r = "<r< + class t_node{ + + public: + + unsigned int id; + + //lists of edge indices for capillaries + //the "in" and "out" just indicate how the geometry is defined: + // edges in the "in" list are geometrically oriented such that the terminal node is last + // edges in the "out" list are geometrically oriented such that the terminal node is first + std::list< fiber_i > in; //edge indices for incoming capillaries + std::list< fiber_i > out; //edge indices for outgoing capillaries + + std::string str(){ + + std::stringstream ss; + + ss<::iterator f; + + for(f = in.begin(); f != in.end(); f++){ + + if(f != in.begin()) + ss<<", "; + ss<<(*f)->n[0]->id; + } + + //if there are nodes in both lists, separate them by a comma + if(out.size() > 0 && in.size() > 0) + ss<<", "; + + for(f = out.begin(); f != out.end(); f++){ + + if(f != out.begin()) + ss<<", "; + ss<<(*f)->n[1]->id; + } + + + return ss.str(); + + + + } + }; + + +//---------------NETWORK CLASS----------------------------- + +protected: + + //list of terminal nodes + std::list N; + + //list of fibers + std::list F; + + /// Sets a unique ID for each terminal node and fiber + void set_names(){ + + unsigned int i; + + i = 0; + for(t_node_i ti = N.begin(); ti != N.end(); ti++) + ti->id = i++; + + i = 0; + for(fiber_i fi = F.begin(); fi != F.end(); fi++) + fi->id = i++; + } + +public: + + std::string str(){ + + + //assign names to elements of the network + set_names(); + + //create a stringstream for output + std::stringstream ss; + + //output the nodes + ss<<"Nodes-----------------------------------------------"<str()<radius(length); + + //output the IDs of the terminal nodes + ss<n[0]->id<<" -- "<n[1]->id<<": length = "< 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< E; //list of pointers to edges. - std::vector V; //list of nodes. - std::vector< stim::vec > allVerticesList; // all nodes before sampling - std::vector > allVerticesAfterSampling ; //all vertices after sampling - - ///Returns the number of edges in the network. - unsigned int - sizeE() - { - return E.size(); - } - - ///Returns the number of nodes in the network. - unsigned int - sizeV() - { - return V.size(); - } - unsigned int - sizeAfterSamplingV() - { - return allVerticesAfterSampling.size(); - } -/* //adds an edge from two std::vectors with positions and radii. - void - addEdge(std::vector< stim::vec > pos, std::vector > radii, int endId) - { - - edge a(pos,radii, endId); - E.push_back(a); - } */ - - ///A complicated method that adds an edge to the network. - ///Most of this functionality will be moved into fiber. - void - addEdge(std::vector< stim::vec > pos, std::vector > radii, int startId, int endId) - { - // - if(startId == -1 && endId == -1) - { - //case where the edge is neither starting nor ending in a fiber. - //i. first fiber. - - //Add two nodes to the nodes vector - V.push_back(node(pos[pos.size()-1])); - V.push_back(node(pos[0])); - - //the edge will be connected to the nodes - edge *a = new edge(pos,radii,(V.size()-2), (V.size()-1)); - - //add fiber to fiber list. - E.push_back(a); - - //The last two nodes added to V will "own" the last edge added to E. - V[V.size()-1].addEdge(E.size()-1); - V[V.size()-2].addEdge(E.size()-1); - - } - - else if(startId != -1 && endId == -1) - { - //case where the edge is starting with a fiber, but not ending in one. - - //split the fiber behind you into two. - unsigned int point = E[startId]->nearest_idx(pos[0]); - - //split the hit fiber at point two parts temp[0], temp[1] - std::vector < stim::fiber > temp = E[startId]->split(point); - if(temp.size() > 1) - { - //add the nearest point in the behind fiber into the hitting fiber. - pos.insert(pos.begin(), E[startId]->nearest(pos[0])); - stim::vec v(E[startId]->radius(point), E[startId]->radius(point)); - radii.insert(radii.begin(), v); - - //reset the fiber at the endId to be a shorter fiber(temp[0]). - std::vector > ce = temp[0].centerline(); - std::vector > cm = temp[0].centerlinemag(); - - //remake the edge, such that the starting point of this edge - //is the same the split point, but the end edge is the old end. - V.push_back(node(ce[ce.size()-1])); - int tempNodeId = E[startId]->Nodes[1]; - E[startId] = new edge(ce, cm, (V.size()-1), E[startId]->Nodes[2]); - V[V.size()-1].addEdge(startId); - - - //add the other part of the fiber (temp[1]) - ce = temp[1].centerline(); - cm = temp[1].centerlinemag(); - E.push_back(new edge(ce, cm,tempNodeId ,(V.size()-1))); - V[V.size()-1].addEdge(E.size()-1); - - V[tempNodeId].removeEdge(startId); - V[tempNodeId].addEdge(E.size()-1); -// V[V.size()-1].removeEdge(startId); - - //make the new hitting fiber.. - V.push_back(node(pos[pos.size()-1])); - edge *a = new edge(pos, radii, (V.size()-2), (V.size()-1)); - E.push_back(a); - V[V.size()-1].addEdge(E.size()-1); - V[V.size()-2].addEdge(E.size()-1); - } else { - stim::vec pz = E[startId]->nearest(pos[0]); - if((V[E[startId]->Nodes[1]].getPosition() - pz).len() < - (V[E[startId]->Nodes[2]].getPosition() - pz).len()) - { - unsigned int point = E[startId]->nearest_idx(pos[0]); - pos.insert(pos.begin(), E[startId]->nearest(pos[0])); - stim::vec v(E[startId]->radius(point), E[startId]->radius(point)); - radii.insert(radii.begin(), v); - V.push_back(node(pos[pos.size()-1])); - edge *a = new edge(pos, radii, E[startId]->Nodes[1], (V.size()-1)); - E.push_back(a); - - V[V.size()-1].addEdge(E.size()-1); - V[ E[startId]->Nodes[1]].addEdge(E.size()-1); - - } - else - { - unsigned int point = E[startId]->nearest_idx(pos[0]); - pos.insert(pos.begin(), E[startId]->nearest(pos[0])); - stim::vec v(E[startId]->radius(point), E[startId]->radius(point)); - radii.insert(radii.begin(), v); - V.push_back(node(pos[pos.size()-1])); - edge *a = new edge(pos, radii, E[startId]->Nodes[2], (V.size()-1)); - E.push_back(a); - - V[V.size()-1].addEdge(E.size()-1); - V[ E[startId]->Nodes[2]].addEdge(E.size()-1); - } - } - } - - //case where the edge is starting at a seedpoint but ends in a fiber. - if(startId == -1 && endId != -1 && endId < sizeE()) - { - //split the hit fiber into two. - unsigned int point = E[endId]->nearest_idx(pos[pos.size() -1]); - - //split the hit fiber at point into two parts temp[0], temp[1] - std::vector < stim::fiber > temp - = E[endId]->split(point); - if(temp.size() > 1) - { - //add the nearest point in the hit fiber into the hitting fiber. - pos.push_back(E[endId]->nearest(pos[pos.size()-1])); - // pos.insert(pos.end(), E[endId].nearest(pos[pos.size()-1])); - stim::vec v(E[endId]->radius(point), E[endId]->radius(point)); - radii.push_back(v); - - //split the hit fiber at point into two parts temp[0], temp[1] - std::vector < stim::fiber > temp - = E[endId]->split(point); - - //reset the fiber at endId to be a shorter fiber (temp[0]). - std::vector > ce = temp[0].centerline(); - std::vector > cm = temp[0].centerlinemag(); - - //remake the edge, such that the ending point of this edge - //is the same as before, but the starting edge is new. - V.push_back(node(ce[ce.size()-1])); //split point. - int tempNodeId = E[endId]->Nodes[2]; - E[endId] = new edge(ce, cm, E[endId]->Nodes[1], (V.size()-1)); - V[V.size()-1].addEdge(endId); - - //add that other part of the fiber (temp[1]) - //such that it begins with the middle node, and ends with - //the last node of the previous fiber. - ce = temp[1].centerline(); - cm = temp[1].centerlinemag(); - E.push_back(new edge(ce, cm, (V.size()-1), tempNodeId)); - V[V.size()-1].addEdge(E.size()-1); - // V[V.size()-1].removeEdge(endId); - - - - //make the new hitting fiber.. - V.push_back(pos[0]); - edge *a = new edge(pos,radii,(V.size()-1), (V.size()-2)); - E.push_back(a); - V[V.size()-1].addEdge(E.size()-1); - V[V.size()-2].addEdge(E.size()-1); - - //in the end we have added two new nodes and two new edges. - } - else { - stim::vec pz = E[endId]->nearest(pos[0]); - if((V[ E[endId]->Nodes[1]].getPosition() - pz).len() < - (V[E[endId]->Nodes[2]].getPosition() - pz).len()) - { - unsigned int point = E[endId]->nearest_idx(pos[0]); - pos.insert(pos.begin(), E[endId]->nearest(pos[0])); - stim::vec v(E[endId]->radius(point), E[endId]->radius(point)); - radii.insert(radii.begin(), v); - V.push_back(node(pos[pos.size()-1])); - edge *a = new edge(pos, radii, E[endId]->Nodes[1], (V.size()-1)); - E.push_back(a); - - V[V.size()-1].addEdge(E.size()-1); - V[ E[endId]->Nodes[1]].addEdge(E.size()-1); - - } - else - { - unsigned int point = E[endId]->nearest_idx(pos[0]); - pos.insert(pos.begin(), E[endId]->nearest(pos[0])); - stim::vec v(E[endId]->radius(point), E[endId]->radius(point)); - radii.insert(radii.begin(), v); - V.push_back(node(pos[pos.size()-1])); - edge *a = new edge(pos, radii, E[endId]->Nodes[2], (V.size()-1)); - E.push_back(a); - - V[V.size()-1].addEdge(E.size()-1); - V[ E[endId]->Nodes[2]].addEdge(E.size()-1); - } - } - } - - if(startId != -1 && endId != -1 && endId < sizeE()) - { - //case where the edge is starting with a fiber, and ends in one. - - //split the fiber behind you into two. - unsigned int point = E[startId]->nearest_idx(pos[0]); -// std::cout << "in merge to both" << std::endl; - - //split the hit fiber at point two parts temp[0], temp[1] - std::vector < stim::fiber > temp = E[startId]->split(point); - if(temp.size() > 1) - { - //extend the previous fiber (guaranteed to be added last) by one - //and add the - pos = E[E.size()-1]->centerline(); - radii = E[E.size()-1]->centerlinemag(); - pos.insert(pos.begin(), E[startId]->nearest(pos[0])); - stim::vec v(E[startId]->radius(point), E[startId]->radius(point)); - radii.insert(radii.begin(), v); - V.erase(V.end()); - V.push_back(node(pos[0])); - - - //something weird is going on here. - E[E.size()-1] = new edge(pos, radii, (V.size()-2), (V.size()-1)); - V[V.size()-1].addEdge(E.size()-1); - - //reset the fiber at the endId to be a shorter fiber(temp[0]). - std::vector > ce = temp[0].centerline(); - std::vector > cm = temp[0].centerlinemag(); - -// std::cout << ce.size() << std::endl; - //remake the edge, such that the starting point of this edge - //is the same as before, but the end edge is new. - int tempNodeId = E[startId]->Nodes[1]; - E[startId] = new edge(ce, cm, (V.size()-1), E[startId]->Nodes[2]); - V[V.size()-1].addEdge(startId); - - //add the other part of the fiber (temp[1]) - ce = temp[1].centerline(); - cm = temp[1].centerlinemag(); - E.push_back(new edge(ce, cm,tempNodeId, (V.size()-1))); - V[V.size()-1].addEdge(E.size()-1); - V[tempNodeId].removeEdge(startId); - V[tempNodeId].addEdge(E.size()-1); -// V[V.size()-1].removeEdge(startId); - } - else { - stim::vec pz = E[endId]->nearest(pos[0]); - if((V[ E[endId]->Nodes[1]].getPosition() - pz).len() < - (V[E[endId]->Nodes[2]].getPosition() - pz).len()) - { - unsigned int point = E[endId]->nearest_idx(pos[0]); - pos.insert(pos.begin(), E[endId]->nearest(pos[0])); - stim::vec v(E[endId]->radius(point), E[endId]->radius(point)); - radii.insert(radii.begin(), v); - V.push_back(node(pos[pos.size()-1])); - edge *a = new edge(pos, radii, E[endId]->Nodes[1], (V.size()-1)); - E.push_back(a); - - V[V.size()-1].addEdge(E.size()-1); - V[ E[endId]->Nodes[1]].addEdge(E.size()-1); - - } - else - { - unsigned int point = E[endId]->nearest_idx(pos[0]); - pos.insert(pos.begin(), E[endId]->nearest(pos[0])); - stim::vec v(E[endId]->radius(point), E[endId]->radius(point)); - radii.insert(radii.begin(), v); - V.push_back(node(pos[pos.size()-1])); - edge *a = new edge(pos, radii, E[endId]->Nodes[2], (V.size()-1)); - E.push_back(a); - - V[V.size()-1].addEdge(E.size()-1); - V[ E[endId]->Nodes[2]].addEdge(E.size()-1); - } - } - } - - } - - ///@param pos: std::vector of stim vecs with the positions of the point in the fiber. - ///@param radii: std::vector of floats with the radii of the fiber at positions in pos. - ///returns the forest as a std::string. For testing only. - std::string - edges_to_str() - { - std::stringstream ss; - for(unsigned int i = 0; i < E.size(); i++) - { - ss << i << ": " << E[i]->str() << std::endl; - } - return(ss.str()); - } - // total number of points on all edges!=fibers in the network - int - totalEdges() - { - int totalEdgeVertices=0;int N=0; - for (unsigned int i=0; i < sizeE(); i ++) - {// FOR N points on the fiber N-1 edges are possible - N = E[i]->n_pts(); - totalEdgeVertices = totalEdgeVertices + N- 1; - } - return totalEdgeVertices; - } - // sum of all the fiber lengths - float - lengthOfNetwork() - { - float networkLength=0;float N=0; - for (unsigned int i=0; i < sizeE(); i ++) - { - N = E[i]->length(); - networkLength = networkLength + N; - } - return networkLength; - } - ///@param i: index of the required fiber. - ///Returns an std::vector with the centerline of the ith fiber in the edgelist. - std::vector< stim::vec > - getEdgeCenterLine(int i) - { - std::vector < stim::vec > a; - return E[i]->centerline(); - } - - ///@param i: index of the required fiber. - ///Returns an std::vector with the centerline radii of the ith fiber in the edgelist. - std::vector< stim::vec > - getEdgeCenterLineMag(int i) - { - std::vector < stim::vec > a; - return E[i]->centerlinemag(); - } - - ///@param i: index of the required fibers nodes.. - ///Returns an std::string with the start and end points of this node.. - std::string - nodes_to_str(int i) - { - std::stringstream ss; - ss << "[from Node " << E[i] -> Nodes[1] << " to " << E[i] -> Nodes[2] << "]"; - return ss.str(); - } - //load an obj file into a network class - stim::network - objToNetwork(stim::obj objInput) - { - stim::network nwc; - //network network2; - // function to convert an obj file loaded using stim/visualization/obj.h - // to a 3D network class using network.h methods. - std::vector< stim::vec > fiberPositions; //initialize positions on the fiber to be added to network - objInput.getLine(1, fiberPositions); int numDim = fiberPositions[0].size();//find dimensions of the vertices in obj file - // to verify if the nodes are already pushed into node list - std::vector validList; - validList.assign(objInput.numV(), false); - // go through each of the lines "l followed by indices in obj and add all start and end vertices as the nodes - // using addNode function for adding nodes - // and the line as an edge(belongs to fiber class) using addEdge function - std::vector vertexIndices(objInput.numV()); - std::vector< stim::vec > vertices = objInput.getAllV(vertexIndices); - nwc.addVertices(vertices); - for (unsigned i =1; i< objInput.numL() + 1; i++) - { - // edges/fibers could be added to the network class - std::vector< stim::vec > fiberPositions; - - objInput.getLine(i, fiberPositions); - // finding size to allocate radii - int numPointsOnFiber = fiberPositions.size(); - // to extend it to a 3D postion if it is a 1D/2D vertex in space - std::vector< stim::vec > fiberPositions1(numPointsOnFiber); - // 2D case append and assign the last dimension to zero - if (numDim == 2) - {// 2D case append and assign the last dimension to zero repeating it for all the points on fiber - for (int j = 0; j < numPointsOnFiber; j ++) - { - fiberPositions1[j][numDim - 2] = fiberPositions[j][0]; - fiberPositions1[j][numDim -1 ] = fiberPositions[j][1]; - fiberPositions1[j][numDim] = 0; - } - } - // 1D case append and assign last two dimensions to zero - else if (numDim == 1) - { - for (int j = 0; j < numPointsOnFiber; j ++) - { - fiberPositions1[j][numDim - 2] = fiberPositions[j][0]; - fiberPositions1[j][numDim -1 ] = 0; - fiberPositions1[j][numDim] = 0; - } - } - else - { - fiberPositions1 = fiberPositions; - } - // then add edge - //edge* a = new edge(fiberPositions1,radii); - //std::cout<<"here"< > newPointList; - newPointList = Resample(fiberPositions1); - int numPointsOnnewFiber = newPointList.size(); - nwc.addVerticesAfterSamplimg(newPointList); - std::vector radii(numPointsOnnewFiber); // allocating radii to zero - nwc.addEdge(newPointList, radii); - // add nodes - stim::vec v0(3);stim::vec vN(3); - int endId = numPointsOnnewFiber -1; - v0[0] = newPointList[0][0];v0[1] = newPointList[0][1];v0[2] = newPointList[0][2]; - vN[0] = newPointList[endId][0];vN[1] = newPointList[endId][1];vN[2] = newPointList[endId][2]; - // VISITED INDEXES OF the nodes are set to true - if(!validList[objInput.getIndex(vertices, v0)]) - { - //V.push_back(node(v0)); - nwc.addNode(v0); - validList[objInput.getIndex(vertices, v0)] = true; - } - if(!validList[objInput.getIndex(vertices, vN)]) - { - //V.push_back(node(vN)); - nwc.addNode(vN); - validList[objInput.getIndex(vertices, vN)] = true; - } - } - return nwc; - } - // convert ground and test case to network ,kd trees - boost::tuple< ANNkd_tree*, ANNkd_tree*, stim::network, stim::network > - LoadNetworks(std::string gold_filename, std::string test_filename) - { - ANNkd_tree* kdTree1;ANNkd_tree* kdTree2; - using namespace stim; - network network1;network network2; - stim::obj objFile1(gold_filename); - network1 = objToNetwork(objFile1); - kdTree1 = generateKdTreeFromNetwork(network1); - std::cout<<"Gold Standard:network and kdtree generated"< objFile2(test_filename); - network2 = objToNetwork(objFile2); - kdTree2 = generateKdTreeFromNetwork(network2); - std::cout<<"Test Network:network and kdtree generated"< p0, std::vector p1) - { - - double sum = 0; - for(unsigned int d = 0; d < 3; d++) - sum += p0[d] * p1[d]; - return sqrt(sum); - } - // sum of elements in a vector - double sum(stim::vec metricList) - { - float sumMetricList = 0; - for (unsigned int count=0; count network) //kd-tree stores all points in the fiber for fast searching - { - std::cout<<"kd trees generated"< node; - for(unsigned int i = 0; i < n_data; i++) - { - //node = network.V[i].getPosition(); - //convert_to_double - for(unsigned int d = 0; d < 3; d++){ //for each dimension - c[i][d] = double(network.allVerticesAfterSampling[i][d]); //copy the coordinate - } - r[i] = r[i]; //copy the radius - } - ANNpointArray pts = (ANNpointArray)c; //create an array of data points - kdt = new ANNkd_tree(pts, n_data, 3); //build a KD tree - return kdt; - } - - ///@param pos: std::vector of stim vecs with the positions of the point in the fiber. - ///@param radii: std::vector of floats with the radii of the fiber at positions in pos. - ///adds an edge from two std::vectors with positions and radii. - void - addEdge(std::vector< stim::vec > pos, std::vector radii) - { - edge *a = new edge(pos,radii); - E.push_back(a); - } - void - addNode(stim::vec nodes) - { - V.push_back(nodes); - } - void - addVertices(std::vector< stim::vec > vertices) - { - for (unsigned int i=0; i < vertices.size(); i ++) - { - allVerticesList.push_back(vertices[i]); - } - } - void - addVerticesAfterSamplimg(std::vector< stim::vec > vertices) - { - for (unsigned int i=0; i < vertices.size(); i ++) - { - allVerticesAfterSampling.push_back(vertices[i]); - } - } - // gaussian function - float gaussianFunction(float x, float std=25) - { - float evaluate = 1.0f - ((exp(-x/(2*std*std)))); - return evaluate; - } - // compares point on a network to a kd tree for two skeletons - boost::tuple< float, float > - compareSkeletons(boost::tuple< ANNkd_tree*, ANNkd_tree*, stim::network, stim::network > networkKdtree) - { - float gFPR, gFNR; - gFPR = CompareNetKD(networkKdtree.get<0>(), networkKdtree.get<3>()); - gFNR = CompareNetKD(networkKdtree.get<1>(), networkKdtree.get<2>()); - return boost::make_tuple(gFPR, gFNR); - } - // gaussian distance of points on network to Kdtree - float - CompareNetKD(ANNkd_tree* kdTreeGroundtruth, stim::network networkTruth) - { - std::vector< stim::vec > fiberPoints; - float netmetsMetric=0; - double eps = 0; // error bound - ANNdistArray dists1;dists1 = new ANNdist[1]; // near neighbor distances - ANNdistArray dists2; dists2 = new ANNdist[1];// near neighbor distances - ANNidxArray nnIdx1; nnIdx1 = new ANNidx[1]; // near neighbor indices // allocate near neigh indices - ANNidxArray nnIdx2; nnIdx2 = new ANNidx[1]; // near neighbor indices // allocate near neigh indices - int N; int numQueryPoints = networkTruth.totalEdges(); - float totalNetworkLength = networkTruth.lengthOfNetwork(); - stim::vec fiberMetric(networkTruth.sizeE()); - //for each fiber - for (unsigned int i=0; i < networkTruth.sizeE(); i ++) - { - std::vector p1(3); std::vector p2(3);//temporary variables to store point positions - fiberPoints = networkTruth.E[i]->centerline(); - N = networkTruth.E[i]->n_pts(); - stim::vec segmentMetric(N-1); - // for each segment in the fiber - for (unsigned int j = 0; j < N - 1; j++) - { - ANNpoint queryPt1; queryPt1 = annAllocPt(3); - ANNpoint queryPt2; queryPt2 = annAllocPt(3); - //for each dimension of the point - for(unsigned int d = 0; d < 3; d++) - { - queryPt1[d] = double(fiberPoints[j][d]); - p1[d] = double(fiberPoints[j][d]); - queryPt2[d] = double(fiberPoints[j + 1][d]); - p2[d] = double(fiberPoints[j + 1][d]);// for each dimension - } - kdTreeGroundtruth->annkSearch( queryPt1, 1, nnIdx1, dists1, eps); // error bound - kdTreeGroundtruth->annkSearch( queryPt2, 1, nnIdx2, dists2, eps); // error bound - std::cout< > Resample(std::vector< stim::vec > fiberPositions, float spacing=25.0) - { - std::vector p1(3), p2(3), v(3); - stim::vec p(3); - std::vector > newPointList; - for(unsigned int f=0; f= spacing ) - { for (int dim=0; dim<3;dim++) //find the direction of travel - {v[dim] = p2[dim] - p1[dim];} - //set the step size to the voxel size - T step; - for(step=0.0; step v) - { - T sum=0; - for (int i=0; istrObj(strArray, num); - //removeCharsFromString(str,"0"); - ofs << "l " << str << "\n"; - } - ofs.close(); - } - ///exports the graph. - void - to_csv() - { - std::ofstream ofs; - ofs.open("Graph.csv", std::ofstream::out | std::ofstream::app); - std::cout<<"here"<name VARCHAR\n"; - for(int i = 0; i < V.size(); i++) - { - ofs << i << "\n"; - } - ofs << "edgedef>Nodes[1] VARCHAR, Nodes[2] VARCHAR, weight INT, length FLOAT, av_radius FLOAT \n"; - for(int i = 0; i < E.size(); i++) - { - ofs << E[i]->Nodes[1] << "," << E[i]->Nodes[2] << "," <n_pts() - << ","<< E[i]->length() << "," << E[i]->average_radius() << "\n"; - } - ofs.close(); - } - -}; -}; -#endif -- libgit2 0.21.4