codebase / stimlib

Commit f48ef9324d608f99f44a0d6c2063889b638946ad

Authored by Jiaming Guo
1 parent fdfdeda0

make netmets works for new cylinder

Showing 2 changed files with 68 additions and 60 deletions   Show diff stats
Inline Side-by-side
stim/structures/kdtree.cuh
  Show/Hide comments   View file @f48ef93
... ... @@ -347,6 +347,10 @@ namespace stim {
347 347 }
348 348 }
349 349  
  350 + /// create a KD-tree given a pointer to an array of reference points and the number of reference points
  351 + /// @param h_reference_points is a host array containing the reference points in (x0, y0, z0, ...., ) order
  352 + /// @param reference_count is the number of reference point in the array
  353 + /// @param max_levels is the deepest number of tree levels allowed
350 354 void create(T *h_reference_points, size_t reference_count, size_t max_levels) {
351 355 #ifdef __CUDACC__
352 356 if (max_levels > 10) {
... ... @@ -552,6 +556,11 @@ namespace stim {
552 556 }
553 557 }
554 558  
  559 + /// search the KD tree for nearest neighbors to a set of specified query points
  560 + /// @param h_query_points an array of query points in (x0, y0, z0, ...) order
  561 + /// @param query_count is the number of query points
  562 + /// @param indices are the indices to the nearest reference point for each query points
  563 + /// @param distances is an array containing the distance between each query point and the nearest reference point
555 564 void search(T *h_query_points, size_t query_count, size_t *indices, T *distances) {
556 565 #ifdef __CUDACC__
557 566 std::vector < typename cpu_kdtree::point<T, D> > query_points(query_count);
... ...
stim/visualization/cylinder.h
  Show/Hide comments   View file @f48ef93
... ... @@ -13,7 +13,7 @@ class cylinder : public centerline&lt;T&gt; {
13 13 protected:
14 14  
15 15 std::vector< stim::vec3<T> > U; //stores the array of U vectors defining the Frenet frame
16   - std::vector< std::vector<T> > M; //stores a list of magnitudes for each point in the centerline (assuming mags[0] is the radius)
  16 + std::vector< T > R; //stores a list of magnitudes for each point in the centerline (assuming mags[0] is the radius)
17 17  
18 18 using stim::centerline<T>::findIdx;
19 19  
... ... @@ -21,6 +21,7 @@ protected:
21 21 // this function assumes that the centerline has already been set
22 22 void init() {
23 23 U.resize(size()); //allocate space for the frenet frame vectors
  24 + R.resize(size());
24 25  
25 26 stim::circle<T> c; //create a circle
26 27 stim::vec3<T> d0, d1;
... ... @@ -44,22 +45,22 @@ public:
44 45  
45 46 cylinder(centerline& c, T r) : centerline(c) {
46 47 init();
47   - add_mag(r);
  48 + //add_mag(r);
48 49 }
49 50  
50 51 //initialize a cylinder with a list of points and magnitude values
51 52 cylinder(centerline& c, std::vector<T> r) : centerline(c){
52 53 init();
53   - add_mag(r);
  54 + //add_mag(r);
54 55 }
55 56  
56 57 ///Returns magnitude i at the given length into the fiber (based on the pvalue).
57 58 ///Interpolates the position along the line.
58 59 ///@param l: the location of the in the cylinder.
59 60 ///@param idx: integer location of the point closest to l but prior to it.
60   - T m(T l, int idx, size_t i) {
  61 + T r(T l, int idx) {
61 62 T a = (l - L[idx]) / (L[idx + 1] - L[idx]);
62   - T v2 = (M[idx][i] + (M[idx + 1][i] - M[idx][i])*a);
  63 + T v2 = (R[idx] + (R[idx + 1] - R[idx])*a);
63 64  
64 65 return v2;
65 66 }
... ... @@ -67,53 +68,53 @@ public:
67 68 ///Returns the ith magnitude at the given p-value (p value ranges from 0 to 1).
68 69 ///interpolates the radius along the line.
69 70 ///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1).
70   - T m(T pvalue, unsigned i = 0) {
71   - if (pvalue <= 0.0) return M[0][i];
72   - if (pvalue >= 1.0) return M[size() - 1][i];
  71 + T rl(T pvalue) {
  72 + if (pvalue <= 0.0) return R[0];
  73 + if (pvalue >= 1.0) return R[size() - 1];
73 74  
74 75 T l = pvalue*L[L.size() - 1];
75 76 int idx = findIdx(l);
76   - return m(l, idx, i);
  77 + return r(l, idx);
77 78 }
78 79  
79 80 /// Returns the magnitude at the given index
80 81 /// @param i is the index of the desired point
81   - /// @param m is the index of the magnitude value
82   - T m(unsigned i, unsigned m) {
83   - return M[i][m];
  82 + /// @param r is the index of the magnitude value
  83 + T r(unsigned i) {
  84 + return R[i];
84 85 }
85 86  
86 87 ///adds a magnitude to each point in the cylinder
87   - void add_mag(T val = 0) {
  88 + /*void add_mag(V val = 0) {
88 89 if (M.size() == 0) M.resize(size()); //if the magnitude vector isn't initialized, resize it to match the centerline
89 90 for (size_t i = 0; i < size(); i++) //for each point
90   - M[i].push_back(val); //add this value to the magnitude vector at each point
91   - }
  91 + R[i].push_back(val); //add this value to the magnitude vector at each point
  92 + }*/
92 93  
93 94 //adds a magnitude based on a list of magnitudes for each point
94   - void add_mag(std::vector<T> val) {
  95 + /*void add_mag(std::vector<T> val) {
95 96 if (M.size() == 0) M.resize(size()); //if the magnitude vector isn't initialized, resize it to match the centerline
96 97 for (size_t i = 0; i < size(); i++) //for each point
97   - M[i].push_back(val[i]); //add this value to the magnitude vector at each point
98   - }
  98 + R[i].push_back(val[i]); //add this value to the magnitude vector at each point
  99 + }*/
99 100  
100 101 //sets the value of magnitude m at point i
101   - void set_mag(size_t m, size_t i, T v) {
102   - M[i][m] = v;
  102 + void set_r(size_t i, T r) {
  103 + R[i] = r;
103 104 }
104 105  
105   - size_t nmags() {
  106 + /*size_t nmags() {
106 107 if (M.size() == 0) return 0;
107   - else return M[0].size();
108   - }
  108 + else return R[0].size();
  109 + }*/
109 110  
110 111 ///Returns a circle representing the cylinder cross section at point i
111   - stim::circle<T> circ(size_t i, size_t m = 0) {
112   - return stim::circle<T>(at(i), M[i][m], d(i), U[i]);
  112 + stim::circle<T> circ(size_t i) {
  113 + return stim::circle<T>(at(i), R[i], d(i), U[i]);
113 114 }
114 115  
115 116 ///Returns an OBJ object representing the cylinder with a radial tesselation value of N using magnitude m
116   - stim::obj<T> OBJ(size_t N, size_t m = 0) {
  117 + stim::obj<T> OBJ(size_t N) {
117 118 stim::obj<T> out; //create an OBJ object
118 119 stim::circle<T> c0, c1;
119 120 std::vector< stim::vec3<T> > p0, p1; //allocate space for the point sets representing the circles bounding each cylinder segment
... ... @@ -165,7 +166,7 @@ public:
165 166 size_t N = std::vector< stim::vec3<T> >::size();
166 167 ss << "---------[" << N << "]---------" << std::endl;
167 168 for (size_t i = 0; i < N; i++)
168   - ss << std::vector< stim::vec3<T> >::at(i) << " r = " << M[i][0] << " u = " << U[i] << std::endl;
  169 + ss << std::vector< stim::vec3<T> >::at(i) << " r = " << R[i] << " u = " << U[i] << std::endl;
169 170 ss << "--------------------" << std::endl;
170 171  
171 172 return ss.str();
... ... @@ -173,19 +174,19 @@ public:
173 174  
174 175 /// Integrates a magnitude value along the cylinder.
175 176 /// @param m is the magnitude value to be integrated (this is usually the radius)
176   - T integrate(size_t m = 0) {
  177 + T integrate() {
177 178  
178 179 T sum = 0; //initialize the integral to zero
179 180 T m0, m1; //allocate space for both magnitudes in a single segment
180   - m0 = M[0][m]; //initialize the first point and magnitude to the first point in the cylinder
181   - T len = L[1]; //allocate space for the segment length
  181 + m0 = R[0]; //initialize the first point and magnitude to the first point in the cylinder
  182 + T len = L[0]; //allocate space for the segment length
182 183  
183 184  
184 185 for (unsigned p = 1; p < size(); p++) { //for every consecutive point in the cylinder
185   - m1 = M[p][m];
186   - if (p > 1) len = (L[p] - L[p - 1]); //calculate the segment length using the L array
  186 + m1 = R[p];
  187 + if (p > 1) len = (L[p - 1] - L[p - 2]); //calculate the segment length using the L array
187 188 sum += (m0 + m1) / (T)2.0 * len; //add the average magnitude, weighted by the segment length
188   - m0 = m1; //move to the next segment by shifting points
  189 + m0 = m1; //move to the next segment by shifting points
189 190 }
190 191 return sum; //return the integral
191 192 }
... ... @@ -196,27 +197,25 @@ public:
196 197 cylinder<T> resample(T spacing) {
197 198 cylinder<T> c = stim::centerline<T>::resample(spacing); //resample the centerline and use it to create a new cylinder
198 199  
199   - size_t nm = nmags(); //get the number of magnitude values in the current cylinder
200   - if (nm > 0) { //if there are magnitude values
201   - std::vector<T> magvec(nm, 0); //create a magnitude vector for a single point
202   - c.M.resize(c.size(), magvec); //allocate space for a magnitude vector at each point of the new cylinder
203   - }
  200 + //size_t nm = nmags(); //get the number of magnitude values in the current cylinder
  201 + //if (nm > 0) { //if there are magnitude values
  202 + // std::vector<T> magvec(nm, 0); //create a magnitude vector for a single point
  203 + // c.M.resize(c.size(), magvec); //allocate space for a magnitude vector at each point of the new cylinder
  204 + //}
204 205  
205 206 T l, t;
206 207 for (size_t i = 0; i < c.size(); i++) { //for each point in the new cylinder
207 208 l = c.L[i]; //get the length along the new cylinder
208 209 t = l / length(); //calculate the parameter value along the new cylinder
209   - for (size_t mag = 0; mag < nm; mag++) { //for each magnitude value
210   - c.M[i][mag] = m(t, mag); //retrieve the interpolated magnitude from the current cylinder and store it in the new one
211   - }
  210 + //for (size_t mag = 0; mag < nm; mag++) { //for each magnitude value
  211 + c.R[i] = r(t); //retrieve the interpolated magnitude from the current cylinder and store it in the new one
  212 + //}
212 213 }
213   -
214   -
215 214 return c;
216 215 }
217 216  
218   - std::vector< stim::cylinder<T> > split(unsigned int idx){
219   -
  217 + std::vector< stim::cylinder<T> > split(unsigned int idx) {
  218 +
220 219 unsigned N = size();
221 220 std::vector< stim::centerline<T> > LL;
222 221 LL.resize(2);
... ... @@ -224,21 +223,21 @@ public:
224 223 std::vector< stim::cylinder<T> > C(LL.size());
225 224 unsigned i = 0;
226 225 C[0] = LL[0];
227   - C[0].M.resize(idx);
228   - for( ; i < idx; i++){
  226 + //C[0].R.resize(idx);
  227 + for (; i < idx; i++) {
229 228 //for(unsigned d = 0; d < 3; d++)
230   - //C[0][i][d] = LL[0].c[i][d];
231   - C[0].M[i] = M[i];
232   - C[0].M[i].resize(1);
  229 + //C[0][i][d] = LL[0].c[i][d];
  230 + C[0].R[i] = R[i];
  231 + //C[0].R[i].resize(1);
233 232 }
234   - if(C.size() == 2){
  233 + if (C.size() == 2) {
235 234 C[1] = LL[1];
236   - C[1].M.resize(N - idx);
237   - for( ; i < N; i++){
  235 + //C[1].M.resize(N - idx);
  236 + for (; i < N; i++) {
238 237 //for(unsigned d = 0; d < 3; d++)
239   - //C[1][i - idx][d] = LL[1].c[i - idx][d];
240   - C[1].M[i - idx] = M[i];
241   - C[1].M[i - idx].resize(1);
  238 + //C[1][i - idx][d] = LL[1].c[i - idx][d];
  239 + C[1].R[i - idx] = R[i];
  240 + //C[1].M[i - idx].resize(1);
242 241 }
243 242 }
244 243  
... ... @@ -273,7 +272,7 @@ public:
273 272 }
274 273  
275 274 stim::vec3<T> dr = (inP[1] - inP[0]).norm();
276   - s = stim::circle<T>(inP[0], inM[0][0], dr, stim::vec3<T>(1,0,0));
  275 + s = stim::circle<T>(inP[0], inR[0][0], dr, stim::vec3<T>(1,0,0));
277 276 norms[0] = s.N;
278 277 e[0] = s;
279 278 Us[0] = e[0].U;
... ... @@ -287,7 +286,7 @@ public:
287 286  
288 287 norms[i] = s.N;
289 288  
290   - s.scale(inM[i][0]/inM[i-1][0]);
  289 + s.scale(inR[i][0]/inR[i-1][0]);
291 290 e[i] = s;
292 291 Us[i] = e[i].U;
293 292 }
... ... @@ -299,7 +298,7 @@ public:
299 298  
300 299 norms[j] = s.N;
301 300  
302   - s.scale(inM[j][0]/inM[j-1][0]);
  301 + s.scale(inR[j][0]/inR[j-1][0]);
303 302 e[j] = s;
304 303 Us[j] = e[j].U;
305 304 }
... ... @@ -426,7 +425,7 @@ public:
426 425 stim::vec<T> zero(0.0,0.0);
427 426 temp.resize(inM.size(), zero);
428 427 for(int i = 0; i < inM.size(); i++)
429   - temp[i][0] = inM[i];
  428 + temp[i][0] = inR[i];
430 429 init(inP, temp);
431 430 }
432 431  
... ...