Blame view

stim/biomodels/network.h 33.8 KB
df036f4c   David Mayerich   added the network...
1
2
3
  #ifndef STIM_NETWORK_H
  #define STIM_NETWORK_H
  
7f27eafa   David Mayerich   simplified the st...
4
  #include <stdlib.h>
9c97e126   David Mayerich   added an axis-ali...
5
  #include <assert.h>
7f27eafa   David Mayerich   simplified the st...
6
7
8
9
10
  #include <sstream>
  #include <fstream>
  #include <algorithm>
  #include <string.h>
  #include <math.h>
308a743c   David Mayerich   fixed class compa...
11
  #include <stim/math/vec3.h>
5687cf1b   Cherub P. Harder   Added #include <s...
12
  #include <stim/visualization/obj.h>
58ee2b21   David Mayerich   incorporated stim...
13
  #include <stim/visualization/cylinder.h>
fbbc07be   Jiaming Guo   replace ANN by kd...
14
  #include <stim/structures/kdtree.cuh>
41100e51   Jiaming Guo   add cpu kdtree ve...
15
  #include <stim/cuda/cudatools/timer.h>
7f27eafa   David Mayerich   simplified the st...
16
  
df036f4c   David Mayerich   added the network...
17
  
c72184d1   Jiaming Guo   add many function...
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
  #ifdef __CUDACC__
  //device gaussian function
  __device__ float gaussianFunction(float x, float std){ return expf(-x/(2*std*std));}
  
  //compute metric in parallel
  template <typename T>
  __global__ void d_metric(T* M, size_t n, T* D, float sigma){
  	size_t x = blockDim.x * blockIdx.x + threadIdx.x;
  	if(x >= n) return;
  	M[x] = 1.0f - gaussianFunction(D[x], sigma);
  }
  
  //find the corresponding edge index from array index
  __global__ void d_findedge(size_t* I, size_t n, unsigned* R, size_t* E, size_t ne){
  	size_t x = blockDim.x * blockIdx.x + threadIdx.x;
  	if(x >= n) return;
  	unsigned i = 0;
  	size_t N = 0;
  	for(unsigned e = 0; e < ne; e++){
  		N += E[e];
  		if(I[x] < N){
  			R[x] = i;
  			break;
  		}
  		i++;
  	}
  }
  #endif
  
  //hard-coded factor
  int threshold_fac = 10;
  float metric_fac = 0.6f;							// might be related to the distance field
  
