Mapowanie tekstury na torus

Z związku z pojawieniem się powierzchni wyższego rzędu przypominających torus dodałem osobną klasę mapowania na torus, której kod bardzo przypomina ten z torusa.

Różnice są na początku. Ponieważ z założenie mapować możemy nie tylko na torus, ale na dowolny kształt nie możemy zakładać, że punkty będą leżeć dokładnie na torusie. Dodatkowo do policzenia mapowania potrzebujemy punkt w lokalnych współrzędnych torusa (tak jak liczyliśmy to dla torusa) zaś tutaj dostajemy go we współrzędnych UVW. Pierwszych kilka linijek liczy więc lokalną współrzędną i mały promień związany z tym punktem na lub poza torusem.

public class TorodidalUVMapper : UVMapper
{
    private double m_local_small_radius;
    private double m_local_big_radius;
    private double m_big_radius;
    private double m_small_radius;

    [YAXNode]
    public double BigRadius
    {
        get
        {
            return m_big_radius;
        }
        set
        {
            m_big_radius = value;
            UpdateLocalRadiuses();
        }
    }

    [YAXNode]
    public double SmallRadius
    {
        get
        {
            return m_small_radius;
        }
        set
        {
            m_small_radius = value;
            UpdateLocalRadiuses();
        }
    }

    [YAXOnDeserialized]
    private void UpdateLocalRadiuses()
    {
        m_local_small_radius = SmallRadius / (SmallRadius + BigRadius);
        m_local_big_radius = BigRadius / (SmallRadius + BigRadius);
    }

    public TorodidalUVMapper(double a_big_radius = 0, double a_small_radius = 0)
    {
        BigRadius = a_big_radius;
        SmallRadius = a_small_radius;
    }

    public override Vector2 Map(Vector3 a_uvw)
    {
        Vector3 local_pos = (a_uvw - Vector3.HALF) / new Vector3(0.5, 0.5 / 
            m_local_small_radius, 0.5);
        Vector3 p1 = new Vector3(local_pos.X, 0, local_pos.Z).Normalized * 
            m_local_big_radius;
        double small_radius = (local_pos - p1).Length;

        double ca = local_pos.Y / small_radius;

        Debug.Assert(ca.IsAlmostRelativeGreaterOrEqualThen(-1));
        Debug.Assert(ca.IsAlmostRelativeLessOrEqualThen(1));
        ca = DoubleExtensions.Limit(ca, -1, 1);

        double v = Math.Acos(ca) / Constants.PI / 2;

        if (local_pos.X * local_pos.X + local_pos.Z * local_pos.Z > 
            m_local_big_radius * m_local_big_radius)
        {
            v = v + 0.25;
        }
        else
        {
            if (local_pos.Y > 0)
                v = 0.25 - v;
            else
                v = 1.25 - v;
        }

        double u = Math.Atan2(local_pos.X, local_pos.Z) / (2 * Constants.PI);

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

        return new Vector2(u, v);
    }

    public override UVMappingMethod Method
    {
        get
        {
            return UVMappingMethod.Toroidal;
        }
    }

    public override string ToString()
    {
        return String.Format("Toroidal uv mapper: R: {0}, r: {1}", BigRadius, 
             SmallRadius);
    }
}

Coraz częste korzystanie z serializacji stanu obiektu

Przypuśćmy, że musimy zserializować stan obiektu, który może mieć różne mapowania UV. Może być ono toroidalne albo stożkowe. Oba wymagają po dwa różne parametry. Wykorzystanie enumeracji do opisania metody mapowania jest niewystarczające. Zamiast tego serializuje cały obiekt. W tym wypadku to jedyne najlepsze rozwiązanie.

Ale np. przez bardzo długi czas lista postprocesorów była budowana w runtime z parametrów opcji renderowania. Wraz z dodaniem postprocesora skalującego, wraz z pojawieniem się filtru jako jego parametru, a filtry mogą mieć różne parametry serializowanie całej listy stało się koniecznością. I tak się dzieje w bardzo wielu miejscach.

Ma to swoją wadę - mniejsza przenośność kodu. Ale na razie póki nie muszę wczytywać innych formatów opisu sceny to rozwiązanie działa bardzo dobrze.

Mapowanie tekstury na stożek

Nazwijmy to mapowanie stożkowym.

Jeśli UVW.Y = 1 to mapowanie przeprowadzamy dla postawy: UV = [UVW.X, UVW.Y]. Przeciwnie mapowanie przeprowadzamy dla powierzchni bocznej stożka.

Jeśli boczną powierzchnię stożka rozwiniemy to otrzymamy wycinek koła. Jeśli wysokość stożka to H, zaś promień jego podstawy to R, to okrąg ten ma promień równy długości tworzącej stożek: $L = \sqrt{R^2 + H^2}$. Obwód stożka u podstawy to $2 \pi R$ i jest to tym samym długość łuku wycinka.

Pierwsze co potrzebujemy to miarę kąta pomiędzy płaszczyzną przechodzącą przez oś stożka i wektor Forward, a płaszczyzną wyznaczoną przez oś stożka i punkt na nim. W moim przypadku wartość tego kąta: $\alpha = \displaystyle\arctan \left ( {\frac{UVW.X-0.5}{-(UVW.Z-0.5)}} \right )$, gdzie $\alpha \subset \left\langle 0,2 \pi\right\rangle$.

Kąt wycinka okręgu, który tworzy bok stożka ma wartość:

$\displaystyle{\frac{\beta}{2 \pi L}=\frac{2 \pi}{2 \pi L}\Rightarrow \beta=\frac{2 \pi R}{L}=\frac{2 \pi R}{\sqrt{R^2 + H^2}}}$

Musimy przeskalować wartość kątu $\alpha$ do zakresu $\left\langle 0, \beta \right\rangle$. Przeskalowana wartość kąta $\alpha$ ma wartość:

