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