path: root/dep/g3dlite/include/G3D/Ray.h

/**
 @file Ray.h
 
 Ray class
 
 @maintainer Morgan McGuire, http://graphics.cs.williams.edu
 
 @created 2002-07-12
 @edited  2009-06-29
 */

#ifndef G3D_Ray_h
#define G3D_Ray_h

#include "G3D/platform.h"
#include "G3D/Vector3.h"
#include "G3D/Triangle.h"

namespace G3D {

/**
 A 3D Ray.
 */
class Ray {
private:
    friend class Intersect;

    Point3  m_origin;

    /** Unit length */
    Vector3 m_direction;

    /** 1.0 / direction */
    Vector3 m_invDirection;

    
    /** The following are for the "ray slope" optimization from
      "Fast Ray / Axis-Aligned Bounding Box Overlap Tests using Ray Slopes" 
      by Martin Eisemann, Thorsten Grosch, Stefan M��ller and Marcus Magnor
      Computer Graphics Lab, TU Braunschweig, Germany and
      University of Koblenz-Landau, Germany */
    enum Classification {MMM, MMP, MPM, MPP, PMM, PMP, PPM, PPP, POO, MOO, OPO, OMO, OOP, OOM, OMM, OMP, OPM, OPP, MOM, MOP, POM, POP, MMO, MPO, PMO, PPO};    

    Classification classification;

    /** ray slope */
    float ibyj, jbyi, kbyj, jbyk, ibyk, kbyi;

    /** Precomputed components */
    float c_xy, c_xz, c_yx, c_yz, c_zx, c_zy;

public:
    /** \param direction Assumed to have unit length */
    void set(const Point3& origin, const Vector3& direction);

    const Point3& origin() const {
        return m_origin;
    }
    
    /** Unit direction vector. */
    const Vector3& direction() const {
        return m_direction;
    }

    /** Component-wise inverse of direction vector.  May have inf() components */
    const Vector3& invDirection() const {
        return m_invDirection;
    }
    
    Ray() {
        set(Point3::zero(), Vector3::unitX());
    }

    /** \param direction Assumed to have unit length */
    Ray(const Point3& origin, const Vector3& direction) {
        set(origin, direction);
    }

    Ray(class BinaryInput& b);
    
    void serialize(class BinaryOutput& b) const;
    void deserialize(class BinaryInput& b);
    
    /**
       Creates a Ray from a origin and a (nonzero) unit direction.
    */
    static Ray fromOriginAndDirection(const Point3& point, const Vector3& direction) {
        return Ray(point, direction);
    }
    
    /** Returns a new ray which has the same direction but an origin
        advanced along direction by @a distance */
    Ray bumpedRay(float distance) const {
        return Ray(m_origin + m_direction * distance, m_direction);
    }

    /** Returns a new ray which has the same direction but an origin
        advanced by \a distance * \a bumpDirection */
    Ray bumpedRay(float distance, const Vector3& bumpDirection) const {
        return Ray(m_origin + bumpDirection * distance, m_direction);
    }

    /**
       \brief Returns the closest point on the Ray to point.
    */
    Point3 closestPoint(const Point3& point) const {
        float t = m_direction.dot(point - m_origin);
        if (t < 0) {
            return m_origin;
        } else {
            return m_origin + m_direction * t;
        }
    }

    /**
     Returns the closest distance between point and the Ray
     */
    float distance(const Point3& point) const {
        return (closestPoint(point) - point).magnitude();
    }

    /**
     Returns the point where the Ray and plane intersect.  If there
     is no intersection, returns a point at infinity.

      Planes are considered one-sided, so the ray will not intersect
      a plane where the normal faces in the traveling direction.
    */
    Point3 intersection(const class Plane& plane) const;

    /**
     Returns the distance until intersection with the sphere or the (solid) ball bounded by the sphere.
     Will be 0 if inside the sphere, inf if there is no intersection.

     The ray direction is <B>not</B> normalized.  If the ray direction
     has unit length, the distance from the origin to intersection
     is equal to the time.  If the direction does not have unit length,
     the distance = time * direction.length().

     See also the G3D::CollisionDetection "movingPoint" methods,
     which give more information about the intersection.

     \param solid If true, rays inside the sphere immediately intersect (good for collision detection).  If false, they hit the opposite side of the sphere (good for ray tracing).
     */
    float intersectionTime(const class Sphere& sphere, bool solid = false) const;

    float intersectionTime(const class Plane& plane) const;

    float intersectionTime(const class Box& box) const;

    float intersectionTime(const class AABox& box) const;

    /**
     The three extra arguments are the weights of vertices 0, 1, and 2
     at the intersection point; they are useful for texture mapping
     and interpolated normals.
     */
    float intersectionTime(
        const Vector3& v0, const Vector3& v1, const Vector3& v2,
        const Vector3& edge01, const Vector3& edge02,
        float& w0, float& w1, float& w2) const;

     /**
     Ray-triangle intersection for a 1-sided triangle.  Fastest version.
       @cite http://www.acm.org/jgt/papers/MollerTrumbore97/
       http://www.graphics.cornell.edu/pubs/1997/MT97.html
     */
    float intersectionTime(
        const Point3& vert0,
        const Point3& vert1,
        const Point3& vert2,
        const Vector3& edge01,
        const Vector3& edge02) const;


