A SIMPLE RAYTRACER

<pragma:references('System.Drawing.dll')>
<pragma:target('exe')>

using System, System.Reflection;
using System.Drawing;

namespace Freya.Demos.RayTracing;

public

    // GEOMETRY ----------------------------------------

    Tolerance = static class
    public
        const Epsilon: Double = 0.000001;
    end;

    Vector = record
    public
        X, Y, Z: Double;
        constructor(X, Y, Z: Double);

        property Length: Double =>
            (X.Sqr + Y.Sqr + Z.Sqr).Sqrt;

        property Sqr: Double =>
            X.Sqr + Y.Sqr + Z.Sqr;

        method Dot(V: Vector): Double;
        method Mirror(Axis: Vector): Vector;
        method Normalize: Vector;
        method ToString: String; override;

    public
        static Null: Vector = new Vector(0, 0, 0); readonly;
        static XRay: Vector = new Vector(1, 0, 0); readonly;
        static YRay: Vector = new Vector(0, 1, 0); readonly;
        static ZRay: Vector = new Vector(0, 0, 1); readonly;

        method+(V1, V2: Vector): Vector;
        method-(V1, V2: Vector): Vector;
        method*(V1, V2: Vector): Double;
        method^(V1, V2: Vector): Vector;
        method*(Factor: Double; V: Vector): Vector;
    end;

    Matrix = record
    public
        A11, A12, A13: Double;
        A21, A22, A23: Double;
        A31, A32, A33: Double;

        constructor(A11, A12, A13, A21, A22, A23, A31, A32, A33: Double);
        constructor(Eye, LookAt, Sky: Vector);

        method Rotate(V: Vector): Vector;
        method Rotate(X, Y, Z: Double): Vector;
        method RotateNormalize(X, Y, Z: Double): Vector;
        method Transpose;
    end;

    <DefaultMember('Items')>
    Ray = sealed class
    public
        Origin, Direction: Vector;

        property Items[Time: Double]: Vector;
        begin
            Result.X := Origin.X + Time * Direction.X;
            Result.Y := Origin.Y + Time * Direction.Y;
            Result.Z := Origin.Z + Time * Direction.Z;
        end;
    end;

    HitInfo = record
    public
        HitPoint, Normal: Vector;
        Material: IMaterial;
        Time: Double;
    end;

    // COLOR MODEL -------------------------------------

    Pixel = record
    public
        R, G, B: Single;
        constructor(R, G, B: Single);
        constructor(R, G, B: Double);
        constructor(Color: System.Drawing.Color);
        constructor(Intensity: Double);

        method Add(F: Single; P: Pixel);
        method Clip;
        method Interpolate(Amount: Single; Other: Pixel): Pixel;
        method Multiply(F: Single; P: Pixel): Pixel;
        method ToColor: System.Drawing.Color;
        method ToString: String; override;

    public
        static Black: Pixel = new Pixel; readonly;
        static White: Pixel = new Pixel(1, 1, 1); readonly;
        static Red:   Pixel = new Pixel(1, 0, 0); readonly;
        static Green: Pixel = new Pixel(0, 1, 0); readonly;
        static Blue:  Pixel = new Pixel(0, 0, 1); readonly;

        method+(P1, P2: Pixel): Pixel;
        method*(P1, P2: Pixel): Pixel;
        method*(Factor: Single; P: Pixel): Pixel;
    end;

    <DefaultMember('Items')>
    PixelMap = sealed class
    public
        constructor(Height, Width: Integer);
        constructor(Colors: Array@2[Pixel]);

        property Width, Height: Integer; readonly;
        property Items[Row, Col: Integer]: Pixel;

        method ToBitmap: System.Drawing.Bitmap;
    end;

    // BASIC INTERFACES --------------------------------

    IShape = interface
        method ShadowTest(R: Ray): Boolean;
        method HitTest(R: Ray; MaxTime: Double; var Info: HitInfo): Boolean;
    end;

    ISampler = interface
        method Render(Scene: Scene): PixelMap;
        method Render(Scene: Scene; Row, Col: Integer): Pixel;
    end;

    ICamera = interface
        property PrimaryRay: Ray; readonly;
        property Location, Target, Up: Vector;
        property Height, Width: Integer;

        method GetRay(Row, Col: Integer);
    end;

    ILight = interface
        property Color: Pixel; readonly;

        method GetDirection(L: Vector): Vector;
        method Initialize(Scene: Scene);
        method Intensity(L: Vector): Single;
    end;

    IMaterial = interface
        method GetColor(Location: Vector): Pixel;
        method Radiance(
            Location, Normal, Reflection: Vector;
            Light: ILight; Color: Pixel): Pixel;
        method Reflection(
            var Location: Vector; Cosine: Double; out Filter: Pixel): Single;
    end;

    Scene = sealed class
    public
        constructor(
            Sampler: ISampler; Camera: ICamera; Root: IShape;
            Lights: Array[ILight]);

        property Sampler: ISampler; readonly;
        property Camera: ICamera; readonly;
        property Root: IShape; readonly;
        property Lights: Array[ILight]; readonly;

        method Render: PixelMap;
        method Render(Row, Col: Integer): Pixel;

        static method Render(Sampler: ISampler; Camera: ICamera; Root: IShape;
            Lights: Array[ILight]): PixelMap;
    end;

    // SCENE IMPLEMENTATIONS ----------------------------

    BasicSampler = sealed class(ISampler)
    public
        constructor(Bounces: Integer; MinWeight: Double);

        property Bounces: Integer; readonly;
        property MinWeight: Double; readonly;

        method Render(Scene: Scene): PixelMap;
        method Render(Scene: Scene; Row, Col: Integer): Pixel;
    protected
        method Trace(R: Ray): Pixel; virtual;
    end;

    PerspectiveCamera = sealed class(ICamera)
    public
        constructor(
            Location, Target, Up: Vector;
            Angle: Double; Height, Width: Integer);
        constructor(
            Location, Target: Vector;
            Angle: Double; Height, Width: Integer);
        constructor(
            Location, Target: Vector;
            Height, Width: Integer);

        property PrimaryRay: Ray; readonly;
        property Location, Target, Up: Vector;
        property Height, Width: Integer;

        method GetRay(Row, Col: Integer);
    end;

    PointLight = sealed class(ILight)
    public
        constructor(Location: Vector; C: Pixel);

        property Color: Pixel; readonly;

        method GetDirection(L: Vector): Vector;
        method Initialize(Scene: Scene);
        method Intensity(L: Vector): Single;
    end;

    // SHAPES IMPLEMENTATIONS ---------------------------

    Sphere = sealed class(IShape)
    public
        constructor(Center: Vector; Radius: Double; Material: IMaterial);
    end;

    Cylinder = sealed class(IShape)
    public
        constructor(Bottom, Top: Vector; Radius: Double; Material: IMaterial);
    end;

    Union = sealed class(IShape)
    public
        constructor(Shapes: Array[IShape]);
        constructor(Shape1, Shape2: IShape);
        constructor(Shape1, Shape2, Shape3: IShape);
        constructor(Shape1, Shape2, Shape3, Shape4: IShape);
    end;

    // MATERIALS ----------------------------------

    Plastic = sealed class(IMaterial)
    public
        constructor(Color: Pixel; Reflection, PhongAmount, PhongSize: Double);
        constructor(Color: Pixel; Reflection: Double);
        constructor(Color: Pixel);
    end;

    Metal = sealed class(IMaterial)
    public
        constructor(Color: Pixel;
            MinReflection, MaxReflection, Diffuse: Double;
            PhongAmount, PhongSize: Double);
        constructor(Color: Pixel;
            MinReflection, MaxReflection, Diffuse: Double);
        constructor(Color: Pixel;
            MinReflection, MaxReflection: Double);
        constructor(Color: Pixel; MinReflection: Double);
    end;

