plane.h 4.01 KB
``````#ifndef STIM_PLANE_H
#define STIM_PLANE_H

#include <iostream>
#include <stim/math/vec3.h>
#include <stim/cuda/cudatools/callable.h>
#include <stim/math/quaternion.h>

namespace stim
{
template<typename T> class plane;
}

template<typename T>
CUDA_CALLABLE stim::plane<T> operator-(stim::plane<T> v);

namespace stim
{

template <typename T>
class plane
{
protected:
stim::vec3<T> P;
stim::vec3<T> N;
stim::vec3<T> U;

///Initializes the plane with standard coordinates.
///
CUDA_CALLABLE void init()
{
P = stim::vec3<T>(0, 0, 0);
N = stim::vec3<T>(0, 0, 1);
U = stim::vec3<T>(1, 0, 0);
}

public:

CUDA_CALLABLE plane()
{
init();
}

CUDA_CALLABLE plane(vec3<T> n, vec3<T> p = vec3<T>(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(vec3<T> a, vec3<T> b, vec3<T> c)
{
init();
P = c;
stim::vec3<T> 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<U>()
{
plane<U> result(N, P);
return result;

}

CUDA_CALLABLE vec3<T> n()
{
return N;
}

CUDA_CALLABLE vec3<T> p()
{
return P;
}

CUDA_CALLABLE vec3<T> u()
{
return U;
}

///flip the plane front-to-back
CUDA_CALLABLE plane<T> flip(){
plane<T> 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(vec3<T> 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(vec3<T> p){

vec3<T> 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 vec3<T> perpendicular(vec3<T> v){
return N * v.dot(N);
}

//compute the projection of v in the plane
CUDA_CALLABLE vec3<T> parallel(vec3<T> v){
return v - perpendicular(v);
}

CUDA_CALLABLE void setU(vec3<T> v)
{
U = (parallel(v.norm())).norm();
}

CUDA_CALLABLE void decompose(vec3<T> v, vec3<T>& para, vec3<T>& 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(vec3<T> v, vec3<T> &v_par, vec3<T> &v_perp){

v_perp = v.dot(N);
v_par = v - v_perp;
}

//compute the reflection of v off of the plane
CUDA_CALLABLE vec3<T> reflect(vec3<T> v){

//compute the reflection using N_prime as the plane normal
vec3<T> par = parallel(v);
vec3<T> r = (-v) + par * 2;
return r;

}

CUDA_CALLABLE stim::plane<T> operator-()
{
stim::plane<T> p = *this;

//negate the normal vector
p.N = -p.N;
return p;
}

//output a string
std::string str(){
std::stringstream ss;
ss<<"P: "<<P<<std::endl;
ss<<"N: "<<N<<std::endl;
ss<<"U: "<<U;
return ss.str();
}

CUDA_CALLABLE void rotate(vec3<T> n)
{
quaternion<T> q;
q.CreateRotation(N, n);
matrix_sq<T, 3> M = q.toMatrix3();
N = M * N;
U = M * U;

}

CUDA_CALLABLE void rotate(vec3<T> n, vec3<T> &Y)
{
quaternion<T> q;
q.CreateRotation(N, n);

N = q.toMatrix3() * N;
U = q.toMatrix3() * U;
Y = q.toMatrix3() * Y;

}

CUDA_CALLABLE void rotate(vec3<T> n, vec3<T> &X, vec3<T> &Y)
{
quaternion<T> q;
q.CreateRotation(N, n);

N = q.toMatrix3() * N;
U = q.toMatrix3() * U;
X = q.toMatrix3() * X;
Y = q.toMatrix3() * Y;

}

};

}
#endif
``````