TerrainGen coverage

terraingen/
Coverage:
low: ≥ 0%
medium: ≥ 75.0%
high: ≥ 90.0%
Lines:
245 / 247
99.2%
Branches:
121 / 128
94.5%
Functions:
24 / 24
100.0%
util/surface.cpp
Line Branch Exec Source
1 #include "terraingen/util/surface.h"
2
3 #include "terraingen/util/checks.h"
4
5 #include <algorithm>
6 #include <cmath>
7 #include <memory>
8 #include <numbers>
9 #include <stdexcept>
10 #include <utility>
11
12 namespace terraingen {
13
14 namespace {
15
16 /// The east and north unit vectors at longitude `lon` and latitude `lat` (radians).
17 1978597 void LocalAxes(double lon, double lat, Vec3d& east, Vec3d& north) {
18 1978597 const double sinLon = std::sin(lon);
19 1978597 const double cosLon = std::cos(lon);
20 1978597 const double sinLat = std::sin(lat);
21 1978597 east = {-sinLon, cosLon, 0.0};
22 1978597 north = {-sinLat * cosLon, -sinLat * sinLon, std::cos(lat)};
23 1978597 }
24
25 } // namespace
26
27 15016452 Vec3d UnitVector(double lonDeg, double latDeg) {
28 15016452 const double lon = lonDeg * kDegree;
29 15016452 const double lat = latDeg * kDegree;
30 15016452 const double cosLat = std::cos(lat);
31 15016452 return {cosLat * std::cos(lon), cosLat * std::sin(lon), std::sin(lat)};
32 }
33
34 539312 Surface Surface::Sphere(double radius) {
35
1/1
✓ Branch 2 taken 539309 times.
 539312 RequirePositive(radius, "A sphere's radius");
36 539309 Surface surface;
37 539309 surface.kind_ = Kind::Body;
38 539309 surface.radius_ = radius;
39 539309 surface.unitLength_ = radius * kDegree;
40 539309 surface.radii_ = {radius, radius, radius};
41 539309 return surface;
42 }
43
44 94 Surface Surface::Ellipsoid(double a, double b, double c) {
45
1/1
✓ Branch 2 taken 93 times.
 94 RequirePositive(a, "An ellipsoid's radius a");
46
1/1
✓ Branch 2 taken 91 times.
 93 RequirePositive(b, "An ellipsoid's radius b");
47
1/1
✓ Branch 2 taken 90 times.
 91 RequirePositive(c, "An ellipsoid's radius c");
48
4/4
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 84 times.
✓ Branch 2 taken 3 times.
✓ Branch 3 taken 3 times.
 90 if (a == b && b == c) {
49
1/1
✓ Branch 1 taken 3 times.
 3 return Sphere(a);
50 }
51 87 Surface surface;
52 87 surface.kind_ = Kind::Body;
53 87 surface.radius_ = std::cbrt(a * b * c);
54 87 surface.unitLength_ = surface.radius_ * kDegree;
55 87 surface.radii_ = {a, b, c};
56 87 return surface;
57 87 }
58
59 75 Surface Surface::ShapeModel(std::shared_ptr<const RadialShape> shape) {
60
2/2
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 74 times.
 75 if (!shape) {
61
1/1
✓ Branch 2 taken 1 time.
 1 throw std::invalid_argument("A shape model surface needs a radial shape");
62 }
63 74 constexpr int kSteps = 180;
64 74 double cubed = 0.0;
65 74 double weight = 0.0;
66
2/2
✓ Branch 0 taken 13320 times.
✓ Branch 1 taken 74 times.
 13394 for (int j = 0; j < kSteps; ++j) {
67 13320 const double lat = -90.0 + (j + 0.5) * 180.0 / kSteps;
68 13320 const double cosLat = std::cos(lat * kDegree);
69
2/2
✓ Branch 0 taken 4795200 times.
✓ Branch 1 taken 13320 times.
 4808520 for (int i = 0; i < 2 * kSteps; ++i) {
70
1/1
✓ Branch 3 taken 4795200 times.
 4795200 const double r = shape->At(UnitVector(-180.0 + (i + 0.5) * 180.0 / kSteps, lat));
71 4795200 cubed += r * r * r * cosLat;
72 4795200 weight += cosLat;
73 }
74 }
75 74 Surface surface;
76 74 surface.kind_ = Kind::Body;
77 74 surface.radius_ = std::cbrt(cubed / weight);
78 74 surface.unitLength_ = surface.radius_ * kDegree;
79 74 surface.shape_ = std::move(shape);
80 74 return surface;
81 }
82
83 7416 double Surface::MinRadius() const {
84
2/2
✓ Branch 1 taken 9 times.
✓ Branch 2 taken 7407 times.
 7416 if (shape_) {
85 9 return shape_->MinRadius();
86 }
87
1/1
✓ Branch 1 taken 7407 times.
 7407 return std::min({radii_.x, radii_.y, radii_.z});
88 }
89
90 4 double Surface::MaxRadius() const {
91
2/2
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 1 time.
 4 if (shape_) {
92 3 return shape_->MaxRadius();
93 }
94
1/1
✓ Branch 1 taken 1 time.
 1 return std::max({radii_.x, radii_.y, radii_.z});
95 }
96
97 2759907 double Surface::RadiusAt(const Vec2d& point) const {
98
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 2759906 times.
 2759907 if (kind_ == Kind::Plane) {
99 1 return 0.0;
100 }
101
2/2
✓ Branch 1 taken 178339 times.
✓ Branch 2 taken 2581567 times.
 2759906 if (IsSphere()) {
102 178339 return radius_;
103 }
104 2581567 const Vec3d u = UnitVector(point.x, point.y);
105
2/2
✓ Branch 1 taken 617818 times.
✓ Branch 2 taken 1963749 times.
 2581567 if (shape_) {
106
1/1
✓ Branch 2 taken 617818 times.
 617818 return shape_->At(u);
107 }
108 1963749 const double ex = u.x / radii_.x;
109 1963749 const double ey = u.y / radii_.y;
110 1963749 const double ez = u.z / radii_.z;
111 1963749 return 1.0 / std::sqrt(ex * ex + ey * ey + ez * ez);
112 }
113
114 2139625 double Surface::AreaPerSteradian(const Vec2d& point) const {
115
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 2139624 times.
 2139625 if (kind_ == Kind::Plane) {
116 1 return 1.0;
117 }
118
2/2
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 2139623 times.
 2139624 if (IsSphere()) {
119 1 return radius_ * radius_;
120 }
121
1/1
✓ Branch 1 taken 2139623 times.
 2139623 const Vec3d p = Position(point);
122 2139623 const double r = Length(p);
123
1/1
✓ Branch 1 taken 2139623 times.
 2139623 const double cosine = Dot(p, Normal(point, p)) / r;
124 2139623 return r * r / cosine;
125 }
126
127 2185762 Vec3d Surface::Normal(const Vec2d& point, const Vec3d& p) const {
128 2185762 Vec3d n;
129
2/2
✓ Branch 1 taken 536205 times.
✓ Branch 2 taken 1649557 times.
 2185762 if (shape_) {
130 536205 constexpr double kStep = 1e-6;
131 536205 const Vec3d u = UnitVector(point.x, point.y);
132 536205 Vec3d east;
133 536205 Vec3d north;
134 536205 LocalAxes(point.x * kDegree, point.y * kDegree, east, north);
135 2144820 const auto at = [&](const Vec3d& axis, double sign) {
136 2144820 const Vec3d d{u.x + sign * kStep * axis.x, u.y + sign * kStep * axis.y,
137 2144820 u.z + sign * kStep * axis.z};
138 2144820 const double length = Length(d);
139 2144820 const Vec3d unit{d.x / length, d.y / length, d.z / length};
140
1/1
✓ Branch 2 taken 2144820 times.
 4289640 return unit * shape_->At(unit);
141 536205 };
142
4/4
✓ Branch 1 taken 536205 times.
✓ Branch 4 taken 536205 times.
✓ Branch 8 taken 536205 times.
✓ Branch 11 taken 536205 times.
 536205 n = Cross(at(east, 1.0) - at(east, -1.0), at(north, 1.0) - at(north, -1.0));
143 } else {
144 1649557 n = {p.x / (radii_.x * radii_.x), p.y / (radii_.y * radii_.y), p.z / (radii_.z * radii_.z)};
145 }
146 2185762 const double length = Length(n);
147 2185762 return {n.x / length, n.y / length, n.z / length};
148 }
149
150 45133 Vec2d Surface::TangentOffset(const Vec2d& centre, const Vec2d& point) const {
151
1/1
✓ Branch 1 taken 45133 times.
 45133 const Vec3d c = Position(centre);
152
1/1
✓ Branch 1 taken 45133 times.
 45133 const Vec3d p = Position(point);
153
1/1
✓ Branch 1 taken 45133 times.
 45133 const LocalFrame frame = Frame(centre);
154 45133 const Vec3d d{p.x - c.x, p.y - c.y, p.z - c.z};
155 90266 return {Dot(d, frame.east), Dot(d, frame.north)};
156 }
157
158 402 Vec2d Surface::ShapeLocate(const Vec2d& centre, const Vec2d& offset) const {
159
4/4
✓ Branch 0 taken 20 times.
✓ Branch 1 taken 382 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 18 times.
 402 if (offset.x == 0.0 && offset.y == 0.0) {
160 2 return centre;
161 }
162
1/1
✓ Branch 1 taken 400 times.
 400 const Vec3d c = Position(centre);
163
1/1
✓ Branch 1 taken 400 times.
 400 const LocalFrame frame = Frame(centre);
164 400 const Vec3d q = c + frame.east * offset.x + frame.north * offset.y;
165 41211 const auto along = [&](double s) { return q - frame.up * s; };
166 40811 const auto excess = [&](double s) {
167 40811 const Vec3d x = along(s);
168 40811 const double r = Length(x);
169
1/1
✓ Branch 2 taken 40811 times.
 40811 return r - shape_->At({x.x / r, x.y / r, x.z / r});
170 400 };
171 400 const double step = shape_->MaxRadius() / 64.0;
172 400 double inside = 0.0;
173 400 double outside = 0.0;
174
3/3
✓ Branch 1 taken 400 times.
✓ Branch 3 taken 374 times.
✓ Branch 4 taken 26 times.
 400 if (excess(0.0) >= 0.0) {
175 374 inside = step;
176
6/7
✓ Branch 1 taken 385 times.
✓ Branch 3 taken 11 times.
✓ Branch 4 taken 374 times.
✓ Branch 7 taken 11 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 11 times.
✓ Branch 10 taken 374 times.
 385 while (excess(inside) >= 0.0 && inside < 4.0 * shape_->MaxRadius()) {
177 11 outside = inside;
178 11 inside += step;
179 }
180 } else {
181 26 outside = -step;
182
3/7
✓ Branch 1 taken 26 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 26 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 26 times.
 26 while (excess(outside) < 0.0 && outside > -4.0 * shape_->MaxRadius()) {
183 ✗ inside = outside;
184 ✗ outside -= step;
185 }
186 }
187
2/2
✓ Branch 0 taken 40000 times.
✓ Branch 1 taken 400 times.
 40400 for (int iteration = 0; iteration < 100; ++iteration) {
188 40000 const double middle = 0.5 * (inside + outside);
189
3/3
✓ Branch 1 taken 40000 times.
✓ Branch 3 taken 21005 times.
✓ Branch 4 taken 18995 times.
 40000 (excess(middle) >= 0.0 ? outside : inside) = middle;
190 }
191 400 return PointOf(along(0.5 * (inside + outside)));
192 }
193
194 602 Vec2d Surface::EllipsoidLocate(const Vec2d& centre, const Vec2d& offset) const {
195
4/4
✓ Branch 0 taken 20 times.
✓ Branch 1 taken 582 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 18 times.
 602 if (offset.x == 0.0 && offset.y == 0.0) {
196 2 return centre;
197 }
198
1/1
✓ Branch 1 taken 600 times.
 600 const Vec3d c = Position(centre);
199
1/1
✓ Branch 1 taken 600 times.
 600 const LocalFrame frame = Frame(centre);
200 600 const Vec3d q{c.x + offset.x * frame.east.x + offset.y * frame.north.x,
201 600 c.y + offset.x * frame.east.y + offset.y * frame.north.y,
202 600 c.z + offset.x * frame.east.z + offset.y * frame.north.z};
203 600 const Vec3d& n = frame.up;
204 600 const Vec3d inverse{1.0 / (radii_.x * radii_.x), 1.0 / (radii_.y * radii_.y),
205 600 1.0 / (radii_.z * radii_.z)};
206 600 const double a = n.x * n.x * inverse.x + n.y * n.y * inverse.y + n.z * n.z * inverse.z;
207 600 const double b = q.x * n.x * inverse.x + q.y * n.y * inverse.y + q.z * n.z * inverse.z;
208 600 const double cq = q.x * q.x * inverse.x + q.y * q.y * inverse.y + q.z * q.z * inverse.z;
209 600 const double discriminant = b * b - a * (cq - 1.0);
210 600 const double t = (b - std::sqrt(std::max(discriminant, 0.0))) / a;
211 600 return PointOf({q.x - t * n.x, q.y - t * n.y, q.z - t * n.z});
212 }
213
214 1345729 Vec2d Surface::SphereOffset(const Vec2d& centre, const Vec2d& point) const {
215 1345729 const Vec3d c = UnitVector(centre.x, centre.y);
216 1345729 const Vec3d p = UnitVector(point.x, point.y);
217 1345729 Vec3d east;
218 1345729 Vec3d north;
219 1345729 LocalAxes(centre.x * kDegree, centre.y * kDegree, east, north);
220 1345729 const double angle = std::atan2(Length(Cross(c, p)), Dot(c, p));
221 1345729 const double e = Dot(p, east);
222 1345729 const double n = Dot(p, north);
223 1345729 const double across = std::hypot(e, n);
224
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 1345723 times.
 1345729 if (across == 0.0) {
225
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
 6 return {0.0, angle == 0.0 ? 0.0 : radius_ * angle};
226 }
227 1345723 const double scale = radius_ * angle / across;
228 1345723 return {e * scale, n * scale};
229 }
230
231 50406 Vec2d Surface::SphereLocate(const Vec2d& centre, const Vec2d& offset) const {
232 50406 const double distance = std::hypot(offset.x, offset.y);
233
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 50396 times.
 50406 if (distance == 0.0) {
234 10 return centre;
235 }
236 50396 const Vec3d c = UnitVector(centre.x, centre.y);
237 50396 Vec3d east;
238 50396 Vec3d north;
239 50396 LocalAxes(centre.x * kDegree, centre.y * kDegree, east, north);
240 50396 const double angle = distance / radius_;
241 50396 const double along = std::sin(angle) / distance;
242 50396 const double cosAngle = std::cos(angle);
243 50396 const Vec3d p{cosAngle * c.x + along * (offset.x * east.x + offset.y * north.x),
244 50396 cosAngle * c.y + along * (offset.x * east.y + offset.y * north.y),
245 50396 cosAngle * c.z + along * (offset.x * east.z + offset.y * north.z)};
246 50396 return PointOf(p);
247 }
248
249 1316487 double Surface::Distance(const Vec2d& a, const Vec2d& b) const {
250
2/2
✓ Branch 0 taken 1213096 times.
✓ Branch 1 taken 103391 times.
 1316487 if (kind_ == Kind::Plane) {
251 1213096 return std::hypot(b.x - a.x, b.y - a.y);
252 }
253
2/2
✓ Branch 1 taken 100401 times.
✓ Branch 2 taken 2990 times.
 103391 if (IsNonSphere()) {
254
1/1
✓ Branch 1 taken 100401 times.
 100401 const Vec3d p = Position(a);
255
1/1
✓ Branch 1 taken 100401 times.
 100401 const Vec3d q = Position(b);
256 100401 return Length(Vec3d{q.x - p.x, q.y - p.y, q.z - p.z});
257 }
258 2990 const Vec3d u = UnitVector(a.x, a.y);
259 2990 const Vec3d v = UnitVector(b.x, b.y);
260 2990 return radius_ * std::atan2(Length(Cross(u, v)), Dot(u, v));
261 }
262
263 3884841 Vec3d Surface::Position(const Vec2d& point) const {
264
2/2
✓ Branch 0 taken 88916 times.
✓ Branch 1 taken 3795925 times.
 3884841 if (kind_ == Kind::Plane) {
265 88916 return {point.x, point.y, 0.0};
266 }
267 3795925 const Vec3d unit = UnitVector(point.x, point.y);
268
3/3
✓ Branch 1 taken 2528942 times.
✓ Branch 2 taken 1266983 times.
✓ Branch 4 taken 2528942 times.
 3795925 const double r = IsNonSphere() ? RadiusAt(point) : radius_;
269 3795925 return {r * unit.x, r * unit.y, r * unit.z};
270 }
271
272 46274 LocalFrame Surface::Frame(const Vec2d& point) const {
273
2/2
✓ Branch 0 taken 7 times.
✓ Branch 1 taken 46267 times.
 46274 if (kind_ == Kind::Plane) {
274 7 return {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}};
275 }
276 46267 LocalFrame frame;
277 46267 LocalAxes(point.x * kDegree, point.y * kDegree, frame.east, frame.north);
278 46267 frame.up = UnitVector(point.x, point.y);
279
2/2
✓ Branch 1 taken 46139 times.
✓ Branch 2 taken 128 times.
 46267 if (IsNonSphere()) {
280
2/2
✓ Branch 1 taken 46139 times.
✓ Branch 4 taken 46139 times.
 46139 const Vec3d n = Normal(point, Position(point));
281 46139 const double along = Dot(frame.east, n);
282 46139 const Vec3d east{frame.east.x - along * n.x, frame.east.y - along * n.y,
283 46139 frame.east.z - along * n.z};
284 46139 const double length = Length(east);
285 46139 frame.up = n;
286 46139 frame.east = {east.x / length, east.y / length, east.z / length};
287 46139 frame.north = Cross(n, frame.east);
288 }
289 46267 return frame;
290 }
291
292 973068 Vec2d Surface::PointOf(const Vec3d& direction) const {
293
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 973067 times.
 973068 if (kind_ == Kind::Plane) {
294 1 return {direction.x, direction.y};
295 }
296 973067 const double lon = std::atan2(direction.y, direction.x) / kDegree;
297 973067 const double lat = std::atan2(direction.z, std::hypot(direction.x, direction.y)) / kDegree;
298
2/2
✓ Branch 0 taken 972784 times.
✓ Branch 1 taken 283 times.
 973067 return {lon == -180.0 ? 180.0 : lon, lat};
299 }
300
301 89594 BucketIndex Surface::Bucket(const Vec2d& point, double size) const {
302
1/1
✓ Branch 1 taken 89594 times.
 89594 const Vec3d p = Position(point);
303 89594 return {static_cast<long long>(std::floor(p.x / size)),
304 89594 static_cast<long long>(std::floor(p.y / size)),
305 179188 static_cast<long long>(std::floor(p.z / size))};
306 }
307
308 10894196 Bounds2d Surface::Expand(const Bounds2d& region, double metres) const {
309
2/2
✓ Branch 0 taken 10870919 times.
✓ Branch 1 taken 23277 times.
 10894196 if (kind_ == Kind::Plane) {
310 10870919 return region.Expanded(metres);
311 }
312 23277 double dLat = metres / unitLength_;
313 23277 double angle = metres / radius_;
314
2/2
✓ Branch 1 taken 7358 times.
✓ Branch 2 taken 15919 times.
 23277 if (IsNonSphere()) {
315
1/1
✓ Branch 1 taken 7358 times.
 7358 const double shortest = MinRadius();
316
2/2
✓ Branch 0 taken 7357 times.
✓ Branch 1 taken 1 time.
 7358 angle = metres < shortest ? std::asin(metres / shortest) : std::numbers::pi;
317 7358 dLat = angle / kDegree;
318 }
319 69831 Bounds2d expanded{region.minX, std::max(-90.0, region.minY - dLat), region.maxX,
320 23277 std::min(90.0, region.maxY + dLat)};
321 23277 const double widest = std::max(std::fabs(region.minY), std::fabs(region.maxY));
322 23277 const double reach = std::sin(std::min(angle, std::numbers::pi / 2.0));
323 23277 const double cosLat = std::cos(std::min(widest, 90.0) * kDegree);
324
2/2
✓ Branch 0 taken 21836 times.
✓ Branch 1 taken 1439 times.
 23275 const double dLon = angle < std::numbers::pi / 2.0 && reach < cosLat
325
2/2
✓ Branch 0 taken 23275 times.
✓ Branch 1 taken 2 times.
 46552 ? std::asin(reach / cosLat) / kDegree
326 23277 : 360.0;
327
5/6
✓ Branch 0 taken 22549 times.
✓ Branch 1 taken 728 times.
✓ Branch 2 taken 21737 times.
✓ Branch 3 taken 812 times.
✓ Branch 4 taken 21737 times.
✗ Branch 5 not taken.
 23277 if (expanded.minY <= -90.0 || expanded.maxY >= 90.0 || dLon >= 180.0 ||
328
2/2
✓ Branch 0 taken 3738 times.
✓ Branch 1 taken 17999 times.
 21737 region.maxX - region.minX + 2.0 * dLon >= 360.0) {
329 5278 expanded.minX = -180.0;
330 5278 expanded.maxX = 180.0;
331 } else {
332 17999 expanded.minX = region.minX - dLon;
333 17999 expanded.maxX = region.maxX + dLon;
334 }
335 23277 return expanded;
336 }
337
338 } // namespace terraingen
339