TerrainGen coverage

terraingen/
Coverage:
low: ≥ 0%
medium: ≥ 75.0%
high: ≥ 90.0%
Lines:
314 / 318
98.7%
Branches:
369 / 410
90.0%
Functions:
33 / 33
100.0%
layers/boulders/boulder_options.cpp
Line Branch Exec Source
1 #include "terraingen/layers/boulders/boulder_options.h"
2
3 #include "terraingen/distribution/golombek_rapp.h"
4 #include "terraingen/distribution/power_law.h"
5 #include "terraingen/distribution/summed.h"
6 #include "terraingen/noise/worley.h"
7 #include "terraingen/util/checks.h"
8
9 #include <algorithm>
10 #include <cmath>
11 #include <numbers>
12 #include <stdexcept>
13 #include <string>
14 #include <utility>
15 #include <vector>
16
17 namespace terraingen {
18
19 namespace {
20
21 /// Whether `distribution` may be evaluated from `low` to `high` metres.
22 71 bool Covers(const SizeFrequencyDistribution& distribution, double low, double high) {
23
2/2
✓ Branch 1 taken 36 times.
✓ Branch 2 taken 35 times.
 107 return distribution.AllowsExtrapolation() ||
24
3/4
✓ Branch 1 taken 34 times.
✓ Branch 2 taken 2 times.
✓ Branch 4 taken 34 times.
✗ Branch 5 not taken.
 107 (distribution.ValidRange().min <= low && distribution.ValidRange().max >= high);
25 }
26
27 /// Throws unless `distribution` is a density valid over the boulder sizes.
28 73 void CheckDistribution(const SizeFrequencyDistribution& distribution, const BoulderSizes& sizes,
29 const char* what) {
30
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 71 times.
 73 if (distribution.Kind() != SfdKind::Density) {
31
3/3
✓ Branch 2 taken 2 times.
✓ Branch 5 taken 2 times.
✓ Branch 8 taken 2 times.
 6 throw std::invalid_argument(std::string(what) + " must be a density, not a rate");
32 }
33
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 69 times.
 71 if (!Covers(distribution, sizes.minDiameter, sizes.maxDiameter)) {
34
2/2
✓ Branch 1 taken 2 times.
✓ Branch 4 taken 2 times.
 4 throw std::invalid_argument(std::string(what) +
35 " must be valid from minDiameter to maxDiameter, or allow "
36
1/1
✓ Branch 2 taken 2 times.
 6 "extrapolation");
37 }
38 69 }
39
40 /// Checks the settings that only matter when craters throw blocks.
41 69 void CheckThrownBlocks(const BoulderEjecta& ejecta, const BoulderSizes& sizes) {
42
2/2
✓ Branch 1 taken 21 times.
✓ Branch 2 taken 48 times.
 69 if (ejecta.distribution) {
43 21 CheckDistribution(*ejecta.distribution, sizes, "The boulder ejecta distribution");
44
3/4
✓ Branch 0 taken 48 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
✓ Branch 3 taken 47 times.
 48 } else if (!(ejecta.coverage > 0.0 && ejecta.coverage < 1.0)) {
45
1/1
✓ Branch 2 taken 1 time.
 1 throw std::invalid_argument("Boulder ejecta coverage must be between 0 and 1");
46 }
47
5/6
✓ Branch 1 taken 66 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 time.
✓ Branch 4 taken 65 times.
✓ Branch 5 taken 1 time.
✓ Branch 6 taken 65 times.
 66 if (!(std::isfinite(ejecta.reach) && ejecta.reach > 1.0)) {
48
1/1
✓ Branch 2 taken 1 time.
 1 throw std::invalid_argument("Boulder ejecta need a reach over 1 radius");
49 }
50
5/6
✓ Branch 1 taken 65 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 time.
✓ Branch 5 taken 64 times.
✓ Branch 6 taken 1 time.
✓ Branch 7 taken 64 times.
 65 if (!(IsPositive(ejecta.decayRadii) && IsInUnitInterval(ejecta.floorShare))) {
51 throw std::invalid_argument(
52
1/1
✓ Branch 2 taken 1 time.
 1 "Boulder ejecta need a positive decayRadii and a floorShare in [0, 1]");
53 }
54
4/6
✓ Branch 0 taken 64 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 64 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 time.
✓ Branch 5 taken 63 times.
 128 if (!(ejecta.minFreshness >= 0.0 && ejecta.minFreshness < 1.0 &&
55
3/4
✓ Branch 1 taken 63 times.
✓ Branch 2 taken 1 time.
✗ Branch 4 not taken.
✓ Branch 5 taken 63 times.
 64 IsPositive(ejecta.lifetimeGyr) && IsNonNegative(ejecta.agingExponent))) {
56 throw std::invalid_argument("Boulder ejecta need minFreshness in [0, 1), a positive "
57
1/1
✓ Branch 2 taken 1 time.
 1 "lifetimeGyr and a non-negative agingExponent");
58 }
59
6/6
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 60 times.
✓ Branch 5 taken 1 time.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 1 time.
✓ Branch 8 taken 62 times.
 63 if (ejecta.gravity && !IsPositive(*ejecta.gravity)) {
60
1/1
✓ Branch 2 taken 1 time.
 1 throw std::invalid_argument("Boulder ejecta gravity must be positive");
61 }
62 62 }
63
64 } // namespace
65
66 418043 bool BoulderModelThreshold::IsModel(double diameter, double spacing) const {
67
6/6
✓ Branch 1 taken 192600 times.
✓ Branch 2 taken 225443 times.
✓ Branch 4 taken 186355 times.
✓ Branch 5 taken 6245 times.
✓ Branch 6 taken 6 times.
✓ Branch 7 taken 411792 times.
 829847 const bool large = (minDiameter && diameter >= *minDiameter) ||
68
4/4
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 2 times.
✓ Branch 4 taken 3 times.
✓ Branch 5 taken 1 time.
 411804 (minSamples && spacing > 0.0 && diameter >= *minSamples * spacing);
69
6/6
✓ Branch 0 taken 6248 times.
✓ Branch 1 taken 411795 times.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 6246 times.
✓ Branch 6 taken 1 time.
✓ Branch 7 taken 1 time.
 424290 return large && !(maxDiameter && diameter >= *maxDiameter) &&
70
5/6
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 6245 times.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 time.
✓ Branch 7 taken 1 time.
 424290 !(maxSamples && spacing > 0.0 && diameter >= *maxSamples * spacing);
71 }
72
73 1010 double BoulderEjecta::MinCraterDiameter() const {
74
2/2
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 994 times.
 1010 const double excavated = regolithThickness ? *regolithThickness / excavationDepthRatio : 0.0;
75 1010 return std::fmax(minCraterDiameter, excavated);
76 }
77
78 9 double BoulderThrowDistance(double craterDiameter, double regolithThickness, double gravity) {
79
2/2
✓ Branch 2 taken 25 times.
✓ Branch 3 taken 7 times.
 32 for (const double value : {craterDiameter, regolithThickness, gravity}) {
80
1/1
✓ Branch 2 taken 23 times.
 25 RequirePositive(value, "Crater diameter, regolith thickness and gravity");
81 }
82 7 return 138.5 / gravity * std::pow(craterDiameter / regolithThickness, 0.7086);
83 }
84
85 400 void ValidateBoulderSizes(const BoulderSizes& sizes) {
86
5/6
✓ Branch 1 taken 399 times.
✓ Branch 2 taken 1 time.
✓ Branch 4 taken 399 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 398 times.
 799 if (!(IsPositive(sizes.minDiameter) && std::isfinite(sizes.maxDiameter) &&
87
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 398 times.
 399 sizes.maxDiameter > sizes.minDiameter)) {
88
1/1
✓ Branch 2 taken 2 times.
 2 throw std::invalid_argument("Boulder diameters must satisfy 0 < minDiameter < maxDiameter");
89 }
90
2/2
✓ Branch 1 taken 52 times.
✓ Branch 2 taken 346 times.
 398 if (sizes.distribution) {
91 52 CheckDistribution(*sizes.distribution, sizes, "The boulder distribution");
92
4/4
✓ Branch 0 taken 344 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 342 times.
 346 } else if (!(sizes.coverage > 0.0 && sizes.coverage < 1.0)) {
93
1/1
✓ Branch 2 taken 4 times.
 4 throw std::invalid_argument("Boulder coverage must be between 0 and 1");
94 }
95 392 }
96
97 393 void ValidateBoulderPatches(const BoulderPatches& patches) {
98
7/8
✓ Branch 1 taken 393 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 391 times.
✓ Branch 4 taken 2 times.
✓ Branch 6 taken 1 time.
✓ Branch 7 taken 390 times.
✓ Branch 8 taken 3 times.
✓ Branch 9 taken 390 times.
 393 if (!(std::isfinite(patches.spread) && patches.spread >= 1.0 && IsPositive(patches.size))) {
99
1/1
✓ Branch 2 taken 3 times.
 3 throw std::invalid_argument("Boulder patches need spread >= 1 and a positive size");
100 }
101 390 }
102
103 381 void ValidateBoulderEjecta(const BoulderEjecta& ejecta, const BoulderSizes& sizes) {
104
2/2
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 380 times.
 381 if (!IsNonNegative(ejecta.minCraterDiameter)) {
105
1/1
✓ Branch 2 taken 1 time.
 1 throw std::invalid_argument("Boulder ejecta minCraterDiameter must not be negative");
106 }
107
6/6
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 374 times.
✓ Branch 5 taken 1 time.
✓ Branch 6 taken 5 times.
✓ Branch 7 taken 1 time.
✓ Branch 8 taken 379 times.
 380 if (ejecta.regolithThickness && !IsPositive(*ejecta.regolithThickness)) {
108
1/1
✓ Branch 2 taken 1 time.
 1 throw std::invalid_argument("Boulder ejecta regolithThickness must be positive");
109 }
110
3/4
✓ Branch 0 taken 378 times.
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
✓ Branch 3 taken 378 times.
 379 if (!(ejecta.excavationDepthRatio > 0.0 && ejecta.excavationDepthRatio <= 1.0)) {
111
1/1
✓ Branch 2 taken 1 time.
 1 throw std::invalid_argument("Boulder ejecta excavationDepthRatio must be in (0, 1]");
112 }
113
2/2
✓ Branch 1 taken 69 times.
✓ Branch 2 taken 309 times.
 378 if (ejecta.MinCraterDiameter() > 0.0) {
114 69 CheckThrownBlocks(ejecta, sizes);
115 }
116 371 }
117
118 388 void ValidateBoulderBedrock(const BoulderBedrock& bedrock) {
119
1/1
✓ Branch 2 taken 386 times.
 388 RequirePositive(bedrock.abundance, "Boulder bedrockAbundance");
120
2/2
✓ Branch 1 taken 28 times.
✓ Branch 2 taken 358 times.
 386 if (bedrock.exposure) {
121 28 ValidateExposureParams(*bedrock.exposure);
122 }
123
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 383 times.
 385 if (!IsPositive(bedrock.exposureSpacing)) {
124
1/1
✓ Branch 2 taken 2 times.
 2 throw std::invalid_argument("Boulder exposureSpacing must be positive");
125 }
126
1/1
✓ Branch 2 taken 383 times.
 383 RequirePositive(bedrock.meanSize, "Boulder bedrock meanSize");
127 383 }
128
129 373 void ValidateBoulderPose(const BoulderPose& pose) {
130
6/6
✓ Branch 0 taken 372 times.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 371 times.
✓ Branch 3 taken 1 time.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 371 times.
 744 if (!(pose.minHeightRatio > 0.0 && pose.maxHeightRatio >= pose.minHeightRatio &&
131
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 371 times.
 371 std::isfinite(pose.maxHeightRatio))) {
132
1/1
✓ Branch 2 taken 2 times.
 2 throw std::invalid_argument("Boulder height ratios must satisfy 0 < min <= max");
133 }
134
5/6
✓ Branch 0 taken 371 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 370 times.
✓ Branch 3 taken 1 time.
✓ Branch 4 taken 1 time.
✓ Branch 5 taken 369 times.
 371 if (!(pose.minBurial >= 0.0 && pose.maxBurial >= pose.minBurial && pose.maxBurial < 1.0)) {
135
1/1
✓ Branch 2 taken 2 times.
 2 throw std::invalid_argument("Boulder burial must satisfy 0 <= min <= max < 1");
136 }
137
2/2
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 368 times.
 369 if (!IsNonNegative(pose.maxLean)) {
138
1/1
✓ Branch 2 taken 1 time.
 1 throw std::invalid_argument("Boulder maxLean must not be negative");
139 }
140
3/4
✓ Branch 0 taken 368 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
✓ Branch 3 taken 367 times.
 368 if (!(pose.burialByWear >= 0.0 && pose.burialByWear <= 1.0)) {
141
1/1
✓ Branch 2 taken 1 time.
 1 throw std::invalid_argument("Boulder burialByWear must be in [0, 1]");
142 }
143
4/6
✓ Branch 0 taken 367 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 367 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 time.
✓ Branch 5 taken 366 times.
 367 if (!(pose.sideShare >= 0.0 && pose.endShare >= 0.0 && pose.sideShare + pose.endShare <= 1.0)) {
144 throw std::invalid_argument(
145
1/1
✓ Branch 2 taken 1 time.
 1 "Boulder sideShare and endShare must not be negative and must sum to at most 1");
146 }
147
2/4
✓ Branch 0 taken 366 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 366 times.
 366 if (!(pose.maxTilt >= 0.0 && pose.maxTilt <= std::numbers::pi)) {
148
0/1
✗ Branch 2 not taken.
 ✗ throw std::invalid_argument("Boulder maxTilt must be in [0, pi]");
149 }
150
8/8
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 350 times.
✓ Branch 4 taken 15 times.
✓ Branch 5 taken 1 time.
✓ Branch 7 taken 1 time.
✓ Branch 8 taken 14 times.
✓ Branch 9 taken 2 times.
✓ Branch 10 taken 364 times.
 366 if (pose.maxSlope && !(*pose.maxSlope > 0.0 && *pose.maxSlope < 0.5 * std::numbers::pi)) {
151
1/1
✓ Branch 2 taken 2 times.
 2 throw std::invalid_argument("Boulder maxSlope must be in (0, pi / 2)");
152 }
153
4/4
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 356 times.
✓ Branch 2 taken 1 time.
✓ Branch 3 taken 363 times.
 372 if (pose.rollDistance &&
154
4/6
✓ Branch 2 taken 7 times.
✓ Branch 3 taken 1 time.
✓ Branch 5 taken 7 times.
✗ Branch 6 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 7 times.
 372 !(pose.maxSlope && *pose.rollDistance > 0.0 && std::isfinite(*pose.rollDistance))) {
155 throw std::invalid_argument("Boulder rollDistance must be positive and finite, with a "
156
1/1
✓ Branch 2 taken 1 time.
 1 "maxSlope");
157 }
158 363 }
159
160 365 void ValidateBoulderFillet(const BoulderFillet& fillet) {
161
3/4
✓ Branch 0 taken 365 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 363 times.
 365 if (!(fillet.heightShare >= 0.0 && fillet.heightShare <= 1.0)) {
162
1/1
✓ Branch 2 taken 2 times.
 2 throw std::invalid_argument("Boulder fillet heightShare must be in [0, 1]");
163 }
164
2/2
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 362 times.
 363 if (!IsPositive(fillet.widthRatio)) {
165
1/1
✓ Branch 2 taken 1 time.
 1 throw std::invalid_argument("Boulder fillet widthRatio must be positive");
166 }
167 362 }
168
169 391 void ValidateBoulderClusters(const BoulderClusters& clusters) {
170
2/4
✓ Branch 0 taken 391 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 391 times.
 391 if (!(clusters.share >= 0.0 && clusters.share <= 1.0)) {
171
0/1
✗ Branch 2 not taken.
 ✗ throw std::invalid_argument("Boulder cluster share must be in [0, 1]");
172 }
173
2/2
✓ Branch 0 taken 11 times.
✓ Branch 1 taken 380 times.
 391 if (clusters.model == BoulderClusterModel::Thomas) {
174
6/6
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 9 times.
✓ Branch 3 taken 1 time.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 9 times.
 20 if (!(clusters.radius > 0.0 && clusters.radius <= 2.0 * clusters.spacing &&
175
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 9 times.
 9 std::isfinite(clusters.spacing))) {
176 throw std::invalid_argument(
177
1/1
✓ Branch 2 taken 2 times.
 2 "Thomas boulder clusters must satisfy 0 < radius <= 2 spacing");
178 }
179 9 return;
180 }
181
5/6
✓ Branch 0 taken 380 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 378 times.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 378 times.
 758 if (!(clusters.radius > 0.0 && clusters.radius < clusters.spacing &&
182
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 378 times.
 378 std::isfinite(clusters.spacing))) {
183
1/1
✓ Branch 2 taken 2 times.
 2 throw std::invalid_argument("Boulder clusters must satisfy 0 < radius < spacing");
184 }
185 }
186
187 2054165 double ClusterBump(double distance, const BoulderClusters& clusters) {
188 2054165 const double r = distance / clusters.radius;
189 2054165 return std::exp(-0.5 * r * r);
190 }
191
192 28 double MeanClusterBump(const BoulderClusters& clusters) {
193 28 constexpr int kSamples = 256;
194 28 constexpr double kCells = 16.0;
195 28 double sum = 0.0;
196
2/2
✓ Branch 0 taken 7168 times.
✓ Branch 1 taken 28 times.
 7196 for (int j = 0; j < kSamples; ++j) {
197
2/2
✓ Branch 0 taken 1835008 times.
✓ Branch 1 taken 7168 times.
 1842176 for (int i = 0; i < kSamples; ++i) {
198 1835008 const double u = (i + 0.5) * kCells / kSamples;
199 1835008 const double v = (j + 0.5) * kCells / kSamples;
200 const WorleyDistances d =
201 1835008 Worley2(static_cast<float>(u), static_cast<float>(v), 0x434C5553u);
202 1835008 sum += ClusterBump(static_cast<double>(d.f1) * clusters.spacing, clusters);
203 }
204 }
205 28 return sum / (kSamples * kSamples);
206 }
207
208 26 double MeanThomasField(const BoulderClusters& clusters, bool body) {
209 26 const double sigma = clusters.radius;
210 26 const double reach = kThomasReachRadii;
211 26 const double perArea = 1.0 / (clusters.spacing * clusters.spacing);
212
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 23 times.
 26 if (body) {
213 3 const double perVolume = perArea / (std::sqrt(2.0 * std::numbers::pi) * sigma);
214 3 const double inside =
215 std::erf(reach / std::numbers::sqrt2) -
216 std::sqrt(2.0 / std::numbers::pi) * reach * std::exp(-0.5 * reach * reach);
217 3 return perVolume * std::pow(2.0 * std::numbers::pi * sigma * sigma, 1.5) * inside;
218 }
219 23 return perArea * 2.0 * std::numbers::pi * sigma * sigma *
220 23 (1.0 - std::exp(-0.5 * reach * reach));
221 }
222
223 819 double MaxClusterAbundance(const BoulderClusters& clusters) {
224
2/2
✓ Branch 0 taken 778 times.
✓ Branch 1 taken 41 times.
 819 if (clusters.share <= 0.0) {
225 778 return 1.0;
226 }
227
2/2
✓ Branch 0 taken 23 times.
✓ Branch 1 taken 18 times.
 41 if (clusters.model == BoulderClusterModel::Nearest) {
228 23 return 1.0 - clusters.share + clusters.share / MeanClusterBump(clusters);
229 }
230 18 constexpr double kChance = 1e-6;
231 18 const double ratio = clusters.radius / clusters.spacing;
232 18 const double mean = 9.0 * std::numbers::pi * ratio * ratio;
233 18 double term = std::exp(-mean);
234 18 double above = 1.0 - term;
235 18 int crowd = 0;
236
2/2
✓ Branch 0 taken 522 times.
✓ Branch 1 taken 18 times.
 540 while (above >= kChance) {
237 522 ++crowd;
238 522 term *= mean / crowd;
239 522 above -= term;
240 }
241 18 return 1.0 - clusters.share +
242 18 clusters.share * std::fmax(1.0, crowd) / MeanThomasField(clusters);
243 }
244
245 389 void ValidateBoulderFragmentation(const BoulderFragmentation& fragmentation) {
246
1/1
✓ Branch 2 taken 389 times.
 389 RequirePositive(fragmentation.minDiameter, "Boulder fragmentation minDiameter");
247
3/4
✓ Branch 0 taken 389 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 388 times.
✓ Branch 3 taken 1 time.
 389 if (!(fragmentation.splitWear >= 0.0 && fragmentation.rubbleWear >= fragmentation.splitWear &&
248
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 388 times.
 388 fragmentation.rubbleWear <= 1.0)) {
249 throw std::invalid_argument(
250
1/1
✓ Branch 2 taken 1 time.
 1 "Boulder fragmentation must satisfy 0 <= splitWear <= rubbleWear <= 1");
251 }
252
3/4
✓ Branch 0 taken 386 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 386 times.
✗ Branch 3 not taken.
 388 if (!(fragmentation.minPieces >= 2 && fragmentation.maxPieces >= fragmentation.minPieces &&
253
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 386 times.
 386 fragmentation.maxPieces <= 12)) {
254
1/1
✓ Branch 2 taken 2 times.
 2 throw std::invalid_argument("Boulder fragments must satisfy 2 <= min <= max <= 12 pieces");
255 }
256
1/1
✓ Branch 2 taken 386 times.
 386 RequireNonNegative(fragmentation.gapRatio, "Boulder fragmentation gapRatio");
257
1/1
✓ Branch 2 taken 386 times.
 386 RequireNonNegative(fragmentation.spreadRatio, "Boulder fragmentation spreadRatio");
258
3/4
✓ Branch 0 taken 386 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
✓ Branch 3 taken 385 times.
 386 if (!(fragmentation.slabShare >= 0.0 && fragmentation.slabShare <= 1.0)) {
259
1/1
✓ Branch 2 taken 1 time.
 1 throw std::invalid_argument("Boulder fragmentation slabShare must be in [0, 1]");
260 }
261 385 }
262
263 371 void ValidateBoulderVentifacts(const BoulderVentifacts& ventifacts) {
264
4/4
✓ Branch 0 taken 739 times.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 737 times.
✓ Branch 3 taken 2 times.
 740 const auto unit = [](double v) { return v >= 0.0 && v <= 1.0; };
265
6/6
✓ Branch 1 taken 369 times.
✓ Branch 2 taken 2 times.
✓ Branch 4 taken 1 time.
✓ Branch 5 taken 368 times.
✓ Branch 6 taken 3 times.
✓ Branch 7 taken 368 times.
 371 if (!(unit(ventifacts.share) && unit(ventifacts.keelShare))) {
266
1/1
✓ Branch 2 taken 3 times.
 3 throw std::invalid_argument("Boulder ventifact share and keelShare must be in [0, 1]");
267 }
268
4/4
✓ Branch 0 taken 367 times.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 365 times.
✓ Branch 3 taken 2 times.
 368 if (!(ventifacts.minAngle > 0.0 && ventifacts.minAngle <= ventifacts.maxAngle &&
269
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 364 times.
 365 ventifacts.maxAngle < 0.5 * std::numbers::pi)) {
270 throw std::invalid_argument(
271
1/1
✓ Branch 2 taken 4 times.
 4 "Boulder ventifact facet angles must satisfy 0 < minAngle <= maxAngle < 90 deg");
272 }
273
6/6
✓ Branch 1 taken 363 times.
✓ Branch 2 taken 1 time.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 355 times.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 362 times.
 735 if (!(std::isfinite(ventifacts.windAzimuth) &&
274
2/2
✓ Branch 3 taken 1 time.
✓ Branch 4 taken 7 times.
 371 (!ventifacts.secondWindAzimuth || std::isfinite(*ventifacts.secondWindAzimuth)))) {
275
1/1
✓ Branch 2 taken 2 times.
 2 throw std::invalid_argument("Boulder ventifact wind azimuths must be finite");
276 }
277 362 }
278
279 2 BoulderVentifacts GaleVentifacts() {
280 2 BoulderVentifacts ventifacts;
281 2 ventifacts.share = 0.5 * 0.78;
282 2 ventifacts.keelShare = 12.0 / 39.0;
283 2 return ventifacts;
284 }
285
286 20854 double FilletReach(const Boulder& boulder, const BoulderFillet& fillet) {
287 20854 return (0.5 + fillet.widthRatio) * boulder.diameter;
288 }
289
290 100500 double FilletHeight(const Boulder& boulder, double distance, const BoulderFillet& fillet) {
291 100500 const double radius = 0.5 * boulder.diameter;
292 100500 const double top = fillet.heightShare * boulder.ExposedHeight();
293
2/2
✓ Branch 0 taken 25504 times.
✓ Branch 1 taken 74996 times.
 100500 if (distance <= radius) {
294 25504 return top;
295 }
296 74996 const double t = (distance - radius) / (fillet.widthRatio * boulder.diameter);
297
2/2
✓ Branch 0 taken 35263 times.
✓ Branch 1 taken 39733 times.
 74996 return t >= 1.0 ? 0.0 : top * (1.0 - t) * (1.0 - t);
298 }
299
300 362 void ValidateBoulderRendering(const BoulderRendering& rendering) {
301 362 const BoulderModelThreshold& models = rendering.models;
302 1798 for (const auto& threshold :
303
2/2
✓ Branch 2 taken 1440 times.
✓ Branch 3 taken 358 times.
 2160 {models.minDiameter, models.minSamples, models.maxDiameter, models.maxSamples}) {
304
2/2
✓ Branch 1 taken 73 times.
✓ Branch 2 taken 1367 times.
 1440 if (threshold) {
305
1/1
✓ Branch 3 taken 69 times.
 73 RequirePositive(*threshold, "Boulder model thresholds");
306 }
307 }
308 358 const BoulderReliefRange& range = rendering.reliefRange;
309
1/1
✓ Branch 2 taken 358 times.
 358 RequireNonNegative(range.minDiameter, "Boulder relief minDiameter");
310
6/6
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 356 times.
✓ Branch 4 taken 1 time.
✓ Branch 5 taken 1 time.
✓ Branch 6 taken 1 time.
✓ Branch 7 taken 357 times.
 358 if (range.maxDiameter && !(*range.maxDiameter > range.minDiameter)) {
311
1/1
✓ Branch 2 taken 1 time.
 1 throw std::invalid_argument("Boulder relief maxDiameter must exceed its minDiameter");
312 }
313 357 ValidateBoulderShapeOptions(rendering.shapes);
314 357 }
315
316 396 void ValidateBoulderLayerOptions(const BoulderLayerOptions& options) {
317 396 ValidateBoulderSizes(options.sizes);
318
2/2
✓ Branch 1 taken 12 times.
✓ Branch 2 taken 379 times.
 391 if (options.densityMap) {
319 12 ValidateDensityMap(*options.densityMap);
320 }
321
2/2
✓ Branch 1 taken 13 times.
✓ Branch 2 taken 378 times.
 391 if (options.slopeMap) {
322 13 ValidateSizeSlopeMap(*options.slopeMap);
323
2/2
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 12 times.
 13 if (!options.sizes.distribution) {
324 throw std::invalid_argument(
325
1/1
✓ Branch 2 taken 1 time.
 1 "A boulder size slope map needs a size distribution, not Golombek-Rapp sizes");
326 }
327 }
328 390 ValidateBoulderPatches(options.patches);
329 389 ValidateBoulderClusters(options.clusters);
330 386 ValidateBoulderFragmentation(options.fragmentation);
331 385 ValidateBoulderBedrock(options.bedrock);
332
4/4
✓ Branch 1 taken 332 times.
✓ Branch 2 taken 50 times.
✓ Branch 3 taken 3 times.
✓ Branch 4 taken 379 times.
 714 if (!options.sizes.distribution &&
333
2/2
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 329 times.
 332 options.sizes.coverage * MaxBackgroundAbundance(options) >= 1.0) {
334 throw std::invalid_argument("Boulder coverage times the patch spread, cluster and bedrock "
335
1/1
✓ Branch 2 taken 3 times.
 3 "abundance must stay below 1");
336 }
337 379 ValidateBoulderEjecta(options.ejecta, options.sizes);
338 370 ValidateBoulderPose(options.pose);
339 362 ValidateBoulderFillet(options.fillet);
340 362 ValidateBoulderVentifacts(options.ventifacts);
341 361 ValidateResolutionFade(options.fade);
342 360 ValidateBoulderRendering(options.rendering);
343 356 }
344
345 namespace {
346
347 /// Coverages the Golombek-Rapp measures are tabulated at, log-spaced.
348 constexpr int kCoverageSteps = 512;
349 constexpr double kLeastCoverage = 1e-6;
350 constexpr double kMostCoverage = 0.999;
351
352 /// What `map` measures of `distribution`, with the areal fraction of Golombek-Rapp boulders in
353 /// closed form so that it needs no largest diameter.
354 6150 double BoulderMeasure(const DensityMap& map, const SizeFrequencyDistribution& distribution) {
355
1/2
✓ Branch 0 taken 6150 times.
✗ Branch 1 not taken.
 6150 const auto* rapp = dynamic_cast<const GolombekRappSfd*>(&distribution);
356
4/4
✓ Branch 0 taken 6144 times.
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 5632 times.
✓ Branch 3 taken 512 times.
 6150 if (rapp != nullptr && map.measure == DensityMeasure::ArealFraction) {
357 const double beyond =
358
1/2
✓ Branch 1 taken 5632 times.
✗ Branch 2 not taken.
 5632 std::isinf(map.maxDiameter) ? 0.0 : rapp->FractionalArea(map.maxDiameter);
359 5632 return rapp->FractionalArea(map.minDiameter) - beyond;
360 }
361 518 return MapMeasure(map, distribution);
362 }
363
364 } // namespace
365
366 39 BoulderMapAbundance::BoulderMapAbundance(const BoulderLayerOptions& options) {
367
2/2
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 38 times.
 39 if (!options.densityMap) {
368
1/1
✓ Branch 2 taken 1 time.
 1 throw std::invalid_argument("Boulder map abundance needs a density map");
369 }
370 38 map_ = *options.densityMap;
371
1/1
✓ Branch 1 taken 38 times.
 38 ValidateDensityMap(map_);
372 38 const BoulderSizes& sizes = options.sizes;
373
6/6
✓ Branch 0 taken 37 times.
✓ Branch 1 taken 1 time.
✓ Branch 3 taken 1 time.
✓ Branch 4 taken 36 times.
✓ Branch 5 taken 1 time.
✓ Branch 6 taken 37 times.
 38 if (map_.total && options.ejecta.MinCraterDiameter() > 0.0) {
374
2/2
✓ Branch 1 taken 1 time.
✓ Branch 5 taken 1 time.
 1 ejectaMeasure_ = BoulderMeasure(map_, *EjectaSizes(options));
375 }
376
6/6
✓ Branch 1 taken 33 times.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 21 times.
✓ Branch 4 taken 12 times.
✓ Branch 5 taken 26 times.
✓ Branch 6 taken 12 times.
 38 if (sizes.distribution || map_.measure == DensityMeasure::Scale) {
377
3/3
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 21 times.
✓ Branch 5 taken 5 times.
 26 measure_ = sizes.distribution ? BoulderMeasure(map_, *sizes.distribution) : 1.0;
378
3/6
✓ Branch 0 taken 26 times.
✗ Branch 1 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 26 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 26 times.
 26 if (!(measure_ > 0.0 && std::isfinite(measure_))) {
379
0/1
✗ Branch 2 not taken.
 ✗ throw std::invalid_argument("A boulder density map must measure some of the sizes");
380 }
381 } else {
382 12 coverage_ = sizes.coverage;
383
2/2
✓ Branch 0 taken 6144 times.
✓ Branch 1 taken 12 times.
 6156 for (int i = 0; i < kCoverageSteps; ++i) {
384 const double k =
385 6144 kLeastCoverage * std::pow(kMostCoverage / kLeastCoverage,
386 6144 static_cast<double>(i) / (kCoverageSteps - 1));
387
1/1
✓ Branch 1 taken 6144 times.
 6144 coverages_.push_back(k);
388
1/1
✓ Branch 1 taken 6144 times.
 6144 measures_.push_back(BoulderMeasure(
389
2/2
✓ Branch 2 taken 6144 times.
✓ Branch 5 taken 6144 times.
 12288 map_, GolombekRappSfd(k, sizes.qConstant, sizes.qCoverageCoefficient)));
390 }
391
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 12 times.
 12 if (!(measures_.back() > 0.0)) {
392
0/1
✗ Branch 2 not taken.
 ✗ throw std::invalid_argument("A boulder density map must measure some of the sizes");
393 }
394 }
395
1/1
✓ Branch 2 taken 38 times.
 38 const double top = map_.field->Max();
396
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 35 times.
 38 const double largest = map_.measure == DensityMeasure::Count ? top / map_.areaUnit : top;
397
1/1
✓ Branch 1 taken 38 times.
 38 max_ = std::fmax(1.0, At(largest, 0.0));
398 41 }
399
400 2777 double BoulderMapAbundance::Coverage(double target) const {
401
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2775 times.
 2777 if (target <= measures_.front()) {
402 2 return coverages_.front() * target / measures_.front();
403 }
404
2/2
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 2774 times.
 2775 if (target >= measures_.back()) {
405 1 return coverages_.back();
406 }
407
1/1
✓ Branch 3 taken 2774 times.
 2774 const auto above = std::lower_bound(measures_.begin(), measures_.end(), target);
408 2774 const auto i = static_cast<std::size_t>(above - measures_.begin());
409 2774 const double t = (target - measures_[i - 1]) / (measures_[i] - measures_[i - 1]);
410 2774 return coverages_[i - 1] + t * (coverages_[i] - coverages_[i - 1]);
411 }
412
413 1083759 double BoulderMapAbundance::At(double value, double ejecta) const {
414
2/2
✓ Branch 1 taken 30001 times.
✓ Branch 2 taken 1053758 times.
 1083759 if (std::isnan(value)) {
415 30001 return 1.0;
416 }
417
2/2
✓ Branch 0 taken 1050975 times.
✓ Branch 1 taken 2783 times.
 1053758 if (map_.measure == DensityMeasure::Scale) {
418 1050975 return std::fmax(0.0, value);
419 }
420 2783 const double target = value - ejecta * ejectaMeasure_;
421
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 2782 times.
 2783 if (!(target > 0.0)) {
422 1 return 0.0;
423 }
424
2/2
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 2777 times.
 2782 return coverages_.empty() ? target / measure_ : Coverage(target) / coverage_;
425 }
426
427 36 double MeanPatchAbundance(const BoulderPatches& patches) {
428
2/2
✓ Branch 0 taken 30 times.
✓ Branch 1 taken 6 times.
 36 if (!(patches.spread > 1.0)) {
429 30 return 1.0;
430 }
431 6 constexpr int kSamples = 256;
432 6 constexpr double kPatches = 32.0;
433 6 double sum = 0.0;
434
2/2
✓ Branch 0 taken 1536 times.
✓ Branch 1 taken 6 times.
 1542 for (int j = 0; j < kSamples; ++j) {
435
2/2
✓ Branch 0 taken 393216 times.
✓ Branch 1 taken 1536 times.
 394752 for (int i = 0; i < kSamples; ++i) {
436 393216 const double x = (i + 0.5) * kPatches / kSamples;
437 393216 const double y = (j + 0.5) * kPatches / kSamples;
438
1/1
✓ Branch 2 taken 393216 times.
 393216 sum += std::pow(patches.spread, PatchNoise({}, x, y, 1.0, 0x50415443u));
439 }
440 }
441 6 return sum / (kSamples * kSamples);
442 }
443
444 815 double MaxBackgroundAbundance(const BoulderLayerOptions& options) {
445 815 const double a = options.bedrock.abundance;
446 815 double patches = options.patches.spread;
447
2/2
✓ Branch 1 taken 64 times.
✓ Branch 2 taken 751 times.
 815 double bedrock = options.bedrock.exposure ? std::fmax(1.0, a) : 1.0;
448 815 double map = 1.0;
449
2/2
✓ Branch 1 taken 26 times.
✓ Branch 2 taken 789 times.
 815 if (options.densityMap) {
450 26 patches /= MeanPatchAbundance(options.patches);
451
2/2
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 20 times.
 26 bedrock = options.bedrock.exposure ? std::fmax(1.0, a) / std::fmin(1.0, a) : 1.0;
452
1/1
✓ Branch 1 taken 26 times.
 26 map = BoulderMapAbundance(options).Max();
453 }
454 815 return patches * bedrock * MaxClusterAbundance(options.clusters) * map;
455 }
456
457 76 std::shared_ptr<const SizeFrequencyDistribution> EjectaSizes(const BoulderLayerOptions& options) {
458 76 const BoulderEjecta& ejecta = options.ejecta;
459
2/2
✓ Branch 1 taken 21 times.
✓ Branch 2 taken 55 times.
 76 if (ejecta.distribution) {
460 21 return ejecta.distribution;
461 }
462 55 return std::make_shared<GolombekRappSfd>(ejecta.coverage, options.sizes.qConstant,
463
1/1
✓ Branch 1 taken 55 times.
 55 options.sizes.qCoverageCoefficient);
464 }
465
466 std::shared_ptr<const SizeFrequencyDistribution>
467 306 BoulderCandidates(const BoulderLayerOptions& options) {
468
1/1
✓ Branch 1 taken 291 times.
 306 ValidateBoulderLayerOptions(options);
469 291 const BoulderSizes& sizes = options.sizes;
470
1/1
✓ Branch 1 taken 291 times.
 291 const double richest = MaxBackgroundAbundance(options);
471 291 std::vector<WeightedSfd> terms;
472
2/2
✓ Branch 1 taken 38 times.
✓ Branch 2 taken 253 times.
 291 if (sizes.distribution) {
473
1/1
✓ Branch 1 taken 38 times.
 38 terms.push_back({richest, sizes.distribution});
474
2/2
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 28 times.
 38 if (options.slopeMap) {
475 10 const SizeSlopeMap& map = *options.slopeMap;
476 10 const double d0 = map.referenceDiameter;
477
1/1
✓ Branch 2 taken 10 times.
 10 const double count = richest * sizes.distribution->Cumulative(d0);
478 10 const double low = map.minExponent;
479 10 const double high = map.maxExponent;
480
1/1
✓ Branch 1 taken 10 times.
 10 terms.push_back(
481 10 {1.0, std::make_shared<PowerLawSfd>(count * std::pow(d0, low) * high / low, low)});
482
1/1
✓ Branch 1 taken 10 times.
 10 terms.push_back({1.0, std::make_shared<PowerLawSfd>(count * std::pow(d0, high), high)});
483 }
484 } else {
485
1/1
✓ Branch 1 taken 253 times.
 253 terms.push_back(
486 253 {1.0, std::make_shared<GolombekRappSfd>(sizes.coverage * richest, sizes.qConstant,
487 253 sizes.qCoverageCoefficient)});
488 }
489
2/2
✓ Branch 1 taken 45 times.
✓ Branch 2 taken 246 times.
 291 if (options.ejecta.MinCraterDiameter() > 0.0) {
490
1/1
✓ Branch 1 taken 45 times.
 45 terms.push_back({1.0, EjectaSizes(options)});
491 }
492
1/1
✓ Branch 2 taken 291 times.
 582 return std::make_shared<SummedSfd>(std::move(terms));
493
4/4
✓ Branch 2 taken 10 times.
✓ Branch 6 taken 10 times.
✓ Branch 10 taken 253 times.
✓ Branch 14 taken 45 times.
 647 }
494
495 } // namespace terraingen
496