#ifndef RITKCUDARAYCASTERCOMMON_H__
#define RITKCUDARAYCASTERCOMMON_H__

#ifndef INFINITY
	#define INFINITY 1e16
#endif

#define DOT3(a,b) (a.x*b.x + a.y*b.y + a.z*b.z)
#define DOT4(a,b) (a.x*b.x + a.y*b.y + a.z*b.z + a.w*b.w)

#include "ThirdParty/cutil_math.h"

// addition
inline __host__ __device__ float3 operator*(float3 a, int3 b)
{
    return make_float3(a.x * b.x, a.y * b.y, a.z * b.z);
}


/// Struct that defines a ray
/// =========================
struct RayType
{
	float OpticalCenter[3];
	float Direction[3];
	float Alpha0;
	float Alpha1;
};


/// Compute a ray
/// ==================
__forceinline__ __device__
void
ComputeRayDirection(float i, float j, const float *InverseProjectionMatrix, float3 *RayDirection)
{
	// Compute ray direction
	RayDirection->x = InverseProjectionMatrix[0]*i + InverseProjectionMatrix[1]*j + InverseProjectionMatrix[2];
	RayDirection->y = InverseProjectionMatrix[3]*i + InverseProjectionMatrix[4]*j + InverseProjectionMatrix[5];
	RayDirection->z = InverseProjectionMatrix[6]*i + InverseProjectionMatrix[7]*j + InverseProjectionMatrix[8];

	// Normalize and revert ray direction
	float Mag = (RayDirection->x)*(RayDirection->x) + (RayDirection->y)*(RayDirection->y) + (RayDirection->z)*(RayDirection->z);
	Mag = -1.f/sqrtf(Mag);
	RayDirection->x *= Mag;
	RayDirection->y *= Mag;
	RayDirection->z *= Mag;
}


/// Intersect a volume with a ray
/// =============================
__forceinline__ __device__
void
IntersectVolume(const float3 *OpticalCenter, const float3 *RayDirection, const ritkCudaRayCasterVolumeParameters *VP, float2 *Alpha)
{
	// compute intersection of ray with all six bbox planes
	float3 invD = make_float3(1.f)/(*RayDirection);
	float3 tbot = invD * (VP->Origin - (*OpticalCenter));
	float3 ttop = invD * (VP->Origin + VP->Spacing*VP->Size - (*OpticalCenter));
    
	// re-order intersections to find smallest and largest on each axis
	float3 tmin = fminf(ttop, tbot);
	float3 tmax = fmaxf(ttop, tbot);

	// find the largest tmin and the smallest tmax
	float largest_tmin = fmaxf(0.f,fmaxf(fmaxf(tmin.x, tmin.y), fmaxf(tmin.x, tmin.z)));
	float smallest_tmax = fminf(fminf(tmax.x, tmax.y), fminf(tmax.x, tmax.z));

	Alpha->x = largest_tmin;
	Alpha->y = smallest_tmax;
}


#endif // RITKCUDARAYCASTERCOMMON_H__