RayTriangleIntersect.h

Go to the documentation of this file.

/*  Copyright (C) 2003-2021 Free Electron Organization
    Any use of this software requires a license.  If a valid license
    was not distributed with this file, visit freeelectron.org. */

/** @file */

/* This file is copyrighted only to the extent of adapting it to FE.
   Substantial source from the original example remains, for which
   no copyright was indicated. */

#ifndef __solve_RayTriangleIntersect_h__
#define __solve_RayTriangleIntersect_h__

#define FE_RTI_DEBUG            FALSE

#define FE_MOLLER_EPSILON       (1e-6)
#define FE_MOLLER_TEST_CULL     FALSE   //* ignore backside using implied normal

namespace fe
{
namespace ext
{

/**************************************************************************//**
    @brief Find intersection between ray and triangle

    @ingroup solve

    http://jgt.akpeters.com/papers/MollerTrumbore97
*//***************************************************************************/
template <typename T>
class RayTriangleIntersect
{
    public:
static  T       solve(const Vector<3,T>& vert0,
                        const Vector<3,T>& vert1,
                        const Vector<3,T>& vert2,
                        const Vector<3,T>& origin,
                        const Vector<3,T>& direction,
                        Barycenter<T>& barycenter);
static  void    resolveContact(
                        const Vector<3,T>& vert0,
                        const Vector<3,T>& vert1,
                        const Vector<3,T>& vert2,
                        const Vector<3,T>& norm0,
                        const Vector<3,T>& norm1,
                        const Vector<3,T>& norm2,
                        const Vector<3,T>& origin,
                        const Vector<3,T>& direction,
                        const Barycenter<T>& barycenter,
                        const T range,
                        Vector<3,T>& intersection,
                        Vector<3,T>& normal);
};

// direction must be normalized
template <typename T>
inline T RayTriangleIntersect<T>::solve(const Vector<3,T>& vert0,
        const Vector<3,T>& vert1,const Vector<3,T>& vert2,
        const Vector<3,T>& origin,const Vector<3,T>& direction,
        Barycenter<T>& barycenter)
{
#if FE_RTI_DEBUG
    feLog("vert %s %s %s\n",c_print(vert0),c_print(vert1),c_print(vert2));
    feLog("origin %s\n",c_print(origin));
    feLog("direction %s\n",c_print(direction));
#endif

    /* find vectors for two edges sharing vert0 */
    Vector<3,T> edge1=vert1-vert0;
    Vector<3,T> edge2=vert2-vert0;

    /* begin calculating determinant - also used to calculate U parameter */
    Vector<3,T> pvec;
    cross3(pvec,direction,edge2);

    /* if determinant is near zero, ray lies in plane of triangle */
    T det=dot(edge1, pvec);

/*****************************************************************************/
#if FE_MOLLER_TEST_CULL     /* if culling is desired */
    if(det<T(FE_MOLLER_EPSILON))
    {
#if FE_RTI_DEBUG
        feLog(" miss: det=%.6G\n",det);
#endif
        return T(-1);
    }

    /* calculate distance from vert0 to ray origin */
    Vector<3,T> tvec=origin-vert0;

    /* calculate U parameter and test bounds */
    T u=dot(tvec,pvec);
    if(u<T(0) || u>det)
    {
#if FE_RTI_DEBUG
        feLog(" miss: u=%.6G\n",u);
#endif
        return T(-1);
    }

    /* prepare to test V parameter */
    Vector<3,T> qvec;
    cross3(qvec,tvec,edge1);

    /* calculate V parameter and test bounds */
    T v=dot(direction,qvec);
    if(v<T(0) || u+v>det)
    {
#if FE_RTI_DEBUG
        feLog(" miss: u=%.6G v=%.6G\n",u,v);
#endif
        return T(-1);
    }

    /* calculate t, scale parameters, ray intersects triangle */
    T range=dot(edge2,qvec);
    T inv_det=T(1)/det;
    range*=inv_det;
    u*=inv_det;
    v*=inv_det;

/*****************************************************************************/
#else                           /* the non-culling branch */
    if(det>-T(FE_MOLLER_EPSILON) && det<T(FE_MOLLER_EPSILON))
    {
#if FE_RTI_DEBUG
        feLog(" miss: det=%.6G\n",det);
#endif
        return T(-1);
    }
    T inv_det=T(1)/det;

    /* calculate distance from vert0 to ray origin */
    Vector<3,T> tvec=origin-vert0;

    /* calculate U parameter and test bounds */
    T u=dot(tvec,pvec)*inv_det;
    if(u<T(0) || u>T(1))
    {
#if FE_RTI_DEBUG
        feLog(" miss: u=%.6G\n",u);
#endif
        return T(-1);
    }

    /* prepare to test V parameter */
    Vector<3,T> qvec;
    cross3(qvec,tvec,edge1);

    /* calculate V parameter and test bounds */
    T v=dot(direction,qvec)*inv_det;
    if(v<T(0) || u+v>T(1))
    {
#if FE_RTI_DEBUG
        feLog(" miss: u=%.6G v=%.6G\n",u,v);
#endif
        return T(-1);
    }

    /* calculate t, ray intersects triangle */
    T range=dot(edge2,qvec)*inv_det;
#endif
/*****************************************************************************/

    barycenter.setUV(u,v);

#if FE_RTI_DEBUG
        feLog(" hit: u=%.6G v=%.6G range=%.6G\n",u,v,range);
#endif

    return range;
}

// direction must be normalized
template <typename T>
inline void RayTriangleIntersect<T>::resolveContact(
        const Vector<3,T>& vert0,
        const Vector<3,T>& vert1,
        const Vector<3,T>& vert2,
        const Vector<3,T>& norm0,
        const Vector<3,T>& norm1,
        const Vector<3,T>& norm2,
        const Vector<3,T>& origin,
        const Vector<3,T>& direction,
        const Barycenter<T>& barycenter,
        const T range,
        Vector<3,T>& intersection,
        Vector<3,T>& normal)
{
    intersection=origin+direction*range;

#if FALSE
    Vector<3,T> edge1=vert1-vert0;
    Vector<3,T> edge2=vert2-vert0;
    cross(normal,edge1,edge2);
#else
    normal=location(barycenter,norm0,norm1,norm2);
#endif
    normalizeSafe(normal);
}

} /* namespace ext */
} /* namespace fe */

#endif /* __solve_RayTriangleIntersect_h__ */