/ray/src/lib/numerics/3d/triangle.cc

Go to the documentation of this file.

#include <lib/numerics/3d/vector3.h>


namespace NUM  // numerics
{


template<typename T,typename V>int IntersectTriangle(
    const Vector3<T> &orig,const Vector3<T> &dir,
    const Vector3<V> &vert0,const Vector3<V> &vert1,const Vector3<V> &vert2,
    T *t,T *u,T *v)
{
static const V EPSILON=0.000001;

    /* find vectors for two edges sharing vert0 */
    Vector3<V> edge1=vert1-vert0;
    Vector3<V> edge2=vert2-vert0;
    
    /* begin calculating determinant - also used to calculate U parameter */
    Vector3<V> pvec=cross(dir,edge2);
    
    /* if determinant is near zero, ray lies in plane of triangle */
    V det = dot(edge1,pvec);
    
#ifdef TEST_CULL           /* define TEST_CULL if culling is desired */
    if(det < EPSILON)  return(0);
    
    /* calculate distance from vert0 to ray origin */
    Vector3<T> tvec=orig-vert0;
    
    /* calculate U parameter and test bounds */
    *u = dot<T>(tvec,pvec);
    if(*u < 0 || *u > det)  return(0);
    
    /* prepare to test V parameter */
    Vector3<V> qvec=cross(tvec,edge1);
    
    /* calculate V parameter and test bounds */
    *v = dot<T>(dir,qvec);
    if(*v < 0 || *u + *v > det)  return(0);
    
    /* calculate t, scale parameters, ray intersects triangle */
    *t = dot<T>(edge2,qvec);
    T inv_det = T(1)/det;
    *t *= inv_det;
    *u *= inv_det;
    *v *= inv_det;
#else                    /* the non-culling branch */
    if(fabs(det) < EPSILON)  return(0);
    
    T inv_det = T(1)/det;
    
    /* calculate distance from vert0 to ray origin */
    Vector3<T> tvec=orig-vert0;
    
    /* calculate U parameter and test bounds */
    *u = dot<T>(tvec,pvec) * inv_det;
    if(*u < 0 || *u > 1)  return(0);
    
    /* prepare to test V parameter */
    Vector3<V> qvec=cross(tvec,edge1);
    
    /* calculate V parameter and test bounds */
    *v = DOT<T>(dir,qvec) * inv_det;
    if(*v < 0 || *u + *v > 1)  return(0);
    
    /* calculate t, ray intersects triangle */
    *t = dot<T>(edge2,qvec) * inv_det;
#endif
    
    return(1);
}


#if 0

template<typename T,typename V>int IntersectTriangle(
    const Vector3<T> &orig,const Vector3<T> &dir,
    const Vector3<V> &vert0,const Vector3<V> &vert1,const Vector3<V> &vert2,
    const Plane3d<V> &planeq,int i1,int i2,
    T *t,T *u,T *v)
{
static const T EPSILON=0.000001;

    /* is ray parallel to plane? */
    T dotval = dot<T>(planeq.normal(),dir);
    if(fabs(dotval) < EPSILON)      /* or use culling check */
        return(0);
    
    /* find distance to plane and intersection point */
    // T dot2 = dot<T>(planeq.normal(),orig);
    // *t = -( planeq.dist() + dot2 ) / dotval;
    *t = - distance(orig) / dotval;
    
    T point[2];
    point[0] = orig[i1] + dir[i1] * *t;
    point[1] = orig[i2] + dir[i2] * *t;
    
    /* begin barycentric intersection algorithm */
    T u0 = point[0] - vert0[i1];
    T v0 = point[1] - vert0[i2];
    V u1 = vert1[i1] - vert0[i1];
    V u2 = vert2[i1] - vert0[i1];
    V v1 = vert1[i2] - vert0[i2];
    V v2 = vert2[i2] - vert0[i2];

    /* calculate and compare barycentric coordinates */
    if(u1 == 0)
    {           /* uncommon case */
        *v = u0 / u2;
        if(*v < 0 || *v > 1)  return(0);
        *u = (v0 - *v * v2) / v1;
    }
    else
    {           /* common case, used for this analysis */
        *v = (v0 * u1 - u0 * v1) / (v2 * u1 - u2 * v1);
        if(*v < 0 || *v > 1)  return(0);
        *u = (u0 - *v * u2) / u1;
    }
    
    if (*u < 0 || *u + *v > 1)  return(0);
    
    return(1);
}
#endif

}  // end of namespace NUM

---