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 280672 : Segments(const Segment &first) : sz_(1), segs_()
58 : {
59 70168 : segs_[0] = first;
60 70168 : }
61 :
62 5548 : Segments(const Segment &first, const Segment &second) : sz_(2), segs_()
63 : {
64 1387 : segs_[0] = first;
65 1387 : segs_[1] = second;
66 1387 : }
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 144497 : std::size_t size() const
78 : {
79 144497 : return sz_;
80 : }
81 :
82 72948 : const Segment &operator[](std::size_t idx) const
83 : {
84 72948 : assert(idx < sz_);
85 72948 : return segs_[idx];
86 : }
87 :
88 : private:
89 : const std::size_t sz_;
90 : /* const */ Segment segs_[3];
91 : };
92 :
93 632985 : Square(const ValuedPoint &upperLeft_, const ValuedPoint &upperRight_,
94 : const ValuedPoint &lowerLeft_, const ValuedPoint &lowerRight_,
95 : uint8_t borders_ = NO_BORDER, bool split_ = false)
96 632985 : : upperLeft(upperLeft_), lowerLeft(lowerLeft_), lowerRight(lowerRight_),
97 : upperRight(upperRight_),
98 632985 : nanCount((std::isnan(upperLeft.value) ? 1 : 0) +
99 632985 : (std::isnan(upperRight.value) ? 1 : 0) +
100 632985 : (std::isnan(lowerLeft.value) ? 1 : 0) +
101 632985 : (std::isnan(lowerRight.value) ? 1 : 0)),
102 632985 : borders(borders_), split(split_)
103 : {
104 632985 : assert(upperLeft.y == upperRight.y);
105 632985 : assert(lowerLeft.y == lowerRight.y);
106 632985 : assert(lowerLeft.x == upperLeft.x);
107 632985 : assert(lowerRight.x == upperRight.x);
108 632985 : assert(!split || nanCount == 0);
109 632985 : }
110 :
111 7997 : Square upperLeftSquare() const
112 : {
113 7997 : assert(!std::isnan(upperLeft.value));
114 : return Square(
115 15994 : upperLeft, upperCenter(), leftCenter(), center(),
116 15994 : (std::isnan(upperRight.value) ? RIGHT_BORDER : NO_BORDER) |
117 7997 : (std::isnan(lowerLeft.value) ? LOWER_BORDER : NO_BORDER),
118 15994 : true);
119 : }
120 :
121 7996 : Square lowerLeftSquare() const
122 : {
123 7996 : assert(!std::isnan(lowerLeft.value));
124 : return Square(
125 15992 : leftCenter(), center(), lowerLeft, lowerCenter(),
126 15992 : (std::isnan(lowerRight.value) ? RIGHT_BORDER : NO_BORDER) |
127 7996 : (std::isnan(upperLeft.value) ? UPPER_BORDER : NO_BORDER),
128 15992 : true);
129 : }
130 :
131 7998 : Square lowerRightSquare() const
132 : {
133 7998 : assert(!std::isnan(lowerRight.value));
134 : return Square(
135 15996 : center(), rightCenter(), lowerCenter(), lowerRight,
136 15996 : (std::isnan(lowerLeft.value) ? LEFT_BORDER : NO_BORDER) |
137 7998 : (std::isnan(upperRight.value) ? UPPER_BORDER : NO_BORDER),
138 15996 : true);
139 : }
140 :
141 7994 : Square upperRightSquare() const
142 : {
143 7994 : assert(!std::isnan(upperRight.value));
144 : return Square(
145 15988 : upperCenter(), upperRight, center(), rightCenter(),
146 15988 : (std::isnan(lowerRight.value) ? LOWER_BORDER : NO_BORDER) |
147 7994 : (std::isnan(upperLeft.value) ? LEFT_BORDER : NO_BORDER),
148 15988 : true);
149 : }
150 :
151 616890 : double maxValue() const
152 : {
153 616890 : assert(nanCount == 0);
154 616890 : return std::max(std::max(upperLeft.value, upperRight.value),
155 616890 : std::max(lowerLeft.value, lowerRight.value));
156 : }
157 :
158 616890 : double minValue() const
159 : {
160 616890 : assert(nanCount == 0);
161 616890 : return std::min(std::min(upperLeft.value, upperRight.value),
162 616890 : std::min(lowerLeft.value, lowerRight.value));
163 : }
164 :
165 16905 : ValuedSegment segment(uint8_t border) const
166 : {
167 16905 : switch (border)
168 : {
169 4225 : case LEFT_BORDER:
170 4225 : return ValuedSegment(upperLeft, lowerLeft);
171 4228 : case LOWER_BORDER:
172 4228 : return ValuedSegment(lowerLeft, lowerRight);
173 4223 : case RIGHT_BORDER:
174 4223 : return ValuedSegment(lowerRight, upperRight);
175 4229 : case UPPER_BORDER:
176 4229 : 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 71561 : Segments segments(double level, double minLevel) const
192 : {
193 71561 : 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 7447 : case (UPPER_LEFT):
202 : // debug("UPPER_LEFT");
203 : return Segments(
204 7447 : Segment(interpolate(UPPER_BORDER, level, minLevel),
205 22341 : interpolate(LEFT_BORDER, level, minLevel)));
206 4249 : case (LOWER_LEFT):
207 : // debug("LOWER_LEFT");
208 : return Segments(
209 4249 : Segment(interpolate(LEFT_BORDER, level, minLevel),
210 12747 : interpolate(LOWER_BORDER, level, minLevel)));
211 2234 : case (LOWER_RIGHT):
212 : // debug("LOWER_RIGHT");
213 : return Segments(
214 2234 : Segment(interpolate(LOWER_BORDER, level, minLevel),
215 6702 : interpolate(RIGHT_BORDER, level, minLevel)));
216 3099 : case (UPPER_RIGHT):
217 : // debug("UPPER_RIGHT");
218 : return Segments(
219 3099 : Segment(interpolate(RIGHT_BORDER, level, minLevel),
220 9297 : interpolate(UPPER_BORDER, level, minLevel)));
221 7624 : case (UPPER_LEFT | LOWER_LEFT):
222 : // debug("UPPER_LEFT | LOWER_LEFT");
223 : return Segments(
224 7624 : Segment(interpolate(UPPER_BORDER, level, minLevel),
225 22872 : interpolate(LOWER_BORDER, level, minLevel)));
226 10575 : case (LOWER_LEFT | LOWER_RIGHT):
227 : // debug("LOWER_LEFT | LOWER_RIGHT");
228 : return Segments(
229 10575 : Segment(interpolate(LEFT_BORDER, level, minLevel),
230 31725 : interpolate(RIGHT_BORDER, level, minLevel)));
231 4599 : case (LOWER_RIGHT | UPPER_RIGHT):
232 : // debug("LOWER_RIGHT | UPPER_RIGHT");
233 : return Segments(
234 4599 : Segment(interpolate(LOWER_BORDER, level, minLevel),
235 13797 : interpolate(UPPER_BORDER, level, minLevel)));
236 11030 : case (UPPER_RIGHT | UPPER_LEFT):
237 : // debug("UPPER_RIGHT | UPPER_LEFT");
238 : return Segments(
239 11030 : Segment(interpolate(RIGHT_BORDER, level, minLevel),
240 33090 : interpolate(LEFT_BORDER, level, minLevel)));
241 2892 : case (ALL_HIGH & ~UPPER_LEFT):
242 : // debug("ALL_HIGH & ~UPPER_LEFT");
243 : return Segments(
244 2892 : Segment(interpolate(LEFT_BORDER, level, minLevel),
245 8676 : interpolate(UPPER_BORDER, level, minLevel)));
246 3728 : case (ALL_HIGH & ~LOWER_LEFT):
247 : // debug("ALL_HIGH & ~LOWER_LEFT");
248 : return Segments(
249 3728 : Segment(interpolate(LOWER_BORDER, level, minLevel),
250 11184 : interpolate(LEFT_BORDER, level, minLevel)));
251 7841 : case (ALL_HIGH & ~LOWER_RIGHT):
252 : // debug("ALL_HIGH & ~LOWER_RIGHT");
253 : return Segments(
254 7841 : Segment(interpolate(RIGHT_BORDER, level, minLevel),
255 23523 : interpolate(LOWER_BORDER, level, minLevel)));
256 4850 : case (ALL_HIGH & ~UPPER_RIGHT):
257 : // debug("ALL_HIGH & ~UPPER_RIGHT");
258 : return Segments(
259 4850 : Segment(interpolate(UPPER_BORDER, level, minLevel),
260 14550 : interpolate(RIGHT_BORDER, level, minLevel)));
261 1387 : 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 1387 : Segment(interpolate(LEFT_BORDER, level, minLevel),
277 1387 : interpolate(LOWER_BORDER, level, minLevel)),
278 1387 : Segment(interpolate(RIGHT_BORDER, level, minLevel),
279 4161 : interpolate(UPPER_BORDER, level, minLevel)));
280 : }
281 0 : assert(false);
282 : return Segments();
283 : }
284 :
285 : template <typename Writer, typename LevelGenerator>
286 632973 : void process(const LevelGenerator &levelGenerator, Writer &writer) const
287 : {
288 632973 : if (nanCount == 4) // nothing to do
289 16083 : return;
290 :
291 632973 : if (nanCount) // split in 4
292 : {
293 16083 : if (!std::isnan(upperLeft.value))
294 7994 : upperLeftSquare().process(levelGenerator, writer);
295 16083 : if (!std::isnan(upperRight.value))
296 7994 : upperRightSquare().process(levelGenerator, writer);
297 16083 : if (!std::isnan(lowerLeft.value))
298 7994 : lowerLeftSquare().process(levelGenerator, writer);
299 16083 : if (!std::isnan(lowerRight.value))
300 7997 : lowerRightSquare().process(levelGenerator, writer);
301 16083 : return;
302 : }
303 :
304 616890 : if (writer.polygonize && borders)
305 : {
306 83930 : 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 67144 : if ((border & borders) == 0)
311 50239 : 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 16905 : const ValuedSegment s(segment(border));
318 :
319 16905 : Point lastPoint(s.first.x, s.first.y);
320 16905 : Point endPoint(s.second.x, s.second.y);
321 16905 : if (s.first.value > s.second.value)
322 1586 : std::swap(lastPoint, endPoint);
323 16905 : bool reverse =
324 18049 : (s.first.value > s.second.value) &&
325 1144 : ((border == UPPER_BORDER) || (border == LEFT_BORDER));
326 :
327 16905 : auto levelIt =
328 16905 : levelGenerator.range(s.first.value, s.second.value);
329 :
330 16905 : auto it = levelIt.begin(); // reused after the for
331 17007 : for (; it != levelIt.end(); ++it)
332 : {
333 102 : const int levelIdx = (*it).first;
334 102 : const double level = (*it).second;
335 :
336 102 : const Point nextPoint(
337 : interpolate(border, level, levelGenerator.minLevel()));
338 102 : if (reverse)
339 28 : writer.addBorderSegment(levelIdx, nextPoint, lastPoint);
340 : else
341 74 : writer.addBorderSegment(levelIdx, lastPoint, nextPoint);
342 102 : lastPoint = nextPoint;
343 : }
344 : // last level (past the end)
345 16905 : if (reverse)
346 996 : writer.addBorderSegment((*it).first, endPoint, lastPoint);
347 : else
348 15909 : writer.addBorderSegment((*it).first, lastPoint, endPoint);
349 : }
350 : }
351 :
352 616890 : auto range = levelGenerator.range(minValue(), maxValue());
353 616890 : auto it = range.begin();
354 616890 : auto itEnd = range.end();
355 616890 : auto next = range.begin();
356 616890 : ++next;
357 :
358 688439 : for (; it != itEnd; ++it, ++next)
359 : {
360 71549 : const int levelIdx = (*it).first;
361 71549 : const double level = (*it).second;
362 :
363 71549 : const Segments segments_ =
364 : segments(level, levelGenerator.minLevel());
365 :
366 144485 : for (std::size_t i = 0; i < segments_.size(); i++)
367 : {
368 72936 : const Segment &s = segments_[i];
369 72936 : writer.addSegment(levelIdx, s.first, s.second);
370 :
371 72936 : 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 14587 : 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 31985 : ValuedPoint center() const
395 : {
396 : return ValuedPoint(
397 31985 : .5 * (upperLeft.x + lowerRight.x),
398 31985 : .5 * (upperLeft.y + lowerRight.y),
399 31985 : ((std::isnan(lowerLeft.value) ? 0 : lowerLeft.value) +
400 31985 : (std::isnan(upperLeft.value) ? 0 : upperLeft.value) +
401 31985 : (std::isnan(lowerRight.value) ? 0 : lowerRight.value) +
402 31985 : (std::isnan(upperRight.value) ? 0 : upperRight.value)) /
403 31985 : (4 - nanCount));
404 : }
405 :
406 15993 : ValuedPoint leftCenter() const
407 : {
408 : return ValuedPoint(
409 15993 : upperLeft.x, .5 * (upperLeft.y + lowerLeft.y),
410 15993 : std::isnan(upperLeft.value)
411 : ? lowerLeft.value
412 11972 : : (std::isnan(lowerLeft.value)
413 11972 : ? upperLeft.value
414 27965 : : .5 * (upperLeft.value + lowerLeft.value)));
415 : }
416 :
417 15994 : ValuedPoint lowerCenter() const
418 : {
419 : return ValuedPoint(
420 15994 : .5 * (lowerLeft.x + lowerRight.x), lowerLeft.y,
421 15994 : std::isnan(lowerRight.value)
422 : ? lowerLeft.value
423 11974 : : (std::isnan(lowerLeft.value)
424 11974 : ? lowerRight.value
425 27968 : : .5 * (lowerRight.value + lowerLeft.value)));
426 : }
427 :
428 15992 : ValuedPoint rightCenter() const
429 : {
430 : return ValuedPoint(
431 15992 : upperRight.x, .5 * (upperRight.y + lowerRight.y),
432 15992 : std::isnan(lowerRight.value)
433 : ? upperRight.value
434 11972 : : (std::isnan(upperRight.value)
435 11972 : ? lowerRight.value
436 27964 : : .5 * (lowerRight.value + upperRight.value)));
437 : }
438 :
439 15991 : ValuedPoint upperCenter() const
440 : {
441 : return ValuedPoint(
442 15991 : .5 * (upperLeft.x + upperRight.x), upperLeft.y,
443 15991 : std::isnan(upperLeft.value)
444 : ? upperRight.value
445 11971 : : (std::isnan(upperRight.value)
446 11971 : ? upperLeft.value
447 27962 : : .5 * (upperLeft.value + upperRight.value)));
448 : }
449 :
450 71561 : uint8_t marchingCase(double level, double minLevel) const
451 : {
452 71561 : return (level < fudge(upperLeft.value, minLevel, level) ? UPPER_LEFT
453 143122 : : ALL_LOW) |
454 71561 : (level < fudge(lowerLeft.value, minLevel, level) ? LOWER_LEFT
455 71561 : : ALL_LOW) |
456 71561 : (level < fudge(lowerRight.value, minLevel, level) ? LOWER_RIGHT
457 71561 : : ALL_LOW) |
458 71561 : (level < fudge(upperRight.value, minLevel, level) ? UPPER_RIGHT
459 71561 : : ALL_LOW);
460 : }
461 :
462 145986 : static double interpolate_(double level, double x1, double x2, double y1,
463 : double y2, bool need_split, double minLevel)
464 : {
465 145986 : 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 143436 : const double xm = .5 * (x1 + x2);
476 143436 : const double ym = .5 * (y1 + y2);
477 143436 : const double fy1 = fudge(y1, minLevel, level);
478 143436 : const double fym = fudge(ym, minLevel, level);
479 143436 : if ((fy1 < level && level < fym) || (fy1 > level && level > fym))
480 : {
481 70150 : x2 = xm;
482 70150 : y2 = ym;
483 : }
484 : else
485 : {
486 73286 : x1 = xm;
487 73286 : y1 = ym;
488 : }
489 : }
490 145986 : const double fy1 = fudge(y1, minLevel, level);
491 145986 : const double ratio = (level - fy1) / (fudge(y2, minLevel, level) - fy1);
492 145986 : return x1 * (1. - ratio) + x2 * ratio;
493 : }
494 :
495 145986 : Point interpolate(uint8_t border, double level, double minLevel) const
496 : {
497 145986 : switch (border)
498 : {
499 41344 : case LEFT_BORDER:
500 41344 : return Point(upperLeft.x,
501 41344 : interpolate_(level, lowerLeft.y, upperLeft.y,
502 41344 : lowerLeft.value, upperLeft.value,
503 82688 : !split, minLevel));
504 31672 : case LOWER_BORDER:
505 31672 : return Point(interpolate_(level, lowerLeft.x, lowerRight.x,
506 31672 : lowerLeft.value, lowerRight.value,
507 31672 : !split, minLevel),
508 63344 : lowerLeft.y);
509 41056 : case RIGHT_BORDER:
510 41056 : return Point(upperRight.x,
511 41056 : interpolate_(level, lowerRight.y, upperRight.y,
512 41056 : lowerRight.value, upperRight.value,
513 82112 : !split, minLevel));
514 31914 : case UPPER_BORDER:
515 31914 : return Point(interpolate_(level, upperLeft.x, upperRight.x,
516 31914 : upperLeft.value, upperRight.value,
517 31914 : !split, minLevel),
518 63828 : upperLeft.y);
519 : }
520 0 : assert(false);
521 : return Point();
522 : }
523 : };
524 : } // namespace marching_squares
525 : #endif
|
