Celestia

1
#include <algorithm>
2
#include <array>
3
#include <cstddef>
4
#include <cstdint>
5
#include <numeric>
6

7
#include "randutils.h"
8

9
namespace celestia::math
10
{
11
namespace
12
{
13
std::mt19937 createRNG()
14
{
15
    std::random_device rd;
16
    std::uniform_int_distribution<std::uint_least32_t> dist{0, UINT32_C(0xffffffff)};
17
    constexpr std::size_t seedSize = std::mt19937::state_size * (std::mt19937::word_size / 32);
18
    std::array<std::uint_least32_t, seedSize> seedData;
19
    std::generate(seedData.begin(), seedData.end(), [&] { return dist(rd); });
20
    std::seed_seq rngSeed(seedData.begin(), seedData.end());
21
    return std::mt19937{ rngSeed };
22
}
23

24
// utility functions for Perlin noise
25
struct PerlinData
26
{
27
    static constexpr int TableSize = (1 << 8);
28
    static constexpr int TableMask = TableSize - 1;
29

30
    std::array<int, TableSize * 2> permutation{};
31
    std::array<float, TableSize> gradients1D{};
32
    std::array<Eigen::Vector2f, TableSize> gradients2D{};
33
    std::array<Eigen::Vector3f, TableSize> gradients3D{};
34

35
    PerlinData()
36
    {
37
        auto& rng = getRNG();
38
        auto permutationMid = permutation.begin() + TableSize;
39
        std::iota(permutation.begin(), permutationMid, 0);
40
        std::shuffle(permutation.begin(), permutationMid, rng);
41
        std::copy(permutation.begin(), permutationMid, permutationMid);
42

43
        std::generate(gradients1D.begin(), gradients1D.end(),
44
                      [&]() { return RealDists<float>::SignedUnit(rng); });
45

46
        std::generate(gradients2D.begin(), gradients2D.end(),
47
                      [&]() { return randomOnCircle<float>(rng); });
48

49
        std::generate(gradients3D.begin(), gradients3D.end(),
50
                      [&]() { return randomOnSphere<float>(rng); });
51
    }
52

53
    inline int permuteIndices(int x, int y) const
54
    {
55
        return permutation[permutation[x] + y];
56
    }
57

58
    inline int permuteIndices(int x, int y, int z) const
59
    {
60
        return permutation[permutation[permutation[x] + y] + z];
61
    }
62

63
    inline float gradientAt(int x) const
64
    {
65
        return gradients1D[permutation[x]];
66
    }
67

68
    inline const Eigen::Vector2f& gradientAt(int x, int y) const
69
    {
70
        return gradients2D[permutation[permutation[x] + y]];
71
    }
72

73
    inline const Eigen::Vector3f& gradientAt(int x, int y, int z) const
74
    {
75
        return gradients3D[permutation[permutation[permutation[x] + y] + z]];
76
    }
77
};
78

79
const PerlinData& getPerlinData()
80
{
81
    static PerlinData perlinData;
82
    return perlinData;
83
}
84

85
constexpr inline float smooth(float t)
86
{
87
    return t * t * (3.0f - 2.0f * t);
88
}
89

90
} // end of anonymous namespace
91

92
float noise(float arg)
93
{
94
    const auto& perlinData = getPerlinData();
95
    int x0 = static_cast<int>(std::floor(arg)) & PerlinData::TableMask;
96
    int x1 = (x0 + 1) & PerlinData::TableMask;
97
    float dx0 = arg - std::floor(arg);
98
    float dx1 = dx0 - 1.0f;
99
    float g0 = perlinData.gradientAt(x0);
100
    float g1 = perlinData.gradientAt(x1);
101
    return lerp(smooth(dx0), dx0 * g0, dx1 * g1);
102
}
103

104
float noise(const Eigen::Vector2f& arg)
105
{
106
    const auto& perlinData = getPerlinData();
107
    int x0[2] {static_cast<int>(std::floor(arg.x())) & PerlinData::TableMask,
108
               static_cast<int>(std::floor(arg.y())) & PerlinData::TableMask};
109
    int x1[2] {(x0[0] + 1) & PerlinData::TableMask,
110
               (x0[1] + 1) & PerlinData::TableMask};
111
    float dx0[2] {arg.x() - std::floor(arg.x()),
112
                  arg.y() - std::floor(arg.y())};
113
    float dx1[2] {dx0[0] - 1.0f, dx0[1] - 1.0f};
114
    const Eigen::Vector2f& g00 = perlinData.gradientAt(x0[0], x0[1]);
115
    const Eigen::Vector2f& g10 = perlinData.gradientAt(x1[0], x0[1]);
116
    const Eigen::Vector2f& g01 = perlinData.gradientAt(x0[0], x1[1]);
117
    const Eigen::Vector2f& g11 = perlinData.gradientAt(x1[0], x1[1]);
118
    Eigen::Vector2f v00{dx0[0], dx0[1]};
119
    Eigen::Vector2f v10{dx1[0], dx0[1]};
120
    Eigen::Vector2f v01{dx0[0], dx1[1]};
121
    Eigen::Vector2f v11{dx1[0], dx1[1]};
122
    float t[2] {smooth(dx0[0]), smooth(dx0[1])};
123
    float nx[2] {lerp(t[0], g00.dot(v00), g10.dot(v10)),
124
                 lerp(t[0], g01.dot(v01), g11.dot(v11))};
125
    return lerp(t[1], nx[0], nx[1]);
126
}
127

128
float noise(const Eigen::Vector3f& arg)
129
{
130
    const auto& perlinData = getPerlinData();
131
    int x0[3] {static_cast<int>(std::floor(arg.x())) & PerlinData::TableMask,
132
               static_cast<int>(std::floor(arg.y())) & PerlinData::TableMask,
133
               static_cast<int>(std::floor(arg.z())) & PerlinData::TableMask};
134
    int x1[3] {(x0[0] + 1) & PerlinData::TableMask,
135
               (x0[1] + 1) & PerlinData::TableMask,
136
               (x0[2] + 1) & PerlinData::TableMask};
137
    float dx0[3] {arg.x() - std::floor(arg.x()),
138
                  arg.y() - std::floor(arg.y()),
139
                  arg.z() - std::floor(arg.z())};
140
    float dx1[3] {dx0[0] - 1.0f, dx0[1] - 1.0f, dx0[2] - 1.0f};
141
    const Eigen::Vector3f& g000 = perlinData.gradientAt(x0[0], x0[1], x0[2]);
142
    const Eigen::Vector3f& g100 = perlinData.gradientAt(x1[0], x0[1], x0[2]);
143
    const Eigen::Vector3f& g010 = perlinData.gradientAt(x0[0], x1[1], x0[2]);
144
    const Eigen::Vector3f& g110 = perlinData.gradientAt(x1[0], x1[1], x0[2]);
145
    const Eigen::Vector3f& g001 = perlinData.gradientAt(x0[0], x0[1], x1[2]);
146
    const Eigen::Vector3f& g101 = perlinData.gradientAt(x1[0], x0[1], x1[2]);
147
    const Eigen::Vector3f& g011 = perlinData.gradientAt(x0[0], x1[1], x1[2]);
148
    const Eigen::Vector3f& g111 = perlinData.gradientAt(x1[0], x1[1], x1[2]);
149
    Eigen::Vector3f v000{dx0[0], dx0[1], dx0[2]};
150
    Eigen::Vector3f v100{dx1[0], dx0[1], dx0[2]};
151
    Eigen::Vector3f v010{dx0[0], dx1[1], dx0[2]};
152
    Eigen::Vector3f v110{dx1[0], dx1[1], dx0[2]};
153
    Eigen::Vector3f v001{dx0[0], dx0[1], dx1[2]};
154
    Eigen::Vector3f v101{dx1[0], dx0[1], dx1[2]};
155
    Eigen::Vector3f v011{dx0[0], dx1[1], dx1[2]};
156
    Eigen::Vector3f v111{dx1[0], dx1[1], dx1[2]};
157
    float t[3] {smooth(dx0[0]), smooth(dx0[1]), smooth(dx0[2])};
158
    float nx[4] {lerp(t[0], g000.dot(v000), g100.dot(v100)),
159
                 lerp(t[0], g010.dot(v010), g110.dot(v110)),
160
                 lerp(t[0], g001.dot(v001), g101.dot(v101)),
161
                 lerp(t[0], g011.dot(v011), g111.dot(v111))};
162
    float ny[2] {lerp(t[1], nx[0], nx[1]),
163
                 lerp(t[1], nx[2], nx[3])};
164
    return lerp(t[2], ny[0], ny[1]);
165
}
166

167
float turbulence(const Eigen::Vector2f& p, float freq)
168
{
169
    float t = 0.0f;
170
    while (freq >= 1.0f)
171
    {
172
        Eigen::Vector2f vec = freq * p;
173
        t += std::abs(noise(vec)) / freq;
174
        freq *= 0.5f;
175
    }
176

177
    return t;
178
}
179

180

181
float turbulence(const Eigen::Vector3f& p, float freq)
182
{
183
    float t = 0.0f;
184
    while (freq >= 1.0f)
185
    {
186
        Eigen::Vector3f vec = freq * p;
187
        t += std::abs(noise(vec)) / freq;
188
        freq *= 0.5f;
189
    }
190

191
    return t;
192
}
193

194

195
float fractalsum(const Eigen::Vector2f& p, float freq)
196
{
197
    float t = 0.0f;
198
    while (freq >= 1.0f)
199
    {
200
        Eigen::Vector2f vec = freq * p;
201
        t += noise(vec) / freq;
202
        freq *= 0.5f;
203
    }
204

205
    return t;
206
}
207

208

209
float fractalsum(const Eigen::Vector3f& p, float freq)
210
{
211
    float t = 0.0f;
212
    while (freq >= 1.0f)
213
    {
214
        Eigen::Vector3f vec = freq * p;
215
        t += noise(vec) / freq;
216
        freq *= 0.5f;
217
    }
218

219
    return t;
220
}
221

222
std::mt19937& getRNG()
223
{
224
    static std::mt19937 rng = createRNG();
225
    return rng;
226
}
227
} // end namespace celestia::math
228