Stożek

Utwórzmy w lokalnym układzie współrzędnych walec którego osią symetri jest Y:

$X^2+Y^2=R^2$

Uzależnijmy jego promień od wysokości Y:

$X^2+Z^2=Y^2$

W ten sposób otrzymujemy równanie powierzchni, której częścią jest stożek. Jego środek znajduje się w środku układu współrzędnych. Jego podstawa niech znajduje się po dodatniej osi Y w odległości H od środka układu współrzędnych.

Dokonując transformacji promienia do lokalnego układu współrzędnych w którym leży nasz stożek skalujemy go także odpowiednio tak by $R=1$ i $H=1$.

Środek naszego stożka znajduje się w połowie jego wysokości.

Transformacja z lokalnego układu współrzędnych do układu współrzędnych świata ma postać:

protected override Vector3 LocalScale
{
    get
    {
        return new Vector3(Radius, Height, Radius);
    }
}

protected override Matrix4 GetLocalToWorldMatrix()
{
    return Matrix4.CreateTranslation(Pos) *
            new Matrix4(Right, Up, Forward) *
            Matrix4.CreateTranslation(0, Height / 2, 0) * 
            Matrix4.CreateRotateX(180) * 
            Matrix4.CreateScale(Scale.Scale(LocalScale));
}

W kolejności najpierw skalujemy stożek, później obracamy go o 180 stopni, tak by czubek znajdował się po dodatniej stronie osi Y. Następnie przesuwamy go tak by jego podstawa znajdowała się na osi XZ. I na końcu transformujemy go do układu współrzędnych świata (rotacja + translacja). Dzięki temu do obliczeń mamy stożek w postaci najprostszej i najszybszej do obliczeń (czubek w środku układu współrzędnych), zaś w układzie świata w postaci intuicyjnej leżącego podstawą na płaszczyźnie XZ.

Podstawiając do równania powierzchni tworzącej parametryczne równanie linii (DIR to wektor znormalizowany):

$P = START + DIR \cdot t$
$x = x_s + x_d t$
$y = y_s + y_d t$
$z = z_s + z_d t$

,otrzymyjemy równanie kwadratowe $At^2 + Bt + C = 0$ o współczynnikach:

$A = x_d^2 + z_d^2 - y_d^2 = DIR \cdot DIR - 2 y_d^2 = 1 - 2 y_d^2$
$B = 2(x_d x_s + z_d z_s - y_d y_s) = 2 \cdot START \cdot DIR - 4 y_d y_s$
$C = x_s^2 + z_s^2 - y_s^2 = START \cdot START - 2 y_s^2$

Jego rozwiązanie wyznaczy nam dystans do boków tworzących stożek. Interesuje nas minimalne rozwiąnie większe od Constants.MINIMAL_DISTANT dla którego współrzędna Y punktu intersekcji spełnia warunek $Constants.MINIMAL\_DISTANT < Y < 1$.

Następnie znajdujemy odległość do płaszczyzny tworzącej stożek korzystając z wzoru na intersekcje promienia z płaszczyzną:

$D = \displaystyle\frac{1-START.Y}{DIR.Y}$

Jeśli odległość ta jest miejsza od poprzednio wyznaczonego dystansu do tworzących stożek możemy wyznaczyć punkt intersekcji i sprawdzić czy znajduje się on w promieniu podstawy:

$POS.X^2 + POS.Z^2 < 1$

Punkt intersekcji został wyznaczony. Teraz trzeba sprawdzić czy uderzenie nastąpiło od zewnątrz czy wewnątrz.

Ponieważ stożek jest bryłą zamkniętą i wypukłą, możemy mieć pewność że promień odbity od wnętrza uderzy we wnętrze stożka. Promień odbity od zewnątrz stożka na pewno w niego nie uderzy. Promień refrakcyjny wychodzący na zewnątrz także w stożek nie uderzy. Promień refrakcyjny wchodzący do stożka na pewno w niego uderzy od środka.

W pozostałych przypadkach, gdy punktem wyjścia nie jest poprzednie uderzenie w stożek musimy liczyć. Jeśli uderzenie nastąpiło w podstawę stożka to wystarczy sprawdzić wartość współrzędnej START.Y. Jeśli uderzenie nastąpiło w bok stożka to uderzenie nastąpiło wewnątrz jeśli $POS.X^2 + POS.Z^2 < POS.Z^2$. Tylko pod warunkiem, że punkt startu promienia znajdował się w płaszczyźnie stożka: $0 \le START.y \le 1$