$\displaystyle\frac{\alpha}{2 \pi L}=\frac{\alpha_1}{2 \pi R} \Rightarrow \alpha_1=\alpha \frac{R}{L}$

Wartość $\displaystyle S=\frac{R}{L}=\frac{R}{\sqrt{R^2 + H^2}}$ jest stała i można ją wyliczyć na starcie renderowania.

W moim przypadku odjąłem od kąta $\alpha_1$ wartość $\displaystyle\frac{\pi}{2}$ i dodałem $\displaystyle\frac{\beta}{2}=S\pi$ ustawiając początek startu nawijania na [u,v]=[0, 0.5]. Nawijanie odbywa się w prawo.

Mając tak przeliczony kąt:

$\displaystyle\alpha_2=\alpha_1 - \frac{\pi}{2} + S\pi$.

Teraz możemy wyliczyć koordynaty tekstur:

$u = 0.5 + 0.5 \cdot UVW.Y \cdot \cos\left(\alpha_2\right)$
$v = 0.5 + 0.5 \cdot UVW.Y \cdot \sin\left(\alpha_2\right)$

UVW.Y oznacza w tym wypadku odległość punku na powierzchni bocznej stożka od jego wierzchołka i zawiera się ona w granicach od zera do jeden.

Kod:

public class ConicalUVMapper : UVMapper
{
    public readonly double S;

    public ConicalUVMapper(double a_radius, double a_height)
    {
        S = a_radius / Math.Sqrt(a_radius * a_radius + a_height * a_height);
    }

    public override Vector2 Map(Vector3 a_uvw)
    {
        if (a_uvw.Y.AlmostEquals(1))
            return new Vector2(a_uvw.X, a_uvw.Z);
        else
        {
            double alpha = (Constants.PI - Math.Atan2(
               a_uvw.X - 0.5, 0.5 - a_uvw.Z)) * S;
            return new Vector2(0.5 - 0.5 * a_uvw.Y * Math.Cos(alpha), 
                                0.5 - 0.5 * a_uvw.Y * Math.Sin(alpha));
        }
    }
}

Mapowanie UV jako osobne klasy

Postanowiłem główne metody mapowania UVW na UV: sferyczne, cylindryczne, stożkowe i toroidalne przenieść do osobnych klas. Zrobiłem to głównie ze względu by stożek mógł być mapowany zarówno sferycznie jak i stożkowo.

To ostatnie to dosyć głupia nazwa.

Mapowanie można wybrać dla większości obiektów. Niektóre z nich jak stożkowe i toroidalne wymagają dodatkowych parametrów. Niektóre obiekty potrafią uzupełnić te parametry automatycznie: torus dla mapowania toroidalnego, stożek dla stożkowego.

Kodu tych klas nie będę podawać. Jest to wyekstraktowany kod metod GetUV()

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;
        }
    }
}

Cylinder

Wszystko jest bardzo podobne do stożka. Bierzemy kod dla stożka i modyfikujemy go. Uwzględniamy inne równanie powierzchni, normalnej. To, że mamy dwie podstawy. To, że w lokalnym układzie współrzędnych cylinder nie jest obrócony o 180 stopni.

Równanie cylindra w lokalnym układzie współrzędnych:

$X^2+Z^2=R^2$

W naszym przypadku skalujemy stożek odpowiednio tak by H=1 i R = 1.

Podstawiając do $X^2+Z^2=1$ parametryczne równanie linii:

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

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

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

Jego rozwiązanie wyznaczy nam dystans do powierzchni bocznej walca. Interesuje nas minimalne rozwiązanie większe od Constants.MINIMAL_DISTANT dla którego współrzędna Y punktu intersekcji spełnia warunek $Constants.MINIMAL\_DISTANT < Y < 1$.

Normalna walca wyliczona z gradientu ma postać: $[x, 0, z]$

Pełny kod klasy cylindra:

public class CylinderObject : RenderableObject
{
    [YAXNode("Radius")]
    private double m_radius;

    [YAXNode("Height")]
    private double m_height;

    public CylinderObject(Vector3 a_right, Vector3 a_up) :
        base(a_right, a_up)
    {
        Update(UpdateFlags.All);
        Closed = true;
    }

    public Material BottomMaterial
    {
        get
        {
            return Materials[1];
        }
        set
        {
            Materials[1] = value;
        }
    }

    public Material TopMaterial
    {
        get
        {
            return Materials[2];
        }
        set
        {
            Materials[2] = value;
        }
    }

    public double Radius
    {
        get
        {
            return m_radius;
        }
        set
        {
            m_radius = value;
            Update(UpdateFlags.BoundBox | UpdateFlags.Matrices);
        }
    }

    public double Height
    {
        get
        {
            return m_height;
        }
        set
        {
            m_height = value;
            Update(UpdateFlags.BoundBox | UpdateFlags.Matrices);
        }
    }

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

    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;

        Vector3 local_pos = default(Vector3);
        double dist = Double.PositiveInfinity;

        double top_dist = (1 - local_start.Y) / local_dir.Y;
        bool top_hit = false;
        bool bottom_hit = false;

        if (top_dist > Constants.MINIMAL_DISTANT)
        {
            local_pos = local_start + local_dir * top_dist;

            if (local_pos.X * local_pos.X + local_pos.Z * local_pos.Z < 1)
            {
                top_hit = true;
                dist = top_dist;
            }
        }

        double bottom_dist = -local_start.Y / local_dir.Y;

        if ((bottom_dist > Constants.MINIMAL_DISTANT) && (bottom_dist < dist))
        {
            Vector3 bottom_pos = local_start + local_dir * bottom_dist;

            if (bottom_pos.X * bottom_pos.X + bottom_pos.Z * bottom_pos.Z < 1)
            {
                top_hit = false;
                bottom_hit = true;
                dist = bottom_dist;
                local_pos = bottom_pos;
            }
        }

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

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

