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 : #include <memory>
22 : #include <algorithm>
23 :
24 : #include "cpl_quad_tree.h"
25 :
26 : #include "point.h"
27 : #include "ogr_api.h"
28 : #include "ogr_geometry.h"
29 :
30 : namespace marching_squares
31 : {
32 :
33 : // Receive rings of different levels and organize them
34 : // into multi-polygons with possible interior rings when requested.
35 : template <typename PolygonWriter> class PolygonRingAppender
36 : {
37 : private:
38 0 : struct Ring
39 : {
40 847 : Ring() : points(), bbox(), interiorRings()
41 : {
42 847 : }
43 :
44 1753 : Ring(const Ring &other) = default;
45 : Ring &operator=(const Ring &other) = default;
46 : // Declaring the copy operations above suppresses the
47 : // implicit move operations, so vector reshuffles and reallocations
48 : // deep-copied entire ring subtrees. Restore them.
49 1686 : Ring(Ring &&other) = default;
50 : Ring &operator=(Ring &&other) = default;
51 :
52 : LineString points;
53 :
54 : // Bounding box, computed once when the ring is complete;
55 : // gives isIn() an O(1) reject so parent search stops being
56 : // O(rings * vertices) per insertion.
57 : OGREnvelope bbox;
58 :
59 847 : void computeBBox()
60 : {
61 847 : bbox = OGREnvelope();
62 489934 : for (const auto &pt : points)
63 489087 : bbox.Merge(pt.x, pt.y);
64 847 : }
65 :
66 : mutable std::vector<Ring> interiorRings;
67 :
68 : const Ring *closestExterior = nullptr;
69 :
70 411 : bool isIn(const Ring &other) const
71 : {
72 : // Check if this is inside other using the winding number algorithm
73 411 : auto checkPoint = this->points.front();
74 : // A point outside the candidate ring's bounding box
75 : // cannot be inside the ring.
76 411 : if (checkPoint.x < other.bbox.MinX ||
77 411 : checkPoint.x > other.bbox.MaxX ||
78 411 : checkPoint.y < other.bbox.MinY ||
79 411 : checkPoint.y > other.bbox.MaxY)
80 : {
81 0 : return false;
82 : }
83 411 : int windingNum = 0;
84 411 : auto otherIter = other.points.begin();
85 : // p1 and p2 define each segment of the ring other that will be
86 : // tested
87 411 : auto p1 = *otherIter;
88 6138680 : while (true)
89 : {
90 6139093 : otherIter++;
91 6139093 : if (otherIter == other.points.end())
92 : {
93 411 : break;
94 : }
95 6138680 : auto p2 = *otherIter;
96 6138680 : if (p1.y <= checkPoint.y)
97 : {
98 5396957 : if (p2.y > checkPoint.y)
99 : {
100 623 : if (isLeft(p1, p2, checkPoint))
101 : {
102 253 : ++windingNum;
103 : }
104 : }
105 : }
106 : else
107 : {
108 741725 : if (p2.y <= checkPoint.y)
109 : {
110 623 : if (!isLeft(p1, p2, checkPoint))
111 : {
112 434 : --windingNum;
113 : }
114 : }
115 : }
116 6138680 : p1 = p2;
117 : }
118 411 : return windingNum != 0;
119 : }
120 :
121 : #ifdef DEBUG
122 : size_t id() const
123 : {
124 : return size_t(static_cast<const void *>(this)) & 0xffff;
125 : }
126 :
127 : void print(std::ostream &ostr) const
128 : {
129 : ostr << id() << ":";
130 : for (const auto &pt : points)
131 : {
132 : ostr << pt.x << "," << pt.y << " ";
133 : }
134 : }
135 : #endif
136 : };
137 :
138 976 : void processTree(const std::vector<Ring> &tree, int level)
139 : {
140 976 : if (level % 2 == 0)
141 : {
142 976 : for (auto &r : tree)
143 : {
144 466 : writer_.addPart(r.points);
145 847 : for (auto &innerRing : r.interiorRings)
146 : {
147 381 : writer_.addInteriorRing(innerRing.points);
148 : }
149 : }
150 : }
151 1823 : for (auto &r : tree)
152 : {
153 847 : processTree(r.interiorRings, level + 1);
154 : }
155 976 : }
156 :
157 : // level -> rings
158 : std::map<double, std::vector<Ring>> rings_;
159 :
160 : // Point-in-polygon accelerator for one target ring: an
161 : // OGRPreparedGeometry (GEOS indexed point-in-area locator) over the
162 : // ring, built lazily when a ring turns out to capture many candidates.
163 : // The pathological case is a domain-spanning ring with millions of
164 : // vertices capturing tens of thousands of earlier rings; testing each
165 : // candidate against the raw ring is O(candidates * vertices). In builds
166 : // without GEOS support the capture step falls back to the linear
167 : // winding test.
168 : struct PreparedRing
169 : {
170 839 : PreparedRing() : poly(), prep()
171 : {
172 839 : }
173 :
174 : OGRPolygon poly;
175 : OGRPreparedGeometryUniquePtr prep;
176 :
177 12 : bool build(const Ring &r)
178 : {
179 12 : poly.empty();
180 12 : prep.reset();
181 24 : auto ring = std::make_unique<OGRLinearRing>();
182 12 : ring->setNumPoints(static_cast<int>(r.points.size()));
183 12 : int i = 0;
184 371468 : for (const auto &pt : r.points)
185 371456 : ring->setPoint(i++, pt.x, pt.y);
186 12 : poly.addRingDirectly(ring.release());
187 12 : poly.closeRings();
188 12 : prep.reset(OGRCreatePreparedGeometry(OGRGeometry::ToHandle(&poly)));
189 24 : return prep != nullptr;
190 : }
191 :
192 252 : bool contains(const Point &p) const
193 : {
194 504 : OGRPoint pt(p.x, p.y);
195 252 : return CPL_TO_BOOL(OGRPreparedGeometryContains(
196 504 : prep.get(), OGRGeometry::ToHandle(&pt)));
197 : }
198 : };
199 :
200 : // Per-level spatial index over TOP-LEVEL rings: a CPLQuadTree over ring
201 : // bounding boxes. Each stored feature points at the ring's slot index in
202 : // the level's ring vector, held in a std::deque so the pointer survives
203 : // growth; vector reallocation of the rings themselves is harmless. Rings
204 : // captured as interior rings of a later ring are removed from the tree
205 : // and their slot tombstoned (points cleared) rather than erased, keeping
206 : // the remaining indices stable.
207 : struct QuadTreeDestroyer
208 : {
209 129 : void operator()(CPLQuadTree *t) const
210 : {
211 129 : CPLQuadTreeDestroy(t);
212 129 : }
213 : };
214 :
215 : std::map<double, std::unique_ptr<CPLQuadTree, QuadTreeDestroyer>> index_;
216 : std::map<double, std::deque<std::size_t>> slots_;
217 : CPLRectObj domain_;
218 :
219 699 : static std::size_t featureSlot(const void *f)
220 : {
221 699 : return *static_cast<const std::size_t *>(f);
222 : }
223 :
224 2149 : static CPLRectObj ringRect(const Ring &r)
225 : {
226 2149 : return CPLRectObj{r.bbox.MinX, r.bbox.MinY, r.bbox.MaxX, r.bbox.MaxY};
227 : }
228 :
229 : PolygonWriter &writer_;
230 :
231 : public:
232 : const bool polygonize = true;
233 :
234 37 : PolygonRingAppender(PolygonWriter &writer, double minX, double minY,
235 : double maxX, double maxY)
236 : : rings_(), index_(), slots_(), domain_{minX, minY, maxX, maxY},
237 37 : writer_(writer)
238 : {
239 37 : }
240 :
241 879 : void addLine(double level, LineString &ls, bool)
242 : {
243 879 : auto &levelRings = rings_[level];
244 879 : auto &levelTree = index_[level];
245 879 : auto &levelSlots = slots_[level];
246 879 : if (!levelTree)
247 129 : levelTree.reset(CPLQuadTreeCreate(&domain_, nullptr));
248 879 : if (ls.empty())
249 : {
250 32 : return;
251 : }
252 : // Create a new ring from the LineString
253 1694 : Ring newRing;
254 847 : newRing.points.swap(ls);
255 847 : newRing.computeBBox();
256 : // Find the top-level parent (if any) through the index instead of
257 : // scanning every top-level ring, then descend the (short) nested
258 : // sibling lists exactly as before.
259 847 : Ring *parentRing = nullptr;
260 : {
261 847 : Ring *top = nullptr;
262 847 : const auto &fp0 = newRing.points.front();
263 847 : CPLRectObj aoi{fp0.x, fp0.y, fp0.x, fp0.y};
264 847 : int nHits = 0;
265 847 : void **hits = CPLQuadTreeSearch(levelTree.get(), &aoi, &nHits);
266 867 : for (int h = 0; h < nHits && top == nullptr; h++)
267 : {
268 20 : Ring &cand = levelRings[featureSlot(hits[h])];
269 20 : if (!cand.points.empty() && newRing.isIn(cand))
270 8 : top = &cand;
271 : }
272 847 : CPLFree(hits);
273 847 : if (top != nullptr)
274 : {
275 8 : parentRing = top;
276 : // This queue holds the rings to be checked
277 16 : std::deque<Ring *> queue;
278 8 : std::transform(
279 : top->interiorRings.begin(), top->interiorRings.end(),
280 0 : std::back_inserter(queue), [](Ring &r) { return &r; });
281 8 : while (!queue.empty())
282 : {
283 0 : Ring *curRing = queue.front();
284 0 : queue.pop_front();
285 0 : if (newRing.isIn(*curRing))
286 : {
287 : // We know that there should only be one ring per
288 : // level that we should fit in, so we can discard the
289 : // rest of the queue and try again with the children
290 : // of this ring
291 0 : parentRing = curRing;
292 0 : queue.clear();
293 0 : std::transform(curRing->interiorRings.begin(),
294 : curRing->interiorRings.end(),
295 : std::back_inserter(queue),
296 0 : [](Ring &r) { return &r; });
297 : }
298 : }
299 : }
300 : }
301 847 : if (parentRing == nullptr)
302 : {
303 : // Top-level insertion: capture existing top-level rings that lie
304 : // inside the new ring, via the index. Build a per-target PIP
305 : // index lazily so a huge ring capturing many candidates costs
306 : // O(V + R * V/B), not O(R * V).
307 1678 : std::vector<std::size_t> captured;
308 1678 : PreparedRing pip;
309 839 : bool pipTried = false;
310 839 : bool pipBuilt = false;
311 839 : std::size_t nCandidates = 0;
312 : {
313 839 : CPLRectObj aoi = ringRect(newRing);
314 839 : int nHits = 0;
315 839 : void **hits = CPLQuadTreeSearch(levelTree.get(), &aoi, &nHits);
316 1518 : for (int h = 0; h < nHits; h++)
317 : {
318 679 : const std::size_t idx = featureSlot(hits[h]);
319 679 : Ring &cand = levelRings[idx];
320 679 : if (cand.points.empty())
321 36 : continue;
322 679 : const auto &fp = cand.points.front();
323 679 : if (fp.x < newRing.bbox.MinX || fp.x > newRing.bbox.MaxX ||
324 643 : fp.y < newRing.bbox.MinY || fp.y > newRing.bbox.MaxY)
325 36 : continue;
326 655 : if (!pipTried && ++nCandidates > 16 &&
327 12 : newRing.points.size() > 512)
328 : {
329 12 : pipTried = true;
330 12 : pipBuilt = pip.build(newRing);
331 : }
332 643 : const bool inside =
333 643 : pipBuilt ? pip.contains(fp) : cand.isIn(newRing);
334 643 : if (inside)
335 471 : captured.push_back(idx);
336 : }
337 839 : CPLFree(hits);
338 : }
339 : // Sorting by slot restores insertion order, so captured rings
340 : // nest in the same order the original linear scan produced.
341 839 : std::sort(captured.begin(), captured.end());
342 839 : captured.erase(std::unique(captured.begin(), captured.end()),
343 839 : captured.end());
344 1310 : for (std::size_t idx : captured)
345 : {
346 471 : CPLRectObj rb = ringRect(levelRings[idx]);
347 471 : CPLQuadTreeRemove(levelTree.get(), &levelSlots[idx], &rb);
348 471 : newRing.interiorRings.push_back(std::move(levelRings[idx]));
349 471 : levelRings[idx].points.clear(); // tombstone the slot
350 : }
351 839 : levelRings.push_back(std::move(newRing));
352 839 : levelSlots.push_back(levelRings.size() - 1);
353 839 : CPLRectObj nb = ringRect(levelRings.back());
354 839 : CPLQuadTreeInsertWithBounds(levelTree.get(), &levelSlots.back(),
355 : &nb);
356 : }
357 : else
358 : {
359 : // Get a pointer to the list we need to check for rings to include
360 : // in this ring
361 8 : std::vector<Ring> *parentRingList = &(parentRing->interiorRings);
362 : // We found a valid parent, so we need to:
363 : // 1. Find all the inner rings of the parent that are inside the new
364 : // ring
365 8 : auto trueGroupIt = std::partition(
366 : parentRingList->begin(), parentRingList->end(),
367 0 : [&newRing](Ring &pRing) { return !pRing.isIn(newRing); });
368 : // 2. Move those rings out of the parent and into the new ring's
369 : // interior rings
370 8 : std::move(trueGroupIt, parentRingList->end(),
371 : std::back_inserter(newRing.interiorRings));
372 : // 3. Get rid of the moved-from elements in the parent's interior
373 : // rings
374 8 : parentRingList->erase(trueGroupIt, parentRingList->end());
375 : // 4. Add the new ring to the parent's interior rings
376 8 : parentRingList->push_back(std::move(newRing));
377 : }
378 : }
379 :
380 37 : ~PolygonRingAppender()
381 : {
382 : // If there's no rings, nothing to do here
383 37 : if (rings_.size() == 0)
384 0 : return;
385 :
386 : // Traverse tree of rings
387 166 : for (auto &r : rings_)
388 : {
389 : // Drop tombstoned slots (rings captured as interior
390 : // rings of later-arriving parents) before traversal.
391 258 : std::vector<Ring> live;
392 129 : live.reserve(r.second.size());
393 968 : for (auto &ring : r.second)
394 839 : if (!ring.points.empty())
395 368 : live.push_back(std::move(ring));
396 : // For each level, create a multipolygon by traversing the tree of
397 : // rings and adding a part for every other level
398 129 : writer_.startPolygon(r.first);
399 129 : processTree(live, 0);
400 129 : writer_.endPolygon();
401 : }
402 37 : }
403 : };
404 :
405 : } // namespace marching_squares
406 :
407 : #endif
|
