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 | 46841 | cells.erase(std::unique(cells.begin(), cells.end()), cells.end()); | |
| 36 | 46841 | } | |
| 37 | |||
| 38 | } // namespace | ||
| 39 | |||
| 40 | 11 | CubeSphereCells::CubeSphereCells(double radius) | |
| 41 | 11 | : radius_(radius), surface_(Surface::Sphere(radius)) {} | |
| 42 | |||
| 43 | 52 | CubeSphereCells::CubeSphereCells(const Surface& surface) : radius_(0.0), surface_(surface) { | |
| 44 | 52 | if (!surface.IsBody()) { | |
| 45 | 1 | throw std::invalid_argument("Cube-sphere cells need a body, not a plane"); | |
| 46 | } | ||
| 47 | 51 | radius_ = surface.MinRadius(); | |
| 48 | 51 | if (surface.IsSphere()) { | |
| 49 | 42 | return; | |
| 50 | } | ||
| 51 | 9 | if (surface.Shape() != nullptr) { | |
| 52 | 2 | weightCap_ = surface.Shape()->MaxAreaPerSteradian(); | |
| 53 | 2 | return; | |
| 54 | } | ||
| 55 | 98 | for (std::size_t face = 0; face < detail::kCubeFaces.size(); ++face) { | |
| 56 | 2730 | for (std::int32_t i = 0; i < kCapCells; ++i) { | |
| 57 | 174720 | for (std::int32_t j = 0; j < kCapCells; ++j) { | |
| 58 | 344064 | weightCap_ = std::max(weightCap_, | |
| 59 | 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 | 1299612 | RequirePositive(cellSize, "Cell size"); | |
| 78 | 1299611 | const double n = std::ceil(radius_ * 2.0 * kQuarter / cellSize); | |
| 79 | 1299611 | if (!(n < static_cast<double>(std::numeric_limits<std::int32_t>::max() / 6))) { | |
| 80 | 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 | 511241 | const std::int32_t n = CellsPerEdge(cellSize); | |
| 87 | 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 | 1022482 | return std::clamp<std::int32_t>(k, 0, n - 1); | |
| 95 | 511241 | }; | |
| 96 | 511241 | return {static_cast<std::int32_t>(face) * n + index(detail::kCubeFaces[face].u), | |
| 97 | 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 | 179404 | for (std::int32_t dj = -rings; dj <= rings; ++dj) { | |
| 109 | 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 | 529193 | if (c >= 0 && c < n && r >= 0 && r < n) { | |
| 113 | 430320 | cells.push_back({face * n + c, r}); | |
| 114 | 430320 | continue; | |
| 115 | } | ||
| 116 | 395492 | for (const double su : {0.1, 0.5, 0.9}) { | |
| 117 | 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 | 889857 | if (!(std::fabs(alpha) < 2.0 * kQuarter - step && | |
| 121 | 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 | 121628 | const Vec2d point = Surface::Sphere(radius_).PointOf(p); | |
| 126 | 121628 | cells.push_back(CellOf(point.x, point.y, cellSize)); | |
| 127 | } | ||
| 128 | } | ||
| 129 | } | ||
| 130 | } | ||
| 131 | 39393 | SortUnique(cells); | |
| 132 | 39393 | return cells; | |
| 133 | } | ||
| 134 | |||
| 135 | 3096 | std::vector<CellIndex> CubeSphereCells::CellsCovering(const Bounds2d& region, | |
| 136 | double cellSize) const { | ||
| 137 | 3096 | RequirePositive(cellSize, "Cell size"); | |
| 138 | 6191 | if (!(std::isfinite(region.minX) && std::isfinite(region.minY) && std::isfinite(region.maxX) && | |
| 139 | 3095 | std::isfinite(region.maxY) && region.maxX >= region.minX && region.maxY >= region.minY)) { | |
| 140 | 2 | throw std::invalid_argument("Region must be finite and not inverted"); | |
| 141 | } | ||
| 142 | 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 | 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 | 221762 | for (std::int64_t c = 0; c <= columns; ++c) { | |
| 154 | 212454 | const double lon = columns == 0 ? box.minX | |
| 155 | 212453 | : box.minX + width * static_cast<double>(c) / | |
| 156 | 212453 | static_cast<double>(columns); | |
| 157 | 212454 | seeds.push_back(CellOf(lon, lat, cellSize)); | |
| 158 | } | ||
| 159 | } | ||
| 160 | 3094 | SortUnique(seeds); | |
| 161 | 3094 | std::vector<CellIndex> cells; | |
| 162 | 31569 | for (const CellIndex& seed : seeds) { | |
| 163 | 28475 | const std::vector<CellIndex> around = Neighbours(seed, cellSize, 1); | |
| 164 | 28475 | cells.insert(cells.end(), around.begin(), around.end()); | |
| 165 | 28475 | } | |
| 166 | 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 | 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 | 46682 | for (const CellIndex& cell : cells) { | |
| 174 | 43588 | if (std::binary_search(seeds.begin(), seeds.end(), cell)) { | |
| 175 | 28475 | continue; | |
| 176 | } | ||
| 177 | 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 | 45339 | for (const double a : {a1, a2}) { | |
| 183 | 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 | 15113 | widest = std::max(widest, std::acos(std::clamp(nearest, -1.0, 1.0))); | |
| 189 | 15113 | ring.push_back(cell); | |
| 190 | 15113 | centres.push_back(sphere.PointOf(centre)); | |
| 191 | } | ||
| 192 | 3094 | if (ring.empty()) { | |
| 193 | 1834 | return cells; | |
| 194 | } | ||
| 195 | 1260 | const Bounds2d reach = sphere.Expand(box, radius_ * (widest * (1.0 + 1e-9) + 1e-12)); | |
| 196 | 1260 | std::vector<CellIndex> kept = seeds; | |
| 197 | 16373 | for (std::size_t k = 0; k < ring.size(); ++k) { | |
| 198 | 15113 | if (sphere.Contains(reach, centres[k])) { | |
| 199 | 2769 | kept.push_back(ring[k]); | |
| 200 | } | ||
| 201 | } | ||
| 202 | 1260 | SortUnique(kept); | |
| 203 | 1260 | return kept; | |
| 204 | 3094 | } | |
| 205 | |||
| 206 | 309165 | double CubeSphereCells::CellArea(const CellIndex& cell, double cellSize) const { | |
| 207 | 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 | 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 | 5593 | detail::kCubeFaces[static_cast<std::size_t>(cell.i / CellsPerEdge(cellSize))]; | |
| 219 | 5593 | double area = 0.0; | |
| 220 | 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 | 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 | 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 | 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 | 172897 | const auto [aLo, aHi] = absolute(a1, a2); | |
| 243 | 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 | 691588 | for (std::size_t k = 0; k < 3; ++k) { | |
| 257 | 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 | 12401 | if (weightCap_ == 0.0) { | |
| 267 | 11439 | return {CellArea(cell, cellSize), 0.0}; | |
| 268 | } | ||
| 269 | 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 | 962 | const auto face = static_cast<std::size_t>(cell.i / CellsPerEdge(cellSize)); | |
| 272 | 962 | const double bound = surface_.Shape() != nullptr | |
| 273 | 1827 | ? weightCap_ | |
| 274 | 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 | 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 | 26914875 | for (int iteration = 0; iteration < 64; ++iteration) { | |
| 294 | 26500800 | const double mid = 0.5 * (lo + hi); | |
| 295 | 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 | 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 | 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 |