Processing math: 100%

Lecture 9: Triangles

For spheres and planes, we plugged the ray equation p(t)=po+t∗→d into the equation of the shape in question to solve for intersection. But what is the equation for a triangle? We have the plane equation, and we have line equations for each edge. We could find the plane intersection point, then check against each line. In fact, this is how you can support arbitrary polygons, but there are more convenient representations for triangles.

Review: Barycentric Coordinates

Barycentric coordinates of point p in terms of a, b, and c are the numbers α, β, and γ such that:

p=α∗a+β∗b+γ∗c

With the constraint:

α+β+γ=1

This must be true if the three values are barycentric coordinates.

However, the values don’t have to be positive.

A point p is inside the the triangle if and only if:

0≤α≤1 0≤β≤1 0≤γ≤1

What does α=0 mean? Any point opposite the vertex a, or rather along the line from b to c.

What about α,β=0? In this case we know γ must be 1, this is the c vertex.

The center of the triangle is α=β=γ=13

One interesting thing to note is that we’re really defining a coordinate system with, e.g., origin a and basis vectors b−a and c−a

p=a+β∗(b−a)+γ∗(c−a)

Intersection

Let’s plug in our standard ray equation:

P0+t∗→d=a+β∗(b−a)+γ∗(c−a)

Since this is a vector form it is really three different equations:

P0x+t∗xd=xa+β∗(xb−xa)+γ∗(xc−xa) P0y+t∗yd=ya+β∗(yb−ya)+γ∗(yc−ya) P0z+t∗zd=za+β∗(zb−za)+γ∗(zc−za)

Here we have three equations and three unknowns - t, β, and γ. It’s a linear system we can solve using matrix algebra. First let’s move all the unknowns to one side:

(xa−xb)∗β+(xa−xc)∗γ+xd∗t=xa−P0x (ya−yb)∗β+(ya−yc)∗γ+yd∗t=ya−P0y (za−zb)∗β+(za−zc)∗γ+zd∗t=za−P0z

Linear System

Let’s rewrite in matrix form:

[xa−xbxa−xcxdya−ybya−ycydza−zbza−zczd][βγt]=[xa−P0xya−P0yza−P0z]

This is a linear system of the form:

Ax=b

This 3x3 linear system can be solved by finding the inverse of the matrix.

Ax=b A−1Ax=A−1b x=A−1b

However, we don’t need to compute the inverse explicitly - there is a faster option.

Cramer’s Rule

Given a linear system of the form

Ax=b

with

x=[x1x2x3] xi=det(Ai)det(A)

where Ai is the matrix formed by replacing the ith column of A by the column vector b.

The denominator det(A) gives us an obvious way to test for intersection - if this is 0, there is no intersection. Cramer’s rule is also convenient because it permits solving for individual elements of the solution vector. So, we can solve for t and check if it’s positive. Them, we can solve for a and check 0≤α≤1, etc.

Determinants

Recall, to compute the 2x2 determinant:

det|abcd|=ad−bc

And for a 3x3 determinant:

det|abcdefghi| =a∗det|efhi|−b∗det|dfgi|+c∗det|degh| =aei−afh−bdi+bfg+cdh−ceg

Note the pattern of (+) (-) (+) for the three 2x2 determinants. We can compute this “expansion” of the determinant using any row or column, not just the first row (as in this example). However, if you use the middle row or middle column, you flip the signs to (-) (+) (-)

Solve

Below are the results of applying Cramer’s rule. Note the vertical bars around a matrix mean to compute the determinant.

β=|xa−P0xxa−xcxdya−P0yya−ycydza−P0zza−zczd||A| γ=|xa−xbxa−P0xxdya−ybya−P0yydza−zbza−P0zzd||A| t=|xa−xbxa−xcxa−P0xya−ybya−ycya−P0yza−zbza−zcza−P0z||A|

Implementation

Common Subterms

3x3 matrices have common subterms that can be exploited by carefully naming variables. Lets look at generic form with dummy variables.

[adgbehcfi][βγt]=[jkl]

Solving for determinants expands to (here we expand the determinant on the third column purposefully):

β=|jdgkehlfi||adgbehcfi| β=j(ei−hf)+k(gf−di)+l(dh−eg)a(ei−hf)+b(gf−di)+c(dh−eg)

The expression (ei−hf) is found in two places - instead of computing this twice, you can create a variable for it, compute once, and use it twice. Of course, also save the denominator value - it is the same for all three variables.

Similarly, we can purposefully expand both γ and t computations to share expressions.

γ=|ajgbkhcli|det(A) γ=i(ak−jb)+h(jc−al)+g(bl−kc)det(A) t=|adjbekcfl|det(A) t=−d(bl−kc)+e(al−jc)−f(ak−jb)det(A)

Which we can then manipulate to create common terms with the γ solution.

t=−f(ak−jb)−e(jc−al)−d(bl−kc)det(A)

Arguably, this also makes cleaner and easier to read/verify code.

Pseudocode

hit_tri(ray r, vert a, vert b, vert c, t0, t1)

    compute t
    if (t < t0) || (t > t1) then
        return false

    compute gamma
    if (gamma < 0) || (gamma > 1) then
        return false

    compute beta
    if (beta < 0) || (beta > 1 - gamma) then
        return false

    return true