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 <cassert> 32 #include <cstdint> 33 #include <math.h> 34 #include "utility.h" 35 #include "point.h" 36 37 namespace marching_squares { 38 39 struct Square 40 { 41 // Bit flags to determine borders around pixel 42 static const uint8_t NO_BORDER = 0; // 0000 0000 43 static const uint8_t LEFT_BORDER = 1 << 0; // 0000 0001 44 static const uint8_t LOWER_BORDER = 1 << 1; // 0000 0010 45 static const uint8_t RIGHT_BORDER = 1 << 2; // 0000 0100 46 static const uint8_t UPPER_BORDER = 1 << 3; // 0000 1000 47 48 // Bit flags for marching square case 49 static const uint8_t ALL_LOW = 0; // 0000 0000 50 static const uint8_t UPPER_LEFT = 1 << 0; // 0000 0001 51 static const uint8_t LOWER_LEFT = 1 << 1; // 0000 0010 52 static const uint8_t LOWER_RIGHT = 1 << 2; // 0000 0100 53 static const uint8_t UPPER_RIGHT = 1 << 3; // 0000 1000 54 static const uint8_t ALL_HIGH = UPPER_LEFT | LOWER_LEFT | LOWER_RIGHT | UPPER_RIGHT; // 0000 1111 55 static const uint8_t SADDLE_NW = UPPER_LEFT | LOWER_RIGHT; // 0000 0101 56 static const uint8_t SADDLE_NE = UPPER_RIGHT | LOWER_LEFT; // 0000 1010 57 58 typedef std::pair< Point, Point > Segment; 59 typedef std::pair< ValuedPoint, ValuedPoint > ValuedSegment; 60 61 // 62 // An array of segments, at most 3 segments 63 struct Segments 64 { SegmentsSquare::Segments65 Segments(): sz_(0) {} SegmentsSquare::Segments66 Segments(const Segment &first): sz_(1), segs_() 67 { 68 segs_[0] = first; 69 } SegmentsSquare::Segments70 Segments(const Segment &first, const Segment &second): sz_(2), segs_() 71 { 72 segs_[0] = first; 73 segs_[1] = second; 74 } SegmentsSquare::Segments75 Segments(const Segment &first, const Segment &second, const Segment &third): sz_(3), segs_() 76 { 77 segs_[0] = first; 78 segs_[1] = second; 79 segs_[2] = third; 80 } 81 sizeSquare::Segments82 std::size_t size() const {return sz_;} 83 84 const Segment & operator[](std::size_t idx) const 85 { 86 assert(idx < sz_); 87 return segs_[idx]; 88 } 89 private: 90 const std::size_t sz_; 91 /* const */ Segment segs_[3]; 92 }; 93 94 95 Square(const ValuedPoint & upperLeft_, const ValuedPoint & upperRight_, 96 const ValuedPoint & lowerLeft_, const ValuedPoint & lowerRight_, uint8_t borders_ = NO_BORDER, bool split_ = false) upperLeftSquare97 : upperLeft(upperLeft_) 98 , lowerLeft(lowerLeft_) 99 , lowerRight(lowerRight_) 100 , upperRight(upperRight_) 101 , nanCount((std::isnan(upperLeft.value) ? 1 : 0) + (std::isnan(upperRight.value) ? 1 : 0) 102 + (std::isnan(lowerLeft.value) ? 1 : 0) + (std::isnan(lowerRight.value) ? 1 : 0)) 103 , borders(borders_) 104 , split(split_) 105 { 106 assert(upperLeft.y == upperRight.y); 107 assert(lowerLeft.y == lowerRight.y); 108 assert(lowerLeft.x == upperLeft.x); 109 assert(lowerRight.x == upperRight.x); 110 assert(!split || nanCount == 0); 111 } 112 upperLeftSquareSquare113 Square upperLeftSquare() const 114 { 115 assert(!std::isnan(upperLeft.value)); 116 return Square( 117 upperLeft, upperCenter(), 118 leftCenter(), center(), 119 (std::isnan(upperRight.value) ? RIGHT_BORDER: NO_BORDER) 120 | (std::isnan(lowerLeft.value) ? LOWER_BORDER : NO_BORDER), true); 121 } 122 lowerLeftSquareSquare123 Square lowerLeftSquare() const 124 { 125 assert(!std::isnan(lowerLeft.value)); 126 return Square( 127 leftCenter(), center(), 128 lowerLeft, lowerCenter(), 129 (std::isnan(lowerRight.value) ? RIGHT_BORDER: NO_BORDER) 130 | (std::isnan(upperLeft.value) ? UPPER_BORDER : NO_BORDER), true); 131 } 132 lowerRightSquareSquare133 Square lowerRightSquare() const 134 { 135 assert(!std::isnan(lowerRight.value)); 136 return Square( 137 center(), rightCenter(), 138 lowerCenter(), lowerRight, 139 (std::isnan(lowerLeft.value) ? LEFT_BORDER: NO_BORDER) 140 | (std::isnan(upperRight.value) ? UPPER_BORDER : NO_BORDER), true); 141 } 142 upperRightSquareSquare143 Square upperRightSquare() const 144 { 145 assert(!std::isnan(upperRight.value)); 146 return Square( 147 upperCenter(), upperRight, 148 center(), rightCenter(), 149 (std::isnan(lowerRight.value) ? LOWER_BORDER: NO_BORDER) 150 | (std::isnan(upperLeft.value) ? LEFT_BORDER : NO_BORDER), true); 151 } 152 maxValueSquare153 double maxValue() const 154 { 155 assert(nanCount==0); 156 return std::max(std::max(upperLeft.value, upperRight.value), 157 std::max(lowerLeft.value, lowerRight.value)); 158 } 159 minValueSquare160 double minValue() const 161 { 162 assert(nanCount==0); 163 return std::min(std::min(upperLeft.value, upperRight.value), 164 std::min(lowerLeft.value, lowerRight.value)); 165 } 166 segmentSquare167 ValuedSegment segment(uint8_t border) const 168 { 169 switch(border) 170 { 171 case LEFT_BORDER: return ValuedSegment(upperLeft, lowerLeft); 172 case LOWER_BORDER: return ValuedSegment(lowerLeft, lowerRight); 173 case RIGHT_BORDER: return ValuedSegment(lowerRight, upperRight); 174 case UPPER_BORDER: return ValuedSegment(upperRight, upperLeft); 175 } 176 assert(false); 177 return ValuedSegment(upperLeft, upperLeft); 178 } 179 180 // returns segments of contour 181 // 182 // segments are oriented: 183 // - they form a vector from their first point to their second point. 184 // - when looking at the vector upward, values greater than the level are on the right 185 // 186 // ^ 187 // - | + segmentsSquare188 Segments segments(double level) const 189 { 190 switch (marchingCase(level)) 191 { 192 case (ALL_LOW): 193 //debug("ALL_LOW"); 194 return Segments(); 195 case (ALL_HIGH): 196 //debug("ALL_HIGH"); 197 return Segments(); 198 case (UPPER_LEFT): 199 //debug("UPPER_LEFT"); 200 return Segments(Segment(interpolate(UPPER_BORDER, level), interpolate(LEFT_BORDER, level))); 201 case (LOWER_LEFT): 202 //debug("LOWER_LEFT"); 203 return Segments(Segment(interpolate(LEFT_BORDER, level), interpolate(LOWER_BORDER, level))); 204 case (LOWER_RIGHT): 205 //debug("LOWER_RIGHT"); 206 return Segments(Segment(interpolate(LOWER_BORDER, level), interpolate(RIGHT_BORDER, level))); 207 case (UPPER_RIGHT): 208 //debug("UPPER_RIGHT"); 209 return Segments(Segment(interpolate(RIGHT_BORDER, level), interpolate(UPPER_BORDER, level))); 210 case (UPPER_LEFT | LOWER_LEFT): 211 //debug("UPPER_LEFT | LOWER_LEFT"); 212 return Segments(Segment(interpolate(UPPER_BORDER, level), interpolate(LOWER_BORDER, level))); 213 case (LOWER_LEFT | LOWER_RIGHT): 214 //debug("LOWER_LEFT | LOWER_RIGHT"); 215 return Segments(Segment(interpolate(LEFT_BORDER, level), interpolate(RIGHT_BORDER, level))); 216 case (LOWER_RIGHT | UPPER_RIGHT): 217 //debug("LOWER_RIGHT | UPPER_RIGHT"); 218 return Segments(Segment(interpolate(LOWER_BORDER, level), interpolate(UPPER_BORDER, level))); 219 case (UPPER_RIGHT | UPPER_LEFT): 220 //debug("UPPER_RIGHT | UPPER_LEFT"); 221 return Segments(Segment(interpolate(RIGHT_BORDER, level), interpolate(LEFT_BORDER, level))); 222 case (ALL_HIGH & ~UPPER_LEFT): 223 //debug("ALL_HIGH & ~UPPER_LEFT"); 224 return Segments(Segment(interpolate(LEFT_BORDER, level), interpolate(UPPER_BORDER, level))); 225 case (ALL_HIGH & ~LOWER_LEFT): 226 //debug("ALL_HIGH & ~LOWER_LEFT"); 227 return Segments(Segment(interpolate(LOWER_BORDER, level), interpolate(LEFT_BORDER, level))); 228 case (ALL_HIGH & ~LOWER_RIGHT): 229 //debug("ALL_HIGH & ~LOWER_RIGHT"); 230 return Segments(Segment(interpolate(RIGHT_BORDER, level), interpolate(LOWER_BORDER, level))); 231 case (ALL_HIGH & ~UPPER_RIGHT): 232 //debug("ALL_HIGH & ~UPPER_RIGHT"); 233 return Segments(Segment(interpolate(UPPER_BORDER, level), interpolate(RIGHT_BORDER, level))); 234 case (SADDLE_NE): 235 case (SADDLE_NW): 236 // From the two possible saddle configurations, we always return the same one. 237 // 238 // The classical marching square algorithm says the ambiguity should be resolved between the two 239 // possible configurations by looking at the value of the center point. 240 // But in certain cases, this may lead to line contours from different levels that cross 241 // each other and then gives invalid polygons. 242 // 243 // Arbitrarily choosing one of the two possible configurations is not really that worse than 244 // deciding based on the center point. 245 return Segments( 246 Segment(interpolate(LEFT_BORDER, level), interpolate(LOWER_BORDER, level)), 247 Segment(interpolate(RIGHT_BORDER, level), interpolate(UPPER_BORDER, level))); 248 } 249 assert(false); 250 return Segments(); 251 } 252 253 template <typename Writer, typename LevelGenerator> processSquare254 void process(const LevelGenerator &levelGenerator, Writer & writer) const 255 { 256 if (nanCount == 4) // nothing to do 257 return; 258 259 260 if (nanCount) // split in 4 261 { 262 if (!std::isnan(upperLeft.value)) 263 upperLeftSquare().process(levelGenerator, writer); 264 if (!std::isnan(upperRight.value)) 265 upperRightSquare().process(levelGenerator, writer); 266 if (!std::isnan(lowerLeft.value)) 267 lowerLeftSquare().process(levelGenerator, writer); 268 if (!std::isnan(lowerRight.value)) 269 lowerRightSquare().process(levelGenerator, writer); 270 return; 271 } 272 273 if ( writer.polygonize && borders ) 274 { 275 for ( uint8_t border : {UPPER_BORDER, LEFT_BORDER, RIGHT_BORDER, LOWER_BORDER} ) 276 { 277 // bitwise AND to test which borders we have on the square 278 if ( ( border & borders ) == 0 ) 279 continue; 280 281 // convention: for a level = L, store borders for the previous level up to 282 // (and including) L in the border of level "L". 283 // For fixed sets of level, this means there is an "Inf" slot for borders of the highest level 284 const ValuedSegment s(segment( border )); 285 286 Point lastPoint(s.first.x, s.first.y); 287 Point endPoint(s.second.x, s.second.y); 288 if ( s.first.value > s.second.value ) 289 std::swap( lastPoint, endPoint ); 290 bool reverse = (s.first.value > s.second.value) && 291 ( (border == UPPER_BORDER) || (border == LEFT_BORDER) ); 292 293 auto levelIt = levelGenerator.range(s.first.value, s.second.value); 294 295 auto it = levelIt.begin(); // reused after the for 296 for ( ; it != levelIt.end(); ++it ) 297 { 298 const int levelIdx = (*it).first; 299 const double level = (*it).second; 300 301 const Point nextPoint(interpolate(border, level)); 302 if ( reverse ) 303 writer.addBorderSegment( levelIdx, nextPoint, lastPoint ); 304 else 305 writer.addBorderSegment( levelIdx, lastPoint, nextPoint ); 306 lastPoint = nextPoint; 307 } 308 // last level (past the end) 309 if ( reverse ) 310 writer.addBorderSegment( (*it).first, endPoint, lastPoint ); 311 else 312 writer.addBorderSegment( (*it).first, lastPoint, endPoint ); 313 } 314 } 315 316 auto range = levelGenerator.range( minValue(), maxValue() ); 317 auto it = range.begin(); 318 auto itEnd = range.end(); 319 auto next = range.begin(); ++next; 320 321 for ( ; it != itEnd; ++it, ++next ) { 322 const int levelIdx = (*it).first; 323 const double level = (*it).second; 324 325 const Segments segments_ = segments( level ); 326 327 for (std::size_t i=0; i<segments_.size(); i++) 328 { 329 const Segment &s = segments_[i]; 330 writer.addSegment( levelIdx, s.first, s.second ); 331 332 if ( writer.polygonize ) { 333 // the contour is used in the polygon of higher level as well 334 // 335 // TODO: copying the segment to the higher level is easy, 336 // but it involves too much memory. We should reuse segment 337 // contours when constructing polygon rings. 338 writer.addSegment( (*next).first, s.first, s.second ); 339 } 340 } 341 } 342 } 343 344 const ValuedPoint upperLeft; 345 const ValuedPoint lowerLeft; 346 const ValuedPoint lowerRight; 347 const ValuedPoint upperRight; 348 const int nanCount; 349 const uint8_t borders; 350 const bool split; 351 352 private: 353 centerSquare354 ValuedPoint center() const 355 { 356 return ValuedPoint( 357 .5*(upperLeft.x + lowerRight.x), 358 .5*(upperLeft.y + lowerRight.y), 359 ( (std::isnan(lowerLeft.value) ? 0 : lowerLeft.value) 360 + (std::isnan(upperLeft.value) ? 0 : upperLeft.value) 361 + (std::isnan(lowerRight.value) ? 0 : lowerRight.value) 362 + (std::isnan(upperRight.value) ? 0 : upperRight.value) 363 ) / (4-nanCount)); 364 } 365 leftCenterSquare366 ValuedPoint leftCenter() const 367 { 368 return ValuedPoint( 369 upperLeft.x, 370 .5*(upperLeft.y + lowerLeft.y), 371 std::isnan(upperLeft.value) 372 ? lowerLeft.value 373 : (std::isnan(lowerLeft.value) ? upperLeft.value : .5*(upperLeft.value + lowerLeft.value))); 374 } 375 lowerCenterSquare376 ValuedPoint lowerCenter() const 377 { 378 return ValuedPoint( 379 .5*(lowerLeft.x + lowerRight.x), 380 lowerLeft.y, 381 std::isnan(lowerRight.value) 382 ? lowerLeft.value 383 : (std::isnan(lowerLeft.value) ? lowerRight.value : .5*(lowerRight.value + lowerLeft.value))); 384 } 385 rightCenterSquare386 ValuedPoint rightCenter() const 387 { 388 return ValuedPoint( 389 upperRight.x, 390 .5*(upperRight.y + lowerRight.y), 391 std::isnan(lowerRight.value) 392 ? upperRight.value 393 : (std::isnan(upperRight.value) ? lowerRight.value : .5*(lowerRight.value + upperRight.value))); 394 } 395 upperCenterSquare396 ValuedPoint upperCenter() const 397 { 398 return ValuedPoint( 399 .5*(upperLeft.x + upperRight.x), 400 upperLeft.y, 401 std::isnan(upperLeft.value) 402 ? upperRight.value 403 : (std::isnan(upperRight.value) ? upperLeft.value : .5*(upperLeft.value + upperRight.value))); 404 } 405 marchingCaseSquare406 uint8_t marchingCase(double level) const 407 { 408 return (level < fudge(level, upperLeft.value) ? UPPER_LEFT : ALL_LOW) 409 | (level < fudge(level, lowerLeft.value) ? LOWER_LEFT : ALL_LOW) 410 | (level < fudge(level, lowerRight.value) ? LOWER_RIGHT : ALL_LOW) 411 | (level < fudge(level, upperRight.value) ? UPPER_RIGHT : ALL_LOW); 412 413 } 414 interpolate_Square415 static double interpolate_(double level, double x1, double x2, double y1, double y2, bool need_split) 416 { 417 if (need_split) 418 { 419 // The two cases are here to avoid numerical roundup errors, for two points, we always compute 420 // the same interpolation. This condition is ensured by the order left->right bottom->top in interpole calls 421 // 422 // To obtain the same value for border (split) and non-border element, we take the 423 // middle value and interpolate from this to the end 424 const double xm = .5*(x1 + x2); 425 const double ym = .5*(y1 + y2); 426 const double fy1 = fudge( level, y1 ); 427 const double fym = fudge( level, ym ); 428 if ( (fy1 < level && level < fym) || (fy1 > level && level > fym) ) 429 { 430 x2 = xm; 431 y2 = ym; 432 } 433 else 434 { 435 x1 = xm; 436 y1 = ym; 437 } 438 } 439 const double fy1 = fudge(level, y1); 440 const double ratio = (level - fy1) / (fudge(level, y2) - fy1); 441 return x1 * (1. - ratio) + x2 * ratio; 442 } 443 interpolateSquare444 Point interpolate(uint8_t border, double level) const 445 { 446 switch (border) 447 { 448 case LEFT_BORDER: 449 return Point(upperLeft.x, interpolate_(level, lowerLeft.y, upperLeft.y, 450 lowerLeft.value, upperLeft.value, !split)); 451 case LOWER_BORDER: 452 return Point(interpolate_(level, lowerLeft.x, lowerRight.x, 453 lowerLeft.value, lowerRight.value, !split), lowerLeft.y); 454 case RIGHT_BORDER: 455 return Point(upperRight.x, interpolate_(level, lowerRight.y, upperRight.y, 456 lowerRight.value, upperRight.value, !split)); 457 case UPPER_BORDER: 458 return Point(interpolate_(level, upperLeft.x, upperRight.x, 459 upperLeft.value, upperRight.value, !split), upperLeft.y); 460 } 461 assert(false); 462 return Point(); 463 } 464 465 }; 466 } 467 #endif 468