Commit e0edbe1317a40526ff0bc4e06ad8fd38e8c62b50

Authored by Laila Saadatifard
2 parents 03428452 9f5c0d4a

Merge branch 'master' of git.stim.ee.uh.edu:codebase/stimlib

stim/biomodels/network.h
@@ -33,7 +33,10 @@ class network{ @@ -33,7 +33,10 @@ class network{
33 public: 33 public:
34 unsigned v[2]; //unique id's designating the starting and ending 34 unsigned v[2]; //unique id's designating the starting and ending
35 // default constructor 35 // default constructor
36 - edge() : cylinder<T>(){v[1] = -1; v[0] = -1;} 36 + edge() : cylinder<T>()
  37 + {
  38 + v[1] = -1; v[0] = -1;
  39 + }
37 /// Constructor - creates an edge from a list of points by calling the stim::fiber constructor 40 /// Constructor - creates an edge from a list of points by calling the stim::fiber constructor
38 41
39 ///@param p is an array of positions in space 42 ///@param p is an array of positions in space
@@ -54,7 +57,7 @@ class network{ @@ -54,7 +57,7 @@ class network{
54 /// Output the edge information as a string 57 /// Output the edge information as a string
55 std::string str(){ 58 std::string str(){
56 std::stringstream ss; 59 std::stringstream ss;
57 - ss<<"("<<cylinder<T>::size()<<")\tl = "<<length()<<"\t"<<v[0]<<"----"<<v[1]; 60 + ss<<"("<<cylinder<T>::size()<<")\tl = "<<this.length()<<"\t"<<v[0]<<"----"<<v[1];
58 return ss.str(); 61 return ss.str();
59 } 62 }
60 63
@@ -64,9 +67,9 @@ class network{ @@ -64,9 +67,9 @@ class network{
64 class vertex : public stim::vec3<T> 67 class vertex : public stim::vec3<T>
65 { 68 {
66 public: 69 public:
67 - //std::vector<unsigned int> edges; //indices of edges connected to this node.  
68 - std::vector<unsigned int> e[2]; //indices of edges going out (e[0]) and coming in (e[1])  
69 - //stim::vec3<T> p; //position of this node in physical space. 70 + //std::vector<unsigned int> edges; //indices of edges connected to this node.
  71 + std::vector<unsigned int> e[2]; //indices of edges going out (e[0]) and coming in (e[1])
  72 + //stim::vec3<T> p; //position of this node in physical space.
70 73
71 //constructor takes a stim::vec 74 //constructor takes a stim::vec
72 vertex(stim::vec3<T> p) : stim::vec3<T>(p){} 75 vertex(stim::vec3<T> p) : stim::vec3<T>(p){}
@@ -90,6 +93,7 @@ class network{ @@ -90,6 +93,7 @@ class network{
90 return ss.str(); 93 return ss.str();
91 } 94 }
92 95
  96 +
93 }; 97 };
94 98
95 protected: 99 protected:
@@ -99,6 +103,18 @@ protected: @@ -99,6 +103,18 @@ protected:
99 103
100 public: 104 public:
101 105
  106 + ///default constructor
  107 + network()
  108 + {
  109 +
  110 + }
  111 +
  112 + ///constructor with a file to load.
  113 + network(std::string fileLocation)
  114 + {
  115 + load_obj(fileLocation);
  116 + }
  117 +
102 ///Returns the number of edges in the network. 118 ///Returns the number of edges in the network.
103 unsigned int edges(){ 119 unsigned int edges(){
104 return E.size(); 120 return E.size();
@@ -109,72 +125,219 @@ public: @@ -109,72 +125,219 @@ public:
109 return V.size(); 125 return V.size();
110 } 126 }
111 127
  128 + std::vector<vertex> operator*(T s){
  129 + for (unsigned i=0; i< vertices; i ++ ){
  130 + V[i] = V[i] * s;
  131 + }
  132 + return V;
  133 + }
  134 +
  135 + std::vector<vertex> operator*(vec<T> s){
  136 + for (unsigned i=0; i< vertices; i ++ ){
  137 + for (unsigned dim = 0 ; dim< 3; dim ++){
  138 + V[i][dim] = V[i][dim] * s[dim];
  139 + }
  140 + }
  141 + return V;
  142 + }
  143 +
  144 + // Returns an average of branching index in the network
  145 +
  146 + double BranchingIndex(){
  147 + double B=0;
  148 + for(unsigned v=0; v < V.size(); v ++){
  149 + B += ((V[v].e[0].size()) + (V[v].e[1].size()));
  150 + }
  151 + B = B / V.size();
  152 + return B;
  153 +
  154 + }
  155 +
  156 + // Returns number of branch points in thenetwork
  157 +
  158 + unsigned int BranchP(){
  159 + unsigned int B=0;
  160 + unsigned int c;
  161 + for(unsigned v=0; v < V.size(); v ++){
  162 + c = ((V[v].e[0].size()) + (V[v].e[1].size()));
  163 + if (c > 2){
  164 + B += 1;}
  165 + }
  166 + return B;
  167 +
  168 + }
  169 +
  170 + // Returns number of end points (tips) in thenetwork
  171 +
  172 + unsigned int EndP(){
  173 + unsigned int B=0;
  174 + unsigned int c;
  175 + for(unsigned v=0; v < V.size(); v ++){
  176 + c = ((V[v].e[0].size()) + (V[v].e[1].size()));
  177 + if (c == 1){
  178 + B += 1;}
  179 + }
  180 + return B;
  181 +
  182 + }
  183 +
  184 + //// Returns a dictionary with the key as the vertex
  185 + //std::map<std::vector<vertex>,unsigned int> DegreeDict(){
  186 + // std::map<std::vector<vertex>,unsigned int> dd;
  187 + // unsigned int c = 0;
  188 + // for(unsigned v=0; v < V.size(); v ++){
  189 + // c = ((V[v].e[0].size()) + (V[v].e[1].size()));
  190 + // dd[V[v]] = c;
  191 + // }
  192 + // return dd;
  193 + //}
  194 +
  195 + //// Return number of branching stems
  196 + //unsigned int Stems(){
  197 + // unsigned int s = 0;
  198 + // std::map<std::vector<vertex>,unsigned int> dd;
  199 + // dd = DegreeDict();
  200 + // //for(unsigned v=0; v < V.size(); v ++){
  201 + // // V[v].e[0].
  202 + // return s;
  203 + //}
  204 +
  205 +
  206 + // Returns an average of fiber/edge lengths in the network
  207 + double Lengths(){
  208 + stim::vec<T> L;double sumLength = 0;
  209 + for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network
  210 + L.push_back(E[e].length()); //append the edge length
  211 + sumLength = sumLength + E[e].length();
  212 + }
  213 + double avg = sumLength / E.size();
  214 + return avg;
  215 + }
  216 +
  217 +
  218 + // Returns an average of tortuosities in the network
  219 + double Tortuosities(){
  220 + stim::vec<T> t;
  221 + stim::vec<T> id1, id2; // starting and ending vertices of the edge
  222 + double distance;double tortuosity;double sumTortuosity = 0;
  223 + for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network
  224 + id1 = E[e][0]; //get the edge starting point
  225 + id2 = E[e][E[e].size() - 1]; //get the edge ending point
  226 + distance = (id1 - id2).len(); //displacement between the starting and ending points
  227 + if(distance > 0){
  228 + tortuosity = E[e].length()/ distance ; // tortuoisty = edge length / edge displacement
  229 + }
  230 + else{
  231 + tortuosity = 0;}
  232 + t.push_back(tortuosity);
  233 + sumTortuosity += tortuosity;
  234 + }
  235 + double avg = sumTortuosity / E.size();
  236 + return avg;
  237 + }
  238 +
  239 + // Returns average contraction of the network
  240 + double Contractions(){
  241 + stim::vec<T> t;
  242 + stim::vec<T> id1, id2; // starting and ending vertices of the edge
  243 + double distance;double contraction;double sumContraction = 0;
  244 + for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network
  245 + id1 = E[e][0]; //get the edge starting point
  246 + id2 = E[e][E[e].size() - 1]; //get the edge ending point
  247 + distance = (id1 - id2).len(); //displacement between the starting and ending points
  248 + contraction = distance / E[e].length(); // tortuoisty = edge length / edge displacement
  249 + t.push_back(contraction);
  250 + sumContraction += contraction;
  251 + }
  252 + double avg = sumContraction / E.size();
  253 + return avg;
  254 + }
  255 +
  256 + // returns average fractal dimension of the branches of the network
  257 + double FractalDimensions(){
  258 + stim::vec<T> t;
  259 + stim::vec<T> id1, id2; // starting and ending vertices of the edge
  260 + double distance;double fract;double sumFractDim = 0;
  261 + for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network
  262 + id1 = E[e][0]; //get the edge starting point
  263 + id2 = E[e][E[e].size() - 1]; //get the edge ending point
  264 + distance = (id1 - id2).len(); //displacement between the starting and ending points
  265 + fract = std::log(distance) / std::log(E[e].length()); // tortuoisty = edge length / edge displacement
  266 + t.push_back(sumFractDim);
  267 + sumFractDim += fract;
  268 + }
  269 + double avg = sumFractDim / E.size();
  270 + return avg;
  271 + }
112 stim::cylinder<T> get_cylinder(unsigned f){ 272 stim::cylinder<T> get_cylinder(unsigned f){
113 - return E[f]; //return the specified edge (casting it to a fiber) 273 + return E[f]; //return the specified edge (casting it to a fiber)
114 } 274 }
115 275
116 //load a network from an OBJ file 276 //load a network from an OBJ file
117 void load_obj(std::string filename){ 277 void load_obj(std::string filename){
118 278
119 - stim::obj<T> O; //create an OBJ object  
120 - O.load(filename); //load the OBJ file as an object 279 + stim::obj<T> O; //create an OBJ object
  280 + O.load(filename); //load the OBJ file as an object
121 281
122 - std::vector<unsigned> id2vert; //this list stores the OBJ vertex ID associated with each network vertex 282 + std::vector<unsigned> id2vert; //this list stores the OBJ vertex ID associated with each network vertex
123 283
124 - unsigned i[2]; //temporary, IDs associated with the first and last points in an OBJ line 284 + unsigned i[2]; //temporary, IDs associated with the first and last points in an OBJ line
125 285
126 //for each line in the OBJ object 286 //for each line in the OBJ object
127 for(unsigned int l = 1; l <= O.numL(); l++){ 287 for(unsigned int l = 1; l <= O.numL(); l++){
128 288
129 - std::vector< stim::vec<T> > c; //allocate an array of points for the vessel centerline 289 + std::vector< stim::vec<T> > c; //allocate an array of points for the vessel centerline
130 O.getLine(l, c); //get the fiber centerline 290 O.getLine(l, c); //get the fiber centerline
131 291
132 std::vector< stim::vec3<T> > c3(c.size()); 292 std::vector< stim::vec3<T> > c3(c.size());
133 for(size_t j = 0; j < c.size(); j++) 293 for(size_t j = 0; j < c.size(); j++)
134 c3[j] = c[j]; 294 c3[j] = c[j];
135 295
136 - edge new_edge = c3; //create an edge from the given centerline  
137 - unsigned int I = new_edge.size(); //calculate the number of points on the centerline 296 + // edge new_edge = c3; ///This is dangerous.
  297 + edge new_edge(c3);
  298 +
  299 + //create an edge from the given centerline
  300 + unsigned int I = new_edge.size(); //calculate the number of points on the centerline
138 301
139 //get the first and last vertex IDs for the line 302 //get the first and last vertex IDs for the line
140 - std::vector< unsigned > id; //create an array to store the centerline point IDs 303 + std::vector< unsigned > id; //create an array to store the centerline point IDs
141 O.getLinei(l, id); //get the list of point IDs for the line 304 O.getLinei(l, id); //get the list of point IDs for the line
142 i[0] = id.front(); //get the OBJ ID for the first element of the line 305 i[0] = id.front(); //get the OBJ ID for the first element of the line
143 i[1] = id.back(); //get the OBJ ID for the last element of the line 306 i[1] = id.back(); //get the OBJ ID for the last element of the line
144 307
145 - std::vector<unsigned>::iterator it; //create an iterator for searching the id2vert array 308 + std::vector<unsigned>::iterator it; //create an iterator for searching the id2vert array
146 unsigned it_idx; //create an integer for the id2vert entry index 309 unsigned it_idx; //create an integer for the id2vert entry index
147 310
148 //find out if the nodes for this fiber have already been created 311 //find out if the nodes for this fiber have already been created
149 it = find(id2vert.begin(), id2vert.end(), i[0]); //look for the first node 312 it = find(id2vert.begin(), id2vert.end(), i[0]); //look for the first node
150 - it_idx = std::distance(id2vert.begin(), it);  
151 if(it == id2vert.end()){ //if i[0] hasn't already been used 313 if(it == id2vert.end()){ //if i[0] hasn't already been used
152 vertex new_vertex = new_edge[0]; //create a new vertex, assign it a position 314 vertex new_vertex = new_edge[0]; //create a new vertex, assign it a position
153 - new_vertex.e[0].push_back(E.size()); //add the current edge as outgoing  
154 - new_edge.v[0] = V.size(); //add the new vertex to the edge  
155 - V.push_back(new_vertex); //add the new vertex to the vertex list  
156 - id2vert.push_back(i[0]); //add the ID to the ID->vertex conversion list 315 + new_vertex.e[0].push_back(E.size()); //add the current edge as outgoing
  316 + new_edge.v[0] = V.size(); //add the new edge to the edge
  317 + V.push_back(new_vertex); //add the new vertex to the vertex list
  318 + id2vert.push_back(i[0]); //add the ID to the ID->vertex conversion list
157 } 319 }
158 - else{ //if the vertex already exists 320 + else{ //if the vertex already exists
  321 + it_idx = std::distance(id2vert.begin(), it);
159 V[it_idx].e[0].push_back(E.size()); //add the current edge as outgoing 322 V[it_idx].e[0].push_back(E.size()); //add the current edge as outgoing
160 new_edge.v[0] = it_idx; 323 new_edge.v[0] = it_idx;
161 } 324 }
162 325
163 - it = find(id2vert.begin(), id2vert.end(), i[1]); //look for the second ID  
164 - it_idx = std::distance(id2vert.begin(), it);  
165 - if(it == id2vert.end()){ //if i[1] hasn't already been used  
166 - vertex new_vertex = new_edge[I-1]; //create a new vertex, assign it a position  
167 - new_vertex.e[1].push_back(E.size()); //add the current edge as incoming  
168 - new_edge.v[1] = V.size();  
169 - V.push_back(new_vertex); //add the new vertex to the vertex list  
170 - id2vert.push_back(i[1]); //add the ID to the ID->vertex conversion list 326 + it = find(id2vert.begin(), id2vert.end(), i[1]); //look for the second ID
  327 + if(it == id2vert.end()){ //if i[1] hasn't already been used
  328 + vertex new_vertex = new_edge[I-1]; //create a new vertex, assign it a position
  329 + new_vertex.e[1].push_back(E.size()); //add the current edge as incoming
  330 + new_edge.v[1] = V.size(); //add the new vertex to the edge
  331 + V.push_back(new_vertex); //add the new vertex to the vertex list
  332 + id2vert.push_back(i[1]); //add the ID to the ID->vertex conversion list
171 } 333 }
172 - else{ //if the vertex already exists 334 + else{ //if the vertex already exists
  335 + it_idx = std::distance(id2vert.begin(), it);
173 V[it_idx].e[1].push_back(E.size()); //add the current edge as incoming 336 V[it_idx].e[1].push_back(E.size()); //add the current edge as incoming
174 new_edge.v[1] = it_idx; 337 new_edge.v[1] = it_idx;
175 } 338 }
176 339
177 - E.push_back(new_edge); //push the edge to the list 340 + E.push_back(new_edge); //push the edge to the list
178 341
179 } 342 }
180 } 343 }
@@ -199,17 +362,17 @@ public: @@ -199,17 +362,17 @@ public:
199 /// This function resamples all fibers in a network given a desired minimum spacing 362 /// This function resamples all fibers in a network given a desired minimum spacing
200 /// @param spacing is the minimum distance between two points on the network 363 /// @param spacing is the minimum distance between two points on the network
201 stim::network<T> resample(T spacing){ 364 stim::network<T> resample(T spacing){
202 - stim::network<T> n; //create a new network that will be an exact copy, with resampled fibers  
203 - n.V = V; //copy all vertices 365 + stim::network<T> n; //create a new network that will be an exact copy, with resampled fibers
  366 + n.V = V; //copy all vertices
204 367
205 - n.E.resize(edges()); //allocate space for the edge list 368 + n.E.resize(edges()); //allocate space for the edge list
206 369
207 //copy all fibers, resampling them in the process 370 //copy all fibers, resampling them in the process
208 - for(unsigned e = 0; e < edges(); e++){ //for each edge in the edge list  
209 - n.E[e] = E[e].resample(spacing); //resample the edge and copy it to the new network 371 + for(unsigned e = 0; e < edges(); e++){ //for each edge in the edge list
  372 + n.E[e] = E[e].resample(spacing); //resample the edge and copy it to the new network
210 } 373 }
211 374
212 - return n; //return the resampled network 375 + return n; //return the resampled network
213 } 376 }
214 377
215 378
@@ -236,14 +399,14 @@ public: @@ -236,14 +399,14 @@ public:
236 /// @param m is the magnitude value to use. The default is 0 (usually radius). 399 /// @param m is the magnitude value to use. The default is 0 (usually radius).
237 T average(unsigned m = 0){ 400 T average(unsigned m = 0){
238 401
239 - T M, L; //allocate space for the total magnitude and length  
240 - M = L = 0; //initialize both the initial magnitude and length to zero  
241 - for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network 402 + T M, L; //allocate space for the total magnitude and length
  403 + M = L = 0; //initialize both the initial magnitude and length to zero
  404 + for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network
242 M += E[e].integrate(m); //get the integrated magnitude 405 M += E[e].integrate(m); //get the integrated magnitude
243 - L += E[e].length(); //get the edge length 406 + L += E[e].length(); //get the edge length
244 } 407 }
245 408
246 - return M / L; //return the average magnitude 409 + return M / L; //return the average magnitude
247 } 410 }
248 411
249 /// This function compares two networks and returns the percentage of the current network that is missing from A. 412 /// This function compares two networks and returns the percentage of the current network that is missing from A.
@@ -256,17 +419,17 @@ public: @@ -256,17 +419,17 @@ public:
256 R = (*this); //initialize the result with the current network 419 R = (*this); //initialize the result with the current network
257 420
258 //generate a KD-tree for network A 421 //generate a KD-tree for network A
259 - float metric = 0.0; // initialize metric to be returned after comparing the networks  
260 - ANNkd_tree* kdt; // initialize a pointer to a kd tree  
261 - double **c; // centerline (array of double pointers) - points on kdtree must be double  
262 - unsigned int n_data = A.total_points(); // set the number of points  
263 - c = (double**) malloc(sizeof(double*) * n_data); // allocate the array pointer  
264 - for(unsigned int i = 0; i < n_data; i++) // allocate space for each point of 3 dimensions 422 + float metric = 0.0; // initialize metric to be returned after comparing the networks
  423 + ANNkd_tree* kdt; // initialize a pointer to a kd tree
  424 + double **c; // centerline (array of double pointers) - points on kdtree must be double
  425 + unsigned int n_data = A.total_points(); // set the number of points
  426 + c = (double**) malloc(sizeof(double*) * n_data); // allocate the array pointer
  427 + for(unsigned int i = 0; i < n_data; i++) // allocate space for each point of 3 dimensions
