tribox3.h

Go to the documentation of this file.

/** @file tribox3.h
 */
/********************************************************/
/* AABB-triangle overlap test code                      */
/* by Tomas Akenine-Möller                              */
/* Function: int triBoxOverlap(float boxcenter[3],      */
/*          float boxhalfsize[3],float triverts[3][3]); */
/* History:                                             */
/*   2001-03-05: released the code in its first version */
/*   2001-06-18: changed the order of the tests, faster */
/*                                                      */
/* Acknowledgement: Many thanks to Pierre Terdiman for  */
/* suggestions and discussions on how to optimize code. */
/* Thanks to David Hunt for finding a ">="-bug!         */
/* See also: https://doi.org/10.1145/1198555.1198747    */
/* http://fileadmin.cs.lth.se/cs/personal/tomas_akenine-moller/code/ */
/********************************************************/
#include <math.h>
#include <stdio.h>
 
#define X 0
#define Y 1
#define Z 2
 
#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]; 
 
#define FINDMINMAX(x0,x1,x2,min,max) \
  min = max = x0;   \
  if(x1<min) min=x1;\
  if(x1>max) max=x1;\
  if(x2<min) min=x2;\
  if(x2>max) max=x2;
 
int planeBoxOverlap(float normal[3], float vert[3], float maxbox[3])  // -NJMP-
{
  int q;
  float vmin[3],vmax[3],v;
  for(q=X;q<=Z;q++)
  {
    v=vert[q];          // -NJMP-
    if(normal[q]>0.0f)
    {
      vmin[q]=-maxbox[q] - v; // -NJMP-
      vmax[q]= maxbox[q] - v; // -NJMP-
    }
    else
    {
      vmin[q]= maxbox[q] - v; // -NJMP-
      vmax[q]=-maxbox[q] - v; // -NJMP-
    }
  }
  if(DOT(normal,vmin)>0.0f) return 0; // -NJMP-
  if(DOT(normal,vmax)>=0.0f) return 1;  // -NJMP-
  
  return 0;
}
 
 
/*======================== X-tests ========================*/
#define AXISTEST_X01(a, b, fa, fb)         \
  p0 = a*v0[Y] - b*v0[Z];                \
  p2 = a*v2[Y] - b*v2[Z];                \
        if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;} \
  rad = fa * boxhalfsize[Y] + fb * boxhalfsize[Z];   \
  if(min>rad || max<-rad) return 0;
 
#define AXISTEST_X2(a, b, fa, fb)        \
  p0 = a*v0[Y] - b*v0[Z];                \
  p1 = a*v1[Y] - b*v1[Z];                \
        if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;} \
  rad = fa * boxhalfsize[Y] + fb * boxhalfsize[Z];   \
  if(min>rad || max<-rad) return 0;
 
/*======================== Y-tests ========================*/
#define AXISTEST_Y02(a, b, fa, fb)         \
  p0 = -a*v0[X] + b*v0[Z];               \
  p2 = -a*v2[X] + b*v2[Z];                     \
        if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;} \
  rad = fa * boxhalfsize[X] + fb * boxhalfsize[Z];   \
  if(min>rad || max<-rad) return 0;
 
#define AXISTEST_Y1(a, b, fa, fb)        \
  p0 = -a*v0[X] + b*v0[Z];               \
  p1 = -a*v1[X] + b*v1[Z];                   \
        if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;} \
  rad = fa * boxhalfsize[X] + fb * boxhalfsize[Z];   \
  if(min>rad || max<-rad) return 0;
 
/*======================== Z-tests ========================*/
 
#define AXISTEST_Z12(a, b, fa, fb)         \
  p1 = a*v1[X] - b*v1[Y];                \
  p2 = a*v2[X] - b*v2[Y];                \
        if(p2<p1) {min=p2; max=p1;} else {min=p1; max=p2;} \
  rad = fa * boxhalfsize[X] + fb * boxhalfsize[Y];   \
  if(min>rad || max<-rad) return 0;
 
#define AXISTEST_Z0(a, b, fa, fb)        \
  p0 = a*v0[X] - b*v0[Y];          \
  p1 = a*v1[X] - b*v1[Y];                \
        if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;} \
  rad = fa * boxhalfsize[X] + fb * boxhalfsize[Y];   \
  if(min>rad || max<-rad) return 0;
 
