codebase / stimlib

Commit 7f27eafa7dfb882cd7c240ce1cb480e5da10c024

Authored by David Mayerich
1 parent 1ab6c061

simplified the stim::network class, added the ability to load a network from an OBJ file

Showing 4 changed files with 636 additions and 1306 deletions   Show diff stats
Inline Side-by-side
stim/visualization/fiber.h renamed to stim/biomodels/fiber.h
  Show/Hide comments   View file @7f27eaf
... ... @@ -13,6 +13,7 @@ namespace stim{
13 13 template<typename T>
14 14 class fiber{
15 15  
  16 +protected:
16 17 unsigned int N; //number of points in the fiber
17 18 double **c; //centerline (array of double pointers)
18 19  
... ... @@ -70,11 +71,16 @@ class fiber{
70 71 kdt = new ANNkd_tree(pts, n_data, 3); //build a KD tree
71 72 }
72 73  
  74 +
73 75 double dist(double* p0, double* p1){
74 76  
75 77 double sum = 0;
76   - for(unsigned int d = 0; d < 3; d++)
77   - sum += p0[d] * p1[d];
  78 + float v;
  79 + for(unsigned int d = 0; d < 3; d++){
  80 + v = p1[d] - p0[d];
  81 + sum +=v * v;
  82 +
  83 + }
78 84  
79 85 return sqrt(sum);
80 86 }
... ... @@ -92,6 +98,19 @@ class fiber{
92 98 return *idx;
93 99 }
94 100  
  101 + /// Returns a stim::vec representing the point at index i
  102 +
  103 + /// @param i is an index of the desired centerline point
  104 + stim::vec<T> get_vec(unsigned i){
  105 + stim::vec<T> r;
  106 + r.resize(3);
  107 + r[0] = c[i][0];
  108 + r[1] = c[i][1];
  109 + r[2] = c[i][2];
  110 +
  111 + return r;
  112 + }
  113 +
95 114  
96 115  
97 116  
... ... @@ -114,6 +133,25 @@ public:
114 133 init(n);
115 134 }
116 135  
  136 + /// Constructor takes a list of stim::vec points, the radius at each point is set to zero
  137 + fiber(std::vector< stim::vec<T> > p){
  138 + init(p.size()); //initialize the fiber
  139 +
  140 + //for each point, set the centerline position and radius
  141 + for(unsigned int i = 0; i < N; i++){
  142 +
  143 + //set the centerline position
  144 + for(unsigned int d = 0; d < 3; d++)
  145 + c[i][d] = (double) p[i][d];
  146 +
  147 + //set the radius
  148 + r[i] = 0;
  149 + }
  150 +
  151 + //generate a kd tree
  152 + gen_kdtree();
  153 + }
  154 +
117 155 /// constructor takes a list of points and radii
118 156 fiber(std::vector< stim::vec< T > > pos, std::vector< T > radii){
119 157 init(pos.size()); //initialize the fiber
... ... @@ -440,6 +478,19 @@ public:
440 478 unsigned int n_pts(){
441 479 return N;
442 480 }
  481 +
  482 + /// Bracket operator returns the element at index i
  483 +
  484 + /// @param i is the index of the element to be returned as a stim::vec
  485 + stim::vec<T> operator[](unsigned i){
  486 + return get_vec(i);
  487 + }
  488 +
  489 + /// Back method returns the last point in the fiber
  490 + stim::vec<T> back(){
  491 + return get_vec(N-1);
  492 + }
  493 +
443 494 };
444 495  
445 496  
... ...
stim/biomodels/network.h
  Show/Hide comments   View file @7f27eaf
1 1 #ifndef STIM_NETWORK_H
2 2 #define STIM_NETWORK_H
3 3  
  4 +#include <list>
  5 +#include <stdlib.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/biomodels/fiber.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){
  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 fiber<T>
  32 + {
  33 + public:
  34 + unsigned v[2]; //unique id's designating the starting and ending
83 35  
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   -
137   - }
  36 + /// Constructor - creates an edge from a list of points by calling the stim::fiber constructor
138 37  
  38 + ///@param p is a position in space
  39 + edge(std::vector< stim::vec<T> > p) : fiber<T>(p){}
  40 +
  41 + /// Output the edge information as a string
139 42 std::string str(){
140 43 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   -
  44 + ss<<"("<<fiber<T>::N<<")\tl = "<<length()<<"\t"<<v[0]<<"----"<<v[1];
148 45 return ss.str();
149 46 }
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   -
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 47  
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<<", ";
186   -
187   - for(f = out.begin(); f != out.end(); f++){
  48 + };
  49 +
  50 + ///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.
  51 + class vertex : public stim::vec<T>
  52 + {
  53 + public:
  54 + //std::vector<unsigned int> edges; //indices of edges connected to this node.
  55 + std::vector<unsigned int> e[2]; //indices of edges going out (e[0]) and coming in (e[1])
  56 + //stim::vec<T> p; //position of this node in physical space.
  57 +
  58 + //constructor takes a stim::vec
  59 + vertex(stim::vec<T> p) : stim::vec<T>(p){}
  60 +
  61 + /// Output the vertex information as a string
  62 + std::string str(){
  63 + std::stringstream ss;
  64 + ss<<"\t(x, y, z) = "<<stim::vec<T>::str();
  65 +
  66 + if(e[0].size() > 0){
  67 + ss<<"\t> ";
  68 + for(unsigned int o = 0; o < e[0].size(); o++)
  69 + ss<<e[0][o]<<" ";
  70 + }
  71 + if(e[1].size() > 0){
  72 + ss<<"\t< ";
  73 + for(unsigned int i = 0; i < e[1].size(); i++)
  74 + ss<<e[1][i]<<" ";
  75 + }
188 76  
189   - if(f != out.begin())
190   - ss<<", ";
191   - ss<<(*f)->n[1]->id;
  77 + return ss.str();
192 78 }
193   -
194   -
195   - return ss.str();
196   -
197   -
198   -
199   - }
  79 +
200 80 };
201 81  
  82 + private:
202 83  
203   -//---------------NETWORK CLASS-----------------------------
204   -
205   -protected:
206   -
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++;
  84 + std::vector<edge> E; //list of edges
  85 + std::vector<vertex> V; //list of vertices.
  86 +
  87 + public:
