codebase / stimlib

Commit 4166e97375f81aae2bdd5b74c039317b6f1f930a

Authored by Pavel Govyadinov
1 parent 47fe6300

Added debug tags to everything that compile only when the debug option is done, … 

…fixed the instability as well as some other issue that came out later
Showing 2 changed files with 117 additions and 82 deletions   Show diff stats
Inline Side-by-side
stim/cuda/testKernel.cuh
  Show/Hide comments   View file @4166e97
... ... @@ -8,7 +8,7 @@
8 8 #include <stim/cuda/cudatools/devices.h>
9 9 #include <stim/cuda/cudatools/threads.h>
10 10 #include <stim/cuda/cuda_texture.cuh>
11   - stim::cuda::cuda_texture tx; //texture object.
  11 +
12 12 float* print;
13 13  
14 14 ///Initialization function, allocates the memory and passes the necessary
... ... @@ -26,16 +26,16 @@
26 26 cudaFree(print); ///temporary
27 27 }
28 28  
29   - __device__
30   - float templ(int x)
31   - {
32   - if(x < 32/6 || x > 32*5/6 || (x > 32*2/6 && x < 32*4/6)){
33   - return 1.0;
34   - }else{
35   - return 0.0;
36   - }
37   -
38   - }
  29 + __device__
  30 + float templ(int x, int max_x)
  31 + {
  32 + if(x < max_x/6 || x > max_x*5/6 || (x > max_x*2/6 && x < max_x*4/6))
  33 + {
  34 + return 1.0;
  35 + }else{
  36 + return 0.0;
  37 + }
  38 + }
39 39  
40 40 ///Find the difference of the given set of samples and the template
41 41 ///using cuda acceleration.
... ... @@ -44,33 +44,24 @@
44 44 ///@param float* result --a pointer to the memory that stores the result.
45 45 __global__
46 46 //void get_diff (float *result)
47   - void get_diff (cudaTextureObject_t texIn, float *print)
  47 + void get_diff (cudaTextureObject_t texIn, float *print, int dx)
48 48 {
49 49 int x = threadIdx.x + blockIdx.x * blockDim.x;
50 50 int y = threadIdx.y + blockIdx.y * blockDim.y;
51   -// int idx = y*64+x;
52   - int idx = y*32+x;
  51 + int idx = y*dx+x;
53 52 // int idx = y*16+x;
54 53  
55 54 float valIn = tex2D<unsigned char>(texIn, x, y);
56   - float templa = templ(x);
57   - print[idx] = valIn; ///temporary
  55 + float templa = templ(x, 32)*255.0;
  56 + print[idx] = abs(valIn-templa); ///temporary
58 57 //print[idx] = abs(templa); ///temporary
59 58  
60 59 }
61 60  
62   -
63   - ///External access-point to the cuda function
64   - ///@param GLuint texbufferID --GLtexture (most be contained in a framebuffer object)
65   - /// that holds the data that will be handed to cuda.
66   - ///@param GLenum texType --either GL_TEXTURE_1D, GL_TEXTURE_2D or GL_TEXTURE_3D
67   - /// may work with other gl texture types, but untested.
68   - ///@param DIM_Y, the number of samples in the template.
69   - void test(GLint texbufferID, GLenum texType, int x, int y)
  61 + void test(cudaTextureObject_t tObj, int x, int y, std::string nam)