        if ((roots >= 1) && (dist1 > Constants.MINIMAL_DISTANT) && (dist1 < dist))
        {
            Vector3 local_pos2 = local_start + local_dir * dist1;

            if ((local_pos2.Y >= Constants.MINIMAL_DISTANT) && (local_pos2.Y < 1))
            {
                top_hit = false;
                bottom_hit = false;
                dist = dist1;
                local_pos = local_pos2;
            }
        }

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

            if ((local_pos2.Y >= Constants.MINIMAL_DISTANT) && (local_pos2.Y < 1))
            {
                top_hit = false;
                bottom_hit = false;
                dist = dist2;
                local_pos = local_pos2;
            }
        }

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

        bool back_hit;

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

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

        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 = top_hit ? TopMaterial : bottom_hit ? BottomMaterial : 
                DefaultMaterial
        };
    }

    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.LocalPos.Y.AlmostEquals(0))
        {
            if (a_intersection.BackHit)
                return Up;
            else
                return -Up;
        }
        else
        {
            if (a_intersection.BackHit)
            {
                return (LocalToWorldNormal *
                    new Vector3(-a_intersection.LocalPos.X, 0, 
                        -a_intersection.LocalPos.Z)).Normalized;
            }
            else
            {
                return (LocalToWorldNormal *
                    new Vector3(a_intersection.LocalPos.X, 0, 
                        a_intersection.LocalPos.Z)).Normalized;
            }
        }
    }

    public override string ToString()
    {
        return String.Format("Cone: {0}", Name);
    }

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

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

    public override Vector2 GetUV(Intersection a_intersection)
    {
        if (a_intersection.UVW.Y.AlmostEquals(0) || 
            a_intersection.UVW.Y.AlmostEquals(1))
        {
            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);
        }
    }

    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_intersection.LocalPos.Y.AlmostEquals(0)))
        {
            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 = -Up;
            }
            else
            {
                a_tangent_x = Vector3.CrossProduct(Up, 
                    a_intersection.Normal).Normalized;
                a_tangent_y = Up;
            }
        }
    }

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

    protected override AABB GetLocalBoundBox()
    {
        return new AABB(new Vector3(-Radius, 0, -Radius),
                        new Vector3(Radius, Height, Radius));
    }

    public override void ScaleAbsolute(double a_scale)
    {
        Radius *= a_scale;
        Height *= a_scale;

        base.ScaleAbsolute(a_scale);
    }
}

Więcej materiałów

Przy okazji implementacji stożka zacząłem się zastanawiać jak mając jeden materiał na obiekt poradzić sobie z teksturowaniem go. Prawdopodobnie inaczej będziemy chcieli pokolorować podstawę, a inaczej boki. Ten sam problem pojawi się przy okazji walca. I co ważniejsze przy okazji mesha. Postanowiłem więc zastąpić pojedynczy materiał ich lista, w bardzo wielu przypadkach lista będzie miała jeden materiał. Dodatkowo całe to zarządzanie materiałami wyrzuciłem do osobnego obiektu RenderableObjectMaterials. Wprowadziłem coś takiego jak RenderableObject.DefaultMaterial, który jest pierwszym materiałem na liście. Dodatkowo np. materiały dla stożka

Torus - mapowanie wypukłości

Procedura jest bardzo podobna do tej dla sfery. Wektor X może się wyzerować jeśli r>R i Pos znajduje się w dołku. Możemy taki punkt dla przyjętego mapowania nazwać osobliwym i nie przejmować się mapowaniem w nim.

public override void GetTangents(
    Intersection a_intersection, out Vector3 a_tangent_x,
    out Vector3 a_tangent_y)
{
    a_tangent_x = Vector3.CrossProduct(Up, a_intersection.Pos).Normalized;

    if (a_intersection.BackHit)
    {
        a_tangent_y = Vector3.CrossProduct(a_intersection.Normal, 
            a_tangent_x).Normalized;
    }
    else
    {
        a_tangent_y = Vector3.CrossProduct(a_tangent_x, 
            a_intersection.Normal).Normalized;
    }
}

Tłumienie światła w ośrodku

Promień światła rozchodzący się w danym ośrodku traci swoją energię poprzez pochłanianie i rozproszenie. Możliwa jest także reemisja którą tutaj się nie zajmujemy. Rozproszeniem w sensie pojawianiem się promieni rozproszenia i oświetlaniem przez nie sceny też tutaj się nie zajmujemy. Spróbujemy zasymulować utratę energii przez promień biegnący przez dany ośrodek poprzez jego pochłanianie.

Jeśli założymy że energia promienia na każdy przebyty odcinek dystansu maleje zawsze o tyle samo, niezależnie od wartości tej energii to energia maleje w sposób ekspotencjalny. Inaczej mówiąc poniższe równanie jest rozwiązaniem równania różniczkowego $\displaystyle\frac{dC_{out}}{dx}=a$.

$C_{OUT} = C*e^{-ax}$, gdzie x to dystans, a to współczynnik tłumienia ośrodka.

Latwo zauważyć, że kolor dla $x=0$ nie jest zmieniamy i jego energia maleje do zera wraz z przebytym dystansem. Dla RGB współczynnik a jest trójką liczb.

Inne podejście zastosowano w OpenGL. Tam wartość tłumienia określa się jako:

$T=1/(Ax^2+Bx+C)$

Jeśli chcemy by dla $x=1$ T wynosiło 1 to $C=1$. Jeśli $C>1$ to definiujemy coś co możemy nazwać energią pochłoniętą zaraz na dzień dobry. W sytuacjach odbicia albo refleksji możemy uznać to za energie pochłoniętą podczas tych zjawisk. Ja generalnie energię absorbowaną podczas odbicia albo załamania określam gdzie indziej. Dalej chcemy by $T<1$ dla $X>0$, czyli $A>0$. Poza tym raczej nie chcemy by dla $X>0$ T w jakimkolwiek przedziale malało. Czyli biorąc pod uwagę przesuniętą o -C parabole $Ax^2+Bx$, jedno miejsce zerowe wypadnie w $X=0$, drugie zaś powinno być po ujemnej stronie X, czyli $B>0$. Takie warunki zapewniają nam, że energia z biegiem promienia nigdy nie rośnie.