implementation for Vector is

    constructor(X, Y, Z: Double);
    begin
        Self.X := X;
        Self.Y := Y;
        Self.Z := Z;
    end;

    method Dot(V: Vector): Double =>
        X * V.X + Y * V.Y + Z * V.Z;

    method Mirror(Axis: Vector): Vector =>
        (Self - (2.0 * Dot(Axis)) * Axis).Normalize;

    method Normalize: Vector =>
        using len := (X.Sqr + Y.Sqr + Z.Sqr).Sqrt do
            if len = 0 then Self
            else new Vector(X / len, Y / len, Z / len);

    method ToString: String =>
        String.Format('[{0},{1},{2}]', X, Y, Z);

    method+(V1, V2: Vector): Vector;
    begin
        Result.X := V1.X + V2.X;
        Result.Y := V1.Y + V2.Y;
        Result.Z := V1.Z + V2.Z;
    end;

    method-(V1, V2: Vector): Vector;
    begin
        Result.X := V1.X - V2.X;
        Result.Y := V1.Y - V2.Y;
        Result.Z := V1.Z - V2.Z;
    end;

    method*(V1, V2: Vector): Double =>
        V1.X * V2.X + V1.Y * V2.Y + V1.Z * V2.Z;

    method^(V1, V2: Vector): Vector =>
        new Vector(
            V1.Y * V2.Z - V1.Z * V2.Y,
            V1.Z * V2.X - V1.X * V2.Z,
            V1.X * V2.Y - V1.Y * V2.X);

    method*(Factor: Double; V: Vector): Vector;
    begin
        Result.X := Factor * V.X;
        Result.Y := Factor * V.Y;
        Result.Z := Factor * V.Z;
    end;