70 62 {
71 63  
72 64 //Bind the Texture in GL and allow access to cuda.
73   - tx.MapCudaTexture(texbufferID, texType);
74 65  
75 66 //initialize the return arrays.
76 67  
... ... @@ -85,44 +76,13 @@
85 76  
86 77  
87 78 // get_diff <<< blocks, threads >>> (tx.getTexture(), print);
88   - get_diff <<< numBlocks, threadsPerBlock >>> (tx.getTexture(), print);
89   -
90   - cudaDeviceSynchronize();
91   - stringstream name; //for debugging
92   - name << "FromTex.bmp";
93   - stim::gpu2image<float>(print, name.str(),x,y,0,255);
94   -
95   - tx.UnmapCudaTexture();
96   - cleanUP();
97   - }
98   -
99   - void test(GLint texbufferID, GLenum texType, int x, int y, std::string nam)
100   - {
101   -
102   - //Bind the Texture in GL and allow access to cuda.
103   - tx.MapCudaTexture(texbufferID, texType);
104   -
105   - //initialize the return arrays.
106   -
107   - initArray(x,y);
108   - dim3 numBlocks(1, y);
109   - dim3 threadsPerBlock(x, 1);
110   - int max_threads = stim::maxThreadsPerBlock();
111   - //dim3 threads(max_threads, 1);
112   - //dim3 blocks(x / threads.x + 1, y);
113   - //dim3 numBlocks(2, 2);
114   - //dim3 threadsPerBlock(8, 108);
115   -
116   -
117   -// get_diff <<< blocks, threads >>> (tx.getTexture(), print);
118   - get_diff <<< numBlocks, threadsPerBlock >>> (tx.getTexture(), print);
  79 + get_diff <<< numBlocks, threadsPerBlock >>> (tObj, print, x);
119 80  
120 81 cudaDeviceSynchronize();
121 82 stringstream name; //for debugging
122 83 name << nam.c_str();
123   - //stim::gpu2image<float>(print, name.str(),x,y,0,255);
  84 + stim::gpu2image<float>(print, name.str(),x,y,0,255);
124 85  
125   - tx.UnmapCudaTexture();
126 86 cleanUP();
127 87 }
128 88  
... ...
stim/gl/gl_spider.h
  Show/Hide comments   View file @4166e97
... ... @@ -134,7 +134,13 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
134 134 stim::cuda::cuda_texture t_mag; //cuda_texture object used as an interface between OpenGL and cuda for size vectors.
135 135  
136 136  
137   - int iter;
  137 + #ifdef DEBUG
  138 + stringstream name;
  139 + int iter;
  140 + int iter_pos;
  141 + int iter_dir;
  142 + int iter_siz;
  143 + #endif
138 144  
139 145 //--------------------------------------------------------------------------//
140 146 //-------------------------------PRIVATE METHODS----------------------------//
... ... @@ -155,11 +161,14 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
155 161 #ifdef TIMING
156 162 direction_time += gpuStopTimer(); //profiling
157 163 #endif
158   - #ifdef TESTING
159   -// test(texbufferID, GL_TEXTURE_2D,2*n_pixels,numSamples*n_pixels, "Final_Cost_Direction.bmp");
160   - #endif
161 164  
162 165 int best = getCost(t_dir.getTexture(), t_dir.getAuxArray() ,numSamples); //find min cost.
  166 + #ifdef DEBUG
  167 + name.str("");
  168 + name << "Final_Cost_Direction_fiber_"<< iter << "_" << iter_dir << ".bmp";
  169 + test(t_dir.getTexture(), n_pixels*2.0, numSamples*n_pixels, name.str());
  170 + iter_dir++;
  171 + #endif
