layers/boulders/boulder_models.cpp
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | #include "terraingen/layers/boulders/boulder_models.h" | ||
| 2 | |||
| 3 | #include <algorithm> | ||
| 4 | #include <array> | ||
| 5 | #include <cmath> | ||
| 6 | #include <limits> | ||
| 7 | #include <numbers> | ||
| 8 | #include <utility> | ||
| 9 | #include <vector> | ||
| 10 | |||
| 11 | namespace terraingen { | ||
| 12 | |||
| 13 | namespace { | ||
| 14 | |||
| 15 | /// Axis `index` of `frame`: 0 east, 1 north, 2 up. | ||
| 16 | 1026 | const Vec3d& Axis(const LocalFrame& frame, int index) { | |
| 17 | 1026 | return index == 0 ? frame.east : (index == 1 ? frame.north : frame.up); | |
| 18 | } | ||
| 19 | |||
| 20 | /// The rotation turning the shape's axis that `rest` names to the vertical. | ||
| 21 | 11083 | RigidTransform RestRotation(BoulderRest rest) { | |
| 22 | 11083 | switch (rest) { | |
| 23 | 68 | case BoulderRest::Side: | |
| 24 | 68 | return RigidTransform::RotationX(0.5 * std::numbers::pi); | |
| 25 | 41 | case BoulderRest::End: | |
| 26 | 41 | return RigidTransform::RotationY(-0.5 * std::numbers::pi); | |
| 27 | 10974 | case BoulderRest::Face: | |
| 28 | 10974 | break; | |
| 29 | } | ||
| 30 | 10974 | return {}; | |
| 31 | } | ||
| 32 | |||
| 33 | /// The centre of mass of a closed mesh of uniform density, from its signed tetrahedra. | ||
| 34 | 11016 | Vec3d MassCentre(const std::vector<Vec3d>& points, | |
| 35 | const std::vector<std::array<std::uint32_t, 3>>& triangles) { | ||
| 36 | 11016 | Vec3d sum{}; | |
| 37 | 11016 | double volume = 0.0; | |
| 38 | 10510962 | for (const auto& t : triangles) { | |
| 39 | 10499946 | const Vec3d& a = points[t[0]]; | |
| 40 | 10499946 | const Vec3d& b = points[t[1]]; | |
| 41 | 10499946 | const Vec3d& c = points[t[2]]; | |
| 42 | 10499946 | const double v = Dot(a, Cross(b, c)); | |
| 43 | 10499946 | sum = sum + (a + b + c) * v; | |
| 44 | 10499946 | volume += v; | |
| 45 | } | ||
| 46 | 11016 | return sum * (1.0 / (4.0 * volume)); | |
| 47 | } | ||
| 48 | |||
| 49 | /// The convex hull of points in the plane, counter-clockwise, without collinear points. | ||
| 50 | 79435 | std::vector<Vec2d> Hull2d(std::vector<Vec2d> points) { | |
| 51 | 79435 | std::sort(points.begin(), points.end(), [](const Vec2d& a, const Vec2d& b) { | |
| 52 | 186931 | return a.x < b.x || (a.x == b.x && a.y < b.y); | |
| 53 | }); | ||
| 54 | 79435 | if (points.size() < 3) { | |
| 55 | 36188 | return points; | |
| 56 | } | ||
| 57 | 86570 | const auto turn = [](const Vec2d& o, const Vec2d& a, const Vec2d& b) { | |
| 58 | 86570 | return (a.x - o.x) * (b.y - o.y) - (a.y - o.y) * (b.x - o.x); | |
| 59 | }; | ||
| 60 | 43247 | std::vector<Vec2d> hull(2 * points.size()); | |
| 61 | 43247 | std::size_t k = 0; | |
| 62 | 173010 | for (std::size_t i = 0; i < points.size(); ++i) { | |
| 63 | 151290 | while (k >= 2 && turn(hull[k - 2], hull[k - 1], points[i]) <= 0.0) { | |
| 64 | 21527 | --k; | |
| 65 | } | ||
| 66 | 129763 | hull[k++] = points[i]; | |
| 67 | } | ||
| 68 | 129763 | for (std::size_t i = points.size() - 1, lower = k + 1; i-- > 0;) { | |
| 69 | 108258 | while (k >= lower && turn(hull[k - 2], hull[k - 1], points[i]) <= 0.0) { | |
| 70 | 21742 | --k; | |
| 71 | } | ||
| 72 | 86516 | hull[k++] = points[i]; | |
| 73 | } | ||
| 74 | 43247 | hull.resize(k - 1); | |
| 75 | 43247 | return hull; | |
| 76 | 43247 | } | |
| 77 | |||
| 78 | /// The point of `hull`'s boundary nearest `p`, and whether `p` lies inside a hull of 3 or more | ||
| 79 | /// points. | ||
| 80 | 79435 | std::pair<Vec2d, bool> NearestOnHull(const std::vector<Vec2d>& hull, const Vec2d& p) { | |
| 81 | 79435 | if (hull.size() == 1) { | |
| 82 | 10994 | return {hull[0], false}; | |
| 83 | } | ||
| 84 | 68441 | bool inside = hull.size() >= 3; | |
| 85 | 68441 | Vec2d best = hull[0]; | |
| 86 | 68441 | double bestDistance = std::numeric_limits<double>::infinity(); | |
| 87 | 248592 | for (std::size_t i = 0; i < hull.size(); ++i) { | |
| 88 | 180151 | const Vec2d& a = hull[i]; | |
| 89 | 180151 | const Vec2d& b = hull[(i + 1) % hull.size()]; | |
| 90 | 180151 | const double ex = b.x - a.x; | |
| 91 | 180151 | const double ey = b.y - a.y; | |
| 92 | 180151 | if (ex * (p.y - a.y) - ey * (p.x - a.x) < 0.0) { | |
| 93 | 64929 | inside = false; | |
| 94 | } | ||
| 95 | const double t = | ||
| 96 | 180151 | std::clamp(((p.x - a.x) * ex + (p.y - a.y) * ey) / (ex * ex + ey * ey), 0.0, 1.0); | |
| 97 | 180151 | const Vec2d q{a.x + t * ex, a.y + t * ey}; | |
| 98 | 180151 | const double d = std::hypot(p.x - q.x, p.y - q.y); | |
| 99 | 180151 | if (d < bestDistance) { | |
| 100 | 105416 | bestDistance = d; | |
| 101 | 105416 | best = q; | |
| 102 | } | ||
| 103 | } | ||
| 104 | 68441 | return {best, inside}; | |
| 105 | } | ||
| 106 | |||
| 107 | /// The rotation that tips `points` over on a level plane until they rest: each step pivots about | ||
| 108 | /// the edge or corner of the ground contact nearest the centre of mass, toward it, until another | ||
| 109 | /// point touches, and stops once the centre of mass stands over the contact. | ||
| 110 | 11016 | std::array<double, 9> SettleRotation(std::vector<Vec3d> points, Vec3d centre) { | |
| 111 | 11016 | std::array<double, 9> total{1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}; | |
| 112 | 11016 | double extent = 0.0; | |
| 113 | 5283021 | for (const Vec3d& p : points) { | |
| 114 | extent = | ||
| 115 | 5272005 | std::fmax(extent, std::fmax(std::fabs(p.x), std::fmax(std::fabs(p.y), std::fabs(p.z)))); | |
| 116 | } | ||
| 117 | 11016 | const double tolerance = 1e-7 * extent; | |
| 118 | 79435 | for (int step = 0; step < 64; ++step) { | |
| 119 | 79435 | double low = std::numeric_limits<double>::infinity(); | |
| 120 | 41820565 | for (const Vec3d& p : points) { | |
| 121 | 41741130 | low = std::fmin(low, p.z); | |
| 122 | } | ||
| 123 | 79435 | std::vector<Vec2d> contact; | |
| 124 | 41820565 | for (const Vec3d& p : points) { | |
| 125 | 41741130 | if (p.z <= low + tolerance) { | |
| 126 | 191145 | contact.push_back({p.x, p.y}); | |
| 127 | } | ||
| 128 | } | ||
| 129 | 79435 | const Vec2d c{centre.x, centre.y}; | |
| 130 | 79435 | const auto [pivot, inside] = NearestOnHull(Hull2d(contact), c); | |
| 131 | 79435 | const double distance = std::hypot(c.x - pivot.x, c.y - pivot.y); | |
| 132 | 79435 | if (inside || distance <= tolerance) { | |
| 133 | break; | ||
| 134 | } | ||
| 135 | 68419 | const Vec3d d{(c.x - pivot.x) / distance, (c.y - pivot.y) / distance, 0.0}; | |
| 136 | 68419 | const Vec3d a{-d.y, d.x, 0.0}; | |
| 137 | 68419 | const Vec3d up{0.0, 0.0, 1.0}; | |
| 138 | 68419 | const Vec3d o{pivot.x, pivot.y, low}; | |
| 139 | 68419 | double rise = 0.0; | |
| 140 | 68419 | double run = 0.0; | |
| 141 | 36537544 | for (const Vec3d& p : points) { | |
| 142 | 36469125 | const double u = Dot(p - o, d); | |
| 143 | 36469125 | if (u > tolerance && (run == 0.0 || (p.z - low) * run < rise * u)) { | |
| 144 | 407406 | rise = p.z - low; | |
| 145 | 407406 | run = u; | |
| 146 | } | ||
| 147 | } | ||
| 148 | 68419 | if (run == 0.0) { | |
| 149 | ✗ | break; | |
| 150 | } | ||
| 151 | 68419 | const double angle = std::atan2(rise, run); | |
| 152 | 68419 | const double cs = std::cos(angle); | |
| 153 | 68419 | const double sn = std::sin(angle); | |
| 154 | 36537544 | const auto turn = [&](const Vec3d& p) { | |
| 155 | 36537544 | const Vec3d v = p - o; | |
| 156 | 36537544 | const double u = Dot(v, d); | |
| 157 | 36537544 | const double w = v.z; | |
| 158 | 36537544 | return o + a * Dot(v, a) + d * (u * cs + w * sn) + up * (w * cs - u * sn); | |
| 159 | 68419 | }; | |
| 160 | 36537544 | for (Vec3d& p : points) { | |
| 161 | 36469125 | p = turn(p); | |
| 162 | } | ||
| 163 | 68419 | centre = turn(centre); | |
| 164 | 68419 | std::array<double, 9> m{}; | |
| 165 | 68419 | const std::array<double, 3> av{a.x, a.y, a.z}; | |
| 166 | 68419 | const std::array<double, 3> dv{d.x, d.y, d.z}; | |
| 167 | 68419 | const std::array<double, 3> zv{0.0, 0.0, 1.0}; | |
| 168 | 273676 | for (std::size_t r = 0; r < 3; ++r) { | |
| 169 | 821028 | for (std::size_t k = 0; k < 3; ++k) { | |
| 170 | 1231542 | m[3 * r + k] = av[r] * av[k] + cs * (dv[r] * dv[k] + zv[r] * zv[k]) + | |
| 171 | 615771 | sn * (dv[r] * zv[k] - zv[r] * dv[k]); | |
| 172 | } | ||
| 173 | } | ||
| 174 | 68419 | std::array<double, 9> next{}; | |
| 175 | 273676 | for (std::size_t r = 0; r < 3; ++r) { | |
| 176 | 821028 | for (std::size_t k = 0; k < 3; ++k) { | |
| 177 | 2463084 | for (std::size_t j = 0; j < 3; ++j) { | |
| 178 | 1847313 | next[3 * r + k] += m[3 * r + j] * total[3 * j + k]; | |
| 179 | } | ||
| 180 | } | ||
| 181 | } | ||
| 182 | 68419 | total = next; | |
| 183 | 79435 | } | |
| 184 | 11016 | return total; | |
| 185 | } | ||
| 186 | |||
| 187 | } // namespace | ||
| 188 | |||
| 189 | 1617 | RigidTransform TiltRotation(double tilt, double azimuth) { | |
| 190 | 1617 | const double ax = -std::sin(azimuth); | |
| 191 | 1617 | const double ay = std::cos(azimuth); | |
| 192 | 1617 | const double c = std::cos(tilt); | |
| 193 | 1617 | const double s = std::sin(tilt); | |
| 194 | 1617 | const double t = 1.0 - c; | |
| 195 | 1617 | RigidTransform turn; | |
| 196 | 1617 | turn.rotation = {c + t * ax * ax, t * ax * ay, s * ay, t * ax * ay, c + t * ay * ay, -s * ax, | |
| 197 | 1617 | -s * ay, s * ax, c}; | |
| 198 | 1617 | return turn; | |
| 199 | } | ||
| 200 | |||
| 201 | 504 | BoulderPlacement PlacementAtSite(const BoulderPlacement& placement, const Vec2d& boulder, | |
| 202 | const Vec2d& site, const Surface& surface) { | ||
| 203 | 504 | if (!surface.IsBody()) { | |
| 204 | 447 | return placement; | |
| 205 | } | ||
| 206 | 57 | const LocalFrame from = surface.Frame(boulder); | |
| 207 | 57 | const LocalFrame to = surface.Frame(site); | |
| 208 | 57 | std::array<double, 9> local{}; | |
| 209 | 228 | for (int r = 0; r < 3; ++r) { | |
| 210 | 684 | for (int c = 0; c < 3; ++c) { | |
| 211 | 513 | local[static_cast<std::size_t>(3 * r + c)] = Dot(Axis(to, r), Axis(from, c)); | |
| 212 | } | ||
| 213 | } | ||
| 214 | 57 | const Vec3d pr = surface.Position(boulder); | |
| 215 | 57 | const Vec3d ps = surface.Position(site); | |
| 216 | 57 | const Vec3d d{pr.x - ps.x, pr.y - ps.y, pr.z - ps.z}; | |
| 217 | 57 | RigidTransform toSite; | |
| 218 | 57 | toSite.rotation = local; | |
| 219 | 57 | toSite.translation = {Dot(to.east, d), Dot(to.north, d), Dot(to.up, d)}; | |
| 220 | 57 | BoulderPlacement moved = placement; | |
| 221 | 57 | moved.pose = toSite.After(placement.pose); | |
| 222 | 57 | return moved; | |
| 223 | } | ||
| 224 | |||
| 225 | 11083 | BoulderPlacement PlaceBoulderModel(const Boulder& boulder, const BoulderShape& shape, | |
| 226 | double groundHeight, bool settle) { | ||
| 227 | 11083 | BoulderPlacement placement; | |
| 228 | 11083 | placement.scale = {boulder.diameter, boulder.diameter, boulder.height / shape.shortAxis}; | |
| 229 | 11083 | placement.pose = RigidTransform::RotationZ(boulder.yaw) | |
| 230 | 11083 | .After(RigidTransform::RotationX(boulder.lean)) | |
| 231 | 11083 | .After(RestRotation(boulder.rest)); | |
| 232 | 11083 | if (settle && !shape.mesh.triangles.empty()) { | |
| 233 | 11016 | std::vector<Vec3d> points; | |
| 234 | 11016 | points.reserve(shape.mesh.positions.size()); | |
| 235 | 5283021 | for (const Vec3d& p : shape.mesh.positions) { | |
| 236 | 5272005 | points.push_back(placement.Apply(p)); | |
| 237 | } | ||
| 238 | 11016 | RigidTransform rest; | |
| 239 | 11016 | rest.rotation = SettleRotation(points, MassCentre(points, shape.mesh.triangles)); | |
| 240 | 11016 | placement.pose = rest.After(placement.pose); | |
| 241 | 11016 | } | |
| 242 | 11083 | if (boulder.tilt != 0.0) { | |
| 243 | 100 | placement.pose = TiltRotation(boulder.tilt, boulder.tiltAzimuth).After(placement.pose); | |
| 244 | } | ||
| 245 | 11083 | double low = std::numeric_limits<double>::infinity(); | |
| 246 | 11083 | double high = -low; | |
| 247 | 5317078 | for (const Vec3d& p : shape.mesh.positions) { | |
| 248 | 5305995 | const double z = placement.Apply(p).z; | |
| 249 | 5305995 | low = std::fmin(low, z); | |
| 250 | 5305995 | high = std::fmax(high, z); | |
| 251 | } | ||
| 252 | 11083 | placement.pose.translation = {boulder.x, boulder.y, | |
| 253 | 11083 | groundHeight - boulder.burial * (high - low) - low}; | |
| 254 | 11083 | return placement; | |
| 255 | } | ||
| 256 | |||
| 257 | 190183 | std::optional<double> TopSurface(const TriangleMesh& placed, double x, double y) { | |
| 258 | 190183 | std::optional<double> top; | |
| 259 | 189739647 | for (const auto& t : placed.triangles) { | |
| 260 | 189549464 | const Vec3d& a = placed.positions[t[0]]; | |
| 261 | 189549464 | const Vec3d& b = placed.positions[t[1]]; | |
| 262 | 189549464 | const Vec3d& c = placed.positions[t[2]]; | |
| 263 | 189549464 | const double area = (b.x - a.x) * (c.y - a.y) - (c.x - a.x) * (b.y - a.y); | |
| 264 | 189549464 | if (!(area > 0.0)) { | |
| 265 | 189499045 | continue; | |
| 266 | } | ||
| 267 | 99624970 | const double u = ((b.x - x) * (c.y - y) - (c.x - x) * (b.y - y)) / area; | |
| 268 | 99624970 | const double v = ((c.x - x) * (a.y - y) - (a.x - x) * (c.y - y)) / area; | |
| 269 | 99624970 | const double w = 1.0 - u - v; | |
| 270 | 99624970 | if (u < 0.0 || v < 0.0 || w < 0.0) { | |
| 271 | 99574551 | continue; | |
| 272 | } | ||
| 273 | 50419 | const double z = u * a.z + v * b.z + w * c.z; | |
| 274 | 50419 | if (!top || z > *top) { | |
| 275 | 50383 | top = z; | |
| 276 | } | ||
| 277 | } | ||
| 278 | 190183 | return top; | |
| 279 | } | ||
| 280 | |||
| 281 | 5 | TriangleMesh BoulderMesh(const Boulder& boulder, const BoulderShapeOptions& options, | |
| 282 | double groundHeight, std::uint32_t prototypes) { | ||
| 283 | const BoulderShape shape = | ||
| 284 | 5 | GenerateBoulderShape(BoulderShapeSeed(boulder.id, prototypes), options); | |
| 285 | 5 | const BoulderPlacement placement = PlaceBoulderModel(boulder, shape, groundHeight); | |
| 286 | 5 | TriangleMesh mesh = shape.mesh; | |
| 287 | 2511 | for (Vec3d& p : mesh.positions) { | |
| 288 | 2506 | p = placement.Apply(p); | |
| 289 | } | ||
| 290 | 10 | return mesh; | |
| 291 | 5 | } | |
| 292 | |||
| 293 | } // namespace terraingen | ||
| 294 |