simplex noise

Sources/Voxelotl/CMakeLists.txt

@@ -23,7 +23,9 @@ add_executable(Voxelotl MACOSX_BUNDLE
   Random/Arc4Random.swift
   Random/PCG32Random.swift
   Random/Xoroshiro128.swift
+  Random/CoherentNoise.swift
   Random/PerlinNoiseGenerator.swift
+  Random/SimplexNoise.swift
   Random/RandomCollectionExtensions.swift
 
   # Resource classes

Sources/Voxelotl/Math/VectorExtensions.swift

@@ -1,7 +1,11 @@
 import simd
 
-extension SIMD3 {
+public extension SIMD2 where Scalar: Numeric & AdditiveArithmetic {
-  var xy: SIMD2<Scalar> {
+  @inline(__always) func dot(_ b: Self) -> Scalar { self.x * b.x + self.y * b.y }
+}
+
+public extension SIMD3 {
+  @inline(__always) var xy: SIMD2<Scalar> {
     get { .init(self.x, self.y) }
     set {
       self.x = newValue.x
@@ -9,7 +13,7 @@ extension SIMD3 {
     }
   }
 
-  var xz: SIMD2<Scalar> {
+  @inline(__always) var xz: SIMD2<Scalar> {
     get { .init(self.x, self.z) }
     set {
       self.x = newValue.x
@@ -18,7 +22,7 @@ extension SIMD3 {
   }
 }
 
-extension SIMD3 where Scalar: FloatingPoint {
+public extension SIMD3 where Scalar: FloatingPoint {
   @inline(__always) static var X: Self      { Self(1, 0, 0) }
   @inline(__always) static var Y: Self      { Self(0, 1, 0) }
   @inline(__always) static var Z: Self      { Self(0, 0, 1) }
@@ -31,6 +35,10 @@ extension SIMD3 where Scalar: FloatingPoint {
   @inline(__always) static var back: Self    {  Z }
 }
 
+public extension SIMD3 where Scalar: Numeric & AdditiveArithmetic {
+  @inline(__always) func dot(_ b: Self) -> Scalar { self.x * b.x + self.y * b.y + self.z * b.z }
+}
+
 extension SIMD3 where Scalar == Float {
   static func * (q: simd_quatf, v: Self) -> Self {
 #if true
@@ -51,8 +59,8 @@ extension SIMD3 where Scalar == Float {
   }
 }
 
-extension SIMD4 {
+public extension SIMD4 {
-  var xy: SIMD2<Scalar> {
+  @inline(__always) var xy: SIMD2<Scalar> {
     get { .init(self.x, self.y) }
     set {
       self.x = newValue.x
@@ -60,7 +68,7 @@ extension SIMD4 {
     }
   }
 
-  var xyz: SIMD3<Scalar> {
+  @inline(__always) var xyz: SIMD3<Scalar> {
     get { .init(self.x, self.y, self.z) }
     set {
       self.x = newValue.x
@@ -69,3 +77,7 @@ extension SIMD4 {
     }
   }
 }
+
+public extension SIMD4 where Scalar: Numeric & AdditiveArithmetic {
+  @inline(__always) func dot(_ b: Self) -> Scalar { self.x * b.x + self.y * b.y + self.z * b.z + self.w * b.w }
+}

Sources/Voxelotl/Random/CoherentNoise.swift

@@ -0,0 +1,17 @@
+public protocol CoherentNoise2D {
+  associatedtype Scalar: FloatingPoint & SIMDScalar
+
+  func get(_ point: SIMD2<Scalar>) -> Scalar
+}
+
+public protocol CoherentNoise3D {
+  associatedtype Scalar: FloatingPoint & SIMDScalar
+
+  func get(_ point: SIMD3<Scalar>) -> Scalar
+}
+
+public protocol CoherentNoise4D {
+  associatedtype Scalar: FloatingPoint & SIMDScalar
+
+  func get(_ point: SIMD4<Scalar>) -> Scalar
+}

Sources/Voxelotl/Random/PerlinNoiseGenerator.swift

@@ -1,10 +1,15 @@
 import Foundation
 
-public struct ImprovedPerlin<T: BinaryFloatingPoint & SIMDScalar> {
+public struct ImprovedPerlin<T: BinaryFloatingPoint & SIMDScalar>: CoherentNoise3D {
-  private let p: [Int32]
+  private let p: [Int16]
+
+  public init(permutation: [Int16]) {
+    assert(permutation.count == 0x100)
+    self.p = permutation
+  }
 
   public init(random: inout any RandomProvider) {
-    self.p = (0..<0x100).map { Int32($0) }.shuffled(using: &random)
+    self.p = (0..<0x100).map { Int16($0) }.shuffled(using: &random)
   }
 
   public func get(_ point: SIMD3<T>) -> T {

Sources/Voxelotl/Random/SimplexNoise.swift

@@ -0,0 +1,214 @@
+import Foundation
+
+public struct SimplexNoise<T: BinaryFloatingPoint & SIMDScalar>: CoherentNoise2D, CoherentNoise3D, CoherentNoise4D {
+  private let p: [Int16], pMod12: [Int16]
+
+  private let grad3: [SIMD3<T>] = [
+    .init(1, 1, 0), .init(-1,  1, 0), .init(1, -1,  0), .init(-1, -1,  0),
+    .init(1, 0, 1), .init(-1,  0, 1), .init(1,  0, -1), .init(-1,  0, -1),
+    .init(0, 1, 1), .init( 0, -1, 1), .init(0,  1, -1), .init( 0, -1, -1)
+  ]
+  private let grad4: [SIMD4<T>] = [
+    .init( 0,  1, 1, 1), .init( 0,  1,  1, -1), .init( 0,  1, -1, 1), .init( 0,  1, -1, -1),
+    .init( 0, -1, 1, 1), .init( 0, -1,  1, -1), .init( 0, -1, -1, 1), .init( 0, -1, -1, -1),
+    .init( 1,  0, 1, 1), .init( 1,  0,  1, -1), .init( 1,  0, -1, 1), .init( 1,  0, -1, -1),
+    .init(-1,  0, 1, 1), .init(-1,  0,  1, -1), .init(-1,  0, -1, 1), .init(-1,  0, -1, -1),
+    .init( 1,  1, 0, 1), .init( 1,  1,  0, -1), .init( 1, -1,  0, 1), .init( 1, -1,  0, -1),
+    .init(-1,  1, 0, 1), .init(-1,  1,  0, -1), .init(-1, -1,  0, 1), .init(-1, -1,  0, -1),
+    .init( 1,  1, 1, 0), .init( 1,  1, -1,  0), .init( 1, -1,  1, 0), .init( 1, -1, -1,  0),
+    .init(-1,  1, 1, 0), .init(-1,  1, -1,  0), .init(-1, -1,  1, 0), .init(-1, -1, -1,  0)
+  ]
+
+  public init(permutation: [Int16]) {
+    assert(permutation.count == 0x100)
+    self.p = permutation
+    self.pMod12 = self.p.map { $0 % 12 }
+  }
+
+  public init(random: inout any RandomProvider) {
+    self.p = (0..<0x100).map { Int16($0) }.shuffled(using: &random)
+    self.pMod12 = self.p.map { $0 % 12 }
+  }
+
+  public func get(_ point: SIMD2<T>) -> T {
+    // Skew space into rhobuses to find which simplex cell we're in
+    let f2 = 0.5 * (T(3).squareRoot() - 1)
+    let g2 = (3 - T(3).squareRoot()) / 6
+    let skewFactor = point.sum() * f2
+    let cellID = SIMD2(floor(point.x + skewFactor), floor(point.y + skewFactor))
+    let cellOrigin = cellID - (cellID.sum() * g2)
+    let corner0 = point - cellOrigin
+
+    // For the 2d case, the simplex shape is an equilateral triangle
+    // Determine which side of the rhombus on to find the simplex
+    let cornerOfs1: SIMD2<Int> = corner0.x > corner0.y ? .init(1, 0) : .init(0, 1)
+    let corner1 = corner0 - SIMD2<T>(cornerOfs1) + g2
+    let corner2 = corner0 - 1 + 2 * g2
+
+    // Compute the hashed gradient indices of the three simplex corners
+    let cellHash = SIMD2<Int>(cellID) & 0xFF
+    let gradIndex0 = permMod12(cellHash.x + perm(cellHash.y))
+    let gradIndex1 = permMod12(cellHash.x + cornerOfs1.x + perm(cellHash.y + cornerOfs1.y))
+    let gradIndex2 = permMod12(cellHash.x + 1 + perm(cellHash.y + 1))
+
+    // Calculate the contribution from the three corners
+    @inline(__always) func cornerContribution(_ corner: SIMD2<T>, _ gradID: Int) -> T {
+      var t = 0.5 - corner.x * corner.x - corner.y * corner.y
+      if t < 0 {
+        return 0
+      } else {
+        t *= t
+        return t * t * self.grad3[gradID].xy.dot(corner)
+      }
+    }
+    let noise0 = cornerContribution(corner0, gradIndex0)
+    let noise1 = cornerContribution(corner1, gradIndex1)
+    let noise2 = cornerContribution(corner2, gradIndex2)
+
+    return 70 * (noise0 + noise1 + noise2)
+  }
+
+  public func get(_ point: SIMD3<T>) -> T {
+    // Skew space into rhombohedrons to find which simplex cell we're in
+    let g3 = 1 / T(6), f3 = 1 / T(3)
+    let skewFactor = point.sum() * f3
+    let cellID = SIMD3(floor(point.x + skewFactor), floor(point.y + skewFactor), floor(point.z + skewFactor))
+    let cellOrigin = cellID - (cellID.sum() * g3)
+    let corner0 = point - cellOrigin
+
+    // For the 3D case, the simplex shape is a slightly irregular tetrahedron
+    // Compute the offsets for the second & third corners of the simplex
+    let (corner1ID, corner2ID): (SIMD3<Int>, SIMD3<Int>) = if corner0.x >= corner0.y {
+      if corner0.y >= corner0.z {
+        (.init(1, 0, 0), .init(1, 1, 0))  // X Y Z
+      } else if corner0.x >= corner0.z {
+        (.init(1, 0, 0), .init(1, 0, 1))  // X Z Y
+      } else {
+        (.init(0, 0, 1), .init(1, 0 ,1))  // Z X Y
+      }
+    } else {
+      if corner0.y < corner0.z {
+        (.init(0, 0, 1), .init(0, 1, 1))  // Z Y X
+      } else if corner0.x < corner0.z {
+        (.init(0, 1, 0), .init(0, 1, 1))  // Y Z X
+      } else {
+        (.init(0, 1, 0), .init(1, 1, 0))  // Y X Z
+      }
+    }
+    let corner1 = corner0 - SIMD3<T>(corner1ID) + g3
+    let corner2 = corner0 - SIMD3<T>(corner2ID) + 2 * g3
+    let corner3 = corner0 - 1 + 3 * g3
+
+    // Compute the hashed gradient indices of the four simplex corners
+    let cellHash = SIMD3<Int>(cellID) & 0xFF
+    let gradCorner0 = permMod12(
+      cellHash.x + perm(
+        cellHash.y + perm(
+          cellHash.z)))
+    let gradCorner1 = permMod12(
+      cellHash.x + corner1ID.x + perm(
+        cellHash.y + corner1ID.y + perm(
+          cellHash.z + corner1ID.z)))
+    let gradCorner2 = permMod12(
+      cellHash.x + corner2ID.x + perm(
+        cellHash.y + corner2ID.y + perm(
+          cellHash.z + corner2ID.z)))
+    let gradCorner3 = permMod12(
+      cellHash.x + 1 + perm(
+        cellHash.y + 1 + perm(
+          cellHash.z + 1)))
+
+    // Calculate the contribution from the four corners
+    @inline(__always) func cornerContribution(_ corner: SIMD3<T>, _ gradID: Int) -> T {
+      var t = 0.6 - corner.x * corner.x - corner.y * corner.y - corner.z * corner.z
+      if t < 0 {
+        return 0
+      } else {
+        t *= t
+        return t * t * self.grad3[gradID].dot(corner)
+      }
+    }
+    let noise0 = cornerContribution(corner0, gradCorner0)
+    let noise1 = cornerContribution(corner1, gradCorner1)
+    let noise2 = cornerContribution(corner2, gradCorner2)
+    let noise3 = cornerContribution(corner3, gradCorner3)
+
+    return 32 * (noise0 + noise1 + noise2 + noise3)
+  }
+
+  public func get(_ point: SIMD4<T>) -> T {
+    let g4 = (5 - T(5).squareRoot()) / 20, f4 = (T(5).squareRoot() - 1) / 4
+    let skewFactor = point.sum() * f4
+    let cellID = SIMD4(floor(point.x + skewFactor), floor(point.y + skewFactor), floor(point.z + skewFactor), floor(point.w + skewFactor))
+    let cellOrigin = cellID - (cellID.sum() * g4)
+    let corner0 = point - cellOrigin
+
+    // Determine which of the 24 simplices we're in
+    // Find the magnitude ordering
+    var rank = SIMD4<Int>.zero
+    if corner0.x > corner0.y { rank.x += 1 } else { rank.y += 1 }
+    if corner0.x > corner0.z { rank.x += 1 } else { rank.z += 1 }
+    if corner0.x > corner0.w { rank.x += 1 } else { rank.w += 1 }
+    if corner0.y > corner0.z { rank.y += 1 } else { rank.z += 1 }
+    if corner0.y > corner0.w { rank.y += 1 } else { rank.w += 1 }
+    if corner0.z > corner0.w { rank.z += 1 } else { rank.w += 1 }
+
+    // Compute 4D corners
+    let cornerOfs1 = SIMD4<Int>.zero.replacing(with: .one, where: rank .>= 3)
+    let cornerOfs2 = SIMD4<Int>.zero.replacing(with: .one, where: rank .>= 2)
+    let cornerOfs3 = SIMD4<Int>.zero.replacing(with: .one, where: rank .>= 1)
+    let corner1 = corner0 - SIMD4<T>(cornerOfs1) + g4
+    let corner2 = corner0 - SIMD4<T>(cornerOfs2) + 2 * g4
+    let corner3 = corner0 - SIMD4<T>(cornerOfs3) + 3 * g4
+    let corner4 = corner0 - 1 + 4 * g4
+
+    // Compute the hashed gradient indices of the five simplex corners
+    let cellHash = SIMD4<Int>(cellID) & 0xFF
+    let gradIndex0 = Int(perm(
+      cellHash.x + perm(
+        cellHash.y + perm(
+          cellHash.z + perm(
+            cellHash.w))))) & 0x1F
+    let gradIndex1 = Int(perm(
+      cellHash.x + cornerOfs1.x + perm(
+        cellHash.y + cornerOfs1.y + perm(
+          cellHash.z + cornerOfs1.z + perm(
+            cellHash.w + cornerOfs1.w))))) & 0x1F
+    let gradIndex2 = Int(perm(
+      cellHash.x + cornerOfs2.x + perm(
+        cellHash.y + cornerOfs2.y + perm(
+          cellHash.z + cornerOfs2.z + perm(
+            cellHash.w + cornerOfs2.w))))) & 0x1F
+    let gradIndex3 = Int(perm(
+      cellHash.x + cornerOfs3.x + perm(
+        cellHash.y + cornerOfs3.y + perm(
+          cellHash.z + cornerOfs3.z + perm(
+            cellHash.w + cornerOfs3.w))))) & 0x1F
+    let gradIndex4 = Int(perm(
+      cellHash.x + 1 + perm(
+        cellHash.y + 1 + perm(
+          cellHash.z + 1 + perm(
+            cellHash.w + 1))))) & 0x1F
+
+    // Calculate the contribution from the five corners
+    @inline(__always) func cornerContribution(_ corner: SIMD4<T>, _ gradID: Int) -> T {
+      var t = corner.indices.reduce(0.6) { accum, i in accum - corner[i] * corner[i] }
+      if t < 0 {
+        return 0
+      } else {
+        t *= t
+        return t * t * self.grad4[gradID].dot(corner)
+      }
+    }
+    let noise0 = cornerContribution(corner0, gradIndex0)
+    let noise1 = cornerContribution(corner1, gradIndex1)
+    let noise2 = cornerContribution(corner2, gradIndex2)
+    let noise3 = cornerContribution(corner3, gradIndex3)
+    let noise4 = cornerContribution(corner4, gradIndex4)
+
+    return 27 * (noise0 + noise1 + noise2 + noise3 + noise4)
+  }
+
+  @inline(__always) fileprivate func perm(_ idx: Int) -> Int { Int(self.p[idx & 0xFF]) }
+  @inline(__always) fileprivate func permMod12(_ idx: Int) -> Int { Int(self.pMod12[idx & 0xFF]) }
+}