163 172 stim::vec<float> next( ///calculate the next vector.
164 173 dV[best][0]*S[0]*R[0],
165 174 dV[best][1]*S[1]*R[1],
... ... @@ -189,10 +198,13 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
189 198 position_time += gpuStopTimer(); ///timer for profiling
190 199 #endif
191 200  
192   - #ifdef TESTING
193   -// test(texbufferID, GL_TEXTURE_2D,2*n_pixels,numSamples*n_pixels, "Final_Cost_Direction.bmp");
194   - #endif
195 201 int best = getCost(t_pos.getTexture(), t_pos.getAuxArray(), numSamplesPos); //find min cost.
  202 + #ifdef DEBUG
  203 + name.str("");
  204 + name << "Final_Cost_Position_" << iter << "_" << iter_pos << ".bmp";
  205 + test(t_pos.getTexture(), n_pixels*2.0, numSamplesPos*n_pixels, name.str());
  206 + iter_pos++;
  207 + #endif
196 208 stim::vec<float> next( //find next position.
197 209 pV[best][0],
198 210 pV[best][1],
... ... @@ -221,9 +233,13 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
221 233 #ifdef TIMING
222 234 size_time += gpuStopTimer();
223 235 #endif
224   - #ifdef TESTING
225   - #endif
226 236 int best = getCost(t_mag.getTexture(), t_mag.getAuxArray(), numSamplesMag); //get best cost.
  237 + #ifdef DEBUG
  238 + name.str("");
  239 + name << "Final_Cost_Size_" << iter << "_" << iter_siz << ".bmp";
  240 + test(t_mag.getTexture(), n_pixels*2.0, numSamplesMag*n_pixels, name.str());
  241 + iter_siz++;
  242 + #endif
227 243 setMagnitude(m*mV[best]); //adjust the magnitude.
228 244 }
229 245  
... ... @@ -242,7 +258,6 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
242 258 if(cL.size() < 4){} ///if the size of the fiber is less then 4 we do nothing.
243 259 else{
244 260 setMatrix(1); ///finds the current transformation matrix
245   -
246 261 DrawLongCylinder(n, l_template, l_square); ///Draw the cylinder.
247 262 stim::cylinder<float> cyl(cL, cM);
248 263 std::vector< stim::vec<float> > result = find_branch(cylinder_texID, GL_TEXTURE_2D, n*l_square, (cL.size()-1)*l_template); ///find all the centers in cuda
... ... @@ -269,8 +284,7 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
269 284 }
270 285 stim::vec3<float> v = cyl.surf(pval, result[i][0]); ///find the coordinates of the point at pval on the surface in tissue space.
271 286 stim::vec3<float> di = cyl.p(pval); ///find the coord of v in tissue space projected on the centerline.
272   - float rad = cyl.r(pval); ///find the radius at the pvalue's location
273   - std::cout << v << rad << std::endl;
  287 + float rad = cyl.r(pval)/2; ///find the radius at the pvalue's location
274 288 if(
275 289 !(v[0] > size[0] || v[1] > size[1]
276 290 || v[2] > size[2] || v[0] < 0
... ... @@ -355,15 +369,23 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
355 369 stim::vec3<float> sph(1, theta, phi); ///combine into a vector in spherical coordinates.
356 370 stim::vec3<float> cart = sph.sph2cart();///convert to cartesian.
357 371 dV.push_back(cart); ///save the generated vector for further use.
  372 + #ifdef DEBUG
  373 +// std::cout << cart << std::endl;
  374 + #endif
358 375 if(cos(Y.dot(cart)) < 0.087) ///make sure that the Y is not parallel to the new vector.
359 376 {
360 377 Y[0] = 0.0; Y[1] = 1.0;
361 378 }else{
362 379 Y[0] = 1.0; Y[1] = 0.0;
363 380 }
  381 +
364 382 hor = stim::rect<float>(mag, ///generate a rectangle with the new vectro as a normal.
365 383 pos, cart,
366 384 ((Y.cross(cart)).cross(cart)).norm());
  385 +
  386 + #ifdef DEBUG
  387 + // std::cout << hor.n() << std::endl;
  388 + #endif
367 389 ver = stim::rect<float>(mag, ///generate another rectangle that's perpendicular the first but parallel to the cart vector.
368 390 pos, cart,
369 391 hor.n());
... ... @@ -389,10 +411,19 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
389 411  
390 412 //Set up the variable necessary for vector creation.
391 413 glNewList(dList+1, GL_COMPILE); ///generate a new display list.
392   - for(int i = 0; i < numSamplesPos; i++) ///for the number of position samples
  414 + pV.push_back(pos);
  415 + hor = stim::rect<float>(mag, ///generate a rec tangle with the new vector as a normal.
  416 + pos, dir,
  417 + ((Y.cross(d)).cross(d))
  418 + .norm());
  419 + ver = stim::rect<float>(mag, ///generate anoth er rectangle that's perpendicular the first but parallel to the cart vector.
  420 + pos, dir,
  421 + hor.n());
  422 + UpdateBuffer(0.0, 0.0+0*n_pixels);
  423 + for(int i = 1; i < numSamplesPos; i++) ///for the number of position samples
393 424 {
394 425 stim::vec3<float> temp = uniformRandVector(); ///generate a random point on a plane.
395   - temp = temp*delta*2.0 - delta/2.0; ///scale the point
  426 + temp = temp*delta*2.0 - delta; ///scale the point and center it around 0.0
396 427 temp[2] = 0.0;
397 428 pV.push_back(temp); ///save the point for further use.
398 429 hor = stim::rect<float>(mag, ///generate a rectangle with the new vector as a normal.
... ... @@ -405,6 +436,10 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
405 436 UpdateBuffer(0.0, 0.0+i*n_pixels); ///sample the necessary points and put them into a display list.
406 437 }
407 438 glEndList(); ///finilize the display list.
  439 + #ifdef DEBUG
  440 + for(int i = 0; i < numSamplesPos; i++)
  441 + std::cout << pV[i] << std::endl;
  442 + #endif
408 443 }
409 444  
410 445 ///@param float delta, How much the rectangles are allowed to expand.
... ... @@ -428,7 +463,6 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
428 463 float max = 1.0+delta; ///largers size.
429 464 float step = (max-min)/(numSamplesMag-1); ///the size variation from one rect to the next.
430 465 float factor;
431   -
432 466 glNewList(dList+2, GL_COMPILE);
433 467 for(int i = 0; i < numSamplesMag; i++){ ///for the number of position samples
434 468 //Create linear index
... ... @@ -545,8 +579,8 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
545 579 glBindTexture(GL_TEXTURE_2D, textureID);
546 580  
547 581 //Textures repeat and use linear interpolation, luminance format.
548   - glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT); ///Set up the texture to repeat at edges.
549   - glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_REPEAT);
  582 + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_BORDER); ///Set up the texture to repeat at edges.
  583 + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_BORDER);
