src/rand/noise.cpp

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