### Blame view

stim/math/plane.h 4.01 KB
 1 2 #ifndef STIM_PLANE_H #define STIM_PLANE_H 3 4 #include 5 #include 6 7 #include #include 8 9 10 11 12 namespace stim { template class plane; 13 14 } 15 16 17 18 19 20 21 22 23 24 template CUDA_CALLABLE stim::plane operator-(stim::plane v); namespace stim { template class plane { protected: 25 26 27 stim::vec3 P; stim::vec3 N; stim::vec3 U; 28 29 30 31 32 ///Initializes the plane with standard coordinates. /// CUDA_CALLABLE void init() { 33 34 35 P = stim::vec3(0, 0, 0); N = stim::vec3(0, 0, 1); U = stim::vec3(1, 0, 0); 36 37 38 } public: 39 40 41 42 43 44 CUDA_CALLABLE plane() { init(); } 45 CUDA_CALLABLE plane(vec3 n, vec3 p = vec3(0, 0, 0)) 46 47 48 49 50 51 52 53 54 55 56 57 58 { 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) 59 CUDA_CALLABLE plane(vec3 a, vec3 b, vec3 c) 60 61 62 { init(); P = c; 63 stim::vec3 n = (c - a).cross(b - a); 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 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; } 87 CUDA_CALLABLE vec3 n() 88 89 90 91 { return N; } 92 CUDA_CALLABLE vec3 p() 93 94 95 96 { return P; } 97 CUDA_CALLABLE vec3 u() 98 99 100 101 102 103 104 105 106 107 108 109 { 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) 110 CUDA_CALLABLE int face(vec3 v){ 111 112 113 114 115 116 117 118 119 120 121 122 123 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) 124 CUDA_CALLABLE int side(vec3 p){ 125 126 vec3 v = p - P; //get the vector from P to the query point p 127 128 129 130 131 return face(v); } //compute the component of v that is perpendicular to the plane 132 CUDA_CALLABLE vec3 perpendicular(vec3 v){ 133 134 135 136 return N * v.dot(N); } //compute the projection of v in the plane 137 CUDA_CALLABLE vec3 parallel(vec3 v){ 138 139 140 return v - perpendicular(v); } 141 CUDA_CALLABLE void setU(vec3 v) 142 143 144 145 { U = (parallel(v.norm())).norm(); } 146 CUDA_CALLABLE void decompose(vec3 v, vec3& para, vec3& perp){ 147 148 149 150 151 perp = N * v.dot(N); para = v - perp; } //get both the parallel and perpendicular components of a vector v w.r.t. the plane 152 CUDA_CALLABLE void project(vec3 v, vec3 &v_par, vec3 &v_perp){ 153 154 155 156 157 158 v_perp = v.dot(N); v_par = v - v_perp; } //compute the reflection of v off of the plane 159 CUDA_CALLABLE vec3 reflect(vec3 v){ 160 161 //compute the reflection using N_prime as the plane normal 162 163 vec3 par = parallel(v); vec3 r = (-v) + par * 2; 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 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) 188 189 190 { quaternion q; q.CreateRotation(N, n); 191 192 193 matrix_sq M = q.toMatrix3(); N = M * N; U = M * U; 194 195 196 } 197 CUDA_CALLABLE void rotate(vec3 n, vec3 &Y) 198 199 200 201 202 203 204 205 206 207 { quaternion q; q.CreateRotation(N, n); N = q.toMatrix3() * N; U = q.toMatrix3() * U; Y = q.toMatrix3() * Y; } 208 CUDA_CALLABLE void rotate(vec3 n, vec3 &X, vec3 &Y) 209 210 211 212 213 214 215 216 217 218 { quaternion q; q.CreateRotation(N, n); N = q.toMatrix3() * N; U = q.toMatrix3() * U; X = q.toMatrix3() * X; Y = q.toMatrix3() * Y; } 219 220 }; 221 222 223 } 224 #endif