Commit 3e5d3ad3f2058215f2576aa0f0cf67dec1368064

Authored by Pavel Govyadinov
2 parents eb5dfb2b 27194b56

merged the changes from Master_Clone_W_Branch with into master

stim/cuda/branch_detection.cuh
@@ -5,167 +5,48 @@ @@ -5,167 +5,48 @@
5 //#include <math.h> 5 //#include <math.h>
6 #include <stim/visualization/colormap.h> 6 #include <stim/visualization/colormap.h>
7 #include <stim/cuda/cuda_texture.cuh> 7 #include <stim/cuda/cuda_texture.cuh>
8 -#include <stim/cuda/templates/gradient.cuh>  
9 -#include <stim/cuda/templates/gaussian_blur.cuh>  
10 -#include <stim/cuda/arraymath.cuh>  
11 -#include <stim/cuda/ivote.cuh>  
12 -#include <stim/cuda/testKernel.cuh> 8 +#include <stim/cuda/filter.cuh>
13 typedef unsigned int uint; 9 typedef unsigned int uint;
14 typedef unsigned int uchar; 10 typedef unsigned int uchar;
15 11
16 -stim::cuda::cuda_texture t;  
17 -float* gpuTable;  
18 -float* gpuGrad;  
19 -float* gpuVote;  
20 -float* gpuI;  
21 -float* gpuCenters;  
22 -  
23 -void atan_2d(float* cpuTable, unsigned int rmax)  
24 -{  
25 - //initialize the width and height of the window which atan2 are computed in.  
26 - int xsize = 2*rmax +1;  
27 - int ysize = 2*rmax +1;  
28 -  
29 - // assign the center coordinates of the atan2 window to yi and xi  
30 - int yi = rmax;  
31 - int xi = rmax;  
32 -  
33 -  
34 - for (int xt = 0; xt < xsize; xt++){  
35 -  
36 - for(int yt = 0; yt < ysize; yt++){  
37 -  
38 - //convert the current 2D coordinates to 1D  
39 - int id = yt * xsize + xt;  
40 - // calculate the distance between the pixel and the center of the atan2 window  
41 - float xd = xi - xt;  
42 - float yd = yi - yt;  
43 -  
44 - // calculate the angle between the pixel and the center of the atan2 window and store the result.  
45 - float atan_2d_vote = atan2(yd, xd);  
46 - cpuTable[id] = atan_2d_vote;  
47 - }  
48 - }  
49 -  
50 -}  
51 -  
52 -void initCuda(unsigned int bytes_table, unsigned int bytes_ds)  
53 -{  
54 - HANDLE_ERROR(  
55 - cudaMalloc((void**) &gpuTable, bytes_table)  
56 - );  
57 - HANDLE_ERROR(  
58 - cudaMalloc((void**) &gpuI, bytes_ds)  
59 - );  
60 - HANDLE_ERROR(  
61 - cudaMalloc((void**) &gpuGrad, bytes_ds*2)  
62 - );  
63 - HANDLE_ERROR(  
64 - cudaMalloc((void**) &gpuVote, bytes_ds)  
65 - );  
66 - HANDLE_ERROR(  
67 - cudaMalloc((void**) &gpuCenters, bytes_ds)  
68 - );  
69 -}  
70 -  
71 -void cleanCuda()  
72 -{  
73 - HANDLE_ERROR(  
74 - cudaFree(gpuTable)  
75 - );  
76 - HANDLE_ERROR(  
77 - cudaFree(gpuGrad)  
78 - );  
79 - HANDLE_ERROR(  
80 - cudaFree(gpuVote)  
81 - );  
82 - HANDLE_ERROR(  
83 - cudaFree(gpuCenters)  
84 - );  
85 - HANDLE_ERROR(  
86 - cudaFree(gpuI)  
87 - );  
88 -}  
89 12
90 std::vector< stim::vec<float> > 13 std::vector< stim::vec<float> >
91 find_branch(GLint texbufferID, GLenum texType, unsigned int x, unsigned int y) 14 find_branch(GLint texbufferID, GLenum texType, unsigned int x, unsigned int y)
92 { 15 {
93 - float phi = 15.1*M_PI/180;  
94 - int iter = 5;  
95 - float dphi = phi/iter;  
96 - float rmax = 10;  
97 - float sigma = 4;  
98 - unsigned int pixels = x * y;  
99 - unsigned int bytes = sizeof(float) * pixels;  
100 - unsigned int bytes_table = sizeof(float) * (2*rmax + 1) * (2*rmax + 1);  
101 - unsigned int x_ds = (x + (x % 1 == 0 ? 0:1));  
102 - unsigned int y_ds = (y + (x % 1 == 0 ? 0:1));  
103 - unsigned int bytes_ds = sizeof(float) * x_ds * y_ds;  
104 - unsigned int conn = 5;  
105 - float final_t = 200.0;  
106 - float* cpuTable = (float*) malloc(bytes_table);  
107 - float* cpuCenters = (float*) malloc(bytes_ds); 16 + float sigma = 2.0;
  17 + unsigned int conn = 7;
  18 + float threshold = 40.0;
  19 + float* cpuCenters = (float*) malloc(x*y*sizeof(float));
  20 + int sizek = 7;
108 21
109 stringstream name; 22 stringstream name;
110 23
111 24
112 -  
113 -  
114 - std::vector<stim::vec<float> > output;  
115 - initCuda(bytes_table, bytes_ds);  
116 -  
117 - atan_2d(cpuTable, rmax);  
118 - cudaMemcpy(gpuTable, cpuTable, bytes_table, cudaMemcpyHostToDevice);  
119 -  
120 -  
121 - t.MapCudaTexture(texbufferID, texType);  
122 - cudaDeviceSynchronize();  
123 - stim::cuda::tex_gaussian_blur2<float>(  
124 - gpuI, sigma, x, y, t.getTexture(), t.getArray()  
125 - ); 25 + cpuCenters = stim::cuda::get_centers(texbufferID, texType, x, y, sizek, sigma, conn, threshold);
126 cudaDeviceSynchronize(); 26 cudaDeviceSynchronize();
127 27
128 28
129 - stim::cuda::gpu_gradient_2d<float>(  
130 - gpuGrad, gpuI, x, y  
131 - );  
132 - cudaDeviceSynchronize();  
133 -  
134 - stim::cuda::gpu_cart2polar<float>(gpuGrad, x, y);  
135 - cudaDeviceSynchronize();  
136 29
137 - cudaDeviceSynchronize();  
138 - for (int i = 0; i < iter; i++)  
139 - {  
140 - stim::cuda::gpu_vote<float>(gpuVote, gpuGrad, gpuTable, phi, rmax, x, y);  
141 - cudaDeviceSynchronize();  
142 - stim::cuda::gpu_update_dir<float>(gpuVote, gpuGrad, gpuTable, phi, rmax, x, y);  
143 - cudaDeviceSynchronize();  
144 - phi = phi - dphi;  
145 - } 30 + std::vector<stim::vec<float> > output;
146 31
147 cudaDeviceSynchronize(); 32 cudaDeviceSynchronize();
148 - stim::cuda::gpu_local_max<float>(gpuCenters, gpuVote, final_t, conn, x, y);  
149 - cudaMemcpy(cpuCenters, gpuCenters, bytes_ds, cudaMemcpyDeviceToHost);  
150 - for(int i = 0; i < pixels; i++) 33 +
  34 + for(int i = 0; i < x; i++)
151 { 35 {
152 - int ix = (i % x);  
153 - int iy = (i / x);  
154 - if((cpuCenters[i] == 1) && (ix > 4) && (ix < x-4)) 36 + for(int j = 0; j < y; j++)
155 { 37 {
156 -  
157 - float x_v = (float) ix;  
158 - float y_v = (float) iy;  
159 - output.push_back(stim::vec<float>((x_v/(float)x),  
160 - (y_v/(float)y), 0.0));  
161 -  
162 - } 38 + int idx = x*j+i;
  39 + if(cpuCenters[idx] != 0)
  40 + {
  41 + float x_v = (float) i;
  42 + float y_v = (float) j;
  43 + output.push_back(stim::vec<float>((x_v/(float)x*360.0),
  44 + (y_v), y_v/8));
  45 + }
  46 +
  47 + }
163 } 48 }
164 -  
165 -  
166 - t.UnmapCudaTexture();  
167 - cleanCuda();  
168 - free(cpuTable); 49 +
169 free(cpuCenters); 50 free(cpuCenters);
170 return output; 51 return output;
171 } 52 }
stim/cuda/branch_detection2.cuh deleted
1 -#include <stim/cuda/templates/gaussian_blur.cuh>  
2 -#include <stim/cuda/templates/gradient.cuh>  
3 -#include <stim/cuda/arraymath.cuh>  
4 -#include <stim/cuda/ivote.cuh>  
5 -  
6 -  
7 -  
8 -  
9 -  
10 -  
11 -  
12 -  
13 -  
14 -  
15 -void atan_2(float* cpuTable, unsigned int rmax){  
16 -  
17 - //initialize the width and height of the window which atan2 are computed in.  
18 - int xsize = 2*rmax +1;  
19 - int ysize = 2*rmax +1;  
20 -  
21 - // assign the center coordinates of the atan2 window to yi and xi  
22 - int yi = rmax;  
23 - int xi = rmax;  
24 -  
25 -  
26 - for (int xt = 0; xt < xsize; xt++){  
27 -  
28 - for(int yt = 0; yt < ysize; yt++){  
29 -  
30 - //convert the current 2D coordinates to 1D  
31 - int id = yt * xsize + xt;  
32 - // calculate the distance between the pixel and the center of the atan2 window  
33 - float xd = xi - xt;  
34 - float yd = yi - yt;  
35 -  
36 - // calculate the angle between the pixel and the center of the atan2 window and store the result.  
37 - float atan_2d_vote = atan2(yd, xd);  
38 - cpuTable[id] = atan_2d_vote;  
39 - }  
40 - }  
41 -  
42 -}  
43 -std::vector<stim::vec<float> >  
44 -find_branch(GLint texbufferID, GLenum texType, unsigned int x, unsigned int y)  
45 -{  
46 -  
47 - float* cpuTable = (float  
48 -  
49 - unsigned int pixels = x * y;  
50 - unsigned int bytes = sizeof(float) * pixels;  
51 -  
52 - //calculate the number of bytes in the atan2 table  
53 -  
54 - unsigned int bytes_table = (2*rmax+1) * (2*rmax+1) * sizeof(float);  
55 -  
56 -  
57 -  
58 - //allocate space on the GPU for the atan2 table  
59 -  
60 - float* gpuTable;  
61 -  
62 - cudaMalloc(&gpuTable, bytes_table);  
63 -  
64 -  
65 -  
66 - cudaMemcpy(gpuTable, cpuTable, bytes_table, cudaMemcpyHostToDevice);  
67 -  
68 - unsigned int sigma_ds = 1/resize;  
69 - unsigned int x_ds = (x/sigma_ds + (x %sigma_ds == 0 ? 0:1));  
70 - unsigned int y_ds = (y/sigma_ds + (y %sigma_ds == 0 ? 0:1));  
71 - unsigned int bytes_ds = sizeof(float) * x_ds * y_ds;  
72 -  
73 -  
74 - float* gpuI;  
75 - cudaMalloc(&gpuI, bytes_ds);  
76 -  
77 -  
78 - float* gpuGrad;  
79 - cudaMalloc(&gpuGrad, bytes_ds*2);  
80 -  
81 - float* gpuVote;  
82 - cudaMalloc(&gpuVote, bytes_ds);  
83 -  
84 - // allocate space on the GPU for the detected cell centes  
85 -  
86 - float* gpuCenters;  
87 -  
88 - cudaMalloc(&gpuCenters, bytes_ds);  
89 -  
90 -  
91 - stim::cuda::gpu_down_sample<float>(gpuI, gpuI0, resize, x , y);  
92 - cudaMemcpy(cpuResize, gpuI, bytes_ds, cudaMemcpyDeviceToHost);  
93 -  
94 -x = x_ds;  
95 - y = y_ds;  
96 - t = t * resize;  
97 - //sigma = sigma * resize;  
98 -  
99 - cudaDeviceSynchronize();  
100 - stim::cuda::gpu_gaussian_blur2<float>(gpuI,sigma, x, y);  
101 - cudaDeviceSynchronize();  
102 - cudaMemcpy(cpuBlur, gpuI, bytes_ds, cudaMemcpyDeviceToHost);  
103 - cudaDeviceSynchronize();  
104 -  
105 - stim::cuda::gpu_gradient_2d<float>(gpuGrad, gpuI, x, y);  
106 - cudaDeviceSynchronize();  
107 - cudaMemcpy(cpuGradient, gpuGrad, bytes_ds*2, cudaMemcpyDeviceToHost);  
108 -  
109 - stim::cuda::gpu_cart2polar<float>(gpuGrad, x, y);  
110 - cudaDeviceSynchronize();  
111 - cudaMemcpy(cpuCart2Polar, gpuGrad, bytes_ds*2, cudaMemcpyDeviceToHost);  
112 -  
113 -  
114 - //multiply the gradient by a constant and calculate the absolute value (to save an image)  
115 -  
116 - stim::cuda::cpu_multiply<float>(cpuCart2Polar, 40, x * y * 2);  
117 -  
118 - cudaDeviceSynchronize();  
119 -  
120 - stim::cuda::cpu_abs<float>(cpuCart2Polar, x * y * 2);  
121 -  
122 - cudaDeviceSynchronize();  
123 -  
124 -  
125 - for (int i =0; i<iter; i++){  
126 -  
127 - stim::cuda::gpu_vote<float>(gpuVote, gpuGrad, gpuTable, phi, rmax, x, y);  
128 - cudaDeviceSynchronize();  
129 - stim::cuda::gpu_update_dir<float>(gpuVote, gpuGrad, gpuTable, phi, rmax, x, y);  
130 - cudaDeviceSynchronize();  
131 - switch (i){  
132 - case 0 : cudaMemcpy(cpuVote1, gpuVote, bytes_ds, cudaMemcpyDeviceToHost);  
133 - break;  
134 - case 1 : cudaMemcpy(cpuVote2, gpuVote, bytes_ds, cudaMemcpyDeviceToHost);  
135 - break;  
136 - case 2 : cudaMemcpy(cpuVote3, gpuVote, bytes_ds, cudaMemcpyDeviceToHost);  
137 - break;  
138 - case 3 : cudaMemcpy(cpuVote4, gpuVote, bytes_ds, cudaMemcpyDeviceToHost);  
139 - break;  
140 - case 4 : cudaMemcpy(cpuVote5, gpuVote, bytes_ds, cudaMemcpyDeviceToHost);  
141 - break;  
142 - default : cudaMemcpy(cpuVote5, gpuVote, bytes_ds, cudaMemcpyDeviceToHost);  
143 - break;  
144 - }  
145 - phi = phi - dphi;  
146 - }  
147 -  
148 - stim::cuda::gpu_local_max<float>(gpuCenters, gpuVote, t, conn, x, y);  
149 - cudaMemcpy(cpuCenters, gpuCenters, bytes_ds, cudaMemcpyDeviceToHost);  
150 -  
151 -}  
stim/cuda/cuda_texture.cuh
@@ -41,6 +41,46 @@ namespace stim @@ -41,6 +41,46 @@ namespace stim
41 texDesc.normalizedCoords = 0; 41 texDesc.normalizedCoords = 0;
42 } 42 }
43 43
  44 +
  45 + ///Enable the nromalized texture coordinates.
  46 + ///@param bool, 1 for on, 0 for off
  47 + void
  48 + SetTextureCoordinates(bool val)
  49 + {
  50 + if(val)
  51 + texDesc.normalizedCoords = 1;
  52 + else
  53 + texDesc.normalizedCoords = 0;
  54 + }
  55 +
  56 + ///sets the dimension dim to used the mode at the borders of the texture.
  57 + ///@param dim : 0-x, 1-y, 2-z
  58 + ///@param mode: cudaAddressModeWrap = 0,
  59 + /// cudaAddressModeClamp = 1,
  60 + /// cudaAddressMNodeMirror = 2,
  61 + /// cudaAddressModeBorder = 3,
  62 + void
  63 + SetAddressMode(int dim, int mode)
  64 + {
  65 + switch(mode)
  66 + {
  67 + case 0:
  68 + texDesc.addressMode[dim] = cudaAddressModeWrap;
  69 + break;
  70 + case 1:
  71 + texDesc.addressMode[dim] = cudaAddressModeClamp;
  72 + break;
  73 + case 2:
  74 + texDesc.addressMode[dim] = cudaAddressModeMirror;
  75 + break;
  76 + case 3:
  77 + texDesc.addressMode[dim] = cudaAddressModeBorder;
  78 + break;
  79 + default:
  80 + break;
  81 + }
  82 + }
  83 +
