packages feed

raytrace 0.1.0.0 → 0.2.0.0

raw patch · 8 files changed

+572/−187 lines, 8 files

Files

CHANGELOG.md view
@@ -3,3 +3,16 @@ ## 0.1.0.0 — 2025-07-18  * First version. Released on an unsuspecting world.++## 0.2.0.0 — 2025-07-31++* New material model allowing ray redirection+* Added camera setting `cs_redirectTargets`+* Added triangles and meshes+* Added `planeShape` function+* Added time parameter to `Geometry`, allowing for motion blur+* Different type signatures for `boxJoin`, `boxHull`, and `bvhTree`+* Removed `autoTree` function+* New materials: `anisotropic` and `lommelSeeliger`+* Volumes can be non-convex+* All rays are normalized
README.md view
@@ -1,39 +1,41 @@ # raytrace -A Haskell ray tracing library. Largely based on the books [Ray Tracing in One Weekend](https://raytracing.github.io/books/RayTracingInOneWeekend.html) and [Ray Tracing: The Next Week](https://raytracing.github.io/books/RayTracingTheNextWeek.html) by Peter Shirley.+A Haskell library for rendering images using path tracing. Inspired by the book series [Ray Tracing in One Weekend](https://raytracing.github.io/) by Peter Shirley and Steve Hollasch. +[View on Hackage](https://hackage.haskell.org/package/raytrace)+ Features:-* Spheres, parallelograms, and boxes (you can define your own shapes in addition to these)-* Volumes (fog and the like)-* A variety of materials, with behaviors including light emission and refraction+* Spheres, parallelograms, and boxes+* Triangle meshes (with a parser for .obj files)+* Volumes (fog, subsurface scattering)+* A variety of materials, with behaviors including light emission and refraction (you can define additional materials in Haskell) * Texture mapping * Perlin noise for procedurally generated textures-* Parallel computation of pixels+* Parallel computation of pixels (if you compile your program with the flags `-threaded -rtsopts -with-rtsopts=-N`) * Bounding volume hierarchies-* Affine transformations (rotation, translation, and so on)-* Optional defocusing to simulate a real camera lens--Possible future additions:+* The ability to prioritize light sources when choosing a scattered ray direction, leading to faster convergence without sacrificing accuracy (see the section "Ray Redirection" below)+* Transformed instances * Motion blur-* Triangular meshes-* Less noise in scenes with small light sources+* Optional defocusing to simulate a real camera lens  ![Example](https://raw.githubusercontent.com/UnaryPlus/raytrace/refs/heads/main/demo1.png) -The image above, with 405 million top-level rays, was generated on my laptop in about 8 minutes. The blurriness in the foreground and background is due to defocusing, wherein only a single plane is in focus. The following image demonstrates texture mapping, light sources, and fog:+The image above, with 405 million top-level rays, was generated on my laptop in about 8 minutes. The blurriness in the foreground and background is due to defocusing, wherein only a single plane is in focus. The following image demonstrates texture mapping, light sources, motion blur, and fog. (The code for both of these images was based on code in the aforementioned books.)  ![Example](https://raw.githubusercontent.com/UnaryPlus/raytrace/refs/heads/main/demo2.png) -(The code for both of these images was based on code in the aforementioned books.)+Here is an image demonstrating polygonal meshes and subsurface scattering: -## Example Usage+![Example](https://raw.githubusercontent.com/UnaryPlus/raytrace/refs/heads/main/pawn_demo.png) +## Example Code+ ```haskell module Main where  import Graphics.Ray import Linear (V3(V3))-import System.Random (newStdGen)+import System.Random (mkStdGen) import Data.Functor.Identity (Identity)  world :: Geometry Identity Material@@ -53,11 +55,23 @@   }  main :: IO ()-main = do-  seed <- newStdGen+main = +  let seed = mkStdGen 100 in   writeImage "example_image.png" (raytrace settings world seed) ```  This produces the following image:  ![Example](https://raw.githubusercontent.com/UnaryPlus/raytrace/refs/heads/main/example_image.png)++## Ray Redirection++For scenes with small light sources, the basic path tracing algorithm can result in very noisy images, as most rays never reach a light source. To reduce the noise, we can alter the distribution of scattered ray directions, increasing the probability of sending a ray toward a light source and decreasing the probability of other directions. This sounds like it would introduce bias and make surfaces appear brighter than they should, but we can correct for that by weighting the resulting colors, downweighting directions that we have made more common and upweighting directions that we have made less common. The mathematical basis for this is known as [Monte Carlo integration](https://en.wikipedia.org/wiki/Monte_Carlo_integration), specifically the "importance sampling" method.++Here is an example using the [Cornell Box](https://en.wikipedia.org/wiki/Cornell_box). The first image was generated with `cs_redirectTargets = []` (the default), and the second with `cs_redirectTargets = [ (0.25, V3 343 554 332, V3 (-130) 0 0, V3 0 0 (-105)) ]`. (This means that 25% of scattered rays are sent toward a random point on the parallelogram defined by the three vectors, which happens to be the light source.) Both images were generated with the same number of samples per pixel (200), but the second image is much less noisy. It also took half as long to render due to rays terminating at the light sooner.++![Cornell box without redirection](https://raw.githubusercontent.com/UnaryPlus/raytrace/refs/heads/main/cornell_box_noisy.png)++![Cornell box with redirection](https://raw.githubusercontent.com/UnaryPlus/raytrace/refs/heads/main/cornell_box_redirect.png)++Of course, the second image still has some noise, which could be reduced by increasing the number of samples per pixel.
raytrace.cabal view
@@ -1,13 +1,12 @@ cabal-version:      3.0 name:               raytrace-version:            0.1.0.0+version:            0.2.0.0  category:           Graphics synopsis:           Ray tracing library description:        -    A Haskell ray tracing library. -    Largely based on the books [Ray Tracing in One Weekend](https://raytracing.github.io/books/RayTracingInOneWeekend.html) -    and [Ray Tracing: The Next Week](https://raytracing.github.io/books/RayTracingTheNextWeek.html) by Peter Shirley.+    A Haskell library for rendering images using path tracing. +    Inspired by the book series [Ray Tracing in One Weekend](https://raytracing.github.io/) by Peter Shirley and Steve Hollasch.          See README.md for more information. 
src/Graphics/Ray.hs view
@@ -25,7 +25,7 @@ import Graphics.Ray.Texture import Graphics.Ray.Noise -import Linear (V2(V2), V3(V3), (*^), (^*), normalize, cross, (^/), zero)+import Linear (V2(V2), V3(V3), (*^), (^*), normalize, cross, (^/), zero, dot) import System.Random (StdGen, random, splitGen) import Data.Massiv.Array (B, D, S, U, Ix2((:.)), (!)) import qualified Data.Massiv.Array as A@@ -35,20 +35,36 @@ import Control.Monad.State (State, state, evalState) import Control.Monad (replicateM) import Data.Functor.Identity (Identity, runIdentity)+import Data.Maybe (listToMaybe)  data CameraSettings = CameraSettings-  { cs_center :: Point3 -- ^ Camera position -  , cs_lookAt :: Point3 -- ^ Point for the camera to look at -  , cs_up :: Vec3 -- ^ Camera \"up\" vector -  , cs_vfov :: Double -- ^ Vertical field of view (in radians) -  , cs_aspectRatio :: Double -- ^ Width-to-height ratio of image -  , cs_imageWidth :: Int -- ^ Image width in pixels -  , cs_samplesPerPixel :: Int -- ^ Number of top-level rays created per pixel -  , cs_maxRecursionDepth :: Int -- ^ Number of times a ray can reflect before recursion stops -  , cs_background :: Ray -> Color -- ^ Background color (which can depend on direction) -  , cs_defocusAngle :: Double -- ^ If this is positive, the image will be somewhat blurry in the foreground and background, -                              -- with only a single plane in focus (like an image produced by a real camera)-  , cs_focusDist :: Double -- ^ Distance from the camera to the plane of focus (only matters if defocus angle is nonzero) +  { cs_center :: Point3 +    -- ^ Camera position +  , cs_lookAt :: Point3 +    -- ^ Point for the camera to look at +  , cs_up :: Vec3 +    -- ^ Camera \"up\" vector +  , cs_vfov :: Double +    -- ^ Vertical field of view (in radians) +  , cs_aspectRatio :: Double +    -- ^ Width-to-height ratio of image +  , cs_imageWidth :: Int +    -- ^ Image width in pixels +  , cs_samplesPerPixel :: Int +    -- ^ Number of top-level rays created per pixel +  , cs_maxRecursionDepth :: Int +    -- ^ Number of times a ray can reflect before recursion stops +  , cs_background :: Ray -> Color +    -- ^ Background color (which can depend on direction) +  , cs_defocusAngle :: Double +    -- ^ If this is positive, the image will be somewhat blurry in the foreground and background, +    -- with only a single plane in focus (like an image produced by a real camera)+  , cs_focusDist :: Double +    -- ^ Distance from the camera to the plane of focus (only matters if defocus angle is nonzero) +  , cs_redirectTargets :: [(Double, Point3, Vec3, Vec3)] +    -- ^ List of parallelograms to send scattered rays toward (if the material allows it) along with +    -- probabilities, which should add up to less than 1. Sending more rays toward lights in +    -- scenes with small light sources can lead to faster convergence (i.e. less noise).   }  -- | By default, the camera is positioned at the origin looking in the negative z direction, with the positive y direction being upward@@ -63,6 +79,7 @@ -- cs_background = const (V3 1 1 1) -- cs_defocusAngle = 0.0 -- cs_focusDist = 10.0+-- cs_redirectTargets = [] -- @ defaultCameraSettings :: CameraSettings defaultCameraSettings = CameraSettings@@ -77,6 +94,7 @@   , cs_background = const (V3 1 1 1)   , cs_defocusAngle = 0.0   , cs_focusDist = 10.0+  , cs_redirectTargets = []   }  class ToRandom m where@@ -90,6 +108,15 @@   toRandom :: State StdGen a -> State StdGen a   toRandom = id +-- [private]+data RedirectTarget = RedirectTarget+  { rt_origin :: Point3+  , rt_U :: Vec3+  , rt_V :: Vec3+  , rt_cross :: Vec3+  , rt_hit :: Ray -> Maybe Double+  }+ -- | Produce an image from the given camera settings, world, and seed. raytrace :: ToRandom m => CameraSettings -> Geometry m Material -> StdGen -> A.Matrix D Color raytrace (CameraSettings {..}) (Geometry _ hitWorld) seed = let@@ -108,6 +135,21 @@   pixelU = across ^/ fromIntegral cs_imageWidth   pixelV = down ^/ fromIntegral imageHeight +  redirectTargets = flip map cs_redirectTargets $ \(_, q, s0, s1) -> RedirectTarget+    { rt_origin = q+    , rt_U = s0+    , rt_V = s1+    , rt_cross = cross s0 s1+    , rt_hit = +        let Geometry _ hit = parallelogram q s0 s1 in+        \ray -> fmap (hr_t . fst) (runIdentity (hit undefined ray (0, infinity)))+    }++  probs = map (\(p, _, _, _) -> p) cs_redirectTargets+  remProb = 1 - sum probs+  thresholds = zip redirectTargets (scanl1 (+) probs)+  getTarget r = fmap fst (listToMaybe (dropWhile ((r >=) . snd) thresholds))+   defocusRadius = cs_focusDist * tan (cs_defocusAngle / 2)   defocusDiskU = u ^* defocusRadius   defocusDiskV = v ^* defocusRadius@@ -127,24 +169,66 @@   getRay i j = do     origin <- sampleDefocusDisk     target <- samplePixel i j-    pure (Ray origin (target - origin))+    pure (Ray origin (normalize (target - origin))) -  rayColor :: Int -> Ray -> State StdGen Color-  rayColor depth ray+  rayColor :: Int -> Double -> Ray -> State StdGen Color+  rayColor depth time ray@(Ray _ rayDir)     | depth <= 0 = pure zero     | otherwise =-    toRandom (hitWorld ray (0.0001, infinity)) >>= \case+    toRandom (hitWorld time ray (0.0001, infinity)) >>= \case       Nothing -> pure (cs_background ray)-      Just (hit, Material mat) -> -        mat ray hit >>= \case-          (emitted, Nothing) -> pure emitted-          (emitted, Just (attenuation, ray')) -> do-            c <- rayColor (depth - 1) ray'-            pure (emitted + attenuation * c)+      Just (hit, Material mat) -> do+        let (emitted, genRes) = mat rayDir hit+        res <- genRes+        (emitted +) <$> case res of+          Absorb -> pure zero+          Scatter attenuation dir -> (attenuation *) <$> rayColor (depth - 1) time (Ray (hr_point hit) dir)++          HemisphereF matF -> do+            choice <- getTarget <$> state random+            dir <- case choice of+              Nothing -> do+                uu <- randomUnitVector+                pure (normalize (hr_normal hit + uu))+              Just RedirectTarget {..} -> do+                (i, j) <- state random+                let lightPt = rt_origin + i *^ rt_U + j *^ rt_V+                pure (normalize (lightPt - hr_point hit))+            let pdf1 = dot dir (hr_normal hit) / pi+            if pdf1 <= 0 then pure zero else do+              let pdfs = flip map redirectTargets $ \RedirectTarget {..} ->+                    case rt_hit (Ray (hr_point hit) dir) of+                      Nothing -> 0+                      Just t -> t * t / abs (dot rt_cross dir)+              -- The pdf from which 'dir' was generated+              let pdf = remProb * pdf1 + sum (zipWith (*) probs pdfs)+              c <- rayColor (depth - 1) time (Ray (hr_point hit) dir)+              pure (matF dir * c ^* (pdf1 / pdf))++          SphereF matF -> do+            choice <- getTarget <$> state random+            dir <- case choice of+              Nothing -> randomUnitVector+              Just RedirectTarget {..} -> do+                (i, j) <- state random+                let lightPt = rt_origin + i *^ rt_U + j *^ rt_V+                pure (normalize (lightPt - hr_point hit))+            let pdf1 = 0.25 / pi+            let pdfs = flip map redirectTargets $ \RedirectTarget {..} ->+                  case rt_hit (Ray (hr_point hit) dir) of+                    Nothing -> 0+                    Just t -> t * t / abs (dot rt_cross dir)+            -- The pdf from which 'dir' was generated+            let pdf = remProb * pdf1 + sum (zipWith (*) probs pdfs) +            c <- rayColor (depth - 1) time (Ray (hr_point hit) dir)+            pure (matF dir * c ^* (pdf1 / pdf))      pixelColor :: Int -> Int -> State StdGen Color   pixelColor i j = do-    colors <- replicateM cs_samplesPerPixel (getRay i j >>= rayColor cs_maxRecursionDepth)+    colors <- replicateM cs_samplesPerPixel $ do+      time <- state random+      ray <- getRay i j +      rayColor cs_maxRecursionDepth time ray     pure (sum colors ^/ fromIntegral cs_samplesPerPixel)    -- array of random seeds for each pixel (constructed using splitGen)
src/Graphics/Ray/Core.hs view
@@ -6,7 +6,7 @@     -- * Intervals   , Interval, inInterval, midpoint, padInterval     -- * Boxes-  , Box, fromCorners, boxJoin, boxHull, allCorners, padBox, longestDim, overlapsBox+  , Box, fromCorners, boxJoin, boxHull, allCorners, padBox, shiftBox, longestDim, overlapsBox     -- * Hit Records   , HitRecord(..)   ) where@@ -16,7 +16,7 @@ import Control.Monad.State (MonadState, state) import Control.Applicative (liftA2) import Data.Maybe (isJust)-import Data.Foldable (foldl')+import Data.List (foldl1')  -- | Floating-point infinity. infinity :: Double@@ -46,7 +46,6 @@   | otherwise = if y > z then Y else Z  -- | If @n@ is the normal vector of a mirror, and @v@ is an incoming light ray, then @reflect n v@ is the outgoing light ray.--- The first argument should be a unit vector, but the second need not be. reflect :: Vec3 -> Vec3 -> Vec3 reflect normal v =    v - 2 * dot normal v *^ normal@@ -68,20 +67,19 @@     then pure vec     else randomInUnitDisk --- | A ray with an origin and a direction. Points on the ray @Ray orig dir@ are parametrized by @orig + t *^ dir@.--- There is no expectation that the direction be a unit vector.+-- | A ray with an origin and a direction (which is always a unit vector). +-- Points on the ray @Ray orig dir@ are parametrized by @orig + t *^ dir@. data Ray = Ray Point3 Vec3   deriving (Show) --- | An interval with a lower bound and an upper bound. --- Variously interpreted as a closed interval or an open interval; it doesn't really matter.+-- | An interval with a lower bound and an upper bound. type Interval = (Double, Double)  -- [private] size :: Interval -> Double size (a, b) = b - a --- | Test whether a number is in the interval.+-- | Test whether a number is in the interval (interpreted as an open interval). inInterval :: Interval -> Double -> Bool inInterval (tmin, tmax) t = tmin < t && t < tmax @@ -114,16 +112,20 @@ fromCorners :: Point3 -> Point3 -> Box fromCorners = liftA2 (\x y -> if x < y then (x, y) else (y, x)) --- | The smallest box containing two boxes.-boxJoin :: Box -> Box -> Box-boxJoin = liftA2 (\(min1, max1) (min2, max2) -> (min min1 min2, max max1 max2))- -- | The smallest box containing all of the boxes.-boxHull :: [Box] -> Box-boxHull = foldl' boxJoin (V3 (infinity, -infinity) (infinity, -infinity) (infinity, -infinity))+boxJoin :: [Box] -> Box+boxJoin = foldl1' (liftA2 (\(min1, max1) (min2, max2) -> (min min1 min2, max max1 max2)))  +-- | The smallest box containing all of the points.+boxHull :: [Point3] -> Box+boxHull pts = let+  coords = sequence pts+  mins = fmap minimum coords+  maxs = fmap maximum coords+  in liftA2 (,) mins maxs+ -- | Get a list of all eight corners of a box.-allCorners :: Box -> [ Point3 ]+allCorners :: Box -> [Point3] allCorners (V3 i1 i2 i3) =    [ V3 (f1 i1) (f2 i2) (f3 i3)    | f1 <- [ fst, snd ], f2 <- [ fst, snd ], f3 <- [ fst, snd ]@@ -133,6 +135,10 @@ padBox :: Double -> Box -> Box padBox padding = fmap (padInterval padding) +-- | Translate the box by the given vector.+shiftBox :: Vec3 -> Box -> Box+shiftBox = liftA2 (\x (a, b) -> (a + x, b + x))+ -- | The dimension in which the box is the longest. longestDim :: Box -> Dim longestDim = argMax . fmap size@@ -153,3 +159,4 @@   , hr_frontSide :: Bool -- ^ Whether the ray hit the "front side" (for a closed surface, the outside).   , hr_uv :: V2 Double -- ^ The texture coordinates of the intersection.   }+  deriving (Show)
src/Graphics/Ray/Geometry.hs view
@@ -1,21 +1,26 @@ {-# LANGUAGE InstanceSigs #-} {-# LANGUAGE LambdaCase #-} {-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE ViewPatterns #-} module Graphics.Ray.Geometry    ( -- * Geometry-    Geometry(Geometry), pureGeometry, boundingBox+    Geometry(Geometry), boundingBox, pureGeometry, transform, moving     -- * Surfaces and Volumes-  , sphere, parallelogram, cuboid, constantMedium+  , sphere, planeShape, parallelogram, cuboid, triangle, triangleMesh, constantMedium+    -- * Polygonal Meshes+  , Mesh(Mesh), transformVertices, parseObj, readObj     -- * Groups-  , group, bvhNode, Tree(Leaf, Node), bvhTree, autoTree+  , group, bvhNode, bvhTree     -- * Transformations-  , transform, translate, rotateX, rotateY, rotateZ+  , translate, rotateX, rotateY, rotateZ, scale   ) where  import Graphics.Ray.Core  import Linear (V2(V2), V3(V3), dot, quadrance, (*^), (^/), cross, norm, M44, inv44, (!*), V4(V4)) import qualified Linear.V4 as V4+import Data.Massiv.Array (U, (!))+import qualified Data.Massiv.Array as A import System.Random (StdGen, random) import Control.Monad.State (State, state) import Control.Monad (guard, foldM)@@ -24,22 +29,26 @@ import Data.Bifunctor (first, second) import Data.Functor.Identity (Identity(Identity), runIdentity) import Data.Functor ((<&>))+import Text.Read (readMaybe)+import Data.Char (isDigit)+import GHC.Base (inline)  -- | A @'Geometry' m a@ has a bounding box (used in the implementation of bounding volume hierarchies),--- as well as a function that takes a ray and an interval, and in the @m@ monad, produces either @Nothing@+-- as well as a function that takes a time, a ray, and an open interval, and in the @m@ monad, produces either @Nothing@ -- (if the ray does not intersect the shape within that interval) or a tuple consisting of a 'HitRecord' and a value of type @a@.--- Typically, @m@ is either 'Identity' or @'State' 'StdGen'@, and @a@ is either @()@ or 'Geometry.Material.Material'. Use the '(<$)' operator--- to add a material to a geometry.-data Geometry m a = Geometry Box (Ray -> Interval -> m (Maybe (HitRecord, a)))+-- The time parameter is between 0 and 1 and allows for motion blur effects.+-- Typically, @m@ is either 'Identity' or @'State' 'StdGen'@, and @a@ is either @()@ or 'Geometry.Material.Material'. +-- Use the '(<$)' operator to add a material to a geometry.+data Geometry m a = Geometry Box (Double -> Ray -> Interval -> m (Maybe (HitRecord, a)))  instance Functor m => Functor (Geometry m) where   {-# SPECIALISE fmap :: (a -> b) -> Geometry Identity a -> Geometry Identity b #-}   fmap :: (a -> b) -> Geometry m a -> Geometry m b-  fmap f (Geometry bbox hit) = Geometry bbox (\ray ival -> fmap (fmap (second f)) (hit ray ival))+  fmap f (Geometry bbox hit) = Geometry bbox (\ray time ival -> fmap (fmap (second f)) (hit ray time ival))  -- | Promote a pure geometry to a monadic one. pureGeometry :: Applicative m => Geometry Identity a -> Geometry m a-pureGeometry (Geometry bbox f) = Geometry bbox (\ray ival -> pure (runIdentity (f ray ival)))+pureGeometry (Geometry bbox f) = Geometry bbox (\ray time ival -> pure (runIdentity (f ray time ival)))  -- | Get a geometry's bounding box. boundingBox :: Geometry m a -> Box@@ -51,18 +60,17 @@   diag = V3 radius radius radius   bbox = fromCorners (center - diag) (center + diag) -  hitSphere (Ray orig dir) bounds = Identity $ do+  hitSphere _ (Ray orig dir) bounds = Identity $ do     let oc = center - orig-    let a = quadrance dir     let h = dot dir oc      let c = quadrance oc - radius*radius -    let discriminant = h*h - a*c+    let discriminant = h*h - c     guard (discriminant >= 0)          let sqrtd = sqrt discriminant-    let root1 = (h - sqrtd) / a-    let root2 = (h + sqrtd) / a+    let root1 = h - sqrtd+    let root2 = h + sqrtd      t <-        if inInterval bounds root1 @@ -95,30 +103,33 @@     u = atan2 x z / (2 * pi) + 0.5     v = acos (-y) / pi  --- | Construct a parallelogram from a corner point and two edge vectors.+-- | Construct a subset of a plane. (See 'parallelogram' and 'triangle' for two instances of this.) -- Which side is the \"front side\" is determined by the right hand rule.-parallelogram :: Point3 -> Vec3 -> Vec3 -> Geometry Identity ()-parallelogram q u v = let+{-# INLINABLE planeShape #-}+planeShape +  :: Point3 -- ^ (0, 0) point on the plane+  -> Vec3 -- ^ First basis vector+  -> Vec3 -- ^ Second basis vector+  -> (Double -> Double -> Bool) -- ^ Whether a point on the plane is in the shape+  -> (Double -> Double -> V2 Double) -- ^ Texture coordinates+  -> Box -- ^ Bounding box (this is padded by a small amount to ensure all dimensions are positive)+  -> Geometry Identity ()+planeShape q u v test getUV bbox = let   cp = cross u v-  area = norm cp-  normal = cp ^/ area-  normalS = normal ^/ area-  n_dot_q = dot normal q+  norm_cp = norm cp+  normal = cp ^/ norm_cp+  normalS = normal ^/ norm_cp -  box1 = fromCorners q (q + u + v)-  box2 = fromCorners (q + u) (q + v)-  bbox = padBox 0.0001 (boxJoin box1 box2)-  -  hitParallelogram (Ray orig dir) bounds = Identity $ do+  hitShape _ (Ray orig dir) bounds = Identity $ do     let denom = dot normal dir     guard (abs denom > 1e-8)-    let t = (n_dot_q - dot normal orig) / denom+    let t = dot normal (q - orig) / denom     guard (inInterval bounds t)     let p = orig + t *^ dir     let p_rel = p - q     let a = normalS `dot` (p_rel `cross` v)     let b = normalS `dot` (u `cross` p_rel)-    guard (0 <= a && a <= 1 && 0 <= b && b <= 1)+    guard (test a b)     let frontSide = denom < 0      let hit = HitRecord@@ -126,12 +137,19 @@           , hr_point = p           , hr_normal = if frontSide then normal else -normal           , hr_frontSide = frontSide-          , hr_uv = V2 a b+          , hr_uv = getUV a b           }     Just (hit, ()) -  in Geometry bbox hitParallelogram +  in Geometry (padBox 0.0001 bbox) hitShape +-- | Construct a parallelogram from a corner point and two edge vectors.+parallelogram :: Point3 -> Vec3 -> Vec3 -> Geometry Identity ()+parallelogram q u v = let+  bbox = boxHull [ q, q + u, q + v, q + u + v ] +  test a b = 0 <= a && a <= 1 && 0 <= b && b <= 1+  in inline planeShape q u v test V2 bbox+ -- | Construct an axis-aligned rectangular cuboid (implemented as a 'group' of parallelograms). cuboid :: Box -> Geometry Identity () cuboid (V3 (xmin, xmax) (ymin, ymax) (zmin, zmax)) = let@@ -147,36 +165,165 @@     , parallelogram (V3 xmin ymin zmin) dx dz -- bottom     ] --- | Construct a constant-density medium (like fog or smoke). --- Typical materials are 'Graphics.Material.isotropic' and 'Graphics.Material.pitchBlack'.+-- | Construct a triangle from three corner points and their texture coordinates.+triangle :: (Point3, V2 Double) -> (Point3, V2 Double) -> (Point3, V2 Double) -> Geometry Identity ()+triangle (p0, uv0) (p1, uv1) (p2, uv2) = let+  s1 = p1 - p0+  s2 = p2 - p0+  bbox = boxHull [ p0, p1, p2 ]+  test a b = a >= 0 && b >= 0 && a + b <= 1+  getUV a b = (1 - a - b) *^ uv0 + a *^ uv1 + b *^ uv2+  in inline planeShape p0 s1 s2 test getUV bbox++-- | A collection of triangles.+data Mesh = Mesh +  (A.Vector U Point3) -- ^ Array of vertex locations.+  (A.Vector U (V2 Double)) -- ^ Array of texture coordinates.+  [V3 (Int, Maybe Int)] -- ^ List of triangles. Each triangle has three vertices, whose locations and (optional) texture coordinates+                        -- are defined by indexing into the two arrays.+  deriving (Show)++-- | Apply an affine transformation (represented as a 4 by 4 matrix whose bottom row is 0 0 0 1) to the vertices of a mesh.+transformVertices :: M44 Double -> Mesh -> Mesh+transformVertices m (Mesh vs vts fs) = +  let vs' = A.compute (A.map (dropLast . (m !*) . V4.point) vs) in+  Mesh vs' vts fs++-- TODO: use IO exceptions instead of Either?+-- | Parse the .obj file at the given location.+readObj :: FilePath -> IO (Either String Mesh) +readObj path = first ((path ++ ", ") ++) . parseObj <$> readFile path++-- | Parse a Wavefront .obj file.+--+-- * Comments beginning with @#@ are ignored, as are all lines that do not+--   begin with @v @, @vt @, or @f @. +-- * Faces with more than three vertices are allowed, but are triangulated before+--   adding them to the 'Mesh' object.+-- * Indices can be either positive or negative. For example, if there are 20+--   @v@ statements in the file, the vertex indices of a face must be in the+--   range [1, 20] or [-20, -1], with -1 meaning the last vertex. In either case,+--   they are converted into 0-based indices.+parseObj :: String -> Either String Mesh+parseObj file = do+  let (vLines, vtLines, fLines) = partitionLines (removeComments file)+  vs <- mapM parseV vLines+  vts <- mapM parseVT vtLines+  fs <- concat <$> mapM (parseF (length vs) (length vts)) fLines+  Right (Mesh (A.fromList A.Seq vs) (A.fromList A.Seq vts) fs)+  +  where+    removeComments :: String -> [String]+    removeComments = map (takeWhile (/= '#')) . lines++    partitionLines :: [String] -> ([(Int, String)], [(Int, String)], [(Int, String)])+    partitionLines = foldr addLine ([], [], []) . zip [1..]++    addLine (k, line) (vs, vts, fs) =+      case line of+        'v':' ':rest -> ((k, rest) : vs, vts, fs)+        'v':'t':' ':rest -> (vs, (k, rest) : vts, fs)+        'f':' ':rest -> (vs, vts, (k, rest) : fs)+        _ -> (vs, vts, fs)++    withLine :: Int -> String -> String+    withLine k err = "line " ++ show k ++ ": " ++ err+    +    -- a 'v' statement must begin with three decimal numbers+    parseV (k, line) = +      case words line of+        (readMaybe -> Just x) : (readMaybe -> Just y) : (readMaybe -> Just z) : _ -> Right (V3 x y z)+        _ -> Left (withLine k "invalid 'v' statement")+    +    -- a 'vt' statement must begin with two decimal numbers (or consist of a single decimal number, in which case v defaults to 0)+    parseVT (k, line) =+      case words line of+        [readMaybe -> Just u] -> Right (V2 u 0)+        (readMaybe -> Just u) : (readMaybe -> Just v) : _ -> Right (V2 u v)+        _ -> Left (withLine k "invalid 'vt' statement")++    parseF numVs numVTs (k, line) = +      case mapM (getIndices numVs numVTs) (words line) of+          Left err -> Left (withLine k err)+          Right (i:is) | length is >= 2 -> Right (map (uncurry (V3 i)) (pairs is))+          Right _ -> Left (withLine k "invalid 'f' statement (fewer than 3 vertices)")++    pairs :: [a] -> [(a, a)]+    pairs = \case+      [] -> []+      [_] -> []+      x:xs@(y:_) -> (x, y) : pairs xs++    processIx len i +      | 1 <= i && i <= len = Right (i - 1)+      | -len <= i && i <= -1 = Right (i + len)+      | otherwise = Left ("index out of bounds: " ++ show i)++    getIndices :: Int -> Int -> String -> Either String (Int, Maybe Int)+    getIndices numVs numVTs str = do+      (i, rest) <- extractInt str+      i' <- processIx numVs i+      case rest of+        "" -> Right (i', Nothing)+        '/':'/':_ -> Right (i', Nothing)+        '/':str' -> do+          (j, _) <- extractInt str'+          j' <- processIx numVTs j+          Right (i', Just j')+        c:_ -> Left ("unexpected character '" ++ c : "'")++    extractInt :: String -> Either String (Int, String)+    extractInt = \case+      '-':str -> first negate <$> extractNat str+      str -> extractNat str+    +    extractNat :: String -> Either String (Int, String)+    extractNat str = +      let (ds, rest) = span isDigit str in+      case readMaybe ds of+        Nothing -> Left "expected number"+        Just i -> Right (i, rest)++-- | Realize a 'Mesh' as a 'Geometry' (implemented as a 'bvhTree' of triangles).+triangleMesh :: Mesh -> Geometry Identity ()+triangleMesh (Mesh verts uvs tris) = +  bvhTree $ flip map tris $ \(V3 (i0, j0) (i1, j1) (i2, j2)) -> let+    uv0 = maybe (V2 0 0) (uvs !) j0+    uv1 = maybe (V2 1 0) (uvs !) j1+    uv2 = maybe (V2 0 1) (uvs !) j2+    in triangle (verts ! i0, uv0) (verts ! i1, uv1) (verts ! i2, uv2)++-- | Construct a constant-density medium; useful for subsurface scattering and fog effects.+-- Typical materials are 'Graphics.Material.isotropic', 'Graphics.Material.anisotropic', and 'Graphics.Material.pitchBlack'. constantMedium   :: Double -- ^ Density -  -> Geometry Identity () -- ^ Surface (assumed to be convex in current implementation)+  -> Geometry Identity () -- ^ Surface (assumed to be a closed surface with the \"front side\" facing outwards)   -> Geometry (State StdGen) () constantMedium density (Geometry bbox hitObj) = let   negInvDensity = -(1 / density) -  hitMedium :: Ray -> Interval -> State StdGen (Maybe (HitRecord, ()))-  hitMedium ray@(Ray orig dir) (tmin, tmax) = -    case do (hit1, ()) <- runIdentity (hitObj ray (-infinity, infinity))-            (hit2, ()) <- runIdentity (hitObj ray (hr_t hit1, infinity))-            let t1 = max tmin (hr_t hit1)-            let t2 = min tmax (hr_t hit2)-            guard (t1 < t2)-            Just (t1, t2, hit1) of+  hitMedium :: Double -> Ray -> Interval -> State StdGen (Maybe (HitRecord, ()))+  hitMedium time ray@(Ray orig dir) (tmin, tmax) = +    case do (hit1, ()) <- runIdentity (hitObj time ray (tmin, infinity)) +            if hr_frontSide hit1 +              then do+                guard (hr_t hit1 < tmax)+                (hit2, ()) <- runIdentity (hitObj time ray (hr_t hit1, infinity))+                Just (hr_t hit1, min tmax (hr_t hit2))+              else Just (tmin, min tmax (hr_t hit1))+    of       Nothing -> pure Nothing -- ray is never in fog within interval-      Just (t1, t2, hit1) -> state random <&> \rand ->-         do let rayScale = norm dir-            let inDist = (t2 - t1) * rayScale+      Just (t1, t2) -> state random <&> \rand ->+         do let inDist = t2 - t1             let hitDist = negInvDensity * log rand             guard (hitDist < inDist)-            let t = t1 + hitDist / rayScale+            let t = t1 + hitDist             let hit = HitRecord                   { hr_t = t                   , hr_point = orig + t *^ dir-                  , hr_normal = hr_normal hit1-                  , hr_frontSide = hr_frontSide hit1-                  , hr_uv = hr_uv hit1+                  , hr_normal = -dir -- arbitrary+                  , hr_frontSide = True+                  , hr_uv = V2 0 0 -- arbitrary                   }             Just (hit, ()) @@ -188,11 +335,11 @@ {-# SPECIALISE group :: [Geometry Identity a] -> Geometry Identity a #-} group :: Monad m => [Geometry m a] -> Geometry m a group obs = let-  bbox = boxHull (map boundingBox obs)-  -  hitGroup ray (tmin, tmax) =+  bbox = boxJoin (map boundingBox obs)++  hitGroup ray time (tmin, tmax) =     let try (tmax', knownHit) (Geometry _ hitObj) =-          hitObj ray (tmin, tmax') <&> \case+          hitObj ray time (tmin, tmax') <&> \case             Nothing -> (tmax', knownHit)             Just (hit, mat) -> (hr_t hit, Just (hit, mat))     in snd <$> foldM try (tmax, Nothing) obs@@ -204,49 +351,43 @@ {-# SPECIALISE bvhNode :: Geometry Identity a -> Geometry Identity a -> Geometry Identity a #-} bvhNode :: Monad m => Geometry m a -> Geometry m a -> Geometry m a bvhNode (Geometry bboxLeft hitLeft) (Geometry bboxRight hitRight) = let-  bbox = boxJoin bboxLeft bboxRight+  bbox = boxJoin [bboxLeft, bboxRight] -  hitBvhNode ray (tmin, tmax)+  hitBvhNode time ray (tmin, tmax)     | overlapsBox bbox ray (tmin, tmax) = -      hitLeft ray (tmin, tmax) >>= \case-        Nothing -> hitRight ray (tmin, tmax)-        res@(Just (hit, _)) -> fmap (<|> res) (hitRight ray (tmin, hr_t hit))+      hitLeft time ray (tmin, tmax) >>= \case+        Nothing -> hitRight time ray (tmin, tmax)+        res@(Just (hit, _)) -> fmap (<|> res) (hitRight time ray (tmin, hr_t hit))     | otherwise = pure Nothing      in Geometry bbox hitBvhNode -data Tree a = Leaf a | Node (Tree a) (Tree a)---- | Group multiple geometric objects (organized as a tree) into a single object. A bounding box is created for every subtree of the --- given tree; if a ray does not intersect the bounding box, it cannot hit any of the child objects, so none of--- them need to be tested further.-{-# SPECIALISE bvhTree :: Tree (Geometry Identity a) -> Geometry Identity a #-}-bvhTree :: Monad m => Tree (Geometry m a) -> Geometry m a+-- | Group multiple geometric objects into a single object. +-- The objects are organized into a tree based on their positions, and then a 'bvhNode'+-- is created for each node in the tree.+{-# SPECIALISE bvhTree :: [Geometry Identity a] -> Geometry Identity a #-}+bvhTree :: Monad m => [Geometry m a] -> Geometry m a bvhTree = \case-  Leaf a -> a-  Node left right -> bvhNode (bvhTree left) (bvhTree right)---- | Organize the geometric objects into a tree based on their positions.-autoTree :: [Geometry m a] -> Tree (Geometry m a)-autoTree = \case-  [] -> error "autoTree: empty list"-  [obj] -> Leaf obj+  [] -> error "bvhTree: empty list"+  [obj] -> obj   obs -> let-    d = longestDim (boxHull (map boundingBox obs))+    d = longestDim (boxJoin (map boundingBox obs))     obs' = sortOn (midpoint . component d . boundingBox) obs     (left, right) = splitAt (length obs `div` 2) obs'-    in Node (autoTree left) (autoTree right)+    in bvhNode (bvhTree left) (bvhTree right)  -- | Apply an affine transformation (represented as a 4 by 4 matrix whose bottom row is 0 0 0 1) to a geometric object.+-- The transformation should be Euclidean (translation, rotation, reflection, or a composition thereof); otherwise, the+-- normal vectors of the resulting geometry will be incorrect. transform :: Functor m => M44 Double -> Geometry m a -> Geometry m a transform m (Geometry bbox hitObj) = let   m34 = dropLast m   inv_m = dropLast (inv44 m)-  cornerCoords = mapM ((m34 !*) . V4.point) (allCorners bbox) :: V3 [Double]-  bbox' = fromCorners (fmap minimum cornerCoords) (fmap maximum cornerCoords)-  in Geometry bbox' $ \(Ray orig dir) ival ->+  corners = map ((m34 !*) . V4.point) (allCorners bbox)+  bbox' = boxHull corners+  in Geometry bbox' $ \time (Ray orig dir) ival ->     let ray' = Ray (inv_m !* V4.point orig) (inv_m !* V4.vector dir) in-    flip (fmap . fmap . first) (hitObj ray' ival) $ \hit@(HitRecord {..}) ->+    flip (fmap . fmap . first) (hitObj time ray' ival) $ \hit@(HitRecord {..}) ->       hit { hr_point = m34 !* V4.point hr_point, hr_normal = m34 !* V4.vector hr_normal }  -- | Translation.@@ -290,6 +431,26 @@     c = cos angle     s = sin angle +-- | Scaling centered at the origin. This should not be used with 'transform'.+scale :: Double -> M44 Double+scale a = V4+  (V4 a 0 0 0)+  (V4 0 a 0 0)+  (V4 0 0 a 0)+  (V4 0 0 0 1)+ -- [private] dropLast :: V4 a -> V3 a dropLast (V4 x y z _) = V3 x y z++-- TODO: this produces incorrect results for objects with solid textures+-- | Create a motion-blurred object that is translated by the first vector at time 0+-- and by the second vector at time 1.+moving :: Functor m => Vec3 -> Vec3 -> Geometry m a -> Geometry m a+moving v0 v1 (Geometry bbox hitObj) =+  let bbox' = boxJoin [shiftBox v0 bbox, shiftBox v1 bbox] in+  Geometry bbox' $ \time (Ray orig dir) ival -> let+    shift = (1 - time) *^ v0 + time *^ v1+    ray' = Ray (orig - shift) dir+    in flip (fmap . fmap . first) (hitObj time ray' ival) $ \hit -> +      hit { hr_point = hr_point hit + shift }
src/Graphics/Ray/Material.hs view
@@ -1,70 +1,96 @@ {-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE TupleSections #-} module Graphics.Ray.Material -  ( Material(Material)-  , lightSource, pitchBlack, lambertian, mirror, metal, dielectric, transparent, isotropic+  ( Material(..), MaterialResult(..)+  , lightSource, pitchBlack, lambertian, lommelSeeliger, mirror, metal, dielectric, transparent, isotropic, anisotropic   ) where  import Graphics.Ray.Core import Graphics.Ray.Texture -import Linear (V3(V3), zero, normalize, dot, quadrance, (*^))+import Linear (V3(V3), zero, dot, quadrance, (*^), normalize) import System.Random (StdGen, random) import Control.Monad.State (State, state) --- | A material is a function that, given the details of a ray-surface intersection, produces an emitted color--- (which is @0@ for all of the materials in this module save 'lightSource') and potentially a reflected ray--- along with a color to scale the result of tracing said ray by.-newtype Material = Material (Ray -> HitRecord -> State StdGen (Color, Maybe (Color, Ray)))+-- | A material is a function that takes the details of a ray-surface intersection (specifically, the direction of the+-- ray, along with a 'HitRecord' object) and produces an emitted color as well as a 'MaterialResult' (see below).+newtype Material = Material (Vec3 -> HitRecord -> (Color, State StdGen MaterialResult)) --- | A material that emits light and does not reflect rays.+-- | The ray tracer uses the 'MaterialResult' to compute the color from scattering, which is added to the color from+-- emission.+data MaterialResult+  = Absorb +    -- ^ Do not create a new ray.+  | Scatter Color Vec3 +    -- ^ Create a new ray in the given direction (which must be a unit vector)+    -- and multiply the result of recursively ray tracing the ray by the given color.+  | HemisphereF (Vec3 -> Color) +    -- ^ Leave the task of computing a direction to the ray tracer (allowing it to send more rays toward light sources).+    -- The function is the BRDF of the material times \(\pi\). Its argument (the new ray direction) is guaranteed to be a unit+    -- vector and to have a positive dot product with 'hr_normal'.+  | SphereF (Vec3 -> Color) +    -- ^ Similar to 'HemisphereF', but the new ray direction can be anywhere on the sphere. (Such materials are typically used+    -- with volumes like 'Graphics.Geometry.constantMedium'.) The function is the albedo times the material's phase function+    -- times \(4\pi\). Its argument is guaranteed to be a unit vector.++-- NOTE: I could generalize 'HemisphereF' and 'SphereF' by allowing the material to specify a unit vector generator (and pdf)+-- to use in the case of no redirection. 'anisotropic' and 'lommelSeeliger' could benefit from this. One issue is that it would require+-- converting a single vector into an orthonormal basis.++-- | A material that emits light equally in all directions and does not scatter rays. lightSource :: Texture -> Material-lightSource (Texture tex) = Material $-  \_ (HitRecord {..}) -> pure (tex hr_point hr_uv, Nothing)+lightSource (Texture tex) = Material $ \ _ HitRecord{..} -> (tex hr_point hr_uv, pure Absorb) --- | A material that absorbs all light. Identical in principle to @'lambertian' ('constantTexture' 0)@, @'isotropic' ('constantTexture' 0)@,--- and so on, but avoids unecessary computation.+-- | A material that absorbs all light. This is the same as @'lightSource' ('constantTexture' 0)@, but the name reflects+-- the fact that it is not actually a light source. pitchBlack :: Material-pitchBlack = Material $ \ _ _ -> pure (zero, Nothing)+pitchBlack = Material $ \ _ _ -> (zero, pure Absorb)  -- | A material that exhibits Lambertian reflectance. The direction of the reflected ray is independent of the direction -- of the incoming ray, resulting in a diffuse (non-shiny) appearance. lambertian :: Texture -> Material lambertian (Texture tex) = Material $-  \_ (HitRecord {..}) -> do-    u <- randomUnitVector-    pure (zero, Just (tex hr_point hr_uv, Ray hr_point (hr_normal + 0.9999 *^ u)))+  \_ (HitRecord {..}) -> (zero, pure $ HemisphereF $ const (tex hr_point hr_uv)) +-- | A material that exhibits Lommel-Seeliger reflectance.+lommelSeeliger :: Texture -> Material+lommelSeeliger (Texture tex) = Material $+  \inDir (HitRecord {..}) -> (zero,) $ pure $ HemisphereF $ \outDir -> let+    mu0 = -(dot inDir hr_normal)+    mu1 = dot outDir hr_normal+    in 0.25 / (mu0 + mu1) *^ tex hr_point hr_uv + -- | A colored mirror. (For no color, use @'constantTexture' 1@.) mirror :: Texture -> Material mirror (Texture tex) = Material $-  \(Ray _ dir) (HitRecord {..}) -> -    pure (zero, Just (tex hr_point hr_uv, Ray hr_point (reflect hr_normal dir)))+  \dir (HitRecord {..}) -> +    (zero, pure $ Scatter (tex hr_point hr_uv) (reflect hr_normal dir)) +-- TODO: make this more physically-based? -- | A metallic-looking material that reflects rays inexactly. The larger the first argument is, the less shiny -- the material. @'metal' 0@ behaves the same as 'mirror'. metal :: Double -> Texture -> Material metal fuzz (Texture tex) = Material $-  \(Ray _ dir) (HitRecord {..}) -> do+  \dir (HitRecord {..}) -> (zero,) $ do     u <- randomUnitVector-    let dir' = normalize (reflect hr_normal dir) + (fuzz *^ u)+    let dir' = reflect hr_normal dir + fuzz *^ u     let scatter = dot dir' hr_normal > 0-    pure (zero, if scatter then Just (tex hr_point hr_uv, Ray hr_point dir') else Nothing)+    pure (if scatter then Scatter (tex hr_point hr_uv) (normalize dir') else Absorb)  -- [private] refract :: Double -> Double -> Vec3 -> Vec3 -> Vec3  refract iorRatio cosTheta normal u = let   perp = iorRatio *^ (u + cosTheta *^ normal) -  para = -(sqrt (abs (1 - quadrance perp)) *^ normal)-  in perp + para+  para = -(sqrt (abs (1 - quadrance perp)) *^ normal) -- NOTE: quadrance perp = (iorRatio * sinTheta)^2+  in perp + para -- unit vector  -- | A material that either reflects or refracts all incoming rays, like clear glass. -- The argument is the index of refraction relative to the surrounding medium. dielectric :: Double -> Material dielectric ior = Material $-  \(Ray _ dir) (HitRecord {..}) -> do+  \dir (HitRecord {..}) -> (zero,) $ do     let iorRatio = if hr_frontSide then 1/ior else ior-    let u = normalize dir-    let cosTheta = min 1 (dot hr_normal (-u))+    let cosTheta = min 1 (dot hr_normal (-dir))     let sinTheta = sqrt (1 - cosTheta * cosTheta)     let cannotRefract = iorRatio * sinTheta > 1 @@ -74,21 +100,30 @@     x <- state random      let dir' = if cannotRefract || x < reflectance-          then reflect hr_normal u-          else refract iorRatio cosTheta hr_normal u+          then reflect hr_normal dir+          else refract iorRatio cosTheta hr_normal dir     -    pure (zero, Just (V3 1 1 1, Ray hr_point dir'))    +    pure (Scatter (V3 1 1 1) dir') --- | A material that lets all light through, with the given tint.+-- | A material that lets light through, with the given tint. transparent :: Texture -> Material transparent (Texture tex) = Material $-  \(Ray _ dir) (HitRecord {..}) ->-    pure (zero, Just (tex hr_point hr_uv, Ray hr_point dir))+  \dir (HitRecord {..}) ->+    (zero, pure $ Scatter (tex hr_point hr_uv) dir)  -- | A material that scatters an incoming ray in a direction chosen uniformly at random from the unit sphere.  -- (Typically used with 'Graphics.Geometry.constantMedium'.) isotropic :: Texture -> Material isotropic (Texture tex) = Material $-  \_ (HitRecord {..}) -> do-    dir <- randomUnitVector-    pure (zero, Just (tex hr_point hr_uv, Ray hr_point dir))+  \_ (HitRecord {..}) -> (zero, pure $ SphereF $ const (tex hr_point hr_uv))++-- | A material that scatters an incoming ray according to the Henyey-Greenstein distribution. The first parameter+-- should be in the range (-1, 1); negative values result in more backward scattering and positive values result in+-- more forward scattering.+-- (Typically used with 'Graphics.Geometry.constantMedium'.)+anisotropic :: Double -> Texture -> Material+anisotropic g (Texture tex) = Material $+  \inDir (HitRecord {..}) -> (zero,) $ pure $ SphereF $ \outDir -> let+    mu = dot inDir outDir+    hg = (1 - g*g) / (1 + g*g - 2*g*mu)**1.5 +    in hg *^ tex hr_point hr_uv
test/Main.hs view
@@ -1,23 +1,32 @@ {-# LANGUAGE DataKinds #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE PartialTypeSignatures #-}-{-# LANGUAGE ViewPatterns #-}+{-# LANGUAGE LambdaCase #-}+{-# LANGUAGE TupleSections #-} module Main where+-- TODO: add unit tests?  import Graphics.Ray -import Linear (V3(V3), (*^), normalize, norm, (!*!))+import Linear (V3(V3), (*^), norm, (!*!)) import System.Random (StdGen, newStdGen, randomR, random, mkStdGen) import Control.Monad.State (State, runState, state) import Control.Monad (forM, replicateM) import Control.Applicative (liftA2) import Data.Functor.Identity (Identity)+import Data.Either (fromRight)  sky :: Ray -> Color-sky (Ray _ (normalize -> V3 _ y _)) = +sky (Ray _ (V3 _ y _)) =    let a = 0.5 * (y + 1) in   (1 - a) *^ V3 1 1 1 + a *^ V3 0.5 0.7 1 +grayFade :: Ray -> Color+grayFade (Ray _ dir) = let+  y = component Y dir+  t = (y + 1) * 0.5+  in t *^ 1+ metalTest :: IO () metalTest = let   materialGround = lambertian (constantTexture (V3 0.8 0.8 0.0))@@ -140,7 +149,7 @@    genWorld :: State StdGen (Geometry Identity Material)   genWorld = do-    fmap (bvhTree . autoTree . (bigSpheres ++) . concat) $ forM (liftA2 (,) [-11..10] [-11..10]) $ \(a, b) -> do+    fmap (bvhTree . (bigSpheres ++) . concat) $ forM (liftA2 (,) [-11..10] [-11..10]) $ \(a, b) -> do       offsetX <- state (randomR (0, 0.9))       offsetZ <- state (randomR (0, 0.9))       let center = V3 (a + offsetX) 0.2 (b + offsetZ)@@ -203,9 +212,10 @@     , cs_vfov = degrees 40     , cs_center = V3 278 278 (-800)     , cs_lookAt = V3 278 278 0+    , cs_redirectTargets = [ (0.25, V3 343 554 332, V3 (-130) 0 0, V3 0 0 (-105)) ]     } -  in writeImageSqrt "cornell_box.png" . raytrace settings world =<< newStdGen+  in writeImageSqrt "cornell_box.png" $ raytrace settings world (mkStdGen 234)  cornellSmoke :: IO () cornellSmoke = let@@ -241,6 +251,7 @@     , cs_vfov = degrees 40     , cs_center = V3 278 278 (-800)     , cs_lookAt = V3 278 278 0+    , cs_redirectTargets = [ (0.25, V3 113 554 127, V3 330 0 0, V3 0 0 305) ]     }    in writeImageSqrt "cornell_smoke.png" . raytrace settings world =<< newStdGen@@ -252,7 +263,7 @@    generateBoxes :: State StdGen (Geometry Identity Material)   generateBoxes = -    fmap ((ground <$) . bvhTree . autoTree) $ +    fmap ((ground <$) . bvhTree) $      forM (liftA2 (,) [0..19] [0..19]) $ \(i, j) -> do       let x0 = -1000 + i * 100       let z0 = -1000 + j * 100@@ -264,16 +275,17 @@      generateBalls :: State StdGen (Geometry Identity Material)   generateBalls =-    fmap ((white <$) . transform (translate (V3 (-100) 270 395) !*! rotateY (degrees 15)) . bvhTree . autoTree) $+    fmap ((white <$) . transform (translate (V3 (-100) 270 395) !*! rotateY (degrees 15)) . bvhTree) $     replicateM 1000 $ do       p <- state (randomR (0, 165))        pure (sphere p 10)      boundary = sphere (V3 360 150 145) 70 +  light f = f (V3 123 554 147) (V3 300 0 0) (V3 0 0 265)      largeObjects earth =-    [ lightSource (constantTexture (V3 7 7 7)) <$ parallelogram (V3 123 554 147) (V3 300 0 0) (V3 0 0 265)-    , lambertian (constantTexture (V3 0.7 0.3 0.1)) <$ sphere (V3 415 400 200) 50+    [ lightSource (constantTexture (V3 7 7 7)) <$ light parallelogram+    , lambertian (constantTexture (V3 0.7 0.3 0.1)) <$ moving 0 (V3 30 0 0) (sphere (V3 400 400 200) 50)     , dielectric 1.5 <$ sphere (V3 260 150 45) 50     , dielectric 1.5 <$ boundary     , metal 1.0 (constantTexture (V3 0.8 0.8 0.9)) <$ sphere (V3 0 150 145) 50@@ -299,21 +311,81 @@     , cs_samplesPerPixel = samplesPerPixel     , cs_maxRecursionDepth = maxRecursionDepth     , cs_background = const 0+    , cs_redirectTargets = [ light (0.25,,,) ]     }    in do     earth <- readImage "images/earthmap.jpg"-    seed <- newStdGen+    let seed = mkStdGen 1234     let (world, seed') = runState (generateWorld earth) seed     writeImageSqrt path (raytrace settings world seed') --- This should take less than 110 seconds+pawnTest :: IO ()+pawnTest = let+  world mesh = +    let pawn = triangleMesh mesh in group+      [ pureGeometry (dielectric 1.5 <$ pawn)+      , isotropic (constantTexture (V3 1 0 0)) <$ constantMedium 5 pawn+      ]+  settings = defaultCameraSettings +    { cs_center = V3 0 3.75 5+    , cs_lookAt = V3 0 2.75 0+    , cs_imageWidth = 500+    , cs_vfov = degrees 80+    , cs_samplesPerPixel = 400+    , cs_maxRecursionDepth = 20+    , cs_background = grayFade+    }+  in+  readObj "images/pawn.obj" >>= \case+    Left err -> putStrLn err+    Right mesh -> +      let mesh' = transformVertices (scale 100) mesh in+      writeImage "pawn_demo.png" (raytrace settings (world mesh') (mkStdGen 55))++lommelSeeligerTest :: IO ()+lommelSeeligerTest = let+  world = group+    [ lommelSeeliger (constantTexture 1) <$ sphere (V3 0 0 (-2)) 1+    , lightSource (constantTexture 160) <$ sphere (V3 0 0 22) 1+    ]++  settings = defaultCameraSettings+    { cs_imageWidth = 500+    , cs_samplesPerPixel = 500+    , cs_background = const 0+    , cs_redirectTargets = [ (0.5, V3 (-1) (-1) 21, V3 2 0 0, V3 0 2 0) ]+    }++  in writeImage "test_image.png" (raytrace settings world (mkStdGen 55))++bunnyTest :: IO ()+bunnyTest = let+  world mesh = lambertian (constantTexture (V3 0.3 0.3 1)) <$ triangleMesh mesh++  settings = defaultCameraSettings +    { cs_center = V3 0 0.5 2+    , cs_lookAt = V3 0 0 0+    , cs_imageWidth = 600+    , cs_samplesPerPixel = 100 +    , cs_background = grayFade+    }++  in do +    mesh <- fromRight undefined <$> readObj "images/bunny.obj"+    let center = fmap midpoint (boundingBox (triangleMesh mesh))+    let mesh' = transformVertices (rotateY (degrees 30) !*! scale 12 !*! translate (-center)) mesh+    writeImage "test_image.png" (raytrace settings (world mesh') (mkStdGen 55))++-- Version 0.1.0.0: ~1:45 (without redirection)+-- Version 0.2.0.0: ~1:00 with redirection, ~2:00 without cornellTest :: IO () cornellTest = cornellBox 200 50 --- This should take less than 70 seconds+-- Version 0.1.0.0: ~1:05 (without redirection)+-- Version 0.2.0.0: ~1:05 with redirection, ~1:10 without demoTest :: IO () demoTest = demo2 "test_image.png" 400 250 4  main :: IO ()-main = noiseTest+main = demoTest