diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,5 @@ +# Revision history for raytrace + +## 0.1.0.0 — 2025-07-18 + +* First version. Released on an unsuspecting world. diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -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) diff --git a/raytrace.cabal b/raytrace.cabal new file mode 100644 --- /dev/null +++ b/raytrace.cabal @@ -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 diff --git a/src/Graphics/Ray.hs b/src/Graphics/Ray.hs new file mode 100644 --- /dev/null +++ b/src/Graphics/Ray.hs @@ -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)) diff --git a/src/Graphics/Ray/Core.hs b/src/Graphics/Ray/Core.hs new file mode 100644 --- /dev/null +++ b/src/Graphics/Ray/Core.hs @@ -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. + } diff --git a/src/Graphics/Ray/Geometry.hs b/src/Graphics/Ray/Geometry.hs new file mode 100644 --- /dev/null +++ b/src/Graphics/Ray/Geometry.hs @@ -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 diff --git a/src/Graphics/Ray/Material.hs b/src/Graphics/Ray/Material.hs new file mode 100644 --- /dev/null +++ b/src/Graphics/Ray/Material.hs @@ -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)) diff --git a/src/Graphics/Ray/Noise.hs b/src/Graphics/Ray/Noise.hs new file mode 100644 --- /dev/null +++ b/src/Graphics/Ray/Noise.hs @@ -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 diff --git a/src/Graphics/Ray/Texture.hs b/src/Graphics/Ray/Texture.hs new file mode 100644 --- /dev/null +++ b/src/Graphics/Ray/Texture.hs @@ -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)) diff --git a/test/Main.hs b/test/Main.hs new file mode 100644 --- /dev/null +++ b/test/Main.hs @@ -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