Powyższym współczynnikom także możemy nadać sens fizyczny. A to energia tracona przez promień na skutek zmniejszającego się strumienia energii wraz ze wzrostem promienia sfery, której środek to punkt emisji. Energia zsumowana po takiej sferze w najlepszym wypadku powinna pozostać bez zmian. B to liniowa utrata energii przez promień, np. na skutek rozproszenia, pochłaniania. Zaś możliwe znaczenie C zostało podane wcześniej.

Ważna cechą obu metod jest to, że jeśli światło przechodzi przez wiele ośrodków, to nie ma znaczenia od której strony liczymy końcową wartość energii jaka dociera do kamery. Innymi słowy możemy liczyć bieg promienia od kamery do światła i po drodze kalkulować efektywny współczynnik tłumienia.

W przypadku gdy nasza scena może znacznie różnić się wymiarami, tzn może być bardzo duża albo bardzo mała we znaki daje się kwestia błędów numerycznych. Najprostszym rozwiązaniem jest normalizacja sceny, a także wyliczanie intersekcji z obiektami w lokalnym znormalizowanym układzie współrzędnych. Kiedy normalizujemy scenę zmieniamy jej wymiary, zmienią się więc odległości przebywane przez promienie i ich kolory - jeśli podlegają one tłumieniu. Musimy wprowadzić dodatkowy czynnik kompensujący. Najlepiej nazwać go prędkością światła. Wraz z zmniejszeniem sceny, chcemy by światło biegło wolniej - wtedy w tej samej jednostce czasu światło przebędzie taką samą drogę jak w scenie powiększonej. Czyli prędkość światła jest proporcjonalna do skali. W naszych równaniach tłumienia X zastąpmy czasem $t=X/V$. Jeśli nasza scena zmalała dwukrotnie, prędkość światła podobnie widzimy, że czas przebycia danego dystansu jest niezmienny.

Takie równania tłumienia są tylko pewnym przybliżeniem rzeczywistości, ale pozwalają dość realistycznie oddać efekt przechodzenia światła przez kolorowe ośrodki, rozpraszanie energii na dalekich odległościach w powietrzu.

Przykład implementacji:

public abstract class Attenuation
{
    public abstract void CalculateAttenuation(float a_light_speed);
    public abstract ColorFloat Attenuate(ColorFloat a_light, float a_dist);
}

public class PolynomialAttenuation : Attenuation
{
    public ColorFloat A;
    public ColorFloat B;
    public ColorFloat C;
    private ColorFloat m_real_a;
    private ColorFloat m_real_b;

    public override ColorFloat Attenuate(ColorFloat a_light, float a_dist)
    {
        return a_light / (m_real_a * a_dist * a_dist + m_real_b * a_dist + C);
    }

    public override void CalculateAttenuation(float a_light_speed)
    {
        Debug.Assert(A >= ColorFloat.Black);
        Debug.Assert(B >= ColorFloat.Black);
        Debug.Assert(C >= ColorFloat.White);

        m_real_a = A / a_light_speed / a_light_speed;
        m_real_b = B / a_light_speed;
    }
}

public class ExponentalAttenuation : Attenuation
{
    private ColorFloat m_real_alpha;
    public ColorFloat Alpha;

    public override void CalculateAttenuation(float a_light_speed)
    {
        m_real_alpha = Alpha / a_light_speed;
    }

    public override ColorFloat Attenuate(ColorFloat a_light, float a_dist)
    {
        return new ColorFloat(
            a_light.R * (float)Math.Exp(-a_dist * m_real_alpha.R),
            a_light.G * (float)Math.Exp(-a_dist * m_real_alpha.G),
            a_light.B * (float)Math.Exp(-a_dist * m_real_alpha.B));
    }
}

Porównywanie liczb z marginesem błędu

W obliczeniach na liczbach zmiennoprzecinkowych w wyniku błędów zaokrąglania w wielu miejscach nie możemy porównywać liczb w sposób absolutny. Np. przyrównywać ich do np. 0 podczas wyliczania wyznaczników wielomianów. W takich sytuacjach pozostaje nam porównanie liczb z pewnym marginesem. Liczby możemy porównywać względnie i bezwzględnie. Np. porównywanie dwóch liczb z dokładnością bezwzględną 1e-12, kiedy różnica pomiędzy nimi może się wahać od 1 do miliona nie ma sensu. Porównywanie liczb bezwzględne ma tylko sens jeśli porównywane liczby mieszczą się w pewnym zakresie - są w jakiś sposób normowane. Oczywiście takie porównywanie jest szybsze niż porównywanie względne.

Oto aktualnie przeze mnie wykorzystywany kod na porównywanie liczb w sposób względny i bezwzględny:

public static bool AlmostEquals(
    this double a_d1, double a_d2, double a_precision)
{
    return Math.Abs(a_d1 - a_d2) < a_precision;
}

public static bool AlmostRelativeEquals(
    this double a_d1, double a_d2, double a_precision)
{
    double mid = Math.Max(Math.Abs(a_d1), Math.Abs(a_d2));
            
 if (Double.IsInfinity(mid))
        return false;
    
    if (mid > a_precision)
        return Math.Abs(a_d1 - a_d2) <= a_precision * mid;
    else
        return a_d1 < a_precision;
}

Dysk

Czyli intersekcja sfery z płaszczyzną. Takie kółko w przestrzeni 3D. Bardzo przydatne do zatykania walców i stożków.

Bardzo prosty do utworzenia element. Zbudowany z niewielkimi modyfikacjami na kodzie płaszczyzny.