implementation for Matrix is

    constructor(A11, A12, A13, A21, A22, A23, A31, A32, A33: Double);
    begin
        Self.A11 := A11; Self.A12 := A12; Self.A13 := A13;
        Self.A21 := A21; Self.A22 := A22; Self.A23 := A23;
        Self.A31 := A31; Self.A32 := A32; Self.A33 := A33;
    end;

    constructor(Eye, LookAt, Sky: Vector);
    begin
        A13 := LookAt.X - Eye.X;
        A23 := LookAt.Y - Eye.Y;
        A33 := LookAt.Z - Eye.Z;
        var dist := 1.0 / (A13.Sqr + A23.Sqr + A33.Sqr).Sqrt;
        A13 *= dist; A23 *= dist; A33 *= dist;
        A11 := Sky.Y * A33 - Sky.Z * A23;
        A21 := Sky.Z * A13 - Sky.X * A33;
        A31 := Sky.X * A23 - Sky.Y * A13;
        dist := 1.0 / (A11.Sqr + A21.Sqr + A31.Sqr).Sqrt;
        A11 *= dist; A21 *= dist; A31 *= dist;
        A12 := A23 * A31 - A33 * A21;
        A22 := A33 * A11 - A13 * A31;
        A32 := A13 * A21 - A23 * A11;
    end;

    method Rotate(V: Vector): Vector =>
        new Vector(
            A11 * V.X + A12 * V.Y + A13 * V.Z,
            A21 * V.X + A22 * V.Y + A23 * V.Z,
            A31 * V.X + A32 * V.Y + A33 * V.Z);

    method Rotate(X, Y, Z: Double): Vector =>
        new Vector(
            A11 * X + A12 * Y + A13 * Z,
            A21 * X + A22 * Y + A23 * Z,
            A31 * X + A32 * Y + A33 * Z);

    method RotateNormalize(X, Y, Z: Double): Vector =>
        (new Vector(
            A11 * X + A12 * Y + A13 * Z,
            A21 * X + A22 * Y + A23 * Z,
            A31 * X + A32 * Y + A33 * Z)).Normalize;

    method Transpose;
    begin
        var T := A12; A12 := A21; A21 := T;
        T := A13; A13 := A31; A31 := T;
        T := A23; A23 := A32; A32 := T;
    end;

implementation for Pixel is

    constructor(R, G, B: Single);
    begin
        Self.R := Math.Min(R, 1.0F);
        Self.G := Math.Min(G, 1.0F);
        Self.B := Math.Min(B, 1.0F);
    end;

    constructor(R, G, B: Double);
    begin
        Self.R := Single(Math.Min(R, 1.0));
        Self.G := Single(Math.Min(G, 1.0));
        Self.B := Single(Math.Min(B, 1.0));
    end;

    constructor(Color: System.Drawing.Color);
    begin
        Self.R := Color.R / 255F;
        Self.G := Color.G / 255F;
        Self.B := Color.B / 255F;
    end;

    constructor(Intensity: Double);
    begin
        Self.R := Single(Math.Min(Intensity, 1.0));
        Self.G := Self.R;
        Self.B := Self.R;
    end;

    method Add(F: Single; P: Pixel);
    begin
        R += F * P.R;
        G += F * P.G;
        B += F * P.B;
    end;

    method Clip;
    begin
        if R < 0.0F then R := 0.0F else if R > 1.0F then R := 1.0F;
        if G < 0.0F then G := 0.0F else if G > 1.0F then G := 1.0F;
        if B < 0.0F then B := 0.0F else if B > 1.0F then B := 1.0F;
    end;

    method Interpolate(Amount: Single; Other: Pixel): Pixel =>
        new Pixel(
            R + (Other.R - R) * Amount,
            G + (Other.G - G) * Amount,
            B + (Other.B - B) * Amount);

    method Multiply(F: Single; P: Pixel): Pixel;
    begin
        Result.R := F * R * P.R;
        Result.G := F * G * P.G;
        Result.B := F * B * P.B;
    end;

    method ToColor: System.Drawing.Color =>
        System.Drawing.Color.FromArgb(255,
            Integer(R * 255F), Integer(G * 255F), Integer(B * 255F));

    method ToString: String =>
        String.Format('<{0}/{1}/{2}>', R, G, B);

    method+(P1, P2: Pixel): Pixel;
    begin
        Result.R := P1.R + P2.R;
        Result.G := P1.G + P2.G;
        Result.B := P1.B + P2.B;
    end;

    method*(P1, P2: Pixel): Pixel;
    begin
        Result.R := P1.R * P2.R;
        Result.G := P1.G * P2.G;
        Result.B := P1.B * P2.B;
    end;

    method*(Factor: Single; P: Pixel): Pixel;
    begin
        Result.R := Factor * P.R;
        Result.G := Factor * P.G;
        Result.B := Factor * P.B;
    end;

