layers/roughness.cpp - TerrainGen coverage

layers/roughness.cpp

Line	Branch	Exec	Source
1			#include "terraingen/layers/roughness.h"
2
3			#include "terraingen/noise/hash.h"
4			#include "terraingen/noise/perlin.h"
5			#include "terraingen/raster/cell_loops.h"
6			#include "terraingen/raster/grid.h"
7			#include "terraingen/util/checks.h"
8			#include "terraingen/util/log.h"
9
10			#include <algorithm>
11			#include <array>
12			#include <cmath>
13			#include <numbers>
14			#include <stdexcept>
15
16			namespace terraingen {
17
18			namespace {
19
20			constexpr int kMaxOctaves = 40;
21
22			/// Salts of the hashes giving each octave its offsets and seed, its 3D rotation, and its 3D
23			/// offset.
24			constexpr std::uint32_t kOctaveSalt = 0x524F5547u;
25			constexpr std::uint32_t kRotationSalt = 0x524F5433u;
26			constexpr std::uint32_t kOffsetZSalt = 7u;
27
28			/// The golden angle, in radians, by which each octave is turned from the last.
29			constexpr double kGoldenAngle = 2.39996322972865332;
30
31			/// A uniform in [0, 1024) from the top 16 bits of a hash, for lattice offsets.
32		418	double LatticeOffset(std::uint32_t hash) {
33		418	return static_cast<double>(hash >> 16) / 65536.0 * 1024.0;
34			}
35
36			/// The rotation matrix, row-major, of a uniform random rotation from three uniform hashes, by
37			/// Shoemake's (1992) unit quaternion.
38		16	std::array<double, 9> UniformRotation(const UVec3& turn) {
39		16	const double u1 = HashToUnitFloat(turn.x);
40		16	const double u2 = 2.0 * std::numbers::pi * static_cast<double>(HashToUnitFloat(turn.y));
41		16	const double u3 = 2.0 * std::numbers::pi * static_cast<double>(HashToUnitFloat(turn.z));
42		16	const double a = std::sqrt(1.0 - u1);
43		16	const double b = std::sqrt(u1);
44		16	const double qx = a * std::sin(u2);
45		16	const double qy = a * std::cos(u2);
46		16	const double qz = b * std::sin(u3);
47		16	const double qw = b * std::cos(u3);
48		16	return {1.0 - 2.0 * (qy * qy + qz * qz), 2.0 * (qx * qy - qz * qw),
49		16	2.0 * (qx * qz + qy * qw), 2.0 * (qx * qy + qz * qw),
50		16	1.0 - 2.0 * (qx * qx + qz * qz), 2.0 * (qy * qz - qx * qw),
51		16	2.0 * (qx * qz - qy * qw), 2.0 * (qy * qz + qx * qw),
52		16	1.0 - 2.0 * (qx * qx + qy * qy)};
53			}
54
55			/// `PerlinStructureFunction(r) / r^2` at r = 2^(k / 4) for k from `kFirstQuarter`, estimated from
56			/// 2^26 pairs each.
57			constexpr int kFirstQuarter = -24;
58			constexpr std::array<double, 33> kStructureOverSquare{
59			0.772366, 0.772241, 0.772129, 0.771929, 0.771444, 0.770975, 0.770202, 0.769434, 0.767832,
60			0.766070, 0.763032, 0.759311, 0.753675, 0.745950, 0.735267, 0.720564, 0.700391, 0.673043,
61			0.636947, 0.589660, 0.530574, 0.459215, 0.377706, 0.291690, 0.209524, 0.141536, 0.094237,
62			0.065236, 0.046153, 0.032397, 0.022947, 0.016231, 0.011443};
63			constexpr double kLargeLag = 0.183120;
64
65			/// `Perlin3StructureFunction(r) / r^2` on the same lags, estimated from 2^24 pairs each.
66			constexpr std::array<double, 33> kStructure3OverSquare{
67			0.494082, 0.494230, 0.494034, 0.493864, 0.493805, 0.493342, 0.493074, 0.492376, 0.491455,
68			0.490546, 0.488813, 0.486484, 0.483336, 0.478718, 0.472398, 0.463827, 0.452007, 0.435827,
69			0.414506, 0.386695, 0.351648, 0.308521, 0.258949, 0.205715, 0.153652, 0.108458, 0.074433,
70			0.051596, 0.036333, 0.025663, 0.018156, 0.012839, 0.009067};
71			constexpr double kLargeLag3 = 0.145198;
72
73			/// A structure function at lag `r` from its `table` of value over r^2, `large` beyond the table.
74		3951	double Interpolated(const std::array<double, 33>& table, double large, double r) {
75	2/2 ✓ Branch 0 taken 2 times. ✓ Branch 1 taken 3949 times.	3951	if (!(r > 0.0)) {
76		2	return 0.0;
77			}
78		3949	const double quarters = 4.0 * std::log2(r) - kFirstQuarter;
79	2/2 ✓ Branch 0 taken 562 times. ✓ Branch 1 taken 3387 times.	3949	if (quarters <= 0.0) {
80		562	return table.front() * r * r;
81			}
82		3387	const auto last = static_cast<double>(table.size() - 1);
83	2/2 ✓ Branch 0 taken 1880 times. ✓ Branch 1 taken 1507 times.	3387	if (quarters >= last) {
84		1880	return large;
85			}
86		1507	const auto k = static_cast<std::size_t>(quarters);
87		1507	const double t = quarters - static_cast<double>(k);
88		1507	return ((1.0 - t) * table[k] + t * table[k + 1]) * r * r;
89			}
90
91			} // namespace
92
93		54	void ValidateRoughnessOptions(const RoughnessOptions& options) {
94	1/1 ✓ Branch 2 taken 52 times.	54	RequirePositive(options.rmsSlope, "Roughness rmsSlope");
95	1/1 ✓ Branch 2 taken 50 times.	52	RequirePositive(options.baseline, "Roughness baseline");
96	4/4 ✓ Branch 0 taken 48 times. ✓ Branch 1 taken 2 times. ✓ Branch 2 taken 3 times. ✓ Branch 3 taken 45 times.	50	if (!(options.hurst > 0.0 && options.hurst <= 1.0)) {
97	1/1 ✓ Branch 2 taken 5 times.	5	throw std::invalid_argument("Roughness hurst must be in (0, 1]");
98			}
99	6/6 ✓ Branch 1 taken 43 times. ✓ Branch 2 taken 2 times. ✓ Branch 3 taken 41 times. ✓ Branch 4 taken 2 times. ✓ Branch 5 taken 6 times. ✓ Branch 6 taken 39 times.	86	if (!(std::isfinite(options.maxWavelength) && options.minWavelength > 0.0 &&
100	2/2 ✓ Branch 0 taken 2 times. ✓ Branch 1 taken 39 times.	41	options.minWavelength <= options.maxWavelength)) {
101			throw std::invalid_argument(
102	1/1 ✓ Branch 2 taken 6 times.	6	"Roughness wavelengths must satisfy 0 < minWavelength <= maxWavelength");
103			}
104	2/2 ✓ Branch 0 taken 2 times. ✓ Branch 1 taken 37 times.	39	if (std::log2(options.maxWavelength / options.minWavelength) >= kMaxOctaves) {
105	1/1 ✓ Branch 1 taken 2 times.	4	throw std::invalid_argument("Roughness wavelengths may span at most " +
106	2/2 ✓ Branch 3 taken 2 times. ✓ Branch 6 taken 2 times.	6	std::to_string(kMaxOctaves) + " octaves");
107			}
108		37	ValidateResolutionFade(options.fade);
109		35	}
110
111		2331	double PerlinStructureFunction(double r) {
112		2331	return Interpolated(kStructureOverSquare, kLargeLag, r);
113			}
114
115		1620	double Perlin3StructureFunction(double r) {
116		1620	return Interpolated(kStructure3OverSquare, kLargeLag3, r);
117			}
118
119		36	RoughnessLayer::RoughnessLayer(RoughnessOptions options, Surface surface)
120		36	: FormulaLayer(surface), options_(options) {
121	1/1 ✓ Branch 1 taken 27 times.	36	ValidateRoughnessOptions(options_);
122		27	const double shortest = options_.minWavelength * (1.0 - 1e-12);
123		27	double wavelength = options_.maxWavelength;
124	2/2 ✓ Branch 0 taken 201 times. ✓ Branch 1 taken 27 times.	228	for (std::uint32_t i = 0; wavelength >= shortest; ++i, wavelength *= 0.5) {
125	1/1 ✓ Branch 2 taken 201 times.	201	octaves_.push_back(MakeOctave(i, wavelength));
126			}
127		27	const double scale = options_.rmsSlope / RmsSlope(options_.baseline);
128	2/2 ✓ Branch 4 taken 201 times. ✓ Branch 5 taken 27 times.	228	for (Octave& octave : octaves_) {
129		201	octave.amplitude *= scale;
130			}
131		45	}
132
133		201	RoughnessLayer::Octave RoughnessLayer::MakeOctave(std::uint32_t i, double wavelength) const {
134		201	const UVec3 hash = Pcg3dHash({options_.seed, i, kOctaveSalt});
135		201	const double angle = kGoldenAngle * static_cast<double>(i);
136		201	Octave octave{wavelength,
137		201	std::pow(wavelength, options_.hurst),
138		201	std::cos(angle),
139		201	std::sin(angle),
140		201	LatticeOffset(hash.x),
141		201	LatticeOffset(hash.y),
142		201	hash.z};
143	2/2 ✓ Branch 1 taken 16 times. ✓ Branch 2 taken 185 times.	201	if (surface_.IsBody()) {
144		16	octave.rotation = UniformRotation(Pcg3dHash({options_.seed, i, kRotationSalt}));
145		16	octave.offsetZ = LatticeOffset(Pcg3dHash({hash.z, i, kOffsetZSalt}).x);
146			}
147		201	return octave;
148			}
149
150		30872304	double RoughnessLayer::OctaveNoise(const Octave& octave, double x, double y) const {
151	2/2 ✓ Branch 1 taken 94260 times. ✓ Branch 2 taken 30778044 times.	30872304	if (surface_.IsBody()) {
152	1/1 ✓ Branch 1 taken 94260 times.	94260	const Vec3d p = surface_.Position({x, y});
153		94260	const std::array<double, 9>& m = octave.rotation;
154		94260	const double scale = 1.0 / octave.wavelength;
155		377040	return static_cast<double>(PerlinNoise3Precise(
156	1/1 ✓ Branch 3 taken 94260 times.	94260	scale * (m[0] * p.x + m[1] * p.y + m[2] * p.z) + octave.offsetX,
157		94260	scale * (m[3] * p.x + m[4] * p.y + m[5] * p.z) + octave.offsetY,
158		188520	scale * (m[6] * p.x + m[7] * p.y + m[8] * p.z) + octave.offsetZ, octave.seed));
159			}
160		30778044	const double u = (octave.cos * x + octave.sin * y) / octave.wavelength + octave.offsetX;
161		30778044	const double v = (octave.cos * y - octave.sin * x) / octave.wavelength + octave.offsetY;
162		30778044	return static_cast<double>(PerlinNoise2Precise(u, v, octave.seed));
163			}
164
165		67	double RoughnessLayer::RmsSlope(double baseline) const {
166		67	double structure = 0.0;
167	2/2 ✓ Branch 4 taken 481 times. ✓ Branch 5 taken 67 times.	548	for (const Octave& octave : octaves_) {
168		481	const double lag = baseline / octave.wavelength;
169		481	structure +=
170		481	octave.amplitude * octave.amplitude *
171	2/2 ✓ Branch 1 taken 16 times. ✓ Branch 2 taken 465 times.	481	(surface_.IsBody() ? Perlin3StructureFunction(lag) : PerlinStructureFunction(lag));
172			}
173		67	return std::sqrt(structure) / baseline;
174			}
175
176		4	std::vector<double> RoughnessLayer::Wavelengths() const {
177		4	std::vector<double> wavelengths;
178	2/2 ✓ Branch 5 taken 24 times. ✓ Branch 6 taken 4 times.	28	for (const Octave& octave : octaves_) {
179	1/1 ✓ Branch 1 taken 24 times.	24	wavelengths.push_back(octave.wavelength);
180			}
181		4	return wavelengths;
182			}
183
184		6428524	float RoughnessLayer::HeightAt(float height, double x, double y, double spacing,
185			const Ground* /ground/) const {
186		6428524	double sum = 0.0;
187	2/2 ✓ Branch 4 taken 32788328 times. ✓ Branch 5 taken 4512500 times.	37300828	for (const Octave& octave : octaves_) {
188	1/1 ✓ Branch 1 taken 32788328 times.	32788328	const double weight = ResolutionWeight(octave.wavelength, spacing, options_.fade);
189	2/2 ✓ Branch 0 taken 1916024 times. ✓ Branch 1 taken 30872304 times.	32788328	if (weight <= 0.0) {
190		1916024	break;
191			}
192	1/1 ✓ Branch 1 taken 30872304 times.	30872304	sum += weight * octave.amplitude * OctaveNoise(octave, x, y);
193			}
194		6428524	return height + static_cast<float>(sum);
195			}
196
197			} // namespace terraingen
198