Intersekcja promienia z trójkątem

Równanie promienia w postaci parametrycznej:

$P = S + \vec{d}t$

Współrzędne punkt trójkąta 3D, mając jego współrzędne barycentryczne liczymy:

$P=P_0+u(P_1-P_0)+v(P_2-P_0)$

Z porównania:

$S + \vec{d}t = P_0+u(P_1-P_0)+v(P_2-P_0)$

Porządkując względem niewiadomych t, u, v:

$S + \vec{d}t = P_0 + uP_1 - uP_0 + vP_2 - vP_0$

$t\vec{d} - u(P_1-P_0) - v(P_2-P_0)=P_0-S$

$-t\vec{d} + u(P_1-P_0) + v(P_2-P_0)=S-P_0$

$\left [ -\vec{d}, P_1-P_0, P_2-P_0 \right ] \begin{bmatrix}t \\ u \\ v \end{bmatrix}=S-P_0$

Oznaczamy $E_1=P_1-P_0$, $E_2=P_2-P_0$:

$\left [ -\vec{d}, E_1, E_2 \right ] \begin{bmatrix}t \\ u \\ v \end{bmatrix}=S-P_0$

Korzystając ze wzorów Cramera rozwiązanie ma postać:

$\begin{bmatrix}t \\ u \\ v \end{bmatrix}=
\displaystyle\frac{1}{\left | -\vec{d}, E_1, E_2 \right |}
\begin{bmatrix}
\left | T, E_1, E_2 \right | \\
\left | -\vec{d}, T, E_2 \right | \\
\left | -\vec{d}, E_1, T \right |
\end{bmatrix}$

Gdzie wyrażenia typu $\left | T, E_1, E_2 \right |$ oznaczają wyznacznik macierzy w której trzy wiersze definiują wektory $T, E_1, E_2$.

Korzystając z właściwości iloczynu mieszanego i iloczynu wektorowego możemy zapisać:

$\begin{bmatrix}t \\ u \\ v \end{bmatrix}=\displaystyle\frac{1}{E_1 \cdot (\vec{d} \times E_2)}\begin{bmatrix}
E_2 \cdot (T \times E_1) \\
T \cdot (\vec{d} \times E_2) \\
\vec{d} \cdot (T \times E_1)
\end{bmatrix}$

Oznaczamy $P=\vec{d}\times E_2$ i $Q=T\times E_1$:

$\begin{bmatrix}t \\ u \\ v \end{bmatrix}=\displaystyle\frac{1}{E_1 \cdot P}\begin{bmatrix}
E_2 \cdot Q \\
T \cdot P \\
\vec{d} \cdot Q
\end{bmatrix}$

Sprawdzając znak t upewniamy się czy przecięcie nie nastąpiło przed punktem startu promienia S.

Za pomocą u i v sprawdzamy czy punkt przecięcia nastąpił we wnętrzu trójkąta. Otrzymane niejako przy okazji współrzędne barycentryczne (u,v) przydają się podczas np. teksturowania, warto je sobie gdzieś zapamiętać, by nie musieć ich liczyć ponownie.

Jeśli $E_1 \cdot P=0$ układ nie ma rozwiązania. Promień biegnie równolegle do płaszczyzny w której zawiera się trójkąt. Wartość tego wyrażenia jest długością rzutu wektora kierunkowego prostej na normalną trójkąta. Badając znak $E_1 \cdot P=E_1$ ustalamy z której strony trójkąta uderzył w niego promień.

Przykład implementacji w C#:

public override Intersection GetIntersection(Ray a_ray)
{
    if (a_ray.PrevHit == this)
        return Scene.NoIntersection;

    if (OneSide)
    {
        if ((a_ray.Dir * Normal) > 0)
            return Scene.NoIntersection;
    }

    var E1 = V2 - V1;
    var E2 = V3 - V1;
       
    Vector3 P = Vector3.CrossProduct(a_ray.Dir, E2);
    double a = E1 * P;
            
    if (a == 0)
        return Scene.NoIntersection;
            
    Vector3 T = a_ray.Start - V1;
    double u = (P * T) / a;

    if (u < 0)
    {
        if (u < -INTERSECTION_ERROR) 
            return Scene.NoIntersection;
        u = 0;
    }

    if (u > 1)
    {
        if (u > 1 + INTERSECTION_ERROR)
            return Scene.NoIntersection;
        u = 1;
    }
            
    Vector3 Q = Vector3.CrossProduct(T, E1);
    double v = (a_ray.Dir * Q) / a;
            
    if (v < 0)
    {
        if (v < -INTERSECTION_ERROR) 
            return Scene.NoIntersection;
        v = 0;
    }
    if (u + v > 1.0)
        return Scene.NoIntersection;
            
    double dist = (Q * E2) / a;
            
    if (dist > 0)
    {
        bool backHit = a < 0;
            
        if (backHit && OneSide)
            return Scene.NoIntersection;
        else
        {
            return new Intersection()
            {
                SceneObject = this,
                SourceRay = a_ray,
                Dist = dist,
                Scene = Scene,
                BackHit = backHit,
                Pos = a_ray.HitPoint(dist),
                BaricentricUV = new Vector2(u, v),
                UV = BaricentricToUV(u, v)
            };
        }
    }
    else
        return Scene.NoIntersection;
}

Kontrola wartości (u,v) została tutaj zmodyfikowana o uwzględnienie błędów obliczeniowych. Dzięki czemu dwóm przystającym trójkątom nie zdarza się mieć prześwitów na przystającej krawędzi. Zastosowana w kodzie metoda na pewno zwiększa pole trójkąta, ale na tyle nieznacznie, że nie powinno mieć to znaczenia. Poza tym wymiana kolejności testowania ich intersekcji z promieniem spowoduje wygenerowanie (nieznacznie) różnej sceny. Stała INTERSECTION_ERROR została ustalona przeze mnie w kodzie testowym. Na chwile obecną jej wartość to 1e-13. I jest sto razy większa niż wyszło to w kodzie testowym. Na kod testowy składa się scena złożona z niebieskiego kwadratu, zaraz za nimi jest komplet żółtych kwadratów dokładnie pokrywający przekątną. Całą sceną losowo kręcimy lub przemieszczamy. Za każdym razem bierzemy nową scenę by uniknąć błędów obliczeniowych. Położenie kwadratu względem kamery nie zmienia. Próbkujemy więc na sztywno piksele wokół przekątnej i sprawdzamy czy są niebieskie. I tak z kilka tysięcy razy dla pewności.