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         629 :         Ring() : points(), interiorRings()
      36             :         {
      37         629 :         }
      38             : 
      39        3397 :         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       12760 :         bool isIn(const Ring &other) const
      49             :         {
      50             :             // Check if this is inside other using the winding number algorithm
      51       12760 :             auto checkPoint = this->points.front();
      52       12760 :             int windingNum = 0;
      53       12760 :             auto otherIter = other.points.begin();
      54             :             // p1 and p2 define each segment of the ring other that will be
      55             :             // tested
      56       12760 :             auto p1 = *otherIter;
      57      512354 :             while (true)
      58             :             {
      59      525114 :                 otherIter++;
      60      525114 :                 if (otherIter == other.points.end())
      61             :                 {
      62       12760 :                     break;
      63             :                 }
      64      512354 :                 auto p2 = *otherIter;
      65      512354 :                 if (p1.y <= checkPoint.y)
      66             :                 {
      67      276656 :                     if (p2.y > checkPoint.y)
      68             :                     {
      69        1069 :                         if (isLeft(p1, p2, checkPoint))
      70             :                         {
      71         417 :                             ++windingNum;
      72             :                         }
      73             :                     }
      74             :                 }
      75             :                 else
      76             :                 {
      77      235698 :                     if (p2.y <= checkPoint.y)
      78             :                     {
      79        1069 :                         if (!isLeft(p1, p2, checkPoint))
      80             :                         {
      81         624 :                             --windingNum;
      82             :                         }
      83             :                     }
      84             :                 }
      85      512354 :                 p1 = p2;
      86             :             }
      87       12760 :             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         735 :     void processTree(const std::vector<Ring> &tree, int level)
     108             :     {
     109         735 :         if (level % 2 == 0)
     110             :         {
     111         735 :             for (auto &r : tree)
     112             :             {
     113         348 :                 writer_.addPart(r.points);
     114         629 :                 for (auto &innerRing : r.interiorRings)
     115             :                 {
     116         281 :                     writer_.addInteriorRing(innerRing.points);
     117             :                 }
     118             :             }
     119             :         }
     120        1364 :         for (auto &r : tree)
     121             :         {
     122         629 :             processTree(r.interiorRings, level + 1);
     123             :         }
     124         735 :     }
     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          29 :     PolygonRingAppender(PolygonWriter &writer) : rings_(), writer_(writer)
     135             :     {
     136          29 :     }
     137             : 
     138         654 :     void addLine(double level, LineString &ls, bool)
     139             :     {
     140         654 :         auto &levelRings = rings_[level];
     141         654 :         if (ls.empty())
     142             :         {
     143          25 :             return;
     144             :         }
     145             :         // Create a new ring from the LineString
     146        1258 :         Ring newRing;
     147         629 :         newRing.points.swap(ls);
     148             :         // This queue holds the rings to be checked
     149        1258 :         std::deque<Ring *> queue;
     150         629 :         std::transform(levelRings.begin(), levelRings.end(),
     151        6384 :                        std::back_inserter(queue), [](Ring &r) { return &r; });
     152         629 :         Ring *parentRing = nullptr;
     153        7013 :         while (!queue.empty())
     154             :         {
     155        6384 :             Ring *curRing = queue.front();
     156        6384 :             queue.pop_front();
     157        6384 :             if (newRing.isIn(*curRing))
     158             :             {
     159             :                 // We know that there should only be one ring per level that we
     160             :                 // should fit in, so we can discard the rest of the queue and
     161             :                 // try again with the children of this ring
     162           8 :                 parentRing = curRing;
     163           8 :                 queue.clear();
     164           8 :                 std::transform(curRing->interiorRings.begin(),
     165             :                                curRing->interiorRings.end(),
     166             :                                std::back_inserter(queue),
     167           0 :                                [](Ring &r) { return &r; });
     168             :             }
     169             :         }
     170             :         // Get a pointer to the list we need to check for rings to include in
     171             :         // this ring
     172             :         std::vector<Ring> *parentRingList;
     173         629 :         if (parentRing == nullptr)
     174             :         {
     175         621 :             parentRingList = &levelRings;
     176             :         }
     177             :         else
     178             :         {
     179           8 :             parentRingList = &(parentRing->interiorRings);
     180             :         }
     181             :         // We found a valid parent, so we need to:
     182             :         // 1. Find all the inner rings of the parent that are inside the new
     183             :         // ring
     184         629 :         auto trueGroupIt = std::partition(
     185             :             parentRingList->begin(), parentRingList->end(),
     186        6376 :             [newRing](Ring &pRing) { return !pRing.isIn(newRing); });
     187             :         // 2. Move those rings out of the parent and into the new ring's
     188             :         // interior rings
     189         629 :         std::move(trueGroupIt, parentRingList->end(),
     190             :                   std::back_inserter(newRing.interiorRings));
     191             :         // 3. Get rid of the moved-from elements in the parent's interior rings
     192         629 :         parentRingList->erase(trueGroupIt, parentRingList->end());
     193             :         // 4. Add the new ring to the parent's interior rings
     194         629 :         parentRingList->push_back(newRing);
     195             :     }
     196             : 
     197          29 :     ~PolygonRingAppender()
     198             :     {
     199             :         // If there's no rings, nothing to do here
     200          29 :         if (rings_.size() == 0)
     201           0 :             return;
     202             : 
     203             :         // Traverse tree of rings
     204         135 :         for (auto &r : rings_)
     205             :         {
     206             :             // For each level, create a multipolygon by traversing the tree of
     207             :             // rings and adding a part for every other level
     208         106 :             writer_.startPolygon(r.first);
     209         106 :             processTree(r.second, 0);
     210         106 :             writer_.endPolygon();
     211             :         }
     212          29 :     }
     213             : };
     214             : 
     215             : }  // namespace marching_squares
     216             : 
     217             : #endif