// *********************************************************
// Ray Class
// Author : Tamy Boubekeur (boubek@gmail.com).
// Copyright (C) 2010 Tamy Boubekeur.
// All rights reserved.
// *********************************************************

#include "Ray.h"

using namespace std;

static const unsigned int NUMDIM = 3, RIGHT = 0, LEFT = 1, MIDDLE = 2;
static const float EPSILON = 0.00001f;

bool Ray::intersect (const BoundingBox & bbox, Vec3Df & intersectionPoint) const
{
	//*
    const Vec3Df& min = bbox.getMin();
    const Vec3Df& max = bbox.getMax();
    Vec3Df dir = direction;
    for(unsigned int j = 0; j < 3; j++)
		if(dir[j] > -EPSILON && dir[j] < EPSILON)
            dir[j] = EPSILON;
    const Vec3Df& invdir = Vec3Df(1,1,1) / dir;
    float tmin=0, tmax=0;
    float tm = 0, tM = 0;

    for(unsigned int j = 0; j < 3; j++)
    {
        tm = (min[j] - origin[j]) * invdir[j];
        tM = (max[j] - origin[j]) * invdir[j];
        tmin = (j == 0 ? std::min(tm,tM) : std::max(std::min(tm,tM), tmin));
        tmax = (j == 0 ? std::max(tm,tM) : std::min(std::max(tm,tM), tmax));
    }

	if(tmax < 0 || tmin > tmax)
		return false;

	intersectionPoint = origin + (tmin < 0.0f ? tmax : tmin) * direction;
    return true;

    /*/
    const Vec3Df & minBb = bbox.getMin ();
    const Vec3Df & maxBb = bbox.getMax ();
    bool inside = true;
    unsigned int  quadrant[NUMDIM];
    register unsigned int i;
    unsigned int whichPlane;
    Vec3Df maxT;
    Vec3Df candidatePlane;

    for (i=0; i<NUMDIM; i++)
        if (origin[i] < minBb[i]) {
            quadrant[i] = LEFT;
            candidatePlane[i] = minBb[i];
            inside = false;
        } else if (origin[i] > maxBb[i]) {
            quadrant[i] = RIGHT;
            candidatePlane[i] = maxBb[i];
            inside = false;
        } else	{
            quadrant[i] = MIDDLE;
        }

    if (inside)	{
        intersectionPoint = origin;
        return (true);
    }

    for (i = 0; i < NUMDIM; i++)
        if (quadrant[i] != MIDDLE && direction[i] !=0.)
            maxT[i] = (candidatePlane[i]-origin[i]) / direction[i];
        else
            maxT[i] = -1.;

    whichPlane = 0;
    for (i = 1; i < NUMDIM; i++)
        if (maxT[whichPlane] < maxT[i])
            whichPlane = i;

    if (maxT[whichPlane] < 0.) return (false);
    for (i = 0; i < NUMDIM; i++)
        if (whichPlane != i) {
            intersectionPoint[i] = origin[i] + maxT[whichPlane] *direction[i];
            if (intersectionPoint[i] < minBb[i] || intersectionPoint[i] > maxBb[i])
                return (false);
        } else {
            intersectionPoint[i] = candidatePlane[i];
        }
    return (true);
    //*/
}

bool Ray::getRange(const BoundingBox& bbox, float& tmin, float& tmax)
{
	const Vec3Df& min = bbox.getMin();
	const Vec3Df& max = bbox.getMax();
	Vec3Df dir = direction;
	for(unsigned int j = 0; j < 3; j++)
		if(dir[j] > -EPSILON && dir[j] < EPSILON)
			dir[j] = EPSILON;
	const Vec3Df& invdir = Vec3Df(1,1,1) / dir;
	float tMin=0, tMax=0;
	float tm = 0, tM = 0;

	for(unsigned int j = 0; j < 3; j++)
	{
		tm = (min[j] - origin[j]) * invdir[j];
		tM = (max[j] - origin[j]) * invdir[j];
		tMin = (j == 0 ? std::min(tm,tM) : std::max(std::min(tm,tM), tMin));
		tMax = (j == 0 ? std::max(tm,tM) : std::min(std::max(tm,tM), tMax));
	}

	if(tMax < 0 || tMin > tMax)
		return false;

	tmin = std::max(0.0f, tMin);
	tmax = tMax;
	return true;
}

bool Ray::intersect(const Vec3Df& tri0, const Vec3Df& tri1, const Vec3Df& tri2, Vec3Df& ipt, float& u, float& v, float& t) const
{
    static const float EPSILON = 0.0000001f;
    Vec3Df edge1, edge2, tvec, pvec, qvec;
	float det, inv_det;

    /* find vectors for two edges sharing vert0 */
    edge1 = tri1 - tri0;
    edge2 = tri2 - tri0;

    /* begin calculating determinant - also used to calculate U parameter */
    pvec = Vec3Df::crossProduct(direction, edge2);

    /* if determinant is near zero, ray lies in plane of triangle */
    det = Vec3Df::dotProduct(edge1, pvec);
    if (det > -EPSILON && det < EPSILON)
        return false;
    inv_det = 1.0f / det;

    /* calculate distance from vert0 to ray origin */
    tvec = origin - tri0;

    /* calculate U parameter and test bounds */
    u = Vec3Df::dotProduct(tvec, pvec) * inv_det;
    if(u < 0.0f || u > 1.0f)
        return false;

    /* prepare to test V parameter */
    qvec = Vec3Df::crossProduct(tvec, edge1);

    /* calculate V parameter and test bounds */
    v = Vec3Df::dotProduct(direction, qvec) * inv_det;
    if(v < 0.0f || u + v > 1.0f)
        return false;

    /* calculate t, ray intersects triangle */
    t = Vec3Df::dotProduct(edge2, qvec) * inv_det;
    ipt = origin + t * direction;

    return true;
}

bool Ray::intersect(const Vec3Df& center, float radius, Vec3Df& ipoint, float& t) const
{
    const float a = direction.dot(direction);
    const float b = 2 * (origin - center).dot(direction);
    const float c = (origin - center).dot(origin - center) - radius * radius;
    const float det = b * b - 4 * a * c;

    if(det < 0.0f)
        return false;
    else
    {
        const float sDet = sqrt(det);
        const float t0 = (-b - sDet) / (2 * a);
        const float t1 = (-b + sDet) / (2 * a);
        const float tmin = std::min(t0,t1);
        const float tmax = std::max(t0,t1);

        if(tmax < 0.0f)
            return false;

        t = (tmin < 0.0f ? tmax : tmin);
        ipoint = origin + t * direction;
        return true;
    }
}

float Ray::getT(const Vec3Df& pt) const
{
    static const float minVal = 0.0000001f;
    for(unsigned int j = 0; j < 3; j++)
        if(direction[j] <= -minVal || direction[j] >= minVal)
            return (pt[j] - origin[j]) / direction[j];

    return (pt[0] - origin[0]) / minVal;
}