550 584 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR); ///Set up the texture to use Linear interpolation
551 585 glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
552 586 glTexImage2D(GL_TEXTURE_2D, 0, GL_LUMINANCE, ///Create the texture with no data.
... ... @@ -746,7 +780,12 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
746 780 numSamples = samples;
747 781 numSamplesPos = samplespos;
748 782 numSamplesMag = samplesmag;
  783 + #ifdef DEBUG
749 784 iter = 0;
  785 + iter_pos = 0;
  786 + iter_dir = 0;
  787 + iter_siz = 0;
  788 + #endif
750 789 }
751 790  
752 791 ///Position constructor: floats.
... ... @@ -770,7 +809,12 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
770 809 numSamples = numsamples;
771 810 numSamplesPos = numsamplespos;
772 811 numSamplesMag = numsamplesmag;
  812 + #ifdef DEBUG
773 813 iter = 0;
  814 + iter_pos = 0;
  815 + iter_dir = 0;
  816 + iter_siz = 0;
  817 + #endif
774 818 }
775 819  
776 820 ///Position constructor: vecs of floats.
... ... @@ -789,7 +833,12 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
789 833 numSamples = samples;
790 834 numSamplesPos = samplesPos;
791 835 numSamplesMag = samplesMag;
  836 + #ifdef DEBUG
792 837 iter = 0;
  838 + iter_pos = 0;
  839 + iter_dir = 0;
  840 + iter_siz = 0;
  841 + #endif
793 842 }
794 843  
795 844 ///destructor
... ... @@ -825,6 +874,9 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
825 874 #endif
826 875 #ifdef DEBUG
827 876 iter = 0;
  877 + iter_pos = 0;
  878 + iter_dir = 0;
  879 + iter_siz = 0;
