src/rand/noise.cpp

   1 #include "GaloisLFSR.hpp"
   2 #include "SimplexNoise.hpp"
   3 #include "WorleyNoise.hpp"
   4
   5 #include <cmath>
   6
   7
   8 namespace {
   9
  10 constexpr float one_third = 1.0f/3.0f;
  11 constexpr float one_sixth = 1.0f/6.0f;
  12
  13 }
  14
  15 namespace blank {
  16
  17 SimplexNoise::SimplexNoise(std::uint64_t seed) noexcept
  18 : grad({
  19         {  1.0f,  1.0f,  0.0f },
  20         { -1.0f,  1.0f,  0.0f },
  21         {  1.0f, -1.0f,  0.0f },
  22         { -1.0f, -1.0f,  0.0f },
  23         {  1.0f,  0.0f,  1.0f },
  24         { -1.0f,  0.0f,  1.0f },
  25         {  1.0f,  0.0f, -1.0f },
  26         { -1.0f,  0.0f, -1.0f },
  27         {  0.0f,  1.0f,  1.0f },
  28         {  0.0f, -1.0f,  1.0f },
  29         {  0.0f,  1.0f, -1.0f },
  30         {  0.0f, -1.0f, -1.0f },
  31 }) {
  32         GaloisLFSR random(seed ^ 0x0123456789ACBDEF);
  33         unsigned char value;
  34         for (size_t i = 0; i < 256; ++i) {
  35                 perm[i] = random(value);
  36                 perm[i] &= 0xFF;
  37                 perm[i + 256] = perm[i];
  38                 perm12[i] = perm[i] % 12;
  39                 perm12[i + 256] = perm12[i];
  40         }
  41 }
  42
  43
  44 float SimplexNoise::operator ()(const glm::vec3 &in) const noexcept {
  45         float skew = (in.x + in.y + in.z) * one_third;
  46
  47         glm::vec3 skewed(glm::floor(in + skew));
  48         float tr = (skewed.x + skewed.y + skewed.z) * one_sixth;
  49
  50         glm::vec3 unskewed(skewed - tr);
  51         glm::vec3 relative(in - unskewed);
  52
  53         glm::vec3 second, third;
  54
  55         if (relative.x >= relative.y) {
  56                 if (relative.y >= relative.z) {
  57                         second = { 1, 0, 0 };
  58                         third = { 1, 1, 0 };
  59                 } else if (relative.x >= relative.z) {
  60                         second = { 1, 0, 0 };
  61                         third = { 1, 0, 1 };
  62                 } else {
  63                         second = { 0, 0, 1 };
  64                         third = { 1, 0, 1 };
  65                 }
  66         } else if (relative.y < relative.z) {
  67                 second = { 0, 0, 1 };
  68                 third = { 0, 1, 1 };
  69         } else if (relative.x < relative.z) {
  70                 second = { 0, 1, 0 };
  71                 third = { 0, 1, 1 };
  72         } else {
  73                 second = { 0, 1, 0 };
  74                 third = { 1, 1, 0 };
  75         }
  76
  77         glm::vec3 offset[4] = {
  78                 in - unskewed,
  79                 relative - second + one_sixth,
  80                 relative - third + one_third,
  81                 relative - 0.5f,
  82         };
  83
  84         int index[3] = {
  85                 (int)(skewed.x) & 0xFF,
  86                 (int)(skewed.y) & 0xFF,
  87                 (int)(skewed.z) & 0xFF,
  88         };
  89
  90         float n = 0.0f;
  91
  92         // 0
  93         float t = 0.6f - dot(offset[0], offset[0]);
  94         if (t > 0.0f) {
  95                 t *= t;
  96                 int corner = Perm12(index[0] + Perm(index[1] + Perm(index[2])));
  97                 n += t * t * dot(Grad(corner), offset[0]);
  98         }
  99
 100         // 1
 101         t = 0.6f - dot(offset[1], offset[1]);
 102         if (t > 0.0f) {
 103                 t *= t;
 104                 int corner = Perm12(index[0] + int(second.x) + Perm(index[1] + int(second.y) + Perm(index[2] + int(second.z))));
 105                 n += t * t * dot(Grad(corner), offset[1]);
 106         }
 107
 108         // 2
 109         t = 0.6f - dot(offset[2], offset[2]);
 110         if (t > 0.0f) {
 111                 t *= t;
 112                 int corner = Perm12(index[0] + int(third.x) + Perm(index[1] + int(third.y) + Perm(index[2] + int(third.z))));
 113                 n += t * t * dot(Grad(corner), offset[2]);
 114         }
 115
 116         // 3
 117         t = 0.6f - dot(offset[3], offset[3]);
 118         if (t > 0.0f) {
 119                 t *= t;
 120                 int corner = Perm12(index[0] + 1 + Perm(index[1] + 1 + Perm(index[2] + 1)));
 121                 n += t * t * dot(Grad(corner), offset[3]);
 122         }
 123
 124         return 32.0f * n;
 125 }
 126
 127
 128 int SimplexNoise::Perm(int idx) const noexcept {
 129         return perm[idx];
 130 }
 131
 132 int SimplexNoise::Perm12(int idx) const noexcept {
 133         return perm12[idx];
 134 }
 135
 136 const glm::vec3 &SimplexNoise::Grad(int idx) const noexcept {
 137         return grad[idx];
 138 }
 139
 140
 141 WorleyNoise::WorleyNoise(unsigned int seed) noexcept
 142 : seed(seed)
 143 , num_points(8) {
 144
 145 }
 146
 147 float WorleyNoise::operator ()(const glm::vec3 &in) const noexcept {
 148         glm::vec3 center = floor(in);
 149
 150         float closest = 1.0f;  // cannot be farther away than 1.0
 151
 152         for (int z = -1; z <= 1; ++z) {
 153                 for (int y = -1; y <= 1; ++y) {
 154                         for (int x = -1; x <= 1; ++x) {
 155                                 glm::vec3 cube(center.x + x, center.y + y, center.z + z);
 156                                 unsigned int cube_rand =
 157                                         (unsigned(cube.x) * 130223) ^
 158                                         (unsigned(cube.y) * 159899) ^
 159                                         (unsigned(cube.z) * 190717) ^
 160                                         seed;
 161
 162                                 for (int i = 0; i < num_points; ++i) {
 163                                         glm::vec3 point(cube);
 164                                         cube_rand = 190667 * cube_rand + 109807;
 165                                         point.x += float(cube_rand % 262144) / 262144.0f;
 166                                         cube_rand = 135899 * cube_rand + 189169;
 167                                         point.y += float(cube_rand % 262144) / 262144.0f;
 168                                         cube_rand = 159739 * cube_rand + 112139;
 169                                         point.z += float(cube_rand % 262144) / 262144.0f;
 170
 171                                         glm::vec3 diff(in - point);
 172                                         float distance = sqrt(dot(diff, diff));
 173                                         if (distance < closest) {
 174                                                 closest = distance;
 175                                         }
 176                                 }
 177                         }
 178                 }
 179         }
 180
 181         // closest ranges (0, 1), so normalizing to (-1,1) is trivial
 182         // though heavily biased towards lower numbers
 183         return 2.0f * closest - 1.0f;
 184 }
 185
 186 }