Jeśli chodzi o intersekcję promienia z dyskiem, wyliczamy najpierw intersekcję promienia z płaszczyzną w której leży dysk. Następnie liczymy odległość do środka dysku i porównujemy ją z promieniem dysku.

W stosunku do kodu płaszczyzny modyfikujemy jeszcze kod zwracający współrzędne UVW i UV i gotowe.

Przykładowy kod:

public override Intersection GetIntersection(
    Intersection a_source_ray_intersection, Ray a_ray)
{
    Vector3 local_start;
    Vector3 local_dir;
    TransformToLocal(a_ray, out local_dir, out local_start);

    if (OneSide)
    {
        if (local_dir.Y > 0)
            return Scene.NoIntersection;
    }

    if ((a_source_ray_intersection != null) && 
        (a_source_ray_intersection.SceneObject == this))
    {
        return Scene.NoIntersection;
    }

    if (local_dir.Y.AlmostRelativeEquals(0))
        return Scene.NoIntersection;

    double dist = local_start.Y / -local_dir.Y;

    if (dist <= 0)
        return Scene.NoIntersection;
    else
    {
        bool back_hit = (local_dir.Y > 0);

        if (back_hit && OneSide)
            return Scene.NoIntersection;
        else
        {
            Vector3 local_pos = local_start + local_dir * dist;

            if (local_pos.X * local_pos.X + local_pos.Z * local_pos.Z > 1)
                return Scene.NoIntersection;

            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
            };
        }
    }  
}

public override Vector3 GetNormal(Intersection a_intersection)
{
    if (a_intersection.BackHit)
        return -Normal;
    else
        return Normal;
}
public override Vector3 GetUVW(Intersection a_intersection)
{
    return base.GetUVW(a_intersection) * 0.5 + new Vector3(0.5, 0, 0.5);
}

public override Vector2 GetUV(Intersection a_intersection)
{
    Debug.Assert(a_intersection.UVW.Y.AlmostRelativeEquals(0));
    return new Vector2(a_intersection.UVW.X, a_intersection.UVW.Z);
}

public override void GetTangents(Intersection a_intersection,
    out Vector3 a_tangent_x, out Vector3 a_tangent_y)
{
    a_tangent_x = Right;
    a_tangent_y = Forward;
}

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

protected override AABB GetLocalBoundBox()
{
    return new AABB(new Vector3(-Radius, -Constants.MINIMAL_DISTANT, -Radius),
                    new Vector3(Radius, Constants.MINIMAL_DISTANT, Radius));
}

AABB

Czyli Axis-Aligned Bound Box. Czyli prostopadłościan umieszczony w przestrzeni 3D o bokach równoległych do osi układu współrzędnych. Zanim przystąpimy do znajdowania intersekcji dla dowolnie skomplikowanego obiektu, bardziej opłacalne czasowo może być dla nas sprawdzanie intersekcji właśnie z jakąś bryłą w której ten obiekt zawiera się, w tym wypadku AABB. AABB są dosyć często używane gdyż badanie intersekcji z nimi jest dosyć szybkie. Są też naturalnym sojusznikiem dla podziału przestrzeni KD-Tree (zgodnie z płaszczyznami układu współrzędnych). Z uwagi na to, że AABB z reguły niezbyt ciasno opinają obiekt możliwa jest duża liczba fałszywych sprawdzeń. Wyliczanie AABB jest wolne, dla każdej transformacji obiektu trzeba je liczyć na nowo. W przeciwieństwie do np. OBB gdzie wystarczy OBB transformować wraz z obiektem. OOB zreguły opina nasz obiekt bardziej niż AABB. Zamiast mówić opina możemy wprowadzić miarę stosunku objętości bryły wpisanej w AABB do objętości AABB.

Do ich reprezentacji wykorzystujemy dwa punkty w przestrzeni: Minimum i Maksimum, gdzie Min <= Max dla każdej ze współrzędnych X, Y, Z.

Intersekcje z promieniem sprawdzamy osobno dla współrzędnych X, Y, Z. Zacznijmy od X. Sprawdzamy dla jakich t promień dany w postaci kierunkowej przecina płaszczyzny wyznaczone przez Min.X i Max.X.

$\begin{align*}
& x = START.X + DIR.X \cdot t \\
\\
& Min.X = START.X + DIR.X \cdot t \\
\\
& t_{MIN} = \frac{Min.X - START.X}{DIR.X} \\
\\
& t_{MAX} = \frac{Max.X - START.X}{DIR.X} \\
\end{align*}$

Teraz wyznaczamy dwie wartości które będą nam towarzyszyć poprzez cały test:

$\begin{align*}
& Far = Max(t_{MIN}, t_{MAX}) \\
& Near = Min(t_{MIN}, t_{MAX}) \\
\end{align*}$

Jeśli $Far<0$ to przecięcie nastąpiło przed punktem startu promienia. Tym samym intersekcji brak. Zauważmy, że Near może być ujemne.

Jeśli współrzędna kierunkowa $DIR.X=0$ to wtedy intersekcja jest możliwa jeśli $Min.X \le START.X \le Max.X$.

Przechodzimy do następnej płaszczyzny Y. Jeśli $DIR.Y=0$ to sprawdzamy czy START.Y leży w zakresie $Min.Y \le START.Y \le Max.Y$. Przeciwnie wyznaczamy podobnie jak poprzednio $t_{MIN}$ i $t_{MAX}$. Wyznaczamy:

$Far_1 = Max(t_{MIN}, t_{MAX})$
$Near_1 = Min(t_{MIN}, t_{MAX})$

I teraz wyliczamy nowa Near i Far:

$Far = Min(Far, Far_1)$
$Near = Max(Near, Near_1)$

