packages feed

raytrace (empty) → 0.1.0.0

raw patch · 11 files changed

+1369/−0 lines, 11 filesdep +Colordep +basedep +linear

Dependencies added: Color, base, linear, massiv, massiv-io, mtl, random, raytrace

Files

+ CHANGELOG.md view
@@ -0,0 +1,5 @@+# Revision history for raytrace++## 0.1.0.0 — 2025-07-18++* First version. Released on an unsuspecting world.
+ LICENSE view
@@ -0,0 +1,29 @@+Copyright (c) 2025, Owen Bechtel+++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++    * Redistributions of source code must retain the above copyright+      notice, this list of conditions and the following disclaimer.++    * Redistributions in binary form must reproduce the above+      copyright notice, this list of conditions and the following+      disclaimer in the documentation and/or other materials provided+      with the distribution.++    * Neither the name of the copyright holder nor the names of its+      contributors may be used to endorse or promote products derived+      from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ README.md view
@@ -0,0 +1,63 @@+# 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.++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+* Texture mapping+* Perlin noise for procedurally generated textures+* Parallel computation of pixels+* Bounding volume hierarchies+* Affine transformations (rotation, translation, and so on)+* Optional defocusing to simulate a real camera lens++Possible future additions:+* Motion blur+* Triangular meshes+* Less noise in scenes with small light sources++![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:++![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.)++## Example Usage++```haskell+module Main where++import Graphics.Ray+import Linear (V3(V3))+import System.Random (newStdGen)+import Data.Functor.Identity (Identity)++world :: Geometry Identity Material+world = group+  [ lambertian (checkerTexture 20 10 0.2 0.8) <$ sphere (V3 0 0 0) 1+  , lambertian (constantTexture (V3 0 0.2 0.5)) <$ sphere (V3 0 (-1000) 0) 999+  , mirror (constantTexture 0.8) <$ parallelogram (V3 (-3.25) (-1) (-0.75)) (V3 1.25 0 (-1.25)) (V3 0 2 0)+  ]++settings :: CameraSettings+settings = defaultCameraSettings+  { cs_center = V3 (-0.75) 0 2+  , cs_lookAt = V3 0 0 (-1)+  , cs_aspectRatio = 16 / 9+  , cs_imageWidth = 600+  , cs_samplesPerPixel = 50+  }++main :: IO ()+main = do+  seed <- newStdGen+  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)
+ raytrace.cabal view
@@ -0,0 +1,62 @@+cabal-version:      3.0+name:               raytrace+version:            0.1.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.+    +    See README.md for more information.++license:            BSD-3-Clause+license-file:       LICENSE++author:             Owen Bechtel+maintainer:         ombspring@gmail.com++build-type:         Simple++extra-doc-files:    CHANGELOG.md, README.md++library+    hs-source-dirs:   src+    default-language: Haskell2010+    ghc-options:      -Wall+    +    exposed-modules:+        Graphics.Ray, +        Graphics.Ray.Core, +        Graphics.Ray.Geometry, +        Graphics.Ray.Texture, +        Graphics.Ray.Material, +        Graphics.Ray.Noise++    build-depends:    +        base >= 4.17.2.1 && < 5, +        linear >= 1.23.2 && < 1.24,+        massiv >= 1.0.5 && < 1.1,+        massiv-io >= 1.0.0 && < 1.1,+        Color >= 0.4.0 && < 0.5, +        random >= 1.3.0 && < 1.4, +        mtl >= 2.2.2 && < 2.3++test-suite test+    type:             exitcode-stdio-1.0+    hs-source-dirs:   test+    default-language: Haskell2010+    ghc-options:      -Wall -threaded -rtsopts -with-rtsopts=-N++    main-is:          Main.hs++    build-depends:    +        base >= 4.17.2.1 && < 5, +        linear >= 1.23.2 && < 1.24,+        massiv >= 1.0.5 && < 1.1,+        massiv-io >= 1.0.0 && < 1.1,+        Color >= 0.4.0 && < 0.5, +        random >= 1.3.0 && < 1.4, +        mtl >= 2.2.2 && < 2.3, +        raytrace
+ src/Graphics/Ray.hs view
@@ -0,0 +1,176 @@+{-# LANGUAGE LambdaCase #-}+{-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE InstanceSigs #-}+{-# LANGUAGE FlexibleInstances #-}+module Graphics.Ray +  ( -- * Camera+    CameraSettings(..), defaultCameraSettings+    -- * Ray Tracing+  , ToRandom(toRandom), raytrace+    -- * Image IO+  , readImage, writeImage, writeImageSqrt+    -- * Re-exports+  , module Graphics.Ray.Core+  , module Graphics.Ray.Geometry+  , module Graphics.Ray.Material+  , module Graphics.Ray.Texture+  , module Graphics.Ray.Noise+  ) where++import Graphics.Ray.Core+import Graphics.Ray.Geometry+import Graphics.Ray.Material+import Graphics.Ray.Texture+import Graphics.Ray.Noise++import Linear (V2(V2), V3(V3), (*^), (^*), normalize, cross, (^/), zero)+import System.Random (StdGen, random, splitGen)+import Data.Massiv.Array (B, D, S, U, Ix2((:.)), (!))+import qualified Data.Massiv.Array as A+import qualified Data.Massiv.Array.IO as I+import Graphics.Pixel.ColorSpace (SRGB, Linearity(Linear, NonLinear))+import qualified Graphics.Pixel.ColorSpace as C+import Control.Monad.State (State, state, evalState)+import Control.Monad (replicateM)+import Data.Functor.Identity (Identity, runIdentity)++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) +  }++-- | By default, the camera is positioned at the origin looking in the negative z direction, with the positive y direction being upward+-- (and the positive x direction being rightward). The remaining attributes are as follows:+-- +-- @+-- cs_vfov = pi / 2+-- cs_aspectRatio = 1.0+-- cs_imageWidth = 100+-- cs_samplesPerPixel = 10+-- cs_maxRecursionDepth = 10+-- cs_background = const (V3 1 1 1)+-- cs_defocusAngle = 0.0+-- cs_focusDist = 10.0+-- @+defaultCameraSettings :: CameraSettings+defaultCameraSettings = CameraSettings+  { cs_center = V3 0 0 0+  , cs_lookAt = V3 0 0 (-1)+  , cs_up = V3 0 1 0+  , cs_vfov = pi / 2+  , cs_aspectRatio = 1.0+  , cs_imageWidth = 100+  , cs_samplesPerPixel = 10+  , cs_maxRecursionDepth = 10+  , cs_background = const (V3 1 1 1)+  , cs_defocusAngle = 0.0+  , cs_focusDist = 10.0+  }++class ToRandom m where+  toRandom :: m a -> State StdGen a++instance ToRandom Identity where+  toRandom :: Identity a -> State StdGen a+  toRandom = pure . runIdentity++instance ToRandom (State StdGen) where+  toRandom :: State StdGen a -> State StdGen a+  toRandom = id++-- | 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+  imageHeight = round (fromIntegral cs_imageWidth / cs_aspectRatio)+  viewportHeight = cs_focusDist * tan (cs_vfov / 2) * 2+  viewportWidth = viewportHeight * fromIntegral cs_imageWidth / fromIntegral imageHeight++  w = normalize (cs_center - cs_lookAt)+  u = normalize (cross cs_up w) +  v = cross w u++  across = viewportWidth *^ u+  down = -(viewportHeight *^ v)++  topLeft = cs_center - w ^* cs_focusDist - across ^/ 2 - down ^/ 2+  pixelU = across ^/ fromIntegral cs_imageWidth+  pixelV = down ^/ fromIntegral imageHeight++  defocusRadius = cs_focusDist * tan (cs_defocusAngle / 2)+  defocusDiskU = u ^* defocusRadius+  defocusDiskV = v ^* defocusRadius+  +  sampleDefocusDisk :: State StdGen Point3+  sampleDefocusDisk = do+    V2 x y <- randomInUnitDisk+    pure (cs_center + x *^ defocusDiskU + y *^ defocusDiskV)++  samplePixel :: Int -> Int -> State StdGen Point3+  samplePixel i j = do+    x <- state random+    y <- state random+    pure (topLeft + (fromIntegral i + x) *^ pixelU + (fromIntegral j + y) *^ pixelV)++  getRay :: Int -> Int -> State StdGen Ray+  getRay i j = do+    origin <- sampleDefocusDisk+    target <- samplePixel i j+    pure (Ray origin (target - origin))++  rayColor :: Int -> Ray -> State StdGen Color+  rayColor depth ray+    | depth <= 0 = pure zero+    | otherwise =+    toRandom (hitWorld 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)+  +  pixelColor :: Int -> Int -> State StdGen Color+  pixelColor i j = do+    colors <- replicateM cs_samplesPerPixel (getRay i j >>= rayColor cs_maxRecursionDepth)+    pure (sum colors ^/ fromIntegral cs_samplesPerPixel)++  -- array of random seeds for each pixel (constructed using splitGen)+  seeds :: A.Matrix B StdGen+  (_, seeds) = A.randomArrayS seed (A.Sz (imageHeight :. cs_imageWidth)) splitGen++  in A.makeArray A.Par (A.Sz (imageHeight :. cs_imageWidth)) (\ix@(j :. i) -> evalState (pixelColor i j) (seeds ! ix))++-- | Read an image file, converting each pixel to linear RGB color space.+readImage :: FilePath -> IO (A.Matrix U Color)+readImage path = A.compute . A.map fromPixel <$> (I.readImageAuto path :: IO (A.Matrix S (C.Pixel (SRGB 'Linear) Double)))+  where+    fromPixel :: C.Pixel (SRGB 'Linear) Double -> Color+    fromPixel (C.Pixel (C.ColorSRGB r g b)) = V3 r g b++-- | Write an array of linear RGB colors to an image file.+writeImage :: (A.Source r Color) => FilePath -> A.Matrix r Color -> IO ()+writeImage path m = I.writeImageAuto path (A.map toPixel m)+  where+    toPixel :: Color -> C.Pixel (SRGB 'Linear) Double+    toPixel (V3 r g b) = C.Pixel (C.ColorSRGB r g b)++-- | Write an array to an image file, using a slightly incorrect color space conversion function.+-- This function exists for testing purposes.+writeImageSqrt :: (A.Source r Color) => FilePath -> A.Matrix r Color -> IO ()+writeImageSqrt path m = I.writeImageAuto path (A.map toPixel m)+  where+    toPixel :: Color -> C.Pixel (SRGB 'NonLinear) Double+    toPixel (V3 r g b) = C.Pixel (C.ColorSRGB (sqrt r) (sqrt g) (sqrt b))
+ src/Graphics/Ray/Core.hs view
@@ -0,0 +1,155 @@+module Graphics.Ray.Core +  ( -- * Numbers+    infinity, degrees+    -- * Vectors and Rays+  , Vec3, Point3, Color, Dim(..), component, reflect, randomUnitVector, randomInUnitDisk, Ray(Ray)+    -- * Intervals+  , Interval, inInterval, midpoint, padInterval+    -- * Boxes+  , Box, fromCorners, boxJoin, boxHull, allCorners, padBox, longestDim, overlapsBox+    -- * Hit Records+  , HitRecord(..)+  ) where++import Linear (V2, V3(V3), dot, quadrance, (*^), (^/))+import System.Random (RandomGen, randomR)+import Control.Monad.State (MonadState, state)+import Control.Applicative (liftA2)+import Data.Maybe (isJust)+import Data.Foldable (foldl')++-- | Floating-point infinity.+infinity :: Double+infinity = 1/0++-- | Convert an angle from degrees to radians. @degrees x@ means @x@ degrees.+degrees :: Double -> Double+degrees x = x * pi / 180++type Vec3 = V3 Double+type Point3 = V3 Double+type Color = V3 Double++data Dim = X | Y | Z+  deriving (Eq, Ord, Enum, Bounded, Show)++-- | Get the X, Y, or Z component of a vector.+component :: Dim -> V3 a -> a+component X (V3 x _ _) = x+component Y (V3 _ y _) = y+component Z (V3 _ _ z) = z++-- [private]+argMax :: Ord a => V3 a -> Dim+argMax (V3 x y z)+  | x > y     = if x > z then X else Z+  | 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++-- | Get a random 3-dimensional vector of norm 1.+randomUnitVector :: (RandomGen g, MonadState g m) => m Vec3+randomUnitVector = do+  vec <- state (randomR (-1, 1))+  let q = quadrance vec+  if 1e-8 <= q && q <= 1+    then pure (vec ^/ sqrt q)+    else randomUnitVector++-- | Get a random 2-dimensional vector of norm less than or equal to 1.+randomInUnitDisk :: (RandomGen g, MonadState g m) => m (V2 Double)+randomInUnitDisk = do+  vec <- state (randomR (-1, 1))+  if quadrance vec <= 1+    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.+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.+type Interval = (Double, Double)++-- [private]+size :: Interval -> Double+size (a, b) = b - a++-- | Test whether a number is in the interval.+inInterval :: Interval -> Double -> Bool+inInterval (tmin, tmax) t = tmin < t && t < tmax++-- | The midpoint of @(a, b)@ is @(a + b) / 2@.+midpoint :: Interval -> Double+midpoint (a, b) = (a + b) / 2++-- | Extend both the lower bound and the upper bound of the interval by the first argument.+padInterval :: Double -> Interval -> Interval+padInterval padding (tmin, tmax) = (tmin - padding, tmax + padding)++-- [private]+isect :: Interval -> Interval -> Maybe Interval+isect (a, b) (c, d) = let+  imin = max a c+  imax = min b d+  in if imin > imax then Nothing else Just (imin, imax)++-- [private]+overlapsInterval :: Interval -> Double -> Double -> Interval+overlapsInterval (tmin, tmax) x d = let+  t0 = (tmin - x) / d+  t1 = (tmax - x) / d+  in if t0 < t1 then (t0, t1) else (t1, t0)++-- | An axis-aligned box; the product of three intervals.+type Box = V3 Interval++-- | Create a box from two opposite corners.+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))++-- | Get a list of all eight corners of a box.+allCorners :: Box -> [ Point3 ]+allCorners (V3 i1 i2 i3) = +  [ V3 (f1 i1) (f2 i2) (f3 i3) +  | f1 <- [ fst, snd ], f2 <- [ fst, snd ], f3 <- [ fst, snd ]+  ]++-- | Extend the box by the first argument in all six directions.+padBox :: Double -> Box -> Box+padBox padding = fmap (padInterval padding)++-- | The dimension in which the box is the longest.+longestDim :: Box -> Dim+longestDim = argMax . fmap size++-- | Test whether any part of the ray, when restricted to the interval, is within the box.+overlapsBox :: Box -> Ray -> Interval -> Bool+overlapsBox (V3 ix iy iz) (Ray (V3 ox oy oz) (V3 dx dy dz)) (tmin, tmax) =+  isJust $ do+    (tmin', tmax') <- isect (tmin, tmax) (overlapsInterval ix ox dx)+    (tmin'', tmax'') <- isect (tmin', tmax') (overlapsInterval iy oy dy)+    isect (tmin'', tmax'') (overlapsInterval iz oz dz)++-- | A record of an intersection of a ray with a surface.+data HitRecord = HitRecord+  { hr_t :: Double -- ^ The @t@ parameter of the intersection.+  , hr_point :: Point3 -- ^ The location of the intersection.+  , hr_normal :: Vec3 -- ^ The normal vector of the surface in the direction that the ray came from. Should always be a unit vector.+  , 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.+  }
+ src/Graphics/Ray/Geometry.hs view
@@ -0,0 +1,295 @@+{-# LANGUAGE InstanceSigs #-}+{-# LANGUAGE LambdaCase #-}+{-# LANGUAGE RecordWildCards #-}+module Graphics.Ray.Geometry +  ( -- * Geometry+    Geometry(Geometry), pureGeometry, boundingBox+    -- * Surfaces and Volumes+  , sphere, parallelogram, cuboid, constantMedium+    -- * Groups+  , group, bvhNode, Tree(Leaf, Node), bvhTree, autoTree+    -- * Transformations+  , transform, translate, rotateX, rotateY, rotateZ+  ) 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 System.Random (StdGen, random)+import Control.Monad.State (State, state)+import Control.Monad (guard, foldM)+import Control.Applicative ((<|>))+import Data.List (sortOn)+import Data.Bifunctor (first, second)+import Data.Functor.Identity (Identity(Identity), runIdentity)+import Data.Functor ((<&>))++-- | 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@+-- (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)))++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))++-- | 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)))++-- | Get a geometry's bounding box.+boundingBox :: Geometry m a -> Box+boundingBox (Geometry b _) = b++-- | Construct a sphere with the given center and radius.+sphere :: Point3 -> Double -> Geometry Identity ()+sphere center radius = let+  diag = V3 radius radius radius+  bbox = fromCorners (center - diag) (center + diag)++  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+    guard (discriminant >= 0)+    +    let sqrtd = sqrt discriminant+    let root1 = (h - sqrtd) / a+    let root2 = (h + sqrtd) / a++    t <- +      if inInterval bounds root1 +        then Just root1 +      else if inInterval bounds root2+        then Just root2+      else Nothing+    +    let point = orig + t *^ dir+    let outwardNormal = (point - center) ^/ radius+    let frontSide = dot dir outwardNormal <= 0+    let hit = HitRecord+          { hr_t = t+          , hr_point = point+          , hr_normal = if frontSide then outwardNormal else -outwardNormal+          , hr_frontSide = frontSide+          , hr_uv = sphereUV outwardNormal -- only computed when necessary thanks to laziness+          }+    Just (hit, ())+  +  in Geometry bbox hitSphere++-- [private]+-- With default camera settings (-z direction is forward, +y direction is up),+-- texture images will be wrapped around the sphere starting and ending on the+-- far side of the sphere.+sphereUV :: Vec3 -> V2 Double+sphereUV (V3 x y z) = V2 u v+  where+    u = atan2 x z / (2 * pi) + 0.5+    v = acos (-y) / pi ++-- | Construct a parallelogram from a corner point and two edge vectors.+-- Which side is the \"front side\" is determined by the right hand rule.+parallelogram :: Point3 -> Vec3 -> Vec3 -> Geometry Identity ()+parallelogram q u v = let+  cp = cross u v+  area = norm cp+  normal = cp ^/ area+  normalS = normal ^/ area+  n_dot_q = dot normal q++  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+    let denom = dot normal dir+    guard (abs denom > 1e-8)+    let t = (n_dot_q - dot normal 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)+    let frontSide = denom < 0++    let hit = HitRecord+          { hr_t = t+          , hr_point = p+          , hr_normal = if frontSide then normal else -normal+          , hr_frontSide = frontSide+          , hr_uv = V2 a b+          }+    Just (hit, ())++  in Geometry bbox hitParallelogram ++-- | 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+  dx = V3 (xmax - xmin) 0 0+  dy = V3 0 (ymax - ymin) 0+  dz = V3 0 0 (zmax - zmin)+  in group +    [ parallelogram (V3 xmin ymin zmax) dx dy -- front+    , parallelogram (V3 xmax ymin zmin) (-dx) dy -- back+    , parallelogram (V3 xmin ymin zmin) dz dy -- left+    , parallelogram (V3 xmax ymin zmax) (-dz) dy -- right+    , parallelogram (V3 xmin ymax zmax) dx (-dz) -- top+    , 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'.+constantMedium+  :: Double -- ^ Density +  -> Geometry Identity () -- ^ Surface (assumed to be convex in current implementation)+  -> 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+      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+            let hitDist = negInvDensity * log rand+            guard (hitDist < inDist)+            let t = t1 + hitDist / rayScale+            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+                  }+            Just (hit, ())++  in Geometry bbox hitMedium++-- | Group multiple geometric objects into a single object. When testing if a ray hits a group, +-- every constituent of the group is tested without regard to its position. With a large number of objects,+-- use 'bvhTree' for greater efficiency.+{-# 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) =+    let try (tmax', knownHit) (Geometry _ hitObj) =+          hitObj ray (tmin, tmax') <&> \case+            Nothing -> (tmax', knownHit)+            Just (hit, mat) -> (hr_t hit, Just (hit, mat))+    in snd <$> foldM try (tmax, Nothing) obs+  +  in Geometry bbox hitGroup++-- | A single node in a bounding volume hierarchy. Before testing whether a ray hits each child,+-- it tests whether the ray hits a bounding box containing the two children.+{-# 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++  hitBvhNode 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))+    | 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+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+  obs -> let+    d = longestDim (boxHull (map boundingBox obs))+    obs' = sortOn (midpoint . component d . boundingBox) obs+    (left, right) = splitAt (length obs `div` 2) obs'+    in Node (autoTree left) (autoTree right)++-- | Apply an affine transformation (represented as a 4 by 4 matrix whose bottom row is 0 0 0 1) to a geometric object.+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 ->+    let ray' = Ray (inv_m !* V4.point orig) (inv_m !* V4.vector dir) in+    flip (fmap . fmap . first) (hitObj ray' ival) $ \hit@(HitRecord {..}) ->+      hit { hr_point = m34 !* V4.point hr_point, hr_normal = m34 !* V4.vector hr_normal }++-- | Translation.+translate :: Vec3 -> M44 Double+translate (V3 x y z) = V4+  (V4 1 0 0 x)+  (V4 0 1 0 y)+  (V4 0 0 1 z)+  (V4 0 0 0 1)++-- | Rotation about the X axis.+rotateX :: Double -> M44 Double+rotateX angle = V4+  (V4 1 0 0 0)+  (V4 0 c (-s) 0)+  (V4 0 s c 0)+  (V4 0 0 0 1)+  where+    c = cos angle+    s = sin angle++-- | Rotation about the Y axis.+rotateY :: Double -> M44 Double+rotateY angle = V4+  (V4 c 0 s 0)+  (V4 0 1 0 0)+  (V4 (-s) 0 c 0)+  (V4 0 0 0 1)+  where +    c = cos angle+    s = sin angle++-- | Rotation about the Z axis.+rotateZ :: Double -> M44 Double+rotateZ angle = V4+  (V4 c (-s) 0 0)+  (V4 s c 0 0)+  (V4 0 0 1 0)+  (V4 0 0 0 1)+  where+    c = cos angle+    s = sin angle++-- [private]+dropLast :: V4 a -> V3 a+dropLast (V4 x y z _) = V3 x y z
+ src/Graphics/Ray/Material.hs view
@@ -0,0 +1,94 @@+{-# LANGUAGE RecordWildCards #-}+module Graphics.Ray.Material +  ( Material(Material)+  , lightSource, pitchBlack, lambertian, mirror, metal, dielectric, transparent, isotropic+  ) where++import Graphics.Ray.Core+import Graphics.Ray.Texture++import Linear (V3(V3), zero, normalize, dot, quadrance, (*^))+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 that emits light and does not reflect rays.+lightSource :: Texture -> Material+lightSource (Texture tex) = Material $+  \_ (HitRecord {..}) -> pure (tex hr_point hr_uv, Nothing)++-- | A material that absorbs all light. Identical in principle to @'lambertian' ('constantTexture' 0)@, @'isotropic' ('constantTexture' 0)@,+-- and so on, but avoids unecessary computation.+pitchBlack :: Material+pitchBlack = Material $ \ _ _ -> pure (zero, Nothing)++-- | 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)))++-- | 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)))++-- | 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+    u <- randomUnitVector+    let dir' = normalize (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)++-- [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++-- | 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+    let iorRatio = if hr_frontSide then 1/ior else ior+    let u = normalize dir+    let cosTheta = min 1 (dot hr_normal (-u))+    let sinTheta = sqrt (1 - cosTheta * cosTheta)+    let cannotRefract = iorRatio * sinTheta > 1++    let r0 = (1 - iorRatio) / (1 + iorRatio)+    let r0' = r0 * r0+    let reflectance = r0' + (1 - r0') * (1 - cosTheta)**5 -- Schlick approximation+    x <- state random++    let dir' = if cannotRefract || x < reflectance+          then reflect hr_normal u+          else refract iorRatio cosTheta hr_normal u+    +    pure (zero, Just (V3 1 1 1, Ray hr_point dir'))    ++-- | A material that lets all 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))++-- | 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))
+ src/Graphics/Ray/Noise.hs view
@@ -0,0 +1,92 @@+module Graphics.Ray.Noise+  ( perlinNoise, fractalNoise, turbulence+  ) where++import Graphics.Ray.Core++import Linear (V3(V3), dot)+import System.Random (StdGen, mkStdGen)+import qualified Data.Massiv.Array as A+import Data.Massiv.Array (U, (!))+import Control.Monad.State (State, evalState)+import Control.Monad (replicateM)+import Data.Bits ((.&.), xor)++smoothstep :: Double -> Double+smoothstep x = x*x * (3 - 2*x)++-- | 3-dimensional Perlin noise.+-- Repeats with period 256 in all three dimensions.+-- Return values are always within the range \([-\sqrt{3}/2, +\sqrt{3}/2]\).+perlinNoise :: V3 Double -> Double+perlinNoise (V3 x y z) = let+  ix = floor x+  iy = floor y+  iz = floor z+  fx = x - fromIntegral ix+  fy = y - fromIntegral iy+  fz = z - fromIntegral iz+  in sum $ do+    i <- [ 0, 1 ]+    j <- [ 0, 1 ]+    k <- [ 0, 1 ]+    let i' = fromIntegral i+    let j' = fromIntegral j +    let k' = fromIntegral k+    let grad = gradients ! ((permX ! ((ix + i) .&. 255)) `xor` (permY ! ((iy + j) .&. 255)) `xor` (permZ ! ((iz + k) .&. 255)))+    let rel = V3 (fx - i') (fy - j') (fz - k')+    let coef = smoothstep (i' * fx + (1 - i') * (1 - fx)) * smoothstep (j' * fy + (1 - j') * (1 - fy)) * smoothstep (k' * fz + (1 - k') * (1 - fz))+    [ coef * (grad `dot` rel) ]++-- | The first argument is the number of perlin noise terms. Each term has half the weight of the previous term, and double the noise frequency.+fractalNoise :: Int -> V3 Double -> Double+fractalNoise depth p = let+  coefs = iterate (/2) 1+  points = iterate (*2) p+  in sum (take depth (zipWith (*) coefs (map perlinNoise points)))++-- | The absolute value of 'fractalNoise'.+turbulence :: Int -> V3 Double -> Double+turbulence depth p = abs (fractalNoise depth p)++permX :: A.Vector U Int+permX = A.fromList A.Seq +  [ 179, 60, 35, 16, 220, 94, 67, 236, 106, 112, 65, 166, 83, 101, 246, 140, 219, 186, 113, 88, 153, 70, 34, 63, 157, 210, 212, 188, 54, 74, 23, 161+  , 28, 137, 126, 107, 183, 58, 134, 127, 211, 225, 17, 123, 150, 243, 160, 68, 75, 239, 173, 221, 89, 109, 61, 72, 159, 80, 154, 18, 214, 144, 197, 24+  , 105, 32, 84, 226, 136, 29, 139, 97, 230, 167, 165, 238, 27, 14, 50, 193, 46, 253, 240, 111, 69, 196, 130, 102, 104, 118, 204, 12, 169, 202, 142, 25+  , 245, 215, 149, 138, 185, 48, 223, 247, 47, 98, 143, 26, 87, 251, 103, 52, 234, 232, 218, 205, 92, 228, 162, 85, 122, 191, 242, 164, 129, 192, 255, 231+  , 147, 91, 178, 213, 176, 36, 120, 155, 241, 222, 177, 20, 152, 141, 51, 171, 250, 95, 71, 119, 254, 172, 53, 146, 135, 124, 125, 163, 235, 99, 7, 22+  , 100, 229, 93, 174, 3, 189, 116, 66, 217, 158, 237, 55, 151, 0, 148, 45, 86, 64, 216, 43, 252, 121, 200, 115, 39, 184, 82, 56, 9, 181, 62, 2+  , 81, 209, 44, 79, 19, 110, 41, 10, 194, 15, 132, 224, 249, 96, 233, 117, 49, 203, 5, 37, 11, 59, 168, 114, 90, 131, 31, 145, 40, 206, 13, 187+  , 133, 207, 4, 199, 170, 78, 30, 182, 248, 21, 6, 227, 57, 180, 73, 42, 128, 175, 108, 33, 244, 201, 198, 77, 195, 8, 38, 190, 76, 156, 208, 1+  ]++permY :: A.Vector U Int+permY = A.fromList A.Seq +  [ 252, 123, 131, 151, 243, 143, 12, 247, 196, 179, 99, 0, 178, 109, 71, 160, 93, 205, 127, 38, 142, 117, 152, 124, 166, 95, 200, 121, 15, 17, 10, 190+  , 116, 158, 173, 75, 248, 191, 197, 58, 70, 184, 226, 146, 6, 239, 165, 113, 218, 34, 83, 77, 74, 5, 176, 85, 112, 147, 59, 66, 14, 31, 2, 21+  , 198, 108, 255, 11, 36, 156, 96, 73, 1, 189, 126, 50, 220, 16, 249, 23, 139, 135, 141, 92, 159, 4, 119, 174, 171, 253, 86, 227, 251, 172, 7, 149+  , 207, 212, 224, 44, 187, 91, 175, 84, 28, 211, 180, 195, 52, 98, 244, 125, 138, 54, 210, 201, 209, 219, 133, 240, 9, 202, 199, 42, 51, 48, 104, 154+  , 88, 53, 64, 105, 242, 63, 223, 222, 67, 238, 32, 134, 72, 62, 101, 150, 94, 161, 19, 236, 215, 97, 206, 22, 188, 230, 170, 18, 13, 162, 129, 90+  , 246, 130, 35, 46, 43, 213, 29, 76, 241, 61, 30, 136, 235, 87, 114, 68, 183, 177, 250, 132, 3, 122, 110, 145, 81, 78, 232, 69, 60, 80, 216, 65+  , 164, 237, 157, 203, 25, 221, 254, 181, 8, 182, 214, 24, 40, 102, 37, 106, 228, 39, 229, 163, 111, 186, 245, 55, 41, 225, 234, 49, 45, 231, 107, 103+  , 168, 153, 47, 56, 118, 120, 137, 155, 148, 82, 33, 204, 79, 169, 144, 20, 27, 128, 192, 217, 100, 208, 115, 140, 233, 89, 193, 185, 167, 57, 26, 194+  ]++permZ :: A.Vector U Int+permZ = A.fromList A.Seq +  [153, 90, 163, 138, 20, 136, 79, 100, 93, 38, 185, 31, 193, 43, 161, 2, 30, 37, 87, 6, 127, 207, 96, 51, 27, 227, 203, 215, 155, 190, 106, 94+  , 65, 183, 114, 71, 74, 219, 245, 39, 140, 216, 195, 191, 3, 214, 13, 23, 168, 179, 218, 133, 102, 53, 194, 124, 166, 108, 116, 246, 35, 109, 220, 121+  , 205, 110, 83, 242, 252, 231, 25, 128, 61, 57, 187, 222, 189, 228, 148, 101, 239, 162, 48, 150, 55, 174, 178, 42, 200, 160, 58, 206, 11, 92, 204, 67+  , 8, 113, 34, 172, 181, 177, 66, 16, 243, 159, 197, 135, 249, 241, 28, 147, 158, 139, 176, 88, 29, 149, 254, 226, 192, 201, 202, 186, 69, 45, 72, 199+  , 107, 99, 209, 10, 230, 217, 247, 89, 50, 129, 15, 64, 248, 167, 180, 146, 210, 169, 4, 81, 97, 238, 5, 237, 125, 236, 95, 9, 104, 221, 32, 86+  , 165, 134, 224, 182, 198, 234, 253, 59, 164, 52, 60, 208, 103, 22, 46, 188, 54, 12, 1, 112, 144, 137, 47, 156, 98, 78, 18, 82, 175, 225, 19, 184+  , 244, 49, 120, 157, 130, 70, 80, 251, 212, 250, 119, 118, 196, 223, 105, 68, 84, 21, 145, 173, 56, 235, 91, 17, 126, 76, 131, 117, 152, 111, 33, 36+  , 233, 141, 122, 85, 7, 123, 44, 62, 115, 26, 0, 40, 229, 211, 213, 63, 143, 77, 14, 24, 142, 154, 73, 132, 171, 41, 170, 240, 75, 151, 232, 255+  ]++gradients :: A.Vector U (V3 Double)+gradients = let+  action = replicateM 256 randomUnitVector :: State StdGen [V3 Double] +  vecs = evalState action (mkStdGen 666)+  in A.fromList A.Seq vecs
+ src/Graphics/Ray/Texture.hs view
@@ -0,0 +1,79 @@+{-# LANGUAGE FlexibleContexts #-}+module Graphics.Ray.Texture +  ( Texture(Texture)+  , constantTexture, solidTexture, uvTexture, imageTexture, checkerTexture, noiseTexture, marbleTexture+  ) where++import Graphics.Ray.Core+import Graphics.Ray.Noise++import Linear (V2(V2), V3, (^*), (*^), dot)+import qualified Data.Massiv.Array as A+import Data.Massiv.Array (Ix2((:.)), (!))+import Data.Bits ((.&.))++newtype Texture = Texture (Point3 -> V2 Double -> Color)++-- | Texture that is the same color everywhere.+constantTexture :: Color -> Texture+constantTexture color = Texture (\_ _ -> color)++-- | Texture assigning a color to each point in 3D space.+solidTexture :: (Point3 -> Color) -> Texture+solidTexture f = Texture (\p _ -> f p)++-- | Texture whose color depends on two parameters (u, v) that vary across a surface.+uvTexture :: (V2 Double -> Color) -> Texture+uvTexture f = Texture (const f)++-- | Convert an image into a texture. (u, v) = (0, 0) maps to the bottom left of the image,+-- (u, v) = (1, 1) maps to the top right, and coordinates outside of this range wrap around.+imageTexture :: (A.Manifest r Color) => A.Matrix r Color -> Texture+imageTexture image = let +  A.Sz (h :. w) = A.size image +  w' = fromIntegral w+  h' = fromIntegral h+  in+  uvTexture $ \(V2 u v) -> let+    i = floor (u * w') `mod` w+    j = floor ((1 - v) * h') `mod` h+    in+    image ! (j :. i)++-- | UV texture with two colors alternating in a checkerboard pattern. The first two arguments+-- are the dimensions of the checkerboard.+checkerTexture :: Int -> Int -> Color -> Color -> Texture+checkerTexture n_u n_v c0 c1 = let+  n_u' = fromIntegral n_u+  n_v' = fromIntegral n_v+  in+  uvTexture $ \(V2 u v) -> let+    i = floor (u * n_u')+    j = floor (v * n_v')+    in if (i + j) .&. 1 == (0 :: Int) then c0 else c1++-- | Perlin noise texture.+noiseTexture +  :: Int -- ^ Number of layers of noise (see 'fractalNoise')+  -> Double -- ^ Noise frequency+  -> V3 Double -- ^ Shift applied before calling noise function+  -> Color -- ^ Color 1+  -> Color -- ^ Color 2+  -> Texture+noiseTexture k freq shift color0 color1 = let+  scale = 0.5 / 0.8+  getNoise p = fractalNoise k (p ^* freq + shift) * scale + 0.5 +  diff = color1 - color0+  in solidTexture $ \p -> color0 + diff ^* getNoise p++-- | Texture with noisy black and white stripes, resulting in a marble-like appearance. +marbleTexture +  :: Vec3 -- ^ Direction of stripes+  -> Double -- ^ Frequency+  -> V3 Double -- ^ Shift applied before calling noise function+  -> Texture +marbleTexture dir freq shift =+  solidTexture $ \p -> let+    sinArg = freq * dot dir p+    noise = 10 * turbulence 7 (0.25 * freq *^ p + shift)+    in 1 ^* (0.5 + 0.5 * sin (sinArg + noise))
+ test/Main.hs view
@@ -0,0 +1,319 @@+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE PartialTypeSignatures #-}+{-# LANGUAGE ViewPatterns #-}+module Main where++import Graphics.Ray++import Linear (V3(V3), (*^), normalize, 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)++sky :: Ray -> Color+sky (Ray _ (normalize -> V3 _ y _)) = +  let a = 0.5 * (y + 1) in+  (1 - a) *^ V3 1 1 1 + a *^ V3 0.5 0.7 1++metalTest :: IO ()+metalTest = let+  materialGround = lambertian (constantTexture (V3 0.8 0.8 0.0))+  materialCenter = lambertian (constantTexture (V3 0.1 0.2 0.5))+  materialLeft = dielectric 1.5+  materialBubble = dielectric (1 / 1.5)+  materialRight = metal 1.0 (constantTexture (V3 0.8 0.6 0.2))++  world = group +    [ materialGround <$ sphere (V3 0 (-100.5) (-1)) 100+    , materialCenter <$ sphere (V3 0 0 (-1.2)) 0.5+    , materialLeft <$ sphere (V3 (-1) 0 (-1)) 0.5 +    , materialBubble <$ sphere (V3 (-1) 0 (-1)) 0.4+    , materialRight <$ sphere (V3 1 0 (-1)) 0.5+    ]++  settings = defaultCameraSettings +    { cs_aspectRatio = 16/9+    , cs_imageWidth = 400+    , cs_samplesPerPixel = 100+    , cs_maxRecursionDepth = 50+    , cs_background = sky+    , cs_center = V3 (-2) 2 1+    , cs_lookAt = V3 0 0 (-1)+    , cs_vfov = degrees 20+    , cs_defocusAngle = degrees 10+    , cs_focusDist = 3.4+    }++  in do+  seed <- newStdGen+  writeImageSqrt "test_image.png" $ raytrace settings world seed++noiseTest :: IO ()+noiseTest = let+  groundMaterial = lambertian (noiseTexture 2 2.0 (V3 10 0 0) 0 1)+  ballMaterial = lambertian (marbleTexture (V3 0 0 1) 4 0)++  world = group +    [ groundMaterial <$ sphere (V3 0 (-1000) 0) 1000+    , ballMaterial <$ sphere (V3 0 2 0) 2+    ]+  +  settings = defaultCameraSettings+    { cs_aspectRatio = 16 / 9+    , cs_imageWidth = 400+    , cs_samplesPerPixel = 100+    , cs_maxRecursionDepth = 50+    , cs_background = sky+    , cs_vfov = degrees 20+    , cs_center = V3 13 2 3+    , cs_lookAt = V3 0 0 0+    }+    +  in do+  seed <- newStdGen+  writeImageSqrt "noise_test.png" $ raytrace settings world seed++quadTest :: IO ()+quadTest = let+  red = lambertian (constantTexture (V3 1.0 0.2 0.2))+  green = lambertian (constantTexture (V3 0.2 1.0 0.2))+  blue = lambertian (constantTexture (V3 0.2 0.2 1.0))+  orange = lambertian (constantTexture (V3 1.0 0.5 0.0))+  teal = lambertian (constantTexture (V3 0.2 0.8 0.8))++  world = group+    [ red <$ parallelogram (V3 (-3) (-2) 5) (V3 0 0 (-4)) (V3 0 4 0) +    , green <$ parallelogram (V3 (-2) (-2) 0) (V3 4 0 0) (V3 0 4 0)+    , blue <$ parallelogram (V3 3 (-2) 1) (V3 0 0 4) (V3 0 4 0)+    , orange <$ parallelogram (V3 (-2) 3 1) (V3 4 0 0) (V3 0 0 4)+    , teal <$ parallelogram (V3 (-2) (-3) 5) (V3 4 0 0) (V3 0 0 (-4))+    ]+  +  settings = defaultCameraSettings+    { cs_aspectRatio = 1+    , cs_imageWidth = 400+    , cs_samplesPerPixel = 100+    , cs_maxRecursionDepth = 50+    , cs_background = sky+    , cs_vfov = degrees 80+    , cs_center = V3 0 0 9+    , cs_lookAt = V3 0 0 0+    }++  in do+  seed <- newStdGen+  writeImageSqrt "test_image.png" $ raytrace settings world seed++cuboidTest :: IO ()+cuboidTest = do+  globe <- readImage "images/earthmap.jpg"+  let globeMaterial = lambertian (imageTexture globe)+  let object = globeMaterial <$ cuboid (fromCorners (-V3 1 2 0.5) (V3 1 2 0.5))+  let world = transform (translate (V3 0 0 (-3)) !*! rotateX (degrees 60)) object+  let settings = defaultCameraSettings { cs_imageWidth = 300 }+  writeImage "test_image.png" . raytrace settings world =<< newStdGen++sphereUVTest :: IO ()+sphereUVTest = do+  globe <- readImage "images/earthmap.jpg"+  let globeMaterial = lambertian (imageTexture globe)+  let world = globeMaterial <$ group [ sphere (V3 0 0 (-2)) 0.4, sphere (V3 0 0 (-1)) 0.4 ]+  let settings = defaultCameraSettings { cs_imageWidth = 1, cs_samplesPerPixel = 1, cs_vfov = 0.0001}+  writeImage "test_image.png" (raytrace settings world (mkStdGen 12))++demo1 :: IO ()+demo1 = let+  materialGround = lambertian (constantTexture (V3 0.5 0.5 0.5))+  materialGlass = dielectric 1.5+  materialDiffuse = lambertian (constantTexture (V3 0.4 0.2 0.1))+  materialMirror = mirror (constantTexture (V3 0.7 0.6 0.5))++  bigSpheres =+    [ materialGround <$ sphere (V3 0 (-1000) 0) 1000+    , materialGlass <$ sphere (V3 0 1 0) 1+    , materialDiffuse <$ sphere (V3 (-4) 1 0) 1+    , materialMirror <$ sphere (V3 4 1 0) 1+    ]++  genWorld :: State StdGen (Geometry Identity Material)+  genWorld = do+    fmap (bvhTree . autoTree . (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)++      if norm (center - V3 4 0.2 0) <= 0.9 then pure [] else do+        chooseMat <- state random+        mat <- +          if (chooseMat :: Double) < 0.8 then do+            color <- liftA2 (*) (state random) (state random)+            pure (lambertian (constantTexture color))+          else if chooseMat < 0.95 then do+            fuzz <- state (randomR (0, 0.5))+            color <- state (randomR (0.5, 1))+            pure (metal fuzz (constantTexture color))+          else pure materialGlass+        pure [ mat <$ sphere center 0.2 ]+  +  settings = defaultCameraSettings+    { cs_aspectRatio = 16 / 9+    , cs_imageWidth = 1200+    , cs_samplesPerPixel = 500+    , cs_maxRecursionDepth = 50+    , cs_vfov = degrees 20+    , cs_center = V3 13 2 3+    , cs_lookAt = V3 0 0 0+    , cs_defocusAngle = degrees 0.6+    , cs_focusDist = 10+    , cs_background = sky+    }++  in do+  seed <- newStdGen+  let (world, seed') = runState genWorld seed+  writeImageSqrt "test_image.png" $ raytrace settings world seed'++cornellBox :: Int -> Int -> IO ()+cornellBox samplesPerPixel maxRecurionDepth = let+  red = lambertian (constantTexture (V3 0.65 0.05 0.05))+  white = lambertian (constantTexture (V3 0.73 0.73 0.73))+  green = lambertian (constantTexture (V3 0.12 0.45 0.15))+  light = lightSource (constantTexture (V3 15 15 15))++  world = group+    [ green <$ parallelogram (V3 555 0 0) (V3 0 555 0) (V3 0 0 555)+    , red <$ parallelogram (V3 0 0 0) (V3 0 555 0) (V3 0 0 555)+    , light <$ parallelogram (V3 343 554 332) (V3 (-130) 0 0) (V3 0 0 (-105))+    , white <$ parallelogram (V3 0 0 0) (V3 555 0 0) (V3 0 0 555)+    , white <$ parallelogram (V3 555 555 555) (V3 (-555) 0 0) (V3 0 0 (-555))+    , white <$ parallelogram (V3 0 0 555) (V3 555 0 0) (V3 0 555 0)+    , transform (translate (V3 265 0 295) !*! rotateY (degrees 15)) $ white <$ cuboid (fromCorners (V3 0 0 0) (V3 165 330 165))+    , transform (translate (V3 130 0 65) !*! rotateY (degrees (-18))) $ white <$ cuboid (fromCorners (V3 0 0 0) (V3 165 165 165))+    ]++  settings = defaultCameraSettings+    { cs_aspectRatio = 1.0+    , cs_imageWidth = 600+    , cs_samplesPerPixel = samplesPerPixel+    , cs_maxRecursionDepth = maxRecurionDepth+    , cs_background = const (V3 0 0 0)+    , cs_vfov = degrees 40+    , cs_center = V3 278 278 (-800)+    , cs_lookAt = V3 278 278 0+    }++  in writeImageSqrt "cornell_box.png" . raytrace settings world =<< newStdGen++cornellSmoke :: IO ()+cornellSmoke = let+  red = lambertian (constantTexture (V3 0.65 0.05 0.05))+  white = lambertian (constantTexture (V3 0.73 0.73 0.73))+  green = lambertian (constantTexture (V3 0.12 0.45 0.15))+  light = lightSource (constantTexture (V3 7 7 7))++  surfaces = group+    [ green <$ parallelogram (V3 555 0 0) (V3 0 555 0) (V3 0 0 555)+    , red <$ parallelogram (V3 0 0 0) (V3 0 555 0) (V3 0 0 555)+    , light <$ parallelogram (V3 113 554 127) (V3 330 0 0) (V3 0 0 305)+    , white <$ parallelogram (V3 0 0 0) (V3 555 0 0) (V3 0 0 555)+    , white <$ parallelogram (V3 555 555 555) (V3 (-555) 0 0) (V3 0 0 (-555))+    , white <$ parallelogram (V3 0 0 555) (V3 555 0 0) (V3 0 555 0)+    ]+  +  cube1 = transform (translate (V3 265 0 295) !*! rotateY (degrees 15)) $ cuboid (fromCorners (V3 0 0 0) (V3 165 330 165))+  cube2 = transform (translate (V3 130 0 65) !*! rotateY (degrees (-18))) $ cuboid (fromCorners (V3 0 0 0) (V3 165 165 165))+  +  world = group+    [ pureGeometry surfaces+    , pitchBlack <$ constantMedium 0.01 cube1+    , isotropic (constantTexture 1) <$ constantMedium 0.01 cube2+    ]++  settings = defaultCameraSettings+    { cs_aspectRatio = 1.0+    , cs_imageWidth = 600+    , cs_samplesPerPixel = 200+    , cs_maxRecursionDepth = 50+    , cs_background = const (V3 0 0 0)+    , cs_vfov = degrees 40+    , cs_center = V3 278 278 (-800)+    , cs_lookAt = V3 278 278 0+    }++  in writeImageSqrt "cornell_smoke.png" . raytrace settings world =<< newStdGen++demo2 :: FilePath -> Int -> Int -> Int -> IO ()+demo2 path imageWidth samplesPerPixel maxRecursionDepth = let+  ground = lambertian (constantTexture (V3 0.48 0.83 0.53))+  white = lambertian (constantTexture (V3 0.73 0.73 0.73))++  generateBoxes :: State StdGen (Geometry Identity Material)+  generateBoxes = +    fmap ((ground <$) . bvhTree . autoTree) $ +    forM (liftA2 (,) [0..19] [0..19]) $ \(i, j) -> do+      let x0 = -1000 + i * 100+      let z0 = -1000 + j * 100+      let x1 = x0 + 100+      let z1 = z0 + 100+      let y0 = 0+      y1 <- state (randomR (1, 101))+      pure (cuboid (fromCorners (V3 x0 y0 z0) (V3 x1 y1 z1)))+  +  generateBalls :: State StdGen (Geometry Identity Material)+  generateBalls =+    fmap ((white <$) . transform (translate (V3 (-100) 270 395) !*! rotateY (degrees 15)) . bvhTree . autoTree) $+    replicateM 1000 $ do+      p <- state (randomR (0, 165)) +      pure (sphere p 10)+  +  boundary = sphere (V3 360 150 145) 70 +  +  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+    , 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+    , lambertian (imageTexture earth) <$ transform (translate (V3 400 0 400) !*! rotateY (pi/2)) (sphere (V3 0 200 0) 100)+    , lambertian (marbleTexture (V3 0 0 0.05) 4 0) <$ sphere (V3 220 280 300) 80+    ]+  +  generateWorld earth = do+    boxes <- generateBoxes+    balls <- generateBalls+    pure $ group +      [ pureGeometry (group (boxes : balls : largeObjects earth))+      , isotropic (constantTexture 1) <$ constantMedium 0.0001 (sphere (V3 0 0 0) 5000)+      , isotropic (constantTexture (V3 0.2 0.4 0.9)) <$ constantMedium 0.2 boundary+      ]+  +  settings = defaultCameraSettings+    { cs_center = V3 478 278 (-600)+    , cs_lookAt = V3 278 278 0+    , cs_vfov = degrees 40+    , cs_aspectRatio = 1.0+    , cs_imageWidth = imageWidth+    , cs_samplesPerPixel = samplesPerPixel+    , cs_maxRecursionDepth = maxRecursionDepth+    , cs_background = const 0+    }++  in do+    earth <- readImage "images/earthmap.jpg"+    seed <- newStdGen+    let (world, seed') = runState (generateWorld earth) seed+    writeImageSqrt path (raytrace settings world seed')++-- This should take less than 110 seconds+cornellTest :: IO ()+cornellTest = cornellBox 200 50++-- This should take less than 70 seconds+demoTest :: IO ()+demoTest = demo2 "test_image.png" 400 250 4++main :: IO ()+main = noiseTest