TerrainGen coverage

terraingen/
Coverage:
low: ≥ 0%
medium: ≥ 75.0%
high: ≥ 90.0%
Lines:
154 / 157
98.1%
Branches:
111 / 122
91.0%
Functions:
22 / 22
100.0%
layers/craters/secondary_model.cpp
Line Branch Exec Source
1 #include "terraingen/layers/craters/secondary_model.h"
2
3 #include "terraingen/layers/fade.h"
4 #include "terraingen/noise/value.h"
5 #include "terraingen/util/checks.h"
6
7 #include <algorithm>
8 #include <cmath>
9 #include <limits>
10 #include <numbers>
11 #include <stdexcept>
12
13 namespace terraingen {
14
15 namespace {
16
17 /// Singer et al. (2020, Table 3): d = a R^-b, a = kSingerScale D^kSingerSizeExponent, b =
18 /// kSingerRangeExponent ln D, in km.
19 constexpr double kSingerScale = 0.022;
20 constexpr double kSingerSizeExponent = 1.96;
21 constexpr double kSingerRangeExponent = 0.14;
22
23 /// The outer range of rays as a multiple of their median length (Elliott et al. 2018).
24 constexpr double kLongestRayOverMedian = 1.6;
25
26 /// How strongly secondaries clump at `u` clump wavelengths along a ray keyed `key`: 0 between
27 /// clumps, 1 in them.
28 30677469 double ClumpBump(double u, std::uint32_t key) {
29 30677469 const double n = static_cast<double>(ValueNoise2(static_cast<float>(u), 0.5f, key));
30
2/2
✓ Branch 1 taken 30677469 times.
✓ Branch 4 taken 30677469 times.
 30677469 return Smoothstep(std::clamp(0.5 + 1.5 * n, 0.0, 1.0));
31 }
32
33 /// The mean of `ClumpBump` along a ray, found once.
34 1030 double MeanClumpBump() {
35 4 static const double mean = [] {
36 4 constexpr int kSamples = 1 << 16;
37 4 double sum = 0.0;
38
2/2
✓ Branch 0 taken 262144 times.
✓ Branch 1 taken 4 times.
 262148 for (int i = 0; i < kSamples; ++i) {
39 262144 sum += ClumpBump(0.0371 * i, 0x434c554du);
40 }
41 4 return sum / kSamples;
42
4/7
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1026 times.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 1030 }();
43 1030 return mean;
44 }
45
46 } // namespace
47
48 3111001 double SingerMaxSecondaryDiameter(double diameter, double range) {
49 3111001 const double d = diameter / 1000.0;
50 3111001 const double r = range / 1000.0;
51 3111001 return 1000.0 * kSingerScale * std::pow(d, kSingerSizeExponent) *
52 3111001 std::pow(r, -kSingerRangeExponent * std::log(d));
53 }
54
55 3225 double SingerRangeForSize(double diameter, double size) {
56 3225 const double d = diameter / 1000.0;
57 3225 const double logD = std::log(d);
58
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 3223 times.
 3225 if (!(logD > 0.0)) {
59 2 return std::numeric_limits<double>::infinity();
60 }
61 3223 const double a = kSingerScale * std::pow(d, kSingerSizeExponent);
62 3223 return 1000.0 * std::pow(a / (size / 1000.0), 1.0 / (kSingerRangeExponent * logD));
63 }
64
65 482165 double SecondaryOnsetRange(double diameter) {
66 482165 const double small = 3000.0;
67 482165 const double large = 31000.0;
68 482165 double rim = 1.8;
69
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 482160 times.
 482165 if (diameter >= large) {
70 5 rim = 1.3;
71
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 482154 times.
 482160 } else if (diameter > small) {
72 6 const double t = std::log(diameter / small) / std::log(large / small);
73 6 rim = 1.8 + t * (1.3 - 1.8);
74 }
75 482165 return (0.5 + rim) * diameter;
76 }
77
78 952099 double ElliottMedianRayLength(double diameter, double coefficient, double exponent) {
79 952099 return 1000.0 * coefficient * std::pow(0.0005 * diameter, exponent);
80 }
81
82 74 void ValidateSecondaryCraterOptions(const SecondaryCraterOptions& options) {
83
6/6
✓ Branch 1 taken 73 times.
✓ Branch 2 taken 1 time.
✓ Branch 5 taken 1 time.
✓ Branch 6 taken 72 times.
✓ Branch 7 taken 2 times.
✓ Branch 8 taken 72 times.
 74 if (!(options.onsetCount && IsPositive(*options.onsetCount))) {
84 throw std::invalid_argument("Secondary craters need an onset count: the sources give their "
85
1/1
✓ Branch 2 taken 2 times.
 2 "sizes, not their number");
86 }
87
6/6
✓ Branch 1 taken 71 times.
✓ Branch 2 taken 1 time.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 69 times.
✓ Branch 6 taken 3 times.
✓ Branch 7 taken 69 times.
 72 if (!(IsPositive(options.countFraction) && IsPositive(options.sizeSlope))) {
88
1/1
✓ Branch 2 taken 3 times.
 3 throw std::invalid_argument("The secondary count fraction and size slope must be positive");
89 }
90
6/6
✓ Branch 1 taken 68 times.
✓ Branch 2 taken 1 time.
✓ Branch 3 taken 1 time.
✓ Branch 4 taken 67 times.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 67 times.
 69 if (!(IsPositive(options.minDiameter) && options.maxDiameter > options.minDiameter)) {
91
1/1
✓ Branch 2 taken 2 times.
 2 throw std::invalid_argument("Secondary diameters must satisfy 0 < min < max");
92 }
93
6/6
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 65 times.
✓ Branch 5 taken 1 time.
✓ Branch 6 taken 1 time.
✓ Branch 7 taken 1 time.
✓ Branch 8 taken 66 times.
 67 if (options.outerRange && !IsPositive(*options.outerRange)) {
94
1/1
✓ Branch 2 taken 1 time.
 1 throw std::invalid_argument("The secondary outer range must be positive");
95 }
96
1/1
✓ Branch 2 taken 64 times.
 66 ValidateCraterSelection(options.selection, "Primary");
97
5/6
✓ Branch 1 taken 63 times.
✓ Branch 2 taken 1 time.
✗ Branch 4 not taken.
✓ Branch 5 taken 63 times.
✓ Branch 6 taken 1 time.
✓ Branch 7 taken 63 times.
 64 if (!(IsPositive(options.cellSizeFactor) && IsPositive(options.maxExpected))) {
98
1/1
✓ Branch 2 taken 1 time.
 1 throw std::invalid_argument("The cell size factor and expected limit must be positive");
99 }
100 63 const SecondaryRays& rays = options.rays;
101 63 ValidateRayLayout(rays.layout);
102
6/6
✓ Branch 1 taken 61 times.
✓ Branch 2 taken 1 time.
✓ Branch 4 taken 60 times.
✓ Branch 5 taken 1 time.
✓ Branch 6 taken 4 times.
✓ Branch 7 taken 58 times.
 122 if (!(IsInUnitInterval(rays.share) && IsInUnitInterval(rays.clustering) &&
103
3/4
✓ Branch 0 taken 59 times.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 59 times.
✗ Branch 3 not taken.
 60 rays.minHalfAngleDegrees > 0.0 && rays.minHalfAngleDegrees < 90.0 &&
104
2/2
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 58 times.
 59 IsPositive(rays.clusterSpacing))) {
105 throw std::invalid_argument("Secondary rays need a share and clustering in [0, 1], a floor "
106
1/1
✓ Branch 2 taken 4 times.
 4 "angle in (0, 90) degrees and a positive cluster spacing");
107 }
108 58 const SecondaryShapes& shapes = options.shapes;
109 58 ValidateValueDistribution(shapes.ellipticity);
110
6/6
✓ Branch 1 taken 56 times.
✓ Branch 2 taken 1 time.
✓ Branch 4 taken 55 times.
✓ Branch 5 taken 1 time.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 55 times.
 112 if (!(IsInUnitInterval(shapes.elongatedShare) && shapes.ellipticity.Lowest() >= 1.0 &&
111
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 55 times.
 55 std::isfinite(shapes.ellipticity.Highest()))) {
112 throw std::invalid_argument("Secondary shapes need a share in [0, 1] and a bounded "
113
1/1
✓ Branch 2 taken 2 times.
 2 "ellipticity of at least 1");
114 }
115 55 }
116
117 481133 double SecondaryOuterRange(double diameter, const SecondaryCraterOptions& options) {
118
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 481131 times.
 481133 if (options.outerRange) {
119 2 return *options.outerRange * diameter;
120 }
121 481131 return kLongestRayOverMedian * ElliottMedianRayLength(diameter);
122 }
123
124 1080 std::optional<SecondaryRayField> SecondaryRayField::Make(const Crater& primary,
125 const SecondaryRays& rays) {
126
2/2
✓ Branch 0 taken 49 times.
✓ Branch 1 taken 1031 times.
 1080 if (!rays.enabled) {
127 49 return std::nullopt;
128 }
129 1031 SecondaryRayField field;
130
1/1
✓ Branch 1 taken 1031 times.
 1031 field.rays = RayDirections(primary, rays.layout);
131
2/2
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 1030 times.
 1031 if (field.rays.empty()) {
132 1 return std::nullopt;
133 }
134 1030 field.diameter = primary.diameter;
135 1030 field.onset = SecondaryOnsetRange(primary.diameter);
136 1030 field.share = rays.share;
137 1030 field.sinFloor = std::sin(rays.minHalfAngleDegrees * std::numbers::pi / 180.0);
138 1030 field.clustering = rays.clustering;
139 1030 field.spacing = rays.clusterSpacing;
140
1/1
✓ Branch 1 taken 1030 times.
 1030 field.meanClump = MeanClumpBump();
141 1030 return field;
142 1031 }
143
144 45984311 double SecondaryRayField::Width(const RayDirection& ray, double range) const {
145 45984311 return std::fmax(diameter * ray.sinHalfWidth, range * sinFloor);
146 }
147
148 32687417 double SecondaryRayField::Profile(const RayDirection& ray, double range, double azimuth) const {
149 32687417 const double delta = std::remainder(azimuth - ray.azimuth, 2.0 * std::numbers::pi);
150
2/2
✓ Branch 0 taken 19390523 times.
✓ Branch 1 taken 13296894 times.
 32687417 if (std::cos(delta) <= 0.0) {
151 19390523 return 0.0;
152 }
153 13296894 const double across = range * std::sin(std::fabs(delta)) / Width(ray, range);
154 13296894 return std::exp(-across * across);
155 }
156
157 30415325 double SecondaryRayField::Clump(const RayDirection& ray, double range) const {
158 30415325 const double u = std::log(range / onset) / spacing;
159 30415325 return (1.0 - clustering) + clustering * ClumpBump(u, ray.key) / meanClump;
160 }
161
162 2934190 double SecondaryRayField::MeanProfile(double range) const {
163 2934190 double sum = 0.0;
164
2/2
✓ Branch 4 taken 32687417 times.
✓ Branch 5 taken 2934190 times.
 35621607 for (const RayDirection& ray : rays) {
165 32687417 sum += std::sqrt(std::numbers::pi) * Width(ray, range);
166 }
167 2934190 return sum / (2.0 * std::numbers::pi * range);
168 }
169
170 2727802 double SecondaryRayField::Density(double range, double azimuth) const {
171 2727802 double sum = 0.0;
172
2/2
✓ Branch 4 taken 30415325 times.
✓ Branch 5 taken 2727802 times.
 33143127 for (const RayDirection& ray : rays) {
173
1/1
✓ Branch 2 taken 30415325 times.
 30415325 sum += Profile(ray, range, azimuth) * Clump(ray, range);
174 }
175 2727802 return (1.0 - share) + share * sum / MeanProfile(range);
176 }
177
178 6251 double SecondaryRayField::Peak(double r0, double r1) const {
179 6251 double peak = 0.0;
180
2/2
✓ Branch 2 taken 18753 times.
✓ Branch 3 taken 6251 times.
 25004 for (const double range : {r0, 0.5 * (r0 + r1), r1}) {
181
2/2
✓ Branch 4 taken 206388 times.
✓ Branch 5 taken 18753 times.
 225141 for (const RayDirection& at : rays) {
182 206388 double sum = 0.0;
183
2/2
✓ Branch 4 taken 2272092 times.
✓ Branch 5 taken 206388 times.
 2478480 for (const RayDirection& ray : rays) {
184 2272092 sum += Profile(ray, range, at.azimuth);
185 }
186 206388 peak = std::fmax(peak, sum / MeanProfile(range));
187 }
188 }
189 6251 const double clump = (1.0 - clustering) + clustering / meanClump;
190 6251 return (1.0 - share) + share * 1.1 * peak * clump;
191 }
192
193 8 SecondaryDensity::SecondaryDensity(const Crater& primary, const SecondaryCraterOptions& options)
194 8 : sizeSlope_(options.sizeSlope), minSize_(options.minDiameter), maxSize_(options.maxDiameter),
195 8 diameter_(primary.diameter) {
196
6/7
✓ Branch 1 taken 8 times.
✓ Branch 3 taken 7 times.
✓ Branch 4 taken 1 time.
✗ Branch 6 not taken.
✓ Branch 7 taken 7 times.
✓ Branch 8 taken 1 time.
✓ Branch 9 taken 7 times.
 8 if (!options.selection.Includes(primary) || !options.onsetCount) {
197 1 return;
198 }
199 7 onset_ = SecondaryOnsetRange(diameter_);
200 7 const double outer = SecondaryOuterRange(diameter_, options);
201
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 7 times.
 7 if (!(outer > onset_)) {
202 ✗ return;
203 }
204 7 outer_ = outer;
205 7 density0_ = *options.onsetCount / (0.25 * std::numbers::pi * diameter_ * diameter_);
206 7 reference_ = options.countFraction * diameter_;
207
1/1
✓ Branch 1 taken 7 times.
 7 rays_ = SecondaryRayField::Make(primary, options.rays);
208 }
209
210 503568 double SecondaryDensity::Measure(const Vec2d& offset, double dMin, double dMax, bool cover) const {
211 503568 const double range = std::hypot(offset.x, offset.y);
212
3/4
✓ Branch 0 taken 441638 times.
✓ Branch 1 taken 61930 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 441638 times.
 503568 if (!(range >= onset_ && range < outer_)) {
213 61930 return 0.0;
214 }
215 441638 const double low = std::fmax(dMin, minSize_);
216 const double high =
217 441638 std::fmin(std::fmin(dMax, maxSize_), SingerMaxSecondaryDiameter(diameter_, range));
218
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 441638 times.
 441638 if (!(high > low)) {
219 ✗ return 0.0;
220 }
221 441638 const double r = 1.0 + (range - onset_) / diameter_;
222 441638 const double s = sizeSlope_;
223 441638 double sizes = std::pow(low / reference_, -s) - std::pow(high / reference_, -s);
224
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 441636 times.
 441638 if (cover) {
225 2 const double scale = 0.25 * std::numbers::pi * s * std::pow(reference_, s);
226
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 1 time.
 3 sizes = s == 2.0 ? scale * std::log(high / low)
227 1 : scale * (std::pow(high, 2.0 - s) - std::pow(low, 2.0 - s)) / (2.0 - s);
228 }
229 441638 const double density = density0_ / (r * r) * sizes;
230
2/2
✓ Branch 1 taken 90000 times.
✓ Branch 2 taken 351638 times.
 441638 return rays_ ? density * rays_->Density(range, std::atan2(offset.y, offset.x)) : density;
231 }
232
233 503566 double SecondaryDensity::At(const Vec2d& offset, double dMin, double dMax) const {
234 503566 return Measure(offset, dMin, dMax, false);
235 }
236
237 2 double SecondaryDensity::CoverAt(const Vec2d& offset, double dMin, double dMax) const {
238 2 return Measure(offset, dMin, dMax, true);
239 }
240
241 2 double SecondaryDensity::Largest() const {
242
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if (!(outer_ > 0.0)) {
243 ✗ return 0.0;
244 }
245 2 return std::fmin(maxSize_, std::fmax(SingerMaxSecondaryDiameter(diameter_, onset_),
246 4 SingerMaxSecondaryDiameter(diameter_, outer_)));
247 }
248
249 202 double LargestSecondary(double diameter, const SecondaryCraterOptions& options) {
250 202 const double onset = SecondaryOnsetRange(diameter);
251 202 const double outer = SecondaryOuterRange(diameter, options);
252
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 200 times.
 202 if (!(outer > onset)) {
253 2 return 0.0;
254 }
255 200 return std::fmin(options.maxDiameter, std::fmax(SingerMaxSecondaryDiameter(diameter, onset),
256 200 SingerMaxSecondaryDiameter(diameter, outer)));
257 }
258
259 } // namespace terraingen
260