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 352820 : Segments(const Segment &first) : sz_(1), segs_()
58 : {
59 88205 : segs_[0] = first;
60 88205 : }
61 :
62 11716 : Segments(const Segment &first, const Segment &second) : sz_(2), segs_()
63 : {
64 2929 : segs_[0] = first;
65 2929 : segs_[1] = second;
66 2929 : }
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 185197 : std::size_t size() const
78 : {
79 185197 : return sz_;
80 : }
81 :
82 94069 : const Segment &operator[](std::size_t idx) const
83 : {
84 94069 : assert(idx < sz_);
85 94069 : return segs_[idx];
86 : }
87 :
88 : private:
89 : const std::size_t sz_;
90 : /* const */ Segment segs_[3];
91 : };
92 :
93 826457 : Square(const ValuedPoint &upperLeft_, const ValuedPoint &upperRight_,
94 : const ValuedPoint &lowerLeft_, const ValuedPoint &lowerRight_,
95 : uint8_t borders_ = NO_BORDER, bool split_ = false)
96 826457 : : upperLeft(upperLeft_), lowerLeft(lowerLeft_), lowerRight(lowerRight_),
97 : upperRight(upperRight_),
98 826457 : nanCount((std::isnan(upperLeft.value) ? 1 : 0) +
99 826457 : (std::isnan(upperRight.value) ? 1 : 0) +
100 826457 : (std::isnan(lowerLeft.value) ? 1 : 0) +
101 826457 : (std::isnan(lowerRight.value) ? 1 : 0)),
102 826457 : borders(borders_), split(split_)
103 : {
104 826457 : assert(upperLeft.y == upperRight.y);
105 826457 : assert(lowerLeft.y == lowerRight.y);
106 826457 : assert(lowerLeft.x == upperLeft.x);
107 826457 : assert(lowerRight.x == upperRight.x);
108 826457 : assert(!split || nanCount == 0);
109 826457 : }
110 :
111 10485 : Square upperLeftSquare() const
112 : {
113 10485 : assert(!std::isnan(upperLeft.value));
114 : return Square(
115 20970 : upperLeft, upperCenter(), leftCenter(), center(),
116 20970 : (std::isnan(upperRight.value) ? RIGHT_BORDER : NO_BORDER) |
117 10485 : (std::isnan(lowerLeft.value) ? LOWER_BORDER : NO_BORDER),
118 20970 : true);
119 : }
120 :
121 10484 : Square lowerLeftSquare() const
122 : {
123 10484 : assert(!std::isnan(lowerLeft.value));
124 : return Square(
125 20968 : leftCenter(), center(), lowerLeft, lowerCenter(),
126 20968 : (std::isnan(lowerRight.value) ? RIGHT_BORDER : NO_BORDER) |
127 10484 : (std::isnan(upperLeft.value) ? UPPER_BORDER : NO_BORDER),
128 20968 : true);
129 : }
130 :
131 10486 : Square lowerRightSquare() const
132 : {
133 10486 : assert(!std::isnan(lowerRight.value));
134 : return Square(
135 20972 : center(), rightCenter(), lowerCenter(), lowerRight,
136 20972 : (std::isnan(lowerLeft.value) ? LEFT_BORDER : NO_BORDER) |
137 10486 : (std::isnan(upperRight.value) ? UPPER_BORDER : NO_BORDER),
138 20972 : true);
139 : }
140 :
141 10482 : Square upperRightSquare() const
142 : {
143 10482 : assert(!std::isnan(upperRight.value));
144 : return Square(
145 20964 : upperCenter(), upperRight, center(), rightCenter(),
146 20964 : (std::isnan(lowerRight.value) ? LOWER_BORDER : NO_BORDER) |
147 10482 : (std::isnan(upperLeft.value) ? LEFT_BORDER : NO_BORDER),
148 20964 : true);
149 : }
150 :
151 805305 : double maxValue() const
152 : {
153 805305 : assert(nanCount == 0);
154 805305 : return std::max(std::max(upperLeft.value, upperRight.value),
155 805305 : std::max(lowerLeft.value, lowerRight.value));
156 : }
157 :
158 805305 : double minValue() const
159 : {
160 805305 : assert(nanCount == 0);
161 805305 : return std::min(std::min(upperLeft.value, upperRight.value),
162 805305 : 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 91140 : Segments segments(double level) const
192 : {
193 91140 : switch (marchingCase(level))
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 9414 : case (UPPER_LEFT):
202 : // debug("UPPER_LEFT");
203 9414 : return Segments(Segment(interpolate(UPPER_BORDER, level),
204 28242 : interpolate(LEFT_BORDER, level)));
205 5935 : case (LOWER_LEFT):
206 : // debug("LOWER_LEFT");
207 5935 : return Segments(Segment(interpolate(LEFT_BORDER, level),
208 17805 : interpolate(LOWER_BORDER, level)));
209 3881 : case (LOWER_RIGHT):
210 : // debug("LOWER_RIGHT");
211 3881 : return Segments(Segment(interpolate(LOWER_BORDER, level),
212 11643 : interpolate(RIGHT_BORDER, level)));
213 4513 : case (UPPER_RIGHT):
214 : // debug("UPPER_RIGHT");
215 4513 : return Segments(Segment(interpolate(RIGHT_BORDER, level),
216 13539 : interpolate(UPPER_BORDER, level)));
217 9088 : case (UPPER_LEFT | LOWER_LEFT):
218 : // debug("UPPER_LEFT | LOWER_LEFT");
219 9088 : return Segments(Segment(interpolate(UPPER_BORDER, level),
220 27264 : interpolate(LOWER_BORDER, level)));
221 12366 : case (LOWER_LEFT | LOWER_RIGHT):
222 : // debug("LOWER_LEFT | LOWER_RIGHT");
223 12366 : return Segments(Segment(interpolate(LEFT_BORDER, level),
224 37098 : interpolate(RIGHT_BORDER, level)));
225 5940 : case (LOWER_RIGHT | UPPER_RIGHT):
226 : // debug("LOWER_RIGHT | UPPER_RIGHT");
227 5940 : return Segments(Segment(interpolate(LOWER_BORDER, level),
228 17820 : interpolate(UPPER_BORDER, level)));
229 12426 : case (UPPER_RIGHT | UPPER_LEFT):
230 : // debug("UPPER_RIGHT | UPPER_LEFT");
231 12426 : return Segments(Segment(interpolate(RIGHT_BORDER, level),
232 37278 : interpolate(LEFT_BORDER, level)));
233 4235 : case (ALL_HIGH & ~UPPER_LEFT):
234 : // debug("ALL_HIGH & ~UPPER_LEFT");
235 4235 : return Segments(Segment(interpolate(LEFT_BORDER, level),
236 12705 : interpolate(UPPER_BORDER, level)));
237 4812 : case (ALL_HIGH & ~LOWER_LEFT):
238 : // debug("ALL_HIGH & ~LOWER_LEFT");
239 4812 : return Segments(Segment(interpolate(LOWER_BORDER, level),
240 14436 : interpolate(LEFT_BORDER, level)));
241 9420 : case (ALL_HIGH & ~LOWER_RIGHT):
242 : // debug("ALL_HIGH & ~LOWER_RIGHT");
243 9420 : return Segments(Segment(interpolate(RIGHT_BORDER, level),
244 28260 : interpolate(LOWER_BORDER, level)));
245 6175 : case (ALL_HIGH & ~UPPER_RIGHT):
246 : // debug("ALL_HIGH & ~UPPER_RIGHT");
247 6175 : return Segments(Segment(interpolate(UPPER_BORDER, level),
248 18525 : interpolate(RIGHT_BORDER, level)));
249 2929 : case (SADDLE_NE):
250 : case (SADDLE_NW):
251 : // From the two possible saddle configurations, we always return
252 : // the same one.
253 : //
254 : // The classical marching square algorithm says the ambiguity
255 : // should be resolved between the two possible configurations by
256 : // looking at the value of the center point. But in certain
257 : // cases, this may lead to line contours from different levels
258 : // that cross each other and then gives invalid polygons.
259 : //
260 : // Arbitrarily choosing one of the two possible configurations
261 : // is not really that worse than deciding based on the center
262 : // point.
263 2929 : return Segments(Segment(interpolate(LEFT_BORDER, level),
264 2929 : interpolate(LOWER_BORDER, level)),
265 2929 : Segment(interpolate(RIGHT_BORDER, level),
266 8787 : interpolate(UPPER_BORDER, level)));
267 : }
268 0 : assert(false);
269 : return Segments();
270 : }
271 :
272 : template <typename Writer, typename LevelGenerator>
273 826445 : void process(const LevelGenerator &levelGenerator, Writer &writer) const
274 : {
275 826445 : if (nanCount == 4) // nothing to do
276 21140 : return;
277 :
278 826444 : if (nanCount) // split in 4
279 : {
280 21139 : if (!std::isnan(upperLeft.value))
281 10482 : upperLeftSquare().process(levelGenerator, writer);
282 21139 : if (!std::isnan(upperRight.value))
283 10482 : upperRightSquare().process(levelGenerator, writer);
284 21139 : if (!std::isnan(lowerLeft.value))
285 10482 : lowerLeftSquare().process(levelGenerator, writer);
286 21139 : if (!std::isnan(lowerRight.value))
287 10485 : lowerRightSquare().process(levelGenerator, writer);
288 21139 : return;
289 : }
290 :
291 805305 : if (writer.polygonize && borders)
292 : {
293 129130 : for (uint8_t border :
294 : {UPPER_BORDER, LEFT_BORDER, RIGHT_BORDER, LOWER_BORDER})
295 : {
296 : // bitwise AND to test which borders we have on the square
297 103304 : if ((border & borders) == 0)
298 77295 : continue;
299 :
300 : // convention: for a level = L, store borders for the previous
301 : // level up to (and including) L in the border of level "L". For
302 : // fixed sets of level, this means there is an "Inf" slot for
303 : // borders of the highest level
304 26009 : const ValuedSegment s(segment(border));
305 :
306 26009 : Point lastPoint(s.first.x, s.first.y);
307 26009 : Point endPoint(s.second.x, s.second.y);
308 26009 : if (s.first.value > s.second.value)
309 1622 : std::swap(lastPoint, endPoint);
310 26009 : bool reverse =
311 27171 : (s.first.value > s.second.value) &&
312 1162 : ((border == UPPER_BORDER) || (border == LEFT_BORDER));
313 :
314 26009 : auto levelIt =
315 26009 : levelGenerator.range(s.first.value, s.second.value);
316 :
317 26009 : auto it = levelIt.begin(); // reused after the for
318 26173 : for (; it != levelIt.end(); ++it)
319 : {
320 164 : const int levelIdx = (*it).first;
321 164 : const double level = (*it).second;
322 :
323 164 : const Point nextPoint(interpolate(border, level));
324 164 : if (reverse)
325 42 : writer.addBorderSegment(levelIdx, nextPoint, lastPoint);
326 : else
327 122 : writer.addBorderSegment(levelIdx, lastPoint, nextPoint);
328 164 : lastPoint = nextPoint;
329 : }
330 : // last level (past the end)
331 26009 : if (reverse)
332 1016 : writer.addBorderSegment((*it).first, endPoint, lastPoint);
333 : else
334 24993 : writer.addBorderSegment((*it).first, lastPoint, endPoint);
335 : }
336 : }
337 :
338 805305 : auto range = levelGenerator.range(minValue(), maxValue());
339 805305 : auto it = range.begin();
340 805305 : auto itEnd = range.end();
341 805305 : auto next = range.begin();
342 805305 : ++next;
343 :
344 896433 : for (; it != itEnd; ++it, ++next)
345 : {
346 91128 : const int levelIdx = (*it).first;
347 91128 : const double level = (*it).second;
348 :
349 91128 : const Segments segments_ = segments(level);
350 :
351 185185 : for (std::size_t i = 0; i < segments_.size(); i++)
352 : {
353 94057 : const Segment &s = segments_[i];
354 94057 : writer.addSegment(levelIdx, s.first, s.second);
355 :
356 94057 : if (writer.polygonize)
357 : {
358 : // the contour is used in the polygon of higher level as
359 : // well
360 : //
361 : // TODO: copying the segment to the higher level is easy,
362 : // but it involves too much memory. We should reuse segment
363 : // contours when constructing polygon rings.
364 16124 : writer.addSegment((*next).first, s.first, s.second);
365 : }
366 : }
367 : }
368 : }
369 :
370 : const ValuedPoint upperLeft;
371 : const ValuedPoint lowerLeft;
372 : const ValuedPoint lowerRight;
373 : const ValuedPoint upperRight;
374 : const int nanCount;
375 : const uint8_t borders;
376 : const bool split;
377 :
378 : private:
379 41937 : ValuedPoint center() const
380 : {
381 : return ValuedPoint(
382 41937 : .5 * (upperLeft.x + lowerRight.x),
383 41937 : .5 * (upperLeft.y + lowerRight.y),
384 41937 : ((std::isnan(lowerLeft.value) ? 0 : lowerLeft.value) +
385 41937 : (std::isnan(upperLeft.value) ? 0 : upperLeft.value) +
386 41937 : (std::isnan(lowerRight.value) ? 0 : lowerRight.value) +
387 41937 : (std::isnan(upperRight.value) ? 0 : upperRight.value)) /
388 41937 : (4 - nanCount));
389 : }
390 :
391 20969 : ValuedPoint leftCenter() const
392 : {
393 : return ValuedPoint(
394 20969 : upperLeft.x, .5 * (upperLeft.y + lowerLeft.y),
395 20969 : std::isnan(upperLeft.value)
396 : ? lowerLeft.value
397 15684 : : (std::isnan(lowerLeft.value)
398 15684 : ? upperLeft.value
399 36653 : : .5 * (upperLeft.value + lowerLeft.value)));
400 : }
401 :
402 20970 : ValuedPoint lowerCenter() const
403 : {
404 : return ValuedPoint(
405 20970 : .5 * (lowerLeft.x + lowerRight.x), lowerLeft.y,
406 20970 : std::isnan(lowerRight.value)
407 : ? lowerLeft.value
408 15686 : : (std::isnan(lowerLeft.value)
409 15686 : ? lowerRight.value
410 36656 : : .5 * (lowerRight.value + lowerLeft.value)));
411 : }
412 :
413 20968 : ValuedPoint rightCenter() const
414 : {
415 : return ValuedPoint(
416 20968 : upperRight.x, .5 * (upperRight.y + lowerRight.y),
417 20968 : std::isnan(lowerRight.value)
418 : ? upperRight.value
419 15684 : : (std::isnan(upperRight.value)
420 15684 : ? lowerRight.value
421 36652 : : .5 * (lowerRight.value + upperRight.value)));
422 : }
423 :
424 20967 : ValuedPoint upperCenter() const
425 : {
426 : return ValuedPoint(
427 20967 : .5 * (upperLeft.x + upperRight.x), upperLeft.y,
428 20967 : std::isnan(upperLeft.value)
429 : ? upperRight.value
430 15683 : : (std::isnan(upperRight.value)
431 15683 : ? upperLeft.value
432 36650 : : .5 * (upperLeft.value + upperRight.value)));
433 : }
434 :
435 91140 : uint8_t marchingCase(double level) const
436 : {
437 182280 : return (level < fudge(upperLeft.value, level) ? UPPER_LEFT : ALL_LOW) |
438 182280 : (level < fudge(lowerLeft.value, level) ? LOWER_LEFT : ALL_LOW) |
439 91140 : (level < fudge(lowerRight.value, level) ? LOWER_RIGHT
440 91140 : : ALL_LOW) |
441 182280 : (level < fudge(upperRight.value, level) ? UPPER_RIGHT : ALL_LOW);
442 : }
443 :
444 188290 : static double interpolate_(double level, double x1, double x2, double y1,
445 : double y2, bool need_split)
446 : {
447 188290 : if (need_split)
448 : {
449 : // The two cases are here to avoid numerical roundup errors, for two
450 : // points, we always compute the same interpolation. This condition
451 : // is ensured by the order left->right bottom->top in interpole
452 : // calls
453 : //
454 : // To obtain the same value for border (split) and non-border
455 : // element, we take the middle value and interpolate from this to
456 : // the end
457 182022 : const double xm = .5 * (x1 + x2);
458 182022 : const double ym = .5 * (y1 + y2);
459 182022 : const double fy1 = fudge(y1, level);
460 182022 : const double fym = fudge(ym, level);
461 182022 : if ((fy1 < level && level < fym) || (fy1 > level && level > fym))
462 : {
463 90476 : x2 = xm;
464 90476 : y2 = ym;
465 : }
466 : else
467 : {
468 91546 : x1 = xm;
469 91546 : y1 = ym;
470 : }
471 : }
472 188290 : const double fy1 = fudge(y1, level);
473 188290 : const double ratio = (level - fy1) / (fudge(y2, level) - fy1);
474 188290 : return x1 * (1. - ratio) + x2 * ratio;
475 : }
476 :
477 188290 : Point interpolate(uint8_t border, double level) const
478 : {
479 188290 : switch (border)
480 : {
481 52176 : case LEFT_BORDER:
482 52176 : return Point(upperLeft.x,
483 52176 : interpolate_(level, lowerLeft.y, upperLeft.y,
484 52176 : lowerLeft.value, upperLeft.value,
485 104352 : !split));
486 42028 : case LOWER_BORDER:
487 42028 : return Point(interpolate_(level, lowerLeft.x, lowerRight.x,
488 42028 : lowerLeft.value, lowerRight.value,
489 42028 : !split),
490 84056 : lowerLeft.y);
491 51767 : case RIGHT_BORDER:
492 51767 : return Point(upperRight.x,
493 51767 : interpolate_(level, lowerRight.y, upperRight.y,
494 51767 : lowerRight.value, upperRight.value,
495 103534 : !split));
496 42319 : case UPPER_BORDER:
497 42319 : return Point(interpolate_(level, upperLeft.x, upperRight.x,
498 42319 : upperLeft.value, upperRight.value,
499 42319 : !split),
500 84638 : upperLeft.y);
501 : }
502 0 : assert(false);
503 : return Point();
504 : }
505 : };
506 : } // namespace marching_squares
507 : #endif
|
