codebase / stimlib

Commit 3800c27f696781bd2349f1121cfa39d57bc56aff

Authored by David Mayerich
1 parent 7c43b7f9

added code to subdivide networks

Showing 2 changed files with 60 additions and 36 deletions   Show diff stats
Inline Side-by-side
stim/biomodels/fiber.h
  Show/Hide comments   View file @3800c27
... ... @@ -5,7 +5,7 @@
5 5 #include <ANN/ANN.h>
6 6  
7 7 namespace stim{
8   -
  8 +
9 9 /** This class stores information about a single fiber represented as a set of geometric points
10 10 * between two branch or end points. This class is used as a fundamental component of the stim::network
11 11 * class to describe an interconnected (often biological) network.
... ... @@ -39,7 +39,7 @@ protected:
39 39  
40 40 for(unsigned int i = 0; i < N; i++) //allocate space for each point
41 41 c[i] = (double*) malloc(sizeof(double) * 3);
42   -
  42 +
43 43 r = (T*) malloc(sizeof(T) * N); //allocate space for the radii
44 44 }
45 45  
... ... @@ -62,7 +62,7 @@ protected:
62 62 gen_kdtree(); //generate the kd tree for the new fiber
63 63 }
64 64  
65   - /// generate a KD tree for points on fiber
  65 + /// generate a KD tree for points on fiber
66 66 void gen_kdtree()
67 67 {
68 68 int n_data = N; //create an array of data points
... ... @@ -123,7 +123,7 @@ public:
123 123 copy(obj);
124 124  
125 125 }
126   -
  126 +
127 127 //temp constructor for graph visualization
128 128 fiber(int n)
129 129 {
... ... @@ -218,7 +218,7 @@ public:
218 218  
219 219 return l; //return the length
220 220 }
221   -
  221 +
222 222 /// Calculates the length and average radius of the fiber
223 223  
224 224 /// @param length is filled with the fiber length
... ... @@ -267,7 +267,7 @@ public:
267 267 T r_sum = 0.;
268 268 for(unsigned int i = 0; i < N; i++)
269 269 {
270   - r_sum = r_sum + r[i];
  270 + r_sum = r_sum + r[i];
271 271 }
272 272 return r_sum/((T) N);
273 273 }
... ... @@ -282,16 +282,16 @@ public:
282 282 T radius(int idx){
283 283 return r[idx];
284 284 }
285   - /// get index of a node on a fiber
  285 + /// get index of a node on a fiber
286 286 // by matching the node on fiber to already set vertices (both strings)
287 287 // used in obj file conversion
288   - int
289   - getIndexes(std::string* input, std::string searched, int sizeV)
  288 + int
  289 + getIndexes(std::string* input, std::string searched, int sizeV)
290 290 {
291 291 int result = 0;
292 292 for (int i = 0; i < sizeV; i++)
293 293 {
294   - if (input[i] == searched)
  294 + if (input[i] == searched)
295 295 {
296 296 result = i + 1;
297 297 }
... ... @@ -453,7 +453,7 @@ public:
453 453  
454 454 /// @param i is the index of the element to be returned as a stim::vec
455 455 stim::vec<T> operator[](unsigned i){
456   - return get_vec(i);
  456 + return get_vec(i);
457 457 }
458 458  
459 459 /// Back method returns the last point in the fiber
... ... @@ -461,14 +461,14 @@ public:
461 461 return get_vec(N-1);
462 462 }
463 463 ////resample a fiber in the network
464   - stim::fiber<T> Resample(T spacing=25.0)
  464 + stim::fiber<T> resample(T spacing=25.0)
465 465 {
466 466  
467 467  
468 468 std::vector<T> v(3); //v-direction vector of the segment
469 469 stim::vec<T> p(3); //- intermediate point to be added
470   - std::vector<T> p1(3); // p1 - starting point of an segment on the fiber,
471   - std::vector<T> p2(3); // p2 - ending point,
  470 + std::vector<T> p1(3); // p1 - starting point of an segment on the fiber,
  471 + std::vector<T> p2(3); // p2 - ending point,
472 472 double sum=0; //distance summation
473 473 std::vector<stim::vec<T> > fiberPositions = centerline();
474 474 std::vector<stim::vec<T> > newPointList; // initialize list of new resampled points on the fiber
... ... @@ -483,7 +483,7 @@ public:
483 483 v[d] = p2[d] - p1[d]; //direction vector
484 484 sum +=v[d] * v[d]; //length of segment-distance between starting and ending point
485 485 }
486   - newPointList.push_back(fiberPositions[f]); //always push the starting point
  486 + //newPointList.push_back(fiberPositions[f]); //always push the starting point
487 487  
488 488 T lengthSegment = sqrt(sum); //find Length of the segment as distance between the starting and ending points of the segment
489 489  
... ...
stim/biomodels/network.h
  Show/Hide comments   View file @3800c27
... ... @@ -37,8 +37,20 @@ class network{
37 37  
38 38 /// Constructor - creates an edge from a list of points by calling the stim::fiber constructor
39 39  
40   - ///@param p is a position in space
  40 + ///@param p is an array of positions in space
41 41 edge(std::vector< stim::vec<T> > p) : fiber<T>(p){}
  42 +
  43 + /// Copy constructor creates an edge from a fiber
  44 + edge(stim::fiber<T> f) : fiber<T>(f) {}
  45 +
  46 + /// Resamples an edge by calling the fiber resampling function
  47 + edge resample(T spacing){
  48 + edge e(fiber<T>::resample()); //call the fiber->edge constructor
  49 + e.v[0] = v[0]; //copy the vertex data
  50 + e.v[1] = v[1];
  51 +
  52 + return e; //return the new edge
  53 + }
42 54  
43 55 /// Output the edge information as a string
44 56 std::string str(){
... ... @@ -98,15 +110,16 @@ class network{
98 110 return V.size();
99 111 }
100 112  
  113 + stim::fiber<T> get_fiber(unsigned f){
  114 + return E[f]; //return the specified edge (casting it to a fiber)
  115 + }
  116 +
101 117 //load a network from an OBJ file
102 118 void load_obj(std::string filename){
103 119  
104 120 stim::obj<T> O; //create an OBJ object
105 121 O.load(filename); //load the OBJ file as an object
106 122  
107   - //prints each line in the obj file - vertex positions and fibers
108   - std::cout<<O.str()<<std::endl;
109   -
110 123 std::vector<unsigned> id2vert; //this list stores the OBJ vertex ID associated with each network vertex
111 124  
112 125 unsigned i[2]; //temporary, IDs associated with the first and last points in an OBJ line
... ... @@ -117,7 +130,7 @@ class network{
117 130 std::vector< stim::vec<T> > c; //allocate an array of points for the vessel centerline
118 131 O.getLine(l, c); //get the fiber centerline
119 132  
120   - edge e = c; //create an edge from the given centerline
  133 + edge new_edge = c; //create an edge from the given centerline
121 134  
122 135 //get the first and last vertex IDs for the line
123 136 std::vector< unsigned > id; //create an array to store the centerline point IDs
... ... @@ -132,36 +145,32 @@ class network{
132 145 it = find(id2vert.begin(), id2vert.end(), i[0]); //look for the first node
133 146 it_idx = std::distance(id2vert.begin(), it);
134 147 if(it == id2vert.end()){ //if i[0] hasn't already been used
135   - vertex v = e[0]; //create a new vertex, assign it a position
136   - v.e[0].push_back(E.size()); //add the current edge as outgoing
137   - e.v[0] = V.size(); //add the new vertex to the edge
138   - V.push_back(v); //add the new vertex to the vertex list
  148 + vertex new_vertex = new_edge[0]; //create a new vertex, assign it a position
  149 + new_vertex.e[0].push_back(E.size()); //add the current edge as outgoing
  150 + new_edge.v[0] = V.size(); //add the new vertex to the edge
  151 + V.push_back(new_vertex); //add the new vertex to the vertex list
139 152 id2vert.push_back(i[0]); //add the ID to the ID->vertex conversion list
140 153 }
141 154 else{ //if the vertex already exists
142 155 V[it_idx].e[0].push_back(E.size()); //add the current edge as outgoing
143   - e.v[0] = it_idx;
  156 + new_edge.v[0] = it_idx;
144 157 }
145 158  
146 159 it = find(id2vert.begin(), id2vert.end(), i[1]); //look for the second ID
147 160 it_idx = std::distance(id2vert.begin(), it);
148 161 if(it == id2vert.end()){ //if i[1] hasn't already been used
149   - vertex v = e.back(); //create a new vertex, assign it a position
150   - v.e[1].push_back(E.size()); //add the current edge as incoming
151   - e.v[1] = V.size();
152   - V.push_back(v); //add the new vertex to the vertex list
  162 + vertex new_vertex = new_edge.back(); //create a new vertex, assign it a position
  163 + new_vertex.e[1].push_back(E.size()); //add the current edge as incoming
  164 + new_edge.v[1] = V.size();
  165 + V.push_back(new_vertex); //add the new vertex to the vertex list
153 166 id2vert.push_back(i[1]); //add the ID to the ID->vertex conversion list
154 167 }
155 168 else{ //if the vertex already exists
156 169 V[it_idx].e[1].push_back(E.size()); //add the current edge as incoming
157   - e.v[1] = it_idx;
  170 + new_edge.v[1] = it_idx;
158 171 }
159   - stim::fiber<T> f = e.Resample(25.0);
160   - std::cout<<"resampled fiber length"<<std::endl;
161   - std::cout<<f.length()<<std::endl;
162   - std::cout<<"resampled fiber--"<<std::endl;
163   - std::cout<<f.str()<<std::endl;
164   - E.push_back(e); //push the edge to the list
  172 +
  173 + E.push_back(new_edge); //push the edge to the list
165 174  
166 175 }
167 176 }
... ... @@ -182,6 +191,21 @@ class network{
182 191  
183 192 return ss.str();
184 193 }
  194 +
  195 + /// This function resamples all fibers in a network given a desired minimum spacing
  196 + stim::network<T> resample(T spacing){
  197 + stim::network<T> n; //create a new network that will be an exact copy, with resampled fibers
  198 + n.V = V; //copy all vertices
  199 +
  200 + n.E.resize(E.size()); //allocate space for the edge list
  201 +
  202 + //copy all fibers, resampling them in the process
  203 + for(unsigned e = 0; e < E.size(); e++){ //for each edge in the edge list
  204 + n.E[e] = E[e].resample(spacing); //resample the edge and copy it to the new network
  205 + }
  206 +
  207 + return n; //return the resampled network
  208 + }
185 209 }; //end stim::network class
186 210 }; //end stim namespace
187 211 #endif
... ...