stim/math/rect.h

#ifndef RTS_RECT_H
#define RTS_RECT_H

//enable CUDA_CALLABLE macro
#include "../cuda/callable.h"
#include "../math/vector.h"
#include "../math/triangle.h"
#include "../math/quaternion.h"
#include <iostream>
#include <iomanip>
#include <algorithm>

namespace stim{

//template for a rectangle class in ND space
template <class T, int N = 3>
struct rect
{
	/*
		^                   O
		|                   
		|                   
		Y         C         
		|                   
		|                   
		O---------X--------->
	*/

private:

	stim::vec<T, N> C;
	stim::vec<T, N> X;
	stim::vec<T, N> Y;

	CUDA_CALLABLE void scale(T factor){
		X *= factor;
		Y *= factor;
	}
	

	CUDA_CALLABLE void normal(vec<T, N> n){		//orient the rectangle along the specified normal

		n = n.norm();								//normalize, just in case
		vec<T, N> n_current = X.cross(Y).norm();	//compute the current normal
		quaternion<T> q;							//create a quaternion
		q.CreateRotation(n_current, n);				//initialize a rotation from n_current to n

		//apply the quaternion to the vectors and position
		X = q.toMatrix3() * X;
		Y = q.toMatrix3() * Y;
	}

	CUDA_CALLABLE void init(){
		C = vec<T, N>(0, 0, 0);
		X = vec<T, N>(1, 0, 0);
		Y = vec<T, N>(0, 1, 0);
	}

public:

	CUDA_CALLABLE rect(){
		init();
	}

	//create a rectangle given a size and position
	CUDA_CALLABLE rect(T size, T z_pos = (T)0){
		init();			//use the default setup
		scale(size);	//scale the rectangle
		C[2] = z_pos;
	}

	//create a rectangle from a center point, normal, and size
	CUDA_CALLABLE rect(T size, vec<T, N> c, vec<T, N> n = vec<T, N>(0, 0, 1)){
		init();			//start with the default setting
		C = c;
		scale(size);	//scale the rectangle
		normal(n);		//orient
	}
	
	//create a rectangle from a center point, normal, and size
	CUDA_CALLABLE rect(vec<T,2> size, vec<T, N> c, vec<T, N> n = vec<T, N>(0, 0, 1)){
		init();			//start with the default setting
		C = c;
		scale(size);	//scale the rectangle
		normal(n);		//orient
	}

	CUDA_CALLABLE rect(vec<T, N> center, vec<T, N> directionX, vec<T, N> directionY )
	{
		C = center;
		X = directionX;
		Y = directionY;
	}

	CUDA_CALLABLE rect(vec<T,N> mag, vec<T, N> center, vec<T, N> directionX, vec<T, N> directionY )
	{	
		C = center;
		X = directionX;
		Y = directionY;
		scale(mag[0], mag[1]);
	}
	
	CUDA_CALLABLE void scale(T factor1, T factor2){
		X *= factor1;
		Y *= factor2;
	}

	//boolean comparison
	bool operator==(const rect<T, N> & rhs)
	{
		if(C == rhs.C && X == rhs.X && Y == rhs.Y)
			return true;
		else
			return false;
	}

	/*******************************************
	Return the normal for the rect
	*******************************************/
	CUDA_CALLABLE stim::vec<T, N> n()
	{
        return (X.cross(Y)).norm();
	}

	//get the world space value given the planar coordinates a, b in [0, 1]
	CUDA_CALLABLE stim::vec<T, N> p(T a, T b)
	{
		stim::vec<T, N> result;
		//given the two parameters a, b = [0 1], returns the position in world space
		vec<T, N> A = C - X * (T)0.5 - Y * (T)0.5;
		result = A + X * a + Y * b;

		return result;
	}

	//parenthesis operator returns the world space given rectangular coordinates a and b in [0 1]
	CUDA_CALLABLE stim::vec<T, N> operator()(T a, T b)
	{
		return p(a, b);
	}

	std::string str()
	{
		std::stringstream ss;
		vec<T, N> A = C - X * (T)0.5 - Y * (T)0.5;
		ss<<std::left<<"B="<<std::setfill('-')<<std::setw(20)<<A + Y<<">"<<"C="<<A + Y + X<<std::endl;
		ss<<std::setfill(' ')<<std::setw(23)<<"|"<<"|"<<std::endl<<std::setw(23)<<"|"<<"|"<<std::endl;
		ss<<std::left<<"A="<<std::setfill('-')<<std::setw(20)<<A<<">"<<"D="<<A + X;

        return ss.str();

	}

	//scales the rectangle by a value rhs
	CUDA_CALLABLE rect<T, N> operator*(T rhs)
	{
		//scales the plane by a scalar value

		//create the new rectangle
		rect<T, N> result = *this;
		result.scale(rhs);

		return result;

	}

	//computes the distance between the specified point and this rectangle
	CUDA_CALLABLE T dist(vec<T, N> p)
	{
        //compute the distance between a point and this rect

		vec<T, N> A = C - X * (T)0.5 - Y * (T)0.5;

        //first break the rect up into two triangles
        triangle<T, N> T0(A, A+X, A+Y);
        triangle<T, N> T1(A+X+Y, A+X, A+Y);


        T d0 = T0.dist(p);
        T d1 = T1.dist(p);

        if(d0 < d1)
            return d0;
        else
            return d1;
	}

	CUDA_CALLABLE T dist_max(vec<T, N> p)
	{
		vec<T, N> A = C - X * (T)0.5 - Y * (T)0.5;
        T da = (A - p).len();
        T db = (A+X - p).len();
        T dc = (A+Y - p).len();
        T dd = (A+X+Y - p).len();

        return std::max( da, std::max(db, std::max(dc, dd) ) );
	}
};

}	//end namespace rts

template <typename T, int N>
std::ostream& operator<<(std::ostream& os, stim::rect<T, N> R)
{
    os<<R.str();
    return os;
}


#endif