diff --git a/stim/cuda/branch_detection.cuh b/stim/cuda/branch_detection.cuh index ec7d6e9..abd9b44 100644 --- a/stim/cuda/branch_detection.cuh +++ b/stim/cuda/branch_detection.cuh @@ -152,16 +152,17 @@ find_branch(GLint texbufferID, GLenum texType, unsigned int x, unsigned int y) stim::gpu2image(gpuCenters, "Centers.jpg", x, y, 0, 1); cudaMemcpy(cpuCenters, gpuCenters, bytes_ds, cudaMemcpyDeviceToHost); stim::cpu2image(cpuCenters, "CentersXPU.jpg", x, y, 0, 1); + std::cout << pixels << " " << x << " " << y << std::endl; for(int i = 0; i < pixels; i++) { int ix = (i % x); int iy = (i / x); - if((cpuCenters[i] != 0)) + if((cpuCenters[i] == 1)) { float x_v = (float) ix; float y_v = (float) iy; - output.push_back(stim::vec((x_v/(float)x), + output.push_back(stim::vec((x_v/(float)x*360), (y_v/(float)y), 0.0)); } diff --git a/stim/gl/gl_spider.h b/stim/gl/gl_spider.h index ae54d45..09f5975 100644 --- a/stim/gl/gl_spider.h +++ b/stim/gl/gl_spider.h @@ -229,21 +229,21 @@ class gl_spider : public virtual gl_texture DrawLongCylinder(n, l_template, l_square); stim::cylinder cyl(cL, cM); std::vector< stim::vec > result = find_branch(btexbufferID, GL_TEXTURE_2D, n*l_square, (cL.size()-1)*l_template); - std::cerr << "the number of points is " << result.size() << std::endl; +// std::cerr << "the number of points is " << result.size() << std::endl; if(!result.empty()) { for(int i = 0; i < result.size(); i++) { - std::cout << "Testing "<< i << ": " << result[i][0] << ", " << result[i][1]*360.0 << std::endl; - stim::vec v = cyl.surf(result[i][0], result[i][1]*360.0); - std::cout << v[0] << " ," << v[1] << " ," << v[2] << std::endl; - stim::vec di = cyl.p(result[i][0]); + std::cout << "Testing "<< i << ": " << result[i][0] << ", " << result[i][1] << std::endl; + stim::vec v = cyl.surf(result[i][1], result[i][0]); +// std::cout << v[0] << " ," << v[1] << " ," << v[2] << std::endl; + stim::vec di = cyl.p(result[i][1]); std::cout << di[0] << " ," << di[1] << " ," << di[2] << std::endl; - float rad = cyl.r(result[i][0]); - std::cout << rad << std::endl; + float rad = cyl.r(result[i][1]); +// std::cout << rad << std::endl; setSeed(v); setSeedVec(stim::vec(-di[0]+v[0], -di[1]+v[1], -di[2]+v[2]).norm()); - setSeedMag(cyl.r(result[i][0])); + setSeedMag(cyl.r(result[i][1])); } } std::cout << "I ran the new branch detection" << std::endl; diff --git a/stim/math/circle.h b/stim/math/circle.h new file mode 100644 index 0000000..f7acde9 --- /dev/null +++ b/stim/math/circle.h @@ -0,0 +1,179 @@ +#ifndef STIM_CIRCLE_H +#define STIM_CIRCLE_H + + +#include +#include +#include +#include +#include +#include +#include + +namespace stim{ + +template +class circle : plane +{ + +private: + + stim::vec Y; + + CUDA_CALLABLE void + init() + { + Y = U.cross(N).norm(); + } + +public: + using stim::plane::n; + using stim::plane::P; + using stim::plane::N; + using stim::plane::U; + using stim::plane::rotate; + using stim::plane::setU; + + ///base constructor + ///@param th value of the angle of the starting point from 0 to 360. + CUDA_CALLABLE + circle() : plane() + { + init(); + } + + ///create a rectangle given a size and position in Z space. + ///@param size: size of the rectangle in ND space. + ///@param z_pos z coordinate of the rectangle. + CUDA_CALLABLE + circle(T size, T z_pos = (T)0) : plane() + { + init(); + center(stim::vec(0,0,z_pos)); + scale(size); + } + + ///create a rectangle from a center point, normal + ///@param c: x,y,z location of the center. + ///@param n: x,y,z direction of the normal. + CUDA_CALLABLE + circle(vec c, vec n = vec(0,0,1)) : plane() + { + center(c); + normal(n); + init(); + } + + ///create a rectangle from a center point, normal, and size + ///@param c: x,y,z location of the center. + ///@param s: size of the rectangle. + ///@param n: x,y,z direction of the normal. + CUDA_CALLABLE + circle(vec c, T s, vec n = vec(0,0,1)) : plane() + { + init(); + center(c); + rotate(n, U, Y); + scale(s); + } + + ///create a rectangle from a center point, normal, and size + ///@param c: x,y,z location of the center. + ///@param s: size of the rectangle. + ///@param n: x,y,z direction of the normal. + ///@param u: x,y,z direction for the zero vector (from where the rotation starts) + CUDA_CALLABLE + circle(vec c, T s, vec n = vec(0,0,1), vec u = vec(1, 0, 0)) : plane() + { + init(); + setU(u); + center(c); + normal(n); + scale(s); + } + + ///scales the circle by a certain factor + ///@param factor: the factor by which the dimensions of the shape are scaled. + CUDA_CALLABLE + void scale(T factor) + { + U *= factor; + Y *= factor; + } + + ///sets the normal for the cirlce + ///@param n: x,y,z direction of the normal. + CUDA_CALLABLE void + normal(vec n) + { + rotate(n, Y); + } + + ///sets the center of the circle. + ///@param n: x,y,z location of the center. + CUDA_CALLABLE T + center(vec p) + { + this->P = p; + } + + ///boolean comparison + bool + operator==(const circle & rhs) + { + if(P == rhs.P && U == rhs.U && Y == rhs.Y) + return true; + else + return false; + } + + ///get the world space value given the planar coordinates a, b in [0, 1] + CUDA_CALLABLE stim::vec p(T a, T b) + { + stim::vec result; + + vec A = this->P - this->U * (T)0.5 - Y * (T)0.5; + result = A + this->U * a + Y * b; + return result; + } + + ///parenthesis operator returns the world space given rectangular coordinates a and b in [0 1] + CUDA_CALLABLE stim::vec operator()(T a, T b) + { + return p(a,b); + } + + ///returns a vector with the points on the initialized circle. + ///connecting the points results in a circle. + ///@param n: integer for the number of points representing the circle. + std::vector > + getPoints(int n) + { + std::vector > result; + stim::vec point; + T x,y; + float step = 360.0/(float) n; + for(float j = 0; j <= 360.0; j += step) + { + y = 0.5*cos(j*2.0*M_PI/360.0)+0.5; + x = 0.5*sin(j*2.0*M_PI/360.0)+0.5; + result.push_back(p(x,y)); + } + return result; + } + + ///returns a vector with the points on the initialized circle. + ///connecting the points results in a circle. + ///@param n: integer for the number of points representing the circle. + stim::vec + p(T theta) + { + T x,y; + y = 0.5*cos(theta*2.0*M_PI/360.0)+0.5; + x = 0.5*sin(theta*2.0*M_PI/360.0)+0.5; + return p(x,y); + } + +}; +} +#endif diff --git a/stim/math/plane.h b/stim/math/plane.h index baa4683..570830a 100644 --- a/stim/math/plane.h +++ b/stim/math/plane.h @@ -194,6 +194,17 @@ class plane } + CUDA_CALLABLE void rotate(vec n, vec &Y) + { + quaternion q; + q.CreateRotation(N, n); + + N = q.toMatrix3() * N; + U = q.toMatrix3() * U; + Y = q.toMatrix3() * Y; + + } + CUDA_CALLABLE void rotate(vec n, vec &X, vec &Y) { quaternion q; diff --git a/stim/visualization/cylinder.h b/stim/visualization/cylinder.h new file mode 100644 index 0000000..c2f7e11 --- /dev/null +++ b/stim/visualization/cylinder.h @@ -0,0 +1,245 @@ +#ifndef STIM_CYLINDER_H +#define STIM_CYLINDER_H +#include +#include +#include + + +namespace stim +{ +template +class cylinder +{ + private: + stim::circle s; //an arbitrary circle + std::vector< stim::vec > pos; //positions of the cylinder. + std::vector< stim::vec > mags; //radii at each position + std::vector< T > L; //length of the cylinder at each position. + + ///default init + void + init() + { + + } + + ///inits the cylinder from a list of points (inP) and radii (inM) + void + init(std::vector > inP, std::vector > inM) + { + pos = inP; + mags = inM; + + //calculate each L. + L.resize(pos.size()); + T temp = (T)0; + L[0] = 0; + for(int i = 1; i < L.size(); i++) + { + temp += (pos[i-1] - pos[i]).len(); + L[i] = temp; + } + } + + ///returns the direction vector at point idx. + stim::vec + d(int idx) + { + if(idx == 0) + { + return (pos[idx+1] - pos[idx]).norm(); + } + else if(idx == pos.size()-1) + { + return (pos[idx] - pos[idx-1]).norm(); + } + else + { + stim::vec v1 = (pos[idx]-pos[idx-1]).norm(); + stim::vec v2 = (pos[idx+1]-pos[idx]).norm(); + return (v1+v2).norm(); + } + + } + + ///returns the total length of the line at index j. + T + getl(int j) + { + T temp = (T) 0; + for(int i = 0; i < j; ++i) + { + temp += (pos[i] - pos[i+1]).len(); + L[i] = temp; + } + } + + ///finds the index of the point closest to the length l on the lower bound. + ///binary search. + int + findIdx(T l) + { + unsigned int i = L.size()/2; + unsigned int max = L.size()-1; + unsigned int min = 0; +// std::cerr << "Index initially: " << i << std::endl; +// std::cerr << "l initially: " << l << std::endl; + while(i > 0 && i < L.size()-1) + { +// std::cerr << "L[i] = " << L[i] << " i= " << i << " maxval= " << L.size()-1 << std::endl; + if(l < L[i]) + { + max = i; +// std::cerr << l << " < " << L[i] << "so max is set to " << i << " and the new i is " << min+(max-min)/2 << std::endl; + i = min+(max-min)/2; + } + else if(L[i] <= l && L[i+1] >= l) + { +// std::cerr << L[i] << " < " << l << " < " << L[i+1] << std::endl; + break; + } + else + { + min = i; +// std::cerr << l << " > " << L[i] << "so min is set to " << i << " and the new i is " << min+(max-min)/2 << std::endl; + i = min+(max-min)/2; + } + } + return i; + } + + public: + ///default constructor + cylinder() + { + + } + + ///constructor to create a cylinder from a set of points, radii, and the number of sides for the cylinder. + ///@param inP: Vector of stim vecs composing the points of the centerline. + ///@param inM: Vector of stim vecs composing the radii of the centerline. + cylinder(std::vector > inP, std::vector > inM) + { + init(inP, inM); + } + + + ///Returns a position vector at the given p-value (p value ranges from 0 to 1). + ///interpolates the position along the line. + ///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1). + stim::vec + p(T pvalue) + { + if(pvalue < 0.0 || pvalue > 1.0) + { + return stim::vec(-1,-1,-1); + } + T l = pvalue*L[L.size()-1]; + int idx = findIdx(l); + T rat = (l-L[idx])/(L[idx+1]-L[idx]); + return( pos[idx] + (pos[idx+1]-pos[idx])*rat); + } + + ///Returns a position vector at the given length into the fiber (based on the pvalue). + ///Interpolates the radius along the line. + ///@param l: the location of the in the cylinder. + ///@param idx: integer location of the point closest to l but prior to it. + stim::vec + p(T l, int idx) + { + T rat = (l-L[idx])/(L[idx+1]-L[idx]); + return( pos[idx] + (pos[idx+1]-pos[idx])*rat); +// return( +// return (pos[idx] + (pos[idx+1]-pos[idx])*((l-L[idx])/(L[idx+1]- L[idx]))); + } + + ///Returns a radius at the given p-value (p value ranges from 0 to 1). + ///interpolates the radius along the line. + ///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1). + T + r(T pvalue) + { + if(pvalue < 0.0 || pvalue > 1.0) + return; + T l = pvalue*L[L.size()-1]; + int idx = findIdx(l); + return (mags[idx] + (mags[idx+1]-mags[idx])*((l-L[idx])/(L[idx+1]- L[idx])))[0]; + } + + ///Returns a radius at the given length into the fiber (based on the pvalue). + ///Interpolates the position along the line. + ///@param l: the location of the in the cylinder. + ///@param idx: integer location of the point closest to l but prior to it. + T + r(T l, int idx) + { + T rat = (l-L[idx])/(L[idx+1]-L[idx]); + return( mags[idx] + (mags[idx+1]-mags[idx])*rat)[0]; +// return (mags[idx] + (mags[idx+1]-mags[idx])*((l-L[idx])/(L[idx+1]- L[idx])))[0]; + } + + + ///returns the position of the point with a given pvalue and theta on the surface + ///in x, y, z coordinates. Theta is in degrees from 0 to 360. + ///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1). + ///@param theta: the angle to the point of a circle. + stim::vec + surf(T pvalue, T theta) + { + if(pvalue < 0.0 || pvalue > 1.0) + { + return stim::vec(-1,-1,-1); + } else { +// std::cout << "AAAAAAAAAAAAAAAAAAAAAAA " << pvalue*L[L.size()-1] << " " << L[L.size()-1] << std::endl; + T l = pvalue*L[L.size()-1]; +// std::cerr << "The l value is: " << l << std::endl; + int idx = findIdx(l); +// std::cerr << "Index: " << idx << std::endl; + stim::vec ps = p(l, idx); + T m = r(l, idx); + stim::vec dr = d(idx); + s = stim::circle(ps, m, dr, stim::vec(1,0,0)); + return(s.p(theta)); + } + } + + ///returns a vector of points necessary to create a circle at every position in the fiber. + ///@param sides: the number of sides of each circle. + std::vector > > + getPoints(int sides) + { + if(pos.size() < 2) + { + return; + } else { + std::vector > > points; + points.resize(pos.size()); + stim::vec dr = (pos[1] - pos[0]).norm(); + s = stim::circle(pos[0], mags[0][0], dr, stim::vec(1,0,0)); + points[0] = s.getPoints(sides); + for(int i = 1; i < pos.size(); i++) + { + dr = d(i); +// dr = (pos[i-1] - pos[i]).norm(); + +// s = stim::circle(pos[i], mags[i][0], dr, stim::vec(1,0,0)); + s.center(pos[i]); + s.normal(dr); +// s.scale(mags[i][0]/mags[i-1][0], mags[i][0]/mags[i-1][0]); + s.scale(mags[i][0]/mags[i-1][0]); + points[i] = s.getPoints(sides); + } + return points; + } + } + + void + print(int idx) + { + std::cout << d(idx) << std::endl; + } + +}; + +} +#endif -- libgit2 0.21.4