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 351064 : Segments(const Segment &first) : sz_(1), segs_()
58 : {
59 87766 : segs_[0] = first;
60 87766 : }
61 :
62 11196 : Segments(const Segment &first, const Segment &second) : sz_(2), segs_()
63 : {
64 2799 : segs_[0] = first;
65 2799 : segs_[1] = second;
66 2799 : }
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 183929 : std::size_t size() const
78 : {
79 183929 : return sz_;
80 : }
81 :
82 93370 : const Segment &operator[](std::size_t idx) const
83 : {
84 93370 : assert(idx < sz_);
85 93370 : return segs_[idx];
86 : }
87 :
88 : private:
89 : const std::size_t sz_;
90 : /* const */ Segment segs_[3];
91 : };
92 :
93 825701 : Square(const ValuedPoint &upperLeft_, const ValuedPoint &upperRight_,
94 : const ValuedPoint &lowerLeft_, const ValuedPoint &lowerRight_,
95 : uint8_t borders_ = NO_BORDER, bool split_ = false)
96 825701 : : upperLeft(upperLeft_), lowerLeft(lowerLeft_), lowerRight(lowerRight_),
97 : upperRight(upperRight_),
98 825701 : nanCount((std::isnan(upperLeft.value) ? 1 : 0) +
99 825701 : (std::isnan(upperRight.value) ? 1 : 0) +
100 825701 : (std::isnan(lowerLeft.value) ? 1 : 0) +
101 825701 : (std::isnan(lowerRight.value) ? 1 : 0)),
102 825701 : borders(borders_), split(split_)
103 : {
104 825701 : assert(upperLeft.y == upperRight.y);
105 825701 : assert(lowerLeft.y == lowerRight.y);
106 825701 : assert(lowerLeft.x == upperLeft.x);
107 825701 : assert(lowerRight.x == upperRight.x);
108 825701 : assert(!split || nanCount == 0);
109 825701 : }
110 :
111 10425 : Square upperLeftSquare() const
112 : {
113 10425 : assert(!std::isnan(upperLeft.value));
114 : return Square(
115 20850 : upperLeft, upperCenter(), leftCenter(), center(),
116 20850 : (std::isnan(upperRight.value) ? RIGHT_BORDER : NO_BORDER) |
117 10425 : (std::isnan(lowerLeft.value) ? LOWER_BORDER : NO_BORDER),
118 20850 : true);
119 : }
120 :
121 10424 : Square lowerLeftSquare() const
122 : {
123 10424 : assert(!std::isnan(lowerLeft.value));
124 : return Square(
125 20848 : leftCenter(), center(), lowerLeft, lowerCenter(),
126 20848 : (std::isnan(lowerRight.value) ? RIGHT_BORDER : NO_BORDER) |
127 10424 : (std::isnan(upperLeft.value) ? UPPER_BORDER : NO_BORDER),
128 20848 : true);
129 : }
130 :
131 10426 : Square lowerRightSquare() const
132 : {
133 10426 : assert(!std::isnan(lowerRight.value));
134 : return Square(
135 20852 : center(), rightCenter(), lowerCenter(), lowerRight,
136 20852 : (std::isnan(lowerLeft.value) ? LEFT_BORDER : NO_BORDER) |
137 10426 : (std::isnan(upperRight.value) ? UPPER_BORDER : NO_BORDER),
138 20852 : true);
139 : }
140 :
141 10422 : Square upperRightSquare() const
142 : {
143 10422 : assert(!std::isnan(upperRight.value));
144 : return Square(
145 20844 : upperCenter(), upperRight, center(), rightCenter(),
146 20844 : (std::isnan(lowerRight.value) ? LOWER_BORDER : NO_BORDER) |
147 10422 : (std::isnan(upperLeft.value) ? LEFT_BORDER : NO_BORDER),
148 20844 : true);
149 : }
150 :
151 804677 : double maxValue() const
152 : {
153 804677 : assert(nanCount == 0);
154 804677 : return std::max(std::max(upperLeft.value, upperRight.value),
155 804677 : std::max(lowerLeft.value, lowerRight.value));
156 : }
157 :
158 804677 : double minValue() const
159 : {
160 804677 : assert(nanCount == 0);
161 804677 : return std::min(std::min(upperLeft.value, upperRight.value),
162 804677 : std::min(lowerLeft.value, lowerRight.value));
163 : }
164 :
165 26009 : ValuedSegment segment(uint8_t border) const
166 : {
167 26009 : switch (border)
168 : {
169 6501 : case LEFT_BORDER:
170 6501 : return ValuedSegment(upperLeft, lowerLeft);
171 6504 : case LOWER_BORDER:
172 6504 : return ValuedSegment(lowerLeft, lowerRight);
173 6499 : case RIGHT_BORDER:
174 6499 : return ValuedSegment(lowerRight, upperRight);
175 6505 : case UPPER_BORDER:
176 6505 : 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 90571 : Segments segments(double level, double minLevel) const
192 : {
193 90571 : 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 9268 : case (UPPER_LEFT):
202 : // debug("UPPER_LEFT");
203 : return Segments(
204 9268 : Segment(interpolate(UPPER_BORDER, level, minLevel),
205 27804 : interpolate(LEFT_BORDER, level, minLevel)));
206 5811 : case (LOWER_LEFT):
207 : // debug("LOWER_LEFT");
208 : return Segments(
209 5811 : Segment(interpolate(LEFT_BORDER, level, minLevel),
210 17433 : interpolate(LOWER_BORDER, level, minLevel)));
211 3758 : case (LOWER_RIGHT):
212 : // debug("LOWER_RIGHT");
213 : return Segments(
214 3758 : Segment(interpolate(LOWER_BORDER, level, minLevel),
215 11274 : interpolate(RIGHT_BORDER, level, minLevel)));
216 4408 : case (UPPER_RIGHT):
217 : // debug("UPPER_RIGHT");
218 : return Segments(
219 4408 : Segment(interpolate(RIGHT_BORDER, level, minLevel),
220 13224 : interpolate(UPPER_BORDER, level, minLevel)));
221 9212 : case (UPPER_LEFT | LOWER_LEFT):
222 : // debug("UPPER_LEFT | LOWER_LEFT");
223 : return Segments(
224 9212 : Segment(interpolate(UPPER_BORDER, level, minLevel),
225 27636 : interpolate(LOWER_BORDER, level, minLevel)));
226 12476 : case (LOWER_LEFT | LOWER_RIGHT):
227 : // debug("LOWER_LEFT | LOWER_RIGHT");
228 : return Segments(
229 12476 : Segment(interpolate(LEFT_BORDER, level, minLevel),
230 37428 : interpolate(RIGHT_BORDER, level, minLevel)));
231 6076 : case (LOWER_RIGHT | UPPER_RIGHT):
232 : // debug("LOWER_RIGHT | UPPER_RIGHT");
233 : return Segments(
234 6076 : Segment(interpolate(LOWER_BORDER, level, minLevel),
235 18228 : interpolate(UPPER_BORDER, level, minLevel)));
236 12559 : case (UPPER_RIGHT | UPPER_LEFT):
237 : // debug("UPPER_RIGHT | UPPER_LEFT");
238 : return Segments(
239 12559 : Segment(interpolate(RIGHT_BORDER, level, minLevel),
240 37677 : interpolate(LEFT_BORDER, level, minLevel)));
241 4124 : case (ALL_HIGH & ~UPPER_LEFT):
242 : // debug("ALL_HIGH & ~UPPER_LEFT");
243 : return Segments(
244 4124 : Segment(interpolate(LEFT_BORDER, level, minLevel),
245 12372 : interpolate(UPPER_BORDER, level, minLevel)));
246 4722 : case (ALL_HIGH & ~LOWER_LEFT):
247 : // debug("ALL_HIGH & ~LOWER_LEFT");
248 : return Segments(
249 4722 : Segment(interpolate(LOWER_BORDER, level, minLevel),
250 14166 : interpolate(LEFT_BORDER, level, minLevel)));
251 9288 : case (ALL_HIGH & ~LOWER_RIGHT):
252 : // debug("ALL_HIGH & ~LOWER_RIGHT");
253 : return Segments(
254 9288 : Segment(interpolate(RIGHT_BORDER, level, minLevel),
255 27864 : interpolate(LOWER_BORDER, level, minLevel)));
256 6064 : case (ALL_HIGH & ~UPPER_RIGHT):
257 : // debug("ALL_HIGH & ~UPPER_RIGHT");
258 : return Segments(
259 6064 : Segment(interpolate(UPPER_BORDER, level, minLevel),
260 18192 : interpolate(RIGHT_BORDER, level, minLevel)));
261 2799 : 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 2799 : Segment(interpolate(LEFT_BORDER, level, minLevel),
277 2799 : interpolate(LOWER_BORDER, level, minLevel)),
278 2799 : Segment(interpolate(RIGHT_BORDER, level, minLevel),
279 8397 : interpolate(UPPER_BORDER, level, minLevel)));
280 : }
281 0 : assert(false);
282 : return Segments();
283 : }
284 :
285 : template <typename Writer, typename LevelGenerator>
286 825689 : void process(const LevelGenerator &levelGenerator, Writer &writer) const
287 : {
288 825689 : if (nanCount == 4) // nothing to do
289 21012 : return;
290 :
291 825688 : if (nanCount) // split in 4
292 : {
293 21011 : if (!std::isnan(upperLeft.value))
294 10422 : upperLeftSquare().process(levelGenerator, writer);
295 21011 : if (!std::isnan(upperRight.value))
296 10422 : upperRightSquare().process(levelGenerator, writer);
297 21011 : if (!std::isnan(lowerLeft.value))
298 10422 : lowerLeftSquare().process(levelGenerator, writer);
299 21011 : if (!std::isnan(lowerRight.value))
300 10425 : lowerRightSquare().process(levelGenerator, writer);
301 21011 : return;
302 : }
303 :
304 804677 : if (writer.polygonize && borders)
305 : {
306 129130 : 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 103304 : if ((border & borders) == 0)
311 77295 : 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 26009 : const ValuedSegment s(segment(border));
318 :
319 26009 : Point lastPoint(s.first.x, s.first.y);
320 26009 : Point endPoint(s.second.x, s.second.y);
321 26009 : if (s.first.value > s.second.value)
322 1622 : std::swap(lastPoint, endPoint);
323 26009 : bool reverse =
324 27171 : (s.first.value > s.second.value) &&
325 1162 : ((border == UPPER_BORDER) || (border == LEFT_BORDER));
326 :
327 26009 : auto levelIt =
328 26009 : levelGenerator.range(s.first.value, s.second.value);
329 :
330 26009 : auto it = levelIt.begin(); // reused after the for
331 26173 : for (; it != levelIt.end(); ++it)
332 : {
333 164 : const int levelIdx = (*it).first;
334 164 : const double level = (*it).second;
335 :
336 164 : const Point nextPoint(
337 : interpolate(border, level, levelGenerator.minLevel()));
338 164 : if (reverse)
339 42 : writer.addBorderSegment(levelIdx, nextPoint, lastPoint);
340 : else
341 122 : writer.addBorderSegment(levelIdx, lastPoint, nextPoint);
342 164 : lastPoint = nextPoint;
343 : }
344 : // last level (past the end)
345 26009 : if (reverse)
346 1016 : writer.addBorderSegment((*it).first, endPoint, lastPoint);
347 : else
348 24993 : writer.addBorderSegment((*it).first, lastPoint, endPoint);
349 : }
350 : }
351 :
352 804677 : auto range = levelGenerator.range(minValue(), maxValue());
353 804677 : auto it = range.begin();
354 804677 : auto itEnd = range.end();
355 804677 : auto next = range.begin();
356 804677 : ++next;
357 :
358 895236 : for (; it != itEnd; ++it, ++next)
359 : {
360 90559 : const int levelIdx = (*it).first;
361 90559 : const double level = (*it).second;
362 :
363 90559 : const Segments segments_ =
364 : segments(level, levelGenerator.minLevel());
365 :
366 183917 : for (std::size_t i = 0; i < segments_.size(); i++)
367 : {
368 93358 : const Segment &s = segments_[i];
369 93358 : writer.addSegment(levelIdx, s.first, s.second);
370 :
371 93358 : 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 16124 : 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 41697 : ValuedPoint center() const
395 : {
396 : return ValuedPoint(
397 41697 : .5 * (upperLeft.x + lowerRight.x),
398 41697 : .5 * (upperLeft.y + lowerRight.y),
399 41697 : ((std::isnan(lowerLeft.value) ? 0 : lowerLeft.value) +
400 41697 : (std::isnan(upperLeft.value) ? 0 : upperLeft.value) +
401 41697 : (std::isnan(lowerRight.value) ? 0 : lowerRight.value) +
402 41697 : (std::isnan(upperRight.value) ? 0 : upperRight.value)) /
403 41697 : (4 - nanCount));
404 : }
405 :
406 20849 : ValuedPoint leftCenter() const
407 : {
408 : return ValuedPoint(
409 20849 : upperLeft.x, .5 * (upperLeft.y + lowerLeft.y),
410 20849 : std::isnan(upperLeft.value)
411 : ? lowerLeft.value
412 15596 : : (std::isnan(lowerLeft.value)
413 15596 : ? upperLeft.value
414 36445 : : .5 * (upperLeft.value + lowerLeft.value)));
415 : }
416 :
417 20850 : ValuedPoint lowerCenter() const
418 : {
419 : return ValuedPoint(
420 20850 : .5 * (lowerLeft.x + lowerRight.x), lowerLeft.y,
421 20850 : std::isnan(lowerRight.value)
422 : ? lowerLeft.value
423 15598 : : (std::isnan(lowerLeft.value)
424 15598 : ? lowerRight.value
425 36448 : : .5 * (lowerRight.value + lowerLeft.value)));
426 : }
427 :
428 20848 : ValuedPoint rightCenter() const
429 : {
430 : return ValuedPoint(
431 20848 : upperRight.x, .5 * (upperRight.y + lowerRight.y),
432 20848 : std::isnan(lowerRight.value)
433 : ? upperRight.value
434 15596 : : (std::isnan(upperRight.value)
435 15596 : ? lowerRight.value
436 36444 : : .5 * (lowerRight.value + upperRight.value)));
437 : }
438 :
439 20847 : ValuedPoint upperCenter() const
440 : {
441 : return ValuedPoint(
442 20847 : .5 * (upperLeft.x + upperRight.x), upperLeft.y,
443 20847 : std::isnan(upperLeft.value)
444 : ? upperRight.value
445 15595 : : (std::isnan(upperRight.value)
446 15595 : ? upperLeft.value
447 36442 : : .5 * (upperLeft.value + upperRight.value)));
448 : }
449 :
450 90571 : uint8_t marchingCase(double level, double minLevel) const
451 : {
452 90571 : return (level < fudge(upperLeft.value, minLevel, level) ? UPPER_LEFT
453 181142 : : ALL_LOW) |
454 90571 : (level < fudge(lowerLeft.value, minLevel, level) ? LOWER_LEFT
455 90571 : : ALL_LOW) |
456 90571 : (level < fudge(lowerRight.value, minLevel, level) ? LOWER_RIGHT
457 90571 : : ALL_LOW) |
458 90571 : (level < fudge(upperRight.value, minLevel, level) ? UPPER_RIGHT
459 90571 : : ALL_LOW);
460 : }
461 :
462 186892 : static double interpolate_(double level, double x1, double x2, double y1,
463 : double y2, bool need_split, double minLevel)
464 : {
465 186892 : 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 180908 : const double xm = .5 * (x1 + x2);
476 180908 : const double ym = .5 * (y1 + y2);
477 180908 : const double fy1 = fudge(y1, minLevel, level);
478 180908 : const double fym = fudge(ym, minLevel, level);
479 180908 : if ((fy1 < level && level < fym) || (fy1 > level && level > fym))
480 : {
481 88929 : x2 = xm;
482 88929 : y2 = ym;
483 : }
484 : else
485 : {
486 91979 : x1 = xm;
487 91979 : y1 = ym;
488 : }
489 : }
490 186892 : const double fy1 = fudge(y1, minLevel, level);
491 186892 : const double ratio = (level - fy1) / (fudge(y2, minLevel, level) - fy1);
492 186892 : return x1 * (1. - ratio) + x2 * ratio;
493 : }
494 :
495 186892 : Point interpolate(uint8_t border, double level, double minLevel) const
496 : {
497 186892 : switch (border)
498 : {
499 51818 : case LEFT_BORDER:
500 51818 : return Point(upperLeft.x,
501 51818 : interpolate_(level, lowerLeft.y, upperLeft.y,
502 51818 : lowerLeft.value, upperLeft.value,
503 103636 : !split, minLevel));
504 41689 : case LOWER_BORDER:
505 41689 : return Point(interpolate_(level, lowerLeft.x, lowerRight.x,
506 41689 : lowerLeft.value, lowerRight.value,
507 41689 : !split, minLevel),
508 83378 : lowerLeft.y);
509 51409 : case RIGHT_BORDER:
510 51409 : return Point(upperRight.x,
511 51409 : interpolate_(level, lowerRight.y, upperRight.y,
512 51409 : lowerRight.value, upperRight.value,
513 102818 : !split, minLevel));
514 41976 : case UPPER_BORDER:
515 41976 : return Point(interpolate_(level, upperLeft.x, upperRight.x,
516 41976 : upperLeft.value, upperRight.value,
517 41976 : !split, minLevel),
518 83952 : upperLeft.y);
519 : }
520 0 : assert(false);
521 : return Point();
522 : }
523 : };
524 : } // namespace marching_squares
525 : #endif
|
