layers/sand_shadow.cpp
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | #include "terraingen/layers/sand_shadow.h" | ||
| 2 | |||
| 3 | #include "terraingen/raster/cell_loops.h" | ||
| 4 | |||
| 5 | #include <algorithm> | ||
| 6 | #include <cmath> | ||
| 7 | #include <numbers> | ||
| 8 | #include <stdexcept> | ||
| 9 | #include <utility> | ||
| 10 | |||
| 11 | namespace terraingen { | ||
| 12 | |||
| 13 | 24 | void ValidateSandShadowOptions(const SandShadowOptions& options) { | |
| 14 | 24 | if (!std::isfinite(options.windAzimuth)) { | |
| 15 | 1 | throw std::invalid_argument("Sand shadow windAzimuth must be finite"); | |
| 16 | } | ||
| 17 | 23 | if (!(options.lengthRatio > 0.0 && std::isfinite(options.lengthRatio))) { | |
| 18 | 3 | throw std::invalid_argument("Sand shadow lengthRatio must be positive and finite"); | |
| 19 | } | ||
| 20 | 20 | if (!(options.flankAngle > 0.0 && options.flankAngle < 0.5 * std::numbers::pi)) { | |
| 21 | 2 | throw std::invalid_argument("Sand shadow flankAngle must be in (0, 90) deg"); | |
| 22 | } | ||
| 23 | 18 | if (!(options.minDiameter >= 0.0 && std::isfinite(options.minDiameter))) { | |
| 24 | 1 | throw std::invalid_argument("Sand shadow minDiameter must be finite and not negative"); | |
| 25 | } | ||
| 26 | 17 | if (!(options.coverThickness > 0.0 && std::isfinite(options.coverThickness))) { | |
| 27 | 1 | throw std::invalid_argument("Sand shadow coverThickness must be positive and finite"); | |
| 28 | } | ||
| 29 | 16 | ValidateResolutionFade(options.fade); | |
| 30 | 15 | } | |
| 31 | |||
| 32 | 21 | SandShadowLayer::SandShadowLayer(std::shared_ptr<const BoulderLayer> boulders, | |
| 33 | 21 | SandShadowOptions options, Surface surface) | |
| 34 | 21 | : boulders_(std::move(boulders)), options_(options), surface_(surface) { | |
| 35 | 21 | ValidateSandShadowOptions(options_); | |
| 36 | 13 | if (!boulders_) { | |
| 37 | 1 | throw std::invalid_argument("Sand shadows need a boulder layer"); | |
| 38 | } | ||
| 39 | 12 | downwind_ = {std::cos(options_.windAzimuth), std::sin(options_.windAzimuth)}; | |
| 40 | 12 | flankSlope_ = std::tan(options_.flankAngle); | |
| 41 | 48 | } | |
| 42 | |||
| 43 | std::vector<SandShadowLayer::Shadow> | ||
| 44 | 62232 | SandShadowLayer::ShadowsOver(const Bounds2d& region, double spacing, const Ground* ground) const { | |
| 45 | 62232 | const double longest = boulders_->MaxFeatureDiameter() * (0.5 + options_.lengthRatio); | |
| 46 | 62232 | std::vector<Shadow> shadows; | |
| 47 | 62232 | for (const Boulder& boulder : | |
| 48 | 294909 | boulders_->Boulders(surface_.Expand(region, longest), spacing, ground)) { | |
| 49 | 170445 | if (boulder.diameter < options_.minDiameter) { | |
| 50 | 242 | continue; | |
| 51 | } | ||
| 52 | 170203 | const double weight = ResolutionWeight(boulder.diameter, spacing, options_.fade); | |
| 53 | 170203 | const double cap = boulder.ExposedHeight(); | |
| 54 | 170203 | if (weight <= 0.0 || cap <= 0.0) { | |
| 55 | 1 | continue; | |
| 56 | } | ||
| 57 | 170202 | const double reach = boulder.diameter * (0.5 + options_.lengthRatio); | |
| 58 | 170202 | shadows.push_back({{boulder.x, boulder.y}, | |
| 59 | 170202 | boulder.diameter, | |
| 60 | cap, | ||
| 61 | weight, | ||
| 62 | 340404 | surface_.Expand({boulder.x, boulder.y, boulder.x, boulder.y}, reach)}); | |
| 63 | 62232 | } | |
| 64 | 62232 | return shadows; | |
| 65 | } | ||
| 66 | |||
| 67 | 132575 | double SandShadowLayer::DepthOf(const Shadow& shadow, double x, double y) const { | |
| 68 | 132575 | const Vec2d offset = surface_.Offset(shadow.centre, {x, y}); | |
| 69 | 132575 | const double along = offset.x * downwind_.x + offset.y * downwind_.y - 0.5 * shadow.width; | |
| 70 | 132575 | const double length = options_.lengthRatio * shadow.width; | |
| 71 | 132575 | if (!(along >= 0.0 && along < length)) { | |
| 72 | 82418 | return 0.0; | |
| 73 | } | ||
| 74 | 50157 | const double across = std::fabs(offset.y * downwind_.x - offset.x * downwind_.y); | |
| 75 | 50157 | const double taper = 1.0 - along / length; | |
| 76 | 50157 | const double half = 0.5 * shadow.width * taper; | |
| 77 | 50157 | if (across >= half) { | |
| 78 | 44025 | return 0.0; | |
| 79 | } | ||
| 80 | 6132 | const double crest = std::min(0.5 * shadow.width * flankSlope_, shadow.cap) * taper; | |
| 81 | 6132 | return shadow.weight * crest * (1.0 - across / half); | |
| 82 | } | ||
| 83 | |||
| 84 | 60690 | double SandShadowLayer::ThicknessAt(double x, double y, double spacing, | |
| 85 | const Ground* ground) const { | ||
| 86 | 60690 | double depth = 0.0; | |
| 87 | 229174 | for (const Shadow& shadow : ShadowsOver({x, y, x, y}, spacing, ground)) { | |
| 88 | 168484 | if (surface_.Contains(shadow.box, {x, y})) { | |
| 89 | 6805 | depth = std::max(depth, DepthOf(shadow, x, y)); | |
| 90 | } | ||
| 91 | 60690 | } | |
| 92 | 60690 | return depth; | |
| 93 | } | ||
| 94 | |||
| 95 | 60674 | double SandShadowLayer::MaterialCoverAt(double x, double y, double spacing, | |
| 96 | const Ground* ground) const { | ||
| 97 | 60674 | return Smoothstep(ThicknessAt(x, y, spacing, ground) / options_.coverThickness); | |
| 98 | } | ||
| 99 | |||
| 100 | 3 | Raster<double> SandShadowLayer::MaterialCoverOver(const RasterGrid& grid, double spacing, | |
| 101 | const Ground* ground, unsigned threads) const { | ||
| 102 | 3 | Raster<double> cover(grid.width, grid.height); | |
| 103 | 3 | cover.SetTransform(grid.transform); | |
| 104 | 3 | const std::vector<Shadow> shadows = ShadowsOver(grid.Bounds(), spacing, ground); | |
| 105 | 3 | std::vector<std::vector<PixelWindow>> windows; | |
| 106 | 318 | for (const Shadow& shadow : shadows) { | |
| 107 | 315 | windows.push_back(PixelWindowsIn(grid, shadow.box, surface_)); | |
| 108 | } | ||
| 109 | 3 | ForEachFeatureCell( | |
| 110 | cover, threads, shadows.size(), | ||
| 111 | 1937 | [&](std::size_t n) -> const std::vector<PixelWindow>& { return windows[n]; }, | |
| 112 | 3 | [&](double& cell, const FeatureCell& at) { | |
| 113 | 66293 | cell = std::max(cell, DepthOf(shadows[at.feature], at.centre.x, at.centre.y)); | |
| 114 | 66293 | }); | |
| 115 | 548 | for (std::size_t row = 0; row < grid.height; ++row) { | |
| 116 | 127266 | for (std::size_t col = 0; col < grid.width; ++col) { | |
| 117 | 126721 | cover(col, row) = Smoothstep(cover(col, row) / options_.coverThickness); | |
| 118 | } | ||
| 119 | } | ||
| 120 | 6 | return cover; | |
| 121 | 3 | } | |
| 122 | |||
| 123 | 1539 | void SandShadowLayer::ApplySamples(SampleBlock& samples, const LayerContext& context) const { | |
| 124 | 1539 | Raster<float>* heights = samples.Grid(); | |
| 125 | const double spacing = | ||
| 126 | heights != nullptr | ||
| 127 | 1539 | ? std::max(context.spacing, RasterGrid::FromRaster(*heights).Spacing(surface_)) | |
| 128 | 1539 | : context.spacing; | |
| 129 | 1539 | const std::vector<Shadow> shadows = ShadowsOver(samples.Bounds(), spacing, context.ground); | |
| 130 | 1539 | if (heights != nullptr) { | |
| 131 | 2 | const std::span<const float> values = std::as_const(*heights).Values(); | |
| 132 | 2 | const std::vector<float> beneath(values.begin(), values.end()); | |
| 133 | 2 | const std::size_t width = heights->Width(); | |
| 134 | 2 | std::vector<std::vector<PixelWindow>> windows; | |
| 135 | 188 | for (const Shadow& shadow : shadows) { | |
| 136 | 186 | windows.push_back(PixelWindowsIn(*heights, shadow.box, surface_)); | |
| 137 | } | ||
| 138 | 4 | ForEachFeatureCell( | |
| 139 | 2 | *heights, context.threads, shadows.size(), | |
| 140 | 1666 | [&](std::size_t n) -> const std::vector<PixelWindow>& { return windows[n]; }, | |
| 141 | 2 | [&](float& cell, const FeatureCell& at) { | |
| 142 | 59153 | const double depth = DepthOf(shadows[at.feature], at.centre.x, at.centre.y); | |
| 143 | 59153 | const double base = static_cast<double>(beneath[at.row * width + at.col]); | |
| 144 | 59153 | cell = static_cast<float>(std::max(static_cast<double>(cell), base + depth)); | |
| 145 | 59153 | }); | |
| 146 | 2 | return; | |
| 147 | 2 | } | |
| 148 | 3074 | for (std::size_t i = 0; i < samples.Size(); ++i) { | |
| 149 | 1537 | if (samples.IsNoData(i)) { | |
| 150 | 1 | continue; | |
| 151 | } | ||
| 152 | 1536 | const Vec2d p = samples.Position(i); | |
| 153 | 1536 | double depth = 0.0; | |
| 154 | 2752 | for (const Shadow& shadow : shadows) { | |
| 155 | 1216 | if (surface_.Contains(shadow.box, p)) { | |
| 156 | 324 | depth = std::max(depth, DepthOf(shadow, p.x, p.y)); | |
| 157 | } | ||
| 158 | } | ||
| 159 | 1536 | samples.Height(i) = static_cast<float>(static_cast<double>(samples.Height(i)) + depth); | |
| 160 | } | ||
| 161 | 1539 | } | |
| 162 | |||
| 163 | } // namespace terraingen | ||
| 164 |