Pamiętajmy, że Far i Near to tak naprawdę parametr t promienia. Wyznaczamy więc dwukrotnie wartości parametru t dla którego prosta przecina dane płaszczyzny. Mamy dwie wartości Near i Far z testu dla X i teraz mamy dwie nowe wartości z testy dla Y. Jeśli prosta przetnie AABB to dokona tego dla konkretnych dwóch wartości t. Podstawiając je do wzorów na Near i Far dla konkretnych płaszczyzn powinniśmy otrzymać współrzędne Min i Max. W zakresie t pomiędzy punktem wejścia w AABB, a wyjścia współrzędne punktów prostej leżą pomiędzy Min i Max. Powyższe wyliczenia znajdują nowy wspólny zbiór parametrów t dla których prosta znajduje się jednocześnie pomiędzy Min.X i Max.X oraz pomiędzy Min.Y, a Max.X.

Jeśli $Far<0$ to intersekcji brak. Zbiór t dla której są one pomiędzy Min.Y, a Max.Y leży przed punktem startu promienia.

Powyższy rysunek obrazuje różne kombinacje możliwych intersekcji:

Dla promienia 2 współrzędna kierunkowa DIR.Y. Ponieważ nie zachodzi $Min.Y \le START.Y \le Max.Y$ prosta ta nigdy nie przetnie AABB. W przeciwieństwie do promienia 5. Dla promienia 1 zbiory $<Near.X, Far.X>$ oraz $<Near.Y, Far.Y>$ są rozłączne. Promień ten nigdy nie przetnie AABB. W przeciwieństwie do promienia 4, gdzie zbiory te mają część wspólną. Promień 3 jest przykładem promienia dla którego $Near<0$ utrzymuje się przez cały czas obliczeń. I tak ujemna wartość zostanie zwrócona jako odległość do miejsca intersekcji. Wartość oznaczająca brak intersekcji to Double.PositiveInfinity. Wartość ujemna oznacza, że punkt startu promienia znajduje się w AABB. Taki wynik oznacza, że powinniśmy także przetestować na intersekcje wszystkie obiekty w tym AABB.

Jeśli $Far<Near$ to intersekcji brak. Zbiory t są rozłączne, innymi słowy prosta nigdy jednocześnie nie znajduje się pomiędzy Min.X i Max.X oraz Min.Y i Max.Y i tym samym nie może przeciąć AABB.

Podobnie jak dla Y czynimy dla Z.

Ostatecznie otrzymuje segment promienia który przecina AABB, z reguły interesuje nas tylko współrzędna wejścia do AABB.

Kod wyznaczania dystansu do punktu pierwszego punktu intersekcji:

public double GetDist(Ray a_ray)
{
    Debug.Assert(Minimum <= Maximum);

    double n = Double.MinValue;
    double f = Double.MaxValue;

    if (a_ray.Dir.X.AlmostRelativeEquals(0))
    {
        if ((a_ray.Start.X <= Minimum.X) || (a_ray.Start.X >= Maximum.X))
            return Double.PositiveInfinity;
    }
    else
    {
        double d1 = (Minimum.X - a_ray.Start.X) / a_ray.Dir.X;
        double d2 = (Maximum.X - a_ray.Start.X) / a_ray.Dir.X;

        f = Math.Max(d1, d2);

        if (f < 0)
            return double.PositiveInfinity;

        n = Math.Min(d1, d2);
    }

    if (a_ray.Dir.Y.AlmostRelativeEquals(0))
    {
        if ((a_ray.Start.Y <= Minimum.Y) || (a_ray.Start.Y >= Maximum.Y))
            return Double.PositiveInfinity;
    }
    else
    {
        double d1 = (Minimum.Y - a_ray.Start.Y) / a_ray.Dir.Y;
        double d2 = (Maximum.Y - a_ray.Start.Y) / a_ray.Dir.Y;

        f = Math.Min(f, Math.Max(d1, d2));

        if (f < 0)
            return double.PositiveInfinity;

        n = Math.Max(n, Math.Min(d1, d2));

        if (n > f)
            return Double.PositiveInfinity;
    }

    if (a_ray.Dir.Z.AlmostRelativeEquals(0))
    {
        if ((a_ray.Start.Z <= Minimum.Z) || (a_ray.Start.Z >= Maximum.Z))
            return Double.PositiveInfinity;
    }
    else
    {
        double d1 = (Minimum.Z - a_ray.Start.Z) / a_ray.Dir.Z;
        double d2 = (Maximum.Z - a_ray.Start.Z) / a_ray.Dir.Z;

        f = Math.Min(f, Math.Max(d1, d2));

        if (f < 0)
            return double.PositiveInfinity;

        n = Math.Max(n, Math.Min(d1, d2));

        if (n > f)
            return Double.PositiveInfinity;
    }

    return n;
}

Torus - mapowanie tekstury

Mamy punkt uderzenia promienia w torus $P(X,Y,Z)$ w lokalnym układzie współrzędnych, w którym torus leży w środku układu współrzędnych w płaszczyźnie XZ.

Współrzędna u jest liczona całkiem podobnie jak dla sfery $u=\arctan(^X/_Z)/2\pi$. Teraz dla drugiej połowy okręgu współrzędną trzeba zwiększyć o jeden. W moim przypadku przyjąłem, że współrzędna u=0 leży na dodatniej osi X.

Ze współrzędną v jest trochę trudniej. Zauważmy, że niezależnie od położenia punktu na torusie w płaszczyźnie go przecinającej przechodzącej przez oś Y kątowa wartość współrzędnej v ma wartość: $\arccos(^Y/_r)$. Teraz tylko musimy ją odpowiednio przekształcić w zależności w której ćwiartce znalazł się punkt uderzenia. Ja przyjąłem, że współrzędna $v=0$ jest wewnątrz torusa, $v=0.25$ na górze (góra w kierunki osi Y) itd.

Znajdująca się w kodzie korekta wartości argumentu dla arccos to korekta błędów numerycznych. Oczywiście jeśli wartość jest większa od jeden albo mniejsza od -1 o więcej niż wynosi dopuszczalny błąd to zdecydowanie jest coś nie tak.

