Cramer's Rule for the calculation of alpha and beta in the 2nd book, chapter 6.5

[Bowen951209]

I want to propose another way to calculate $\alpha$ and $\beta$ in quad intersection in chapter 6.5 in the 2nd book.
In the book, we calculate $\alpha$ and $\beta$ using this formula:

$$ \alpha = \vec{w} \cdot (\vec{p} \times \vec{v}) $$ $$ \beta = \vec{w} \cdot (\vec{p} \times \vec{p}) $$

where $\vec{w} = \frac{\vec{n}}{\vec{n}\cdot\vec{n}}$.

However, I've never seen the trick described in the book to calculate this. Instead, I would use Cramer's Rule. For a quick reminder, if we want to solve:

$$ \begin{aligned} a_1 x + b_1 y &= c_1 \\ a_2 x + b_2 y &= c_2 \end{aligned} $$

The solutions for $x$ and $y$ are:

$$ x = \frac{\Delta_x}{\Delta}, \quad y = \frac{\Delta_y}{\Delta} $$ , where $$ \Delta = \begin{vmatrix} a_1 & b_1 \\ a_2 & b_2 \end{vmatrix} = a_1 b_2 - a_2 b_1\\ $$

$$ \Delta x = \begin{vmatrix} c_1 & b_1 \\ c_2 & b_2 \end{vmatrix} = c_1 b_2 - c_2 b_1 $$

$$ \Delta y = \begin{vmatrix} a_1 & c_1 \\ a_2 & c_2 \end{vmatrix} = a_1 c_2 - a_2 c_1 $$

If we assume that $\vec{p}$ can be the linear combination of $\vec{u}$ and $\vec{v}$, then there must be solutions for $\alpha$ and $\beta$.
Set $\vec{p} = (p_{x}, p_{y}, p_{z})$, $\vec{u} = (u_{x}, u_{y}, u_{z})$, $\vec{v} = (v_{x}, v_{y}, v_{z})$. And we want to solve the below equation for $\alpha$ and $\beta$:

$$ (p_{x}, p_{y}, p_{z}) = \alpha (u_{x}, u_{y}, u_{z}) + \beta (v_{x}, v_{y}, v_{z}) $$

In a clearer form:

$$ \begin{aligned} \alpha u_{x} + \beta v_{x} = x \\ \alpha u_{y} + \beta v_{y} = y \\ \alpha u_{z} + \beta v_{z} = z \end{aligned} $$

To solve $\alpha$ and $\beta$, you might think only two of the above equations are needed to get $\Delta$, $\Delta \alpha$, and $\Delta \beta$ in Cramer's Rule, but it's actually not the case. We can only use the components that ensures $\Delta$ is not $0$, and this how I check it in GLSL:

float delta, alpha, beta;

if((delta = u.x * v.y - u.y * v.x) != 0) {
    alpha = (planar_hitpt_vector.x * v.y - planar_hitpt_vector.y * v.x) / delta;
    beta = (planar_hitpt_vector.y * u.x - planar_hitpt_vector.x * u.y) / delta;
} else if((delta = u.x * v.z - u.z * v.x) != 0) {
    alpha = (planar_hitpt_vector.x * v.z - planar_hitpt_vector.z * v.x) / delta;
    beta = (planar_hitpt_vector.z * u.x - planar_hitpt_vector.x * u.z) / delta;
} else {
    delta = u.y * v.z - u.z * v.y;
    alpha = (planar_hitpt_vector.y * v.z - planar_hitpt_vector.z * v.y) / delta;
    beta = (planar_hitpt_vector.z * u.y - planar_hitpt_vector.y * u.z) / delta;
}

We simply find a non-zero $\Delta$ (and there must be at least one) and calculate the respective $\Delta \alpha, \Delta \beta$, and solve $\alpha = \frac{\Delta \alpha}{\Delta}$, $\beta = \frac{\Delta \beta}{\Delta}$.

In my implementation in OpenGL, I get a 10% boost in the 5 quads scene at max ray depth of 50 and 200 samples.

Book's approach:

My approach: