Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Marching square algorithm
4 : * Purpose: Core algorithm implementation for contour line generation.
5 : * Author: Oslandia <infos at oslandia dot com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2018, Oslandia <infos at oslandia dot com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 : #ifndef MARCHING_SQUARE_SQUARE_H
13 : #define MARCHING_SQUARE_SQUARE_H
14 :
15 : #include <algorithm>
16 : #include <cassert>
17 : #include <cmath>
18 : #include <cstdint>
19 : #include <math.h>
20 : #include "utility.h"
21 : #include "point.h"
22 :
23 : namespace marching_squares
24 : {
25 :
26 : struct Square
27 : {
28 : // Bit flags to determine borders around pixel
29 : static const uint8_t NO_BORDER = 0; // 0000 0000
30 : static const uint8_t LEFT_BORDER = 1 << 0; // 0000 0001
31 : static const uint8_t LOWER_BORDER = 1 << 1; // 0000 0010
32 : static const uint8_t RIGHT_BORDER = 1 << 2; // 0000 0100
33 : static const uint8_t UPPER_BORDER = 1 << 3; // 0000 1000
34 :
35 : // Bit flags for marching square case
36 : static const uint8_t ALL_LOW = 0; // 0000 0000
37 : static const uint8_t UPPER_LEFT = 1 << 0; // 0000 0001
38 : static const uint8_t LOWER_LEFT = 1 << 1; // 0000 0010
39 : static const uint8_t LOWER_RIGHT = 1 << 2; // 0000 0100
40 : static const uint8_t UPPER_RIGHT = 1 << 3; // 0000 1000
41 : static const uint8_t ALL_HIGH =
42 : UPPER_LEFT | LOWER_LEFT | LOWER_RIGHT | UPPER_RIGHT; // 0000 1111
43 : static const uint8_t SADDLE_NW = UPPER_LEFT | LOWER_RIGHT; // 0000 0101
44 : static const uint8_t SADDLE_NE = UPPER_RIGHT | LOWER_LEFT; // 0000 1010
45 :
46 : typedef std::pair<Point, Point> Segment;
47 : typedef std::pair<ValuedPoint, ValuedPoint> ValuedSegment;
48 :
49 : //
50 : // An array of segments, at most 3 segments
51 : struct Segments
52 : {
53 24 : Segments() : sz_(0)
54 : {
55 6 : }
56 :
57 339120 : Segments(const Segment &first) : sz_(1), segs_()
58 : {
59 84780 : segs_[0] = first;
60 84780 : }
61 :
62 10156 : Segments(const Segment &first, const Segment &second) : sz_(2), segs_()
63 : {
64 2539 : segs_[0] = first;
65 2539 : segs_[1] = second;
66 2539 : }
67 :
68 : Segments(const Segment &first, const Segment &second,
69 : const Segment &third)
70 : : sz_(3), segs_()
71 : {
72 : segs_[0] = first;
73 : segs_[1] = second;
74 : segs_[2] = third;
75 : }
76 :
77 177177 : std::size_t size() const
78 : {
79 177177 : return sz_;
80 : }
81 :
82 89864 : const Segment &operator[](std::size_t idx) const
83 : {
84 89864 : assert(idx < sz_);
85 89864 : return segs_[idx];
86 : }
87 :
88 : private:
89 : const std::size_t sz_;
90 : /* const */ Segment segs_[3];
91 : };
92 :
93 824024 : Square(const ValuedPoint &upperLeft_, const ValuedPoint &upperRight_,
94 : const ValuedPoint &lowerLeft_, const ValuedPoint &lowerRight_,
95 : uint8_t borders_ = NO_BORDER, bool split_ = false)
96 824024 : : upperLeft(upperLeft_), lowerLeft(lowerLeft_), lowerRight(lowerRight_),
97 : upperRight(upperRight_),
98 824024 : nanCount((std::isnan(upperLeft.value) ? 1 : 0) +
99 824024 : (std::isnan(upperRight.value) ? 1 : 0) +
100 824024 : (std::isnan(lowerLeft.value) ? 1 : 0) +
101 824024 : (std::isnan(lowerRight.value) ? 1 : 0)),
102 824024 : borders(borders_), split(split_)
103 : {
104 824024 : assert(upperLeft.y == upperRight.y);
105 824024 : assert(lowerLeft.y == lowerRight.y);
106 824024 : assert(lowerLeft.x == upperLeft.x);
107 824024 : assert(lowerRight.x == upperRight.x);
108 824024 : assert(!split || nanCount == 0);
109 824024 : }
110 :
111 10278 : Square upperLeftSquare() const
112 : {
113 10278 : assert(!std::isnan(upperLeft.value));
114 : return Square(
115 20556 : upperLeft, upperCenter(), leftCenter(), center(),
116 20556 : (std::isnan(upperRight.value) ? RIGHT_BORDER : NO_BORDER) |
117 10278 : (std::isnan(lowerLeft.value) ? LOWER_BORDER : NO_BORDER),
118 20556 : true);
119 : }
120 :
121 10277 : Square lowerLeftSquare() const
122 : {
123 10277 : assert(!std::isnan(lowerLeft.value));
124 : return Square(
125 20554 : leftCenter(), center(), lowerLeft, lowerCenter(),
126 20554 : (std::isnan(lowerRight.value) ? RIGHT_BORDER : NO_BORDER) |
127 10277 : (std::isnan(upperLeft.value) ? UPPER_BORDER : NO_BORDER),
128 20554 : true);
129 : }
130 :
131 10279 : Square lowerRightSquare() const
132 : {
133 10279 : assert(!std::isnan(lowerRight.value));
134 : return Square(
135 20558 : center(), rightCenter(), lowerCenter(), lowerRight,
136 20558 : (std::isnan(lowerLeft.value) ? LEFT_BORDER : NO_BORDER) |
137 10279 : (std::isnan(upperRight.value) ? UPPER_BORDER : NO_BORDER),
138 20558 : true);
139 : }
140 :
141 10275 : Square upperRightSquare() const
142 : {
143 10275 : assert(!std::isnan(upperRight.value));
144 : return Square(
145 20550 : upperCenter(), upperRight, center(), rightCenter(),
146 20550 : (std::isnan(lowerRight.value) ? LOWER_BORDER : NO_BORDER) |
147 10275 : (std::isnan(upperLeft.value) ? LEFT_BORDER : NO_BORDER),
148 20550 : true);
149 : }
150 :
151 803345 : double maxValue() const
152 : {
153 803345 : assert(nanCount == 0);
154 803345 : return std::max(std::max(upperLeft.value, upperRight.value),
155 803345 : std::max(lowerLeft.value, lowerRight.value));
156 : }
157 :
158 803345 : double minValue() const
159 : {
160 803345 : assert(nanCount == 0);
161 803345 : return std::min(std::min(upperLeft.value, upperRight.value),
162 803345 : std::min(lowerLeft.value, lowerRight.value));
163 : }
164 :
165 25913 : ValuedSegment segment(uint8_t border) const
166 : {
167 25913 : switch (border)
168 : {
169 6477 : case LEFT_BORDER:
170 6477 : return ValuedSegment(upperLeft, lowerLeft);
171 6480 : case LOWER_BORDER:
172 6480 : return ValuedSegment(lowerLeft, lowerRight);
173 6475 : case RIGHT_BORDER:
174 6475 : return ValuedSegment(lowerRight, upperRight);
175 6481 : case UPPER_BORDER:
176 6481 : return ValuedSegment(upperRight, upperLeft);
177 : }
178 0 : assert(false);
179 : return ValuedSegment(upperLeft, upperLeft);
180 : }
181 :
182 : // returns segments of contour
183 : //
184 : // segments are oriented:
185 : // - they form a vector from their first point to their second point.
186 : // - when looking at the vector upward, values greater than the level are
187 : // on the right
188 : //
189 : // ^
190 : // - | +
191 87325 : Segments segments(double level, double minLevel) const
192 : {
193 87325 : switch (marchingCase(level, minLevel))
194 : {
195 1 : case (ALL_LOW):
196 : // debug("ALL_LOW");
197 1 : return Segments();
198 5 : case (ALL_HIGH):
199 : // debug("ALL_HIGH");
200 5 : return Segments();
201 8970 : case (UPPER_LEFT):
202 : // debug("UPPER_LEFT");
203 : return Segments(
204 8970 : Segment(interpolate(UPPER_BORDER, level, minLevel),
205 26910 : interpolate(LEFT_BORDER, level, minLevel)));
206 5557 : case (LOWER_LEFT):
207 : // debug("LOWER_LEFT");
208 : return Segments(
209 5557 : Segment(interpolate(LEFT_BORDER, level, minLevel),
210 16671 : interpolate(LOWER_BORDER, level, minLevel)));
211 3491 : case (LOWER_RIGHT):
212 : // debug("LOWER_RIGHT");
213 : return Segments(
214 3491 : Segment(interpolate(LOWER_BORDER, level, minLevel),
215 10473 : interpolate(RIGHT_BORDER, level, minLevel)));
216 4192 : case (UPPER_RIGHT):
217 : // debug("UPPER_RIGHT");
218 : return Segments(
219 4192 : Segment(interpolate(RIGHT_BORDER, level, minLevel),
220 12576 : interpolate(UPPER_BORDER, level, minLevel)));
221 8986 : case (UPPER_LEFT | LOWER_LEFT):
222 : // debug("UPPER_LEFT | LOWER_LEFT");
223 : return Segments(
224 8986 : Segment(interpolate(UPPER_BORDER, level, minLevel),
225 26958 : interpolate(LOWER_BORDER, level, minLevel)));
226 12124 : case (LOWER_LEFT | LOWER_RIGHT):
227 : // debug("LOWER_LEFT | LOWER_RIGHT");
228 : return Segments(
229 12124 : Segment(interpolate(LEFT_BORDER, level, minLevel),
230 36372 : interpolate(RIGHT_BORDER, level, minLevel)));
231 5830 : case (LOWER_RIGHT | UPPER_RIGHT):
232 : // debug("LOWER_RIGHT | UPPER_RIGHT");
233 : return Segments(
234 5830 : Segment(interpolate(LOWER_BORDER, level, minLevel),
235 17490 : interpolate(UPPER_BORDER, level, minLevel)));
236 12351 : case (UPPER_RIGHT | UPPER_LEFT):
237 : // debug("UPPER_RIGHT | UPPER_LEFT");
238 : return Segments(
239 12351 : Segment(interpolate(RIGHT_BORDER, level, minLevel),
240 37053 : interpolate(LEFT_BORDER, level, minLevel)));
241 3877 : case (ALL_HIGH & ~UPPER_LEFT):
242 : // debug("ALL_HIGH & ~UPPER_LEFT");
243 : return Segments(
244 3877 : Segment(interpolate(LEFT_BORDER, level, minLevel),
245 11631 : interpolate(UPPER_BORDER, level, minLevel)));
246 4540 : case (ALL_HIGH & ~LOWER_LEFT):
247 : // debug("ALL_HIGH & ~LOWER_LEFT");
248 : return Segments(
249 4540 : Segment(interpolate(LOWER_BORDER, level, minLevel),
250 13620 : interpolate(LEFT_BORDER, level, minLevel)));
251 9022 : case (ALL_HIGH & ~LOWER_RIGHT):
252 : // debug("ALL_HIGH & ~LOWER_RIGHT");
253 : return Segments(
254 9022 : Segment(interpolate(RIGHT_BORDER, level, minLevel),
255 27066 : interpolate(LOWER_BORDER, level, minLevel)));
256 5840 : case (ALL_HIGH & ~UPPER_RIGHT):
257 : // debug("ALL_HIGH & ~UPPER_RIGHT");
258 : return Segments(
259 5840 : Segment(interpolate(UPPER_BORDER, level, minLevel),
260 17520 : interpolate(RIGHT_BORDER, level, minLevel)));
261 2539 : case (SADDLE_NE):
262 : case (SADDLE_NW):
263 : // From the two possible saddle configurations, we always return
264 : // the same one.
265 : //
266 : // The classical marching square algorithm says the ambiguity
267 : // should be resolved between the two possible configurations by
268 : // looking at the value of the center point. But in certain
269 : // cases, this may lead to line contours from different levels
270 : // that cross each other and then gives invalid polygons.
271 : //
272 : // Arbitrarily choosing one of the two possible configurations
273 : // is not really that worse than deciding based on the center
274 : // point.
275 : return Segments(
276 2539 : Segment(interpolate(LEFT_BORDER, level, minLevel),
277 2539 : interpolate(LOWER_BORDER, level, minLevel)),
278 2539 : Segment(interpolate(RIGHT_BORDER, level, minLevel),
279 7617 : interpolate(UPPER_BORDER, level, minLevel)));
280 : }
281 0 : assert(false);
282 : return Segments();
283 : }
284 :
285 : template <typename Writer, typename LevelGenerator>
286 824012 : void process(const LevelGenerator &levelGenerator, Writer &writer) const
287 : {
288 824012 : if (nanCount == 4) // nothing to do
289 20667 : return;
290 :
291 824012 : if (nanCount) // split in 4
292 : {
293 20667 : if (!std::isnan(upperLeft.value))
294 10275 : upperLeftSquare().process(levelGenerator, writer);
295 20667 : if (!std::isnan(upperRight.value))
296 10275 : upperRightSquare().process(levelGenerator, writer);
297 20667 : if (!std::isnan(lowerLeft.value))
298 10275 : lowerLeftSquare().process(levelGenerator, writer);
299 20667 : if (!std::isnan(lowerRight.value))
300 10278 : lowerRightSquare().process(levelGenerator, writer);
301 20667 : return;
302 : }
303 :
304 803345 : if (writer.polygonize && borders)
305 : {
306 128770 : for (uint8_t border :
307 : {UPPER_BORDER, LEFT_BORDER, RIGHT_BORDER, LOWER_BORDER})
308 : {
309 : // bitwise AND to test which borders we have on the square
310 103016 : if ((border & borders) == 0)
311 77103 : continue;
312 :
313 : // convention: for a level = L, store borders for the previous
314 : // level up to (and including) L in the border of level "L". For
315 : // fixed sets of level, this means there is an "Inf" slot for
316 : // borders of the highest level
317 25913 : const ValuedSegment s(segment(border));
318 :
319 25913 : Point lastPoint(s.first.x, s.first.y);
320 25913 : Point endPoint(s.second.x, s.second.y);
321 25913 : if (s.first.value > s.second.value)
322 1598 : std::swap(lastPoint, endPoint);
323 25913 : bool reverse =
324 27063 : (s.first.value > s.second.value) &&
325 1150 : ((border == UPPER_BORDER) || (border == LEFT_BORDER));
326 :
327 25913 : auto levelIt =
328 25913 : levelGenerator.range(s.first.value, s.second.value);
329 :
330 25913 : auto it = levelIt.begin(); // reused after the for
331 26043 : for (; it != levelIt.end(); ++it)
332 : {
333 130 : const int levelIdx = (*it).first;
334 130 : const double level = (*it).second;
335 :
336 130 : const Point nextPoint(
337 : interpolate(border, level, levelGenerator.minLevel()));
338 130 : if (reverse)
339 36 : writer.addBorderSegment(levelIdx, nextPoint, lastPoint);
340 : else
341 94 : writer.addBorderSegment(levelIdx, lastPoint, nextPoint);
342 130 : lastPoint = nextPoint;
343 : }
344 : // last level (past the end)
345 25913 : if (reverse)
346 1004 : writer.addBorderSegment((*it).first, endPoint, lastPoint);
347 : else
348 24909 : writer.addBorderSegment((*it).first, lastPoint, endPoint);
349 : }
350 : }
351 :
352 803345 : auto range = levelGenerator.range(minValue(), maxValue());
353 803345 : auto it = range.begin();
354 803345 : auto itEnd = range.end();
355 803345 : auto next = range.begin();
356 803345 : ++next;
357 :
358 890658 : for (; it != itEnd; ++it, ++next)
359 : {
360 87313 : const int levelIdx = (*it).first;
361 87313 : const double level = (*it).second;
362 :
363 87313 : const Segments segments_ =
364 : segments(level, levelGenerator.minLevel());
365 :
366 177165 : for (std::size_t i = 0; i < segments_.size(); i++)
367 : {
368 89852 : const Segment &s = segments_[i];
369 89852 : writer.addSegment(levelIdx, s.first, s.second);
370 :
371 89852 : if (writer.polygonize)
372 : {
373 : // the contour is used in the polygon of higher level as
374 : // well
375 : //
376 : // TODO: copying the segment to the higher level is easy,
377 : // but it involves too much memory. We should reuse segment
378 : // contours when constructing polygon rings.
379 16073 : writer.addSegment((*next).first, s.first, s.second);
380 : }
381 : }
382 : }
383 : }
384 :
385 : const ValuedPoint upperLeft;
386 : const ValuedPoint lowerLeft;
387 : const ValuedPoint lowerRight;
388 : const ValuedPoint upperRight;
389 : const int nanCount;
390 : const uint8_t borders;
391 : const bool split;
392 :
393 : private:
394 41109 : ValuedPoint center() const
395 : {
396 : return ValuedPoint(
397 41109 : .5 * (upperLeft.x + lowerRight.x),
398 41109 : .5 * (upperLeft.y + lowerRight.y),
399 41109 : ((std::isnan(lowerLeft.value) ? 0 : lowerLeft.value) +
400 41109 : (std::isnan(upperLeft.value) ? 0 : upperLeft.value) +
401 41109 : (std::isnan(lowerRight.value) ? 0 : lowerRight.value) +
402 41109 : (std::isnan(upperRight.value) ? 0 : upperRight.value)) /
403 41109 : (4 - nanCount));
404 : }
405 :
406 20555 : ValuedPoint leftCenter() const
407 : {
408 : return ValuedPoint(
409 20555 : upperLeft.x, .5 * (upperLeft.y + lowerLeft.y),
410 20555 : std::isnan(upperLeft.value)
411 : ? lowerLeft.value
412 15388 : : (std::isnan(lowerLeft.value)
413 15388 : ? upperLeft.value
414 35943 : : .5 * (upperLeft.value + lowerLeft.value)));
415 : }
416 :
417 20556 : ValuedPoint lowerCenter() const
418 : {
419 : return ValuedPoint(
420 20556 : .5 * (lowerLeft.x + lowerRight.x), lowerLeft.y,
421 20556 : std::isnan(lowerRight.value)
422 : ? lowerLeft.value
423 15390 : : (std::isnan(lowerLeft.value)
424 15390 : ? lowerRight.value
425 35946 : : .5 * (lowerRight.value + lowerLeft.value)));
426 : }
427 :
428 20554 : ValuedPoint rightCenter() const
429 : {
430 : return ValuedPoint(
431 20554 : upperRight.x, .5 * (upperRight.y + lowerRight.y),
432 20554 : std::isnan(lowerRight.value)
433 : ? upperRight.value
434 15388 : : (std::isnan(upperRight.value)
435 15388 : ? lowerRight.value
436 35942 : : .5 * (lowerRight.value + upperRight.value)));
437 : }
438 :
439 20553 : ValuedPoint upperCenter() const
440 : {
441 : return ValuedPoint(
442 20553 : .5 * (upperLeft.x + upperRight.x), upperLeft.y,
443 20553 : std::isnan(upperLeft.value)
444 : ? upperRight.value
445 15387 : : (std::isnan(upperRight.value)
446 15387 : ? upperLeft.value
447 35940 : : .5 * (upperLeft.value + upperRight.value)));
448 : }
449 :
450 87325 : uint8_t marchingCase(double level, double minLevel) const
451 : {
452 87325 : return (level < fudge(upperLeft.value, minLevel, level) ? UPPER_LEFT
453 174650 : : ALL_LOW) |
454 87325 : (level < fudge(lowerLeft.value, minLevel, level) ? LOWER_LEFT
455 87325 : : ALL_LOW) |
456 87325 : (level < fudge(lowerRight.value, minLevel, level) ? LOWER_RIGHT
457 87325 : : ALL_LOW) |
458 87325 : (level < fudge(upperRight.value, minLevel, level) ? UPPER_RIGHT
459 87325 : : ALL_LOW);
460 : }
461 :
462 179846 : static double interpolate_(double level, double x1, double x2, double y1,
463 : double y2, bool need_split, double minLevel)
464 : {
465 179846 : if (need_split)
466 : {
467 : // The two cases are here to avoid numerical roundup errors, for two
468 : // points, we always compute the same interpolation. This condition
469 : // is ensured by the order left->right bottom->top in interpole
470 : // calls
471 : //
472 : // To obtain the same value for border (split) and non-border
473 : // element, we take the middle value and interpolate from this to
474 : // the end
475 174708 : const double xm = .5 * (x1 + x2);
476 174708 : const double ym = .5 * (y1 + y2);
477 174708 : const double fy1 = fudge(y1, minLevel, level);
478 174708 : const double fym = fudge(ym, minLevel, level);
479 174708 : if ((fy1 < level && level < fym) || (fy1 > level && level > fym))
480 : {
481 85785 : x2 = xm;
482 85785 : y2 = ym;
483 : }
484 : else
485 : {
486 88923 : x1 = xm;
487 88923 : y1 = ym;
488 : }
489 : }
490 179846 : const double fy1 = fudge(y1, minLevel, level);
491 179846 : const double ratio = (level - fy1) / (fudge(y2, minLevel, level) - fy1);
492 179846 : return x1 * (1. - ratio) + x2 * ratio;
493 : }
494 :
495 179846 : Point interpolate(uint8_t border, double level, double minLevel) const
496 : {
497 179846 : switch (border)
498 : {
499 50008 : case LEFT_BORDER:
500 50008 : return Point(upperLeft.x,
501 50008 : interpolate_(level, lowerLeft.y, upperLeft.y,
502 50008 : lowerLeft.value, upperLeft.value,
503 100016 : !split, minLevel));
504 39980 : case LOWER_BORDER:
505 39980 : return Point(interpolate_(level, lowerLeft.x, lowerRight.x,
506 39980 : lowerLeft.value, lowerRight.value,
507 39980 : !split, minLevel),
508 79960 : lowerLeft.y);
509 49605 : case RIGHT_BORDER:
510 49605 : return Point(upperRight.x,
511 49605 : interpolate_(level, lowerRight.y, upperRight.y,
512 49605 : lowerRight.value, upperRight.value,
513 99210 : !split, minLevel));
514 40253 : case UPPER_BORDER:
515 40253 : return Point(interpolate_(level, upperLeft.x, upperRight.x,
516 40253 : upperLeft.value, upperRight.value,
517 40253 : !split, minLevel),
518 80506 : upperLeft.y);
519 : }
520 0 : assert(false);
521 : return Point();
522 : }
523 : };
524 : } // namespace marching_squares
525 : #endif
|