Kod:

public override Vector2 GetUV(Intersection a_intersection)
{
    Vector3 local_pos = a_intersection.LocalPos;

    float ca = local_pos.Y / m_local_small_radius;
    Debug.Assert(ca.IsAlmostRelativeGreaterOrEqualThen(-1));
    Debug.Assert(ca.IsAlmostRelativeLessOrEqualThen(1));
    ca = SingleExtensions.Limit(ca, -1, 1);

    float v = (float)Math.Acos(ca) / Constants.PI / 2;

    if (local_pos.X * local_pos.X + local_pos.Z * local_pos.Z > 
        m_local_big_radius * m_local_big_radius)
    {
        v = v + 0.25f;
    }
    else
    {
        if (local_pos.Y > 0)
            v = 0.25f - v;
        else
            v = 1.25f - v;
    }

    float u = (float)Math.Atan2(local_pos.X, local_pos.Z) / (2 * Constants.PI);

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

    return new Vector2(u, v);
}

Torus - normalna

Normalną możemy policzyć korzystając z gradientu funkcji lub poprzez podejście geometryczne. Podejście geometryczne daje poprawne wyniki kiedy r>R i kiedy R=0.

Równanie torusa leżącego w płaszczyźnie XZ:

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

Normalna to gradient funkcji:

$f(x,y,z)=(\sqrt{x^2+y^2}-R)^2+y^2-r^2$
$\nabla f=\left[2x-\frac{2xR}{\sqrt{x^2+z^2}},2y,2z-\frac{2zR}{\sqrt{x^2+z^2}}\right]$

Możemy także najpierw lekko uporządkować funkcje:

$(\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)$
$f(x,y,z)=(x^2+y^2+z^2+R^2-r^2)^2-4R^2(x^2+z^2)$
$\begin{align*}
\nabla f=[& 4x(x^2+y^2+z^2+R^2-r^2)^2-8xR^2, \\& 4y(x^2+y^2+z^2+R^2-r^2)^2, \\& 4z(x^2+y^2+z^2+R^2-r^2)^2-8zR^2]
\end{align*}$

W podejściu geometrycznym mamy punkt w który uderzył promień należący do torusa $P(x,y,z)$. Torus leży w płaszczyźnie $XZ$. Rzut wektora $\vec P$ (od środka torusa do punktu P) na płaszczyznę $XZ$ to $\vec A=[x,y,0]$. Wektor o długości $R$ o takim samym kierunku to $\vec B = R\hat A$. Wektor $\hat{\vec P-\vec B}$ to normalna. Jeśli mamy uderzenie od środka torusa to należy ją wziąć ze znakiem minus.

Kod wszystkich trzech wersji:

public override Vector3 GetNormal(
   Intersection a_intersection)
{
    Vector3 local_pos = a_intersection.LocalPos;
            
    //float R2 = m_local_big_radius * m_local_big_radius;
    //float n = local_pos.SqrLen + R2 - m_local_small_radius * m_local_small_radius;
    //        
    //Vector3 normal = new Vector3(
    //    local_pos.X * n - 2 * R2 * local_pos.X,
    //    local_pos.Y * n,
    //    local_pos.Z * n - 2 * R2 * local_pos.Z);
    //
    //if (a_intersection.BackHit)
    //    return (LocalToWorldNormal * -normal).Normalized;
    //else
    //    return (LocalToWorldNormal * normal).Normalized;

    if (a_intersection.BackHit)
    {
        return (LocalToWorldNormal * (new Vector3(local_pos.X, 0, 
                local_pos.Z).Normalized * m_local_big_radius - 
                local_pos)).Normalized;
    }
    else
    {
        return (LocalToWorldNormal * (local_pos - new Vector3(local_pos.X, 0, 
                local_pos.Z).Normalized * m_local_big_radius)).Normalized;
    }

    //var nn = LocalToWorldNormal * new Vector3(
    //    local_pos.X - local_pos.X * m_local_big_radius / 
    //    (float)Math.Sqrt(local_pos.X * local_pos.X + local_pos.Z * local_pos.Z),
    //    local_pos.Y,
    //    local_pos.Z - local_pos.Z * m_local_big_radius / 
    //    (float)Math.Sqrt(local_pos.X * local_pos.X + local_pos.Z * local_pos.Z)
    //    ).Normalized;
    //if (a_intersection.BackHit)
    //    return -nn;
    //else
    //    return nn;
}

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;
}

Przeźroczyste cienie

Wraz z dodaniem materiałów przeźroczystych pojawiła się potrzeba dodania czegoś co można by nazwać półprzeźroczystymi cieniami. Na takim materiale promień światła ulega załamaniu. W obecnym stanie raytracera, kiedy nie ma globalnego modelu oświetlenia nie możemy tego efektu uwzględnić podczas obliczania czy punkt jest w cieniu światła. Musimy więc ucieć się do triku. Wypuszczamy z miejsca uderzenia promień w kierunku źródła światła i badamy wszystkie intersekcje po drodze. W przypadku natrafienia na materiał nieprzepuszczalny punkt jest w cieniu dla tego żródła światła. Na każdym materiale półprzeźroczystym dokonujemy modyfikacji koloru światła. Lolor światła modyfikujemy także w zależności od tłumienia ośrodka. Obie modyfikacje koloru światła są niezależne od siebie oraz są niezależne od kolejności w jakich ich dokonujemy. Możemy więc poruszać się od miejsca uderzenia do źródła światła (jak i w druga stronę). Promień na materiale półprzeźroczystym nie ulega załamaniu.

Co możemy zrobić w przypadku uderzenia w materiał półprzeźroczysty. Możemy pomniejszyć kolor światła źródła o transmitancje materiału. Możemy pomniejszyć go o wartość kosinusa kąta pomiędzy normalną, a kierunkiem promienia symulując tym samym ubytek energii promienia w wyniku rozproszenia. Możemy wreszcie uwzględnić efekt Fresnella.

