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:

- Triangles are precalculated according to the authors. Planes in general are defined as:
- 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:

- 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*:

- 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*) - 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.