Compare View
Commits (11)
Showing
7 changed files
Show diff stats
python/structen.py
... | ... | @@ -226,7 +226,7 @@ def st2amira(filename, T): |
226 | 226 | A[..., 2] = T[..., 3] |
227 | 227 | |
228 | 228 | #roll the tensor axis so that it is the leading component |
229 | - A = numpy.rollaxis(A, A.ndim - 1) | |
229 | + #A = numpy.rollaxis(A, A.ndim - 1) | |
230 | 230 | A.tofile(filename) |
231 | 231 | print("\n", A.shape) |
232 | 232 | ... | ... |
stim/biomodels/network.h
... | ... | @@ -14,22 +14,23 @@ |
14 | 14 | #include <stim/visualization/cylinder.h> |
15 | 15 | #include <stim/structures/kdtree.cuh> |
16 | 16 | #include <stim/cuda/cudatools/timer.h> |
17 | +#include <stim/cuda/cudatools/callable.h> | |
17 | 18 | |
19 | +//********************help function******************** | |
20 | +// gaussian_function | |
21 | +CUDA_CALLABLE float gaussianFunction(float x, float std = 25) { return exp(-x / (2 * std*std)); } // std default sigma value is 25 | |
18 | 22 | |
23 | +// compute metric in parallel | |
19 | 24 | #ifdef __CUDACC__ |
20 | -//device gaussian function | |
21 | -__device__ float gaussianFunction(float x, float std){ return expf(-x/(2*std*std));} | |
22 | - | |
23 | -//compute metric in parallel | |
24 | 25 | template <typename T> |
25 | -__global__ void d_metric(T* M, size_t n, T* D, float sigma){ | |
26 | +__global__ void find_metric_parallel(T* M, size_t n, T* D, float sigma){ | |
26 | 27 | size_t x = blockDim.x * blockIdx.x + threadIdx.x; |
27 | 28 | if(x >= n) return; |
28 | 29 | M[x] = 1.0f - gaussianFunction(D[x], sigma); |
29 | 30 | } |
30 | 31 | |
31 | 32 | //find the corresponding edge index from array index |
32 | -__global__ void d_findedge(size_t* I, size_t n, unsigned* R, size_t* E, size_t ne){ | |
33 | +__global__ void find_edge_index_parallel(size_t* I, size_t n, unsigned* R, size_t* E, size_t ne){ | |
33 | 34 | size_t x = blockDim.x * blockIdx.x + threadIdx.x; |
34 | 35 | if(x >= n) return; |
35 | 36 | unsigned i = 0; |
... | ... | @@ -46,13 +47,13 @@ __global__ void d_findedge(size_t* I, size_t n, unsigned* R, size_t* E, size_t n |
46 | 47 | #endif |
47 | 48 | |
48 | 49 | //hard-coded factor |
49 | -int threshold_fac = 0; | |
50 | +int threshold_fac; | |
50 | 51 | |
51 | 52 | namespace stim{ |
52 | 53 | /** This is the a class that interfaces with gl_spider in order to store the currently |
53 | 54 | * segmented network. The following data is stored and can be extracted: |
54 | 55 | * 1)Network geometry and centerline. |
55 | - * 2)Network connectivity (a graph of nodes and edges), reconstructed using ANN library. | |
56 | + * 2)Network connectivity (a graph of nodes and edges), reconstructed using kdtree. | |
56 | 57 | */ |
57 | 58 | |
58 | 59 | template<typename T> |
... | ... | @@ -136,15 +137,12 @@ class network{ |
136 | 137 | |
137 | 138 | return ss.str(); |
138 | 139 | } |
139 | - | |
140 | - | |
141 | 140 | }; |
142 | 141 | |
143 | 142 | protected: |
144 | 143 | |
145 | 144 | std::vector<edge> E; //list of edges |
146 | 145 | std::vector<vertex> V; //list of vertices. |
147 | - //std::vector<int> NT; //stores a list of neuronal type for each point in the centerline (will set value only when centerline is built from swc file) | |
148 | 146 | |
149 | 147 | public: |
150 | 148 | |
... | ... | @@ -171,12 +169,12 @@ public: |
171 | 169 | } |
172 | 170 | |
173 | 171 | ///Returns the radius at specific point in the edge |
174 | - T get_R(unsigned e, unsigned i) { | |
172 | + T get_r(unsigned e, unsigned i) { | |
175 | 173 | return E[e].r(i); |
176 | 174 | } |
177 | 175 | |
178 | 176 | ///Returns the average radius of specific edge |
179 | - T get_average_R(unsigned e) { | |
177 | + T get_average_r(unsigned e) { | |
180 | 178 | T result = 0.0; |
181 | 179 | unsigned n = E[e].size(); |
182 | 180 | for (unsigned p = 0; p < n; p++) |
... | ... | @@ -186,7 +184,7 @@ public: |
186 | 184 | } |
187 | 185 | |
188 | 186 | ///Returns the length of current edge |
189 | - T get_L(unsigned e) { | |
187 | + T get_l(unsigned e) { | |
190 | 188 | return E[e].length(); |
191 | 189 | } |
192 | 190 | |
... | ... | @@ -219,12 +217,13 @@ public: |
219 | 217 | } |
220 | 218 | |
221 | 219 | ///Set radius |
222 | - void set_R(unsigned e, std::vector<T> radius) { | |
220 | + void set_r(unsigned e, std::vector<T> radius) { | |
223 | 221 | E[e].cylinder<T>::copy_r(radius); |
224 | 222 | } |
225 | 223 | |
226 | - void set_R(unsigned e, unsigned i, T radius) { | |
227 | - E[e].cylinder<T>::set_r(i, radius); | |
224 | + void set_r(unsigned e, T radius) { | |
225 | + for (size_t i = 0; i < E[e].size(); i++) | |
226 | + E[e].cylinder<T>::set_r(i, radius); | |
228 | 227 | } |
229 | 228 | //scale the network by some constant value |
230 | 229 | // I don't think these work?????? |
... | ... | @@ -580,9 +579,6 @@ public: |
580 | 579 | return n; //return the resampled network |
581 | 580 | } |
582 | 581 | |
583 | - // host gaussian function | |
584 | - __host__ float gaussianFunction(float x, float std = 25){ return exp(-x/(2*std*std));} // std default value is 25 | |
585 | - | |
586 | 582 | /// Calculate the total number of points on all edges. |
587 | 583 | unsigned total_points(){ |
588 | 584 | unsigned n = 0; |
... | ... | @@ -624,6 +620,19 @@ public: |
624 | 620 | } |
625 | 621 | } |
626 | 622 | |
623 | + // get list of metric | |
624 | + std::vector<T> metric() { | |
625 | + std::vector<T> result; | |
626 | + T m; | |
627 | + for (size_t e = 0; e < E.size(); e++) { | |
628 | + for (size_t p = 0; p < E[e].size(); p++) { | |
629 | + m = E[e].r(p); | |
630 | + result.push_back(m); | |
631 | + } | |
632 | + } | |
633 | + return result; | |
634 | + } | |
635 | + | |
627 | 636 | /// Calculate the average magnitude across the entire network. |
628 | 637 | /// @param m is the magnitude value to use. The default is 0 (usually radius). |
629 | 638 | T average(unsigned m = 0){ |
... | ... | @@ -694,7 +703,7 @@ public: |
694 | 703 | size_t threads = (1024>n)?n:1024; |
695 | 704 | size_t blocks = n/threads + (n%threads)?1:0; |
696 | 705 | |
697 | - d_metric<<<blocks, threads>>>(d_m1, n, d_dists, sigma); //calculate the metric value based on the distance | |
706 | + find_metric_parallel<<<blocks, threads>>>(d_m1, n, d_dists, sigma); //calculate the metric value based on the distance | |
698 | 707 | |
699 | 708 | cudaMemcpy(m1, d_m1, n * sizeof(T), cudaMemcpyDeviceToHost); |
700 | 709 | |
... | ... | @@ -717,9 +726,9 @@ public: |
717 | 726 | T* dists = new T[n]; |
718 | 727 | size* nnIdx = new size_t[n]; |
719 | 728 | |
720 | - edge2array(queryPt, R.E[e]); | |
729 | + edge2array(query, R.E[e]); | |
721 | 730 | |
722 | - kdt.cpu_search(queryPt, n, nnIdx, dists); //find the distance between A and the current network | |
731 | + kdt.cpu_search(query, n, nnIdx, dists); //find the distance between A and the current network | |
723 | 732 | |
724 | 733 | for(unsigned p = 0; p < R.E[e].size(); p++){ |
725 | 734 | m1[p] = 1.0f - gaussianFunction((T)dists[p], sigma); //calculate the metric value based on the distance |
... | ... | @@ -795,7 +804,7 @@ public: |
795 | 804 | size_t threads = (1024>n)?n:1024; // test to see whether the number of point in current edge is more than 1024 |
796 | 805 | size_t blocks = n/threads + (n%threads)?1:0; |
797 | 806 | |
798 | - d_metric<<<blocks, threads>>>(d_m1, n, d_dists, sigma); // calculate the metrics in parallel | |
807 | + find_metric_parallel<<<blocks, threads>>>(d_m1, n, d_dists, sigma); // calculate the metrics in parallel | |
799 | 808 | |
800 | 809 | cudaMemcpy(m1, d_m1, n * sizeof(T), cudaMemcpyDeviceToHost); |
801 | 810 | |
... | ... | @@ -816,7 +825,7 @@ public: |
816 | 825 | cudaMalloc((void**)&d_edge_point_num, NB * sizeof(size_t)); |
817 | 826 | cudaMemcpy(d_edge_point_num, &edge_point_num[0], NB * sizeof(size_t), cudaMemcpyHostToDevice); |
818 | 827 | |
819 | - d_findedge<<<blocks, threads>>>(d_nnIdx, n, d_relation, d_edge_point_num, NB); // find the edge corresponding to the array index in parallel | |
828 | + find_edge_index_parallel<<<blocks, threads>>>(d_nnIdx, n, d_relation, d_edge_point_num, NB); // find the edge corresponding to the array index in parallel | |
820 | 829 | |
821 | 830 | cudaMemcpy(&relation[e][0], d_relation, n * sizeof(unsigned), cudaMemcpyDeviceToHost); //copy relationship from device to host |
822 | 831 | } | ... | ... |
stim/cuda/cudatools/callable.h
stim/iVote/ivote2/local_max.cuh
... | ... | @@ -10,30 +10,22 @@ namespace stim{ |
10 | 10 | |
11 | 11 | // this kernel calculates the local maximum for finding the cell centers |
12 | 12 | template<typename T> |
13 | - __global__ void cuda_local_max(T* gpuCenters, T* gpuVote, int conn, size_t x, size_t y){ | |
13 | + __global__ void cuda_local_max(T* gpuCenters, T* gpuVote, int conn, int x, int y){ | |
14 | 14 | |
15 | 15 | // calculate the 2D coordinates for this current thread. |
16 | - size_t xi = blockIdx.x * blockDim.x + threadIdx.x; | |
17 | - size_t yi = blockIdx.y * blockDim.y + threadIdx.y; | |
16 | + int xi = blockIdx.x * blockDim.x + threadIdx.x; | |
17 | + int yi = blockIdx.y * blockDim.y + threadIdx.y; | |
18 | 18 | |
19 | 19 | if(xi >= x || yi >= y) |
20 | 20 | return; |
21 | 21 | |
22 | 22 | // convert 2D coordinates to 1D |
23 | - size_t i = yi * x + xi; | |
23 | + int i = yi * x + xi; | |
24 | 24 | |
25 | 25 | gpuCenters[i] = 0; //initialize the value at this location to zero |
26 | 26 | |
27 | 27 | T val = gpuVote[i]; |
28 | 28 | |
29 | - //compare to the threshold | |
30 | - //if(val < final_t) return; | |
31 | - | |
32 | - //define an array to store indices with same vote value | |
33 | - /*int * IdxEq; | |
34 | - IdxEq = new int [2*conn]; | |
35 | - int n = 0;*/ | |
36 | - | |
37 | 29 | for(int xl = xi - conn; xl < xi + conn; xl++){ |
38 | 30 | for(int yl = yi - conn; yl < yi + conn; yl++){ |
39 | 31 | if(xl >= 0 && xl < x && yl >= 0 && yl < y){ |
... | ... | @@ -42,8 +34,7 @@ namespace stim{ |
42 | 34 | return; |
43 | 35 | } |
44 | 36 | if (gpuVote[il] == val){ |
45 | - /*IdxEq[n] = il; | |
46 | - n = n+1;*/ | |
37 | + | |
47 | 38 | if( il > i){ |
48 | 39 | return; |
49 | 40 | } |
... | ... | @@ -51,12 +42,7 @@ namespace stim{ |
51 | 42 | } |
52 | 43 | } |
53 | 44 | } |
54 | - /*if (n!=0){ | |
55 | - if(IdxEq[n/2] !=i){ | |
56 | - return; | |
57 | - } | |
58 | - } */ | |
59 | - //gpuCenters[i] = 1; | |
45 | + | |
60 | 46 | gpuCenters[i] = gpuVote[i]; |
61 | 47 | } |
62 | 48 | ... | ... |
stim/math/vec3.h
... | ... | @@ -5,6 +5,8 @@ |
5 | 5 | #include <stim/cuda/cudatools/callable.h> |
6 | 6 | #include <cmath> |
7 | 7 | |
8 | +#include <sstream> | |
9 | + | |
8 | 10 | |
9 | 11 | namespace stim{ |
10 | 12 | |
... | ... | @@ -243,9 +245,9 @@ public: |
243 | 245 | return false; |
244 | 246 | } |
245 | 247 | |
246 | -#ifndef __NVCC__ | |
248 | +//#ifndef __CUDACC__ | |
247 | 249 | /// Outputs the vector as a string |
248 | - std::string str() const{ | |
250 | +std::string str() const{ | |
249 | 251 | std::stringstream ss; |
250 | 252 | |
251 | 253 | const size_t N = 3; |
... | ... | @@ -261,7 +263,7 @@ public: |
261 | 263 | |
262 | 264 | return ss.str(); |
263 | 265 | } |
264 | -#endif | |
266 | +//#endif | |
265 | 267 | |
266 | 268 | size_t size(){ return 3; } |
267 | 269 | ... | ... |
stim/visualization/colormap.h
... | ... | @@ -166,14 +166,12 @@ static void gpu2gpu(T* gpuSource, unsigned char* gpuDest, unsigned int nVals, T |
166 | 166 | gridX = 65535; |
167 | 167 | } |
168 | 168 | dim3 dimGrid(gridX, gridY); |
169 | - //int gridDim = (nVals + blockDim - 1)/blockDim; | |
170 | 169 | if(cm == cmGrayscale) |
171 | 170 | applyGrayscale<<<dimGrid, blockDim>>>(gpuSource, gpuDest, nVals, minVal, maxVal); |
172 | 171 | else if(cm == cmBrewer) |
173 | 172 | { |
174 | 173 | initBrewer(); |
175 | 174 | applyBrewer<<<dimGrid, blockDim>>>(gpuSource, gpuDest, nVals, minVal, maxVal); |
176 | - //HANDLE_ERROR(cudaMemset(gpuDest, 0, sizeof(unsigned char) * nVals * 3)); | |
177 | 175 | destroyBrewer(); |
178 | 176 | } |
179 | 177 | |
... | ... | @@ -190,13 +188,9 @@ static void gpu2cpu(T* gpuSource, unsigned char* cpuDest, unsigned int nVals, T |
190 | 188 | unsigned char* gpuDest; |
191 | 189 | HANDLE_ERROR(cudaMalloc( (void**)&gpuDest, sizeof(unsigned char) * nVals * 3 )); |
192 | 190 | |
193 | - //HANDLE_ERROR(cudaMemset(gpuSource, 0, sizeof(T) * nVals)); | |
194 | - | |
195 | 191 | //create the image on the gpu |
196 | 192 | gpu2gpu(gpuSource, gpuDest, nVals, minVal, maxVal, cm); |
197 | - | |
198 | - //HANDLE_ERROR(cudaMemset(gpuDest, 0, sizeof(unsigned char) * nVals * 3)); | |
199 | - | |
193 | + | |
200 | 194 | //copy the image from the GPU to the CPU |
201 | 195 | HANDLE_ERROR(cudaMemcpy(cpuDest, gpuDest, sizeof(unsigned char) * nVals * 3, cudaMemcpyDeviceToHost)); |
202 | 196 | ... | ... |
stim/visualization/gl_network.h
... | ... | @@ -158,17 +158,18 @@ public: |
158 | 158 | ///render the network cylinder as a series of tubes(when only one network loaded) |
159 | 159 | void glCylinder0() { |
160 | 160 | |
161 | + float r1, r2; | |
161 | 162 | if (!glIsList(dlist)) { // if dlist isn't a display list, create it |
162 | 163 | dlist = glGenLists(1); // generate a display list |
163 | 164 | glNewList(dlist, GL_COMPILE); // start a new display list |
164 | - T radius = 0.0; | |
165 | 165 | for (unsigned e = 0; e < E.size(); e++) { // for each edge in the network |
166 | - radius = get_average_R(e) / get_average_R(0); // calculate the radius | |
167 | 166 | for (unsigned p = 1; p < E[e].size(); p++) { // for each point on that edge |
168 | 167 | stim::circle<T> C1 = E[e].circ(p - 1); |
169 | 168 | stim::circle<T> C2 = E[e].circ(p); |
170 | - C1.set_R(2 * radius); // re-scale the circle to the same | |
171 | - C2.set_R(2 * radius); | |
169 | + r1 = E[e].r(p - 1); | |
170 | + r2 = E[e].r(p); | |
171 | + C1.set_R(2 * r1); // re-scale the circle to the same | |
172 | + C2.set_R(2 * r2); | |
172 | 173 | std::vector< stim::vec3<T> > Cp1 = C1.glpoints(20);// get 20 points on the circle plane |
173 | 174 | std::vector< stim::vec3<T> > Cp2 = C2.glpoints(20); |
174 | 175 | glBegin(GL_QUAD_STRIP); |
... | ... | @@ -180,8 +181,17 @@ public: |
180 | 181 | } // set the texture coordinate based on the specified magnitude index |
181 | 182 | } |
182 | 183 | for (unsigned n = 0; n < V.size(); n++) { |
183 | - size_t num = V[n].e[0].size(); // store the number of outgoing edge of that vertex | |
184 | - renderBall(V[n][0], V[n][1], V[n][2], 2 * radius, 20); | |
184 | + for (unsigned i = 0; i < E.size(); i++) { | |
185 | + if (E[i].v[0] == n) { | |
186 | + r1 = E[i].r(0); | |
187 | + break; | |
188 | + } | |
189 | + else if (E[i].v[1] == n) { | |
190 | + r1 = E[i].r(E[i].size() - 1); | |
191 | + break; | |
192 | + } | |
193 | + } | |
194 | + renderBall(V[n][0], V[n][1], V[n][2], 2 * r1, 20); | |
185 | 195 | } |
186 | 196 | glEndList(); // end the display list |
187 | 197 | } | ... | ... |