Wszystkie wspomniane metody zostały zaimplementowane jako osobne klasy wywodzące się z ShadowTester. Jest to specjalna klasa do ustalania koloru światła w punkcie uderzenia.

Pierwszy najprostszy tester SolidShadowTester nie implementuje przeźroczystych cieni. Drugi AttenuationShadowTester uwzględnia tłumienie ośrodka. Trzeci AttenuationTransmittanceShadowTester dodatkowo uwzględnia trasmitancje materiału. Czwarty AttenuationTransmittanceAngleShadowTester bierze pod uwagę kąt padania promienia do normalnej. Przy czym dla bardziej realistycznego efektu podnoszony jest on do potęgi Constants.TRANSFLUENT_SHADOWS_ANGLE_POWER_EXPONENT=0.25. Ponieważ wartości kąta zawierają się od 0 do 1, wyginamy więc krzywą tak, że $cos^{0.25}>cos$. Cienie są więc o wiele bardziej jasne, i w ten sposób są bardziej realistyczne. Ostatni tester AttenuationTransmittanceAngleFresnellShadowTester w warunkach lokalnego oświetlenia okazał się najgorszy. Potrafi wygenerować całkowicie ciemne cienie w wyniku całkowitego wewnętrznego odbicia. Ale zostawiłem go mimo wszystko dla porównania.

Klasy promieni

Wraz z rozwojem raytracera kod promienia musiał się zmienić i to znacznie. No dobra może nie musiał aż tak.

Po pierwsze intersekcje zacząłem sprawdzać w lokalnym układzie współrzędnych.

Po drugie testowanie cieni półprzeźroczystych wymagało wielu testów dla promienia podążającego w tym samym kierunku, tak że zmienia się tylko punkt startu dla kolejnych powierzchni półprzeźroczystych. Pokusa manipulowania danymi promienia okazała się duża. Z drugiej strony promień jest wykorzystywany w wielu miejscach. Po takiej manipulacji i przed powrotem z funkcji trzeba go przywrócić do stanu pierwotnego, nie wspominając o innych możliwych niepożądanych interakcjach, kiedy zmodyfikujemy promień z klasy intersekcji, a później chcemy coś z niej liczyć - wiele elementów tej klasy jest liczona w sposób opóźniony na żądanie. Tak więc modyfikacje promieni musiały zostać zakazane, a promień musiał stać się klasą immutable. Tzn. poprzednio nie był immutable, ale raczej go nie modyfikowałem. Teraz zostało to zakazane jawnie.

Wyleciała z niego przede wszystkim właściwość wskazująca na źródłową intersekcję. Została ona dodana jako parametr wielu funkcji. Wraz z tym z promienia poleciało wiele funkcji, których wartość wymagała do obliczenia intersekcji. Okazało się, że kod dało się bez problemu przekształcić, tak by się bez nich obejść, ale w kodzie promienia musiała się pojawić jedna dodatkowa informacja mówiąca nam o tym z jakim promieniem mamy do czynienia: odbitym czy załamanym. Bez tej informacji wiele funkcji nie może działać.

Stary kod promienia:

public class Ray
{
    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
    private Vector3 m_dir;

    public int Depth = 0;
    public Vector3 Start;
    public Intersection SourceIntersection;
    public bool Refracted;

    public RenderableObject PrevHit
    {
        get
        {
            if (SourceIntersection != null)
                return SourceIntersection.SceneObject;
            else
                return null;
        }
    }

    public bool PrevBackHit
    {
        get
        {
            if (SourceIntersection != null)
                return SourceIntersection.BackHit;
            else
                return false;
        }
    }

    public Vector3 Dir
    {
        get
        {
            return m_dir;
        }
        set
        {
            Debug.Assert(value.IsNormalized);
            m_dir = value;
        }
    }

    public Vector3 HitPoint(float a_distance)
    {
        return Start + Dir * a_distance;
    }

    public override string ToString()
    {
        return String.Format("start: {0}; direction: {1}; depth: {2}",
            Start, Dir, Depth);
    }

    public void Transform(Matrix4 a_matrix, out Vector3 a_local_dir, out Vector3 a_local_start)
    {
        a_local_dir = (a_matrix.Rotation * Dir).Normalized;
        a_local_start = a_matrix * Start;
    }
}

Nowy kod:

public enum RaySurfaceSide
{
    /// 

    /// i.e. Camera ray, LightTestRay.
    /// 
    NonDetermined, 

    /// 

    /// i.e. RefractedRay. ShadowRay.
    /// 
    RefractedSide,

    /// 

    /// i.e. ReflectedRay, ShadowRay.
    /// 
    ReflectedSide
}

public abstract class Ray
{
    [DebuggerBrowsable(DebuggerBrowsableState.Never)]
    private Vector3 m_dir;

    public int Depth { get; protected set; }
    public Vector3 Start { get; protected set; }
    public RaySurfaceSide RaySurfaceSide { get; protected set; }

    protected Ray(Scene a_scene)
    {
        a_scene.Statistics.RaysCreated++;
    }

    public Vector3 Dir
    {
        get
        {
            return m_dir;
        }
        protected set
        {
            Debug.Assert(value.IsNormalized);
            m_dir = value;
        }
    }

    public Vector3 HitPoint(float a_distance)
    {
        return Start + Dir * a_distance;
    }

    public override string ToString()
    {
        return String.Format("{3}, Start: {0}; Dir: {1}; Depth: {2}, Type: {4}",
            Start, Dir, Depth, GetType().Name, RaySurfaceSide);
    }

    public void Transform(Matrix4 a_matrix, out Vector3 a_local_dir, out Vector3 a_local_start)
    {
        a_local_dir = (a_matrix.Rotation * Dir).Normalized;
        a_local_start = a_matrix * Start;
    }
}