TerrainGen coverage


terraingen/
Coverage:
low: ≥ 0%
medium: ≥ 75.0%
high: ≥ 90.0%
Lines:
212 / 212
100.0%
Branches:
143 / 147
97.3%
Functions:
20 / 20
100.0%

util/cube_sphere_cells.cpp
Line Branch Exec Source
1 #include "terraingen/util/cells.h"
2
3 #include "terraingen/util/checks.h"
4 #include "terraingen/util/cube_faces.h"
5 #include "terraingen/util/surface.h"
6
7 #include <algorithm>
8 #include <array>
9 #include <cmath>
10 #include <limits>
11 #include <numbers>
12 #include <stdexcept>
13
14 namespace terraingen {
15
16 namespace {
17
18 constexpr double kQuarter = std::numbers::pi / 4.0;
19
20 /// Cells per face edge over which an ellipsoid's global bound of the area per steradian is taken.
21 constexpr std::int32_t kCapCells = 64;
22
23 /// The solid angle of the gnomonic rectangle [0, `a`] x [0, `b`] on a face, signed.
24 108900008 double Corner(double a, double b) {
25 108900008 return std::atan(a * b / std::sqrt(1.0 + a * a + b * b));
26 }
27
28 /// The gnomonic coordinate of the edge `k` of `n` cells along a face.
29 3645388 double Edge(std::int32_t k, std::int32_t n) {
30 3645388 return std::tan(-kQuarter + 2.0 * kQuarter * static_cast<double>(k) / static_cast<double>(n));
31 }
32
33 46841 void SortUnique(std::vector<CellIndex>& cells) {
34 46841 std::sort(cells.begin(), cells.end());
35
2/2
✓ Branch 5 taken 46841 times.
✓ Branch 9 taken 46841 times.
46841 cells.erase(std::unique(cells.begin(), cells.end()), cells.end());
36 46841 }
37
38 } // namespace
39
40 11 CubeSphereCells::CubeSphereCells(double radius)
41
1/1
✓ Branch 2 taken 10 times.
11 : radius_(radius), surface_(Surface::Sphere(radius)) {}
42
43 52 CubeSphereCells::CubeSphereCells(const Surface& surface) : radius_(0.0), surface_(surface) {
44
2/2
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 51 times.
52 if (!surface.IsBody()) {
45
1/1
✓ Branch 2 taken 1 time.
1 throw std::invalid_argument("Cube-sphere cells need a body, not a plane");
46 }
47
1/1
✓ Branch 1 taken 51 times.
51 radius_ = surface.MinRadius();
48
2/2
✓ Branch 1 taken 42 times.
✓ Branch 2 taken 9 times.
51 if (surface.IsSphere()) {
49 42 return;
50 }
51
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 7 times.
9 if (surface.Shape() != nullptr) {
52 2 weightCap_ = surface.Shape()->MaxAreaPerSteradian();
53 2 return;
54 }
55
2/2
✓ Branch 0 taken 42 times.
✓ Branch 1 taken 7 times.
98 for (std::size_t face = 0; face < detail::kCubeFaces.size(); ++face) {
56
2/2
✓ Branch 0 taken 2688 times.
✓ Branch 1 taken 42 times.
2730 for (std::int32_t i = 0; i < kCapCells; ++i) {
57
2/2
✓ Branch 0 taken 172032 times.
✓ Branch 1 taken 2688 times.
174720 for (std::int32_t j = 0; j < kCapCells; ++j) {
58 344064 weightCap_ = std::max(weightCap_,
59
1/1
✓ Branch 5 taken 172032 times.
172032 WeightBound(face, Edge(i, kCapCells), Edge(i + 1, kCapCells),
60 Edge(j, kCapCells), Edge(j + 1, kCapCells)));
61 }
62 }
63 }
64 2 }
65
66 325240 std::array<double, 4> CubeSphereCells::Bounds(const CellIndex& cell, double cellSize) const {
67 325240 const std::int32_t n = CellsPerEdge(cellSize);
68 325240 const std::int32_t col = cell.i % n;
69 325240 return {Edge(col, n), Edge(col + 1, n), Edge(cell.j, n), Edge(cell.j + 1, n)};
70 }
71
72 251403 double CubeSphereCells::Weight(const Vec2d& point) const {
73 251403 return surface_.AreaPerSteradian(point);
74 }
75
76 1299612 std::int32_t CubeSphereCells::CellsPerEdge(double cellSize) const {
77
1/1
✓ Branch 2 taken 1299611 times.
1299612 RequirePositive(cellSize, "Cell size");
78 1299611 const double n = std::ceil(radius_ * 2.0 * kQuarter / cellSize);
79
2/2
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 1299610 times.
1299611 if (!(n < static_cast<double>(std::numeric_limits<std::int32_t>::max() / 6))) {
80
1/1
✓ Branch 2 taken 1 time.
1 throw std::invalid_argument("Cells are too small for the sphere");
81 }
82 1299610 return std::max<std::int32_t>(1, static_cast<std::int32_t>(n));
83 }
84
85 511241 CellIndex CubeSphereCells::CellOf(double lonDeg, double latDeg, double cellSize) const {
86
1/1
✓ Branch 1 taken 511241 times.
511241 const std::int32_t n = CellsPerEdge(cellSize);
87
1/1
✓ Branch 1 taken 511241 times.
511241 const Vec3d p = UnitVector(lonDeg, latDeg);
88 511241 const std::size_t face = detail::CubeFaceOf(detail::kCubeFaces, p);
89 511241 const double best = Dot(p, detail::kCubeFaces[face].normal);
90 1022482 const auto index = [&](const Vec3d& axis) {
91 1022482 const double angle = std::atan(Dot(p, axis) / best);
92 1022482 const auto k = static_cast<std::int32_t>(
93 1022482 std::floor((angle + kQuarter) / (2.0 * kQuarter) * static_cast<double>(n)));
94
1/1
✓ Branch 1 taken 1022482 times.
1022482 return std::clamp<std::int32_t>(k, 0, n - 1);
95 511241 };
96
1/1
✓ Branch 2 taken 511241 times.
511241 return {static_cast<std::int32_t>(face) * n + index(detail::kCubeFaces[face].u),
97
1/1
✓ Branch 2 taken 511241 times.
1022482 index(detail::kCubeFaces[face].v)};
98 }
99
100 39393 std::vector<CellIndex> CubeSphereCells::Neighbours(const CellIndex& cell, double cellSize,
101 int rings) const {
102 39393 const std::int32_t n = CellsPerEdge(cellSize);
103 39393 const std::int32_t face = cell.i / n;
104 39393 const detail::CubeFace& f = detail::kCubeFaces[static_cast<std::size_t>(face)];
105 39393 const std::int32_t col = cell.i % n;
106 39393 const double step = 2.0 * kQuarter / static_cast<double>(n);
107 39393 std::vector<CellIndex> cells;
108
2/2
✓ Branch 0 taken 140011 times.
✓ Branch 1 taken 39393 times.
179404 for (std::int32_t dj = -rings; dj <= rings; ++dj) {
109
2/2
✓ Branch 0 taken 529193 times.
✓ Branch 1 taken 140011 times.
669204 for (std::int32_t di = -rings; di <= rings; ++di) {
110 529193 const std::int32_t c = col + di;
111 529193 const std::int32_t r = cell.j + dj;
112
8/8
✓ Branch 0 taken 495976 times.
✓ Branch 1 taken 33217 times.
✓ Branch 2 taken 462595 times.
✓ Branch 3 taken 33381 times.
✓ Branch 4 taken 446748 times.
✓ Branch 5 taken 15847 times.
✓ Branch 6 taken 430320 times.
✓ Branch 7 taken 16428 times.
529193 if (c >= 0 && c < n && r >= 0 && r < n) {
113
1/1
✓ Branch 1 taken 430320 times.
430320 cells.push_back({face * n + c, r});
114 430320 continue;
115 }
116
2/2
✓ Branch 2 taken 296619 times.
✓ Branch 3 taken 98873 times.
395492 for (const double su : {0.1, 0.5, 0.9}) {
117
2/2
✓ Branch 2 taken 889857 times.
✓ Branch 3 taken 296619 times.
1186476 for (const double sv : {0.1, 0.5, 0.9}) {
118 889857 const double alpha = -kQuarter + (static_cast<double>(c) + su) * step;
119 889857 const double beta = -kQuarter + (static_cast<double>(r) + sv) * step;
120
2/2
✓ Branch 0 taken 236853 times.
✓ Branch 1 taken 653004 times.
889857 if (!(std::fabs(alpha) < 2.0 * kQuarter - step &&
121
2/2
✓ Branch 0 taken 115225 times.
✓ Branch 1 taken 121628 times.
236853 std::fabs(beta) < 2.0 * kQuarter - step)) {
122 768229 continue;
123 }
124 121628 const Vec3d p = detail::CubeFacePoint(f, std::tan(alpha), std::tan(beta));
125
2/2
✓ Branch 1 taken 121628 times.
✓ Branch 4 taken 121628 times.
121628 const Vec2d point = Surface::Sphere(radius_).PointOf(p);
126
2/2
✓ Branch 1 taken 121628 times.
✓ Branch 4 taken 121628 times.
121628 cells.push_back(CellOf(point.x, point.y, cellSize));
127 }
128 }
129 }
130 }
131
1/1
✓ Branch 1 taken 39393 times.
39393 SortUnique(cells);
132 39393 return cells;
133 }
134
135 3096 std::vector<CellIndex> CubeSphereCells::CellsCovering(const Bounds2d& region,
136 double cellSize) const {
137
1/1
✓ Branch 2 taken 3096 times.
3096 RequirePositive(cellSize, "Cell size");
138
6/8
✓ Branch 1 taken 3096 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3096 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3095 times.
✓ Branch 8 taken 1 time.
✓ Branch 9 taken 2 times.
✓ Branch 10 taken 3094 times.
6191 if (!(std::isfinite(region.minX) && std::isfinite(region.minY) && std::isfinite(region.maxX) &&
139
4/6
✓ Branch 1 taken 3095 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3094 times.
✓ Branch 4 taken 1 time.
✗ Branch 5 not taken.
✓ Branch 6 taken 3094 times.
3095 std::isfinite(region.maxY) && region.maxX >= region.minX && region.maxY >= region.minY)) {
140
1/1
✓ Branch 2 taken 2 times.
2 throw std::invalid_argument("Region must be finite and not inverted");
141 }
142
1/1
✓ Branch 1 taken 3094 times.
3094 const Surface sphere = Surface::Sphere(radius_);
143 9282 const Bounds2d box{region.minX, std::max(-90.0, region.minY), region.maxX,
144 3094 std::min(90.0, region.maxY)};
145 3094 const double step = cellSize / 4.0 / sphere.UnitLength();
146 3094 const auto rows = static_cast<std::int64_t>(std::ceil((box.maxY - box.minY) / step));
147 3094 std::vector<CellIndex> seeds;
148
2/2
✓ Branch 0 taken 9308 times.
✓ Branch 1 taken 3094 times.
12402 for (std::int64_t r = 0; r <= rows; ++r) {
149 9308 const double lat = std::min(box.maxY, box.minY + static_cast<double>(r) * step);
150 9308 const double width = std::min(box.maxX - box.minX, 360.0);
151 9308 const double cosLat = std::cos(lat * kDegree);
152 9308 const auto columns = static_cast<std::int64_t>(std::ceil(width * cosLat / step));
153
2/2
✓ Branch 0 taken 212454 times.
✓ Branch 1 taken 9308 times.
221762 for (std::int64_t c = 0; c <= columns; ++c) {
154
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 212453 times.
212454 const double lon = columns == 0 ? box.minX
155 212453 : box.minX + width * static_cast<double>(c) /
156 212453 static_cast<double>(columns);
157
2/2
✓ Branch 1 taken 212454 times.
✓ Branch 4 taken 212454 times.
212454 seeds.push_back(CellOf(lon, lat, cellSize));
158 }
159 }
160
1/1
✓ Branch 1 taken 3094 times.
3094 SortUnique(seeds);
161 3094 std::vector<CellIndex> cells;
162
2/2
✓ Branch 5 taken 28475 times.
✓ Branch 6 taken 3094 times.
31569 for (const CellIndex& seed : seeds) {
163
1/1
✓ Branch 1 taken 28475 times.
28475 const std::vector<CellIndex> around = Neighbours(seed, cellSize, 1);
164
1/1
✓ Branch 5 taken 28475 times.
28475 cells.insert(cells.end(), around.begin(), around.end());
165 28475 }
166
1/1
✓ Branch 1 taken 3094 times.
3094 SortUnique(cells);
167 // A ring cell reaches the region only if its centre lies within its farthest corner's angle of
168 // it, since cell edges are great circles; one expansion by the largest such angle tests them.
169
1/1
✓ Branch 1 taken 3094 times.
3094 const std::int32_t n = CellsPerEdge(cellSize);
170 3094 std::vector<CellIndex> ring;
171 3094 std::vector<Vec2d> centres;
172 3094 double widest = 0.0;
173
2/2
✓ Branch 5 taken 43588 times.
✓ Branch 6 taken 3094 times.
46682 for (const CellIndex& cell : cells) {
174
3/3
✓ Branch 3 taken 43588 times.
✓ Branch 5 taken 28475 times.
✓ Branch 6 taken 15113 times.
43588 if (std::binary_search(seeds.begin(), seeds.end(), cell)) {
175 28475 continue;
176 }
177
1/1
✓ Branch 1 taken 15113 times.
15113 const auto [a1, a2, b1, b2] = Bounds(cell, cellSize);
178 15113 const detail::CubeFace& face = detail::kCubeFaces[static_cast<std::size_t>(cell.i / n)];
179 const Vec3d centre =
180 15113 Normalized(detail::CubeFacePoint(face, 0.5 * (a1 + a2), 0.5 * (b1 + b2)));
181 15113 double nearest = 1.0;
182
2/2
✓ Branch 2 taken 30226 times.
✓ Branch 3 taken 15113 times.
45339 for (const double a : {a1, a2}) {
183
2/2
✓ Branch 2 taken 60452 times.
✓ Branch 3 taken 30226 times.
90678 for (const double b : {b1, b2}) {
184 60452 nearest =
185 60452 std::min(nearest, Dot(centre, Normalized(detail::CubeFacePoint(face, a, b))));
186 }
187 }
188
1/1
✓ Branch 1 taken 15113 times.
15113 widest = std::max(widest, std::acos(std::clamp(nearest, -1.0, 1.0)));
189
1/1
✓ Branch 1 taken 15113 times.
15113 ring.push_back(cell);
190
2/2
✓ Branch 1 taken 15113 times.
✓ Branch 4 taken 15113 times.
15113 centres.push_back(sphere.PointOf(centre));
191 }
192
2/2
✓ Branch 1 taken 1834 times.
✓ Branch 2 taken 1260 times.
3094 if (ring.empty()) {
193 1834 return cells;
194 }
195
1/1
✓ Branch 1 taken 1260 times.
1260 const Bounds2d reach = sphere.Expand(box, radius_ * (widest * (1.0 + 1e-9) + 1e-12));
196
1/1
✓ Branch 1 taken 1260 times.
1260 std::vector<CellIndex> kept = seeds;
197
2/2
✓ Branch 1 taken 15113 times.
✓ Branch 2 taken 1260 times.
16373 for (std::size_t k = 0; k < ring.size(); ++k) {
198
2/2
✓ Branch 2 taken 2769 times.
✓ Branch 3 taken 12344 times.
15113 if (sphere.Contains(reach, centres[k])) {
199
1/1
✓ Branch 2 taken 2769 times.
2769 kept.push_back(ring[k]);
200 }
201 }
202
1/1
✓ Branch 1 taken 1260 times.
1260 SortUnique(kept);
203 1260 return kept;
204 3094 }
205
206 309165 double CubeSphereCells::CellArea(const CellIndex& cell, double cellSize) const {
207
1/1
✓ Branch 1 taken 309165 times.
309165 const auto [a1, a2, b1, b2] = Bounds(cell, cellSize);
208 309165 const double solid = Corner(a2, b2) - Corner(a1, b2) - Corner(a2, b1) + Corner(a1, b1);
209
2/2
✓ Branch 1 taken 303572 times.
✓ Branch 2 taken 5593 times.
309165 if (surface_.IsSphere()) {
210 303572 return radius_ * radius_ * solid;
211 }
212 5593 constexpr std::array<double, 5> kNodes{-0.906179845938664, -0.538469310105683, 0.0,
213 0.538469310105683, 0.906179845938664};
214 5593 constexpr std::array<double, 5> kWeights{0.236926885056189, 0.478628670499366,
215 0.568888888888889, 0.478628670499366,
216 0.236926885056189};
217 const detail::CubeFace& face =
218
1/1
✓ Branch 1 taken 5593 times.
5593 detail::kCubeFaces[static_cast<std::size_t>(cell.i / CellsPerEdge(cellSize))];
219 5593 double area = 0.0;
220
2/2
✓ Branch 0 taken 27965 times.
✓ Branch 1 taken 5593 times.
67116 for (std::size_t i = 0; i < kNodes.size(); ++i) {
221 27965 const double a = 0.5 * (a1 + a2) + 0.5 * (a2 - a1) * kNodes[i];
222
2/2
✓ Branch 0 taken 139825 times.
✓ Branch 1 taken 27965 times.
195755 for (std::size_t j = 0; j < kNodes.size(); ++j) {
223 139825 const double b = 0.5 * (b1 + b2) + 0.5 * (b2 - b1) * kNodes[j];
224 139825 const double jacobian = std::pow(1.0 + a * a + b * b, -1.5);
225 const double weight =
226
2/2
✓ Branch 2 taken 139825 times.
✓ Branch 5 taken 139825 times.
139825 surface_.AreaPerSteradian(surface_.PointOf(detail::CubeFacePoint(face, a, b)));
227 139825 area += kWeights[i] * kWeights[j] * weight * jacobian;
228 }
229 }
230 5593 return area * 0.25 * (a2 - a1) * (b2 - b1);
231 }
232
233 172897 double CubeSphereCells::WeightBound(std::size_t face, double a1, double a2, double b1,
234 double b2) const {
235 345794 const auto absolute = [](double lo, double hi) {
236
4/4
✓ Branch 0 taken 178682 times.
✓ Branch 1 taken 167112 times.
✓ Branch 2 taken 11568 times.
✓ Branch 3 taken 167114 times.
345794 if (lo <= 0.0 && hi >= 0.0) {
237 11568 return std::array<double, 2>{0.0, std::max(-lo, hi)};
238 }
239 334226 return std::array<double, 2>{std::min(std::fabs(lo), std::fabs(hi)),
240 334226 std::max(std::fabs(lo), std::fabs(hi))};
241 };
242
1/1
✓ Branch 1 taken 172897 times.
172897 const auto [aLo, aHi] = absolute(a1, a2);
243
1/1
✓ Branch 1 taken 172897 times.
172897 const auto [bLo, bHi] = absolute(b1, b2);
244 172897 const std::array<double, 3> largest{1.0 / (1.0 + aLo * aLo + bLo * bLo),
245 172897 aHi * aHi / (1.0 + aHi * aHi + bLo * bLo),
246 172897 bHi * bHi / (1.0 + aLo * aLo + bHi * bHi)};
247 172897 const std::array<double, 3> smallest{1.0 / (1.0 + aHi * aHi + bHi * bHi),
248 172897 aLo * aLo / (1.0 + aLo * aLo + bHi * bHi),
249 172897 bLo * bLo / (1.0 + aHi * aHi + bLo * bLo)};
250 172897 const detail::CubeFace& f = detail::kCubeFaces[face];
251 172897 const std::array<Vec3d, 3> axes{f.normal, f.u, f.v};
252 172897 const Vec3d radii = surface_.Radii();
253 172897 const std::array<double, 3> r{radii.x, radii.y, radii.z};
254 172897 double numerator = 0.0;
255 172897 double denominator = 0.0;
256
2/2
✓ Branch 0 taken 518691 times.
✓ Branch 1 taken 172897 times.
691588 for (std::size_t k = 0; k < 3; ++k) {
257
4/4
✓ Branch 1 taken 345794 times.
✓ Branch 2 taken 172897 times.
✓ Branch 4 taken 172897 times.
✓ Branch 5 taken 172897 times.
518691 const std::size_t axis = axes[k].x != 0.0 ? 0 : (axes[k].y != 0.0 ? 1 : 2);
258 518691 const double inverse = 1.0 / (r[axis] * r[axis]);
259 518691 numerator += inverse * inverse * largest[k];
260 518691 denominator += inverse * smallest[k];
261 }
262 345794 return std::sqrt(numerator) / (denominator * denominator);
263 }
264
265 12401 SurfaceCells::Sampling CubeSphereCells::CellSampling(const CellIndex& cell, double cellSize) const {
266
2/2
✓ Branch 0 taken 11439 times.
✓ Branch 1 taken 962 times.
12401 if (weightCap_ == 0.0) {
267
1/1
✓ Branch 1 taken 11439 times.
11439 return {CellArea(cell, cellSize), 0.0};
268 }
269
1/1
✓ Branch 1 taken 962 times.
962 const auto [a1, a2, b1, b2] = Bounds(cell, cellSize);
270 962 const double solid = Corner(a2, b2) - Corner(a1, b2) - Corner(a2, b1) + Corner(a1, b1);
271
1/1
✓ Branch 1 taken 962 times.
962 const auto face = static_cast<std::size_t>(cell.i / CellsPerEdge(cellSize));
272 962 const double bound = surface_.Shape() != nullptr
273
2/2
✓ Branch 0 taken 97 times.
✓ Branch 1 taken 865 times.
1827 ? weightCap_
274
1/1
✓ Branch 1 taken 865 times.
865 : std::min(WeightBound(face, a1, a2, b1, b2), weightCap_);
275 962 return {solid * bound, bound};
276 }
277
278 414075 Vec2d CubeSphereCells::SamplePoint(double u, double v, const CellIndex& cell,
279 double cellSize) const {
280
1/1
✓ Branch 1 taken 414075 times.
414075 const std::int32_t n = CellsPerEdge(cellSize);
281 414075 const detail::CubeFace& face = detail::kCubeFaces[static_cast<std::size_t>(cell.i / n)];
282 414075 const std::int32_t col = cell.i % n;
283 414075 const double a1 = Edge(col, n);
284 414075 const double a2 = Edge(col + 1, n);
285 414075 const double b1 = Edge(cell.j, n);
286 414075 const double b2 = Edge(cell.j + 1, n);
287 26914875 const auto strip = [&](double a) {
288 26914875 return Corner(a, b2) - Corner(a, b1) - Corner(a1, b2) + Corner(a1, b1);
289 414075 };
290 414075 const double total = strip(a2);
291 414075 double lo = a1;
292 414075 double hi = a2;
293
2/2
✓ Branch 0 taken 26500800 times.
✓ Branch 1 taken 414075 times.
26914875 for (int iteration = 0; iteration < 64; ++iteration) {
294 26500800 const double mid = 0.5 * (lo + hi);
295
2/2
✓ Branch 1 taken 13336332 times.
✓ Branch 2 taken 13164468 times.
26500800 (strip(mid) < u * total ? lo : hi) = mid;
296 }
297 414075 const double a = 0.5 * (lo + hi);
298 1242225 const auto ratio = [&](double b) { return b / std::sqrt(1.0 + a * a + b * b); };
299 414075 const double k = ratio(b1) + v * (ratio(b2) - ratio(b1));
300 414075 const double b = k * std::sqrt(1.0 + a * a) / std::sqrt(1.0 - k * k);
301
2/2
✓ Branch 1 taken 414075 times.
✓ Branch 5 taken 414075 times.
414075 return Surface::Sphere(radius_).PointOf(detail::CubeFacePoint(face, a, b));
302 }
303
304 203984 Vec3d CubeSphereCells::Direction(int face, double alpha, double beta) {
305 const detail::CubeFace& f =
306
1/1
✓ Branch 1 taken 203984 times.
203984 detail::kCubeFaces[static_cast<std::size_t>(std::clamp(face, 0, 5))];
307 203984 const Vec3d p = detail::CubeFacePoint(f, std::tan(alpha), std::tan(beta));
308 203984 const double length = std::sqrt(p.x * p.x + p.y * p.y + p.z * p.z);
309 203984 return {p.x / length, p.y / length, p.z / length};
310 }
311
312 } // namespace terraingen
313