physics/cloud.cpp
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // SPDX-License-Identifier: BSD-3-Clause | ||
| 2 | /// @file cloud.cpp | ||
| 3 | /// @brief Cloud droplet phase tables, Monte Carlo slab transport, and tileable noise textures. | ||
| 4 | // | ||
| 5 | // The noise textures' recipes follow three-geospatial's clouds (MIT, Copyright (c) 2024 Shota | ||
| 6 | // Matsuda) and Hillaire's TileableVolumeNoise (MIT, Copyright (c) 2017 Sebastien Hillaire); see | ||
| 7 | // cloud.h. | ||
| 8 | #include "ernie/physics/cloud.h" | ||
| 9 | |||
| 10 | #include "ernie/physics/scattering.h" | ||
| 11 | #include "ernie/util/disk_cache.h" | ||
| 12 | #include "ernie/util/math.h" | ||
| 13 | #include "ernie/util/noise.h" | ||
| 14 | #include "ernie/util/parallel.h" | ||
| 15 | #include "ernie/util/random.h" | ||
| 16 | |||
| 17 | #include <glm/glm.hpp> | ||
| 18 | |||
| 19 | #include <algorithm> | ||
| 20 | #include <cmath> | ||
| 21 | #include <complex> | ||
| 22 | #include <cstddef> | ||
| 23 | #include <filesystem> | ||
| 24 | #include <future> | ||
| 25 | #include <numbers> | ||
| 26 | #include <optional> | ||
| 27 | #include <span> | ||
| 28 | #include <string> | ||
| 29 | #include <utility> | ||
| 30 | |||
| 31 | namespace ernie { | ||
| 32 | |||
| 33 | namespace { | ||
| 34 | |||
| 35 | using std::numbers::pi; | ||
| 36 | |||
| 37 | /// The integral of f(theta) sin(theta) over [theta0, theta1] by the trapezoid rule; times 2 pi it | ||
| 38 | /// is the integral over that band of the sphere. | ||
| 39 | 6435 | double SphereTrapezoid(double f0, double theta0, double f1, double theta1) { | |
| 40 | 6435 | return 0.5 * (f0 * std::sin(theta0) + f1 * std::sin(theta1)) * (theta1 - theta0); | |
| 41 | } | ||
| 42 | |||
| 43 | // ---------------------------------------------------------------------------------------------- | ||
| 44 | // Droplet phase functions | ||
| 45 | |||
| 46 | /// Liquid water's refractive index across the visible, its absorption negligible there. | ||
| 47 | constexpr std::complex<double> kWater(1.333, 0.0); | ||
| 48 | |||
| 49 | /// One band's phase function with its forward peak's share of the light, which it is divided by | ||
| 50 | /// 1 minus. | ||
| 51 | 48 | struct BandPhase { | |
| 52 | std::vector<double> phase; | ||
| 53 | double peak = 0.0; | ||
| 54 | }; | ||
| 55 | |||
| 56 | /// The share of light within `cut` of forward, by the trapezoid rule over the sphere. | ||
| 57 | 12 | double ForwardPeak(const std::vector<double>& phase, const std::vector<double>& thetas, | |
| 58 | double cut) { | ||
| 59 | 12 | double peak = 0.0; | |
| 60 | 1332 | for (std::size_t i = 0; i + 1 < thetas.size() && thetas[i] < cut; ++i) { | |
| 61 | 1320 | const double to = std::min(thetas[i + 1], cut); | |
| 62 | 1320 | const double along = (to - thetas[i]) / (thetas[i + 1] - thetas[i]); | |
| 63 | 1320 | const double end = phase[i] + (phase[i + 1] - phase[i]) * along; | |
| 64 | 1320 | peak += SphereTrapezoid(phase[i], thetas[i], end, to) * 2.0 * pi; | |
| 65 | } | ||
| 66 | 12 | return std::clamp(peak, 0.0, 0.95); | |
| 67 | } | ||
| 68 | |||
| 69 | /// One band's Mie phase function against the extinction less its forward peak: the peak kept for | ||
| 70 | /// single scattering (its silver lining about the sun), folded away from the rest. | ||
| 71 | 12 | BandPhase DropletBandPhase(const LogNormalMode& mode, double wavelengthNm, | |
| 72 | const std::vector<double>& thetas, const std::vector<double>& cosines, | ||
| 73 | double cut) { | ||
| 74 | 12 | ParticleOptics optics = MieParticles(std::span(&mode, 1), kWater, wavelengthNm, cosines, 48); | |
| 75 | 12 | BandPhase band{std::move(optics.phase), 0.0}; | |
| 76 | 12 | band.peak = ForwardPeak(band.phase, thetas, cut); | |
| 77 | 12300 | for (double& value : band.phase) { | |
| 78 | 12288 | value /= 1.0 - band.peak; | |
| 79 | } | ||
| 80 | 12 | return band; | |
| 81 | 12 | } | |
| 82 | |||
| 83 | /// The bands' mean beyond the peak, held at the peak's edge value within it so light scattered | ||
| 84 | /// many times shows no hole about the sun, and normalised to integrate to 1 again. | ||
| 85 | 3 | std::vector<double> FoldedMean(std::vector<double> mean, const std::vector<double>& thetas, | |
| 86 | double cut) { | ||
| 87 | 3 | const std::size_t n = thetas.size(); | |
| 88 | 3 | std::size_t edge = 0; | |
| 89 | 333 | while (edge + 1 < n && thetas[edge] < cut) { | |
| 90 | ++edge; | ||
| 91 | } | ||
| 92 | 3 | std::fill(mean.begin(), mean.begin() + static_cast<std::ptrdiff_t>(edge), mean[edge]); | |
| 93 | 3 | double total = 0.0; | |
| 94 | 3072 | for (std::size_t i = 0; i + 1 < n; ++i) { | |
| 95 | 3069 | total += SphereTrapezoid(mean[i], thetas[i], mean[i + 1], thetas[i + 1]) * 2.0 * pi; | |
| 96 | } | ||
| 97 | 3075 | for (double& value : mean) { | |
| 98 | 3072 | value /= total; | |
| 99 | } | ||
| 100 | 3 | return mean; | |
| 101 | } | ||
| 102 | |||
| 103 | // ---------------------------------------------------------------------------------------------- | ||
| 104 | // Transport in slabs | ||
| 105 | |||
| 106 | /// Scattering angles drawn from a phase function's table, by inverting its cumulative distribution | ||
| 107 | /// over the sphere. | ||
| 108 | 2 | class PhaseSampler { | |
| 109 | public: | ||
| 110 | 2 | explicit PhaseSampler(const DropletPhaseTable& phase) { | |
| 111 | 2 | const std::size_t n = phase.angles; | |
| 112 | 2 | const float* mean = phase.values.data() + phase.values.size() - n; | |
| 113 | 2 | thetas_.resize(n); | |
| 114 | 2 | cdf_.assign(n, 0.0); | |
| 115 | 2050 | for (std::size_t i = 0; i < n; ++i) { | |
| 116 | 2048 | thetas_[i] = phase.Angle(i); | |
| 117 | 2048 | if (i > 0) { | |
| 118 | 2046 | cdf_[i] = cdf_[i - 1] + SphereTrapezoid(double(mean[i - 1]), thetas_[i - 1], | |
| 119 | 2046 | double(mean[i]), thetas_[i]); | |
| 120 | } | ||
| 121 | } | ||
| 122 | 2050 | for (double& value : cdf_) { | |
| 123 | 2048 | value /= cdf_.back(); | |
| 124 | } | ||
| 125 | 2 | } | |
| 126 | |||
| 127 | /// The scattering angle at cumulative probability `xi`. | ||
| 128 | 67850525 | double Theta(double xi) const { | |
| 129 | 67850525 | return Interpolate<double>(cdf_, thetas_, xi); | |
| 130 | } | ||
| 131 | |||
| 132 | private: | ||
| 133 | std::vector<double> thetas_; | ||
| 134 | std::vector<double> cdf_; | ||
| 135 | }; | ||
| 136 | |||
| 137 | /// `direction` turned by scattering angle `theta` about azimuth `phi`. | ||
| 138 | 67850525 | glm::dvec3 Scatter(const glm::dvec3& direction, double theta, double phi) { | |
| 139 | 67850525 | const double cosT = std::cos(theta); | |
| 140 | 67850525 | const double sinT = std::sin(theta); | |
| 141 | 67850525 | const double cosP = std::cos(phi); | |
| 142 | 67850525 | const double sinP = std::sin(phi); | |
| 143 | 67850525 | const double s = std::sqrt(std::max(0.0, 1.0 - direction.z * direction.z)); | |
| 144 | 67850525 | if (s < 1e-6) { | |
| 145 | 1844190 | return {sinT * cosP, sinT * sinP, direction.z < 0.0 ? -cosT : cosT}; | |
| 146 | } | ||
| 147 | 66006335 | return {sinT * (direction.x * direction.z * cosP - direction.y * sinP) / s + direction.x * cosT, | |
| 148 | 66006335 | sinT * (direction.y * direction.z * cosP + direction.x * sinP) / s + direction.y * cosT, | |
| 149 | 66006335 | -sinT * cosP * s + direction.z * cosT}; | |
| 150 | } | ||
| 151 | |||
| 152 | /// One slab's tallies: per tally set (light scattered once, then the order sets), per depth, view | ||
| 153 | /// cosine and azimuth, the crossing weights and their entry offsets' sums. | ||
| 154 | 72 | class SlabTallies { | |
| 155 | public: | ||
| 156 | 72 | SlabTallies(const CloudTransportShape& shape, double thickness) | |
| 157 | 72 | : shape_(shape), perDepth_(shape.viewCosines * shape.azimuths), | |
| 158 | 72 | perSet_(shape.depths * perDepth_), values_((kCloudOrderSets + 1) * perSet_ * 4, 0.0) { | |
| 159 | 720 | for (std::size_t k = 0; k < shape.depths; ++k) { | |
| 160 | 648 | planes_.push_back(shape.DepthFraction(k) * thickness); | |
| 161 | } | ||
| 162 | 72 | } | |
| 163 | |||
| 164 | /// Adds a photon's crossings of the depth planes between `from` and `to`, travelling along | ||
| 165 | /// `direction` after scattering `order` times; its entry point was the origin. | ||
| 166 | 99307805 | void AddSegment(std::size_t order, const glm::dvec3& from, const glm::dvec3& to, | |
| 167 | const glm::dvec3& direction) { | ||
| 168 | 99307805 | if (order == 0 || direction.z == 0.0) { | |
| 169 | return; | ||
| 170 | } | ||
| 171 | 67850525 | const std::size_t set = order == 1 ? 0 : 1 + CloudOrderSet(order); | |
| 172 | 67850525 | const double lo = std::min(from.z, to.z); | |
| 173 | 67850525 | const double hi = std::max(from.z, to.z); | |
| 174 | 67850525 | const double weight = 1.0 / std::abs(direction.z); | |
| 175 | 67850525 | const std::size_t cosine = | |
| 176 | 135701050 | std::min(static_cast<std::size_t>( | |
| 177 | 67850525 | std::max(0.0, (direction.z + 1.0) * 0.5 * double(shape_.viewCosines))), | |
| 178 | 67850525 | shape_.viewCosines - 1); | |
| 179 | 67850525 | const std::size_t azimuth = | |
| 180 | 135701050 | std::min(static_cast<std::size_t>(std::atan2(std::abs(direction.y), direction.x) / pi * | |
| 181 | 67850525 | double(shape_.azimuths)), | |
| 182 | 67850525 | shape_.azimuths - 1); | |
| 183 | 67850525 | const double side = direction.y < 0.0 ? -1.0 : 1.0; | |
| 184 | 678505250 | for (std::size_t k = 0; k < shape_.depths; ++k) { | |
| 185 | 610654725 | if (!(planes_[k] > lo && planes_[k] <= hi)) { | |
| 186 | 497985394 | continue; | |
| 187 | } | ||
| 188 | 112669331 | const double along = (planes_[k] - from.z) / direction.z; | |
| 189 | 112669331 | const double u = -(from.x + direction.x * along); | |
| 190 | 112669331 | const double v = -(from.y + direction.y * along) * side; | |
| 191 | 112669331 | double* tally = | |
| 192 | 112669331 | values_.data() + | |
| 193 | 112669331 | (set * perSet_ + k * perDepth_ + cosine * shape_.azimuths + azimuth) * 4; | |
| 194 | 112669331 | tally[0] += weight; | |
| 195 | 112669331 | tally[1] += weight * u; | |
| 196 | 112669331 | tally[2] += weight * v; | |
| 197 | 112669331 | tally[3] += weight * (u * u + v * v); | |
| 198 | } | ||
| 199 | } | ||
| 200 | |||
| 201 | /// The four sums of tally set `set`, cell `i` within the slab. | ||
| 202 | 746496 | const double* At(std::size_t set, std::size_t i) const { | |
| 203 | 746496 | return values_.data() + (set * perSet_ + i) * 4; | |
| 204 | } | ||
| 205 | |||
| 206 | private: | ||
| 207 | const CloudTransportShape& shape_; | ||
| 208 | std::size_t perDepth_; | ||
| 209 | std::size_t perSet_; | ||
| 210 | std::vector<double> planes_; | ||
| 211 | std::vector<double> values_; | ||
| 212 | }; | ||
| 213 | |||
| 214 | /// Traces a slab's photons, each entering the lit face at the origin and scattering without | ||
| 215 | /// absorption until it leaves. | ||
| 216 | 72 | SlabTallies SimulateSlab(const CloudTransportShape& shape, const PhaseSampler& sampler, | |
| 217 | std::size_t light, std::size_t thickness, std::size_t photons, | ||
| 218 | std::uint64_t seed) { | ||
| 219 | 72 | const double muL = shape.LightCosine(light); | |
| 220 | 72 | const double tau = shape.Thickness(thickness); | |
| 221 | 72 | SlabTallies tallies(shape, tau); | |
| 222 | 72 | SplitMix64 random(seed); | |
| 223 | 31457352 | for (std::size_t p = 0; p < photons; ++p) { | |
| 224 | 31457280 | glm::dvec3 position(0.0); | |
| 225 | 31457280 | glm::dvec3 direction(std::sqrt(1.0 - muL * muL), 0.0, muL); | |
| 226 | 31457280 | for (std::size_t order = 0;; ++order) { | |
| 227 | 99307805 | const double step = -std::log(1.0 - random.NextUniform()); | |
| 228 | 99307805 | const glm::dvec3 next = position + direction * step; | |
| 229 | 99307805 | tallies.AddSegment(order, position, next, direction); | |
| 230 | 99307805 | if (next.z < 0.0 || next.z > tau) { | |
| 231 | break; | ||
| 232 | } | ||
| 233 | 67850525 | position = next; | |
| 234 | 135701050 | direction = Scatter(direction, sampler.Theta(random.NextUniform()), | |
| 235 | 67850525 | 2.0 * pi * random.NextUniform()); | |
| 236 | 67850525 | } | |
| 237 | } | ||
| 238 | 72 | return tallies; | |
| 239 | } | ||
| 240 | |||
| 241 | /// Writes a slab's tallies into the table as radiance and collectors. Each photon carries the | ||
| 242 | /// light's flux through the lit face (its cosine) over the photons, over each bin's solid angle | ||
| 243 | /// (both sides of the light's plane). | ||
| 244 | 72 | void StoreSlab(const SlabTallies& tallies, std::size_t light, std::size_t thickness, | |
| 245 | std::size_t photons, CloudTransportTable& table) { | ||
| 246 | 72 | const CloudTransportShape& shape = table.shape; | |
| 247 | 72 | const std::size_t perSlab = shape.depths * shape.viewCosines * shape.azimuths; | |
| 248 | 72 | const double solidAngle = 2.0 / double(shape.viewCosines) * 2.0 * pi / double(shape.azimuths); | |
| 249 | 72 | const double scale = shape.LightCosine(light) / (static_cast<double>(photons) * solidAngle); | |
| 250 | 72 | const std::size_t first = shape.Cell(light, thickness, 0, 0, 0); | |
| 251 | 720 | for (std::size_t set = 0; set <= kCloudOrderSets; ++set) { | |
| 252 | 747144 | for (std::size_t i = 0; i < perSlab; ++i) { | |
| 253 | 746496 | const double* tally = tallies.At(set, i); | |
| 254 | 746496 | const double weight = tally[0]; | |
| 255 | 746496 | const auto radiance = static_cast<float>(weight * scale); | |
| 256 | 746496 | if (set == 0) { | |
| 257 | 82944 | table.single[first + i] = radiance; | |
| 258 | 82944 | continue; | |
| 259 | } | ||
| 260 | 663552 | float* cell = table.values.data() + ((set - 1) * shape.Cells() + first + i) * 4; | |
| 261 | 663552 | cell[0] = radiance; | |
| 262 | 663552 | if (weight > 0.0) { | |
| 263 | 406887 | const double u = tally[1] / weight; | |
| 264 | 406887 | const double v = tally[2] / weight; | |
| 265 | 406887 | cell[1] = static_cast<float>(u); | |
| 266 | 406887 | cell[2] = static_cast<float>(v); | |
| 267 | 813774 | cell[3] = | |
| 268 | 793357 | static_cast<float>(std::sqrt(std::max(0.0, tally[3] / weight - u * u - v * v))); | |
| 269 | } | ||
| 270 | } | ||
| 271 | } | ||
| 272 | 72 | } | |
| 273 | |||
| 274 | // ---------------------------------------------------------------------------------------------- | ||
| 275 | // Noise textures | ||
| 276 | |||
| 277 | /// The weather map's Worley sum: four octaves, each doubling in frequency and losing 5%. | ||
| 278 | 32768 | double WeatherWorley(const glm::dvec3& p, double frequency) { | |
| 279 | 32768 | double noise = 0.0; | |
| 280 | 32768 | double amplitude = 0.4; | |
| 281 | 163840 | for (int i = 0; i < 4; ++i) { | |
| 282 | 131072 | noise += amplitude * (1.0 - TileableWorley(p, frequency)); | |
| 283 | 131072 | frequency *= 2.0; | |
| 284 | 131072 | amplitude *= 0.95; | |
| 285 | } | ||
| 286 | 32768 | return noise; | |
| 287 | } | ||
| 288 | |||
| 289 | /// The centre of texel (i, j, k) of an n-texel texture, in [0, 1]. | ||
| 290 | 54272 | glm::dvec3 TexelCentre(std::size_t i, std::size_t j, std::size_t k, std::size_t n) { | |
| 291 | 54272 | return (glm::dvec3(i, j, k) + 0.5) / double(n); | |
| 292 | } | ||
| 293 | |||
| 294 | /// The shape: Perlin-Worley noise, remapped by a Worley fractal sum. | ||
| 295 | 1 | std::vector<float> MakeShapeNoise(std::size_t n) { | |
| 296 | 1 | std::vector<float> shape(n * n * n, 0.0f); | |
| 297 | 1 | ParallelLoop(n, [&](std::size_t k) { | |
| 298 | 1056 | for (std::size_t j = 0; j < n; ++j) { | |
| 299 | 33792 | for (std::size_t i = 0; i < n; ++i) { | |
| 300 | 32768 | const glm::dvec3 p = TexelCentre(i, j, k, n); | |
| 301 | 32768 | const double perlin = std::clamp(TileablePerlin(p, glm::dvec3(8.0), 3), 0.0, 1.0); | |
| 302 | 32768 | const double worley = 0.625 * (1.0 - TileableWorley(p, 8.0)) + | |
| 303 | 32768 | 0.25 * (1.0 - TileableWorley(p, 32.0)) + | |
| 304 | 32768 | 0.125 * (1.0 - TileableWorley(p, 56.0)); | |
| 305 | 32768 | const double perlinWorley = worley + perlin * (1.0 - worley); | |
| 306 | 32768 | const double fbm = WorleyFbm(p, 4.0); | |
| 307 | 32768 | shape[(k * n + j) * n + i] = | |
| 308 | 32768 | static_cast<float>((perlinWorley - (fbm - 1.0)) / (2.0 - fbm)); | |
| 309 | } | ||
| 310 | } | ||
| 311 | 32 | }); | |
| 312 | 1 | return shape; | |
| 313 | ✗ | } | |
| 314 | |||
| 315 | /// The detail: a Worley fractal sum. | ||
| 316 | 1 | std::vector<float> MakeDetailNoise(std::size_t n) { | |
| 317 | 1 | std::vector<float> detail(n * n * n, 0.0f); | |
| 318 | 1 | ParallelLoop(n, [&](std::size_t k) { | |
| 319 | 272 | for (std::size_t j = 0; j < n; ++j) { | |
| 320 | 4352 | for (std::size_t i = 0; i < n; ++i) { | |
| 321 | 8192 | detail[(k * n + j) * n + i] = | |
| 322 | 4096 | static_cast<float>(WorleyFbm(TexelCentre(i, j, k, n), 1.0)); | |
| 323 | } | ||
| 324 | } | ||
| 325 | 16 | }); | |
| 326 | 1 | return detail; | |
| 327 | ✗ | } | |
| 328 | |||
| 329 | /// The local weather: the low and mid layers' Worley cells, the high layer's streaky Perlin. | ||
| 330 | 1 | std::vector<float> MakeWeatherMap(std::size_t n) { | |
| 331 | 1 | std::vector<float> weather(n * n * 4, 0.0f); | |
| 332 | 1 | ParallelLoop(n, [&](std::size_t j) { | |
| 333 | 16512 | for (std::size_t i = 0; i < n; ++i) { | |
| 334 | 16384 | const glm::dvec3 p = TexelCentre(i, j, 0, n) * glm::dvec3(1.0, 1.0, 0.0); | |
| 335 | 16384 | const double mid = SmoothStep(1.0, 1.4, WeatherWorley(p + glm::dvec3(0.5), 8.0)); | |
| 336 | 16384 | const double low = | |
| 337 | 16384 | std::clamp(SmoothStep(0.8, 1.4, WeatherWorley(p, 16.0)) - mid, 0.0, 1.0); | |
| 338 | 16384 | const double high = | |
| 339 | 16384 | SmoothStep(-0.5, 0.5, TileablePerlin(p, glm::dvec3(6.0, 12.0, 1.0), 8)); | |
| 340 | 16384 | float* texel = weather.data() + (j * n + i) * 4; | |
| 341 | 16384 | texel[0] = static_cast<float>(low); | |
| 342 | 16384 | texel[1] = static_cast<float>(mid); | |
| 343 | 16384 | texel[2] = static_cast<float>(high); | |
| 344 | 16384 | texel[3] = 1.0f; | |
| 345 | } | ||
| 346 | 128 | }); | |
| 347 | 1 | return weather; | |
| 348 | ✗ | } | |
| 349 | |||
| 350 | /// The turbulence: the direction of Perlin noise's curl, mapped to [0, 1]. | ||
| 351 | 1 | std::vector<float> MakeTurbulenceMap(std::size_t n) { | |
| 352 | 1 | std::vector<float> turbulence(n * n * 3, 0.0f); | |
| 353 | 1 | ParallelLoop(n, [&](std::size_t j) { | |
| 354 | 32 | constexpr double kDelta = 0.1; | |
| 355 | 1056 | for (std::size_t i = 0; i < n; ++i) { | |
| 356 | 1024 | const glm::dvec3 p = TexelCentre(i, j, 0, n) * glm::dvec3(1.0, 1.0, 0.0); | |
| 357 | 1024 | const glm::dvec3 x0 = TileablePerlin3(p - glm::dvec3(kDelta, 0.0, 0.0)); | |
| 358 | 1024 | const glm::dvec3 x1 = TileablePerlin3(p + glm::dvec3(kDelta, 0.0, 0.0)); | |
| 359 | 1024 | const glm::dvec3 y0 = TileablePerlin3(p - glm::dvec3(0.0, kDelta, 0.0)); | |
| 360 | 1024 | const glm::dvec3 y1 = TileablePerlin3(p + glm::dvec3(0.0, kDelta, 0.0)); | |
| 361 | 1024 | const glm::dvec3 z0 = TileablePerlin3(p - glm::dvec3(0.0, 0.0, kDelta)); | |
| 362 | 1024 | const glm::dvec3 z1 = TileablePerlin3(p + glm::dvec3(0.0, 0.0, kDelta)); | |
| 363 | 1024 | glm::dvec3 curl(y1.z - y0.z - z1.y + z0.y, z1.x - z0.x - x1.z + x0.z, | |
| 364 | 1024 | x1.y - x0.y - y1.x + y0.x); | |
| 365 | 1024 | const double length = glm::length(curl); | |
| 366 | 1024 | curl = length > 0.0 ? curl / length : glm::dvec3(0.0); | |
| 367 | 1024 | float* texel = turbulence.data() + (j * n + i) * 3; | |
| 368 | 4096 | for (int c = 0; c < 3; ++c) { | |
| 369 | 3072 | texel[c] = static_cast<float>(0.5 * curl[c] + 0.5); | |
| 370 | } | ||
| 371 | } | ||
| 372 | 32 | }); | |
| 373 | 1 | return turbulence; | |
| 374 | ✗ | } | |
| 375 | |||
| 376 | /// The path of a cloud cache file for a key, or nothing when caching is off. | ||
| 377 | ✗ | std::optional<std::filesystem::path> CloudCachePath(const CacheKey& key) { | |
| 378 | ✗ | const std::optional<std::filesystem::path> directory = | |
| 379 | ✗ | CacheDirectory("ERNIE_CLOUD_CACHE", "clouds"); | |
| 380 | ✗ | if (!directory) { | |
| 381 | ✗ | return std::nullopt; | |
| 382 | } | ||
| 383 | ✗ | return *directory / (key.Hex() + ".bin"); | |
| 384 | ✗ | } | |
| 385 | |||
| 386 | } // namespace | ||
| 387 | |||
| 388 | // ---------------------------------------------------------------------------------------------- | ||
| 389 | // Droplet phase functions | ||
| 390 | |||
| 391 | 3 | DropletPhaseTable DropletPhase(double effectiveRadiusUm, double spread, const SpectralBands& bands, | |
| 392 | std::size_t angles) { | ||
| 393 | 3 | DropletPhaseTable table; | |
| 394 | 3 | table.angles = std::max<std::size_t>(angles, 2); | |
| 395 | 3 | const std::size_t n = table.angles; | |
| 396 | 3 | const auto count = static_cast<std::size_t>(std::max(bands.count, 0)); | |
| 397 | 3 | std::vector<double> thetas(n); | |
| 398 | 3 | std::vector<double> cosines(n); | |
| 399 | 3075 | for (std::size_t i = 0; i < n; ++i) { | |
| 400 | 3072 | thetas[i] = table.Angle(i); | |
| 401 | 3072 | cosines[i] = std::cos(thetas[i]); | |
| 402 | } | ||
| 403 | 3 | const LogNormalMode mode = LogNormalMode::FromEffective(std::max(effectiveRadiusUm, 0.5) * 1e-6, | |
| 404 | 3 | 1.0 / (std::max(spread, 0.1) + 2.0)); | |
| 405 | 3 | const double cut = kDropletForwardDegrees * pi / 180.0; | |
| 406 | 3 | std::vector<std::future<BandPhase>> jobs; | |
| 407 | 15 | for (int b = 0; b < bands.count; ++b) { | |
| 408 | 24 | jobs.push_back(std::async(std::launch::async, [&, b] { | |
| 409 | 12 | return DropletBandPhase(mode, bands.Centre(b), thetas, cosines, cut); | |
| 410 | })); | ||
| 411 | } | ||
| 412 | // Each band's row, then the bands' mean beyond the peak, summed in doubles. | ||
| 413 | 3 | table.values.assign((count + 1) * n, 0.0f); | |
| 414 | 3 | std::vector<double> mean(n, 0.0); | |
| 415 | 15 | for (std::size_t b = 0; b < count; ++b) { | |
| 416 | 12 | const BandPhase band = jobs[b].get(); | |
| 417 | 12 | table.forwardFraction += band.peak / static_cast<double>(count); | |
| 418 | 12300 | for (std::size_t i = 0; i < n; ++i) { | |
| 419 | 12288 | table.values[b * n + i] = static_cast<float>(band.phase[i]); | |
| 420 | 12288 | if (thetas[i] >= cut) { | |
| 421 | 10968 | mean[i] += band.phase[i] / static_cast<double>(count); | |
| 422 | } | ||
| 423 | } | ||
| 424 | 12 | } | |
| 425 | 3 | const std::vector<double> folded = FoldedMean(std::move(mean), thetas, cut); | |
| 426 | 3075 | for (std::size_t i = 0; i < n; ++i) { | |
| 427 | 3072 | table.values[count * n + i] = static_cast<float>(folded[i]); | |
| 428 | } | ||
| 429 | 3 | return table; | |
| 430 | 3 | } | |
| 431 | |||
| 432 | // ---------------------------------------------------------------------------------------------- | ||
| 433 | // Transport in slabs | ||
| 434 | |||
| 435 | 46871731 | std::size_t CloudOrderSet(std::size_t order) { | |
| 436 | 46871731 | if (order <= 2) { | |
| 437 | return 0; | ||
| 438 | } | ||
| 439 | 34691167 | if (order <= 8) { | |
| 440 | 21938716 | return (order - 1) / 2; | |
| 441 | } | ||
| 442 | 12752451 | if (order <= 12) { | |
| 443 | return 4; | ||
| 444 | } | ||
| 445 | 8504746 | if (order <= 18) { | |
| 446 | return 5; | ||
| 447 | } | ||
| 448 | 5065312 | return order <= 30 ? 6 : 7; | |
| 449 | } | ||
| 450 | |||
| 451 | 193 | double CloudTransportShape::LightCosine(std::size_t i) const { | |
| 452 | 193 | return static_cast<double>(i + 1) / double(lightCosines); | |
| 453 | } | ||
| 454 | |||
| 455 | 121 | double CloudTransportShape::Thickness(std::size_t i) const { | |
| 456 | 121 | return 0.25 * std::exp2(static_cast<double>(i)); | |
| 457 | } | ||
| 458 | |||
| 459 | 648 | double CloudTransportShape::DepthFraction(std::size_t i) const { | |
| 460 | 648 | return CosineSpacing(i, depths); | |
| 461 | } | ||
| 462 | |||
| 463 | 6912 | double CloudTransportShape::ViewCosine(std::size_t i) const { | |
| 464 | 6912 | return -1.0 + (static_cast<double>(i) + 0.5) * 2.0 / double(viewCosines); | |
| 465 | } | ||
| 466 | |||
| 467 | ✗ | double CloudTransportShape::Azimuth(std::size_t i) const { | |
| 468 | ✗ | return (static_cast<double>(i) + 0.5) * pi / double(azimuths); | |
| 469 | } | ||
| 470 | |||
| 471 | 712733 | std::size_t CloudTransportShape::Cells() const { | |
| 472 | 712733 | return lightCosines * thicknesses * depths * viewCosines * azimuths; | |
| 473 | } | ||
| 474 | |||
| 475 | 6273 | std::size_t CloudTransportShape::Cell(std::size_t light, std::size_t thickness, std::size_t depth, | |
| 476 | std::size_t viewCosine, std::size_t azimuth) const { | ||
| 477 | 6273 | return (((light * thicknesses + thickness) * depths + depth) * viewCosines + viewCosine) * | |
| 478 | 6273 | azimuths + | |
| 479 | 6273 | azimuth; | |
| 480 | } | ||
| 481 | |||
| 482 | 2 | CloudTransportTable SimulateCloudTransport(const DropletPhaseTable& phase, | |
| 483 | const CloudTransportOptions& options) { | ||
| 484 | 2 | const CloudTransportShape& shape = options.shape; | |
| 485 | 2 | const PhaseSampler sampler(phase); | |
| 486 | 2 | CloudTransportTable table; | |
| 487 | 2 | table.shape = shape; | |
| 488 | 2 | table.values.assign(kCloudOrderSets * shape.Cells() * 4, 0.0f); | |
| 489 | 2 | table.single.assign(shape.Cells(), 0.0f); | |
| 490 | // Each light and thickness's slab on its own; the thickest first, as they take longest. | ||
| 491 | 2 | ParallelLoop(shape.lightCosines * shape.thicknesses, [&](std::size_t job) { | |
| 492 | 72 | const std::size_t light = job % shape.lightCosines; | |
| 493 | 72 | const std::size_t thickness = shape.thicknesses - 1 - job / shape.lightCosines; | |
| 494 | 72 | const std::uint64_t seed = HashCombine(options.seed, job); | |
| 495 | 72 | const SlabTallies tallies = | |
| 496 | 72 | SimulateSlab(shape, sampler, light, thickness, options.photons, seed); | |
| 497 | 72 | StoreSlab(tallies, light, thickness, options.photons, table); | |
| 498 | 72 | }); | |
| 499 | 2 | return table; | |
| 500 | 2 | } | |
| 501 | |||
| 502 | ✗ | CloudTransportTable CachedCloudTransport(const CloudTransportOptions& options) { | |
| 503 | ✗ | const CloudTransportShape& shape = options.shape; | |
| 504 | ✗ | CacheKey key; | |
| 505 | ✗ | key.Add("ernie-cloud-transport-1"); | |
| 506 | ✗ | for (const std::uint64_t value : | |
| 507 | ✗ | {std::uint64_t{shape.lightCosines}, std::uint64_t{shape.thicknesses}, | |
| 508 | ✗ | std::uint64_t{shape.depths}, std::uint64_t{shape.viewCosines}, | |
| 509 | ✗ | std::uint64_t{shape.azimuths}, std::uint64_t{options.photons}, options.seed}) { | |
| 510 | ✗ | key.Add(std::to_string(value)); | |
| 511 | } | ||
| 512 | ✗ | const std::optional<std::filesystem::path> path = CloudCachePath(key); | |
| 513 | ✗ | if (path) { | |
| 514 | ✗ | if (std::optional<std::vector<float>> values = | |
| 515 | ✗ | ReadCachedFloats(*path, kCloudOrderSets * shape.Cells() * 4)) { | |
| 516 | ✗ | return {shape, std::move(*values), {}}; | |
| 517 | ✗ | } | |
| 518 | } | ||
| 519 | ✗ | CloudTransportTable table = | |
| 520 | ✗ | SimulateCloudTransport(DropletPhase(7.0, 2.0, {400.0, 700.0, 8}, 1024), options); | |
| 521 | ✗ | table.single.clear(); | |
| 522 | ✗ | if (path) { | |
| 523 | ✗ | WriteCachedFloats(*path, table.values); | |
| 524 | } | ||
| 525 | ✗ | return table; | |
| 526 | ✗ | } | |
| 527 | |||
| 528 | // ---------------------------------------------------------------------------------------------- | ||
| 529 | // Noise textures | ||
| 530 | |||
| 531 | 1 | CloudNoise MakeCloudNoise(std::size_t shapeSize, std::size_t detailSize, std::size_t weatherSize, | |
| 532 | std::size_t turbulenceSize) { | ||
| 533 | 1 | CloudNoise noise{shapeSize, {}, detailSize, {}, weatherSize, {}, turbulenceSize, {}}; | |
| 534 | 1 | noise.shape = MakeShapeNoise(noise.shapeSize); | |
| 535 | 1 | noise.detail = MakeDetailNoise(noise.detailSize); | |
| 536 | 1 | noise.weather = MakeWeatherMap(noise.weatherSize); | |
| 537 | 1 | noise.turbulence = MakeTurbulenceMap(noise.turbulenceSize); | |
| 538 | 1 | return noise; | |
| 539 | ✗ | } | |
| 540 | |||
| 541 | ✗ | CloudNoise CachedCloudNoise() { | |
| 542 | ✗ | CloudNoise noise; | |
| 543 | ✗ | std::vector<float>* parts[] = {&noise.shape, &noise.detail, &noise.weather, &noise.turbulence}; | |
| 544 | ✗ | const std::size_t sizes[] = {noise.shapeSize * noise.shapeSize * noise.shapeSize, | |
| 545 | noise.detailSize * noise.detailSize * noise.detailSize, | ||
| 546 | noise.weatherSize * noise.weatherSize * 4, | ||
| 547 | ✗ | noise.turbulenceSize * noise.turbulenceSize * 3}; | |
| 548 | ✗ | std::size_t total = 0; | |
| 549 | ✗ | for (const std::size_t size : sizes) { | |
| 550 | ✗ | total += size; | |
| 551 | } | ||
| 552 | ✗ | CacheKey key; | |
| 553 | ✗ | key.Add("ernie-cloud-noise-1"); | |
| 554 | ✗ | const std::optional<std::filesystem::path> path = CloudCachePath(key); | |
| 555 | ✗ | if (path) { | |
| 556 | // The four textures one after another in one file. | ||
| 557 | ✗ | if (const std::optional<std::vector<float>> values = ReadCachedFloats(*path, total)) { | |
| 558 | std::size_t offset = 0; | ||
| 559 | ✗ | for (std::size_t p = 0; p < 4; ++p) { | |
| 560 | ✗ | parts[p]->assign(values->begin() + static_cast<std::ptrdiff_t>(offset), | |
| 561 | ✗ | values->begin() + static_cast<std::ptrdiff_t>(offset + sizes[p])); | |
| 562 | ✗ | offset += sizes[p]; | |
| 563 | } | ||
| 564 | ✗ | return noise; | |
| 565 | ✗ | } | |
| 566 | } | ||
| 567 | ✗ | noise = MakeCloudNoise(); | |
| 568 | ✗ | if (path) { | |
| 569 | ✗ | std::vector<float> values; | |
| 570 | ✗ | values.reserve(total); | |
| 571 | ✗ | for (const std::vector<float>* part : parts) { | |
| 572 | ✗ | values.insert(values.end(), part->begin(), part->end()); | |
| 573 | } | ||
| 574 | ✗ | WriteCachedFloats(*path, values); | |
| 575 | ✗ | } | |
| 576 | return noise; | ||
| 577 | ✗ | } | |
| 578 | |||
| 579 | } // namespace ernie | ||
| 580 |