TerrainGen coverage

terraingen/
Coverage:
low: ≥ 0%
medium: ≥ 75.0%
high: ≥ 90.0%
Lines:
111 / 111
100.0%
Branches:
107 / 119
89.9%
Functions:
14 / 14
100.0%
sources/geotiff_source.cpp
Line Branch Exec Source
1 #include "terraingen/sources/geotiff_source.h"
2
3 #include "terraingen/raster/gaussian_blur.h"
4 #include "terraingen/util/checks.h"
5
6 #include <cmath>
7 #include <limits>
8 #include <stdexcept>
9 #include <string>
10
11 namespace terraingen {
12
13 namespace {
14
15 /// `value` as an integer when it is one to within rounding, as offsets between lattices are.
16 64 std::optional<long long> Whole(double value) {
17 64 const double rounded = std::round(value);
18
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 60 times.
 64 if (std::fabs(value - rounded) > 1e-6) {
19 4 return std::nullopt;
20 }
21 60 return static_cast<long long>(rounded);
22 }
23
24 } // namespace
25
26 46 GeoTiffSource::GeoTiffSource(std::filesystem::path path, GeoTiffSourceOptions options)
27
1/1
✓ Branch 4 taken 45 times.
 46 : path_(std::move(path)), options_(options), info_(ReadGeoTiffInfo(path_)) {
28 45 const std::size_t band = options_.band;
29
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 44 times.
 45 if (band >= info_.bands) {
30
4/4
✓ Branch 1 taken 1 time.
✓ Branch 4 taken 1 time.
✓ Branch 7 taken 1 time.
✓ Branch 10 taken 1 time.
 2 throw std::invalid_argument("Band " + std::to_string(band) + " is missing; the file has " +
31
2/2
✓ Branch 1 taken 1 time.
✓ Branch 4 taken 1 time.
 4 std::to_string(info_.bands) +
32
2/3
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 time.
 3 (info_.bands == 1 ? " band" : " bands"));
33 }
34 44 const std::optional<double>& scale = options_.scale;
35 44 const std::optional<double>& offset = options_.offset;
36
10/10
✓ Branch 1 taken 7 times.
✓ Branch 2 taken 37 times.
✓ Branch 5 taken 6 times.
✓ Branch 6 taken 1 time.
✓ Branch 8 taken 3 times.
✓ Branch 9 taken 40 times.
✓ Branch 12 taken 1 time.
✓ Branch 13 taken 2 times.
✓ Branch 14 taken 2 times.
✓ Branch 15 taken 42 times.
 44 if ((scale && !std::isfinite(*scale)) || (offset && !std::isfinite(*offset))) {
37
1/1
✓ Branch 2 taken 2 times.
 2 throw std::invalid_argument("Scale and offset must be finite");
38 }
39
1/1
✓ Branch 2 taken 41 times.
 42 RequireNonNegative(options_.lowPass, "Low-pass deviation");
40 41 scaling_ = {scale.value_or(info_.bandScaling[band].scale),
41 41 offset.value_or(info_.bandScaling[band].offset)};
42
3/3
✓ Branch 1 taken 41 times.
✓ Branch 3 taken 1 time.
✓ Branch 4 taken 40 times.
 41 if (info_.geoKeys.IsGeographic()) {
43 throw std::invalid_argument(
44 "The raster's CRS is geographic (degrees); craters are sized in metres, so reproject "
45
1/1
✓ Branch 2 taken 1 time.
 1 "it first, for example with gdalwarp -t_srs");
46 }
47 40 grid_ = {info_.width, info_.height, info_.transform};
48
2/2
✓ Branch 1 taken 9 times.
✓ Branch 2 taken 31 times.
 40 if (options_.pixelSize) {
49
1/1
✓ Branch 2 taken 7 times.
 9 grid_ = ResampledGrid(grid_, *options_.pixelSize);
50 }
51 75 }
52
53 5 Vec2d GeoTiffSource::ReferencePoint() const {
54 5 return {grid_.transform.originX, grid_.transform.originY};
55 }
56
57 4958 void GeoTiffSource::Load() const {
58
1/1
✓ Branch 1 taken 4958 times.
 4958 std::call_once(loaded_, [this] {
59
1/1
✓ Branch 1 taken 33 times.
 33 GeoTiff tiff = ReadGeoTiff(path_, options_.band);
60
2/2
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 25 times.
 33 if (!scaling_.Identity()) {
61 8 Raster<float>& raster = tiff.raster;
62 8 constexpr float kMissing = std::numeric_limits<float>::quiet_NaN();
63
2/2
✓ Branch 5 taken 4683 times.
✓ Branch 6 taken 8 times.
 4691 for (float& value : raster.Values()) {
64 9366 value = raster.IsNoData(value)
65
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 4678 times.
 4683 ? kMissing
66 4678 : static_cast<float>(static_cast<double>(value) * scaling_.scale +
67 4678 scaling_.offset);
68 }
69
2/2
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 3 times.
 8 if (raster.NoData()) {
70 5 raster.SetNoData(kMissing);
71 }
72 }
73
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 28 times.
 33 if (options_.lowPass > 0.0) {
74 5 const GeoTransform& transform = tiff.raster.Transform();
75 tiff.raster =
76 5 GaussianBlur(tiff.raster, options_.lowPass / std::fabs(transform.pixelWidth),
77
1/1
✓ Branch 1 taken 5 times.
 5 options_.lowPass / std::fabs(transform.pixelHeight));
78 }
79
1/1
✓ Branch 2 taken 33 times.
 33 native_ = std::make_shared<const Raster<float>>(std::move(tiff.raster));
80
6/8
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 28 times.
✓ Branch 4 taken 5 times.
✓ Branch 7 taken 5 times.
✓ Branch 10 taken 5 times.
✓ Branch 11 taken 28 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 71 heights_ = options_.pixelSize ? std::make_shared<const Raster<float>>(Resample(
81 5 *native_, *options_.pixelSize, options_.resampling))
82 33 : native_;
83 33 });
84 4958 }
85
86 4811 const Raster<float>& GeoTiffSource::Heights() const {
87 4811 Load();
88 4811 return *heights_;
89 }
90
91 147 const Raster<float>& GeoTiffSource::Native() const {
92 147 Load();
93 147 return *native_;
94 }
95
96 std::optional<GeoTiffSource::LatticeOffset>
97 32 GeoTiffSource::OffsetOnLattice(const GeoTransform& target) const {
98 32 const GeoTransform& source = grid_.transform;
99 62 const auto same = [](double a, double b) { return std::fabs(a - b) <= 1e-12 * std::fabs(b); };
100 const std::optional<long long> col =
101 32 Whole((target.originX - source.originX) / source.pixelWidth);
102 const std::optional<long long> row =
103 32 Whole((target.originY - source.originY) / source.pixelHeight);
104 32 if (!same(target.pixelWidth, source.pixelWidth) ||
105
3/4
✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 29 times.
✓ Branch 4 taken 1 time.
 30 !same(target.pixelHeight, source.pixelHeight) || target.rowRotation != source.rowRotation ||
106
8/10
✓ Branch 0 taken 30 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 29 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 28 times.
✓ Branch 6 taken 1 time.
✗ Branch 8 not taken.
✓ Branch 9 taken 28 times.
✓ Branch 10 taken 4 times.
✓ Branch 11 taken 28 times.
 62 target.colRotation != source.colRotation || !col || !row) {
107 4 return std::nullopt;
108 }
109 28 return LatticeOffset{*col, *row};
110 }
111
112 28 void GeoTiffSource::CopyFromLattice(Raster<float>& heights, LatticeOffset offset) const {
113 28 const Raster<float>& data = Heights();
114 28 const float fill = data.NoData().value_or(std::numeric_limits<float>::quiet_NaN());
115 28 bool outside = false;
116
2/2
✓ Branch 1 taken 2543 times.
✓ Branch 2 taken 28 times.
 2571 for (std::size_t row = 0; row < heights.Height(); ++row) {
117 2543 const long long sourceRow = offset.row + static_cast<long long>(row);
118
2/2
✓ Branch 1 taken 918058 times.
✓ Branch 2 taken 2543 times.
 920601 for (std::size_t col = 0; col < heights.Width(); ++col) {
119 918058 const long long sourceCol = offset.col + static_cast<long long>(col);
120
1/2
✓ Branch 0 taken 918058 times.
✗ Branch 1 not taken.
 918058 if (sourceRow < 0 || sourceCol < 0 ||
121
5/6
✓ Branch 0 taken 918058 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 918054 times.
✓ Branch 4 taken 4 times.
✓ Branch 5 taken 8 times.
✓ Branch 6 taken 918050 times.
 2754170 sourceRow >= static_cast<long long>(data.Height()) ||
122
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 918050 times.
 918054 sourceCol >= static_cast<long long>(data.Width())) {
123 8 heights(col, row) = fill;
124 8 outside = true;
125 8 continue;
126 }
127 918050 heights(col, row) =
128 918050 data(static_cast<std::size_t>(sourceCol), static_cast<std::size_t>(sourceRow));
129 }
130 }
131 28 heights.SetNoData(FilledNoData(data.NoData(), outside));
132 28 }
133
134 4 void GeoTiffSource::FillBySampling(Raster<float>& heights) const {
135 4 const GeoTransform& target = heights.Transform();
136
3/4
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
 4 if (target.rowRotation != 0.0 || target.colRotation != 0.0) {
137
2/2
✓ Branch 1 taken 1 time.
✓ Branch 4 taken 1 time.
 2 throw std::invalid_argument("A rotated grid off the lattice of " + path_.string() +
138
2/2
✓ Branch 2 taken 1 time.
✓ Branch 5 taken 1 time.
 3 " cannot be filled");
139 }
140
1/1
✓ Branch 2 taken 3 times.
 3 SampleAtCentres(heights, Heights().NoData(),
141
1/1
✓ Branch 1 taken 3 times.
 100 [this](double x, double y) { return HeightAt(x, y); });
142 3 }
143
144 32 void GeoTiffSource::Fill(Raster<float>& heights) const {
145
2/2
✓ Branch 3 taken 28 times.
✓ Branch 4 taken 4 times.
 32 if (const auto offset = OffsetOnLattice(heights.Transform())) {
146
1/1
✓ Branch 2 taken 28 times.
 28 CopyFromLattice(heights, *offset);
147 } else {
148
1/1
✓ Branch 1 taken 3 times.
 4 FillBySampling(heights);
149 }
150 31 }
151
152 4780 float GeoTiffSource::HeightAt(double x, double y) const {
153
1/1
✓ Branch 1 taken 4780 times.
 4780 const Raster<float>& data = Heights();
154
1/1
✓ Branch 1 taken 4780 times.
 4780 const Vec2d pixel = grid_.transform.WorldToPixel(x, y);
155 4780 const double col = std::floor(pixel.x);
156 4780 const double row = std::floor(pixel.y);
157 4780 const bool centred =
158
3/4
✓ Branch 0 taken 4633 times.
✓ Branch 1 taken 147 times.
✓ Branch 2 taken 4633 times.
✗ Branch 3 not taken.
 4780 std::fabs(pixel.x - col - 0.5) < 1e-9 && std::fabs(pixel.y - row - 0.5) < 1e-9;
159
2/2
✓ Branch 0 taken 147 times.
✓ Branch 1 taken 4633 times.
 4780 if (!centred) {
160
2/2
✓ Branch 1 taken 147 times.
✓ Branch 4 taken 147 times.
 147 return SampleRaster(Native(), x, y, options_.resampling);
161 }
162
7/8
✓ Branch 0 taken 4623 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 4615 times.
✓ Branch 3 taken 8 times.
✓ Branch 5 taken 4615 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 18 times.
✓ Branch 8 taken 4615 times.
 9248 if (!(col >= 0.0 && row >= 0.0 && col < static_cast<double>(data.Width()) &&
163
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 4615 times.
 4615 row < static_cast<double>(data.Height()))) {
164 18 return std::numeric_limits<float>::quiet_NaN();
165 }
166 4615 const float value = data(static_cast<std::size_t>(col), static_cast<std::size_t>(row));
167
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 4613 times.
 4615 return data.IsNoData(value) ? std::numeric_limits<float>::quiet_NaN() : value;
168 }
169
170 } // namespace terraingen
171