TerrainGen coverage

terraingen/
Coverage:
low: ≥ 0%
medium: ≥ 75.0%
high: ≥ 90.0%
Lines:
399 / 406
98.3%
Branches:
264 / 307
86.0%
Functions:
30 / 30
100.0%
layers/dunes.cpp
Line Branch Exec Source
1 #include "terraingen/layers/dunes.h"
2
3 #include "terraingen/noise/perlin.h"
4 #include "terraingen/random/stream.h"
5 #include "terraingen/util/checks.h"
6 #include "terraingen/util/surface.h"
7
8 #include <algorithm>
9 #include <cmath>
10 #include <numbers>
11 #include <stdexcept>
12
13 namespace terraingen {
14
15 namespace {
16
17 constexpr std::uint32_t kDuneStream = 0x44554e45u;
18 constexpr double kHalfPi = 0.5 * std::numbers::pi;
19 /// Wavelengths per lattice unit of the crest wander and of the defects along crests.
20 constexpr double kWanderWavelengths = 4.0;
21 constexpr double kDefectWavelengths = 3.0;
22 /// Phase by which barchanoid crest centres trail their links, in wavelengths, and the length of
23 /// a crest segment from link to link, in wavelengths.
24 constexpr double kBarchanoidBend = 0.15;
25 constexpr double kBarchanoidSegment = 1.5;
26 /// Newton steps for the closest point on a barchan's brink.
27 constexpr int kBrinkSteps = 6;
28 /// Share of the dune height over which sand ripples come in at a dune's edge.
29 constexpr double kSandEdge = 0.05;
30 /// Share of the ripple height kept on a slip face, where avalanches leave them faint.
31 constexpr double kSlipFaceRipples = 0.2;
32 /// Ripple crest wander and kinks: wavelengths per lattice unit of each, and the kinks' share of
33 /// the wander.
34 constexpr double kRippleWander = 5.0;
35 constexpr double kRippleKinkScale = 1.5;
36 constexpr double kRippleKinks = 0.25;
37 /// Share of the ripple height that fades along crests, and its scale in wavelengths.
38 constexpr double kRippleFading = 0.6;
39 constexpr double kRippleFadeScale = 4.0;
40 /// Forked crests: the smallest side in wavelengths of the cells their pairs are drawn in (cells
41 /// otherwise hold one pair on average), the separation of a pair in wavelengths, and the radius in
42 /// wavelengths within which its phase is whole and beyond which it is gone.
43 constexpr std::uint32_t kForkStream = 0x464f524bu;
44 constexpr double kForkMinCell = 0.5;
45 constexpr double kForkMin = 1.5;
46 constexpr double kForkMax = 3.0;
47 constexpr double kForkInner = 3.2;
48 constexpr double kForkOuter = 4.0;
49 /// Distance in wavelengths over which the ripples fade in from a fork's end.
50 constexpr double kForkCore = 0.5;
51 /// Crest breaks: the stream they are drawn from, and the distance in wavelengths over which a
52 /// crest falls to nothing at a break.
53 constexpr std::uint32_t kCrestBreakStream = 0x42524b53u;
54 constexpr double kBreakTaper = 1.0;
55 /// Grainflows: the share of a scarp's length its theatre head takes to reach full depth, the share
56 /// of a lobe's length it takes to rise, and the share after which its toe rounds off.
57 constexpr std::uint32_t kGrainflowStream = 0x47464c57u;
58 constexpr double kGrainflowHead = 0.15;
59 constexpr double kGrainflowRise = 0.25;
60 constexpr double kGrainflowToe = 0.6;
61
62 /// A cos^alpha bump falling from 1 at s = 0 to 0 at s = 1.
63 13605102 double Bump(double s, double alpha) {
64
2/2
✓ Branch 0 taken 6851610 times.
✓ Branch 1 taken 6753492 times.
 13605102 return s >= 1.0 ? 0.0 : std::pow(std::cos(kHalfPi * s), alpha);
65 }
66
67 /// The fractional part of `x`, in [0, 1).
68 11297102 double Fract(double x) {
69 11297102 return x - std::floor(x);
70 }
71
72 /// One isolated dune: its height, width and how far it is a barchan rather than a dome.
73 struct Dune {
74 double height = 0.0;
75 double width = 0.0;
76 double barchan = 1.0; ///< 1 for a barchan, 0 for a dome, between for the transition.
77 };
78
79 /// The height of a barchan at (x, y) in its wind frame: x downwind from the brink's apex, y across.
80 31740372 double BarchanAt(double x, double y, const Dune& dune, const DuneOptions& options) {
81 31740372 const double half = 0.5 * dune.width;
82
2/2
✓ Branch 0 taken 25615282 times.
✓ Branch 1 taken 6125090 times.
 31740372 if (std::fabs(y) >= half) {
83 25615282 return 0.0;
84 }
85 6125090 const double eta = y / half;
86 6125090 const double crest = dune.height * (1.0 - eta * eta);
87 6125090 const double bend = options.hornLength * dune.height / (half * half);
88 6125090 const double brink = bend * y * y;
89
2/2
✓ Branch 0 taken 5199478 times.
✓ Branch 1 taken 925612 times.
 6125090 if (x <= brink) {
90 5199478 const double stoss = dune.height / std::tan(options.stossSlope);
91 5199478 return crest * Bump((brink - x) / stoss, options.stossExponent);
92 }
93 925612 double t = y;
94
2/2
✓ Branch 0 taken 5553672 times.
✓ Branch 1 taken 925612 times.
 6479284 for (int step = 0; step < kBrinkSteps; ++step) {
95 5553672 const double f = 2.0 * bend * bend * t * t * t + (1.0 - 2.0 * bend * x) * t - y;
96 5553672 const double df = 6.0 * bend * bend * t * t + 1.0 - 2.0 * bend * x;
97
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5553672 times.
 5553672 if (df == 0.0) {
98 ✗ break;
99 }
100
1/1
✓ Branch 1 taken 5553672 times.
 5553672 t = std::clamp(t - f / df, -half, half);
101 }
102 925612 const double tau = t / half;
103 925612 const double distance = std::hypot(x - bend * t * t, y - t);
104 1851224 return std::max(0.0,
105 925612 dune.height * (1.0 - tau * tau) - std::tan(options.slipFaceAngle) * distance);
106 }
107
108 /// The height of a dome at (x, y) about its centre in its wind frame.
109 3733380 double DomeAt(double x, double y, const Dune& dune, const DuneOptions& options) {
110 3733380 const double along = x / (0.5 * options.domeLength * dune.width);
111 3733380 const double across = y / (0.5 * dune.width);
112 3733380 return dune.height * Bump(std::hypot(along, across), options.domeExponent);
113 }
114
115 /// The farthest a dune of `width` reaches from its centre (m).
116 31 double Reach(double width, const DuneOptions& options) {
117 31 const double height = options.heightRatio * width;
118 ✗ const double barchan = std::max(
119 31 {height / std::tan(options.stossSlope),
120
1/1
✓ Branch 1 taken 31 times.
 31 options.hornLength * height + height / std::tan(options.slipFaceAngle), 0.5 * width});
121 31 const double dome = 0.5 * width * std::max(1.0, options.domeLength);
122 62 return std::max(barchan, dome);
123 }
124
125 /// The ripple profile at phase `t` in [0, 1): a straight rise to a sharp crest, then a straight
126 /// fall over `lee`.
127 277936 double RippleProfile(double t, double lee) {
128 277936 const double upwind = 1.0 - lee;
129
2/2
✓ Branch 0 taken 231517 times.
✓ Branch 1 taken 46419 times.
 277936 return t < upwind ? t / upwind : 1.0 - (t - upwind) / lee;
130 }
131
132 } // namespace
133
134 61 void ValidateDuneOptions(const DuneOptions& options) {
135 61 const double right = kHalfPi;
136
4/4
✓ Branch 0 taken 60 times.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 59 times.
✓ Branch 3 taken 1 time.
 61 if (!(options.slipFaceAngle > 0.0 && options.slipFaceAngle < right &&
137
3/4
✓ Branch 0 taken 59 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 56 times.
✓ Branch 3 taken 3 times.
 59 options.stossSlope > 0.0 && options.stossSlope < options.slipFaceAngle &&
138
2/4
✓ Branch 0 taken 56 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 56 times.
 56 options.flankSlope > 0.0 && options.flankSlope < right)) {
139 throw std::invalid_argument(
140
1/1
✓ Branch 2 taken 5 times.
 5 "Dune slopes must be positive and below 90 deg, the stoss below the slip face");
141 }
142
4/6
✓ Branch 1 taken 56 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 56 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 54 times.
 112 if (!(IsPositive(options.stossExponent) && IsPositive(options.domeExponent) &&
143
4/4
✓ Branch 1 taken 55 times.
✓ Branch 2 taken 1 time.
✓ Branch 4 taken 54 times.
✓ Branch 5 taken 1 time.
 56 IsPositive(options.spacing) && IsPositive(options.width) &&
144
2/4
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 54 times.
✗ Branch 5 not taken.
 54 IsNonNegative(options.minWidth) && IsPositive(options.heightRatio) &&
145
2/4
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 54 times.
✗ Branch 5 not taken.
 54 IsNonNegative(options.hornLength) && IsPositive(options.domeLength) &&
146
2/4
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 54 times.
 54 IsPositive(options.wavelength) && IsNonNegative(options.height))) {
147
1/1
✓ Branch 2 taken 2 times.
 2 throw std::invalid_argument("Dune sizes and exponents must be positive");
148 }
149
5/6
✓ Branch 0 taken 54 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 53 times.
✓ Branch 3 taken 1 time.
✓ Branch 4 taken 52 times.
✓ Branch 5 taken 1 time.
 54 if (!(options.density >= 0.0 && options.density <= 1.0 && options.defects >= 0.0 &&
150
4/6
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 52 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 time.
✓ Branch 5 taken 51 times.
 52 options.defects <= 1.0 && options.widthSpread >= 0.0 && options.widthSpread < 1.0)) {
151 throw std::invalid_argument(
152
1/1
✓ Branch 2 taken 3 times.
 3 "Dune density and defects must be in [0, 1], and the width spread in [0, 1)");
153 }
154
5/6
✓ Branch 1 taken 50 times.
✓ Branch 2 taken 1 time.
✓ Branch 4 taken 50 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 49 times.
 101 if (!(IsNonNegative(options.sinuosity) && IsNonNegative(options.azimuthSpread) &&
155
3/4
✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 49 times.
✓ Branch 5 taken 1 time.
 50 std::isfinite(options.windAzimuth) && std::isfinite(options.anchor.x) &&
156
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 49 times.
 49 std::isfinite(options.anchor.y))) {
157 throw std::invalid_argument(
158
1/1
✓ Branch 2 taken 2 times.
 2 "Dune sinuosity and azimuth spread must be finite and not negative, the anchor finite");
159 }
160 49 const DuneRipples& ripples = options.ripples;
161
4/6
✓ Branch 1 taken 49 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 49 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 5 times.
✓ Branch 7 taken 44 times.
 98 if (!(IsNonNegative(ripples.wavelength) && IsNonNegative(ripples.height) &&
162
3/4
✓ Branch 0 taken 49 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 48 times.
✓ Branch 3 taken 1 time.
 49 ripples.leeFraction > 0.0 && ripples.leeFraction < 1.0 &&
163
2/4
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 48 times.
✗ Branch 5 not taken.
 48 IsNonNegative(ripples.sinuosity) && IsNonNegative(ripples.slopeCoupling) &&
164
3/4
✓ Branch 0 taken 47 times.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 47 times.
✗ Branch 3 not taken.
 48 ripples.interduneCover >= 0.0 && ripples.interduneCover <= 1.0 &&
165
4/4
✓ Branch 1 taken 46 times.
✓ Branch 2 taken 1 time.
✓ Branch 4 taken 45 times.
✓ Branch 5 taken 1 time.
 47 IsNonNegative(ripples.defectDensity) && IsNonNegative(options.defectDensity) &&
166
3/4
✓ Branch 0 taken 45 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
✓ Branch 3 taken 44 times.
 45 ripples.flankStraightening >= 0.0 && ripples.flankStraightening <= 1.0)) {
167 throw std::invalid_argument(
168 "Dune ripples need sizes and coupling not negative, a lee fraction in (0, 1) and an "
169 "interdune cover and flank straightening in [0, 1], and defect densities not "
170
1/1
✓ Branch 2 taken 5 times.
 5 "negative");
171 }
172
1/1
✓ Branch 1 taken 43 times.
 44 ValidateResolutionFade(options.fade);
173
1/1
✓ Branch 1 taken 43 times.
 43 ValidateResolutionFade(ripples.fade);
174 43 const DuneGrainflows& flows = options.grainflows;
175 161 const auto range = [](double low, double high) {
176
5/6
✓ Branch 1 taken 160 times.
✓ Branch 2 taken 1 time.
✓ Branch 3 taken 159 times.
✓ Branch 4 taken 1 time.
✓ Branch 6 taken 159 times.
✗ Branch 7 not taken.
 161 return IsPositive(low) && high >= low && std::isfinite(high);
177 };
178
5/6
✓ Branch 0 taken 43 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 41 times.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 5 times.
✓ Branch 5 taken 38 times.
 84 if (!(flows.cover >= 0.0 && flows.cover <= 1.0 &&
179
2/2
✓ Branch 1 taken 40 times.
✓ Branch 2 taken 1 time.
 41 range(flows.minScarpWidth, flows.maxScarpWidth) &&
180
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 range(flows.minScarpLength, flows.maxScarpLength) &&
181
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 range(flows.minLobeWidth, flows.maxLobeWidth) &&
182
4/4
✓ Branch 1 taken 39 times.
✓ Branch 2 taken 1 time.
✓ Branch 4 taken 1 time.
✓ Branch 5 taken 38 times.
 40 range(flows.minLobeLength, flows.maxLobeLength) && IsNonNegative(flows.thickness))) {
183 throw std::invalid_argument(
184 "Dune grainflows need a cover in [0, 1], positive sizes with each minimum at most its "
185
1/1
✓ Branch 2 taken 5 times.
 5 "maximum, and a thickness not negative");
186 }
187
1/1
✓ Branch 1 taken 37 times.
 38 ValidateResolutionFade(flows.fade);
188 37 }
189
190 27 DuneRipples MarsDuneRipples() {
191 27 DuneRipples ripples;
192 27 ripples.wavelength = 2.1;
193 27 ripples.height = 0.2;
194 27 const double leeSlope = std::tan(30.0 * kDegree);
195 27 ripples.leeFraction = ripples.height / (ripples.wavelength * leeSlope);
196 27 ripples.sinuosity = 0.8;
197 27 ripples.defectDensity = 0.099;
198 27 ripples.flankStraightening = 1.0;
199 27 return ripples;
200 }
201
202 12 DuneGrainflows MarsDuneGrainflows() {
203 12 DuneGrainflows flows;
204 12 flows.cover = 1.0;
205 12 flows.minScarpWidth = 0.54;
206 12 flows.maxScarpWidth = 0.67;
207 12 flows.minScarpLength = 1.0;
208 12 flows.maxScarpLength = 1.6;
209 12 flows.minLobeWidth = 0.57;
210 12 flows.maxLobeWidth = 0.66;
211 12 flows.minLobeLength = 2.9;
212 12 flows.maxLobeLength = 6.0;
213 12 flows.thickness = 0.015;
214 12 return flows;
215 }
216
217 8 DuneOptions MarsBarchans() {
218 8 DuneOptions dunes;
219 8 dunes.shape = DuneShape::Barchan;
220 8 dunes.slipFaceAngle = 29.0 * kDegree;
221 8 dunes.stossSlope = 8.0 * kDegree;
222 8 dunes.width = 370.0;
223 8 dunes.minWidth = 200.0;
224 8 dunes.domeLength = 2.0;
225 8 dunes.spacing = 900.0;
226 8 dunes.density = 0.5;
227 8 dunes.ripples = MarsDuneRipples();
228 8 return dunes;
229 }
230
231 10 DuneOptions MarsTransverseDunes() {
232 10 DuneOptions dunes;
233 10 dunes.shape = DuneShape::Transverse;
234 10 dunes.slipFaceAngle = 29.0 * kDegree;
235 10 dunes.stossSlope = 8.0 * kDegree;
236 10 dunes.wavelength = 248.0;
237 10 dunes.height = 28.0;
238 10 dunes.defectDensity = 0.0027;
239 10 dunes.ripples = MarsDuneRipples();
240 10 return dunes;
241 }
242
243 53 DuneLayer::DuneLayer(DuneOptions options, Surface surface)
244 53 : FormulaLayer(surface), options_(options) {
245
1/1
✓ Branch 1 taken 31 times.
 53 ValidateDuneOptions(options_);
246 31 wind_ = {std::cos(options_.windAzimuth), std::sin(options_.windAzimuth)};
247 31 const double widest = options_.width * (1.0 + options_.widthSpread);
248
1/1
✓ Branch 1 taken 31 times.
 31 rings_ = static_cast<std::int64_t>(std::ceil(Reach(widest, options_) / options_.spacing));
249
2/2
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 30 times.
 31 if (surface_.IsBody()) {
250
1/1
✓ Branch 1 taken 1 time.
 1 anchorDirection_ = UnitVector(options_.anchor.x, options_.anchor.y);
251 }
252 75 }
253
254 9992116 bool DuneLayer::Planar(double x, double y, Vec2d& planar) const {
255
2/2
✓ Branch 1 taken 9992113 times.
✓ Branch 2 taken 3 times.
 9992116 if (!surface_.IsBody()) {
256 9992113 planar = {x, y};
257 9992113 return true;
258 }
259
3/3
✓ Branch 1 taken 3 times.
✓ Branch 4 taken 1 time.
✓ Branch 5 taken 2 times.
 3 if (Dot(UnitVector(x, y), anchorDirection_) <= 0.0) {
260 1 return false;
261 }
262
1/1
✓ Branch 1 taken 2 times.
 2 planar = surface_.Offset(options_.anchor, {x, y});
263 2 return true;
264 }
265
266 3941528 double DuneLayer::Barchans(const Vec2d& q) const {
267 3941528 const auto ix = static_cast<std::int64_t>(std::floor(q.x / options_.spacing));
268 3941528 const auto iy = static_cast<std::int64_t>(std::floor(q.y / options_.spacing));
269 3941528 double highest = 0.0;
270
2/2
✓ Branch 0 taken 11824584 times.
✓ Branch 1 taken 3941528 times.
 15766112 for (std::int64_t gy = iy - rings_; gy <= iy + rings_; ++gy) {
271
2/2
✓ Branch 0 taken 35473752 times.
✓ Branch 1 taken 11824584 times.
 47298336 for (std::int64_t gx = ix - rings_; gx <= ix + rings_; ++gx) {
272 35473752 RandomStream random(options_.seed, kDuneStream, static_cast<std::uint32_t>(gx),
273 static_cast<std::uint32_t>(gy),
274
1/1
✓ Branch 1 taken 35473752 times.
 35473752 static_cast<std::uint32_t>((gx >> 32) ^ ((gy >> 32) << 16)));
275
2/3
✓ Branch 1 taken 35473752 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 35473752 times.
 35473752 if (random.NextDouble() >= options_.density) {
276 ✗ continue;
277 }
278
1/1
✓ Branch 1 taken 35473752 times.
 35473752 const double cx = (static_cast<double>(gx) + random.NextDouble()) * options_.spacing;
279
1/1
✓ Branch 1 taken 35473752 times.
 35473752 const double cy = (static_cast<double>(gy) + random.NextDouble()) * options_.spacing;
280 35473752 Dune dune;
281 35473752 dune.width =
282
1/1
✓ Branch 1 taken 35473752 times.
 35473752 options_.width * (1.0 + options_.widthSpread * (2.0 * random.NextDouble() - 1.0));
283 35473752 dune.height = options_.heightRatio * dune.width;
284 const double azimuth =
285
1/1
✓ Branch 1 taken 35473752 times.
 35473752 options_.windAzimuth + options_.azimuthSpread * (2.0 * random.NextDouble() - 1.0);
286 70947504 dune.barchan = options_.shape == DuneShape::Dome
287
2/2
✓ Branch 0 taken 32811075 times.
✓ Branch 1 taken 2662677 times.
 35473752 ? 0.0
288 32811075 : Smoothstep((dune.width - options_.minWidth) /
289
1/1
✓ Branch 2 taken 32811075 times.
 32811075 (0.5 * std::max(options_.minWidth, 1e-9)));
290 35473752 const double dx = q.x - cx;
291 35473752 const double dy = q.y - cy;
292 35473752 const double c = std::cos(azimuth);
293 35473752 const double s = std::sin(azimuth);
294 35473752 const double x = dx * c + dy * s;
295 35473752 const double y = dy * c - dx * s;
296 35473752 double h = 0.0;
297
2/2
✓ Branch 0 taken 31740372 times.
✓ Branch 1 taken 3733380 times.
 35473752 if (dune.barchan > 0.0) {
298
1/1
✓ Branch 1 taken 31740372 times.
 31740372 h += dune.barchan * BarchanAt(x, y, dune, options_);
299 }
300
2/2
✓ Branch 0 taken 3733380 times.
✓ Branch 1 taken 31740372 times.
 35473752 if (dune.barchan < 1.0) {
301 3733380 h += (1.0 - dune.barchan) * DomeAt(x, y, dune, options_);
302 }
303 35473752 highest = std::max(highest, h);
304 }
305 }
306 3941528 return highest;
307 }
308
309 11019166 double DuneLayer::Crests(const Vec2d& q) const {
310 11019166 const double lambda = options_.wavelength;
311 11019166 const Vec2d across{-wind_.y, wind_.x};
312 11019166 const bool linear = options_.shape == DuneShape::Linear;
313
2/2
✓ Branch 0 taken 1445906 times.
✓ Branch 1 taken 9573260 times.
 11019166 const Vec2d normal = linear ? across : wind_;
314 8 const auto segment = [&](const Vec2d& p) {
315 8 return std::fabs(std::cos(std::numbers::pi * (p.x * across.x + p.y * across.y) /
316 8 (kBarchanoidSegment * lambda)));
317 11019166 };
318 33057498 const auto phase = [&](const Vec2d& p) {
319 double value =
320 33057498 (p.x * normal.x + p.y * normal.y) / lambda +
321 66114996 options_.sinuosity * static_cast<double>(PerlinNoise2Precise(
322 33057498 p.x / (kWanderWavelengths * lambda),
323 33057498 p.y / (kWanderWavelengths * lambda), options_.seed));
324
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 33057492 times.
 33057498 if (options_.shape == DuneShape::Barchanoid) {
325 6 value += kBarchanoidBend * segment(p);
326 }
327 33057498 return value;
328 11019166 };
329 11019166 const double psi = phase(q);
330 11019166 const double step = 0.01 * lambda;
331 11019166 const double gx = (phase({q.x + step, q.y}) - psi) / step;
332 11019166 const double gy = (phase({q.x, q.y + step}) - psi) / step;
333 11019166 const double stretch = std::max(1e-6, std::hypot(gx, gy) * lambda);
334 double amplitude =
335 22038332 1.0 - options_.defects *
336 22038332 (0.5 + 0.5 * static_cast<double>(PerlinNoise2Precise(
337 11019166 q.x / (kDefectWavelengths * lambda),
338 11019166 q.y / (kDefectWavelengths * lambda), options_.seed + 1u)));
339
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 11019164 times.
 11019166 if (options_.shape == DuneShape::Barchanoid) {
340 2 amplitude *= 0.6 + 0.4 * segment(q);
341 }
342
2/2
✓ Branch 0 taken 1512000 times.
✓ Branch 1 taken 9507166 times.
 11019166 if (options_.defectDensity > 0.0) {
343
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1512000 times.
 1512000 const Vec2d along = linear ? wind_ : across;
344 1512000 amplitude *=
345
1/1
✓ Branch 1 taken 1512000 times.
 1512000 CrestBreaks(static_cast<std::int64_t>(std::floor(psi)), q.x * along.x + q.y * along.y);
346 }
347 11019166 const double height = amplitude * options_.height;
348
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 11019166 times.
 11019166 if (height <= 0.0) {
349 ✗ return 0.0;
350 }
351 11019166 const double u = Fract(psi);
352
2/2
✓ Branch 0 taken 1445906 times.
✓ Branch 1 taken 9573260 times.
 11019166 if (linear) {
353 const double flank =
354 1445906 std::min(0.5, height * stretch / (std::tan(options_.flankSlope) * lambda));
355 1445906 const double offset = std::fabs(u - 0.5);
356
2/2
✓ Branch 0 taken 1059845 times.
✓ Branch 1 taken 386061 times.
 1445906 return offset >= flank ? 0.0 : height * Bump(offset / flank, 2.0);
357 }
358 9573260 const double leeSlope = std::tan(options_.slipFaceAngle);
359 9573260 const double stossSlope = std::tan(options_.stossSlope);
360 9573260 const double tallest = options_.height * stretch / lambda;
361 9573260 const double brink = std::min(1.0 - std::min(1.0, tallest / leeSlope), tallest / stossSlope);
362 9573260 const double stoss = std::min(brink, height * stretch / (stossSlope * lambda));
363
2/2
✓ Branch 0 taken 4286183 times.
✓ Branch 1 taken 5287077 times.
 9573260 if (u <= brink) {
364
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4286183 times.
 4286183 return stoss <= 0.0 ? 0.0 : height * Bump((brink - u) / stoss, options_.stossExponent);
365 }
366 5287077 const double down = (u - brink) * lambda / stretch;
367 5287077 return std::max(0.0, height - leeSlope * down);
368 }
369
370 14960694 double DuneLayer::Height(const Vec2d& q) const {
371
2/3
✓ Branch 0 taken 3941528 times.
✓ Branch 1 taken 11019166 times.
✗ Branch 2 not taken.
 14960694 switch (options_.shape) {
372 3941528 case DuneShape::Barchan:
373 case DuneShape::Dome:
374 3941528 return Barchans(q);
375 11019166 case DuneShape::Transverse:
376 case DuneShape::Barchanoid:
377 case DuneShape::Linear:
378 11019166 return Crests(q);
379 }
380 ✗ return 0.0;
381 }
382
383 1512000 double DuneLayer::CrestBreaks(std::int64_t crest, double along) const {
384 1512000 const double taper = kBreakTaper * options_.wavelength;
385 1512000 const double length = 1.0 / options_.defectDensity + taper;
386 1512000 const auto cell = static_cast<std::int64_t>(std::floor(along / length));
387 1512000 double weight = 1.0;
388
2/2
✓ Branch 0 taken 4536000 times.
✓ Branch 1 taken 1512000 times.
 6048000 for (std::int64_t c = cell - 1; c <= cell + 1; ++c) {
389 4536000 RandomStream random(options_.seed, kCrestBreakStream, static_cast<std::uint32_t>(crest),
390 static_cast<std::uint32_t>(c),
391
1/1
✓ Branch 1 taken 4536000 times.
 4536000 static_cast<std::uint32_t>((crest >> 32) ^ ((c >> 32) << 16)));
392
1/1
✓ Branch 1 taken 4536000 times.
 4536000 const double at = (static_cast<double>(c) + random.NextDouble()) * length;
393
1/1
✓ Branch 1 taken 4536000 times.
 4536000 weight = std::min(weight, Smoothstep(std::fabs(along - at) / taper));
394 }
395 1512000 return weight;
396 }
397
398 1242145 Vec2d DuneLayer::Gradient(const Vec2d& q, double step) const {
399
2/2
✓ Branch 1 taken 1242145 times.
✓ Branch 4 taken 1242145 times.
 1242145 return {(Height({q.x + step, q.y}) - Height({q.x - step, q.y})) / (2.0 * step),
400
2/2
✓ Branch 1 taken 1242145 times.
✓ Branch 4 taken 1242145 times.
 1242145 (Height({q.x, q.y + step}) - Height({q.x, q.y - step})) / (2.0 * step)};
401 }
402
403 964209 double DuneLayer::SlipFaceWeight(const Vec2d& g) const {
404 964209 const double grade = std::hypot(g.x, g.y);
405
2/2
✓ Branch 0 taken 294672 times.
✓ Branch 1 taken 669537 times.
 964209 if (grade <= 0.0) {
406 294672 return 0.0;
407 }
408 669537 const double facing = -(g.x * wind_.x + g.y * wind_.y) / grade;
409 669537 return Smoothstep(2.0 * facing) *
410 669537 Smoothstep(2.0 * grade / std::tan(options_.slipFaceAngle) - 1.0);
411 }
412
413 964209 double DuneLayer::Grainflows(const Vec2d& q, double face) const {
414 964209 const DuneGrainflows& flows = options_.grainflows;
415
2/2
✓ Branch 0 taken 851609 times.
✓ Branch 1 taken 112600 times.
 964209 if (face <= 0.0) {
416 851609 return 0.0;
417 }
418 112600 const Vec2d across{-wind_.y, wind_.x};
419 112600 const double u = q.x * wind_.x + q.y * wind_.y;
420 112600 const double v = q.x * across.x + q.y * across.y;
421 112600 const double cellAlong = flows.maxScarpLength + flows.maxLobeLength;
422 112600 const double cellAcross = std::max(flows.maxScarpWidth, flows.maxLobeWidth);
423 112600 const auto iu = static_cast<std::int64_t>(std::floor(u / cellAlong));
424 112600 const auto iv = static_cast<std::int64_t>(std::floor(v / cellAcross));
425 3896748 const auto between = [](double low, double high, double t) { return low + (high - low) * t; };
426 112600 double relief = 0.0;
427
2/2
✓ Branch 0 taken 337800 times.
✓ Branch 1 taken 112600 times.
 450400 for (std::int64_t cu = iu - 1; cu <= iu + 1; ++cu) {
428
2/2
✓ Branch 0 taken 1013400 times.
✓ Branch 1 taken 337800 times.
 1351200 for (std::int64_t cv = iv - 1; cv <= iv + 1; ++cv) {
429 1013400 RandomStream random(options_.seed, kGrainflowStream, static_cast<std::uint32_t>(cu),
430 static_cast<std::uint32_t>(cv),
431
1/1
✓ Branch 1 taken 1013400 times.
 1013400 static_cast<std::uint32_t>((cu >> 32) ^ ((cv >> 32) << 16)));
432
3/3
✓ Branch 1 taken 1013400 times.
✓ Branch 3 taken 39213 times.
✓ Branch 4 taken 974187 times.
 1013400 if (random.NextDouble() >= flows.cover) {
433 39213 continue;
434 }
435
1/1
✓ Branch 1 taken 974187 times.
 974187 const double head = (static_cast<double>(cu) + random.NextDouble()) * cellAlong;
436
1/1
✓ Branch 1 taken 974187 times.
 974187 const double middle = (static_cast<double>(cv) + random.NextDouble()) * cellAcross;
437 const double scarpWidth =
438
1/1
✓ Branch 1 taken 974187 times.
 974187 between(flows.minScarpWidth, flows.maxScarpWidth, random.NextDouble());
439 const double scarpLength =
440
1/1
✓ Branch 1 taken 974187 times.
 974187 between(flows.minScarpLength, flows.maxScarpLength, random.NextDouble());
441 const double lobeWidth =
442
1/1
✓ Branch 1 taken 974187 times.
 974187 between(flows.minLobeWidth, flows.maxLobeWidth, random.NextDouble());
443 const double lobeLength =
444
1/1
✓ Branch 1 taken 974187 times.
 974187 between(flows.minLobeLength, flows.maxLobeLength, random.NextDouble());
445 974187 const double s = u - head;
446 974187 const double n = std::fabs(v - middle);
447
4/4
✓ Branch 0 taken 488587 times.
✓ Branch 1 taken 485600 times.
✓ Branch 2 taken 57004 times.
✓ Branch 3 taken 431583 times.
 974187 if (s >= 0.0 && s < scarpLength) {
448 57004 const double half = 0.5 * scarpWidth;
449 const double profile =
450
1/1
✓ Branch 1 taken 57004 times.
 57004 Smoothstep(s / (kGrainflowHead * scarpLength)) * (1.0 - s / scarpLength);
451 57004 relief -= profile * std::max(0.0, 1.0 - (n / half) * (n / half));
452 }
453 974187 const double t = s - scarpLength;
454
4/4
✓ Branch 0 taken 431583 times.
✓ Branch 1 taken 542604 times.
✓ Branch 2 taken 196362 times.
✓ Branch 3 taken 235221 times.
 974187 if (t >= 0.0 && t < lobeLength) {
455 392724 const double toe = std::max(0.0, (t - kGrainflowToe * lobeLength) /
456 196362 ((1.0 - kGrainflowToe) * lobeLength));
457 196362 const double half = 0.5 * lobeWidth * std::sqrt(std::max(0.0, 1.0 - toe * toe));
458
2/2
✓ Branch 0 taken 54631 times.
✓ Branch 1 taken 141731 times.
 196362 if (n < half) {
459
1/1
✓ Branch 1 taken 54631 times.
 54631 const double profile = Smoothstep(t / (kGrainflowRise * lobeLength));
460 54631 relief += profile * (1.0 - (n / half) * (n / half));
461 }
462 }
463 }
464 }
465 112600 return face * flows.thickness * relief;
466 }
467
468 277936 double DuneLayer::Ripples(const Vec2d& q, double dune) const {
469 277936 const DuneRipples& ripples = options_.ripples;
470 277936 const double lambda = ripples.wavelength;
471
1/1
✓ Branch 1 taken 277936 times.
 277936 const Vec2d g = Gradient(q, 0.25 * lambda);
472 277936 const double gx = g.x;
473 277936 const double gy = g.y;
474 277936 const double kx = wind_.x - ripples.slopeCoupling * gx;
475 277936 const double ky = wind_.y - ripples.slopeCoupling * gy;
476 277936 const double tighten = std::max(1e-6, std::hypot(kx, ky));
477 277936 double sinuosity = ripples.sinuosity;
478
2/2
✓ Branch 0 taken 276751 times.
✓ Branch 1 taken 1185 times.
 277936 if (ripples.flankStraightening > 0.0) {
479 276751 const double grade = std::hypot(gx, gy);
480 const double stoss =
481
2/2
✓ Branch 0 taken 276184 times.
✓ Branch 1 taken 567 times.
 276751 grade > 0.0 ? std::max(0.0, (gx * wind_.x + gy * wind_.y) / grade) : 1.0;
482 553502 sinuosity *= 1.0 - ripples.flankStraightening *
483
1/1
✓ Branch 1 taken 276751 times.
 276751 Smoothstep(grade / std::tan(options_.stossSlope)) * (1.0 - stoss);
484 }
485 277936 double core = 1.0;
486 const double psi =
487 555872 (q.x * wind_.x + q.y * wind_.y - ripples.slopeCoupling * dune) / lambda +
488 833808 sinuosity * (static_cast<double>(PerlinNoise2Precise(q.x / (kRippleWander * lambda),
489 277936 q.y / (kRippleWander * lambda),
490 277936 options_.seed + 2u)) +
491 555872 kRippleKinks * static_cast<double>(PerlinNoise2Precise(
492 277936 q.x / (kRippleKinkScale * lambda),
493 277936 q.y / (kRippleKinkScale * lambda), options_.seed + 3u))) +
494
1/1
✓ Branch 1 taken 277936 times.
 277936 Forks(q, lambda, ripples.defectDensity, core);
495 const double fading =
496 1.0 -
497 555872 kRippleFading * (0.5 + 0.5 * static_cast<double>(PerlinNoise2Precise(
498 277936 q.x / (kRippleFadeScale * lambda),
499 277936 q.y / (kRippleFadeScale * lambda), options_.seed + 4u)));
500
1/1
✓ Branch 1 taken 277936 times.
 277936 const double sand = SandOf(dune);
501 277936 const double steep = std::hypot(gx, gy) / std::tan(options_.slipFaceAngle);
502
1/1
✓ Branch 1 taken 277936 times.
 277936 const double avalanche = 1.0 - (1.0 - kSlipFaceRipples) * Smoothstep(2.0 * steep - 1.0);
503 277936 const double share =
504 277936 (ripples.interduneCover + (1.0 - ripples.interduneCover) * sand) * avalanche;
505 555872 return share * fading * core * ripples.height / std::max(1.0, tighten) *
506 555872 RippleProfile(Fract(psi), ripples.leeFraction);
507 }
508
509 277936 double DuneLayer::Forks(const Vec2d& q, double wavelength, double density, double& core) const {
510
2/2
✓ Branch 0 taken 2370 times.
✓ Branch 1 taken 275566 times.
 277936 if (density <= 0.0) {
511 2370 return 0.0;
512 }
513 275566 const double cell = std::max(std::sqrt(wavelength / density), kForkMinCell * wavelength);
514 275566 const double none = std::exp(-density * cell * cell / wavelength);
515 275566 const auto rings = static_cast<std::int64_t>(std::ceil(kForkOuter * wavelength / cell));
516 275566 const double inner = kForkInner * wavelength;
517 275566 const double outer = kForkOuter * wavelength;
518 275566 const auto ix = static_cast<std::int64_t>(std::floor(q.x / cell));
519 275566 const auto iy = static_cast<std::int64_t>(std::floor(q.y / cell));
520 275566 double phase = 0.0;
521
2/2
✓ Branch 0 taken 1389936 times.
✓ Branch 1 taken 275566 times.
 1665502 for (std::int64_t gy = iy - rings; gy <= iy + rings; ++gy) {
522
2/2
✓ Branch 0 taken 7034422 times.
✓ Branch 1 taken 1389936 times.
 8424358 for (std::int64_t gx = ix - rings; gx <= ix + rings; ++gx) {
523 7034422 RandomStream random(options_.seed, kForkStream, static_cast<std::uint32_t>(gx),
524 static_cast<std::uint32_t>(gy),
525
1/1
✓ Branch 1 taken 7034422 times.
 7034422 static_cast<std::uint32_t>((gx >> 32) ^ ((gy >> 32) << 16)));
526
1/1
✓ Branch 1 taken 7034422 times.
 7034422 double draw = random.NextDouble();
527
2/2
✓ Branch 0 taken 6824610 times.
✓ Branch 1 taken 7034422 times.
 13859032 while (draw > none) {
528
1/1
✓ Branch 1 taken 6824610 times.
 6824610 draw *= random.NextDouble();
529
1/1
✓ Branch 1 taken 6824610 times.
 6824610 const Vec2d centre{(static_cast<double>(gx) + random.NextDouble()) * cell,
530
1/1
✓ Branch 1 taken 6824610 times.
 6824610 (static_cast<double>(gy) + random.NextDouble()) * cell};
531
1/1
✓ Branch 1 taken 6824610 times.
 6824610 const double turn = 2.0 * std::numbers::pi * random.NextDouble();
532 const double half =
533
1/1
✓ Branch 1 taken 6824610 times.
 6824610 0.5 * wavelength * (kForkMin + (kForkMax - kForkMin) * random.NextDouble());
534
3/3
✓ Branch 1 taken 6824610 times.
✓ Branch 3 taken 3362440 times.
✓ Branch 4 taken 3462170 times.
 6824610 const double sign = random.NextDouble() < 0.5 ? -1.0 : 1.0;
535 6824610 const double r = std::hypot(q.x - centre.x, q.y - centre.y);
536
1/1
✓ Branch 1 taken 6824610 times.
 6824610 const double window = 1.0 - Smoothstep((r - inner) / (outer - inner));
537
2/2
✓ Branch 0 taken 3979005 times.
✓ Branch 1 taken 2845605 times.
 6824610 if (window <= 0.0) {
538 3979005 continue;
539 }
540 2845605 const Vec2d a{centre.x - half * std::cos(turn), centre.y - half * std::sin(turn)};
541 2845605 const Vec2d b{centre.x + half * std::cos(turn), centre.y + half * std::sin(turn)};
542 double subtended =
543 2845605 std::atan2(q.y - a.y, q.x - a.x) - std::atan2(q.y - b.y, q.x - b.x);
544 2845605 subtended -= 2.0 * std::numbers::pi *
545 2845605 std::floor((subtended + std::numbers::pi) / (2.0 * std::numbers::pi));
546 2845605 phase += sign * window * subtended / (2.0 * std::numbers::pi);
547 const double nearest =
548 2845605 std::min(std::hypot(q.x - a.x, q.y - a.y), std::hypot(q.x - b.x, q.y - b.y));
549
1/1
✓ Branch 1 taken 2845605 times.
 2845605 core = std::min(core, Smoothstep(nearest / (kForkCore * wavelength)));
550 }
551 }
552 }
553 275566 return phase;
554 }
555
556 279755 double DuneLayer::SandOf(double dune) const {
557 279755 const double reference =
558
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4389 times.
 4389 options_.shape == DuneShape::Barchan || options_.shape == DuneShape::Dome
559
2/2
✓ Branch 0 taken 4389 times.
✓ Branch 1 taken 275366 times.
 284144 ? options_.heightRatio * options_.width
560 : options_.height;
561
1/1
✓ Branch 2 taken 279755 times.
 279755 return Smoothstep(dune / std::max(1e-9, kSandEdge * reference));
562 }
563
564 1819 double DuneLayer::MaterialCoverAt(double x, double y, double /*spacing*/,
565 const Ground* /*ground*/) const {
566 1819 Vec2d q;
567
2/3
✓ Branch 1 taken 1819 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1819 times.
 1819 if (!Planar(x, y, q)) {
568 ✗ return 0.0;
569 }
570
2/2
✓ Branch 1 taken 1819 times.
✓ Branch 4 taken 1819 times.
 1819 const double sand = SandOf(Height(q));
571 1819 const DuneRipples& ripples = options_.ripples;
572
3/4
✓ Branch 0 taken 405 times.
✓ Branch 1 taken 1414 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 405 times.
 1819 if (ripples.wavelength <= 0.0 || ripples.height <= 0.0) {
573 1414 return sand;
574 }
575 405 return ripples.interduneCover + (1.0 - ripples.interduneCover) * sand;
576 }
577
578 7254345 double DuneLayer::DuneHeightAt(double x, double y) const {
579 7254345 Vec2d q;
580
4/4
✓ Branch 1 taken 7254345 times.
✓ Branch 3 taken 7254344 times.
✓ Branch 4 taken 1 time.
✓ Branch 6 taken 7254344 times.
 14508690 return Planar(x, y, q) ? Height(q) : 0.0;
581 }
582
583 2735952 float DuneLayer::HeightAt(float height, double x, double y, double spacing,
584 const Ground* /*ground*/) const {
585 2735952 Vec2d q;
586
2/3
✓ Branch 1 taken 2735952 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2735952 times.
 2735952 if (!Planar(x, y, q)) {
587 ✗ return height;
588 }
589
3/4
✓ Branch 0 taken 1399960 times.
✓ Branch 1 taken 1335992 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1399960 times.
 2735952 const bool isolated = options_.shape == DuneShape::Barchan || options_.shape == DuneShape::Dome;
590
2/2
✓ Branch 0 taken 1335992 times.
✓ Branch 1 taken 1399960 times.
 2735952 const double size = isolated ? options_.width : options_.wavelength;
591
1/1
✓ Branch 1 taken 2735952 times.
 2735952 const double duneWeight = ResolutionWeight(size, spacing, options_.fade);
592 2735952 const DuneRipples& ripples = options_.ripples;
593
1/2
✓ Branch 0 taken 1328562 times.
✗ Branch 1 not taken.
 1328562 const double rippleWeight = ripples.wavelength > 0.0 && ripples.height > 0.0
594
3/3
✓ Branch 0 taken 1328562 times.
✓ Branch 1 taken 1407390 times.
✓ Branch 3 taken 1328562 times.
 4064514 ? ResolutionWeight(ripples.wavelength, spacing, ripples.fade)
595 2735952 : 0.0;
596 2735952 const DuneGrainflows& flows = options_.grainflows;
597
1/2
✓ Branch 0 taken 2014889 times.
✗ Branch 1 not taken.
 2014889 const double flowWeight = flows.cover > 0.0 && flows.thickness > 0.0
598
3/3
✓ Branch 0 taken 2014889 times.
✓ Branch 1 taken 721063 times.
✓ Branch 3 taken 2014889 times.
 4750841 ? ResolutionWeight(flows.minLobeWidth, spacing, flows.fade)
599 2735952 : 0.0;
600
4/6
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 2735951 times.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 time.
✗ Branch 5 not taken.
 2735952 if (duneWeight <= 0.0 && rippleWeight <= 0.0 && flowWeight <= 0.0) {
601 1 return height;
602 }
603
1/1
✓ Branch 1 taken 2735951 times.
 2735951 const double dune = Height(q);
604 2735951 double relief = duneWeight * dune;
605
2/2
✓ Branch 0 taken 277936 times.
✓ Branch 1 taken 2458015 times.
 2735951 if (rippleWeight > 0.0) {
606
1/1
✓ Branch 1 taken 277936 times.
 277936 relief += rippleWeight * Ripples(q, dune);
607 }
608
2/2
✓ Branch 0 taken 964209 times.
✓ Branch 1 taken 1771742 times.
 2735951 if (flowWeight > 0.0) {
609
1/1
✓ Branch 1 taken 964209 times.
 964209 const Vec2d g = Gradient(q, 0.25 * flows.minScarpLength);
610
2/2
✓ Branch 1 taken 964209 times.
✓ Branch 4 taken 964209 times.
 964209 relief += flowWeight * Grainflows(q, SlipFaceWeight(g));
611 }
612 2735951 return static_cast<float>(static_cast<double>(height) + relief);
613 }
614
615 } // namespace terraingen
616