implementation for PixelMap is

    Colors: Array@2[Pixel];

    constructor(Height, Width: Integer);
    begin
        Self.Height := Height;
        Self.Width := Width;
        Self.Colors := new Array@2[Pixel](Height, Width);
    end;

    constructor(Colors: Array@2[Pixel]);
    begin
        Self.Height := Colors.GetUpperBound(0) + 1;
        Self.Width := Colors.GetUpperBound(1) + 1;
        Self.Colors := Colors;
    end;

    property Items(Row, Col: Integer): Pixel;
    begin
        Result := Colors[Row, Col];
    end;

    property Items(Row, Col: Integer; Value: Pixel);
    begin
        Colors[Row, Col] := Value;
    end;

    method ToBitmap: System.Drawing.Bitmap;
    begin
        Result := new System.Drawing.Bitmap(Width, Height);
        try
            for row := 0 to Height - 1 do
                for col := 0 to Width - 1 do
                    Result.SetPixel(col, row,
                        Colors[Height - row - 1, col].ToColor);
        except
            Result.Dispose;
            raise;
        end;
    end;

implementation for Scene is

    constructor(
        Sampler: ISampler; Camera: ICamera; Root: IShape;
        Lights: Array[ILight]);
    begin
        Self.Sampler := Sampler;
        Self.Camera := Camera;
        Self.Root := Root;
        Self.Lights := Lights ?? new Array[ILight](0);
    end;

    method Render: PixelMap => Sampler.Render(Self);

    method Render(Sampler: ISampler; Camera: ICamera; Root: IShape;
        Lights: Array[ILight]): PixelMap =>
        (new Scene(Sampler, Camera, Root, Lights)).Render;

    method Render(Row, Col: Integer): Pixel =>
        Sampler.Render(Self, Row, Col);

implementation for BasicSampler is

    Root: IShape;
    Lights: Array[ILight];

    constructor(Bounces: Integer; MinWeight: Double);
    begin
        Self.Bounces := Bounces;
        Self.MinWeight := MinWeight;
    end;

    method Initialize(Scene: Scene);
    begin
        Self.Root := Scene.Root;
        Self.Lights := Scene.Lights;
        for light in Self.Lights do
            light.Initialize(Scene);
    end;

    method Render(Scene: Scene): PixelMap;
    begin
        Initialize(Scene);
        var Camera := Scene.Camera;
        var Wdth := Camera.Width;
        var Pixels := new Array@2[Pixel](Camera.Height, Wdth);
        for Row: Integer := 0 to Camera.Height - 1 do
        begin
            for Col: Integer := Wdth - 1 downto 0 do
            begin
                Camera.GetRay(Row, Col);
                Pixels[Row, Col] := Trace(Camera.PrimaryRay);
            end;
        end;
        Result := new PixelMap(Pixels);
    end;

    method Render(Scene: Scene; Row, Col: Integer): Pixel;
    begin
        Initialize(Scene);
        Scene.Camera.GetRay(Row, Col);
        Result := Trace(Scene.Camera.PrimaryRay);
    end;

    method Trace(R: Ray): Pixel;
    var
        info: HitInfo;
        lightFilter: Pixel;
        reflection: Vector;
    begin
        var bncs := Self.Bounces;
        var weight := 1.0F;
        var filter := Pixel.White;
        while Root.HitTest(R, Double.MaxValue, var info) do
        begin
            var location := R[info.Time];
            var cosine := R.Direction.Dot(info.Normal);
            reflection.X := R.Direction.X - 2 * cosine * info.Normal.X;
            reflection.Y := R.Direction.Y - 2 * cosine * info.Normal.Y;
            reflection.Z := R.Direction.Z - 2 * cosine * info.Normal.Z;
            var mat := info.Material;
            var surface := mat.GetColor(location);
            // Multiply the surface color by the ambient light factor.
            var color := 0.2F * surface;
            for light in Self.Lights do
            begin
                var factor: Single := light.Intensity(location);
                if factor > 0.0F then
                    color.Add(factor, mat.Radiance(
                        location, info.Normal, reflection, light, surface));
            end;
            Result += color.Multiply(weight, filter);
            bncs--;
            if bncs <= 0 then Break;
            weight *= mat.Reflection(var location, cosine, out lightFilter);
            if weight < MinWeight then Break;
            filter *= lightFilter;
            R.Origin := location;
            R.Direction := reflection;
        end;
        Result.Clip;
    end;