44 //-------------------------------------------------------------------------// 84 //-------------------------------------------------------------------------//
45 //-------------------------------CUDA_MAPPING------------------------------// 85 //-------------------------------CUDA_MAPPING------------------------------//
46 //-------------------------------------------------------------------------// 86 //-------------------------------------------------------------------------//
stim/cuda/filter.cuh 0 → 100644
  1 +#ifndef STIM_FILTER_H
  2 +#define STIM_FILTER_H
  3 +
  4 +#include <assert.h>
  5 +#include <cuda.h>
  6 +#include <cuda_runtime.h>
  7 +#include <stdio.h>
  8 +#include <stim/visualization/colormap.h>
  9 +#include <sstream>
  10 +#include <stim/cuda/cudatools/devices.h>
  11 +#include <stim/cuda/cudatools/threads.h>
  12 +#include <stim/cuda/cuda_texture.cuh>
  13 +#include <stim/cuda/ivote.cuh>
  14 +#include <stim/cuda/arraymath.cuh>
  15 +
  16 +#define IMAD(a,b,c) ( __mul24((a), (b)) + (c) )
  17 +#define M_PI 3.141592654f
  18 +
  19 +
  20 +namespace stim
  21 +{
  22 + namespace cuda
  23 + {
  24 +
  25 + float* gpuLoG;
  26 + float* LoG;
  27 + float* res;
  28 + float* centers;
  29 + stim::cuda::cuda_texture tx;
  30 +
  31 +
  32 +
  33 + void initArray(int DIM_X, int DIM_Y, int kl)
  34 + {
  35 +
  36 + LoG = (float*) malloc(kl*kl*sizeof(float));
  37 + HANDLE_ERROR(
  38 + cudaMalloc( (void**) &gpuLoG, kl*kl*sizeof(float))
  39 + );
  40 + // checkCUDAerrors("Memory Allocation, LoG");
  41 + HANDLE_ERROR(
  42 + cudaMalloc( (void**) &res, DIM_Y*DIM_X*sizeof(float))
  43 + );
  44 + HANDLE_ERROR(
  45 + cudaMalloc( (void**) &centers, DIM_Y*DIM_X*sizeof(float))
  46 + );
  47 + // checkCUDAerrors("Memory Allocation, Result");
  48 + }
  49 +
  50 + void cleanUp(cudaGraphicsResource_t src)
  51 + {
  52 + HANDLE_ERROR(
  53 + cudaFree(gpuLoG)
  54 + );
  55 + HANDLE_ERROR(
  56 + cudaFree(res)
  57 + );
  58 + HANDLE_ERROR(
  59 + cudaFree(centers)
  60 + );
  61 + free(LoG);
  62 + }
  63 +
  64 + void
  65 + filterKernel (float kl, float sigma, float *LoG)
  66 + {
  67 + float t = 0.0;
  68 + float kr = kl/2;
  69 + float x, y;
  70 + int idx;
  71 + for(int i = 0; i < kl; i++){
  72 + for(int j = 0; j < kl; j++){
  73 + idx = j*kl+i;
  74 + x = i - 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))
  77 + *expf(-(powf(x,2)+powf(y,2))/2/powf(sigma,2));
  78 + t +=LoG[idx];
  79 + }
  80 + }
  81 +
  82 + for(int i = 0; i < kl*kl; i++)
  83 + {
  84 + LoG[i] = LoG[i]/t;
  85 + }
  86 +
  87 + }
  88 +
  89 + //Shared memory would be better.
  90 + __global__
  91 + void
  92 + applyFilter(cudaTextureObject_t texIn, unsigned int DIM_X, unsigned int DIM_Y, int kr, int kl, float *res, float* gpuLoG){
  93 + //R = floor(size/2)
  94 + //THIS IS A NAIVE WAY TO DO IT, and there is a better way)
  95 +
  96 + __shared__ float shared[7][7];
  97 + int x = blockIdx.x;
  98 + int y = blockIdx.y;
  99 + int xi = threadIdx.x;
  100 + int yi = threadIdx.y;
  101 + float val = 0;
  102 + float tu = (x-kr+xi)/(float)DIM_X;
  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));
  105 + __syncthreads();
  106 +
  107 +
  108 + //x = max(0,x);
  109 + //x = min(x, width-1);
  110 + //y = max(y, 0);
  111 + //y = min(y, height - 1);
  112 +
  113 + int idx = y*DIM_X+x;
  114 + int k_idx;
  115 + for(unsigned int step = blockDim.x/2; step >= 1; step >>= 1)
  116 + {
  117 + __syncthreads();
  118 + if (xi < step)
  119 + {
  120 + shared[xi][yi] += shared[xi + step][yi];
  121 + }
  122 + __syncthreads();
  123 + }
  124 + __syncthreads();
  125 +
  126 + for(unsigned int step = blockDim.y/2; step >= 1; step >>= 1)
  127 + {
  128 + __syncthreads();
  129 + if(yi < step)
  130 + {
  131 + shared[xi][yi] += shared[xi][yi + step];
  132 + }
  133 + __syncthreads();
  134 + }
  135 + __syncthreads();
  136 + if(xi == 0 && yi == 0)
  137 + res[idx] = shared[0][0];
  138 + }
  139 +
  140 + extern "C"
  141 + float *
  142 + get_centers(GLint texbufferID, GLenum texType, int DIM_X, int DIM_Y, int sizeK, float sigma, float conn, float threshold)
  143 + {
  144 + tx.SetTextureCoordinates(1);
  145 + tx.SetAddressMode(1, 3);
  146 + tx.MapCudaTexture(texbufferID, texType);
  147 + float* result = (float*) malloc(DIM_X*DIM_Y*sizeof(float));
  148 +
  149 + initArray(DIM_X, DIM_Y, sizeK);
  150 +
  151 + filterKernel(sizeK, sigma, LoG);
  152 + cudaMemcpy(gpuLoG, LoG, sizeK*sizeK*sizeof(float), cudaMemcpyHostToDevice);
  153 + dim3 numBlocks(DIM_X, DIM_Y);
  154 + dim3 threadsPerBlock(sizeK, sizeK);
  155 +
  156 + applyFilter <<< numBlocks, threadsPerBlock >>> (tx.getTexture(), DIM_X, DIM_Y, floor(sizeK/2), sizeK, res, gpuLoG);
  157 +
  158 +
  159 + stim::cuda::gpu_local_max<float>(centers, res, threshold, conn, DIM_X, DIM_Y);
  160 +
  161 + cudaDeviceSynchronize();
  162 +
  163 +
  164 + cudaMemcpy(result, centers, DIM_X*DIM_Y*sizeof(float), cudaMemcpyDeviceToHost);
  165 +
  166 + tx.UnmapCudaTexture();
  167 + cleanUP();
  168 + return result;
  169 + }
  170 +
  171 + }
  172 +}
  173 +#endif
