codebase / stimlib

Commit 8761b649ae2c4f4d0921fff6bff6a1a4aa465508

Authored by Pavel Govyadinov
1 parent 8c4f5d84

moved a lot of the tracing functionality into the spider class. From the outside… 

… just attach spider and run trace()
Showing 4 changed files with 395 additions and 475 deletions   Show diff stats
Inline Side-by-side
stim/biomodels/network.h
  Show/Hide comments   View file @8761b64
1 1 #ifndef STIM_NETWORK_H
2 2 #define STIM_NETWORK_H
3 3  
  4 +#include <stdlib.h>
  5 +#include <assert.h>
  6 +#include <sstream>
  7 +#include <fstream>
  8 +#include <algorithm>
  9 +#include <string.h>
  10 +#include <math.h>
4 11 #include <stim/math/vector.h>
5 12 #include <stim/visualization/obj.h>
6   -#include <list>
  13 +#include <stim/visualization/cylinder.h>
7 14 #include <ANN/ANN.h>
  15 +#include <boost/tuple/tuple.hpp>
  16 +
8 17  
9 18 namespace stim{
  19 +/** This is the a class that interfaces with gl_spider in order to store the currently
  20 + * segmented network. The following data is stored and can be extracted:
  21 + * 1)Network geometry and centerline.
  22 + * 2)Network connectivity (a graph of nodes and edges), reconstructed using ANN library.
  23 +*/
10 24  
11   -/** This class provides an interface for dealing with biological networks.
12   - * It takes the following aspects into account:
13   - * 1) Network geometry and centerlines
14   - * 2) Network connectivity (a graph structure can be extracted)
15   - * 3) Network surface structure (the surface is represented as a triangular mesh and referenced to the centerline)
16   - */
17 25  
18 26 template<typename T>
19 27 class network{
20 28  
21   - //-------------------HELPER CLASSES-----------------------
22   - /// Stores information about a geometric point on the network centerline (including point position and radius)
23   - //template<typename T>
24   - class point : public stim::vec<T>{
25   -
26   - public:
27   - T r;
28   -
29   - point() : stim::vec<T>(){}
30   -
31   - //casting constructor
32   - point(stim::vec<T> rhs) : stim::vec<T>(rhs){}
33   - };
34   -
35   - //template<typename T>
36   - class t_node;
37   - class fiber;
38   -
39   - //create typedefs for the iterators to simplify the network code
40   - typedef typename std::list< fiber >::iterator fiber_i;
41   - typedef typename std::list< t_node >::iterator t_node_i;
42   -
43   - /// Stores information about a single capillary (a length of vessel between two branch or end points)
44   - //template<typename T>
45   - class fiber : public std::list< point >{
46   -
47   - using std::list< point >::begin;
48   - using std::list< point >::end;
49   - using std::list< point >::size;
50   -
51   -
52   - public:
53   - //std::list< point > P; //geometric point positions
54   -
55   - typename std::list< t_node >::iterator n[2]; //indices to terminal nodes
56   - unsigned int id;
57   -
58   - public:
59   -
60   - /// Calculate the length of the fiber and return it.
61   - T length(){
62   -
63   - point p0, p1;
64   - T l = 0; //initialize the length to zero
65   -
66   - //for each point
67   - typename std::list< point >::iterator i; //create a point iterator
68   - for(i = begin(); i != end(); i++){ //for each point in the fiber
69   -
70   - if(i == begin()) //if this is the first point, just store it
71   - p1 = *i;
72   - else{ //if this is any other point
73   - p0 = p1; //shift p1->p0
74   - p1 = *i; //set p1 to the new point
75   - l += (p1 - p0).len(); //add the length of p1 - p0 to the running sum
76   - }
77   - }
78   -
79   - return l; //return the length
80   - }
81   -
82   - T radius(T& length){
83   -
84   - point p0, p1; //temporary variables to store point positions
85   - T r0, r1; //temporary variables to store radii at points
86   - T l, r; //temporary variable to store the length and average radius of a fiber segment
87   - T length_sum = 0; //initialize the length to zero
88   - T radius_sum = 0; //initialize the radius sum to zero
89   -
90   - //for each point
91   - typename std::list< point >::iterator i; //create a point iterator
92   - for(i = begin(); i != end(); i++){ //for each point in the fiber
93   -
94   - if(i == begin()){ //if this is the first point, just store it
95   - p1 = *i;
96   - r1 = i->r;
97   - }
98   - else{ //if this is any other point
99   - p0 = p1; //shift p1->p0 and r1->r0
100   - r0 = r1;
101   - p1 = *i; //set p1 to the new point
102   - r1 = i->r; //and r1
103   -
104   - l = (p1 - p0).len(); //calculate the length of the p0-p1 segment
105   - r = (r0 + r1) / 2; //calculate the average radius of the segment
106   -
107   - radius_sum += r * l; //add the radius scaled by the length to a running sum
108   - length_sum += l; //add the length of p1 - p0 to the running sum
109   - }
110   - }
111   -
112   - length = length_sum; //store the total length
113   -
114   - //if the total length is zero, store a radius of zero
115   - if(length == 0)
116   - return 0;
117   - else
118   - return radius_sum / length; //return the average radius of the fiber
119   - }
120   -
121   - std::vector< stim::vec<T> > geometry(){
122   -
123   - std::vector< stim::vec<T> > result; //create an array to store the fiber geometry
124   - result.resize( size() ); //pre-allocate the array
125   -
126   - typename std::list< point >::iterator p; //create a list iterator
127   - unsigned int pi = 0; //create an index into the result array
128   -
129   - //for each geometric point on the fiber centerline
130   - for(p = begin(); p != end(); p++){
131   - result[pi] = *p;
132   - pi++;
133   - }
134   -
135   - return result; //return the geometry array
136   -
  29 + ///Each edge is a fiber with two nodes.
  30 + ///Each node is an in index to the endpoint of the fiber in the nodes array.
  31 + class edge : public cylinder<T>
  32 + {
  33 + public:
  34 + unsigned v[2]; //unique id's designating the starting and ending
  35 + // default constructor
  36 + edge() : cylinder<T>(){v[1] = -1; v[0] = -1;}
  37 + /// Constructor - creates an edge from a list of points by calling the stim::fiber constructor
  38 +
  39 + ///@param p is an array of positions in space
  40 + edge(std::vector< stim::vec<T> > p) : cylinder<T>(p){}
  41 +
  42 + /// Copy constructor creates an edge from a fiber
  43 + edge(stim::cylinder<T> f) : cylinder<T>(f) {}
  44 +
  45 + /// Resamples an edge by calling the fiber resampling function
  46 + edge resample(T spacing){
  47 + edge e(cylinder<T>::resample(spacing)); //call the fiber->edge constructor
  48 + e.v[0] = v[0]; //copy the vertex data
  49 + e.v[1] = v[1];
  50 +
  51 + return e; //return the new edge
137 52 }
138 53  
  54 + /// Output the edge information as a string
139 55 std::string str(){
140 56 std::stringstream ss;
141   -
142   - //create an iterator for the point list
143   - typename std::list<point>::iterator i;
144   - for(i = begin(); i != end(); i++){
145   - ss<<i->str()<<" r = "<<i->r<<std::endl;
146   - }
147   -
  57 + ss<<"("<<cylinder<T>::size()<<")\tl = "<<length()<<"\t"<<v[0]<<"----"<<v[1];
148 58 return ss.str();
149 59 }
150   - };
151   -
152   - /// Terminal node for a capillary. This is analogous to a graph vertex and contains a list of edge indices.
153   - //template<typename T>
154   - class t_node{
155 60  
156   - public:
157   -
158   - unsigned int id;
159   -
160   - //lists of edge indices for capillaries
161   - //the "in" and "out" just indicate how the geometry is defined:
162   - // edges in the "in" list are geometrically oriented such that the terminal node is last
163   - // edges in the "out" list are geometrically oriented such that the terminal node is first
164   - std::list< fiber_i > in; //edge indices for incoming capillaries
165   - std::list< fiber_i > out; //edge indices for outgoing capillaries
166   -
167   - std::string str(){
168   -
169   - std::stringstream ss;
170   -
171   - ss<<id<<": "; //output the node ID
172   -
173   - //output the IDs for both lists
174   - typename std::list< fiber_i >::iterator f;
175   -
176   - for(f = in.begin(); f != in.end(); f++){
177   -
178   - if(f != in.begin())
179   - ss<<", ";
180   - ss<<(*f)->n[0]->id;
181   - }
182   -
183   - //if there are nodes in both lists, separate them by a comma
184   - if(out.size() > 0 && in.size() > 0)
185   - ss<<", ";
  61 + };
186 62  
187   - for(f = out.begin(); f != out.end(); f++){
  63 + ///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.
  64 + class vertex : public stim::vec<T>
  65 + {
  66 + public:
  67 + //std::vector<unsigned int> edges; //indices of edges connected to this node.
  68 + std::vector<unsigned int> e[2]; //indices of edges going out (e[0]) and coming in (e[1])
  69 + //stim::vec<T> p; //position of this node in physical space.
  70 +
  71 + //constructor takes a stim::vec
  72 + vertex(stim::vec<T> p) : stim::vec<T>(p){}
  73 +
  74 + /// Output the vertex information as a string
  75 + std::string str(){
  76 + std::stringstream ss;
  77 + ss<<"\t(x, y, z) = "<<stim::vec<T>::str();
  78 +
  79 + if(e[0].size() > 0){
  80 + ss<<"\t> ";
  81 + for(unsigned int o = 0; o < e[0].size(); o++)
  82 + ss<<e[0][o]<<" ";
  83 + }
  84 + if(e[1].size() > 0){
  85 + ss<<"\t< ";
  86 + for(unsigned int i = 0; i < e[1].size(); i++)
  87 + ss<<e[1][i]<<" ";
  88 + }
188 89  
189   - if(f != out.begin())
190   - ss<<", ";
191   - ss<<(*f)->n[1]->id;
  90 + return ss.str();
192 91 }
193 92  
194   -
195   - return ss.str();
196   -
197   -
198   -
199   - }
200 93 };
201 94  
202   -
203   -//---------------NETWORK CLASS-----------------------------
204   -
205 95 protected:
206 96  
207   - //list of terminal nodes
208   - std::list<t_node> N;
209   -
210   - //list of fibers
211   - std::list<fiber> F;
212   -
213   - /// Sets a unique ID for each terminal node and fiber
214   - void set_names(){
215   -
216   - unsigned int i;
217   -
218   - i = 0;
219   - for(t_node_i ti = N.begin(); ti != N.end(); ti++)
220   - ti->id = i++;
221   -
222   - i = 0;
223   - for(fiber_i fi = F.begin(); fi != F.end(); fi++)
224   - fi->id = i++;
225   - }
  97 + std::vector<edge> E; //list of edges
  98 + std::vector<vertex> V; //list of vertices.
226 99  
227 100 public:
228 101  
229   - std::string str(){
230   -
231   -
232   - //assign names to elements of the network
233   - set_names();
234   -
235   - //create a stringstream for output
236   - std::stringstream ss;
237   -
238   - //output the nodes
239   - ss<<"Nodes-----------------------------------------------"<<std::endl;
240   - for(t_node_i i = N.begin(); i != N.end(); i++){
241   - ss<<i->str()<<std::endl;
242   - }
243   -
244   - //output the fibers
245   - ss<<std::endl<<"Fibers----------------------------------------------"<<std::endl;
246   -
247   - T length, radius;
248   - //output every fiber
249   - for(fiber_i f = F.begin(); f != F.end(); f++){
250   -
251   - //calculate the length and average radius
252   - radius = f->radius(length);
253   -
254   - //output the IDs of the terminal nodes
255   - ss<<f->n[0]->id<<" -- "<<f->n[1]->id<<": length = "<<length<<", average radius = "<<radius<<std::endl;
256   - }
  102 + ///Returns the number of edges in the network.
  103 + unsigned int edges(){
  104 + return E.size();
  105 + }
257 106  
258   - return ss.str();
  107 + ///Returns the number of nodes in the network.
  108 + unsigned int vertices(){
  109 + return V.size();
259 110 }
260 111  
261   - /// Load a network from an OBJ object
  112 + stim::cylinder<T> get_cylinder(unsigned f){
  113 + return E[f]; //return the specified edge (casting it to a fiber)
  114 + }
262 115  
263   - /// @param object is the object file to be used as the basis for the network
264   - void load( stim::obj<T> object){
  116 + //load a network from an OBJ file
  117 + void load_obj(std::string filename){
265 118  
266   - //get the number of vertices in the object
267   - unsigned int nV = object.numV();
  119 + stim::obj<T> O; //create an OBJ object
  120 + O.load(filename); //load the OBJ file as an object
268 121  
269   - //allocate an array of pointers to nodes, which will be used to preserve connectivity
270   - //initiate all values to T.end()
271   - std::vector< t_node_i > node_hash(nV, N.end());
  122 + std::vector<unsigned> id2vert; //this list stores the OBJ vertex ID associated with each network vertex
272 123  
273   - unsigned int nL = object.numL(); //get the number of lines in the OBJ
  124 + unsigned i[2]; //temporary, IDs associated with the first and last points in an OBJ line
274 125  
275   - //for each line in the OBJ structure
276   - for(unsigned int li = 0; li < nL; li++){
  126 + //for each line in the OBJ object
  127 + for(unsigned int l = 1; l <= O.numL(); l++){
277 128  
278   - F.push_back(fiber()); //push a new fiber onto the fiber list
  129 + std::vector< stim::vec<T> > c; //allocate an array of points for the vessel centerline
  130 + O.getLine(l, c); //get the fiber centerline
279 131  
280   - fiber_i f = --(F.end()); //get an iterator to the new fiber
  132 + edge new_edge = c; //create an edge from the given centerline
  133 + unsigned int I = new_edge.size(); //calculate the number of points on the centerline
281 134  
282   - //----------Handle the terminating nodes for the fiber
  135 + //get the first and last vertex IDs for the line
  136 + std::vector< unsigned > id; //create an array to store the centerline point IDs
  137 + O.getLinei(l, id); //get the list of point IDs for the line
  138 + i[0] = id.front(); //get the OBJ ID for the first element of the line
  139 + i[1] = id.back(); //get the OBJ ID for the last element of the line
283 140  
284   - //get the indices of the line vertices
285   - std::vector< unsigned int > Li = object.getL_Vi(li);
286   - unsigned int i0 = Li.front() - 1;
287   - unsigned int i1 = Li.back() - 1;
  141 + std::vector<unsigned>::iterator it; //create an iterator for searching the id2vert array
  142 + unsigned it_idx; //create an integer for the id2vert entry index
288 143  
289   - //deal with the first end point of the capillary
290   - if(node_hash[i0] != N.end()){ //if the node has been used before
291   - (*f).n[0] = node_hash[i0]; //assign the node to the new capillary
292   - (*node_hash[i0]).out.push_back(f); //add an out pointer to the existing node
  144 + //find out if the nodes for this fiber have already been created
  145 + it = find(id2vert.begin(), id2vert.end(), i[0]); //look for the first node
  146 + it_idx = std::distance(id2vert.begin(), it);
  147 + if(it == id2vert.end()){ //if i[0] hasn't already been used
  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
  152 + id2vert.push_back(i[0]); //add the ID to the ID->vertex conversion list
293 153 }
294   - else{ //otherwise
295   - N.push_back(t_node()); //create a new node and add it to the node list
296   - t_node_i t = --(N.end()); //get an iterator to the new node
297   - node_hash[i0] = t; //add a pointer to the new node to the hash list
298   - (*f).n[0] = t; //add a pointer to the new node to the capillary
299   - (*t).out.push_back(f); //add a pointer to the capillary to the new node
  154 + else{ //if the vertex already exists
  155 + V[it_idx].e[0].push_back(E.size()); //add the current edge as outgoing
  156 + new_edge.v[0] = it_idx;
300 157 }
301 158  
302   - //deal with the last end point of the capillary
303   - if(node_hash[i1] != N.end()){
304   - (*f).n[1] = node_hash[i1];
305   - (*node_hash[i1]).in.push_back(f);
  159 + it = find(id2vert.begin(), id2vert.end(), i[1]); //look for the second ID
  160 + it_idx = std::distance(id2vert.begin(), it);
  161 + if(it == id2vert.end()){ //if i[1] hasn't already been used
  162 + vertex new_vertex = new_edge[I-1]; //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
  166 + id2vert.push_back(i[1]); //add the ID to the ID->vertex conversion list
306 167 }
307   - else{
308   - N.push_back(t_node());
309   - t_node_i t = --(N.end());
310   - node_hash[i1] = t; //add the new node to the hash list
311   - (*f).n[1] = t;
312   - (*t).in.push_back(f);
  168 + else{ //if the vertex already exists
  169 + V[it_idx].e[1].push_back(E.size()); //add the current edge as incoming
  170 + new_edge.v[1] = it_idx;
313 171 }
314 172  
315   - //-------------Handle the geometric points for the fiber
316   - std::vector< vec<T> > L = object.getL_V(li);
317   - std::vector< vec<T> > R = object.getL_VT(li);
  173 + E.push_back(new_edge); //push the edge to the list
318 174  
319   - unsigned int nP = L.size(); //get the number of geometric points in the fiber
320   - //for each vertex in the fiber
321   - for(unsigned int pi = 0; pi < nP; pi++){
322   - point p = (point)L[pi]; //move the geometric coordinates into a point structure
323   - p.r = R[pi][0]; //store the radius
324   - f->push_back(p); //push the point onto the current fiber
325   - }
326 175 }
  176 + }
327 177  
328   - } //end load()
329   -
330   - /// Returns an array of node positions
331   - std::vector< stim::vec<T> > get_node_positions(){
332   -
333   - std::vector< stim::vec<T> > result; //create an array to store the result
334   - result.resize(N.size()); //set the array size
335   -
336   - t_node_i ni; //create a terminal node iterator
337   - unsigned int vi = 0; //vertex index into the result array
338   -
339   - //for every terminal node
340   - for(ni = N.begin(); ni != N.end(); ni++){
341   -
342   - //create a vector based on the node position
343   -
344   - //if the number of outgoing nodes is nonzero
345   - if(ni->out.size() != 0)
346   - result[vi] = ni->out.front()->front();
347   - else if(ni->in.size() != 0)
348   - result[vi] = ni->in.front()->back();
  178 + /// Output the network as a string
  179 + std::string str(){
349 180  
350   - vi++; //increment the array index
  181 + std::stringstream ss;
  182 + ss<<"Nodes ("<<V.size()<<")--------"<<std::endl;
  183 + for(unsigned int v = 0; v < V.size(); v++){
  184 + ss<<"\t"<<v<<V[v].str()<<std::endl;
351 185 }
352 186  
353   - //return the resulting array
354   - return result;
355   - }
  187 + ss<<"Edges ("<<E.size()<<")--------"<<std::endl;
  188 + for(unsigned e = 0; e < E.size(); e++){
  189 + ss<<"\t"<<e<<E[e].str()<<std::endl;
  190 + }
356 191  
357   - std::vector< stim::vec<T> > get_fiber_geometry( fiber_i f ){
358   - return f->geometry();
  192 + return ss.str();
359 193 }
360 194  
361   - /// Generate an OBJ file from the network
362   -
363   - stim::obj<T> obj(){
364   -
365   - //create an OBJ object
366   - stim::obj<T> object;
367   -
368   - //name the nodes
369   - set_names();
  195 + /// This function resamples all fibers in a network given a desired minimum spacing
  196 + /// @param spacing is the minimum distance between two points on the network
  197 + stim::network<T> resample(T spacing){
  198 + stim::network<T> n; //create a new network that will be an exact copy, with resampled fibers
  199 + n.V = V; //copy all vertices
370 200  
371   - //retrieve a list of terminal node positions
372   - std::vector< stim::vec<T> > node_pos = get_node_positions();
  201 + n.E.resize(edges()); //allocate space for the edge list
373 202  
374   - //add the nodes to the obj file
375   - object.addV(node_pos);
376   -
377   - //counter for vertex indices in the object class
378   - unsigned int nP;
379   -
380   - //for each fiber
381   - fiber_i fi; //create a fiber iterator
382   - for(fi = F.begin(); fi != F.end(); fi++){
  203 + //copy all fibers, resampling them in the process
  204 + for(unsigned e = 0; e < edges(); e++){ //for each edge in the edge list
  205 + n.E[e] = E[e].resample(spacing); //resample the edge and copy it to the new network
  206 + }
383 207  
384   - //get an array of fiber points
385   - std::vector< stim::vec<T> > fiber_p = get_fiber_geometry(fi);
  208 + return n; //return the resampled network
  209 + }
386 210  
387   - //create a subset of this array
388   - typename std::vector< stim::vec<T> >::iterator start = fiber_p.begin() + 1;
389   - typename std::vector< stim::vec<T> >::iterator end = fiber_p.end() - 1;
390   - typename std::vector< stim::vec<T> > fiber_subset(start, end);
391 211  
392   - //add this subset to the geometry object
393   - nP = object.addV(fiber_subset);
394 212  
395   - //create an array to hold vertex indices for a line
396   - std::vector<unsigned int> line;
397   - line.resize(fiber_p.size());
  213 + /// Calculate the total number of points on all edges.
  214 + unsigned total_points(){
  215 + unsigned n = 0;
  216 + for(unsigned e = 0; e < E.size(); e++)
  217 + n += E[e].size();
  218 + return n;
  219 + }
398 220  
399   - //add the terminal nodes to the line list (make sure to add 1 to make them compatible with the OBJ)
400   - line[0] = fi->n[0]->id + 1;
401   - line[line.size() - 1] = fi->n[1]->id + 1;
  221 + // gaussian function
  222 + float gaussianFunction(float x, float std=25){ return exp(-x/(2*std*std));} // by default std = 25
402 223  
403   - //add the intermediate vertex indices to the line array
404   - for(unsigned int i = 0; i < fiber_subset.size(); i++){
405   - line[1 + i] = nP + i;
406   - }
  224 + // stim 3d vector to annpoint of 3 dimensions
  225 + void stim2ann(ANNpoint &a, stim::vec<T> b){
  226 + a[0] = b[0];
  227 + a[1] = b[1];
  228 + a[2] = b[2];
  229 + }
407 230  
408   - //add the line list to the object class
409   - object.addLine(line);
  231 + /// Calculate the average magnitude across the entire network.
  232 + /// @param m is the magnitude value to use. The default is 0 (usually radius).
  233 + T average(unsigned m = 0){
410 234  
  235 + T M, L; //allocate space for the total magnitude and length
  236 + M = L = 0; //initialize both the initial magnitude and length to zero
  237 + for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network
  238 + M += E[e].integrate(m); //get the integrated magnitude
  239 + L += E[e].length(); //get the edge length
411 240 }
412 241  
413   - return object;
  242 + return M / L; //return the average magnitude
414 243 }
415 244  
416   - /// This function returns the information necessary for a simple graph-based physical (ex. fluid) simulation.
417   -
418   - /// @param n0 is a array which will contain the list of source nodes
419   - /// @param n1 is a array which will contain the list of destination nodes
420   - /// @param length is a array containing the lengths of fibers in the network
421   - /// @param radius is a array containing the average radii of fibers in the network
422   - void build_simgraph(std::vector<unsigned int>& n0, std::vector<unsigned int>& n1, std::vector<T>& length, std::vector<T>& radius){
  245 + /// This function compares two networks and returns the percentage of the current network that is missing from A.
423 246  
424   - //determine the number of fibers in the network
425   - unsigned int nF = F.size();
  247 + /// @param A is the network to compare to - the field is generated for A
  248 + /// @param sigma is the user-defined tolerance value - smaller values provide a stricter comparison
  249 + stim::network<T> compare(stim::network<T> A, float sigma){
426 250  
427   - //allocate the necessary space to store the fiber information
428   - n0.resize(nF);
429   - n1.resize(nF);
430   - length.resize(nF);
431   - radius.resize(nF);
  251 + stim::network<T> R; //generate a network storing the result of the comparison
  252 + R = (*this); //initialize the result with the current network
432 253  
433   - //assign names (identifiers) to the network components
434   - set_names();
  254 + //generate a KD-tree for network A
  255 + float metric = 0.0; // initialize metric to be returned after comparing the networks
  256 + ANNkd_tree* kdt; // initialize a pointer to a kd tree
  257 + double **c; // centerline (array of double pointers) - points on kdtree must be double
  258 + unsigned int n_data = A.total_points(); // set the number of points
  259 + c = (double**) malloc(sizeof(double*) * n_data); // allocate the array pointer
  260 + for(unsigned int i = 0; i < n_data; i++) // allocate space for each point of 3 dimensions
  261 + c[i] = (double*) malloc(sizeof(double) * 3);
435 262  
436   - //fill the arrays
437   - unsigned int i = 0;
438   - T l, r;
439   - for(fiber_i f = F.begin(); f != F.end(); f++){
440   - n0[i] = f->n[0]->id; //get the identifiers for the first and second nodes for the current fiber
441   - n1[i] = f->n[1]->id;
  263 + unsigned t = 0;
  264 + for(unsigned e = 0; e < A.E.size(); e++){ //for each edge in the network
  265 + for(unsigned p = 0; p < A.E[e].size(); p++){ //for each point in the edge
  266 + for(unsigned d = 0; d < 3; d++){ //for each coordinate
442 267  
443   - r = f->radius(l); //get the length and radius of the capillary (calculated at the same time)
  268 + c[t][d] = A.E[e][p][d];
  269 + }
  270 + t++;
  271 + }
  272 + }
444 273  
445   - radius[i] = r; //store the radius in the output array
446   - length[i] = l; //store the length in the output array
  274 + //compare each point in the current network to the field produced by A
  275 + ANNpointArray pts = (ANNpointArray)c; // create an array of data points of type double
  276 + kdt = new ANNkd_tree(pts, n_data, 3); // build a KD tree using the annpointarray
  277 + double eps = 0; // error bound
  278 + ANNdistArray dists = new ANNdist[1]; // near neighbor distances
  279 + ANNidxArray nnIdx = new ANNidx[1]; // near neighbor indices // allocate near neigh indices
  280 +
  281 + stim::vec<T> p0, p1;
  282 + float m0, m1;
  283 + float M = 0; //stores the total metric value
  284 + float l; //stores the segment length
  285 + float L = 0; //stores the total network length
  286 + ANNpoint queryPt = annAllocPt(3);
  287 + for(unsigned e = 0; e < R.E.size(); e++){ //for each edge in A
  288 + R.E[e].add_mag(0); //add a new magnitude for the metric
  289 +
  290 + for(unsigned p = 0; p < R.E[e].size(); p++){ //for each point in the edge
  291 +
  292 + p1 = R.E[e][p]; //get the next point in the edge
  293 + stim2ann(queryPt, p1);
  294 + kdt->annkSearch( queryPt, 1, nnIdx, dists, eps); //find the distance between A and the current network
  295 + m1 = 1.0f - gaussianFunction(dists[0], sigma); //calculate the metric value based on the distance
  296 + R.E[e].set_mag(m1, p, 1); //set the error for the second point in the segment
447 297  
448   - i++; //increment the array index
  298 + }
449 299 }
450 300  
451   -
  301 + return R; //return the resulting network
452 302 }
453 303  
454   -};
455   -
456   -}; //end namespace stim
  304 + /// Returns the number of magnitude values stored in each edge. This should be uniform across the network.
  305 + unsigned nmags(){
  306 + return E[0].nmags();
  307 + }
457 308  
458 309  
  310 +}; //end stim::network class
  311 +}; //end stim namespace
459 312 #endif
... ...
stim/gl/gl_spider.h
  Show/Hide comments   View file @8761b64
... ... @@ -96,7 +96,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
96 96 std::vector< stim::vec<float> > cD; //Direction of line currently being traced.
97 97 std::vector< stim::vec<float> > cM; //Magnitude of line currently being traced.
98 98  
99   - stim::glObj<float> sk; //object to store the skeleton.
  99 +// stim::glObj<float> sk; //object to store the skeleton.
100 100 stim::glnetwork<float> nt; //object for storing the network.
101 101  
102 102 stim::vec<float> rev; //reverse vector;
... ... @@ -1157,7 +1157,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1157 1157 {
1158 1158 float x, y, z, u, v, w, m;
1159 1159 myfile >> x >> y >> z >> u >> v >> w >> m;
1160   - setSeed(x, y , z);
  1160 + setSeed(x, y, z);
1161 1161 setSeedVec(u, v, w);
1162 1162 setSeedMag(m);
1163 1163 }
... ... @@ -1171,14 +1171,28 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1171 1171 void
1172 1172 saveNetwork(std::string name)
1173 1173 {
  1174 + stim::glObj<float> sk;
  1175 + for(int i = 0; i < nt.sizeE(); i++)
  1176 + {
  1177 + std::vector<stim::vec< float > > cm = nt.getEdgeCenterLineMag(i);
  1178 + std::vector<stim::vec< float > > ce = nt.getEdgeCenterLine(i);
  1179 + sk.Begin(stim::OBJ_LINE);
  1180 + for(int j = 0; j < ce.size(); j++)
  1181 + {
  1182 + sk.TexCoord(cm[j][0]);
  1183 + sk.Vertex(ce[j][0], ce[j][1], ce[j][2]);
  1184 + }
  1185 + sk.End();
  1186 + }
1174 1187 sk.save(name);
1175 1188 }
1176 1189  
  1190 + ///Depreciated, but might be reused later()
1177 1191 ///returns a COPY of the entire stim::glObj object.
1178 1192 stim::glObj<float>
1179 1193 getNetwork()
1180 1194 {
1181   - return sk;
  1195 +// return sk;
1182 1196 }
1183 1197  
1184 1198 ///returns a COPY of the entire stim::glnetwork object.
... ... @@ -1326,114 +1340,35 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1326 1340 void
1327 1341 trace(int min_cost)
1328 1342 {
1329   - Bind();
1330   - rev = stim::vec<float>(0.0,0.0,1.0);
  1343 +// rev = stim::vec<float>(0.0,0.0,1.0);
1331 1344 bool sEmpty = true;
1332 1345 float lastmag = 16.0;;
1333   - while(!seeds.empty())
  1346 + stim::vec<float> curSeed;
  1347 + stim::vec<float> curSeedVec;
  1348 + float curSeedMag;
  1349 + while(!Empty())
1334 1350 {
1335 1351 //clear the currently traced line and start a new one.
1336 1352 cL.clear();
1337 1353 cM.clear();
1338   - sk.Begin(stim::OBJ_LINE);
1339   - stim::vec<float> curSeed = seeds.top();
1340   -// std::cout << "The current seeds is " << curSeed << std::endl;
1341   - stim::vec<float> curSeedVec = seedsvecs.top();
1342   - float curSeedMag = seedsmags.top();
  1354 + cD.clear();
  1355 + curSeed = seeds.top();
  1356 + curSeedVec = seedsvecs.top();
  1357 + curSeedMag = seedsmags.top();
1343 1358 seeds.pop();
1344 1359 seedsvecs.pop();
1345 1360 seedsmags.pop();
1346 1361 // std::cout << "The current seed Vector is " << curSeedVec << std::endl;
1347 1362 setPosition(curSeed);
1348 1363 setDirection(curSeedVec);
1349   - cL.push_back(curSeed);
1350   - cM.push_back(curSeedMag);
1351   - sk.createFromSelf(GL_SELECT);
1352   - traceLine(min_cost);
1353   -
1354   - sk.rev();
1355   - // std::cout << "reversed" << std::endl;
1356   - std::reverse(cL.begin(), cL.end());
1357   - std::reverse(cM.begin(), cM.end());
1358   - setPosition(curSeed);
1359   - setDirection(-rev);
1360   - setMagnitude(16.0);
1361   - sk.createFromSelf(GL_SELECT);
1362   - traceLine(min_cost);
1363   - sk.End();
  1364 + setMagnitude(curSeedMag);
  1365 +// cL.push_back(curSeed);
  1366 +// cM.push_back(curSeedMag);
  1367 +// cD.push_back(curSeedMag);
  1368 + pair<stim::fiber<float>, int> a = traceLine(p, m, min_cost);
1364 1369 }
1365   - Unbind();
1366 1370 }
1367 1371  
1368   - ///@param min_cost the cost value used for tracing
1369   - ///traces the seedpoint passed to completion in one directions.
1370   - void
1371   - traceLine(int min_cost)
1372   - {
1373   - stim::vec<float> pos;
1374   - stim::vec<float> mag;
1375   - int h;
1376   - bool started = false;
1377   - bool running = true;
1378   - stim::vec<float> size(S[0]*R[0], S[1]*R[1], S[2]*R[2]);
1379   - while(running)
1380   - {
1381   - int cost = Step();
1382   - if (cost > min_cost){
1383   - running = false;
1384   - break;
1385   - } else {
1386   - //Have we found an edge?
1387   - pos = getPosition();
1388   - if(pos[0] > size[0] || pos[1] > size[1]
1389   - || pos[2] > size[2] || pos[0] < 0
1390   - || pos[1] < 0 || pos[2] < 0)
1391   - {
1392   -// std::cout << "Found Edge" << std::endl;
1393   - running = false;
1394   - break;
1395   - }
1396   - //If this is the first step in the trace,
1397   - // save the direction
1398   - //(to be used later to trace the fiber in the opposite direction)
1399   - if(started == false){
1400   - rev = -getDirection();
1401   - started = true;
1402   - }
1403   -// std::cout << i << p << std::endl;
1404   - m = getMagnitude();
1405   - //Has the template size gotten unreasonable?
1406   - if(m[0] > 75 || m[0] < 1){
1407   -// std::cout << "Magnitude Limit" << std::endl;
1408   - running = false;
1409   - break;
1410   - }
1411   - else
1412   - {
1413   - h = selectObject(pos, getDirection(), m[0]);
1414   - //Have we hit something previously traced?
1415   - if(h != -1){
1416   - std::cout << "I hit a line" << h << std::endl;
1417   - running = false;
1418   - break;
1419   - }
1420   - else {
1421   - cL.push_back(stim::vec<float>(p[0], p[1],p[2]));//
1422   - sk.TexCoord(m[0]);
1423   - sk.Vertex(p[0], p[1], p[2]);
1424   - Bind(btexbufferID, bfboID, 27);
1425   - CHECK_OPENGL_ERROR
1426   - branchDetection();
1427   - CHECK_OPENGL_ERROR
1428   - Unbind();
1429   - CHECK_OPENGL_ERROR
1430   - }
1431   - }
1432   - }
1433   - }
1434   - }
1435   -
1436   -
1437 1372 int
1438 1373 selectObject(stim::vec<float> loc, stim::vec<float> dir, float mag)
1439 1374 {
... ... @@ -1597,7 +1532,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1597 1532 stim::vec<float> sdir = getDirection();
1598 1533  
1599 1534 Bind();
1600   - sk.Begin(stim::OBJ_LINE);
  1535 +// sk.Begin(stim::OBJ_LINE);
1601 1536  
1602 1537  
1603 1538 // sk.createFromSelf(GL_SELECT);
... ... @@ -1618,7 +1553,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1618 1553 if (cost > min_cost){
1619 1554 std::cout << "Cost Limit" << std::endl;
1620 1555 running = false;
1621   - sk.End();
  1556 +// sk.End();
1622 1557 branchDetection2();
1623 1558 pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), -1);
1624 1559 addToNetwork(a, spos, smag, sdir);
... ... @@ -1634,7 +1569,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1634 1569 std::cout << "Edge Limit" << std::endl;
1635 1570 // std::cout << "Found Edge" << std::endl;
1636 1571 running = false;
1637   - sk.End();
  1572 +// sk.End();
1638 1573 branchDetection2();
1639 1574 pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), -1);
1640 1575 addToNetwork(a, spos, smag, sdir);
... ... @@ -1654,7 +1589,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1654 1589 if(mag[0] > 75 || mag[0] < 1){
1655 1590 std::cout << "Magnitude Limit" << std::endl;
1656 1591 running = false;
1657   - sk.End();
  1592 +// sk.End();
1658 1593 branchDetection2();
1659 1594 pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), -1);
1660 1595 addToNetwork(a, spos, smag, sdir);
... ... @@ -1668,7 +1603,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1668 1603 if(h != -1){
1669 1604 std::cout << "Hit Limit" << std::endl;
1670 1605 running = false;
1671   - sk.End();
  1606 +// sk.End();
1672 1607 branchDetection2();
1673 1608 pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), h);
1674 1609 addToNetwork(a, spos, smag, sdir);
... ... @@ -1680,8 +1615,8 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1680 1615 cM.push_back(stim::vec<float>(m[0], m[0]));
1681 1616 // cM.push_back(m[0]);
1682 1617  
1683   - sk.TexCoord(m[0]);
1684   - sk.Vertex(p[0], p[1], p[2]);
  1618 +// sk.TexCoord(m[0]);
  1619 +// sk.Vertex(p[0], p[1], p[2]);
1685 1620 Bind(btexbufferID, bfboID, 27);
1686 1621 CHECK_OPENGL_ERROR
1687 1622 // branchDetection();
... ...
stim/visualization/aaboundingbox.h 0 → 100644
  Show/Hide comments   View file @8761b64
  1 +#ifndef STIM_AABB
  2 +#define STIM_AABB
  3 +
  4 +namespace stim{
  5 +
  6 +
  7 +/// This class describes a structure for an axis-aligned bounding box
  8 +template< typename T >
  9 +class aaboundingbox{
  10 +
  11 +public:
  12 + bool set; //has the bounding box been set to include any points?
  13 + stim::vec<T> A; //minimum point in the bounding box
  14 + stim::vec<T> B; //maximum point in the bounding box
  15 +
  16 + aaboundingbox(){ //constructor generates an empty bounding box
  17 + set = false;
  18 + }
  19 +
  20 +
  21 + /// Test if a point is inside of the bounding box and returns true if it is.
  22 +
  23 + /// @param p is the point to be tested
  24 + bool test(stim::vec<T> p){
  25 +
  26 + for(unsigned d = 0; d < p.size(); p++){ //for each dimension
  27 + if(p[d] < A[d]) return false; //if the point is less than the minimum bound, return false
  28 + if(p[d] > B[d]) return false; //if the point is greater than the max bound, return false
  29 + }
  30 + return true;
  31 + }
  32 +
  33 + /// Expand the bounding box to include the specified point.
  34 +
  35 + /// @param p is the point to be included
  36 + void expand(stim::vec<T> p){
  37 +
  38 + if(!set){ //if the bounding box is empty, fill it with the current point
  39 + A = B = p;
  40 + set = true;
  41 + }
  42 +
  43 + for(unsigned d = 0; d < p.size(); d++){ //for each dimension
  44 + if(p[d] < A[d]) A[d] = p[d]; //expand the bounding box as necessary
  45 + if(p[d] > B[d]) B[d] = p[d];
  46 + }
  47 + }
  48 +
  49 + /// Return the center point of the bounding box as a stim::vec
  50 + stim::vec<T> center(){
  51 + return (B + A) * 0.5;
  52 + }
  53 +
  54 + /// Return the size of the bounding box as a stim::vec
  55 + stim::vec<T> size(){
  56 + return (B - A);
  57 + }
  58 +
  59 + /// Generate a string for the bounding box
  60 + std::string str(){
  61 + std::stringstream ss;
  62 + ss<<A.str()<<"----->"<<B.str();
  63 + return ss.str();
  64 + }
  65 +
  66 +
  67 +}; //end stim::aabb
  68 +
  69 +
  70 +}; //end namespace stim
  71 +
  72 +#endif
0 73 \ No newline at end of file
... ...
stim/visualization/gl_aaboundingbox.h 0 → 100644
  Show/Hide comments   View file @8761b64
  1 +#ifndef STIM_GL_AABB
  2 +#define STIM_GL_AABB
  3 +
  4 +#include <stim/visualization/aaboundingbox.h>
  5 +#include <GL/gl.h>
  6 +
  7 +namespace stim{
  8 +
  9 +template <typename T>
  10 +class gl_aaboundingbox : public aaboundingbox<T>{
  11 +
  12 +public:
  13 +
  14 + //default constructor
  15 + gl_aaboundingbox() : stim::aaboundingbox<T>(){}
  16 +
  17 + //constructor takes an AABB
  18 + gl_aaboundingbox(stim::aaboundingbox<T> b) : stim::aaboundingbox<T>(b){}
  19 +
  20 +
  21 + /// Specifies vertices of the bounding box using CW winding. Use GL_LINE_LOOP for wireframe or GL_QUADS for a solid.
  22 + void glWire(){
  23 +
  24 + //front plane (in A[2])
  25 + glBegin(GL_LINE_LOOP);
  26 + glVertex3f(A[0], A[1], A[2]);
  27 + glVertex3f(A[0], B[1], A[2]);
  28 + glVertex3f(B[0], B[1], A[2]);
  29 + glVertex3f(B[0], A[1], A[2]);
  30 + glEnd();
  31 +
  32 + //back plane (in B[2])
  33 + glBegin(GL_LINE_LOOP);
  34 + glVertex3f(B[0], B[1], B[2]);
  35 + glVertex3f(A[0], B[1], B[2]);
  36 + glVertex3f(A[0], A[1], B[2]);
  37 + glVertex3f(B[0], A[1], B[2]);
  38 + glEnd();
  39 +
  40 + //fill out the rest of the lines to connect the two faces
  41 + glBegin(GL_LINES);
  42 + glVertex3f(A[0], B[1], A[2]);
  43 + glVertex3f(A[0], B[1], B[2]);
  44 + glVertex3f(B[0], B[1], B[2]);
  45 + glVertex3f(B[0], B[1], A[2]);
  46 + glVertex3f(B[0], A[1], A[2]);
  47 + glVertex3f(B[0], A[1], B[2]);
  48 + glVertex3f(A[0], A[1], B[2]);
  49 + glVertex3f(A[0], A[1], A[2]);
  50 + glEnd();
  51 +
  52 + }
  53 +
  54 +
  55 +}; //end stim::gl_aabb
  56 +
  57 +
  58 +}; //end namespace stim
  59 +
  60 +#endif
0 61 \ No newline at end of file
... ...