implementation for PerspectiveCamera is

    Transform: Matrix;
    Angle, Distance: Double;
    A, BRow, BCol: Double;

    constructor(
        Location, Target, Up: Vector;
        Angle: Double; Height, Width: Integer);
    begin
        Self.Location := Location;
        Self.Target := Target;
        var diff := Target - Location;
        if (Up ^ Diff).Length < Tolerance.Epsilon then
            if (Vector.ZRay ^ Diff).Length < Tolerance.Epsilon then
                Up := Vector.YRay
            else
                Up := Vector.ZRay;
        Self.Up := Up;
        Self.Angle := Angle;
        Self.Height := Height;
        Self.Width := Width;
        Self.PrimaryRay := new Ray;
        Self.Transform := new Matrix(Location, Target, Up);
        Self.Distance := (Location - Target).Length;
        var pixelSize := Distance * (Angle * Math.PI / 360).Tan;
        if Width >= Height then
            pixelSize /= Width
        else
            pixelSize /= Height;
        Self.A := 2 * pixelSize;
        Self.BRow := pixelSize * (1 - Height);
        Self.BCol := pixelSize * (1 - Width);
    end;

    constructor(Location, Target: Vector;
        Angle: Double; Height, Width: Integer)
        : Self(Location, Target, Vector.YRay, Angle, Height, Width);

    constructor(Location, Target: Vector; Height, Width: Integer)
        : Self(Location, Target, Vector.YRay, 60.0, Height, Width);

    method GetRay(Row, Col: Integer);
    begin
        PrimaryRay.Origin := Location;
        PrimaryRay.Direction := Transform.RotateNormalize(
            A * Col + BCol, A * Row + BRow, Distance);
    end;

implementation for PointLight is

    Location: Vector;
    TestRay: Ray;
    Root: IShape;

    constructor(Location: Vector; C: Pixel);
    begin
        Self.Location := Location;
        Self.Color := C;
        Self.TestRay := new Ray;
    end;

    method GetDirection(L: Vector): Vector =>
        (Location - L).Normalize;

    method Initialize(Scene: Scene);
    begin
        Self.Root := Scene.Root;
    end;

    method Intensity(L: Vector): Single;
    begin
        TestRay.Origin := Self.Location;
        TestRay.Direction := Self.Location - L;
        Result := if Root.ShadowTest(TestRay) then 0.0F else 1.0F;
    end;

implementation for Union is

    Shapes: Array[IShape];

    constructor(Shapes: Array[IShape]);
    begin
        Self.Shapes := Shapes;
    end;

    constructor(Shape1, Shape2: IShape);
    begin
        Self.Shapes := [Shape1, Shape2];
    end;

    constructor(Shape1, Shape2, Shape3: IShape);
    begin
        Self.Shapes := [Shape1, Shape2, Shape3];
    end;

    constructor(Shape1, Shape2, Shape3, Shape4: IShape);
    begin
        Self.Shapes := [Shape1, Shape2, Shape3, Shape4];
    end;

    method IShape.ShadowTest(R: Ray): Boolean;
    begin
        for shape in Shapes do
            if shape.ShadowTest(R) then
            begin
                Result := true;
                Exit;
            end;
        Result := false;
    end;

    method IShape.HitTest(
        R: Ray; MaxTime: Double; var Info: HitInfo): Boolean;
    begin
        for I: Integer := 0 to Shapes.Length - 1 do
            if Shapes[I].HitTest(R, MaxTime, var Info) then
            begin
                for J: Integer := I + 1 to Shapes.Length - 1 do
                    Shapes[J].HitTest(R, Info.Time, var Info);
                Result := true;
                Exit;
            end;
        Result := false;
    end;

