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