raytri.cpp

Go to the documentation of this file.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include <math.h>
 
#include "raytri.h"
 
namespace ConvexDecomposition
{
 
 
/* a = b - c */
#define vector(a,b,c) \
        (a)[0] = (b)[0] - (c)[0];       \
        (a)[1] = (b)[1] - (c)[1];       \
        (a)[2] = (b)[2] - (c)[2];
 
 
 
#define innerProduct(v,q) \
                ((v)[0] * (q)[0] + \
                (v)[1] * (q)[1] + \
                (v)[2] * (q)[2])
 
#define crossProduct(a,b,c) \
        (a)[0] = (b)[1] * (c)[2] - (c)[1] * (b)[2]; \
        (a)[1] = (b)[2] * (c)[0] - (c)[2] * (b)[0]; \
        (a)[2] = (b)[0] * (c)[1] - (c)[0] * (b)[1];
 
bool rayIntersectsTriangle(const double *p,const double *d,const double *v0,const double *v1,const double *v2,double &t)
{
 
        double e1[3],e2[3],h[3],s[3],q[3];
        double a,f,u,v;
 
        vector(e1,v1,v0);
        vector(e2,v2,v0);
        crossProduct(h,d,e2);
        a = innerProduct(e1,h);
 
        if (a > -0.00001 && a < 0.00001)
                return(false);
 
        f = 1/a;
        vector(s,p,v0);
        u = f * (innerProduct(s,h));
 
        if (u < 0.0 || u > 1.0)
                return(false);
 
        crossProduct(q,s,e1);
        v = f * innerProduct(d,q);
        if (v < 0.0 || u + v > 1.0)
                return(false);
        // at this stage we can compute t to find out where
        // the intersection point is on the line
        t = f * innerProduct(e2,q);
        if (t > 0) // ray intersection
                return(true);
        else // this means that there is a line intersection
                 // but not a ray intersection
                 return (false);
}
 
 
bool lineIntersectsTriangle(const double *rayStart,const double *rayEnd,const double *p1,const double *p2,const double *p3,double *sect)
{
        double dir[3];
 
  dir[0] = rayEnd[0] - rayStart[0];
  dir[1] = rayEnd[1] - rayStart[1];
  dir[2] = rayEnd[2] - rayStart[2];
 
  double d = sqrt(dir[0]*dir[0] + dir[1]*dir[1] + dir[2]*dir[2]);
  double r = 1.0f / d;
 
  dir[0]*=r;
  dir[1]*=r;
  dir[2]*=r;
 
 
  double t;
 
        bool ret = rayIntersectsTriangle(rayStart, dir, p1, p2, p3, t );
 
        if ( ret )
        {
                if ( t > d )
                {
                        sect[0] = rayStart[0] + dir[0]*t;
                        sect[1] = rayStart[1] + dir[1]*t;
                        sect[2] = rayStart[2] + dir[2]*t;
                }
                else
                {
                        ret = false;
                }
        }
 
  return ret;
}
 
}; // end of namespace