df036f4c   David Mayerich   added the network...
51
  namespace stim{
7f27eafa   David Mayerich   simplified the st...
52
53
54
55
56
  /** This is the a class that interfaces with gl_spider in order to store the currently
   *   segmented network. The following data is stored and can be extracted:
   *   1)Network geometry and centerline.
   *   2)Network connectivity (a graph of nodes and edges), reconstructed using ANN library.
  */
df036f4c   David Mayerich   added the network...
57
  
df036f4c   David Mayerich   added the network...
58
59
60
  template<typename T>
  class network{
  
7f27eafa   David Mayerich   simplified the st...
61
62
  	///Each edge is a fiber with two nodes.
  	///Each node is an in index to the endpoint of the fiber in the nodes array.
58ee2b21   David Mayerich   incorporated stim...
63
  	class edge : public cylinder<T>
7f27eafa   David Mayerich   simplified the st...
64
65
66
  	{
  		public:
  		unsigned v[2];		//unique id's designating the starting and ending
7c43b7f9   pranathivemuri   return resampled ...
67
  		// default constructor
c72184d1   Jiaming Guo   add many function...
68
  		edge() : cylinder<T>() {
e8c9a82b   David Mayerich   fixed GCC errors ...
69
  			v[1] = (unsigned)(-1); v[0] = (unsigned)(-1);
2e0f052a   Pavel Govyadinov   minor changes to ...
70
  		}
7f27eafa   David Mayerich   simplified the st...
71
  		/// Constructor - creates an edge from a list of points by calling the stim::fiber constructor
e376212f   David Mayerich   added support for...
72
  
3800c27f   David Mayerich   added code to sub...
73
  		///@param p is an array of positions in space
c72184d1   Jiaming Guo   add many function...
74
  		edge(centerline p) : cylinder<T>(p){}
3800c27f   David Mayerich   added code to sub...
75
  
c72184d1   Jiaming Guo   add many function...
76
  		/// Copy constructor creates an edge from a cylinder
58ee2b21   David Mayerich   incorporated stim...
77
  		edge(stim::cylinder<T> f) : cylinder<T>(f) {}
3800c27f   David Mayerich   added code to sub...
78
79
80
  
  		/// Resamples an edge by calling the fiber resampling function
  		edge resample(T spacing){
58ee2b21   David Mayerich   incorporated stim...
81
  			edge e(cylinder<T>::resample(spacing));	//call the fiber->edge constructor
3800c27f   David Mayerich   added code to sub...
82
83
84
85
86
  			e.v[0] = v[0];					//copy the vertex data
  			e.v[1] = v[1];
  
  			return e;						//return the new edge
  		}
9c97e126   David Mayerich   added an axis-ali...
87
  
7f27eafa   David Mayerich   simplified the st...
88
  		/// Output the edge information as a string
df036f4c   David Mayerich   added the network...
89
90
  		std::string str(){
  			std::stringstream ss;
fbbc07be   Jiaming Guo   replace ANN by kd...
91
  			ss<<"("<<cylinder<T>::size()<<")\tl = "<<this->length()<<"\t"<<v[0]<<"----"<<v[1];
df036f4c   David Mayerich   added the network...
92
93
  			return ss.str();
  		}
df036f4c   David Mayerich   added the network...
94
  
c72184d1   Jiaming Guo   add many function...
95
96
97
98
99
100
101
102
103
104
105
106
107
  		std::vector<edge> split(unsigned int idx){
  		
  			std::vector< stim::cylinder<T> > C;
  			C.resize(2);
  			C =	(*this).cylinder<T>::split(idx);
  			std::vector<edge> E(C.size());
  
  			for(unsigned e = 0; e < E.size(); e++){
  				E[e] = C[e];
  			}
  			return E;
  		}
  
9c97e126   David Mayerich   added an axis-ali...
108
109
  	};
  
7f27eafa   David Mayerich   simplified the st...
110
  	///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.
308a743c   David Mayerich   fixed class compa...
111
  	class vertex : public stim::vec3<T>
7f27eafa   David Mayerich   simplified the st...
112
113
  	{
  		public:
2e0f052a   Pavel Govyadinov   minor changes to ...
114
115
116
  			//std::vector<unsigned int> edges;					//indices of edges connected to this node.
  			std::vector<unsigned int> e[2];						//indices of edges going out (e[0]) and coming in (e[1])
  			//stim::vec3<T> p;							//position of this node in physical space.
7f27eafa   David Mayerich   simplified the st...
117
118
  
  			//constructor takes a stim::vec
308a743c   David Mayerich   fixed class compa...
119
  			vertex(stim::vec3<T> p) : stim::vec3<T>(p){}
7f27eafa   David Mayerich   simplified the st...
120
121
122
123
  
  			/// Output the vertex information as a string
  			std::string	str(){
  				std::stringstream ss;
308a743c   David Mayerich   fixed class compa...
124
  				ss<<"\t(x, y, z) = "<<stim::vec3<T>::str();
7f27eafa   David Mayerich   simplified the st...
125
126
127
128
129
130
131
132
133
134
135
  
  				if(e[0].size() > 0){
  					ss<<"\t> ";
  					for(unsigned int o = 0; o < e[0].size(); o++)
  						ss<<e[0][o]<<" ";
  				}
  				if(e[1].size() > 0){
  					ss<<"\t< ";
  					for(unsigned int i = 0; i < e[1].size(); i++)
  						ss<<e[1][i]<<" ";
  				}
df036f4c   David Mayerich   added the network...
136
  
7f27eafa   David Mayerich   simplified the st...
137
  				return ss.str();
df036f4c   David Mayerich   added the network...
138
  			}
9c97e126   David Mayerich   added an axis-ali...
139
  
fe0c7539   pranathivemuri   added functions t...
140
  
df036f4c   David Mayerich   added the network...
141
142
  	};
  
9c97e126   David Mayerich   added an axis-ali...
143
  protected:
df036f4c   David Mayerich   added the network...
144
  
7f27eafa   David Mayerich   simplified the st...
145
146
  	std::vector<edge> E;       //list of edges
  	std::vector<vertex> V;	    //list of vertices.
9c97e126   David Mayerich   added an axis-ali...
147
148
  
  public:
df036f4c   David Mayerich   added the network...
149
  
2e0f052a   Pavel Govyadinov   minor changes to ...
150
151
152
153
154
155
156
157
158
159
160
161
  	///default constructor
  	network()
  	{
  		
  	}
  
  	///constructor with a file to load.
  	network(std::string fileLocation)
  	{
  		load_obj(fileLocation);
  	}
  
7f27eafa   David Mayerich   simplified the st...
162
163
164
  	///Returns the number of edges in the network.
  	unsigned int edges(){
  		return E.size();
df036f4c   David Mayerich   added the network...
165
166
  	}
  
7f27eafa   David Mayerich   simplified the st...
167
168
169
  	///Returns the number of nodes in the network.
  	unsigned int vertices(){
  		return V.size();
df036f4c   David Mayerich   added the network...
170
171
  	}
  
87d0a1d2   David Mayerich   minor changes for...
172
173
174
  	//scale the network by some constant value
  	//	I don't think these work??????
  	/*std::vector<vertex> operator*(T s){
fe0c7539   pranathivemuri   added functions t...
175
176
177
178
179
180
181
182
183
184
185
186
187
  		for (unsigned i=0; i< vertices; i ++ ){
  			V[i] = V[i] * s;
  		}
  		return V;
  	}
  
  	std::vector<vertex> operator*(vec<T> s){
  		for (unsigned i=0; i< vertices; i ++ ){
  			for (unsigned dim = 0 ; dim< 3; dim ++){
  				V[i][dim] = V[i][dim] * s[dim];
  			}
  		}
  		return V;
87d0a1d2   David Mayerich   minor changes for...
188
  	}*/
fe0c7539   pranathivemuri   added functions t...
189
  
abce2443   pranathivemuri   average branching...
190
  	// Returns an average of branching index in the network
abce2443   pranathivemuri   average branching...
191
192
193
194
195
196
197
198
199
200
  	double BranchingIndex(){
  		double B=0;
  		for(unsigned v=0; v < V.size(); v ++){
  			B += ((V[v].e[0].size()) + (V[v].e[1].size()));
  		}
  		B = B / V.size();
  		return B;
  
  	}
  
15a53122   pranathivemuri   Added other simpl...
201
  	// Returns number of branch points in thenetwork
15a53122   pranathivemuri   Added other simpl...
202
203
204
205
206
207
208
209
210
211
212
213
214
  	unsigned int BranchP(){
  		unsigned int B=0;
  		unsigned int c;
  		for(unsigned v=0; v < V.size(); v ++){
  			c = ((V[v].e[0].size()) + (V[v].e[1].size()));
  			if (c > 2){
  			B += 1;}
  		}		
  		return B;
  
  	}
  
  	// Returns number of end points (tips) in thenetwork
15a53122   pranathivemuri   Added other simpl...
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
  	unsigned int EndP(){
  		unsigned int B=0;
  		unsigned int c;
  		for(unsigned v=0; v < V.size(); v ++){
  			c = ((V[v].e[0].size()) + (V[v].e[1].size()));
  			if (c == 1){
  			B += 1;}
  		}		
  		return B;
  
  	}
  
  	//// Returns a dictionary with the key as the vertex
  	//std::map<std::vector<vertex>,unsigned int> DegreeDict(){
  	//	std::map<std::vector<vertex>,unsigned int> dd;
  	//	unsigned int c = 0;
  	//	for(unsigned v=0; v < V.size(); v ++){
  	//		c = ((V[v].e[0].size()) + (V[v].e[1].size()));
  	//		dd[V[v]] = c;
  	//	}
  	//	return dd;
  	//}
  
  	//// Return number of branching stems
  	//unsigned int Stems(){
  	//	unsigned int s = 0;
  	//	std::map<std::vector<vertex>,unsigned int> dd;
  	//	dd = DegreeDict();
  	//	//for(unsigned v=0; v < V.size(); v ++){
  	//	//	V[v].e[0].
  	//	return s;
  	//}
  
87d0a1d2   David Mayerich   minor changes for...
248
  	//Calculate Metrics---------------------------------------------------
fe0c7539   pranathivemuri   added functions t...
249
250
  	// Returns an average of fiber/edge lengths in the network
  	double Lengths(){
1b2858a2   David Mayerich   fixed some commen...
251
252
  		stim::vec<T> L;
  		double sumLength = 0;
fe0c7539   pranathivemuri   added functions t...
253
254
255
256
257
258
259
260
261
262
263
264
265
  		for(unsigned e = 0; e < E.size(); e++){				//for each edge in the network
  			L.push_back(E[e].length());						//append the edge length
  			sumLength = sumLength + E[e].length();
  		}
  		double avg = sumLength / E.size();
  		return avg;
  	}
  
  
  	// Returns an average of tortuosities in the network
  	double Tortuosities(){
  		stim::vec<T> t;
  		stim::vec<T> id1, id2;                        // starting and ending vertices of the edge
fe0c7539   pranathivemuri   added functions t...
266
267
268
269
270
  		double distance;double tortuosity;double sumTortuosity = 0;
  		for(unsigned e = 0; e < E.size(); e++){				//for each edge in the network
  			id1 = E[e][0];									//get the edge starting point
  			id2 = E[e][E[e].size() - 1];					//get the edge ending point
  			distance = (id1 - id2).len();                   //displacement between the starting and ending points
fe0c7539   pranathivemuri   added functions t...
271
272
273
274
275
276
277
278
279
  			if(distance > 0){
  				tortuosity = E[e].length()/	distance	;		// tortuoisty = edge length / edge displacement
  			}
  			else{
  				tortuosity = 0;}
  			t.push_back(tortuosity);
  			sumTortuosity += tortuosity;
  		}
  		double avg = sumTortuosity / E.size();
fe0c7539   pranathivemuri   added functions t...
280
281
282
  		return avg;
  	}
  
15a53122   pranathivemuri   Added other simpl...
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
  	// Returns average contraction of the network
  	double Contractions(){
  	stim::vec<T> t;
  	stim::vec<T> id1, id2;                        // starting and ending vertices of the edge
  	double distance;double contraction;double sumContraction = 0;
  	for(unsigned e = 0; e < E.size(); e++){				//for each edge in the network
  		id1 = E[e][0];									//get the edge starting point
  		id2 = E[e][E[e].size() - 1];					//get the edge ending point
  		distance = (id1 - id2).len();                   //displacement between the starting and ending points
  		contraction = distance / E[e].length();		// tortuoisty = edge length / edge displacement
  		t.push_back(contraction);
  		sumContraction += contraction;
  	}
  	double avg = sumContraction / E.size();
  	return avg;
  	}
fe0c7539   pranathivemuri   added functions t...
299
  
15a53122   pranathivemuri   Added other simpl...
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
  	// returns average fractal dimension of the branches of the network
  	double FractalDimensions(){
  	stim::vec<T> t;
  	stim::vec<T> id1, id2;                        // starting and ending vertices of the edge
  	double distance;double fract;double sumFractDim = 0;
  	for(unsigned e = 0; e < E.size(); e++){				//for each edge in the network
  		id1 = E[e][0];									//get the edge starting point
  		id2 = E[e][E[e].size() - 1];					//get the edge ending point
  		distance = (id1 - id2).len();                   //displacement between the starting and ending points
  		fract = std::log(distance) / std::log(E[e].length());		// tortuoisty = edge length / edge displacement
  		t.push_back(sumFractDim);
  		sumFractDim += fract;
  	}
  	double avg = sumFractDim / E.size();
  	return avg;
  	}
87d0a1d2   David Mayerich   minor changes for...
316
317
318
319
  
  	//returns a cylinder represented a given fiber (based on edge index)
  	stim::cylinder<T> get_cylinder(unsigned e){
  		return E[e];									//return the specified edge (casting it to a fiber)
3800c27f   David Mayerich   added code to sub...
320
321
  	}
  
7f27eafa   David Mayerich   simplified the st...
322
323
  	//load a network from an OBJ file
  	void load_obj(std::string filename){
e376212f   David Mayerich   added support for...
324
  
2e0f052a   Pavel Govyadinov   minor changes to ...
325
326
  		stim::obj<T> O;									//create an OBJ object
  		O.load(filename);								//load the OBJ file as an object
df036f4c   David Mayerich   added the network...
327
  
2e0f052a   Pavel Govyadinov   minor changes to ...
328
  		std::vector<unsigned> id2vert;							//this list stores the OBJ vertex ID associated with each network vertex
df036f4c   David Mayerich   added the network...
329
  
2e0f052a   Pavel Govyadinov   minor changes to ...
330
  		unsigned i[2];									//temporary, IDs associated with the first and last points in an OBJ line
df036f4c   David Mayerich   added the network...
331
  
7f27eafa   David Mayerich   simplified the st...
332
333
  		//for each line in the OBJ object
  		for(unsigned int l = 1; l <= O.numL(); l++){
df036f4c   David Mayerich   added the network...
334
  
2e0f052a   Pavel Govyadinov   minor changes to ...
335
  			std::vector< stim::vec<T> > c;						//allocate an array of points for the vessel centerline
7f27eafa   David Mayerich   simplified the st...
336
  			O.getLine(l, c);							//get the fiber centerline
df036f4c   David Mayerich   added the network...
337
  
c72184d1   Jiaming Guo   add many function...
338
  			stim::centerline<T> c3(c.size());
308a743c   David Mayerich   fixed class compa...
339
340
341
  			for(size_t j = 0; j < c.size(); j++)
  				c3[j] = c[j];
  
2e0f052a   Pavel Govyadinov   minor changes to ...
342
343
344
345
346
  	//		edge new_edge = c3;		///This is dangerous.
  			edge new_edge(c3);
  					
  			//create an edge from the given centerline
  			unsigned int I = new_edge.size();					//calculate the number of points on the centerline
df036f4c   David Mayerich   added the network...
347
  
7f27eafa   David Mayerich   simplified the st...
348
  			//get the first and last vertex IDs for the line
2e0f052a   Pavel Govyadinov   minor changes to ...
349
  			std::vector< unsigned > id;						//create an array to store the centerline point IDs
9c97e126   David Mayerich   added an axis-ali...
350
  			O.getLinei(l, id);							//get the list of point IDs for the line
7f27eafa   David Mayerich   simplified the st...
351
352
  			i[0] = id.front();							//get the OBJ ID for the first element of the line
  			i[1] = id.back();							//get the OBJ ID for the last element of the line
df036f4c   David Mayerich   added the network...
353
  
2e0f052a   Pavel Govyadinov   minor changes to ...
354
  			std::vector<unsigned>::iterator it;					//create an iterator for searching the id2vert array
7f27eafa   David Mayerich   simplified the st...
355
  			unsigned it_idx;							//create an integer for the id2vert entry index
df036f4c   David Mayerich   added the network...
356
  
7f27eafa   David Mayerich   simplified the st...
357
  			//find out if the nodes for this fiber have already been created
7f27eafa   David Mayerich   simplified the st...
358
  			it = find(id2vert.begin(), id2vert.end(), i[0]);	//look for the first node
7f27eafa   David Mayerich   simplified the st...
359
  			if(it == id2vert.end()){							//if i[0] hasn't already been used
4df9ec0e   pranathivemuri   compare networks ...
360
  				vertex new_vertex = new_edge[0];				//create a new vertex, assign it a position
2e0f052a   Pavel Govyadinov   minor changes to ...
361
362
363
364
  				new_vertex.e[0].push_back(E.size());				//add the current edge as outgoing
  				new_edge.v[0] = V.size();					//add the new edge to the edge
  				V.push_back(new_vertex);					//add the new vertex to the vertex list
  				id2vert.push_back(i[0]);					//add the ID to the ID->vertex conversion list
df036f4c   David Mayerich   added the network...
365
  			}
2e0f052a   Pavel Govyadinov   minor changes to ...
366
  			else{									//if the vertex already exists
15a53122   pranathivemuri   Added other simpl...
367
  				it_idx = std::distance(id2vert.begin(), it);
df480014   pranathivemuri   added in comments...
368
  				V[it_idx].e[0].push_back(E.size());				//add the current edge as outgoing
3800c27f   David Mayerich   added code to sub...
369
  				new_edge.v[0] = it_idx;
df036f4c   David Mayerich   added the network...
370
371
  			}
  
2e0f052a   Pavel Govyadinov   minor changes to ...
372
  			it = find(id2vert.begin(), id2vert.end(), i[1]);			//look for the second ID
2e0f052a   Pavel Govyadinov   minor changes to ...
373
374
375
  			if(it == id2vert.end()){						//if i[1] hasn't already been used
  				vertex new_vertex = new_edge[I-1];				//create a new vertex, assign it a position
  				new_vertex.e[1].push_back(E.size());				//add the current edge as incoming
b2f48bca   Pavel Govyadinov   fixed merge confl...
376
  				new_edge.v[1] = V.size();                                  	//add the new vertex to the edge
2e0f052a   Pavel Govyadinov   minor changes to ...
377
378
  				V.push_back(new_vertex);					//add the new vertex to the vertex list
  				id2vert.push_back(i[1]);					//add the ID to the ID->vertex conversion list
df036f4c   David Mayerich   added the network...
379
  			}
2e0f052a   Pavel Govyadinov   minor changes to ...
380
  			else{									//if the vertex already exists
15a53122   pranathivemuri   Added other simpl...
381
  				it_idx = std::distance(id2vert.begin(), it);
df480014   pranathivemuri   added in comments...
382
  				V[it_idx].e[1].push_back(E.size());				//add the current edge as incoming
3800c27f   David Mayerich   added code to sub...
383
  				new_edge.v[1] = it_idx;
df036f4c   David Mayerich   added the network...
384
  			}
9c97e126   David Mayerich   added an axis-ali...
385
  
2e0f052a   Pavel Govyadinov   minor changes to ...
386
  			E.push_back(new_edge);							//push the edge to the list
df036f4c   David Mayerich   added the network...
387
  
df036f4c   David Mayerich   added the network...
388
  		}
e376212f   David Mayerich   added support for...
389
390
  	}
  
7f27eafa   David Mayerich   simplified the st...
391
392
  	/// Output the network as a string
  	std::string str(){
e376212f   David Mayerich   added support for...
393
  
7f27eafa   David Mayerich   simplified the st...
394
395
396
397
  		std::stringstream ss;
  		ss<<"Nodes ("<<V.size()<<")--------"<<std::endl;
  		for(unsigned int v = 0; v < V.size(); v++){
  			ss<<"\t"<<v<<V[v].str()<<std::endl;
e376212f   David Mayerich   added support for...
398
399
  		}
  
7f27eafa   David Mayerich   simplified the st...
400
401
402
  		ss<<"Edges ("<<E.size()<<")--------"<<std::endl;
  		for(unsigned e = 0; e < E.size(); e++){
  			ss<<"\t"<<e<<E[e].str()<<std::endl;
df036f4c   David Mayerich   added the network...
403
404
  		}
  
7f27eafa   David Mayerich   simplified the st...
405
  		return ss.str();
df036f4c   David Mayerich   added the network...
406
  	}
3800c27f   David Mayerich   added code to sub...
407
408
  
  	/// This function resamples all fibers in a network given a desired minimum spacing
9c97e126   David Mayerich   added an axis-ali...
409
  	/// @param spacing is the minimum distance between two points on the network
3800c27f   David Mayerich   added code to sub...
410
  	stim::network<T> resample(T spacing){
2e0f052a   Pavel Govyadinov   minor changes to ...
411
412
  		stim::network<T> n;								//create a new network that will be an exact copy, with resampled fibers
  		n.V = V;									//copy all vertices
3800c27f   David Mayerich   added code to sub...
413
  
2e0f052a   Pavel Govyadinov   minor changes to ...
414
  		n.E.resize(edges());								//allocate space for the edge list
3800c27f   David Mayerich   added code to sub...
415
416
  
  		//copy all fibers, resampling them in the process
2e0f052a   Pavel Govyadinov   minor changes to ...
417
418
  		for(unsigned e = 0; e < edges(); e++){						//for each edge in the edge list
  			n.E[e] = E[e].resample(spacing);					//resample the edge and copy it to the new network
3800c27f   David Mayerich   added code to sub...
419
420
  		}
  
2e0f052a   Pavel Govyadinov   minor changes to ...
421
  		return n;							              	//return the resampled network
4df9ec0e   pranathivemuri   compare networks ...
422
423
  	}
  
c72184d1   Jiaming Guo   add many function...
424
425
  	// host gaussian function
  	__host__ float gaussianFunction(float x, float std = 25){ return exp(-x/(2*std*std));}  // std default value is 25
4df9ec0e   pranathivemuri   compare networks ...
426
  
9c97e126   David Mayerich   added an axis-ali...
427
  	/// Calculate the total number of points on all edges.
797f8ab0   David Mayerich   cleaned up the st...
428
429
430
431
  	unsigned total_points(){
  		unsigned n = 0;
  		for(unsigned e = 0; e < E.size(); e++)
  			n += E[e].size();
4df9ec0e   pranathivemuri   compare networks ...
432
  		return n;
797f8ab0   David Mayerich   cleaned up the st...
433
  	}
4df9ec0e   pranathivemuri   compare networks ...
434
  
fcd2eb7c   David Mayerich   added support for...
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
  	//Copy the point cloud representing the centerline for the network into an array
  	void centerline_cloud(T* dst) {
  		size_t p;										//stores the current edge point
  		size_t P;										//stores the number of points in an edge
  		size_t i = 0;									//index into the output array of points
  		for (size_t e = 0; e < E.size(); e++) {			//for each edge in the network
  			P = E[e].size();							//get the number of points in this edge
  			for (p = 0; p < P; p++) {
  				dst[i * 3 + 0] = E[e][p][0];		
  				dst[i * 3 + 1] = E[e][p][1];
  				dst[i * 3 + 2] = E[e][p][2];
  				i++;
  			}
  		}
  	}
  
fbbc07be   Jiaming Guo   replace ANN by kd...
451
452
      // convert vec3 to array
  	void stim2array(float *a, stim::vec3<T> b){
797f8ab0   David Mayerich   cleaned up the st...
453
454
455
456
457
  		a[0] = b[0];
  		a[1] = b[1];
  		a[2] = b[2];
  	}
  
c72184d1   Jiaming Guo   add many function...
458
459
460
461
462
463
464
465
466
467
  	// convert vec3 to array in bunch
  	void edge2array(T* a, std::vector< typename stim::vec3<T> > b){
  		size_t n = b.size();
  		for(size_t i = 0; i < n; i++){
  			a[i * 3 + 0] = b[i][0];
  			a[i * 3 + 1] = b[i][1];
  			a[i * 3 + 2] = b[i][2];	 
  		}
  	}
  
9c97e126   David Mayerich   added an axis-ali...
468
469
470
471
  	/// Calculate the average magnitude across the entire network.
  	/// @param m is the magnitude value to use. The default is 0 (usually radius).
  	T average(unsigned m = 0){
  
2e0f052a   Pavel Govyadinov   minor changes to ...
472
473
474
  		T M, L;										//allocate space for the total magnitude and length
  		M = L = 0;									//initialize both the initial magnitude and length to zero
  		for(unsigned e = 0; e < E.size(); e++){						//for each edge in the network
c72184d1   Jiaming Guo   add many function...
475
  			M += E[e].integrate();							//get the integrated magnitude
2e0f052a   Pavel Govyadinov   minor changes to ...
476
  			L += E[e].length();							//get the edge length
9c97e126   David Mayerich   added an axis-ali...
477
  		}
4df9ec0e   pranathivemuri   compare networks ...
478
  
2e0f052a   Pavel Govyadinov   minor changes to ...
479
  		return M / L;									//return the average magnitude
9c97e126   David Mayerich   added an axis-ali...
480
481
482
  	}
  
  	/// This function compares two networks and returns the percentage of the current network that is missing from A.
9c97e126   David Mayerich   added an axis-ali...
483
484
  	/// @param A is the network to compare to - the field is generated for A
  	/// @param sigma is the user-defined tolerance value - smaller values provide a stricter comparison
c72184d1   Jiaming Guo   add many function...
485
  	/// @param device is the device that user want to use
41100e51   Jiaming Guo   add cpu kdtree ve...
486
  	stim::network<T> compare(stim::network<T> A, float sigma, int device){
9c97e126   David Mayerich   added an axis-ali...
487
  
fcd2eb7c   David Mayerich   added support for...
488
489
  		stim::network<T> R;										//generate a network storing the result of the comparison
  		R = (*this);											//initialize the result with the current network
9c97e126   David Mayerich   added an axis-ali...
490
  
fcd2eb7c   David Mayerich   added support for...
491
  		T *c;						                 			// centerline (array of double pointers) - points on kdtree must be double
41100e51   Jiaming Guo   add cpu kdtree ve...
492
  		size_t n_data = A.total_points();          				// set the number of points
fcd2eb7c   David Mayerich   added support for...
493
  		c = (T*) malloc(sizeof(T) * n_data * 3);				//allocate an array to store all points in the data set				
9c97e126   David Mayerich   added an axis-ali...
494
  
797f8ab0   David Mayerich   cleaned up the st...
495
  		unsigned t = 0;
fcd2eb7c   David Mayerich   added support for...
496
497
  		for(unsigned e = 0; e < A.E.size(); e++){				//for each edge in the network
  			for(unsigned p = 0; p < A.E[e].size(); p++){		//for each point in the edge
797f8ab0   David Mayerich   cleaned up the st...
498
499
  				for(unsigned d = 0; d < 3; d++){				//for each coordinate
  
fcd2eb7c   David Mayerich   added support for...
500
  					c[t * 3 + d] = A.E[e][p][d];				//copy the point into the array c
4df9ec0e   pranathivemuri   compare networks ...
501
  				}
797f8ab0   David Mayerich   cleaned up the st...
502
  				t++;
4df9ec0e   pranathivemuri   compare networks ...
503
  			}
797f8ab0   David Mayerich   cleaned up the st...
504
505
  		}
  
41100e51   Jiaming Guo   add cpu kdtree ve...
506
  		//generate a KD-tree for network A
fcd2eb7c   David Mayerich   added support for...
507
  		size_t MaxTreeLevels = 3;								// max tree level
41100e51   Jiaming Guo   add cpu kdtree ve...
508
509
510
  		
  #ifdef __CUDACC__
  		cudaSetDevice(device);
c72184d1   Jiaming Guo   add many function...
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
  		stim::kdtree<T, 3> kdt;								// initialize a pointer to a kd tree
  
  		kdt.create(c, n_data, MaxTreeLevels);				// build a KD tree
  
  		for(unsigned e = 0; e < R.E.size(); e++){					//for each edge in A
  			R.E[e].add_mag(0);							//add a new magnitude for the metric
  			size_t errormag_id = R.E[e].nmags() - 1;		//get the id for the new magnitude
  			
  			size_t n = R.E[e].size();						// the number of points in current edge
  			T* queryPt = new T[n];
  			T* m1 = new T[n];
  			T* dists = new T[n];
  			size_t* nnIdx = new size_t[n];
  
  			T* d_dists;										
  			T* d_m1;										
  			cudaMalloc((void**)&d_dists, n * sizeof(T));
  			cudaMalloc((void**)&d_m1, n * sizeof(T));
  
  			edge2array(queryPt, R.E[e]);
  			kdt.search(queryPt, n, nnIdx, dists);		
  
  			cudaMemcpy(d_dists, dists, n * sizeof(T), cudaMemcpyHostToDevice);					// copy dists from host to device
  
  			// configuration parameters
  			size_t threads = (1024>n)?n:1024;
  			size_t blocks = n/threads + (n%threads)?1:0;
  
  			d_metric<<<blocks, threads>>>(d_m1, n, d_dists, sigma);					//calculate the metric value based on the distance
  
  			cudaMemcpy(m1, d_m1, n * sizeof(T), cudaMemcpyDeviceToHost);
  
  			for(unsigned p = 0; p < n; p++){
  				R.E[e].set_mag(errormag_id, p, m1[p]);
  			}
  
  			//d_set_mag<<<blocks, threads>>>(R.E[e].M, errormag_id, n, m1);
  		}
  
  #else
  		stim::kdtree<T, 3> kdt;
  		kdt.create(c, n_data, MaxTreeLevels);
  	
  		for(unsigned e = 0; e < R.E.size(); e++){			//for each edge in A
  			R.E[e].add_mag(0);								//add a new magnitude for the metric
  			size_t errormag_id = R.E[e].nmags() - 1;
  
  			size_t n = R.E[e].size();						// the number of points in current edge
  			T* query = new T[n];
  			T* m1 = new T[n];
  			T* dists = new T[n];
  			size* nnIdx = new size_t[n];
  
  			edge2array(queryPt, R.E[e]);
  
  			kdt.cpu_search(queryPt, n, nnIdx, dists);			//find the distance between A and the current network
  
  			for(unsigned p = 0; p < R.E[e].size(); p++){
  				m1[p] = 1.0f - gaussianFunction((T)dists[p], sigma);		//calculate the metric value based on the distance
  				R.E[e].set_mag(errormag_id, p, m1[p]);					//set the error for the second point in the segment
  			}
  		}
  #endif
  		return R;		//return the resulting network
  	}
  
  	/// This function compares two networks and split the current one according to the nearest neighbor of each point in each edge
  	/// @param A is the network to split
  	/// @param B is the corresponding mapping network
  	/// @param sigma is the user-defined tolerance value - smaller values provide a stricter comparison
  	/// @param device is the device that user want to use
  	void split(stim::network<T> A, stim::network<T> B, float sigma, int device){
  
  		T *c;						                 	
  		size_t n_data = B.total_points();          				
  		c = (T*) malloc(sizeof(T) * n_data * 3); 				
  
  		size_t NB = B.E.size();								// the number of edges in B
  		unsigned t = 0;
  		for(unsigned e = 0; e < NB; e++){					// for every edge in B			
  			for(unsigned p = 0; p < B.E[e].size(); p++){	// for every points in B.E[e]
  				for(unsigned d = 0; d < 3; d++){				
  
  					c[t * 3 + d] = B.E[e][p][d];			// convert to array
  				}
  				t++;
  			}
  		}
  		size_t MaxTreeLevels = 3;							// max tree level
  
  #ifdef __CUDACC__
  		cudaSetDevice(device);
  		stim::kdtree<T, 3> kdt;								// initialize a pointer to a kd tree
41100e51   Jiaming Guo   add cpu kdtree ve...
604
  	
9c97e126   David Mayerich   added an axis-ali...
605
  		//compare each point in the current network to the field produced by A
41100e51   Jiaming Guo   add cpu kdtree ve...
606
  		kdt.create(c, n_data, MaxTreeLevels);				// build a KD tree
797f8ab0   David Mayerich   cleaned up the st...
607
  
c72184d1   Jiaming Guo   add many function...
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
  		std::vector<std::vector<unsigned>> relation;		// the relationship between GT and T corresponding to NN
  		relation.resize(A.E.size());										
  
  		for(unsigned e = 0; e < A.E.size(); e++){			//for each edge in A
  			A.E[e].add_mag(0);								//add a new magnitude for the metric
  			size_t errormag_id = A.E[e].nmags() - 1;
  			
  			size_t n = A.E[e].size();						// the number of edges in A
  
  			T* queryPt = new T[n];							// set of all the points in current edge
  			T* m1 = new T[n];								// array of metrics for every point in current edge
  			T* dists = new T[n];							// store the distances for every point in current edge
  			size_t* nnIdx = new size_t[n];					// store the indices for every point in current edge
  			
  			// define pointers in device
  			T* d_dists;														
  			T* d_m1;
  			size_t* d_nnIdx;
  
  			// allocate memory for defined pointers
  			cudaMalloc((void**)&d_dists, n * sizeof(T));
  			cudaMalloc((void**)&d_m1, n * sizeof(T));
  			cudaMalloc((void**)&d_nnIdx, n * sizeof(size_t));
  
  			edge2array(queryPt, A.E[e]);						// convert edge to array
  			kdt.search(queryPt, n, nnIdx, dists);				// search the tree to find the NN for every point in current edge
  
  			cudaMemcpy(d_dists, dists, n * sizeof(T), cudaMemcpyHostToDevice);					// copy dists from host to device
  			cudaMemcpy(d_nnIdx, nnIdx, n * sizeof(size_t), cudaMemcpyHostToDevice);				// copy Idx from host to device
  
  			// configuration parameters
  			size_t threads = (1024>n)?n:1024;													// test to see whether the number of point in current edge is more than 1024
  			size_t blocks = n/threads + (n%threads)?1:0;
  
  			d_metric<<<blocks, threads>>>(d_m1, n, d_dists, sigma);								// calculate the metrics in parallel
  
  			cudaMemcpy(m1, d_m1, n * sizeof(T), cudaMemcpyDeviceToHost);
  
  			for(unsigned p = 0; p < n; p++){
  				A.E[e].set_mag(errormag_id, p, m1[p]);											// set the error(radius) value to every point in current edge
  			}
  
  			relation[e].resize(n);																// resize every edge relation size
  
  			unsigned* d_relation;
  			cudaMalloc((void**)&d_relation, n * sizeof(unsigned));								// allocate memory
  
  			std::vector<size_t> edge_point_num(NB);												// %th element is the number of points that %th edge has
  			for(unsigned ee = 0; ee < NB; ee++)
  				edge_point_num[ee] = B.E[ee].size();
  
  			size_t* d_edge_point_num;
  			cudaMalloc((void**)&d_edge_point_num, NB * sizeof(size_t));
  			cudaMemcpy(d_edge_point_num, &edge_point_num[0], NB * sizeof(size_t), cudaMemcpyHostToDevice);
  
  			d_findedge<<<blocks, threads>>>(d_nnIdx, n, d_relation, d_edge_point_num, NB);			// find the edge corresponding to the array index in parallel
  
  			cudaMemcpy(&relation[e][0], d_relation, n * sizeof(unsigned), cudaMemcpyDeviceToHost);	//copy relationship from device to host
  		}
  #else
  		stim::kdtree<T, 3> kdt;
  		kdt.create(c, n_data, MaxTreeLevels);
  	
  		std::vector<std::vector<unsigned>> relation;		// the mapping relationship between two networks
  		relation.resize(A.E.size());										
  		for(unsigned i = 0; i < A.E.size(); i++)
  			relation[i].resize(A.E[i].size());
  
  		std::vector<size_t> edge_point_num(NB);				//%th element is the number of points that %th edge has
  		for(unsigned ee = 0; ee < NB; ee++)
  			edge_point_num[ee] = B.E[ee].size();
  
41100e51   Jiaming Guo   add cpu kdtree ve...
680
  		T* queryPt = new T[3];
c72184d1   Jiaming Guo   add many function...
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
  		for(unsigned e = 0; e < A.E.size(); e++){			//for each edge in A
  			A.E[e].add_mag(0);								//add a new magnitude for the metric
  			size_t errormag_id = A.E[e].nmags() - 1;
  
  			size_t n = A.E[e].size();						// the number of edges in A
  
  			T* queryPt = new T[n];
  			T* m1 = new T[n];
  			T* dists = new T[n];							//store the distances
  			size_t* nnIdx = new size_t[n];					//store the indices
  			
  			edge2array(queryPt, A.E[e]);
  			kdt.search(queryPt, n, nnIdx, dists);		
  
  			for(unsigned p = 0; p < A.E[e].size(); p++){
  				m1[p] = 1.0f - gaussianFunction((T)dists[p], sigma);	//calculate the metric value based on the distance
  				A.E[e].set_mag(errormag_id, p, m1[p]);					//set the error for the second point in the segment
  				
  				unsigned id = 0;																	//mapping edge's idx
  				size_t num = 0;																		//total number of points before #th edge
  				for(unsigned i = 0; i < NB; i++){
  					num += B.E[i].size();
  					if(nnIdx[p] < num){																//find the edge it belongs to
  						relation[e][p] = id;
  						break;
  					}
  					id++;																			//current edge won't be the one, move to next edge
  				}
  			}
  		}
  #endif
  		E = A.E;
  		V = A.V;
  
  		unsigned int id = 0;									// split value								
  		for(unsigned e = 0; e < E.size(); e++){					// for every edge
  			for(unsigned p = 0; p < E[e].size() - 1; p++){		// for every point in each edge
  				if(relation[e][p] != relation[e][p + 1]){		// find the nearest edge changing point
  					id = p + 1;									// if there is no change in NN
  					if(id < E[e].size()/threshold_fac || (E[e].size() - id) < E[e].size()/threshold_fac)				// set tolerance to 10, tentatively--should find out the mathematical threshold!
  						id = E[e].size() - 1;																			// extreme situation is not acceptable
  					else
  						break;
  				}
  				if(p == E[e].size() - 2)																// if there is no splitting index, set the id to the last point index of current edge
  					id = p + 1;
  			}
  			unsigned errormag_id = E[e].nmags() - 1;
  			T G = 0;																					// test to see whether it has its nearest neighbor
  			for(unsigned i = 0; i < E[e].size(); i++)
  				G += E[e].m(i, errormag_id);															// won't split special edges
  			if(G / E[e].size() > metric_fac)															// should based on the color map
  				id = E[e].size() - 1;																	// set split idx to outgoing direction vertex
  
  			std::vector<edge> tmpe;
  			tmpe.resize(2);
  			tmpe = E[e].split(id);
  			vertex tmpv = stim::vec3<T>(-1, -1, 0);														// store the split point as vertex
  			if(tmpe.size() == 2){
  				relation.resize(relation.size() + 1);
  				for(unsigned d = id; d < E[e].size(); d++)
  					relation[relation.size() - 1].push_back(relation[e][d]);
  				tmpe[0].v[0] = E[e].v[0];																// begining vertex of first half edge -> original begining vertex
  				tmpe[1].v[1] = E[e].v[1];																// ending vertex of second half edge -> original ending vertex
  				tmpv = E[e][id];
  				V.push_back(tmpv);
  				tmpe[0].v[1] = (unsigned)V.size() - 1;													// ending vertex of first half edge -> new vertex
  				tmpe[1].v[0] = (unsigned)V.size() - 1;													// begining vertex of second half edge -> new vertex
  				edge tmp(E[e]);
  				E[e] = tmpe[0];																			// replace original edge by first half edge
  				E.push_back(tmpe[1]);																	// push second half edge to the last
  				V[V.size() - 1].e[1].push_back(e);														// push first half edge to the incoming of new vertex
  				V[V.size() - 1].e[0].push_back((unsigned)E.size() - 1);									// push second half edge to the outgoing of new vertex
  				for(unsigned i = 0; i < V[tmp.v[1]].e[1].size(); i++)									// find the incoming edge of original ending vertex
  					if(V[tmp.v[1]].e[1][i] == e)
  						V[tmp.v[1]].e[1][i] = (unsigned)V.size() - 1;
  			}
  		}
  	}
797f8ab0   David Mayerich   cleaned up the st...
760
  
c72184d1   Jiaming Guo   add many function...
761
762
763
764
765
766
767
  	/// This function compares two splitted networks and yields a mapping relationship between them according to NN
  	/// @param B is the network that the current network is going to map to
  	/// @param C is the mapping relationship: C[e1] = _e1 means e1 edge in current network is mapping the _e1 edge in B
  	/// @param device is the device that user want to use
  	void mapping(stim::network<T> B, std::vector<unsigned> &C, int device){
  		stim::network<T> A;								//generate a network storing the result of the comparison
  		A = (*this);
797f8ab0   David Mayerich   cleaned up the st...
768
  
c72184d1   Jiaming Guo   add many function...
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
  		size_t n = A.E.size();							// the number of edges in A
  		size_t NB = B.E.size();							// the number of edges in B
  
  		C.resize(A.E.size());	
  
  		T *c;						                 	// centerline (array of double pointers) - points on kdtree must be double
  		size_t n_data = B.total_points();          		// set the number of points
  		c = (T*) malloc(sizeof(T) * n_data * 3); 				
  
  		unsigned t = 0;
  		for(unsigned e = 0; e < NB; e++){					// for each edge in the network
  			for(unsigned p = 0; p < B.E[e].size(); p++){	// for each point in the edge
  				for(unsigned d = 0; d < 3; d++){			// for each coordinate
  
  					c[t * 3 + d] = B.E[e][p][d];
  				}
  				t++;
  			}
  		}
fbbc07be   Jiaming Guo   replace ANN by kd...
788
  
c72184d1   Jiaming Guo   add many function...
789
790
791
792
793
794
795
796
797
  		//generate a KD-tree for network A
  		//float metric = 0.0;                               		// initialize metric to be returned after comparing the network
  		size_t MaxTreeLevels = 3;									// max tree level
  		
  #ifdef __CUDACC__
  		cudaSetDevice(device);
  		stim::kdtree<T, 3> kdt;								// initialize a pointer to a kd tree
  	
  		kdt.create(c, n_data, MaxTreeLevels);				// build a KD tree
797f8ab0   David Mayerich   cleaned up the st...
798
  
c72184d1   Jiaming Guo   add many function...
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
  		for(unsigned e = 0; e < n; e++){					//for each edge in A
  			size_t errormag_id = A.E[e].nmags() - 1;		//get the id for the new magnitude
  			
  			//pre-judge to get rid of impossibly mapping edges
  			T M = 0;
  			for(unsigned p = 0; p < A.E[e].size(); p++)
  				M += A.E[e].m(p, errormag_id);
  			M = M / A.E[e].size();
  			if(M > metric_fac)
  				C[e] = (unsigned)-1;						//set the nearest edge of impossibly mapping edges to maximum of unsigned
  			else{
  				T* queryPt = new T[1];
  				T* dists = new T[1];
  				size_t* nnIdx = new size_t[1];
  
  				stim2array(queryPt, A.E[e][A.E[e].size()/2]);
  				kdt.search(queryPt, 1, nnIdx, dists);
  				
  				unsigned id = 0;							//mapping edge's idx
  				size_t num = 0;								//total number of points before #th edge
  				for(unsigned i = 0; i < NB; i++){
  					num += B.E[i].size();
  					if(nnIdx[0] < num){
  						C[e] = id;
  						break;
  					}
  					id++;
  				}
797f8ab0   David Mayerich   cleaned up the st...
827
828
  			}
  		}
41100e51   Jiaming Guo   add cpu kdtree ve...
829
  #else
c72184d1   Jiaming Guo   add many function...
830
  		stim::kdtree<T, 3> kdt;
41100e51   Jiaming Guo   add cpu kdtree ve...
831
832
  		kdt.create(c, n_data, MaxTreeLevels);
  		T *dists = new T[1];								// near neighbor distances
fcd2eb7c   David Mayerich   added support for...
833
  		size_t *nnIdx = new size_t[1];						// near neighbor indices // allocate near neigh indices
797f8ab0   David Mayerich   cleaned up the st...
834
  
41100e51   Jiaming Guo   add cpu kdtree ve...
835
  		stim::vec3<T> p0, p1;
41100e51   Jiaming Guo   add cpu kdtree ve...
836
  		T* queryPt = new T[3];
41100e51   Jiaming Guo   add cpu kdtree ve...
837
  
c72184d1   Jiaming Guo   add many function...
838
839
840
841
842
843
844
845
846
847
  		for(unsigned e = 0; e < R.E.size(); e++){			// for each edge in A
  			T M;											// the sum of metrics of current edge
  			unsigned errormag_id = R.E[e].nmags() - 1;
  			for(unsigned p = 0; p < R.E[e].size(); p++)
  				M += A.E[e].m(p, errormag_id);
  			M = M / A.E[e].size();
  			if(M > metric_fac)								// if the sum is larger than the metric_fac, we assume that it doesn't has corresponding edge in B
  				C[e] = (unsigned)-1;
  			else{											// if it should have corresponding edge in B, then...
  				p1 = R.E[e][R.E[e].size()/2];							
41100e51   Jiaming Guo   add cpu kdtree ve...
848
  				stim2array(queryPt, p1);
c72184d1   Jiaming Guo   add many function...
849
850
851
852
853
854
855
856
857
858
859
860
  				kdt.cpu_search(queryPt, 1, nnIdx, dists);	// search the tree		
  				
  				unsigned id = 0;							//mapping edge's idx
  				size_t num = 0;								//total number of points before #th edge
  				for(unsigned i = 0; i < NB; i++){
  					num += B.E[i].size();
  					if(nnIdx[0] < num){
  						C[e] = id;
  						break;
  					}
  					id++;
  				}
41100e51   Jiaming Guo   add cpu kdtree ve...
861
862
863
  			}
  		}
  #endif
3800c27f   David Mayerich   added code to sub...
864
  	}
9c97e126   David Mayerich   added an axis-ali...
865
  
b3a38641   David Mayerich   added the ability...
866
867
868
869
  	/// Returns the number of magnitude values stored in each edge. This should be uniform across the network.
  	unsigned nmags(){
  		return E[0].nmags();
  	}
10b95d6e   pranathivemuri   load text file to...
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
  	// split a string in text by the character sep
  	stim::vec<T> split(std::string &text, char sep) 
  	{
  		stim::vec<T> tokens;
  		std::size_t start = 0, end = 0;
  		while ((end = text.find(sep, start)) != std::string::npos) {
  		tokens.push_back(atof(text.substr(start, end - start).c_str()));
  		start = end + 1;
  		}
  		tokens.push_back(atof(text.substr(start).c_str()));
  		return tokens;
  	}
  	// load a network in text file to a network class
  	void load_txt(std::string filename)
  	{
  		std::vector <std::string> file_contents;
e8c9a82b   David Mayerich   fixed GCC errors ...
886
  		std::ifstream file(filename.c_str());
10b95d6e   pranathivemuri   load text file to...
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
  		std::string line;
  		std::vector<unsigned> id2vert;	//this list stores the vertex ID associated with each network vertex
  		//for each line in the text file, store them as strings in file_contents
  		while (std::getline(file, line))
  		{
  			std::stringstream ss(line);
  			file_contents.push_back(ss.str());
  		}
  		unsigned int numEdges = atoi(file_contents[0].c_str()); //number of edges in the network
  		unsigned int I = atoi(file_contents[1].c_str()) ;		//calculate the number of points3d on the first edge
  		unsigned int count = 1; unsigned int k = 2; // count is global counter through the file contents, k is for the vertices on the edges
  		// for each edge in the network.
  		for (unsigned int i = 0; i < numEdges; i ++ )
  		{
  			// pre allocate a position vector p with number of points3d on the edge p
  			std::vector< stim::vec<T> > p(0, I);
  			// for each point on the nth edge
  		  for (unsigned int j = k; j < I + k; j++)
  		  {
  			 // split the points3d of floats with separator space and form a float3 position vector out of them
  			  p.push_back(split(file_contents[j], ' '));
  		  }
  			count += p.size() + 1; // increment count to point at the next edge in the network
  			I = atoi(file_contents[count].c_str()); // read in the points3d at the next edge and convert it to an integer
  			k = count + 1;
  			edge new_edge = p; // create an edge with a vector of points3d  on the edge
  			E.push_back(new_edge); // push the edge into the network
  		}
  		unsigned int numVertices = atoi(file_contents[count].c_str()); // this line in the text file gives the number of distinct vertices
  		count = count + 1; // this line of text file gives the first verrtex
  		// push each vertex into V
  		for (unsigned int i = 0; i < numVertices; i ++)
  		{
  			vertex new_vertex = split(file_contents[count], ' ');
  			V.push_back(new_vertex);
  			count += atoi(file_contents[count + 1].c_str()) + 2; // Skip number of edge ids + 2 to point to the next vertex
  		}
  	} // end load_txt function
b3a38641   David Mayerich   added the ability...
925
  
abce2443   pranathivemuri   average branching...
926
927
928
929
930
931
932
933
934
  	// strTxt returns a string of edges
  	std::string
  	strTxt(std::vector< stim::vec<T> > p)
  	{
  		std::stringstream ss;
  		std::stringstream oss;
  		for(unsigned int i = 0; i < p.size(); i++){
  			ss.str(std::string());
  			for(unsigned int d = 0; d < 3; d++){
9f5c0d4a   Pavel Govyadinov   edits to make net...
935
  				ss<<p[i][d];
abce2443   pranathivemuri   average branching...
936
  			}
8b7d6f9a   David Mayerich   light modificatio...
937
  			ss << "\n";
abce2443   pranathivemuri   average branching...
938
939
940
941
942
943
944
945
946
947
948
949
950
  		}
  		return ss.str();
  	}
  	// removes specified character from string
  	void removeCharsFromString(std::string &str, char* charsToRemove ) {
  	   for ( unsigned int i = 0; i < strlen(charsToRemove); ++i ) {
  		  str.erase( remove(str.begin(), str.end(), charsToRemove[i]), str.end() );
  	   }
  	}
  	//exports network to txt file
  	void
  	to_txt(std::string filename)
  	{
e8c9a82b   David Mayerich   fixed GCC errors ...
951
952
  		std::ofstream ofs(filename.c_str(), std::ofstream::out | std::ofstream::app);
  		//int num;
abce2443   pranathivemuri   average branching...
953
954
955
956
957
958
959
960
961
962
963
964
  		ofs << (E.size()).str() << "\n";
  		for(unsigned int i = 0; i < E.size(); i++)
  		{
  			 std::string str;
  			 ofs << (E[i].size()).str() << "\n";
  			 str = E[i].strTxt();
               ofs << str << "\n";
  		} 
  		for(int i = 0; i < V.size(); i++)
  		{
  			std::string str;
  			str = V[i].str();
e8c9a82b   David Mayerich   fixed GCC errors ...
965
966
  			char temp[4] = "[],";
  			removeCharsFromString(str, temp);
abce2443   pranathivemuri   average branching...
967
968
969
970
  			ofs << str << "\n";
  		}
  		ofs.close();
  	}
7f27eafa   David Mayerich   simplified the st...
971
972
  };		//end stim::network class
  };		//end stim namespace
c72184d1   Jiaming Guo   add many function...
973
  #endif