io/raster.h
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // SPDX-License-Identifier: BSD-3-Clause | ||
| 2 | /// @file raster.h | ||
| 3 | /// @brief Rasters read a rectangle at a time: TIFF/GeoTIFF, PDS3 and decoded images, mosaics of | ||
| 4 | /// them, planetary DEMs and the local heightfields drawn of them. | ||
| 5 | #pragma once | ||
| 6 | |||
| 7 | #include "ernie/export.h" | ||
| 8 | |||
| 9 | #include "ernie/util/file.h" | ||
| 10 | #include "ernie/util/tensor.h" | ||
| 11 | |||
| 12 | #include "ernie/util/glm.h" | ||
| 13 | |||
| 14 | #include <array> | ||
| 15 | #include <cmath> | ||
| 16 | #include <cstddef> | ||
| 17 | #include <cstdint> | ||
| 18 | #include <filesystem> | ||
| 19 | #include <memory> | ||
| 20 | #include <optional> | ||
| 21 | #include <span> | ||
| 22 | #include <string> | ||
| 23 | #include <string_view> | ||
| 24 | #include <vector> | ||
| 25 | |||
| 26 | namespace ernie { | ||
| 27 | |||
| 28 | // ---------------------------------------------------------------------------------------------- | ||
| 29 | // Rasters | ||
| 30 | |||
| 31 | /// @brief How a raster's samples are stored. | ||
| 32 | enum class GridSampleType { | ||
| 33 | UInt8, | ||
| 34 | UInt16, | ||
| 35 | Int16, | ||
| 36 | UInt32, | ||
| 37 | Int32, | ||
| 38 | Float32, | ||
| 39 | }; | ||
| 40 | |||
| 41 | /// @brief The size of a sample type. | ||
| 42 | /// @param type The type. | ||
| 43 | /// @return Bytes per sample. | ||
| 44 | ERNIE_API std::uint32_t ByteSize(GridSampleType type); | ||
| 45 | |||
| 46 | /// @brief Decodes one sample in host byte order. | ||
| 47 | /// @param sample The sample's bytes. | ||
| 48 | /// @param type How it is stored. | ||
| 49 | /// @return Its value. | ||
| 50 | ERNIE_API double DecodeSample(const std::byte* sample, GridSampleType type); | ||
| 51 | |||
| 52 | /// @brief Where an equirectangular raster lies. | ||
| 53 | struct LatLongGrid { | ||
| 54 | double west = -180.0; ///< Longitude of the west edge, in degrees east. | ||
| 55 | double north = 90.0; ///< Latitude of the north edge, in degrees north. | ||
| 56 | double degreesAcross = 0.0; ///< Degrees per pixel across; positive. | ||
| 57 | double degreesDown = 0.0; ///< Degrees per pixel down; positive. | ||
| 58 | }; | ||
| 59 | |||
| 60 | /// @brief A rectangle of a raster's pixels. | ||
| 61 | struct GridRect { | ||
| 62 | std::uint32_t x = 0; ///< Left column. | ||
| 63 | std::uint32_t y = 0; ///< Top row. | ||
| 64 | std::uint32_t width = 0; ///< Width, in pixels. | ||
| 65 | std::uint32_t height = 0; ///< Height, in pixels. | ||
| 66 | }; | ||
| 67 | |||
| 68 | /// @brief A raster image read a rectangle at a time. | ||
| 69 | /// Lets images larger than memory be baked, resampled and streamed in parts. Read by one thread at | ||
| 70 | /// a time. | ||
| 71 | 56 | class ERNIE_API GridSource { | |
| 72 | public: | ||
| 73 | 5 | virtual ~GridSource() = default; | |
| 74 | |||
| 75 | /// @return Width, in pixels. | ||
| 76 | virtual std::uint32_t Width() const = 0; | ||
| 77 | /// @return Height, in pixels. | ||
| 78 | virtual std::uint32_t Height() const = 0; | ||
| 79 | /// @return Samples per pixel (1 grey, 3 RGB, 4 RGBA). | ||
| 80 | virtual std::uint32_t Channels() const = 0; | ||
| 81 | /// @return How samples are stored. | ||
| 82 | virtual GridSampleType Type() const = 0; | ||
| 83 | |||
| 84 | /// @return Bytes per sample. | ||
| 85 | 428 | std::uint32_t SampleBytes() const { | |
| 86 | 428 | return ByteSize(Type()); | |
| 87 | } | ||
| 88 | |||
| 89 | /// @brief Reads a rectangle of samples in host byte order. | ||
| 90 | /// @param x The left column. | ||
| 91 | /// @param y The top row. | ||
| 92 | /// @param w The width, in pixels. | ||
| 93 | /// @param h The height, in pixels. | ||
| 94 | /// @param out w * h * Channels() * SampleBytes() bytes, rows top to bottom, interleaved. | ||
| 95 | /// @throws std::out_of_range If the rectangle passes the image or `out` is too small. | ||
| 96 | /// @throws std::runtime_error If the data is damaged. | ||
| 97 | virtual void Read(std::uint32_t x, std::uint32_t y, std::uint32_t w, std::uint32_t h, | ||
| 98 | std::span<std::byte> out) const = 0; | ||
| 99 | |||
| 100 | /// @return The scale in sample * Scale() + Offset(), the physical value (PDS SCALING_FACTOR; | ||
| 101 | /// usually 1). | ||
| 102 | ✗ | virtual double Scale() const { | |
| 103 | ✗ | return 1.0; | |
| 104 | } | ||
| 105 | /// @return The offset to physical values (usually 0). | ||
| 106 | ✗ | virtual double Offset() const { | |
| 107 | ✗ | return 0.0; | |
| 108 | } | ||
| 109 | /// @return The sample marking no data, if any. | ||
| 110 | ✗ | virtual std::optional<double> NoData() const { | |
| 111 | ✗ | return std::nullopt; | |
| 112 | } | ||
| 113 | /// @return Its latitude-longitude placement, if the file says. | ||
| 114 | ✗ | virtual std::optional<LatLongGrid> Grid() const { | |
| 115 | ✗ | return std::nullopt; | |
| 116 | } | ||
| 117 | /// @return Rectangles best read whole one after another, such as a mosaic's tiles; the whole | ||
| 118 | /// image by default. | ||
| 119 | 12 | virtual std::vector<GridRect> Parts() const { | |
| 120 | 12 | return {GridRect{0, 0, Width(), Height()}}; | |
| 121 | } | ||
| 122 | /// @brief Frees memory held for reading, such as a decoded image; it is reloaded if needed. | ||
| 123 | 1 | virtual void Release() const {} | |
| 124 | |||
| 125 | /// @return Every sample, as `Read` of the whole image. | ||
| 126 | std::vector<std::byte> ReadAll() const; | ||
| 127 | |||
| 128 | protected: | ||
| 129 | /// @brief Checks a `Read`'s rectangle and room. | ||
| 130 | /// @throws std::out_of_range If the rectangle passes the image or `out` is too small. | ||
| 131 | void CheckRead(std::uint32_t x, std::uint32_t y, std::uint32_t w, std::uint32_t h, | ||
| 132 | std::span<const std::byte> out) const; | ||
| 133 | }; | ||
| 134 | |||
| 135 | /// @brief A grid of like tiles read as one raster. | ||
| 136 | /// For global mosaics too large for one file, such as Blue Marble's eight 21600 x 21600 tiles. | ||
| 137 | /// Parts are tiles, decoded one at a time; the grid is the first tile's, spread over the rest. | ||
| 138 | class ERNIE_API MosaicRaster : public GridSource { | ||
| 139 | public: | ||
| 140 | /// @brief Opens the tiles. | ||
| 141 | /// @param tiles The tile files, row by row, top to bottom. | ||
| 142 | /// @param columns Tiles per row. | ||
| 143 | /// @throws std::invalid_argument If the tiles do not fill whole rows or differ; else as | ||
| 144 | /// `OpenRaster`. | ||
| 145 | MosaicRaster(const std::vector<std::filesystem::path>& tiles, std::uint32_t columns); | ||
| 146 | MosaicRaster(const MosaicRaster&) = delete; | ||
| 147 | MosaicRaster& operator=(const MosaicRaster&) = delete; | ||
| 148 | |||
| 149 | 42 | std::uint32_t Width() const override { | |
| 150 | 42 | return columns_ * tileWidth_; | |
| 151 | } | ||
| 152 | 42 | std::uint32_t Height() const override { | |
| 153 | 42 | return rows_ * tileHeight_; | |
| 154 | } | ||
| 155 | 56 | std::uint32_t Channels() const override { | |
| 156 | 56 | return tiles_.front()->Channels(); | |
| 157 | } | ||
| 158 | 60 | GridSampleType Type() const override { | |
| 159 | 60 | return tiles_.front()->Type(); | |
| 160 | } | ||
| 161 | ✗ | double Scale() const override { | |
| 162 | ✗ | return tiles_.front()->Scale(); | |
| 163 | } | ||
| 164 | ✗ | double Offset() const override { | |
| 165 | ✗ | return tiles_.front()->Offset(); | |
| 166 | } | ||
| 167 | ✗ | std::optional<double> NoData() const override { | |
| 168 | ✗ | return tiles_.front()->NoData(); | |
| 169 | } | ||
| 170 | ✗ | std::optional<LatLongGrid> Grid() const override { | |
| 171 | ✗ | return tiles_.front()->Grid(); | |
| 172 | } | ||
| 173 | std::vector<GridRect> Parts() const override; | ||
| 174 | void Release() const override; | ||
| 175 | void Read(std::uint32_t x, std::uint32_t y, std::uint32_t w, std::uint32_t h, | ||
| 176 | std::span<std::byte> out) const override; | ||
| 177 | |||
| 178 | private: | ||
| 179 | std::vector<std::unique_ptr<GridSource>> tiles_; | ||
| 180 | std::uint32_t columns_ = 0; | ||
| 181 | std::uint32_t rows_ = 0; | ||
| 182 | std::uint32_t tileWidth_ = 0; | ||
| 183 | std::uint32_t tileHeight_ = 0; | ||
| 184 | /// The tile read last, released as another is read. | ||
| 185 | mutable std::size_t current_ = 0; | ||
| 186 | }; | ||
| 187 | |||
| 188 | /// @brief Opens the mosaic a `.mosaic` file describes. | ||
| 189 | /// Each line is a row of tile paths relative to the file, top to bottom, separated by spaces; `#` | ||
| 190 | /// starts a comment. | ||
| 191 | /// @param path The `.mosaic` file. | ||
| 192 | /// @return The mosaic. | ||
| 193 | /// @throws std::runtime_error If it cannot be read. | ||
| 194 | /// @throws std::invalid_argument For rows of different lengths, or as `MosaicRaster`. | ||
| 195 | ERNIE_API std::unique_ptr<MosaicRaster> OpenMosaic(const std::filesystem::path& path); | ||
| 196 | |||
| 197 | /// @brief Opens a raster by extension. | ||
| 198 | /// TIFF/GeoTIFF and PDS3 (.img with .lbl, or the .lbl) read in place; 8/16-bit PNG and JPEG decode | ||
| 199 | /// whole on first read; `.mosaic` as `OpenMosaic`. | ||
| 200 | /// @param path The file. | ||
| 201 | /// @return The raster. | ||
| 202 | /// @throws std::runtime_error If it cannot be read. | ||
| 203 | ERNIE_API std::unique_ptr<GridSource> OpenRaster(const std::filesystem::path& path); | ||
| 204 | |||
| 205 | // ---------------------------------------------------------------------------------------------- | ||
| 206 | // Grid views | ||
| 207 | |||
| 208 | /// @brief A non-owning, north-up grid of `float` samples. | ||
| 209 | /// Ernie takes heightfields and albedo maps in this form. Row 0 is the northern edge; the world | ||
| 210 | /// frame is x east, y north, z up, in camera units. | ||
| 211 | struct GridView { | ||
| 212 | const float* data = nullptr; ///< The first sample. | ||
| 213 | std::size_t width = 0; ///< Columns. | ||
| 214 | std::size_t height = 0; ///< Rows. | ||
| 215 | std::size_t rowStride = 0; ///< Elements from one row to the next; 0 means `width`. | ||
| 216 | double originX = 0.0; ///< West edge of the first cell, in world units. | ||
| 217 | double originY = 0.0; ///< North edge of the first cell, in world units. | ||
| 218 | double cellSize = 1.0; ///< Side of a cell in world units; positive. | ||
| 219 | std::optional<float> noData; ///< Value marking cells without data, if any. | ||
| 220 | |||
| 221 | /// @return The row stride in elements. | ||
| 222 | ✗ | std::size_t Stride() const { | |
| 223 | ✗ | return rowStride == 0 ? width : rowStride; | |
| 224 | } | ||
| 225 | /// @brief Reads one cell. | ||
| 226 | /// @param col The column. | ||
| 227 | /// @param row The row, 0 at the north. | ||
| 228 | /// @return The sample. | ||
| 229 | float operator()(std::size_t col, std::size_t row) const { | ||
| 230 | return data[row * Stride() + col]; | ||
| 231 | } | ||
| 232 | /// @brief Tests a sample for data. | ||
| 233 | /// @param value The sample. | ||
| 234 | /// @return Whether it is finite and not the nodata value. | ||
| 235 | bool IsValid(float value) const { | ||
| 236 | return std::isfinite(value) && !(noData && value == *noData); | ||
| 237 | } | ||
| 238 | }; | ||
| 239 | |||
| 240 | /// @brief Checks that a grid view is usable. | ||
| 241 | /// @param view The view to check. | ||
| 242 | /// @param what Its name, for the error message. | ||
| 243 | /// @param minSide The fewest cells allowed along each side. | ||
| 244 | /// @throws std::invalid_argument If it has no data, too few cells, a row stride below its width, a | ||
| 245 | /// non-finite origin, or a cell size that is not positive and finite. | ||
| 246 | ERNIE_API void CheckGridView(const GridView& view, const std::string& what, std::size_t minSide); | ||
| 247 | |||
| 248 | // ---------------------------------------------------------------------------------------------- | ||
| 249 | // PDS3 | ||
| 250 | |||
| 251 | /// @brief A PDS3 label's keywords, in the order the label gives them. | ||
| 252 | 20 | struct ERNIE_API Pds3Label { | |
| 253 | /// @brief One `KEY = value` line. | ||
| 254 | struct Keyword { | ||
| 255 | /// The OBJECTs holding it, outermost first, each followed by a dot (`IMAGE.`); empty at the | ||
| 256 | /// top level. | ||
| 257 | std::string object; | ||
| 258 | /// Its key, as the label writes it. | ||
| 259 | std::string key; | ||
| 260 | /// Its value as written: comments dropped and continued lines joined, units (`<km>`) and | ||
| 261 | /// quotes kept. | ||
| 262 | std::string value; | ||
| 263 | }; | ||
| 264 | |||
| 265 | /// The keywords; OBJECT and END_OBJECT lines are not among them. | ||
| 266 | std::vector<Keyword> keywords; | ||
| 267 | |||
| 268 | /// @brief A keyword by its object path and key. | ||
| 269 | /// @param name Its objects and key, such as `IMAGE.LINES`; case matters. | ||
| 270 | /// @return Its value without units or quotes, the label's last of that name, or none. | ||
| 271 | std::optional<std::string> Get(std::string_view name) const; | ||
| 272 | /// @brief A keyword by its object path and key, as a number. | ||
| 273 | /// @param name Its objects and key, such as `IMAGE.LINES`. | ||
| 274 | /// @return Its value. | ||
| 275 | /// @throws std::invalid_argument If it is missing or not a number. | ||
| 276 | double Number(std::string_view name) const; | ||
| 277 | /// @brief A keyword by its object path and key. | ||
| 278 | /// @param name Its objects and key, such as `IMAGE.LINES`. | ||
| 279 | /// @return Its value without units or quotes. | ||
| 280 | /// @throws std::invalid_argument If it is missing. | ||
| 281 | std::string Text(std::string_view name) const; | ||
| 282 | /// @brief A keyword by its key alone, in any object. | ||
| 283 | /// @param key Its key, in any case. | ||
| 284 | /// @return Its value as written, the label's first with that key, or none. | ||
| 285 | std::optional<std::string> Find(std::string_view key) const; | ||
| 286 | }; | ||
| 287 | |||
| 288 | /// @brief A PDS3 value's text without its units and outer quotes, as `Pds3Label::Get` gives it. | ||
| 289 | /// @param value The value as written; a sequence keeps its units. | ||
| 290 | /// @return Its text. | ||
| 291 | ERNIE_API std::string CleanPds3Value(std::string_view value); | ||
| 292 | |||
| 293 | /// @brief Parses a PDS3 label's text up to its END line. | ||
| 294 | /// @param text The label: `KEY = value` lines, OBJECT and END_OBJECT nesting, comments (`/* */`) | ||
| 295 | /// and values continued over lines inside quotes, parentheses or braces. | ||
| 296 | /// @return Its keywords. | ||
| 297 | /// | ||
| 298 | /// References: | ||
| 299 | /// - NASA PDS, "Planetary Data System Standards Reference" 3.8, 2009 (the label's syntax, | ||
| 300 | /// chapter 12). | ||
| 301 | ERNIE_API Pds3Label ParsePds3Label(std::string_view text); | ||
| 302 | |||
| 303 | /// @brief A PDS3 image read in place from its mapped file. | ||
| 304 | /// Supports one band of 8/16/32-bit integers or floats in either byte order, SCALING_FACTOR, | ||
| 305 | /// OFFSET, MISSING_CONSTANT and simple cylindrical grids. The label is attached or in a .lbl | ||
| 306 | /// beside it. | ||
| 307 | /// | ||
| 308 | /// References: | ||
| 309 | /// - NASA PDS, "Planetary Data System Standards Reference" 3.8, 2009 (the label's syntax, chapter | ||
| 310 | /// 12; IMAGE objects, appendix A.21; IMAGE_MAP_PROJECTION, appendix A.22). | ||
| 311 | /// - LOLA team, "LRO LOLA Gridded Data Record" (LDEM_*.LBL). | ||
| 312 | class ERNIE_API PdsImage : public GridSource { | ||
| 313 | public: | ||
| 314 | /// @brief Opens an image. | ||
| 315 | /// @param path The label or data file; the label's ^IMAGE names the data beside it. | ||
| 316 | /// @throws std::runtime_error If the label cannot be found or read, or describes unsupported | ||
| 317 | /// samples. | ||
| 318 | explicit PdsImage(const std::filesystem::path& path); | ||
| 319 | |||
| 320 | 29 | std::uint32_t Width() const override { | |
| 321 | 29 | return width_; | |
| 322 | } | ||
| 323 | 27 | std::uint32_t Height() const override { | |
| 324 | 27 | return height_; | |
| 325 | } | ||
| 326 | 11 | std::uint32_t Channels() const override { | |
| 327 | 11 | return 1; | |
| 328 | } | ||
| 329 | 35 | GridSampleType Type() const override { | |
| 330 | 35 | return type_; | |
| 331 | } | ||
| 332 | 5 | double Scale() const override { | |
| 333 | 5 | return scale_; | |
| 334 | } | ||
| 335 | 5 | double Offset() const override { | |
| 336 | 5 | return offset_; | |
| 337 | } | ||
| 338 | 5 | std::optional<double> NoData() const override { | |
| 339 | 5 | return noData_; | |
| 340 | } | ||
| 341 | 5 | std::optional<LatLongGrid> Grid() const override { | |
| 342 | 5 | return grid_; | |
| 343 | } | ||
| 344 | void Read(std::uint32_t x, std::uint32_t y, std::uint32_t w, std::uint32_t h, | ||
| 345 | std::span<std::byte> out) const override; | ||
| 346 | |||
| 347 | /// @brief Reads a label value, without quotes or units. | ||
| 348 | /// @param key The key; the first object naming it wins. | ||
| 349 | /// @return The value, or none if absent. | ||
| 350 | std::optional<std::string> Label(std::string_view key) const; | ||
| 351 | |||
| 352 | /// @return The whole label, for keywords by their object path. | ||
| 353 | 12 | const Pds3Label& Keywords() const { | |
| 354 | 12 | return label_; | |
| 355 | } | ||
| 356 | |||
| 357 | private: | ||
| 358 | Pds3Label label_; | ||
| 359 | std::shared_ptr<const MappedFile> file_; | ||
| 360 | std::uint64_t start_ = 0; ///< The first sample's byte offset. | ||
| 361 | bool bigEndian_ = false; | ||
| 362 | std::uint32_t width_ = 0; | ||
| 363 | std::uint32_t height_ = 0; | ||
| 364 | GridSampleType type_ = GridSampleType::UInt8; | ||
| 365 | double scale_ = 1.0; | ||
| 366 | double offset_ = 0.0; | ||
| 367 | std::optional<double> noData_; | ||
| 368 | std::optional<LatLongGrid> grid_; | ||
| 369 | }; | ||
| 370 | |||
| 371 | // ---------------------------------------------------------------------------------------------- | ||
| 372 | // TIFF | ||
| 373 | |||
| 374 | /// @brief A TIFF or GeoTIFF image decoded on demand from its mapped file, so images larger than | ||
| 375 | /// memory read in parts. | ||
| 376 | /// Classic and BigTIFF in either byte order; strips or tiles; 8/16/32-bit integers and 32-bit | ||
| 377 | /// floats; none, LZW or Deflate with predictors 2 and 3. GeoTIFF pixel scale, tie point and GDAL | ||
| 378 | /// no-data. | ||
| 379 | /// | ||
| 380 | /// References: | ||
| 381 | /// - Adobe, "TIFF Revision 6.0", 1992 (the layout, LZW, predictor 2). | ||
| 382 | /// - Adobe, "Adobe Photoshop TIFF Technical Note 3", 2005 (Deflate, the floating-point predictor). | ||
| 383 | /// - Aware Systems, "BigTIFF" (64-bit offsets). | ||
| 384 | /// https://www.awaresystems.be/imaging/tiff/bigtiff.html | ||
| 385 | /// - OGC, "GeoTIFF Standard" 1.1, 2019 (ModelPixelScale, ModelTiepoint); GDAL's GDAL_NODATA tag | ||
| 386 | /// (42113). | ||
| 387 | class ERNIE_API TiffImage : public GridSource { | ||
| 388 | public: | ||
| 389 | /// @brief Opens an image and reads its layout. | ||
| 390 | /// @param path The file. | ||
| 391 | /// @throws std::runtime_error If the file cannot be read, is not a TIFF, or has an unsupported | ||
| 392 | /// layout or encoding. | ||
| 393 | explicit TiffImage(const std::filesystem::path& path); | ||
| 394 | |||
| 395 | 53 | std::uint32_t Width() const override { | |
| 396 | 53 | return width_; | |
| 397 | } | ||
| 398 | 51 | std::uint32_t Height() const override { | |
| 399 | 51 | return height_; | |
| 400 | } | ||
| 401 | 33 | std::uint32_t Channels() const override { | |
| 402 | 33 | return channels_; | |
| 403 | } | ||
| 404 | 295 | GridSampleType Type() const override { | |
| 405 | 295 | return type_; | |
| 406 | } | ||
| 407 | /// @return GeoTIFF's ModelPixelScale (x, y, z), if present. | ||
| 408 | 3 | std::optional<std::array<double, 3>> PixelScale() const { | |
| 409 | 3 | return pixelScale_; | |
| 410 | } | ||
| 411 | /// @return GeoTIFF's ModelTiepoint (raster i, j, k to model x, y, z), if present. | ||
| 412 | 3 | std::optional<std::array<double, 6>> Tiepoint() const { | |
| 413 | 3 | return tiepoint_; | |
| 414 | } | ||
| 415 | /// @return GDAL's no-data value, if present. | ||
| 416 | 2 | std::optional<double> NoData() const override { | |
| 417 | 2 | return noData_; | |
| 418 | } | ||
| 419 | /// @return The grid from the tie point and pixel scale, when both exist and the model space is | ||
| 420 | /// in degrees. | ||
| 421 | std::optional<LatLongGrid> Grid() const override; | ||
| 422 | void Read(std::uint32_t x, std::uint32_t y, std::uint32_t w, std::uint32_t h, | ||
| 423 | std::span<std::byte> out) const override; | ||
| 424 | |||
| 425 | private: | ||
| 426 | /// Decodes a strip or tile to host order, predictor undone. The last one decoded per plane is | ||
| 427 | /// kept for the next `Read`. | ||
| 428 | const std::vector<std::byte>& Chunk(std::size_t chunk) const; | ||
| 429 | |||
| 430 | /// Samples per pixel in a chunk: all channels, or one when planes are stored apart. | ||
| 431 | 247 | std::uint32_t ChunkChannels() const { | |
| 432 | 247 | return planes_ > 1 ? 1u : channels_; | |
| 433 | } | ||
| 434 | |||
| 435 | std::shared_ptr<const MappedFile> file_; | ||
| 436 | bool bigEndian_ = false; | ||
| 437 | std::uint32_t width_ = 0; | ||
| 438 | std::uint32_t height_ = 0; | ||
| 439 | std::uint32_t channels_ = 1; | ||
| 440 | std::uint32_t bits_ = 8; | ||
| 441 | GridSampleType type_ = GridSampleType::UInt8; | ||
| 442 | std::uint32_t compression_ = 1; | ||
| 443 | std::uint32_t predictor_ = 1; | ||
| 444 | /// Strips are `width_` wide and `chunkHeight_` rows; tiles `chunkWidth_` by `chunkHeight_`, | ||
| 445 | /// `chunksAcross_` to a row. | ||
| 446 | bool tiled_ = false; | ||
| 447 | std::uint32_t chunkWidth_ = 0; | ||
| 448 | std::uint32_t chunkHeight_ = 0; | ||
| 449 | std::uint32_t chunksAcross_ = 1; | ||
| 450 | /// Chunks of one plane; with the planes stored apart (PlanarConfiguration 2), `planes_` of | ||
| 451 | /// them one after another. | ||
| 452 | std::size_t chunksPerPlane_ = 0; | ||
| 453 | std::uint32_t planes_ = 1; | ||
| 454 | std::vector<std::uint64_t> offsets_; | ||
| 455 | std::vector<std::uint64_t> byteCounts_; | ||
| 456 | std::optional<std::array<double, 3>> pixelScale_; | ||
| 457 | std::optional<std::array<double, 6>> tiepoint_; | ||
| 458 | std::optional<double> noData_; | ||
| 459 | mutable std::vector<std::size_t> cachedChunks_; | ||
| 460 | mutable std::vector<std::vector<std::byte>> cached_; | ||
| 461 | }; | ||
| 462 | |||
| 463 | // ---------------------------------------------------------------------------------------------- | ||
| 464 | // Planetary DEMs | ||
| 465 | // | ||
| 466 | // Digital elevation models as the PDS archives them, and the local heightfields Ernie's Terrain | ||
| 467 | // draws of them. | ||
| 468 | // | ||
| 469 | // References: | ||
| 470 | // - G. A. Neumann, "Lunar Reconnaissance Orbiter Lunar Orbiter Laser Altimeter Reduced Data Record | ||
| 471 | // and Derived Products Software Interface Specification", version 2.42, 2010 (the LDEM products, | ||
| 472 | // their labels and simple cylindrical grids). | ||
| 473 | // - J. P. Snyder, "Map Projections: A Working Manual", USGS Professional Paper 1395, 1987 (the | ||
| 474 | // equirectangular projection; the polar stereographic, eq. 21-33 and 21-34). | ||
| 475 | |||
| 476 | /// @brief A planet's radius on an equirectangular grid, as LOLA's LDEM products (simple | ||
| 477 | /// cylindrical, true at the equator) and LROC's NAC DTMs (true at a standard parallel) hold it. | ||
| 478 | /// `heights` are metres above `referenceRadius` (NaN where the product has no data), row by row | ||
| 479 | /// from the northern, each row from the western sample eastwards, `pixelsPerDegree` along the | ||
| 480 | /// meridians and that times cos(standardParallel) along the parallels. Sample s of line l is | ||
| 481 | /// centred at latitude maxLatitude - (l + 0.5) / ppd and east longitude | ||
| 482 | /// westLongitude + (s + 0.5) / (ppd cos(standardParallel)). | ||
| 483 | 3 | struct ERNIE_API CylindricalDem { | |
| 484 | std::uint32_t lines = 0; | ||
| 485 | std::uint32_t samples = 0; | ||
| 486 | double pixelsPerDegree = 1.0; | ||
| 487 | double maxLatitudeDeg = 90.0; | ||
| 488 | double westLongitudeDeg = 0.0; | ||
| 489 | double standardParallelDeg = 0.0; | ||
| 490 | double referenceRadius = 0.0; ///< Metres. | ||
| 491 | std::vector<float> heights; | ||
| 492 | |||
| 493 | /// @brief The height at a point, bilinear between the cells' centres, clamped to the grid's | ||
| 494 | /// edges. | ||
| 495 | /// @param latitudeDeg The latitude. | ||
| 496 | /// @param longitudeDeg The east longitude, any multiple of 360 apart. | ||
| 497 | /// @return The height, or NaN beside a cell without data. | ||
| 498 | double HeightAt(double latitudeDeg, double longitudeDeg) const; | ||
| 499 | }; | ||
| 500 | |||
| 501 | /// @brief Reads the DEM a PDS3 label describes, on its simple cylindrical or equirectangular grid | ||
| 502 | /// (true at CENTER_LATITUDE). | ||
| 503 | /// Heights are samples times SCALING_FACTOR, plus OFFSET unless OFFSET is the reference sphere's | ||
| 504 | /// radius (as LOLA's are); PDS's real null (below -1e30), MISSING_CONSTANT and non-finite samples | ||
| 505 | /// are no data. | ||
| 506 | /// @param labelPath The label (or the data, with its label attached or beside it). | ||
| 507 | /// @return The DEM. | ||
| 508 | /// @throws std::runtime_error If a file cannot be read. | ||
| 509 | /// @throws std::invalid_argument For another projection or a label it cannot follow. | ||
| 510 | ERNIE_API CylindricalDem ReadCylindricalDem(const std::filesystem::path& labelPath); | ||
| 511 | |||
| 512 | /// @brief A planet's radius on a polar stereographic grid (spherical, true at the pole), as | ||
| 513 | /// LOLA's polar LDEM products hold it. | ||
| 514 | /// `heights` are metres above `referenceRadius`, row by row, `mapScale` metres a pixel at the | ||
| 515 | /// pole. A point at polar distance rho = 2 R tan(45 + lat / 2) (south) or 2 R tan(45 - lat / 2) | ||
| 516 | /// (north) and longitude lon from the centre longitude lies at x = rho sin(lon) and y = | ||
| 517 | /// rho cos(lon) (south) or -rho cos(lon) (north), the cell centre at line | ||
| 518 | /// lineProjectionOffset - y / mapScale and sample sampleProjectionOffset + x / mapScale. | ||
| 519 | 1 | struct ERNIE_API PolarStereographicDem { | |
| 520 | std::uint32_t lines = 0; | ||
| 521 | std::uint32_t samples = 0; | ||
| 522 | double mapScale = 1.0; ///< Metres a pixel. | ||
| 523 | double lineProjectionOffset = 0.0; | ||
| 524 | double sampleProjectionOffset = 0.0; | ||
| 525 | double centerLatitudeDeg = -90.0; ///< -90 or 90. | ||
| 526 | double centerLongitudeDeg = 0.0; | ||
| 527 | double referenceRadius = 0.0; ///< Metres. | ||
| 528 | std::vector<float> heights; | ||
| 529 | |||
| 530 | /// @brief The fractional (line, sample) of a point. | ||
| 531 | /// @param latitudeDeg The latitude. | ||
| 532 | /// @param longitudeDeg The east longitude. | ||
| 533 | /// @return Its line and sample. | ||
| 534 | glm::dvec2 ToGrid(double latitudeDeg, double longitudeDeg) const; | ||
| 535 | /// @brief The height at a point, bilinear between the cells' centres, clamped to the grid's | ||
| 536 | /// edges. | ||
| 537 | /// @param latitudeDeg The latitude. | ||
| 538 | /// @param longitudeDeg The east longitude. | ||
| 539 | /// @return The height. | ||
| 540 | double HeightAt(double latitudeDeg, double longitudeDeg) const; | ||
| 541 | }; | ||
| 542 | |||
| 543 | /// @brief Reads a DEM as `ReadCylindricalDem` does, but on a polar stereographic grid. | ||
| 544 | /// @param labelPath The label. | ||
| 545 | /// @return The DEM. | ||
| 546 | /// @throws std::invalid_argument As `ReadCylindricalDem`, and for a projection not centred on a | ||
| 547 | /// pole. | ||
| 548 | ERNIE_API PolarStereographicDem ReadPolarStereographicDem(const std::filesystem::path& labelPath); | ||
| 549 | |||
| 550 | /// @brief A heightfield of a DEM about one point, in a local frame there: x east, y north and z | ||
| 551 | /// up, its origin on the reference sphere at the point. | ||
| 552 | /// Each cell holds the height above that tangent plane of the surface straight below or above it, | ||
| 553 | /// so the planet's curvature is in the heights (a drop of about x^2 / 2R). Where the DEM has no | ||
| 554 | /// data the surface is the reference sphere. | ||
| 555 | 3 | struct LocalTerrain { | |
| 556 | Tensor heights; ///< float32 (rows, cols), the first northern. | ||
| 557 | glm::dvec2 origin; ///< x and y of the north-western cell. | ||
| 558 | double cellSize = 1.0; | ||
| 559 | /// Body-fixed coordinates (metres, the planet's centre at the origin) to the local frame. | ||
| 560 | glm::dmat4 bodyToLocal{1.0}; | ||
| 561 | }; | ||
| 562 | |||
| 563 | /// @brief Makes the heightfield of a DEM about a point. | ||
| 564 | /// @param dem The DEM. | ||
| 565 | /// @param latitudeDeg The point's latitude. | ||
| 566 | /// @param longitudeDeg The point's east longitude. | ||
| 567 | /// @param cells Cells along each side, at least 2. | ||
| 568 | /// @param cellSize Metres between cells; positive. | ||
| 569 | /// @return The `cells` x `cells` heightfield. | ||
| 570 | /// @throws std::invalid_argument For fewer than 2 cells or a non-positive cell size. | ||
| 571 | ERNIE_API LocalTerrain MakeLocalTerrain(const CylindricalDem& dem, double latitudeDeg, | ||
| 572 | double longitudeDeg, std::uint32_t cells, double cellSize); | ||
| 573 | /// @copydoc MakeLocalTerrain(const CylindricalDem&, double, double, std::uint32_t, double) | ||
| 574 | ERNIE_API LocalTerrain MakeLocalTerrain(const PolarStereographicDem& dem, double latitudeDeg, | ||
| 575 | double longitudeDeg, std::uint32_t cells, double cellSize); | ||
| 576 | |||
| 577 | /// @brief Makes a heightfield as `MakeLocalTerrain` does, taking `fallback`'s heights wherever | ||
| 578 | /// `dem` has none: a high-resolution DTM inside a regional DEM. Both share `dem`'s reference | ||
| 579 | /// radius. | ||
| 580 | /// @param dem The DEM. | ||
| 581 | /// @param fallback The DEM used where `dem` has no data. | ||
| 582 | /// @param latitudeDeg The point's latitude. | ||
| 583 | /// @param longitudeDeg The point's east longitude. | ||
| 584 | /// @param cells Cells along each side, at least 2. | ||
| 585 | /// @param cellSize Metres between cells; positive. | ||
| 586 | /// @return The heightfield. | ||
| 587 | ERNIE_API LocalTerrain MakeLocalTerrain(const CylindricalDem& dem, const CylindricalDem& fallback, | ||
| 588 | double latitudeDeg, double longitudeDeg, | ||
| 589 | std::uint32_t cells, double cellSize); | ||
| 590 | |||
| 591 | /// @brief The body-fixed unit vector of a planetocentric latitude and east longitude. | ||
| 592 | /// @param latitudeDeg The latitude. | ||
| 593 | /// @param longitudeDeg The east longitude. | ||
| 594 | /// @return The direction. | ||
| 595 | ERNIE_API glm::dvec3 BodyDirection(double latitudeDeg, double longitudeDeg); | ||
| 596 | |||
| 597 | } // namespace ernie | ||
| 598 |