public override Intersection GetIntersection(
    Intersection a_source_ray_intersection, Ray a_ray)
{
    bool back_hit_expected = false;

    if ((a_source_ray_intersection != null) && 
        a_source_ray_intersection.SceneObject == this)
    {
        Debug.Assert((a_ray.RaySurfaceSide == RaySurfaceSide.RefractedSide) ||
                        (a_ray.RaySurfaceSide == RaySurfaceSide.ReflectedSide));
        if (a_source_ray_intersection.BackHit ^ (a_ray.RaySurfaceSide ==  
            RaySurfaceSide.ReflectedSide))
        {
            return Scene.NoIntersection;
        }
        else
        {
            if (OneSide)
                return Scene.NoIntersection;
            else
                back_hit_expected = true;
        }
    }

    Vector3 local_dir = default(Vector3);
    Vector3 local_start = default(Vector3);

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

    double dist = Double.PositiveInfinity;

    PolynomialSolverAlgebraic solver = new PolynomialSolverAlgebraic(
        new Polynomial(
            1 - 2 * local_dir.Y * local_dir.Y,
            2 * (local_dir.X * local_start.X + local_dir.Z * local_start.Z - 
                local_dir.Y * local_start.Y),
            local_start.X * local_start.X + local_start.Z * local_start.Z - 
                local_start.Y * local_start.Y));

    double dist1 = default(Double), dist2 = default(Double);
    int roots = solver.SolveAlgebraicQuadratic(ref dist1, ref dist2);

    if (roots == 0)
        return Scene.NoIntersection;

    if ((roots >= 1) && (dist1 > Constants.MINIMAL_DISTANT))
    {
        double ly = local_start.Y + local_dir.Y * dist1;

        if ((ly >= Constants.MINIMAL_DISTANT) && (ly < 1))
            dist = dist1;
    }

    if ((roots == 2) && (dist2 > Constants.MINIMAL_DISTANT) && (dist2 < dist))
    {
        double ly = local_start.Y + local_dir.Y * dist2;

        if ((ly >= Constants.MINIMAL_DISTANT) && (ly < 1))
            dist = dist2;
    }

    double base_dist = (1 - local_start.Y) / local_dir.Y;
    bool base_hit = false;

    if ((base_dist > Constants.MINIMAL_DISTANT) && (base_dist < dist))
    {
        double lx = local_start.X + local_dir.X * base_dist;
        double lz = local_start.Z + local_dir.Z * base_dist;

        if (lx * lx + lz * lz < 1)
        {
            base_hit = true;
            dist = base_dist;
        }
    }

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

    bool back_hit = false;

    if (back_hit_expected)
        back_hit = true;
    else
    {
        if (base_hit)
            back_hit = local_start.Y < 1;
        else
        {
            if ((local_start.Y > 0) && (local_start.Y < 1))
            {
                back_hit = local_start.X * local_start.X + local_start.Z * 
                    local_start.Z < local_start.Y * local_start.Y;
            }
        }
    }

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

    Vector3 local_pos = local_start + local_dir * dist;
    Vector3 world_pos = LocalToWorld * local_pos;
    dist = (world_pos - a_ray.Start).Length;

    return new Intersection()
    {
        PrevIntersection = a_source_ray_intersection,
        SceneObject = this,
        SourceRay = a_ray,
        Dist = dist,
        Scene = Scene,
        BackHit = back_hit,
        Pos = world_pos,
        LocalPos = local_pos,
        Material = base_hit ? BaseMaterial : DefaultMaterial
    };
}

Normalną w miejscu uderzenia wyznaczymy z gradientu funkcji:

$f(x,y,z)=X^2+Z^2-Y^2$
$\nabla f=[2X,2Z,-2Y]=[X,Z,-Y]$

Jeśli uderzenie nastąpiło od środka wartość normalnej bierzemy ze znakiem ujemnym.

public override Vector3 GetNormal(Intersection a_intersection)
{
    if (a_intersection.LocalPos.Y.AlmostEquals(1))
    {
        if (a_intersection.BackHit)
            return Up;
        else
            return -Up;
    }
    else
    {
        if (a_intersection.BackHit)
        {
            return (LocalToWorldNormal * new Vector3(-a_intersection.LocalPos.X,
                a_intersection.LocalPos.Y, -a_intersection.LocalPos.Z)).Normalized;
        }
        else
        {
            return (LocalToWorldNormal * new Vector3(a_intersection.LocalPos.X,
                -a_intersection.LocalPos.Y, a_intersection.LocalPos.Z)).Normalized;
        }
    }
}

Kod na wyliczenie współrzędnych UVW (0..1):

public override Vector3 GetUVW(Intersection a_intersection)
{
    Vector3 uvw = base.GetUVW(a_intersection);

    return new Vector3(uvw.X / 2 + 0.5, uvw.Y, -uvw.Z / 2 + 0.5);
}

Współrzędna Y jest odwracana tak by góra tesktury (Y=0) znajdowała się na czubku stożka.

Dla UV jeśli uderzenie nastąpiło w podstawę stożka mapujemy współrzędną wprost. Dla boku stożka nawijamy na niego teksturę. Tutaj kod jest podobny do kodu dla sfery.

public override Vector2 GetUV(Intersection a_intersection)
{
    if (a_intersection.UVW.Y.AlmostEquals(0))
        return new Vector2(a_intersection.UVW.X, a_intersection.UVW.Z);
    else
    {
        double u = Math.Atan2(a_intersection.LocalPos.X, 
            -a_intersection.LocalPos.Z) / 2 / Constants.PI;

        if (a_intersection.LocalPos.X < 0)
            u = u + 1;

        return new Vector2(u, a_intersection.UVW.Y);
    }
}

Tekstura wokół osi stożka nawija się od dodatniej osi Z w prawo. No i pozostało nam wyliczanie stycznych do powierzchni potrzebnych do mapowania wypukłości.

public override void GetTangents(
    Intersection a_intersection,
    out Vector3 a_tangent_x, out Vector3 a_tangent_y)
{
    if (a_intersection.LocalPos.Y.AlmostEquals(1))
    {
        a_tangent_x = Right;
        a_tangent_y = Forward;
    }
    else
    {
        if (a_intersection.BackHit)
        {
            a_tangent_x = Vector3.CrossProduct(a_intersection.Normal, 
               Up).Normalized;
            a_tangent_y = Vector3.CrossProduct(a_intersection.Normal, 
                a_tangent_x).Normalized;
        }
        else
        {
            a_tangent_x = Vector3.CrossProduct(Up, 
                a_intersection.Normal).Normalized;
            a_tangent_y = Vector3.CrossProduct(a_tangent_x, 
                a_intersection.Normal).Normalized;
        }
    }
}