265 c[i] = (double*) malloc(sizeof(double) * 3); 428 c[i] = (double*) malloc(sizeof(double) * 3);
266 429
267 unsigned t = 0; 430 unsigned t = 0;
268 for(unsigned e = 0; e < A.E.size(); e++){ //for each edge in the network 431 for(unsigned e = 0; e < A.E.size(); e++){ //for each edge in the network
269 - for(unsigned p = 0; p < A.E[e].size(); p++){ //for each point in the edge 432 + for(unsigned p = 0; p < A.E[e].size(); p++){ //for each point in the edge
270 for(unsigned d = 0; d < 3; d++){ //for each coordinate 433 for(unsigned d = 0; d < 3; d++){ //for each coordinate
271 434
272 c[t][d] = A.E[e][p][d]; 435 c[t][d] = A.E[e][p][d];
@@ -276,27 +439,27 @@ public: @@ -276,27 +439,27 @@ public:
276 } 439 }
277 440
278 //compare each point in the current network to the field produced by A 441 //compare each point in the current network to the field produced by A
279 - ANNpointArray pts = (ANNpointArray)c; // create an array of data points of type double  
280 - kdt = new ANNkd_tree(pts, n_data, 3); // build a KD tree using the annpointarray 442 + ANNpointArray pts = (ANNpointArray)c; // create an array of data points of type double
  443 + kdt = new ANNkd_tree(pts, n_data, 3); // build a KD tree using the annpointarray
281 double eps = 0; // error bound 444 double eps = 0; // error bound
282 - ANNdistArray dists = new ANNdist[1]; // near neighbor distances  
283 - ANNidxArray nnIdx = new ANNidx[1]; // near neighbor indices // allocate near neigh indices 445 + ANNdistArray dists = new ANNdist[1]; // near neighbor distances
  446 + ANNidxArray nnIdx = new ANNidx[1]; // near neighbor indices // allocate near neigh indices
284 447
285 stim::vec3<T> p0, p1; 448 stim::vec3<T> p0, p1;
286 float m1; 449 float m1;
287 - float M = 0; //stores the total metric value  
288 - float L = 0; //stores the total network length 450 + float M = 0; //stores the total metric value
  451 + float L = 0; //stores the total network length
289 ANNpoint queryPt = annAllocPt(3); 452 ANNpoint queryPt = annAllocPt(3);
290 for(unsigned e = 0; e < R.E.size(); e++){ //for each edge in A 453 for(unsigned e = 0; e < R.E.size(); e++){ //for each edge in A
291 - R.E[e].add_mag(0); //add a new magnitude for the metric 454 + R.E[e].add_mag(0); //add a new magnitude for the metric
292 455
293 - for(unsigned p = 0; p < R.E[e].size(); p++){ //for each point in the edge 456 + for(unsigned p = 0; p < R.E[e].size(); p++){ //for each point in the edge
294 457
295 - p1 = R.E[e][p]; //get the next point in the edge 458 + p1 = R.E[e][p]; //get the next point in the edge
296 stim2ann(queryPt, p1); 459 stim2ann(queryPt, p1);
297 - kdt->annkSearch( queryPt, 1, nnIdx, dists, eps); //find the distance between A and the current network 460 + kdt->annkSearch( queryPt, 1, nnIdx, dists, eps); //find the distance between A and the current network
298 m1 = 1.0f - gaussianFunction((float)dists[0], sigma); //calculate the metric value based on the distance 461 m1 = 1.0f - gaussianFunction((float)dists[0], sigma); //calculate the metric value based on the distance
299 - R.E[e].set_mag(m1, p, 1); //set the error for the second point in the segment 462 + R.E[e].set_mag(m1, p, 1); //set the error for the second point in the segment
300 463
301 } 464 }
302 } 465 }
@@ -308,8 +471,106 @@ public: @@ -308,8 +471,106 @@ public:
308 unsigned nmags(){ 471 unsigned nmags(){
309 return E[0].nmags(); 472 return E[0].nmags();
310 } 473 }
  474 + // split a string in text by the character sep
  475 + stim::vec<T> split(std::string &text, char sep)
  476 + {
  477 + stim::vec<T> tokens;
  478 + std::size_t start = 0, end = 0;
  479 + while ((end = text.find(sep, start)) != std::string::npos) {
  480 + tokens.push_back(atof(text.substr(start, end - start).c_str()));
  481 + start = end + 1;
  482 + }
  483 + tokens.push_back(atof(text.substr(start).c_str()));
  484 + return tokens;
  485 + }
  486 + // load a network in text file to a network class
  487 + void load_txt(std::string filename)
  488 + {
  489 + std::vector <std::string> file_contents;
  490 + std::ifstream file(filename);
  491 + std::string line;
  492 + std::vector<unsigned> id2vert; //this list stores the vertex ID associated with each network vertex
  493 + //for each line in the text file, store them as strings in file_contents
  494 + while (std::getline(file, line))
  495 + {
  496 + std::stringstream ss(line);
  497 + file_contents.push_back(ss.str());
  498 + }
  499 + unsigned int numEdges = atoi(file_contents[0].c_str()); //number of edges in the network
  500 + unsigned int I = atoi(file_contents[1].c_str()) ; //calculate the number of points3d on the first edge
  501 + unsigned int count = 1; unsigned int k = 2; // count is global counter through the file contents, k is for the vertices on the edges
  502 + // for each edge in the network.
  503 + for (unsigned int i = 0; i < numEdges; i ++ )
  504 + {
  505 + // pre allocate a position vector p with number of points3d on the edge p
  506 + std::vector< stim::vec<T> > p(0, I);
  507 + // for each point on the nth edge
  508 + for (unsigned int j = k; j < I + k; j++)
  509 + {
  510 + // split the points3d of floats with separator space and form a float3 position vector out of them
  511 + p.push_back(split(file_contents[j], ' '));
  512 + }
  513 + count += p.size() + 1; // increment count to point at the next edge in the network
  514 + I = atoi(file_contents[count].c_str()); // read in the points3d at the next edge and convert it to an integer
  515 + k = count + 1;
  516 + edge new_edge = p; // create an edge with a vector of points3d on the edge
  517 + E.push_back(new_edge); // push the edge into the network
  518 + }
  519 + unsigned int numVertices = atoi(file_contents[count].c_str()); // this line in the text file gives the number of distinct vertices
  520 + count = count + 1; // this line of text file gives the first verrtex
  521 + // push each vertex into V
  522 + for (unsigned int i = 0; i < numVertices; i ++)
  523 + {
  524 + vertex new_vertex = split(file_contents[count], ' ');
  525 + V.push_back(new_vertex);
  526 + count += atoi(file_contents[count + 1].c_str()) + 2; // Skip number of edge ids + 2 to point to the next vertex
  527 + }
  528 + } // end load_txt function