221 88  
222   - i = 0;
223   - for(fiber_i fi = F.begin(); fi != F.end(); fi++)
224   - fi->id = i++;
  89 + ///Returns the number of edges in the network.
  90 + unsigned int edges(){
  91 + return E.size();
225 92 }
226 93  
227   -public:
228   -
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   - }
257   -
258   - return ss.str();
  94 + ///Returns the number of nodes in the network.
  95 + unsigned int vertices(){
  96 + return V.size();
259 97 }
260 98  
261   - /// Load a network from an OBJ object
  99 + //load a network from an OBJ file
  100 + void load_obj(std::string filename){
262 101  
263   - /// @param object is the object file to be used as the basis for the network
264   - void load( stim::obj<T> object){
  102 + stim::obj<T> O; //create an OBJ object
  103 + O.load(filename); //load the OBJ file as an object
265 104  
266   - //get the number of vertices in the object
267   - unsigned int nV = object.numV();
  105 + //grab each line in the OBJ object - these will be fibers
  106 + std::cout<<O.str()<<std::endl;
268 107  
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());
  108 + std::vector<unsigned> id2vert; //this list stores the OBJ vertex ID associated with each network vertex
272 109  
273   - unsigned int nL = object.numL(); //get the number of lines in the OBJ
  110 + unsigned i[2]; //temporary, IDs associated with the first and last points in an OBJ line
274 111  
275   - //for each line in the OBJ structure
276   - for(unsigned int li = 0; li < nL; li++){
  112 + //for each line in the OBJ object
  113 + for(unsigned int l = 1; l <= O.numL(); l++){
277 114  
278   - F.push_back(fiber()); //push a new fiber onto the fiber list
  115 + std::vector< stim::vec<T> > c; //allocate an array of points for the vessel centerline
  116 + O.getLine(l, c); //get the fiber centerline
279 117  
280   - fiber_i f = --(F.end()); //get an iterator to the new fiber
  118 + edge e(c); //create an edge from the given centerline
281 119  
282   - //----------Handle the terminating nodes for the fiber
  120 + //get the first and last vertex IDs for the line
  121 + std::vector< unsigned > id; //create an array to store the centerline point IDs
  122 + O.getLinei(l, id); //get the list of point IDs for the line
  123 + i[0] = id.front(); //get the OBJ ID for the first element of the line
  124 + i[1] = id.back(); //get the OBJ ID for the last element of the line
283 125  
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;
  126 + std::vector<unsigned>::iterator it; //create an iterator for searching the id2vert array
  127 + unsigned it_idx; //create an integer for the id2vert entry index
288 128  
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
  129 + //find out if the nodes for this fiber have already been created
  130 + it = find(id2vert.begin(), id2vert.end(), i[0]); //look for the first node
  131 + it_idx = std::distance(id2vert.begin(), it);
  132 + if(it == id2vert.end()){ //if i[0] hasn't already been used
  133 + vertex v = e[0]; //create a new vertex, assign it a position
  134 + v.e[0].push_back(E.size()); //add the current edge as outgoing
  135 + e.v[0] = V.size(); //add the new vertex to the edge
  136 + V.push_back(v); //add the new vertex to the vertex list
  137 + id2vert.push_back(i[0]); //add the ID to the ID->vertex conversion list
293 138 }
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
  139 + else{ //if the vertex already exists
  140 + V[it_idx].e[0].push_back(E.size()); //add the current edge as outgoing
  141 + e.v[0] = it_idx;
300 142 }
301 143  
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);
  144 + it = find(id2vert.begin(), id2vert.end(), i[1]); //look for the second ID
  145 + it_idx = std::distance(id2vert.begin(), it);
  146 + if(it == id2vert.end()){ //if i[1] hasn't already been used
  147 + vertex v = e.back(); //create a new vertex, assign it a position
  148 + v.e[1].push_back(E.size()); //add the current edge as incoming
  149 + e.v[1] = V.size();
  150 + V.push_back(v); //add the new vertex to the vertex list
  151 + id2vert.push_back(i[1]); //add the ID to the ID->vertex conversion list
306 152 }
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);
  153 + else{ //if the vertex already exists
  154 + V[it_idx].e[1].push_back(E.size()); //add the current edge as incoming
  155 + e.v[1] = it_idx;
313 156 }
314 157  
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);
  158 + E.push_back(e); //push the edge to the list
318 159  
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 160 }
327   -
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();
349   -
350   - vi++; //increment the array index
351   - }
352   -
353   - //return the resulting array
354   - return result;
355 161 }
356 162  
357   - std::vector< stim::vec<T> > get_fiber_geometry( fiber_i f ){
358   - return f->geometry();
359   - }
360   -
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();
370   -
371   - //retrieve a list of terminal node positions
372   - std::vector< stim::vec<T> > node_pos = get_node_positions();
373   -
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++){
383   -
384   - //get an array of fiber points
385   - std::vector< stim::vec<T> > fiber_p = get_fiber_geometry(fi);
386   -
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   -
392   - //add this subset to the geometry object
393   - nP = object.addV(fiber_subset);
394   -
395   - //create an array to hold vertex indices for a line
396   - std::vector<unsigned int> line;
397   - line.resize(fiber_p.size());
398   -
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;
402   -
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   - }
407   -
408   - //add the line list to the object class
409   - object.addLine(line);
  163 + /// Output the network as a string
  164 + std::string str(){
410 165  
  166 + std::stringstream ss;
  167 + ss<<"Nodes ("<<V.size()<<")--------"<<std::endl;
  168 + for(unsigned int v = 0; v < V.size(); v++){
  169 + ss<<"\t"<<v<<V[v].str()<<std::endl;
411 170 }
412 171  
413   - return object;
414   - }
415   -
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){
423   -
424   - //determine the number of fibers in the network
425   - unsigned int nF = F.size();
426   -
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);
432   -
433   - //assign names (identifiers) to the network components
434   - set_names();
435   -
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;
442   -
443   - r = f->radius(l); //get the length and radius of the capillary (calculated at the same time)
444   -
445   - radius[i] = r; //store the radius in the output array
446   - length[i] = l; //store the length in the output array
447   -
448   - i++; //increment the array index
  172 + ss<<"Edges ("<<E.size()<<")--------"<<std::endl;
  173 + for(unsigned e = 0; e < E.size(); e++){
  174 + ss<<"\t"<<e<<E[e].str()<<std::endl;
449 175 }
450 176  
451   -
  177 + return ss.str();
452 178 }
453   -
454   -};
455   -
456   -}; //end namespace stim
457   -
458   -
  179 +}; //end stim::network class
  180 +}; //end stim namespace
