Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL algorithms
4 : * Purpose: Tests for the marching squares algorithm
5 : * Author: Hugo Mercier, <hugo dot mercier at oslandia dot com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2018, Hugo Mercier, <hugo dot mercier at oslandia dot com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "gdal_unit_test.h"
14 :
15 : #include "gdal_alg.h"
16 :
17 : #include "marching_squares/level_generator.h"
18 : #include "marching_squares/segment_merger.h"
19 : #include "marching_squares/contour_generator.h"
20 :
21 : #include <limits>
22 :
23 : #include "gtest_include.h"
24 :
25 : namespace marching_squares
26 : {
27 : class TestRingAppender
28 : {
29 : public:
30 : struct Point
31 : {
32 172 : Point(double xx, double yy) : x(xx), y(yy)
33 : {
34 172 : }
35 :
36 : double x;
37 : double y;
38 :
39 : bool operator<(const Point &b) const
40 : {
41 : return x == b.x ? y < b.y : x < b.x;
42 : }
43 :
44 169 : bool operator==(const Point &b) const
45 : {
46 169 : return std::fabs(x - b.x) < 0.001 && std::fabs(y - b.y) < 0.001;
47 : }
48 :
49 169 : bool operator!=(const Point &b) const
50 : {
51 169 : return !(*this == b);
52 : }
53 : };
54 :
55 8 : void addLine(double level, LineString &ls, bool /* closed */)
56 : {
57 8 : auto &v = points_[level];
58 16 : std::vector<Point> ring;
59 98 : for (const auto &pt : ls)
60 : {
61 90 : ring.push_back(Point(pt.x, pt.y));
62 : }
63 8 : v.push_back(ring);
64 8 : }
65 :
66 8 : bool hasRing(double level, const std::vector<Point> &other) const
67 : {
68 8 : auto it = points_.find(level);
69 8 : if (it == points_.end())
70 : {
71 0 : return false;
72 : }
73 :
74 8 : const auto &rings = it->second;
75 9 : for (const auto &ring : rings)
76 : {
77 9 : if (ringEquals_(ring, other))
78 : {
79 8 : return true;
80 : }
81 : else
82 : {
83 : // test also the reverse ring
84 5 : auto rev = other;
85 5 : std::reverse(rev.begin(), rev.end());
86 5 : if (ringEquals_(ring, rev))
87 : {
88 4 : return true;
89 : }
90 : }
91 : }
92 0 : return false;
93 : }
94 :
95 : void out(std::ostream &o, double level)
96 : {
97 : for (const auto &p : points_[level])
98 : {
99 : out_(o, p);
100 : }
101 : }
102 :
103 : private:
104 : // level -> vector of rings
105 : std::map<double, std::vector<std::vector<Point>>> points_;
106 :
107 14 : bool ringEquals_(const std::vector<Point> &aRing,
108 : const std::vector<Point> &bRing) const
109 : {
110 14 : if (aRing.size() - 1 != bRing.size())
111 : {
112 0 : return false;
113 : }
114 :
115 : // rings do not really have a "first" point, but since
116 : // we represent them with a vector, we need to find a common "first"
117 : // point
118 14 : Point pfirst = aRing[0];
119 14 : size_t offset = 0;
120 73 : while (offset < bRing.size() && pfirst != bRing[offset])
121 59 : offset++;
122 14 : if (offset >= bRing.size())
123 : {
124 : // can't find a common point
125 2 : return false;
126 : }
127 : // now compare each point of the two rings
128 106 : for (size_t i = 0; i < aRing.size(); i++)
129 : {
130 98 : const Point &p2 = bRing[(i + offset) % bRing.size()];
131 98 : if (aRing[i] != p2)
132 : {
133 4 : return false;
134 : }
135 : }
136 8 : return true;
137 : }
138 :
139 : void out_(std::ostream &o, const std::vector<Point> &points) const
140 : {
141 : o << "{ ";
142 : for (const auto &pt : points)
143 : {
144 : o << "{" << pt.x << "," << pt.y << "}, ";
145 : }
146 : o << "}, ";
147 : }
148 : };
149 : } // namespace marching_squares
150 :
151 : namespace
152 : {
153 : using namespace marching_squares;
154 :
155 : // Common fixture with test data
156 : struct test_ms_contour : public ::testing::Test
157 : {
158 : };
159 :
160 : // Dummy test
161 4 : TEST_F(test_ms_contour, dummy)
162 : {
163 : // one pixel
164 2 : std::vector<double> data = {2.0};
165 2 : TestRingAppender w;
166 : {
167 : IntervalLevelRangeIterator levels(
168 1 : 0.0, 10.0, -std::numeric_limits<double>::infinity());
169 : SegmentMerger<TestRingAppender, IntervalLevelRangeIterator> writer(
170 2 : w, levels, /* polygonize */ true);
171 : ContourGenerator<decltype(writer), IntervalLevelRangeIterator> cg(
172 2 : 1, 1, /* hasNoData */ false, NaN, writer, levels);
173 1 : cg.feedLine(&data[0]);
174 :
175 : // "Polygon ring"
176 1 : EXPECT_TRUE(w.hasRing(10.0, {{0.0, 0.0},
177 : {0.5, 0.0},
178 : {1.0, 0.0},
179 : {1.0, 0.5},
180 : {1.0, 1.0},
181 : {0.5, 1.0},
182 : {0.0, 1.0},
183 : {0.0, 0.5}}));
184 : }
185 1 : }
186 :
187 4 : TEST_F(test_ms_contour, two_pixels)
188 : {
189 : // two pixels
190 : // 10 7
191 : // levels = 8
192 2 : std::vector<double> data = {10.0, 7.0};
193 2 : TestRingAppender w;
194 :
195 : {
196 : IntervalLevelRangeIterator levels(
197 1 : 8.0, 10.0, -std::numeric_limits<double>::infinity());
198 : SegmentMerger<TestRingAppender, IntervalLevelRangeIterator> writer(
199 2 : w, levels, /* polygonize */ true);
200 : ContourGenerator<decltype(writer), IntervalLevelRangeIterator> cg(
201 2 : 2, 1, /* hasNoData */ false, NaN, writer, levels);
202 1 : cg.feedLine(&data[0]);
203 :
204 : // "Polygon #0"
205 1 : EXPECT_TRUE(w.hasRing(8.0, {{1.166, 0.0},
206 : {1.5, 0.0},
207 : {2.0, 0.0},
208 : {2.0, 0.5},
209 : {2.0, 1.0},
210 : {1.5, 1.0},
211 : {1.166, 1.0},
212 : {1.166, 0.5}}));
213 : // "Polygon #1"
214 1 : EXPECT_TRUE(w.hasRing(18.0, {{1.166, 0.0},
215 : {1.0, 0.0},
216 : {0.5, 0.0},
217 : {0.0, 0.0},
218 : {0.0, 0.5},
219 : {0.0, 1.0},
220 : {0.5, 1.0},
221 : {1.0, 1.0},
222 : {1.166, 1.0},
223 : {1.166, 0.5}}));
224 : }
225 1 : }
226 :
227 4 : TEST_F(test_ms_contour, four_pixels)
228 : {
229 : // four pixels
230 : // 10 7
231 : // 4 5
232 : // levels = 8
233 : // pixels
234 : // +-----+-----+-----+-----+
235 : // | | | | |
236 : // | NaN | NaN | NaN | NaN |
237 : // | | | | |
238 : // +-----+-----+-----+-----+
239 : // | | | | |
240 : // | NaN | 10 | 7 | NaN |
241 : // | | | | |
242 : // +-----+-----+-----+-----+
243 : // | | | | |
244 : // | NaN | 4 | 5 | NaN |
245 : // | | | | |
246 : // +-----+-----+-----+-----+
247 : // | | | | |
248 : // | NaN | NaN | NaN | NaN |
249 : // | | | | |
250 : // +-----+-----+-----+-----+
251 : //
252 : // squares
253 : // +-----+-----+-----+-----+
254 : // |NaN | NaN | NaN | NaN |
255 : // | +.....+.....+.....+ |
256 : // | : | : | : | : |
257 : // +--:--+--:--+--:--+--:--+
258 : // | : |10: | 7: |NaN |
259 : // NaN+.....+.....+.....+ |
260 : // | : | : | : | : |
261 : // +--:--+--:--+--:--+--:--+
262 : // | : | 4: | 5: |NaN |
263 : // NaN+.....+.....+.....+ |
264 : // | : | : | : | : |
265 : // +--:--+--:--+--:--+--:--+
266 : // | : | : | : | : |
267 : // | +.....+.....+.....+ |
268 : // | NaN | NaN | NaN | NaN |
269 : // +-----+-----+-----+-----+
270 : //
271 : // subsquares
272 : // legend:
273 : // : contour
274 : // = border (level 8)
275 : // # border (level 18)
276 : //
277 : // NaN NaN NaN
278 : // +------------------+------------------+------------------+
279 : // | | | |
280 : // | (0,0) | (1,0) | (2,0) |
281 : // | 10 10| 8.5 7| 7 |
282 : // | +#########+########+###o=====+========++ |
283 : // | # | | : | || |
284 : // | # | | : | || |
285 : // | # | | : | || |
286 : // +--------+---------+--------+---o-----+--------++--------+
287 : // |NaN 10# 10| ........: 7| 7 || NaN|
288 : // | o.........o..: | || |
289 : // | || | | || |
290 : // | 7++---------+ 7 6 +--------++ |
291 : // | || | | || |
292 : // | || | | || |
293 : // | || | 4.5 | || |
294 : // +-------++---------+--------+---------+--------++--------+
295 : // |NaN 4|| 4 | | 5| 5 || NaN|
296 : // | || | | | || |
297 : // | || | | | || |
298 : // | ++=========+========+=========+========++ |
299 : // | 4 4 | 4.5 5| 5 |
300 : // | (0,2) | (1,2) | (2,2) |
301 : // | | | |
302 : // +------------------+------------------+------------------+
303 : // NaN NaN NaN NaN
304 :
305 2 : std::vector<double> data = {10.0, 7.0, 4.0, 5.0};
306 2 : TestRingAppender w;
307 :
308 : {
309 : IntervalLevelRangeIterator levels(
310 1 : 8.0, 10.0, -std::numeric_limits<double>::infinity());
311 : SegmentMerger<TestRingAppender, IntervalLevelRangeIterator> writer(
312 2 : w, levels, /* polygonize */ true);
313 : ContourGenerator<decltype(writer), IntervalLevelRangeIterator> cg(
314 2 : 2, 2, /* hasNoData */ false, NaN, writer, levels);
315 1 : cg.feedLine(&data[0]);
316 1 : cg.feedLine(&data[2]);
317 :
318 : // "Polygon #0"
319 1 : EXPECT_TRUE(w.hasRing(8.0, {{2.0, 0.0},
320 : {2.0, 0.5},
321 : {2.0, 1.0},
322 : {2.0, 1.5},
323 : {2.0, 2.0},
324 : {1.5, 2.0},
325 : {1.0, 2.0},
326 : {0.5, 2.0},
327 : {0.0, 2.0},
328 : {0.0, 1.5},
329 : {0.0, 1.0},
330 : {0.0, 0.833},
331 : {0.5, 0.833},
332 : {1.167, 0.5},
333 : {1.167, 0.0},
334 : {1.5, 0.0}}));
335 : // "Polygon #1"
336 1 : EXPECT_TRUE(w.hasRing(18.0, {{0.0, 0.0},
337 : {0.5, 0.0},
338 : {1.0, 0.0},
339 : {1.167, 0.0},
340 : {1.167, 0.5},
341 : {0.5, 0.833},
342 : {0, 0.833},
343 : {0.0, 0.5}}));
344 : }
345 1 : }
346 :
347 4 : TEST_F(test_ms_contour, saddle_point)
348 : {
349 : // four pixels
350 : // two rings
351 : // with a saddle point
352 : // 5 10
353 : // 10 5
354 : // levels = 8
355 : // pixels
356 : // +-----+-----+-----+-----+
357 : // | | | | |
358 : // | NaN | NaN | NaN | NaN |
359 : // | | | | |
360 : // +-----+-----+-----+-----+
361 : // | | | | |
362 : // | NaN | 5 | 10 | NaN |
363 : // | | | | |
364 : // +-----+-----+-----+-----+
365 : // | | | | |
366 : // | NaN | 10 | 5 | NaN |
367 : // | | | | |
368 : // +-----+-----+-----+-----+
369 : // | | | | |
370 : // | NaN | NaN | NaN | NaN |
371 : // | | | | |
372 : // +-----+-----+-----+-----+
373 : //
374 : // squares
375 : // +-----+-----+-----+-----+
376 : // |NaN | NaN | NaN | NaN |
377 : // | +.....+.....+.....+ |
378 : // | : | : | : | : |
379 : // +--:--+--:--+--:--+--:--+
380 : // | : | 5: |10: |NaN |
381 : // NaN+.....+.....+.....+ |
382 : // | : | : | : | : |
383 : // +--:--+--:--+--:--+--:--+
384 : // | : |10: | 5: |NaN |
385 : // NaN+.....+.....+.....+ |
386 : // | : | : | : | : |
387 : // +--:--+--:--+--:--+--:--+
388 : // | : | : | : | : |
389 : // | +.....+.....+.....+ |
390 : // | NaN | NaN | NaN | NaN |
391 : // +-----+-----+-----+-----+
392 : //
393 : // subsquares
394 : // legend:
395 : // : contour
396 : // # border (level 8)
397 : // = border (level 18)
398 : //
399 : // NaN NaN NaN
400 : // +------------------+------------------+------------------+
401 : // | | | |
402 : // | (0,0) | (1,0) | (2,0) |
403 : // | 5 5| 7.5 10| 10 |
404 : // | +#########+########+###o=====+========++ |
405 : // | # | | : | || |
406 : // | # | | : | || |
407 : // | # | | : | || |
408 : // +--------+---------+--------+---o-----+--------++--------+
409 : // |NaN 5 # 5| \ 10| 10|| NaN|
410 : // | # | \___o........o |
411 : // | # | | # |
412 : // | 7.5++---------+7.5 7.5+--------+ |
413 : // | # | | # |
414 : // | o.........o\_ | # |
415 : // | || | \_ 7.5 | # |
416 : // +-------++---------+----\o--+---------+--------+---------+
417 : // |NaN 10|| 10| : | 5| 5 # NaN|
418 : // | || | : | | # |
419 : // | || | : | | # |
420 : // | ++=========+=====o##+#########+########+ |
421 : // | 10 10| 7.5 5| 5 |
422 : // | (0,2) | (1,2) | (2,2) |
423 : // | | | |
424 : // +------------------+------------------+------------------+
425 : // NaN NaN NaN NaN
426 :
427 2 : std::vector<double> data = {5.0, 10.0, 10.0, 5.0};
428 2 : TestRingAppender w;
429 :
430 : {
431 : IntervalLevelRangeIterator levels(
432 1 : 8.0, 10.0, -std::numeric_limits<double>::infinity());
433 : SegmentMerger<TestRingAppender, IntervalLevelRangeIterator> writer(
434 2 : w, levels, /* polygonize */ true);
435 : ContourGenerator<decltype(writer), IntervalLevelRangeIterator> cg(
436 2 : 2, 2, /* hasNoData */ false, NaN, writer, levels);
437 1 : cg.feedLine(&data[0]);
438 1 : cg.feedLine(&data[2]);
439 :
440 : // "Polygon #0"
441 1 : EXPECT_TRUE(w.hasRing(8.0, {{1.5, 2},
442 : {2, 2},
443 : {2, 1.5},
444 : {2, 1},
445 : {2, 0.9},
446 : {1.5, 0.9},
447 : {1.1, 0.5},
448 : {1.1, 0},
449 : {1, 0},
450 : {0.5, 0},
451 : {0, 0},
452 : {0, 0.5},
453 : {0, 1},
454 : {0, 1.1},
455 : {0.5, 1.1},
456 : {0.9, 1.5},
457 : {0.9, 2},
458 : {1, 2}}));
459 : // "Polygon #1, Ring #0"
460 1 : EXPECT_TRUE(w.hasRing(18.0, {{2, 0.9},
461 : {2, 0.5},
462 : {2, 0},
463 : {1.5, 0},
464 : {1.1, 0},
465 : {1.1, 0.5},
466 : {1.5, 0.9}}));
467 : // "Polygon #1, Ring #1"
468 1 : EXPECT_TRUE(w.hasRing(18.0, {{0.9, 1.5},
469 : {0.5, 1.1},
470 : {0, 1.1},
471 : {0, 1.5},
472 : {0, 2},
473 : {0.5, 2},
474 : {0.9, 2}}));
475 : }
476 1 : }
477 : } // namespace
|