implementation for Sphere is

    Center: Vector;
    Radius, R2, RInv: Double;
    Material: IMaterial;

    constructor(Center: Vector; Radius: Double; Material: IMaterial);
    begin
        Self.Center := Center;
        Self.Radius := Radius;
        Self.R2 := Radius.Sqr;
        Self.RInv := 1.0 / Radius;
        Self.Material := Material;
    end;

    method IShape.ShadowTest(R: Ray): Boolean;
    begin
        var delta := Center - R.Origin;
        var ray2 := 1.0 / R.Direction.Sqr;
        var beta := ray2 * R.Direction.Dot(delta);
        var discr := beta.Sqr - (delta.Sqr - R2) * ray2;
        if discr < 0.0 then
            Result := false
        else
        begin
            discr := discr.Sqrt;
            var t := beta - discr;
            if t < Tolerance.Epsilon then
            begin
                t := beta + discr;
                Result := t >= Tolerance.Epsilon and t <= 1.0;
            end
            else
                Result := t <= 1.0;
        end;
    end;

    method IShape.HitTest(
        R: Ray; MaxTime: Double; var Info: HitInfo): Boolean;
    var
        delta: Vector;
        beta, discr, t: Double;
    begin
        delta := Self.Center - R.Origin;
        beta := R.Direction.Dot(delta);
        discr := beta.Sqr - (delta.Sqr - R2);
        if discr < 0.0 then
            Result := false
        else
        begin
            discr := discr.Sqrt;
            t := beta - discr;
            if t < Tolerance.Epsilon then
            begin
                t := beta + discr;
                if t < Tolerance.Epsilon then
                begin
                    Result := false;
                    Exit;
                end;
            end;
            if t > MaxTime then
            begin
                Result := false;
                Exit;
            end;
            Info.Time := t;
            Info.HitPoint := R[t];
            Info.Normal := RInv * (Info.HitPoint - Center);
            Info.Material := Material;
            Result := true;
        end;
    end;

implementation for Cylinder is

    Bottom, Top, BNormal, TNormal: Vector;
    Radius, Height, R2, RInv: Double;
    Transform, Inverse: Matrix;
    Material: IMaterial;

    constructor(Bottom, Top: Vector; Radius: Double; Material: IMaterial);
    begin
        Self.Bottom := Bottom;
        Self.Top := Top;
        Self.Radius := Radius;
        Self.Material := Material;
        Self.R2 := Radius * Radius;
        Self.RInv := 1.0 / Radius;
        Self.Height := (Top - Bottom).Length;
        var v := (Top - Bottom).Normalize;
        var proj := (1.0 - v.Y.Sqr).Sqrt;
        if proj.Abs > Tolerance.Epsilon then
            Self.Transform := new Matrix(
                v.Z / proj, v.X, v.X * v.Y / proj,
                0.0, v.Y, -proj,
                -v.X / proj, v.Z, v.Y * v.Z / proj)
        else if v.Y > 0.0 then
            Self.Transform := new Matrix(+1, 0, 0, 0, +1, 0, 0, 0, +1)
        else
            Self.Transform := new Matrix(+1, 0, 0, 0, -1, 0, 0, 0, +1);
        Inverse := Transform;
        Self.BNormal := new Vector(-Inverse.A12, -Inverse.A22, -Inverse.A32);
        Self.TNormal := new Vector(Inverse.A12, Inverse.A22, Inverse.A32);
        Inverse.Transpose;
    end;

    method IShape.ShadowTest(R: Ray): Boolean;
    begin
        var Org := Inverse.Rotate(R.Origin - Bottom);
        var Dir := Inverse.Rotate(R.Direction);
        var tt0 := Tolerance.Epsilon;
        var tt1 := 1.0;
        var temp := 1.0 / Dir.Y;
        var t0 := -Org.Y * temp;
        var t1 := (Height - Org.Y) * temp;
        if temp >= 0.0 then
        begin
            if t0 > tt0 then tt0 := t0;
            if t1 < 1.0 then tt1 := t1;
        end
        else
        begin
            if t1 > tt0 then tt0 := t1;
            if t0 < 1.0 then tt1 := t0;
        end;
        if tt0 > tt1 then
            Exit;
        temp := 1.0 / (Dir.X.Sqr + Dir.Z.Sqr);
        var beta := -(Dir.X.Sqr + Dir.Z.Sqr) * temp;
        temp := beta.Sqr - (Org.X.Sqr + Org.Z.Sqr - R2) * temp;
        if temp < 0.0 then
            Exit;
        temp := temp.Sqrt;
        t0 := beta - temp;
        t1 := beta + temp;
        if t0 > tt0 then tt0 := t0;
        if t1 < tt1 then tt1 := t1;
        Result := tt0 <= tt1;
    end;

    method IShape.HitTest(
        R: Ray; MaxTime: Double; var Info: HitInfo): Boolean;
    var
        tt0, tt1: Double;
    begin
        var Org := Inverse.Rotate(R.Origin - Bottom);
        var Dir := Inverse.Rotate(R.Direction);
        // Compute time bounds for the Y axis.
        var temp := 1.0 / Dir.Y;
        if temp >= 0 then
        begin
           tt0 := -Org.Y * temp;
           tt1 := (Height - Org.Y) * temp;
        end
        else
        begin
            tt1 := -Org.Y * temp;
            tt0 := (Height - Org.Y) * temp;
        end;
        // Compute time bounds for the XZ plane.
        temp := 1.0 / (1.0 - Dir.Y.Sqr);
        var beta := -(Dir.X * Org.X + Dir.Z * Org.Z) * temp;
        temp := beta * beta - (Org.X.Sqr + Org.Z.Sqr - R2) * temp;
        if temp < 0 then
            Exit;
        temp := temp.Sqrt;
        var t0 := beta - temp;
        var t1 := beta + temp;
        if t0 > tt0 then tt0 := t0;
        if t1 < tt1 then tt1 := t1;
        if tt0 > tt1 then
            Exit;
        // Check the first intersection.
        if tt0 < Tolerance.Epsilon then
        begin
            tt0 := tt1;
            if tt0 < Tolerance.Epsilon then
                Exit;
        end;
        if tt0 > MaxTime then
            Exit;
        // Fill the hit information record.
        Info.Time := tt0;
        beta := Org.Y + tt0 * Dir.Y;
        if beta < Tolerance.Epsilon then
            Info.Normal := BNormal
        else if Height - beta < Tolerance.Epsilon then
            Info.Normal := TNormal
        else
            Info.Normal := Transform.Rotate(
                (Org.X + tt0 * Dir.X) * RInv,
                0.0,
                (Org.Z + tt0 * Dir.Z) * RInv);
        Info.HitPoint := R[tt0];
        Info.Material := material;
        Result := true;
    end;

