TerrainGen coverage

terraingen/
Coverage:
low: ≥ 0%
medium: ≥ 75.0%
high: ≥ 90.0%
Lines:
50 / 50
100.0%
Branches:
51 / 60
85.0%
Functions:
14 / 14, 1 excl
100.0%
util/surface.h
Line Branch Exec Source
1 /// @file
2 /// @brief `Surface`, the plane, sphere or triaxial ellipsoid a terrain lies on, with local frames
3 /// and distances.
4 #pragma once
5
6 #include "terraingen/util/geometry.h"
7 #include "terraingen/util/radial_shape.h"
8
9 #include <memory>
10
11 namespace terraingen {
12
13 /// The east, north and up unit vectors at a surface point.
14 struct LocalFrame {
15 Vec3d east; ///< Unit vector east.
16 Vec3d north; ///< Unit vector north.
17 Vec3d up; ///< Unit vector up: the surface normal, away from the centre of a body.
18 };
19
20 /// The integer coordinates of a cube of a 3D grid.
21 struct BucketIndex {
22 1 long long i = 0;
23 1 long long j = 0;
24 1 long long k = 0;
25
26 /// Componentwise equality.
27
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 time.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 time.
 1 bool operator==(const BucketIndex&) const = default;
28 };
29
30 /// The shape a terrain lies on, a plane or a body (a sphere, a triaxial ellipsoid of which the
31 /// sphere is the special case of equal radii, or a star-shaped shape model), measuring distances
32 /// and directions between surface coordinates.
33 ///
34 /// On a plane, surface coordinates are metres in a projected frame. On a body they are
35 /// planetocentric longitude and latitude in degrees (east and north positive): the direction from
36 /// the centre, with heights along it. On a sphere of radius R the local frame of a point is its
37 /// azimuthal equidistant projection in metres; on an ellipsoid or shape model it is the tangent
38 /// plane at the point, onto which offsets are projected along the normal, exact for features small
39 /// next to the local curvature radius.
40 class Surface {
41 public:
42 /// The kind of surface.
43 enum class Kind {
44 Plane, ///< A plane.
45 Body, ///< A body: a sphere, a triaxial ellipsoid or a shape model.
46 };
47
48 /// A plane, the default.
49 1120853 Surface() = default;
50
51 /// A plane.
52 490 static Surface Plane() {
53 490 return {};
54 }
55
56 /// A sphere of `radius` metres.
57 /// @throws std::invalid_argument unless `radius` is positive and finite.
58 static Surface Sphere(double radius);
59
60 /// A triaxial ellipsoid with semi-axes `a`, `b` and `c` metres along x (longitude 0), y
61 /// (longitude 90) and z (the north pole); equal radii give `Sphere`.
62 /// @throws std::invalid_argument unless every radius is positive and finite.
63 static Surface Ellipsoid(double a, double b, double c);
64
65 /// A star-shaped body whose radius by direction is `shape`, such as an asteroid from a shape
66 /// model; up is the normal of the interpolated surface, from finite differences.
67 /// @throws std::invalid_argument for a null `shape`.
68 static Surface ShapeModel(std::shared_ptr<const RadialShape> shape);
69
70 /// The kind of surface.
71 Kind GetKind() const {
72 return kind_;
73 }
74
75 /// Whether the surface is a body (a sphere, an ellipsoid or a shape model) rather than a plane.
76 136012211 bool IsBody() const {
77 136012211 return kind_ == Kind::Body;
78 }
79
80 /// Whether the surface is exactly a sphere, of one radius.
81 10481846 bool IsSphere() const {
82
8/8
✓ Branch 0 taken 10481845 times.
✓ Branch 1 taken 1 time.
✓ Branch 3 taken 8669703 times.
✓ Branch 4 taken 1812142 times.
✓ Branch 5 taken 6190636 times.
✓ Branch 6 taken 2479067 times.
✓ Branch 7 taken 3026472 times.
✓ Branch 8 taken 3164164 times.
 10481846 return kind_ == Kind::Body && !shape_ && radii_.x == radii_.y && radii_.y == radii_.z;
83 }
84
85 /// The body's shape model; null for a plane, a sphere or an ellipsoid.
86 1284 const RadialShape* Shape() const {
87 1284 return shape_.get();
88 }
89
90 /// The smallest distance in metres from the centre to the surface; 0 on a plane.
91 double MinRadius() const;
92
93 /// The largest distance in metres from the centre to the surface; 0 on a plane.
94 double MaxRadius() const;
95
96 /// The radius in metres of a sphere, or the volumetric mean radius of an ellipsoid
97 /// ((abc)^(1/3)) or shape model; 0 for a plane.
98 382 double Radius() const {
99 382 return radius_;
100 }
101
102 /// The semi-axes in metres along x, y and z; all equal on a sphere, 0 on a plane or a shape
103 /// model.
104 172902 Vec3d Radii() const {
105 172902 return radii_;
106 }
107
108 /// The distance in metres from the centre to the surface at `point`: the sphere's radius, or
109 /// the ellipsoid's radius in the point's direction; 0 on a plane.
110 double RadiusAt(const Vec2d& point) const;
111
112 /// The surface area per steradian of direction from the centre at `point`, r^2 / (u . n) in
113 /// m^2: R^2 on a sphere; 1 on a plane, where it is the area per square metre.
114 double AreaPerSteradian(const Vec2d& point) const;
115
116 /// Metres per surface unit along a meridian: 1 on a plane, R pi / 180 on a sphere, and that of
117 /// the mean radius `Radius` on an ellipsoid, where the true scale varies with position.
118 19076 double UnitLength() const {
119
2/2
✓ Branch 0 taken 4345 times.
✓ Branch 1 taken 14731 times.
 19076 return kind_ == Kind::Plane ? 1.0 : unitLength_;
120 }
121
122 /// The east and north offsets in metres of `point` in the local frame of `centre` (both in
123 /// surface coordinates); exactly `point` - `centre` on a plane.
124 36596549 Vec2d Offset(const Vec2d& centre, const Vec2d& point) const {
125
2/2
✓ Branch 0 taken 35205687 times.
✓ Branch 1 taken 1390862 times.
 36596549 if (kind_ == Kind::Plane) {
126 35205687 return {point.x - centre.x, point.y - centre.y};
127 }
128
2/2
✓ Branch 1 taken 45133 times.
✓ Branch 2 taken 1345729 times.
 1390862 return IsNonSphere() ? TangentOffset(centre, point) : SphereOffset(centre, point);
129 }
130
131 /// The point, in surface coordinates, at east and north `offset` metres from `centre`: the
132 /// inverse of `Offset`. On a sphere the longitude is in (-180, 180].
133 41199409 Vec2d Locate(const Vec2d& centre, const Vec2d& offset) const {
134
2/2
✓ Branch 0 taken 41148000 times.
✓ Branch 1 taken 51409 times.
 41199409 if (kind_ == Kind::Plane) {
135 41148000 return {centre.x + offset.x, centre.y + offset.y};
136 }
137
2/2
✓ Branch 1 taken 402 times.
✓ Branch 2 taken 51007 times.
 51409 if (shape_) {
138 402 return ShapeLocate(centre, offset);
139 }
140
2/2
✓ Branch 1 taken 602 times.
✓ Branch 2 taken 50405 times.
 51007 return IsNonSphere() ? EllipsoidLocate(centre, offset) : SphereLocate(centre, offset);
141 }
142
143 /// The distance in metres between two points: straight on a plane, along a great circle on
144 /// a sphere, and the chord on an ellipsoid, which matches the surface distance to second order
145 /// for nearby points and never exceeds it.
146 double Distance(const Vec2d& a, const Vec2d& b) const;
147
148 /// A surface point in 3D, in metres: (x, y, 0) on a plane; `RadiusAt` times `UnitVector` on a
149 /// body.
150 Vec3d Position(const Vec2d& point) const;
151
152 /// The east, north and up unit vectors at a point in the frame of `Position`: the x, y
153 /// and z axes on a plane; the local frame of `Offset` on a body, with up along the normal.
154 LocalFrame Frame(const Vec2d& point) const;
155
156 /// The surface point of a nonzero `direction` from the body's centre; (x, y) on a plane.
157 Vec2d PointOf(const Vec3d& direction) const;
158
159 /// A box containing every point within a distance of a region; on a body latitudes
160 /// are clamped to [-90, 90] and a box reaching a pole or wider than a full turn takes every
161 /// longitude, [-180, 180]. Longitudes are not wrapped.
162 /// @param metres The distance in metres.
163 Bounds2d Expand(const Bounds2d& region, double metres) const;
164
165 /// Whether a region holds a point, half-open as `Bounds2d::Contains`; on a sphere a
166 /// longitude also counts 360 degrees either side, and a box reaching latitude 90 holds the
167 /// north pole.
168 1524650273 bool Contains(const Bounds2d& region, const Vec2d& point) const {
169
2/2
✓ Branch 1 taken 22371801 times.
✓ Branch 2 taken 1502278472 times.
 1524650273 if (region.Contains(point.x, point.y)) {
170 22371801 return true;
171 }
172
2/2
✓ Branch 0 taken 1501362512 times.
✓ Branch 1 taken 915960 times.
 1502278472 if (kind_ == Kind::Plane) {
173 1501362512 return false;
174 }
175
4/4
✓ Branch 0 taken 96 times.
✓ Branch 1 taken 915864 times.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 91 times.
 915960 const Vec2d pole{point.x, point.y >= 90.0 && region.maxY >= 90.0 ? region.minY : point.y};
176
6/6
✓ Branch 1 taken 915955 times.
✓ Branch 2 taken 5 times.
✓ Branch 4 taken 914980 times.
✓ Branch 5 taken 975 times.
✓ Branch 6 taken 1406 times.
✓ Branch 7 taken 913574 times.
 1830940 return region.Contains(pole.x, pole.y) || region.Contains(pole.x - 360.0, pole.y) ||
177 1830940 region.Contains(pole.x + 360.0, pole.y);
178 }
179
180 /// The 3D bucket holding `Position(point)`, for hashing nearby features; on a plane k is
181 /// always 0.
182 /// @param size The bucket side in metres.
183 BucketIndex Bucket(const Vec2d& point, double size) const;
184
185 /// Whether two surfaces are the same kind and radii, and share one shape model or have none.
186
3/4
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 101 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 101 times.
 205 bool operator==(const Surface&) const = default;
187
188 private:
189 /// Whether the surface is a body the sphere formulas do not apply to: an ellipsoid of unequal
190 /// radii or a shape model.
191 5298545 bool IsNonSphere() const {
192
3/4
✓ Branch 0 taken 5298545 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 2728575 times.
✓ Branch 4 taken 2569970 times.
 5298545 return kind_ == Kind::Body && !IsSphere();
193 }
194
195 Vec2d SphereOffset(const Vec2d& centre, const Vec2d& point) const;
196 Vec2d SphereLocate(const Vec2d& centre, const Vec2d& offset) const;
197 Vec2d TangentOffset(const Vec2d& centre, const Vec2d& point) const;
198 Vec2d EllipsoidLocate(const Vec2d& centre, const Vec2d& offset) const;
199 Vec2d ShapeLocate(const Vec2d& centre, const Vec2d& offset) const;
200 /// The outward unit normal at surface coordinates `point`, whose position is `p`, on an
201 /// ellipsoid or shape model.
202 Vec3d Normal(const Vec2d& point, const Vec3d& p) const;
203
204
3/4
terraingen::Surface::IsNonSphere() const:
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 58 times.
terraingen::Surface::operator==(terraingen::Surface const&) const:
✗ Branch 0 not taken.
✓ Branch 1 taken 45 times.
 104 Kind kind_ = Kind::Plane;
205
3/4
terraingen::Surface::IsNonSphere() const:
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 57 times.
terraingen::Surface::operator==(terraingen::Surface const&) const:
✗ Branch 0 not taken.
✓ Branch 1 taken 45 times.
 103 double radius_ = 0.0;
206
2/4
terraingen::Surface::IsNonSphere() const:
✗ Branch 0 not taken.
✓ Branch 1 taken 57 times.
terraingen::Surface::operator==(terraingen::Surface const&) const:
✗ Branch 0 not taken.
✓ Branch 1 taken 45 times.
 102 double unitLength_ = 1.0;
207 102 Vec3d radii_{};
208 101 std::shared_ptr<const RadialShape> shape_;
209 };
210
211 /// The unit vector of a longitude and latitude: x towards longitude 0 on the equator, y
212 /// towards longitude 90 and z towards the north pole.
213 /// Both angles are in degrees.
214 Vec3d UnitVector(double lonDeg, double latDeg);
215
216 } // namespace terraingen
217