layers/craters/shaping/apply.cpp
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | #include "terraingen/layers/craters/shaping/apply.h" | ||
| 2 | |||
| 3 | #include "terraingen/gpu/kernel.h" | ||
| 4 | #include "terraingen/layers/craters/shaping/crater_kernel.h" | ||
| 5 | #include "terraingen/noise/fractal.h" | ||
| 6 | #include "terraingen/noise/hash.h" | ||
| 7 | #include "terraingen/raster/cell_loops.h" | ||
| 8 | #include "terraingen/raster/grid.h" | ||
| 9 | #include "terraingen/util/checks.h" | ||
| 10 | #include "terraingen/util/log.h" | ||
| 11 | #include "terraingen/util/parallel.h" | ||
| 12 | |||
| 13 | #include <algorithm> | ||
| 14 | #include <array> | ||
| 15 | #include <cmath> | ||
| 16 | #include <cstdint> | ||
| 17 | #include <limits> | ||
| 18 | #include <map> | ||
| 19 | #include <numbers> | ||
| 20 | #include <optional> | ||
| 21 | #include <stdexcept> | ||
| 22 | #include <utility> | ||
| 23 | #include <vector> | ||
| 24 | |||
| 25 | namespace terraingen { | ||
| 26 | |||
| 27 | namespace { | ||
| 28 | |||
| 29 | /// Nodes and weights of 9-point Gauss-Hermite quadrature (the node at 0 and the positive half): | ||
| 30 | /// the mean of f over a normal distribution of mean m and deviation s is | ||
| 31 | /// sum(w f(m + sqrt(2) s x)) / sqrt(pi). | ||
| 32 | constexpr std::array<double, 5> kHermiteNodes{0.0, 0.723551018752838, 1.46855328921667, | ||
| 33 | 2.26658058453184, 3.19099320178153}; | ||
| 34 | constexpr std::array<double, 5> kHermiteWeights{0.720235215606051, 0.432651559002556, | ||
| 35 | 0.0884745273943766, 0.00494362427553695, | ||
| 36 | 0.0000396069772632645}; | ||
| 37 | |||
| 38 | /// The profile's height at `distance`, its rim rounded by `rim`: blended by `rim.Weight` toward | ||
| 39 | /// the profile averaged over a normal spread of distances of deviation `sigmaRadii` radii. | ||
| 40 | 790043 | double ProfileHeight(const CraterProfile& profile, double distance, double diameter, | |
| 41 | const CraterState& state, const RimSmoothing& rim) { | ||
| 42 | 790043 | const double height = profile.Height(distance, diameter, state); | |
| 43 | 790043 | const double radius = 0.5 * diameter; | |
| 44 | 790043 | const double weight = rim.Weight(distance / radius); | |
| 45 | 790043 | if (weight == 0.0) { | |
| 46 | 775174 | return height; | |
| 47 | } | ||
| 48 | 14869 | const double step = std::numbers::sqrt2 * rim.sigmaRadii * radius; | |
| 49 | 14869 | double smoothed = kHermiteWeights[0] * height; | |
| 50 | 148690 | for (std::size_t i = 1; i < kHermiteNodes.size(); ++i) { | |
| 51 | 178428 | for (const double sign : {-1.0, 1.0}) { | |
| 52 | 118952 | const double at = std::fabs(distance + sign * step * kHermiteNodes[i]); | |
| 53 | 118952 | smoothed += kHermiteWeights[i] * profile.Height(at, diameter, state); | |
| 54 | } | ||
| 55 | } | ||
| 56 | 14869 | smoothed /= std::sqrt(std::numbers::pi); | |
| 57 | 14869 | return height + weight * (smoothed - height); | |
| 58 | } | ||
| 59 | |||
| 60 | /// The noise `noise` adds to `crater` at `offset` metres from its centre and `radii` crater radii | ||
| 61 | /// out: footprint noise fading to nothing from the rim to `reachRadii`, and rim noise in its band. | ||
| 62 | 790043 | double NoiseHeight(const Crater& crater, const Vec2d& offset, double radii, double reachRadii, | |
| 63 | const CraterNoise& noise) { | ||
| 64 | 790043 | if (noise.amplitude == 0.0 && noise.rimAmplitude == 0.0) { | |
| 65 | 757775 | return 0.0; | |
| 66 | } | ||
| 67 | 32268 | const double radius = 0.5 * crater.diameter; | |
| 68 | const std::uint32_t key = | ||
| 69 | 32268 | PcgHash(static_cast<std::uint32_t>(crater.id) ^ | |
| 70 | 32268 | PcgHash(static_cast<std::uint32_t>(crater.id >> 32) ^ noise.seed)); | |
| 71 | 37668 | const auto fractal = [&](double wavelengthRadii, std::uint32_t seed) { | |
| 72 | 37668 | FractalParams params; | |
| 73 | 37668 | params.seed = seed; | |
| 74 | 37668 | params.wavelength = static_cast<float>(wavelengthRadii * radius); | |
| 75 | 37668 | params.amplitude = 1.0f; | |
| 76 | 37668 | params.octaves = noise.octaves; | |
| 77 | 37668 | params.gain = static_cast<float>(noise.gain); | |
| 78 | return static_cast<double>( | ||
| 79 | 75336 | Fractal2(static_cast<float>(offset.x), static_cast<float>(offset.y), params)); | |
| 80 | 32268 | }; | |
| 81 | const double footprint = | ||
| 82 | 32268 | radii <= 1.0 ? 1.0 : Falloff(radii, 1.0, reachRadii, BlendShape::Smootherstep); | |
| 83 | const double rim = | ||
| 84 | radii <= 1.0 | ||
| 85 | 32268 | ? Falloff(2.0 - radii, 1.0, 2.0 - noise.rimInnerRadii, BlendShape::Smootherstep) | |
| 86 | 18768 | : Falloff(radii, 1.0, noise.rimOuterRadii, BlendShape::Smootherstep); | |
| 87 | 32268 | double height = 0.0; | |
| 88 | 32268 | if (noise.amplitude > 0.0 && footprint > 0.0) { | |
| 89 | 32268 | height += noise.amplitude * footprint * fractal(noise.wavelengthRadii, key); | |
| 90 | } | ||
| 91 | 32268 | if (noise.rimAmplitude > 0.0 && rim > 0.0) { | |
| 92 | 5400 | height += noise.rimAmplitude * rim * fractal(noise.rimWavelengthRadii, key + 1u); | |
| 93 | } | ||
| 94 | 32268 | return crater.diameter * height; | |
| 95 | } | ||
| 96 | |||
| 97 | /// Fractal noise of `octaves` octaves streaked radially about a crater: `across` radii wide around | ||
| 98 | /// it at the rim and `along` radii long outward, at `radii` radii and `angle` radians from its | ||
| 99 | /// centre. | ||
| 100 | 37152 | double RadialStreaks(double radii, double angle, double across, double along, int octaves, | |
| 101 | std::uint32_t seed) { | ||
| 102 | 37152 | FractalParams params; | |
| 103 | 37152 | params.seed = seed; | |
| 104 | 37152 | params.wavelength = 1.0f; | |
| 105 | 37152 | params.amplitude = 1.0f; | |
| 106 | 37152 | params.octaves = octaves; | |
| 107 | 37152 | const double ring = 1.0 / across; | |
| 108 | return static_cast<double>( | ||
| 109 | 74304 | Fractal3(ring * std::cos(angle), ring * std::sin(angle), radii / along, params)); | |
| 110 | } | ||
| 111 | |||
| 112 | /// The texture `texture` adds to `crater` in `state` at `offset` metres from its centre and | ||
| 113 | /// `radii` crater radii out: gullies down the wall and ejecta texture, faded with wear. | ||
| 114 | 790043 | double TextureHeight(const Crater& crater, const CraterState& state, const Vec2d& offset, | |
| 115 | double radii, const CraterTexture& texture) { | ||
| 116 | 790043 | if (texture.wallStreaks == 0.0 && texture.ejecta == 0.0) { | |
| 117 | 749543 | return 0.0; | |
| 118 | } | ||
| 119 | 40500 | const double wear = (state.Freshness() - texture.minFreshness) / (1.0 - texture.minFreshness); | |
| 120 | 40500 | const double fade = Smoothstep(std::clamp(wear, 0.0, 1.0)); | |
| 121 | 40500 | if (!(fade > 0.0)) { | |
| 122 | ✗ | return 0.0; | |
| 123 | } | ||
| 124 | const std::uint32_t key = | ||
| 125 | 40500 | PcgHash(static_cast<std::uint32_t>(crater.id) ^ | |
| 126 | 40500 | PcgHash(static_cast<std::uint32_t>(crater.id >> 32) ^ texture.seed ^ 0x54455854u)); | |
| 127 | 40500 | const double angle = std::atan2(offset.y, offset.x); | |
| 128 | 40500 | double height = 0.0; | |
| 129 | 40500 | if (texture.wallStreaks > 0.0 && radii > 0.25 && radii < 1.0) { | |
| 130 | 10152 | const double band = Smoothstep(std::clamp((radii - 0.25) / 0.35, 0.0, 1.0)) * | |
| 131 | 10152 | Smoothstep(std::clamp((1.0 - radii) / 0.1, 0.0, 1.0)); | |
| 132 | 20304 | height += texture.wallStreaks * crater.diameter * band * | |
| 133 | 10152 | RadialStreaks(radii, angle, texture.wallStreakWidthRadii, 0.5, 2, key); | |
| 134 | } | ||
| 135 | 40500 | if (texture.ejecta > 0.0 && radii > 1.0 && radii < texture.ejectaReachRadii) { | |
| 136 | 27000 | const double radius = 0.5 * crater.diameter; | |
| 137 | const double thickness = | ||
| 138 | 27000 | 0.14 * std::pow(radius, 0.74) * EjectaThicknessShare(radii, texture.ejectaReachRadii); | |
| 139 | 27000 | const double rise = Smoothstep(std::clamp((radii - 1.0) / 0.1, 0.0, 1.0)); | |
| 140 | 54000 | height += texture.ejecta * thickness * rise * | |
| 141 | 27000 | RadialStreaks(radii, angle, texture.ejectaLineationRadii, | |
| 142 | 27000 | 2.0 * texture.ejectaLineationRadii, 3, key + 1u); | |
| 143 | } | ||
| 144 | 40500 | return fade * height; | |
| 145 | } | ||
| 146 | |||
| 147 | /// Height of the raster at world point (x, y) by bilinear interpolation of the cell centres, or | ||
| 148 | /// nothing if a contributing cell is outside the raster or holds no data. | ||
| 149 | 151893 | std::optional<double> Sample(const Raster<float>& raster, double x, double y) { | |
| 150 | 151893 | const Vec2d pixel = raster.Transform().WorldToPixel(x, y); | |
| 151 | 151893 | const double u = pixel.x - 0.5; | |
| 152 | 151893 | const double v = pixel.y - 0.5; | |
| 153 | 151893 | const double col0 = std::floor(u); | |
| 154 | 151893 | const double row0 = std::floor(v); | |
| 155 | 151893 | const double maxCol = static_cast<double>(raster.Width()) - 1.0; | |
| 156 | 151893 | const double maxRow = static_cast<double>(raster.Height()) - 1.0; | |
| 157 | 151893 | if (col0 < 0.0 || row0 < 0.0 || col0 > maxCol || row0 > maxRow) { | |
| 158 | 39589 | return std::nullopt; | |
| 159 | } | ||
| 160 | 112304 | const auto c0 = static_cast<std::size_t>(col0); | |
| 161 | 112304 | const auto r0 = static_cast<std::size_t>(row0); | |
| 162 | 112304 | const std::size_t c1 = std::min(c0 + 1, raster.Width() - 1); | |
| 163 | 112304 | const std::size_t r1 = std::min(r0 + 1, raster.Height() - 1); | |
| 164 | 112304 | const double fu = u - col0; | |
| 165 | 112304 | const double fv = v - row0; | |
| 166 | 112304 | const std::array<float, 4> cells{raster(c0, r0), raster(c1, r0), raster(c0, r1), | |
| 167 | 112304 | raster(c1, r1)}; | |
| 168 | 112304 | std::array<double, 4> values{}; | |
| 169 | 673824 | for (std::size_t i = 0; i < cells.size(); ++i) { | |
| 170 | 449216 | if (raster.IsNoData(cells[i])) { | |
| 171 | ✗ | return std::nullopt; | |
| 172 | } | ||
| 173 | 449216 | values[i] = cells[i]; | |
| 174 | } | ||
| 175 | 112304 | return (1.0 - fv) * ((1.0 - fu) * values[0] + fu * values[1]) + | |
| 176 | 112304 | fv * ((1.0 - fu) * values[2] + fu * values[3]); | |
| 177 | } | ||
| 178 | |||
| 179 | /// The box a crater reaching `reachRadii` radii can change, in surface coordinates. | ||
| 180 | 5742809 | Bounds2d Footprint(const PreparedCrater& prepared, double reachRadii, const Surface& surface) { | |
| 181 | 5742809 | const Crater& crater = prepared.crater; | |
| 182 | 5742809 | const double reach = prepared.outline.MaxScale() * reachRadii * (0.5 * crater.diameter); | |
| 183 | 5742809 | return surface.Expand({crater.x, crater.y, crater.x, crater.y}, reach); | |
| 184 | } | ||
| 185 | |||
| 186 | /// Craters indexed for `KernelDatumRelief`, with the kernel that reads them. | ||
| 187 | /// | ||
| 188 | /// Records hold each crater's position as the float nearest it and `lows` what that float drops, | ||
| 189 | /// so the kernel takes differences between nearby craters exactly; bands and cells are at absolute | ||
| 190 | /// positions. The relief summed at a point therefore depends only on the craters reaching it, not | ||
| 191 | /// on which others an index holds, and grids and points fit the same datums. | ||
| 192 | struct DatumKernelIndex { | ||
| 193 | CraterDatumKernel kernel; | ||
| 194 | CraterKernelBatch batch; | ||
| 195 | std::vector<float> lows; | ||
| 196 | /// Per size band: its cell size, which is its largest reach, and its first cell's column and | ||
| 197 | /// row. | ||
| 198 | std::vector<float> bandFloats; | ||
| 199 | /// Per size band: its cells across and down and its first cell. | ||
| 200 | std::vector<std::uint32_t> bandCells; | ||
| 201 | /// Per cell: its first entry and its entry count. | ||
| 202 | std::vector<std::uint32_t> cells; | ||
| 203 | /// Indices into the records, ascending within each cell. | ||
| 204 | std::vector<std::uint32_t> entries; | ||
| 205 | /// The record of each prepared crater, for the craters indexed. | ||
| 206 | std::vector<std::uint32_t> slot; | ||
| 207 | }; | ||
| 208 | |||
| 209 | /// Most cells `IndexDatumRelief` lays out, and the largest cell column or row a float holds | ||
| 210 | /// exactly enough, before datums are left to the double path. | ||
| 211 | constexpr std::size_t kMaxReliefCells = std::size_t{1} << 24; | ||
| 212 | constexpr float kMaxReliefCellIndex = 8388608.0f; | ||
| 213 | |||
| 214 | /// The size band of `crater`: its diameter stretched by its outline, as a power of two. | ||
| 215 | 1434035 | int ReliefBand(const PreparedCrater& crater) { | |
| 216 | 1434035 | return std::ilogb(crater.crater.diameter * crater.outline.MaxScale()); | |
| 217 | } | ||
| 218 | |||
| 219 | /// Indexes the craters of `prepared` that can reach the datum footprint of any of `fitted`, in | ||
| 220 | /// size bands of cells as wide as each band's largest reach; nothing for craters spread too far. | ||
| 221 | 2722 | std::optional<DatumKernelIndex> IndexDatumRelief(const std::vector<PreparedCrater>& prepared, | |
| 222 | const std::vector<std::size_t>& fitted, | ||
| 223 | const CraterProfile& profile, | ||
| 224 | const CraterShapeKernel& shape, | ||
| 225 | const CraterApplyOptions& options) { | ||
| 226 | 2722 | constexpr double kInfinity = std::numeric_limits<double>::infinity(); | |
| 227 | 2722 | Bounds2d region{kInfinity, kInfinity, -kInfinity, -kInfinity}; | |
| 228 | 2722 | std::vector<bool> isFitted(prepared.size(), false); | |
| 229 | 133233 | for (const std::size_t k : fitted) { | |
| 230 | 130511 | const Crater& c = prepared[k].crater; | |
| 231 | 130511 | const double fitRadius = options.datumFitRadii * 0.5 * c.diameter; | |
| 232 | 130511 | region = {std::fmin(region.minX, c.x - fitRadius), std::fmin(region.minY, c.y - fitRadius), | |
| 233 | 130511 | std::fmax(region.maxX, c.x + fitRadius), std::fmax(region.maxY, c.y + fitRadius)}; | |
| 234 | 130511 | isFitted[k] = true; | |
| 235 | } | ||
| 236 | 2722 | const std::size_t last = *std::max_element(fitted.begin(), fitted.end()); | |
| 237 | 2722 | DatumKernelIndex index; | |
| 238 | 2722 | index.slot.assign(prepared.size(), 0); | |
| 239 | 2722 | std::vector<std::size_t> kept; | |
| 240 | 4384816 | for (std::size_t n = 0; n <= last; ++n) { | |
| 241 | 4382094 | const PreparedCrater& p = prepared[n]; | |
| 242 | 4382094 | const double reach = p.outline.MaxScale() * profile.ReachRadii() * 0.5 * p.crater.diameter; | |
| 243 | 11137317 | if (isFitted[n] || (p.crater.x + reach > region.minX && p.crater.x - reach < region.maxX && | |
| 244 | 6755223 | p.crater.y + reach > region.minY && p.crater.y - reach < region.maxY)) { | |
| 245 | 1434035 | index.slot[n] = static_cast<std::uint32_t>(kept.size()); | |
| 246 | 1434035 | kept.push_back(n); | |
| 247 | } | ||
| 248 | } | ||
| 249 | 2722 | index.batch = MakeCraterKernelBatch(shape, profile, prepared, kept, options); | |
| 250 | 2722 | if (index.batch.data.empty()) { | |
| 251 | 94 | index.batch.data.push_back(0.0f); | |
| 252 | } | ||
| 253 | 2722 | index.lows.resize(2 * kept.size()); | |
| 254 | 2722 | std::map<int, std::vector<std::uint32_t>> bands; | |
| 255 | 1436757 | for (std::size_t m = 0; m < kept.size(); ++m) { | |
| 256 | 1434035 | const Crater& c = prepared[kept[m]].crater; | |
| 257 | 1434035 | CraterRecord& record = index.batch.records[m]; | |
| 258 | 1434035 | record.x = static_cast<float>(c.x); | |
| 259 | 1434035 | record.y = static_cast<float>(c.y); | |
| 260 | 1434035 | index.lows[2 * m] = static_cast<float>(c.x - static_cast<double>(record.x)); | |
| 261 | 1434035 | index.lows[2 * m + 1] = static_cast<float>(c.y - static_cast<double>(record.y)); | |
| 262 | 1434035 | bands[ReliefBand(prepared[kept[m]])].push_back(static_cast<std::uint32_t>(m)); | |
| 263 | } | ||
| 264 | 2722 | std::vector<std::uint32_t> cellOf; | |
| 265 | 2722 | std::vector<std::uint32_t> filled; | |
| 266 | 15354 | for (const auto& [band, members] : bands) { | |
| 267 | 12632 | const auto size = static_cast<float>(1.0001 * profile.ReachRadii() * std::ldexp(1.0, band)); | |
| 268 | 12632 | float minI = std::numeric_limits<float>::infinity(); | |
| 269 | 12632 | float minJ = minI; | |
| 270 | 12632 | float maxI = -minI; | |
| 271 | 12632 | float maxJ = -minI; | |
| 272 | 1446667 | for (const std::uint32_t m : members) { | |
| 273 | 1434035 | const float i = std::floor(index.batch.records[m].x / size); | |
| 274 | 1434035 | const float j = std::floor(index.batch.records[m].y / size); | |
| 275 | 1434035 | minI = std::fmin(minI, i); | |
| 276 | 1434035 | minJ = std::fmin(minJ, j); | |
| 277 | 1434035 | maxI = std::fmax(maxI, i); | |
| 278 | 1434035 | maxJ = std::fmax(maxJ, j); | |
| 279 | } | ||
| 280 | 12632 | const std::size_t first = index.cells.size() / 2; | |
| 281 | 12632 | const auto across = static_cast<std::size_t>(maxI - minI) + 1; | |
| 282 | 12632 | const auto down = static_cast<std::size_t>(maxJ - minJ) + 1; | |
| 283 | 25264 | if (std::fmax(std::fabs(minI), std::fabs(maxI)) >= kMaxReliefCellIndex || | |
| 284 | 25264 | std::fmax(std::fabs(minJ), std::fabs(maxJ)) >= kMaxReliefCellIndex || | |
| 285 | 12632 | first + across * down > kMaxReliefCells) { | |
| 286 | ✗ | return std::nullopt; | |
| 287 | } | ||
| 288 | 12632 | index.bandFloats.insert(index.bandFloats.end(), {size, minI, minJ}); | |
| 289 | 12632 | index.bandCells.insert(index.bandCells.end(), {static_cast<std::uint32_t>(across), | |
| 290 | static_cast<std::uint32_t>(down), | ||
| 291 | static_cast<std::uint32_t>(first)}); | ||
| 292 | 12632 | index.cells.resize(2 * (first + across * down), 0); | |
| 293 | 12632 | cellOf.resize(members.size()); | |
| 294 | 1446667 | for (std::size_t m = 0; m < members.size(); ++m) { | |
| 295 | 1434035 | const CraterRecord& r = index.batch.records[members[m]]; | |
| 296 | 1434035 | const auto i = static_cast<std::size_t>(std::floor(r.x / size) - minI); | |
| 297 | 1434035 | const auto j = static_cast<std::size_t>(std::floor(r.y / size) - minJ); | |
| 298 | 1434035 | cellOf[m] = static_cast<std::uint32_t>(first + j * across + i); | |
| 299 | 1434035 | ++index.cells[2 * cellOf[m] + 1]; | |
| 300 | } | ||
| 301 | 12632 | auto start = static_cast<std::uint32_t>(index.entries.size()); | |
| 302 | 6125702 | for (std::size_t cell = first; cell < first + across * down; ++cell) { | |
| 303 | 6113070 | index.cells[2 * cell] = start; | |
| 304 | 6113070 | start += index.cells[2 * cell + 1]; | |
| 305 | } | ||
| 306 | 12632 | index.entries.resize(start); | |
| 307 | 12632 | filled.assign(across * down, 0); | |
| 308 | 1446667 | for (std::size_t m = 0; m < members.size(); ++m) { | |
| 309 | 1434035 | const std::uint32_t cell = cellOf[m]; | |
| 310 | 1434035 | index.entries[index.cells[2 * cell] + filled[cell - first]++] = members[m]; | |
| 311 | } | ||
| 312 | } | ||
| 313 | 2722 | index.kernel.shape = shape; | |
| 314 | 2722 | index.kernel.fitRadii = static_cast<float>(options.datumFitRadii); | |
| 315 | 2722 | index.kernel.latticeSide = static_cast<std::uint32_t>(kDatumLatticeSide); | |
| 316 | 2722 | index.kernel.bandCount = static_cast<std::uint32_t>(bands.size()); | |
| 317 | 2722 | return index; | |
| 318 | 2722 | } | |
| 319 | |||
| 320 | /// Fits the datums of craters `fitted` of `prepared`, all indexed in `index`: the ground under | ||
| 321 | /// each lattice point is read, the relief of the craters before each crater added by | ||
| 322 | /// `KernelDatumRelief` on GPU `backend` when it can take it and on CPU threads otherwise, and the | ||
| 323 | /// fit solved in double. | ||
| 324 | 17217 | void FitDatumsWithKernel(std::vector<PreparedCrater>& prepared, | |
| 325 | const std::vector<std::size_t>& fitted, const HeightSampler& ground, | ||
| 326 | const DatumKernelIndex& index, const CraterApplyOptions& options, | ||
| 327 | Backend backend, unsigned threads) { | ||
| 328 | 17217 | constexpr auto kPoints = static_cast<std::size_t>(kDatumLatticeSide * kDatumLatticeSide); | |
| 329 | 17217 | Raster<float> lattice(kPoints, fitted.size(), 0.0f); | |
| 330 | 17217 | std::vector<std::uint8_t> valid(lattice.Size(), 0); | |
| 331 | 76958 | const auto forEachPoint = [&](std::size_t n, const auto& visit) { | |
| 332 | 76958 | const double fitRadius = options.datumFitRadii * 0.5 * prepared[fitted[n]].crater.diameter; | |
| 333 | 615664 | for (int j = 0; j < kDatumLatticeSide; ++j) { | |
| 334 | 538706 | const double dy = DatumLatticeOffset(j, fitRadius); | |
| 335 | 4309648 | for (int i = 0; i < kDatumLatticeSide; ++i) { | |
| 336 | 3770942 | const double dx = DatumLatticeOffset(i, fitRadius); | |
| 337 | 3770942 | if (dx * dx + dy * dy < fitRadius * fitRadius) { | |
| 338 | 2847446 | visit(n * kPoints + static_cast<std::size_t>(j * kDatumLatticeSide + i), dx, | |
| 339 | dy); | ||
| 340 | } | ||
| 341 | } | ||
| 342 | } | ||
| 343 | 76958 | }; | |
| 344 | 17217 | ParallelFor(fitted.size(), threads, [&](std::size_t n) { | |
| 345 | 38479 | const Crater& c = prepared[fitted[n]].crater; | |
| 346 | 38479 | forEachPoint(n, [&](std::size_t at, double dx, double dy) { | |
| 347 | 1423723 | const float height = ground(c.x + dx, c.y + dy); | |
| 348 | 1423723 | if (!std::isnan(height)) { | |
| 349 | 1384142 | lattice.Values()[at] = height; | |
| 350 | 1384142 | valid[at] = 1; | |
| 351 | } | ||
| 352 | 1423723 | }); | |
| 353 | 38479 | }); | |
| 354 | 17217 | std::vector<std::uint32_t> targets(fitted.size()); | |
| 355 | 55696 | for (std::size_t n = 0; n < fitted.size(); ++n) { | |
| 356 | 38479 | targets[n] = index.slot[fitted[n]]; | |
| 357 | } | ||
| 358 | 17217 | const CraterKernelBatch& batch = index.batch; | |
| 359 | const std::array<KernelBuffer, 8> extra{ | ||
| 360 | 17217 | KernelBuffer{std::as_bytes(std::span(batch.records)), batch.records.size()}, | |
| 361 | 17217 | KernelBuffer{std::as_bytes(std::span(batch.data)), batch.data.size()}, | |
| 362 | 17217 | KernelBuffer{std::as_bytes(std::span(index.lows)), index.lows.size()}, | |
| 363 | 17217 | KernelBuffer{std::as_bytes(std::span(index.bandFloats)), index.bandFloats.size()}, | |
| 364 | 17217 | KernelBuffer{std::as_bytes(std::span(index.bandCells)), index.bandCells.size()}, | |
| 365 | 17217 | KernelBuffer{std::as_bytes(std::span(index.cells)), index.cells.size()}, | |
| 366 | 17217 | KernelBuffer{std::as_bytes(std::span(index.entries)), index.entries.size()}, | |
| 367 | 120519 | KernelBuffer{std::as_bytes(std::span(targets)), targets.size()}}; | |
| 368 | 17217 | Backend ran = backend; | |
| 369 | 17217 | if (!RunKernel(backend, CompiledKernel<CraterDatumKernel>(), | |
| 370 | std::as_bytes(std::span(&index.kernel, 1)), lattice, {0.0, 0.0}, extra)) { | ||
| 371 | 17213 | ran = Backend::Cpu; | |
| 372 | 17213 | ParallelFor(fitted.size(), threads, [&](std::size_t n) { | |
| 373 | 1827950 | for (std::size_t point = 0; point < kPoints; ++point) { | |
| 374 | 1791391 | const std::size_t at = n * kPoints + point; | |
| 375 | 1791391 | if (valid[at] != 0) { | |
| 376 | 2664906 | lattice.Values()[at] += KernelDatumRelief( | |
| 377 | 2664906 | index.kernel, batch.Data(), batch.records.data(), index.lows.data(), | |
| 378 | 1332453 | index.bandFloats.data(), index.bandCells.data(), index.cells.data(), | |
| 379 | 1332453 | index.entries.data(), targets[n], static_cast<std::uint32_t>(point)); | |
| 380 | } | ||
| 381 | } | ||
| 382 | 36559 | }); | |
| 383 | } | ||
| 384 | 17217 | ParallelFor(fitted.size(), threads, [&](std::size_t n) { | |
| 385 | 38479 | PlaneFit fit; | |
| 386 | 38479 | forEachPoint(n, [&](std::size_t at, double dx, double dy) { | |
| 387 | 1423723 | if (valid[at] != 0) { | |
| 388 | 1384142 | fit.Add(dx, dy, lattice.Values()[at]); | |
| 389 | } | ||
| 390 | 1423723 | }); | |
| 391 | 38479 | std::optional<Plane>& datum = prepared[fitted[n]].datum; | |
| 392 | 38479 | if (fit.Empty()) { | |
| 393 | 554 | datum.reset(); | |
| 394 | } else { | ||
| 395 | 37925 | datum = options.datum == CraterDatum::Plane ? fit.Solve() : Plane{fit.Mean(), 0.0, 0.0}; | |
| 396 | } | ||
| 397 | 38479 | }); | |
| 398 | 17217 | LogDebug("craters", "Fitted ", fitted.size(), " crater datums on ", BackendName(ran)); | |
| 399 | 17217 | } | |
| 400 | |||
| 401 | /// Gives every crater of `prepared` its shape (`SlopeShape`) and wear (`SlopeWearRate`) on the | ||
| 402 | /// gradient of `ground`, the ground beneath the layer, across its datum footprint. Done before any | ||
| 403 | /// datum is fitted, so the relief datums are fitted to is the relief applied, and each crater's | ||
| 404 | /// slope depends on no other crater of its layer. | ||
| 405 | 3393 | void ShapeOnSlopes(std::vector<PreparedCrater>& prepared, const HeightSampler& ground, | |
| 406 | const DegradationModel& degradation, const CraterApplyOptions& options) { | ||
| 407 | 3393 | const CraterSlopeResponse& response = options.slope; | |
| 408 | 3393 | if (!response.enabled && !response.criticalGradient) { | |
| 409 | 3386 | return; | |
| 410 | } | ||
| 411 | 1385 | for (PreparedCrater& crater : prepared) { | |
| 412 | 1378 | const Crater& c = crater.crater; | |
| 413 | 1378 | const double half = options.datumFitRadii * 0.5 * c.diameter; | |
| 414 | 5512 | const auto at = [&](double east, double north) { | |
| 415 | 5512 | const Vec2d p = options.surface.Locate({c.x, c.y}, {east, north}); | |
| 416 | 11024 | return static_cast<double>(ground(p.x, p.y)); | |
| 417 | 1378 | }; | |
| 418 | 1378 | const double east = (at(half, 0.0) - at(-half, 0.0)) / (2.0 * half); | |
| 419 | 1378 | const double north = (at(0.0, half) - at(0.0, -half)) / (2.0 * half); | |
| 420 | 1378 | if (!std::isfinite(east) || !std::isfinite(north)) { | |
| 421 | 2 | continue; | |
| 422 | } | ||
| 423 | 1376 | crater.slope = SlopeShape(c, east, north, response); | |
| 424 | 1376 | const double rate = SlopeWearRate(response, east, north); | |
| 425 | 1376 | if (rate > 1.0) { | |
| 426 | 597 | crater.state = degradation.WornState(c, rate); | |
| 427 | } | ||
| 428 | } | ||
| 429 | } | ||
| 430 | |||
| 431 | /// Fits the datums of craters `fitted` of `prepared` to `ground`, or to the raster as given when | ||
| 432 | /// it is empty; with the crater datum kernel when `shape` is set and it can take them. | ||
| 433 | 51 | void FitDatums(std::vector<PreparedCrater>& prepared, const std::vector<std::size_t>& fitted, | |
| 434 | const Raster<float>& raster, const CraterProfile& profile, | ||
| 435 | const std::optional<CraterShapeKernel>& shape, const CraterApplyOptions& options, | ||
| 436 | const CraterSampling& sampling) { | ||
| 437 | 51 | std::optional<Raster<float>> input; | |
| 438 | 51 | HeightSampler ground = sampling.ground; | |
| 439 | 51 | if (!ground) { | |
| 440 | 29 | input.emplace(raster); | |
| 441 | 29 | ground = [&input](double x, double y) { | |
| 442 | return static_cast<float>( | ||
| 443 | 149481 | Sample(*input, x, y).value_or(std::numeric_limits<double>::quiet_NaN())); | |
| 444 | 29 | }; | |
| 445 | } | ||
| 446 | 51 | if (shape && options.datum != CraterDatum::CentreHeight) { | |
| 447 | 47 | if (const std::optional<DatumKernelIndex> index = | |
| 448 | 47 | IndexDatumRelief(prepared, fitted, profile, *shape, options)) { | |
| 449 | 47 | FitDatumsWithKernel(prepared, fitted, ground, *index, options, sampling.backend, | |
| 450 | 47 | sampling.threads); | |
| 451 | 47 | return; | |
| 452 | 47 | } | |
| 453 | } | ||
| 454 | 4 | const CraterReliefIndex index(prepared, profile, MaxOutlineScale(options.irregularity), | |
| 455 | 4 | options.surface); | |
| 456 | 4 | ParallelFor(fitted.size(), sampling.threads, [&](std::size_t n) { | |
| 457 | 108 | const std::size_t k = fitted[n]; | |
| 458 | 108 | prepared[k].datum = | |
| 459 | 216 | FitCraterDatum(prepared[k].crater, options, DatumGround(ground, index, k)); | |
| 460 | 108 | }); | |
| 461 | 98 | } | |
| 462 | |||
| 463 | /// A prepared crater with the box its profile can reach. | ||
| 464 | struct ReachingCrater { | ||
| 465 | PreparedCrater* prepared; | ||
| 466 | Bounds2d footprint; | ||
| 467 | double reach; | ||
| 468 | |||
| 469 | 1367427277 | bool Reaches(const Vec2d& p, const Surface& surface) const { | |
| 470 | 1368191852 | return surface.Contains(footprint, p) && | |
| 471 | 764575 | CraterProfileDistance(prepared->crater, prepared->outline, p.x, p.y, surface) < | |
| 472 | 1368191852 | reach; | |
| 473 | } | ||
| 474 | }; | ||
| 475 | |||
| 476 | /// The craters of `prepared` reaching `reachRadii` radii that reach any of `samples`, in order. | ||
| 477 | 129391 | std::vector<ReachingCrater> CratersReaching(const SampleBlock& samples, | |
| 478 | std::vector<PreparedCrater>& prepared, | ||
| 479 | double reachRadii, const Surface& surface) { | ||
| 480 | 129391 | std::vector<ReachingCrater> reaching; | |
| 481 | 5728245 | for (PreparedCrater& crater : prepared) { | |
| 482 | 5598854 | const ReachingCrater candidate{&crater, Footprint(crater, reachRadii, surface), | |
| 483 | 5598854 | reachRadii * 0.5 * crater.crater.diameter}; | |
| 484 | 1372585764 | for (std::size_t i = 0; i < samples.Size(); ++i) { | |
| 485 | 1367250099 | if (candidate.Reaches(samples.Position(i), surface)) { | |
| 486 | 263189 | reaching.push_back(candidate); | |
| 487 | 263189 | break; | |
| 488 | } | ||
| 489 | } | ||
| 490 | } | ||
| 491 | 129391 | return reaching; | |
| 492 | } | ||
| 493 | |||
| 494 | /// Side of the tiles craters are binned into for the kernel, in cells. | ||
| 495 | constexpr std::size_t kTileCells = 32; | ||
| 496 | |||
| 497 | /// Shapes the craters `onRaster` of `prepared` into `raster` with the kernel: each crater's | ||
| 498 | /// `windows` are cut into the tiles they cover, the tiles run on `threads` threads, and a tile | ||
| 499 | /// shapes its craters in order over their own cells, so every cell sees them in order. | ||
| 500 | 366 | void ShapeTiles(Raster<float>& raster, unsigned threads, const CraterShapeKernel& kernel, | |
| 501 | const CraterKernelBatch& batch, const std::vector<PreparedCrater>& prepared, | ||
| 502 | const std::vector<std::size_t>& onRaster, | ||
| 503 | const std::vector<std::vector<PixelWindow>>& windows) { | ||
| 504 | 366 | const std::size_t tilesX = (raster.Width() + kTileCells - 1) / kTileCells; | |
| 505 | 366 | const std::size_t tilesY = (raster.Height() + kTileCells - 1) / kTileCells; | |
| 506 | struct Entry { | ||
| 507 | std::uint32_t crater; | ||
| 508 | PixelWindow cells; | ||
| 509 | }; | ||
| 510 | 366 | std::vector<std::vector<Entry>> tiles(tilesX * tilesY); | |
| 511 | 108945 | for (std::size_t n = 0; n < onRaster.size(); ++n) { | |
| 512 | 217158 | for (const PixelWindow& window : windows[n]) { | |
| 513 | 271137 | for (std::size_t ty = window.rowBegin / kTileCells; | |
| 514 | 271137 | ty <= (window.rowEnd - 1) / kTileCells; ++ty) { | |
| 515 | 162558 | for (std::size_t tx = window.colBegin / kTileCells; | |
| 516 | 441221 | tx <= (window.colEnd - 1) / kTileCells; ++tx) { | |
| 517 | 278663 | PixelWindow cells; | |
| 518 | 278663 | cells.colBegin = std::max(window.colBegin, tx * kTileCells); | |
| 519 | 278663 | cells.colEnd = std::min(window.colEnd, (tx + 1) * kTileCells); | |
| 520 | 278663 | cells.rowBegin = std::max(window.rowBegin, ty * kTileCells); | |
| 521 | 278663 | cells.rowEnd = std::min(window.rowEnd, (ty + 1) * kTileCells); | |
| 522 | 278663 | tiles[ty * tilesX + tx].push_back({static_cast<std::uint32_t>(n), cells}); | |
| 523 | } | ||
| 524 | } | ||
| 525 | } | ||
| 526 | } | ||
| 527 | 366 | const CraterData data = batch.Data(); | |
| 528 | 366 | ParallelFor(tiles.size(), threads, [&](std::size_t tile) { | |
| 529 | 326670 | for (const Entry& entry : tiles[tile]) { | |
| 530 | 278663 | const Crater& c = prepared[onRaster[entry.crater]].crater; | |
| 531 | 278663 | const CraterRecord& record = batch.records[entry.crater]; | |
| 532 | 4203735 | for (std::size_t row = entry.cells.rowBegin; row < entry.cells.rowEnd; ++row) { | |
| 533 | 73424549 | for (std::size_t col = entry.cells.colBegin; col < entry.cells.colEnd; ++col) { | |
| 534 | 69499477 | float& cell = raster(col, row); | |
| 535 | 69499477 | if (!raster.IsNoData(cell)) { | |
| 536 | 69499476 | const Vec2d centre = raster.CellCenter(col, row); | |
| 537 | 69499476 | cell = ApplyCraterKernel(kernel, data, record, cell, | |
| 538 | 69499476 | static_cast<float>(centre.x - c.x), | |
| 539 | 69499476 | static_cast<float>(centre.y - c.y)); | |
| 540 | } | ||
| 541 | } | ||
| 542 | } | ||
| 543 | } | ||
| 544 | 48007 | }); | |
| 545 | 366 | } | |
| 546 | |||
| 547 | /// Side of the tiles craters are binned into on a GPU, in cells. | ||
| 548 | constexpr std::size_t kGpuTileCells = 16; | ||
| 549 | |||
| 550 | /// Shapes the craters `onRaster` of `prepared` into `raster` on GPU `backend`, each cell running | ||
| 551 | /// the craters of its tile in order; false, leaving `raster` unchanged, when the backend is the | ||
| 552 | /// CPU or cannot take the grid (`RunKernel`). | ||
| 553 | 378 | bool ShapeTilesOnGpu(Raster<float>& raster, Backend backend, const CraterShapeKernel& kernel, | |
| 554 | CraterKernelBatch& batch, const std::vector<PreparedCrater>& prepared, | ||
| 555 | const std::vector<std::size_t>& onRaster, | ||
| 556 | const std::vector<std::vector<PixelWindow>>& windows) { | ||
| 557 | 378 | if (backend == Backend::Cpu || batch.records.empty()) { | |
| 558 | 366 | return false; | |
| 559 | } | ||
| 560 | 12 | const std::size_t tilesX = (raster.Width() + kGpuTileCells - 1) / kGpuTileCells; | |
| 561 | 12 | const std::size_t tilesY = (raster.Height() + kGpuTileCells - 1) / kGpuTileCells; | |
| 562 | 12 | std::vector<std::vector<std::uint32_t>> lists(tilesX * tilesY); | |
| 563 | 2148 | for (std::size_t n = 0; n < onRaster.size(); ++n) { | |
| 564 | 4272 | for (const PixelWindow& window : windows[n]) { | |
| 565 | 12451 | for (std::size_t ty = window.rowBegin / kGpuTileCells; | |
| 566 | 12451 | ty <= (window.rowEnd - 1) / kGpuTileCells; ++ty) { | |
| 567 | 73485 | for (std::size_t tx = window.colBegin / kGpuTileCells; | |
| 568 | 73485 | tx <= (window.colEnd - 1) / kGpuTileCells; ++tx) { | |
| 569 | 63170 | std::vector<std::uint32_t>& list = lists[ty * tilesX + tx]; | |
| 570 | 63170 | if (list.empty() || list.back() != n) { | |
| 571 | 63170 | list.push_back(static_cast<std::uint32_t>(n)); | |
| 572 | } | ||
| 573 | } | ||
| 574 | } | ||
| 575 | } | ||
| 576 | } | ||
| 577 | std::vector<std::uint32_t> tiles{static_cast<std::uint32_t>(tilesX), | ||
| 578 | 24 | static_cast<std::uint32_t>(kGpuTileCells)}; | |
| 579 | 12 | std::vector<std::uint32_t> entries; | |
| 580 | 1353 | for (const std::vector<std::uint32_t>& list : lists) { | |
| 581 | 1341 | tiles.push_back(static_cast<std::uint32_t>(entries.size())); | |
| 582 | 1341 | tiles.push_back(static_cast<std::uint32_t>(list.size())); | |
| 583 | 1341 | entries.insert(entries.end(), list.begin(), list.end()); | |
| 584 | } | ||
| 585 | 12 | if (entries.empty()) { | |
| 586 | ✗ | entries.push_back(0); | |
| 587 | } | ||
| 588 | 12 | const Bounds2d bounds = RasterGrid::FromRaster(raster).Bounds(); | |
| 589 | 12 | const Vec2d origin{0.5 * (bounds.minX + bounds.maxX), 0.5 * (bounds.minY + bounds.maxY)}; | |
| 590 | 2148 | for (std::size_t n = 0; n < onRaster.size(); ++n) { | |
| 591 | 2136 | const Crater& c = prepared[onRaster[n]].crater; | |
| 592 | 2136 | batch.records[n].x = static_cast<float>(c.x - origin.x); | |
| 593 | 2136 | batch.records[n].y = static_cast<float>(c.y - origin.y); | |
| 594 | } | ||
| 595 | 12 | if (batch.data.empty()) { | |
| 596 | 9 | batch.data.push_back(0.0f); | |
| 597 | } | ||
| 598 | const std::array<KernelBuffer, 4> extra{ | ||
| 599 | 12 | KernelBuffer{std::as_bytes(std::span(batch.records)), batch.records.size()}, | |
| 600 | 12 | KernelBuffer{std::as_bytes(std::span(batch.data)), batch.data.size()}, | |
| 601 | 12 | KernelBuffer{std::as_bytes(std::span(tiles)), tiles.size()}, | |
| 602 | 36 | KernelBuffer{std::as_bytes(std::span(entries)), entries.size()}}; | |
| 603 | 12 | if (!RunKernel(backend, CompiledKernel<CraterShapeKernel>(), | |
| 604 | std::as_bytes(std::span(&kernel, 1)), raster, origin, extra)) { | ||
| 605 | ✗ | return false; | |
| 606 | } | ||
| 607 | 12 | LogDebug("craters", "Shaped ", onRaster.size(), " craters on ", BackendName(backend)); | |
| 608 | 12 | return true; | |
| 609 | 12 | } | |
| 610 | |||
| 611 | } // namespace | ||
| 612 | |||
| 613 | 1773812 | double CraterReachRadii(const CraterProfile& profile, const CraterApplyOptions& options) { | |
| 614 | 1773812 | const double reach = profile.ReachRadii(); | |
| 615 | 1773812 | return options.texture.ReachesBeyondRim() ? std::fmax(reach, options.texture.ejectaReachRadii) | |
| 616 | 1773812 | : reach; | |
| 617 | } | ||
| 618 | |||
| 619 | 1369491 | float ApplyCraterTo(float height, double x, double y, const PreparedCrater& crater, | |
| 620 | const CraterProfile& profile, const CraterApplyOptions& options) { | ||
| 621 | 1369491 | if (const std::optional<CraterShapeKernel> kernel = MakeCraterShapeKernel(profile, options)) { | |
| 622 | const CraterKernelBatch batch = | ||
| 623 | 1585237 | MakeCraterKernelBatch(*kernel, profile, {crater}, {0}, options); | |
| 624 | 317046 | return ApplyCraterKernel(*kernel, batch.Data(), batch.records[0], height, | |
| 625 | 317046 | static_cast<float>(x - crater.crater.x), | |
| 626 | 634092 | static_cast<float>(y - crater.crater.y)); | |
| 627 | 317046 | } | |
| 628 | 1052444 | const double diameter = crater.crater.diameter; | |
| 629 | 1052444 | const double radius = 0.5 * diameter; | |
| 630 | 1052444 | const Vec2d offset = options.surface.Offset({crater.crater.x, crater.crater.y}, {x, y}); | |
| 631 | 1052444 | const CraterSlopeShape& slope = crater.slope; | |
| 632 | const double distance = | ||
| 633 | 1052444 | CraterProfileDistance(crater.crater, crater.outline, | |
| 634 | 1052444 | slope.Level() ? offset : SlopeShiftedOffset(slope, offset, radius)); | |
| 635 | 1052444 | if (distance >= CraterReachRadii(profile, options) * radius) { | |
| 636 | 262401 | return height; | |
| 637 | } | ||
| 638 | 790043 | const double radii = distance / radius; | |
| 639 | 790043 | const double current = height; | |
| 640 | 790043 | double replaced = 0.0; | |
| 641 | 790043 | if (crater.datum) { | |
| 642 | 236217 | double weight = options.interiorReplacement * options.blend.Weight(radii); | |
| 643 | 236217 | if (options.texture.mantling > 0.0 && radii > 1.0) { | |
| 644 | 27000 | weight = std::fmax(weight, | |
| 645 | 27000 | options.texture.mantling * | |
| 646 | 27000 | EjectaThicknessShare(radii, options.texture.ejectaReachRadii)); | |
| 647 | } | ||
| 648 | 236217 | replaced = weight * (crater.datum->At(offset.x, offset.y) - current); | |
| 649 | } | ||
| 650 | 790043 | double bare = ProfileHeight(profile, distance, diameter, crater.state, options.rimSmoothing); | |
| 651 | 790043 | if (bare > 0.0 && !slope.Level()) { | |
| 652 | 15708 | bare *= SlopeRaisedShare(slope, offset); | |
| 653 | } | ||
| 654 | const double shaped = | ||
| 655 | 790043 | bare + NoiseHeight(crater.crater, offset, radii, profile.ReachRadii(), options.noise) + | |
| 656 | 790043 | TextureHeight(crater.crater, crater.state, offset, radii, options.texture); | |
| 657 | 790043 | return static_cast<float>(current + crater.resolution * shaped + crater.resolution * replaced); | |
| 658 | 317048 | } | |
| 659 | |||
| 660 | 512 | void ApplyCraters(Raster<float>& raster, const CraterList& craters, const CraterShapeModels& models, | |
| 661 | const CraterApplyOptions& options, const CraterSampling& sampling) { | ||
| 662 | 512 | ValidateCraterApplyOptions(options); | |
| 663 | 504 | RequireNonNegative(sampling.spacing, "Sample spacing"); | |
| 664 | 503 | if (raster.Empty() || craters.Empty()) { | |
| 665 | 435 | return; | |
| 666 | } | ||
| 667 | 446 | const Stopwatch stopwatch; | |
| 668 | 446 | const double resolved = ResolvedSpacing(raster, sampling.spacing, options.surface); | |
| 669 | 445 | const CraterProfile& profile = *models.profile; | |
| 670 | std::vector<PreparedCrater> prepared = | ||
| 671 | 445 | PrepareCraters(craters, options, resolved, *models.degradation, sampling.reshape); | |
| 672 | 445 | std::vector<std::size_t> onRaster; | |
| 673 | 445 | std::vector<std::vector<PixelWindow>> windows; | |
| 674 | 144400 | for (std::size_t k = 0; k < prepared.size(); ++k) { | |
| 675 | std::vector<PixelWindow> around = PixelWindowsIn( | ||
| 676 | 143955 | raster, Footprint(prepared[k], CraterReachRadii(profile, options), options.surface), | |
| 677 | 287910 | options.surface); | |
| 678 | 143955 | if (!around.empty()) { | |
| 679 | 112435 | onRaster.push_back(k); | |
| 680 | 112435 | windows.push_back(std::move(around)); | |
| 681 | } | ||
| 682 | 143955 | } | |
| 683 | 445 | const std::optional<CraterShapeKernel> kernel = MakeCraterShapeKernel(profile, options); | |
| 684 | 445 | if (options.datum != CraterDatum::None && !onRaster.empty()) { | |
| 685 | const HeightSampler beneath = | ||
| 686 | 51 | sampling.ground ? sampling.ground : HeightSampler([&raster](double x, double y) { | |
| 687 | return static_cast<float>( | ||
| 688 | 2412 | Sample(raster, x, y).value_or(std::numeric_limits<double>::quiet_NaN())); | |
| 689 | 51 | }); | |
| 690 | 51 | ShapeOnSlopes(prepared, beneath, *models.degradation, options); | |
| 691 | 51 | FitDatums(prepared, onRaster, raster, profile, kernel, options, sampling); | |
| 692 | 51 | } | |
| 693 | 445 | if (kernel) { | |
| 694 | CraterKernelBatch batch = | ||
| 695 | 378 | MakeCraterKernelBatch(*kernel, profile, prepared, onRaster, options); | |
| 696 | 378 | if (!ShapeTilesOnGpu(raster, sampling.backend, *kernel, batch, prepared, onRaster, | |
| 697 | windows)) { | ||
| 698 | 366 | ShapeTiles(raster, sampling.threads, *kernel, batch, prepared, onRaster, windows); | |
| 699 | } | ||
| 700 | 378 | LogInfo("craters", "Applied ", onRaster.size(), " of ", craters.Size(), " craters (", | |
| 701 | 378 | craters.Size() - prepared.size(), " unresolved at ", resolved, " m, ", | |
| 702 | 378 | prepared.size() - onRaster.size(), " off the raster) in ", | |
| 703 | 756 | FormatDuration(stopwatch.Milliseconds())); | |
| 704 | 378 | return; | |
| 705 | 378 | } | |
| 706 | 134 | ForEachFeatureCell( | |
| 707 | 67 | raster, sampling.threads, onRaster.size(), | |
| 708 | 10344 | [&](std::size_t n) -> const std::vector<PixelWindow>& { return windows[n]; }, | |
| 709 | 67 | [&](float& cell, const FeatureCell& at) { | |
| 710 | 734295 | cell = ApplyCraterTo(cell, at.centre.x, at.centre.y, prepared[onRaster[at.feature]], | |
| 711 | profile, options); | ||
| 712 | 734295 | }); | |
| 713 | 67 | LogInfo("craters", "Applied ", onRaster.size(), " of ", craters.Size(), " craters (", | |
| 714 | 67 | craters.Size() - prepared.size(), " unresolved at ", resolved, " m, ", | |
| 715 | 67 | prepared.size() - onRaster.size(), " off the raster) in ", | |
| 716 | 134 | FormatDuration(stopwatch.Milliseconds())); | |
| 717 | 1201 | } | |
| 718 | |||
| 719 | 245651 | void ApplyCraters(SampleBlock& samples, const CraterList& craters, const CraterShapeModels& models, | |
| 720 | const CraterApplyOptions& options, const CraterSampling& sampling) { | ||
| 721 | 245651 | if (Raster<float>* grid = samples.Grid()) { | |
| 722 | 438 | ApplyCraters(*grid, craters, models, options, sampling); | |
| 723 | 245079 | return; | |
| 724 | } | ||
| 725 | 245213 | ValidateCraterApplyOptions(options); | |
| 726 | 245213 | RequirePositive(sampling.spacing, "Sample spacing"); | |
| 727 | 245213 | if (options.datum != CraterDatum::None && !sampling.ground) { | |
| 728 | ✗ | throw std::invalid_argument("Craters with a datum at points need the ground"); | |
| 729 | } | ||
| 730 | 245213 | if (samples.Size() == 0 || craters.Empty()) { | |
| 731 | 115822 | return; | |
| 732 | } | ||
| 733 | 129391 | const CraterProfile& profile = *models.profile; | |
| 734 | std::vector<PreparedCrater> prepared = | ||
| 735 | 129391 | PrepareCraters(craters, options, sampling.spacing, *models.degradation, sampling.reshape); | |
| 736 | const std::vector<ReachingCrater> reaching = | ||
| 737 | 129391 | CratersReaching(samples, prepared, CraterReachRadii(profile, options), options.surface); | |
| 738 | 129391 | const std::optional<CraterShapeKernel> kernel = MakeCraterShapeKernel(profile, options); | |
| 739 | 129391 | std::vector<std::size_t> indices; | |
| 740 | 392580 | for (const ReachingCrater& crater : reaching) { | |
| 741 | 263189 | indices.push_back(static_cast<std::size_t>(crater.prepared - prepared.data())); | |
| 742 | } | ||
| 743 | 129391 | if (options.datum != CraterDatum::None && !reaching.empty()) { | |
| 744 | 3342 | ShapeOnSlopes(prepared, sampling.ground, *models.degradation, options); | |
| 745 | 3342 | std::optional<DatumKernelIndex> kernelIndex; | |
| 746 | 3342 | bool indexed = false; | |
| 747 | 3342 | std::optional<CraterReliefIndex> index; | |
| 748 | 17241 | const auto fit = [&](std::size_t k) -> std::optional<Plane> { | |
| 749 | 17241 | if (kernel && options.datum != CraterDatum::CentreHeight) { | |
| 750 | 17170 | if (!indexed) { | |
| 751 | 2675 | kernelIndex = IndexDatumRelief(prepared, indices, profile, *kernel, options); | |
| 752 | 2675 | indexed = true; | |
| 753 | } | ||
| 754 | 17170 | if (kernelIndex) { | |
| 755 | 34340 | FitDatumsWithKernel(prepared, {k}, sampling.ground, *kernelIndex, options, | |
| 756 | Backend::Cpu, 1); | ||
| 757 | 17170 | return prepared[k].datum; | |
| 758 | } | ||
| 759 | } | ||
| 760 | 71 | if (!index) { | |
| 761 | 62 | index.emplace(prepared, profile, MaxOutlineScale(options.irregularity), | |
| 762 | 62 | options.surface); | |
| 763 | } | ||
| 764 | 71 | return FitCraterDatum(prepared[k].crater, options, | |
| 765 | 142 | DatumGround(sampling.ground, *index, k)); | |
| 766 | 3342 | }; | |
| 767 | 136030 | for (const std::size_t k : indices) { | |
| 768 | 132688 | if (sampling.cache != nullptr) { | |
| 769 | 132226 | prepared[k].datum = *sampling.cache->Get<std::optional<Plane>>( | |
| 770 | 132226 | sampling.owner, prepared[k].crater.id, sampling.spacing, | |
| 771 | 413457 | [&] { return fit(k); }); | |
| 772 | } else { | ||
| 773 | 462 | prepared[k].datum = fit(k); | |
| 774 | } | ||
| 775 | } | ||
| 776 | 3342 | } | |
| 777 | 129391 | if (kernel) { | |
| 778 | const CraterKernelBatch batch = | ||
| 779 | 128819 | MakeCraterKernelBatch(*kernel, profile, prepared, indices, options); | |
| 780 | 128819 | const CraterData data = batch.Data(); | |
| 781 | 1036034 | for (std::size_t i = 0; i < samples.Size(); ++i) { | |
| 782 | 907215 | if (samples.IsNoData(i)) { | |
| 783 | ✗ | continue; | |
| 784 | } | ||
| 785 | 907215 | const Vec2d p = samples.Position(i); | |
| 786 | 907215 | float height = samples.Height(i); | |
| 787 | 35131287 | for (std::size_t j = 0; j < reaching.size(); ++j) { | |
| 788 | 34224072 | if (options.surface.Contains(reaching[j].footprint, p)) { | |
| 789 | 3191454 | const Crater& c = reaching[j].prepared->crater; | |
| 790 | 3191454 | height = ApplyCraterKernel(*kernel, data, batch.records[j], height, | |
| 791 | 3191454 | static_cast<float>(p.x - c.x), | |
| 792 | 3191454 | static_cast<float>(p.y - c.y)); | |
| 793 | } | ||
| 794 | } | ||
| 795 | 907215 | samples.Height(i) = height; | |
| 796 | } | ||
| 797 | 128819 | return; | |
| 798 | 128819 | } | |
| 799 | 141409 | for (std::size_t i = 0; i < samples.Size(); ++i) { | |
| 800 | 140837 | const Vec2d p = samples.Position(i); | |
| 801 | 318015 | for (const ReachingCrater& crater : reaching) { | |
| 802 | 177178 | if (!samples.IsNoData(i) && crater.Reaches(p, options.surface)) { | |
| 803 | 11749 | samples.Height(i) = | |
| 804 | 11749 | ApplyCraterTo(samples.Height(i), p.x, p.y, *crater.prepared, profile, options); | |
| 805 | } | ||
| 806 | } | ||
| 807 | } | ||
| 808 | 387029 | } | |
| 809 | |||
| 810 | } // namespace terraingen | ||
| 811 |