implementation for Plastic is

    Color: Pixel;
    Reflection: Single;
    PhongAmount, PhongSize: Double;

    constructor(Color: Pixel; Reflection, PhongAmount, PhongSize: Double);
    begin
        Self.Color := Color;
        Self.Reflection := Math.Max(Math.Min(1.0F, Single(Reflection)), 0.0F);
        Self.PhongAmount := PhongAmount;
        Self.PhongSize := PhongSize + 1.0;
    end;

    constructor(Color: Pixel; Reflection: Double)
        : Self(Color, Reflection, 0.0, 1.0);

    constructor(Color: Pixel) : Self(Color, 0.0, 0.0, 1.0);

    method IMaterial.GetColor(Location: Vector): Pixel => Color;

    method IMaterial.Radiance(
        Location, Normal, Reflection: Vector;
        Light: ILight; Color: Pixel): Pixel;
    begin
        var lightRay := Light.GetDirection(Location);
        var cosine := Single(lightRay.Dot(Normal));
        if cosine > 0.0F then
            Result := Light.Color.Multiply(cosine, Color);
        if PhongAmount > 0.0 then
        begin
            var f := lightRay.Dot(Reflection);
            if f > 0.0 then
                Result.Add(
                    Single(PhongAmount * Math.Pow(f, PhongSize)),
                    Light.Color);
        end;
    end;

    method IMaterial.Reflection(
        var Location: Vector; Cosine: Double; out Filter: Pixel): Single;
    begin
        Filter := Pixel.White;
        Result := Reflection;
    end;

