network.h
15 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
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
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
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
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
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
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
#ifndef STIM_NETWORK_H
#define STIM_NETWORK_H
#include <stdlib.h>
#include <assert.h>
#include <sstream>
#include <fstream>
#include <algorithm>
#include <string.h>
#include <math.h>
#include <stim/math/vector.h>
#include <stim/visualization/obj.h>
#include <stim/biomodels/fiber.h>
#include <ANN/ANN.h>
#include <boost/tuple/tuple.hpp>
namespace stim{
/** 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.
*/
template<typename T>
class network{
///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
// default constructor
edge() : fiber<T>(){v[1] = -1; v[0] = -1;}
/// Constructor - creates an edge from a list of points by calling the stim::fiber constructor
///@param p is an array of positions in space
edge(std::vector< stim::vec<T> > p) : fiber<T>(p){}
/// 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){
edge e(fiber<T>::resample(spacing)); //call the fiber->edge constructor
e.v[0] = v[0]; //copy the vertex data
e.v[1] = v[1];
return e; //return the new edge
}
/// Output the edge information as a string
std::string str(){
std::stringstream ss;
ss<<"("<<fiber<T>::N<<")\tl = "<<length()<<"\t"<<v[0]<<"----"<<v[1];
return ss.str();
}
};
///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])
//stim::vec<T> p; //position of this node in physical space.
//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]<<" ";
}
return ss.str();
}
};
private:
std::vector<edge> E; //list of edges
std::vector<vertex> V; //list of vertices.
public:
///Returns the number of edges in the network.
unsigned int edges(){
return E.size();
}
///Returns the number of nodes in the network.
unsigned int vertices(){
return V.size();
}
stim::fiber<T> get_fiber(unsigned f){
return E[f]; //return the specified edge (casting it to a fiber)
}
//load a network from an OBJ file
void load_obj(std::string filename){
stim::obj<T> O; //create an OBJ object
O.load(filename); //load the OBJ file as an object
std::vector<unsigned> id2vert; //this list stores the OBJ vertex ID associated with each network vertex
unsigned i[2]; //temporary, IDs associated with the first and last points in an OBJ line
//for each line in the OBJ object
for(unsigned int l = 1; l <= O.numL(); l++){
std::vector< stim::vec<T> > c; //allocate an array of points for the vessel centerline
O.getLine(l, c); //get the fiber centerline
edge new_edge = c; //create an edge from the given centerline
//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
std::vector<unsigned>::iterator it; //create an iterator for searching the id2vert array
unsigned it_idx; //create an integer for the id2vert entry index
//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
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
id2vert.push_back(i[0]); //add the ID to the ID->vertex conversion list
}
else{ //if the vertex already exists
V[it_idx].e[0].push_back(E.size()); //add the current edge as outgoing
new_edge.v[0] = it_idx;
}
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
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
id2vert.push_back(i[1]); //add the ID to the ID->vertex conversion list
}
else{ //if the vertex already exists
V[it_idx].e[1].push_back(E.size()); //add the current edge as incoming
new_edge.v[1] = it_idx;
}
E.push_back(new_edge); //push the edge to the list
}
}
/// Output the network as a string
std::string str(){
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;
}
ss<<"Edges ("<<E.size()<<")--------"<<std::endl;
for(unsigned e = 0; e < E.size(); e++){
ss<<"\t"<<e<<E[e].str()<<std::endl;
}
return ss.str();
}
/// 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
n.E.resize(edges()); //allocate space for the edge list
//copy all fibers, resampling them in the process
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
}
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
unsigned total_points(){
unsigned n = 0;
for(unsigned e = 0; e < E.size(); e++)
n += E[e].size();
//unsigned int n = points_afterResampling().size();
return n;
}
// 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);
}
void stim2ann(ANNpoint &a, stim::vec<T> b){
a[0] = b[0];
a[1] = b[1];
a[2] = b[2];
}
/// This function compares two networks and returns a metric
float compare(stim::network<T> A, float sigma){
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
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];
}
t++;
}
}
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
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
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
for (unsigned int i=0; i < A.edges(); i ++) // loop through all the edges/fibers in the network
{
FIBER = A.get_fiber(i); // Get the fiber corresponding to the index i
std::vector< stim::vec<T> > fiberPoints = FIBER.centerline(); // Get the points on the fiber
N = FIBER.size(); // number of points on the fiber
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) ;
}
fiberMetric[i] = sum(segmentMetric);
}
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
return metric;*/
}
}; //end stim::network class
}; //end stim namespace
#endif