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