Commit 5068402b75ea90260c5640544c339e565034671c
1 parent
72175dcf
add function to access protected parameters outside the class
Showing
2 changed files
with
50 additions
and
31 deletions
Show diff stats
stim/biomodels/network.h
@@ -14,22 +14,23 @@ | @@ -14,22 +14,23 @@ | ||
14 | #include <stim/visualization/cylinder.h> | 14 | #include <stim/visualization/cylinder.h> |
15 | #include <stim/structures/kdtree.cuh> | 15 | #include <stim/structures/kdtree.cuh> |
16 | #include <stim/cuda/cudatools/timer.h> | 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 | #ifdef __CUDACC__ | 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 | template <typename T> | 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 | size_t x = blockDim.x * blockIdx.x + threadIdx.x; | 27 | size_t x = blockDim.x * blockIdx.x + threadIdx.x; |
27 | if(x >= n) return; | 28 | if(x >= n) return; |
28 | M[x] = 1.0f - gaussianFunction(D[x], sigma); | 29 | M[x] = 1.0f - gaussianFunction(D[x], sigma); |
29 | } | 30 | } |
30 | 31 | ||
31 | //find the corresponding edge index from array index | 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 | size_t x = blockDim.x * blockIdx.x + threadIdx.x; | 34 | size_t x = blockDim.x * blockIdx.x + threadIdx.x; |
34 | if(x >= n) return; | 35 | if(x >= n) return; |
35 | unsigned i = 0; | 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,13 +47,13 @@ __global__ void d_findedge(size_t* I, size_t n, unsigned* R, size_t* E, size_t n | ||
46 | #endif | 47 | #endif |
47 | 48 | ||
48 | //hard-coded factor | 49 | //hard-coded factor |
49 | -int threshold_fac = 0; | 50 | +int threshold_fac; |
50 | 51 | ||
51 | namespace stim{ | 52 | namespace stim{ |
52 | /** This is the a class that interfaces with gl_spider in order to store the currently | 53 | /** This is the a class that interfaces with gl_spider in order to store the currently |
53 | * segmented network. The following data is stored and can be extracted: | 54 | * segmented network. The following data is stored and can be extracted: |
54 | * 1)Network geometry and centerline. | 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 | template<typename T> | 59 | template<typename T> |
@@ -136,15 +137,12 @@ class network{ | @@ -136,15 +137,12 @@ class network{ | ||
136 | 137 | ||
137 | return ss.str(); | 138 | return ss.str(); |
138 | } | 139 | } |
139 | - | ||
140 | - | ||
141 | }; | 140 | }; |
142 | 141 | ||
143 | protected: | 142 | protected: |
144 | 143 | ||
145 | std::vector<edge> E; //list of edges | 144 | std::vector<edge> E; //list of edges |
146 | std::vector<vertex> V; //list of vertices. | 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 | public: | 147 | public: |
150 | 148 | ||
@@ -171,12 +169,12 @@ public: | @@ -171,12 +169,12 @@ public: | ||
171 | } | 169 | } |
172 | 170 | ||
173 | ///Returns the radius at specific point in the edge | 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 | return E[e].r(i); | 173 | return E[e].r(i); |
176 | } | 174 | } |
177 | 175 | ||
178 | ///Returns the average radius of specific edge | 176 | ///Returns the average radius of specific edge |
179 | - T get_average_R(unsigned e) { | 177 | + T get_average_r(unsigned e) { |
180 | T result = 0.0; | 178 | T result = 0.0; |
181 | unsigned n = E[e].size(); | 179 | unsigned n = E[e].size(); |
182 | for (unsigned p = 0; p < n; p++) | 180 | for (unsigned p = 0; p < n; p++) |
@@ -186,7 +184,7 @@ public: | @@ -186,7 +184,7 @@ public: | ||
186 | } | 184 | } |
187 | 185 | ||
188 | ///Returns the length of current edge | 186 | ///Returns the length of current edge |
189 | - T get_L(unsigned e) { | 187 | + T get_l(unsigned e) { |
190 | return E[e].length(); | 188 | return E[e].length(); |
191 | } | 189 | } |
192 | 190 | ||
@@ -219,12 +217,13 @@ public: | @@ -219,12 +217,13 @@ public: | ||
219 | } | 217 | } |
220 | 218 | ||
221 | ///Set radius | 219 | ///Set radius |
222 | - void set_R(unsigned e, std::vector<T> radius) { | 220 | + void set_r(unsigned e, std::vector<T> radius) { |
223 | E[e].cylinder<T>::copy_r(radius); | 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 | //scale the network by some constant value | 228 | //scale the network by some constant value |
230 | // I don't think these work?????? | 229 | // I don't think these work?????? |
@@ -580,9 +579,6 @@ public: | @@ -580,9 +579,6 @@ public: | ||
580 | return n; //return the resampled network | 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 | /// Calculate the total number of points on all edges. | 582 | /// Calculate the total number of points on all edges. |
587 | unsigned total_points(){ | 583 | unsigned total_points(){ |
588 | unsigned n = 0; | 584 | unsigned n = 0; |
@@ -624,6 +620,19 @@ public: | @@ -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 | /// Calculate the average magnitude across the entire network. | 636 | /// Calculate the average magnitude across the entire network. |
628 | /// @param m is the magnitude value to use. The default is 0 (usually radius). | 637 | /// @param m is the magnitude value to use. The default is 0 (usually radius). |
629 | T average(unsigned m = 0){ | 638 | T average(unsigned m = 0){ |
@@ -694,7 +703,7 @@ public: | @@ -694,7 +703,7 @@ public: | ||
694 | size_t threads = (1024>n)?n:1024; | 703 | size_t threads = (1024>n)?n:1024; |
695 | size_t blocks = n/threads + (n%threads)?1:0; | 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 | cudaMemcpy(m1, d_m1, n * sizeof(T), cudaMemcpyDeviceToHost); | 708 | cudaMemcpy(m1, d_m1, n * sizeof(T), cudaMemcpyDeviceToHost); |
700 | 709 | ||
@@ -717,9 +726,9 @@ public: | @@ -717,9 +726,9 @@ public: | ||
717 | T* dists = new T[n]; | 726 | T* dists = new T[n]; |
718 | size* nnIdx = new size_t[n]; | 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 | for(unsigned p = 0; p < R.E[e].size(); p++){ | 733 | for(unsigned p = 0; p < R.E[e].size(); p++){ |
725 | m1[p] = 1.0f - gaussianFunction((T)dists[p], sigma); //calculate the metric value based on the distance | 734 | m1[p] = 1.0f - gaussianFunction((T)dists[p], sigma); //calculate the metric value based on the distance |
@@ -795,7 +804,7 @@ public: | @@ -795,7 +804,7 @@ public: | ||
795 | size_t threads = (1024>n)?n:1024; // test to see whether the number of point in current edge is more than 1024 | 804 | size_t threads = (1024>n)?n:1024; // test to see whether the number of point in current edge is more than 1024 |
796 | size_t blocks = n/threads + (n%threads)?1:0; | 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 | cudaMemcpy(m1, d_m1, n * sizeof(T), cudaMemcpyDeviceToHost); | 809 | cudaMemcpy(m1, d_m1, n * sizeof(T), cudaMemcpyDeviceToHost); |
801 | 810 | ||
@@ -816,7 +825,7 @@ public: | @@ -816,7 +825,7 @@ public: | ||
816 | cudaMalloc((void**)&d_edge_point_num, NB * sizeof(size_t)); | 825 | cudaMalloc((void**)&d_edge_point_num, NB * sizeof(size_t)); |
817 | cudaMemcpy(d_edge_point_num, &edge_point_num[0], NB * sizeof(size_t), cudaMemcpyHostToDevice); | 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 | cudaMemcpy(&relation[e][0], d_relation, n * sizeof(unsigned), cudaMemcpyDeviceToHost); //copy relationship from device to host | 830 | cudaMemcpy(&relation[e][0], d_relation, n * sizeof(unsigned), cudaMemcpyDeviceToHost); //copy relationship from device to host |
822 | } | 831 | } |
stim/visualization/gl_network.h
@@ -158,17 +158,18 @@ public: | @@ -158,17 +158,18 @@ public: | ||
158 | ///render the network cylinder as a series of tubes(when only one network loaded) | 158 | ///render the network cylinder as a series of tubes(when only one network loaded) |
159 | void glCylinder0() { | 159 | void glCylinder0() { |
160 | 160 | ||
161 | + float r1, r2; | ||
161 | if (!glIsList(dlist)) { // if dlist isn't a display list, create it | 162 | if (!glIsList(dlist)) { // if dlist isn't a display list, create it |
162 | dlist = glGenLists(1); // generate a display list | 163 | dlist = glGenLists(1); // generate a display list |
163 | glNewList(dlist, GL_COMPILE); // start a new display list | 164 | glNewList(dlist, GL_COMPILE); // start a new display list |
164 | - T radius = 0.0; | ||
165 | for (unsigned e = 0; e < E.size(); e++) { // for each edge in the network | 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 | for (unsigned p = 1; p < E[e].size(); p++) { // for each point on that edge | 166 | for (unsigned p = 1; p < E[e].size(); p++) { // for each point on that edge |
168 | stim::circle<T> C1 = E[e].circ(p - 1); | 167 | stim::circle<T> C1 = E[e].circ(p - 1); |
169 | stim::circle<T> C2 = E[e].circ(p); | 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 | std::vector< stim::vec3<T> > Cp1 = C1.glpoints(20);// get 20 points on the circle plane | 173 | std::vector< stim::vec3<T> > Cp1 = C1.glpoints(20);// get 20 points on the circle plane |
173 | std::vector< stim::vec3<T> > Cp2 = C2.glpoints(20); | 174 | std::vector< stim::vec3<T> > Cp2 = C2.glpoints(20); |
174 | glBegin(GL_QUAD_STRIP); | 175 | glBegin(GL_QUAD_STRIP); |
@@ -180,8 +181,17 @@ public: | @@ -180,8 +181,17 @@ public: | ||
180 | } // set the texture coordinate based on the specified magnitude index | 181 | } // set the texture coordinate based on the specified magnitude index |
181 | } | 182 | } |
182 | for (unsigned n = 0; n < V.size(); n++) { | 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 | glEndList(); // end the display list | 196 | glEndList(); // end the display list |
187 | } | 197 | } |