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_POLYGON_RING_APPENDER_H
      13             : #define MARCHING_SQUARE_POLYGON_RING_APPENDER_H
      14             : 
      15             : #include <vector>
      16             : #include <list>
      17             : #include <map>
      18             : #include <deque>
      19             : #include <cassert>
      20             : #include <iterator>
      21             : 
      22             : #include "point.h"
      23             : #include "ogr_geometry.h"
      24             : 
      25             : namespace marching_squares
      26             : {
      27             : 
      28             : // Receive rings of different levels and organize them
      29             : // into multi-polygons with possible interior rings when requested.
      30             : template <typename PolygonWriter> class PolygonRingAppender
      31             : {
      32             :   private:
      33         144 :     struct Ring
      34             :     {
      35         610 :         Ring() : points(), interiorRings()
      36             :         {
      37         610 :         }
      38             : 
      39        3338 :         Ring(const Ring &other) = default;
      40             :         Ring &operator=(const Ring &other) = default;
      41             : 
      42             :         LineString points;
      43             : 
      44             :         mutable std::vector<Ring> interiorRings;
      45             : 
      46             :         const Ring *closestExterior = nullptr;
      47             : 
      48       12756 :         bool isIn(const Ring &other) const
      49             :         {
      50             :             // Check if this is inside other using the winding number algorithm
      51       12756 :             auto checkPoint = this->points.front();
      52       12756 :             int windingNum = 0;
      53       12756 :             auto otherIter = other.points.begin();
      54             :             // p1 and p2 define each segment of the ring other that will be
      55             :             // tested
      56       12756 :             auto p1 = *otherIter;
      57      506898 :             while (true)
      58             :             {
      59      519654 :                 otherIter++;
      60      519654 :                 if (otherIter == other.points.end())
      61             :                 {
      62       12756 :                     break;
      63             :                 }
      64      506898 :                 auto p2 = *otherIter;
      65      506898 :                 if (p1.y <= checkPoint.y)
      66             :                 {
      67      273540 :                     if (p2.y > checkPoint.y)
      68             :                     {
      69        1068 :                         if (isLeft(p1, p2, checkPoint))
      70             :                         {
      71         416 :                             ++windingNum;
      72             :                         }
      73             :                     }
      74             :                 }
      75             :                 else
      76             :                 {
      77      233358 :                     if (p2.y <= checkPoint.y)
      78             :                     {
      79        1068 :                         if (!isLeft(p1, p2, checkPoint))
      80             :                         {
      81         620 :                             --windingNum;
      82             :                         }
      83             :                     }
      84             :                 }
      85      506898 :                 p1 = p2;
      86             :             }
      87       12756 :             return windingNum != 0;
      88             :         }
      89             : 
      90             : #ifdef DEBUG
      91             :         size_t id() const
      92             :         {
      93             :             return size_t(static_cast<const void *>(this)) & 0xffff;
      94             :         }
      95             : 
      96             :         void print(std::ostream &ostr) const
      97             :         {
      98             :             ostr << id() << ":";
      99             :             for (const auto &pt : points)
     100             :             {
     101             :                 ostr << pt.x << "," << pt.y << " ";
     102             :             }
     103             :         }
     104             : #endif
     105             :     };
     106             : 
     107         688 :     void processTree(const std::vector<Ring> &tree, int level)
     108             :     {
     109         688 :         if (level % 2 == 0)
     110             :         {
     111         688 :             for (auto &r : tree)
     112             :             {
     113         334 :                 writer_.addPart(r.points);
     114         610 :                 for (auto &innerRing : r.interiorRings)
     115             :                 {
     116         276 :                     writer_.addInteriorRing(innerRing.points);
     117             :                 }
     118             :             }
     119             :         }
     120        1298 :         for (auto &r : tree)
     121             :         {
     122         610 :             processTree(r.interiorRings, level + 1);
     123             :         }
     124         688 :     }
     125             : 
     126             :     // level -> rings
     127             :     std::map<double, std::vector<Ring>> rings_;
     128             : 
     129             :     PolygonWriter &writer_;
     130             : 
     131             :   public:
     132             :     const bool polygonize = true;
     133             : 
     134          19 :     PolygonRingAppender(PolygonWriter &writer) : rings_(), writer_(writer)
     135             :     {
     136          19 :     }
     137             : 
     138         610 :     void addLine(double level, LineString &ls, bool)
     139             :     {
     140             :         // Create a new ring from the LineString
     141        1220 :         Ring newRing;
     142         610 :         newRing.points.swap(ls);
     143         610 :         auto &levelRings = rings_[level];
     144             :         // This queue holds the rings to be checked
     145        1220 :         std::deque<Ring *> queue;
     146         610 :         std::transform(levelRings.begin(), levelRings.end(),
     147        6381 :                        std::back_inserter(queue), [](Ring &r) { return &r; });
     148         610 :         Ring *parentRing = nullptr;
     149        6991 :         while (!queue.empty())
     150             :         {
     151        6381 :             Ring *curRing = queue.front();
     152        6381 :             queue.pop_front();
     153        6381 :             if (newRing.isIn(*curRing))
     154             :             {
     155             :                 // We know that there should only be one ring per level that we
     156             :                 // should fit in, so we can discard the rest of the queue and
     157             :                 // try again with the children of this ring
     158           6 :                 parentRing = curRing;
     159           6 :                 queue.clear();
     160           6 :                 std::transform(curRing->interiorRings.begin(),
     161             :                                curRing->interiorRings.end(),
     162             :                                std::back_inserter(queue),
     163           0 :                                [](Ring &r) { return &r; });
     164             :             }
     165             :         }
     166             :         // Get a pointer to the list we need to check for rings to include in
     167             :         // this ring
     168             :         std::vector<Ring> *parentRingList;
     169         610 :         if (parentRing == nullptr)
     170             :         {
     171         604 :             parentRingList = &levelRings;
     172             :         }
     173             :         else
     174             :         {
     175           6 :             parentRingList = &(parentRing->interiorRings);
     176             :         }
     177             :         // We found a valid parent, so we need to:
     178             :         // 1. Find all the inner rings of the parent that are inside the new
     179             :         // ring
     180         610 :         auto trueGroupIt = std::partition(
     181             :             parentRingList->begin(), parentRingList->end(),
     182        6375 :             [newRing](Ring &pRing) { return !pRing.isIn(newRing); });
     183             :         // 2. Move those rings out of the parent and into the new ring's
     184             :         // interior rings
     185         610 :         std::move(trueGroupIt, parentRingList->end(),
     186             :                   std::back_inserter(newRing.interiorRings));
     187             :         // 3. Get rid of the moved-from elements in the parent's interior rings
     188         610 :         parentRingList->erase(trueGroupIt, parentRingList->end());
     189             :         // 4. Add the new ring to the parent's interior rings
     190         610 :         parentRingList->push_back(newRing);
     191         610 :     }
     192             : 
     193          19 :     ~PolygonRingAppender()
     194             :     {
     195             :         // If there's no rings, nothing to do here
     196          19 :         if (rings_.size() == 0)
     197           0 :             return;
     198             : 
     199             :         // Traverse tree of rings
     200          97 :         for (auto &r : rings_)
     201             :         {
     202             :             // For each level, create a multipolygon by traversing the tree of
     203             :             // rings and adding a part for every other level
     204          78 :             writer_.startPolygon(r.first);
     205          78 :             processTree(r.second, 0);
     206          78 :             writer_.endPolygon();
     207             :         }
     208          19 :     }
     209             : };
     210             : 
     211             : }  // namespace marching_squares
     212             : 
     213             : #endif