TerrainGen coverage

terraingen/
Coverage:
low: ≥ 0%
medium: ≥ 75.0%
high: ≥ 90.0%
Lines:
135 / 135
100.0%
Branches:
125 / 133
94.0%
Functions:
12 / 12
100.0%
albedo/raster_albedo.cpp
Line Branch Exec Source
1 #include "terraingen/albedo/raster_albedo.h"
2
3 #include "terraingen/layers/roughness.h"
4 #include "terraingen/noise/hash.h"
5 #include "terraingen/noise/perlin.h"
6 #include "terraingen/util/checks.h"
7
8 #include <algorithm>
9 #include <array>
10 #include <cmath>
11 #include <stdexcept>
12 #include <utility>
13
14 namespace terraingen {
15
16 namespace {
17
18 /// Seed of the fixed sample the RMS of `PerlinNoise2` is measured on.
19 constexpr std::uint32_t kRmsSeed = 0x5eed1234u;
20 /// Step between the seeds of successive detail octaves (2^32 over the golden ratio).
21 constexpr std::uint32_t kOctaveSeedStep = 0x9e3779b9u;
22
23 /// The lags, in cells, the map's structure function is measured at.
24 constexpr std::array<std::size_t, 4> kLags{1, 2, 4, 8};
25
26 /// The RMS of `PerlinNoise2` over a fixed set of points, so octaves can be given an RMS.
27 152035 double PerlinRms() {
28 5 static const double rms = [] {
29 5 double sum = 0.0;
30 5 constexpr int kSide = 256;
31
2/2
✓ Branch 0 taken 1280 times.
✓ Branch 1 taken 5 times.
 1285 for (int j = 0; j < kSide; ++j) {
32
2/2
✓ Branch 0 taken 327680 times.
✓ Branch 1 taken 1280 times.
 328960 for (int i = 0; i < kSide; ++i) {
33 const float v =
34 327680 PerlinNoise2(i / 8, j / 8, static_cast<float>(i % 8) / 8.0f + 1.0f / 16.0f,
35 327680 static_cast<float>(j % 8) / 8.0f + 1.0f / 16.0f, kRmsSeed);
36 327680 sum += static_cast<double>(v) * static_cast<double>(v);
37 }
38 }
39 5 return std::sqrt(sum / (kSide * kSide));
40
4/4
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 152019 times.
✓ Branch 3 taken 5 times.
✓ Branch 4 taken 11 times.
 152035 }();
41 152035 return rms;
42 }
43
44 /// The logarithm of the mean of the bands at a cell, or NaN where it has no data.
45 14526 double LogBrightness(const std::vector<Raster<float>>& bands, std::size_t col, std::size_t row) {
46 14526 double sum = 0.0;
47
2/2
✓ Branch 4 taken 15043 times.
✓ Branch 5 taken 14524 times.
 29567 for (const Raster<float>& band : bands) {
48 15043 const float value = band(col, row);
49
6/8
✓ Branch 1 taken 15041 times.
✓ Branch 2 taken 2 times.
✓ Branch 4 taken 15041 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 15041 times.
✓ Branch 8 taken 2 times.
✓ Branch 9 taken 15041 times.
 15043 if (band.IsNoData(value) || !std::isfinite(value) || value <= 0.0f) {
50 2 return NAN;
51 }
52 15041 sum += static_cast<double>(value);
53 }
54 14524 return std::log(sum / static_cast<double>(bands.size()));
55 }
56
57 } // namespace
58
59 19 std::vector<double> LogBrightnessGrid(const std::vector<Raster<float>>& bands) {
60 19 const std::size_t width = bands.front().Width();
61 19 const std::size_t height = bands.front().Height();
62
1/1
✓ Branch 1 taken 19 times.
 19 std::vector<double> logs(width * height);
63
2/2
✓ Branch 0 taken 436 times.
✓ Branch 1 taken 19 times.
 455 for (std::size_t row = 0; row < height; ++row) {
64
2/2
✓ Branch 0 taken 14526 times.
✓ Branch 1 taken 436 times.
 14962 for (std::size_t col = 0; col < width; ++col) {
65
1/1
✓ Branch 1 taken 14526 times.
 14526 logs[row * width + col] = LogBrightness(bands, col, row);
66 }
67 }
68 19 return logs;
69 }
70
71 75 std::optional<double> StructureFunction(const std::vector<double>& values, std::size_t width,
72 std::size_t height, std::size_t lag) {
73 75 double sum = 0.0;
74 75 std::size_t pairs = 0;
75 103158 const auto add = [&](double a, double b) {
76
6/6
✓ Branch 1 taken 103150 times.
✓ Branch 2 taken 8 times.
✓ Branch 4 taken 103146 times.
✓ Branch 5 taken 4 times.
✓ Branch 6 taken 103146 times.
✓ Branch 7 taken 12 times.
 103158 if (std::isfinite(a) && std::isfinite(b)) {
77 103146 sum += (a - b) * (a - b);
78 103146 ++pairs;
79 }
80 103233 };
81
2/2
✓ Branch 0 taken 1742 times.
✓ Branch 1 taken 75 times.
 1817 for (std::size_t row = 0; row < height; ++row) {
82
2/2
✓ Branch 0 taken 58098 times.
✓ Branch 1 taken 1742 times.
 59840 for (std::size_t col = 0; col < width; ++col) {
83 58098 const double here = values[row * width + col];
84
2/2
✓ Branch 0 taken 51610 times.
✓ Branch 1 taken 6488 times.
 58098 if (col + lag < width) {
85 51610 add(here, values[row * width + col + lag]);
86 }
87
2/2
✓ Branch 0 taken 51548 times.
✓ Branch 1 taken 6550 times.
 58098 if (row + lag < height) {
88 51548 add(here, values[(row + lag) * width + col]);
89 }
90 }
91 }
92
4/4
✓ Branch 0 taken 68 times.
✓ Branch 1 taken 7 times.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 64 times.
 75 if (!(pairs > 0 && sum > 0.0)) {
93 11 return std::nullopt;
94 }
95 64 return sum / static_cast<double>(pairs);
96 }
97
98 17 double FitSlope(const std::vector<double>& x, const std::vector<double>& y) {
99 17 const auto n = static_cast<double>(x.size());
100 17 double mx = 0.0;
101 17 double my = 0.0;
102
2/2
✓ Branch 1 taken 65 times.
✓ Branch 2 taken 17 times.
 82 for (std::size_t i = 0; i < x.size(); ++i) {
103 65 mx += x[i] / n;
104 65 my += y[i] / n;
105 }
106 17 double sxy = 0.0;
107 17 double sxx = 0.0;
108
2/2
✓ Branch 1 taken 65 times.
✓ Branch 2 taken 17 times.
 82 for (std::size_t i = 0; i < x.size(); ++i) {
109 65 sxy += (x[i] - mx) * (y[i] - my);
110 65 sxx += (x[i] - mx) * (x[i] - mx);
111 }
112 17 return sxy / sxx;
113 }
114
115 29 void ValidateAlbedoDetail(const AlbedoDetail& detail) {
116
3/4
✓ Branch 0 taken 26 times.
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 26 times.
 29 if (detail.octaves < 1 || detail.octaves > 16) {
117
1/1
✓ Branch 2 taken 3 times.
 3 throw std::invalid_argument("Albedo detail needs 1 to 16 octaves");
118 }
119
2/2
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 25 times.
 26 if (!IsNonNegative(detail.scale)) {
120
1/1
✓ Branch 2 taken 1 time.
 1 throw std::invalid_argument("The albedo detail scale must not be negative");
121 }
122
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 23 times.
 25 if (detail.hurst) {
123
1/1
✓ Branch 3 taken 1 time.
 2 RequireUnitInterval(*detail.hurst, "The albedo detail Hurst exponent");
124 }
125 24 ValidateResolutionFade(detail.fade);
126 23 }
127
128 26 RasterAlbedo::RasterAlbedo(std::vector<Raster<float>> bands, RasterAlbedoOptions options,
129 26 Surface surface)
130 26 : bands_(std::move(bands)), options_(options), surface_(surface) {
131
6/6
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 23 times.
✓ Branch 4 taken 1 time.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 1 time.
✓ Branch 7 taken 25 times.
 26 if (bands_.size() != 1 && bands_.size() != kMaxAlbedoBands) {
132
1/1
✓ Branch 2 taken 1 time.
 1 throw std::invalid_argument("An albedo map must have 1 or 3 bands");
133 }
134 25 const Raster<float>& first = bands_.front();
135
6/6
✓ Branch 1 taken 24 times.
✓ Branch 2 taken 1 time.
✓ Branch 4 taken 23 times.
✓ Branch 5 taken 1 time.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 23 times.
 48 if (first.Empty() || first.Transform().rowRotation != 0.0 ||
136
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 23 times.
 23 first.Transform().colRotation != 0.0) {
137
1/1
✓ Branch 2 taken 2 times.
 2 throw std::invalid_argument("An albedo map must be neither empty nor rotated");
138 }
139
2/2
✓ Branch 5 taken 26 times.
✓ Branch 6 taken 22 times.
 48 for (const Raster<float>& band : bands_) {
140
5/6
✓ Branch 2 taken 25 times.
✓ Branch 3 taken 1 time.
✓ Branch 6 taken 25 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 time.
✓ Branch 9 taken 25 times.
 51 if (band.Width() != first.Width() || band.Height() != first.Height() ||
141
1/2
✗ Branch 3 not taken.
✓ Branch 4 taken 25 times.
 25 !(band.Transform() == first.Transform())) {
142
1/1
✓ Branch 2 taken 1 time.
 1 throw std::invalid_argument("The bands of an albedo map must share one grid");
143 }
144 }
145
1/1
✓ Branch 1 taken 18 times.
 22 ValidateAlbedoDetail(options_.detail);
146
2/2
✓ Branch 1 taken 18 times.
✓ Branch 4 taken 18 times.
 18 cell_ = RasterGrid::FromRaster(first).Spacing(surface_);
147
1/1
✓ Branch 1 taken 18 times.
 18 MeasureRoughness();
148
2/2
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 17 times.
 18 if (options_.detail.hurst) {
149 1 hurst_ = *options_.detail.hurst;
150 }
151 42 }
152
153 18 void RasterAlbedo::MeasureRoughness() {
154 18 const std::size_t width = bands_.front().Width();
155 18 const std::size_t height = bands_.front().Height();
156
1/1
✓ Branch 1 taken 18 times.
 18 const std::vector<double> logs = LogBrightnessGrid(bands_);
157 18 std::vector<double> logLags;
158 18 std::vector<double> logStructure;
159
2/2
✓ Branch 0 taken 72 times.
✓ Branch 1 taken 18 times.
 90 for (const std::size_t lag : kLags) {
160
2/2
✓ Branch 2 taken 63 times.
✓ Branch 3 taken 9 times.
 72 if (const std::optional<double> structure = StructureFunction(logs, width, height, lag)) {
161
1/1
✓ Branch 1 taken 63 times.
 63 logLags.push_back(std::log(static_cast<double>(lag)));
162
1/1
✓ Branch 2 taken 63 times.
 63 logStructure.push_back(std::log(*structure));
163 }
164 }
165
5/6
✓ Branch 1 taken 17 times.
✓ Branch 2 taken 1 time.
✓ Branch 4 taken 17 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 17 times.
✓ Branch 7 taken 1 time.
 18 if (!logLags.empty() && logLags.front() == 0.0) {
166 17 roughness_ = std::sqrt(0.5 * std::exp(logStructure.front()));
167 }
168
2/2
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 2 times.
 18 if (logLags.size() >= 2) {
169
1/1
✓ Branch 2 taken 16 times.
 16 hurst_ = std::clamp(0.5 * FitSlope(logLags, logStructure), 0.0, 1.0);
170 }
171 18 }
172
173 153638 double RasterAlbedo::DetailSum(double x, double y, double spacing) const {
174 153638 const AlbedoDetail& detail = options_.detail;
175 153638 const Vec2d origin = bands_.front().CellCenter(0, 0);
176 153638 double sum = 0.0;
177
2/2
✓ Branch 0 taken 288963 times.
✓ Branch 1 taken 66591 times.
 355554 for (int octave = 0; octave < detail.octaves; ++octave) {
178 288963 const double wavelength = cell_ / std::exp2(octave + 1);
179
1/1
✓ Branch 1 taken 288963 times.
 288963 const double weight = ResolutionWeight(wavelength, spacing, detail.fade);
180
2/2
✓ Branch 0 taken 87047 times.
✓ Branch 1 taken 201916 times.
 288963 if (weight == 0.0) {
181 87047 break;
182 }
183 const std::uint32_t seed =
184 201916 PcgHash(detail.seed + kOctaveSeedStep * static_cast<std::uint32_t>(octave));
185 201916 float noise = 0.0f;
186
2/2
✓ Branch 1 taken 1605 times.
✓ Branch 2 taken 200311 times.
 201916 if (surface_.IsBody()) {
187
1/1
✓ Branch 1 taken 1605 times.
 1605 const Vec3d p = surface_.Position({x, y});
188
1/1
✓ Branch 1 taken 1605 times.
 1605 noise = PerlinNoise3Precise(p.x / wavelength, p.y / wavelength, p.z / wavelength, seed);
189 } else {
190 noise =
191 200311 PerlinNoise2Precise((x - origin.x) / wavelength, (y - origin.y) / wavelength, seed);
192 }
193 201916 sum += weight * std::exp2(-(octave + 1) * hurst_) * static_cast<double>(noise);
194 }
195 153638 return sum;
196 }
197
198 156218 AlbedoSample RasterAlbedo::ValueAt(double x, double y, const AlbedoContext& context) const {
199 156218 AlbedoSample value{};
200
2/2
✓ Branch 1 taken 158017 times.
✓ Branch 2 taken 155385 times.
 313402 for (std::size_t band = 0; band < bands_.size(); ++band) {
201
1/1
✓ Branch 2 taken 158017 times.
 158017 value[band] = SampleRaster(bands_[band], x, y, options_.method);
202
2/2
✓ Branch 2 taken 833 times.
✓ Branch 3 taken 157184 times.
 158017 if (std::isnan(value[band])) {
203
1/1
✓ Branch 1 taken 833 times.
 833 value.fill(NAN);
204 833 return value;
205 }
206 }
207 155385 const AlbedoDetail& detail = options_.detail;
208
5/6
✓ Branch 0 taken 153639 times.
✓ Branch 1 taken 1746 times.
✓ Branch 2 taken 153638 times.
✓ Branch 3 taken 1 time.
✗ Branch 4 not taken.
✓ Branch 5 taken 153638 times.
 155385 if (!detail.enabled || detail.scale == 0.0 || roughness_ == 0.0) {
209 1747 return value;
210 }
211
1/1
✓ Branch 1 taken 153638 times.
 153638 const double sum = DetailSum(x, y, context.spacing);
212 const double rms =
213
3/3
✓ Branch 1 taken 1603 times.
✓ Branch 2 taken 152035 times.
✓ Branch 4 taken 1603 times.
 153638 surface_.IsBody() ? std::sqrt(0.5 * Perlin3StructureFunction(1e6)) : PerlinRms();
214 153638 const double factor = std::exp(detail.scale * roughness_ * sum / rms);
215
2/2
✓ Branch 1 taken 155180 times.
✓ Branch 2 taken 153638 times.
 308818 for (std::size_t band = 0; band < bands_.size(); ++band) {
216 155180 value[band] = static_cast<float>(static_cast<double>(value[band]) * factor);
217 }
218 153638 return value;
219 }
220
221 } // namespace terraingen
222