459 181 #endif
... ...
stim/biomodels/network_dep.h 0 → 100644
  Show/Hide comments   View file @7f27eaf
  1 +#ifndef STIM_NETWORK_H
  2 +#define STIM_NETWORK_H
  3 +
  4 +#include <stim/math/vector.h>
  5 +#include <stim/visualization/obj.h>
  6 +#include <list>
  7 +#include <ANN/ANN.h>
  8 +
  9 +namespace stim{
  10 +
  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 +
  18 +template<typename T>
  19 +class network{
  20 +
  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 +
  137 + }
  138 +
  139 + std::string str(){
  140 + 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 +
  148 + return ss.str();
  149 + }
  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 +
  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<<", ";
  186 +
  187 + for(f = out.begin(); f != out.end(); f++){
  188 +
  189 + if(f != out.begin())
  190 + ss<<", ";
  191 + ss<<(*f)->n[1]->id;
  192 + }
  193 +
  194 +
  195 + return ss.str();
  196 +
  197 +
  198 +
  199 + }
  200 + };
  201 +
  202 +
  203 +//---------------NETWORK CLASS-----------------------------
  204 +
  205 +protected:
  206 +
  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 + }
  226 +
  227 +public:
  228 +
  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 + }
  257 +
  258 + return ss.str();
  259 + }
  260 +
  261 + /// Load a network from an OBJ object
  262 +
  263 + /// @param object is the object file to be used as the basis for the network
  264 + void load( stim::obj<T> object){
  265 +
  266 + //get the number of vertices in the object
  267 + unsigned int nV = object.numV();
  268 +
  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());
  272 +
  273 + unsigned int nL = object.numL(); //get the number of lines in the OBJ
  274 +
  275 + //for each line in the OBJ structure
  276 + for(unsigned int li = 0; li < nL; li++){
  277 +
  278 + F.push_back(fiber()); //push a new fiber onto the fiber list
  279 +
  280 + fiber_i f = --(F.end()); //get an iterator to the new fiber
  281 +
  282 + //----------Handle the terminating nodes for the fiber
  283 +
  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;
  288 +
  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
  293 + }
  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
  300 + }
  301 +
  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);
  306 + }
  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);
  313 + }
  314 +
  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);
  318 +
  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 + }
  327 +
  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();
  349 +
  350 + vi++; //increment the array index
  351 + }
  352 +
  353 + //return the resulting array
  354 + return result;
  355 + }
  356 +
  357 + std::vector< stim::vec<T> > get_fiber_geometry( fiber_i f ){
  358 + return f->geometry();
  359 + }
  360 +
  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();
  370 +
  371 + //retrieve a list of terminal node positions
  372 + std::vector< stim::vec<T> > node_pos = get_node_positions();
  373 +
  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++){
  383 +
  384 + //get an array of fiber points
  385 + std::vector< stim::vec<T> > fiber_p = get_fiber_geometry(fi);
  386 +
  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 +
  392 + //add this subset to the geometry object
  393 + nP = object.addV(fiber_subset);
  394 +
  395 + //create an array to hold vertex indices for a line
  396 + std::vector<unsigned int> line;
  397 + line.resize(fiber_p.size());
  398 +
  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;
  402 +
  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 + }
  407 +
  408 + //add the line list to the object class
  409 + object.addLine(line);
  410 +
  411 + }
  412 +
  413 + return object;
  414 + }
  415 +
  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){
  423 +
  424 + //determine the number of fibers in the network
  425 + unsigned int nF = F.size();
  426 +
  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);
  432 +
  433 + //assign names (identifiers) to the network components
  434 + set_names();
  435 +
  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;
  442 +
  443 + r = f->radius(l); //get the length and radius of the capillary (calculated at the same time)
  444 +
  445 + radius[i] = r; //store the radius in the output array
  446 + length[i] = l; //store the length in the output array
  447 +
  448 + i++; //increment the array index
  449 + }
  450 +
  451 +
  452 + }
  453 +
  454 +};
  455 +
  456 +}; //end namespace stim
  457 +
  458 +
  459 +#endif