311 529
312 - 530 + // strTxt returns a string of edges
  531 + std::string
  532 + strTxt(std::vector< stim::vec<T> > p)
  533 + {
  534 + std::stringstream ss;
  535 + std::stringstream oss;
  536 + for(unsigned int i = 0; i < p.size(); i++){
  537 + ss.str(std::string());
  538 + for(unsigned int d = 0; d < 3; d++){
  539 + ss<<p[i][d];
  540 + }
  541 + ss < "\n";
  542 + }
  543 + return ss.str();
  544 + }
  545 + // removes specified character from string
  546 + void removeCharsFromString(std::string &str, char* charsToRemove ) {
  547 + for ( unsigned int i = 0; i < strlen(charsToRemove); ++i ) {
  548 + str.erase( remove(str.begin(), str.end(), charsToRemove[i]), str.end() );
  549 + }
  550 + }
  551 + //exports network to txt file
  552 + void
  553 + to_txt(std::string filename)
  554 + {
  555 + std::ofstream ofs(filename, std::ofstream::out | std::ofstream::app);
  556 + int num;
  557 + ofs << (E.size()).str() << "\n";
  558 + for(unsigned int i = 0; i < E.size(); i++)
  559 + {
  560 + std::string str;
  561 + ofs << (E[i].size()).str() << "\n";
  562 + str = E[i].strTxt();
  563 + ofs << str << "\n";
  564 + }
  565 + for(int i = 0; i < V.size(); i++)
  566 + {
  567 + std::string str;
  568 + str = V[i].str();
  569 + removeCharsFromString(str, "[],");
  570 + ofs << str << "\n";
  571 + }
  572 + ofs.close();
  573 + }