    float intersectionTime(
        const Point3& vert0,
        const Point3& vert1,
        const Point3& vert2) const {

        return intersectionTime(vert0, vert1, vert2, vert1 - vert0, vert2 - vert0);
    }


    float intersectionTime(
        const Point3&  vert0,
        const Point3&  vert1,
        const Point3&  vert2,
        float&         w0,
        float&         w1,
        float&         w2) const {

        return intersectionTime(vert0, vert1, vert2, vert1 - vert0, vert2 - vert0, w0, w1, w2);
    }


    /* One-sided triangle 
       */
    float intersectionTime(const Triangle& triangle) const {
        return intersectionTime(
            triangle.vertex(0), triangle.vertex(1), triangle.vertex(2),
            triangle.edge01(), triangle.edge02());
    }


    float intersectionTime
    (const Triangle& triangle,
     float&         w0,
     float&         w1,
     float&         w2) const {
        return intersectionTime(triangle.vertex(0), triangle.vertex(1), triangle.vertex(2),
            triangle.edge01(), triangle.edge02(), w0, w1, w2);
    }

    /** Refracts about the normal
        using G3D::Vector3::refractionDirection
        and bumps the ray slightly from the newOrigin. */
    Ray refract(
        const Vector3&  newOrigin,
        const Vector3&  normal,
        float           iInside,
        float           iOutside) const;

    /** Reflects about the normal
        using G3D::Vector3::reflectionDirection
        and bumps the ray slightly from
        the newOrigin. */
    Ray reflect(
        const Vector3&  newOrigin,
        const Vector3&  normal) const;
};


#define EPSILON 0.000001
#define CROSS(dest,v1,v2) \
          dest[0]=v1[1]*v2[2]-v1[2]*v2[1]; \
          dest[1]=v1[2]*v2[0]-v1[0]*v2[2]; \
          dest[2]=v1[0]*v2[1]-v1[1]*v2[0];

#define DOT(v1,v2) (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2])

#define SUB(dest,v1,v2) \
          dest[0]=v1[0]-v2[0]; \
          dest[1]=v1[1]-v2[1]; \
          dest[2]=v1[2]-v2[2]; 

inline float Ray::intersectionTime(
    const Point3& vert0,
    const Point3& vert1,
    const Point3& vert2,
    const Vector3& edge1,
    const Vector3& edge2) const {

    (void)vert1;
    (void)vert2;

    // Barycenteric coords
    float u, v;

    float tvec[3], pvec[3], qvec[3];
    
    // begin calculating determinant - also used to calculate U parameter
    CROSS(pvec, m_direction, edge2);
    
    // if determinant is near zero, ray lies in plane of triangle
    const float det = DOT(edge1, pvec);
    
    if (det < EPSILON) {
        return finf();
    }
    
    // calculate distance from vert0 to ray origin
    SUB(tvec, m_origin, vert0);
    
    // calculate U parameter and test bounds
    u = DOT(tvec, pvec);
    if ((u < 0.0f) || (u > det)) {
        // Hit the plane outside the triangle
        return finf();
    }
    
    // prepare to test V parameter
    CROSS(qvec, tvec, edge1);
    
    // calculate V parameter and test bounds
    v = DOT(m_direction, qvec);
    if ((v < 0.0f) || (u + v > det)) {
        // Hit the plane outside the triangle
        return finf();
    }
    

    // Case where we don't need correct (u, v):
    const float t = DOT(edge2, qvec);
    
    if (t >= 0.0f) {
        // Note that det must be positive
        return t / det;
    } else {
        // We had to travel backwards in time to intersect
        return finf();
    }
}


inline float Ray::intersectionTime
(const Point3&  vert0,
 const Point3&  vert1,
 const Point3&  vert2,
 const Vector3&  edge1,
 const Vector3&  edge2,
 float&         w0,
 float&         w1,
 float&         w2) const {

    (void)vert1;
    (void)vert2;

    // Barycenteric coords
    float u, v;

    float tvec[3], pvec[3], qvec[3];

    // begin calculating determinant - also used to calculate U parameter
    CROSS(pvec, m_direction, edge2);
    
    // if determinant is near zero, ray lies in plane of triangle
    const float det = DOT(edge1, pvec);
    
    if (det < EPSILON) {
        return finf();
    }
    
    // calculate distance from vert0 to ray origin
    SUB(tvec, m_origin, vert0);
    
    // calculate U parameter and test bounds
    u = DOT(tvec, pvec);
    if ((u < 0.0f) || (u > det)) {
        // Hit the plane outside the triangle
        return finf();
    }
    
    // prepare to test V parameter
    CROSS(qvec, tvec, edge1);
    
    // calculate V parameter and test bounds
    v = DOT(m_direction, qvec);
    if ((v < 0.0f) || (u + v > det)) {
        // Hit the plane outside the triangle
        return finf();
    }
    
    float t = DOT(edge2, qvec);
    
    if (t >= 0) {
        const float inv_det = 1.0f / det;
        t *= inv_det;
        u *= inv_det;
        v *= inv_det;

        w0 = (1.0f - u - v);
        w1 = u;
        w2 = v;

        return t;
    } else {
        // We had to travel backwards in time to intersect
        return finf();
    }
}

#undef EPSILON
#undef CROSS
#undef DOT
#undef SUB

}// namespace

#endif