stim/cuda/filter.h deleted
1 -#include <assert.h>  
2 -#include <cuda.h>  
3 -#include <cuda_runtime.h>  
4 -#include <stdio.h>  
5 -#include <stim/visualization/colormap.h>  
6 -#include <sstream>  
7 -  
8 -#define IMAD(a,b,c) ( __mul24((a), (b)) + (c) )  
9 -  
10 -int kr;  
11 -int kl;  
12 -float sigma;  
13 -float* LoG;  
14 -float* result;  
15 -cudaArray* srcArray;  
16 -texture<uchar, cudaTextureType2D, cudaReadModeElementType> texIn;  
17 -  
18 -  
19 -__device__ float filterKernel ()  
20 -{  
21 - float t = 0;  
22 - idx = j*kl+i;  
23 - for(int i = 0; i < kl; i++){  
24 - for(int j = 0; j < kl; j++){  
25 - x = i - floor(kl);  
26 - y = j - floor(kl);  
27 - LoG(idx) = (-1/M_PI/sigma^4)* (1 - (x^2+y^2)/2/sigma^2)  
28 - *exp(-(x^2+y^2)/2/sigma^2);  
29 - t +=LoG(idx);  
30 - }  
31 - }  
32 - LoG =/ t;  
33 -}  
34 -  
35 -void initArray(cudaGraphicsResource_t src, int DIM_X, int DIM_Y)  
36 -{  
37 - HANDLE_ERROR(  
38 - cudaGraphicsMapResources(1, &src)  
39 - );  
40 - HANDLE_ERROR(  
41 - cudaGraphicsSubResourceGetMappedArray(&srcArray, src, 0,0)  
42 - );  
43 - HANDLE_ERROR(  
44 - cudaBindTertureToArray(texIn, srcArray)  
45 - );  
46 - cudaMalloc( (void**) &LoG, kl*kl*sizeof(float));  
47 - checkCUDAerrors("Memory Allocation, LoG");  
48 - cudaMalloc( (void**) &result, DIM_Y*DIM_X*sizeof(float));  
49 - checkCUDAerrors("Memory Allocation, Result");  
50 -}  
51 -  
52 -void cleanUp(cudaGraphicsResource_t src);  
53 -{  
54 - HANDLE_ERROR(  
55 - cudaUnbindTexture(texIn)  
56 - );  
57 - HANDLE_ERROR(  
58 - cudaFree(LoG)  
59 - );  
60 - HANDLE_ERROR(  
61 - cudaFree(result)  
62 - );  
63 - HANDLE_ERROR(  
64 - cudaGraphicsUnmapResources(1, &src)  
65 - );  
66 -}  
67 -  
68 -//Shared memory would be better.  
69 -__global__  
70 -void  
71 -applyFilter(unsigned int DIM_X, unsigned int DIM_Y){  
72 -//R = floor(size/2)  
73 -//THIS IS A NAIVE WAY TO DO IT, and there is a better way)  
74 - //__shared__ float shared[(DIM_X+2*R), (DIM_Y+2*R)];  
75 -  
76 - const int x = IMAD(blockDim.x, blockIdx.x, threadIdx.x);  
77 - const int y = IMAD(blockDim.y, blockIdx.y, threadIdx.y);  
78 - float val = 0;  
79 - //x = max(0,x);  
80 - //x = min(x, width-1);  
81 - //y = max(y, 0);  
82 - //y = min(y, height - 1);  
83 -  
84 - int idx = y*DIM_X+x;  
85 - //unsigned int bindex = threadIdx.y * blockDim.y + threadIdx.x;  
86 -  
87 - //float valIn = tex2D(texIn, x, y);  
88 - for (int i = -kr; i <= kr; i++){ //rows  
89 - for (int j = -kr; i <= kr; j++){ //colls  
90 - k_idx = (j+kr)+(i+kr)*kl;  
91 - xi = max(0, x+i);  
92 - xi = min(x+i, DIM_X-1);  
93 - yj = max(y+j, 0);  
94 - yj = min(y+j, DIM_Y-1);  
95 - val += LoG(k_idx)*tex2D(texIn,x+i, y+j);  
96 - }  
97 - }  
98 -  
99 - result[idx] = val;  
100 -}  
stim/cuda/sharedmem.cuh
@@ -35,34 +35,6 @@ namespace stim{ @@ -35,34 +35,6 @@ namespace stim{
35 } 35 }
36 } 36 }
37 37
38 - template<typename T, typename D>  
39 - __device__ void sharedMemcpy_tex2D(T* dest, cudaTextureObject_t src,  
40 - unsigned int x, unsigned int y, unsigned int X, unsigned int Y,  
41 - dim3 threadIdx, dim3 blockDim){  
42 -  
43 - //calculate the number of iterations required for the copy  
44 - unsigned int xI, yI;  
45 - xI = X/blockDim.x + 1; //number of iterations along X  
46 - yI = Y/blockDim.y + 1; //number of iterations along Y  
47 -  
48 - //for each iteration  
49 - for(unsigned int xi = 0; xi < xI; xi++){  
50 - for(unsigned int yi = 0; yi < yI; yi++){  
51 -  
52 - //calculate the index into shared memory  
53 - unsigned int sx = xi * blockDim.x + threadIdx.x;  
54 - unsigned int sy = yi * blockDim.y + threadIdx.y;  
55 -  
56 - //calculate the index into the texture  
57 - unsigned int tx = x + sx;  
58 - unsigned int ty = y + sy;  
59 -  
60 - //perform the copy  
61 - if(sx < X && sy < Y)  
62 - dest[sy * X + sx] = abs(255 - tex2D<D>(src, tx, ty));  
63 - }  
64 - }  
65 - }  
66 38
67 } 39 }
68 } 40 }
stim/cuda/testKernel.cuh
@@ -14,9 +14,9 @@ @@ -14,9 +14,9 @@
14 ///Initialization function, allocates the memory and passes the necessary 14 ///Initialization function, allocates the memory and passes the necessary
15 ///handles from OpenGL and Cuda. 15 ///handles from OpenGL and Cuda.
16 ///@param DIM_Y --integer controlling how much memory to allocate. 16 ///@param DIM_Y --integer controlling how much memory to allocate.
17 - void initArray() 17 + void initArray(int x, int y)
18 { 18 {
19 - cudaMalloc( (void**) &print, 216*16*sizeof(float)); ///temporary 19 + cudaMalloc( (void**) &print, x*y*sizeof(float)); ///temporary
20 } 20 }
21 21
22 ///Deinit function that frees the memery used and releases the texture resource 22 ///Deinit function that frees the memery used and releases the texture resource
@@ -48,12 +48,14 @@ @@ -48,12 +48,14 @@
48 { 48 {
49 int x = threadIdx.x + blockIdx.x * blockDim.x; 49 int x = threadIdx.x + blockIdx.x * blockDim.x;
50 int y = threadIdx.y + blockIdx.y * blockDim.y; 50 int y = threadIdx.y + blockIdx.y * blockDim.y;
51 - int idx = y*16+x; 51 + int idx = y*64+x;
  52 +// int idx = y*32+x;
  53 +// int idx = y*16+x;
52 54
53 float valIn = tex2D<unsigned char>(texIn, x, y); 55 float valIn = tex2D<unsigned char>(texIn, x, y);
54 float templa = templ(x); 56 float templa = templ(x);
55 - //print[idx] = abs(valIn); ///temporary  
56 - print[idx] = abs(templa); ///temporary 57 + print[idx] = valIn; ///temporary
  58 + //print[idx] = abs(templa); ///temporary
57 59
58 } 60 }
59 61
@@ -64,7 +66,7 @@ @@ -64,7 +66,7 @@
64 ///@param GLenum texType --either GL_TEXTURE_1D, GL_TEXTURE_2D or GL_TEXTURE_3D 66 ///@param GLenum texType --either GL_TEXTURE_1D, GL_TEXTURE_2D or GL_TEXTURE_3D
65 /// may work with other gl texture types, but untested. 67 /// may work with other gl texture types, but untested.
66 ///@param DIM_Y, the number of samples in the template. 68 ///@param DIM_Y, the number of samples in the template.
67 - void test(GLint texbufferID, GLenum texType) 69 + void test(GLint texbufferID, GLenum texType, int x, int y)
68 { 70 {
69 71
70 //Bind the Texture in GL and allow access to cuda. 72 //Bind the Texture in GL and allow access to cuda.
@@ -72,16 +74,12 @@ @@ -72,16 +74,12 @@
72 74
73 //initialize the return arrays. 75 //initialize the return arrays.
74 76
75 - initArray();  
76 -  
77 - int x = 16;  
78 - int y = 27*8;  
79 - y = 8* 1089; 77 + initArray(x,y);
  78 + dim3 numBlocks(1, y);
  79 + dim3 threadsPerBlock(x, 1);
80 int max_threads = stim::maxThreadsPerBlock(); 80 int max_threads = stim::maxThreadsPerBlock();
81 //dim3 threads(max_threads, 1); 81 //dim3 threads(max_threads, 1);
82 //dim3 blocks(x / threads.x + 1, y); 82 //dim3 blocks(x / threads.x + 1, y);
83 - dim3 numBlocks(1, 1089);  
84 - dim3 threadsPerBlock(16, 8);  
85 //dim3 numBlocks(2, 2); 83 //dim3 numBlocks(2, 2);
86 //dim3 threadsPerBlock(8, 108); 84 //dim3 threadsPerBlock(8, 108);
87 85
@@ -92,7 +90,7 @@ @@ -92,7 +90,7 @@
92 cudaDeviceSynchronize(); 90 cudaDeviceSynchronize();
93 stringstream name; //for debugging 91 stringstream name; //for debugging
94 name << "FromTex.bmp"; 92 name << "FromTex.bmp";
95 - stim::gpu2image<float>(print, name.str(),16,1089*8,0,1.0); 93 + stim::gpu2image<float>(print, name.str(),x,y,0,255);
96 94
97 tx.UnmapCudaTexture(); 95 tx.UnmapCudaTexture();
98 cleanUP(); 96 cleanUP();
stim/gl/gl_spider.h
@@ -23,6 +23,7 @@ @@ -23,6 +23,7 @@
23 #include <stim/cuda/branch_detection.cuh> 23 #include <stim/cuda/branch_detection.cuh>
24 #include "../../../volume-spider/fiber.h" 24 #include "../../../volume-spider/fiber.h"
25 #include "../../../volume-spider/glnetwork.h" 25 #include "../../../volume-spider/glnetwork.h"
  26 +#include <stim/visualization/cylinder.h>
26 //#include <stim/cuda/testKernel.cuh> 27 //#include <stim/cuda/testKernel.cuh>
27 28
28 //#include <stim/cuda/testKernel.cuh> 29 //#include <stim/cuda/testKernel.cuh>
@@ -81,8 +82,8 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -81,8 +82,8 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
81 int numSamplesPos; 82 int numSamplesPos;
82 int numSamplesMag; 83 int numSamplesMag;
83 84
84 -// float stepsize = 4.0; //Step size.  
85 - float stepsize = 3.0; //Step size. 85 + float stepsize = 5.0; //Step size.
  86 +// float stepsize = 3.0; //Step size.
86 int current_cost; //variable to store the cost of the current step. 87 int current_cost; //variable to store the cost of the current step.
87 88
88 89
@@ -95,7 +96,6 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -95,7 +96,6 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
95 std::vector< stim::vec<float> > cD; //Direction of line currently being traced. 96 std::vector< stim::vec<float> > cD; //Direction of line currently being traced.
96 std::vector< stim::vec<float> > cM; //Magnitude of line currently being traced. 97 std::vector< stim::vec<float> > cM; //Magnitude of line currently being traced.
97 98
98 - stim::glObj<float> sk; //object to store the skeleton.  
99 stim::glnetwork<float> nt; //object for storing the network. 99 stim::glnetwork<float> nt; //object for storing the network.
100 100
101 stim::vec<float> rev; //reverse vector; 101 stim::vec<float> rev; //reverse vector;
@@ -169,6 +169,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -169,6 +169,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
169 169
170 ///subject to change. 170 ///subject to change.
171 ///finds branches. 171 ///finds branches.
  172 + ///depreciated
172 void 173 void
173 branchDetection() 174 branchDetection()
174 { 175 {
@@ -196,14 +197,8 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -196,14 +197,8 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
196 -p[1] + cylp[1]*S[1]*R[1], 197 -p[1] + cylp[1]*S[1]*R[1],
197 -p[2] + cylp[2]*S[2]*R[2]); 198 -p[2] + cylp[2]*S[2]*R[2]);
198 seeddir = seeddir.norm(); 199 seeddir = seeddir.norm();
199 -// float seedm = m[0]/2.0;  
200 float seedm = m[0]; 200 float seedm = m[0];
201 // Uncomment for global run 201 // Uncomment for global run
202 -/* stim::vec<float> lSeed = getLastSeed();  
203 - if(sqrt(pow((lSeed[0] - vec[0]),2)  
204 - + pow((lSeed[1] - vec[1]),2) +  
205 - pow((lSeed[2] - vec[2]),2)) > m[0]/4.0  
206 - && */  
207 if( 202 if(
208 !(vec[0] > size[0] || vec[1] > size[1] 203 !(vec[0] > size[0] || vec[1] > size[1]
209 || vec[2] > size[2] || vec[0] < 0 204 || vec[2] > size[2] || vec[0] < 0
@@ -218,6 +213,56 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -218,6 +213,56 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
218 213
219 } 214 }
220 215
  216 +
  217 + ///finds all the branches in the a given fiber.
  218 + ///using LoG method.
  219 + void
  220 + branchDetection2(int n = 8, int l_template = 8, int l_square = 8)
  221 + {
  222 + if(cL.size() < 4){}
  223 + else{
  224 + setMatrix(1);
  225 + DrawLongCylinder(n, l_template, l_square);
  226 + stim::cylinder<float> cyl(cL, cM);
  227 + std::vector< stim::vec<float> > result = find_branch(btexbufferID, GL_TEXTURE_2D, n*l_square, (cL.size()-1)*l_template);
  228 + stim::vec<float> size(S[0]*R[0], S[1]*R[1], S[2]*R[2]);
  229 + float pval;
  230 + if(!result.empty())
  231 + {
  232 + for(int i = 0; i < result.size(); i++)
  233 + {
  234 + int id = result[i][2];
  235 + if(fmod(result[i][2], id) != 0 && id != 0)
  236 + {
  237 +
  238 + pval = ((cyl.getl(id+1)-cyl.getl(id))*
  239 + (fmod(result[i][2], id))+cyl.getl(id))/cyl.getl(cL.size()-1);
  240 + }
  241 + else if(id == 0)
  242 + {
  243 + pval = (cyl.getl(id+1)*result[i][2])/cyl.getl(cL.size()-1);
  244 + }
  245 + else
  246 + {
  247 + pval = (cyl.getl(id)/cyl.getl(cL.size()-1));
  248 + }
  249 + stim::vec<float> v = cyl.surf(pval, result[i][0]);
  250 + stim::vec<float> di = cyl.p(pval);
  251 + float rad = cyl.r(pval);
  252 + if(
  253 + !(v[0] > size[0] || v[1] > size[1]
  254 + || v[2] > size[2] || v[0] < 0
  255 + || v[1] < 0 || v[2] < 0))
  256 + {
  257 + setSeed(v);
  258 + setSeedVec((v-di).norm());
  259 + setSeedMag(rad);
  260 + }
  261 + }
  262 + }
  263 + }
  264 + }
  265 +
221 266
222 //--------------------------------------------------------------------------// 267 //--------------------------------------------------------------------------//
223 //---------------------TEMPLATE CREATION METHODS----------------------------// 268 //---------------------TEMPLATE CREATION METHODS----------------------------//
@@ -459,6 +504,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -459,6 +504,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
459 void 504 void
460 GenerateFBO(unsigned int width, unsigned int height, GLuint &textureID, GLuint &framebufferID) 505 GenerateFBO(unsigned int width, unsigned int height, GLuint &textureID, GLuint &framebufferID)
461 { 506 {
  507 + glDeleteFramebuffers(1, &framebufferID);
462 glGenFramebuffers(1, &framebufferID); 508 glGenFramebuffers(1, &framebufferID);
463 glBindFramebuffer(GL_FRAMEBUFFER, framebufferID); 509 glBindFramebuffer(GL_FRAMEBUFFER, framebufferID);
464 int numChannels = 1; 510 int numChannels = 1;
@@ -506,45 +552,60 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -506,45 +552,60 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
506 } 552 }
507 553
508 554
  555 + ///IF type == 0
509 ///Method for using the gl manipulation to align templates from 556 ///Method for using the gl manipulation to align templates from
510 ///Template space (-0.5 0.5) to Texture space (0.0, 1.0), 557 ///Template space (-0.5 0.5) to Texture space (0.0, 1.0),
511 ///Based on the p of the spider in real space (arbitrary). 558 ///Based on the p of the spider in real space (arbitrary).
  559 +
  560 + ///IF type == 1
  561 + ///Method for using the gl manipulation to set up a matrix
  562 + ///To transform from tissue space into texture space.
512 ///All transformation happen in glMatrixMode(GL_TEXTURE). 563 ///All transformation happen in glMatrixMode(GL_TEXTURE).
513 - void setMatrix() 564 + ///All transformation happen in glMatrixMode(GL_TEXTURE).
  565 + void setMatrix(int type = 0)
514 { 566 {
515 - float curTrans[16]; //array to store the matrix values.  
516 - stim::vec<float> rot = getRotation(d); //get the rotation parameters for the current direction vector.  
517 - glMatrixMode(GL_TEXTURE);  
518 - glLoadIdentity(); 567 + if(type == 0)
  568 + {
  569 + float curTrans[16]; //array to store the matrix values.
  570 + stim::vec<float> rot = getRotation(d); //get the rotation parameters for the current direction vector.
  571 + glMatrixMode(GL_TEXTURE);
  572 + glLoadIdentity();
519 573
520 - //Scale by the voxel size and number of slices.  
521 - glScalef(1.0/S[0]/R[0], 1.0/S[1]/R[1], 1.0/S[2]/R[2]);  
522 - //translate to the current position of the spider in the texture.  
523 - glTranslatef(p[0],  
524 - p[1],  
525 - p[2]);  
526 - //rotate to the current direction of the spider.  
527 - glRotatef(rot[0], rot[1], rot[2], rot[3]);  
528 - //scale to the magnitude of the spider.  
529 - glScalef(m[0],  
530 - m[0],  
531 - m[0]);  
532 - //get and store the current transformation matrix for later use.  
533 - glGetFloatv(GL_TEXTURE_MATRIX, curTrans);  
534 - cT.set(curTrans);  
535 - // printTransform();  
536 -  
537 - CHECK_OPENGL_ERROR  
538 - //revert back to default gl mode.  
539 - glMatrixMode(GL_MODELVIEW); 574 + //Scale by the voxel size and number of slices.
  575 + glScalef(1.0/S[0]/R[0], 1.0/S[1]/R[1], 1.0/S[2]/R[2]);
  576 + //translate to the current position of the spider in the texture.
  577 + glTranslatef(p[0],
  578 + p[1],
  579 + p[2]);
  580 + //rotate to the current direction of the spider.
  581 + glRotatef(rot[0], rot[1], rot[2], rot[3]);
  582 + //scale to the magnitude of the spider.
  583 + glScalef(m[0],
  584 + m[0],
  585 + m[0]);
  586 + //get and store the current transformation matrix for later use.
  587 + glGetFloatv(GL_TEXTURE_MATRIX, curTrans);
  588 + cT.set(curTrans);
  589 + // printTransform();
  590 +
  591 + CHECK_OPENGL_ERROR
  592 + //revert back to default gl mode.
  593 + glMatrixMode(GL_MODELVIEW);
  594 + }
  595 + else if(type == 1)
  596 + {
  597 + glMatrixMode(GL_TEXTURE);
  598 + glLoadIdentity();
  599 + glScalef(1.0/S[0]/R[0], 1.0/S[1]/R[1], 1.0/S[2]/R[2]);
  600 + glMatrixMode(GL_MODELVIEW);
  601 + }
540 } 602 }
541 603
542 ///Method for controling the buffer and texture binding. 604 ///Method for controling the buffer and texture binding.
543 ///Clears the buffer upon binding. 605 ///Clears the buffer upon binding.
544 void 606 void
545 - Bind() 607 + Bind(float len = 8.0)
546 { 608 {
547 - float len = 8.0;  
548 glBindFramebuffer(GL_FRAMEBUFFER, fboID);//set up GL buffer 609 glBindFramebuffer(GL_FRAMEBUFFER, fboID);//set up GL buffer
549 glFramebufferTexture2D( 610 glFramebufferTexture2D(
550 GL_FRAMEBUFFER, 611 GL_FRAMEBUFFER,
@@ -576,9 +637,8 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -576,9 +637,8 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
576 ///@param GLuint &framebufferID, framebuffer used for storage. 637 ///@param GLuint &framebufferID, framebuffer used for storage.
577 ///@param int nSamples, number of rectanges to create. 638 ///@param int nSamples, number of rectanges to create.
578 void 639 void
579 - Bind(GLuint &textureID, GLuint &framebufferID, int nSamples) 640 + Bind(GLuint &textureID, GLuint &framebufferID, int nSamples, float len = 8.0)
580 { 641 {
581 - float len = 8.0;  
582 glBindFramebuffer(GL_FRAMEBUFFER, framebufferID);//set up GL buffer 642 glBindFramebuffer(GL_FRAMEBUFFER, framebufferID);//set up GL buffer
583 glFramebufferTexture2D( 643 glFramebufferTexture2D(
584 GL_FRAMEBUFFER, 644 GL_FRAMEBUFFER,
@@ -1085,7 +1145,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1085,7 +1145,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1085 { 1145 {
1086 float x, y, z, u, v, w, m; 1146 float x, y, z, u, v, w, m;
1087 myfile >> x >> y >> z >> u >> v >> w >> m; 1147 myfile >> x >> y >> z >> u >> v >> w >> m;
1088 - setSeed(x, y , z); 1148 + setSeed(x, y, z);
1089 setSeedVec(u, v, w); 1149 setSeedVec(u, v, w);
1090 setSeedMag(m); 1150 setSeedMag(m);
1091 } 1151 }
@@ -1099,14 +1159,28 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1099,14 +1159,28 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1099 void 1159 void
1100 saveNetwork(std::string name) 1160 saveNetwork(std::string name)
1101 { 1161 {
  1162 + stim::glObj<float> sk;
  1163 + for(int i = 0; i < nt.sizeE(); i++)
  1164 + {
  1165 + std::vector<stim::vec< float > > cm = nt.getEdgeCenterLineMag(i);
  1166 + std::vector<stim::vec< float > > ce = nt.getEdgeCenterLine(i);
  1167 + sk.Begin(stim::OBJ_LINE);
  1168 + for(int j = 0; j < ce.size(); j++)
  1169 + {
  1170 + sk.TexCoord(cm[j][0]);
  1171 + sk.Vertex(ce[j][0], ce[j][1], ce[j][2]);
  1172 + }
  1173 + sk.End();
  1174 + }
1102 sk.save(name); 1175 sk.save(name);
1103 } 1176 }
1104 1177
  1178 + ///Depreciated, but might be reused later()
1105 ///returns a COPY of the entire stim::glObj object. 1179 ///returns a COPY of the entire stim::glObj object.
1106 stim::glObj<float> 1180 stim::glObj<float>
1107 getNetwork() 1181 getNetwork()
1108 { 1182 {
1109 - return sk; 1183 +// return sk;
1110 } 1184 }
1111 1185
1112 ///returns a COPY of the entire stim::glnetwork object. 1186 ///returns a COPY of the entire stim::glnetwork object.
@@ -1216,6 +1290,37 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1216,6 +1290,37 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1216 glEnd(); 1290 glEnd();
1217 glEndList(); 1291 glEndList();
1218 } 1292 }
  1293 +
  1294 +///need to return the cylinder.
  1295 + void
  1296 + DrawLongCylinder(int n = 8, int l_template = 8,int l_square = 8)
  1297 + {
  1298 + int cylLen = cL.size()-1;
  1299 + GenerateFBO(n*l_square, cylLen*l_template, btexbufferID, bfboID);
  1300 + Bind(btexbufferID, bfboID, cylLen, l_template*l_square/2.0);
  1301 + stim::cylinder<float> cyl(cL, cM);
  1302 + std::vector<std::vector<stim::vec<float> > > p = cyl.getPoints(n);
  1303 + for(int i = 0; i < p.size()-1; i++) ///number of circles
  1304 + {
  1305 + for(int j = 0; j < p[0].size()-1; j++) ///points in the circle
  1306 + {
  1307 + glBegin(GL_QUADS);
  1308 + glTexCoord3f(p[i][j][0], p[i][j][1], p[i][j][2]);
  1309 + glVertex2f(j*l_square, i*(float)l_template);
  1310 +
  1311 + glTexCoord3f(p[i][j+1][0], p[i][j+1][1], p[i][j+1][2]);
  1312 + glVertex2f(j*l_square+l_square, i*(float)l_template);
  1313 +
  1314 + glTexCoord3f(p[i+1][j+1][0], p[i+1][j+1][1], p[i+1][j+1][2]);
  1315 + glVertex2f(j*l_square+l_square, i*(float)l_template+(float)l_template);
  1316 +
  1317 + glTexCoord3f(p[i+1][j][0], p[i+1][j][1], p[i+1][j][2]);
  1318 + glVertex2f(j*l_square,i*(float)l_template+(float)l_template);
  1319 + glEnd();
  1320 + }
  1321 + }
  1322 + Unbind();
  1323 + }
1219 1324
1220 1325
1221 ///@param min_cost the cost value used for tracing 1326 ///@param min_cost the cost value used for tracing
@@ -1223,114 +1328,35 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1223,114 +1328,35 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1223 void 1328 void
1224 trace(int min_cost) 1329 trace(int min_cost)
1225 { 1330 {
1226 - Bind();  
1227 - rev = stim::vec<float>(0.0,0.0,1.0); 1331 +// rev = stim::vec<float>(0.0,0.0,1.0);
1228 bool sEmpty = true; 1332 bool sEmpty = true;
1229 float lastmag = 16.0;; 1333 float lastmag = 16.0;;
1230 - while(!seeds.empty()) 1334 + stim::vec<float> curSeed;
  1335 + stim::vec<float> curSeedVec;
  1336 + float curSeedMag;
  1337 + while(!Empty())
1231 { 1338 {
1232 //clear the currently traced line and start a new one. 1339 //clear the currently traced line and start a new one.
1233 cL.clear(); 1340 cL.clear();
1234 cM.clear(); 1341 cM.clear();
1235 - sk.Begin(stim::OBJ_LINE);  
1236 - stim::vec<float> curSeed = seeds.top();  
1237 -// std::cout << "The current seeds is " << curSeed << std::endl;  
1238 - stim::vec<float> curSeedVec = seedsvecs.top();  
1239 - float curSeedMag = seedsmags.top(); 1342 + cD.clear();
  1343 + curSeed = seeds.top();
  1344 + curSeedVec = seedsvecs.top();
  1345 + curSeedMag = seedsmags.top();
1240 seeds.pop(); 1346 seeds.pop();
1241 seedsvecs.pop(); 1347 seedsvecs.pop();
1242 seedsmags.pop(); 1348 seedsmags.pop();
1243 // std::cout << "The current seed Vector is " << curSeedVec << std::endl; 1349 // std::cout << "The current seed Vector is " << curSeedVec << std::endl;
1244 setPosition(curSeed); 1350 setPosition(curSeed);
1245 setDirection(curSeedVec); 1351 setDirection(curSeedVec);
1246 - cL.push_back(curSeed);  
1247 - cM.push_back(curSeedMag);  
1248 - sk.createFromSelf(GL_SELECT);  
1249 - traceLine(min_cost);  
1250 -  
1251 - sk.rev();  
1252 - // std::cout << "reversed" << std::endl;  
1253 - std::reverse(cL.begin(), cL.end());  
1254 - std::reverse(cM.begin(), cM.end());  
1255 - setPosition(curSeed);  
1256 - setDirection(-rev);  
1257 - setMagnitude(16.0);  
1258 - sk.createFromSelf(GL_SELECT);  
1259 - traceLine(min_cost);  
1260 - sk.End(); 1352 + setMagnitude(curSeedMag);
  1353 +// cL.push_back(curSeed);
  1354 +// cM.push_back(curSeedMag);
  1355 +// cD.push_back(curSeedMag);
  1356 + pair<stim::fiber<float>, int> a = traceLine(p, m, min_cost);
1261 } 1357 }
1262 - Unbind();  
1263 } 1358 }
1264 1359
1265 - ///@param min_cost the cost value used for tracing  
1266 - ///traces the seedpoint passed to completion in one directions.  
1267 - void  
1268 - traceLine(int min_cost)  
1269 - {  
1270 - stim::vec<float> pos;  
1271 - stim::vec<float> mag;  
1272 - int h;  
1273 - bool started = false;  
1274 - bool running = true;  
1275 - stim::vec<float> size(S[0]*R[0], S[1]*R[1], S[2]*R[2]);  
1276 - while(running)  
1277 - {  
1278 - int cost = Step();  
1279 - if (cost > min_cost){  
1280 - running = false;  
1281 - break;  
1282 - } else {  
1283 - //Have we found an edge?  
1284 - pos = getPosition();  
1285 - if(pos[0] > size[0] || pos[1] > size[1]  
1286 - || pos[2] > size[2] || pos[0] < 0  
1287 - || pos[1] < 0 || pos[2] < 0)  
1288 - {  
1289 -// std::cout << "Found Edge" << std::endl;  
1290 - running = false;  
1291 - break;  
1292 - }  
1293 - //If this is the first step in the trace,  
1294 - // save the direction  
1295 - //(to be used later to trace the fiber in the opposite direction)  
1296 - if(started == false){  
1297 - rev = -getDirection();  
1298 - started = true;  
1299 - }  
1300 -// std::cout << i << p << std::endl;  
1301 - m = getMagnitude();  
1302 - //Has the template size gotten unreasonable?  
1303 - if(m[0] > 75 || m[0] < 1){  
1304 -// std::cout << "Magnitude Limit" << std::endl;  
1305 - running = false;  
1306 - break;  
1307 - }  
1308 - else  
1309 - {  
1310 - h = selectObject(pos, getDirection(), m[0]);  
1311 - //Have we hit something previously traced?  
1312 - if(h != -1){  
1313 - std::cout << "I hit a line" << h << std::endl;  
1314 - running = false;  
1315 - break;  
1316 - }  
1317 - else {  
1318 - cL.push_back(stim::vec<float>(p[0], p[1],p[2]));//  
1319 - sk.TexCoord(m[0]);  
1320 - sk.Vertex(p[0], p[1], p[2]);  
1321 - Bind(btexbufferID, bfboID, 27);  
1322 - CHECK_OPENGL_ERROR  
1323 - branchDetection();  
1324 - CHECK_OPENGL_ERROR  
1325 - Unbind();  
1326 - CHECK_OPENGL_ERROR  
1327 - }  
1328 - }  
1329 - }  
1330 - }  
1331 - }  
1332 -  
1333 -  
1334 int 1360 int
1335 selectObject(stim::vec<float> loc, stim::vec<float> dir, float mag) 1361 selectObject(stim::vec<float> loc, stim::vec<float> dir, float mag)
1336 { 1362 {
@@ -1468,10 +1494,10 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1468,10 +1494,10 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1468 spos[2] = spos[2]-sdir[2]*smag[0]/2.; 1494 spos[2] = spos[2]-sdir[2]*smag[0]/2.;
1469 int h = selectObject(spos, -sdir, smag[0]); 1495 int h = selectObject(spos, -sdir, smag[0]);
1470 //did start with a fiber? 1496 //did start with a fiber?
1471 - if(h != -1){ 1497 + if(h != -1 && h != nt.sizeE()){
1472 // std::cout << "got here double" << smag.str() << std::endl; 1498 // std::cout << "got here double" << smag.str() << std::endl;
1473 nt.addEdge(ce,cm, h, in.second); 1499 nt.addEdge(ce,cm, h, in.second);
1474 - } 1500 + } else { nt.addEdge(ce,cm, -1, -1);}
1475 } 1501 }
1476 } 1502 }
1477 } 1503 }
@@ -1494,7 +1520,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1494,7 +1520,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1494 stim::vec<float> sdir = getDirection(); 1520 stim::vec<float> sdir = getDirection();
1495 1521
1496 Bind(); 1522 Bind();
1497 - sk.Begin(stim::OBJ_LINE); 1523 +// sk.Begin(stim::OBJ_LINE);
1498 1524
1499 1525
1500 // sk.createFromSelf(GL_SELECT); 1526 // sk.createFromSelf(GL_SELECT);
@@ -1514,7 +1540,8 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1514,7 +1540,8 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1514 int cost = Step(); 1540 int cost = Step();
1515 if (cost > min_cost){ 1541 if (cost > min_cost){
1516 running = false; 1542 running = false;
1517 - sk.End(); 1543 +// sk.End();
  1544 + branchDetection2();
1518 pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), -1); 1545 pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), -1);
1519 addToNetwork(a, spos, smag, sdir); 1546 addToNetwork(a, spos, smag, sdir);
1520 return a; 1547 return a;
@@ -1526,9 +1553,8 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1526,9 +1553,8 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1526 || pos[2] > size[2] || pos[0] < 0 1553 || pos[2] > size[2] || pos[0] < 0
1527 || pos[1] < 0 || pos[2] < 0) 1554 || pos[1] < 0 || pos[2] < 0)
1528 { 1555 {
1529 -// std::cout << "Found Edge" << std::endl;  
1530 running = false; 1556 running = false;
1531 - sk.End(); 1557 + branchDetection2();
1532 pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), -1); 1558 pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), -1);
1533 addToNetwork(a, spos, smag, sdir); 1559 addToNetwork(a, spos, smag, sdir);
1534 return a; 1560 return a;
@@ -1541,13 +1567,11 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1541,13 +1567,11 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1541 rev = -getDirection(); 1567 rev = -getDirection();
1542 started = true; 1568 started = true;
1543 } 1569 }
1544 -// std::cout << i << p << std::endl;  
1545 //Has the template size gotten unreasonable? 1570 //Has the template size gotten unreasonable?
1546 mag = getMagnitude(); 1571 mag = getMagnitude();
1547 if(mag[0] > 75 || mag[0] < 1){ 1572 if(mag[0] > 75 || mag[0] < 1){
1548 -// std::cout << "Magnitude Limit" << std::endl;  
1549 running = false; 1573 running = false;
1550 - sk.End(); 1574 + branchDetection2();
1551 pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), -1); 1575 pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), -1);
1552 addToNetwork(a, spos, smag, sdir); 1576 addToNetwork(a, spos, smag, sdir);
1553 return a; 1577 return a;
@@ -1559,7 +1583,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1559,7 +1583,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1559 //Have we hit something previously traced? 1583 //Have we hit something previously traced?
1560 if(h != -1){ 1584 if(h != -1){
1561 running = false; 1585 running = false;
1562 - sk.End(); 1586 + branchDetection2();
1563 pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), h); 1587 pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), h);
1564 addToNetwork(a, spos, smag, sdir); 1588 addToNetwork(a, spos, smag, sdir);
1565 return a; 1589 return a;
@@ -1568,14 +1592,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1568,14 +1592,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1568 else { 1592 else {
1569 cL.push_back(stim::vec<float>(p[0], p[1],p[2])); 1593 cL.push_back(stim::vec<float>(p[0], p[1],p[2]));
1570 cM.push_back(stim::vec<float>(m[0], m[0])); 1594 cM.push_back(stim::vec<float>(m[0], m[0]));
1571 -// cM.push_back(m[0]);  
1572 -  
1573 - sk.TexCoord(m[0]);  
1574 - sk.Vertex(p[0], p[1], p[2]);  
1575 - Bind(btexbufferID, bfboID, 27);  
1576 - CHECK_OPENGL_ERROR  
1577 - branchDetection();  
1578 - CHECK_OPENGL_ERROR 1595 +// Bind(btexbufferID, bfboID, 27);
1579 Unbind(); 1596 Unbind();
1580 CHECK_OPENGL_ERROR 1597 CHECK_OPENGL_ERROR
1581 1598
stim/math/circle.h
1 #ifndef STIM_CIRCLE_H 1 #ifndef STIM_CIRCLE_H
2 #define STIM_CIRCLE_H 2 #define STIM_CIRCLE_H
3 3
4 -//enable CUDA_CALLABLE macro  
5 #include <stim/cuda/cudatools/callable.h> 4 #include <stim/cuda/cudatools/callable.h>
  5 +#include <stim/math/plane.h>
6 #include <stim/math/vector.h> 6 #include <stim/math/vector.h>
7 #include <stim/math/triangle.h> 7 #include <stim/math/triangle.h>
8 -#include <stim/math/quaternion.h>  
9 -#include <stim/math/rect.h>  
10 -#include <iostream>  
11 -#include <iomanip> 8 +#include <assert.h>
12 #include <algorithm> 9 #include <algorithm>
  10 +#include <iostream>
13 11
14 -namespace stim  
15 -{ 12 +namespace stim{
16 13
17 -template <class T>  
18 -struct circle : rect<T> 14 +template <typename T>
  15 +class circle : plane<T>
19 { 16 {
20 - private:  
21 - T theta;  
22 -  
23 - public:  
24 -  
25 - using stim::rect<T>::p;  
26 - using stim::rect<T>::normal;  
27 - using stim::rect<T>::center;  
28 - using stim::rect<T>::scale;  
29 - ///base constructor  
30 - ///@param th value of the angle of the starting point from 0 to 360.  
31 - CUDA_CALLABLE circle(float th = 0.0) : rect<T>()  
32 - {  
33 - theta = th;  
34 - }  
35 17
36 - ///create a rectangle given a size and position in Z space.  
37 - ///@param size: size of the rectangle in ND space.  
38 - ///@param z_pos z coordinate of the rectangle.  
39 - ///@param th value of the angle of the starting point from 0 to 360.  
40 - CUDA_CALLABLE circle(T size, T zpos = (T)0, float th = 0.0) : rect<T>(size, zpos)  
41 - {  
42 - theta = th;