313 }; //end stim::network class 574 }; //end stim::network class
314 }; //end stim namespace 575 }; //end stim namespace
315 #endif 576 #endif
stim/cuda/cuda_texture.cuh
@@ -134,15 +134,15 @@ namespace stim @@ -134,15 +134,15 @@ namespace stim
134 void 134 void
135 UnmapCudaTexture() 135 UnmapCudaTexture()
136 { 136 {
137 - HANDLE_ERROR(  
138 - cudaGraphicsUnmapResources(1, &resource)  
139 - );  
140 - HANDLE_ERROR(  
141 - cudaGraphicsUnregisterResource(resource)  
142 - );  
143 - HANDLE_ERROR(  
144 - cudaDestroyTextureObject(tObj)  
145 - ); 137 + // HANDLE_ERROR(
  138 + // cudaGraphicsUnmapResources(1, &resource)
  139 + // );
  140 + // HANDLE_ERROR(
  141 + // cudaGraphicsUnregisterResource(resource)
  142 + // );
  143 + // HANDLE_ERROR(
  144 + // cudaDestroyTextureObject(tObj)
  145 + // );
146 // HANDLE_ERROR( 146 // HANDLE_ERROR(
147 // cudaFreeArray(srcArray) 147 // cudaFreeArray(srcArray)
148 // ); 148 // );
stim/cuda/filter.cuh
@@ -7,6 +7,7 @@ @@ -7,6 +7,7 @@
7 #include <stdio.h> 7 #include <stdio.h>
8 #include <stim/visualization/colormap.h> 8 #include <stim/visualization/colormap.h>
9 #include <sstream> 9 #include <sstream>
  10 +#include <stim/math/constants.h>
10 #include <stim/cuda/cudatools/devices.h> 11 #include <stim/cuda/cudatools/devices.h>
11 #include <stim/cuda/cudatools/threads.h> 12 #include <stim/cuda/cudatools/threads.h>
12 #include <stim/cuda/cuda_texture.cuh> 13 #include <stim/cuda/cuda_texture.cuh>
@@ -14,7 +15,6 @@ @@ -14,7 +15,6 @@
14 #include <stim/cuda/arraymath.cuh> 15 #include <stim/cuda/arraymath.cuh>
15 16
16 #define IMAD(a,b,c) ( __mul24((a), (b)) + (c) ) 17 #define IMAD(a,b,c) ( __mul24((a), (b)) + (c) )
17 -#define M_PI 3.141592654f  
18 18
19 19
20 namespace stim 20 namespace stim
@@ -73,7 +73,7 @@ namespace stim @@ -73,7 +73,7 @@ namespace stim
73 idx = j*kl+i; 73 idx = j*kl+i;
74 x = i - kr - 0.5; 74 x = i - kr - 0.5;
75 y = j - kr - 0.5; 75 y = j - kr - 0.5;
76 - LoG[idx] = (-1.0/M_PI/powf(sigma, 4))* (1 - (powf(x,2)+powf(y,2))/2.0/powf(sigma, 2)) 76 + LoG[idx] = (-1.0/PI/powf(sigma, 4))* (1 - (powf(x,2)+powf(y,2))/2.0/powf(sigma, 2))
77 *expf(-(powf(x,2)+powf(y,2))/2/powf(sigma,2)); 77 *expf(-(powf(x,2)+powf(y,2))/2/powf(sigma,2));
78 t +=LoG[idx]; 78 t +=LoG[idx];
79 } 79 }
@@ -98,7 +98,7 @@ namespace stim @@ -98,7 +98,7 @@ namespace stim
98 int y = blockIdx.y; 98 int y = blockIdx.y;
99 int xi = threadIdx.x; 99 int xi = threadIdx.x;
100 int yi = threadIdx.y; 100 int yi = threadIdx.y;
101 - float val = 0; 101 + // float val = 0;
102 float tu = (x-kr+xi)/(float)DIM_X; 102 float tu = (x-kr+xi)/(float)DIM_X;
103 float tv = (y-kr+yi)/(float)DIM_Y; 103 float tv = (y-kr+yi)/(float)DIM_Y;
104 shared[xi][yi] = gpuLoG[yi*kl+xi]*(255.0-(float)tex2D<unsigned char>(texIn, tu, tv)); 104 shared[xi][yi] = gpuLoG[yi*kl+xi]*(255.0-(float)tex2D<unsigned char>(texIn, tu, tv));
@@ -111,7 +111,7 @@ namespace stim @@ -111,7 +111,7 @@ namespace stim
111 //y = min(y, height - 1); 111 //y = min(y, height - 1);
112 112
113 int idx = y*DIM_X+x; 113 int idx = y*DIM_X+x;
114 - int k_idx; 114 + // int k_idx;
115 for(unsigned int step = blockDim.x/2; step >= 1; step >>= 1) 115 for(unsigned int step = blockDim.x/2; step >= 1; step >>= 1)
116 { 116 {
117 __syncthreads(); 117 __syncthreads();
stim/cuda/ivote/down_sample.cuh
@@ -44,7 +44,7 @@ namespace stim{ @@ -44,7 +44,7 @@ namespace stim{
44 unsigned int y_ds = (y/sigma_ds + (y %sigma_ds == 0 ? 0:1)); 44 unsigned int y_ds = (y/sigma_ds + (y %sigma_ds == 0 ? 0:1));
45 45
46 //get the number of pixels in the image 46 //get the number of pixels in the image
47 - unsigned int pixels_ds = x_ds * y_ds; 47 +// unsigned int pixels_ds = x_ds * y_ds;
48 48
49 unsigned int max_threads = stim::maxThreadsPerBlock(); 49 unsigned int max_threads = stim::maxThreadsPerBlock();
50 dim3 threads(max_threads, 1); 50 dim3 threads(max_threads, 1);
@@ -97,4 +97,4 @@ namespace stim{ @@ -97,4 +97,4 @@ namespace stim{
97 } 97 }
98 } 98 }
99 99
100 -#endif  
101 \ No newline at end of file 100 \ No newline at end of file
  101 +#endif
stim/envi/agilent_binary.h 0 → 100644
  1 +//make sure that this header file is only loaded once
  2 +#ifndef STIM_AGILENT_BINARY_H
  3 +#define STIM_AGILENT_BINARY_H
  4 +
  5 +#include <string>
  6 +#include <fstream>
  7 +
  8 +//CUDA
  9 +#ifdef CUDA_FOUND
  10 + #include <cuda_runtime.h>
  11 + #include "cufft.h"
  12 + #include <stim/cuda/cudatools/error.h>
  13 +#endif
  14 +
  15 +namespace stim{
  16 +
  17 +template<typename T>
  18 +class agilent_binary{
  19 +
  20 +protected:
  21 + std::string fname;
  22 + T* ptr;
  23 + size_t R[3];
  24 + static const size_t header = 1020;
  25 + double Z[2];
  26 +
  27 +public:
  28 + size_t size(){
  29 + return (size_t)R[0] * (size_t)R[1] * (size_t)R[2];
  30 + }
  31 +
  32 + size_t bytes(){
  33 + return size() * sizeof(T);
  34 + }
  35 + void alloc(){
  36 + ptr = (T*) malloc(bytes());
  37 + }
  38 + void alloc(short x, short y, short z){
  39 + R[0] = x;
  40 + R[1] = y;
  41 + R[2] = z;
  42 + alloc();
  43 + }
  44 +
  45 + void deep_copy(agilent_binary<T>* dst, const agilent_binary<T>* src){
  46 + dst->alloc(src->R[0], src->R[1], src->R[2]); //allocate memory
  47 + memcpy(dst->ptr, src->ptr, bytes()); //copy the data
  48 + memcpy(dst->Z, src->Z, sizeof(double) * 2); //copy the data z range
  49 + }
  50 +
  51 + agilent_binary(){
  52 + memset(R, 0, sizeof(short) * 3); //set the resolution to zero
  53 + ptr = NULL;
  54 + }
  55 +
  56 + /// Constructor with resolution
  57 + agilent_binary(short x, short y, short z){
  58 + alloc(x, y, z);
  59 + }
  60 +
  61 + /// Constructor with filename
  62 + agilent_binary(std::string filename){
  63 + ptr = NULL;
  64 + load(filename);
  65 + }
  66 +
  67 + /// Copy constructor
  68 + agilent_binary(const agilent_binary<T> &obj){
  69 + deep_copy(this, &obj);
  70 + }
  71 +
  72 + agilent_binary<T>& operator=(const agilent_binary<T> rhs){
  73 + if(this != &rhs){ //check for self-assignment
  74 + deep_copy(this, &rhs); //make a deep copy
  75 + }
  76 + return *this; //return the result
  77 + }
  78 +
  79 + ~agilent_binary(){
  80 + free(ptr);
  81 + }
  82 +
  83 + void load(std::string filename){
  84 + if(ptr != NULL) free(ptr); //if memory has been allocated, free it
  85 +
  86 + fname = filename; //save the filename
  87 +
  88 + short x, y, z;
  89 +
  90 + std::ifstream infile(fname, std::ios::binary); //open the input file
  91 + infile.seekg(9, std::ios::beg); //seek past 9 bytes from the beginning of the file
  92 +
  93 + infile.read((char*)(&z), 2); //read two bytes of data (the number of samples is stored as a 16-bit integer)
  94 +
  95 + infile.seekg(13, std::ios::cur); //skip another 13 bytes
  96 + infile.read((char*)(&x), 2); //read the X and Y dimensions
  97 + infile.read((char*)(&y), 2);
  98 +
  99 + infile.seekg(header, std::ios::beg); //seek to the start of the data
  100 +
  101 + alloc(x, y, z);
  102 + ptr = (T*) malloc(bytes()); //allocate space for the data
  103 + infile.read((char*)ptr, bytes()); //read the data
  104 + infile.close();
  105 + }
  106 +
  107 + void save(std::string filename){
  108 + std::ofstream outfile(filename, std::ios::binary); //create an output file
  109 +
  110 + char zero = 0;
  111 + for(size_t i = 0; i < 9; i++) outfile.write(&zero, 1); //write 9 zeros
  112 + outfile.write((char*)&R[0], 2);
  113 + for(size_t i = 0; i < 13; i++) outfile.write(&zero, 1); //write 13 zeros
  114 + outfile.write((char*)&R[1], 2);
  115 + outfile.write((char*)&R[2], 2);
  116 + for(size_t i = 0; i < 992; i++) outfile.write(&zero, 1); //write 992 zeros
  117 + //char zerovec[1020];
  118 + //outfile.write((char*)zerovec, 1020);
  119 +
  120 + size_t b = bytes();
  121 + outfile.write((char*)ptr, b); //write the data to the output file
  122 + outfile.close();
  123 + }
  124 +
  125 + stim::envi_header create_header(){
  126 + stim::envi_header header;
  127 + header.samples = R[0];
  128 + header.lines = R[1];
  129 + header.bands = R[2];
  130 +
  131 + double z_delta = (double)(Z[1] - Z[0]) / (double)(R[2] - 1);
  132 + header.wavelength.resize(R[2]);
  133 + for(size_t i = 0; i < R[2]; i++)
  134 + header.wavelength[i] = i * z_delta + Z[0];
  135 +
  136 + return header;
  137 + }
  138 +
  139 + /// Calculate the absorbance spectrum from the transmission spectrum given a background
  140 + void absorbance(stim::agilent_binary<T>* background){
  141 + size_t N = size(); //calculate the number of values to be ratioed
  142 + if(N != background->size()){
  143 + std::cerr<<"ERROR in stim::agilent_binary::absorbance() - transmission image size doesn't match background"<<std::endl;
  144 + exit(1);
  145 + }
  146 + for(size_t i = 0; i < N; i++)
  147 + ptr[i] = -log10(ptr[i] / background->ptr[i]);
  148 + }
  149 +
  150 +#ifdef CUDA_FOUND
  151 + /// Perform an FFT and return a binary file with bands in the specified range
  152 + agilent_binary<T> fft(float band_min, float band_max){
  153 + auto total_start = std::chrono::high_resolution_clock::now();
  154 +
  155 + auto start = std::chrono::high_resolution_clock::now();
  156 + T* cpu_data = (T*) malloc( bytes() ); //allocate space for the transposed data
  157 + for(size_t b = 0; b < R[2]; b++){
  158 + for(size_t x = 0; x < R[0] * R[1]; x++){
  159 + cpu_data[x * R[2] + b] = ptr[b * R[0] * R[1] + x];
  160 + }
  161 + }
  162 + auto end = std::chrono::high_resolution_clock::now();
  163 + std::chrono::duration<double> diff = end-start;
  164 + // std::cout << "Transpose data: " << diff.count() << " s\n";
  165 +
  166 + start = std::chrono::high_resolution_clock::now();
  167 + cufftHandle plan; //allocate space for a cufft plan
  168 + cufftReal* gpu_data; //create a pointer to the data
  169 + size_t batch = R[0] * R[1]; //calculate the batch size (X * Y)
  170 + HANDLE_ERROR(cudaMalloc((void**)&gpu_data, bytes())); //allocate space on the GPU
  171 + HANDLE_ERROR(cudaMemcpy(gpu_data, cpu_data, bytes(), cudaMemcpyHostToDevice)); //copy the data to the GPU
  172 + cufftComplex* gpu_fft;
  173 + HANDLE_ERROR(cudaMalloc((void**)&gpu_fft, R[0] * R[1] * (R[2]/2 + 1) * sizeof(cufftComplex)));
  174 + end = std::chrono::high_resolution_clock::now();
  175 + diff = end-start;
  176 + //std::cout << "Allocate/copy: " << diff.count() << " s\n";
  177 +
  178 + start = std::chrono::high_resolution_clock::now();
  179 + int N[1]; //create an array with the interferogram size (required for cuFFT input)
  180 + N[0] = R[2]; //set the only array value to the interferogram size
  181 + if(cufftPlanMany(&plan, 1, N, NULL, 1, R[2], NULL, 1, R[2], CUFFT_R2C, batch) != CUFFT_SUCCESS){
  182 + std::cout<<"cuFFT Error: unable to create 1D plan."<<std::endl;
  183 + exit(1);
  184 + }
  185 + end = std::chrono::high_resolution_clock::now();
  186 + diff = end-start;
  187 + //std::cout << "Create a plan: " << diff.count() << " s\n";
  188 +
  189 + start = std::chrono::high_resolution_clock::now();
  190 + if (cufftExecR2C(plan, gpu_data, gpu_fft) != CUFFT_SUCCESS){ //execute the (implicitly forward) transform
  191 + std::cout<<"CUFFT error: ExecR2C Forward failed";
  192 + exit(1);
  193 + }
  194 + end = std::chrono::high_resolution_clock::now();
  195 + diff = end-start;
  196 + //std::cout << "Perform FFT: " << diff.count() << " s\n";
  197 +
  198 + start = std::chrono::high_resolution_clock::now();
  199 + std::complex<T>* cpu_fft = (std::complex<T>*) malloc( R[0] * R[1] * (R[2]/2+1) * sizeof(std::complex<T>) );
  200 + HANDLE_ERROR(cudaMemcpy(cpu_fft, gpu_fft, R[0] * R[1] * (R[2]/2+1) * sizeof(cufftComplex), cudaMemcpyDeviceToHost)); //copy data from the host to the device
  201 +
  202 + double int_delta = 0.00012656; //interferogram sample spacing in centimeters
  203 + double int_length = int_delta * R[2]; //interferogram length in centimeters
  204 + double fft_delta = 1/int_length; //spectrum spacing (in inverse centimeters, wavenumber
  205 +
  206 + size_t start_i = std::ceil(band_min / fft_delta); //calculate the first band to store
  207 + size_t size_i = std::floor(band_max / fft_delta) - start_i; //calculate the number of bands to store
  208 + size_t end_i = start_i + size_i; //last band number
  209 + agilent_binary<T> result(R[0], R[1], size_i);
  210 + result.Z[0] = start_i * fft_delta; //set the range for the FFT result
  211 + result.Z[1] = end_i * fft_delta;
  212 +
  213 + for(size_t b = start_i; b < end_i; b++){
  214 + for(size_t x = 0; x < R[0] * R[1]; x++){
  215 + result.ptr[(b - start_i) * R[0] * R[1] + x] = abs(cpu_fft[x * (R[2]/2+1) + b]);
  216 + }
  217 + }
  218 + end = std::chrono::high_resolution_clock::now();
  219 + diff = end-start;
  220 + //std::cout << "Transpose/Crop: " << diff.count() << " s\n";
  221 +
  222 + auto total_end = std::chrono::high_resolution_clock::now();
  223 + diff = total_end-total_start;
  224 +
  225 + cufftDestroy(plan);
  226 + HANDLE_ERROR(cudaFree(gpu_data));
  227 + HANDLE_ERROR(cudaFree(gpu_fft));
  228 + free(cpu_data);
  229 + free(cpu_fft);
  230 +
  231 + return result;
  232 + }
  233 +#endif
  234 +
  235 +};
  236 +
  237 +}
  238 +
  239 +#endif
0 \ No newline at end of file 240 \ No newline at end of file
@@ -235,6 +235,7 @@ public: @@ -235,6 +235,7 @@ public:
235 235
236 for(unsigned long long i = 0; i < B; i++) //for each element of the spectrum 236 for(unsigned long long i = 0; i < B; i++) //for each element of the spectrum
237 p[i] = (double) temp[i]; //cast each element to a double value 237 p[i] = (double) temp[i]; //cast each element to a double value
  238 + free(temp); //free temporary memory
238 return true; 239 return true;
239 } 240 }
240 241
stim/gl/gl_spider.h
@@ -14,6 +14,7 @@ @@ -14,6 +14,7 @@
14 #include <stim/math/vec3.h> 14 #include <stim/math/vec3.h>
15 #include <stim/math/rect.h> 15 #include <stim/math/rect.h>
16 #include <stim/math/matrix.h> 16 #include <stim/math/matrix.h>
  17 +#include <stim/math/constants.h>
