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 347244 : Segments(const Segment &first) : sz_(1), segs_()
58 : {
59 86811 : segs_[0] = first;
60 86811 : }
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 182019 : std::size_t size() const
78 : {
79 182019 : return sz_;
80 : }
81 :
82 92415 : const Segment &operator[](std::size_t idx) const
83 : {
84 92415 : assert(idx < sz_);
85 92415 : return segs_[idx];
86 : }
87 :
88 : private:
89 : const std::size_t sz_;
90 : /* const */ Segment segs_[3];
91 : };
92 :
93 825860 : Square(const ValuedPoint &upperLeft_, const ValuedPoint &upperRight_,
94 : const ValuedPoint &lowerLeft_, const ValuedPoint &lowerRight_,
95 : uint8_t borders_ = NO_BORDER, bool split_ = false)
96 825860 : : upperLeft(upperLeft_), lowerLeft(lowerLeft_), lowerRight(lowerRight_),
97 : upperRight(upperRight_),
98 825860 : nanCount((std::isnan(upperLeft.value) ? 1 : 0) +
99 825860 : (std::isnan(upperRight.value) ? 1 : 0) +
100 825860 : (std::isnan(lowerLeft.value) ? 1 : 0) +
101 825860 : (std::isnan(lowerRight.value) ? 1 : 0)),
102 825860 : borders(borders_), split(split_)
103 : {
104 825860 : assert(upperLeft.y == upperRight.y);
105 825860 : assert(lowerLeft.y == lowerRight.y);
106 825860 : assert(lowerLeft.x == upperLeft.x);
107 825860 : assert(lowerRight.x == upperRight.x);
108 825860 : assert(!split || nanCount == 0);
109 825860 : }
110 :
111 10446 : Square upperLeftSquare() const
112 : {
113 10446 : assert(!std::isnan(upperLeft.value));
114 : return Square(
115 20892 : upperLeft, upperCenter(), leftCenter(), center(),
116 20892 : (std::isnan(upperRight.value) ? RIGHT_BORDER : NO_BORDER) |
117 10446 : (std::isnan(lowerLeft.value) ? LOWER_BORDER : NO_BORDER),
118 20892 : true);
119 : }
120 :
121 10445 : Square lowerLeftSquare() const
122 : {
123 10445 : assert(!std::isnan(lowerLeft.value));
124 : return Square(
125 20890 : leftCenter(), center(), lowerLeft, lowerCenter(),
126 20890 : (std::isnan(lowerRight.value) ? RIGHT_BORDER : NO_BORDER) |
127 10445 : (std::isnan(upperLeft.value) ? UPPER_BORDER : NO_BORDER),
128 20890 : true);
129 : }
130 :
131 10447 : Square lowerRightSquare() const
132 : {
133 10447 : assert(!std::isnan(lowerRight.value));
134 : return Square(
135 20894 : center(), rightCenter(), lowerCenter(), lowerRight,
136 20894 : (std::isnan(lowerLeft.value) ? LEFT_BORDER : NO_BORDER) |
137 10447 : (std::isnan(upperRight.value) ? UPPER_BORDER : NO_BORDER),
138 20894 : true);
139 : }
140 :
141 10443 : Square upperRightSquare() const
142 : {
143 10443 : assert(!std::isnan(upperRight.value));
144 : return Square(
145 20886 : upperCenter(), upperRight, center(), rightCenter(),
146 20886 : (std::isnan(lowerRight.value) ? LOWER_BORDER : NO_BORDER) |
147 10443 : (std::isnan(upperLeft.value) ? LEFT_BORDER : NO_BORDER),
148 20886 : true);
149 : }
150 :
151 804788 : double maxValue() const
152 : {
153 804788 : assert(nanCount == 0);
154 804788 : return std::max(std::max(upperLeft.value, upperRight.value),
155 804788 : std::max(lowerLeft.value, lowerRight.value));
156 : }
157 :
158 804788 : double minValue() const
159 : {
160 804788 : assert(nanCount == 0);
161 804788 : return std::min(std::min(upperLeft.value, upperRight.value),
162 804788 : 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 89616 : Segments segments(double level) const
192 : {
193 89616 : 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 9265 : case (UPPER_LEFT):
202 : // debug("UPPER_LEFT");
203 9265 : return Segments(Segment(interpolate(UPPER_BORDER, level),
204 27795 : interpolate(LEFT_BORDER, level)));
205 5808 : case (LOWER_LEFT):
206 : // debug("LOWER_LEFT");
207 5808 : return Segments(Segment(interpolate(LEFT_BORDER, level),
208 17424 : interpolate(LOWER_BORDER, level)));
209 3758 : case (LOWER_RIGHT):
210 : // debug("LOWER_RIGHT");
211 3758 : return Segments(Segment(interpolate(LOWER_BORDER, level),
212 11274 : interpolate(RIGHT_BORDER, level)));
213 4405 : case (UPPER_RIGHT):
214 : // debug("UPPER_RIGHT");
215 4405 : return Segments(Segment(interpolate(RIGHT_BORDER, level),
216 13215 : interpolate(UPPER_BORDER, level)));
217 8975 : case (UPPER_LEFT | LOWER_LEFT):
218 : // debug("UPPER_LEFT | LOWER_LEFT");
219 8975 : return Segments(Segment(interpolate(UPPER_BORDER, level),
220 26925 : interpolate(LOWER_BORDER, level)));
221 12243 : case (LOWER_LEFT | LOWER_RIGHT):
222 : // debug("LOWER_LEFT | LOWER_RIGHT");
223 12243 : return Segments(Segment(interpolate(LEFT_BORDER, level),
224 36729 : interpolate(RIGHT_BORDER, level)));
225 5840 : case (LOWER_RIGHT | UPPER_RIGHT):
226 : // debug("LOWER_RIGHT | UPPER_RIGHT");
227 5840 : return Segments(Segment(interpolate(LOWER_BORDER, level),
228 17520 : interpolate(UPPER_BORDER, level)));
229 12322 : case (UPPER_RIGHT | UPPER_LEFT):
230 : // debug("UPPER_RIGHT | UPPER_LEFT");
231 12322 : return Segments(Segment(interpolate(RIGHT_BORDER, level),
232 36966 : interpolate(LEFT_BORDER, level)));
233 4124 : case (ALL_HIGH & ~UPPER_LEFT):
234 : // debug("ALL_HIGH & ~UPPER_LEFT");
235 4124 : return Segments(Segment(interpolate(LEFT_BORDER, level),
236 12372 : interpolate(UPPER_BORDER, level)));
237 4721 : case (ALL_HIGH & ~LOWER_LEFT):
238 : // debug("ALL_HIGH & ~LOWER_LEFT");
239 4721 : return Segments(Segment(interpolate(LOWER_BORDER, level),
240 14163 : interpolate(LEFT_BORDER, level)));
241 9287 : case (ALL_HIGH & ~LOWER_RIGHT):
242 : // debug("ALL_HIGH & ~LOWER_RIGHT");
243 9287 : return Segments(Segment(interpolate(RIGHT_BORDER, level),
244 27861 : interpolate(LOWER_BORDER, level)));
245 6063 : case (ALL_HIGH & ~UPPER_RIGHT):
246 : // debug("ALL_HIGH & ~UPPER_RIGHT");
247 6063 : return Segments(Segment(interpolate(UPPER_BORDER, level),
248 18189 : interpolate(RIGHT_BORDER, level)));
249 2799 : 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 2799 : return Segments(Segment(interpolate(LEFT_BORDER, level),
264 2799 : interpolate(LOWER_BORDER, level)),
265 2799 : Segment(interpolate(RIGHT_BORDER, level),
266 8397 : interpolate(UPPER_BORDER, level)));
267 : }
268 0 : assert(false);
269 : return Segments();
270 : }
271 :
272 : template <typename Writer, typename LevelGenerator>
273 825848 : void process(const LevelGenerator &levelGenerator, Writer &writer) const
274 : {
275 825848 : if (nanCount == 4) // nothing to do
276 21060 : return;
277 :
278 825847 : if (nanCount) // split in 4
279 : {
280 21059 : if (!std::isnan(upperLeft.value))
281 10443 : upperLeftSquare().process(levelGenerator, writer);
282 21059 : if (!std::isnan(upperRight.value))
283 10443 : upperRightSquare().process(levelGenerator, writer);
284 21059 : if (!std::isnan(lowerLeft.value))
285 10443 : lowerLeftSquare().process(levelGenerator, writer);
286 21059 : if (!std::isnan(lowerRight.value))
287 10446 : lowerRightSquare().process(levelGenerator, writer);
288 21059 : return;
289 : }
290 :
291 804788 : 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 804788 : auto range = levelGenerator.range(minValue(), maxValue());
339 804788 : auto it = range.begin();
340 804788 : auto itEnd = range.end();
341 804788 : auto next = range.begin();
342 804788 : ++next;
343 :
344 894392 : for (; it != itEnd; ++it, ++next)
345 : {
346 89604 : const int levelIdx = (*it).first;
347 89604 : const double level = (*it).second;
348 :
349 89604 : const Segments segments_ = segments(level);
350 :
351 182007 : for (std::size_t i = 0; i < segments_.size(); i++)
352 : {
353 92403 : const Segment &s = segments_[i];
354 92403 : writer.addSegment(levelIdx, s.first, s.second);
355 :
356 92403 : 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 41781 : ValuedPoint center() const
380 : {
381 : return ValuedPoint(
382 41781 : .5 * (upperLeft.x + lowerRight.x),
383 41781 : .5 * (upperLeft.y + lowerRight.y),
384 41781 : ((std::isnan(lowerLeft.value) ? 0 : lowerLeft.value) +
385 41781 : (std::isnan(upperLeft.value) ? 0 : upperLeft.value) +
386 41781 : (std::isnan(lowerRight.value) ? 0 : lowerRight.value) +
387 41781 : (std::isnan(upperRight.value) ? 0 : upperRight.value)) /
388 41781 : (4 - nanCount));
389 : }
390 :
391 20891 : ValuedPoint leftCenter() const
392 : {
393 : return ValuedPoint(
394 20891 : upperLeft.x, .5 * (upperLeft.y + lowerLeft.y),
395 20891 : std::isnan(upperLeft.value)
396 : ? lowerLeft.value
397 15626 : : (std::isnan(lowerLeft.value)
398 15626 : ? upperLeft.value
399 36517 : : .5 * (upperLeft.value + lowerLeft.value)));
400 : }
401 :
402 20892 : ValuedPoint lowerCenter() const
403 : {
404 : return ValuedPoint(
405 20892 : .5 * (lowerLeft.x + lowerRight.x), lowerLeft.y,
406 20892 : std::isnan(lowerRight.value)
407 : ? lowerLeft.value
408 15628 : : (std::isnan(lowerLeft.value)
409 15628 : ? lowerRight.value
410 36520 : : .5 * (lowerRight.value + lowerLeft.value)));
411 : }
412 :
413 20890 : ValuedPoint rightCenter() const
414 : {
415 : return ValuedPoint(
416 20890 : upperRight.x, .5 * (upperRight.y + lowerRight.y),
417 20890 : std::isnan(lowerRight.value)
418 : ? upperRight.value
419 15626 : : (std::isnan(upperRight.value)
420 15626 : ? lowerRight.value
421 36516 : : .5 * (lowerRight.value + upperRight.value)));
422 : }
423 :
424 20889 : ValuedPoint upperCenter() const
425 : {
426 : return ValuedPoint(
427 20889 : .5 * (upperLeft.x + upperRight.x), upperLeft.y,
428 20889 : std::isnan(upperLeft.value)
429 : ? upperRight.value
430 15625 : : (std::isnan(upperRight.value)
431 15625 : ? upperLeft.value
432 36514 : : .5 * (upperLeft.value + upperRight.value)));
433 : }
434 :
435 89616 : uint8_t marchingCase(double level) const
436 : {
437 179232 : return (level < fudge(upperLeft.value, level) ? UPPER_LEFT : ALL_LOW) |
438 179232 : (level < fudge(lowerLeft.value, level) ? LOWER_LEFT : ALL_LOW) |
439 89616 : (level < fudge(lowerRight.value, level) ? LOWER_RIGHT
440 89616 : : ALL_LOW) |
441 179232 : (level < fudge(upperRight.value, level) ? UPPER_RIGHT : ALL_LOW);
442 : }
443 :
444 184982 : static double interpolate_(double level, double x1, double x2, double y1,
445 : double y2, bool need_split)
446 : {
447 184982 : 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 178986 : const double xm = .5 * (x1 + x2);
458 178986 : const double ym = .5 * (y1 + y2);
459 178986 : const double fy1 = fudge(y1, level);
460 178986 : const double fym = fudge(ym, level);
461 178986 : if ((fy1 < level && level < fym) || (fy1 > level && level > fym))
462 : {
463 88933 : x2 = xm;
464 88933 : y2 = ym;
465 : }
466 : else
467 : {
468 90053 : x1 = xm;
469 90053 : y1 = ym;
470 : }
471 : }
472 184982 : const double fy1 = fudge(y1, level);
473 184982 : const double ratio = (level - fy1) / (fudge(y2, level) - fy1);
474 184982 : return x1 * (1. - ratio) + x2 * ratio;
475 : }
476 :
477 184982 : Point interpolate(uint8_t border, double level) const
478 : {
479 184982 : switch (border)
480 : {
481 51341 : case LEFT_BORDER:
482 51341 : return Point(upperLeft.x,
483 51341 : interpolate_(level, lowerLeft.y, upperLeft.y,
484 51341 : lowerLeft.value, upperLeft.value,
485 102682 : !split));
486 41211 : case LOWER_BORDER:
487 41211 : return Point(interpolate_(level, lowerLeft.x, lowerRight.x,
488 41211 : lowerLeft.value, lowerRight.value,
489 41211 : !split),
490 82422 : lowerLeft.y);
491 50934 : case RIGHT_BORDER:
492 50934 : return Point(upperRight.x,
493 50934 : interpolate_(level, lowerRight.y, upperRight.y,
494 50934 : lowerRight.value, upperRight.value,
495 101868 : !split));
496 41496 : case UPPER_BORDER:
497 41496 : return Point(interpolate_(level, upperLeft.x, upperRight.x,
498 41496 : upperLeft.value, upperRight.value,
499 41496 : !split),
500 82992 : upperLeft.y);
501 : }
502 0 : assert(false);
503 : return Point();
504 : }
505 : };
506 : } // namespace marching_squares
507 : #endif
|
