codebase / stimlib

Commit 59360849ce3cf589292a26067b2297dd7839b772

Authored by Pavel Govyadinov
1 parent 09b24449

test changes to the plane class

Showing 1 changed file with 157 additions and 0 deletions   Show diff stats
Inline Side-by-side
stim/math/plane.h
  Show/Hide comments   View file @5936084
1 1 #ifndef STIM_PLANE_H
2 2 #define STIM_PLANE_H
3 3  
  4 +#include <iostream>
  5 +#include <stim/math/vector.h>
  6 +#include <stim/cuda/cudatools/callable.h>
  7 +#include <stim/math/quaternion.h>
  8 +
4 9 namespace stim
5 10 {
6 11 template <typename T> class plane
... ... @@ -8,11 +13,15 @@ template &lt;typename T&gt; class plane
8 13 private:
9 14 stim::vec<T> P;
10 15 stim::vec<T> N;
  16 + stim::vec<T> U;
11 17  
  18 + ///Initializes the plane with standard coordinates.
  19 + ///
12 20 CUDA_CALLABLE void init()
13 21 {
14 22 P = stim::vec<T>(0, 0, 0);
15 23 N = stim::vec<T>(0, 0, 1);
  24 + U = stim::vec<T>(1, 0, 0);
16 25 }
17 26  
18 27 public:
... ... @@ -21,6 +30,154 @@ template &lt;typename T&gt; class plane
21 30 {
22 31 init();
23 32 }
  33 +
  34 + CUDA_CALLABLE plane(vec<T> n, vec<T> p = vec<T>(0, 0, 0))
  35 + {
  36 + init();
  37 + P = p;
  38 + rotate(N, n.norm());
  39 + }
  40 +
  41 + CUDA_CALLABLE plane(T z_pos)
  42 + {
  43 + init();
  44 + P[2] = z_pos;
  45 + }
  46 +
  47 + //create a plane from three points (a triangle)
  48 + CUDA_CALLABLE plane(vec<T> a, vec<T> b, vec<T> c)
  49 + {
  50 + init();
  51 + P = c;
  52 + stim::vec<T> n = (c - a).cross(b - a);
  53 + try
  54 + {
  55 + if(n.len() != 0)
  56 + {
  57 + rotate(N, n.norm());
  58 + }
  59 + }
  60 + catch
  61 + {
  62 + std::cerr << "No plane can be creates as all points a,b,c lie on a straight line" << std::endl;
  63 + }
  64 + }
  65 +
  66 + template< typename U >
  67 + CUDA_CALLABLE operator plane<U>()
  68 + {
  69 + plane<U> result(N, P);
  70 + return result;
  71 +
  72 + }
  73 +
  74 + CUDA_CALLABLE vec<T> norm()
  75 + {
  76 + return N;
  77 + }
  78 +
  79 + CUDA_CALLABLE vec<T> p()
  80 + {
  81 + return P;
  82 + }
  83 +
  84 + CUDA_CALLABLE vec<T> u()
  85 + {
  86 + return U;
  87 + }
  88 +
  89 + ///flip the plane front-to-back
  90 + CUDA_CALLABLE plane<T> flip(){
  91 + plane<T> result = *this;
  92 + result.N = -result.N;
  93 + return result;
  94 + }
  95 +
  96 + //determines how a vector v intersects the plane (1 = intersects front, 0 = within plane, -1 = intersects back)
  97 + CUDA_CALLABLE int face(vec<T> v){
  98 +
  99 + T dprod = v.dot(N); //get the dot product between v and N
  100 +
  101 + //conditional returns the appropriate value
  102 + if(dprod < 0)
  103 + return 1;
  104 + else if(dprod > 0)
  105 + return -1;
  106 + else
  107 + return 0;
  108 + }
  109 +
  110 + //determine on which side of the plane a point lies (1 = front, 0 = on the plane, -1 = bac k)
  111 + CUDA_CALLABLE int side(vec<T> p){
  112 +
  113 + vec<T> v = p - P; //get the vector from P to the query point p
  114 +
  115 + return face(v);
  116 + }
  117 +
  118 + //compute the component of v that is perpendicular to the plane
  119 + CUDA_CALLABLE vec<T> perpendicular(vec<T> v){
  120 + return N * v.dot(N);
  121 + }
  122 +
  123 + //compute the projection of v in the plane
  124 + CUDA_CALLABLE vec<T> parallel(vec<T> v){
  125 + return v - perpendicular(v);
  126 + }
  127 +
  128 + CUDA_CALLABLE void decompose(vec<T> v, vec<T>& para, vec<T>& perp){
  129 + perp = N * v.dot(N);
  130 + para = v - perp;
  131 + }
  132 +
  133 + //get both the parallel and perpendicular components of a vector v w.r.t. the plane
  134 + CUDA_CALLABLE void project(vec<T> v, vec<T> &v_par, vec<T> &v_perp){
  135 +
  136 + v_perp = v.dot(N);
  137 + v_par = v - v_perp;
  138 + }
  139 +
  140 + //compute the reflection of v off of the plane
  141 + CUDA_CALLABLE vec<T> reflect(vec<T> v){
  142 +
  143 + //compute the reflection using N_prime as the plane normal
  144 + vec<T> par = parallel(v);
  145 + vec<T> r = (-v) + par * 2;
  146 + return r;
  147 +
  148 + }
  149 +
  150 + CUDA_CALLABLE stim::plane<T> operator-()
  151 + {
  152 + rts::plane<T> p = *this;
  153 +
  154 + //negate the normal vector
  155 + p.N = -p.N;
  156 + return p;
  157 + }
  158 +
  159 + //output a string
  160 + std::string str(){
  161 + std::stringstream ss;
  162 + ss<<"P: "<<P<<std::endl;
  163 + ss<<"N: "<<N<<std::endl;
  164 + ss<<"D: "<<D;
  165 + return ss.str();
  166 + }
  167 +
  168 +
  169 + CUDA_CALLABLE void rotate(vec<T> n)
  170 + {
  171 + quaternion<T> q;
  172 + q.CreateRotation(N, n);
  173 +
  174 + N = q.toMatrix3() * N;
  175 + U = q.toMatrix3() * U;
  176 +
  177 + }
  178 +
  179 +};
  180 +
24 181  
25 182 }
26 183 #endif
... ...