codebase / stimlib

  #ifndef RTS_PLANE_H
  #define RTS_PLANE_H
  
  #include <iostream>

  #include <stim/math/vector.h>

  #include "rts/cuda/callable.h"
  
  

  namespace stim{

  template <typename T, int D> class plane;
  }
  
  template <typename T, int D>

  CUDA_CALLABLE stim::plane<T, D> operator-(stim::plane<T, D> v);

  namespace stim{

  
  template <class T, int D = 3>
  class plane{
  
  	//a plane is defined by a point and a normal
  
  private:
  
  	vec<T, D> P;	//point on the plane
  	vec<T, D> N;	//plane normal
  
  	CUDA_CALLABLE void init(){
  		P = vec<T, D>(0, 0, 0);
  		N = vec<T, D>(0, 0, 1);
  	}
  
  
  public:
  
  	//default constructor
  	CUDA_CALLABLE plane(){
  		init();
  	}
  	
  	CUDA_CALLABLE plane(vec<T, D> n, vec<T, D> p = vec<T, D>(0, 0, 0)){
  		P = p;
  		N = 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, D> a, vec<T, D> b, vec<T, D> c){
  		P = c;
  		N = (c - a).cross(b - a);
  		if(N.len() == 0)	//handle the degenerate case when two vectors are the same, N = 0
  			N = 0;
  		else
  			N = N.norm();
  	}
  

  	template< typename U >
  	CUDA_CALLABLE operator plane<U, D>(){
  
  		plane<U, D> result(N, P);
  		return result;
  	}
  

  	CUDA_CALLABLE vec<T, D> norm(){
  		return N;
  	}
  
  	CUDA_CALLABLE vec<T, D> p(){
  		return P;
  	}
  
  	//flip the plane front-to-back
  	CUDA_CALLABLE plane<T, D> flip(){
  		plane<T, D> 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, D> 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 = back)
  	CUDA_CALLABLE int side(vec<T, D> p){
  
  		vec<T, D> 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, D> perpendicular(vec<T, D> v){
  		return N * v.dot(N);
  	}
  
  	//compute the projection of v in the plane
  	CUDA_CALLABLE vec<T, D> parallel(vec<T, D> v){
  		return v - perpendicular(v);
  	}
  

  	CUDA_CALLABLE void decompose(vec<T, D> v, vec<T, D>& para, vec<T, D>& 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, D> v, vec<T, D> &v_par, vec<T, D> &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, D> reflect(vec<T, D> v){
  
  		//compute the reflection using N_prime as the plane normal
  		vec<T, D> par = parallel(v);
  		vec<T, D> r = (-v) + par * 2;
  
  		/*std::cout<<"----------------REFLECT-----------------------------"<<std::endl;
  		std::cout<<str()<<std::endl;
  		std::cout<<"v: "<<v<<std::endl;
  		std::cout<<"r: "<<r<<std::endl;
  		std::cout<<"Perpendicular: "<<perpendicular(v)<<std::endl;
  		std::cout<<"Parallel: "<<par<<std::endl;*/
  		return r;
  
  	}
  

  	CUDA_CALLABLE rts::plane<T, D> operator-()
  	{
  		rts::plane<T, D> 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;
  		return ss.str();
  	}
  
  	///////Friendship

  	//friend CUDA_CALLABLE rts::plane<T, D> operator- <> (rts::plane<T, D> v);

  
  
  
  };
  
  }
  
  //arithmetic operators
  
  //negative operator flips the plane (front to back)

  //template <typename T, int D>

  #endif