albedo/crater_albedo.cpp
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | #include "terraingen/albedo/crater_albedo.h" | ||
| 2 | |||
| 3 | #include "terraingen/layers/craters/secondaries.h" | ||
| 4 | #include "terraingen/noise/fractal.h" | ||
| 5 | #include "terraingen/noise/hash.h" | ||
| 6 | #include "terraingen/noise/value.h" | ||
| 7 | #include "terraingen/random/stream.h" | ||
| 8 | #include "terraingen/raster/grid.h" | ||
| 9 | #include "terraingen/util/checks.h" | ||
| 10 | #include "terraingen/util/log.h" | ||
| 11 | |||
| 12 | #include <algorithm> | ||
| 13 | #include <cmath> | ||
| 14 | #include <numbers> | ||
| 15 | #include <stdexcept> | ||
| 16 | #include <utility> | ||
| 17 | |||
| 18 | namespace terraingen { | ||
| 19 | |||
| 20 | namespace { | ||
| 21 | |||
| 22 | constexpr std::uint32_t kCraterAlbedoStream = 0x414c4244u; | ||
| 23 | constexpr std::uint32_t kObliqueStream = 0x4f424c51u; | ||
| 24 | constexpr std::uint32_t kZoneStream = 0x5a4f4e45u; | ||
| 25 | constexpr std::size_t kBlockRows = 32; | ||
| 26 | |||
| 27 | } // namespace | ||
| 28 | |||
| 29 | 846495 | double RayRetention::Myr(double diameter) const { | |
| 30 | 846495 | double myr = largeMyr; | |
| 31 | 846495 | if (diameter <= smallDiameter) { | |
| 32 | 365831 | myr = smallMyr; | |
| 33 | 480664 | } else if (diameter < midDiameter) { | |
| 34 | 480521 | const double t = std::log(diameter / smallDiameter) / std::log(midDiameter / smallDiameter); | |
| 35 | 480521 | myr = smallMyr * std::pow(midMyr / smallMyr, t); | |
| 36 | 143 | } else if (diameter <= largeDiameter) { | |
| 37 | 142 | myr = midMyr; | |
| 38 | } | ||
| 39 | 846495 | return scale * myr; | |
| 40 | } | ||
| 41 | |||
| 42 | 424 | void ValidateRayRetention(const RayRetention& retention) { | |
| 43 | 848 | if (!(IsPositive(retention.smallDiameter) && IsPositive(retention.midDiameter) && | |
| 44 | 424 | IsPositive(retention.largeDiameter) && retention.smallDiameter < retention.midDiameter && | |
| 45 | 423 | retention.midDiameter <= retention.largeDiameter)) { | |
| 46 | 1 | throw std::invalid_argument("Ray retention diameters must be positive and increasing"); | |
| 47 | } | ||
| 48 | 845 | if (!(IsPositive(retention.smallMyr) && IsPositive(retention.midMyr) && | |
| 49 | 422 | IsPositive(retention.largeMyr) && IsPositive(retention.scale))) { | |
| 50 | 1 | throw std::invalid_argument("Ray retention times and their scale must be positive"); | |
| 51 | } | ||
| 52 | 422 | } | |
| 53 | |||
| 54 | 33 | double ObliqueImpacts::MaxStretch() const { | |
| 55 | 33 | if (!enabled) { | |
| 56 | 1 | return 1.0; | |
| 57 | } | ||
| 58 | 32 | const double ratio = downrangeSpeedRatio.Highest(); | |
| 59 | 32 | const double a = (ratio - 1.0) / (ratio + 1.0); | |
| 60 | 32 | return (1.0 + a) * (1.0 + a) / (1.0 + 0.5 * a * a); | |
| 61 | } | ||
| 62 | |||
| 63 | 431 | void ValidateObliqueImpacts(const ObliqueImpacts& oblique) { | |
| 64 | 431 | ValidateValueDistribution(oblique.downrangeSpeedRatio); | |
| 65 | 859 | if (oblique.downrangeSpeedRatio.Lowest() < 1.0 || | |
| 66 | 429 | !std::isfinite(oblique.downrangeSpeedRatio.Highest())) { | |
| 67 | 1 | throw std::invalid_argument("The downrange speed ratio must be at least 1 and bounded"); | |
| 68 | } | ||
| 69 | 429 | if (!(oblique.butterflyDegrees >= 0.0 && oblique.butterflyDegrees <= oblique.forbiddenDegrees && | |
| 70 | 428 | oblique.forbiddenDegrees <= oblique.asymmetricDegrees && | |
| 71 | 428 | oblique.asymmetricDegrees <= 90.0)) { | |
| 72 | throw std::invalid_argument( | ||
| 73 | 2 | "Oblique impact angles must satisfy 0 <= butterfly <= forbidden <= asymmetric <= 90"); | |
| 74 | } | ||
| 75 | 427 | if (!(oblique.zoneHalfWidthDegrees >= 0.0 && oblique.zoneHalfWidthDegrees <= 180.0)) { | |
| 76 | 1 | throw std::invalid_argument("A forbidden zone's half width must be in [0, 180] degrees"); | |
| 77 | } | ||
| 78 | 426 | if (oblique.heading && !std::isfinite(*oblique.heading)) { | |
| 79 | 1 | throw std::invalid_argument("An impact heading must be finite"); | |
| 80 | } | ||
| 81 | 425 | } | |
| 82 | |||
| 83 | 470965 | double CraterRays::MedianLength(double diameter) const { | |
| 84 | 470965 | return ElliottMedianRayLength(diameter, lengthCoefficient, lengthExponent); | |
| 85 | } | ||
| 86 | |||
| 87 | 433 | void ValidateCraterRays(const CraterRays& rays) { | |
| 88 | 433 | ValidateValueDistribution(rays.count); | |
| 89 | 433 | ValidateValueDistribution(rays.lengthSpread); | |
| 90 | 433 | ValidateValueDistribution(rays.widthDegrees); | |
| 91 | 433 | if (rays.count.Lowest() < 0.0 || rays.count.Highest() > 64.0) { | |
| 92 | 2 | throw std::invalid_argument("Craters must have 0 to 64 rays"); | |
| 93 | } | ||
| 94 | 431 | if (!(rays.widthDegrees.Lowest() > 0.0 && rays.widthDegrees.Highest() < 180.0)) { | |
| 95 | 1 | throw std::invalid_argument("Ray widths must be in (0, 180) degrees"); | |
| 96 | } | ||
| 97 | 430 | if (!(IsPositive(rays.lengthCoefficient) && IsPositive(rays.lengthExponent))) { | |
| 98 | 1 | throw std::invalid_argument("The ray length law needs a positive coefficient and exponent"); | |
| 99 | } | ||
| 100 | 429 | if (!IsNonNegative(rays.excess)) { | |
| 101 | 2 | throw std::invalid_argument("The ray maturity excess must not be negative"); | |
| 102 | } | ||
| 103 | 427 | if (!(IsInUnitInterval(rays.patchiness) && IsPositive(rays.patchRadii))) { | |
| 104 | throw std::invalid_argument( | ||
| 105 | 2 | "Ray patchiness must be in [0, 1] with a positive patch wavelength"); | |
| 106 | } | ||
| 107 | 425 | ValidateObliqueImpacts(rays.oblique); | |
| 108 | 425 | if (rays.layout) { | |
| 109 | 5 | ValidateRayLayout(*rays.layout); | |
| 110 | } | ||
| 111 | 424 | } | |
| 112 | |||
| 113 | 435 | void ValidateCraterAlbedoTexture(const CraterAlbedoTexture& texture) { | |
| 114 | 435 | if (!(IsInUnitInterval(texture.streaks) && IsInUnitInterval(texture.floorShare))) { | |
| 115 | ✗ | throw std::invalid_argument("Crater albedo streaks and floor share must be in [0, 1]"); | |
| 116 | } | ||
| 117 | 435 | if (!(texture.lobes >= 0.0 && texture.lobes <= 0.9 && texture.lobeCount >= 1 && | |
| 118 | 433 | texture.lobeCount <= 64)) { | |
| 119 | throw std::invalid_argument( | ||
| 120 | 2 | "Crater albedo lobes must be in [0, 0.9] with a lobe count of 1 to 64"); | |
| 121 | } | ||
| 122 | 433 | RequirePositive(texture.streakWidthDegrees, "Crater albedo streak width"); | |
| 123 | 433 | RequirePositive(texture.streakLengthRadii, "Crater albedo streak length"); | |
| 124 | 433 | RequireNonNegative(texture.continuousRadii, "Crater albedo continuous ejecta reach"); | |
| 125 | 432 | } | |
| 126 | |||
| 127 | 1843 | void ValidateReflectanceZone(const ReflectanceZone& zone) { | |
| 128 | 1843 | RequireUnitInterval(zone.share, "The share of craters with a reflectance zone"); | |
| 129 | 1841 | ValidateValueDistribution(zone.change); | |
| 130 | 1839 | ValidateValueDistribution(zone.scale); | |
| 131 | 1837 | if (!(zone.change.Lowest() > -1.0 && std::isfinite(zone.change.Highest()))) { | |
| 132 | 2 | throw std::invalid_argument("A reflectance zone's change must be above -1 and bounded"); | |
| 133 | } | ||
| 134 | 3663 | if (!(IsPositive(zone.coefficient) && IsPositive(zone.exponent) && zone.scale.Lowest() > 0.0 && | |
| 135 | 1828 | std::isfinite(zone.scale.Highest()))) { | |
| 136 | throw std::invalid_argument( | ||
| 137 | 7 | "A reflectance zone needs a positive reach law and a positive, bounded reach scale"); | |
| 138 | } | ||
| 139 | 1828 | if (!(zone.raggedness >= 0.0 && zone.raggedness < 1.0 && zone.feather > 0.0 && | |
| 140 | 1824 | zone.feather <= 1.0)) { | |
| 141 | throw std::invalid_argument( | ||
| 142 | 4 | "A reflectance zone's raggedness must be in [0, 1) and its feather in (0, 1]"); | |
| 143 | } | ||
| 144 | 1824 | RequirePositive(zone.fingerWidthDegrees, "A reflectance zone's finger width"); | |
| 145 | 1822 | if (zone.retentionMyr) { | |
| 146 | 1358 | RequirePositive(*zone.retentionMyr, "A reflectance zone's retention time"); | |
| 147 | } | ||
| 148 | 1821 | } | |
| 149 | |||
| 150 | 3001 | double CraterReflectanceZones::MaxReach(double diameter) const { | |
| 151 | 3001 | double reach = 0.0; | |
| 152 | 15005 | for (const ReflectanceZone* zone : {&proximalHigh, &proximalLow, &distalHigh, &distalLow}) { | |
| 153 | 24008 | reach = std::fmax(reach, zone->coefficient * std::pow(diameter, zone->exponent) * | |
| 154 | 12004 | zone->scale.Highest()); | |
| 155 | } | ||
| 156 | 3001 | return reach; | |
| 157 | } | ||
| 158 | |||
| 159 | 467 | void ValidateCraterReflectanceZones(const CraterReflectanceZones& zones) { | |
| 160 | 2288 | for (const ReflectanceZone* zone : | |
| 161 | 2755 | {&zones.proximalHigh, &zones.proximalLow, &zones.distalHigh, &zones.distalLow}) { | |
| 162 | 1842 | ValidateReflectanceZone(*zone); | |
| 163 | } | ||
| 164 | 446 | RequirePositive(zones.maxDiameter, "The largest crater with reflectance zones"); | |
| 165 | 443 | RequirePositive(zones.excavationRatio, "The excavation depth ratio"); | |
| 166 | 441 | RequireNonNegative(zones.matureDepth, "The mature regolith depth"); | |
| 167 | 439 | const ZoneShadows& s = zones.shadows; | |
| 168 | 875 | if (!(s.jetAngleDegrees > 0.0 && s.jetAngleDegrees < 90.0 && IsPositive(s.depth) && | |
| 169 | 436 | s.azimuths >= 8 && s.azimuths <= 8192 && s.radialSamples >= 8 && | |
| 170 | 435 | s.radialSamples <= 8192)) { | |
| 171 | throw std::invalid_argument("Zone shadows need a jet angle in (0, 90) degrees, a positive " | ||
| 172 | 5 | "depth and 8 to 8192 azimuths and samples"); | |
| 173 | } | ||
| 174 | 434 | } | |
| 175 | |||
| 176 | 452 | void ValidateCraterAlbedo(const CraterAlbedoOptions& options) { | |
| 177 | 452 | ValidateCraterReflectanceZones(options.zones); | |
| 178 | 433 | ValidateCraterAlbedoTexture(options.texture); | |
| 179 | 432 | ValidateCraterRays(options.rays); | |
| 180 | 424 | ValidateRayRetention(options.retention); | |
| 181 | 422 | ValidateResolutionFade(options.fade); | |
| 182 | 421 | ValidateCraterSelection(options.selection, "Crater albedo"); | |
| 183 | 419 | RequireUnitInterval(options.share, "The share of craters with bright material"); | |
| 184 | 833 | if (!(IsPositive(options.backgroundMaturity) && std::isfinite(options.maxMaturity) && | |
| 185 | 416 | options.maxMaturity >= options.backgroundMaturity)) { | |
| 186 | throw std::invalid_argument( | ||
| 187 | 2 | "The background maturity must be positive and the cap at least as high"); | |
| 188 | } | ||
| 189 | 415 | if (!(IsNonNegative(options.interiorExcess) && IsNonNegative(options.ejectaExcess))) { | |
| 190 | 1 | throw std::invalid_argument("Crater maturity excesses must not be negative"); | |
| 191 | } | ||
| 192 | 827 | if (!(IsPositive(options.ejectaExponent) && std::isfinite(options.ejectaReachRadii) && | |
| 193 | 413 | options.ejectaReachRadii > 1.0 && IsPositive(options.fadeExponent) && | |
| 194 | 411 | std::isfinite(options.originReflectance))) { | |
| 195 | throw std::invalid_argument("Crater albedo needs a positive ejecta exponent and fade " | ||
| 196 | 3 | "exponent and a reach over 1 radius"); | |
| 197 | } | ||
| 198 | 411 | } | |
| 199 | |||
| 200 | /// The key of `crater`'s texture. | ||
| 201 | 6671622 | std::uint32_t TextureKey(const Crater& crater, const CraterAlbedoTexture& texture) { | |
| 202 | 6671622 | return PcgHash(static_cast<std::uint32_t>(crater.id) ^ | |
| 203 | 13343244 | PcgHash(static_cast<std::uint32_t>(crater.id >> 32) ^ texture.seed)); | |
| 204 | } | ||
| 205 | |||
| 206 | /// How far `crater`'s ejecta edge reaches at `angle` radians around it, over its mean reach. | ||
| 207 | 9655290 | double LobeScale(double angle, std::uint32_t key, const CraterAlbedoTexture& texture) { | |
| 208 | 9655290 | if (texture.lobes == 0.0) { | |
| 209 | 8381190 | return 1.0; | |
| 210 | } | ||
| 211 | 1274100 | FractalParams params; | |
| 212 | 1274100 | params.seed = key ^ 0x4C4F4245u; | |
| 213 | 1274100 | params.wavelength = 1.0f; | |
| 214 | 1274100 | params.amplitude = 1.0f; | |
| 215 | 1274100 | params.octaves = 2; | |
| 216 | 1274100 | const double ring = static_cast<double>(texture.lobeCount) / (2.0 * std::numbers::pi); | |
| 217 | const double n = | ||
| 218 | 1274100 | static_cast<double>(Fractal3(ring * std::cos(angle), ring * std::sin(angle), 0.0, params)); | |
| 219 | 1274100 | return 1.0 + texture.lobes * std::clamp(n, -1.0, 1.0); | |
| 220 | } | ||
| 221 | |||
| 222 | /// The share of the excess the filaments leave at `radii` radii and `angle` radians around a | ||
| 223 | /// crater: 1 along a filament, down to 1 - `streaks` between them. | ||
| 224 | 1115074 | double StreakShare(double radii, double angle, std::uint32_t key, | |
| 225 | const CraterAlbedoTexture& texture) { | ||
| 226 | 1115074 | if (texture.streaks == 0.0) { | |
| 227 | 632951 | return 1.0; | |
| 228 | } | ||
| 229 | 482123 | FractalParams params; | |
| 230 | 482123 | params.seed = key; | |
| 231 | 482123 | params.wavelength = 1.0f; | |
| 232 | 482123 | params.amplitude = 1.0f; | |
| 233 | 482123 | params.octaves = 3; | |
| 234 | 482123 | const double ring = 180.0 / (std::numbers::pi * texture.streakWidthDegrees); | |
| 235 | 964246 | const double n = static_cast<double>(Fractal3(ring * std::cos(angle), ring * std::sin(angle), | |
| 236 | 482123 | radii / texture.streakLengthRadii, params)); | |
| 237 | 482123 | const double filament = Smoothstep(std::clamp(0.5 + 1.5 * n, 0.0, 1.0)); | |
| 238 | 482123 | return 1.0 - texture.streaks * (1.0 - filament); | |
| 239 | } | ||
| 240 | |||
| 241 | /// How far a jet from a crater's rim runs above the ground, on `azimuths` rays of `samples` | ||
| 242 | /// samples around the crater, of which `rays` rays from `first` and their first `used` samples | ||
| 243 | /// are held: those reaching the region evaluated. | ||
| 244 | struct ShadowProfile { | ||
| 245 | std::size_t azimuths = 0; | ||
| 246 | std::size_t samples = 0; | ||
| 247 | std::size_t first = 0; | ||
| 248 | std::size_t rays = 0; | ||
| 249 | std::size_t used = 0; | ||
| 250 | double start = 0.0; | ||
| 251 | double step = 0.0; | ||
| 252 | double depth = 1.0; | ||
| 253 | std::vector<float> gap; | ||
| 254 | |||
| 255 | /// How deep in shadow the ground is at `distance` metres and `angle` radians from the | ||
| 256 | /// crater's centre, from 0 in the open to 1 in full shadow. | ||
| 257 | 5005872 | double Shade(double distance, double angle) const { | |
| 258 | 5005872 | if (distance <= start) { | |
| 259 | 39 | return 0.0; | |
| 260 | } | ||
| 261 | 5005833 | const auto rings = static_cast<double>(azimuths); | |
| 262 | 5005833 | const double r = std::fmin((distance - start) / step, static_cast<double>(samples - 1)); | |
| 263 | 5005833 | double a = angle / (2.0 * std::numbers::pi) * rings; | |
| 264 | 5005833 | a -= std::floor(a / rings) * rings; | |
| 265 | 5005833 | const std::size_t a0 = std::min(static_cast<std::size_t>(a), azimuths - 1); | |
| 266 | 5005833 | const std::size_t a1 = (a0 + 1) % azimuths; | |
| 267 | 5005833 | const std::size_t r0 = std::min(static_cast<std::size_t>(r), samples - 1); | |
| 268 | 5005833 | const std::size_t r1 = std::min(r0 + 1, samples - 1); | |
| 269 | 5005833 | const double fa = a - static_cast<double>(a0); | |
| 270 | 5005833 | const double fr = r - static_cast<double>(r0); | |
| 271 | 20023332 | const auto at = [&](std::size_t i, std::size_t j) { | |
| 272 | 20023332 | const std::size_t k = (i + azimuths - first) % azimuths; | |
| 273 | 20023332 | return k < rays ? static_cast<double>(gap[k * used + std::min(j, used - 1)]) : 0.0; | |
| 274 | 5005833 | }; | |
| 275 | 5005833 | const double g = (1.0 - fa) * ((1.0 - fr) * at(a0, r0) + fr * at(a0, r1)) + | |
| 276 | 5005833 | fa * ((1.0 - fr) * at(a1, r0) + fr * at(a1, r1)); | |
| 277 | 5005833 | return Smoothstep(std::clamp(g / depth, 0.0, 1.0)); | |
| 278 | } | ||
| 279 | }; | ||
| 280 | |||
| 281 | /// The rays, of `azimuths` around `centre`, that reach `region`, with one more on each side, and | ||
| 282 | /// the farthest distance in metres of the region from `centre`. | ||
| 283 | struct RayWindow { | ||
| 284 | std::size_t first = 0; | ||
| 285 | std::size_t rays = 0; | ||
| 286 | double far = 0.0; | ||
| 287 | }; | ||
| 288 | |||
| 289 | 86 | RayWindow RaysOver(const Vec2d& centre, const Bounds2d& region, const Surface& surface, | |
| 290 | std::size_t azimuths) { | ||
| 291 | 86 | RayWindow window; | |
| 292 | 86 | window.rays = azimuths; | |
| 293 | 86 | double minX = INFINITY; | |
| 294 | 86 | double maxX = -INFINITY; | |
| 295 | 86 | double minY = INFINITY; | |
| 296 | 86 | double maxY = -INFINITY; | |
| 297 | 86 | std::vector<Vec2d> corners; | |
| 298 | 258 | for (const double x : {region.minX, region.maxX}) { | |
| 299 | 516 | for (const double y : {region.minY, region.maxY}) { | |
| 300 | 344 | const Vec2d c = surface.Offset(centre, {x, y}); | |
| 301 | 344 | corners.push_back(c); | |
| 302 | 344 | minX = std::fmin(minX, c.x); | |
| 303 | 344 | maxX = std::fmax(maxX, c.x); | |
| 304 | 344 | minY = std::fmin(minY, c.y); | |
| 305 | 344 | maxY = std::fmax(maxY, c.y); | |
| 306 | 344 | window.far = std::fmax(window.far, std::hypot(c.x, c.y)); | |
| 307 | } | ||
| 308 | } | ||
| 309 | 86 | if (minX <= 0.0 && maxX >= 0.0 && minY <= 0.0 && maxY >= 0.0) { | |
| 310 | 4 | return window; | |
| 311 | } | ||
| 312 | 82 | const double mid = std::atan2(0.5 * (minY + maxY), 0.5 * (minX + maxX)); | |
| 313 | 82 | double half = 0.0; | |
| 314 | 410 | for (const Vec2d& c : corners) { | |
| 315 | 328 | half = std::fmax( | |
| 316 | 328 | half, std::fabs(std::remainder(std::atan2(c.y, c.x) - mid, 2.0 * std::numbers::pi))); | |
| 317 | } | ||
| 318 | 82 | const double ray = 2.0 * std::numbers::pi / static_cast<double>(azimuths); | |
| 319 | 82 | const double from = std::floor((mid - half) / ray) - 1.0; | |
| 320 | 82 | const double to = std::ceil((mid + half) / ray) + 1.0; | |
| 321 | 82 | if (to - from + 1.0 >= static_cast<double>(azimuths)) { | |
| 322 | ✗ | return window; | |
| 323 | } | ||
| 324 | 82 | const auto n = static_cast<std::int64_t>(azimuths); | |
| 325 | 82 | window.first = static_cast<std::size_t>(((static_cast<std::int64_t>(from) % n) + n) % n); | |
| 326 | 82 | window.rays = static_cast<std::size_t>(to - from) + 1; | |
| 327 | 82 | return window; | |
| 328 | 86 | } | |
| 329 | |||
| 330 | /// One reflectance zone of a crater, drawn. | ||
| 331 | struct DrawnZone { | ||
| 332 | double change = 0.0; | ||
| 333 | double strength = 0.0; | ||
| 334 | double resolution = 0.0; | ||
| 335 | double reach = 0.0; | ||
| 336 | double raggedness = 0.0; | ||
| 337 | double ring = 0.0; | ||
| 338 | double feather = 1.0; | ||
| 339 | bool streaked = false; | ||
| 340 | bool distal = false; | ||
| 341 | std::uint32_t key = 0; | ||
| 342 | |||
| 343 | /// How much of the change the zone gives at `distance` metres and `angle` radians from the | ||
| 344 | /// crater's centre, from 1 inside to 0 beyond its ragged edge, whose fingers bend and break | ||
| 345 | /// with distance. | ||
| 346 | 17960299 | double Cover(double distance, double angle) const { | |
| 347 | 17960299 | if (distance >= reach) { | |
| 348 | 13216211 | return 0.0; | |
| 349 | } | ||
| 350 | 4744088 | double finger = 1.0; | |
| 351 | 4744088 | if (raggedness > 0.0) { | |
| 352 | 4730767 | FractalParams params; | |
| 353 | 4730767 | params.seed = key; | |
| 354 | 4730767 | params.wavelength = 1.0f; | |
| 355 | 4730767 | params.amplitude = 1.0f; | |
| 356 | 4730767 | params.octaves = 4; | |
| 357 | 9461534 | const double n = static_cast<double>(Fractal3( | |
| 358 | 4730767 | ring * std::cos(angle), ring * std::sin(angle), 3.0 * distance / reach, params)); | |
| 359 | 4730767 | finger = 1.0 - raggedness * std::clamp(0.5 - n, 0.0, 1.0); | |
| 360 | } | ||
| 361 | 4744088 | const double t = distance / (reach * finger); | |
| 362 | 4744088 | return 1.0 - Smoothstep(std::clamp((t - (1.0 - feather)) / feather, 0.0, 1.0)); | |
| 363 | } | ||
| 364 | }; | ||
| 365 | |||
| 366 | /// One ray, in the frame of its crater. | ||
| 367 | struct CraterAlbedo::Ray { | ||
| 368 | double cos = 1.0; | ||
| 369 | double sin = 0.0; | ||
| 370 | double length = 0.0; | ||
| 371 | double sinHalfWidth = 0.0; | ||
| 372 | std::uint32_t seed = 0; | ||
| 373 | |||
| 374 | /// The ray's excess at `along` metres from the crater's centre along it and `lateral` across, | ||
| 375 | /// for a crater of `radius` metres. | ||
| 376 | 84871305 | double ExcessAt(double along, double lateral, double radius, const CraterRays& rays) const { | |
| 377 | 84871305 | if (along <= radius || along >= length) { | |
| 378 | 59323996 | return 0.0; | |
| 379 | } | ||
| 380 | 25547309 | const double half = std::fmin(along, 2.0 * radius) * sinHalfWidth; | |
| 381 | 25547309 | if (lateral >= half) { | |
| 382 | 24860960 | return 0.0; | |
| 383 | } | ||
| 384 | 686349 | const double q = 1.0 - (lateral / half) * (lateral / half); | |
| 385 | 686349 | double strength = rays.excess * q * q * (length - along) / (length - radius); | |
| 386 | 686349 | if (rays.patchiness > 0.0) { | |
| 387 | 684119 | const double wavelength = rays.patchRadii * radius; | |
| 388 | 684119 | const double u = along / wavelength; | |
| 389 | 684119 | const double iu = std::floor(u); | |
| 390 | 684119 | const float noise = ValueNoise2(static_cast<std::int32_t>(iu), 0, | |
| 391 | 684119 | static_cast<float>(u - iu), 0.5f, seed); | |
| 392 | 684119 | strength *= 1.0 - rays.patchiness * (0.5 + 0.5 * static_cast<double>(noise)); | |
| 393 | } | ||
| 394 | 686349 | return strength; | |
| 395 | } | ||
| 396 | }; | ||
| 397 | |||
| 398 | /// A crater's bright material, ready to evaluate at points. | ||
| 399 | struct CraterAlbedo::Prepared { | ||
| 400 | Crater crater; | ||
| 401 | CraterOutline outline; | ||
| 402 | double weight = 0.0; | ||
| 403 | double reach = 0.0; | ||
| 404 | std::vector<Ray> rays; | ||
| 405 | std::vector<DrawnZone> zones; | ||
| 406 | std::shared_ptr<const ShadowProfile> shadow; | ||
| 407 | |||
| 408 | /// The relative reflectance change of the zones at (`x`, `y`), inner zones covering outer as | ||
| 409 | /// long as they last, the bright proximal zone broken by `texture`'s filaments. | ||
| 410 | 9266057 | double Change(const CraterAlbedoTexture& texture, const Surface& surface, double x, | |
| 411 | double y) const { | ||
| 412 | 9266057 | if (zones.empty()) { | |
| 413 | 3974710 | return 0.0; | |
| 414 | } | ||
| 415 | 5291347 | const Vec2d offset = surface.Offset({crater.x, crater.y}, {x, y}); | |
| 416 | 5291347 | const double distance = std::hypot(offset.x, offset.y); | |
| 417 | 5291347 | const double angle = std::atan2(offset.y, offset.x); | |
| 418 | 5291347 | const double shade = shadow ? shadow->Shade(distance, angle) : 0.0; | |
| 419 | 5291347 | double change = 0.0; | |
| 420 | 5291347 | double open = 1.0; | |
| 421 | 23251646 | for (const DrawnZone& zone : zones) { | |
| 422 | 17960299 | const double cover = zone.Cover(distance, angle) * (zone.distal ? 1.0 - shade : 1.0); | |
| 423 | 17960299 | if (cover == 0.0) { | |
| 424 | 15210213 | continue; | |
| 425 | } | ||
| 426 | 2750086 | const double streak = zone.streaked | |
| 427 | 2827345 | ? StreakShare(2.0 * distance / crater.diameter, angle, | |
| 428 | 77259 | TextureKey(crater, texture), texture) | |
| 429 | 2750086 | : 1.0; | |
| 430 | 2750086 | change += open * cover * zone.strength * zone.resolution * streak * zone.change; | |
| 431 | 2750086 | open *= 1.0 - cover * zone.strength; | |
| 432 | } | ||
| 433 | 5291347 | return change; | |
| 434 | } | ||
| 435 | |||
| 436 | /// The excess optical maturity at (`x`, `y`). | ||
| 437 | 9733051 | double Excess(const CraterAlbedoOptions& o, const Surface& surface, double x, double y) const { | |
| 438 | 9733051 | if (weight == 0.0) { | |
| 439 | 1 | return 0.0; | |
| 440 | } | ||
| 441 | 9733050 | const double radius = 0.5 * crater.diameter; | |
| 442 | 9733050 | const double radii = CraterProfileDistance(crater, outline, x, y, surface) / radius; | |
| 443 | 9733050 | const CraterAlbedoTexture& texture = o.texture; | |
| 444 | 9733050 | const bool textured = texture.streaks > 0.0 || texture.lobes > 0.0; | |
| 445 | 9733050 | const Vec2d offset = surface.Offset({crater.x, crater.y}, {x, y}); | |
| 446 | 9733050 | const double angle = textured ? std::atan2(offset.y, offset.x) : 0.0; | |
| 447 | 9733050 | const std::uint32_t key = textured ? TextureKey(crater, texture) : 0u; | |
| 448 | 9733050 | if (radii < 1.0) { | |
| 449 | 77760 | double interior = o.interiorExcess; | |
| 450 | 77760 | const double wall = Smoothstep(std::clamp((radii - 0.3) / 0.6, 0.0, 1.0)); | |
| 451 | 77760 | if (texture.floorShare < 1.0) { | |
| 452 | 38215 | interior *= texture.floorShare + (1.0 - texture.floorShare) * wall; | |
| 453 | } | ||
| 454 | 77760 | const double streak = StreakShare(radii, angle, key, texture); | |
| 455 | 77760 | return weight * interior * (1.0 - (1.0 - streak) * wall); | |
| 456 | } | ||
| 457 | 9655290 | double excess = 0.0; | |
| 458 | 9655290 | const double lobed = radii / LobeScale(angle, key, texture); | |
| 459 | 9655290 | if (lobed < o.ejectaReachRadii) { | |
| 460 | 589742 | const double edge = std::pow(o.ejectaReachRadii, -o.ejectaExponent); | |
| 461 | 589742 | excess = o.ejectaExcess * (std::pow(std::fmax(lobed, 1.0), -o.ejectaExponent) - edge) / | |
| 462 | 589742 | (1.0 - edge); | |
| 463 | } | ||
| 464 | 94526595 | for (const Ray& ray : rays) { | |
| 465 | 84871305 | const double along = offset.x * ray.cos + offset.y * ray.sin; | |
| 466 | 84871305 | const double lateral = std::fabs(offset.y * ray.cos - offset.x * ray.sin); | |
| 467 | 84871305 | excess += ray.ExcessAt(along, lateral, radius, o.rays); | |
| 468 | } | ||
| 469 | 9655290 | if (excess == 0.0) { | |
| 470 | 8695235 | return 0.0; | |
| 471 | } | ||
| 472 | 960055 | const double streak = StreakShare(radii, angle, key, texture); | |
| 473 | 960055 | if (texture.continuousRadii > 0.0) { | |
| 474 | const double ramp = | ||
| 475 | 364020 | Smoothstep(std::clamp((radii - 1.0) / texture.continuousRadii, 0.0, 1.0)); | |
| 476 | 364020 | return weight * excess * (1.0 - (1.0 - streak) * ramp); | |
| 477 | } | ||
| 478 | 596035 | return weight * excess * streak; | |
| 479 | } | ||
| 480 | }; | ||
| 481 | |||
| 482 | 399 | CraterAlbedo::CraterAlbedo(std::shared_ptr<const CraterLayer> craters, CraterAlbedoOptions options) | |
| 483 | 399 | : craters_(std::move(craters)), options_(options) { | |
| 484 | 399 | if (!craters_) { | |
| 485 | 1 | throw std::invalid_argument("Crater albedo needs a crater layer"); | |
| 486 | } | ||
| 487 | 398 | ValidateCraterAlbedo(options_); | |
| 488 | 419 | } | |
| 489 | |||
| 490 | 858018 | bool CraterAlbedo::Eligible(const Crater& crater, RandomStream& random) const { | |
| 491 | 858018 | return options_.selection.Includes(crater) && random.NextDouble() < options_.share; | |
| 492 | } | ||
| 493 | |||
| 494 | 29929238 | CraterAlbedo::Ray CraterAlbedo::DrawRay(RandomStream& random, double median, | |
| 495 | const CraterImpact& impact, | ||
| 496 | const RayDirection* direction) const { | ||
| 497 | 29929238 | const CraterRays& rays = options_.rays; | |
| 498 | 29929238 | const ObliqueImpacts& oblique = rays.oblique; | |
| 499 | 29929238 | const double asymmetry = (impact.speedRatio - 1.0) / (impact.speedRatio + 1.0); | |
| 500 | 29929238 | const double zone = oblique.zoneHalfWidthDegrees * std::numbers::pi / 180.0; | |
| 501 | 29929238 | const double drawn = 2.0 * std::numbers::pi * random.NextDouble(); | |
| 502 | 29929238 | const double azimuth = direction != nullptr ? direction->azimuth : drawn; | |
| 503 | 29929238 | const double l1 = random.NextDouble(); | |
| 504 | 29929238 | const double l2 = random.NextDouble(); | |
| 505 | 29929238 | const double w1 = random.NextDouble(); | |
| 506 | 29929238 | const double w2 = random.NextDouble(); | |
| 507 | 29929238 | Ray ray; | |
| 508 | 29929238 | ray.cos = std::cos(azimuth); | |
| 509 | 29929238 | ray.sin = std::sin(azimuth); | |
| 510 | 29929238 | ray.length = median * std::exp(rays.lengthSpread.Sample(l1, l2)); | |
| 511 | 29929238 | const double fromDownrange = std::remainder(azimuth - impact.heading, 2.0 * std::numbers::pi); | |
| 512 | 29929238 | if (asymmetry > 0.0) { | |
| 513 | 15396597 | const double speed = 1.0 + asymmetry * std::cos(fromDownrange); | |
| 514 | 15396597 | ray.length *= speed * speed / (1.0 + 0.5 * asymmetry * asymmetry); | |
| 515 | } | ||
| 516 | 29929238 | const bool uprange = std::fabs(fromDownrange) > std::numbers::pi - zone; | |
| 517 | 29929238 | const bool downrange = std::fabs(fromDownrange) < zone; | |
| 518 | 29929238 | if ((oblique.enabled && impact.angleDegrees < oblique.forbiddenDegrees && uprange) || | |
| 519 | 28195354 | (oblique.enabled && impact.angleDegrees < oblique.butterflyDegrees && downrange)) { | |
| 520 | 1945808 | ray.length = 0.0; | |
| 521 | } | ||
| 522 | 29929238 | ray.sinHalfWidth = | |
| 523 | direction != nullptr | ||
| 524 | 29929238 | ? direction->sinHalfWidth | |
| 525 | 29923401 | : std::sin(0.5 * rays.widthDegrees.Sample(w1, w2) * std::numbers::pi / 180.0); | |
| 526 | 29929238 | ray.seed = random.NextU32(); | |
| 527 | 29929238 | return ray; | |
| 528 | } | ||
| 529 | |||
| 530 | 470931 | void CraterAlbedo::AddRays(Prepared& prepared, RandomStream& random) const { | |
| 531 | 470931 | const Crater& crater = prepared.crater; | |
| 532 | 470931 | const double c1 = random.NextDouble(); | |
| 533 | 470931 | const double c2 = random.NextDouble(); | |
| 534 | 470931 | const std::vector<RayDirection> directions = options_.rays.layout | |
| 535 | 470931 | ? RayDirections(crater, *options_.rays.layout) | |
| 536 | 470931 | : std::vector<RayDirection>{}; | |
| 537 | 470931 | const auto count = options_.rays.layout | |
| 538 | 470931 | ? static_cast<int>(directions.size()) | |
| 539 | 470202 | : static_cast<int>(std::lround(options_.rays.count.Sample(c1, c2))); | |
| 540 | 470931 | const double median = options_.rays.MedianLength(crater.diameter); | |
| 541 | 470931 | const CraterImpact impact = Impact(crater); | |
| 542 | 30400169 | for (int i = 0; i < count; ++i) { | |
| 543 | const Ray ray = | ||
| 544 | 29935075 | DrawRay(random, median, impact, | |
| 545 | 29935075 | options_.rays.layout ? &directions[static_cast<std::size_t>(i)] : nullptr); | |
| 546 | 29929238 | if (ray.length > 0.5 * crater.diameter) { | |
| 547 | 27982977 | prepared.reach = std::fmax(prepared.reach, ray.length); | |
| 548 | 27982977 | prepared.rays.push_back(ray); | |
| 549 | } | ||
| 550 | } | ||
| 551 | 470931 | } | |
| 552 | |||
| 553 | 12115 | void CraterAlbedo::AddZones(Prepared& prepared, double left, double resolution) const { | |
| 554 | 12115 | const CraterReflectanceZones& z = options_.zones; | |
| 555 | 12115 | const double diameter = prepared.crater.diameter; | |
| 556 | 12115 | const IdWords id = SplitId(prepared.crater.id); | |
| 557 | 12115 | RandomStream random(z.seed, kZoneStream, id.low, id.high); | |
| 558 | 12115 | const bool immature = z.excavationRatio * diameter > z.matureDepth; | |
| 559 | 12115 | std::vector<DrawnZone> distal; | |
| 560 | 60575 | for (const ReflectanceZone* zone : | |
| 561 | 72690 | {&z.proximalHigh, &z.proximalLow, &z.distalHigh, &z.distalLow}) { | |
| 562 | 48460 | const double present = random.NextDouble(); | |
| 563 | 48460 | const double c1 = random.NextDouble(); | |
| 564 | 48460 | const double c2 = random.NextDouble(); | |
| 565 | 48460 | const double s1 = random.NextDouble(); | |
| 566 | 48460 | const double s2 = random.NextDouble(); | |
| 567 | 48460 | DrawnZone drawn; | |
| 568 | 48460 | drawn.key = random.NextU32(); | |
| 569 | const double life = | ||
| 570 | 48460 | zone->retentionMyr | |
| 571 | 48460 | ? std::fmin(left, 1.0 - 1000.0 * prepared.crater.ageGyr / *zone->retentionMyr) | |
| 572 | 48460 | : left; | |
| 573 | 48460 | if (present >= zone->share || (zone == &z.proximalHigh && !immature) || life <= 0.0) { | |
| 574 | 24152 | continue; | |
| 575 | } | ||
| 576 | 24308 | drawn.change = zone->change.Sample(c1, c2); | |
| 577 | 24308 | drawn.strength = std::pow(life, options_.fadeExponent); | |
| 578 | 24308 | drawn.resolution = resolution; | |
| 579 | 24308 | drawn.reach = | |
| 580 | 24308 | zone->coefficient * std::pow(diameter, zone->exponent) * zone->scale.Sample(s1, s2); | |
| 581 | 24308 | drawn.raggedness = zone->raggedness; | |
| 582 | 24308 | drawn.ring = 180.0 / (std::numbers::pi * zone->fingerWidthDegrees); | |
| 583 | 24308 | drawn.feather = zone->feather; | |
| 584 | 24308 | drawn.streaked = zone == &z.proximalHigh; | |
| 585 | 24308 | drawn.distal = zone == &z.distalHigh || zone == &z.distalLow; | |
| 586 | 24308 | prepared.reach = std::fmax(prepared.reach, drawn.reach); | |
| 587 | 24308 | (zone == &z.distalHigh || zone == &z.distalLow ? distal : prepared.zones).push_back(drawn); | |
| 588 | } | ||
| 589 | 12115 | std::sort(distal.begin(), distal.end(), | |
| 590 | 5920 | [](const DrawnZone& a, const DrawnZone& b) { return a.reach < b.reach; }); | |
| 591 | 12115 | prepared.zones.insert(prepared.zones.end(), distal.begin(), distal.end()); | |
| 592 | 12115 | } | |
| 593 | |||
| 594 | 858018 | CraterAlbedo::Prepared CraterAlbedo::Prepare(const Crater& crater, double spacing) const { | |
| 595 | 858018 | Prepared prepared; | |
| 596 | 858018 | prepared.crater = crater; | |
| 597 | 858018 | const CraterAlbedoOptions& o = options_; | |
| 598 | 858018 | const IdWords id = SplitId(crater.id); | |
| 599 | 858018 | RandomStream random(o.seed, kCraterAlbedoStream, id.low, id.high); | |
| 600 | 858018 | if (!Eligible(crater, random)) { | |
| 601 | 11531 | return prepared; | |
| 602 | } | ||
| 603 | 846487 | const double left = 1.0 - 1000.0 * crater.ageGyr / o.retention.Myr(crater.diameter); | |
| 604 | 846487 | const double resolution = ResolutionWeight(crater.diameter, spacing, o.fade); | |
| 605 | 846487 | if (left <= 0.0 || resolution == 0.0) { | |
| 606 | 361766 | return prepared; | |
| 607 | } | ||
| 608 | 484721 | prepared.weight = std::pow(left, o.fadeExponent) * resolution; | |
| 609 | 484721 | prepared.outline = OutlineOf(crater, craters_->Options().apply.irregularity); | |
| 610 | 484721 | prepared.reach = o.ejectaReachRadii * (0.5 * crater.diameter) * prepared.outline.MaxScale() * | |
| 611 | 484721 | (1.0 + o.texture.lobes); | |
| 612 | 484721 | if (o.rays.enabled) { | |
| 613 | 470931 | AddRays(prepared, random); | |
| 614 | } | ||
| 615 | 484721 | if (o.zones.enabled && crater.diameter <= o.zones.maxDiameter) { | |
| 616 | 12115 | AddZones(prepared, left, resolution); | |
| 617 | } | ||
| 618 | 484721 | return prepared; | |
| 619 | } | ||
| 620 | |||
| 621 | 471314 | double CraterAlbedo::Excess(const Crater& crater, double x, double y, double spacing) const { | |
| 622 | 471314 | return Prepare(crater, spacing).Excess(options_, craters_->Options().surface, x, y); | |
| 623 | } | ||
| 624 | |||
| 625 | 4320 | double CraterAlbedo::ReflectanceChange(const Crater& crater, double x, double y, double spacing, | |
| 626 | const Ground* ground) const { | ||
| 627 | 4320 | Prepared prepared = Prepare(crater, spacing); | |
| 628 | 4320 | if (ground != nullptr && options_.zones.shadows.enabled) { | |
| 629 | 11 | AddShadows(prepared, *ground, PointBounds(x, y)); | |
| 630 | } | ||
| 631 | 8640 | return prepared.Change(options_.texture, craters_->Options().surface, x, y); | |
| 632 | 4320 | } | |
| 633 | |||
| 634 | 314 | void CraterAlbedo::AddShadows(Prepared& prepared, const Ground& ground, | |
| 635 | const Bounds2d& region) const { | ||
| 636 | 314 | const ZoneShadows& s = options_.zones.shadows; | |
| 637 | 314 | double reach = 0.0; | |
| 638 | 881 | for (const DrawnZone& zone : prepared.zones) { | |
| 639 | 567 | reach = zone.distal ? std::fmax(reach, zone.reach) : reach; | |
| 640 | } | ||
| 641 | 314 | const Crater& crater = prepared.crater; | |
| 642 | 314 | const double radius = 0.5 * crater.diameter; | |
| 643 | 314 | if (reach <= radius) { | |
| 644 | 228 | return; | |
| 645 | } | ||
| 646 | 86 | auto profile = std::make_shared<ShadowProfile>(); | |
| 647 | 86 | const auto azimuths = static_cast<std::size_t>(s.azimuths); | |
| 648 | 86 | const auto samples = static_cast<std::size_t>(s.radialSamples); | |
| 649 | 86 | profile->azimuths = azimuths; | |
| 650 | 86 | profile->samples = samples; | |
| 651 | 86 | profile->start = radius; | |
| 652 | 86 | profile->step = (reach - radius) / static_cast<double>(samples - 1); | |
| 653 | 86 | profile->depth = s.depth; | |
| 654 | 86 | const Surface& surface = craters_->Options().surface; | |
| 655 | 86 | const RayWindow window = RaysOver({crater.x, crater.y}, region, surface, azimuths); | |
| 656 | 86 | const double used = std::ceil((window.far - radius) / profile->step) + 2.0; | |
| 657 | 86 | profile->first = window.first; | |
| 658 | 86 | profile->rays = window.rays; | |
| 659 | 86 | profile->used = static_cast<std::size_t>(std::clamp(used, 2.0, static_cast<double>(samples))); | |
| 660 | 86 | const std::size_t count = profile->rays * profile->used; | |
| 661 | 86 | std::vector<Vec2d> points; | |
| 662 | 86 | points.reserve(count); | |
| 663 | 3934 | for (std::size_t k = 0; k < profile->rays; ++k) { | |
| 664 | 3848 | const std::size_t i = (window.first + k) % azimuths; | |
| 665 | 3848 | const double angle = | |
| 666 | 3848 | 2.0 * std::numbers::pi * static_cast<double>(i) / static_cast<double>(azimuths); | |
| 667 | 768120 | for (std::size_t j = 0; j < profile->used; ++j) { | |
| 668 | 764272 | const double distance = radius + static_cast<double>(j) * profile->step; | |
| 669 | 764272 | points.push_back(surface.Locate( | |
| 670 | 764272 | {crater.x, crater.y}, {distance * std::cos(angle), distance * std::sin(angle)})); | |
| 671 | } | ||
| 672 | } | ||
| 673 | 86 | std::vector<float> heights(count); | |
| 674 | 86 | ground.HeightsAt(points, heights, profile->step); | |
| 675 | 86 | const double drop = std::tan(s.jetAngleDegrees * std::numbers::pi / 180.0) * profile->step; | |
| 676 | 86 | profile->gap.assign(count, 0.0f); | |
| 677 | 3934 | for (std::size_t k = 0; k < profile->rays; ++k) { | |
| 678 | 3848 | const std::size_t row = k * profile->used; | |
| 679 | 3848 | double jet = static_cast<double>(heights[row]); | |
| 680 | 764272 | for (std::size_t j = 1; j < profile->used; ++j) { | |
| 681 | 760424 | const double h = static_cast<double>(heights[row + j]); | |
| 682 | 760424 | if (std::isnan(h) || std::isnan(jet)) { | |
| 683 | 9 | jet = h; | |
| 684 | 9 | continue; | |
| 685 | } | ||
| 686 | 760415 | jet = std::fmax(h, jet - drop); | |
| 687 | 760415 | profile->gap[row + j] = static_cast<float>(jet - h); | |
| 688 | } | ||
| 689 | } | ||
| 690 | 86 | prepared.shadow = std::move(profile); | |
| 691 | 86 | } | |
| 692 | |||
| 693 | 491233 | CraterImpact CraterAlbedo::Impact(const Crater& crater) const { | |
| 694 | 491233 | const ObliqueImpacts& oblique = options_.rays.oblique; | |
| 695 | 491233 | CraterImpact impact; | |
| 696 | 491233 | if (!oblique.enabled) { | |
| 697 | 2891 | return impact; | |
| 698 | } | ||
| 699 | 488342 | const IdWords id = SplitId(crater.id); | |
| 700 | 488342 | RandomStream random(options_.seed, kObliqueStream, id.low, id.high); | |
| 701 | 488342 | impact.angleDegrees = std::asin(std::sqrt(random.NextDouble())) * 180.0 / std::numbers::pi; | |
| 702 | 488342 | const double h = random.NextDouble(); | |
| 703 | 488342 | impact.heading = oblique.heading ? *oblique.heading : 2.0 * std::numbers::pi * h; | |
| 704 | 488342 | const double r1 = random.NextDouble(); | |
| 705 | 488342 | const double r2 = random.NextDouble(); | |
| 706 | 488342 | if (impact.angleDegrees < oblique.asymmetricDegrees) { | |
| 707 | 251148 | impact.speedRatio = oblique.downrangeSpeedRatio.Sample(r1, r2); | |
| 708 | } | ||
| 709 | 488342 | return impact; | |
| 710 | } | ||
| 711 | |||
| 712 | 10 | double CraterAlbedo::Reach(const Crater& crater) const { | |
| 713 | 10 | return Prepare(crater, 0.0).reach; | |
| 714 | } | ||
| 715 | |||
| 716 | 380 | std::vector<CraterAlbedo::Prepared> CraterAlbedo::Gather(const Bounds2d& region, double spacing, | |
| 717 | const Ground* ground) const { | ||
| 718 | 380 | const double largest = craters_->MaxFeatureDiameter(); | |
| 719 | 380 | const double rayReach = options_.rays.enabled | |
| 720 | 380 | ? options_.rays.MedianLength(largest) * | |
| 721 | 31 | std::exp(options_.rays.lengthSpread.Highest()) * | |
| 722 | 31 | options_.rays.oblique.MaxStretch() | |
| 723 | 380 | : 0.0; | |
| 724 | const double ejectaReach = | ||
| 725 | 760 | std::fmax(rayReach, 0.5 * largest * options_.ejectaReachRadii * | |
| 726 | 380 | MaxOutlineScale(craters_->Options().apply.irregularity) * | |
| 727 | 380 | (1.0 + options_.texture.lobes)); | |
| 728 | 380 | const Surface& surface = craters_->Options().surface; | |
| 729 | 380 | std::vector<Prepared> prepared; | |
| 730 | 1713 | const auto take = [&](const CraterList& list) { | |
| 731 | 384087 | for (std::size_t i = 0; i < list.Size(); ++i) { | |
| 732 | 382374 | Prepared one = Prepare(list.At(i), spacing); | |
| 733 | 382374 | if (one.weight == 0.0) { | |
| 734 | 373294 | continue; | |
| 735 | } | ||
| 736 | 9080 | const Crater& c = one.crater; | |
| 737 | 9080 | if (!surface.Contains(surface.Expand(region, one.reach), {c.x, c.y})) { | |
| 738 | 6544 | continue; | |
| 739 | } | ||
| 740 | 2536 | prepared.push_back(std::move(one)); | |
| 741 | 382374 | } | |
| 742 | 1713 | }; | |
| 743 | 380 | take(craters_->Craters(surface.Expand(region, ejectaReach), spacing)); | |
| 744 | 380 | if (options_.zones.enabled) { | |
| 745 | 333 | for (double high = std::fmin(largest, options_.zones.maxDiameter); | |
| 746 | 1666 | options_.zones.MaxReach(high) > ejectaReach; high *= 0.5) { | |
| 747 | 1333 | take(craters_->Craters(surface.Expand(region, options_.zones.MaxReach(high)), spacing, | |
| 748 | 0.5 * high, high)); | ||
| 749 | } | ||
| 750 | } | ||
| 751 | 380 | std::sort(prepared.begin(), prepared.end(), | |
| 752 | 11652 | [](const Prepared& a, const Prepared& b) { return a.crater.id < b.crater.id; }); | |
| 753 | 380 | prepared.erase(std::unique(prepared.begin(), prepared.end(), | |
| 754 | 2166 | [](const Prepared& a, const Prepared& b) { | |
| 755 | 2166 | return a.crater.id == b.crater.id; | |
| 756 | }), | ||
| 757 | 380 | prepared.end()); | |
| 758 | 380 | if (ground != nullptr && options_.zones.enabled && options_.zones.shadows.enabled) { | |
| 759 | 69 | std::size_t shadowed = 0; | |
| 760 | 372 | for (Prepared& p : prepared) { | |
| 761 | 303 | AddShadows(p, *ground, region); | |
| 762 | 303 | shadowed += p.shadow ? 1u : 0u; | |
| 763 | } | ||
| 764 | 69 | LogDebug("albedo", "Sampled the ground for the zone shadows of ", shadowed, " craters"); | |
| 765 | } | ||
| 766 | 760 | return prepared; | |
| 767 | } | ||
| 768 | |||
| 769 | namespace { | ||
| 770 | |||
| 771 | /// Moves `value` to the reflectance at the background maturity plus `excess`, capped, then | ||
| 772 | /// changes it by the relative reflectance `change`, never below 0. | ||
| 773 | 5964560 | void Brighten(AlbedoSample& value, std::size_t bands, double excess, double change, | |
| 774 | const CraterAlbedoOptions& o) { | ||
| 775 | 5964560 | if (excess == 0.0 && change == 0.0) { | |
| 776 | 3982581 | return; | |
| 777 | } | ||
| 778 | 1981979 | const double maturity = std::fmin(o.backgroundMaturity + excess, o.maxMaturity); | |
| 779 | 1981979 | const double factor = std::fmax(0.0, 1.0 + change); | |
| 780 | 3987969 | for (std::size_t band = 0; band < bands; ++band) { | |
| 781 | 2005990 | if (std::isnan(value[band])) { | |
| 782 | 1 | continue; | |
| 783 | } | ||
| 784 | 2005989 | double reflectance = value[band]; | |
| 785 | 2005989 | if (excess != 0.0) { | |
| 786 | 1023960 | reflectance = ReflectanceAtMaturity(reflectance, o.backgroundMaturity, maturity, | |
| 787 | 1023960 | o.originReflectance); | |
| 788 | } | ||
| 789 | 2005989 | value[band] = static_cast<float>(factor * reflectance); | |
| 790 | } | ||
| 791 | } | ||
| 792 | |||
| 793 | } // namespace | ||
| 794 | |||
| 795 | 351 | void CraterAlbedo::ApplyAt(AlbedoSample& value, std::size_t bands, double x, double y, | |
| 796 | const AlbedoContext& context) const { | ||
| 797 | 351 | ValidateAlbedoContext(context); | |
| 798 | 350 | const Surface& surface = craters_->Options().surface; | |
| 799 | 350 | double excess = 0.0; | |
| 800 | 350 | double change = 0.0; | |
| 801 | 761 | for (const Prepared& p : Gather(PointBounds(x, y), context.spacing, context.ground)) { | |
| 802 | 411 | excess += p.Excess(options_, surface, x, y); | |
| 803 | 411 | change += p.Change(options_.texture, surface, x, y); | |
| 804 | 350 | } | |
| 805 | 350 | Brighten(value, bands, excess, change, options_); | |
| 806 | 350 | } | |
| 807 | |||
| 808 | 31 | void CraterAlbedo::Apply(std::vector<Raster<float>>& bands, const AlbedoContext& context) const { | |
| 809 | 31 | ValidateAlbedoContext(context); | |
| 810 | 31 | if (bands.empty()) { | |
| 811 | 1 | return; | |
| 812 | } | ||
| 813 | 30 | const Surface& surface = craters_->Options().surface; | |
| 814 | 30 | const Raster<float>& first = bands.front(); | |
| 815 | const std::vector<Prepared> prepared = | ||
| 816 | 30 | Gather(RasterGrid::FromRaster(first).Bounds(), context.spacing, context.ground); | |
| 817 | 30 | const std::size_t count = bands.size(); | |
| 818 | 30 | ForEachAlbedoBlock(bands, context.threads, kBlockRows, [&](std::size_t row0, std::size_t row1) { | |
| 819 | 268 | const double ya = first.CellCenter(0, row0).y; | |
| 820 | 268 | const double yb = first.CellCenter(0, row1 - 1).y; | |
| 821 | 268 | const double yMin = std::fmin(ya, yb); | |
| 822 | 268 | const double yMax = std::fmax(ya, yb); | |
| 823 | 268 | std::vector<const Prepared*> near; | |
| 824 | 13284 | for (const Prepared& p : prepared) { | |
| 825 | 13016 | const double reachY = p.reach / surface.UnitLength(); | |
| 826 | 13016 | if (p.crater.y + reachY >= yMin && p.crater.y - reachY <= yMax) { | |
| 827 | 1151 | near.push_back(&p); | |
| 828 | } | ||
| 829 | } | ||
| 830 | 536 | return AlbedoCellVisitor([this, near = std::move(near), &surface, | |
| 831 | 268 | count](AlbedoSample& value, double x, double y) { | |
| 832 | 5964210 | double excess = 0.0; | |
| 833 | 5964210 | double change = 0.0; | |
| 834 | 37682631 | for (const Prepared* p : near) { | |
| 835 | 31718421 | if (!surface.IsBody() && std::fabs(x - p->crater.x) > p->reach) { | |
| 836 | 22457095 | continue; | |
| 837 | } | ||
| 838 | 9261326 | excess += p->Excess(options_, surface, x, y); | |
| 839 | 9261326 | change += p->Change(options_.texture, surface, x, y); | |
| 840 | } | ||
| 841 | 5964210 | Brighten(value, count, excess, change, options_); | |
| 842 | 5964746 | }); | |
| 843 | 268 | }); | |
| 844 | 30 | } | |
| 845 | |||
| 846 | } // namespace terraingen | ||
| 847 |