Dev Log 5: Ray cast! (part2) Bonus!

This part is a bonus as it wasnt going to be until last lines were written in part1.

It is about more precise ray intersection testing that not only test against an AABB but to concrete shapes and return the intersection point and normal (As seen in the last post gif).

I will discuss two shapes, box or oriented bounding box and sphere.

Sphere

To represent the sphere we need position and radius and to save processing costs I'll add radius square.

In vectorial form the sphere equation looks like this:

||P - C||^2 = r^2

Where P is the point to test, C is the center of the sphere and r is the radius. And the above equation can be written as:

dot((P - C), (P - C)) = r^2

Now is time for the intersection. We have 2 equations:

Sphere: dot((P - C), (P - C)) = r^2 Ray: t * D + O

The ray is now the Point (P) we want to test in the sphere so now the equation to solve becomes:

dot((t * D + O - C), (t * D + O - C)) = r^2

After expansion and rearrangement it becomes:

t^2 * dot(D, D) + 2t * dot(D, O - C) + dot(O - C, O - C) - r^2 = 0

As you can see this is in the form of a quadratic equation being t our unknown variable and the constants are:

a = dot(D, D) b = 2 * dot(D, O - C) c = dot(O - C, O - C) - r^2

Solving by Bhaskara formula:

t = (-b ± sqrt(b*b - 4 * a * c) ) / (2 * a)

where (b*b - 4 * a * c) is the discriminant.

• if discriminant < 0 the ray misses the sphere.

• if discriminant = 0 the ray touches the sphere. (one point intersection)

• if discriminant > 0 the ray intersect the sphere. (two points)

4 cases of intersection between ray and sphere

Pseudocode:

calculate_coheficients(a, b, c) calculate_discriminant(discriminant) if discriminant < 0 return false if discriminant == 0 t = -b / (2a) else sqrt_d = sqrt(discriminant) t0 = (-b + sqrt_d) / (2a) t1 = (-b - sqrt_d) / (2a) if t1 > 0 && t0 > 0 t = min(t0, t1) // If inside the sphere else t = t0 < 0 ? t1 : t0

At this point we could return true but first we will calculate the intersection point and the normal at that point.

We know that the point the ray hits is

D * t + O

and we now have all the values.

From the image on the left we can see that we can calculate the normal as:

P = D * t + O normal = (P - C).normalized

and then return true.

Box

To represent the box we need the half extents the position and the rotation (the last 2 I represent them with the transform matrix).

U and V are unit vectors defining box local space

It is very similar to the aabb-ray intersection but since we have more operations per dimension we are gonna enclose them into a for loop for better readability.

We start by getting the local axis and the half extent on that axis.

for the û axis: float half_extent = half_extent[0] vec3 axis = rot_matrix[0] (column-major)

Now we calculate the vector from ray origin to box center and the length on the previous axis.

Yellow: vector from ray origin to box center. Cyan: yellow projected on û

Also we calculate the proportion of the ray direction on the axis and now to early out we test if we didn't intersect and remember if one axis match this conditions that's it no intersection happens. The condition is:

|direction_on_axis| < 0 && (-distance_on_axis - half_extent > 0 || -distance_on_axis + half_extent > 0)

Where:

float direction_on_axis = dot(axis, ray_direction) float distance_on_axis = dot(axis, center - ray_origin)

Then we do the same as in an AABB vs ray and say:

float t0 = (distance_on_axis + half_extent) / direction_on_axis float t1 = (distance_on_axis - half_extent) / direction_on_axis if (t0 > t1) swap(t0, t1) tNear = max(t0, tNear) tFar = min(t1, tFar) if (tNear > tFar || tFar < 0) return false

Where tNear and tFar are outside the loop. And then after the loop we find the solution and the position of the hit:

solution = (tNear > 0) ? tNear : tFar hit_position = ray_origin + ray_direction * solution

Now for the normal, surprisingly you already have all the normals since they are the local axes of the box so what we need is to know wich axis and what sign is.

pseudocode:

int sgn = 1 vec3 axis for every axis: normal_end_point = center + axis * half_extent if abs(dot(hit_point - normal_end_point, axis)) < epsilon: break else: normal_end_point_inv = center - axis * half_extent if abs(dot(hit_point - normal_end_point_inv, axis)) < epsilon: sgn = -1; break normal = sgn * axis return true

where center is the center of the box and hit_pos is the hit position we calculated earlier.

End result

That's all for this post, folk. Bye!