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>
|
4df9ec0e
pranathivemuri
compare networks ...
|
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>
|
c61fd046
David Mayerich
fixed vessel mode...
|
11
|
#include <stim/math/vector.h>
|
5687cf1b
Cherub P. Harder
Added #include <s...
|
12
|
#include <stim/visualization/obj.h>
|
7f27eafa
David Mayerich
simplified the st...
|
13
|
#include <stim/biomodels/fiber.h>
|
df036f4c
David Mayerich
added the network...
|
14
|
#include <ANN/ANN.h>
|
7f27eafa
David Mayerich
simplified the st...
|
15
16
|
#include <boost/tuple/tuple.hpp>
|
df036f4c
David Mayerich
added the network...
|
17
18
|
namespace stim{
|
7f27eafa
David Mayerich
simplified the st...
|
19
20
21
22
23
|
/** 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...
|
24
|
|
df036f4c
David Mayerich
added the network...
|
25
26
27
28
|
template<typename T>
class network{
|
7f27eafa
David Mayerich
simplified the st...
|
29
30
31
32
33
34
|
///Each edge is a fiber with two nodes.
///Each node is an in index to the endpoint of the fiber in the nodes array.
class edge : public fiber<T>
{
public:
unsigned v[2]; //unique id's designating the starting and ending
|
7c43b7f9
pranathivemuri
return resampled ...
|
35
|
// default constructor
|
4df9ec0e
pranathivemuri
compare networks ...
|
36
|
edge() : fiber<T>(){v[1] = -1; v[0] = -1;}
|
7f27eafa
David Mayerich
simplified the st...
|
37
|
/// Constructor - creates an edge from a list of points by calling the stim::fiber constructor
|
e376212f
David Mayerich
added support for...
|
38
|
|
3800c27f
David Mayerich
added code to sub...
|
39
|
///@param p is an array of positions in space
|
7f27eafa
David Mayerich
simplified the st...
|
40
|
edge(std::vector< stim::vec<T> > p) : fiber<T>(p){}
|
3800c27f
David Mayerich
added code to sub...
|
41
42
43
44
45
46
|
/// Copy constructor creates an edge from a fiber
edge(stim::fiber<T> f) : fiber<T>(f) {}
/// Resamples an edge by calling the fiber resampling function
edge resample(T spacing){
|
797f8ab0
David Mayerich
cleaned up the st...
|
47
|
edge e(fiber<T>::resample(spacing)); //call the fiber->edge constructor
|
3800c27f
David Mayerich
added code to sub...
|
48
49
50
51
52
|
e.v[0] = v[0]; //copy the vertex data
e.v[1] = v[1];
return e; //return the new edge
}
|
7f27eafa
David Mayerich
simplified the st...
|
53
54
|
/// Output the edge information as a string
|
df036f4c
David Mayerich
added the network...
|
55
56
|
std::string str(){
std::stringstream ss;
|
7f27eafa
David Mayerich
simplified the st...
|
57
|
ss<<"("<<fiber<T>::N<<")\tl = "<<length()<<"\t"<<v[0]<<"----"<<v[1];
|
df036f4c
David Mayerich
added the network...
|
58
59
|
return ss.str();
}
|
df036f4c
David Mayerich
added the network...
|
60
|
|
7f27eafa
David Mayerich
simplified the st...
|
61
62
63
64
65
66
67
68
|
};
///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.
class vertex : public stim::vec<T>
{
public:
//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])
|
df480014
pranathivemuri
added in comments...
|
69
|
//stim::vec<T> p; //position of this node in physical space.
|
7f27eafa
David Mayerich
simplified the st...
|
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
|
//constructor takes a stim::vec
vertex(stim::vec<T> p) : stim::vec<T>(p){}
/// Output the vertex information as a string
std::string str(){
std::stringstream ss;
ss<<"\t(x, y, z) = "<<stim::vec<T>::str();
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...
|
89
|
|
7f27eafa
David Mayerich
simplified the st...
|
90
|
return ss.str();
|
df036f4c
David Mayerich
added the network...
|
91
|
}
|
7f27eafa
David Mayerich
simplified the st...
|
92
|
|
df036f4c
David Mayerich
added the network...
|
93
94
|
};
|
7f27eafa
David Mayerich
simplified the st...
|
95
|
private:
|
df036f4c
David Mayerich
added the network...
|
96
|
|
7f27eafa
David Mayerich
simplified the st...
|
97
98
99
100
|
std::vector<edge> E; //list of edges
std::vector<vertex> V; //list of vertices.
public:
|
df036f4c
David Mayerich
added the network...
|
101
|
|
7f27eafa
David Mayerich
simplified the st...
|
102
103
104
|
///Returns the number of edges in the network.
unsigned int edges(){
return E.size();
|
df036f4c
David Mayerich
added the network...
|
105
106
|
}
|
7f27eafa
David Mayerich
simplified the st...
|
107
108
109
|
///Returns the number of nodes in the network.
unsigned int vertices(){
return V.size();
|
df036f4c
David Mayerich
added the network...
|
110
111
|
}
|
3800c27f
David Mayerich
added code to sub...
|
112
113
114
115
|
stim::fiber<T> get_fiber(unsigned f){
return E[f]; //return the specified edge (casting it to a fiber)
}
|
7f27eafa
David Mayerich
simplified the st...
|
116
117
|
//load a network from an OBJ file
void load_obj(std::string filename){
|
e376212f
David Mayerich
added support for...
|
118
|
|
7f27eafa
David Mayerich
simplified the st...
|
119
120
|
stim::obj<T> O; //create an OBJ object
O.load(filename); //load the OBJ file as an object
|
df036f4c
David Mayerich
added the network...
|
121
|
|
7f27eafa
David Mayerich
simplified the st...
|
122
|
std::vector<unsigned> id2vert; //this list stores the OBJ vertex ID associated with each network vertex
|
df036f4c
David Mayerich
added the network...
|
123
|
|
4df9ec0e
pranathivemuri
compare networks ...
|
124
|
unsigned i[2]; //temporary, IDs associated with the first and last points in an OBJ line
|
df036f4c
David Mayerich
added the network...
|
125
|
|
7f27eafa
David Mayerich
simplified the st...
|
126
127
|
//for each line in the OBJ object
for(unsigned int l = 1; l <= O.numL(); l++){
|
df036f4c
David Mayerich
added the network...
|
128
|
|
7f27eafa
David Mayerich
simplified the st...
|
129
130
|
std::vector< stim::vec<T> > c; //allocate an array of points for the vessel centerline
O.getLine(l, c); //get the fiber centerline
|
df036f4c
David Mayerich
added the network...
|
131
|
|
3800c27f
David Mayerich
added code to sub...
|
132
|
edge new_edge = c; //create an edge from the given centerline
|
df036f4c
David Mayerich
added the network...
|
133
|
|
7f27eafa
David Mayerich
simplified the st...
|
134
135
136
137
138
|
//get the first and last vertex IDs for the line
std::vector< unsigned > id; //create an array to store the centerline point IDs
O.getLinei(l, id); //get the list of point IDs for the line
i[0] = id.front(); //get the OBJ ID for the first element of the line
i[1] = id.back(); //get the OBJ ID for the last element of the line
|
df036f4c
David Mayerich
added the network...
|
139
|
|
7f27eafa
David Mayerich
simplified the st...
|
140
141
|
std::vector<unsigned>::iterator it; //create an iterator for searching the id2vert array
unsigned it_idx; //create an integer for the id2vert entry index
|
df036f4c
David Mayerich
added the network...
|
142
|
|
7f27eafa
David Mayerich
simplified the st...
|
143
144
145
146
|
//find out if the nodes for this fiber have already been created
it = find(id2vert.begin(), id2vert.end(), i[0]); //look for the first node
it_idx = std::distance(id2vert.begin(), it);
if(it == id2vert.end()){ //if i[0] hasn't already been used
|
4df9ec0e
pranathivemuri
compare networks ...
|
147
148
149
150
|
vertex new_vertex = new_edge[0]; //create a new vertex, assign it a position
new_vertex.e[0].push_back(E.size()); //add the current edge as outgoing
new_edge.v[0] = V.size(); //add the new vertex to the edge
V.push_back(new_vertex); //add the new vertex to the vertex list
|
7f27eafa
David Mayerich
simplified the st...
|
151
|
id2vert.push_back(i[0]); //add the ID to the ID->vertex conversion list
|
df036f4c
David Mayerich
added the network...
|
152
|
}
|
7f27eafa
David Mayerich
simplified the st...
|
153
|
else{ //if the vertex already exists
|
df480014
pranathivemuri
added in comments...
|
154
|
V[it_idx].e[0].push_back(E.size()); //add the current edge as outgoing
|
3800c27f
David Mayerich
added code to sub...
|
155
|
new_edge.v[0] = it_idx;
|
df036f4c
David Mayerich
added the network...
|
156
157
|
}
|
7f27eafa
David Mayerich
simplified the st...
|
158
159
160
|
it = find(id2vert.begin(), id2vert.end(), i[1]); //look for the second ID
it_idx = std::distance(id2vert.begin(), it);
if(it == id2vert.end()){ //if i[1] hasn't already been used
|
3800c27f
David Mayerich
added code to sub...
|
161
162
163
164
|
vertex new_vertex = new_edge.back(); //create a new vertex, assign it a position
new_vertex.e[1].push_back(E.size()); //add the current edge as incoming
new_edge.v[1] = V.size();
V.push_back(new_vertex); //add the new vertex to the vertex list
|
7f27eafa
David Mayerich
simplified the st...
|
165
|
id2vert.push_back(i[1]); //add the ID to the ID->vertex conversion list
|
df036f4c
David Mayerich
added the network...
|
166
|
}
|
7f27eafa
David Mayerich
simplified the st...
|
167
|
else{ //if the vertex already exists
|
df480014
pranathivemuri
added in comments...
|
168
|
V[it_idx].e[1].push_back(E.size()); //add the current edge as incoming
|
3800c27f
David Mayerich
added code to sub...
|
169
|
new_edge.v[1] = it_idx;
|
df036f4c
David Mayerich
added the network...
|
170
|
}
|
3800c27f
David Mayerich
added code to sub...
|
171
|
|
4df9ec0e
pranathivemuri
compare networks ...
|
172
|
E.push_back(new_edge); //push the edge to the list
|
df036f4c
David Mayerich
added the network...
|
173
|
|
df036f4c
David Mayerich
added the network...
|
174
|
}
|
e376212f
David Mayerich
added support for...
|
175
176
|
}
|
7f27eafa
David Mayerich
simplified the st...
|
177
178
|
/// Output the network as a string
std::string str(){
|
e376212f
David Mayerich
added support for...
|
179
|
|
7f27eafa
David Mayerich
simplified the st...
|
180
181
182
183
|
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...
|
184
185
|
}
|
7f27eafa
David Mayerich
simplified the st...
|
186
187
188
|
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...
|
189
190
|
}
|
7f27eafa
David Mayerich
simplified the st...
|
191
|
return ss.str();
|
df036f4c
David Mayerich
added the network...
|
192
|
}
|
3800c27f
David Mayerich
added code to sub...
|
193
194
195
196
197
198
|
/// This function resamples all fibers in a network given a desired minimum spacing
stim::network<T> resample(T spacing){
stim::network<T> n; //create a new network that will be an exact copy, with resampled fibers
n.V = V; //copy all vertices
|
4df9ec0e
pranathivemuri
compare networks ...
|
199
|
n.E.resize(edges()); //allocate space for the edge list
|
3800c27f
David Mayerich
added code to sub...
|
200
201
|
//copy all fibers, resampling them in the process
|
4df9ec0e
pranathivemuri
compare networks ...
|
202
|
for(unsigned e = 0; e < edges(); e++){ //for each edge in the edge list
|
3800c27f
David Mayerich
added code to sub...
|
203
204
205
|
n.E[e] = E[e].resample(spacing); //resample the edge and copy it to the new network
}
|
4df9ec0e
pranathivemuri
compare networks ...
|
206
207
208
209
210
211
212
213
214
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
|
return n; //return the resampled network
}
// this function gives sum of lengths of all the fibers in the network
float lengthOfNetwork(){
stim::fiber<T> FIBER; // initialize a fiber used in looping through all edges casting into fibers in the network
float networkLength=0;float N=0; // initialize variables for finding total length of all the fibers
// for each edge in the network
for (unsigned int i=0; i < E.size(); i ++)
{
FIBER = get_fiber(i); // cast each edge to fiber
N= FIBER.length(); // find length of the fiber
networkLength = networkLength + N; // running sum of fiber lengths
}
return networkLength;
}
// list of all the points after resampling -- function used only on a resampled network
std::vector<stim::vec<T> > points_afterResampling(){
std::vector<stim::vec<T> > pointsVector; // points in the resampled network
stim::fiber<T> FIBER; // initialize a fiber used in looping through all edges casting into fibers in the network
std::vector<stim::vec<T> > pointsFiber; // resampled points on each fiber of the network
for(unsigned e = 0; e < E.size(); e++){ // for each edge in the edge list
FIBER = get_fiber(e); // Cast edge to a fiber
pointsFiber = FIBER.centerline(); // find points on the edge returns a stim vec
for (unsigned v = 0; v < FIBER.n_pts(); v++){ // iterate one point at a time from the stim::vec
pointsVector.push_back(pointsFiber[v]);} //add the points on centerline to the stim::vec points list
}
return pointsVector;
}
// total number of points on all fibers after resampling -- function used only on a resampled network
|
797f8ab0
David Mayerich
cleaned up the st...
|
241
242
243
244
245
|
unsigned total_points(){
unsigned n = 0;
for(unsigned e = 0; e < E.size(); e++)
n += E[e].size();
//unsigned int n = points_afterResampling().size();
|
4df9ec0e
pranathivemuri
compare networks ...
|
246
|
return n;
|
797f8ab0
David Mayerich
cleaned up the st...
|
247
|
}
|
4df9ec0e
pranathivemuri
compare networks ...
|
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
|
// gaussian function
float gaussianFunction(float x, float std=25){ return exp(-x/(2*std*std));} // by default std = 25
// sum of values in a stim vector
float sum(stim::vec<T> metricList){
float sumMetricList = 0; // Initialize variable to calculate sum
for (unsigned int count=0; count<metricList.size(); count++) // for each element in the stim vector
{ sumMetricList += metricList[count];} // running sum of values
return sumMetricList;
}
// distance between two points
double dist(stim::vec<T> p0, stim::vec<T> p1){
double sum = 0; // initialize variables
stim::vec<T> v = p1 - p0;; // direction vector
for(unsigned int d = 0; d < 3; d++){ // for each dimension
sum += v[d] * v[d];} // running sum of modulus of direction vector
return sqrt(sum);
}
|
797f8ab0
David Mayerich
cleaned up the st...
|
270
271
272
273
274
275
|
void stim2ann(ANNpoint &a, stim::vec<T> b){
a[0] = b[0];
a[1] = b[1];
a[2] = b[2];
}
|
4df9ec0e
pranathivemuri
compare networks ...
|
276
277
|
/// This function compares two networks and returns a metric
|
797f8ab0
David Mayerich
cleaned up the st...
|
278
|
float compare(stim::network<T> A, float sigma){
|
4df9ec0e
pranathivemuri
compare networks ...
|
279
280
281
282
283
284
|
float metric = 0.0; // initialize metric to be returned after comparing the networks
ANNkd_tree* kdt; // initialize a pointer to a kd tree
double **c; // centerline (array of double pointers) - points on kdtree must be double
unsigned int n_data = total_points(); // set the number of points
c = (double**) malloc(sizeof(double*) * n_data); // allocate the array pointer
for(unsigned int i = 0; i < n_data; i++) // allocate space for each point of 3 dimensions
|
797f8ab0
David Mayerich
cleaned up the st...
|
285
286
287
288
289
290
291
292
|
c[i] = (double*) malloc(sizeof(double) * 3);
//std::vector<stim::vec<T> > pointsVector = points_afterResampling(); //get vertices in the network after resampling
unsigned t = 0;
for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network
for(unsigned p = 0; p < E[e].size(); p++){ //for each point in the edge
for(unsigned d = 0; d < 3; d++){ //for each coordinate
c[t][d] = E[e][p][d];
|
4df9ec0e
pranathivemuri
compare networks ...
|
293
|
}
|
797f8ab0
David Mayerich
cleaned up the st...
|
294
|
t++;
|
4df9ec0e
pranathivemuri
compare networks ...
|
295
|
}
|
797f8ab0
David Mayerich
cleaned up the st...
|
296
297
|
}
|
4df9ec0e
pranathivemuri
compare networks ...
|
298
299
300
|
ANNpointArray pts = (ANNpointArray)c; // create an array of data points of type double
kdt = new ANNkd_tree(pts, n_data, 3); // build a KD tree using the annpointarray
double eps = 0; // error bound
|
797f8ab0
David Mayerich
cleaned up the st...
|
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
|
ANNdistArray dists = new ANNdist[1]; // near neighbor distances
ANNidxArray nnIdx = new ANNidx[1]; // near neighbor indices // allocate near neigh indices
stim::vec<T> p0, p1;
float m0, m1;
float M = 0; //stores the total metric value
float l; //stores the segment length
float L = 0; //stores the total network length
ANNpoint queryPt = annAllocPt(3);
for(unsigned e = 0; e < A.E.size(); e++){ //for each edge in A
M = 0; //total metric value for the fiber
p0 = A.E[e][0]; //get the first point in the edge
stim2ann(queryPt, p0);
kdt->annkSearch( queryPt, 1, nnIdx, dists, eps); //find the distance between A and the current network
m0 = 1.0f - gaussianFunction(dists[0], sigma); //calculate the metric value based on the distance
for(unsigned p = 1; p < A.E[e].size(); p++){ //for each point in the edge
p1 = A.E[e][p]; //get the next point in the edge
stim2ann(queryPt, p1);
kdt->annkSearch( queryPt, 1, nnIdx, dists, eps); //find the distance between A and the current network
m1 = 1.0f - gaussianFunction(dists[0], sigma); //calculate the metric value based on the distance
//average the metric and scale
l = (p1 - p0).len(); //calculate the length of the segment
L += l; //add the length of this segment to the total network length
M += (m0 + m1)/2 * l; //scale the metric based on segment length
//swap points
p0 = p1;
m0 = m1;
}
}
return M/L;
/*int N; // number of points on the fiber in the second network
float totalNetworkLength = A.lengthOfNetwork();
stim::vec<float> fiberMetric(A.edges()); // allocate space for accumulating fiber metric as each fiber is evaluated
|
4df9ec0e
pranathivemuri
compare networks ...
|
341
342
|
stim::fiber<T> FIBER; // Initialize a fiber will be used in calculating the metric
//for each fiber in the second network, find nearest distance in the kd tree
|
797f8ab0
David Mayerich
cleaned up the st...
|
343
|
for (unsigned int i=0; i < A.edges(); i ++) // loop through all the edges/fibers in the network
|
4df9ec0e
pranathivemuri
compare networks ...
|
344
|
{
|
797f8ab0
David Mayerich
cleaned up the st...
|
345
|
FIBER = A.get_fiber(i); // Get the fiber corresponding to the index i
|
4df9ec0e
pranathivemuri
compare networks ...
|
346
|
std::vector< stim::vec<T> > fiberPoints = FIBER.centerline(); // Get the points on the fiber
|
797f8ab0
David Mayerich
cleaned up the st...
|
347
|
N = FIBER.size(); // number of points on the fiber
|
4df9ec0e
pranathivemuri
compare networks ...
|
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
|
stim::vec<float> segmentMetric(N-1); // allocate space for a vec array that stores metric for each segmen in the fiber
// for each segment in the fiber
for (unsigned j = 0; j < N - 1; j++)
{
stim::vec<T> p1 = fiberPoints[j];stim::vec<T> p2 = fiberPoints[j+1]; // starting and ending points on the segments
ANNpoint queryPt1; queryPt1 = annAllocPt(3); // allocate memory for query points
ANNpoint queryPt2; queryPt2 = annAllocPt(3);
//for each dimension of the points connecting segment
for (unsigned d = 0; d < 3; d++){
queryPt1[d] = double(fiberPoints[j][d]); // starting point on segment in network whose closest distance on KD tree should be found
queryPt2[d] = double(fiberPoints[j + 1][d]); // ending point on segment in network whose closest distance on KD tree should be found
}
kdt->annkSearch( queryPt1, 1, nnIdx1, dists1, eps); // search the nearest point in KD tree to query point(point on network), find the distance
kdt->annkSearch( queryPt2, 1, nnIdx2, dists2, eps); // search the nearest point in KD tree to query point(point on network), find the distance
// find the gaussian of the distance and subtract it from 1 to calculate the error
float error1 = 1.0f - gaussianFunction(float(dists1[0]), sigma);float error2 = 1.0f - gaussianFunction(float(dists2[0]), sigma);
// average the error and scale it with distance between the points
segmentMetric[j] = ((error1 + error2) / 2) * dist(p1, p2) ;
|
4df9ec0e
pranathivemuri
compare networks ...
|
367
368
|
}
fiberMetric[i] = sum(segmentMetric);
|
4df9ec0e
pranathivemuri
compare networks ...
|
369
|
}
|
4df9ec0e
pranathivemuri
compare networks ...
|
370
371
|
metric = sum(fiberMetric)/totalNetworkLength; //normalize the scaled error of all the points with total length of the network
assert (0<=metric <=1); //assert metroc os always less than or equal to one and greater than or equal to zero
|
797f8ab0
David Mayerich
cleaned up the st...
|
372
|
return metric;*/
|
3800c27f
David Mayerich
added code to sub...
|
373
|
}
|
7f27eafa
David Mayerich
simplified the st...
|
374
375
|
}; //end stim::network class
}; //end stim namespace
|
df036f4c
David Mayerich
added the network...
|
376
|
#endif
|