Algoritmo de intersección de Möller-Trumbore

El rayo se define mediante un punto de origen.${\displaystyle O}$y un vector de dirección${\displaystyle D}$Cada punto en el rayo puede expresarse mediante${\displaystyle {\vec {r}}(t)=O+tD}$donde el parámetro${\displaystyle t}$abarca desde menos infinito hasta infinito. El triángulo está definido por tres vértices, llamados${\displaystyle v_{1}}$,${\displaystyle v_{2}}$,${\displaystyle v_{3}}$. El plano sobre el que se encuentra el triángulo, que es necesario para calcular la intersección rayo-triángulo, se define mediante un punto en el plano, como por ejemplo:${\displaystyle v_{1}}$y un vector que es ortogonal a cada punto de ese plano, como el producto vectorial entre el vector de${\displaystyle v_{1}}$a${\displaystyle v_{2}}$y el vector de${\displaystyle v_{1}}$a${\displaystyle v_{3}}$:

${\displaystyle {\vec {n}}\cdot (P_{1}-P_{2})=0}$, dónde${\displaystyle {\vec {n}}=(v_{2}-v_{1})\times (v_{3}-v_{1})}$, y${\displaystyle P_{1}}$y${\displaystyle P_{2}}$¿Hay algún punto en el avión?

Primero, determine si la línea generada por el rayo interseca el plano sobre el que se encuentra el triángulo y, de ser así, encuentre las coordenadas de dicha intersección. La única forma en que la línea no interseca el plano es si el vector director del rayo es paralelo al plano. ^{[ 3 ]} Cuando esto sucede, el producto escalar entre el vector director del rayo y el vector normal del plano será cero. De lo contrario, la línea interseca el plano en algún punto , pero no necesariamente dentro del triángulo.

Utilizando coordenadas baricéntricas , cualquier punto del triángulo puede expresarse como una combinación convexa de los vértices del triángulo:

${\displaystyle P=wv_{1}+uv_{2}+vv_{3}}$

Los coeficientes deben ser no negativos y sumar 1, por lo tanto${\displaystyle w}$puede ser reemplazado por${\displaystyle 1-uv}$:

${\displaystyle {\begin{aligned}P&=(1-uv)v_{1}+uv_{2}+vv_{3}\\P&=v_{1}+u(v_{2}-v_{1})+v(v_{3}-v_{1})\end{aligned}}}$

dónde${\displaystyle P}$es cualquier punto del plano. Observa que${\displaystyle {\vec {e_{1}}}=v_{2}-v_{1}}$y${\displaystyle {\vec {e_{2}}}=v_{3}-v_{1}}$son vectores en el borde del triángulo y, juntos, abarcan un plano (que pasa por el origen). Cada punto en ese plano se puede escribir como${\displaystyle ue_{1}+ve_{2}}$y puede ser traducido por${\displaystyle v_{1}}$"mover" ese punto al plano en el que se encuentra el triángulo.

Para encontrar${\displaystyle u}$y${\displaystyle v}$Para una intersección particular, iguale la expresión del rayo a la expresión del plano y coloque las variables a un lado y las constantes al otro.

${\displaystyle {\begin{aligned}O+tD&=v_{1}+u(v_{2}-v_{1})+v(v_{3}-v_{1})\\O-v_{1}&=-tD+u(v_{2}-v_{1})+v(v_{3}-v_{1})\end{aligned}}}$

Este es un sistema de ecuaciones lineales con tres ecuaciones (una para cada una)${\displaystyle x}$,${\displaystyle y}$,${\displaystyle z}$) y tres incógnitas (${\displaystyle t}$,${\displaystyle u}$, y${\displaystyle v}$), y puede representarse como una multiplicación matriz-vector.

${\displaystyle {\begin{bmatrix}\vert &\vert &\vert \\-D&(v2-v1)&(v3-v1)\\\vert &\vert &\vert \end{bmatrix}}{\begin{bmatrix}t\\u\\v\end{bmatrix}}=O-v_{1}}$

Esta ecuación siempre tendrá una solución cuando la matriz tenga tres vectores columna linealmente independientes en${\displaystyle \mathbb {R} ^{3}}$y por lo tanto es invertible . Esto sucede si y solo si los vértices del triángulo no son colineales y el rayo no es paralelo al plano.

El algoritmo puede utilizar la regla de Cramer para encontrar la${\displaystyle t}$,${\displaystyle u}$, y${\displaystyle v}$valores para una intersección, y si se encuentra dentro del triángulo, las coordenadas exactas de la intersección se pueden encontrar introduciendo${\displaystyle t}$a la ecuación del rayo.

A continuación se muestra una implementación del algoritmo en C++ utilizando la biblioteca glm:

