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