codebase / stimlib

Commit aef72106a381392f5112f7e94c02b513cdf69f90

Authored by David Mayerich
2 parents 963d0676 ef5cebe5

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

Showing 8 changed files with 323 additions and 225 deletions   Show diff stats
Inline Side-by-side
stim/cuda/cuda_texture.cuh
  Show/Hide comments   View file @aef7210
... ... @@ -99,8 +99,8 @@ namespace stim
99 99 &resource,
100 100 tex,
101 101 target,
102   -// cudaGraphicsMapFlagsReadOnly
103   - cudaGraphicsRegisterFlagsNone
  102 + cudaGraphicsMapFlagsReadOnly
  103 +// cudaGraphicsRegisterFlagsNone
104 104 )
105 105 );
106 106  
... ...
stim/cuda/spider_cost.cuh
  Show/Hide comments   View file @aef7210
... ... @@ -8,6 +8,7 @@
8 8 #include <stim/visualization/colormap.h>
9 9 #include <sstream>
10 10 #include <stim/math/vector.h>
  11 +#include <stim/cuda/cudatools/timer.h>
11 12 #include <stim/cuda/cudatools/devices.h>
12 13 #include <stim/cuda/cudatools/threads.h>
13 14 #include <stim/cuda/cuda_texture.cuh>
... ... @@ -40,11 +41,10 @@ namespace stim{
40 41 ///Returns the value of the template pixel.
41 42 ///@param int x --location of a pixel.
42 43 __device__
43   - float Template(int x)
  44 + float Template(int x, int max_x)
44 45 {
45   - if(x < 32/6 || x > 32*5/6 || (x > 32*2/6 && x < 32*4/6)){
46   -// if(x < 2 || x > 13 || (x > 5 && x < 10)){
47   -// if(x < ceilf(16/6) || x > floorf(16*5/6) || (x > floorf(16*2/6) && x < ceilf(16*4/6))){
  46 + if(x < max_x/6 || x > max_x*5/6 || (x > max_x*2/6 && x < max_x*4/6))
  47 + {
48 48 return 1.0;
49 49 }else{
50 50 return 0.0;
... ... @@ -71,7 +71,7 @@ namespace stim{
71 71 int g_idx = blockIdx.y;
72 72  
73 73 float valIn = tex2D<unsigned char>(texIn, x, y)/255.0;
74   - float valTemp = Template(x);
  74 + float valTemp = Template(x, dx);
75 75  
76 76 // print[idx] = abs(valIn); ///temporary
77 77  
... ... @@ -121,9 +121,18 @@ namespace stim{
121 121 {
122 122  
123 123 //Bind the Texture in GL and allow access to cuda.
  124 +// #ifdef TIMING
  125 +// gpuStartTimer();
  126 +// #endif
124 127 t.MapCudaTexture(texbufferID, texType);
  128 +// #ifdef TIMING
  129 +// std::cout << " " << gpuStopTimer();
  130 +// #endif
125 131  
126 132 //initialize the return arrays.
  133 +// #ifdef TIMING
  134 +// gpuStartTimer();
  135 +// #endif
127 136 float* output;
128 137 output = (float* ) malloc(DIM_Y*sizeof(float));
129 138  
... ... @@ -134,14 +143,26 @@ namespace stim{
134 143 //variables for finding the min.
135 144 float mini = 10000000000000000.0;
136 145 int idx = 0;
  146 +// #ifdef TIMING
  147 +// std::cout << " " << gpuStopTimer();
  148 +// #endif
137 149  
138 150 //cuda launch variables.
  151 +// #ifdef TIMING
  152 +// gpuStartTimer();
  153 +// #endif
139 154 dim3 numBlocks(1, DIM_Y);
140 155 dim3 threadsPerBlock(dx, dy);
141 156  
142   -
143 157 get_diff <<< numBlocks, threadsPerBlock, dx*dy*sizeof(float) >>> (t.getTexture(), result, dx, dy);
144   -
  158 + cudaDeviceSynchronize();
  159 +// #ifdef TIMING
  160 +// std::cout << " " << gpuStopTimer();
  161 +// #endif
  162 +
  163 +// #ifdef TIMING
  164 +// gpuStartTimer();
  165 +// #endif
145 166 HANDLE_ERROR(
146 167 cudaMemcpy(output, result, DIM_Y*sizeof(float), cudaMemcpyDeviceToHost)
147 168 );
... ... @@ -149,10 +170,16 @@ namespace stim{
149 170 for( int i = 0; i<DIM_Y; i++){
150 171 if(output[i] < mini){
151 172 mini = output[i];
152   - idx = i;
  173 + idx = i;
153 174 }
154 175 }
  176 +// #ifdef TIMING
  177 +// std::cout << " " << gpuStopTimer();
  178 +// #endif
155 179  
  180 +// #ifdef TIMING
  181 +// gpuStartTimer();
  182 +// #endif
156 183 // stringstream name; //for debugging
157 184 // name << "Test.bmp";
158 185 // stim::gpu2image<float>(print, name.str(),16,218,0,256);
... ... @@ -162,6 +189,9 @@ namespace stim{
162 189 ret[0] = idx; ret[1] = (int) output[idx];
163 190 // std::cout << "The cost is " << output[idx] << std::endl;
164 191 free(output);
  192 +// #ifdef TIMING
  193 +// std::cout << " " << gpuStopTimer() << std::endl;
  194 +// #endif
165 195 return ret;
166 196 }
167 197  
... ...
stim/cuda/testKernel.cuh
  Show/Hide comments   View file @aef7210
... ... @@ -120,7 +120,7 @@
120 120 cudaDeviceSynchronize();
121 121 stringstream name; //for debugging
122 122 name << nam.c_str();
123   - stim::gpu2image<float>(print, name.str(),x,y,0,255);
  123 + //stim::gpu2image<float>(print, name.str(),x,y,0,255);
124 124  
125 125 tx.UnmapCudaTexture();
126 126 cleanUP();
... ...
stim/gl/gl_spider.h
  Show/Hide comments   View file @aef7210
... ... @@ -7,13 +7,14 @@
7 7 #include <cuda_gl_interop.h>
8 8 #include <cudaGL.h>
9 9 #include <math.h>
10   -#include "stim/gl/gl_texture.h"
11   -#include "stim/visualization/camera.h"
12   -#include "stim/gl/error.h"
13   -#include "stim/math/vector.h"
14   -#include "stim/math/rect.h"
15   -#include "stim/math/matrix.h"
16   -#include "stim/cuda/spider_cost.cuh"
  10 +#include <stim/gl/gl_texture.h>
  11 +#include <stim/visualization/camera.h>
  12 +#include <stim/gl/error.h>
  13 +#include <stim/math/vector.h>
  14 +#include <stim/math/vec3.h>
  15 +#include <stim/math/rect.h>
  16 +#include <stim/math/matrix.h>
  17 +#include <stim/cuda/spider_cost.cuh>
17 18 #include <stim/cuda/cudatools/glbind.h>
18 19 #include <stim/cuda/arraymath.cuh>
19 20 #include <stim/cuda/cudatools.h>
... ... @@ -30,6 +31,11 @@
30 31  
31 32 #include <iostream>
32 33 #include <fstream>
  34 +#ifdef TIMING
  35 + #include <ctime>
  36 + #include <cstdio>
  37 +#endif
  38 +
33 39 #ifdef TESTING
34 40 #include <iostream>
35 41 #include <cstdio>
... ... @@ -47,16 +53,26 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
47 53 //doesn't even need the texture iD really.
48 54 private:
49 55  
  56 + #ifdef TIMING
  57 + double branch_time = 0;
  58 + double direction_time = 0;
  59 + double position_time = 0;
  60 + double size_time = 0;
  61 + double cost_time = 0;
  62 + double network_time = 0;
  63 + double hit_time = 0;
  64 + #endif
  65 +
50 66 //
51   - stim::vec<float> p; //vector designating the position of the spider.
52   - stim::vec<float> d; //vector designating the orientation of the spider
  67 + stim::vec3<float> p; //vector designating the position of the spider.
  68 + stim::vec3<float> d; //vector designating the orientation of the spider
53 69 //always a unit vector.
54 70 stim::vec<float> m; //magnitude of the spider vector.
55 71 //mag[0] = length.
56 72 //mag[1] = width.
57   - std::vector<stim::vec<float> > dV; //A list of all the direction vectors.
58   - std::vector<stim::vec<float> > pV; //A list of all the position vectors.
59   - std::vector<stim::vec<float> > mV; //A list of all the size vectors.
  73 + std::vector<stim::vec3<float> > dV; //A list of all the direction vectors.
  74 + std::vector<stim::vec3<float> > pV; //A list of all the position vectors.
  75 + std::vector<stim::vec3<float> > mV; //A list of all the size vectors.
60 76  
61 77 stim::matrix<float, 4> cT; //current Transformation matrix
62 78 //From tissue space to texture space.
... ... @@ -88,21 +104,21 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
88 104  
89 105  
90 106 //Tracing variables.
91   - std::stack< stim::vec<float> > seeds; //seed positions.
92   - std::stack< stim::vec<float> > seedsvecs; //seed directions.
  107 + std::stack< stim::vec3<float> > seeds; //seed positions.
  108 + std::stack< stim::vec3<float> > seedsvecs; //seed directions.
93 109 std::stack< float > seedsmags; //seed magnitudes.
94 110  
95   - std::vector< stim::vec<float> > cL; //Positions of line currently being traced.
96   - std::vector< stim::vec<float> > cD; //Direction of line currently being traced.
  111 + std::vector< stim::vec3<float> > cL; //Positions of line currently being traced.
  112 + std::vector< stim::vec3<float> > cD; //Direction of line currently being traced.
97 113 std::vector< stim::vec<float> > cM; //Magnitude of line currently being traced.
98 114  
99 115 stim::glnetwork<float> nt; //object for storing the network.
100 116  
101 117 stim::vec<float> rev; //reverse vector;
102 118 stim::camera camSel;
103   - stim::vec<float> ps;
104   - stim::vec<float> ups;
105   - stim::vec<float> ds;
  119 + stim::vec3<float> ps;
  120 + stim::vec3<float> ups;
  121 + stim::vec3<float> ds;
106 122  
107 123 static const float t_length = 16.0;
108 124  
... ... @@ -117,11 +133,19 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
117 133 void
118 134 findOptimalDirection()
119 135 {
  136 + #ifdef TIMING
  137 + gpuStartTimer();
  138 + #endif
120 139 setMatrix(); //create the transformation matrix.
121 140 glCallList(dList); //move the templates to p, d, m.
  141 + glFlush();
  142 + #ifdef TIMING
  143 + direction_time += gpuStopTimer();
  144 + #endif
122 145 #ifdef TESTING
123   - test(texbufferID, GL_TEXTURE_2D,2*t_length,numSamples*t_length, "Final_Cost_Direction.bmp");
  146 +// test(texbufferID, GL_TEXTURE_2D,2*t_length,numSamples*t_length, "Final_Cost_Direction.bmp");
124 147 #endif
  148 +
125 149 int best = getCost(texbufferID,numSamples); //find min cost.
126 150 stim::vec<float> next( //find next vector.
127 151 dV[best][0]*S[0]*R[0],
... ... @@ -142,10 +166,18 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
142 166 void
143 167 findOptimalPosition()
144 168 {
  169 + #ifdef TIMING
  170 + gpuStartTimer();
  171 + #endif
145 172 setMatrix(); //create the transformation matrix.
146 173 glCallList(dList+1); //move the templates to p, d, m.
  174 + glFlush();
  175 + #ifdef TIMING
  176 + position_time += gpuStopTimer();
  177 + #endif
  178 +
147 179 #ifdef TESTING
148   - test(ptexbufferID, GL_TEXTURE_2D,2*t_length, numSamplesPos*t_length, "Final_Cost_Position.bmp");
  180 +// test(ptexbufferID, GL_TEXTURE_2D,2*t_length, numSamplesPos*t_length, "Final_Cost_Position.bmp");
149 181 #endif
150 182 int best = getCost(ptexbufferID, numSamplesPos); //find min cost.
151 183 // std::cerr << best << std::endl;
... ... @@ -168,10 +200,17 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
168 200 void
169 201 findOptimalScale()
170 202 {
  203 + #ifdef TIMING
  204 + gpuStartTimer();
  205 + #endif
171 206 setMatrix(); //create the transformation.
172 207 glCallList(dList+2); //move the templates to p, d, m.
  208 + glFlush();
  209 + #ifdef TIMING
  210 + size_time += gpuStopTimer();
  211 + #endif
173 212 #ifdef TESTING
174   - test(mtexbufferID, GL_TEXTURE_2D, 2*t_length, numSamplesMag*t_length, "Final_Cost_Position.bmp");
  213 +// test(mtexbufferID, GL_TEXTURE_2D, 2*t_length, numSamplesMag*t_length, "Final_Cost_Position.bmp");
175 214 #endif
176 215 int best = getCost(mtexbufferID, numSamplesMag); //get best cost.
177 216 setMagnitude(m[0]*mV[best][0]); //adjust the magnitude.
... ... @@ -188,7 +227,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
188 227 glCallList(dList+3);
189 228 std::vector< stim::vec<float> > result = find_branch(
190 229 btexbufferID, GL_TEXTURE_2D, 16, 216);
191   - stim::vec<float> size(S[0]*R[0], S[1]*R[1], S[2]*R[2]);
  230 + stim::vec3<float> size(S[0]*R[0], S[1]*R[1], S[2]*R[2]);
192 231 if(!result.empty())
193 232 {
194 233 for(int i = 1; i < result.size(); i++)
... ... @@ -200,11 +239,11 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
200 239 1.0);
201 240 cylp = cT*cylp;
202 241  
203   - stim::vec<float> vec(
  242 + stim::vec3<float> vec(
204 243 cylp[0]*S[0]*R[0],
205 244 cylp[1]*S[1]*R[1],
206 245 cylp[2]*S[2]*R[2]);
207   - stim::vec<float> seeddir(-p[0] + cylp[0]*S[0]*R[0],
  246 + stim::vec3<float> seeddir(-p[0] + cylp[0]*S[0]*R[0],
208 247 -p[1] + cylp[1]*S[1]*R[1],
209 248 -p[2] + cylp[2]*S[2]*R[2]);
210 249 seeddir = seeddir.norm();
... ... @@ -230,13 +269,17 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
230 269 void
231 270 branchDetection2(int n = 8, int l_template = 8, int l_square = 8)
232 271 {
  272 + #ifdef TIMING
  273 + gpuStartTimer();
  274 + #endif
  275 +
233 276 if(cL.size() < 4){}
234 277 else{
235 278 setMatrix(1);
236 279 DrawLongCylinder(n, l_template, l_square);
237 280 stim::cylinder<float> cyl(cL, cM);
238 281 std::vector< stim::vec<float> > result = find_branch(btexbufferID, GL_TEXTURE_2D, n*l_square, (cL.size()-1)*l_template);
239   - stim::vec<float> size(S[0]*R[0], S[1]*R[1], S[2]*R[2]);
  282 + stim::vec3<float> size(S[0]*R[0], S[1]*R[1], S[2]*R[2]);
240 283 float pval;
241 284 if(!result.empty())
242 285 {
... ... @@ -257,8 +300,8 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
257 300 {
258 301 pval = (cyl.getl(id)/cyl.getl(cL.size()-1));
259 302 }
260   - stim::vec<float> v = cyl.surf(pval, result[i][0]);
261   - stim::vec<float> di = cyl.p(pval);
  303 + stim::vec3<float> v = cyl.surf(pval, result[i][0]);
  304 + stim::vec3<float> di = cyl.p(pval);
262 305 float rad = cyl.r(pval);
263 306 if(
264 307 !(v[0] > size[0] || v[1] > size[1]
... ... @@ -272,6 +315,9 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
272 315 }
273 316 }
274 317 }
  318 + #ifdef TIMING
  319 + branch_time += gpuStopTimer();
  320 + #endif
275 321 }
276 322  
277 323  
... ... @@ -287,15 +333,15 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
287 333 return cos(2.0*atan(1.0)*u2)*sqrt(-1.0*log(u1));
288 334 }
289 335  
290   - stim::vec<float> uniformRandVector()
  336 + stim::vec3<float> uniformRandVector()
291 337 {
292   - stim::vec<float> r(uniformRandom(), uniformRandom(), 1.0);
  338 + stim::vec3<float> r(uniformRandom(), uniformRandom(), 1.0);
293 339 return r;
294 340 }
295 341  
296   - stim::vec<float> normalRandVector()
  342 + stim::vec3<float> normalRandVector()
297 343 {
298   - stim::vec<float> r(normalRandom(), normalRandom(), 1.0);
  344 + stim::vec3<float> r(normalRandom(), normalRandom(), 1.0);
299 345 return r;
300 346 }
301 347  
... ... @@ -304,58 +350,6 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
304 350 //---------------------TEMPLATE CREATION METHODS----------------------------//
305 351 //--------------------------------------------------------------------------//
306 352  
307   - void
308   - MonteCarloDirectionVectors(int nSamples, float solidAngle = 2*M_PI)
309   - {
310   -///equal area projection.
311   -// std::vector<stim::vec<float> > vecsUni;
312   -// for(int i = 0; i < nSamples; i++)
313   -// {
314   -// stim::vec<float> a = uniformRandVector();
315   -// a[1] = a[1]*2*M_PI;
316   -// a[0] = a[0]*(solidAngle);
317   -// a = a.cyl2cart();
318   -// stim::vec<float> b(
319   -// sqrt(1.-(pow(a[0],2)+(pow(a[1],2)))/4)*a[0],
320   -// sqrt(1.-(pow(a[0],2)+(pow(a[1],2)))/4)*a[1],
321   -// 1.-(pow(a[0],2)+pow(a[1],2))/2);
322   -// 2.*a[0]/(1.+pow(a[0],2)+pow(a[1],2)),
323   -// 2.*a[1]/(1.+pow(a[0],2)+pow(a[1],2)),
324   -// (-1.0+pow(a[0],2)+pow(a[1],2))/(1.0+pow(a[0],2)+pow(a[1],2)));
325   -// b = b.norm();
326   -// b[2] = -b[2]; //this doeswn't want to change
327   -// vecsUni.push_back(b);
328   -// }
329   - float PHI[2], Z[2], range;
330   - PHI[0] = asin(solidAngle/2);
331   - PHI[1] = asin(0);
332   -
333   - Z[0] = cos(PHI[0]);
334   - Z[1] = cos(PHI[1]);
335   -
336   -
337   - range = Z[0] - Z[1];
338   -
339   -
340   - float z, theta, phi;
341   -
342   - std::vector<stim::vec<float> > vecsUni;
343   - for(int i = 0; i < numSamplesPos; i++)
344   - {
345   - stim::vec<float> a(uniformRandom()*0.8, uniformRandom()*0.8, 0.0);
346   - a[0] = a[0]-0.4;
347   - a[1] = a[1]-0.4;
348   - vecsUni.push_back(a);
349   - }
350   -
351   - stringstream name;
352   - for(int i = 0; i < numSamplesPos; i++)
353   - name << vecsUni[i].str() << std::endl;
354   -
355   - std::ofstream outFile;
356   - outFile.open("New_Pos_Vectors.txt");
357   - outFile << name.str().c_str();
358   - }
359 353  
360 354  
361 355 ///@param solidAngle, the size of the arc to sample.
... ... @@ -369,10 +363,10 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
369 363 {
370 364  
371 365 //Set up the vectors necessary for Rectangle creation.
372   - vec<float> Y(1.0,0.0,0.0); //orthogonal vec.
373   - vec<float> pos(0.0,0.0,0.0);
374   - vec<float> mag(1.0, 1.0, 1.0);
375   - vec<float> dir(0.0, 0.0, 1.0);
  366 + stim::vec3<float> Y(1.0,0.0,0.0); //orthogonal vec.
  367 + stim::vec3<float> pos(0.0,0.0,0.0);
  368 + stim::vec3<float> mag(1.0, 1.0, 1.0);
  369 + stim::vec3<float> dir(0.0, 0.0, 1.0);
376 370  
377 371 float PHI[2], Z[2], range;
378 372 PHI[0] = solidAngle/2;
... ... @@ -390,8 +384,8 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
390 384 z = uniformRandom()*range + Z[1];
391 385 theta = uniformRandom()*2*M_PI;
392 386 phi = acos(z);
393   - stim::vec<float> sph(1, theta, phi);
394   - stim::vec<float> cart = sph.sph2cart();
  387 + stim::vec3<float> sph(1, theta, phi);
  388 + stim::vec3<float> cart = sph.sph2cart();
395 389 dV.push_back(cart);
396 390 if(cos(Y.dot(cart)) < 0.087)
397 391 {
... ... @@ -420,16 +414,16 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
420 414 genPositionVectors(float delta = 0.4)
421 415 {
422 416 //Set up the vectors necessary for Rectangle creation.
423   - vec<float> Y(1.0,0.0,0.0); //orthogonal vec.
424   - vec<float> pos(0.0,0.0,0.0);
425   - vec<float> mag(1.0, 1.0, 1.0);
426   - vec<float> dir(0.0, 0.0, 1.0);
  417 + stim::vec3<float> Y(1.0,0.0,0.0); //orthogonal vec.
  418 + stim::vec3<float> pos(0.0,0.0,0.0);
  419 + stim::vec3<float> mag(1.0, 1.0, 1.0);
  420 + stim::vec3<float> dir(0.0, 0.0, 1.0);
427 421  
428 422 //Set up the variable necessary for vector creation.
429 423 glNewList(dList+1, GL_COMPILE);
430 424 for(int i = 0; i < numSamplesPos; i++)
431 425 {
432   - stim::vec<float> temp = uniformRandVector();
  426 + stim::vec3<float> temp = uniformRandVector();
433 427 temp = temp*delta*2.0 - delta/2.0;
434 428 temp[2] = 0.0;
435 429 pV.push_back(temp);
... ... @@ -456,10 +450,10 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
456 450 {
457 451  
458 452 //Set up the vectors necessary for Rectangle creation.
459   - vec<float> Y(1.0, 0.0, 0.0); //orthogonal vec.
460   - vec<float> pos(0.0, 0.0, 0.0);
461   - vec<float> mag(1.0, 1.0, 1.0);
462   - vec<float> dir(0.0, 0.0, 1.0);
  453 + stim::vec3<float> Y(1.0, 0.0, 0.0); //orthogonal vec.
  454 + stim::vec3<float> pos(0.0, 0.0, 0.0);
  455 + stim::vec3<float> mag(1.0, 1.0, 1.0);
  456 + stim::vec3<float> dir(0.0, 0.0, 1.0);
463 457  
464 458 //Set up the variable necessary for vector creation.
465 459 int dim = (sqrt(numSamplesMag)-1)/2;
... ... @@ -467,7 +461,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
467 461 float max = 1.0+delta;
468 462 float step = (max-min)/(numSamplesMag-1);
469 463 float factor;
470   - vec<float> temp(0.0,0.0,0.0);
  464 + stim::vec3<float> temp(0.0,0.0,0.0);
471 465  
472 466 glNewList(dList+2, GL_COMPILE);
473 467 for(int i = 0; i < numSamplesMag; i++){
... ... @@ -495,10 +489,10 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
495 489 void
496 490 UpdateBuffer(float v_x, float v_y)
497 491 {
498   - stim::vec<float>p1;
499   - stim::vec<float>p2;
500   - stim::vec<float>p3;
501   - stim::vec<float>p4;
  492 + stim::vec3<float>p1;
  493 + stim::vec3<float>p2;
  494 + stim::vec3<float>p3;
  495 + stim::vec3<float>p4;
502 496 p1 = hor.p(1,1);
503 497 p2 = hor.p(1,0);
504 498 p3 = hor.p(0,0);
... ... @@ -712,10 +706,6 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
712 706 Bind(GLuint &textureID, GLuint &framebufferID, int nSamples, float len = 8.0)
713 707 {
714 708  
715   -// std::cout << glGetIntegerv(GL_MAX_TEXTURE_SIZE) << std::endl;
716   - GLint max;
717   - glGetIntegerv(GL_MAX_TEXTURE_SIZE, &max);
718   -// std::cout << max << std::endl;
719 709 glBindFramebuffer(GL_FRAMEBUFFER, framebufferID);//set up GL buffer
720 710 glFramebufferTexture2D(
721 711 GL_FRAMEBUFFER,
... ... @@ -727,8 +717,6 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
727 717 GLenum DrawBuffers[1] = {GL_COLOR_ATTACHMENT0};
728 718 glDrawBuffers(1, DrawBuffers);
729 719 glBindTexture(GL_TEXTURE_2D, textureID);
730   -// glClearColor(1,1,1,1);
731   -// glClear(GL_COLOR_BUFFER_BIT);
732 720 glMatrixMode(GL_PROJECTION);
733 721 glLoadIdentity();
734 722 glMatrixMode(GL_MODELVIEW);
... ... @@ -764,9 +752,6 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
764 752 {
765 753 Bind();
766 754 CHECK_OPENGL_ERROR
767   - #ifdef TESTING
768   - start = std::clock();
769   - #endif
770 755 findOptimalDirection();
771 756 findOptimalPosition();
772 757 findOptimalScale();
... ... @@ -774,12 +759,6 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
774 759 Bind(btexbufferID, bfboID, 27);
775 760 branchDetection();
776 761 Unbind();
777   -
778   - #ifdef TESTING
779   - duration_sampling = duration_sampling +
780   - (std::clock() - start) / (double) CLOCKS_PER_SEC;
781   - num_sampling = num_sampling + 1.0;
782   - #endif
783 762 return current_cost;
784 763 }
785 764  
... ... @@ -797,17 +776,9 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
797 776 int
798 777 getCost()
799 778 {
800   - #ifdef TESTING
801   - start = std::clock();
802   - #endif
803 779 stim::vec<int> cost =
804 780 stim::cuda::get_cost(texbufferID, GL_TEXTURE_2D, numSamples);
805 781 cudaDeviceSynchronize();
806   - #ifdef TESTING
807   - duration_cuda = duration_cuda +
808   - (std::clock() - start) / (double) CLOCKS_PER_SEC;
809   - num_cuda = num_cuda + 1.0;
810   - #endif
811 782 current_cost = cost[1];
812 783 return cost[0];
813 784 }
... ... @@ -815,16 +786,13 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
815 786 int
816 787 getCost(GLuint tID, int n)
817 788 {
818   - #ifdef TESTING
819   - start = std::clock();
  789 + #ifdef TIMING
  790 + gpuStartTimer();
820 791 #endif
821 792 stim::vec<int> cost =
822 793 stim::cuda::get_cost(tID, GL_TEXTURE_2D, n, 2*t_length, t_length);
823   - cudaDeviceSynchronize();
824   - #ifdef TESTING
825   - duration_cuda = duration_cuda +
826   - (std::clock() - start) / (double) CLOCKS_PER_SEC;
827   - num_cuda = num_cuda + 1.0;
  794 + #ifdef TIMING
  795 + cost_time += gpuStopTimer();
828 796 #endif
829 797 current_cost = cost[1];
830 798 return cost[0];
... ... @@ -837,10 +805,6 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
837 805 initCuda()
838 806 {
839 807 stim::cudaSetDevice();
840   - MonteCarloDirectionVectors(500);
841   - //GLint max;
842   - //glGetIntegerv(GL_MAX_TEXTURE_SIZE, &max);
843   - //std::cout << max << std::endl;
844 808 }
845 809  
846 810 //horizonal rectangle forming the spider.
... ... @@ -848,13 +812,9 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
848 812 //vectical rectangle forming the spider.
849 813 stim::rect<float> ver;
850 814  
851   - //Testing and Timing variables.
  815 + //Timing variables.
852 816 #ifdef TESTING
853 817 std::clock_t start;
854   - double duration_sampling = 0.0;
855   - double duration_cuda = 0.0;
856   - double num_sampling = 0.0;
857   - double num_cuda = 0.0;
858 818 #endif
859 819  
860 820 //--------------------------------------------------------------------------//
... ... @@ -869,11 +829,11 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
869 829 (int samples = 1089, int samplespos = 441,int samplesmag = 144)
870 830 {
871 831 // std::cout << "I ran this constructor" << std::endl;
872   - p = vec<float>(0.0, 0.0, 0.0);
873   - d = vec<float>(0.0, 0.0, 1.0);
874   - m = vec<float>(1.0, 1.0);
875   - S = vec<float>(1.0, 1.0, 1.0);
876   - R = vec<float>(1.0, 1.0, 1.0);
  832 + p = stim::vec3<float>(0.0, 0.0, 0.0);
  833 + d = stim::vec3<float>(0.0, 0.0, 1.0);
  834 + m = stim::vec<float>(1.0, 1.0);
  835 + S = stim::vec3<float>(1.0, 1.0, 1.0);
  836 + R = stim::vec3<float>(1.0, 1.0, 1.0);
877 837 // std::cout << samples << std::endl;
878 838 numSamples = samples;
879 839 // std::cout << numSamples << std::endl;
... ... @@ -894,11 +854,11 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
894 854 (float pos_x, float pos_y, float pos_z, float dir_x, float dir_y, float dir_z,
895 855 float mag_x, int numsamples = 1089, int numsamplespos = 441, int numsamplesmag =144)
896 856 {
897   - p = vec<float>(pos_x, pos_y, pos_z);
898   - d = vec<float>(dir_x, dir_y, dir_z);
899   - m = vec<float>(mag_x, mag_x, mag_x);
900   - S = vec<float>(1.0,1.0,1.0);
901   - R = vec<float>(1.0,1.0,1.0);
  857 + p = stim::vec3<float>(pos_x, pos_y, pos_z);
  858 + d = stim::vec3<float>(dir_x, dir_y, dir_z);
  859 + m = stim::vec<float>(mag_x, mag_x, mag_x);
  860 + S = stim::vec3<float>(1.0,1.0,1.0);
  861 + R = stim::vec3<float>(1.0,1.0,1.0);
902 862 numSamples = numsamples;
903 863 numSamplesPos = numsamplespos;
904 864 numSamplesMag = numsamplesmag;
... ... @@ -910,13 +870,13 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
910 870 ///@param float mag, size of the vector.
911 871 ///@param int samples, number of templates this spider is going to use.
912 872 gl_spider
913   - (stim::vec<float> pos, stim::vec<float> dir, float mag, int samples = 1089, int samplesPos = 441, int samplesMag = 144)
  873 + (stim::vec3<float> pos, stim::vec3<float> dir, float mag, int samples = 1089, int samplesPos = 441, int samplesMag = 144)
914 874 {
915 875 p = pos;
916 876 d = dir;
917 877 m = vec<float>(mag, mag, mag);
918   - S = vec<float>(1.0,1.0,1.0);
919   - R = vec<float>(1.0,1.0,1.0);
  878 + S = vec3<float>(1.0,1.0,1.0);
  879 + R = vec3<float>(1.0,1.0,1.0);
920 880 numSamples = samples;
921 881 numSamplesPos = samplesPos;
922 882 numSamplesMag = samplesMag;
... ... @@ -948,10 +908,13 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
948 908 texID = id;
949 909 //GenerateFBO(16, numSamples*8);
950 910 GenerateFBO(t_length*2, numSamples*t_length, texbufferID, fboID);
  911 + std::cout << numSamples << std::endl;
951 912 CHECK_OPENGL_ERROR
952 913 GenerateFBO(t_length*2, numSamplesPos*t_length, ptexbufferID, pfboID);
  914 + std::cout << numSamplesPos << std::endl;
953 915 CHECK_OPENGL_ERROR
954 916 GenerateFBO(t_length*2, numSamplesMag*t_length, mtexbufferID, mfboID);
  917 + std::cout << numSamplesMag << std::endl;
955 918 CHECK_OPENGL_ERROR
956 919 GenerateFBO(16, 216, btexbufferID, bfboID);
957 920 CHECK_OPENGL_ERROR
... ... @@ -978,21 +941,21 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
978 941 //-----------------------------ACCESS METHODS-------------------------------//
979 942 //--------------------------------------------------------------------------//
980 943 ///Returns the p vector.
981   - vec<float>
  944 + vec3<float>
982 945 getPosition()
983 946 {
984 947 return p;
985 948 }
986 949  
987 950 ///Returns the d vector.
988   - vec<float>
  951 + vec3<float>
989 952 getDirection()
990 953 {
991 954 return d;
992 955 }
993 956  
994 957 ///Returns the m vector.
995   - vec<float>
  958 + stim::vec<float>
996 959 getMagnitude()
997 960 {
998 961 return m;
... ... @@ -1001,7 +964,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1001 964 ///@param stim::vec<float> pos, the new p.
1002 965 ///Sets the p vector to input vector pos.
1003 966 void
1004   - setPosition(vec<float> pos)
  967 + setPosition(stim::vec3<float> pos)
1005 968 {
1006 969 p = pos;
1007 970 }
... ... @@ -1021,7 +984,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1021 984 ///@param stim::vec<float> dir, the new d.
1022 985 ///Sets the d vector to input vector dir.
1023 986 void
1024   - setDirection(vec<float> dir)
  987 + setDirection(stim::vec3<float> dir)
1025 988 {
1026 989 d = dir;
1027 990 }
... ... @@ -1041,7 +1004,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1041 1004 ///@param stim::vec<float> dir, the new d.
1042 1005 ///Sets the m vector to the input vector mag.
1043 1006 void
1044   - setMagnitude(vec<float> mag)
  1007 + setMagnitude(stim::vec<float> mag)
1045 1008 {
1046 1009 m[0] = mag[0];
1047 1010 m[1] = mag[0];
... ... @@ -1071,7 +1034,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1071 1034 ///@param stim::vec<float> Dims, voxel size.
1072 1035 ///Sets the voxel sizes in each direction. necessary for non-standard data.
1073 1036 void
1074   - setDims(stim::vec<float> Dims)
  1037 + setDims(stim::vec3<float> Dims)
1075 1038 {
1076 1039 S = Dims;
1077 1040 }
... ... @@ -1091,7 +1054,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1091 1054 ///@param stim::vec<float> Dims, size of the volume.
1092 1055 ///Sets the data volume sizes in each direction.
1093 1056 void
1094   - setSize(stim::vec<float> Siz)
  1057 + setSize(stim::vec3<float> Siz)
1095 1058 {
1096 1059 S = Siz;
1097 1060 }
... ... @@ -1101,7 +1064,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1101 1064 ///rotates from 0.0, 0.0, 1.0 to dir.
1102 1065 ///return is in degrees for use with glRotatef.
1103 1066 stim::vec<float>
1104   - getRotation(stim::vec<float> dir)
  1067 + getRotation(stim::vec3<float> dir)
1105 1068 {
1106 1069 stim::vec<float> out(0.0,0.0,0.0,0.0);
1107 1070 stim::vec<float> from(0.0,0.0,1.0);
... ... @@ -1113,7 +1076,8 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1113 1076 out[3] = 1.0;
1114 1077 } else {
1115 1078 stim::vec<float> temp(0.0, 0.0, 0.0);;
1116   - temp = (from.cross(dir)).norm();
  1079 + stim::vec<float> dir1(dir[0], dir[1], dir[2]);
  1080 + temp = (from.cross(dir1)).norm();
1117 1081 out[1] = temp[0];
1118 1082 out[2] = temp[1];
1119 1083 out[3] = temp[2];
... ... @@ -1125,7 +1089,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1125 1089 ///Adds a seed to the seed list.
1126 1090 ///Assumes that the coordinates passes are in tissue space.
1127 1091 void
1128   - setSeed(stim::vec<float> pos)
  1092 + setSeed(stim::vec3<float> pos)
1129 1093 {
1130 1094 seeds.push(pos);
1131 1095 }
... ... @@ -1133,7 +1097,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1133 1097 ///@param stim::vec<float> dir, the direction of the seed to be added.
1134 1098 ///Adds a seed to the seed directions list.
1135 1099 void
1136   - setSeedVec(stim::vec<float> dir)
  1100 + setSeedVec(stim::vec3<float> dir)
1137 1101 {
1138 1102 seedsvecs.push(dir);
1139 1103 }
... ... @@ -1173,18 +1137,18 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1173 1137 }
1174 1138  
1175 1139 ///Method to get the top of the seed positions stack.
1176   - stim::vec<float>
  1140 + stim::vec3<float>
1177 1141 getLastSeed()
1178 1142 {
1179   - stim::vec<float> tp = seeds.top();
  1143 + stim::vec3<float> tp = seeds.top();
1180 1144 return tp;
1181 1145 }
1182 1146  
1183 1147 ///Method to get the top of the seed direction stack.
1184   - stim::vec<float>
  1148 + stim::vec3<float>
1185 1149 getLastSeedVec()
1186 1150 {
1187   - stim::vec<float> tp = seedsvecs.top();
  1151 + stim::vec3<float> tp = seedsvecs.top();
1188 1152 return tp;
1189 1153 }
1190 1154  
... ... @@ -1206,12 +1170,53 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1206 1170 }
1207 1171  
1208 1172 ///returns the seeds position stack.
1209   - std::stack<stim::vec<float> >
  1173 + std::stack<stim::vec3<float> >
1210 1174 getSeeds()
1211 1175 {
1212 1176 return seeds;
1213 1177 }
1214 1178  
  1179 + ///sets the number of direction templates.
  1180 + void
  1181 + setNumberDirectionTemplates(int n)
  1182 + {
  1183 + numSamples = n;
  1184 + }
  1185 +
  1186 + ///sets the number of position templates.
  1187 + void
  1188 + setNumberPositionTemplates(int n)
  1189 + {
  1190 + numSamplesPos = n;
  1191 + }
  1192 +
  1193 + ///sets the number of position templates.
  1194 + void
  1195 + setNumberSizeTemplates(int n)
  1196 + {
  1197 + numSamplesMag = n;
  1198 + }
  1199 +
  1200 + #ifdef TIMING
  1201 + ///Returns the timings at the moment the method is called.
  1202 + ///In the following order: Branch, Direction, Position, Size, Cost, Network, Hit_detetion.
  1203 + std::vector<double>
  1204 + getTimings()
  1205 + {
  1206 + std::vector <double> ret;
  1207 + ret.resize(7);
  1208 + ret[0] = branch_time;
  1209 + ret[1] = direction_time;
  1210 + ret[2] = position_time;
  1211 + ret[3] = size_time;
  1212 + ret[4] = cost_time;
  1213 + ret[5] = network_time;
  1214 + ret[6] = hit_time;
  1215 +
  1216 + return ret;
  1217 + }
  1218 + #endif
  1219 +
1215 1220 ///returns true if all seed stacks are empty, else false.
1216 1221 bool
1217 1222 Empty()
... ... @@ -1252,7 +1257,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1252 1257 for(int i = 0; i < nt.sizeE(); i++)
1253 1258 {
1254 1259 std::vector<stim::vec< float > > cm = nt.getEdgeCenterLineMag(i);
1255   - std::vector<stim::vec< float > > ce = nt.getEdgeCenterLine(i);
  1260 + std::vector<stim::vec3< float > > ce = nt.getEdgeCenterLine(i);
1256 1261 sk.Begin(stim::OBJ_LINE);
1257 1262 for(int j = 0; j < ce.size(); j++)
1258 1263 {
... ... @@ -1269,7 +1274,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1269 1274 stim::glObj<float>
1270 1275 getNetwork()
1271 1276 {
1272   -// return sk;
  1277 +
1273 1278 }
1274 1279  
1275 1280 ///returns a COPY of the entire stim::glnetwork object.
... ... @@ -1295,7 +1300,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1295 1300 void
1296 1301 Update()
1297 1302 {
1298   - vec<float> Y(1.0,0.0,0.0);
  1303 + vec3<float> Y(1.0,0.0,0.0);
1299 1304 if(cos(Y.dot(d))< 0.087){
1300 1305 Y[0] = 0.0; Y[1] = 1.0;}
1301 1306 hor = stim::rect<float>(m, p, d.norm(),
... ... @@ -1310,9 +1315,6 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1310 1315 {
1311 1316 Bind(texbufferID, fboID, numSamples, t_length);
1312 1317 CHECK_OPENGL_ERROR
1313   - #ifdef TESTING
1314   - start = std::clock();
1315   - #endif
1316 1318 findOptimalDirection();
1317 1319 Unbind();
1318 1320 Bind(ptexbufferID, pfboID, numSamplesPos, t_length);
... ... @@ -1323,12 +1325,6 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1323 1325 Unbind();
1324 1326 CHECK_OPENGL_ERROR
1325 1327  
1326   - #ifdef TESTING
1327   - duration_sampling = duration_sampling +
1328   - (std::clock() - start) / (double) CLOCKS_PER_SEC;
1329   - num_sampling = num_sampling + 1.0;
1330   - #endif
1331   - std::cout << current_cost << std::endl;
1332 1328 return current_cost;
1333 1329 }
1334 1330  
... ... @@ -1345,6 +1341,39 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1345 1341 //--------------------------------------------------------------------------//
1346 1342  
1347 1343 void
  1344 + MonteCarloDirectionVectors(int nSamples, float solidAngle = 2*M_PI)
  1345 + {
  1346 + float PHI[2], Z[2], range;
  1347 + PHI[0] = asin(solidAngle/2);
  1348 + PHI[1] = asin(0);
  1349 +
  1350 + Z[0] = cos(PHI[0]);
  1351 + Z[1] = cos(PHI[1]);
  1352 +
  1353 +
  1354 + range = Z[0] - Z[1];
  1355 +
  1356 +
  1357 + float z, theta, phi;
  1358 +
  1359 + std::vector<stim::vec3<float> > vecsUni;
  1360 + for(int i = 0; i < numSamplesPos; i++)
  1361 + {
  1362 + stim::vec3<float> a(uniformRandom()*0.8, uniformRandom()*0.8, 0.0);
  1363 + a[0] = a[0]-0.4;
  1364 + a[1] = a[1]-0.4;
  1365 + vecsUni.push_back(a);
  1366 + }
  1367 +
  1368 + stringstream name;
  1369 + for(int i = 0; i < numSamplesPos; i++)
  1370 + name << vecsUni[i].str() << std::endl;
  1371 +
  1372 + std::ofstream outFile;
  1373 + outFile.open("New_Pos_Vectors.txt");
  1374 + outFile << name.str().c_str();
  1375 + }
  1376 + void
1348 1377 DrawCylinder()
1349 1378 {
1350 1379 glNewList(dList+3, GL_COMPILE);
... ... @@ -1389,7 +1418,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1389 1418 GenerateFBO(n*l_square, cylLen*l_template, btexbufferID, bfboID);
1390 1419 Bind(btexbufferID, bfboID, cylLen, l_template*l_square/2.0);
1391 1420 stim::cylinder<float> cyl(cL, cM);
1392   - std::vector<std::vector<stim::vec<float> > > p = cyl.getPoints(n);
  1421 + std::vector<std::vector<stim::vec3<float> > > p = cyl.getPoints(n);
1393 1422 for(int i = 0; i < p.size()-1; i++) ///number of circles
1394 1423 {
1395 1424 for(int j = 0; j < p[0].size()-1; j++) ///points in the circle
... ... @@ -1418,11 +1447,10 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1418 1447 void
1419 1448 trace(int min_cost)
1420 1449 {
1421   -// rev = stim::vec<float>(0.0,0.0,1.0);
1422 1450 bool sEmpty = true;
1423 1451 float lastmag = 16.0;;
1424   - stim::vec<float> curSeed;
1425   - stim::vec<float> curSeedVec;
  1452 + stim::vec3<float> curSeed;
  1453 + stim::vec3<float> curSeedVec;
1426 1454 float curSeedMag;
1427 1455 while(!Empty())
1428 1456 {
... ... @@ -1448,9 +1476,13 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1448 1476 }
1449 1477  
1450 1478 int
1451   - selectObject(stim::vec<float> loc, stim::vec<float> dir, float mag)
  1479 + selectObject(stim::vec3<float> loc, stim::vec3<float> dir, float mag)
1452 1480 {
1453 1481 //Define the varibles and turn on Selection Mode
  1482 +
  1483 + #ifdef TIMING
  1484 + gpuStartTimer();
  1485 + #endif
1454 1486  
1455 1487 float s = 3.0;
1456 1488 GLuint selectBuf[2048];
... ... @@ -1510,6 +1542,10 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1510 1542 // glEnable(GL_CULL_FACE);
1511 1543 hits = glRenderMode(GL_RENDER);
1512 1544 int found_hits = processHits(hits, selectBuf);
  1545 + #ifdef TIMING
  1546 + hit_time += gpuStopTimer();
  1547 + #endif
  1548 +
1513 1549 return found_hits;
1514 1550 }
1515 1551  
... ... @@ -1554,10 +1590,14 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1554 1590 }
1555 1591  
1556 1592 void
1557   - addToNetwork(pair<stim::fiber<float>, int> in, stim::vec<float> spos,
1558   - stim::vec<float> smag, stim::vec<float> sdir)
  1593 + addToNetwork(pair<stim::fiber<float>, int> in, stim::vec3<float> spos,
  1594 + stim::vec<float> smag, stim::vec3<float> sdir)
1559 1595 {
1560   - std::vector<stim::vec<float> > ce = in.first.centerline();
  1596 + #ifdef TIMING
  1597 + double s = std::clock();
  1598 + #endif
  1599 +
  1600 + std::vector<stim::vec3<float> > ce = in.first.centerline();
1561 1601 std::vector<stim::vec<float> > cm = in.first.centerlinemag();
1562 1602 //if the fiber is longer than 2 steps (the number it takes to diverge)
1563 1603 if(ce.size() > 3)
... ... @@ -1590,6 +1630,11 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1590 1630 } else { nt.addEdge(ce,cm, -1, -1);}
1591 1631 }
1592 1632 }
  1633 +
  1634 + #ifdef TIMING
  1635 + double network_time = (std::clock() - s) / (double) CLOCKS_PER_SEC;
  1636 + network_time += network_time * 1000.0;
  1637 + #endif
1593 1638 }
1594 1639  
1595 1640  
... ... @@ -1602,12 +1647,12 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1602 1647 }
1603 1648  
1604 1649 std::pair<stim::fiber<float>, int >
1605   - traceLine(stim::vec<float> pos, stim::vec<float> mag, int min_cost)
  1650 + traceLine(stim::vec3<float> pos, stim::vec<float> mag, int min_cost)
1606 1651 {
1607 1652 //starting (seed) position and magnitude.
1608   - stim::vec<float> spos = getPosition();
  1653 + stim::vec3<float> spos = getPosition();
1609 1654 stim::vec<float> smag = getMagnitude();
1610   - stim::vec<float> sdir = getDirection();
  1655 + stim::vec3<float> sdir = getDirection();
1611 1656  
1612 1657 Bind();
1613 1658 // sk.Begin(stim::OBJ_LINE);
... ... @@ -1624,7 +1669,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1624 1669 int h;
1625 1670 bool started = false;
1626 1671 bool running = true;
1627   - stim::vec<float> size(S[0]*R[0], S[1]*R[1], S[2]*R[2]);
  1672 + stim::vec3<float> size(S[0]*R[0], S[1]*R[1], S[2]*R[2]);
1628 1673 while(running)
1629 1674 {
1630 1675 int cost = Step();
... ... @@ -1684,7 +1729,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1684 1729 break;
1685 1730 }
1686 1731 else {
1687   - cL.push_back(stim::vec<float>(p[0], p[1],p[2]));
  1732 + cL.push_back(stim::vec3<float>(p[0], p[1],p[2]));
1688 1733 cM.push_back(stim::vec<float>(m[0], m[0]));
1689 1734 // Bind(btexbufferID, bfboID, 27);
1690 1735 Unbind();
... ...
stim/image/image.h
  Show/Hide comments   View file @aef7210
... ... @@ -193,6 +193,7 @@ public:
193 193 get_interleaved_bgr(buffer);
194 194 cv::Mat cvImage((int)Y(), (int)X(), cv_type(), buffer);
195 195 cv::imwrite(filename, cvImage);
  196 + free(buffer);
196 197 }
197 198  
198 199 void set_interleaved(T* buffer, size_t width, size_t height, size_t channels){
... ...
stim/math/vec3.h
  Show/Hide comments   View file @aef7210
... ... @@ -37,10 +37,14 @@ public:
37 37 }
38 38  
39 39 //access an element using an index
40   - CUDA_CALLABLE T& operator[](int idx){
  40 + CUDA_CALLABLE T& operator[](size_t idx){
41 41 return ptr[idx];
42 42 }
43 43  
  44 + CUDA_CALLABLE T* data(){
  45 + return ptr;
  46 + }
  47 +
44 48 /// Casting operator. Creates a new vector with a new type U.
45 49 template< typename U >
46 50 CUDA_CALLABLE operator vec3<U>(){
... ... @@ -249,4 +253,4 @@ std::ostream&amp; operator&lt;&lt;(std::ostream&amp; os, stim::vec3&lt;T&gt; const&amp; rhs){
249 253 return os;
250 254 }
251 255  
252   -#endif
253 256 \ No newline at end of file
  257 +#endif
... ...
stim/math/vector.h
  Show/Hide comments   View file @aef7210
... ... @@ -74,6 +74,13 @@ struct vec : public std::vector&lt;T&gt;
74 74 at(i) = other[i];
75 75 }
76 76 }
  77 +
  78 +// vec( vec3<T>& other){
  79 +// resize(3); //resize the current vector to match the copy
  80 +// for(size_t i=0; i<3; i++){ //copy each element
  81 +// at(i) = other[i];
  82 +// }
  83 +// }
77 84  
78 85 //I'm not sure what these were doing here.
79 86 //Keep them now, we'll worry about it later.
... ... @@ -328,15 +335,16 @@ struct vec : public std::vector&lt;T&gt;
328 335 return *this;
329 336 }
330 337  
331   - /// Cast to a vec3
332   - operator stim::vec3<T>(){
333   - stim::vec3<T> r;
334   - size_t N = std::min<size_t>(size(), 3);
335   - for(size_t i = 0; i < N; i++)
336   - r[i] = at(i);
337   - return r;
338   - }
339   -
  338 + /// Cast to a vec3
  339 + operator stim::vec3<T>(){
  340 + stim::vec3<T> r;
  341 + size_t N = std::min<size_t>(size(), 3);
  342 + for(size_t i = 0; i < N; i++)
  343 + r[i] = at(i);
  344 + return r;
  345 + }
  346 +
  347 +
340 348 /// Casting and assignment
341 349 template<typename Y>
342 350 vec<T> & operator=(vec<Y> rhs){
... ... @@ -347,6 +355,16 @@ struct vec : public std::vector&lt;T&gt;
347 355 at(i) = rhs[i];
348 356 return *this;
349 357 }
  358 +
  359 + /// Assign a vec = vec3
  360 + template<typename Y>
  361 + vec<T> & operator=(vec3<Y> rhs)
  362 + {
  363 + resize(3);
  364 + for(size_t i=0; i<3; i++)
  365 + at(i) = rhs[i];
  366 + return *this;
  367 + }
350 368  
351 369 /// Unary minus (returns the negative of the vector)
352 370 vec<T> operator-() const{
... ...
stim/visualization/cylinder.h
  Show/Hide comments   View file @aef7210
... ... @@ -23,7 +23,7 @@ class cylinder
23 23  
24 24 }
25 25  
26   - ///inits the cylinder from a list of points (inP) and radii (inM)
  26 + ///inits the cylinder from a list of points (std::vector of stim::vec3 --inP) and radii (inM)
27 27 void
28 28 init(std::vector<stim::vec3<T> > inP, std::vector<stim::vec<T> > inM)
29 29 {
... ... @@ -142,14 +142,15 @@ class cylinder
142 142 }
143 143  
144 144 ///constructor to create a cylinder from a set of points, radii, and the number of sides for the cylinder.
145   - ///@param inP: Vector of stim vecs composing the points of the centerline.
  145 + ///@param inP: Vector of stim vec3 composing the points of the centerline.
146 146 ///@param inM: Vector of stim vecs composing the radii of the centerline.
147   - cylinder(std::vector<stim::vec3<T> > inP, std::vector<stim::vec3<T> > inM){
  147 + cylinder(std::vector<stim::vec3<T> > inP, std::vector<stim::vec<T> > inM){
148 148 init(inP, inM);
149 149 }
150 150  
  151 +
151 152 ///Constructor defines a cylinder with centerline inP and magnitudes of zero
152   - ///@param inP: Vector of stim vecs composing the points of the centerline
  153 + ///@param inP: Vector of stim vec3 composing the points of the centerline
153 154 cylinder(std::vector< stim::vec3<T> > inP){
154 155 std::vector< stim::vec<T> > inM; //create an array of arbitrary magnitudes
155 156  
... ... @@ -160,7 +161,6 @@ class cylinder
160 161 init(inP, inM);
161 162 }
162 163  
163   -
164 164 ///Returns the number of points on the cylinder centerline
165 165  
166 166 unsigned int size(){
... ...