# Implementing ray-triangle collision algorithm

Hello,

## Introduction:

For my game I am looking for deterministic physics and for that I use fixed-point math I wrote myself.
That is why you will see the following things in my code that I want to explain first:

• The type `Q` is my fixed number which is a `(i32)` type. The math works like normal integer math, it just shifts bits at multiplications and division. The least significant 8 bits are the fractional bits, which means the precision of the number is 1/256.
• The type `Q3` is a tuple of 3 `Q` numbers. I use it for vector math. Note that the implementation of `Q3 * Q3` is the dot-product.
• F64Vec is a `(f64, f64, f64)` type. I only use it to generate the triangle data and plan to replace it in the future, I just created it to test if the algorithm works. That is why you will see nasty `sub(a, b)` methods instead of proper `op`'s.
• Let us assume for now that within my Q-Number system there are no big problems for my project here so I don't have to overhelm you with much more code.

Part of the game physics is of course collision checking. I decided to use the algorithm by Jiri Havel and Adam Herout. You can download the paper about it here:
(Note the code for the intersect algorithm is included in that paper in C Language)

I will qoute some things from that document I think that is relevant. If you want to help me, please also check the document itself in case I missed or misunderstood something.

## Algorithm:

1. Triangles are precalculated according to the authors. Planes in general are defined as:
2. The triangle is made of three planes, first the plane that descibes the surfase of the triangle:
...and second and third which are defined as:
3. The ray that might hit the triangle is defined as this (note that O and P are three-dimensional), with P being a specific point when we fill in the parameter t:
4. Within the algorithm, several values are needed which are:

(variables with ' at them are prefixed with det in my code, so t' becomes det_t)
5. The ray hits the triangle when the following conditions are met:

Note that t_max is the maximum value for t in the ray, which sets the maximum length of the ray.
If the ray did hit the triangle, then the new t is calculated according to 4. and is given back from the function.

## My code:

This is the Triangle struct:

``````pub struct Triangle{
n0: Q3, d0: Q,  //triangle plane
n1: Q3, d1: Q,  //first support plane
n2: Q3, d2: Q,  //second support plane
}
``````

First I calculate the needed data for a triangle from three points. I wrote this function:

``````impl Triangle{
pub fn from_corners(a: &F64Vec, b: &F64Vec, c: &F64Vec) -> Self{
//calculate triangle plane
let ab = F64Vec::sub(&b, &a);
let ac = F64Vec::sub(&c, &a);
let n0_f64 = F64Vec::cross_product(&ab, &ac);
let n0_length_square = n0_f64.0 * n0_f64.0 + n0_f64.1 * n0_f64.1 + n0_f64.2 * n0_f64.2;
let d0_f64 = -F64Vec::dot_product(&n0_f64, &a);
//calculate support planes
let n1_f64 = F64Vec::cross_product(&ac, &n0_f64).div(n0_length_square);
let d1_f64 = -F64Vec::dot_product(&n1_f64, &a);
let n2_f64 = F64Vec::cross_product(&n0_f64, &ab).div(n0_length_square);
let d2_f64 = -F64Vec::dot_product(&n2_f64, &a);
//convert to fixed-point math
Self {
n0: n0_f64.to_q3(), d0: Q::from_f64(d0_f64),
n1: n1_f64.to_q3(), d1: Q::from_f64(d1_f64),
n2: n2_f64.to_q3(), d2: Q::from_f64(d2_f64)
}
}
}
``````

This is the intersect algorithm: (ray_length is t_max from the paper)

``````impl Triangle{
pub fn intersect(&mut self, ray_origin: &Q3, ray_direction: &Q3, ray_length: Q) -> Option<Q>{
let det: Q = self.n0 * ray_direction;
let det_t: Q = self.d0 - (self.n0 * ray_origin);
if Q::ZERO == Q::SIGN_MASK & (det_t ^ (det * ray_length - det_t)){
let det_p: Q3 = *ray_origin * det + *ray_direction * det_t;
let det_u: Q = det_p * self.n1 + det * self.d1;
let det_sub_det_u: Q = det - det_u;
if Q::ZERO == Q::SIGN_MASK & (det_u ^ det_sub_det_u){
let det_v: Q = det_p * self.n2 + det * self.d2;
self.debug_det: Q = det;
self.debug_det_u: Q = det_u;
self.debug_det_v: Q = det_v;
if Q::ZERO == Q::SIGN_MASK & (det_v ^ (det_sub_det_u - det_v)){
return Some(det_t / det);
}
}
}
None
}
}
``````

I also wrote a function that prints all data from an intersection check case, so I can debug more easily:

``````impl Triangle{
pub fn print_intersect_variables(&self, ray_origin: &Q3, ray_direction: &Q3, ray_length: Q){
println!("Triangle plane: n = ({}, {}, {}), d = {}",
(self.n0).0.to_f32(), (self.n0).1.to_f32(), (self.n0).2.to_f32(), self.d0.to_f32());
println!("Support plane 1: n = ({}, {}, {}), d = {}",
(self.n1).0.to_f32(), (self.n1).1.to_f32(), (self.n1).2.to_f32(), self.d1.to_f32());
println!("Support plane 2: n = ({}, {}, {}), d = {}",
(self.n2).0.to_f32(), (self.n2).1.to_f32(), (self.n2).2.to_f32(), self.d2.to_f32());
let det = self.n0 * ray_direction;
let det_t = self.d0 - (self.n0 * ray_origin);
let det_p = *ray_origin * det + *ray_direction * det_t;
let det_u = det_p * self.n1 + det * self.d1;
let det_v = det_p * self.n2 + det * self.d2;
let det_inv = Q::ONE / det;
let t = det_t * det_inv;
let u = det_u * det_inv;
let v = det_v * det_inv;
let eq_det_t_sign = det * ray_length - det_t;
let eq_det_u_sign = det - det_u;
let eq_det_v_sign = det - det_u - det_v;
println!("det: {}, det_t: {}, det_u: {}, det_v: {}", det.to_f32(), det_t.to_f32(), det_u.to_f32(), det_v.to_f32());
println!("t: {}, u: {}, v: {}", t.to_f32(), u.to_f32(), v.to_f32());
println!("1st if: sign({}) == sign({})", det_t.to_f32(), eq_det_t_sign.to_f32());
println!("2nd if: sign({}) == sign({})", det_u.to_f32(), eq_det_u_sign.to_f32());
println!("3rd if: sign({}) == sign({})", det_v.to_f32(), eq_det_v_sign.to_f32());
}
}
``````

The algorithm works with very simple triangles on the x-z plane. But with more complcated normal vectors it fails.

For example this triangle:

``````//A(2,1,1), B(1,2,2), C(3,2,1)
let c0 = F64Vec::new(2f64,1f64,1f64);
let c1 = F64Vec::new(1f64,2f64,2f64);
let c2 = F64Vec::new(3f64,2f64,1f64);
let mut t = Triangle::from_corners(&c0, &c1, &c2);
``````

With that I already know the corners would be points that must result in a collision. So I apply the following ray to the triangle and let it print the algorithm's variables:

``````//attempting to hit point B(1,2,2)
let ray_origin = Q3::new(
Q::ONE,
Q::ONE*4, //ray_origin points to a position that is 2 above B (with y being the vertical compund)
Q::ONE*2
);
let ray_direction = Q3::new(
Q::ZERO,
-Q::ONE, //vector is pointing straight down with the length of 1
Q::ZERO
);
t.print_intersect_variables(&ray_origin, &ray_direction, Q::ONE*10); //t_max = 10 * ray_direction.length
``````

This prints the following text.

``````Triangle plane: n = (-1, 1, -2), d = 3
Support plane 1: n = (-0.33203125, 0.33203125, 0.33203125), d = 0
Support plane 2: n = (0.5, 0.5, 0), d = -1.5
det: -10, det_t: 4, det_u: -29.882813, det_v: -30
t: -0.390625, u: 2.9179688, v: 2.9296875
1st if: sign(4) == sign(-14)
2nd if: sign(-29.882813) == sign(19.882813)
3rd if: sign(-30) == sign(49.882813)
``````

## Problem:

Let's investigate the output. You can see that all three if cases would fail because the signs are different. This is incorrect if my algorithm would work correctly.
Note that I am aware that the conversation from f64 to my fixed point math results in rounding errors. But other values help us to comprehend the situation better:

Even if we barely miss the point because of rounding errors, t, u, and v must be somewhere near their required values. But they aren't as you can see at the output:
`t: -0.390625, u: 2.9179688, v: 2.9296875`

Correct values for this intersection would be
`t: 0.2, u: 1, v: 0`

So what did I do wrong? I know this takes a lot of time and quite some math skills to get into this matter before even being able to find the error, so I highly appreciate anyone who tries to help me. If I left something unclear, please ask.

Thank you.

I did some hand-math to calculate the triangle's `n0` and `d0` and calculated also `det` and `det_t` with these.

I did that to see if I did implement the algorithm incorrectly or if it works as I want it to work.
I came to the same failed first if check, which means so far the algorithm I wrote in Rust is the same I have in my head. But it also means that the algorithm in my head is incorrect.

So this became more a math than a Rust question, as the task now is to find out what I misunderstood in the paper. I am not sure if it is still justified to be an open help topic in this Rust forum, but in case it isn't, it can be closed and I go to a math-orientated forum. @moderators

This topic was automatically closed 90 days after the last reply. We invite you to open a new topic if you have further questions or comments.