828 880 #endif
829 881 stepsize = 3.0;
830 882 n_pixels = 16.0;
... ... @@ -852,7 +904,7 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
852 904 genDirectionVectors(5*stim::PI/4);
853 905 Unbind();
854 906 Bind(position_texID, position_buffID, numSamplesPos, n_pixels);
855   - genPositionVectors();
  907 + genPositionVectors(0.2);
856 908 Unbind();
857 909 Bind(radius_texID, radius_buffID, numSamplesMag, n_pixels);
858 910 genMagnitudeVectors();
... ... @@ -983,17 +1035,20 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
983 1035 getRotation(stim::vec3<float> dir)
984 1036 {
985 1037 stim::vec<float> out(0.0,0.0,0.0,0.0);
986   - stim::vec<float> from(0.0,0.0,1.0);
  1038 + stim::vec3<float> from(0.0,0.0,1.0);
987 1039 out[0] = acos(dir.dot(from))*180/stim::PI;
988   - if(out[0] < 1.0){
  1040 + #ifdef DEBUG
  1041 + std::cout << "out is " << out << std::endl;
  1042 + std::cout << "when rotating from " << from << " to " << dir << std::endl;
  1043 + #endif
  1044 + if(out[0] < 0.01){
989 1045 out[0] = 0.0;
990 1046 out[1] = 0.0;
991 1047 out[2] = 0.0;
992 1048 out[3] = 1.0;
993 1049 } else {
994   - stim::vec<float> temp(0.0, 0.0, 0.0);;
995   - stim::vec<float> dir1(dir[0], dir[1], dir[2]);
996   - temp = (from.cross(dir1)).norm();
  1050 + stim::vec3<float> temp(0.0, 0.0, 0.0);;
  1051 + temp = (from.cross(dir)).norm();
997 1052 out[1] = temp[0];
998 1053 out[2] = temp[1];
999 1054 out[3] = temp[2];
... ... @@ -1226,7 +1281,9 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
1226 1281 int
1227 1282 Step()
1228 1283 {
  1284 + #ifdef DEBUG
1229 1285 std::cerr << "Took a step" << std::endl;
  1286 + #endif
1230 1287 Bind(direction_texID, direction_buffID, numSamples, n_pixels);
1231 1288 CHECK_OPENGL_ERROR
1232 1289 findOptimalDirection();
... ... @@ -1376,7 +1433,19 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
1376 1433 setPosition(curSeed);
1377 1434 setDirection(curSeedVec);
1378 1435 setMagnitude(curSeedMag);
1379   - findOptimalScale();
  1436 +
  1437 + #ifdef DEBUG
  1438 + std::cout << "The new seed " << curSeed << curSeedVec << curSeedMag << std::endl;
  1439 + #endif
  1440 +
  1441 +// Bind(direction_texID, direction_buffID, numSamples, n_pixels);
  1442 +// CHECK_OPENGL_ERROR
  1443 +// findOptimalDirection();
  1444 +// Unbind();
  1445 +
  1446 +// Bind(radius_texID, radius_buffID, numSamplesMag, n_pixels);
  1447 +// findOptimalScale();
  1448 +// Unbind();
1380 1449  
1381 1450 // cL.push_back(curSeed);
1382 1451 // cM.push_back(curSeedMag);
... ... @@ -1552,6 +1621,9 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
1552 1621 sk.End();
1553 1622  
1554 1623 nt.addEdge(L,M);
  1624 + #ifdef DEBUG
  1625 + iter++;
  1626 + #endif
1555 1627 }
1556 1628 }
1557 1629  
... ... @@ -1646,6 +1718,9 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
1646 1718 }
1647 1719 }
1648 1720 }
  1721 + #ifdef DEBUG
  1722 + std::cout << "I broke out!" << std::endl;
  1723 + #endif
1649 1724 }
1650 1725 };
1651 1726 }
... ...