codebase / stimlib

Commit 26aee9ed47c1c0452be707a77bb59e0ed7c62e3f

Authored by Pavel Govyadinov
1 parent d10af691

changed circle class to extend plane.h instead of extending rect.h. Other minor … 

…stability changes. STABLE- PRECYLINDER WORK
Showing 5 changed files with 446 additions and 10 deletions   Show diff stats
Inline Side-by-side
stim/cuda/branch_detection.cuh
  Show/Hide comments   View file @26aee9e
... ... @@ -152,16 +152,17 @@ find_branch(GLint texbufferID, GLenum texType, unsigned int x, unsigned int y)
152 152 stim::gpu2image<float>(gpuCenters, "Centers.jpg", x, y, 0, 1);
153 153 cudaMemcpy(cpuCenters, gpuCenters, bytes_ds, cudaMemcpyDeviceToHost);
154 154 stim::cpu2image<float>(cpuCenters, "CentersXPU.jpg", x, y, 0, 1);
  155 + std::cout << pixels << " " << x << " " << y << std::endl;
155 156 for(int i = 0; i < pixels; i++)
156 157 {
157 158 int ix = (i % x);
158 159 int iy = (i / x);
159   - if((cpuCenters[i] != 0))
  160 + if((cpuCenters[i] == 1))
160 161 {
161 162  
162 163 float x_v = (float) ix;
163 164 float y_v = (float) iy;
164   - output.push_back(stim::vec<float>((x_v/(float)x),
  165 + output.push_back(stim::vec<float>((x_v/(float)x*360),
165 166 (y_v/(float)y), 0.0));
166 167  
167 168 }
... ...
stim/gl/gl_spider.h
  Show/Hide comments   View file @26aee9e
... ... @@ -229,21 +229,21 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
229 229 DrawLongCylinder(n, l_template, l_square);
230 230 stim::cylinder<float> cyl(cL, cM);
231 231 std::vector< stim::vec<float> > result = find_branch(btexbufferID, GL_TEXTURE_2D, n*l_square, (cL.size()-1)*l_template);
232   - std::cerr << "the number of points is " << result.size() << std::endl;
  232 +// std::cerr << "the number of points is " << result.size() << std::endl;
233 233 if(!result.empty())
234 234 {
235 235 for(int i = 0; i < result.size(); i++)
236 236 {
237   - std::cout << "Testing "<< i << ": " << result[i][0] << ", " << result[i][1]*360.0 << std::endl;
238   - stim::vec<float> v = cyl.surf(result[i][0], result[i][1]*360.0);
239   - std::cout << v[0] << " ," << v[1] << " ," << v[2] << std::endl;
240   - stim::vec<float> di = cyl.p(result[i][0]);
  237 + std::cout << "Testing "<< i << ": " << result[i][0] << ", " << result[i][1] << std::endl;
  238 + stim::vec<float> v = cyl.surf(result[i][1], result[i][0]);
  239 +// std::cout << v[0] << " ," << v[1] << " ," << v[2] << std::endl;
  240 + stim::vec<float> di = cyl.p(result[i][1]);
241 241 std::cout << di[0] << " ," << di[1] << " ," << di[2] << std::endl;
242   - float rad = cyl.r(result[i][0]);
243   - std::cout << rad << std::endl;
  242 + float rad = cyl.r(result[i][1]);
  243 +// std::cout << rad << std::endl;
244 244 setSeed(v);
245 245 setSeedVec(stim::vec<float>(-di[0]+v[0], -di[1]+v[1], -di[2]+v[2]).norm());
246   - setSeedMag(cyl.r(result[i][0]));
  246 + setSeedMag(cyl.r(result[i][1]));
247 247 }
248 248 }
249 249 std::cout << "I ran the new branch detection" << std::endl;
... ...
stim/math/circle.h 0 → 100644
  Show/Hide comments   View file @26aee9e
  1 +#ifndef STIM_CIRCLE_H
  2 +#define STIM_CIRCLE_H
  3 +
  4 +
  5 +#include <stim/cuda/cudatools/callable.h>
  6 +#include <stim/math/plane.h>
  7 +#include <stim/math/vector.h>
  8 +#include <stim/math/triangle.h>
  9 +#include <assert.h>
  10 +#include <algorithm>
  11 +#include <iostream>
  12 +
  13 +namespace stim{
  14 +
  15 +template <typename T>
  16 +class circle : plane<T>
  17 +{
  18 +
  19 +private:
  20 +
  21 + stim::vec<T> Y;
  22 +
  23 + CUDA_CALLABLE void
  24 + init()
  25 + {
  26 + Y = U.cross(N).norm();
  27 + }
  28 +
  29 +public:
  30 + using stim::plane<T>::n;
  31 + using stim::plane<T>::P;
  32 + using stim::plane<T>::N;
  33 + using stim::plane<T>::U;
  34 + using stim::plane<T>::rotate;
  35 + using stim::plane<T>::setU;
  36 +
  37 + ///base constructor
  38 + ///@param th value of the angle of the starting point from 0 to 360.
  39 + CUDA_CALLABLE
  40 + circle() : plane<T>()
  41 + {
  42 + init();
  43 + }
  44 +
  45 + ///create a rectangle given a size and position in Z space.
  46 + ///@param size: size of the rectangle in ND space.
  47 + ///@param z_pos z coordinate of the rectangle.
  48 + CUDA_CALLABLE
  49 + circle(T size, T z_pos = (T)0) : plane<T>()
  50 + {
  51 + init();
  52 + center(stim::vec<T>(0,0,z_pos));
  53 + scale(size);
  54 + }
  55 +
  56 + ///create a rectangle from a center point, normal
  57 + ///@param c: x,y,z location of the center.
  58 + ///@param n: x,y,z direction of the normal.
  59 + CUDA_CALLABLE
  60 + circle(vec<T> c, vec<T> n = vec<T>(0,0,1)) : plane<T>()
  61 + {
  62 + center(c);
  63 + normal(n);
  64 + init();
  65 + }
  66 +
  67 + ///create a rectangle from a center point, normal, and size
  68 + ///@param c: x,y,z location of the center.
  69 + ///@param s: size of the rectangle.
  70 + ///@param n: x,y,z direction of the normal.
  71 + CUDA_CALLABLE
  72 + circle(vec<T> c, T s, vec<T> n = vec<T>(0,0,1)) : plane<T>()
  73 + {
  74 + init();
  75 + center(c);
  76 + rotate(n, U, Y);
  77 + scale(s);
  78 + }
  79 +
  80 + ///create a rectangle from a center point, normal, and size
  81 + ///@param c: x,y,z location of the center.
  82 + ///@param s: size of the rectangle.
  83 + ///@param n: x,y,z direction of the normal.
  84 + ///@param u: x,y,z direction for the zero vector (from where the rotation starts)
  85 + CUDA_CALLABLE
  86 + circle(vec<T> c, T s, vec<T> n = vec<T>(0,0,1), vec<T> u = vec<T>(1, 0, 0)) : plane<T>()
  87 + {
  88 + init();
  89 + setU(u);
  90 + center(c);
  91 + normal(n);
  92 + scale(s);
  93 + }
  94 +
  95 + ///scales the circle by a certain factor
  96 + ///@param factor: the factor by which the dimensions of the shape are scaled.
  97 + CUDA_CALLABLE
  98 + void scale(T factor)
  99 + {
  100 + U *= factor;
  101 + Y *= factor;
  102 + }
  103 +
  104 + ///sets the normal for the cirlce
  105 + ///@param n: x,y,z direction of the normal.
  106 + CUDA_CALLABLE void
  107 + normal(vec<T> n)
  108 + {
  109 + rotate(n, Y);
  110 + }
  111 +
  112 + ///sets the center of the circle.
  113 + ///@param n: x,y,z location of the center.
  114 + CUDA_CALLABLE T
  115 + center(vec<T> p)
  116 + {
  117 + this->P = p;
  118 + }
  119 +
  120 + ///boolean comparison
  121 + bool
  122 + operator==(const circle<T> & rhs)
  123 + {
  124 + if(P == rhs.P && U == rhs.U && Y == rhs.Y)
  125 + return true;
  126 + else
  127 + return false;
  128 + }
  129 +
  130 + ///get the world space value given the planar coordinates a, b in [0, 1]
  131 + CUDA_CALLABLE stim::vec<T> p(T a, T b)
  132 + {
  133 + stim::vec<T> result;
  134 +
  135 + vec<T> A = this->P - this->U * (T)0.5 - Y * (T)0.5;
  136 + result = A + this->U * a + Y * b;
  137 + return result;
  138 + }
  139 +
  140 + ///parenthesis operator returns the world space given rectangular coordinates a and b in [0 1]
  141 + CUDA_CALLABLE stim::vec<T> operator()(T a, T b)
  142 + {
  143 + return p(a,b);
  144 + }
  145 +
  146 + ///returns a vector with the points on the initialized circle.
  147 + ///connecting the points results in a circle.
  148 + ///@param n: integer for the number of points representing the circle.
  149 + std::vector<stim::vec<T> >
  150 + getPoints(int n)
  151 + {
  152 + std::vector<stim::vec<T> > result;
  153 + stim::vec<T> point;
  154 + T x,y;
  155 + float step = 360.0/(float) n;
  156 + for(float j = 0; j <= 360.0; j += step)
  157 + {
  158 + y = 0.5*cos(j*2.0*M_PI/360.0)+0.5;
  159 + x = 0.5*sin(j*2.0*M_PI/360.0)+0.5;
  160 + result.push_back(p(x,y));
  161 + }
  162 + return result;
  163 + }
  164 +
  165 + ///returns a vector with the points on the initialized circle.
  166 + ///connecting the points results in a circle.
  167 + ///@param n: integer for the number of points representing the circle.
  168 + stim::vec<T>
  169 + p(T theta)
  170 + {
  171 + T x,y;
  172 + y = 0.5*cos(theta*2.0*M_PI/360.0)+0.5;
  173 + x = 0.5*sin(theta*2.0*M_PI/360.0)+0.5;
  174 + return p(x,y);
  175 + }
  176 +
  177 +};
  178 +}
  179 +#endif
... ...
stim/math/plane.h
  Show/Hide comments   View file @26aee9e
... ... @@ -194,6 +194,17 @@ class plane
194 194  
195 195 }
196 196  
  197 + CUDA_CALLABLE void rotate(vec<T> n, vec<T> &Y)
  198 + {
  199 + quaternion<T> q;
  200 + q.CreateRotation(N, n);
  201 +
  202 + N = q.toMatrix3() * N;
  203 + U = q.toMatrix3() * U;
  204 + Y = q.toMatrix3() * Y;
  205 +
  206 + }
  207 +
197 208 CUDA_CALLABLE void rotate(vec<T> n, vec<T> &X, vec<T> &Y)
198 209 {
199 210 quaternion<T> q;
... ...
stim/visualization/cylinder.h 0 → 100644
  Show/Hide comments   View file @26aee9e
  1 +#ifndef STIM_CYLINDER_H
  2 +#define STIM_CYLINDER_H
  3 +#include <iostream>
  4 +#include <stim/math/circle.h>
  5 +#include <stim/math/vector.h>
  6 +
  7 +
  8 +namespace stim
  9 +{
  10 +template<typename T>
  11 +class cylinder
  12 +{
  13 + private:
  14 + stim::circle<T> s; //an arbitrary circle
  15 + std::vector< stim::vec<T> > pos; //positions of the cylinder.
  16 + std::vector< stim::vec<T> > mags; //radii at each position
  17 + std::vector< T > L; //length of the cylinder at each position.
  18 +
  19 + ///default init
  20 + void
  21 + init()
  22 + {
  23 +
  24 + }
  25 +
  26 + ///inits the cylinder from a list of points (inP) and radii (inM)
  27 + void
  28 + init(std::vector<stim::vec<T> > inP, std::vector<stim::vec<T> > inM)
  29 + {
  30 + pos = inP;
  31 + mags = inM;
  32 +
  33 + //calculate each L.
  34 + L.resize(pos.size());
  35 + T temp = (T)0;
  36 + L[0] = 0;
  37 + for(int i = 1; i < L.size(); i++)
  38 + {
  39 + temp += (pos[i-1] - pos[i]).len();
  40 + L[i] = temp;
  41 + }
  42 + }
  43 +
  44 + ///returns the direction vector at point idx.
  45 + stim::vec<T>
  46 + d(int idx)
  47 + {
  48 + if(idx == 0)
  49 + {
  50 + return (pos[idx+1] - pos[idx]).norm();
  51 + }
  52 + else if(idx == pos.size()-1)
  53 + {
  54 + return (pos[idx] - pos[idx-1]).norm();
  55 + }
  56 + else
  57 + {
  58 + stim::vec<float> v1 = (pos[idx]-pos[idx-1]).norm();
  59 + stim::vec<float> v2 = (pos[idx+1]-pos[idx]).norm();
  60 + return (v1+v2).norm();
  61 + }
  62 +
  63 + }
  64 +
  65 + ///returns the total length of the line at index j.
  66 + T
  67 + getl(int j)
  68 + {
  69 + T temp = (T) 0;
  70 + for(int i = 0; i < j; ++i)
  71 + {
  72 + temp += (pos[i] - pos[i+1]).len();
  73 + L[i] = temp;
  74 + }
  75 + }
  76 +
  77 + ///finds the index of the point closest to the length l on the lower bound.
  78 + ///binary search.
  79 + int
  80 + findIdx(T l)
  81 + {
  82 + unsigned int i = L.size()/2;
  83 + unsigned int max = L.size()-1;
  84 + unsigned int min = 0;
  85 +// std::cerr << "Index initially: " << i << std::endl;
  86 +// std::cerr << "l initially: " << l << std::endl;
  87 + while(i > 0 && i < L.size()-1)
  88 + {
  89 +// std::cerr << "L[i] = " << L[i] << " i= " << i << " maxval= " << L.size()-1 << std::endl;
  90 + if(l < L[i])
  91 + {
  92 + max = i;
  93 +// std::cerr << l << " < " << L[i] << "so max is set to " << i << " and the new i is " << min+(max-min)/2 << std::endl;
  94 + i = min+(max-min)/2;
  95 + }
  96 + else if(L[i] <= l && L[i+1] >= l)
  97 + {
  98 +// std::cerr << L[i] << " < " << l << " < " << L[i+1] << std::endl;
  99 + break;
  100 + }
  101 + else
  102 + {
  103 + min = i;
  104 +// std::cerr << l << " > " << L[i] << "so min is set to " << i << " and the new i is " << min+(max-min)/2 << std::endl;
  105 + i = min+(max-min)/2;
  106 + }
  107 + }
  108 + return i;
  109 + }
  110 +
  111 + public:
  112 + ///default constructor
  113 + cylinder()
  114 + {
  115 +
  116 + }
  117 +
  118 + ///constructor to create a cylinder from a set of points, radii, and the number of sides for the cylinder.
  119 + ///@param inP: Vector of stim vecs composing the points of the centerline.
  120 + ///@param inM: Vector of stim vecs composing the radii of the centerline.
  121 + cylinder(std::vector<stim::vec<T> > inP, std::vector<stim::vec<T> > inM)
  122 + {
  123 + init(inP, inM);
  124 + }
  125 +
  126 +
  127 + ///Returns a position vector at the given p-value (p value ranges from 0 to 1).
  128 + ///interpolates the position along the line.
  129 + ///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1).
  130 + stim::vec<T>
  131 + p(T pvalue)
  132 + {
  133 + if(pvalue < 0.0 || pvalue > 1.0)
  134 + {
  135 + return stim::vec<float>(-1,-1,-1);
  136 + }
  137 + T l = pvalue*L[L.size()-1];
  138 + int idx = findIdx(l);
  139 + T rat = (l-L[idx])/(L[idx+1]-L[idx]);
  140 + return( pos[idx] + (pos[idx+1]-pos[idx])*rat);
  141 + }
  142 +
  143 + ///Returns a position vector at the given length into the fiber (based on the pvalue).
  144 + ///Interpolates the radius along the line.
  145 + ///@param l: the location of the in the cylinder.
  146 + ///@param idx: integer location of the point closest to l but prior to it.
  147 + stim::vec<T>
  148 + p(T l, int idx)
  149 + {
  150 + T rat = (l-L[idx])/(L[idx+1]-L[idx]);
  151 + return( pos[idx] + (pos[idx+1]-pos[idx])*rat);
  152 +// return(
  153 +// return (pos[idx] + (pos[idx+1]-pos[idx])*((l-L[idx])/(L[idx+1]- L[idx])));
  154 + }
  155 +
  156 + ///Returns a radius at the given p-value (p value ranges from 0 to 1).
  157 + ///interpolates the radius along the line.
  158 + ///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1).
  159 + T
  160 + r(T pvalue)
  161 + {
  162 + if(pvalue < 0.0 || pvalue > 1.0)
  163 + return;
  164 + T l = pvalue*L[L.size()-1];
  165 + int idx = findIdx(l);
  166 + return (mags[idx] + (mags[idx+1]-mags[idx])*((l-L[idx])/(L[idx+1]- L[idx])))[0];
  167 + }
  168 +
  169 + ///Returns a radius at the given length into the fiber (based on the pvalue).
  170 + ///Interpolates the position along the line.
  171 + ///@param l: the location of the in the cylinder.
  172 + ///@param idx: integer location of the point closest to l but prior to it.
  173 + T
  174 + r(T l, int idx)
  175 + {
  176 + T rat = (l-L[idx])/(L[idx+1]-L[idx]);
  177 + return( mags[idx] + (mags[idx+1]-mags[idx])*rat)[0];
  178 +// return (mags[idx] + (mags[idx+1]-mags[idx])*((l-L[idx])/(L[idx+1]- L[idx])))[0];
  179 + }
  180 +
  181 +
  182 + ///returns the position of the point with a given pvalue and theta on the surface
  183 + ///in x, y, z coordinates. Theta is in degrees from 0 to 360.
  184 + ///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1).
  185 + ///@param theta: the angle to the point of a circle.
  186 + stim::vec<T>
  187 + surf(T pvalue, T theta)
  188 + {
  189 + if(pvalue < 0.0 || pvalue > 1.0)
  190 + {
  191 + return stim::vec<float>(-1,-1,-1);
  192 + } else {
  193 +// std::cout << "AAAAAAAAAAAAAAAAAAAAAAA " << pvalue*L[L.size()-1] << " " << L[L.size()-1] << std::endl;
  194 + T l = pvalue*L[L.size()-1];
  195 +// std::cerr << "The l value is: " << l << std::endl;
  196 + int idx = findIdx(l);
  197 +// std::cerr << "Index: " << idx << std::endl;
  198 + stim::vec<T> ps = p(l, idx);
  199 + T m = r(l, idx);
  200 + stim::vec<T> dr = d(idx);
  201 + s = stim::circle<T>(ps, m, dr, stim::vec<T>(1,0,0));
  202 + return(s.p(theta));
  203 + }
  204 + }
  205 +
  206 + ///returns a vector of points necessary to create a circle at every position in the fiber.
  207 + ///@param sides: the number of sides of each circle.
  208 + std::vector<std::vector<vec<T> > >
  209 + getPoints(int sides)
  210 + {
  211 + if(pos.size() < 2)
  212 + {
  213 + return;
  214 + } else {
  215 + std::vector<std::vector <vec<T> > > points;
  216 + points.resize(pos.size());
  217 + stim::vec<T> dr = (pos[1] - pos[0]).norm();
  218 + s = stim::circle<T>(pos[0], mags[0][0], dr, stim::vec<T>(1,0,0));
  219 + points[0] = s.getPoints(sides);
  220 + for(int i = 1; i < pos.size(); i++)
  221 + {
  222 + dr = d(i);
  223 +// dr = (pos[i-1] - pos[i]).norm();
  224 +
  225 +// s = stim::circle<T>(pos[i], mags[i][0], dr, stim::vec<T>(1,0,0));
  226 + s.center(pos[i]);
  227 + s.normal(dr);
  228 +// s.scale(mags[i][0]/mags[i-1][0], mags[i][0]/mags[i-1][0]);
  229 + s.scale(mags[i][0]/mags[i-1][0]);
  230 + points[i] = s.getPoints(sides);
  231 + }
  232 + return points;
  233 + }
  234 + }
  235 +
  236 + void
  237 + print(int idx)
  238 + {
  239 + std::cout << d(idx) << std::endl;
  240 + }
  241 +
  242 +};
  243 +
  244 +}
  245 +#endif
... ...