Line data Source code
1 : ///////////////////////////////////////////////////////////////////////////////
2 : //
3 : // Project: C++ Test Suite for GDAL/OGR
4 : // Purpose: Test viewshed algorithm
5 : // Author: Andrew Bell
6 : //
7 : ///////////////////////////////////////////////////////////////////////////////
8 : /*
9 : * SPDX-License-Identifier: MIT
10 : ****************************************************************************/
11 :
12 : #include <algorithm>
13 : #include <array>
14 : #include <utility>
15 :
16 : #include "gdal_unit_test.h"
17 :
18 : #include "gtest_include.h"
19 :
20 : #include "viewshed/viewshed.h"
21 :
22 : namespace gdal
23 : {
24 : namespace viewshed
25 : {
26 :
27 : namespace
28 : {
29 : using Coord = std::pair<int, int>;
30 : using DatasetPtr = std::unique_ptr<GDALDataset>;
31 : using Transform = std::array<double, 6>;
32 : Transform identity{0, 1, 0, 0, 0, 1};
33 :
34 17 : Options stdOptions(int x, int y)
35 : {
36 17 : Options opts;
37 17 : opts.observer.x = x;
38 17 : opts.observer.y = y;
39 17 : opts.outputFilename = "none";
40 17 : opts.outputFormat = "mem";
41 17 : opts.curveCoeff = 0;
42 :
43 17 : return opts;
44 : }
45 :
46 5 : Options stdOptions(const Coord &observer)
47 : {
48 5 : return stdOptions(observer.first, observer.second);
49 : }
50 :
51 22 : DatasetPtr runViewshed(const int8_t *in, int xlen, int ylen,
52 : const Options &opts)
53 : {
54 44 : Viewshed v(opts);
55 :
56 22 : GDALDriver *driver = (GDALDriver *)GDALGetDriverByName("MEM");
57 22 : GDALDataset *dataset = driver->Create("", xlen, ylen, 1, GDT_Int8, nullptr);
58 22 : EXPECT_TRUE(dataset);
59 22 : dataset->SetGeoTransform(identity.data());
60 22 : GDALRasterBand *band = dataset->GetRasterBand(1);
61 22 : EXPECT_TRUE(band);
62 22 : CPLErr err = band->RasterIO(GF_Write, 0, 0, xlen, ylen, (int8_t *)in, xlen,
63 22 : ylen, GDT_Int8, 0, 0, nullptr);
64 22 : EXPECT_EQ(err, CE_None);
65 :
66 22 : EXPECT_TRUE(v.run(band));
67 44 : return v.output();
68 : }
69 :
70 : } // namespace
71 :
72 4 : TEST(Viewshed, all_visible)
73 : {
74 : // clang-format off
75 1 : const int xlen = 3;
76 1 : const int ylen = 3;
77 1 : std::array<int8_t, xlen * ylen> in
78 : {
79 : 1, 2, 3,
80 : 4, 5, 6,
81 : 3, 2, 1
82 : };
83 : // clang-format on
84 :
85 2 : SCOPED_TRACE("all_visible");
86 2 : DatasetPtr output = runViewshed(in.data(), xlen, ylen, stdOptions(1, 1));
87 :
88 : std::array<uint8_t, xlen * ylen> out;
89 1 : GDALRasterBand *band = output->GetRasterBand(1);
90 1 : CPLErr err = band->RasterIO(GF_Read, 0, 0, xlen, ylen, out.data(), xlen,
91 1 : ylen, GDT_Int8, 0, 0, nullptr);
92 1 : EXPECT_EQ(err, CE_None);
93 :
94 10 : for (uint8_t o : out)
95 9 : EXPECT_EQ(o, 127);
96 1 : }
97 :
98 4 : TEST(Viewshed, simple_height)
99 : {
100 : // clang-format off
101 1 : const int xlen = 5;
102 1 : const int ylen = 5;
103 1 : std::array<int8_t, xlen * ylen> in
104 : {
105 : -1, 0, 1, 0, -1,
106 : -1, 2, 0, 4, -1,
107 : -1, 1, 0, -1, -1,
108 : 0, 3, 0, 2, 0,
109 : -1, 0, 0, 3, -1
110 : };
111 :
112 1 : std::array<double, xlen * ylen> observable
113 : {
114 : 4, 2, 0, 4, 8,
115 : 3, 2, 0, 4, 3,
116 : 2, 1, 0, -1, -2,
117 : 4, 3, 0, 2, 1,
118 : 6, 3, 0, 2, 4
119 : };
120 : // clang-format on
121 :
122 : {
123 2 : SCOPED_TRACE("simple_height:normal");
124 :
125 : DatasetPtr output =
126 2 : runViewshed(in.data(), xlen, ylen, stdOptions(2, 2));
127 :
128 : std::array<int8_t, xlen * ylen> out;
129 1 : GDALRasterBand *band = output->GetRasterBand(1);
130 1 : CPLErr err = band->RasterIO(GF_Read, 0, 0, xlen, ylen, out.data(), xlen,
131 1 : ylen, GDT_Int8, 0, 0, nullptr);
132 1 : EXPECT_EQ(err, CE_None);
133 :
134 : // We expect the cell to be observable when the input is higher than the observable
135 : // height.
136 : std::array<int8_t, xlen * ylen> expected;
137 26 : for (size_t i = 0; i < in.size(); ++i)
138 25 : expected[i] = (in[i] >= observable[i] ? 127 : 0);
139 :
140 1 : EXPECT_EQ(expected, out);
141 : }
142 :
143 : {
144 : std::array<double, xlen * ylen> dem;
145 2 : SCOPED_TRACE("simple_height:dem");
146 2 : Options opts = stdOptions(2, 2);
147 1 : opts.outputMode = OutputMode::DEM;
148 :
149 2 : DatasetPtr output = runViewshed(in.data(), xlen, ylen, opts);
150 :
151 1 : GDALRasterBand *band = output->GetRasterBand(1);
152 1 : CPLErr err = band->RasterIO(GF_Read, 0, 0, xlen, ylen, dem.data(), xlen,
153 1 : ylen, GDT_Float64, 0, 0, nullptr);
154 1 : EXPECT_EQ(err, CE_None);
155 :
156 : // DEM values are observable values clamped at 0. Not sure why.
157 1 : std::array<double, xlen *ylen> expected = observable;
158 26 : for (double &d : expected)
159 25 : d = std::max(0.0, d);
160 :
161 : // Double equality is fine here as all the values are small integers.
162 1 : EXPECT_EQ(dem, expected);
163 : }
164 :
165 : {
166 : std::array<double, xlen * ylen> ground;
167 2 : SCOPED_TRACE("simple_height:ground");
168 2 : Options opts = stdOptions(2, 2);
169 1 : opts.outputMode = OutputMode::Ground;
170 2 : DatasetPtr output = runViewshed(in.data(), xlen, ylen, opts);
171 :
172 1 : GDALRasterBand *band = output->GetRasterBand(1);
173 1 : CPLErr err = band->RasterIO(GF_Read, 0, 0, xlen, ylen, ground.data(),
174 1 : xlen, ylen, GDT_Float64, 0, 0, nullptr);
175 1 : EXPECT_EQ(err, CE_None);
176 :
177 : std::array<double, xlen * ylen> expected;
178 26 : for (size_t i = 0; i < expected.size(); ++i)
179 25 : expected[i] = std::max(0.0, observable[i] - in[i]);
180 :
181 : // Double equality is fine here as all the values are small integers.
182 1 : EXPECT_EQ(expected, ground);
183 : }
184 1 : }
185 :
186 : // Addresses cases in #9501
187 4 : TEST(Viewshed, dem_vs_ground)
188 : {
189 : // Run gdal_viewshed on the input 8 x 1 array in both ground and dem mode and
190 : // verify the results are what are expected.
191 5 : auto run = [](const std::array<int8_t, 8> &in, Coord observer,
192 : const std::array<double, 8> &ground,
193 : const std::array<double, 8> &dem)
194 : {
195 5 : const int xlen = 8;
196 5 : const int ylen = 1;
197 :
198 : std::array<double, xlen * ylen> out;
199 10 : Options opts = stdOptions(observer);
200 5 : opts.outputMode = OutputMode::Ground;
201 :
202 : // Verify ground mode.
203 10 : DatasetPtr ds = runViewshed(in.data(), xlen, ylen, opts);
204 5 : GDALRasterBand *band = ds->GetRasterBand(1);
205 5 : CPLErr err = band->RasterIO(GF_Read, 0, 0, xlen, ylen, out.data(), xlen,
206 5 : ylen, GDT_Float64, 0, 0, nullptr);
207 5 : EXPECT_EQ(err, CE_None);
208 45 : for (size_t i = 0; i < ground.size(); ++i)
209 40 : EXPECT_DOUBLE_EQ(out[i], ground[i]);
210 :
211 : // Verify DEM mode.
212 5 : opts.outputMode = OutputMode::DEM;
213 5 : ds = runViewshed(in.data(), xlen, ylen, opts);
214 5 : band = ds->GetRasterBand(1);
215 5 : err = band->RasterIO(GF_Read, 0, 0, xlen, ylen, out.data(), xlen, ylen,
216 : GDT_Float64, 0, 0, nullptr);
217 5 : EXPECT_EQ(err, CE_None);
218 45 : for (size_t i = 0; i < dem.size(); ++i)
219 40 : EXPECT_DOUBLE_EQ(out[i], dem[i]);
220 5 : };
221 :
222 : // Input / Observer / Minimum expected above ground / Minimum expected above zero
223 1 : run({0, 0, 0, 1, 0, 0, 0, 0}, {2, 0}, {0, 0, 0, 0, 2, 3, 4, 5},
224 : {0, 0, 0, 1, 2, 3, 4, 5});
225 1 : run({1, 1, 0, 1, 0, 1, 2, 2}, {3, 0}, {0, 0, 0, 0, 0, 0, 0, 1 / 3.0},
226 : {1, 0, 0, 1, 0, 0, 1, 7 / 3.0});
227 1 : run({0, 0, 0, 1, 1, 0, 0, 0}, {0, 0},
228 : {0, 0, 0, 0, 1 / 3.0, 5 / 3.0, 6 / 3.0, 7 / 3.0},
229 : {0, 0, 0, 0, 4 / 3.0, 5 / 3.0, 6 / 3.0, 7 / 3.0});
230 1 : run({0, 0, 1, 2, 3, 4, 5, 6}, {0, 0}, {0, 0, 0, 0, 0, 0, 0, 0},
231 : {0, 0, 0, 3 / 2.0, 8 / 3.0, 15 / 4.0, 24 / 5.0, 35 / 6.0});
232 1 : run({0, 0, 1, 1, 3, 4, 5, 4}, {0, 0}, {0, 0, 0, .5, 0, 0, 0, 11 / 6.0},
233 : {0, 0, 0, 3 / 2.0, 2, 15 / 4.0, 24 / 5.0, 35 / 6.0});
234 1 : }
235 :
236 : // Test an observer to the right of the raster.
237 4 : TEST(Viewshed, oor_right)
238 : {
239 : // clang-format off
240 1 : const int xlen = 5;
241 1 : const int ylen = 3;
242 1 : std::array<int8_t, xlen * ylen> in
243 : {
244 : 1, 2, 0, 4, 1,
245 : 0, 0, 2, 1, 0,
246 : 1, 0, 0, 3, 3
247 : };
248 : // clang-format on
249 :
250 : {
251 2 : Options opts = stdOptions(6, 1);
252 1 : opts.outputMode = OutputMode::DEM;
253 2 : DatasetPtr ds = runViewshed(in.data(), xlen, ylen, opts);
254 1 : GDALRasterBand *band = ds->GetRasterBand(1);
255 : std::array<double, xlen * ylen> out;
256 1 : CPLErr err = band->RasterIO(GF_Read, 0, 0, xlen, ylen, out.data(), xlen,
257 1 : ylen, GDT_Float64, 0, 0, nullptr);
258 1 : EXPECT_EQ(err, CE_None);
259 :
260 : // clang-format off
261 1 : std::array<double, xlen * ylen> expected
262 : {
263 : 16 / 3.0, 29 / 6.0, 13 / 3.0, 1, 1,
264 : 3, 2.5, 4 / 3.0, 0, 0,
265 : 13 / 3.0, 23 / 6.0, 10 / 3.0, 3, 3
266 : };
267 : // clang-format on
268 :
269 16 : for (size_t i = 0; i < out.size(); ++i)
270 15 : EXPECT_DOUBLE_EQ(out[i], expected[i]);
271 : }
272 :
273 : {
274 2 : Options opts = stdOptions(6, 2);
275 1 : opts.outputMode = OutputMode::DEM;
276 2 : DatasetPtr ds = runViewshed(in.data(), xlen, ylen, opts);
277 1 : GDALRasterBand *band = ds->GetRasterBand(1);
278 : std::array<double, xlen * ylen> out;
279 1 : CPLErr err = band->RasterIO(GF_Read, 0, 0, xlen, ylen, out.data(), xlen,
280 1 : ylen, GDT_Float64, 0, 0, nullptr);
281 1 : EXPECT_EQ(err, CE_None);
282 :
283 : // clang-format off
284 1 : std::array<double, xlen * ylen> expected
285 : {
286 : 26 / 5.0, 17 / 4.0, 11 / 3.0, .5, 1,
287 : 6, 4.5, 3, 1.5, 0,
288 : 9, 7.5, 6, 4.5, 3
289 : };
290 : // clang-format on
291 :
292 16 : for (size_t i = 0; i < out.size(); ++i)
293 15 : EXPECT_DOUBLE_EQ(out[i], expected[i]);
294 : }
295 1 : }
296 :
297 : // Test an observer to the left of the raster.
298 4 : TEST(Viewshed, oor_left)
299 : {
300 : // clang-format off
301 1 : const int xlen = 5;
302 1 : const int ylen = 3;
303 1 : std::array<int8_t, xlen * ylen> in
304 : {
305 : 1, 2, 0, 4, 1,
306 : 0, 0, 2, 1, 0,
307 : 1, 0, 0, 3, 3
308 : };
309 : // clang-format on
310 :
311 : {
312 2 : Options opts = stdOptions(-2, 1);
313 1 : opts.outputMode = OutputMode::DEM;
314 2 : DatasetPtr ds = runViewshed(in.data(), xlen, ylen, opts);
315 1 : GDALRasterBand *band = ds->GetRasterBand(1);
316 : std::array<double, xlen * ylen> out;
317 1 : CPLErr err = band->RasterIO(GF_Read, 0, 0, xlen, ylen, out.data(), xlen,
318 1 : ylen, GDT_Float64, 0, 0, nullptr);
319 1 : EXPECT_EQ(err, CE_None);
320 :
321 : // clang-format off
322 1 : std::array<double, xlen * ylen> expected
323 : {
324 : 1, 1, 2, 2.5, 4.5,
325 : 0, 0, 0, 2.5, 3,
326 : 1, 1, 1, 1.5, 3.5
327 : };
328 : // clang-format on
329 :
330 16 : for (size_t i = 0; i < out.size(); ++i)
331 15 : EXPECT_DOUBLE_EQ(out[i], expected[i]);
332 : }
333 :
334 : {
335 2 : Options opts = stdOptions(-2, 2);
336 1 : opts.outputMode = OutputMode::DEM;
337 2 : DatasetPtr ds = runViewshed(in.data(), xlen, ylen, opts);
338 1 : GDALRasterBand *band = ds->GetRasterBand(1);
339 : std::array<double, xlen * ylen> out;
340 1 : CPLErr err = band->RasterIO(GF_Read, 0, 0, xlen, ylen, out.data(), xlen,
341 1 : ylen, GDT_Float64, 0, 0, nullptr);
342 1 : EXPECT_EQ(err, CE_None);
343 :
344 : // clang-format off
345 1 : std::array<double, xlen * ylen> expected
346 : {
347 : 1, .5, 5 / 3.0, 2.25, 4.2,
348 : 0, .5, 1, 2.5, 3.1,
349 : 1, 1.5, 2, 2.5, 3.6
350 : };
351 : // clang-format on
352 :
353 16 : for (size_t i = 0; i < out.size(); ++i)
354 15 : EXPECT_DOUBLE_EQ(out[i], expected[i]);
355 : }
356 1 : }
357 :
358 : // Test an observer above the raster
359 4 : TEST(Viewshed, oor_above)
360 : {
361 : // clang-format off
362 1 : const int xlen = 5;
363 1 : const int ylen = 3;
364 1 : std::array<int8_t, xlen * ylen> in
365 : {
366 : 1, 2, 0, 4, 1,
367 : 0, 0, 2, 1, 0,
368 : 1, 0, 0, 3, 3
369 : };
370 : // clang-format on
371 :
372 : {
373 2 : Options opts = stdOptions(2, -2);
374 1 : opts.outputMode = OutputMode::DEM;
375 2 : DatasetPtr ds = runViewshed(in.data(), xlen, ylen, opts);
376 1 : GDALRasterBand *band = ds->GetRasterBand(1);
377 : std::array<double, xlen * ylen> out;
378 1 : CPLErr err = band->RasterIO(GF_Read, 0, 0, xlen, ylen, out.data(), xlen,
379 1 : ylen, GDT_Float64, 0, 0, nullptr);
380 :
381 1 : EXPECT_EQ(err, CE_None);
382 :
383 : // clang-format off
384 1 : std::array<double, xlen * ylen> expected
385 : {
386 : 1, 2, 0, 4, 1,
387 : 2.5, 2, 0, 4, 4.5,
388 : 3, 8 / 3.0, 8 / 3.0, 14 / 3.0, 17 / 3.0
389 : };
390 : // clang-format on
391 :
392 16 : for (size_t i = 0; i < out.size(); ++i)
393 15 : EXPECT_DOUBLE_EQ(out[i], expected[i]);
394 : }
395 :
396 : {
397 2 : Options opts = stdOptions(-2, -2);
398 1 : opts.outputMode = OutputMode::DEM;
399 2 : DatasetPtr ds = runViewshed(in.data(), xlen, ylen, opts);
400 1 : GDALRasterBand *band = ds->GetRasterBand(1);
401 : std::array<double, xlen * ylen> out;
402 1 : CPLErr err = band->RasterIO(GF_Read, 0, 0, xlen, ylen, out.data(), xlen,
403 1 : ylen, GDT_Float64, 0, 0, nullptr);
404 1 : EXPECT_EQ(err, CE_None);
405 :
406 : // clang-format off
407 1 : std::array<double, xlen * ylen> expected
408 : {
409 : 1, 2, 0, 4, 1,
410 : 0, 1.5, 2.5, 1.25, 3.15,
411 : 1, 0.5, 2, 3, 2.2
412 : };
413 : // clang-format on
414 :
415 16 : for (size_t i = 0; i < out.size(); ++i)
416 15 : EXPECT_DOUBLE_EQ(out[i], expected[i]);
417 : }
418 1 : }
419 :
420 : // Test an observer below the raster
421 4 : TEST(Viewshed, oor_below)
422 : {
423 : // clang-format off
424 1 : const int xlen = 5;
425 1 : const int ylen = 3;
426 1 : std::array<int8_t, xlen * ylen> in
427 : {
428 : 1, 2, 0, 4, 1,
429 : 0, 0, 2, 1, 0,
430 : 1, 0, 0, 3, 3
431 : };
432 : // clang-format on
433 :
434 : {
435 2 : Options opts = stdOptions(2, 4);
436 1 : opts.outputMode = OutputMode::DEM;
437 2 : DatasetPtr ds = runViewshed(in.data(), xlen, ylen, opts);
438 1 : GDALRasterBand *band = ds->GetRasterBand(1);
439 : std::array<double, xlen * ylen> out;
440 1 : CPLErr err = band->RasterIO(GF_Read, 0, 0, xlen, ylen, out.data(), xlen,
441 1 : ylen, GDT_Float64, 0, 0, nullptr);
442 :
443 1 : EXPECT_EQ(err, CE_None);
444 :
445 : // clang-format off
446 1 : std::array<double, xlen * ylen> expected
447 : {
448 : 1 / 3.0, 2 / 3.0, 8 / 3.0, 11 / 3.0, 5,
449 : 0.5, 0, 0, 3, 4.5,
450 : 1, 0, 0, 3, 3
451 : };
452 : // clang-format on
453 :
454 16 : for (size_t i = 0; i < out.size(); ++i)
455 15 : EXPECT_DOUBLE_EQ(out[i], expected[i]);
456 : }
457 :
458 : {
459 2 : Options opts = stdOptions(6, 4);
460 1 : opts.outputMode = OutputMode::DEM;
461 2 : DatasetPtr ds = runViewshed(in.data(), xlen, ylen, opts);
462 1 : GDALRasterBand *band = ds->GetRasterBand(1);
463 : std::array<double, xlen * ylen> out;
464 1 : CPLErr err = band->RasterIO(GF_Read, 0, 0, xlen, ylen, out.data(), xlen,
465 1 : ylen, GDT_Float64, 0, 0, nullptr);
466 1 : EXPECT_EQ(err, CE_None);
467 :
468 : // clang-format off
469 1 : std::array<double, xlen * ylen> expected
470 : {
471 : 4.2, 6, 6, 1.5, 1,
472 : 1.35, 2.25, 4.5, 4.5, 0,
473 : 1, 0, 0, 3, 3
474 : };
475 : // clang-format on
476 :
477 16 : for (size_t i = 0; i < out.size(); ++i)
478 15 : EXPECT_DOUBLE_EQ(out[i], expected[i]);
479 : }
480 1 : }
481 :
482 : } // namespace viewshed
483 : } // namespace gdal
|
