codebase / stimlib

  #ifndef STIM_PLANE_H
  #define STIM_PLANE_H
  

  #include <iostream>
  #include <stim/math/vector.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::vec<T> P;
  		stim::vec<T> N;

  		stim::vec<T> U;

  		///Initializes the plane with standard coordinates.
  		///

  		CUDA_CALLABLE void init()
  		{
  			P = stim::vec<T>(0, 0, 0);
  			N = stim::vec<T>(0, 0, 1);

  			U = stim::vec<T>(1, 0, 0);

  		}
  
  	public:

  		CUDA_CALLABLE plane()
  		{
  			init();
  		}

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

  		{
  			return N;
  		}
  
  		CUDA_CALLABLE vec<T> p()
  		{
  			return P;
  		}
  
  		CUDA_CALLABLE vec<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(vec<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(vec<T> p){
  
  			vec<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 vec<T> perpendicular(vec<T> v){
  			return N * v.dot(N);
  		}
  
  		//compute the projection of v in the plane
  		CUDA_CALLABLE vec<T> parallel(vec<T> v){
  			return v - perpendicular(v);
  		}
  

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

  		CUDA_CALLABLE void decompose(vec<T> v, vec<T>& para, vec<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(vec<T> v, vec<T> &v_par, vec<T> &v_perp){
  
  			v_perp = v.dot(N);
  			v_par = v - v_perp;
  		}
  
  		//compute the reflection of v off of the plane
  		CUDA_CALLABLE vec<T> reflect(vec<T> v){
  
  			//compute the reflection using N_prime as the plane normal
  			vec<T> par = parallel(v);
  			vec<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(vec<T> n)
  		{
  			quaternion<T> q;
  			q.CreateRotation(N, n);
  			
  			N = q.toMatrix3() * N;
  			U = q.toMatrix3() * U;
  
  		}
  

  		CUDA_CALLABLE void rotate(vec<T> n, vec<T> &X, vec<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