layers/craters/crater_layer.cpp
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | #include "terraingen/layers/craters/crater_layer.h" | ||
| 2 | |||
| 3 | #include "terraingen/distribution/summed.h" | ||
| 4 | #include "terraingen/layers/craters/generation/spacing.h" | ||
| 5 | #include "terraingen/raster/grid.h" | ||
| 6 | #include "terraingen/util/cells.h" | ||
| 7 | |||
| 8 | #include <algorithm> | ||
| 9 | #include <cmath> | ||
| 10 | #include <optional> | ||
| 11 | #include <stdexcept> | ||
| 12 | #include <utility> | ||
| 13 | #include <vector> | ||
| 14 | |||
| 15 | namespace terraingen { | ||
| 16 | |||
| 17 | namespace { | ||
| 18 | |||
| 19 | /// The largest primary (m) searched for the smallest one throwing secondaries a map measures, | ||
| 20 | /// about where Singer et al.'s largest secondary stops growing with the primary. | ||
| 21 | constexpr double kLargestThrower = 5e5; | ||
| 22 | |||
| 23 | } // namespace | ||
| 24 | |||
| 25 | 697 | CraterLayer::CraterLayer(CraterModels models, CraterLayerOptions options) | |
| 26 | 697 | : models_(std::move(models)), options_(std::move(options)) { | |
| 27 | 697 | if (!models_.density || !models_.ages || !models_.degradation || !models_.profile) { | |
| 28 | throw std::invalid_argument( | ||
| 29 | 1 | "A crater layer needs a density, ages, degradation and profile"); | |
| 30 | } | ||
| 31 | 696 | ValidateCraterSpacing(options_.generation.spacing); | |
| 32 | 695 | ValidateCraterApplyOptions(options_.apply); | |
| 33 | 694 | options_.generation.marginRadii = CraterReachRadii(*models_.profile, options_.apply) * | |
| 34 | 694 | MaxOutlineScale(options_.apply.irregularity); | |
| 35 | 694 | options_.generation.surface = options_.surface; | |
| 36 | 694 | options_.apply.surface = options_.surface; | |
| 37 | 694 | if (options_.surface.IsBody()) { | |
| 38 | 32 | cells_ = std::make_shared<CubeSphereCells>(options_.surface); | |
| 39 | } else { | ||
| 40 | 662 | cells_ = std::make_shared<PlanarCells>(); | |
| 41 | } | ||
| 42 | 694 | if (options_.cacheCraters > 0) { | |
| 43 | 693 | cache_ = std::make_shared<CraterCellCache>(options_.cacheCraters); | |
| 44 | } | ||
| 45 | 704 | for (const SecondaryCraterOptions& secondaries : options_.secondaries) { | |
| 46 | 10 | ValidateSecondaryCraterOptions(secondaries); | |
| 47 | } | ||
| 48 | 694 | candidates_ = models_.density; | |
| 49 | 694 | if (!options_.densityMap) { | |
| 50 | 688 | return; | |
| 51 | } | ||
| 52 | 6 | scale_ = std::make_shared<const DensityScale>(*options_.densityMap, *models_.density); | |
| 53 | 6 | if (scale_->Max() > 1.0) { | |
| 54 | candidates_ = | ||
| 55 | 20 | std::make_shared<SummedSfd>(std::vector<WeightedSfd>{{scale_->Max(), models_.density}}); | |
| 56 | } | ||
| 57 | 12 | options_.generation.keepShare = [scale = scale_](const Crater& crater) { | |
| 58 | 17758 | return scale->At(crater.x, crater.y) / scale->Max(); | |
| 59 | 6 | }; | |
| 60 | 6 | const DensityMap& map = *options_.densityMap; | |
| 61 | 6 | if (!map.total || map.measure == DensityMeasure::Scale) { | |
| 62 | 2 | return; | |
| 63 | } | ||
| 64 | 6 | for (const SecondaryCraterOptions& secondaries : options_.secondaries) { | |
| 65 | 2 | const double smallest = std::fmax(map.minDiameter, secondaries.minDiameter); | |
| 66 | 202 | const auto throws = [&](double diameter) { | |
| 67 | 202 | return LargestSecondary(diameter, secondaries) > smallest; | |
| 68 | 2 | }; | |
| 69 | 2 | double low = std::log(1e-3); | |
| 70 | 2 | double high = std::log(kLargestThrower); | |
| 71 | 2 | if (!throws(kLargestThrower)) { | |
| 72 | ✗ | continue; | |
| 73 | } | ||
| 74 | 202 | for (int i = 0; i < 100; ++i) { | |
| 75 | 200 | const double middle = 0.5 * (low + high); | |
| 76 | 200 | (throws(std::exp(middle)) ? high : low) = middle; | |
| 77 | } | ||
| 78 | 2 | throwFrom_ = | |
| 79 | 2 | std::fmin(throwFrom_, std::fmax(std::exp(high), secondaries.selection.minDiameter)); | |
| 80 | } | ||
| 81 | 4 | if (std::isfinite(throwFrom_) && options_.generation.maxCraters > 0) { | |
| 82 | throw std::invalid_argument("Taking secondaries out of a density map needs every crater, " | ||
| 83 | 1 | "not the maxCraters largest"); | |
| 84 | } | ||
| 85 | 3 | options_.generation.spacing.outrankFrom = throwFrom_; | |
| 86 | 42 | } | |
| 87 | |||
| 88 | 3 | CraterKeepShare CraterLayer::KeepShare(const Bounds2d& region) const { | |
| 89 | 3 | CraterGenerationOptions large = options_.generation; | |
| 90 | 3 | large.dMin = std::fmax(large.dMin, throwFrom_); | |
| 91 | 3 | large.minResolvedDiameter = 0.0; | |
| 92 | 3 | double reach = 0.0; | |
| 93 | 6 | for (const SecondaryCraterOptions& secondaries : options_.secondaries) { | |
| 94 | 3 | reach = std::fmax(reach, SecondaryOuterRange(MaxFeatureDiameter(), secondaries)); | |
| 95 | } | ||
| 96 | 3 | const Bounds2d around = options_.surface.Expand(region, reach); | |
| 97 | 3 | CraterList throwers; | |
| 98 | 3 | if (large.dMin < ResolveDMax(large.dMax, *candidates_)) { | |
| 99 | ✗ | if (cache_) { | |
| 100 | large.cellCraters = | ||
| 101 | ✗ | CachedCellCraters(*cache_, *cells_, *candidates_, *models_.ages, options_.seed); | |
| 102 | } | ||
| 103 | throwers = | ||
| 104 | ✗ | GenerateCraters(*cells_, around, *candidates_, *models_.ages, options_.seed, large); | |
| 105 | } else { | ||
| 106 | 3 | AppendManualCraters(large.manual, throwers, [&](const Crater& crater) { | |
| 107 | 3 | return options_.surface.Contains(around, {crater.x, crater.y}); | |
| 108 | }); | ||
| 109 | } | ||
| 110 | struct Thrower { | ||
| 111 | Vec2d centre; | ||
| 112 | SecondaryDensity density; | ||
| 113 | }; | ||
| 114 | 3 | auto thrown = std::make_shared<std::vector<Thrower>>(); | |
| 115 | 5 | for (std::size_t n = 0; n < throwers.Size(); ++n) { | |
| 116 | 2 | const Crater crater = throwers.At(n); | |
| 117 | 4 | for (const SecondaryCraterOptions& secondaries : options_.secondaries) { | |
| 118 | 2 | SecondaryDensity density(crater, secondaries); | |
| 119 | 2 | if (density.Reach() > 0.0) { | |
| 120 | 4 | thrown->push_back({{crater.x, crater.y}, std::move(density)}); | |
| 121 | } | ||
| 122 | 2 | } | |
| 123 | } | ||
| 124 | 3 | const DensityMap& map = *options_.densityMap; | |
| 125 | 3 | const bool cover = map.measure == DensityMeasure::ArealFraction; | |
| 126 | 6 | return [scale = scale_, map, cover, thrown, from = throwFrom_, | |
| 127 | 3 | surface = options_.surface](const Crater& crater) { | |
| 128 | 513392 | const double share = scale->At(crater.x, crater.y) / scale->Max(); | |
| 129 | 513392 | if (crater.diameter >= from) { | |
| 130 | ✗ | return share; | |
| 131 | } | ||
| 132 | 513392 | double value = MapValueAt(map, crater.x, crater.y); | |
| 133 | 513392 | if (std::isnan(value)) { | |
| 134 | ✗ | return share; | |
| 135 | } | ||
| 136 | 1022941 | for (const Thrower& thrower : *thrown) { | |
| 137 | 509549 | const Vec2d offset = surface.Offset(thrower.centre, {crater.x, crater.y}); | |
| 138 | 509549 | if (std::hypot(offset.x, offset.y) < thrower.density.Reach()) { | |
| 139 | 559128 | value -= cover ? thrower.density.CoverAt(offset, map.minDiameter, map.maxDiameter) | |
| 140 | 279564 | : thrower.density.At(offset, map.minDiameter, map.maxDiameter); | |
| 141 | } | ||
| 142 | } | ||
| 143 | 513392 | return scale->Scale(value) / scale->Max(); | |
| 144 | 6 | }; | |
| 145 | 5 | } | |
| 146 | |||
| 147 | 253829 | CraterList CraterLayer::Craters(const Bounds2d& region, double spacing, double minDiameter, | |
| 148 | double maxDiameter) const { | ||
| 149 | 253829 | CraterGenerationOptions generation = options_.generation; | |
| 150 | 253829 | if (options_.generate) { | |
| 151 | 219444 | if (spacing > 0.0) { | |
| 152 | 219423 | generation.minResolvedDiameter = options_.apply.fade.startSamples * spacing; | |
| 153 | } | ||
| 154 | 219444 | generation.minResolvedDiameter = std::fmax(generation.minResolvedDiameter, minDiameter); | |
| 155 | 219444 | generation.maxBandDiameter = std::fmin(generation.maxBandDiameter, maxDiameter); | |
| 156 | 219444 | if (cache_) { | |
| 157 | generation.cellCraters = | ||
| 158 | 219441 | CachedCellCraters(*cache_, *cells_, *candidates_, *models_.ages, options_.seed); | |
| 159 | } | ||
| 160 | 219444 | if (std::isfinite(throwFrom_)) { | |
| 161 | 3 | generation.keepShare = KeepShare(region); | |
| 162 | } | ||
| 163 | 219444 | return GenerateCraters(*cells_, region, *candidates_, *models_.ages, options_.seed, | |
| 164 | 219444 | generation); | |
| 165 | } | ||
| 166 | 34385 | CraterList craters; | |
| 167 | 34388 | AppendManualCraters(generation.manual, craters, [&](const Crater& crater) { | |
| 168 | 82908 | const double margin = 0.5 * generation.marginRadii * crater.diameter; | |
| 169 | 165816 | return options_.surface.Contains(options_.surface.Expand(region, margin), | |
| 170 | 248724 | {crater.x, crater.y}); | |
| 171 | }); | ||
| 172 | 34382 | return craters; | |
| 173 | 253832 | } | |
| 174 | |||
| 175 | 244606 | void CraterLayer::ApplySamples(SampleBlock& samples, const LayerContext& context) const { | |
| 176 | 244606 | const Raster<float>* grid = samples.Grid(); | |
| 177 | 244606 | if (grid == nullptr && options_.apply.datum != CraterDatum::None) { | |
| 178 | 5002 | RequireGround(context.ground, "Point queries of craters with a datum"); | |
| 179 | } | ||
| 180 | 244605 | if (grid == nullptr && options_.generation.maxCraters > 0) { | |
| 181 | throw std::invalid_argument( | ||
| 182 | 1 | "Point queries cannot use maxCraters, which ranks craters over a whole grid"); | |
| 183 | } | ||
| 184 | const double spacing = | ||
| 185 | grid != nullptr | ||
| 186 | 244604 | ? std::max(context.spacing, RasterGrid::FromRaster(*grid).Spacing(options_.surface)) | |
| 187 | 244604 | : context.spacing; | |
| 188 | 244604 | const Bounds2d bounds = samples.Bounds(); | |
| 189 | 244604 | const double margin = DatumMargin(options_.apply, *models_.profile, MaxFeatureDiameter()); | |
| 190 | 244604 | const CraterList craters = Craters(options_.surface.Expand(bounds, margin), spacing); | |
| 191 | const HeightSampler ground = | ||
| 192 | 244604 | context.ground ? context.ground->Sampler(spacing) : HeightSampler{}; | |
| 193 | 244604 | CraterSampling sampling{spacing, context.threads, ground}; | |
| 194 | 244604 | sampling.cache = context.cache; | |
| 195 | 244604 | sampling.owner = this; | |
| 196 | 244604 | sampling.backend = context.backend; | |
| 197 | 244604 | ApplyCraters(samples, craters, models_.Shape(), options_.apply, sampling); | |
| 198 | 244604 | } | |
| 199 | |||
| 200 | 63437 | CraterAttributes CraterLayer::Attributes(const Crater& crater) const { | |
| 201 | 63437 | return ComputeCraterAttributes(crater, models_.degradation->State(crater), *models_.profile); | |
| 202 | } | ||
| 203 | |||
| 204 | 121 | CraterFeatures CraterLayer::CratersOver(const RasterGrid& grid, double spacing) const { | |
| 205 | 121 | CraterFeatures placed{Craters(grid.Bounds(), std::max(spacing, grid.Spacing(options_.surface))), | |
| 206 | 242 | {}}; | |
| 207 | 61372 | for (std::size_t n = 0; n < placed.craters.Size(); ++n) { | |
| 208 | 61251 | placed.attributes.push_back(Attributes(placed.craters.At(n))); | |
| 209 | } | ||
| 210 | 121 | return placed; | |
| 211 | } | ||
| 212 | |||
| 213 | 253431 | double CraterLayer::MaxFeatureDiameter() const { | |
| 214 | 253431 | const CraterGenerationOptions& generation = options_.generation; | |
| 215 | 253431 | double largest = options_.generate ? ResolveDMax(generation.dMax, *models_.density) : 0.0; | |
| 216 | 339003 | for (std::size_t n = 0; n < generation.manual.Size(); ++n) { | |
| 217 | 85572 | largest = std::fmax(largest, generation.manual.At(n).diameter); | |
| 218 | } | ||
| 219 | 253431 | return largest; | |
| 220 | } | ||
| 221 | |||
| 222 | } // namespace terraingen | ||
| 223 |