Torus - równanie i intersekcja z promieniem

Mały promień torusa oznaczmy przez r, duży przez R.

Równanie okręgu: $x^2+y^2=r^2$, którego obrót tworzy torus przesuniętego o R to: $(x-R)^2+y^2=r^2$

Utwórzmy torus poprzez obrót tego okręgu wokół osi Y. Jeśli wyobrazimy sobie taki torus i przetniemy go dowolną płaszczyzną przechodzącą przez oś Y to nasz torus jest definiowany przez podane wyżej równanie 2D. Zauważmy, że x w nim w przypadku 3D oznacza: $\sqrt{(x^2+z^2)}$, czyli równanie torusa w 3D ma postać:

$(\sqrt{x^2+z^2}-R)^2+y^2=r^2$

W zależności od stosunku R do r możemy otrzymać torusy o różnym kształcie. W tym sferę dla R=0 (w zasadzie dwie nakładające się na siebie sfery). Okręg w 3D dla r=0. Torus z dziurką dla R>r i torusy bez dziurki dla r>R. Dla takich przypadków wraz jak r staje się większe od R to nasz torus coraz bardziej upodabnia się do sfery, a wgłębienia zanikają.

Równanie to możemy zapisać jako:

$(\sqrt{x^2+z^2}-R)^2+y^2=r^2$
$x^2+z^2-2R\sqrt{x^2+z^2}+R^2+y^2=r^2$
$x^2+y^2+z^2+R^2-r^2=2R\sqrt{x^2+z^2}$
$(x^2+y^2+z^2+R^2-r^2)^2=4R^2(x^2+z^2)$

Podstawiając do powyższego parametryczne równanie linii o punkcie startu S i kierunku D:

$P = S+ \vec D \cdot t$
$x = x_s + x_d t$
$y = y_s + y_d t$
$z = z_s + z_d t$

, otrzymujemy:

$\begin{align*} (D \cdot D \cdot t^2 + 2 \cdot D \cdot S \cdot t + S \cdot S+ R^2 - r^2)^2= & 4R^2((D \cdot D- y_dy_d) \cdot t^2 + \\ & 2 \cdot (D \cdot S- y_dy_s) \cdot t + \\ & S\cdot S- y_sy_s ) \end{align*}$

, ponieważ wektor $\vec D$ jest znormalizowany to:


$\begin{align*} (t^2 + 2 \cdot D \cdot S \cdot t + S \cdot S+ R^2 - r^2)^2= & 4R^2((1 - y_dy_d) \cdot t^2 + \\ & 2 \cdot (D \cdot S- y_dy_s) \cdot t + \\ & S\cdot S- y_sy_s ) \end{align*}$

, po oznaczeniu:

$\beta = 2 \cdot S \cdot D$
$\gamma = S \cdot S + R^2 - r^2$

, otrzymujemy:

$(t^2 + \beta t + \gamma)^2=4R^2((1 - y_dy_d) \cdot t^2 + (\beta-2y_dy_s) \cdot t + S \cdot S - y_sy_s)$

Stąd już możemy łatwo wyliczyć potrzebne nam współczynniki do rozwiązania wielomianu stopnia czwartego:

$At^4 + Bt^3 + Ct^2 + Dt + E = 0$

$A = 1
$B = 2 \beta$
$C = \beta ^ 2 + 2 \gamma + (2 R \cdot D.Y)^2$
$D = 2 \beta \gamma + 8 R^2 \cdot S.Y \cdot D.Y$
$E = \gamma ^ 2 + (2 R \cdot S.Y)^2 - (2 R r)^2$

Równanie to jest znormalizowane, tzn. $A=1$. Z czterech możliwych rozwiązań rzeczywistych liczy się dla nas to najmniejsze, ale większe od zera, a dokładniej większej od minimalnej odległości pomiędzy obiektami Constants.MINIMAL_DISTANCE.

Teraz musimy jeszcze ustalić w którą stronę torusa uderzył promień. W tym celu najlepiej wyliczyć normalną i sprawdzić kąt pomiędzy nią, a promieniem. O tym jak ją liczyć będzie w innym poście. Tak samo jak o metodach rozwiązywania wielomianów.

Przykładowy kod:

public override Intersection GetIntersection(
    Intersection a_source_ray_intersection, Ray a_ray)
{
    Vector3 local_dir = default(Vector3);
    Vector3 local_start = default(Vector3);
    bool back_hit = false;

    if (!TransformToLocalToBoundBox(a_ray, ref local_dir, ref local_start))
        return Scene.NoIntersection;

    double r2 = m_local_small_radius * m_local_small_radius;
    double R2 = m_local_big_radius * m_local_big_radius;
    double beta = 2 * local_dir * local_start;
    double gamma = local_start.SqrLen - r2 - R2;

    double b = 2 * beta;
    double c = beta * beta + 2 * gamma + 4 * R2 * local_dir.Y * local_dir.Y;
    double d = 2 * beta * gamma + 8 * R2 * local_start.Y * local_dir.Y;
    double e = gamma * gamma + 4 * R2 * local_start.Y * local_start.Y - 4 * R2 * r2;

    double dist = m_solver.SolveQuartic(1, b, c, d, e);

    if (dist == Double.PositiveInfinity)
        return Scene.NoIntersection;

    Vector3 local_pos = local_start + local_dir * dist;

    Vector3 normal = (local_pos - 
        new Vector3(local_pos.X, 0, local_pos.Z).Normalized * m_local_big_radius).Normalized;

    back_hit = normal * local_dir > 0;

    if (back_hit && OneSide)
        return Scene.NoIntersection;

    Vector3 world_pos = LocalToWorld * local_pos;
    double world_dist = (world_pos - a_ray.Start).Length;

    Intersection intersection = new Intersection()
    {
        PrevIntersection = a_source_ray_intersection, 
        SceneObject = this,
        SourceRay = a_ray,
        LocalPos = local_pos,
        Dist = world_dist,
        Scene = Scene,
        BackHit = back_hit,
        Pos = world_pos
    };

    // Numerical errors.
    if (intersection.Normal * intersection.SourceRay.Dir >= 0)
        return Scene.NoIntersection;

    return intersection;
}