... ...
stim/visualization/network.h deleted
View file @1ab6c06
1   -#ifndef STIM_NETWORK_H
2   -#define STIM_NETWORK_H
3   -
4   -#include <list>
5   -#include <stdlib.h>
6   -#include <sstream>
7   -#include <fstream>
8   -#include <algorithm>
9   -#include <string.h>
10   -#include <math.h>
11   -#include <stim/math/vector.h>
12   -#include <stim/visualization/obj.h>
13   -#include <stim/visualization/fiber.h>
14   -#include <ANN/ANN.h>
15   -#include <boost/tuple/tuple.hpp>
16   -
17   -
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   -*/
24   -
25   -
26   -template<typename T>
27   -class network{
28   -
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 fiber<T>
32   - {
33   - public:
34   - int Nodes[2]; //unique id's designating the starting and ending
35   -
36   - ///default constructor
37   - edge() : fiber<T>()
38   - {
39   - Nodes[1] = -1; Nodes[2] = -1;
40   - }
41   -
42   - ///sized costructor with two known nodes.
43   - ///@param startId: int storing idx assigned to Nodes[1].
44   - ///@param endId: int storing idx assigned to Nodes[2].
45   - ///@param n: int for the number of points in the fiber.
46   - edge(int startId, int endId, int n)
47   - :fiber<T>(n)
48   - {
49   - Nodes[1] = startId; Nodes[2] = endId;
50   - }
51   -
52   - ///constructor from a std::vector of stim::vecs of positions and radii.
53   - ///@param pos: Vector of stim vecs storing the positions of the fiber.
54   - ///@param mag: Vector of stim vecs storing the radii of the fiber.
55   - edge(std::vector< stim::vec<T> > pos, std::vector< stim::vec<T> > radii)
56   - : fiber<T>(pos, radii)
57   - {
58   - Nodes[1] = -1; Nodes[2] = -1;
59   - }
60   -
61   - ///constructor from an std::vector of stim::vecs of positions and radii as well as a known starting and ending node.
62   - ///@param pos: Vector of stim vecs storing the positions of the fiber.
63   - ///@param mag: Vector of stim vecs storing the radii of the fiber.
64   - ///@param id1: int storing idx assigned to Nodes[1].
65   - ///@param id2: int storing idx assigned to Nodes[2].
66   - edge(std::vector< stim::vec<T> > pos, std::vector< stim::vec<T> > radii, int id1, int id2)
67   - : fiber<T>(pos, radii)
68   - {
69   - Nodes[1] = id1; Nodes[2] = id2;
70   - }
71   -
72   -
73   - edge(std::vector< stim::vec<T> > pos, std::vector<T> radii)
74   - : fiber<T>(pos, radii)
75   - {
76   - Nodes[1] = -1; Nodes[2] = -1;
77   - }
78   -
79   - ///sets the endpoints to the two int values.
80   - ///@param int id1: index of Nodes[1].
81   - ///@param int id2: index of Nodes[2].
82   - void
83   - setEndpoints(int id1, int id2)
84   - {
85   - Nodes[1] = id1; Nodes[2] = id2;
86   - }
87   -
88   - };
89   -
90   - ///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.
91   - class node
92   - {
93   - private:
94   - std::vector<int> edges; //indices of edges connected to this node.
95   - stim::vec<T> p; //position of this node in physical space.
96   - public:
97   - //no default constructor
98   -
99   - ///@param pos: stim vec with the x, y, z position of the edge.
100   - ///stim::vec constructure with a position but no attached edges.
101   - node(stim::vec<T> pos)
102   - {
103   - p = pos;
104   - }
105   -
106   - ///@param pos: stim vec with the x, y, z position of the edge.
107   - ///@param i: int i storing the index of an attached edge.
108   - //stim::vec constructor with a position and an attached edge.
109   - node(stim::vec<T> pos, int i)
110   - {
111   - p = pos;
112   - edges.push_back(i);
113   - }
114   -
115   - ///@param x: x coordinate of the node..
116   - ///@param y: y coordinate of the node..
117   - ///@param z: z coordinate of the node..
118   - //float value constructor.
119   - node(T x, T y, T z)
120   - {
121   - p = stim::vec<T>(x,y,z);
122   - }
123   - ///@param x: x coordinate of the node..
124   - ///@param y: y coordinate of the node..
125   - ///@param z: z coordinate of the node..
126   - ///@param i: int i storing the index of an attached edge.
127   - //float value constructor with an attached edge.
128   - node(T x, T y, T z, int i)
129   - {
130   - p = stim::vec<T>(x,y,z);
131   - edges.push_back(i);
132   - }
133   -
134   - ///@param id: int index of the fiber to attach to this node.
135   - ///adds the fiber to the rest of the fibers connected to this node.
136   - void
137   - addEdge(int id)
138   - {
139   - edges.push_back(id);
140   - }
141   -
142   - ///returns the position of the node.
143   - stim::vec<T>
144   - getPosition()
145   - {
146   - return p;
147   - }
148   -
149   - ///@param id: int index of the fiber to detach to this node.
150   - ///removes the edge from the list of the edges attached to this node.
151   - void
152   - removeEdge(int id)
153   - {
154   - for(int i = 0; i < edges.size(); i++)
155   - {
156   - if(edges[i] == id)
157   - edges.erase(edges.begin()+i);
158   - }
159   - }
160   -
161   - ///returns and std::string with the position of this node.
162   - std::string
163   - str()
164   - {
165   - return p.str();
166   - }
167   -
168   - ///returns and std::string with the list of edges attached to this node.
169   - std::string
170   - edges_to_str()
171   - {
172   - std::ostringstream ss;
173   - std::cout<<"here"<<std::endl;
174   -// ss << "[";
175   - for(int i = 0; i < edges.size()-1; i++)
176   - {
177   - std::cout<<"here"<<i<<std::endl;
178   - ss << edges[i] << ";";
179   - std::cout<<edges.size()-1<<std::endl;
180   - }
181   - ss << edges[edges.size()-1];
182   -// ss << "]";
183   - return ss.str();
184   - }
185   -
186   - };
187   -
188   - public:
189   -
190   - std::vector<edge*> E; //list of pointers to edges.
191   - std::vector<node> V; //list of nodes.
192   - std::vector< stim::vec<T> > allVerticesList; // all nodes before sampling
193   - std::vector<stim::vec<T> > allVerticesAfterSampling ; //all vertices after sampling
194   -
195   - ///Returns the number of edges in the network.
196   - unsigned int
197   - sizeE()
198   - {
199   - return E.size();
200   - }
201   -
202   - ///Returns the number of nodes in the network.
203   - unsigned int
204   - sizeV()
205   - {
206   - return V.size();
207   - }
208   - unsigned int
209   - sizeAfterSamplingV()
210   - {
211   - return allVerticesAfterSampling.size();
212   - }
213   -/* //adds an edge from two std::vectors with positions and radii.
214   - void
215   - addEdge(std::vector< stim::vec<T> > pos, std::vector<stim::vec<T> > radii, int endId)
216   - {
217   -
218   - edge a(pos,radii, endId);
219   - E.push_back(a);
220   - } */
221   -
222   - ///A complicated method that adds an edge to the network.
223   - ///Most of this functionality will be moved into fiber.
224   - void
225   - addEdge(std::vector< stim::vec<T> > pos, std::vector<stim::vec<T> > radii, int startId, int endId)
226   - {
227   - //
228   - if(startId == -1 && endId == -1)
229   - {
230   - //case where the edge is neither starting nor ending in a fiber.
231   - //i. first fiber.
232   -
233   - //Add two nodes to the nodes vector
234   - V.push_back(node(pos[pos.size()-1]));
235   - V.push_back(node(pos[0]));
236   -
237   - //the edge will be connected to the nodes
238   - edge *a = new edge(pos,radii,(V.size()-2), (V.size()-1));
239   -
240   - //add fiber to fiber list.
241   - E.push_back(a);
242   -
243   - //The last two nodes added to V will "own" the last edge added to E.
244   - V[V.size()-1].addEdge(E.size()-1);
245   - V[V.size()-2].addEdge(E.size()-1);
246   -
247   - }
248   -
249   - else if(startId != -1 && endId == -1)
250   - {
251   - //case where the edge is starting with a fiber, but not ending in one.
252   -
253   - //split the fiber behind you into two.
254   - unsigned int point = E[startId]->nearest_idx(pos[0]);
255   -
256   - //split the hit fiber at point two parts temp[0], temp[1]
257   - std::vector < stim::fiber <T> > temp = E[startId]->split(point);
258   - if(temp.size() > 1)
259   - {
260   - //add the nearest point in the behind fiber into the hitting fiber.
261   - pos.insert(pos.begin(), E[startId]->nearest(pos[0]));
262   - stim::vec<T> v(E[startId]->radius(point), E[startId]->radius(point));
263   - radii.insert(radii.begin(), v);
264   -
265   - //reset the fiber at the endId to be a shorter fiber(temp[0]).
266   - std::vector<stim::vec<T> > ce = temp[0].centerline();
267   - std::vector<stim::vec<T> > cm = temp[0].centerlinemag();
268   -
269   - //remake the edge, such that the starting point of this edge
270   - //is the same the split point, but the end edge is the old end.
271   - V.push_back(node(ce[ce.size()-1]));
272   - int tempNodeId = E[startId]->Nodes[1];
273   - E[startId] = new edge(ce, cm, (V.size()-1), E[startId]->Nodes[2]);
274   - V[V.size()-1].addEdge(startId);
275   -
276   -
277   - //add the other part of the fiber (temp[1])
278   - ce = temp[1].centerline();
279   - cm = temp[1].centerlinemag();
280   - E.push_back(new edge(ce, cm,tempNodeId ,(V.size()-1)));
281   - V[V.size()-1].addEdge(E.size()-1);
282   -
283   - V[tempNodeId].removeEdge(startId);
284   - V[tempNodeId].addEdge(E.size()-1);
285   -// V[V.size()-1].removeEdge(startId);
286   -
287   - //make the new hitting fiber..
288   - V.push_back(node(pos[pos.size()-1]));
289   - edge *a = new edge(pos, radii, (V.size()-2), (V.size()-1));
290   - E.push_back(a);
291   - V[V.size()-1].addEdge(E.size()-1);
292   - V[V.size()-2].addEdge(E.size()-1);
293   - } else {
294   - stim::vec<T> pz = E[startId]->nearest(pos[0]);
295   - if((V[E[startId]->Nodes[1]].getPosition() - pz).len() <
296   - (V[E[startId]->Nodes[2]].getPosition() - pz).len())
297   - {
298   - unsigned int point = E[startId]->nearest_idx(pos[0]);
299   - pos.insert(pos.begin(), E[startId]->nearest(pos[0]));
300   - stim::vec<T> v(E[startId]->radius(point), E[startId]->radius(point));
301   - radii.insert(radii.begin(), v);
302   - V.push_back(node(pos[pos.size()-1]));
303   - edge *a = new edge(pos, radii, E[startId]->Nodes[1], (V.size()-1));
304   - E.push_back(a);
305   -
306   - V[V.size()-1].addEdge(E.size()-1);
307   - V[ E[startId]->Nodes[1]].addEdge(E.size()-1);
308   -
309   - }
310   - else
311   - {
312   - unsigned int point = E[startId]->nearest_idx(pos[0]);
313   - pos.insert(pos.begin(), E[startId]->nearest(pos[0]));
314   - stim::vec<T> v(E[startId]->radius(point), E[startId]->radius(point));
315   - radii.insert(radii.begin(), v);
316   - V.push_back(node(pos[pos.size()-1]));
317   - edge *a = new edge(pos, radii, E[startId]->Nodes[2], (V.size()-1));
318   - E.push_back(a);
319   -
320   - V[V.size()-1].addEdge(E.size()-1);
321   - V[ E[startId]->Nodes[2]].addEdge(E.size()-1);
322   - }
323   - }
324   - }
325   -
326   - //case where the edge is starting at a seedpoint but ends in a fiber.
327   - if(startId == -1 && endId != -1 && endId < sizeE())
328   - {
329   - //split the hit fiber into two.
330   - unsigned int point = E[endId]->nearest_idx(pos[pos.size() -1]);
331   -
332   - //split the hit fiber at point into two parts temp[0], temp[1]
333   - std::vector < stim::fiber <T> > temp
334   - = E[endId]->split(point);
335   - if(temp.size() > 1)
336   - {
337   - //add the nearest point in the hit fiber into the hitting fiber.
338   - pos.push_back(E[endId]->nearest(pos[pos.size()-1]));
339   - // pos.insert(pos.end(), E[endId].nearest(pos[pos.size()-1]));
340   - stim::vec<T> v(E[endId]->radius(point), E[endId]->radius(point));
341   - radii.push_back(v);
342   -
343   - //split the hit fiber at point into two parts temp[0], temp[1]
344   - std::vector < stim::fiber <T> > temp
345   - = E[endId]->split(point);
346   -
347   - //reset the fiber at endId to be a shorter fiber (temp[0]).
348   - std::vector<stim::vec<T> > ce = temp[0].centerline();
349   - std::vector<stim::vec<T> > cm = temp[0].centerlinemag();
350   -
351   - //remake the edge, such that the ending point of this edge
352   - //is the same as before, but the starting edge is new.
353   - V.push_back(node(ce[ce.size()-1])); //split point.
354   - int tempNodeId = E[endId]->Nodes[2];
355   - E[endId] = new edge(ce, cm, E[endId]->Nodes[1], (V.size()-1));
356   - V[V.size()-1].addEdge(endId);
357   -
358   - //add that other part of the fiber (temp[1])
359   - //such that it begins with the middle node, and ends with
360   - //the last node of the previous fiber.
361   - ce = temp[1].centerline();
362   - cm = temp[1].centerlinemag();
363   - E.push_back(new edge(ce, cm, (V.size()-1), tempNodeId));
364   - V[V.size()-1].addEdge(E.size()-1);
365   - // V[V.size()-1].removeEdge(endId);
366   -
367   -
368   -
369   - //make the new hitting fiber..
370   - V.push_back(pos[0]);
371   - edge *a = new edge(pos,radii,(V.size()-1), (V.size()-2));
372   - E.push_back(a);
373   - V[V.size()-1].addEdge(E.size()-1);
374   - V[V.size()-2].addEdge(E.size()-1);
375   -
376   - //in the end we have added two new nodes and two new edges.
377   - }
378   - else {
379   - stim::vec<T> pz = E[endId]->nearest(pos[0]);
380   - if((V[ E[endId]->Nodes[1]].getPosition() - pz).len() <
381   - (V[E[endId]->Nodes[2]].getPosition() - pz).len())
382   - {
383   - unsigned int point = E[endId]->nearest_idx(pos[0]);
384   - pos.insert(pos.begin(), E[endId]->nearest(pos[0]));
385   - stim::vec<T> v(E[endId]->radius(point), E[endId]->radius(point));
386   - radii.insert(radii.begin(), v);
387   - V.push_back(node(pos[pos.size()-1]));
388   - edge *a = new edge(pos, radii, E[endId]->Nodes[1], (V.size()-1));
389   - E.push_back(a);
390   -
391   - V[V.size()-1].addEdge(E.size()-1);
392   - V[ E[endId]->Nodes[1]].addEdge(E.size()-1);
393   -
394   - }
395   - else
396   - {
397   - unsigned int point = E[endId]->nearest_idx(pos[0]);
398   - pos.insert(pos.begin(), E[endId]->nearest(pos[0]));
399   - stim::vec<T> v(E[endId]->radius(point), E[endId]->radius(point));
400   - radii.insert(radii.begin(), v);
401   - V.push_back(node(pos[pos.size()-1]));
402   - edge *a = new edge(pos, radii, E[endId]->Nodes[2], (V.size()-1));
403   - E.push_back(a);
404   -
405   - V[V.size()-1].addEdge(E.size()-1);
406   - V[ E[endId]->Nodes[2]].addEdge(E.size()-1);
407   - }
408   - }
409   - }
410   -
411   - if(startId != -1 && endId != -1 && endId < sizeE())
412   - {
413   - //case where the edge is starting with a fiber, and ends in one.
414   -
415   - //split the fiber behind you into two.
416   - unsigned int point = E[startId]->nearest_idx(pos[0]);
417   -// std::cout << "in merge to both" << std::endl;
418   -
419   - //split the hit fiber at point two parts temp[0], temp[1]
420   - std::vector < stim::fiber <T> > temp = E[startId]->split(point);
421   - if(temp.size() > 1)
422   - {
423   - //extend the previous fiber (guaranteed to be added last) by one
424   - //and add the
425   - pos = E[E.size()-1]->centerline();
426   - radii = E[E.size()-1]->centerlinemag();
427   - pos.insert(pos.begin(), E[startId]->nearest(pos[0]));
428   - stim::vec<T> v(E[startId]->radius(point), E[startId]->radius(point));
429   - radii.insert(radii.begin(), v);
430   - V.erase(V.end());
431   - V.push_back(node(pos[0]));
432   -
433   -
434   - //something weird is going on here.
435   - E[E.size()-1] = new edge(pos, radii, (V.size()-2), (V.size()-1));
436   - V[V.size()-1].addEdge(E.size()-1);
437   -
438   - //reset the fiber at the endId to be a shorter fiber(temp[0]).
439   - std::vector<stim::vec<T> > ce = temp[0].centerline();
440   - std::vector<stim::vec<T> > cm = temp[0].centerlinemag();
441   -
442   -// std::cout << ce.size() << std::endl;
443   - //remake the edge, such that the starting point of this edge
444   - //is the same as before, but the end edge is new.
445   - int tempNodeId = E[startId]->Nodes[1];
446   - E[startId] = new edge(ce, cm, (V.size()-1), E[startId]->Nodes[2]);
447   - V[V.size()-1].addEdge(startId);
448   -
449   - //add the other part of the fiber (temp[1])
450   - ce = temp[1].centerline();
451   - cm = temp[1].centerlinemag();
452   - E.push_back(new edge(ce, cm,tempNodeId, (V.size()-1)));
453   - V[V.size()-1].addEdge(E.size()-1);
454   - V[tempNodeId].removeEdge(startId);
455   - V[tempNodeId].addEdge(E.size()-1);
456   -// V[V.size()-1].removeEdge(startId);
457   - }
458   - else {
459   - stim::vec<T> pz = E[endId]->nearest(pos[0]);
460   - if((V[ E[endId]->Nodes[1]].getPosition() - pz).len() <
461   - (V[E[endId]->Nodes[2]].getPosition() - pz).len())
462   - {
463   - unsigned int point = E[endId]->nearest_idx(pos[0]);
464   - pos.insert(pos.begin(), E[endId]->nearest(pos[0]));
465   - stim::vec<T> v(E[endId]->radius(point), E[endId]->radius(point));
466   - radii.insert(radii.begin(), v);
467   - V.push_back(node(pos[pos.size()-1]));
468   - edge *a = new edge(pos, radii, E[endId]->Nodes[1], (V.size()-1));
469   - E.push_back(a);
470   -
471   - V[V.size()-1].addEdge(E.size()-1);
472   - V[ E[endId]->Nodes[1]].addEdge(E.size()-1);
473   -
474   - }
475   - else
476   - {
477   - unsigned int point = E[endId]->nearest_idx(pos[0]);
478   - pos.insert(pos.begin(), E[endId]->nearest(pos[0]));
479   - stim::vec<T> v(E[endId]->radius(point), E[endId]->radius(point));
480   - radii.insert(radii.begin(), v);
481   - V.push_back(node(pos[pos.size()-1]));
482   - edge *a = new edge(pos, radii, E[endId]->Nodes[2], (V.size()-1));
483   - E.push_back(a);
484   -
485   - V[V.size()-1].addEdge(E.size()-1);
486   - V[ E[endId]->Nodes[2]].addEdge(E.size()-1);
487   - }
488   - }
489   - }
490   -
491   - }
492   -
493   - ///@param pos: std::vector of stim vecs with the positions of the point in the fiber.
494   - ///@param radii: std::vector of floats with the radii of the fiber at positions in pos.
495   - ///returns the forest as a std::string. For testing only.
496   - std::string
497   - edges_to_str()
498   - {
499   - std::stringstream ss;
500   - for(unsigned int i = 0; i < E.size(); i++)
501   - {
502   - ss << i << ": " << E[i]->str() << std::endl;
503   - }
504   - return(ss.str());
505   - }
506   - // total number of points on all edges!=fibers in the network
507   - int
508   - totalEdges()
509   - {
510   - int totalEdgeVertices=0;int N=0;
511   - for (unsigned int i=0; i < sizeE(); i ++)
512   - {// FOR N points on the fiber N-1 edges are possible
513   - N = E[i]->n_pts();
514   - totalEdgeVertices = totalEdgeVertices + N- 1;
515   - }
516   - return totalEdgeVertices;
517   - }
518   - // sum of all the fiber lengths
519   - float
520   - lengthOfNetwork()
521   - {
522   - float networkLength=0;float N=0;
523   - for (unsigned int i=0; i < sizeE(); i ++)
524   - {
525   - N = E[i]->length();
526   - networkLength = networkLength + N;
527   - }
528   - return networkLength;
529   - }
530   - ///@param i: index of the required fiber.
531   - ///Returns an std::vector with the centerline of the ith fiber in the edgelist.
532   - std::vector< stim::vec<T> >
533   - getEdgeCenterLine(int i)
534   - {
535   - std::vector < stim::vec<T> > a;
536   - return E[i]->centerline();
537   - }
538   -
539   - ///@param i: index of the required fiber.
540   - ///Returns an std::vector with the centerline radii of the ith fiber in the edgelist.
541   - std::vector< stim::vec<T> >
542   - getEdgeCenterLineMag(int i)
543   - {
544   - std::vector < stim::vec<T> > a;
545   - return E[i]->centerlinemag();
546   - }
547   -
548   - ///@param i: index of the required fibers nodes..
549   - ///Returns an std::string with the start and end points of this node..
550   - std::string
551   - nodes_to_str(int i)
552   - {
553   - std::stringstream ss;
554   - ss << "[from Node " << E[i] -> Nodes[1] << " to " << E[i] -> Nodes[2] << "]";
555   - return ss.str();
556   - }
557   - //load an obj file into a network class
558   - stim::network<T>
559   - objToNetwork(stim::obj<T> objInput)
560   - {
561   - stim::network<T> nwc;
562   - //network network2;
563   - // function to convert an obj file loaded using stim/visualization/obj.h
564   - // to a 3D network class using network.h methods.
565   - std::vector< stim::vec<T> > fiberPositions; //initialize positions on the fiber to be added to network
566   - objInput.getLine(1, fiberPositions); int numDim = fiberPositions[0].size();//find dimensions of the vertices in obj file
567   - // to verify if the nodes are already pushed into node list
568   - std::vector<bool> validList;
569   - validList.assign(objInput.numV(), false);
570   - // go through each of the lines "l followed by indices in obj and add all start and end vertices as the nodes
571   - // using addNode function for adding nodes
572   - // and the line as an edge(belongs to fiber class) using addEdge function
573   - std::vector<unsigned> vertexIndices(objInput.numV());
574   - std::vector< stim::vec<T> > vertices = objInput.getAllV(vertexIndices);
575   - nwc.addVertices(vertices);
576   - for (unsigned i =1; i< objInput.numL() + 1; i++)
577   - {
578   - // edges/fibers could be added to the network class
579   - std::vector< stim::vec<T> > fiberPositions;
580   -
581   - objInput.getLine(i, fiberPositions);
582   - // finding size to allocate radii
583   - int numPointsOnFiber = fiberPositions.size();
584   - // to extend it to a 3D postion if it is a 1D/2D vertex in space
585   - std::vector< stim::vec<T> > fiberPositions1(numPointsOnFiber);
586   - // 2D case append and assign the last dimension to zero
587   - if (numDim == 2)
588   - {// 2D case append and assign the last dimension to zero repeating it for all the points on fiber
589   - for (int j = 0; j < numPointsOnFiber; j ++)
590   - {
591   - fiberPositions1[j][numDim - 2] = fiberPositions[j][0];
592   - fiberPositions1[j][numDim -1 ] = fiberPositions[j][1];
593   - fiberPositions1[j][numDim] = 0;
594   - }
595   - }
596   - // 1D case append and assign last two dimensions to zero
597   - else if (numDim == 1)
598   - {
599   - for (int j = 0; j < numPointsOnFiber; j ++)
600   - {
601   - fiberPositions1[j][numDim - 2] = fiberPositions[j][0];
602   - fiberPositions1[j][numDim -1 ] = 0;
603   - fiberPositions1[j][numDim] = 0;
604   - }
605   - }
606   - else
607   - {
608   - fiberPositions1 = fiberPositions;
609   - }
610   - // then add edge
611   - //edge* a = new edge(fiberPositions1,radii);
612   - //std::cout<<"here"<<std::endl;
613   - //E.push_back(a);
614   - std::vector<stim::vec<T> > newPointList;
615   - newPointList = Resample(fiberPositions1);
616   - int numPointsOnnewFiber = newPointList.size();
617   - nwc.addVerticesAfterSamplimg(newPointList);
618   - std::vector<T> radii(numPointsOnnewFiber); // allocating radii to zero
619   - nwc.addEdge(newPointList, radii);
620   - // add nodes
621   - stim::vec<T> v0(3);stim::vec<T> vN(3);
622   - int endId = numPointsOnnewFiber -1;
623   - v0[0] = newPointList[0][0];v0[1] = newPointList[0][1];v0[2] = newPointList[0][2];
624   - vN[0] = newPointList[endId][0];vN[1] = newPointList[endId][1];vN[2] = newPointList[endId][2];
625   - // VISITED INDEXES OF the nodes are set to true
626   - if(!validList[objInput.getIndex(vertices, v0)])
627   - {
628   - //V.push_back(node(v0));
629   - nwc.addNode(v0);
630   - validList[objInput.getIndex(vertices, v0)] = true;
631   - }
632   - if(!validList[objInput.getIndex(vertices, vN)])
633   - {
634   - //V.push_back(node(vN));
635   - nwc.addNode(vN);
636   - validList[objInput.getIndex(vertices, vN)] = true;
637   - }
638   - }
639   - return nwc;
640   - }
641   - // convert ground and test case to network ,kd trees
642   - boost::tuple< ANNkd_tree*, ANNkd_tree*, stim::network<T>, stim::network<T> >
643   - LoadNetworks(std::string gold_filename, std::string test_filename)
644   - {
645   - ANNkd_tree* kdTree1;ANNkd_tree* kdTree2;
646   - using namespace stim;
647   - network network1;network network2;
648   - stim::obj<T> objFile1(gold_filename);
649   - network1 = objToNetwork(objFile1);
650   - kdTree1 = generateKdTreeFromNetwork(network1);
651   - std::cout<<"Gold Standard:network and kdtree generated"<<std::endl;
652   - stim::obj<T> objFile2(test_filename);
653   - network2 = objToNetwork(objFile2);
654   - kdTree2 = generateKdTreeFromNetwork(network2);
655   - std::cout<<"Test Network:network and kdtree generated"<<std::endl;
656   - return boost::make_tuple(kdTree1, kdTree2, network1, network2);
657   - std::cout<<"out of this loop"<<std::endl;
658   - }
659   - //distance between two points
660   - double dist(std::vector<T> p0, std::vector<T> p1)
661   - {
662   -
663   - double sum = 0;
664   - for(unsigned int d = 0; d < 3; d++)
665   - sum += p0[d] * p1[d];
666   - return sqrt(sum);
667   - }
668   - // sum of elements in a vector
669   - double sum(stim::vec<T> metricList)
670   - {
671   - float sumMetricList = 0;
672   - for (unsigned int count=0; count<metricList.size(); count++)
673   - { sumMetricList += metricList[count];}
674   - return sumMetricList;
675   - }
676   - //generate a kd tree from network
677   - ANNkd_tree*
678   - generateKdTreeFromNetwork(stim::network<T> network) //kd-tree stores all points in the fiber for fast searching
679   - {
680   - std::cout<<"kd trees generated"<<std::endl;
681   - ANNkd_tree* kdt;
682   - double **c; //centerline (array of double pointers)
683   - float* r; // array of fiber radii
684   - unsigned int n_data = network.sizeAfterSamplingV(); //set the number of points
685   - c = (double**) malloc(sizeof(double*) * n_data); //allocate the array pointer
686   - for(unsigned int i = 0; i < n_data; i++) //allocate space for each point
687   - {c[i] = (double*) malloc(sizeof(double) * 3);}
688   - r = (float*) malloc(sizeof(float) * n_data); //allocate space for the radii
689   - //stim::vec<T> node;
690   - for(unsigned int i = 0; i < n_data; i++)
691   - {
692   - //node = network.V[i].getPosition();
693   - //convert_to_double
694   - for(unsigned int d = 0; d < 3; d++){ //for each dimension
695   - c[i][d] = double(network.allVerticesAfterSampling[i][d]); //copy the coordinate
696   - }
697   - r[i] = r[i]; //copy the radius
698   - }
699   - ANNpointArray pts = (ANNpointArray)c; //create an array of data points
700   - kdt = new ANNkd_tree(pts, n_data, 3); //build a KD tree
701   - return kdt;
702   - }
703   -
704   - ///@param pos: std::vector of stim vecs with the positions of the point in the fiber.
705   - ///@param radii: std::vector of floats with the radii of the fiber at positions in pos.
706   - ///adds an edge from two std::vectors with positions and radii.
707   - void
708   - addEdge(std::vector< stim::vec<T> > pos, std::vector<T> radii)
709   - {
710   - edge *a = new edge(pos,radii);
711   - E.push_back(a);
712   - }
713   - void
714   - addNode(stim::vec<T> nodes)
715   - {
716   - V.push_back(nodes);
717   - }
718   - void
719   - addVertices(std::vector< stim::vec<T> > vertices)
720   - {
721   - for (unsigned int i=0; i < vertices.size(); i ++)
722   - {
723   - allVerticesList.push_back(vertices[i]);
724   - }
725   - }
726   - void
727   - addVerticesAfterSamplimg(std::vector< stim::vec<T> > vertices)
728   - {
729   - for (unsigned int i=0; i < vertices.size(); i ++)
730   - {
731   - allVerticesAfterSampling.push_back(vertices[i]);
732   - }
733   - }
734   - // gaussian function
735   - float gaussianFunction(float x, float std=25)
736   - {
737   - float evaluate = 1.0f - ((exp(-x/(2*std*std))));
738   - return evaluate;
739   - }
740   - // compares point on a network to a kd tree for two skeletons
741   - boost::tuple< float, float >
742   - compareSkeletons(boost::tuple< ANNkd_tree*, ANNkd_tree*, stim::network<float>, stim::network<float> > networkKdtree)
743   - {
744   - float gFPR, gFNR;
745   - gFPR = CompareNetKD(networkKdtree.get<0>(), networkKdtree.get<3>());
746   - gFNR = CompareNetKD(networkKdtree.get<1>(), networkKdtree.get<2>());
747   - return boost::make_tuple(gFPR, gFNR);
748   - }
749   - // gaussian distance of points on network to Kdtree
750   - float
751   - CompareNetKD(ANNkd_tree* kdTreeGroundtruth, stim::network<T> networkTruth)
752   - {
753   - std::vector< stim::vec<T> > fiberPoints;
754   - float netmetsMetric=0;
755   - double eps = 0; // error bound
756   - ANNdistArray dists1;dists1 = new ANNdist[1]; // near neighbor distances
757   - ANNdistArray dists2; dists2 = new ANNdist[1];// near neighbor distances
758   - ANNidxArray nnIdx1; nnIdx1 = new ANNidx[1]; // near neighbor indices // allocate near neigh indices
759   - ANNidxArray nnIdx2; nnIdx2 = new ANNidx[1]; // near neighbor indices // allocate near neigh indices
760   - int N; int numQueryPoints = networkTruth.totalEdges();
761   - float totalNetworkLength = networkTruth.lengthOfNetwork();
762   - stim::vec<float> fiberMetric(networkTruth.sizeE());
763   - //for each fiber
764   - for (unsigned int i=0; i < networkTruth.sizeE(); i ++)
765   - {
766   - std::vector<T> p1(3); std::vector<T> p2(3);//temporary variables to store point positions
767   - fiberPoints = networkTruth.E[i]->centerline();
768   - N = networkTruth.E[i]->n_pts();
769   - stim::vec<float> segmentMetric(N-1);
770   - // for each segment in the fiber
771   - for (unsigned int j = 0; j < N - 1; j++)
772   - {
773   - ANNpoint queryPt1; queryPt1 = annAllocPt(3);
774   - ANNpoint queryPt2; queryPt2 = annAllocPt(3);
775   - //for each dimension of the point
776   - for(unsigned int d = 0; d < 3; d++)
777   - {
778   - queryPt1[d] = double(fiberPoints[j][d]);
779   - p1[d] = double(fiberPoints[j][d]);
780   - queryPt2[d] = double(fiberPoints[j + 1][d]);
781   - p2[d] = double(fiberPoints[j + 1][d]);// for each dimension
782   - }
783   - kdTreeGroundtruth->annkSearch( queryPt1, 1, nnIdx1, dists1, eps); // error bound
784   - kdTreeGroundtruth->annkSearch( queryPt2, 1, nnIdx2, dists2, eps); // error bound
785   - std::cout<<float(dists1[0])<<"dist1-----"<<float(dists2[0])<<"dist2---"<<std::endl;
786   - float dist1 = gaussianFunction(float(dists1[0]));float dist2 = gaussianFunction(float(dists2[0]));
787   - segmentMetric[j] = (((dist1 + dist2) / 2) * dist(p1, p2)) ;
788   - }
789   - fiberMetric[i] = sum(segmentMetric)/;
790   - }
791   - netmetsMetric = sum(fiberMetric)/totalNetworkLength ;
792   - return netmetsMetric;
793   - }
794   - //resample a fiber in the network
795   - std::vector<stim::vec<T> > Resample(std::vector< stim::vec<T> > fiberPositions, float spacing=25.0)
796   - {
797   - std::vector<T> p1(3), p2(3), v(3);
798   - stim::vec<T> p(3);
799   - std::vector<stim::vec<T> > newPointList;
800   - for(unsigned int f=0; f<fiberPositions.size()-1; f++)
801   - {
802   - T lengthSegment = dist(p1,p2);
803   - for(int d=0; d<3;d++)
804   - {
805   - p1[d] = T(fiberPositions[f][d]);
806   - p2[d] = T(fiberPositions[f + 1][d]);
807   - }// for each dimension
808   - if( lengthSegment >= spacing )
809   - { for (int dim=0; dim<3;dim++) //find the direction of travel
810   - {v[dim] = p2[dim] - p1[dim];}
811   - //set the step size to the voxel size
812   - T step;
813   - for(step=0.0; step<lengthSegment; step+=spacing)
814   - {
815   - for(int i=0; i<3;i++)
816   - {p[i] = p1[i] + v[i]*(step/lengthSegment);}
817   - newPointList.push_back(p);
818   - }
819   - }
820   - else
821   - newPointList = fiberPositions;
822   - }
823   - return newPointList;
824   - }
825   - // modulus of a vector
826   - T
827   - Length(std::vector<T> v)
828   - {
829   - T sum=0;
830   - for (int i=0; i<v.size(); i++)
831   - sum += v[i] * v[i];
832   - return sum;
833   - }
834   - // function to remove characters from a string
835   - void removeCharsFromString(std::string &str, char* charsToRemove ) {
836   - for ( unsigned int i = 0; i < strlen(charsToRemove); ++i ) {
837   - str.erase( remove(str.begin(), str.end(), charsToRemove[i]), str.end() );
838   - }
839   - }
840   - ///exports the network graph to obj
841   - void
842   - to_obj()
843   - {
844   - std::ofstream ofs;
845   - ofs.open("Graph.obj", std::ofstream::out | std::ofstream::app);
846   - int num = allVerticesList.size();
847   - std::string *strArray = new std::string[num];
848   - for(unsigned int i = 0; i < allVerticesList.size(); i++)
849   - {
850   - std::string str;
851   - str = allVerticesList[i].str();
852   - removeCharsFromString(str, "[],");
853   - ofs << "v " << str << "\n";
854   - removeCharsFromString(str," ");
855   - strArray[i] = str;
856   - }
857   - for(unsigned int i = 0; i < E.size(); i++)
858   - {
859   - std::string str;
860   - str = E[i]->strObj(strArray, num);
861   - //removeCharsFromString(str,"0");
862   - ofs << "l " << str << "\n";
863   - }
864   - ofs.close();
865   - }
866   - ///exports the graph.
867   - void
868   - to_csv()
869   - {
870   - std::ofstream ofs;
871   - ofs.open("Graph.csv", std::ofstream::out | std::ofstream::app);
872   - std::cout<<"here"<<std::endl;
873   - for(int i = 0; i < V.size(); i++)
874   - {
875   - std::cout<<"here"<<i<<std::endl;
876   - ofs << V[i].edges_to_str() << "\n";
877   - }
878   - ofs.close();
879   - }
880   - ///exports the graph.
881   - void
882   - to_gdf()
883   - {
884   - std::ofstream ofs;
885   - ofs.open("Graph.gdf", std::ofstream::out | std::ofstream::app);
886   - ofs << "nodedef>name VARCHAR\n";
887   - for(int i = 0; i < V.size(); i++)
888   - {
889   - ofs << i << "\n";
890   - }
891   - ofs << "edgedef>Nodes[1] VARCHAR, Nodes[2] VARCHAR, weight INT, length FLOAT, av_radius FLOAT \n";
892   - for(int i = 0; i < E.size(); i++)
893   - {
894   - ofs << E[i]->Nodes[1] << "," << E[i]->Nodes[2] << "," <<E[i]->n_pts()
895   - << ","<< E[i]->length() << "," << E[i]->average_radius() << "\n";
896   - }
897   - ofs.close();
898   - }
899   -
900   -};
901   -};
902   -#endif