implementation for Metal is

    Color, LightFilter: Pixel;
    MinReflection, MaxReflection, Delta: Single;
    Diffuse, PhongAmount, PhongSize: Double;

    constructor(Color: Pixel;
        MinReflection, MaxReflection, Diffuse: Double;
        PhongAmount, PhongSize: Double);
    begin
        Self.Color := Color;
        Self.MinReflection :=
            Math.Max(Math.Min(1.0F, Single(MinReflection)), 0.0F);
        Self.MaxReflection :=
            Math.Max(Math.Min(1.0F, Single(MaxReflection)), 0.0F);
        Self.Diffuse := Diffuse;
        Self.PhongAmount := PhongAmount;
        Self.PhongSize := PhongSize + 1.0;
        Self.Delta := Self.MaxReflection - Self.MinReflection;
        var maxClr := Math.Max(Color.R, Math.Max(Color.G, Color.B));
        if maxClr < 1.0 / 255.0 then
            Self.LightFilter := Pixel.White
        else
            Self.LightFilter := new Pixel(
                1.0F - Color.R / maxClr,
                1.0F - Color.G / maxClr,
                1.0F - Color.B / maxClr);
    end;

    constructor(Color: Pixel; MinReflection, MaxReflection, Diffuse: Double)
        : Self(Color, MinReflection, MaxReflection, Diffuse, 0.0, 1.0);

    constructor(Color: Pixel; MinReflection, MaxReflection: Double)
        : Self(Color, MinReflection, MaxReflection, 1.0, 0.0, 1.0);

    constructor(Color: Pixel; MinReflection: Double)
        : Self(Color, MinReflection, 1.0, 1.0, 0.0, 1.0);

    method IMaterial.GetColor(Location: Vector): Pixel => Color;

    method IMaterial.Radiance(
        Location, Normal, Reflection: Vector;
        Light: ILight; Color: Pixel): Pixel;
    begin
        var lightRay := Light.GetDirection(Location);
        var cosine := Single(Diffuse * lightRay.Dot(Normal));
        if cosine > 0.0F then
            Result := cosine * Color;
        if PhongAmount > 0.0 then
        begin
            var f := lightRay.Dot(Reflection);
            if f > 0.0 then
                Result += Single(PhongAmount * Math.Pow(f, PhongSize)) *
                    Color.Interpolate(cosine * cosine.Sqr.Sqr, Light.Color);
        end;
    end;

    method IMaterial.Reflection(
        var Location: Vector; Cosine: Double; out Filter: Pixel): Single;
    begin
        if Cosine < 0.0 then
        begin
            var c5 := 1.0F + Single(Cosine);
            c5 *= c5.Sqr.Sqr;
            // Interpolate between lightFilter and Pixel.White.
            Filter := new Pixel(
                1.0F - c5 * LightFilter.R,
                1.0F - c5 * LightFilter.G,
                1.0F - c5 * LightFilter.B);
            // Schlick's aproximation for Fresnel coefficient.
            Result := MinReflection + c5 * Delta;
        end
        else
        begin
            Filter := new Pixel(
                1.0F - LightFilter.R, 1.0F - LightFilter.G, 1.0F - LightFilter.B);
            Result := MaxReflection;
        end;
    end;

implementation

    method CreateMaterial1(AColor: System.Drawing.Color): IMaterial =>
        new Metal(new Pixel(AColor), 0.2, 1.0, 1.0, 0.4, 10.0);

    method CreateMaterial2(AColor: System.Drawing.Color): IMaterial =>
        new Plastic(new Pixel(AColor), 0.6);

    method CreateScene1: Scene =>
        new Scene(
            new BasicSampler(10, 0.001),
            new PerspectiveCamera(
                new Vector(0, 0, -10), Vector.Null, 200, 320),
            new Union(
                [
                new Sphere(new Vector( 0.0,0,0), 1.0,
                    CreateMaterial1(System.Drawing.Color.Orchid)),
                new Sphere(new Vector(-1.7,0,0), 0.7,
                    CreateMaterial1(System.Drawing.Color.RoyalBlue)),
                new Sphere(new Vector(+1.7,0,0), 0.7,
                    CreateMaterial1(System.Drawing.Color.RoyalBlue)),
                new Sphere(new Vector(-2.9,0,0), 0.5,
                    CreateMaterial1(System.Drawing.Color.ForestGreen)),
                new Sphere(new Vector(+2.9,0,0), 0.5,
                    CreateMaterial1(System.Drawing.Color.ForestGreen))
                ]),
            [new PointLight(new Vector(10, 10, -10), Pixel.White)]);

    method CreateScene2: Scene =>
        new Scene(
            new BasicSampler(12, 0.001),
            new PerspectiveCamera(
                new Vector(0, 0, -5), Vector.Null, 300, 300),
            new Union(
                new Sphere(new Vector(+1.0,-1.0,0.0), 1.41,
                    CreateMaterial2(System.Drawing.Color.Blue)),
                new Sphere(new Vector(-1.0,+1.0,0.0), 1.41,
                    CreateMaterial2(System.Drawing.Color.Green))),
            [new PointLight(new Vector(-2, 3, -10), new Pixel(0.6))]);

    method Main;
    begin
        try
            var time := Environment.TickCount;
            var map := CreateScene1.Render;
            time := Environment.TickCount - time;
            Console.WriteLine(String.Format(
                'Time ellapsed: {0} milliseconds.', time));
            if map <> nil then
                map.ToBitmap.Save('c:\xsbench.bmp');
        except on e: Exception do
            Console.WriteLine(e.Message);
            Console.WriteLine(e.StackTrace);
        end;
    end;

end.