codebase / stimlib

Commit 09b24449af2408093318ce4c0fcc419d78f040d7

Authored by Pavel Govyadinov
1 parent 3f15dade

changes to plane and additing of the old rts plane class

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