sources/field.cpp
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | #include "terraingen/sources/field.h" | ||
| 2 | |||
| 3 | #include "terraingen/raster/geotiff/geotiff.h" | ||
| 4 | #include "terraingen/util/checks.h" | ||
| 5 | |||
| 6 | #include <algorithm> | ||
| 7 | #include <cmath> | ||
| 8 | #include <stdexcept> | ||
| 9 | #include <tuple> | ||
| 10 | #include <utility> | ||
| 11 | |||
| 12 | namespace terraingen { | ||
| 13 | |||
| 14 | 15 | void ValidateRasterFieldOptions(const RasterFieldOptions& options) { | |
| 15 | 15 | if (options.method == ResampleMethod::Cubic) { | |
| 16 | throw std::invalid_argument( | ||
| 17 | 2 | "A raster field reads nearest or bilinear, which stay within the cells read"); | |
| 18 | } | ||
| 19 | 13 | RequirePositive(options.scale, "Raster field scale"); | |
| 20 | 12 | } | |
| 21 | |||
| 22 | 15 | RasterField::RasterField(Raster<float> raster, RasterFieldOptions options) | |
| 23 | 15 | : raster_(std::move(raster)), options_(options) { | |
| 24 | 15 | ValidateRasterFieldOptions(options_); | |
| 25 | 12 | if (raster_.Empty()) { | |
| 26 | 1 | throw std::invalid_argument("A raster field needs a raster"); | |
| 27 | } | ||
| 28 | 11 | const GeoTransform& transform = raster_.Transform(); | |
| 29 | 11 | if (transform.rowRotation != 0.0 || transform.colRotation != 0.0) { | |
| 30 | 1 | throw std::invalid_argument("A raster field cannot read a rotated raster"); | |
| 31 | } | ||
| 32 | 78 | for (const float value : raster_.Values()) { | |
| 33 | 125 | if (std::isfinite(value) && !raster_.IsNoData(value) && | |
| 34 | 57 | !(static_cast<double>(value) * options_.scale <= max_)) { | |
| 35 | 29 | max_ = static_cast<double>(value) * options_.scale; | |
| 36 | } | ||
| 37 | } | ||
| 38 | 10 | if (std::isfinite(options_.fallback) && !(options_.fallback <= max_)) { | |
| 39 | 1 | max_ = options_.fallback; | |
| 40 | } | ||
| 41 | 20 | } | |
| 42 | |||
| 43 | 40022 | double RasterField::Sample(double x, double y) const { | |
| 44 | 40022 | const Vec2d pixel = raster_.Transform().WorldToPixel(x, y); | |
| 45 | 40022 | const auto width = static_cast<double>(raster_.Width()); | |
| 46 | 40022 | const auto height = static_cast<double>(raster_.Height()); | |
| 47 | 40022 | if (!(pixel.x >= 0.0 && pixel.y >= 0.0 && pixel.x < width && pixel.y < height)) { | |
| 48 | 30008 | return std::numeric_limits<double>::quiet_NaN(); | |
| 49 | } | ||
| 50 | 17524 | const auto valid = [&](std::size_t col, std::size_t row) { | |
| 51 | 17524 | const float value = raster_(col, row); | |
| 52 | 17524 | return std::isfinite(value) && !raster_.IsNoData(value); | |
| 53 | 10014 | }; | |
| 54 | 10014 | if (options_.method == ResampleMethod::Nearest) { | |
| 55 | 3 | const auto col = static_cast<std::size_t>(pixel.x); | |
| 56 | 3 | const auto row = static_cast<std::size_t>(pixel.y); | |
| 57 | 3 | return valid(col, row) ? static_cast<double>(raster_(col, row)) | |
| 58 | 3 | : std::numeric_limits<double>::quiet_NaN(); | |
| 59 | } | ||
| 60 | 10011 | const double u = std::clamp(pixel.x - 0.5, 0.0, width - 1.0); | |
| 61 | 10011 | const double v = std::clamp(pixel.y - 0.5, 0.0, height - 1.0); | |
| 62 | 10011 | const auto col = static_cast<std::size_t>(u); | |
| 63 | 10011 | const auto row = static_cast<std::size_t>(v); | |
| 64 | 10011 | const double tx = u - static_cast<double>(col); | |
| 65 | 10011 | const double ty = v - static_cast<double>(row); | |
| 66 | 10011 | const std::size_t right = std::min(col + 1, raster_.Width() - 1); | |
| 67 | 10011 | const std::size_t below = std::min(row + 1, raster_.Height() - 1); | |
| 68 | 10011 | double sum = 0.0; | |
| 69 | 10011 | double weights = 0.0; | |
| 70 | 90099 | for (const auto& [c, r, w] : | |
| 71 | 30033 | {std::tuple{col, row, (1.0 - tx) * (1.0 - ty)}, std::tuple{right, row, tx * (1.0 - ty)}, | |
| 72 | 60066 | std::tuple{col, below, (1.0 - tx) * ty}, std::tuple{right, below, tx * ty}}) { | |
| 73 | 40044 | if (w > 0.0 && valid(c, r)) { | |
| 74 | 17516 | sum += w * static_cast<double>(raster_(c, r)); | |
| 75 | 17516 | weights += w; | |
| 76 | } | ||
| 77 | } | ||
| 78 | 10011 | return weights > 0.0 ? sum / weights : std::numeric_limits<double>::quiet_NaN(); | |
| 79 | } | ||
| 80 | |||
| 81 | 40019 | double RasterField::ValueAt(double x, double y) const { | |
| 82 | 40019 | double value = Sample(x, y); | |
| 83 | 40019 | if (std::isnan(value) && options_.periodicLongitude) { | |
| 84 | 2 | value = Sample(x + 360.0, y); | |
| 85 | 2 | if (std::isnan(value)) { | |
| 86 | 1 | value = Sample(x - 360.0, y); | |
| 87 | } | ||
| 88 | } | ||
| 89 | 40019 | if (std::isnan(value)) { | |
| 90 | 30007 | return options_.fallback; | |
| 91 | } | ||
| 92 | 10012 | return value * options_.scale; | |
| 93 | } | ||
| 94 | |||
| 95 | 4 | std::shared_ptr<const RasterField> ReadRasterField(const std::filesystem::path& path, | |
| 96 | RasterFieldOptions options, std::size_t band) { | ||
| 97 | 5 | return std::make_shared<RasterField>(ReadGeoTiff(path, band).raster, options); | |
| 98 | } | ||
| 99 | |||
| 100 | } // namespace terraingen | ||
| 101 |