diff --git a/stim/gl/gl_spider.h b/stim/gl/gl_spider.h index 4217c71..8d970af 100644 --- a/stim/gl/gl_spider.h +++ b/stim/gl/gl_spider.h @@ -81,8 +81,8 @@ class gl_spider : public virtual gl_texture int numSamplesPos; int numSamplesMag; -// float stepsize = 4.0; //Step size. - float stepsize = 3.0; //Step size. + float stepsize = 5.0; //Step size. +// float stepsize = 3.0; //Step size. int current_cost; //variable to store the cost of the current step. @@ -140,11 +140,12 @@ class gl_spider : public virtual gl_texture setMatrix(); //create the transformation matrix. glCallList(dList+1); //move the templates to p, d, m. int best = getCost(ptexbufferID, numSamplesPos); //find min cost. + std::cerr << best << std::endl; stim::vec next( //find next position. - pV[best][0], - pV[best][1], - pV[best][2], - 1); + pV[best][0], + pV[best][1], + pV[best][2], + 1); next = cT*next; //find next position. setPosition( next[0]*S[0]*R[0], @@ -727,7 +728,7 @@ class gl_spider : public virtual gl_texture ///Best results if samples is can create a perfect root. ///Default Constructor gl_spider - (int samples = 1089, int samplespos = 400,int samplesmag = 144) + (int samples = 1089, int samplespos = 441,int samplesmag = 144) { p = vec(0.0, 0.0, 0.0); d = vec(0.0, 0.0, 1.0); @@ -750,7 +751,7 @@ class gl_spider : public virtual gl_texture ///@param int samples, number of templates this spider is going to use. gl_spider (float pos_x, float pos_y, float pos_z, float dir_x, float dir_y, float dir_z, - float mag_x, int numsamples = 1089, int numsamplespos = 400, int numsamplesmag =144) + float mag_x, int numsamples = 1089, int numsamplespos = 441, int numsamplesmag =144) { p = vec(pos_x, pos_y, pos_z); d = vec(dir_x, dir_y, dir_z); @@ -768,7 +769,7 @@ class gl_spider : public virtual gl_texture ///@param float mag, size of the vector. ///@param int samples, number of templates this spider is going to use. gl_spider - (stim::vec pos, stim::vec dir, float mag, int samples = 1089, int samplesPos = 400, int samplesMag = 144) + (stim::vec pos, stim::vec dir, float mag, int samples = 1089, int samplesPos = 441, int samplesMag = 144) { p = pos; d = dir; @@ -815,13 +816,13 @@ class gl_spider : public virtual gl_texture glListBase(dList); Bind(texbufferID, fboID, numSamples); genDirectionVectors(5*M_PI/4); -// Unbind(); - Bind(ptexbufferID, bfboID, numSamplesPos); + Unbind(); + Bind(ptexbufferID, pfboID, numSamplesPos); genPositionVectors(); -// Unbind(); + Unbind(); Bind(mtexbufferID, mfboID, numSamplesMag); genMagnitudeVectors(); -// Unbind(); + Unbind(); Bind(btexbufferID, bfboID, 27); DrawCylinder(); Unbind(); @@ -1150,13 +1151,13 @@ class gl_spider : public virtual gl_texture start = std::clock(); #endif findOptimalDirection(); - //Unbind(); + Unbind(); Bind(ptexbufferID, pfboID, numSamplesPos); findOptimalPosition(); - //Unbind(); + Unbind(); Bind(mtexbufferID, mfboID, numSamplesMag); findOptimalScale(); - //Unbind(); + Unbind(); CHECK_OPENGL_ERROR #ifdef TESTING @@ -1566,7 +1567,9 @@ class gl_spider : public virtual gl_texture } else { cL.push_back(stim::vec(p[0], p[1],p[2])); - cM.push_back(m[0]); + cM.push_back(stim::vec(m[0], m[0])); +// cM.push_back(m[0]); + sk.TexCoord(m[0]); sk.Vertex(p[0], p[1], p[2]); Bind(btexbufferID, bfboID, 27); diff --git a/stim/math/plane.h b/stim/math/plane.h index b009cad..baa4683 100644 --- a/stim/math/plane.h +++ b/stim/math/plane.h @@ -1,179 +1,213 @@ -#ifndef RTS_PLANE_H -#define RTS_PLANE_H +#ifndef STIM_PLANE_H +#define STIM_PLANE_H #include #include -#include "rts/cuda/callable.h" +#include +#include -namespace stim{ -template class plane; +namespace stim +{ +template class plane; } -template -CUDA_CALLABLE stim::plane operator-(stim::plane v); - -namespace stim{ - -template -class plane{ - - //a plane is defined by a point and a normal - -private: - - vec P; //point on the plane - vec N; //plane normal - - CUDA_CALLABLE void init(){ - P = vec(0, 0, 0); - N = vec(0, 0, 1); - } - - -public: - - //default constructor - CUDA_CALLABLE plane(){ - init(); - } +template +CUDA_CALLABLE stim::plane operator-(stim::plane v); + +namespace stim +{ + +template +class plane +{ + protected: + stim::vec P; + stim::vec N; + stim::vec U; + + ///Initializes the plane with standard coordinates. + /// + CUDA_CALLABLE void init() + { + P = stim::vec(0, 0, 0); + N = stim::vec(0, 0, 1); + U = stim::vec(1, 0, 0); + } + + public: - CUDA_CALLABLE plane(vec n, vec p = vec(0, 0, 0)){ - P = p; - N = n.norm(); - } - - CUDA_CALLABLE plane(T z_pos){ - init(); - P[2] = z_pos; - } - - //create a plane from three points (a triangle) - CUDA_CALLABLE plane(vec a, vec b, vec c){ - P = c; - N = (c - a).cross(b - a); - if(N.len() == 0) //handle the degenerate case when two vectors are the same, N = 0 - N = 0; - else - N = N.norm(); - } - - template< typename U > - CUDA_CALLABLE operator plane(){ - - plane result(N, P); - return result; - } - - CUDA_CALLABLE vec norm(){ - return N; - } - - CUDA_CALLABLE vec p(){ - return P; - } - - //flip the plane front-to-back - CUDA_CALLABLE plane flip(){ - plane result = *this; - result.N = -result.N; - return result; - } - - //determines how a vector v intersects the plane (1 = intersects front, 0 = within plane, -1 = intersects back) - CUDA_CALLABLE int face(vec v){ - - T dprod = v.dot(N); //get the dot product between v and N - - //conditional returns the appropriate value - if(dprod < 0) - return 1; - else if(dprod > 0) - return -1; - else - return 0; - } - - //determine on which side of the plane a point lies (1 = front, 0 = on the plane, -1 = back) - CUDA_CALLABLE int side(vec p){ - - vec v = p - P; //get the vector from P to the query point p - - return face(v); - } - - //compute the component of v that is perpendicular to the plane - CUDA_CALLABLE vec perpendicular(vec v){ - return N * v.dot(N); - } - - //compute the projection of v in the plane - CUDA_CALLABLE vec parallel(vec v){ - return v - perpendicular(v); - } - - CUDA_CALLABLE void decompose(vec v, vec& para, vec& perp){ - perp = N * v.dot(N); - para = v - perp; - } - - //get both the parallel and perpendicular components of a vector v w.r.t. the plane - CUDA_CALLABLE void project(vec v, vec &v_par, vec &v_perp){ - - v_perp = v.dot(N); - v_par = v - v_perp; - } - - //compute the reflection of v off of the plane - CUDA_CALLABLE vec reflect(vec v){ - - //compute the reflection using N_prime as the plane normal - vec par = parallel(v); - vec r = (-v) + par * 2; - - /*std::cout<<"----------------REFLECT-----------------------------"< operator-() - { - rts::plane p = *this; - - //negate the normal vector - p.N = -p.N; - - return p; - } - - //output a string - std::string str(){ - std::stringstream ss; - ss<<"P: "< operator- <> (rts::plane v); - - + CUDA_CALLABLE plane() + { + init(); + } + + CUDA_CALLABLE plane(vec n, vec p = vec(0, 0, 0)) + { + init(); + P = p; + rotate(n.norm()); + } + + CUDA_CALLABLE plane(T z_pos) + { + init(); + P[2] = z_pos; + } + + //create a plane from three points (a triangle) + CUDA_CALLABLE plane(vec a, vec b, vec c) + { + init(); + P = c; + stim::vec n = (c - a).cross(b - a); + try + { + if(n.len() != 0) + { + rotate(n.norm()); + } else { + throw 42; + } + } + catch(int i) + { + std::cerr << "No plane can be creates as all points a,b,c lie on a straight line" << std::endl; + } + } + + template< typename U > + CUDA_CALLABLE operator plane() + { + plane result(N, P); + return result; + + } + + CUDA_CALLABLE vec n() + { + return N; + } + + CUDA_CALLABLE vec p() + { + return P; + } + + CUDA_CALLABLE vec u() + { + return U; + } + + ///flip the plane front-to-back + CUDA_CALLABLE plane flip(){ + plane result = *this; + result.N = -result.N; + return result; + } + + //determines how a vector v intersects the plane (1 = intersects front, 0 = within plane, -1 = intersects back) + CUDA_CALLABLE int face(vec v){ + + T dprod = v.dot(N); //get the dot product between v and N + + //conditional returns the appropriate value + if(dprod < 0) + return 1; + else if(dprod > 0) + return -1; + else + return 0; + } + + //determine on which side of the plane a point lies (1 = front, 0 = on the plane, -1 = bac k) + CUDA_CALLABLE int side(vec p){ + + vec v = p - P; //get the vector from P to the query point p + + return face(v); + } + + //compute the component of v that is perpendicular to the plane + CUDA_CALLABLE vec perpendicular(vec v){ + return N * v.dot(N); + } + + //compute the projection of v in the plane + CUDA_CALLABLE vec parallel(vec v){ + return v - perpendicular(v); + } + + CUDA_CALLABLE void setU(vec v) + { + U = (parallel(v.norm())).norm(); + } + + CUDA_CALLABLE void decompose(vec v, vec& para, vec& perp){ + perp = N * v.dot(N); + para = v - perp; + } + + //get both the parallel and perpendicular components of a vector v w.r.t. the plane + CUDA_CALLABLE void project(vec v, vec &v_par, vec &v_perp){ + + v_perp = v.dot(N); + v_par = v - v_perp; + } + + //compute the reflection of v off of the plane + CUDA_CALLABLE vec reflect(vec v){ + + //compute the reflection using N_prime as the plane normal + vec par = parallel(v); + vec r = (-v) + par * 2; + return r; + + } + + CUDA_CALLABLE stim::plane operator-() + { + stim::plane p = *this; + + //negate the normal vector + p.N = -p.N; + return p; + } + + //output a string + std::string str(){ + std::stringstream ss; + ss<<"P: "< n) + { + quaternion q; + q.CreateRotation(N, n); + + N = q.toMatrix3() * N; + U = q.toMatrix3() * U; + + } + + CUDA_CALLABLE void rotate(vec n, vec &X, vec &Y) + { + quaternion q; + q.CreateRotation(N, n); + + N = q.toMatrix3() * N; + U = q.toMatrix3() * U; + X = q.toMatrix3() * X; + Y = q.toMatrix3() * Y; + + } }; - + + } - -//arithmetic operators - -//negative operator flips the plane (front to back) -//template - - - - #endif diff --git a/stim/math/rect.h b/stim/math/rect.h index 9767f26..969f88b 100644 --- a/stim/math/rect.h +++ b/stim/math/rect.h @@ -1,26 +1,28 @@ -#ifndef RTS_RECT_H -#define RTS_RECT_H +#ifndef STIM_RECT_H +#define STIM_RECT_H + //enable CUDA_CALLABLE macro #include +#include #include #include -#include #include #include #include +#include namespace stim{ //template for a rectangle class in ND space -template -struct rect +template +class rect : plane { /* ^ O | | - Y C + Y P | | O---------X---------> @@ -28,106 +30,143 @@ struct rect private: - stim::vec C; stim::vec X; stim::vec Y; - CUDA_CALLABLE void scale(T factor){ - X *= factor; - Y *= factor; - } - CUDA_CALLABLE void normal(vec n){ //orient the rectangle along the specified normal - - n = n.norm(); //normalize, just in case - vec n_current = X.cross(Y).norm(); //compute the current normal - quaternion q; //create a quaternion - q.CreateRotation(n_current, n); //initialize a rotation from n_current to n - - //apply the quaternion to the vectors and position - X = q.toMatrix3() * X; - Y = q.toMatrix3() * Y; - } - - CUDA_CALLABLE void init(){ - C = vec(0, 0, 0); - X = vec(1, 0, 0); - Y = vec(0, 1, 0); - } public: - CUDA_CALLABLE rect(){ + using stim::plane::n; + using stim::plane::P; + using stim::plane::N; + using stim::plane::U; + using stim::plane::rotate; + + ///base constructor. + CUDA_CALLABLE rect() + : plane() + { init(); } - //create a rectangle given a size and position - CUDA_CALLABLE rect(T size, T z_pos = (T)0){ + ///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 rect(T size, T z_pos = (T)0) + : plane(z_pos) + { init(); //use the default setup scale(size); //scale the rectangle - C[2] = z_pos; } - //create a rectangle from a center point, normal, and size - CUDA_CALLABLE rect(vec c, vec n = vec(0, 0, 1)){ + ///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 rect(vec c, vec n = vec(0, 0, 1)) + : plane() + { init(); //start with the default setting - C = c; normal(n); //orient } + ///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 rect(vec c, T s, vec n = vec(0, 0, 1)) + : plane() + { + init(); //start with the default setting + scale(s); + center(c); + rotate(n, X, Y); + } + + ///creates a rectangle from a centerpoint and an X and Y direction vectors. + ///@param center: x,y,z location of the center. + ///@param directionX: u,v,w direction of the X vector. + ///@param directionY: u,v,w direction of the Y vector. CUDA_CALLABLE rect(vec center, vec directionX, vec directionY ) + : plane((directionX.cross(directionY)).norm(),center) { - C = center; X = directionX; Y = directionY; } + ///creates a rectangle from a size, centerpoint, X, and Y direction vectors. + ///@param size of the rectangle in ND space. + ///@param center: x,y,z location of the center. + ///@param directionX: u,v,w direction of the X vector. + ///@param directionY: u,v,w direction of the Y vector. CUDA_CALLABLE rect(T size, vec center, vec directionX, vec directionY ) + : plane((directionX.cross(directionY)).norm(),center) { - C = center; X = directionX; Y = directionY; scale(size); } - - CUDA_CALLABLE rect(vec size, vec center, vec directionX, vec directionY ) + + ///creates a rectangle from a size, centerpoint, X, and Y direction vectors. + ///@param size of the rectangle in ND space, size[0] = size in X, size[1] = size in Y. + ///@param center: x,y,z location of the center. + ///@param directionX: u,v,w direction of the X vector. + ///@param directionY: u,v,w direction of the Y vector. + CUDA_CALLABLE rect(vec size, vec center, vec directionX, vec directionY) + : plane((directionX.cross(directionY)).norm(), center) { - C = center; X = directionX; Y = directionY; scale(size[0], size[1]); } - - CUDA_CALLABLE void scale(T factor1, T factor2){ + + CUDA_CALLABLE void scale(T factor){ + X *= factor; + Y *= factor; + } + + ///scales a rectangle in ND space. + ///@param factor1: size of the scale in the X-direction. + ///@param factor2: size of the scale in the Y-direction. + CUDA_CALLABLE void scale(T factor1, T factor2) + { X *= factor1; Y *= factor2; } + ///@param n; vector with the normal. + ///Orients the rectangle along the normal n. + CUDA_CALLABLE void normal(vec n) + { + //orient the rectangle along the specified normal + rotate(n, X, Y); + } + + ///general init method that sets a general rectangle. + CUDA_CALLABLE void init() + { + X = vec(1, 0, 0); + Y = vec(0, 1, 0); + } + //boolean comparison bool operator==(const rect & rhs) { - if(C == rhs.C && X == rhs.X && Y == rhs.Y) + if(P == rhs.P && X == rhs.X && Y == rhs.Y) return true; else return false; } - /******************************************* - Return the normal for the rect - *******************************************/ - CUDA_CALLABLE stim::vec n() - { - return (X.cross(Y)).norm(); - } //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; //given the two parameters a, b = [0 1], returns the position in world space - vec A = C - X * (T)0.5 - Y * (T)0.5; + vec A = this->P - X * (T)0.5 - Y * (T)0.5; result = A + X * a + Y * b; return result; @@ -142,16 +181,16 @@ public: std::string str() { std::stringstream ss; - vec A = C - X * (T)0.5 - Y * (T)0.5; + vec A = P - X * (T)0.5 - Y * (T)0.5; ss<"<<"C="<"<<"D="< operator*(T rhs) { //scales the plane by a scalar value @@ -164,36 +203,44 @@ public: } - //computes the distance between the specified point and this rectangle + ///computes the distance between the specified point and this rectangle. + ///@param p: x, y, z coordinates of the point to calculate distance to. CUDA_CALLABLE T dist(vec p) { //compute the distance between a point and this rect - vec A = C - X * (T)0.5 - Y * (T)0.5; + vec A = P - X * (T)0.5 - Y * (T)0.5; - //first break the rect up into two triangles - triangle T0(A, A+X, A+Y); - triangle T1(A+X+Y, A+X, A+Y); + //first break the rect up into two triangles + triangle T0(A, A+X, A+Y); + triangle T1(A+X+Y, A+X, A+Y); - T d0 = T0.dist(p); - T d1 = T1.dist(p); + T d0 = T0.dist(p); + T d1 = T1.dist(p); - if(d0 < d1) - return d0; - else - return d1; + if(d0 < d1) + return d0; + else + return d1; + } + + CUDA_CALLABLE T center(vec p) + { + this->P = p; } + ///Returns the maximum distance of the rectangle from a point p to the sides of the rectangle. + ///@param p: x, y, z point. CUDA_CALLABLE T dist_max(vec p) { - vec A = C - X * (T)0.5 - Y * (T)0.5; - T da = (A - p).len(); - T db = (A+X - p).len(); - T dc = (A+Y - p).len(); - T dd = (A+X+Y - p).len(); + vec A = P - X * (T)0.5 - Y * (T)0.5; + T da = (A - p).len(); + T db = (A+X - p).len(); + T dc = (A+Y - p).len(); + T dd = (A+X+Y - p).len(); - return std::max( da, std::max(db, std::max(dc, dd) ) ); + return std::max( da, std::max(db, std::max(dc, dd) ) ); } }; -- libgit2 0.21.4