17 #include <stim/cuda/spider_cost.cuh> 18 #include <stim/cuda/spider_cost.cuh>
18 #include <stim/cuda/cudatools/glbind.h> 19 #include <stim/cuda/cudatools/glbind.h>
19 #include <stim/cuda/arraymath.cuh> 20 #include <stim/cuda/arraymath.cuh>
@@ -369,7 +370,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -369,7 +370,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
369 ///Stored in a display list. 370 ///Stored in a display list.
370 ///uses the default d vector <0,0,1> 371 ///uses the default d vector <0,0,1>
371 void 372 void
372 - genDirectionVectors(float solidAngle = 5/M_PI*4) 373 + genDirectionVectors(float solidAngle = M_PI/2)
373 { 374 {
374 375
375 //Set up the vectors necessary for Rectangle creation. 376 //Set up the vectors necessary for Rectangle creation.
@@ -584,7 +585,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -584,7 +585,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
584 glDeleteFramebuffers(1, &framebufferID); 585 glDeleteFramebuffers(1, &framebufferID);
585 glGenFramebuffers(1, &framebufferID); 586 glGenFramebuffers(1, &framebufferID);
586 glBindFramebuffer(GL_FRAMEBUFFER, framebufferID); 587 glBindFramebuffer(GL_FRAMEBUFFER, framebufferID);
587 - int numChannels = 1; 588 +// int numChannels = 1;
588 // unsigned char* texels = new unsigned char[width * height * numChannels]; 589 // unsigned char* texels = new unsigned char[width * height * numChannels];
589 glGenTextures(1, &textureID); 590 glGenTextures(1, &textureID);
590 glBindTexture(GL_TEXTURE_2D, textureID); 591 glBindTexture(GL_TEXTURE_2D, textureID);
@@ -609,7 +610,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -609,7 +610,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
609 { 610 {
610 glGenFramebuffers(1, &fboID); 611 glGenFramebuffers(1, &fboID);
611 glBindFramebuffer(GL_FRAMEBUFFER, fboID); 612 glBindFramebuffer(GL_FRAMEBUFFER, fboID);
612 - int numChannels = 1; 613 +// int numChannels = 1;
613 // unsigned char* texels = new unsigned char[width * height * numChannels]; 614 // unsigned char* texels = new unsigned char[width * height * numChannels];
614 glGenTextures(1, &texbufferID); 615 glGenTextures(1, &texbufferID);
615 glBindTexture(GL_TEXTURE_2D, texbufferID); 616 glBindTexture(GL_TEXTURE_2D, texbufferID);
@@ -938,7 +939,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -938,7 +939,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
938 network_time = 0; 939 network_time = 0;
939 hit_time = 0; 940 hit_time = 0;
940 #endif 941 #endif
941 - stepsize = 3.0; 942 + stepsize = 2.5;
942 t_length = 16.0; 943 t_length = 16.0;
943 944
944 srand(100); 945 srand(100);
@@ -1392,13 +1393,9 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1392,13 +1393,9 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1392 1393
1393 Z[0] = cos(PHI[0]); 1394 Z[0] = cos(PHI[0]);
1394 Z[1] = cos(PHI[1]); 1395 Z[1] = cos(PHI[1]);
1395 -  
1396 1396
1397 range = Z[0] - Z[1]; 1397 range = Z[0] - Z[1];
1398 1398
1399 -  
1400 - float z, theta, phi;  
1401 -  
1402 std::vector<stim::vec3<float> > vecsUni; 1399 std::vector<stim::vec3<float> > vecsUni;
1403 for(int i = 0; i < numSamplesPos; i++) 1400 for(int i = 0; i < numSamplesPos; i++)
1404 { 1401 {
@@ -1490,8 +1487,6 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1490,8 +1487,6 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1490 void 1487 void
1491 trace(int min_cost) 1488 trace(int min_cost)
1492 { 1489 {
1493 - bool sEmpty = true;  
1494 - float lastmag = 16.0;;  
1495 stim::vec3<float> curSeed; 1490 stim::vec3<float> curSeed;
1496 stim::vec3<float> curSeedVec; 1491 stim::vec3<float> curSeedVec;
1497 float curSeedMag; 1492 float curSeedMag;
@@ -1527,7 +1522,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1527,7 +1522,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1527 gpuStartTimer(); 1522 gpuStartTimer();
1528 #endif 1523 #endif
1529 1524
1530 - float s = 3.0; 1525 +// float s = 3.0;
1531 GLuint selectBuf[2048]; 1526 GLuint selectBuf[2048];
1532 GLint hits; 1527 GLint hits;
1533 glSelectBuffer(2048, selectBuf); 1528 glSelectBuffer(2048, selectBuf);
@@ -1541,17 +1536,17 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1541,17 +1536,17 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1541 CHECK_OPENGL_ERROR 1536 CHECK_OPENGL_ERROR
1542 //What would that vessel see in front of it. 1537 //What would that vessel see in front of it.
1543 camSel.setPosition(loc); 1538 camSel.setPosition(loc);
1544 - camSel.setFocalDistance(mag/s);  
1545 - camSel.LookAt((loc[0]+dir[0]*mag/s),  
1546 - (loc[1]+dir[1]*mag/s),  
1547 - (loc[2]+dir[2]*mag/s)); 1539 + camSel.setFocalDistance(mag/stepsize);
  1540 + camSel.LookAt((loc[0]+dir[0]*mag/stepsize),
  1541 + (loc[1]+dir[1]*mag/stepsize),
  1542 + (loc[2]+dir[2]*mag/stepsize));
1548 ps = camSel.getPosition(); 1543 ps = camSel.getPosition();
1549 ups = camSel.getUp(); 1544 ups = camSel.getUp();
1550 ds = camSel.getLookAt(); 1545 ds = camSel.getLookAt();
1551 glMatrixMode(GL_PROJECTION); 1546 glMatrixMode(GL_PROJECTION);
1552 glPushMatrix(); 1547 glPushMatrix();
1553 glLoadIdentity(); 1548 glLoadIdentity();
1554 - glOrtho(-mag/s/2.0, mag/s/2.0, -mag/s/2.0, mag/s/2.0, 0.0, mag/s/2.0); 1549 + glOrtho(-mag/stepsize/2.0, mag/stepsize/2.0, -mag/stepsize/2.0, mag/stepsize/2.0, 0.0, mag/stepsize/2.0);
1555 glMatrixMode(GL_MODELVIEW); 1550 glMatrixMode(GL_MODELVIEW);
1556 glPushMatrix(); 1551 glPushMatrix();
1557 glLoadIdentity(); 1552 glLoadIdentity();
@@ -1597,21 +1592,11 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1597,21 +1592,11 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1597 int 1592 int
1598 processHits(GLint hits, GLuint buffer[]) 1593 processHits(GLint hits, GLuint buffer[])
1599 { 1594 {
1600 - GLuint names, *ptr;  
1601 - //printf("hits = %u\n", hits); 1595 + GLuint *ptr;
1602 ptr = (GLuint *) buffer; 1596 ptr = (GLuint *) buffer;
1603 - // for (int i = 0; i < hits; i++) { /* for each hit */  
1604 - names = *ptr;  
1605 - // printf (" number of names for hit = %u\n", names);  
1606 ptr++; 1597 ptr++;
1607 ptr++; //Skip the minimum depth value. 1598 ptr++; //Skip the minimum depth value.
1608 ptr++; //Skip the maximum depth value. 1599 ptr++; //Skip the maximum depth value.
1609 - // printf (" the name is ");  
1610 - // for (int j = 0; j < names; j++) { /* for each name */  
1611 - // printf ("%u ", *ptr); ptr++;  
1612 - // }  
1613 - // printf ("\n");  
1614 - // }  
1615 1600
1616 1601
1617 if(hits == 0) 1602 if(hits == 0)
@@ -1675,11 +1660,18 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1675,11 +1660,18 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1675 } 1660 }
1676 1661
1677 #ifdef TIMING 1662 #ifdef TIMING
1678 - double network_time = (std::clock() - s) / (double) CLOCKS_PER_SEC;  
1679 - network_time += network_time * 1000.0; 1663 + double nt = (std::clock() - s) / (double) CLOCKS_PER_SEC;
  1664 + network_time += nt * 1000.0;
1680 #endif 1665 #endif
1681 } 1666 }
1682 1667
  1668 +// void
  1669 +// addToNetwork(pair<stim::fiber<float>, int> in, stim::vec3<float> spos,
  1670 +// stim::vec<float> smag, stim::vec3<float> sdir)
  1671 +// {
  1672 +//
  1673 +// }
  1674 +
1683 1675
1684 void 1676 void
1685 printSizes() 1677 printSizes()
@@ -1783,9 +1775,6 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1783,9 +1775,6 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1783 } 1775 }
1784 } 1776 }
1785 } 1777 }
1786 -  
1787 -  
1788 -  
1789 }; 1778 };
1790 } 1779 }
1791 #endif 1780 #endif
stim/grids/image_stack.h
@@ -33,10 +33,10 @@ public: @@ -33,10 +33,10 @@ public:
33 stim::filename file_path(file_mask); 33 stim::filename file_path(file_mask);
34 34
35 //if the file path is relative, update it with the current working directory 35 //if the file path is relative, update it with the current working directory
36 - if(file_path.is_relative()){  
37 - stim::filename wd = stim::filename::cwd();  
38 - file_path = wd.get_relative(file_mask);  
39 - } 36 +// if(file_path.is_relative()){
  37 +// stim::filename wd = stim::filename::cwd();
  38 +// file_path = wd.get_relative(file_mask);
  39 +// }
