Home > Source code of ray collision >

An essay in 3d computation, ray collision

Fun work with Elric, using its OpenGL library we worked together on a collision detection algorithm in order to compute as fast as it would be possible a 3d lighting map. His work was inspired and motivated by the photo-realism rendering of radiosity algorithms.

Download raycol.c

raycol.c

Read raycol.c

  1. /*
  2.  * extract of the Silence OpenGL project
  3.  * (sglMath.cpp part)
  4.  * v1.0
  5.  *
  6.  * 02/05/02
  7.  *
  8.  * by SINTES 'elric' Arnaud / PESCE 'joke' François
  9.  * elric(at)silentbreed.com / joke(at)r00tworld.com
  10.  *
  11.  */
  12.  
  13. /*
  14.  * During the development, we (elric+joke) have to found the
  15.  * _fastest_ algo for Ray/Triangle collision.
  16.  *
  17.  * here is our result.
  18.  *
  19.  */
  20.  
  21. /*
  22.  * NOTES:
  23.  *
  24.  * assumes that it's basic C/OpenGL/SGL syntax.
  25.  *
  26.  * GLfloat is standard C float
  27.  * GLBoolean is standard C bool
  28.  * GL_TRUE is true and GL_FALSE is false
  29.  *
  30.  * SGLpoint3f is a typedef struct covered with 3 GLfloat x, y and z
  31.  *
  32.  * SGL_EPSILUM is equal to 0.0005f
  33.  *
  34.  */
  35.  
  36. /*
  37.  * the sglRayCol function parameters are:
  38.  *
  39.  * p1, p2 : start and end of the ray line
  40.  * pa, pb, pc : points of the triangle to test intersection with
  41.  * function return false if no collision detected,
  42.  * else true and p is set to the collision point
  43.  *
  44.  */
  45.  
  46.  
  47. /*
  48.  * min/max funcs for imediate rejection
  49.  */
  50. static inline GLfloat min2(GLfloat x1, GLfloat x2)
  51. {
  52.     return (x1 <= x2) ? x1 : x2;
  53. }
  54.  
  55. static inline GLfloat max2(GLfloat x1, GLfloat x2)
  56. {
  57.     return (x1 >= x2) ? x1 : x2;
  58. }
  59.  
  60. static inline GLfloat min3(GLfloat x1, GLfloat x2, GLfloat x3)
  61. {
  62.     return (x1 <= x2) ? ((x1 <= x3) ? x1 : x3) : ((x2 <= x3) ? x2 : x3);
  63. }
  64.  
  65. static inline GLfloat max3(GLfloat x1, GLfloat x2, GLfloat x3)
  66. {
  67.     return (x1 >= x2) ? ((x1 >= x3) ? x1 : x3) : ((x2 >= x3) ? x2 : x3);
  68. }
  69.  
  70.  
  71.  
  72. static inline GLboolean
  73. sglImmediateReject(const GLfloat p1x, const GLfloat p1y, const GLfloat p1z,
  74.                    const GLfloat p2x, const GLfloat p2y, const GLfloat p2z,
  75.                    const GLfloat pax, const GLfloat pay, const GLfloat paz,
  76.                    const GLfloat pbx, const GLfloat pby, const GLfloat pbz,
  77.                    const GLfloat pcx, const GLfloat pcy, const GLfloat pcz)
  78. {
  79.     /*
  80.      * According to "Introduction to algorithmic" by Cormen Rivest Leiserson :
  81.      *
  82.      * > 2 step :
  83.      * >     1. Immediate reject : the cube containing the segment must intersect the polygon.
  84.      * >     lp stand for limit point.
  85.      * >     Cube1 = (lp1, lp2);
  86.      * >     Cube2 = (lp3, lp4);
  87.      */
  88.     /* limit boxes */
  89.     GLfloat lpM, lpm;
  90.     /* lp4z lp1z */
  91.     lpM = max3(paz, pbz, pcz);
  92.     lpm = min2(p1z, p2z);
  93.     if (lpM < lpm)
  94.         return GL_TRUE;
  95.     /* lp2z lp3z */
  96.     lpM = max2(p1z, p2z);
  97.     lpm = min3(paz, pbz, pcz);
  98.     if (lpM < lpm)
  99.         return GL_TRUE;
  100.     /* lp2x lp3x */
  101.     lpM = max2(p1x, p2x);
  102.     lpm = min3(pax, pbx, pcx);
  103.     if (lpM < lpm)
  104.         return GL_TRUE;
  105.     /* lp4x lp1x */
  106.     lpM = max3(pax, pbx, pcx);
  107.     lpm = min2(p1x, p2x);
  108.     if (lpM < lpm)
  109.         return GL_TRUE;
  110.     /* lp4y lp1y */
  111.     lpM = max3(pay, pby, pcy);
  112.     lpm = min2(p1y, p2y);
  113.     if (lpM < lpm)
  114.         return GL_TRUE;
  115.     /* lp2y lp3y */
  116.     lpM = max2(p1y, p2y);
  117.     lpm = min3(pay, pby, pcy);
  118.     if (lpM < lpm)
  119.         return GL_TRUE;
  120.     return GL_FALSE;
  121. }
  122.  
  123.  
  124.  
  125. GLboolean
  126. sglRayCol(const SGLpoint3f p1, const SGLpoint3f p2, const SGLpoint3f pa,
  127.           const SGLpoint3f pb, const SGLpoint3f pc, SGLpoint3f * p)
  128. {
  129.     /*
  130.      * Determine whether or not the line segment p1,p2
  131.      * Intersects the 3 vertex facet bounded by pa,pb,pc
  132.      *
  133.      * Return SGL_NULL if no collison or the intersection
  134.      * point p else
  135.      *
  136.      *
  137.      * The equation of the line is p = p1 + mu (p2 - p1)
  138.      * The equation of the plane is a x + b y + c z + d = 0
  139.      *                        n.x x + n.y y + n.z z + d = 0
  140.      *
  141.      * this version is heavily inspired from
  142.      * softSurfer (www.softsurfer.com)
  143.      * and hardly optimised.
  144.      */
  145.     GLfloat p1x, p1y, p1z;
  146.     GLfloat p2x, p2y, p2z;
  147.     GLfloat pax, pay, paz;
  148.     GLfloat pbx, pby, pbz;
  149.     GLfloat pcx, pcy, pcz;
  150.     GLfloat uu, uv, uw;
  151.     GLfloat vu, vv, vw;
  152.     GLfloat nu, nv, nw;
  153.     GLfloat diru, dirv, dirw;
  154.     GLfloat w0u, w0v, w0w;
  155.     GLfloat a, b, r, D, s, t;
  156.     GLfloat uud, uvd, vvd;
  157.     GLfloat wu, wv, ww;
  158.     GLfloat wud, wvd;
  159.     p1x = p1.x;
  160.     p1y = p1.y;
  161.     p1z = p1.z;
  162.     p2x = p2.x;
  163.     p2y = p2.y;
  164.     p2z = p2.z;
  165.     pax = pa.x;
  166.     pay = pa.y;
  167.     paz = pa.z;
  168.     pbx = pb.x;
  169.     pby = pb.y;
  170.     pbz = pb.z;
  171.     pcx = pc.x;
  172.     pcy = pc.y;
  173.     pcz = pc.z;
  174.     /* imediate reject */
  175.     if (sglImmediateReject(p1x, p1y, p1z, p2x, p2y, p2z, pax, pay, paz, pbx, pby, pbz, pcx, pcy, pcz))
  176.         return GL_FALSE;
  177.     /* get triangle edge vectors and plane normal */
  178.     uu = pbx - pax;
  179.     uv = pby - pay;
  180.     uw = pbz - paz;
  181.     vu = pcx - pax;
  182.     vv = pcy - pay;
  183.     vw = pcz - paz;
  184.     nu = uv * vw - uw * vv;
  185.     nv = uw * vu - uu * vw;
  186.     nw = uu * vv - uv * vu;
  187.     /* ray direction vector */
  188.     diru = p2x - p1x;
  189.     dirv = p2y - p1y;
  190.     dirw = p2z - p1z;
  191.     b = nu * diru + nv * dirv + nw * dirw;
  192.     if (b < SGL_EPSILUM)
  193.         /* reduction of  rayCol with Culled triangles
  194.          * consider ray parallel to triangle plane and
  195.          * ray lies in triangle plane refer to cull test
  196.          */
  197.         return GL_FALSE;
  198.     /* old code:
  199.      *
  200.      * if( b < 0 )          b = -b;         // invert culling
  201.      * if( b < SGL_EPSILUM) {               // ray is parallel to triangle plane
  202.      *              return GL_FALSE;
  203.      *      if( a == 0)                     // ray lies in triangle plane
  204.      *              return GL_TRUE;         // in fact, return another value
  205.      *      else
  206.      *              return GL_FALSE;        // ray disjoint from plane
  207.      * }
  208.      */
  209.     /* get intersect point of ray with triangle plane */
  210.     w0u = p1x - pax;
  211.     w0v = p1y - pay;
  212.     w0w = p1z - paz;
  213.     a = nu * w0u + nv * w0v + nw * w0w;
  214.     r = -a / b;
  215.     if (r < 0 || r > 1)         /* ray goes away from triangle */
  216.         return GL_FALSE;        /* => no intersect */
  217.     /* intersect point of ray and plane */
  218.     p->x = p1x + r * diru;
  219.     p->y = p1y + r * dirv;
  220.     p->z = p1z + r * dirw;
  221.     /* is p inside the triangle ? */
  222.     uud = uu * uu + uv * uv + uw * uw;
  223.     uvd = uu * vu + uv * vv + uw * vw;
  224.     vvd = vu * vu + vv * vv + vw * vw;
  225.     wu = p->x - pax;
  226.     wv = p->y - pay;
  227.     ww = p->z - paz;
  228.     wud = wu * uu + wv * uv + ww * uw;
  229.     wvd = wu * vu + wv * vv + ww * vw;
  230.     D = 1 / (uvd * uvd - uud * vvd);
  231.     /* get and test parametric coords */
  232.     s = (uvd * wvd - vvd * wud) * D;
  233.     if (s < 0 || s > 1)         /* p is outside triangle */
  234.         return GL_FALSE;
  235.     t = (uvd * wud - uud * wvd) * D;
  236.     if (t < 0 || (s + t) > 1)   /* p is outside triangle */
  237.         return GL_FALSE;
  238.     return GL_TRUE;             /* p is inside triangle */
  239. }
  240.