A Common Lisp port of the Möller–Trumbore ray-triangle intersection algorithm originally authored in C.

ray-triangle-intersect.lisp

;;; A Common Lisp port of the Möller–Trumbore ray-triangle intersection
;;; algorithm originally authored in C.

;;; Introduction
;;; ============

;;; The Möller–Trumbore ray-triangle intersection algorithm is credited to its
;;; inventors Tomas Möller and Ben Trumbore.

;;; The original C code was accessed from Möller Trumbore's article: `Practical
;;; Analysis of Optimized Ray-Triangle Intersection`. Link:
;;; https://fileadmin.cs.lth.se/cs/Personal/Tomas_Akenine-Moller/raytri/

;; We extend our gratitude and thanks to Tomas Möller, not only for his
;; insightful article but also for making the source code publicly accessible
;; under a permissive license.


;;; Common Lisp Port
;;; ================

;;; `intersect-triangle` below, is a port of `intersect_triangle` from raytri.c
;;; including the original copyright messages that were found.




;; Ray-Triangle Intersection Test Routines
;; Different optimizations of my and Ben Trumbore's
;; code from journals of graphics tools (JGT)
;; http://www.acm.org/jgt/
;; by Tomas Moller, May 2000

;; Copyright 2020 Tomas Akenine-Möller

;; Permission is hereby granted, free of charge, to any person obtaining a copy
;; of this software and associated documentation files (the "Software"), to deal
;; in the Software without restriction, including without limitation the rights
;; to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
;; copies of the Software, and to permit persons to whom the Software is
;; furnished to do so, subject to the following conditions:

;; The above copyright notice and this permission notice shall be included in
;; all copies or substantial portions of the Software.

;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
;; IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
;; FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
;; AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
;; LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
;; OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
;; SOFTWARE.

(defun intersect-triangle (tri ray)
  (let ((epsilon 0.000001))
    (let ((orig (origin.geometry.ray:origin ray))
          (dir (origin.geometry.ray:direction ray))
          (vert0 (origin.geometry.triangle::a tri))
          (vert1 (origin.geometry.triangle::b tri))
          (vert2 (origin.geometry.triangle::c tri)))
      (let* (
             ;; find vectors for two edges sharing vert0
             (edge1 (p:- vert1 vert0))
             (edge2 (p:- vert2 vert0))
             ;; begin calculating determinant - also used to calculate U
             (pvec (p:cross dir edge2))
             ;; if determinant is near zero, ray lies in plane of triangle
             (det (p:dot edge1 pvec)))
        (when (and (> det (- epsilon)) (< det epsilon))
          (return-from intersect-triangle nil))
        (let* ((inv-det (/ 1.0 det))
               ;; calculate distance from vert0 to ray origin
               (tvec (p:- orig vert0))
               ;; calculate U and test bounds
               (u (* (p:dot tvec pvec) inv-det)))
          (when (or (< u 0.0) (> u 1.0))
            (return-from intersect-triangle nil))
          (let* (;;prepare to test V
                 (qvec (p:cross tvec edge1))
                 ;; calculate V parameter and test bounds
                 (v (* (p:dot dir qvec) inv-det)))
            (when (or (< v 0.0) (> (+ u v) 1.0))
              (return-from intersect-triangle nil))
            (let (;; calculate t, ray intersects triangle
                  ;; (using `t_` since `t` clashes in Common Lisp)
                  (t_ (* (p:dot edge2 qvec) inv-det)))
              ;; The original C code returns 1 if intersection occurs and 0
              ;; otherwise, while the actual values computed by the function,
              ;; namely t, u, v are returned indirectly via pointer arguments
              ;; whoose contents are populated by the function.

              ;; Instead, we return `(values t u v)` if intersection occurs and
              ;; nil otherwise.
              (values t_ u v))))))))


;;; The original C code - raytri.c
;;; ==============================

;;; For reference and safekeeping, below is copy of the original author's source
;;; code - raytri.c - taken from:
;;; https://fileadmin.cs.lth.se/cs/Personal/Tomas_Akenine-Moller/raytri/raytri.c

;;; Also see Tomas Möller's insightful article:
;;; https://fileadmin.cs.lth.se/cs/Personal/Tomas_Akenine-Moller/raytri/

;; Again, We extend our gratitude and thanks to Tomas Möller, not only for his
;; insightful article but also for making the source code publicly accessible
;; under a permissive license.




;;; /* Ray-Triangle Intersection Test Routines          */
;;; /* Different optimizations of my and Ben Trumbore's */
;;; /* code from journals of graphics tools (JGT)       */
;;; /* http://www.acm.org/jgt/                          */
;;; /* by Tomas Moller, May 2000                        */
;;; 
;;; /*
;;; Copyright 2020 Tomas Akenine-Möller
;;; 
;;; Permission is hereby granted, free of charge, to any person obtaining a copy
;;; of this software and associated documentation files (the "Software"), to
;;; deal in the Software without restriction, including without limitation the
;;; rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
;;; sell copies of the Software, and to permit persons to whom the Software is
;;; furnished to do so, subject to the following conditions:

;;; The above copyright notice and this permission notice shall be included in
;;; all copies or substantial portions of the Software.

;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
;;; IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
;;; FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
;;; AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
;;; LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
;;; IN THE SOFTWARE.
;;; */
;;
;; #include <math.h>
;;
;; #define EPSILON 0.000001
;; #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];
;;
;; /* the original jgt code */
;; int intersect_triangle(double orig[3], double dir[3],
;; 		       double vert0[3], double vert1[3], double vert2[3],
;; 		       double *t, double *u, double *v)
;; {
;;    double edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
;;    double det,inv_det;
;;
;;    /* find vectors for two edges sharing vert0 */
;;    SUB(edge1, vert1, vert0);
;;    SUB(edge2, vert2, vert0);
;;
;;    /* begin calculating determinant - also used to calculate U parameter */
;;    CROSS(pvec, dir, edge2);
;;
;;    /* if determinant is near zero, ray lies in plane of triangle */
;;    det = DOT(edge1, pvec);
;;
;;    if (det > -EPSILON && det < EPSILON)
;;      return 0;
;;    inv_det = 1.0 / det;
;;
;;    /* calculate distance from vert0 to ray origin */
;;    SUB(tvec, orig, vert0);
;;
;;    /* calculate U parameter and test bounds */
;;    *u = DOT(tvec, pvec) * inv_det;
;;    if (*u < 0.0 || *u > 1.0)
;;      return 0;
;;
;;    /* prepare to test V parameter */
;;    CROSS(qvec, tvec, edge1);
;;
;;    /* calculate V parameter and test bounds */
;;    *v = DOT(dir, qvec) * inv_det;
;;    if (*v < 0.0 || *u + *v > 1.0)
;;      return 0;
;;
;;    /* calculate t, ray intersects triangle */
;;    *t = DOT(edge2, qvec) * inv_det;
;;
;;    return 1;
;; }
;;
;;
;; /* code rewritten to do tests on the sign of the determinant */
;; /* the division is at the end in the code                    */
;; int intersect_triangle1(double orig[3], double dir[3],
;; 			double vert0[3], double vert1[3], double vert2[3],
;; 			double *t, double *u, double *v)
;; {
;;    double edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
;;    double det,inv_det;
;;
;;    /* find vectors for two edges sharing vert0 */
;;    SUB(edge1, vert1, vert0);
;;    SUB(edge2, vert2, vert0);
;;
;;    /* begin calculating determinant - also used to calculate U parameter */
;;    CROSS(pvec, dir, edge2);
;;
;;    /* if determinant is near zero, ray lies in plane of triangle */
;;    det = DOT(edge1, pvec);
;;
;;    if (det > EPSILON)
;;    {
;;       /* calculate distance from vert0 to ray origin */
;;       SUB(tvec, orig, vert0);
;;
;;       /* calculate U parameter and test bounds */
;;       *u = DOT(tvec, pvec);
;;       if (*u < 0.0 || *u > det)
;; 	 return 0;
;;
;;       /* prepare to test V parameter */
;;       CROSS(qvec, tvec, edge1);
;;
;;       /* calculate V parameter and test bounds */
;;       *v = DOT(dir, qvec);
;;       if (*v < 0.0 || *u + *v > det)
;; 	 return 0;
;;
;;    }
;;    else if(det < -EPSILON)
;;    {
;;       /* calculate distance from vert0 to ray origin */
;;       SUB(tvec, orig, vert0);
;;
;;       /* calculate U parameter and test bounds */
;;       *u = DOT(tvec, pvec);
;; /*      printf("*u=%f\n",(float)*u); */
;; /*      printf("det=%f\n",det); */
;;       if (*u > 0.0 || *u < det)
;; 	 return 0;
;;
;;       /* prepare to test V parameter */
;;       CROSS(qvec, tvec, edge1);
;;
;;       /* calculate V parameter and test bounds */
;;       *v = DOT(dir, qvec) ;
;;       if (*v > 0.0 || *u + *v < det)
;; 	 return 0;
;;    }
;;    else return 0;  /* ray is parallell to the plane of the triangle */
;;
;;
;;    inv_det = 1.0 / det;
;;
;;    /* calculate t, ray intersects triangle */
;;    *t = DOT(edge2, qvec) * inv_det;
;;    (*u) *= inv_det;
;;    (*v) *= inv_det;
;;
;;    return 1;
;; }
;;
;; /* code rewritten to do tests on the sign of the determinant */
;; /* the division is before the test of the sign of the det    */
;; int intersect_triangle2(double orig[3], double dir[3],
;; 			double vert0[3], double vert1[3], double vert2[3],
;; 			double *t, double *u, double *v)
;; {
;;    double edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
;;    double det,inv_det;
;;
;;    /* find vectors for two edges sharing vert0 */
;;    SUB(edge1, vert1, vert0);
;;    SUB(edge2, vert2, vert0);
;;
;;    /* begin calculating determinant - also used to calculate U parameter */
;;    CROSS(pvec, dir, edge2);
;;
;;    /* if determinant is near zero, ray lies in plane of triangle */
;;    det = DOT(edge1, pvec);
;;
;;    /* calculate distance from vert0 to ray origin */
;;    SUB(tvec, orig, vert0);
;;    inv_det = 1.0 / det;
;;
;;    if (det > EPSILON)
;;    {
;;       /* calculate U parameter and test bounds */
;;       *u = DOT(tvec, pvec);
;;       if (*u < 0.0 || *u > det)
;; 	 return 0;
;;
;;       /* prepare to test V parameter */
;;       CROSS(qvec, tvec, edge1);
;;
;;       /* calculate V parameter and test bounds */
;;       *v = DOT(dir, qvec);
;;       if (*v < 0.0 || *u + *v > det)
;; 	 return 0;
;;
;;    }
;;    else if(det < -EPSILON)
;;    {
;;       /* calculate U parameter and test bounds */
;;       *u = DOT(tvec, pvec);
;;       if (*u > 0.0 || *u < det)
;; 	 return 0;
;;
;;       /* prepare to test V parameter */
;;       CROSS(qvec, tvec, edge1);
;;
;;       /* calculate V parameter and test bounds */
;;       *v = DOT(dir, qvec) ;
;;       if (*v > 0.0 || *u + *v < det)
;; 	 return 0;
;;    }
;;    else return 0;  /* ray is parallell to the plane of the triangle */
;;
;;    /* calculate t, ray intersects triangle */
;;    *t = DOT(edge2, qvec) * inv_det;
;;    (*u) *= inv_det;
;;    (*v) *= inv_det;
;;
;;    return 1;
;; }
;;
;; /* code rewritten to do tests on the sign of the determinant */
;; /* the division is before the test of the sign of the det    */
;; /* and one CROSS has been moved out from the if-else if-else */
;; int intersect_triangle3(double orig[3], double dir[3],
;; 			double vert0[3], double vert1[3], double vert2[3],
;; 			double *t, double *u, double *v)
;; {
;;    double edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
;;    double det,inv_det;
;;
;;    /* find vectors for two edges sharing vert0 */
;;    SUB(edge1, vert1, vert0);
;;    SUB(edge2, vert2, vert0);
;;
;;    /* begin calculating determinant - also used to calculate U parameter */
;;    CROSS(pvec, dir, edge2);
;;
;;    /* if determinant is near zero, ray lies in plane of triangle */
;;    det = DOT(edge1, pvec);
;;
;;    /* calculate distance from vert0 to ray origin */
;;    SUB(tvec, orig, vert0);
;;    inv_det = 1.0 / det;
;;
;;    CROSS(qvec, tvec, edge1);
;;
;;    if (det > EPSILON)
;;    {
;;       *u = DOT(tvec, pvec);
;;       if (*u < 0.0 || *u > det)
;; 	 return 0;
;;
;;       /* calculate V parameter and test bounds */
;;       *v = DOT(dir, qvec);
;;       if (*v < 0.0 || *u + *v > det)
;; 	 return 0;
;;
;;    }
;;    else if(det < -EPSILON)
;;    {
;;       /* calculate U parameter and test bounds */
;;       *u = DOT(tvec, pvec);
;;       if (*u > 0.0 || *u < det)
;; 	 return 0;
;;
;;       /* calculate V parameter and test bounds */
;;       *v = DOT(dir, qvec) ;
;;       if (*v > 0.0 || *u + *v < det)
;; 	 return 0;
;;    }
;;    else return 0;  /* ray is parallell to the plane of the triangle */
;;
;;    *t = DOT(edge2, qvec) * inv_det;
;;    (*u) *= inv_det;
;;    (*v) *= inv_det;
;;
;;    return 1;
;; }