40 40
41 //get the list of files 41 //get the list of files
42 std::vector<stim::filename> file_list = file_path.get_list(); 42 std::vector<stim::filename> file_list = file_path.get_list();
@@ -134,10 +134,10 @@ public: @@ -134,10 +134,10 @@ public:
134 stim::filename file_path(file_mask); 134 stim::filename file_path(file_mask);
135 135
136 //if the file path is relative, update it with the current working directory 136 //if the file path is relative, update it with the current working directory
137 - if(file_path.is_relative()){  
138 - stim::filename wd = stim::filename::cwd();  
139 - file_path = wd.get_relative(file_mask);  
140 - } 137 +// if(file_path.is_relative()){
  138 +// stim::filename wd = stim::filename::cwd();
  139 +// file_path = wd.get_relative(file_mask);
  140 +// }
141 141
142 //create a list of file names 142 //create a list of file names
143 std::vector<std::string> file_list = stim::wildcards::increment(file_path.str(), 0, R[3]-1, 1); 143 std::vector<std::string> file_list = stim::wildcards::increment(file_path.str(), 0, R[3]-1, 1);
stim/math/matrix.h
@@ -53,13 +53,13 @@ struct matrix @@ -53,13 +53,13 @@ struct matrix
53 53
54 template<typename Y> 54 template<typename Y>
55 vec<Y> operator*(vec<Y> rhs){ 55 vec<Y> operator*(vec<Y> rhs){
56 - unsigned int N = rhs.size(); 56 + unsigned int M = rhs.size();
57 57
58 vec<Y> result; 58 vec<Y> result;
59 - result.resize(N); 59 + result.resize(M);
60 60
61 - for(int r=0; r<N; r++)  
62 - for(int c=0; c<N; c++) 61 + for(int r=0; r<M; r++)
  62 + for(int c=0; c<M; c++)
63 result[r] += (*this)(r, c) * rhs[c]; 63 result[r] += (*this)(r, c) * rhs[c];
64 64
65 return result; 65 return result;
@@ -3,6 +3,7 @@ @@ -3,6 +3,7 @@
3 3
4 4
5 #include <stim/cuda/cudatools/callable.h> 5 #include <stim/cuda/cudatools/callable.h>
  6 +#include <cmath>
6 7
7 8
8 namespace stim{ 9 namespace stim{
@@ -216,7 +217,6 @@ public: @@ -216,7 +217,6 @@ public:
216 return result; 217 return result;
217 } 218 }
218 219
219 -#ifndef __NVCC__  
220 /// Outputs the vector as a string 220 /// Outputs the vector as a string
221 std::string str() const{ 221 std::string str() const{
222 std::stringstream ss; 222 std::stringstream ss;
@@ -234,7 +234,6 @@ public: @@ -234,7 +234,6 @@ public:
234 234
235 return ss.str(); 235 return ss.str();
236 } 236 }
237 -#endif  
238 237
239 size_t size(){ return 3; } 238 size_t size(){ return 3; }
240 239
stim/optics/scalarbeam.h
@@ -121,12 +121,12 @@ CUDA_CALLABLE void lut_lookup(T* lut_values, T* lut, T val, size_t N, T min_val, @@ -121,12 +121,12 @@ CUDA_CALLABLE void lut_lookup(T* lut_values, T* lut, T val, size_t N, T min_val,
121 121
122 template <typename T> 122 template <typename T>
123 CUDA_CALLABLE stim::complex<T> clerp(stim::complex<T> v0, stim::complex<T> v1, T t) { 123 CUDA_CALLABLE stim::complex<T> clerp(stim::complex<T> v0, stim::complex<T> v1, T t) {
124 - return stim::complex<T>( fma(t, v1.r, fma(-t, v0.r, v0.r)), fma(t, v1.i, fma(-t, v0.i, v0.i)) ); 124 + return stim::complex<T>( fmaf(t, v1.r, fmaf(-t, v0.r, v0.r)), fmaf(t, v1.i, fmaf(-t, v0.i, v0.i)) );
125 } 125 }
126 126
127 template <typename T> 127 template <typename T>
128 CUDA_CALLABLE T lerp(T v0, T v1, T t) { 128 CUDA_CALLABLE T lerp(T v0, T v1, T t) {
129 - return fma(t, v1, fma(-t, v0, v0)); 129 + return fmaf(t, v1, fmaf(-t, v0, v0));
130 } 130 }
131 131
132 #ifdef CUDA_FOUND 132 #ifdef CUDA_FOUND
stim/optics/scalarfield.h
@@ -293,6 +293,8 @@ protected: @@ -293,6 +293,8 @@ protected:
293 stim::complex<T>* E; 293 stim::complex<T>* E;
294 size_t R[2]; 294 size_t R[2];
295 locationType loc; 295 locationType loc;
  296 + using rect<T>::X;
  297 + using rect<T>::Y;
296 298
297 T* p[3]; //scalar position for each point in E 299 T* p[3]; //scalar position for each point in E
298 300
stim/parser/filename.h
@@ -2,218 +2,295 @@ @@ -2,218 +2,295 @@
2 #define STIM_FILENAME_H 2 #define STIM_FILENAME_H
3 3
4 #include <stdio.h> /* defines FILENAME_MAX */ 4 #include <stdio.h> /* defines FILENAME_MAX */
  5 +#include <sstream>
  6 +#include <vector>
  7 +#include <algorithm>
  8 +#include <iostream>
  9 +#include <iomanip>
  10 +
  11 +#include <stim/parser/parser.h>
  12 +
  13 +// set OS dependent defines, including:
  14 +// 1) path division character ('/' or '\') used for output only
  15 +// 2) command for getting the current working directory
  16 +// 3) header files for getting the current working directory
  17 +// 4) include BOOST filesystem class if compiling on GCC
  18 +#define STIM_DIV '/'
5 #ifdef _WIN32 19 #ifdef _WIN32
6 #include <windows.h> 20 #include <windows.h>
7 - #include <direct.h>  
8 - #define GetCurrentDir _getcwd  
9 - #define STIM_FILENAME_DIV '/' 21 + #include <direct.h>
  22 + #define GetCurrentDir _getcwd
  23 + #define STIM_FILENAME_DIV '\\'
10 #else 24 #else
11 - #include <unistd.h>  
12 - #define GetCurrentDir getcwd 25 + #ifdef BOOST_PRECOMPILED
  26 + #include <boost/filesystem.hpp>
  27 + #endif
  28 + #include <unistd.h>
  29 + #define GetCurrentDir getcwd
13 #define STIM_FILENAME_DIV '/' 30 #define STIM_FILENAME_DIV '/'
14 #endif 31 #endif
15 32
16 -#include <string>  
17 -#include <sstream>  
18 -#include <iomanip>  
19 -#include <vector>  
20 -#include <stack>  
21 -#include <algorithm>  
22 -#include <iostream>  
23 -  
24 -#include "../parser/parser.h"  
25 -#ifdef BOOST_PRECOMPILED  
26 -#include <boost/filesystem.hpp>  
27 -#endif  
28 namespace stim{ 33 namespace stim{
29 34
30 -//filename class designed to work with both Windows and Unix 35 +class filepath{
  36 +protected:
  37 + std::string _drive; //drive on which the file is located (used for Windows)
  38 + std::vector<std::string> _path; //path for the specified file (list of hierarchical directories)
31 39
32 -class filename{ 40 + /// replaces win32 dividers with the Linux standard (used internally by default)
  41 + std::string unix_div(std::string s) {
  42 + std::replace( s.begin(), s.end(), '\\', '/');
  43 + return s;
  44 + }
33 45
34 -private:  
35 - void init(){ 46 + /// gets the directory hierarchy for the current working directory
  47 + static std::string cwd(){
  48 + char cCurrentPath[FILENAME_MAX];
36 49
37 - absolute = false; 50 + if (!GetCurrentDir(cCurrentPath, sizeof(cCurrentPath))){
  51 + std::cout<<"ERROR getting current working directory."<<std::endl;
  52 + exit(1);
  53 + }
38 54
  55 + std::stringstream ss;
  56 + ss<<cCurrentPath;
  57 + return ss.str();
39 } 58 }
40 59
41 -protected: 60 + /// convert a relative path to an absolute path
  61 + void get_absolute(std::string &drive, std::vector<std::string> &absolute, std::vector<std::string> relative){
42 62
43 - std::string drive; //drive letter (for Windows)  
44 - std::vector<std::string> path; //directory hierarchy  
45 - std::string prefix; //file prefix (without extension)  
46 - std::string ext; //file extension 63 + std::string current = cwd(); //get the current directory as a string
47 64
48 - bool absolute; //filename is an absolute path 65 + std::string current_drive;
  66 + std::vector<std::string> current_dir;
  67 + parse_path(current_drive, current_dir, current); //get the current drive and directories
49 68
50 - //replaces win32 dividers with the Linux standard  
51 - std::string unix_div(std::string s) {  
52 - std::replace( s.begin(), s.end(), '\\', '/');  
53 - return s;  
54 - }  
55 -  
56 - //break up a filename into a prefix and extension  
57 - void parse_name(std::string fname){ 69 + // step through each directory in the relative path, adjusting the current directory
  70 + // index depending on whether the relative path is specified with '.' or '..'
  71 + int current_i = (int)current_dir.size() - 1;
  72 + int relative_i;
  73 + for(relative_i = 0; relative_i < relative.size(); relative_i++){
  74 + if(relative[relative_i] == "..") current_i--;
  75 + else if(relative[relative_i] != ".") break;
  76 + }
  77 + if(current_i < 0){
  78 + std::cerr<<"ERROR stim::filepath - relative path is incompatible with working directory"<<std::endl;
  79 + exit(1);
  80 + }
58 81
59 - //find the extension  
60 - size_t xpos = fname.find_last_of('.'); //find the position of the extension period (if there is one)  
61 - if(xpos != std::string::npos){ //split the prefix and extension  
62 - prefix = fname.substr(0, xpos);  
63 - ext = fname.substr(xpos + 1); 82 + absolute.clear();
  83 + for(size_t i = 0; i <= current_i; i++){
  84 + absolute.push_back(current_dir[i]);
  85 + }
  86 + for(size_t i = relative_i; i < relative.size(); i++){
  87 + absolute.push_back(relative[i]);
64 } 88 }
65 - else  
66 - prefix = fname;  
67 } 89 }
68 90
69 - //parse a file locator string  
70 - void parse(std::string loc){  
71 - if(loc.size() == 0) return;  
72 -  
73 - loc = unix_div(loc);  
74 -  
75 - //determine the drive (if Windows)  
76 - drive = "";  
77 - #ifdef _WIN32  
78 - //if the second character is a colon, the character right before it is the drive  
79 - if(loc[1] == ':'){  
80 - absolute = true; //this is an absolute path  
81 - drive = loc[0]; //copy the drive letter  
82 - loc = loc.substr(3); //remove the drive information from the locator string 91 + /// Parses a directory string into a drive (NULL if not Windows) and list of hierarchical directories
  92 + /// Returns true if the path is relative, false if it is absolute
  93 + void parse_path(std::string &drive, std::vector<std::string> &absolute, std::string dir){
  94 + drive = ""; //initialize the drive to NULL (it will stay that way for Windows)
  95 + std::vector<std::string> path;
  96 + bool relative = true; //the default path is relative
  97 +
  98 + if(dir.length() == 0) return; //return if the file locator is empty
  99 + std::string unix_dir = unix_div(dir); //replace any Windows division characters with Unix
  100 +
  101 + if(unix_dir.length() > 1 && unix_dir[1] == ':'){ //if a drive identifier is given
  102 + if(unix_dir[0] > 64 && unix_dir[0] < 91) //if the drive letter is upper case
  103 + _drive = unix_dir[0] + 32; //convert it to lower case
  104 + else if(unix_dir[0] > 96 && unix_dir[0] < 123) //if the drive character is lower case
  105 + _drive = unix_dir[0]; //store it as-is
  106 + else{ //otherwise throw an error
  107 + std::cerr<<"ERROR stim::filename - drive letter is invalid: "<<unix_dir[0]<<std::endl;
  108 + exit(1);
83 } 109 }
84 - #else  
85 - if(loc[0] == STIM_FILENAME_DIV){  
86 - absolute = true;  
87 - loc = loc.substr(1);  
88 - }  
89 - #endif  
90 -  
91 - //determine the file name  
92 - std::string fname = loc.substr( loc.find_last_of('/') + 1 ); //find the file name (including extension) 110 + unix_dir = unix_dir.substr(2, unix_dir.length()-2); //extract the directory structure
  111 + }
93 112
94 - //parse the file name  
95 - parse_name(fname); 113 + if(unix_dir.front() == '/'){ //if there is a leading slash
  114 + relative = false; //the path is not relative
  115 + unix_dir = unix_dir.substr(1, unix_dir.length() - 1); //remove the slash
  116 + }
  117 + if(unix_dir.back() == '/')
  118 + unix_dir = unix_dir.substr(0, unix_dir.length() - 1);
96 119
97 - //find the directory hierarchy  
98 - std::string dir = loc.substr(0, loc.find_last_of('/') + 1 );  
99 - path = stim::parser::split(dir, '/'); 120 + path = stim::parser::split(unix_dir, '/'); //split up the directory structure
  121 +
  122 + if(relative)
  123 + get_absolute(drive, absolute, path);
  124 + else
  125 + absolute = path;
100 } 126 }
101 127
102 -public: 128 +
103 129
104 - filename(){  
105 - init(); 130 +public:
  131 + filepath(){
  132 + _drive = "";
106 } 133 }
107 134
108 - filename(std::string loc){  
109 - init();  
110 - parse(loc); 135 + filepath(const stim::filepath& p){
  136 + *this = p;
111 } 137 }
112 138
  139 + filepath(const std::string s){
  140 + parse_path(_drive, _path, s);
  141 + }
113 142
114