Implementing ray-triangle collision algorithm

[Urben]

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.

Task:

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:
http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.159.5135&rep=rep1&type=pdf
(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.

[Urben]

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