io/raster.cpp
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // SPDX-License-Identifier: BSD-3-Clause | ||
| 2 | /// @file raster.cpp | ||
| 3 | /// @brief Raster reading: decoded images, mosaics, PDS3 labels and images, TIFF strips and tiles, | ||
| 4 | /// and planetary DEMs. | ||
| 5 | #include "ernie/io/raster.h" | ||
| 6 | |||
| 7 | #include "ernie/io/image.h" | ||
| 8 | |||
| 9 | #include "ernie/util/bytes.h" | ||
| 10 | #include "ernie/util/file.h" | ||
| 11 | #include "ernie/util/profiler.h" | ||
| 12 | #include "ernie/util/text.h" | ||
| 13 | |||
| 14 | #include <stb_image.h> | ||
| 15 | |||
| 16 | #include <algorithm> | ||
| 17 | #include <bit> | ||
| 18 | #include <climits> | ||
| 19 | #include <cstring> | ||
| 20 | #include <fstream> | ||
| 21 | #include <limits> | ||
| 22 | #include <numbers> | ||
| 23 | #include <sstream> | ||
| 24 | #include <stdexcept> | ||
| 25 | #include <utility> | ||
| 26 | |||
| 27 | namespace ernie { | ||
| 28 | |||
| 29 | namespace fs = std::filesystem; | ||
| 30 | |||
| 31 | namespace { | ||
| 32 | |||
| 33 | // ---------------------------------------------------------------------------------------------- | ||
| 34 | // Decoded images | ||
| 35 | |||
| 36 | /// An 8- or 16-bit PNG or JPEG, its header read at once and its pixels decoded whole into memory | ||
| 37 | /// at the first read. | ||
| 38 | class DecodedRaster : public GridSource { | ||
| 39 | public: | ||
| 40 | 34 | explicit DecodedRaster(const fs::path& path) : name_(path.string()) { | |
| 41 | 34 | int width = 0; | |
| 42 | 34 | int height = 0; | |
| 43 | 34 | int channels = 0; | |
| 44 | 34 | if (stbi_info(name_.c_str(), &width, &height, &channels) == 0) { | |
| 45 | ✗ | throw std::runtime_error("Cannot read " + name_ + ": " + stbi_failure_reason()); | |
| 46 | } | ||
| 47 | 34 | if (stbi_is_16_bit(name_.c_str()) != 0) { | |
| 48 | 1 | type_ = GridSampleType::UInt16; | |
| 49 | } | ||
| 50 | 34 | width_ = static_cast<std::uint32_t>(width); | |
| 51 | 34 | height_ = static_cast<std::uint32_t>(height); | |
| 52 | 34 | channels_ = static_cast<std::uint32_t>(channels); | |
| 53 | 34 | } | |
| 54 | |||
| 55 | 149 | std::uint32_t Width() const override { | |
| 56 | 149 | return width_; | |
| 57 | } | ||
| 58 | 145 | std::uint32_t Height() const override { | |
| 59 | 145 | return height_; | |
| 60 | } | ||
| 61 | 160 | std::uint32_t Channels() const override { | |
| 62 | 160 | return channels_; | |
| 63 | } | ||
| 64 | 208 | GridSampleType Type() const override { | |
| 65 | 208 | return type_; | |
| 66 | } | ||
| 67 | |||
| 68 | 37 | void Read(std::uint32_t x, std::uint32_t y, std::uint32_t w, std::uint32_t h, | |
| 69 | std::span<std::byte> out) const override { | ||
| 70 | 37 | CheckRead(x, y, w, h, out); | |
| 71 | 37 | if (pixels_ == nullptr) { | |
| 72 | 28 | Decode(); | |
| 73 | } | ||
| 74 | 37 | const std::size_t pixel = std::size_t{channels_} * SampleBytes(); | |
| 75 | 37 | const auto* bytes = static_cast<const std::byte*>(pixels_.get()); | |
| 76 | 11191 | for (std::uint32_t row = 0; row < h; ++row) { | |
| 77 | 11154 | std::memcpy(out.data() + std::size_t{row} * w * pixel, | |
| 78 | 11154 | bytes + (std::size_t{y + row} * width_ + x) * pixel, | |
| 79 | 11154 | std::size_t{w} * pixel); | |
| 80 | } | ||
| 81 | 37 | } | |
| 82 | |||
| 83 | 38 | void Release() const override { | |
| 84 | 38 | pixels_.reset(); | |
| 85 | 38 | } | |
| 86 | |||
| 87 | private: | ||
| 88 | 28 | void Decode() const { | |
| 89 | 28 | int width = 0; | |
| 90 | 28 | int height = 0; | |
| 91 | 28 | int channels = 0; | |
| 92 | 28 | void* data = | |
| 93 | 28 | type_ == GridSampleType::UInt16 | |
| 94 | 28 | ? static_cast<void*>(stbi_load_16(name_.c_str(), &width, &height, &channels, 0)) | |
| 95 | 28 | : static_cast<void*>(stbi_load(name_.c_str(), &width, &height, &channels, 0)); | |
| 96 | 28 | if (data == nullptr) { | |
| 97 | ✗ | throw std::runtime_error("Cannot read " + name_ + ": " + stbi_failure_reason()); | |
| 98 | } | ||
| 99 | 28 | pixels_.reset(data); | |
| 100 | 28 | if (static_cast<std::uint32_t>(width) != width_ || | |
| 101 | 28 | static_cast<std::uint32_t>(height) != height_ || | |
| 102 | 28 | static_cast<std::uint32_t>(channels) != channels_) { | |
| 103 | ✗ | pixels_.reset(); | |
| 104 | ✗ | throw std::runtime_error(name_ + " changed as it was read"); | |
| 105 | } | ||
| 106 | 28 | } | |
| 107 | |||
| 108 | struct Free { | ||
| 109 | 26 | void operator()(void* data) const { | |
| 110 | 26 | stbi_image_free(data); | |
| 111 | 26 | } | |
| 112 | }; | ||
| 113 | |||
| 114 | std::string name_; | ||
| 115 | std::uint32_t width_ = 0; | ||
| 116 | std::uint32_t height_ = 0; | ||
| 117 | std::uint32_t channels_ = 0; | ||
| 118 | GridSampleType type_ = GridSampleType::UInt8; | ||
| 119 | mutable std::unique_ptr<void, Free> pixels_; | ||
| 120 | }; | ||
| 121 | |||
| 122 | // ---------------------------------------------------------------------------------------------- | ||
| 123 | // PDS3 label syntax | ||
| 124 | |||
| 125 | /// `line` without its comments; `inComment` says whether one is open at its start and at its end. | ||
| 126 | 175 | std::string WithoutComments(const std::string& line, bool& inComment) { | |
| 127 | 175 | std::string text; | |
| 128 | 3813 | for (std::size_t at = 0; at < line.size();) { | |
| 129 | 3644 | if (inComment) { | |
| 130 | 14 | const std::size_t end = line.find("*/", at); | |
| 131 | 14 | if (end == std::string::npos) { | |
| 132 | break; | ||
| 133 | } | ||
| 134 | 8 | inComment = false; | |
| 135 | 8 | at = end + 2; | |
| 136 | 3630 | } else if (line.compare(at, 2, "/*") == 0) { | |
| 137 | 8 | inComment = true; | |
| 138 | 8 | at += 2; | |
| 139 | } else { | ||
| 140 | 3622 | text.push_back(line[at++]); | |
| 141 | } | ||
| 142 | } | ||
| 143 | 175 | return text; | |
| 144 | ✗ | } | |
| 145 | |||
| 146 | /// Whether `value` opens a quote, parenthesis or brace it does not close, so it goes on over the | ||
| 147 | /// next line. | ||
| 148 | 151 | bool Unclosed(const std::string& value) { | |
| 149 | 151 | int depth = 0; | |
| 150 | 151 | bool quoted = false; | |
| 151 | 1477 | for (const char c : value) { | |
| 152 | 1326 | if (c == '"') { | |
| 153 | 40 | quoted = !quoted; | |
| 154 | 1286 | } else if (!quoted && (c == '(' || c == '{')) { | |
| 155 | 3 | ++depth; | |
| 156 | 1010 | } else if (!quoted && (c == ')' || c == '}')) { | |
| 157 | 2 | --depth; | |
| 158 | } | ||
| 159 | } | ||
| 160 | 151 | return quoted || depth > 0; | |
| 161 | } | ||
| 162 | |||
| 163 | // ---------------------------------------------------------------------------------------------- | ||
| 164 | // TIFF layout | ||
| 165 | |||
| 166 | // TIFF tags read. | ||
| 167 | constexpr std::uint16_t kImageWidth = 256; | ||
| 168 | constexpr std::uint16_t kImageLength = 257; | ||
| 169 | constexpr std::uint16_t kBitsPerSample = 258; | ||
| 170 | constexpr std::uint16_t kCompression = 259; | ||
| 171 | constexpr std::uint16_t kStripOffsets = 273; | ||
| 172 | constexpr std::uint16_t kSamplesPerPixel = 277; | ||
| 173 | constexpr std::uint16_t kRowsPerStrip = 278; | ||
| 174 | constexpr std::uint16_t kStripByteCounts = 279; | ||
| 175 | constexpr std::uint16_t kPlanarConfiguration = 284; | ||
| 176 | constexpr std::uint16_t kPredictor = 317; | ||
| 177 | constexpr std::uint16_t kTileWidth = 322; | ||
| 178 | constexpr std::uint16_t kTileLength = 323; | ||
| 179 | constexpr std::uint16_t kTileOffsets = 324; | ||
| 180 | constexpr std::uint16_t kTileByteCounts = 325; | ||
| 181 | constexpr std::uint16_t kSampleFormat = 339; | ||
| 182 | constexpr std::uint16_t kModelPixelScale = 33550; | ||
| 183 | constexpr std::uint16_t kModelTiepoint = 33922; | ||
| 184 | constexpr std::uint16_t kGdalNoData = 42113; | ||
| 185 | |||
| 186 | /// The size of a TIFF field type: BYTE, ASCII, SHORT, LONG, ..., DOUBLE, LONG8; 0 if unknown. | ||
| 187 | 284 | std::size_t TypeBytes(std::uint16_t type) { | |
| 188 | 284 | switch (type) { | |
| 189 | case 1: // BYTE | ||
| 190 | case 2: // ASCII | ||
| 191 | case 6: // SBYTE | ||
| 192 | case 7: // UNDEFINED | ||
| 193 | return 1; | ||
| 194 | 246 | case 3: // SHORT | |
| 195 | 246 | case 8: // SSHORT | |
| 196 | 246 | return 2; | |
| 197 | 20 | case 4: // LONG | |
| 198 | 20 | case 9: // SLONG | |
| 199 | 20 | case 11: // FLOAT | |
| 200 | 20 | return 4; | |
| 201 | 14 | case 5: // RATIONAL | |
| 202 | 14 | case 10: // SRATIONAL | |
| 203 | 14 | case 12: // DOUBLE | |
| 204 | 14 | case 16: // LONG8 | |
| 205 | 14 | case 17: // SLONG8 | |
| 206 | 14 | return 8; | |
| 207 | ✗ | default: | |
| 208 | ✗ | return 0; | |
| 209 | } | ||
| 210 | } | ||
| 211 | |||
| 212 | /// Reads a TIFF's integers in its byte order, checking every access against the file's end. | ||
| 213 | class TiffReader { | ||
| 214 | public: | ||
| 215 | 12 | TiffReader(std::span<const std::byte> bytes, bool bigEndian) | |
| 216 | 12 | : bytes_(bytes), bigEndian_(bigEndian) {} | |
| 217 | |||
| 218 | 1044 | std::uint64_t Unsigned(std::uint64_t at, std::size_t size) const { | |
| 219 | 1044 | if (at > bytes_.size() || size > bytes_.size() - at) { | |
| 220 | ✗ | throw std::runtime_error("TIFF: read past the file's end"); | |
| 221 | } | ||
| 222 | std::uint64_t value = 0; | ||
| 223 | 4198 | for (std::size_t i = 0; i < size; ++i) { | |
| 224 | 3154 | const auto byte = | |
| 225 | 3154 | std::to_integer<std::uint64_t>(bytes_[at + (bigEndian_ ? i : size - 1 - i)]); | |
| 226 | 3154 | value = (value << 8) | byte; | |
| 227 | } | ||
| 228 | 1044 | return value; | |
| 229 | } | ||
| 230 | |||
| 231 | 18 | double Double(std::uint64_t at) const { | |
| 232 | 18 | return std::bit_cast<double>(Unsigned(at, 8)); | |
| 233 | } | ||
| 234 | |||
| 235 | private: | ||
| 236 | std::span<const std::byte> bytes_; | ||
| 237 | bool bigEndian_; | ||
| 238 | }; | ||
| 239 | |||
| 240 | /// One IFD entry: its type, count and where its values are (inline in the entry or at an offset). | ||
| 241 | struct TiffField { | ||
| 242 | std::uint16_t type = 0; | ||
| 243 | std::uint64_t count = 0; | ||
| 244 | std::uint64_t at = 0; | ||
| 245 | }; | ||
| 246 | |||
| 247 | /// TIFF's LZW: 9- to 12-bit codes, most significant bit first, 256 clearing the table and 257 | ||
| 248 | /// ending; widths grow a code early. | ||
| 249 | 26 | std::vector<std::byte> DecodeLzw(std::span<const std::byte> input, std::size_t expected) { | |
| 250 | 26 | constexpr std::uint32_t kClear = 256; | |
| 251 | 26 | constexpr std::uint32_t kEnd = 257; | |
| 252 | 26 | std::vector<std::byte> out; | |
| 253 | 26 | out.reserve(expected); | |
| 254 | 26 | std::array<std::uint16_t, 4096> prefix{}; | |
| 255 | 26 | std::array<std::byte, 4096> suffix{}; | |
| 256 | 26 | std::array<std::byte, 4096> first{}; | |
| 257 | 26 | std::array<std::uint16_t, 4096> length{}; | |
| 258 | 6682 | for (std::uint32_t i = 0; i < 256; ++i) { | |
| 259 | 6656 | suffix[i] = first[i] = static_cast<std::byte>(i); | |
| 260 | 6656 | length[i] = 1; | |
| 261 | } | ||
| 262 | 26 | std::uint64_t bit = 0; | |
| 263 | 26 | const std::uint64_t bits = std::uint64_t{input.size()} * 8; | |
| 264 | 26 | std::uint32_t width = 9; | |
| 265 | 26 | std::uint32_t next = 258; | |
| 266 | 26 | std::int32_t old = -1; | |
| 267 | 1463 | const auto read = [&]() -> std::uint32_t { | |
| 268 | 1437 | if (bit + width > bits) { | |
| 269 | return kEnd; | ||
| 270 | } | ||
| 271 | std::uint32_t code = 0; | ||
| 272 | 14370 | for (std::uint32_t i = 0; i < width; ++i, ++bit) { | |
| 273 | 12933 | const auto byte = std::to_integer<std::uint32_t>(input[bit >> 3]); | |
| 274 | 12933 | code = (code << 1) | ((byte >> (7 - (bit & 7))) & 1); | |
| 275 | } | ||
| 276 | return code; | ||
| 277 | 26 | }; | |
| 278 | 1437 | const auto emit = [&](std::uint32_t code) { | |
| 279 | 1411 | const std::size_t at = out.size(); | |
| 280 | 1411 | out.resize(at + length[code]); | |
| 281 | 14213 | for (std::uint32_t c = code, i = length[code]; i-- > 0;) { | |
| 282 | 12802 | out[at + i] = suffix[c]; | |
| 283 | 12802 | c = prefix[c]; | |
| 284 | } | ||
| 285 | 1437 | }; | |
| 286 | 1411 | const auto add = [&](std::uint32_t from, std::byte byte) { | |
| 287 | 1385 | if (next < 4096) { | |
| 288 | 1385 | prefix[next] = static_cast<std::uint16_t>(from); | |
| 289 | 1385 | suffix[next] = byte; | |
| 290 | 1385 | first[next] = first[from]; | |
| 291 | 1385 | length[next] = static_cast<std::uint16_t>(length[from] + 1); | |
| 292 | 1385 | ++next; | |
| 293 | } | ||
| 294 | 1385 | if (next >= (1u << width) - 1 && width < 12) { | |
| 295 | ✗ | ++width; | |
| 296 | } | ||
| 297 | 1411 | }; | |
| 298 | 1437 | for (std::uint32_t code = read(); code != kEnd; code = read()) { | |
| 299 | 1437 | if (code == kClear) { | |
| 300 | 26 | width = 9; | |
| 301 | 26 | next = 258; | |
| 302 | 26 | old = -1; | |
| 303 | 26 | continue; | |
| 304 | } | ||
| 305 | 1411 | if (old < 0) { | |
| 306 | 26 | if (code > 255) { | |
| 307 | ✗ | throw std::runtime_error("TIFF: damaged LZW data"); | |
| 308 | } | ||
| 309 | 26 | emit(code); | |
| 310 | 1385 | } else if (code < next) { | |
| 311 | 996 | emit(code); | |
| 312 | 996 | add(static_cast<std::uint32_t>(old), first[code]); | |
| 313 | 389 | } else if (code == next) { | |
| 314 | 389 | add(static_cast<std::uint32_t>(old), first[static_cast<std::uint32_t>(old)]); | |
| 315 | 389 | emit(code); | |
| 316 | } else { | ||
| 317 | ✗ | throw std::runtime_error("TIFF: damaged LZW data"); | |
| 318 | } | ||
| 319 | 1411 | old = static_cast<std::int32_t>(code); | |
| 320 | 1411 | if (out.size() >= expected) { | |
| 321 | break; | ||
| 322 | } | ||
| 323 | } | ||
| 324 | 26 | return out; | |
| 325 | ✗ | } | |
| 326 | |||
| 327 | /// zlib-wrapped Deflate, by stb_image's decoder. | ||
| 328 | 77 | std::vector<std::byte> DecodeDeflate(std::span<const std::byte> input, std::size_t expected) { | |
| 329 | 77 | if (expected > static_cast<std::size_t>(INT_MAX) || | |
| 330 | 77 | input.size() > static_cast<std::size_t>(INT_MAX)) { | |
| 331 | ✗ | throw std::runtime_error("TIFF: a Deflate strip over 2 GB"); | |
| 332 | } | ||
| 333 | 77 | std::vector<std::byte> out(expected); | |
| 334 | 77 | const int size = stbi_zlib_decode_buffer( | |
| 335 | 77 | reinterpret_cast<char*>(out.data()), static_cast<int>(expected), | |
| 336 | 77 | reinterpret_cast<const char*>(input.data()), static_cast<int>(input.size())); | |
| 337 | 77 | if (size < 0) { | |
| 338 | ✗ | throw std::runtime_error("TIFF: damaged Deflate data"); | |
| 339 | } | ||
| 340 | 77 | out.resize(static_cast<std::size_t>(size)); | |
| 341 | 77 | return out; | |
| 342 | ✗ | } | |
| 343 | |||
| 344 | /// Undoes predictor 3 in place: the bytes of floats, most significant first, differenced across | ||
| 345 | /// each row. | ||
| 346 | 5 | void UndoFloatingPointPredictor(std::vector<std::byte>& data, std::size_t rowBytes, | |
| 347 | std::uint32_t rows, std::uint32_t channels, | ||
| 348 | std::uint32_t sampleBytes) { | ||
| 349 | 5 | static_assert(std::endian::native == std::endian::little, | |
| 350 | "TIFF floats are reordered for a little-endian host"); | ||
| 351 | 5 | std::vector<std::byte> row(rowBytes); | |
| 352 | 5 | const std::size_t samples = rowBytes / sampleBytes; | |
| 353 | 28 | for (std::uint32_t y = 0; y < rows; ++y) { | |
| 354 | 23 | std::byte* line = data.data() + y * rowBytes; | |
| 355 | 3404 | for (std::size_t i = channels; i < rowBytes; ++i) { | |
| 356 | 3381 | line[i] = static_cast<std::byte>(std::to_integer<unsigned>(line[i]) + | |
| 357 | 3381 | std::to_integer<unsigned>(line[i - channels])); | |
| 358 | } | ||
| 359 | 874 | for (std::size_t i = 0; i < samples; ++i) { | |
| 360 | 4255 | for (std::uint32_t b = 0; b < sampleBytes; ++b) { | |
| 361 | 3404 | row[i * sampleBytes + b] = line[(sampleBytes - 1 - b) * samples + i]; | |
| 362 | } | ||
| 363 | } | ||
| 364 | 23 | std::memcpy(line, row.data(), rowBytes); | |
| 365 | } | ||
| 366 | 5 | } | |
| 367 | |||
| 368 | /// Undoes predictor 2 in place: each sample the difference from its left neighbour's. | ||
| 369 | 26 | void UndoHorizontalPredictor(std::vector<std::byte>& data, std::size_t rowBytes, std::uint32_t rows, | |
| 370 | std::uint32_t channels, std::uint32_t sampleBytes) { | ||
| 371 | 26 | const std::size_t samples = rowBytes / sampleBytes; | |
| 372 | 241 | for (std::uint32_t y = 0; y < rows; ++y) { | |
| 373 | 215 | std::byte* line = data.data() + y * rowBytes; | |
| 374 | 11843 | for (std::size_t i = channels; i < samples; ++i) { | |
| 375 | 11628 | switch (sampleBytes) { | |
| 376 | 10800 | case 1: | |
| 377 | 10800 | line[i] = static_cast<std::byte>(std::to_integer<unsigned>(line[i]) + | |
| 378 | 10800 | std::to_integer<unsigned>(line[i - channels])); | |
| 379 | 10800 | break; | |
| 380 | 828 | case 2: { | |
| 381 | 828 | std::uint16_t a, b; | |
| 382 | 828 | std::memcpy(&a, line + 2 * i, 2); | |
| 383 | 828 | std::memcpy(&b, line + 2 * (i - channels), 2); | |
| 384 | 828 | a = static_cast<std::uint16_t>(a + b); | |
| 385 | 828 | std::memcpy(line + 2 * i, &a, 2); | |
| 386 | 828 | break; | |
| 387 | } | ||
| 388 | ✗ | default: { | |
| 389 | ✗ | std::uint32_t a, b; | |
| 390 | ✗ | std::memcpy(&a, line + 4 * i, 4); | |
| 391 | ✗ | std::memcpy(&b, line + 4 * (i - channels), 4); | |
| 392 | ✗ | a += b; | |
| 393 | ✗ | std::memcpy(line + 4 * i, &a, 4); | |
| 394 | ✗ | break; | |
| 395 | } | ||
| 396 | } | ||
| 397 | } | ||
| 398 | } | ||
| 399 | 26 | } | |
| 400 | |||
| 401 | // ---------------------------------------------------------------------------------------------- | ||
| 402 | // DEMs | ||
| 403 | |||
| 404 | /// A DEM's grid: its size, its heights (metres above the reference sphere) and that sphere's | ||
| 405 | /// radius. | ||
| 406 | ✗ | struct DemGrid { | |
| 407 | std::uint32_t lines = 0; | ||
| 408 | std::uint32_t samples = 0; | ||
| 409 | double referenceRadius = 0.0; | ||
| 410 | std::vector<float> heights; | ||
| 411 | }; | ||
| 412 | |||
| 413 | /// Reads a DEM's heights from its PDS image: samples times the scale, plus the offset unless the | ||
| 414 | /// offset is the reference sphere's radius (as LOLA's is; the projection's A axis is otherwise | ||
| 415 | /// the sphere, in km). PDS's real null, the missing constant and non-finite samples are NaN. | ||
| 416 | 4 | DemGrid ReadDemGrid(const PdsImage& image) { | |
| 417 | 4 | const double offset = image.Offset(); | |
| 418 | 4 | const bool offsetIsRadius = offset > 1e5; | |
| 419 | 4 | DemGrid dem; | |
| 420 | 4 | dem.lines = image.Height(); | |
| 421 | 4 | dem.samples = image.Width(); | |
| 422 | 8 | dem.referenceRadius = offsetIsRadius | |
| 423 | 5 | ? offset | |
| 424 | 1 | : 1e3 * image.Keywords().Number("IMAGE_MAP_PROJECTION.A_AXIS_RADIUS"); | |
| 425 | 4 | const double addend = offsetIsRadius ? 0.0 : offset; | |
| 426 | 4 | const double scale = image.Scale(); | |
| 427 | 4 | const std::optional<double> noData = image.NoData(); | |
| 428 | 4 | const GridSampleType type = image.Type(); | |
| 429 | 4 | const std::uint32_t bytes = image.SampleBytes(); | |
| 430 | 4 | const std::vector<std::byte> raw = image.ReadAll(); | |
| 431 | 4 | dem.heights.resize(std::size_t{dem.lines} * dem.samples); | |
| 432 | 645 | for (std::size_t i = 0; i < dem.heights.size(); ++i) { | |
| 433 | 641 | const double value = DecodeSample(raw.data() + i * bytes, type); | |
| 434 | 641 | const bool missing = | |
| 435 | 641 | (noData && value == *noData) || | |
| 436 | 8 | (type == GridSampleType::Float32 && !(value > -1e30 && std::isfinite(value))); | |
| 437 | 641 | dem.heights[i] = missing ? std::numeric_limits<float>::quiet_NaN() | |
| 438 | 640 | : static_cast<float>(value * scale + addend); | |
| 439 | } | ||
| 440 | 4 | return dem; | |
| 441 | 4 | } | |
| 442 | |||
| 443 | /// The map projection type a DEM's label gives, checked against those accepted. | ||
| 444 | 6 | std::string DemProjection(const PdsImage& image, std::initializer_list<std::string_view> accepted) { | |
| 445 | 6 | const std::string projection = | |
| 446 | 6 | image.Keywords().Text("IMAGE_MAP_PROJECTION.MAP_PROJECTION_TYPE"); | |
| 447 | 6 | if (std::find(accepted.begin(), accepted.end(), projection) == accepted.end()) { | |
| 448 | 6 | throw std::invalid_argument("The DEM's projection '" + projection + "' is not " + | |
| 449 | 6 | std::string(*accepted.begin())); | |
| 450 | } | ||
| 451 | 4 | return projection; | |
| 452 | 2 | } | |
| 453 | |||
| 454 | /// `heights` (`lines` rows of `samples`) at fractional line `line` and sample `sample` (cell | ||
| 455 | /// centres at whole numbers), bilinear, clamped to the grid's edges; 0 for an empty grid. Corners | ||
| 456 | /// of no weight are left out, so a cell without data (NaN) reaches only the points beside it. | ||
| 457 | 1880 | double SampleDem(const std::vector<float>& heights, std::uint32_t lines, std::uint32_t samples, | |
| 458 | double line, double sample) { | ||
| 459 | 1880 | if (lines == 0 || samples == 0) { | |
| 460 | return 0.0; | ||
| 461 | } | ||
| 462 | 1880 | const double y = std::clamp(line, 0.0, double(lines - 1)); | |
| 463 | 1880 | const double x = std::clamp(sample, 0.0, double(samples - 1)); | |
| 464 | 1880 | const auto l0 = | |
| 465 | 1880 | std::min<std::uint32_t>(static_cast<std::uint32_t>(y), lines > 1 ? lines - 2 : 0); | |
| 466 | 1880 | const auto s0 = | |
| 467 | 1880 | std::min<std::uint32_t>(static_cast<std::uint32_t>(x), samples > 1 ? samples - 2 : 0); | |
| 468 | 1880 | const std::uint32_t l1 = std::min(l0 + 1, lines - 1); | |
| 469 | 1880 | const std::uint32_t s1 = std::min(s0 + 1, samples - 1); | |
| 470 | 1880 | const double ty = y - l0; | |
| 471 | 1880 | const double tx = x - s0; | |
| 472 | 1880 | const auto at = [&](std::uint32_t l, std::uint32_t s) { | |
| 473 | 1880 | return static_cast<double>(heights[std::size_t(l) * samples + s]); | |
| 474 | 1880 | }; | |
| 475 | 1880 | const double weights[4] = {(1.0 - tx) * (1.0 - ty), tx * (1.0 - ty), (1.0 - tx) * ty, tx * ty}; | |
| 476 | 1880 | const double values[4] = {at(l0, s0), at(l0, s1), at(l1, s0), at(l1, s1)}; | |
| 477 | 1880 | double sum = 0.0; | |
| 478 | 9400 | for (int i = 0; i < 4; ++i) { | |
| 479 | 7520 | if (weights[i] > 0.0) { | |
| 480 | 7293 | sum += weights[i] * values[i]; | |
| 481 | } | ||
| 482 | } | ||
| 483 | return sum; | ||
| 484 | } | ||
| 485 | |||
| 486 | /// The local terrain about a point of a sphere of `radius`, its surface `heightAt(latitude, | ||
| 487 | /// longitude)` above it (`MakeLocalTerrain`). | ||
| 488 | template <typename HeightAt> | ||
| 489 | 4 | LocalTerrain LocalTerrainAround(double radius, const HeightAt& heightAt, double latitudeDeg, | |
| 490 | double longitudeDeg, std::uint32_t cells, double cellSize) { | ||
| 491 | 4 | if (cells < 2 || !(cellSize > 0.0)) { | |
| 492 | 1 | throw std::invalid_argument("A local terrain needs at least 2 cells a positive size apart"); | |
| 493 | } | ||
| 494 | // The local frame: up along the point's direction, east and north along its parallel and | ||
| 495 | // meridian. | ||
| 496 | 3 | const glm::dvec3 up = BodyDirection(latitudeDeg, longitudeDeg); | |
| 497 | 3 | const double lon = glm::radians(longitudeDeg); | |
| 498 | 3 | const glm::dvec3 east(-std::sin(lon), std::cos(lon), 0.0); | |
| 499 | 3 | const glm::dvec3 north = glm::cross(up, east); | |
| 500 | 3 | const glm::dvec3 centre = radius * up; | |
| 501 | 3 | LocalTerrain terrain; | |
| 502 | 3 | glm::dmat4 localToBody(1.0); | |
| 503 | 3 | localToBody[0] = glm::dvec4(east, 0.0); | |
| 504 | 3 | localToBody[1] = glm::dvec4(north, 0.0); | |
| 505 | 3 | localToBody[2] = glm::dvec4(up, 0.0); | |
| 506 | 3 | localToBody[3] = glm::dvec4(centre, 1.0); | |
| 507 | 3 | terrain.bodyToLocal = glm::inverse(localToBody); | |
| 508 | 3 | terrain.cellSize = cellSize; | |
| 509 | 3 | const double half = 0.5 * (cells - 1) * cellSize; | |
| 510 | 3 | terrain.origin = glm::dvec2(-half, half); | |
| 511 | |||
| 512 | 3 | std::vector<float> heights(std::size_t(cells) * cells); | |
| 513 | 30 | for (std::uint32_t row = 0; row < cells; ++row) { | |
| 514 | 486 | for (std::uint32_t col = 0; col < cells; ++col) { | |
| 515 | 459 | const double x = -half + col * cellSize; | |
| 516 | 459 | const double y = half - row * cellSize; | |
| 517 | // The surface point whose local x and y are these: on the sphere of radius R + h | ||
| 518 | // through them, h the DEM's height where it lands, settled in a few rounds. | ||
| 519 | 459 | double h = 0.0; | |
| 520 | 459 | double z = 0.0; | |
| 521 | 2295 | for (int i = 0; i < 4; ++i) { | |
| 522 | 1836 | const double r = radius + h; | |
| 523 | 1836 | z = std::sqrt(std::max(r * r - x * x - y * y, 0.0)) - radius; | |
| 524 | 1836 | const glm::dvec3 d = glm::normalize(centre + x * east + y * north + z * up); | |
| 525 | 3672 | const double lat = glm::degrees(std::asin(std::clamp(d.z, -1.0, 1.0))); | |
| 526 | 1836 | const double along = glm::degrees(std::atan2(d.y, d.x)); | |
| 527 | 1836 | h = heightAt(lat, along); | |
| 528 | 1836 | if (std::isnan(h)) { | |
| 529 | ✗ | h = 0.0; | |
| 530 | } | ||
| 531 | } | ||
| 532 | 459 | const double r = radius + h; | |
| 533 | 459 | z = std::sqrt(std::max(r * r - x * x - y * y, 0.0)) - radius; | |
| 534 | 459 | heights[std::size_t(row) * cells + col] = static_cast<float>(z); | |
| 535 | } | ||
| 536 | } | ||
| 537 | 3 | terrain.heights = Tensor::Copy<float>(heights, {std::int64_t(cells), std::int64_t(cells)}); | |
| 538 | 3 | return terrain; | |
| 539 | 3 | } | |
| 540 | |||
| 541 | } // namespace | ||
| 542 | |||
| 543 | // ---------------------------------------------------------------------------------------------- | ||
| 544 | // Rasters | ||
| 545 | |||
| 546 | 433 | std::uint32_t ByteSize(GridSampleType type) { | |
| 547 | 433 | switch (type) { | |
| 548 | case GridSampleType::UInt8: | ||
| 549 | return 1; | ||
| 550 | 26 | case GridSampleType::UInt16: | |
| 551 | 26 | case GridSampleType::Int16: | |
| 552 | 26 | return 2; | |
| 553 | case GridSampleType::UInt32: | ||
| 554 | case GridSampleType::Int32: | ||
| 555 | case GridSampleType::Float32: | ||
| 556 | break; | ||
| 557 | } | ||
| 558 | 47 | return 4; | |
| 559 | } | ||
| 560 | |||
| 561 | 643 | double DecodeSample(const std::byte* sample, GridSampleType type) { | |
| 562 | 643 | switch (type) { | |
| 563 | ✗ | case GridSampleType::UInt8: | |
| 564 | ✗ | return LoadAs<std::uint8_t>(sample); | |
| 565 | ✗ | case GridSampleType::UInt16: | |
| 566 | ✗ | return LoadAs<std::uint16_t>(sample); | |
| 567 | 634 | case GridSampleType::Int16: | |
| 568 | 634 | return LoadAs<std::int16_t>(sample); | |
| 569 | ✗ | case GridSampleType::UInt32: | |
| 570 | ✗ | return LoadAs<std::uint32_t>(sample); | |
| 571 | ✗ | case GridSampleType::Int32: | |
| 572 | ✗ | return LoadAs<std::int32_t>(sample); | |
| 573 | case GridSampleType::Float32: | ||
| 574 | break; | ||
| 575 | } | ||
| 576 | 9 | return LoadAs<float>(sample); | |
| 577 | } | ||
| 578 | |||
| 579 | 12 | std::vector<std::byte> GridSource::ReadAll() const { | |
| 580 | 12 | std::vector<std::byte> out(std::size_t{Width()} * Height() * Channels() * SampleBytes()); | |
| 581 | 12 | Read(0, 0, Width(), Height(), out); | |
| 582 | 12 | return out; | |
| 583 | ✗ | } | |
| 584 | |||
| 585 | 77 | void GridSource::CheckRead(std::uint32_t x, std::uint32_t y, std::uint32_t w, std::uint32_t h, | |
| 586 | std::span<const std::byte> out) const { | ||
| 587 | 77 | if (x > Width() || w > Width() - x || y > Height() || h > Height() - y) { | |
| 588 | 2 | throw std::out_of_range("A rectangle past the raster"); | |
| 589 | } | ||
| 590 | 75 | if (out.size() < std::size_t{w} * h * Channels() * SampleBytes()) { | |
| 591 | 1 | throw std::out_of_range("Too little room for the rectangle"); | |
| 592 | } | ||
| 593 | 74 | } | |
| 594 | |||
| 595 | 9 | MosaicRaster::MosaicRaster(const std::vector<fs::path>& tiles, std::uint32_t columns) | |
| 596 | 9 | : columns_(columns) { | |
| 597 | 9 | if (columns == 0 || tiles.empty() || tiles.size() % columns != 0) { | |
| 598 | 2 | throw std::invalid_argument("A mosaic needs whole rows of " + std::to_string(columns) + | |
| 599 | 3 | " tiles"); | |
| 600 | } | ||
| 601 | 8 | rows_ = static_cast<std::uint32_t>(tiles.size() / columns); | |
| 602 | 34 | for (const fs::path& tile : tiles) { | |
| 603 | 56 | tiles_.push_back(OpenRaster(tile)); | |
| 604 | 28 | const GridSource& first = *tiles_.front(); | |
| 605 | 28 | const GridSource& last = *tiles_.back(); | |
| 606 | 80 | if (last.Width() != first.Width() || last.Height() != first.Height() || | |
| 607 | 80 | last.Channels() != first.Channels() || last.Type() != first.Type()) { | |
| 608 | 4 | throw std::invalid_argument("A mosaic's tiles are alike: " + tile.string() + | |
| 609 | 8 | " is not as " + tiles.front().string()); | |
| 610 | } | ||
| 611 | } | ||
| 612 | 6 | tileWidth_ = tiles_.front()->Width(); | |
| 613 | 6 | tileHeight_ = tiles_.front()->Height(); | |
| 614 | 9 | } | |
| 615 | |||
| 616 | 9 | std::vector<GridRect> MosaicRaster::Parts() const { | |
| 617 | 9 | std::vector<GridRect> parts; | |
| 618 | 27 | for (std::uint32_t row = 0; row < rows_; ++row) { | |
| 619 | 56 | for (std::uint32_t column = 0; column < columns_; ++column) { | |
| 620 | 38 | parts.push_back({column * tileWidth_, row * tileHeight_, tileWidth_, tileHeight_}); | |
| 621 | } | ||
| 622 | } | ||
| 623 | 9 | return parts; | |
| 624 | ✗ | } | |
| 625 | |||
| 626 | 4 | void MosaicRaster::Release() const { | |
| 627 | 20 | for (const std::unique_ptr<GridSource>& tile : tiles_) { | |
| 628 | 16 | tile->Release(); | |
| 629 | } | ||
| 630 | 4 | } | |
| 631 | |||
| 632 | 18 | void MosaicRaster::Read(std::uint32_t x, std::uint32_t y, std::uint32_t w, std::uint32_t h, | |
| 633 | std::span<std::byte> out) const { | ||
| 634 | 18 | CheckRead(x, y, w, h, out); | |
| 635 | 18 | const std::size_t pixel = std::size_t{Channels()} * SampleBytes(); | |
| 636 | 18 | std::vector<std::byte> part; | |
| 637 | 37 | for (std::uint32_t row = y / tileHeight_; row < rows_ && row * tileHeight_ < y + h; ++row) { | |
| 638 | 19 | for (std::uint32_t column = x / tileWidth_; | |
| 639 | 43 | column < columns_ && column * tileWidth_ < x + w; ++column) { | |
| 640 | 24 | const std::size_t index = std::size_t{row} * columns_ + column; | |
| 641 | 24 | if (index != current_) { | |
| 642 | 18 | tiles_[current_]->Release(); | |
| 643 | 18 | current_ = index; | |
| 644 | } | ||
| 645 | // The rectangle's overlap with the tile, in the tile. | ||
| 646 | 24 | const std::uint32_t left = std::max(x, column * tileWidth_); | |
| 647 | 24 | const std::uint32_t top = std::max(y, row * tileHeight_); | |
| 648 | 24 | const std::uint32_t right = std::min(x + w, (column + 1) * tileWidth_); | |
| 649 | 24 | const std::uint32_t bottom = std::min(y + h, (row + 1) * tileHeight_); | |
| 650 | 24 | const std::uint32_t across = right - left; | |
| 651 | 24 | part.resize(std::size_t{across} * (bottom - top) * pixel); | |
| 652 | 24 | tiles_[index]->Read(left - column * tileWidth_, top - row * tileHeight_, across, | |
| 653 | bottom - top, part); | ||
| 654 | 668 | for (std::uint32_t r = top; r < bottom; ++r) { | |
| 655 | 644 | std::memcpy(out.data() + (std::size_t{r - y} * w + (left - x)) * pixel, | |
| 656 | 644 | part.data() + std::size_t{r - top} * across * pixel, | |
| 657 | std::size_t{across} * pixel); | ||
| 658 | } | ||
| 659 | } | ||
| 660 | } | ||
| 661 | 18 | } | |
| 662 | |||
| 663 | 4 | std::unique_ptr<MosaicRaster> OpenMosaic(const fs::path& path) { | |
| 664 | 4 | std::ifstream file(path); | |
| 665 | 4 | if (!file) { | |
| 666 | ✗ | throw std::runtime_error("Cannot read " + path.string()); | |
| 667 | } | ||
| 668 | 4 | std::vector<fs::path> tiles; | |
| 669 | 4 | std::size_t columns = 0; | |
| 670 | 4 | std::string line; | |
| 671 | 12 | while (std::getline(file, line)) { | |
| 672 | 10 | if (const std::size_t comment = line.find('#'); comment != std::string::npos) { | |
| 673 | 3 | line.resize(comment); | |
| 674 | } | ||
| 675 | 10 | std::istringstream names(line); | |
| 676 | 10 | std::size_t count = 0; | |
| 677 | 24 | for (std::string name; names >> name; ++count) { | |
| 678 | 70 | tiles.push_back(path.parent_path() / name); | |
| 679 | ✗ | } | |
| 680 | 10 | if (count == 0) { | |
| 681 | 3 | continue; | |
| 682 | } | ||
| 683 | 7 | if (columns != 0 && count != columns) { | |
| 684 | 8 | throw std::invalid_argument(path.string() + ": a row of " + std::to_string(count) + | |
| 685 | 8 | " tiles, not " + std::to_string(columns)); | |
| 686 | } | ||
| 687 | 5 | columns = count; | |
| 688 | 10 | } | |
| 689 | 6 | return std::make_unique<MosaicRaster>(tiles, static_cast<std::uint32_t>(columns)); | |
| 690 | 6 | } | |
| 691 | |||
| 692 | 43 | std::unique_ptr<GridSource> OpenRaster(const fs::path& path) { | |
| 693 | 129 | const std::string extension = ToLower(path.extension().string()); | |
| 694 | 43 | if (extension == ".tif" || extension == ".tiff") { | |
| 695 | 4 | return std::make_unique<TiffImage>(path); | |
| 696 | } | ||
| 697 | 39 | if (extension == ".img" || extension == ".lbl") { | |
| 698 | ✗ | return std::make_unique<PdsImage>(path); | |
| 699 | } | ||
| 700 | 39 | if (extension == ".mosaic") { | |
| 701 | 4 | return OpenMosaic(path); | |
| 702 | } | ||
| 703 | 35 | if (IsPngOrJpeg(path)) { | |
| 704 | 37 | return std::make_unique<DecodedRaster>(path); | |
| 705 | } | ||
| 706 | 1 | throw std::runtime_error( | |
| 707 | 3 | path.string() + " is not a raster read here (TIFF, PDS, PNG, JPEG or a mosaic of them)"); | |
| 708 | 40 | } | |
| 709 | |||
| 710 | // ---------------------------------------------------------------------------------------------- | ||
| 711 | // Grid views | ||
| 712 | |||
| 713 | ✗ | void CheckGridView(const GridView& view, const std::string& what, std::size_t minSide) { | |
| 714 | ✗ | if (view.width < minSide || view.height < minSide) { | |
| 715 | ✗ | throw std::invalid_argument(what + " must be at least " + std::to_string(minSide) + " x " + | |
| 716 | ✗ | std::to_string(minSide) + " cells"); | |
| 717 | } | ||
| 718 | ✗ | if (view.data == nullptr) { | |
| 719 | ✗ | throw std::invalid_argument(what + " has no data"); | |
| 720 | } | ||
| 721 | ✗ | if (view.Stride() < view.width) { | |
| 722 | ✗ | throw std::invalid_argument(what + " row stride is shorter than its width"); | |
| 723 | } | ||
| 724 | ✗ | if (!std::isfinite(view.originX) || !std::isfinite(view.originY) || | |
| 725 | ✗ | !std::isfinite(view.cellSize) || !(view.cellSize > 0.0)) { | |
| 726 | ✗ | throw std::invalid_argument(what + " needs a finite origin and a positive cell size"); | |
| 727 | } | ||
| 728 | ✗ | } | |
| 729 | |||
| 730 | // ---------------------------------------------------------------------------------------------- | ||
| 731 | // PDS3 | ||
| 732 | |||
| 733 | 137 | std::string CleanPds3Value(std::string_view value) { | |
| 734 | 137 | std::string text(value); | |
| 735 | 137 | if (const std::size_t unit = text.find('<'); | |
| 736 | 137 | unit != std::string::npos && text.front() != '"' && text.front() != '(') { | |
| 737 | 46 | text = Trim(text.substr(0, unit)); | |
| 738 | } | ||
| 739 | 137 | if (text.size() >= 2 && (text.front() == '"' || text.front() == '\'') && | |
| 740 | 18 | text.back() == text.front()) { | |
| 741 | 18 | text = text.substr(1, text.size() - 2); | |
| 742 | } | ||
| 743 | 137 | return Trim(text); | |
| 744 | 137 | } | |
| 745 | |||
| 746 | 36 | std::optional<std::string> Pds3Label::Get(std::string_view name) const { | |
| 747 | 136 | for (auto it = keywords.rbegin(); it != keywords.rend(); ++it) { | |
| 748 | 133 | if (name.size() == it->object.size() + it->key.size() && name.starts_with(it->object) && | |
| 749 | 35 | name.ends_with(it->key)) { | |
| 750 | 66 | return CleanPds3Value(it->value); | |
| 751 | } | ||
| 752 | } | ||
| 753 | 3 | return std::nullopt; | |
| 754 | } | ||
| 755 | |||
| 756 | 20 | double Pds3Label::Number(std::string_view name) const { | |
| 757 | 20 | const std::string text = Text(name); | |
| 758 | 19 | try { | |
| 759 | 19 | return std::stod(text); | |
| 760 | 1 | } catch (const std::exception&) { | |
| 761 | 3 | throw std::invalid_argument("PDS3 keyword " + std::string(name) + " = '" + text + | |
| 762 | 3 | "' is not a number"); | |
| 763 | 1 | } | |
| 764 | 18 | } | |
| 765 | |||
| 766 | 32 | std::string Pds3Label::Text(std::string_view name) const { | |
| 767 | 32 | if (std::optional<std::string> value = Get(name)) { | |
| 768 | 31 | return *std::move(value); | |
| 769 | 31 | } | |
| 770 | 2 | throw std::invalid_argument("The PDS3 label has no " + std::string(name)); | |
| 771 | } | ||
| 772 | |||
| 773 | 125 | std::optional<std::string> Pds3Label::Find(std::string_view key) const { | |
| 774 | 125 | const std::string wanted = ToUpper(key); | |
| 775 | 1101 | for (const Keyword& keyword : keywords) { | |
| 776 | 1070 | if (ToUpper(keyword.key) == wanted) { | |
| 777 | 94 | return keyword.value; | |
| 778 | } | ||
| 779 | } | ||
| 780 | 31 | return std::nullopt; | |
| 781 | 125 | } | |
| 782 | |||
| 783 | 10 | Pds3Label ParsePds3Label(std::string_view text) { | |
| 784 | 10 | Pds3Label label; | |
| 785 | 10 | std::vector<std::string> objects; | |
| 786 | 10 | std::istringstream lines{std::string(text)}; | |
| 787 | 10 | std::string line; | |
| 788 | 10 | std::string key; | |
| 789 | 10 | std::string value; | |
| 790 | 10 | bool inComment = false; | |
| 791 | 10 | bool continued = false; | |
| 792 | 168 | const auto finish = [&] { | |
| 793 | 158 | if (key == "OBJECT") { | |
| 794 | 36 | objects.push_back(CleanPds3Value(value)); | |
| 795 | 140 | } else if (key == "END_OBJECT") { | |
| 796 | 17 | if (!objects.empty()) { | |
| 797 | 17 | objects.pop_back(); | |
| 798 | } | ||
| 799 | 123 | } else if (!key.empty()) { | |
| 800 | 113 | std::string object; | |
| 801 | 200 | for (const std::string& name : objects) { | |
| 802 | 174 | object += name + "."; | |
| 803 | } | ||
| 804 | 113 | label.keywords.push_back({std::move(object), std::move(key), Trim(value)}); | |
| 805 | 113 | } | |
| 806 | 158 | key.clear(); | |
| 807 | 158 | value.clear(); | |
| 808 | 158 | continued = false; | |
| 809 | 271 | }; | |
| 810 | 175 | while (std::getline(lines, line)) { | |
| 811 | 175 | line = WithoutComments(line, inComment); | |
| 812 | 175 | if (continued) { | |
| 813 | 6 | value += ' ' + Trim(line); | |
| 814 | } else { | ||
| 815 | 172 | const std::string trimmed = Trim(line); | |
| 816 | 172 | if (trimmed == "END") { | |
| 817 | break; | ||
| 818 | } | ||
| 819 | 162 | const std::size_t equals = trimmed.find('='); | |
| 820 | 162 | if (equals == std::string::npos) { | |
| 821 | 14 | continue; | |
| 822 | } | ||
| 823 | 148 | key = Trim(trimmed.substr(0, equals)); | |
| 824 | 148 | value = Trim(trimmed.substr(equals + 1)); | |
| 825 | 172 | } | |
| 826 | 151 | continued = Unclosed(value); | |
| 827 | 151 | if (!continued) { | |
| 828 | 148 | finish(); | |
| 829 | } | ||
| 830 | } | ||
| 831 | 10 | finish(); | |
| 832 | 20 | return label; | |
| 833 | 10 | } | |
| 834 | |||
| 835 | 9 | PdsImage::PdsImage(const fs::path& path) { | |
| 836 | 9 | const std::string name = path.string(); | |
| 837 | // The label: this file if it is one, else one beside the data, else attached at the data's | ||
| 838 | // start. | ||
| 839 | 9 | fs::path labelPath; | |
| 840 | 27 | if (ToUpper(path.extension().string()) == ".LBL") { | |
| 841 | 7 | labelPath = path; | |
| 842 | 2 | } else if (const fs::path beside = | |
| 843 | 9 | FindFileAnyCase(path.parent_path(), path.stem().string() + ".lbl"); | |
| 844 | 2 | fs::exists(beside)) { | |
| 845 | 1 | labelPath = beside; | |
| 846 | 2 | } | |
| 847 | 9 | std::string text; | |
| 848 | 9 | const bool attached = labelPath.empty(); | |
| 849 | 9 | if (!attached) { | |
| 850 | 8 | text = ReadTextFile(labelPath); | |
| 851 | } else { | ||
| 852 | 1 | std::ifstream file(path, std::ios::binary); | |
| 853 | 1 | text.resize(1 << 16); | |
| 854 | 1 | file.read(text.data(), static_cast<std::streamsize>(text.size())); | |
| 855 | 1 | text.resize(static_cast<std::size_t>(file.gcount())); | |
| 856 | 1 | } | |
| 857 | 8 | if (text.find("PDS_VERSION_ID") == std::string::npos) { | |
| 858 | ✗ | throw std::runtime_error(name + " has no PDS3 label"); | |
| 859 | } | ||
| 860 | 8 | label_ = ParsePds3Label(text); | |
| 861 | 82 | const auto number = [&](std::string_view key, | |
| 862 | std::optional<double> fallback = std::nullopt) -> double { | ||
| 863 | 74 | const std::optional<std::string> value = Label(key); | |
| 864 | 74 | if (!value) { | |
| 865 | 21 | if (fallback) { | |
| 866 | 21 | return *fallback; | |
| 867 | } | ||
| 868 | ✗ | throw std::runtime_error(name + ": the PDS label has no " + std::string(key)); | |
| 869 | } | ||
| 870 | 53 | try { | |
| 871 | 74 | return std::stod(*value); | |
| 872 | ✗ | } catch (const std::exception&) { | |
| 873 | ✗ | throw std::runtime_error(name + ": the PDS label's " + std::string(key) + | |
| 874 | ✗ | " is not a number"); | |
| 875 | ✗ | } | |
| 876 | 74 | }; | |
| 877 | 8 | width_ = static_cast<std::uint32_t>(number("LINE_SAMPLES")); | |
| 878 | 8 | height_ = static_cast<std::uint32_t>(number("LINES")); | |
| 879 | 8 | if (number("BANDS", 1.0) != 1.0) { | |
| 880 | ✗ | throw std::runtime_error(name + ": PDS images of several bands are not read"); | |
| 881 | } | ||
| 882 | 8 | const auto bits = static_cast<std::uint32_t>(number("SAMPLE_BITS")); | |
| 883 | 17 | const std::string sampleType = ToUpper(Label("SAMPLE_TYPE").value_or("")); | |
| 884 | 8 | const bool real = sampleType.find("REAL") != std::string::npos; | |
| 885 | 8 | const bool unsignedType = sampleType.find("UNSIGNED") != std::string::npos; | |
| 886 | 9 | bigEndian_ = sampleType.rfind("LSB", 0) != 0 && sampleType.rfind("PC_", 0) != 0 && | |
| 887 | 1 | sampleType.rfind("VAX", 0) != 0; | |
| 888 | 8 | if (real && bits == 32) { | |
| 889 | 2 | type_ = GridSampleType::Float32; | |
| 890 | 6 | } else if (!real && bits == 8) { | |
| 891 | ✗ | type_ = GridSampleType::UInt8; | |
| 892 | 6 | } else if (!real && bits == 16) { | |
| 893 | 12 | type_ = unsignedType ? GridSampleType::UInt16 : GridSampleType::Int16; | |
| 894 | ✗ | } else if (!real && bits == 32) { | |
| 895 | ✗ | type_ = unsignedType ? GridSampleType::UInt32 : GridSampleType::Int32; | |
| 896 | } else { | ||
| 897 | ✗ | throw std::runtime_error(name + ": PDS samples " + sampleType + " of " + | |
| 898 | ✗ | std::to_string(bits) + " bits are not read"); | |
| 899 | } | ||
| 900 | 8 | scale_ = number("SCALING_FACTOR", 1.0); | |
| 901 | 8 | offset_ = number("OFFSET", 0.0); | |
| 902 | 8 | if (Label("MISSING_CONSTANT")) { | |
| 903 | 1 | noData_ = number("MISSING_CONSTANT"); | |
| 904 | } | ||
| 905 | |||
| 906 | // The data: ^IMAGE names a file beside the label ("NAME.IMG", or ("NAME.IMG", record)), or a | ||
| 907 | // record or byte of this one. | ||
| 908 | 8 | const std::string pointer = label_.Find("^IMAGE").value_or(std::string()); | |
| 909 | 8 | fs::path dataPath = path; | |
| 910 | 8 | std::uint64_t record = 1; | |
| 911 | 8 | bool bytes = false; | |
| 912 | 8 | if (const std::size_t quote = pointer.find('"'); quote != std::string::npos) { | |
| 913 | 7 | const std::size_t end = pointer.find('"', quote + 1); | |
| 914 | 21 | dataPath = FindFileAnyCase(path.parent_path(), pointer.substr(quote + 1, end - quote - 1)); | |
| 915 | 7 | if (const std::size_t comma = pointer.find(',', end); comma != std::string::npos) { | |
| 916 | ✗ | const std::string rest = pointer.substr(comma + 1); | |
| 917 | ✗ | record = std::stoull(rest); | |
| 918 | ✗ | bytes = ToUpper(rest).find("<BYTES>") != std::string::npos; | |
| 919 | ✗ | } | |
| 920 | 1 | } else if (!pointer.empty()) { | |
| 921 | 1 | record = std::stoull(pointer); | |
| 922 | 1 | bytes = ToUpper(pointer).find("<BYTES>") != std::string::npos; | |
| 923 | ✗ | } else if (!attached) { | |
| 924 | ✗ | throw std::runtime_error(name + ": the PDS label has no ^IMAGE"); | |
| 925 | } | ||
| 926 | 8 | const auto recordBytes = static_cast<std::uint64_t>(number("RECORD_BYTES", 0.0)); | |
| 927 | 8 | start_ = bytes ? record - 1 : (record - 1) * recordBytes; | |
| 928 | 8 | file_ = std::make_shared<const MappedFile>(dataPath); | |
| 929 | 8 | if (file_->Bytes().size() < start_ + std::uint64_t{width_} * height_ * SampleBytes()) { | |
| 930 | ✗ | throw std::runtime_error(dataPath.string() + " is shorter than its PDS label says"); | |
| 931 | } | ||
| 932 | |||
| 933 | // Simple cylindrical grids. | ||
| 934 | 16 | const std::string projection = ToUpper(Label("MAP_PROJECTION_TYPE").value_or("")); | |
| 935 | 9 | if ((projection == "SIMPLE CYLINDRICAL" || projection == "EQUIRECTANGULAR") && | |
| 936 | 25 | number("CENTER_LATITUDE", 0.0) == 0.0 && Label("WESTERNMOST_LONGITUDE") && | |
| 937 | 12 | Label("MAXIMUM_LATITUDE")) { | |
| 938 | 4 | LatLongGrid grid; | |
| 939 | 4 | grid.west = number("WESTERNMOST_LONGITUDE"); | |
| 940 | 4 | grid.north = number("MAXIMUM_LATITUDE"); | |
| 941 | 4 | if (Label("MAP_RESOLUTION")) { | |
| 942 | 4 | grid.degreesAcross = grid.degreesDown = 1.0 / number("MAP_RESOLUTION"); | |
| 943 | } else { | ||
| 944 | ✗ | grid.degreesAcross = (number("EASTERNMOST_LONGITUDE") - grid.west) / width_; | |
| 945 | ✗ | grid.degreesDown = (grid.north - number("MINIMUM_LATITUDE")) / height_; | |
| 946 | } | ||
| 947 | 8 | grid_ = grid; | |
| 948 | } | ||
| 949 | 27 | } | |
| 950 | |||
| 951 | 111 | std::optional<std::string> PdsImage::Label(std::string_view key) const { | |
| 952 | 111 | const std::optional<std::string> value = label_.Find(key); | |
| 953 | 111 | if (!value) { | |
| 954 | 29 | return std::nullopt; | |
| 955 | } | ||
| 956 | 164 | return CleanPds3Value(*value); | |
| 957 | 111 | } | |
| 958 | |||
| 959 | 7 | void PdsImage::Read(std::uint32_t x, std::uint32_t y, std::uint32_t w, std::uint32_t h, | |
| 960 | std::span<std::byte> out) const { | ||
| 961 | 7 | CheckRead(x, y, w, h, out); | |
| 962 | 6 | const std::uint32_t bytes = SampleBytes(); | |
| 963 | 6 | const std::span<const std::byte> data = file_->Bytes(); | |
| 964 | 49 | for (std::uint32_t row = 0; row < h; ++row) { | |
| 965 | 43 | const std::span<std::byte> line = | |
| 966 | 43 | out.subspan(std::size_t{row} * w * bytes, std::size_t{w} * bytes); | |
| 967 | 43 | std::memcpy(line.data(), | |
| 968 | 43 | data.data() + start_ + (std::uint64_t{y + row} * width_ + x) * bytes, | |
| 969 | line.size()); | ||
| 970 | 86 | if (NeedsByteSwap(bigEndian_, bytes)) { | |
| 971 | 2 | SwapByteOrder(line, bytes); | |
| 972 | } | ||
| 973 | } | ||
| 974 | 6 | } | |
| 975 | |||
| 976 | // ---------------------------------------------------------------------------------------------- | ||
| 977 | // TIFF | ||
| 978 | |||
| 979 | 13 | TiffImage::TiffImage(const fs::path& path) : file_(std::make_shared<const MappedFile>(path)) { | |
| 980 | 13 | const std::span<const std::byte> bytes = file_->Bytes(); | |
| 981 | 13 | const std::string name = path.string(); | |
| 982 | 13 | if (bytes.size() < 16) { | |
| 983 | ✗ | throw std::runtime_error(name + " is not a TIFF file"); | |
| 984 | } | ||
| 985 | 13 | const auto b0 = std::to_integer<char>(bytes[0]); | |
| 986 | 13 | const auto b1 = std::to_integer<char>(bytes[1]); | |
| 987 | 13 | if (b0 != b1 || (b0 != 'I' && b0 != 'M')) { | |
| 988 | 2 | throw std::runtime_error(name + " is not a TIFF file"); | |
| 989 | } | ||
| 990 | 12 | bigEndian_ = b0 == 'M'; | |
| 991 | 12 | const TiffReader r(bytes, bigEndian_); | |
| 992 | 12 | const std::uint64_t magic = r.Unsigned(2, 2); | |
| 993 | 12 | const bool big = magic == 43; | |
| 994 | 12 | if (magic != 42 && !big) { | |
| 995 | ✗ | throw std::runtime_error(name + " is not a TIFF file"); | |
| 996 | } | ||
| 997 | // The first image's directory. | ||
| 998 | 12 | const std::uint64_t ifd = big ? r.Unsigned(8, 8) : r.Unsigned(4, 4); | |
| 999 | 12 | const std::uint64_t entries = big ? r.Unsigned(ifd, 8) : r.Unsigned(ifd, 2); | |
| 1000 | 12 | const std::size_t entryBytes = big ? 20 : 12; | |
| 1001 | 10 | const std::size_t inlineBytes = big ? 8 : 4; | |
| 1002 | 12 | std::vector<std::pair<std::uint16_t, TiffField>> fields; | |
| 1003 | 166 | for (std::uint64_t i = 0; i < entries; ++i) { | |
| 1004 | 154 | const std::uint64_t at = ifd + (big ? 8 : 2) + i * entryBytes; | |
| 1005 | 154 | const auto tag = static_cast<std::uint16_t>(r.Unsigned(at, 2)); | |
| 1006 | 154 | TiffField field; | |
| 1007 | 154 | field.type = static_cast<std::uint16_t>(r.Unsigned(at + 2, 2)); | |
| 1008 | 154 | field.count = big ? r.Unsigned(at + 4, 8) : r.Unsigned(at + 4, 4); | |
| 1009 | 154 | const std::uint64_t valueAt = at + (big ? 12 : 8); | |
| 1010 | 308 | field.at = TypeBytes(field.type) * field.count <= inlineBytes | |
| 1011 | 154 | ? valueAt | |
| 1012 | 40 | : r.Unsigned(valueAt, inlineBytes); | |
| 1013 | 154 | fields.emplace_back(tag, field); | |
| 1014 | } | ||
| 1015 | 194 | const auto find = [&](std::uint16_t tag) -> const TiffField* { | |
| 1016 | 1536 | for (const auto& [t, field] : fields) { | |
| 1017 | 1492 | if (t == tag) { | |
| 1018 | 138 | return &field; | |
| 1019 | } | ||
| 1020 | } | ||
| 1021 | 44 | return nullptr; | |
| 1022 | 12 | }; | |
| 1023 | 146 | const auto values = [&](std::uint16_t tag) { | |
| 1024 | 134 | std::vector<std::uint64_t> out; | |
| 1025 | 134 | if (const TiffField* field = find(tag)) { | |
| 1026 | 130 | const std::size_t size = TypeBytes(field->type); | |
| 1027 | 130 | if (size == 0) { | |
| 1028 | ✗ | throw std::runtime_error(name + ": a TIFF field of an unknown type"); | |
| 1029 | } | ||
| 1030 | 604 | for (std::uint64_t i = 0; i < field->count; ++i) { | |
| 1031 | 474 | out.push_back(r.Unsigned(field->at + i * size, size)); | |
| 1032 | } | ||
| 1033 | } | ||
| 1034 | 134 | return out; | |
| 1035 | 12 | }; | |
| 1036 | 110 | const auto one = [&](std::uint16_t tag, std::uint64_t fallback) { | |
| 1037 | 98 | const std::vector<std::uint64_t> v = values(tag); | |
| 1038 | 98 | return v.empty() ? fallback : v.front(); | |
| 1039 | 110 | }; | |
| 1040 | 12 | width_ = static_cast<std::uint32_t>(one(kImageWidth, 0)); | |
| 1041 | 12 | height_ = static_cast<std::uint32_t>(one(kImageLength, 0)); | |
| 1042 | 12 | channels_ = static_cast<std::uint32_t>(one(kSamplesPerPixel, 1)); | |
| 1043 | 12 | const std::vector<std::uint64_t> bits = values(kBitsPerSample); | |
| 1044 | 12 | bits_ = bits.empty() ? 1u : static_cast<std::uint32_t>(bits.front()); | |
| 1045 | 12 | if (width_ == 0 || height_ == 0 || channels_ == 0) { | |
| 1046 | ✗ | throw std::runtime_error(name + ": an empty TIFF image"); | |
| 1047 | } | ||
| 1048 | 36 | if (std::any_of(bits.begin(), bits.end(), [&](std::uint64_t b) { return b != bits_; })) { | |
| 1049 | ✗ | throw std::runtime_error(name + ": TIFF samples of mixed sizes"); | |
| 1050 | } | ||
| 1051 | // Planes stored apart: each channel's strips or tiles in turn. | ||
| 1052 | 12 | planes_ = one(kPlanarConfiguration, 1) == 2 ? channels_ : 1u; | |
| 1053 | 12 | const std::uint64_t format = one(kSampleFormat, 1); | |
| 1054 | 12 | if (format == 3 && bits_ == 32) { | |
| 1055 | 4 | type_ = GridSampleType::Float32; | |
| 1056 | 8 | } else if (format == 1 && bits_ == 8) { | |
| 1057 | 6 | type_ = GridSampleType::UInt8; | |
| 1058 | ✗ | } else if (format == 1 && bits_ == 16) { | |
| 1059 | ✗ | type_ = GridSampleType::UInt16; | |
| 1060 | 2 | } else if (format == 2 && bits_ == 16) { | |
| 1061 | 2 | type_ = GridSampleType::Int16; | |
| 1062 | ✗ | } else if (format == 1 && bits_ == 32) { | |
| 1063 | ✗ | type_ = GridSampleType::UInt32; | |
| 1064 | ✗ | } else if (format == 2 && bits_ == 32) { | |
| 1065 | ✗ | type_ = GridSampleType::Int32; | |
| 1066 | } else { | ||
| 1067 | ✗ | throw std::runtime_error(name + ": TIFF samples of " + std::to_string(bits_) + | |
| 1068 | ✗ | " bits, format " + std::to_string(format) + ", are not read"); | |
| 1069 | } | ||
| 1070 | 12 | compression_ = static_cast<std::uint32_t>(one(kCompression, 1)); | |
| 1071 | 12 | if (compression_ != 1 && compression_ != 5 && compression_ != 8 && compression_ != 32946) { | |
| 1072 | ✗ | throw std::runtime_error(name + ": TIFF compression " + std::to_string(compression_) + | |
| 1073 | ✗ | " is not read"); | |
| 1074 | } | ||
| 1075 | 12 | predictor_ = static_cast<std::uint32_t>(one(kPredictor, 1)); | |
| 1076 | 12 | if (predictor_ > 3 || (predictor_ == 3 && type_ != GridSampleType::Float32)) { | |
| 1077 | ✗ | throw std::runtime_error(name + ": TIFF predictor " + std::to_string(predictor_) + | |
| 1078 | ✗ | " is not read"); | |
| 1079 | } | ||
| 1080 | 12 | tiled_ = find(kTileOffsets) != nullptr; | |
| 1081 | 12 | if (tiled_) { | |
| 1082 | 2 | chunkWidth_ = static_cast<std::uint32_t>(one(kTileWidth, 0)); | |
| 1083 | 2 | chunkHeight_ = static_cast<std::uint32_t>(one(kTileLength, 0)); | |
| 1084 | 2 | offsets_ = values(kTileOffsets); | |
| 1085 | 2 | byteCounts_ = values(kTileByteCounts); | |
| 1086 | } else { | ||
| 1087 | 10 | chunkWidth_ = width_; | |
| 1088 | 20 | chunkHeight_ = static_cast<std::uint32_t>( | |
| 1089 | 10 | std::min<std::uint64_t>(one(kRowsPerStrip, height_), height_)); | |
| 1090 | 10 | offsets_ = values(kStripOffsets); | |
| 1091 | 10 | byteCounts_ = values(kStripByteCounts); | |
| 1092 | } | ||
| 1093 | 12 | if (chunkWidth_ == 0 || chunkHeight_ == 0) { | |
| 1094 | ✗ | throw std::runtime_error(name + ": TIFF strips or tiles of no size"); | |
| 1095 | } | ||
| 1096 | 12 | chunksAcross_ = (width_ + chunkWidth_ - 1) / chunkWidth_; | |
| 1097 | 12 | chunksPerPlane_ = std::size_t{chunksAcross_} * ((height_ + chunkHeight_ - 1) / chunkHeight_); | |
| 1098 | 12 | const std::size_t chunks = chunksPerPlane_ * planes_; | |
| 1099 | 12 | cachedChunks_.assign(planes_, ~std::size_t{0}); | |
| 1100 | 12 | cached_.resize(planes_); | |
| 1101 | 12 | if (offsets_.size() < chunks || byteCounts_.size() < chunks) { | |
| 1102 | ✗ | throw std::runtime_error(name + ": TIFF strips or tiles missing"); | |
| 1103 | } | ||
| 1104 | // GeoTIFF's grid, and GDAL's no-data value (ASCII). | ||
| 1105 | 12 | if (const TiffField* field = find(kModelPixelScale); | |
| 1106 | 2 | field != nullptr && field->count >= 3 && field->type == 12) { | |
| 1107 | 2 | pixelScale_ = std::array<double, 3>{r.Double(field->at), r.Double(field->at + 8), | |
| 1108 | 4 | r.Double(field->at + 16)}; | |
| 1109 | } | ||
| 1110 | 12 | if (const TiffField* field = find(kModelTiepoint); | |
| 1111 | 2 | field != nullptr && field->count >= 6 && field->type == 12) { | |
| 1112 | 2 | std::array<double, 6> tie{}; | |
| 1113 | 14 | for (std::size_t i = 0; i < 6; ++i) { | |
| 1114 | 12 | tie[i] = r.Double(field->at + 8 * i); | |
| 1115 | } | ||
| 1116 | 4 | tiepoint_ = tie; | |
| 1117 | } | ||
| 1118 | 12 | if (const TiffField* field = find(kGdalNoData); | |
| 1119 | 2 | field != nullptr && field->type == 2 && field->count > 0) { | |
| 1120 | 2 | std::string text; | |
| 1121 | 14 | for (std::uint64_t i = 0; i < field->count; ++i) { | |
| 1122 | 14 | const auto c = static_cast<char>(r.Unsigned(field->at + i, 1)); | |
| 1123 | 14 | if (c == '\0') { | |
| 1124 | break; | ||
| 1125 | } | ||
| 1126 | 12 | text.push_back(c); | |
| 1127 | } | ||
| 1128 | 2 | try { | |
| 1129 | 2 | noData_ = std::stod(text); | |
| 1130 | ✗ | } catch (const std::exception&) { | |
| 1131 | // Not a number ("nan" parses; anything else is ignored). | ||
| 1132 | ✗ | } | |
| 1133 | 2 | } | |
| 1134 | 15 | } | |
| 1135 | |||
| 1136 | 405 | const std::vector<std::byte>& TiffImage::Chunk(std::size_t chunk) const { | |
| 1137 | 405 | const std::size_t plane = chunk / chunksPerPlane_; | |
| 1138 | 405 | if (chunk == cachedChunks_[plane]) { | |
| 1139 | 171 | return cached_[plane]; | |
| 1140 | } | ||
| 1141 | 234 | ERNIE_ZONE("TIFF chunk"); | |
| 1142 | 234 | const std::uint32_t channels = ChunkChannels(); | |
| 1143 | 234 | const std::uint32_t row0 = | |
| 1144 | 234 | static_cast<std::uint32_t>(chunk % chunksPerPlane_ / chunksAcross_) * chunkHeight_; | |
| 1145 | // A strip at the image's foot holds only the rows left; tiles are always whole. | ||
| 1146 | 243 | const std::uint32_t rows = tiled_ ? chunkHeight_ : std::min(chunkHeight_, height_ - row0); | |
| 1147 | 234 | const std::uint32_t sampleBytes = SampleBytes(); | |
| 1148 | 234 | const std::size_t rowBytes = std::size_t{chunkWidth_} * channels * sampleBytes; | |
| 1149 | 234 | const std::size_t expected = rowBytes * rows; | |
| 1150 | 234 | const std::span<const std::byte> file = file_->Bytes(); | |
| 1151 | 234 | const std::uint64_t offset = offsets_[chunk]; | |
| 1152 | 234 | const std::uint64_t count = byteCounts_[chunk]; | |
| 1153 | 234 | if (offset > file.size() || count > file.size() - offset) { | |
| 1154 | ✗ | throw std::runtime_error("TIFF: a strip or tile past the file's end"); | |
| 1155 | } | ||
| 1156 | 234 | const std::span<const std::byte> input = file.subspan(offset, count); | |
| 1157 | 234 | std::vector<std::byte> data; | |
| 1158 | 234 | switch (compression_) { | |
| 1159 | 26 | case 5: | |
| 1160 | 52 | data = DecodeLzw(input, expected); | |
| 1161 | 26 | break; | |
| 1162 | 77 | case 8: | |
| 1163 | 77 | case 32946: | |
| 1164 | 154 | data = DecodeDeflate(input, expected); | |
| 1165 | 77 | break; | |
| 1166 | 131 | default: | |
| 1167 | 131 | data.assign(input.begin(), input.end()); | |
| 1168 | break; | ||
| 1169 | } | ||
| 1170 | 234 | if (data.size() < expected) { | |
| 1171 | ✗ | throw std::runtime_error("TIFF: a strip or tile cut short"); | |
| 1172 | } | ||
| 1173 | 234 | data.resize(expected); | |
| 1174 | // Predictor 3 reorders the bytes itself; otherwise samples are put in host order before | ||
| 1175 | // predictor 2 is undone. | ||
| 1176 | 234 | if (predictor_ == 3) { | |
| 1177 | 5 | UndoFloatingPointPredictor(data, rowBytes, rows, channels, sampleBytes); | |
| 1178 | 229 | } else if (NeedsByteSwap(bigEndian_, sampleBytes)) { | |
| 1179 | 1 | SwapByteOrder(data, sampleBytes); | |
| 1180 | } | ||
| 1181 | 234 | if (predictor_ == 2) { | |
| 1182 | 26 | UndoHorizontalPredictor(data, rowBytes, rows, channels, sampleBytes); | |
| 1183 | } | ||
| 1184 | 234 | cached_[plane] = std::move(data); | |
| 1185 | 234 | cachedChunks_[plane] = chunk; | |
| 1186 | 234 | return cached_[plane]; | |
| 1187 | 234 | } | |
| 1188 | |||
| 1189 | 15 | void TiffImage::Read(std::uint32_t x, std::uint32_t y, std::uint32_t w, std::uint32_t h, | |
| 1190 | std::span<std::byte> out) const { | ||
| 1191 | 15 | CheckRead(x, y, w, h, out); | |
| 1192 | 13 | const std::size_t pixelBytes = std::size_t{channels_} * SampleBytes(); | |
| 1193 | // Each plane's samples go to their channel of the pixels. | ||
| 1194 | 13 | const std::size_t sampleBytes = SampleBytes(); | |
| 1195 | 13 | const std::size_t chunkPixel = ChunkChannels() * sampleBytes; | |
| 1196 | 34 | for (std::uint32_t plane = 0; plane < planes_; ++plane) { | |
| 1197 | 380 | for (std::uint32_t row = y; row < y + h; ++row) { | |
| 1198 | 359 | const std::uint32_t chunkRow = row / chunkHeight_; | |
| 1199 | 764 | for (std::uint32_t column = x; column < x + w;) { | |
| 1200 | 405 | const std::uint32_t chunkColumn = column / chunkWidth_; | |
| 1201 | 810 | const std::vector<std::byte>& chunk = Chunk( | |
| 1202 | 405 | plane * chunksPerPlane_ + std::size_t{chunkRow} * chunksAcross_ + chunkColumn); | |
| 1203 | 405 | const std::uint32_t left = chunkColumn * chunkWidth_; | |
| 1204 | 405 | const std::uint32_t span = std::min(x + w, left + chunkWidth_) - column; | |
| 1205 | 405 | std::byte* to = out.data() + (std::size_t{row - y} * w + (column - x)) * pixelBytes; | |
| 1206 | 405 | const std::byte* from = | |
| 1207 | 405 | chunk.data() + | |
| 1208 | 405 | (std::size_t{row - chunkRow * chunkHeight_} * chunkWidth_ + (column - left)) * | |
| 1209 | 405 | chunkPixel; | |
| 1210 | 405 | if (planes_ == 1) { | |
| 1211 | 189 | std::memcpy(to, from, span * pixelBytes); | |
| 1212 | } else { | ||
| 1213 | 5712 | for (std::uint32_t i = 0; i < span; ++i) { | |
| 1214 | 5496 | std::memcpy(to + i * pixelBytes + plane * sampleBytes, | |
| 1215 | 5496 | from + i * sampleBytes, sampleBytes); | |
| 1216 | } | ||
| 1217 | } | ||
| 1218 | column += span; | ||
| 1219 | } | ||
| 1220 | } | ||
| 1221 | } | ||
| 1222 | 13 | } | |
| 1223 | |||
| 1224 | 4 | std::optional<LatLongGrid> TiffImage::Grid() const { | |
| 1225 | 4 | if (!pixelScale_ || !tiepoint_) { | |
| 1226 | 1 | return std::nullopt; | |
| 1227 | } | ||
| 1228 | 3 | const std::array<double, 3>& scale = *pixelScale_; | |
| 1229 | 3 | const std::array<double, 6>& tie = *tiepoint_; | |
| 1230 | 3 | const LatLongGrid grid{tie[3] - tie[0] * scale[0], tie[4] + tie[1] * scale[1], scale[0], | |
| 1231 | 3 | scale[1]}; | |
| 1232 | // Degrees: within a turn of longitude and between the poles (metres of a projected grid are far | ||
| 1233 | // larger). | ||
| 1234 | 3 | constexpr double kSlack = 1e-6; | |
| 1235 | 3 | if (!(grid.degreesAcross > 0.0) || !(grid.degreesDown > 0.0) || | |
| 1236 | 3 | grid.degreesAcross * width_ > 360.0 + kSlack || | |
| 1237 | 6 | grid.degreesDown * height_ > 180.0 + kSlack || std::abs(grid.west) > 360.0 || | |
| 1238 | 3 | std::abs(grid.north) > 90.0 + kSlack) { | |
| 1239 | ✗ | return std::nullopt; | |
| 1240 | } | ||
| 1241 | 3 | return grid; | |
| 1242 | } | ||
| 1243 | |||
| 1244 | // ---------------------------------------------------------------------------------------------- | ||
| 1245 | // Planetary DEMs | ||
| 1246 | |||
| 1247 | 1841 | double CylindricalDem::HeightAt(double latitudeDeg, double longitudeDeg) const { | |
| 1248 | 1841 | const double east = std::fmod(std::fmod(longitudeDeg - westLongitudeDeg, 360.0) + 360.0, 360.0); | |
| 1249 | 1841 | return SampleDem(heights, lines, samples, | |
| 1250 | 1841 | (maxLatitudeDeg - latitudeDeg) * pixelsPerDegree - 0.5, | |
| 1251 | 1841 | east * pixelsPerDegree * std::cos(glm::radians(standardParallelDeg)) - 0.5); | |
| 1252 | } | ||
| 1253 | |||
| 1254 | 5 | CylindricalDem ReadCylindricalDem(const fs::path& labelPath) { | |
| 1255 | 5 | const PdsImage image(labelPath); | |
| 1256 | 4 | const std::string projection = DemProjection(image, {"SIMPLE CYLINDRICAL", "EQUIRECTANGULAR"}); | |
| 1257 | 3 | const Pds3Label& label = image.Keywords(); | |
| 1258 | 3 | DemGrid grid = ReadDemGrid(image); | |
| 1259 | 3 | CylindricalDem dem; | |
| 1260 | 3 | dem.lines = grid.lines; | |
| 1261 | 3 | dem.samples = grid.samples; | |
| 1262 | 3 | dem.referenceRadius = grid.referenceRadius; | |
| 1263 | 3 | dem.heights = std::move(grid.heights); | |
| 1264 | 3 | dem.pixelsPerDegree = label.Number("IMAGE_MAP_PROJECTION.MAP_RESOLUTION"); | |
| 1265 | 3 | dem.maxLatitudeDeg = label.Number("IMAGE_MAP_PROJECTION.MAXIMUM_LATITUDE"); | |
| 1266 | 3 | dem.westLongitudeDeg = label.Number("IMAGE_MAP_PROJECTION.WESTERNMOST_LONGITUDE"); | |
| 1267 | 6 | dem.standardParallelDeg = | |
| 1268 | 2 | projection == "EQUIRECTANGULAR" && label.Get("IMAGE_MAP_PROJECTION.CENTER_LATITUDE") | |
| 1269 | 4 | ? label.Number("IMAGE_MAP_PROJECTION.CENTER_LATITUDE") | |
| 1270 | : 0.0; | ||
| 1271 | 6 | return dem; | |
| 1272 | 4 | } | |
| 1273 | |||
| 1274 | 40 | glm::dvec2 PolarStereographicDem::ToGrid(double latitudeDeg, double longitudeDeg) const { | |
| 1275 | 40 | const double lat = glm::radians(latitudeDeg); | |
| 1276 | 40 | const double lon = glm::radians(longitudeDeg - centerLongitudeDeg); | |
| 1277 | 40 | const bool south = centerLatitudeDeg < 0.0; | |
| 1278 | // Snyder 1987, eq. 21-33 and 21-34 for the sphere (k0 = 1, true at the pole): | ||
| 1279 | // rho = 2 R tan(pi/4 -+ phi/2). | ||
| 1280 | 80 | const double rho = | |
| 1281 | 40 | 2.0 * referenceRadius * std::tan(0.25 * std::numbers::pi + (south ? 0.5 : -0.5) * lat); | |
| 1282 | 40 | const double x = rho * std::sin(lon); | |
| 1283 | 40 | const double y = (south ? 1.0 : -1.0) * rho * std::cos(lon); | |
| 1284 | 40 | return {lineProjectionOffset - y / mapScale, sampleProjectionOffset + x / mapScale}; | |
| 1285 | } | ||
| 1286 | |||
| 1287 | 39 | double PolarStereographicDem::HeightAt(double latitudeDeg, double longitudeDeg) const { | |
| 1288 | 39 | const glm::dvec2 at = ToGrid(latitudeDeg, longitudeDeg); | |
| 1289 | 39 | return SampleDem(heights, lines, samples, at.x, at.y); | |
| 1290 | } | ||
| 1291 | |||
| 1292 | 2 | PolarStereographicDem ReadPolarStereographicDem(const fs::path& labelPath) { | |
| 1293 | 2 | const PdsImage image(labelPath); | |
| 1294 | 2 | DemProjection(image, {"POLAR STEREOGRAPHIC"}); | |
| 1295 | 1 | const Pds3Label& label = image.Keywords(); | |
| 1296 | 1 | DemGrid grid = ReadDemGrid(image); | |
| 1297 | 1 | PolarStereographicDem dem; | |
| 1298 | 1 | dem.lines = grid.lines; | |
| 1299 | 1 | dem.samples = grid.samples; | |
| 1300 | 1 | dem.referenceRadius = grid.referenceRadius; | |
| 1301 | 1 | dem.heights = std::move(grid.heights); | |
| 1302 | 1 | dem.centerLatitudeDeg = label.Number("IMAGE_MAP_PROJECTION.CENTER_LATITUDE"); | |
| 1303 | 1 | if (std::abs(std::abs(dem.centerLatitudeDeg) - 90.0) > 1e-9) { | |
| 1304 | ✗ | throw std::invalid_argument( | |
| 1305 | ✗ | "The DEM's polar stereographic projection is not centred on a pole"); | |
| 1306 | } | ||
| 1307 | 1 | dem.centerLongitudeDeg = label.Get("IMAGE_MAP_PROJECTION.CENTER_LONGITUDE") | |
| 1308 | 1 | ? label.Number("IMAGE_MAP_PROJECTION.CENTER_LONGITUDE") | |
| 1309 | : 0.0; | ||
| 1310 | // The scale's units are read from the value as written, since `Number` drops them. | ||
| 1311 | 1 | dem.mapScale = label.Number("IMAGE_MAP_PROJECTION.MAP_SCALE"); | |
| 1312 | 1 | if (ToUpper(label.Find("MAP_SCALE").value_or("")).find("<KM") != std::string::npos) { | |
| 1313 | ✗ | dem.mapScale *= 1e3; | |
| 1314 | } | ||
| 1315 | 1 | dem.lineProjectionOffset = label.Number("IMAGE_MAP_PROJECTION.LINE_PROJECTION_OFFSET"); | |
| 1316 | 1 | dem.sampleProjectionOffset = label.Number("IMAGE_MAP_PROJECTION.SAMPLE_PROJECTION_OFFSET"); | |
| 1317 | 2 | return dem; | |
| 1318 | 2 | } | |
| 1319 | |||
| 1320 | 4 | glm::dvec3 BodyDirection(double latitudeDeg, double longitudeDeg) { | |
| 1321 | 4 | const double lat = glm::radians(latitudeDeg); | |
| 1322 | 4 | const double lon = glm::radians(longitudeDeg); | |
| 1323 | 4 | return {std::cos(lat) * std::cos(lon), std::cos(lat) * std::sin(lon), std::sin(lat)}; | |
| 1324 | } | ||
| 1325 | |||
| 1326 | 2 | LocalTerrain MakeLocalTerrain(const CylindricalDem& dem, double latitudeDeg, double longitudeDeg, | |
| 1327 | std::uint32_t cells, double cellSize) { | ||
| 1328 | 2 | return LocalTerrainAround( | |
| 1329 | 1764 | dem.referenceRadius, [&dem](double lat, double lon) { return dem.HeightAt(lat, lon); }, | |
| 1330 | 2 | latitudeDeg, longitudeDeg, cells, cellSize); | |
| 1331 | } | ||
| 1332 | |||
| 1333 | 1 | LocalTerrain MakeLocalTerrain(const PolarStereographicDem& dem, double latitudeDeg, | |
| 1334 | double longitudeDeg, std::uint32_t cells, double cellSize) { | ||
| 1335 | 1 | return LocalTerrainAround( | |
| 1336 | 36 | dem.referenceRadius, [&dem](double lat, double lon) { return dem.HeightAt(lat, lon); }, | |
| 1337 | 1 | latitudeDeg, longitudeDeg, cells, cellSize); | |
| 1338 | } | ||
| 1339 | |||
| 1340 | 1 | LocalTerrain MakeLocalTerrain(const CylindricalDem& dem, const CylindricalDem& fallback, | |
| 1341 | double latitudeDeg, double longitudeDeg, std::uint32_t cells, | ||
| 1342 | double cellSize) { | ||
| 1343 | 1 | return LocalTerrainAround( | |
| 1344 | 1 | dem.referenceRadius, | |
| 1345 | 38 | [&](double lat, double lon) { | |
| 1346 | 36 | const double h = dem.HeightAt(lat, lon); | |
| 1347 | 36 | return std::isnan(h) ? fallback.HeightAt(lat, lon) : h; | |
| 1348 | }, | ||
| 1349 | 1 | latitudeDeg, longitudeDeg, cells, cellSize); | |
| 1350 | } | ||
| 1351 | |||
| 1352 | } // namespace ernie | ||
| 1353 |