int triBoxOverlap(float boxcenter[3],float boxhalfsize[3],float triverts[3][3])
{
 
  /*    use separating axis theorem to test overlap between triangle and box */
  /*    need to test for overlap in these directions: */
  /*    1) the {x,y,z}-directions (actually, since we use the AABB of the triangle */
  /*       we do not even need to test these) */
  /*    2) normal of the triangle */
  /*    3) crossproduct(edge from tri, {x,y,z}-directin) */
  /*       this gives 3x3=9 more tests */
   float v0[3],v1[3],v2[3];
//   float axis[3];
   float min,max,p0,p1,p2,rad,fex,fey,fez;    // -NJMP- "d" local variable removed
   float normal[3],e0[3],e1[3],e2[3];
 
   /* This is the fastest branch on Sun */
   /* move everything so that the boxcenter is in (0,0,0) */
   SUB(v0,triverts[0],boxcenter);
   SUB(v1,triverts[1],boxcenter);
   SUB(v2,triverts[2],boxcenter);
 
   /* compute triangle edges */
   SUB(e0,v1,v0);      /* tri edge 0 */
   SUB(e1,v2,v1);      /* tri edge 1 */
   SUB(e2,v0,v2);      /* tri edge 2 */
 
   /* Bullet 3:  */
   /*  test the 9 tests first (this was faster) */
   fex = fabsf(e0[X]);
   fey = fabsf(e0[Y]);
   fez = fabsf(e0[Z]);
   AXISTEST_X01(e0[Z], e0[Y], fez, fey);
   AXISTEST_Y02(e0[Z], e0[X], fez, fex);
   AXISTEST_Z12(e0[Y], e0[X], fey, fex);
 
   fex = fabsf(e1[X]);
   fey = fabsf(e1[Y]);
   fez = fabsf(e1[Z]);
   AXISTEST_X01(e1[Z], e1[Y], fez, fey);
   AXISTEST_Y02(e1[Z], e1[X], fez, fex);
   AXISTEST_Z0(e1[Y], e1[X], fey, fex);
 
   fex = fabsf(e2[X]);
   fey = fabsf(e2[Y]);
   fez = fabsf(e2[Z]);
   AXISTEST_X2(e2[Z], e2[Y], fez, fey);
   AXISTEST_Y1(e2[Z], e2[X], fez, fex);
   AXISTEST_Z12(e2[Y], e2[X], fey, fex);
 
   /* Bullet 1: */
   /*  first test overlap in the {x,y,z}-directions */
   /*  find min, max of the triangle each direction, and test for overlap in */
   /*  that direction -- this is equivalent to testing a minimal AABB around */
   /*  the triangle against the AABB */
 
   /* test in X-direction */
   FINDMINMAX(v0[X],v1[X],v2[X],min,max);
   if(min>boxhalfsize[X] || max<-boxhalfsize[X]) return 0;
 
   /* test in Y-direction */
   FINDMINMAX(v0[Y],v1[Y],v2[Y],min,max);
   if(min>boxhalfsize[Y] || max<-boxhalfsize[Y]) return 0;
 
   /* test in Z-direction */
   FINDMINMAX(v0[Z],v1[Z],v2[Z],min,max);
   if(min>boxhalfsize[Z] || max<-boxhalfsize[Z]) return 0;
 
   /* Bullet 2: */
   /*  test if the box intersects the plane of the triangle */
   /*  compute plane equation of triangle: normal*x+d=0 */
   CROSS(normal,e0,e1);
   // -NJMP- (line removed here)
   if(!planeBoxOverlap(normal,v0,boxhalfsize)) return 0;  // -NJMP-
 
   return 1;   /* box and triangle overlaps */
}
 
bool segBoxOverlap (coord * min, coord * max,
        coord * p1, coord * p2)
{
  // Find min and max X for the segment
 
  double minX = p1->x;
  double maxX = p2->x;
 
  if(p1->x > p2->x)
    {
      minX = p2->x;
      maxX = p1->x;
    }
 
  // Find the intersection of the segment's and rectangle's x-projections
 
  if(maxX > max->x)
    {
      maxX = max->x;
    }
 
  if(minX < min->x)
    {
      minX = min->x;
    }
 
  if(minX > maxX) // If their projections do not intersect return false
    {
      return false;
    }
 
  // Find corresponding min and max Y for min and max X we found before
 
  double minY = p1->y;
  double maxY = p2->y;
 
  double dx = p2->x - p1->x;
 
  if(fabs(dx) > 0.0000001)
    {
      double a = (p2->y - p1->y) / dx;
      double b = p1->y - a * p1->x;
      minY = a * minX + b;
      maxY = a * maxX + b;
    }
 
  if(minY > maxY)
    {
      double tmp = maxY;
      maxY = minY;
      minY = tmp;
    }
 
  // Find the intersection of the segment's and rectangle's y-projections
 
  if(maxY > max->y)
    {
      maxY = max->y;
    }
 
  if(minY < min->y)
    {
      minY = min->y;
    }
 
  if(minY > maxY) // If Y-projections do not intersect return false
    {
      return false;
    }
 
  return true;
}