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 : * Permission is hereby granted, free of charge, to any person obtaining a
11 : * copy of this software and associated documentation files (the "Software"),
12 : * to deal in the Software without restriction, including without limitation
13 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
14 : * and/or sell copies of the Software, and to permit persons to whom the
15 : * Software is furnished to do so, subject to the following conditions:
16 : *
17 : * The above copyright notice and this permission notice shall be included
18 : * in all copies or substantial portions of the Software.
19 : *
20 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
21 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
23 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
24 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
25 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
26 : * DEALINGS IN THE SOFTWARE.
27 : ****************************************************************************/
28 : #ifndef MARCHING_SQUARE_SQUARE_H
29 : #define MARCHING_SQUARE_SQUARE_H
30 :
31 : #include <algorithm>
32 : #include <cassert>
33 : #include <cstdint>
34 : #include <math.h>
35 : #include "utility.h"
36 : #include "point.h"
37 :
38 : namespace marching_squares
39 : {
40 :
41 : struct Square
42 : {
43 : // Bit flags to determine borders around pixel
44 : static const uint8_t NO_BORDER = 0; // 0000 0000
45 : static const uint8_t LEFT_BORDER = 1 << 0; // 0000 0001
46 : static const uint8_t LOWER_BORDER = 1 << 1; // 0000 0010
47 : static const uint8_t RIGHT_BORDER = 1 << 2; // 0000 0100
48 : static const uint8_t UPPER_BORDER = 1 << 3; // 0000 1000
49 :
50 : // Bit flags for marching square case
51 : static const uint8_t ALL_LOW = 0; // 0000 0000
52 : static const uint8_t UPPER_LEFT = 1 << 0; // 0000 0001
53 : static const uint8_t LOWER_LEFT = 1 << 1; // 0000 0010
54 : static const uint8_t LOWER_RIGHT = 1 << 2; // 0000 0100
55 : static const uint8_t UPPER_RIGHT = 1 << 3; // 0000 1000
56 : static const uint8_t ALL_HIGH =
57 : UPPER_LEFT | LOWER_LEFT | LOWER_RIGHT | UPPER_RIGHT; // 0000 1111
58 : static const uint8_t SADDLE_NW = UPPER_LEFT | LOWER_RIGHT; // 0000 0101
59 : static const uint8_t SADDLE_NE = UPPER_RIGHT | LOWER_LEFT; // 0000 1010
60 :
61 : typedef std::pair<Point, Point> Segment;
62 : typedef std::pair<ValuedPoint, ValuedPoint> ValuedSegment;
63 :
64 : //
65 : // An array of segments, at most 3 segments
66 : struct Segments
67 : {
68 24 : Segments() : sz_(0)
69 : {
70 6 : }
71 :
72 326280 : Segments(const Segment &first) : sz_(1), segs_()
73 : {
74 81570 : segs_[0] = first;
75 81570 : }
76 :
77 9436 : Segments(const Segment &first, const Segment &second) : sz_(2), segs_()
78 : {
79 2359 : segs_[0] = first;
80 2359 : segs_[1] = second;
81 2359 : }
82 :
83 : Segments(const Segment &first, const Segment &second,
84 : const Segment &third)
85 : : sz_(3), segs_()
86 : {
87 : segs_[0] = first;
88 : segs_[1] = second;
89 : segs_[2] = third;
90 : }
91 :
92 170217 : std::size_t size() const
93 : {
94 170217 : return sz_;
95 : }
96 :
97 86294 : const Segment &operator[](std::size_t idx) const
98 : {
99 86294 : assert(idx < sz_);
100 86294 : return segs_[idx];
101 : }
102 :
103 : private:
104 : const std::size_t sz_;
105 : /* const */ Segment segs_[3];
106 : };
107 :
108 418365 : Square(const ValuedPoint &upperLeft_, const ValuedPoint &upperRight_,
109 : const ValuedPoint &lowerLeft_, const ValuedPoint &lowerRight_,
110 : uint8_t borders_ = NO_BORDER, bool split_ = false)
111 418365 : : upperLeft(upperLeft_), lowerLeft(lowerLeft_), lowerRight(lowerRight_),
112 : upperRight(upperRight_),
113 418365 : nanCount((std::isnan(upperLeft.value) ? 1 : 0) +
114 418365 : (std::isnan(upperRight.value) ? 1 : 0) +
115 418365 : (std::isnan(lowerLeft.value) ? 1 : 0) +
116 418365 : (std::isnan(lowerRight.value) ? 1 : 0)),
117 418365 : borders(borders_), split(split_)
118 : {
119 418365 : assert(upperLeft.y == upperRight.y);
120 418365 : assert(lowerLeft.y == lowerRight.y);
121 418365 : assert(lowerLeft.x == upperLeft.x);
122 418365 : assert(lowerRight.x == upperRight.x);
123 418365 : assert(!split || nanCount == 0);
124 418365 : }
125 :
126 5234 : Square upperLeftSquare() const
127 : {
128 5234 : assert(!std::isnan(upperLeft.value));
129 : return Square(
130 10468 : upperLeft, upperCenter(), leftCenter(), center(),
131 10468 : (std::isnan(upperRight.value) ? RIGHT_BORDER : NO_BORDER) |
132 5234 : (std::isnan(lowerLeft.value) ? LOWER_BORDER : NO_BORDER),
133 10468 : true);
134 : }
135 :
136 5236 : Square lowerLeftSquare() const
137 : {
138 5236 : assert(!std::isnan(lowerLeft.value));
139 : return Square(
140 10472 : leftCenter(), center(), lowerLeft, lowerCenter(),
141 10472 : (std::isnan(lowerRight.value) ? RIGHT_BORDER : NO_BORDER) |
142 5236 : (std::isnan(upperLeft.value) ? UPPER_BORDER : NO_BORDER),
143 10472 : true);
144 : }
145 :
146 5241 : Square lowerRightSquare() const
147 : {
148 5241 : assert(!std::isnan(lowerRight.value));
149 : return Square(
150 10482 : center(), rightCenter(), lowerCenter(), lowerRight,
151 10482 : (std::isnan(lowerLeft.value) ? LEFT_BORDER : NO_BORDER) |
152 5241 : (std::isnan(upperRight.value) ? UPPER_BORDER : NO_BORDER),
153 10482 : true);
154 : }
155 :
156 5231 : Square upperRightSquare() const
157 : {
158 5231 : assert(!std::isnan(upperRight.value));
159 : return Square(
160 10462 : upperCenter(), upperRight, center(), rightCenter(),
161 10462 : (std::isnan(lowerRight.value) ? LOWER_BORDER : NO_BORDER) |
162 5231 : (std::isnan(upperLeft.value) ? LEFT_BORDER : NO_BORDER),
163 10462 : true);
164 : }
165 :
166 407802 : double maxValue() const
167 : {
168 407802 : assert(nanCount == 0);
169 407802 : return std::max(std::max(upperLeft.value, upperRight.value),
170 407802 : std::max(lowerLeft.value, lowerRight.value));
171 : }
172 :
173 407802 : double minValue() const
174 : {
175 407802 : assert(nanCount == 0);
176 407802 : return std::min(std::min(upperLeft.value, upperRight.value),
177 407802 : std::min(lowerLeft.value, lowerRight.value));
178 : }
179 :
180 10473 : ValuedSegment segment(uint8_t border) const
181 : {
182 10473 : switch (border)
183 : {
184 2617 : case LEFT_BORDER:
185 2617 : return ValuedSegment(upperLeft, lowerLeft);
186 2620 : case LOWER_BORDER:
187 2620 : return ValuedSegment(lowerLeft, lowerRight);
188 2615 : case RIGHT_BORDER:
189 2615 : return ValuedSegment(lowerRight, upperRight);
190 2621 : case UPPER_BORDER:
191 2621 : return ValuedSegment(upperRight, upperLeft);
192 : }
193 0 : assert(false);
194 : return ValuedSegment(upperLeft, upperLeft);
195 : }
196 :
197 : // returns segments of contour
198 : //
199 : // segments are oriented:
200 : // - they form a vector from their first point to their second point.
201 : // - when looking at the vector upward, values greater than the level are
202 : // on the right
203 : //
204 : // ^
205 : // - | +
206 83935 : Segments segments(double level) const
207 : {
208 83935 : switch (marchingCase(level))
209 : {
210 1 : case (ALL_LOW):
211 : // debug("ALL_LOW");
212 1 : return Segments();
213 5 : case (ALL_HIGH):
214 : // debug("ALL_HIGH");
215 5 : return Segments();
216 7633 : case (UPPER_LEFT):
217 : // debug("UPPER_LEFT");
218 7633 : return Segments(Segment(interpolate(UPPER_BORDER, level),
219 22899 : interpolate(LEFT_BORDER, level)));
220 4970 : case (LOWER_LEFT):
221 : // debug("LOWER_LEFT");
222 4970 : return Segments(Segment(interpolate(LEFT_BORDER, level),
223 14910 : interpolate(LOWER_BORDER, level)));
224 3089 : case (LOWER_RIGHT):
225 : // debug("LOWER_RIGHT");
226 3089 : return Segments(Segment(interpolate(LOWER_BORDER, level),
227 9267 : interpolate(RIGHT_BORDER, level)));
228 3849 : case (UPPER_RIGHT):
229 : // debug("UPPER_RIGHT");
230 3849 : return Segments(Segment(interpolate(RIGHT_BORDER, level),
231 11547 : interpolate(UPPER_BORDER, level)));
232 20343 : case (UPPER_LEFT | LOWER_LEFT):
233 : // debug("UPPER_LEFT | LOWER_LEFT");
234 20343 : return Segments(Segment(interpolate(UPPER_BORDER, level),
235 61029 : interpolate(LOWER_BORDER, level)));
236 8790 : case (LOWER_LEFT | LOWER_RIGHT):
237 : // debug("LOWER_LEFT | LOWER_RIGHT");
238 8790 : return Segments(Segment(interpolate(LEFT_BORDER, level),
239 26370 : interpolate(RIGHT_BORDER, level)));
240 3706 : case (LOWER_RIGHT | UPPER_RIGHT):
241 : // debug("LOWER_RIGHT | UPPER_RIGHT");
242 3706 : return Segments(Segment(interpolate(LOWER_BORDER, level),
243 11118 : interpolate(UPPER_BORDER, level)));
244 9164 : case (UPPER_RIGHT | UPPER_LEFT):
245 : // debug("UPPER_RIGHT | UPPER_LEFT");
246 9164 : return Segments(Segment(interpolate(RIGHT_BORDER, level),
247 27492 : interpolate(LEFT_BORDER, level)));
248 3308 : case (ALL_HIGH & ~UPPER_LEFT):
249 : // debug("ALL_HIGH & ~UPPER_LEFT");
250 3308 : return Segments(Segment(interpolate(LEFT_BORDER, level),
251 9924 : interpolate(UPPER_BORDER, level)));
252 4073 : case (ALL_HIGH & ~LOWER_LEFT):
253 : // debug("ALL_HIGH & ~LOWER_LEFT");
254 4073 : return Segments(Segment(interpolate(LOWER_BORDER, level),
255 12219 : interpolate(LEFT_BORDER, level)));
256 7521 : case (ALL_HIGH & ~LOWER_RIGHT):
257 : // debug("ALL_HIGH & ~LOWER_RIGHT");
258 7521 : return Segments(Segment(interpolate(RIGHT_BORDER, level),
259 22563 : interpolate(LOWER_BORDER, level)));
260 5124 : case (ALL_HIGH & ~UPPER_RIGHT):
261 : // debug("ALL_HIGH & ~UPPER_RIGHT");
262 5124 : return Segments(Segment(interpolate(UPPER_BORDER, level),
263 15372 : interpolate(RIGHT_BORDER, level)));
264 2359 : case (SADDLE_NE):
265 : case (SADDLE_NW):
266 : // From the two possible saddle configurations, we always return
267 : // the same one.
268 : //
269 : // The classical marching square algorithm says the ambiguity
270 : // should be resolved between the two possible configurations by
271 : // looking at the value of the center point. But in certain
272 : // cases, this may lead to line contours from different levels
273 : // that cross each other and then gives invalid polygons.
274 : //
275 : // Arbitrarily choosing one of the two possible configurations
276 : // is not really that worse than deciding based on the center
277 : // point.
278 2359 : return Segments(Segment(interpolate(LEFT_BORDER, level),
279 2359 : interpolate(LOWER_BORDER, level)),
280 2359 : Segment(interpolate(RIGHT_BORDER, level),
281 7077 : interpolate(UPPER_BORDER, level)));
282 : }
283 0 : assert(false);
284 : return Segments();
285 : }
286 :
287 : template <typename Writer, typename LevelGenerator>
288 418353 : void process(const LevelGenerator &levelGenerator, Writer &writer) const
289 : {
290 418353 : if (nanCount == 4) // nothing to do
291 10547 : return;
292 :
293 418348 : if (nanCount) // split in 4
294 : {
295 10546 : if (!std::isnan(upperLeft.value))
296 5231 : upperLeftSquare().process(levelGenerator, writer);
297 10546 : if (!std::isnan(upperRight.value))
298 5231 : upperRightSquare().process(levelGenerator, writer);
299 10546 : if (!std::isnan(lowerLeft.value))
300 5234 : lowerLeftSquare().process(levelGenerator, writer);
301 10545 : if (!std::isnan(lowerRight.value))
302 5240 : lowerRightSquare().process(levelGenerator, writer);
303 10542 : return;
304 : }
305 :
306 407802 : if (writer.polygonize && borders)
307 : {
308 51890 : for (uint8_t border :
309 : {UPPER_BORDER, LEFT_BORDER, RIGHT_BORDER, LOWER_BORDER})
310 : {
311 : // bitwise AND to test which borders we have on the square
312 41512 : if ((border & borders) == 0)
313 31039 : continue;
314 :
315 : // convention: for a level = L, store borders for the previous
316 : // level up to (and including) L in the border of level "L". For
317 : // fixed sets of level, this means there is an "Inf" slot for
318 : // borders of the highest level
319 10473 : const ValuedSegment s(segment(border));
320 :
321 10473 : Point lastPoint(s.first.x, s.first.y);
322 10473 : Point endPoint(s.second.x, s.second.y);
323 10473 : if (s.first.value > s.second.value)
324 34 : std::swap(lastPoint, endPoint);
325 10473 : bool reverse =
326 10497 : (s.first.value > s.second.value) &&
327 24 : ((border == UPPER_BORDER) || (border == LEFT_BORDER));
328 :
329 10473 : auto levelIt =
330 10473 : levelGenerator.range(s.first.value, s.second.value);
331 :
332 10473 : auto it = levelIt.begin(); // reused after the for
333 10503 : for (; it != levelIt.end(); ++it)
334 : {
335 30 : const int levelIdx = (*it).first;
336 30 : const double level = (*it).second;
337 :
338 30 : const Point nextPoint(interpolate(border, level));
339 30 : if (reverse)
340 8 : writer.addBorderSegment(levelIdx, nextPoint, lastPoint);
341 : else
342 22 : writer.addBorderSegment(levelIdx, lastPoint, nextPoint);
343 30 : lastPoint = nextPoint;
344 : }
345 : // last level (past the end)
346 10473 : if (reverse)
347 20 : writer.addBorderSegment((*it).first, endPoint, lastPoint);
348 : else
349 10453 : writer.addBorderSegment((*it).first, lastPoint, endPoint);
350 : }
351 : }
352 :
353 407802 : auto range = levelGenerator.range(minValue(), maxValue());
354 407798 : auto it = range.begin();
355 407798 : auto itEnd = range.end();
356 407798 : auto next = range.begin();
357 407798 : ++next;
358 :
359 491721 : for (; it != itEnd; ++it, ++next)
360 : {
361 83923 : const int levelIdx = (*it).first;
362 83923 : const double level = (*it).second;
363 :
364 83923 : const Segments segments_ = segments(level);
365 :
366 170205 : for (std::size_t i = 0; i < segments_.size(); i++)
367 : {
368 86282 : const Segment &s = segments_[i];
369 86282 : writer.addSegment(levelIdx, s.first, s.second);
370 :
371 86282 : 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 4463 : 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 20942 : ValuedPoint center() const
395 : {
396 : return ValuedPoint(
397 20942 : .5 * (upperLeft.x + lowerRight.x),
398 20942 : .5 * (upperLeft.y + lowerRight.y),
399 20942 : ((std::isnan(lowerLeft.value) ? 0 : lowerLeft.value) +
400 20942 : (std::isnan(upperLeft.value) ? 0 : upperLeft.value) +
401 20942 : (std::isnan(lowerRight.value) ? 0 : lowerRight.value) +
402 20942 : (std::isnan(upperRight.value) ? 0 : upperRight.value)) /
403 20942 : (4 - nanCount));
404 : }
405 :
406 10470 : ValuedPoint leftCenter() const
407 : {
408 : return ValuedPoint(
409 10470 : upperLeft.x, .5 * (upperLeft.y + lowerLeft.y),
410 10470 : std::isnan(upperLeft.value)
411 : ? lowerLeft.value
412 7832 : : (std::isnan(lowerLeft.value)
413 7832 : ? upperLeft.value
414 18302 : : .5 * (upperLeft.value + lowerLeft.value)));
415 : }
416 :
417 10477 : ValuedPoint lowerCenter() const
418 : {
419 : return ValuedPoint(
420 10477 : .5 * (lowerLeft.x + lowerRight.x), lowerLeft.y,
421 10477 : std::isnan(lowerRight.value)
422 : ? lowerLeft.value
423 7843 : : (std::isnan(lowerLeft.value)
424 7843 : ? lowerRight.value
425 18320 : : .5 * (lowerRight.value + lowerLeft.value)));
426 : }
427 :
428 10472 : ValuedPoint rightCenter() const
429 : {
430 : return ValuedPoint(
431 10472 : upperRight.x, .5 * (upperRight.y + lowerRight.y),
432 10472 : std::isnan(lowerRight.value)
433 : ? upperRight.value
434 7838 : : (std::isnan(upperRight.value)
435 7838 : ? lowerRight.value
436 18310 : : .5 * (lowerRight.value + upperRight.value)));
437 : }
438 :
439 10465 : ValuedPoint upperCenter() const
440 : {
441 : return ValuedPoint(
442 10465 : .5 * (upperLeft.x + upperRight.x), upperLeft.y,
443 10465 : std::isnan(upperLeft.value)
444 : ? upperRight.value
445 7831 : : (std::isnan(upperRight.value)
446 7831 : ? upperLeft.value
447 18296 : : .5 * (upperLeft.value + upperRight.value)));
448 : }
449 :
450 83935 : uint8_t marchingCase(double level) const
451 : {
452 167870 : return (level < fudge(level, upperLeft.value) ? UPPER_LEFT : ALL_LOW) |
453 167870 : (level < fudge(level, lowerLeft.value) ? LOWER_LEFT : ALL_LOW) |
454 83935 : (level < fudge(level, lowerRight.value) ? LOWER_RIGHT
455 83935 : : ALL_LOW) |
456 167870 : (level < fudge(level, upperRight.value) ? UPPER_RIGHT : ALL_LOW);
457 : }
458 :
459 172606 : static double interpolate_(double level, double x1, double x2, double y1,
460 : double y2, bool need_split)
461 : {
462 172606 : if (need_split)
463 : {
464 : // The two cases are here to avoid numerical roundup errors, for two
465 : // points, we always compute the same interpolation. This condition
466 : // is ensured by the order left->right bottom->top in interpole
467 : // calls
468 : //
469 : // To obtain the same value for border (split) and non-border
470 : // element, we take the middle value and interpolate from this to
471 : // the end
472 140088 : const double xm = .5 * (x1 + x2);
473 140088 : const double ym = .5 * (y1 + y2);
474 140088 : const double fy1 = fudge(level, y1);
475 140088 : const double fym = fudge(level, ym);
476 140088 : if ((fy1 < level && level < fym) || (fy1 > level && level > fym))
477 : {
478 69847 : x2 = xm;
479 69847 : y2 = ym;
480 : }
481 : else
482 : {
483 70241 : x1 = xm;
484 70241 : y1 = ym;
485 : }
486 : }
487 172606 : const double fy1 = fudge(level, y1);
488 172606 : const double ratio = (level - fy1) / (fudge(level, y2) - fy1);
489 172606 : return x1 * (1. - ratio) + x2 * ratio;
490 : }
491 :
492 172606 : Point interpolate(uint8_t border, double level) const
493 : {
494 172606 : switch (border)
495 : {
496 40305 : case LEFT_BORDER:
497 40305 : return Point(upperLeft.x,
498 40305 : interpolate_(level, lowerLeft.y, upperLeft.y,
499 40305 : lowerLeft.value, upperLeft.value,
500 80610 : !split));
501 46071 : case LOWER_BORDER:
502 46071 : return Point(interpolate_(level, lowerLeft.x, lowerRight.x,
503 46071 : lowerLeft.value, lowerRight.value,
504 46071 : !split),
505 92142 : lowerLeft.y);
506 39900 : case RIGHT_BORDER:
507 39900 : return Point(upperRight.x,
508 39900 : interpolate_(level, lowerRight.y, upperRight.y,
509 39900 : lowerRight.value, upperRight.value,
510 79800 : !split));
511 46330 : case UPPER_BORDER:
512 46330 : return Point(interpolate_(level, upperLeft.x, upperRight.x,
513 46330 : upperLeft.value, upperRight.value,
514 46330 : !split),
515 92660 : upperLeft.y);
516 : }
517 0 : assert(false);
518 : return Point();
519 : }
520 : };
521 : } // namespace marching_squares
522 : #endif
|
