TerrainGen coverage

terraingen/
Coverage:
low: ≥ 0%
medium: ≥ 75.0%
high: ≥ 90.0%
Lines:
74 / 74
100.0%
Branches:
80 / 82
97.6%
Functions:
9 / 9
100.0%
raster/resample.cpp
Line Branch Exec Source
1 #include "terraingen/raster/resample.h"
2
3 #include <algorithm>
4 #include <array>
5 #include <cmath>
6 #include <limits>
7 #include <optional>
8 #include <stdexcept>
9
10 namespace terraingen {
11
12 namespace {
13
14 /// Keys cubic convolution kernel with a = -0.5.
15 14352332 double CubicWeight(double x) {
16 14352332 constexpr double a = -0.5;
17 14352332 x = std::fabs(x);
18
2/2
✓ Branch 0 taken 7220037 times.
✓ Branch 1 taken 7132295 times.
 14352332 if (x <= 1.0) {
19 7220037 return ((a + 2.0) * x - (a + 3.0)) * x * x + 1.0;
20 }
21
2/2
✓ Branch 0 taken 7088428 times.
✓ Branch 1 taken 43867 times.
 7132295 if (x < 2.0) {
22 7088428 return ((a * x - 5.0 * a) * x + 8.0 * a) * x - 4.0 * a;
23 }
24 43867 return 0.0;
25 }
26
27 14650088 std::size_t Clamp(double index, std::size_t size) {
28
1/1
✓ Branch 1 taken 14650088 times.
 14650088 return static_cast<std::size_t>(std::clamp(index, 0.0, static_cast<double>(size - 1)));
29 }
30
31 /// `raster` interpolated at pixel coordinates (`u`, `v`) measured from the centre of the first
32 /// cell, or nothing where a contributing cell holds no data.
33 767439 std::optional<float> Interpolate(const Raster<float>& raster, double u, double v,
34 ResampleMethod method) {
35
2/2
✓ Branch 0 taken 111 times.
✓ Branch 1 taken 767328 times.
 767439 if (method == ResampleMethod::Nearest) {
36
2/2
✓ Branch 2 taken 111 times.
✓ Branch 6 taken 111 times.
 111 return raster(Clamp(std::round(u), raster.Width()), Clamp(std::round(v), raster.Height()));
37 }
38 767328 const double u0 = std::floor(u);
39 767328 const double v0 = std::floor(v);
40 767328 const double fu = u - u0;
41 767328 const double fv = v - v0;
42
2/2
✓ Branch 0 taken 49595 times.
✓ Branch 1 taken 717733 times.
 767328 const int first = method == ResampleMethod::Bilinear ? 0 : -1;
43
2/2
✓ Branch 0 taken 49595 times.
✓ Branch 1 taken 717733 times.
 767328 const int last = method == ResampleMethod::Bilinear ? 1 : 2;
44 767328 double sum = 0.0;
45
2/2
✓ Branch 0 taken 2969725 times.
✓ Branch 1 taken 767053 times.
 3736778 for (int j = first; j <= last; ++j) {
46 const double wy =
47
4/4
✓ Branch 0 taken 99181 times.
✓ Branch 1 taken 2870544 times.
✓ Branch 2 taken 49595 times.
✓ Branch 3 taken 49586 times.
 2969725 method == ResampleMethod::Bilinear ? (j == 0 ? 1.0 - fv : fv) : CubicWeight(fv - j);
48 2969725 const std::size_t y = Clamp(v0 + j, raster.Height());
49
2/2
✓ Branch 0 taken 11680141 times.
✓ Branch 1 taken 2969450 times.
 14649591 for (int i = first; i <= last; ++i) {
50 const double wx =
51
4/4
✓ Branch 0 taken 198353 times.
✓ Branch 1 taken 11481788 times.
✓ Branch 2 taken 99181 times.
✓ Branch 3 taken 99172 times.
 11680141 method == ResampleMethod::Bilinear ? (i == 0 ? 1.0 - fu : fu) : CubicWeight(fu - i);
52
1/1
✓ Branch 2 taken 11680141 times.
 11680141 const float value = raster(Clamp(u0 + i, raster.Width()), y);
53
6/6
✓ Branch 0 taken 11475035 times.
✓ Branch 1 taken 205106 times.
✓ Branch 3 taken 275 times.
✓ Branch 4 taken 11474760 times.
✓ Branch 5 taken 275 times.
✓ Branch 6 taken 11679866 times.
 11680141 if (wx * wy != 0.0 && raster.IsNoData(value)) {
54 275 return std::nullopt;
55 }
56 11679866 sum += wx * wy * static_cast<double>(value);
57 }
58 }
59 767053 return static_cast<float>(sum);
60 }
61
62 158328 void CheckSource(const Raster<float>& raster) {
63
2/2
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 158327 times.
 158328 if (raster.Empty()) {
64
1/1
✓ Branch 2 taken 1 time.
 1 throw std::invalid_argument("Cannot resample an empty raster");
65 }
66 158327 const GeoTransform& source = raster.Transform();
67
3/4
✓ Branch 0 taken 158326 times.
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
✓ Branch 3 taken 158326 times.
 158327 if (source.rowRotation != 0.0 || source.colRotation != 0.0) {
68
1/1
✓ Branch 2 taken 1 time.
 1 throw std::invalid_argument("Cannot resample a rotated raster");
69 }
70 158326 }
71
72 } // namespace
73
74 158324 float SampleRaster(const Raster<float>& raster, double x, double y, ResampleMethod method) {
75
1/1
✓ Branch 1 taken 158323 times.
 158324 CheckSource(raster);
76
1/1
✓ Branch 2 taken 158323 times.
 158323 const Vec2d pixel = raster.Transform().WorldToPixel(x, y);
77
8/8
✓ Branch 0 taken 157990 times.
✓ Branch 1 taken 333 times.
✓ Branch 2 taken 157723 times.
✓ Branch 3 taken 267 times.
✓ Branch 5 taken 157718 times.
✓ Branch 6 taken 5 times.
✓ Branch 7 taken 609 times.
✓ Branch 8 taken 157714 times.
 316041 if (!(pixel.x >= 0.0 && pixel.y >= 0.0 && pixel.x < static_cast<double>(raster.Width()) &&
78
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 157714 times.
 157718 pixel.y < static_cast<double>(raster.Height()))) {
79 609 return std::numeric_limits<float>::quiet_NaN();
80 }
81
1/1
✓ Branch 1 taken 157714 times.
 157714 const std::optional<float> value = Interpolate(raster, pixel.x - 0.5, pixel.y - 0.5, method);
82
6/6
✓ Branch 1 taken 157456 times.
✓ Branch 2 taken 258 times.
✓ Branch 5 taken 1 time.
✓ Branch 6 taken 157455 times.
✓ Branch 7 taken 259 times.
✓ Branch 8 taken 157455 times.
 157714 if (!value || raster.IsNoData(*value)) {
83 259 return std::numeric_limits<float>::quiet_NaN();
84 }
85 157455 return *value;
86 }
87
88 4 Raster<float> ResampleTo(const Raster<float>& raster, const RasterGrid& grid,
89 ResampleMethod method) {
90 4 CheckSource(raster);
91
3/4
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
 3 if (grid.width == 0 || grid.height == 0) {
92
1/1
✓ Branch 2 taken 1 time.
 1 throw std::invalid_argument("Cannot resample onto an empty grid");
93 }
94 2 Raster<float> result = grid.MakeRaster();
95
1/1
✓ Branch 1 taken 2 times.
 2 SampleAtCentres(result, raster.NoData(),
96 54 [&](double x, double y) { return SampleRaster(raster, x, y, method); });
97 2 return result;
98 }
99
100 38 std::optional<float> FilledNoData(std::optional<float> noData, bool missing) {
101
6/6
✓ Branch 1 taken 31 times.
✓ Branch 2 taken 7 times.
✓ Branch 3 taken 26 times.
✓ Branch 4 taken 5 times.
✓ Branch 5 taken 33 times.
✓ Branch 6 taken 5 times.
 38 if (noData || !missing) {
102 33 return noData;
103 }
104 5 return std::numeric_limits<float>::quiet_NaN();
105 }
106
107 23 Raster<float> Resample(const Raster<float>& raster, double pixelSize, ResampleMethod method) {
108
2/2
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 22 times.
 23 if (raster.Empty()) {
109
1/1
✓ Branch 2 taken 1 time.
 1 throw std::invalid_argument("Cannot resample an empty raster");
110 }
111
3/3
✓ Branch 1 taken 22 times.
✓ Branch 4 taken 18 times.
✓ Branch 7 taken 18 times.
 22 Raster<float> result = ResampledGrid(RasterGrid::FromRaster(raster), pixelSize).MakeRaster();
112 18 result.SetNoData(raster.NoData());
113 18 const double scaleX = std::fabs(raster.Transform().pixelWidth) / pixelSize;
114 18 const double scaleY = std::fabs(raster.Transform().pixelHeight) / pixelSize;
115 18 const float fill = raster.NoData().value_or(std::nanf(""));
116
2/2
✓ Branch 1 taken 1749 times.
✓ Branch 2 taken 18 times.
 1767 for (std::size_t row = 0; row < result.Height(); ++row) {
117 1749 const double v = (static_cast<double>(row) + 0.5) / scaleY - 0.5;
118
2/2
✓ Branch 1 taken 609725 times.
✓ Branch 2 taken 1749 times.
 611474 for (std::size_t col = 0; col < result.Width(); ++col) {
119 609725 const double u = (static_cast<double>(col) + 0.5) / scaleX - 0.5;
120
1/1
✓ Branch 1 taken 609725 times.
 609725 const std::optional<float> value = Interpolate(raster, u, v, method);
121 609725 result(col, row) = value.value_or(fill);
122 }
123 }
124 36 return result;
125 }
126
127 } // namespace terraingen
128