Bounding boxes (#5)

RayTracing/BoundingBox.fs

@@ -0,0 +1,99 @@
+namespace RayTracing
+
+[<Struct ; NoComparison ; NoEquality>]
+type BoundingBox =
+    {
+        Min : Point
+        Max : Point
+    }
+
+[<RequireQualifiedAccess>]
+module BoundingBox =
+
+    let volume (box : BoundingBox) =
+        (Point.coordinate 0 box.Max - Point.coordinate 0 box.Min) *
+        (Point.coordinate 1 box.Max - Point.coordinate 1 box.Min) *
+        (Point.coordinate 2 box.Max - Point.coordinate 2 box.Min)
+
+    let make (min : Point) (max : Point) =
+        {
+            Min = min
+            Max = max
+        }
+
+
+    let inverseDirections (ray : Ray) =
+        struct(1.0 / (Ray.vector ray |> UnitVector.coordinate 0), 1.0 / (Ray.vector ray |> UnitVector.coordinate 1), 1.0 / (Ray.vector ray |> UnitVector.coordinate 2))
+
+    let hits (struct(invX, invY, invZ)) { Ray.Origin = Point (x, y, z) ; Vector = UnitVector (Vector (dx, dy, dz))} (box : BoundingBox) : bool =
+        // The line is (x, y, z) + t (dx, dy, dz)
+        // The line goes through the cuboid iff it passes through the interval in each component:
+        //   there is t such that boxMin.X <= x + t dx <= boxMax.X,
+        // and moreover the acceptable t brackets all overlap.
+        // That is,
+        //  boxMin.X - x <= t dx <= boxMax.X - x
+        let mutable tMin = -infinity
+        let mutable tMax = infinity
+
+        let bailOut =
+            let mutable t0 = (Point.coordinate 0 box.Min - x) * invX
+            let mutable t1 = (Point.coordinate 0 box.Max - x) * invX
+            if invX < 0.0 then
+                let tmp = t1
+                t1 <- t0
+                t0 <- tmp
+
+            tMin <- if t0 > tMin then t0 else tMin
+            tMax <- if t1 < tMax then t1 else tMax
+
+            tMax < tMin || 0.0 >= tMax
+
+        if bailOut then false else
+
+        let bailOut =
+            let mutable t0 = (Point.coordinate 1 box.Min - y) * invY
+            let mutable t1 = (Point.coordinate 1 box.Max - y) * invY
+
+            if invY < 0.0 then
+                let tmp = t1
+                t1 <- t0
+                t0 <- tmp
+
+            tMin <- if t0 > tMin then t0 else tMin
+            tMax <- if t1 < tMax then t1 else tMax
+
+            tMax < tMin || 0.0 >= tMax
+
+        if bailOut then false else
+
+        let mutable t0 = (Point.coordinate 2 box.Min - z) * invZ
+        let mutable t1 = (Point.coordinate 2 box.Max - z) * invZ
+
+        if invZ < 0.0 then
+            let tmp = t1
+            t1 <- t0
+            t0 <- tmp
+
+        tMin <- if t0 > tMin then t0 else tMin
+        tMax <- if t1 < tMax then t1 else tMax
+        tMax >= tMin && tMax >= 0.0
+
+    let mergeTwo (i : BoundingBox) (j : BoundingBox) : BoundingBox =
+        {
+            Min =
+                Point.make
+                    (min (Point.coordinate 0 i.Min) (Point.coordinate 0 j.Min))
+                    (min (Point.coordinate 1 i.Min) (Point.coordinate 1 j.Min))
+                    (min (Point.coordinate 2 i.Min) (Point.coordinate 2 j.Min))
+            Max =
+                Point.make
+                    (max (Point.coordinate 0 i.Max) (Point.coordinate 0 j.Max))
+                    (max (Point.coordinate 1 i.Max) (Point.coordinate 1 j.Max))
+                    (max (Point.coordinate 2 i.Max) (Point.coordinate 2 j.Max))
+        }
+
+    let merge (boxes : BoundingBox []) : BoundingBox option =
+        if boxes.Length = 0 then None else
+        boxes
+        |> Array.reduce mergeTwo
+        |> Some

RayTracing/BoundingBoxTree.fs

@@ -0,0 +1,39 @@
+namespace RayTracing
+
+type BoundingBoxTree =
+    | Leaf of hittable : Hittable * BoundingBox
+    | Branch of left : BoundingBoxTree * right : BoundingBoxTree * all : BoundingBox
+
+[<RequireQualifiedAccess>]
+module BoundingBoxTree =
+    let make (boxes : (Hittable * BoundingBox) array) : BoundingBoxTree option =
+        if boxes.Length = 0 then None else
+
+        let rec go (boxes : (Hittable * BoundingBox) array) =
+            let boundAll =
+                BoundingBox.merge (boxes |> Array.map snd) |> Option.get
+
+            if boxes.Length = 1 then Leaf boxes.[0] else
+            if boxes.Length = 2 then Branch (Leaf boxes.[0], Leaf boxes.[1], boundAll) else
+
+            let choices =
+                Array.init 3 (fun axis ->
+                    let boxes =
+                        boxes
+                        |> Array.sortBy (fun (_, b) -> Point.coordinate axis b.Min)
+                    let leftHalf = boxes.[0..boxes.Length / 2]
+                    let rightHalf = boxes.[(boxes.Length / 2) + 1..]
+                    let leftBound = leftHalf |> Array.map snd |> BoundingBox.merge |> Option.get
+                    let rightBound = rightHalf |> Array.map snd |> BoundingBox.merge |> Option.get
+                    (leftHalf, leftBound), (rightHalf, rightBound)
+                )
+            let (leftHalf, _), (rightHalf, _) =
+                choices
+                |> Array.minBy (fun ((_, leftBound), (_, rightBound)) ->
+                    (BoundingBox.volume leftBound) + (BoundingBox.volume rightBound)
+                )
+
+            Branch (go leftHalf, go rightHalf, boundAll)
+
+        go boxes
+        |> Some

RayTracing/Camera.fs

@@ -23,6 +23,8 @@ type Camera =
         FocalLength : float
         /// How many samples will we take per pixel, for anti-aliasing?
         SamplesPerPixel : int
+        /// How many bounces before we consider ourselves to have lost track of a light ray?
+        BounceDepth : int
     }
 
 [<RequireQualifiedAccess>]
