From 2e0f052a4108ad47246a9b2be824b8f87fc7e1c5 Mon Sep 17 00:00:00 2001 From: pgovyadi Date: Thu, 4 Aug 2016 16:00:42 -0500 Subject: [PATCH] minor changes to network --- stim/biomodels/network.h | 136 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------------------- stim/math/vec3.h | 1 + stim/visualization/gl_aaboundingbox.h | 5 ++++- 3 files changed, 82 insertions(+), 60 deletions(-) diff --git a/stim/biomodels/network.h b/stim/biomodels/network.h index 8cc699e..59db491 100644 --- a/stim/biomodels/network.h +++ b/stim/biomodels/network.h @@ -33,7 +33,10 @@ class network{ public: unsigned v[2]; //unique id's designating the starting and ending // default constructor - edge() : cylinder(){v[1] = -1; v[0] = -1;} + edge() : cylinder() + { + v[1] = -1; v[0] = -1; + } /// Constructor - creates an edge from a list of points by calling the stim::fiber constructor ///@param p is an array of positions in space @@ -54,7 +57,7 @@ class network{ /// Output the edge information as a string std::string str(){ std::stringstream ss; - ss<<"("<::size()<<")\tl = "<::size()<<")\tl = "< { 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::vec3 p; //position of this node in physical space. + //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::vec3 p; //position of this node in physical space. //constructor takes a stim::vec vertex(stim::vec3 p) : stim::vec3(p){} @@ -99,6 +102,18 @@ protected: public: + ///default constructor + network() + { + + } + + ///constructor with a file to load. + network(std::string fileLocation) + { + load_obj(fileLocation); + } + ///Returns the number of edges in the network. unsigned int edges(){ return E.size(); @@ -110,71 +125,74 @@ public: } stim::cylinder get_cylinder(unsigned f){ - return E[f]; //return the specified edge (casting it to a fiber) + return E[f]; //return the specified edge (casting it to a fiber) } //load a network from an OBJ file void load_obj(std::string filename){ - stim::obj O; //create an OBJ object - O.load(filename); //load the OBJ file as an object + stim::obj O; //create an OBJ object + O.load(filename); //load the OBJ file as an object - std::vector id2vert; //this list stores the OBJ vertex ID associated with each network vertex + std::vector id2vert; //this list stores the OBJ vertex ID associated with each network vertex - unsigned i[2]; //temporary, IDs associated with the first and last points in an OBJ line + unsigned i[2]; //temporary, IDs associated with the first and last points in an OBJ line //for each line in the OBJ object for(unsigned int l = 1; l <= O.numL(); l++){ - std::vector< stim::vec > c; //allocate an array of points for the vessel centerline + std::vector< stim::vec > c; //allocate an array of points for the vessel centerline O.getLine(l, c); //get the fiber centerline std::vector< stim::vec3 > c3(c.size()); for(size_t j = 0; j < c.size(); j++) c3[j] = c[j]; - edge new_edge = c3; //create an edge from the given centerline - unsigned int I = new_edge.size(); //calculate the number of points on the centerline + // edge new_edge = c3; ///This is dangerous. + edge new_edge(c3); + + //create an edge from the given centerline + unsigned int I = new_edge.size(); //calculate the number of points on the centerline //get the first and last vertex IDs for the line - std::vector< unsigned > id; //create an array to store the centerline point IDs + std::vector< unsigned > id; //create an array to store the centerline point IDs O.getLinei(l, id); //get the list of point IDs for the line i[0] = id.front(); //get the OBJ ID for the first element of the line i[1] = id.back(); //get the OBJ ID for the last element of the line - std::vector::iterator it; //create an iterator for searching the id2vert array + std::vector::iterator it; //create an iterator for searching the id2vert array unsigned it_idx; //create an integer for the id2vert entry index //find out if the nodes for this fiber have already been created - it = find(id2vert.begin(), id2vert.end(), i[0]); //look for the first node + it = find(id2vert.begin(), id2vert.end(), i[0]); //look for the first node it_idx = std::distance(id2vert.begin(), it); - if(it == id2vert.end()){ //if i[0] hasn't already been used + if(it == id2vert.end()){ //if i[0] hasn't already been used vertex new_vertex = new_edge[0]; //create a new vertex, assign it a position - new_vertex.e[0].push_back(E.size()); //add the current edge as outgoing - new_edge.v[0] = V.size(); //add the new vertex to the edge - V.push_back(new_vertex); //add the new vertex to the vertex list - id2vert.push_back(i[0]); //add the ID to the ID->vertex conversion list + new_vertex.e[0].push_back(E.size()); //add the current edge as outgoing + new_edge.v[0] = V.size(); //add the new edge to the edge + V.push_back(new_vertex); //add the new vertex to the vertex list + id2vert.push_back(i[0]); //add the ID to the ID->vertex conversion list } - else{ //if the vertex already exists + else{ //if the vertex already exists V[it_idx].e[0].push_back(E.size()); //add the current edge as outgoing new_edge.v[0] = it_idx; } - it = find(id2vert.begin(), id2vert.end(), i[1]); //look for the second ID + it = find(id2vert.begin(), id2vert.end(), i[1]); //look for the second ID it_idx = std::distance(id2vert.begin(), it); - if(it == id2vert.end()){ //if i[1] hasn't already been used - vertex new_vertex = new_edge[I-1]; //create a new vertex, assign it a position - new_vertex.e[1].push_back(E.size()); //add the current edge as incoming + if(it == id2vert.end()){ //if i[1] hasn't already been used + vertex new_vertex = new_edge[I-1]; //create a new vertex, assign it a position + new_vertex.e[1].push_back(E.size()); //add the current edge as incoming new_edge.v[1] = V.size(); - V.push_back(new_vertex); //add the new vertex to the vertex list - id2vert.push_back(i[1]); //add the ID to the ID->vertex conversion list + V.push_back(new_vertex); //add the new vertex to the vertex list + id2vert.push_back(i[1]); //add the ID to the ID->vertex conversion list } - else{ //if the vertex already exists + else{ //if the vertex already exists V[it_idx].e[1].push_back(E.size()); //add the current edge as incoming new_edge.v[1] = it_idx; } - E.push_back(new_edge); //push the edge to the list + E.push_back(new_edge); //push the edge to the list } } @@ -199,17 +217,17 @@ public: /// This function resamples all fibers in a network given a desired minimum spacing /// @param spacing is the minimum distance between two points on the network stim::network resample(T spacing){ - stim::network n; //create a new network that will be an exact copy, with resampled fibers - n.V = V; //copy all vertices + stim::network n; //create a new network that will be an exact copy, with resampled fibers + n.V = V; //copy all vertices - n.E.resize(edges()); //allocate space for the edge list + n.E.resize(edges()); //allocate space for the edge list //copy all fibers, resampling them in the process - for(unsigned e = 0; e < edges(); e++){ //for each edge in the edge list - n.E[e] = E[e].resample(spacing); //resample the edge and copy it to the new network + for(unsigned e = 0; e < edges(); e++){ //for each edge in the edge list + n.E[e] = E[e].resample(spacing); //resample the edge and copy it to the new network } - return n; //return the resampled network + return n; //return the resampled network } @@ -236,14 +254,14 @@ public: /// @param m is the magnitude value to use. The default is 0 (usually radius). T average(unsigned m = 0){ - T M, L; //allocate space for the total magnitude and length - M = L = 0; //initialize both the initial magnitude and length to zero - for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network + T M, L; //allocate space for the total magnitude and length + M = L = 0; //initialize both the initial magnitude and length to zero + for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network M += E[e].integrate(m); //get the integrated magnitude - L += E[e].length(); //get the edge length + L += E[e].length(); //get the edge length } - return M / L; //return the average magnitude + return M / L; //return the average magnitude } /// This function compares two networks and returns the percentage of the current network that is missing from A. @@ -256,17 +274,17 @@ public: R = (*this); //initialize the result with the current network //generate a KD-tree for network A - float metric = 0.0; // initialize metric to be returned after comparing the networks - ANNkd_tree* kdt; // initialize a pointer to a kd tree - double **c; // centerline (array of double pointers) - points on kdtree must be double - unsigned int n_data = A.total_points(); // set the number of points - c = (double**) malloc(sizeof(double*) * n_data); // allocate the array pointer - for(unsigned int i = 0; i < n_data; i++) // allocate space for each point of 3 dimensions + float metric = 0.0; // initialize metric to be returned after comparing the networks + ANNkd_tree* kdt; // initialize a pointer to a kd tree + double **c; // centerline (array of double pointers) - points on kdtree must be double + unsigned int n_data = A.total_points(); // set the number of points + c = (double**) malloc(sizeof(double*) * n_data); // allocate the array pointer + for(unsigned int i = 0; i < n_data; i++) // allocate space for each point of 3 dimensions c[i] = (double*) malloc(sizeof(double) * 3); unsigned t = 0; for(unsigned e = 0; e < A.E.size(); e++){ //for each edge in the network - for(unsigned p = 0; p < A.E[e].size(); p++){ //for each point in the edge + for(unsigned p = 0; p < A.E[e].size(); p++){ //for each point in the edge for(unsigned d = 0; d < 3; d++){ //for each coordinate c[t][d] = A.E[e][p][d]; @@ -276,27 +294,27 @@ public: } //compare each point in the current network to the field produced by A - ANNpointArray pts = (ANNpointArray)c; // create an array of data points of type double - kdt = new ANNkd_tree(pts, n_data, 3); // build a KD tree using the annpointarray + ANNpointArray pts = (ANNpointArray)c; // create an array of data points of type double + kdt = new ANNkd_tree(pts, n_data, 3); // build a KD tree using the annpointarray double eps = 0; // error bound - ANNdistArray dists = new ANNdist[1]; // near neighbor distances - ANNidxArray nnIdx = new ANNidx[1]; // near neighbor indices // allocate near neigh indices + ANNdistArray dists = new ANNdist[1]; // near neighbor distances + ANNidxArray nnIdx = new ANNidx[1]; // near neighbor indices // allocate near neigh indices stim::vec3 p0, p1; float m1; - float M = 0; //stores the total metric value - float L = 0; //stores the total network length + float M = 0; //stores the total metric value + float L = 0; //stores the total network length ANNpoint queryPt = annAllocPt(3); for(unsigned e = 0; e < R.E.size(); e++){ //for each edge in A - R.E[e].add_mag(0); //add a new magnitude for the metric + R.E[e].add_mag(0); //add a new magnitude for the metric - for(unsigned p = 0; p < R.E[e].size(); p++){ //for each point in the edge + for(unsigned p = 0; p < R.E[e].size(); p++){ //for each point in the edge - p1 = R.E[e][p]; //get the next point in the edge + p1 = R.E[e][p]; //get the next point in the edge stim2ann(queryPt, p1); - kdt->annkSearch( queryPt, 1, nnIdx, dists, eps); //find the distance between A and the current network + kdt->annkSearch( queryPt, 1, nnIdx, dists, eps); //find the distance between A and the current network m1 = 1.0f - gaussianFunction((float)dists[0], sigma); //calculate the metric value based on the distance - R.E[e].set_mag(m1, p, 1); //set the error for the second point in the segment + R.E[e].set_mag(m1, p, 1); //set the error for the second point in the segment } } diff --git a/stim/math/vec3.h b/stim/math/vec3.h index 80f7257..880535e 100644 --- a/stim/math/vec3.h +++ b/stim/math/vec3.h @@ -3,6 +3,7 @@ #include +#include namespace stim{ diff --git a/stim/visualization/gl_aaboundingbox.h b/stim/visualization/gl_aaboundingbox.h index de81ae8..3e3b0b6 100644 --- a/stim/visualization/gl_aaboundingbox.h +++ b/stim/visualization/gl_aaboundingbox.h @@ -11,6 +11,9 @@ class gl_aaboundingbox : public aaboundingbox{ public: + using stim::aaboundingbox::A; + using stim::aaboundingbox::B; + //default constructor gl_aaboundingbox() : stim::aaboundingbox(){} @@ -57,4 +60,4 @@ public: }; //end namespace stim -#endif \ No newline at end of file +#endif -- libgit2 0.21.4