std :: optional <vec3> ray_intersects_triangle ( const vec3 & ray_origin , const vec3 & ray_vector , const triangle3 & triangle ) { constexpr float epsilon = std :: numeric_limits <float> :: epsilon ( ) ;​vec3 edge1 = triángulo . b - triángulo . a ; vec3 edge2 = triángulo . c - triángulo . a ;// Eliminación de caras posteriores, asumiendo triángulos enrollados en sentido antihorario. const vec3 normal = cross ( edge1 , edge2 ); // No es necesario normalizar if ( dot ( normal , ray_vector ) > 0 ) return {};vec3 ray_cross_e2 = cross ( ray_vector , edge2 ); float det = dot ( edge1 , ray_cross_e2 );if ( abs ( det ) < epsilon ) return {}; // El rayo es paralelo al triángulofloat inv_det = 1.0 / det ; vec3 s = ray_origin - triangle . a ; float u = inv_det * dot ( s , ray_cross_e2 );if ( u < - epsilon || u - 1 > epsilon ) return {}; // El rayo pasa fuera de los límites de edge2vec3 s_cross_e1 = cross ( s , edge1 ); float v = inv_det * dot ( ray_vector , s_cross_e1 );if ( v < - epsilon || u + v - 1 > epsilon ) return {}; // El rayo pasa fuera de los límites de edge1// La línea del rayo interseca con el triángulo. // Calculamos t para encontrar dónde en el rayo está la intersección. float t = inv_det * dot ( edge2 , s_cross_e1 );if ( t > epsilon ) // Intersección de rayos { return vec3 ( ray_origin + ray_vector * t ); } else // Esto significa que hay una intersección de líneas pero no una intersección de rayos. return {}; }

A continuación se muestra una implementación del algoritmo en Rust utilizando la biblioteca glam:

fn moller_trumbore_intersection ( origin : Vec3 , direction : Vec3 , triangle : Triangle ) -> Option < Vec3 > { let e1 = triangle . b - triangle . a ; let e2 = triangle . c - triangle . a ;let ray_cross_e2 = direction.cross ( e2 ) ; let det = e1.dot ( ray_cross_e2 ) ;​​Si det > - f32 :: EPSILON && det < f32 :: EPSILON { return None ; // Este rayo es paralelo a este triángulo. }let inv_det = 1.0 / det ; let s = origin - triangle.a ; let u = inv_det * s.dot ( ray_cross_e2 ) ; if u < 0.0 || u > 1.0 { return None ; }​​​let s_cross_e1 = s . cross ( e1 ); let v = inv_det * direction . dot ( s_cross_e1 ); if v < 0.0 || u + v > 1.0 { return None ; } // En esta etapa podemos calcular t para averiguar dónde está el punto de intersección en la línea. let t = inv_det * e2 . dot ( s_cross_e1 );if t > f32 :: EPSILON { // Intersección de rayos let intersection_point = origin + direction * t ; return Some ( intersection_point ); } else { // Esto significa que hay una intersección de líneas pero no una intersección de rayos. return None ; } }

A continuación se muestra una implementación del algoritmo en Java utilizando javax.vecmathla API de Java 3D :

clase pública MollerTrumbore {privado estático final doble EPSILON = 0.0000001 ;public static boolean rayIntersectsTriangle ( Point3d rayOrigin , Vector3d rayVector , Triangle inTriangle , Point3d outIntersectionPoint ) { Point3d vertex0 = inTriangle . getVertex0 (); Point3d vertex1 = inTriangle . getVertex1 (); Point3d vertex2 = inTriangle . getVertex2 (); Vector3d edge1 = new Vector3d (); Vector3d edge2 = new Vector3d (); Vector3d h = new Vector3d (); Vector3d s = new Vector3d (); Vector3d q = new Vector3d (); double a , f , u , v ; edge1 . sub ( vertex1 , vertex0 ); edge2 . sub ( vertex2 , vertex0 ); h . cross ( rayVector , edge2 ); a = edge1 . dot ( h );if ( a > - EPSILON && a < EPSILON ) { return false ; // Este rayo es paralelo a este triángulo. }f = 1.0 / a ; s . sub ( rayOrigin , vertex0 ); u = f * ( s . dot ( h ));Si ( u < 0.0 || u > 1.0 ) { devolver falso ; }q . cruz ( s , borde1 ); v = f * rayoVector . punto ( q );Si ( v < 0.0 || u + v > 1.0 ) { devolver falso ; }// En esta etapa podemos calcular t para averiguar dónde está el punto de intersección en la línea. double t = f * edge2 . dot ( q ); if ( t > EPSILON ) // intersección de rayos { outIntersectionPoint . set ( 0.0 , 0.0 , 0.0 ); outIntersectionPoint . scaleAdd ( t , rayVector , rayOrigin ); return true ; } else // Esto significa que hay una intersección de línea pero no una intersección de rayos. { return false ; } } }