@@ -53,4 +55,5 @@ module Camera =
             ViewportXAxis = xAxis
             ViewportYAxis = yAxis
             SamplesPerPixel = samplesPerPixel
+            BounceDepth = 150
         }

RayTracing/Hittable.fs

@@ -0,0 +1,38 @@
+namespace RayTracing
+
+type Hittable =
+    | Sphere of Sphere
+    | UnboundedSphere of Sphere
+    | InfinitePlane of InfinitePlane
+
+    member this.Reflection (incoming : LightRay) (strikePoint : Point) =
+        match this with
+        | Sphere s
+        | UnboundedSphere s -> s.Reflection incoming strikePoint
+        | InfinitePlane p -> p.Reflection incoming strikePoint
+
+    member this.BoundingBox : BoundingBox option =
+        match this with
+        | Sphere s -> Sphere.boundingBox s |> Some
+        | UnboundedSphere _
+        | InfinitePlane _ -> None
+
+[<RequireQualifiedAccess>]
+module Hittable =
+
+    let inline boundingBox (h : Hittable) = h.BoundingBox
+
+    /// Returns the distance we must walk along this ray before we first hit an object, the
+    /// colour the resulting light ray is after the interaction, and the new ray.
+    let hits
+        (ray : Ray)
+        (h : Hittable)
+        : float voption
+        =
+        match h with
+        | UnboundedSphere s
+        | Sphere s ->
+            Sphere.firstIntersection s ray
+        | InfinitePlane plane ->
+            InfinitePlane.intersection plane ray
+

RayTracing/Point.fs

@@ -6,12 +6,10 @@ open System.Runtime.CompilerServices
 /// We don't let you compare these for equality, because floats are hard.
 [<NoEquality ; NoComparison ; Struct ; IsReadOnly>]
 type Point =
-    private
     | Point of struct(float * float * float)
 
 [<NoEquality ; NoComparison ; Struct ; IsReadOnly>]
 type Vector =
-    private
     | Vector of struct(float * float * float)
 
 [<Struct ; IsReadOnly ; NoEquality ; NoComparison>]
@@ -80,10 +78,22 @@ module UnitVector =
             Vector (0.0, 0.0, 1.0) |> UnitVector
         |]
 
+    let inline coordinate (i : int) (UnitVector (Vector (a, b, c))) : float =
+        match i with
+        | 0 -> a
+        | 1 -> b
+        | 2 -> c
+        | _ -> failwithf "Bad coordinate: %i" i
+
 [<RequireQualifiedAccess>]
 module Point =
 
