Noise.hs

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