-    let xCoordinate (Point (x, _, _)) = x
+    let inline coordinate (i : int) (Point (x, y, z)) =
+        match i with
+        | 0 -> x
+        | 1 -> y
+        | 2 -> z
+        | _ -> failwithf "Bad coordinate: %i" i
 
     let sum (Point (a, b, c)) (Point (x, y, z)) : Point =
         Point (a + x, b + y, c + z)

RayTracing/Ray.fs

@@ -1,11 +1,10 @@
 namespace RayTracing
 
 type Ray =
-    private
-        {
-            Origin : Point
-            Vector : UnitVector
-        }
+    {
+        Origin : Point
+        Vector : UnitVector
+    }
 
 [<RequireQualifiedAccess>]
 module Ray =
@@ -47,8 +46,8 @@ module Ray =
                 Float.equal t t3
             else false
 
-    let vector r = r.Vector
-    let origin r = r.Origin
+    let inline vector r = r.Vector
+    let inline origin r = r.Origin
 
     let flip (r : Ray) =
         {

RayTracing/Ray.fsi

@@ -1,6 +1,12 @@
 namespace RayTracing
 
-type Ray
+type Ray =
+    {
+        /// For performance reasons, this is public, but please don't use it
+        Origin : Point
+        /// For performance reasons, this is public, but please don't use it
+        Vector : UnitVector
+    }
 
 [<RequireQualifiedAccess>]
 module Ray =
@@ -13,7 +19,7 @@ module Ray =
 
     val liesOn : Point -> Ray -> bool
 
-    val vector : Ray -> UnitVector
-    val origin : Ray -> Point
+    val inline vector : Ray -> UnitVector
+    val inline origin : Ray -> Point
 
     val flip : Ray -> Ray

RayTracing/RayTracing.fsproj

@@ -13,10 +13,13 @@
         <Compile Include="Ray.fs" />
         <Compile Include="Plane.fs" />
         <Compile Include="LightRay.fs" />
+        <Compile Include="BoundingBox.fs" />
         <Compile Include="Sphere.fs" />
         <Compile Include="InfinitePlane.fs" />
         <Compile Include="ImageOutput.fs" />
         <Compile Include="Camera.fs" />
+        <Compile Include="Hittable.fs" />
+        <Compile Include="BoundingBoxTree.fs" />
         <Compile Include="Scene.fs" />
     </ItemGroup>
 

RayTracing/Scene.fs

@@ -2,93 +2,79 @@ namespace RayTracing
 
 open System
 
-type Hittable =
-    | Sphere of Sphere
-    | InfinitePlane of InfinitePlane
-
-    member this.Reflection (incoming : LightRay) (strikePoint : Point) =
-        match this with
-        | Sphere s -> s.Reflection incoming strikePoint
-        | InfinitePlane p -> p.Reflection incoming strikePoint
-
-[<RequireQualifiedAccess>]
-module Hittable =
-    /// Returns the distance we must walk along this ray before we first hit an object, the
-    /// colour the resulting light ray is after the interaction, and the new ray.
-    let hits
-        (ray : Ray)
-        (h : Hittable)
-        : float voption
-        =
-        match h with
-        | Sphere s ->
-            Sphere.firstIntersection s ray
-        | InfinitePlane plane ->
-            InfinitePlane.intersection plane ray
-
 type Scene =
-    {
-        Objects : Hittable array
-    }
+    private
+        {
+            UnboundedObjects : Hittable array
+            BoundingBoxes : BoundingBoxTree option
+        }
 
 [<RequireQualifiedAccess>]
 module Scene =
 
+    let make (objects : Hittable array) =
+        let bounded, unbounded =
+            objects
+            |> Array.map (fun h -> h, Hittable.boundingBox h)
+            |> Array.partition (snd >> Option.isSome)
+        let bounded = bounded |> Array.map (fun (h, box) -> h, Option.get box)
+        let unbounded = unbounded |> Array.map fst
+        let tree =
+            bounded
+            |> BoundingBoxTree.make
+        {
+            UnboundedObjects = unbounded
+            BoundingBoxes = tree
+        }
+
+    let rec bestCandidate (inverseDirections : struct(float * float * float)) (ray : Ray) (bestFloat : float) (bestObject : Hittable) (bestLength : float) (box : BoundingBoxTree) : struct(float * Hittable * float) =
+        match box with
+        | BoundingBoxTree.Leaf (object, box) ->
+            if BoundingBox.hits inverseDirections ray box then
+                match Hittable.hits ray object with
+                | ValueNone -> struct (bestFloat, bestObject, bestLength)
+                | ValueSome point ->
+                    let a = point * point
+                    if a < bestFloat then
+                        struct (a, object, point)
+                    else
+                        struct (bestFloat, bestObject, bestLength)
+            else struct (bestFloat, bestObject, bestLength)
+        | BoundingBoxTree.Branch (left, right, all) ->
+            if BoundingBox.hits inverseDirections ray all then
+                let struct (bestFloat, bestObject, bestLength) = bestCandidate inverseDirections ray bestFloat bestObject bestLength left
+                bestCandidate inverseDirections ray bestFloat bestObject bestLength right
+            else struct (bestFloat, bestObject, bestLength)
+
     let hitObject
         (s : Scene)
         (ray : Ray)
-        : (int * Point) option
+        : (Hittable * Point) option
         =
-        let mutable bestIndex = -1
+        let mutable best = Unchecked.defaultof<_>
         let mutable bestLength = nan
         let mutable bestFloat = infinity
-        for i in 0..s.Objects.Length - 1 do
-            match Hittable.hits ray s.Objects.[i] with
+
+        match s.BoundingBoxes with
+        | None -> ()
+        | Some boundingBoxes ->
+            let struct(f, o, l) = bestCandidate (BoundingBox.inverseDirections ray) ray bestFloat best bestLength boundingBoxes
+            bestFloat <- f
+            best <- o
+            bestLength <- l
+
+        for i in s.UnboundedObjects do
+            match Hittable.hits ray i with
             | ValueNone -> ()
             | ValueSome point ->
                 let a = point * point
-                match Float.compare a bestFloat with
-                | Less ->
+                if Float.compare a bestFloat = Less then
                     bestFloat <- a
-                    bestIndex <- i
+                    best <- i
                     bestLength <- point
-                | _ -> ()
 
         if Double.IsNaN bestLength then None else
-        Some (bestIndex, Ray.walkAlong ray bestLength)
-
-    let internal traceRayPrinting
-        (print : string -> unit)
-        (maxCount : int)
-        (scene : Scene)
-        (ray : LightRay)
-        : Pixel
-        =
-        let rec go (bounces : int) (ray : LightRay) : Pixel =
-            let (Point(x, y, z)) = Ray.origin ray.Ray
-            let (UnitVector (Vector(a, b, c))) = Ray.vector ray.Ray
-            print (sprintf "Ray, colour %i,%i,%i\n  origin (%f, %f, %f)\n  vector (%f, %f, %f)" ray.Colour.Red ray.Colour.Green ray.Colour.Blue x y z a b c)
-            if bounces > maxCount then Colour.HotPink else
-
-            let thingsWeHit = hitObject scene ray.Ray
-            match thingsWeHit with
-            | None ->
-                print ">>> No object collision; black."
-                // Ray goes off into the distance and is never heard from again
-                Colour.Black
-            | Some (objectNumber, strikePoint) ->
-                let (Point(x, y, z)) = strikePoint
-                print (sprintf ">>> collided with object %i at (%f, %f, %f)" objectNumber x y z)
-                let outgoingRay = scene.Objects.[objectNumber].Reflection ray strikePoint
-                match outgoingRay with
-                | Absorbs colour ->
-                    print (sprintf ">>>   surface absorbs, yielding colour %i,%i,%i" colour.Red colour.Green colour.Blue)
-                    colour
-                | Continues outgoingRay ->
-                    print ">>>   continuing tracing."
-                    go (bounces + 1) outgoingRay
-
-        go 0 ray
+        Some (best, Ray.walkAlong ray bestLength)
 
     let internal traceRay
         (maxCount : int)
@@ -106,8 +92,8 @@ module Scene =
             | None ->
                 // Ray goes off into the distance and is never heard from again
                 Colour.Black
-            | Some (objectNumber, strikePoint) ->
-                let outgoingRay = scene.Objects.[objectNumber].Reflection ray strikePoint
+            | Some (object, strikePoint) ->
+                let outgoingRay = object.Reflection ray strikePoint
                 match outgoingRay with
                 | Absorbs colour ->
                     colour
@@ -118,7 +104,7 @@ module Scene =
 
     /// Trace a ray to this one pixel, updating the PixelStats with the result.
     /// n.b. not thread safe
-    let private traceOnce (print : string -> unit) (scene : Scene) (rand : FloatProducer) (camera : Camera) (maxWidthCoord : int) (maxHeightCoord : int) row col stats =
+    let private traceOnce (scene : Scene) (rand : FloatProducer) (camera : Camera) (maxWidthCoord : int) (maxHeightCoord : int) row col stats =
         let struct(rand1, rand2) = rand.GetTwo ()
         let landingPoint =
             ((float col + rand1) * camera.ViewportWidth) / float maxWidthCoord
@@ -134,15 +120,10 @@ module Scene =
             |> Option.get
 
         // Here we've hardcoded that the eye is emitting white light through a medium with refractance 1.
-        let result = traceRay 150 scene { Ray = ray ; Colour = Colour.White }
-        //if result = Colour.HotPink then
-        //    print "hi"
-        //    traceRayPrinting print 150 scene { Ray = ray ; Colour = Colour.White ; Refractance = 1.0<ior> }
-        //    |> ignore
-        //    failwith "Stopping."
+        let result = traceRay camera.BounceDepth scene { Ray = ray ; Colour = Colour.White }
         PixelStats.add result stats
 
-    let renderPixel (print : string -> unit) (scene : Scene) (rand : FloatProducer) (camera : Camera) maxWidthCoord maxHeightCoord row col =
+    let renderPixel (_ : string -> unit) (scene : Scene) (rand : FloatProducer) (camera : Camera) maxWidthCoord maxHeightCoord row col =
         // Where does this pixel correspond to, on the imaginary canvas?
         // For the early prototype, we'll just take the upper right quadrant
         // from the camera.
@@ -151,12 +132,12 @@ module Scene =
         let firstTrial = min 5 (camera.SamplesPerPixel / 2)
 
         for _ in 0..firstTrial do
-            traceOnce print scene rand camera maxWidthCoord maxHeightCoord row col stats
+            traceOnce scene rand camera maxWidthCoord maxHeightCoord row col stats
 
         let oldMean = PixelStats.mean stats
 
         for _ in 1..firstTrial do
-            traceOnce print scene rand camera maxWidthCoord maxHeightCoord row col stats
+            traceOnce scene rand camera maxWidthCoord maxHeightCoord row col stats
 
         let newMean = PixelStats.mean stats
         let difference = Pixel.difference newMean oldMean
@@ -168,7 +149,7 @@ module Scene =
         else
 
             for _ in 1..(camera.SamplesPerPixel - 2 * firstTrial - 1) do
-                traceOnce print scene rand camera maxWidthCoord maxHeightCoord row col stats
+                traceOnce scene rand camera maxWidthCoord maxHeightCoord row col stats
 
             PixelStats.mean stats
 

RayTracing/Sphere.fs

@@ -17,6 +17,7 @@ type Sphere =
             /// what colour ray does it output and in what direction?
             Reflection : LightRay -> Point -> LightDestination
             RadiusSquared : float
+            BoundingBox : BoundingBox
         }
 
 type SphereStyle =
@@ -152,9 +153,9 @@ module Sphere =
         | SphereStyle.LightSource colour ->
             Absorbs (Pixel.combine incomingLight.Colour colour)
         | SphereStyle.LightSourceCap colour ->
-            let circleCentreZCoord = Point.xCoordinate centre
+            let circleCentreZCoord = Point.coordinate 0 centre
             let zCoordLowerBound = circleCentreZCoord + (radius - (radius / 4.0))
-            let strikeZCoord = Point.xCoordinate strikePoint
+            let strikeZCoord = Point.coordinate 0 strikePoint
             let colour =
                 match Float.compare strikeZCoord zCoordLowerBound with
                 | Greater ->
@@ -238,8 +239,11 @@ module Sphere =
             Radius = radius
             Reflection = reflection style centre radius radiusSquared (Float.compare radius 0.0 = Less)
             RadiusSquared = radiusSquared
+            BoundingBox = BoundingBox.make (Point.sum centre (Point.make -radius -radius -radius)) (Point.sum centre (Point.make radius radius radius))
         }
 
+    let boundingBox (s : Sphere) = s.BoundingBox
+
     let liesOn (point : Point) (sphere : Sphere) : bool =
         liesOn' sphere.Centre sphere.Radius point