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_SQUARES_SEGMENT_MERGER_H
13 : #define MARCHING_SQUARES_SEGMENT_MERGER_H
14 :
15 : #include "cpl_error.h"
16 : #include "point.h"
17 :
18 : #include <algorithm>
19 : #include <cassert>
20 : #include <list>
21 : #include <map>
22 : #include <vector>
23 :
24 : #include <iostream>
25 :
26 : namespace marching_squares
27 : {
28 :
29 : // SegmentMerger: join segments into linestrings and possibly into rings of
30 : // polygons
31 : template <typename LineWriter, typename LevelGenerator> struct SegmentMerger
32 : {
33 : struct LineStringEx
34 : {
35 : LineString ls = LineString();
36 : bool isMerged = false;
37 : };
38 :
39 : // a collection of unmerged linestrings
40 : typedef std::list<LineStringEx> Lines;
41 :
42 82 : SegmentMerger(LineWriter &lineWriter, const LevelGenerator &levelGenerator,
43 : bool polygonize_)
44 : : polygonize(polygonize_), lineWriter_(lineWriter), lines_(),
45 82 : endpointIndex_(), levelGenerator_(levelGenerator), m_anSkipLevels()
46 : {
47 82 : }
48 :
49 82 : ~SegmentMerger()
50 : {
51 82 : if (polygonize)
52 : {
53 189 : for (auto it = lines_.begin(); it != lines_.end(); ++it)
54 : {
55 148 : if (!it->second.empty())
56 0 : debug("remaining unclosed contour");
57 : }
58 : }
59 : // write all remaining (non-closed) lines
60 696 : for (auto it = lines_.begin(); it != lines_.end(); ++it)
61 : {
62 614 : const int levelIdx = it->first;
63 :
64 : // Skip levels that should be skipped
65 614 : if (std::find(m_anSkipLevels.begin(), m_anSkipLevels.end(),
66 1228 : levelIdx) != m_anSkipLevels.end())
67 : {
68 29 : continue;
69 : }
70 895 : while (it->second.begin() != it->second.end())
71 : {
72 310 : lineWriter_.addLine(levelGenerator_.level(levelIdx),
73 310 : it->second.begin()->ls, /* closed */ false);
74 310 : it->second.pop_front();
75 : }
76 : }
77 82 : }
78 :
79 476343 : void addSegment(int levelIdx, const Point &start, const Point &end)
80 : {
81 476343 : addSegment_(levelIdx, start, end);
82 476343 : }
83 :
84 95784 : void addBorderSegment(int levelIdx, const Point &start, const Point &end)
85 : {
86 95784 : addSegment_(levelIdx, start, end);
87 95784 : }
88 :
89 5886 : void beginningOfLine()
90 : {
91 5886 : if (polygonize)
92 3810 : return;
93 :
94 : // mark lines as non merged
95 23032 : for (auto &l : lines_)
96 : {
97 48921 : for (auto &ls : l.second)
98 : {
99 27965 : ls.isMerged = false;
100 : }
101 : }
102 : }
103 :
104 5886 : void endOfLine()
105 : {
106 5886 : if (polygonize)
107 3810 : return;
108 :
109 : // At the end of the line, we know that if no segment has been merged to
110 : // an existing line, it means there won't be anything more in the
111 : // future, we can then emit the line (this both speeds up and saves
112 : // memory)
113 :
114 23498 : for (auto &l : lines_)
115 : {
116 21422 : const int levelIdx = l.first;
117 21422 : auto it = l.second.begin();
118 54322 : while (it != l.second.end())
119 : {
120 32900 : if (!it->isMerged)
121 : {
122 : // Note that emitLine_ erases `it` and returns an iterator
123 : // advanced to the next element.
124 4625 : it = emitLine_(levelIdx, it, /* closed */ false);
125 : }
126 : else
127 : {
128 28275 : ++it;
129 : }
130 : }
131 : }
132 : }
133 :
134 : // non copyable
135 : SegmentMerger(const SegmentMerger<LineWriter, LevelGenerator> &) = delete;
136 : SegmentMerger<LineWriter, LevelGenerator> &
137 : operator=(const SegmentMerger<LineWriter, LevelGenerator> &) = delete;
138 :
139 : /**
140 : * @brief setSkipLevels sets the levels that should be skipped
141 : * when polygonize option is set.
142 : * @param anSkipLevels integer 0-based levels to skip.
143 : */
144 32 : void setSkipLevels(const std::vector<int> &anSkipLevels)
145 : {
146 : // Warn if polygonize is not set
147 32 : if (!polygonize)
148 : {
149 0 : CPLError(
150 : CE_Warning, CPLE_NotSupported,
151 : "setSkipLevels is ignored when polygonize option is not set");
152 : }
153 32 : m_anSkipLevels = anSkipLevels;
154 32 : }
155 :
156 : const bool polygonize;
157 :
158 : private:
159 : LineWriter &lineWriter_;
160 : // lines of each level
161 : std::map<int, Lines> lines_;
162 :
163 : // Index of the open lines of each level, keyed by their two endpoints.
164 : // Marching squares emits segments whose endpoints have degree two, so a
165 : // point can only ever be shared by two segment ends; a lookup here
166 : // replaces a linear scan over all open lines, which is quadratic on
167 : // rasters that keep many lines open at once (e.g. many long parallel
168 : // contours crossing each scanline).
169 : struct PointCompare
170 : {
171 36132400 : bool operator()(const Point &a, const Point &b) const
172 : {
173 36132400 : return a.x < b.x || (a.x == b.x && a.y < b.y);
174 : }
175 : };
176 :
177 : // A marching squares point has degree two, so at most two line ends can
178 : // sit on it (the two ends of the same line, when it has closed): a fixed
179 : // two-slot entry avoids a heap allocation per endpoint.
180 : struct EndpointEntry
181 : {
182 : typename Lines::iterator its[2];
183 : int n = 0;
184 : };
185 :
186 : typedef std::map<Point, EndpointEntry, PointCompare> EndpointIndex;
187 : std::map<int, EndpointIndex> endpointIndex_;
188 :
189 : const LevelGenerator &levelGenerator_;
190 :
191 : // Store 0-indexed levels to skip when polygonize option is set
192 : std::vector<int> m_anSkipLevels;
193 :
194 426385 : void registerEndpoint_(EndpointIndex &idx, typename Lines::iterator it,
195 : const Point &p)
196 : {
197 426385 : EndpointEntry &e = idx[p];
198 426385 : assert(e.n < 2);
199 426385 : e.its[e.n++] = it;
200 426385 : }
201 :
202 426385 : void unregisterEndpoint_(EndpointIndex &idx, typename Lines::iterator it,
203 : const Point &p)
204 : {
205 426385 : auto f = idx.find(p);
206 426385 : assert(f != idx.end());
207 426385 : EndpointEntry &e = f->second;
208 426385 : if (e.its[0] == it)
209 410007 : e.its[0] = e.its[1];
210 : else
211 16378 : assert(e.n == 2 && e.its[1] == it);
212 426385 : e.n--;
213 426385 : if (e.n == 0)
214 409880 : idx.erase(f);
215 426385 : }
216 :
217 572127 : void addSegment_(int levelIdx, const Point &start, const Point &end)
218 : {
219 :
220 572127 : Lines &lines = lines_[levelIdx];
221 :
222 572127 : if (start == end)
223 : {
224 0 : debug("degenerate segment (%f %f)", start.x, start.y);
225 0 : return;
226 : }
227 :
228 572127 : auto idxIt = endpointIndex_.find(levelIdx);
229 572127 : if (idxIt == endpointIndex_.end())
230 : {
231 178827 : addSegmentLinear_(levelIdx, lines, start, end);
232 : // The linear scans above are quadratic when many lines stay open
233 : // at once (e.g. long parallel contours crossing each scanline):
234 : // past this size, switch the level to the endpoint index. Below
235 : // it, the scans are cheaper than maintaining the index.
236 178827 : constexpr std::size_t INDEX_THRESHOLD = 100;
237 178827 : if (lines.size() > INDEX_THRESHOLD)
238 : {
239 2 : EndpointIndex &idx = endpointIndex_[levelIdx];
240 204 : for (auto it = lines.begin(); it != lines.end(); ++it)
241 : {
242 202 : registerEndpoint_(idx, it, it->ls.front());
243 202 : registerEndpoint_(idx, it, it->ls.back());
244 : }
245 : }
246 : }
247 : else
248 : {
249 393300 : addSegmentIndexed_(levelIdx, lines, idxIt->second, start, end);
250 : }
251 : }
252 :
253 : // Merge the segment by scanning the open lines: cheap while there are
254 : // few of them, quadratic when there are many.
255 178827 : void addSegmentLinear_(int levelIdx, Lines &lines, const Point &start,
256 : const Point &end)
257 : {
258 : // attempt to merge segment with existing line
259 178827 : auto it = lines.begin();
260 411037 : for (; it != lines.end(); ++it)
261 : {
262 393932 : if (it->ls.back() == end)
263 : {
264 2095 : it->ls.push_back(start);
265 2095 : it->isMerged = true;
266 2095 : break;
267 : }
268 391837 : if (it->ls.front() == end)
269 : {
270 91233 : it->ls.push_front(start);
271 91233 : it->isMerged = true;
272 91233 : break;
273 : }
274 300604 : if (it->ls.back() == start)
275 : {
276 64450 : it->ls.push_back(end);
277 64450 : it->isMerged = true;
278 64450 : break;
279 : }
280 236154 : if (it->ls.front() == start)
281 : {
282 3944 : it->ls.push_front(end);
283 3944 : it->isMerged = true;
284 3944 : break;
285 : }
286 : }
287 :
288 178827 : if (it == lines.end())
289 : {
290 : // new line
291 17105 : lines.push_back(LineStringEx());
292 17105 : lines.back().ls.push_back(start);
293 17105 : lines.back().ls.push_back(end);
294 17105 : lines.back().isMerged = true;
295 : }
296 161722 : else if (polygonize && (it->ls.front() == it->ls.back()))
297 : {
298 : // ring closed
299 760 : emitLine_(levelIdx, it, /* closed */ true);
300 760 : return;
301 : }
302 : else
303 : {
304 : // try to perform linemerge with another line
305 : // since we got out of the previous loop on the first match
306 : // there is no need to test previous elements
307 : // also: a segment merges at most two lines, no need to stall here
308 : // ;)
309 160962 : auto other = it;
310 160962 : ++other;
311 265325 : for (; other != lines.end(); ++other)
312 : {
313 115571 : if (it->ls.back() == other->ls.front())
314 : {
315 7739 : it->ls.pop_back();
316 7739 : it->ls.splice(it->ls.end(), other->ls);
317 7739 : it->isMerged = true;
318 7739 : lines.erase(other);
319 : // if that makes a closed ring, returns it
320 7739 : if (it->ls.front() == it->ls.back())
321 0 : emitLine_(levelIdx, it, /* closed */ true);
322 7739 : break;
323 : }
324 107832 : else if (other->ls.back() == it->ls.front())
325 : {
326 2461 : it->ls.pop_front();
327 2461 : other->ls.splice(other->ls.end(), it->ls);
328 2461 : other->isMerged = true;
329 2461 : lines.erase(it);
330 : // if that makes a closed ring, returns it
331 2461 : if (other->ls.front() == other->ls.back())
332 0 : emitLine_(levelIdx, other, /* closed */ true);
333 2461 : break;
334 : }
335 : // two lists must be merged but one is in the opposite direction
336 105371 : else if (it->ls.back() == other->ls.back())
337 : {
338 65 : it->ls.pop_back();
339 919 : for (auto rit = other->ls.rbegin(); rit != other->ls.rend();
340 854 : ++rit)
341 : {
342 854 : it->ls.push_back(*rit);
343 : }
344 65 : it->isMerged = true;
345 65 : lines.erase(other);
346 : // if that makes a closed ring, returns it
347 65 : if (it->ls.front() == it->ls.back())
348 0 : emitLine_(levelIdx, it, /* closed */ true);
349 65 : break;
350 : }
351 105306 : else if (it->ls.front() == other->ls.front())
352 : {
353 943 : it->ls.pop_front();
354 5437 : for (auto rit = other->ls.begin(); rit != other->ls.end();
355 4494 : ++rit)
356 : {
357 4494 : it->ls.push_front(*rit);
358 : }
359 943 : it->isMerged = true;
360 943 : lines.erase(other);
361 : // if that makes a closed ring, returns it
362 943 : if (it->ls.front() == it->ls.back())
363 0 : emitLine_(levelIdx, it, /* closed */ true);
364 943 : break;
365 : }
366 : }
367 : }
368 : }
369 :
370 : // Merge the segment through the endpoint index. Same merging rules as
371 : // the linear scan (a marching squares point has degree two, so at most
372 : // one line can match each of the segment's endpoints), but O(log n)
373 : // lookups instead of O(n) scans.
374 393300 : void addSegmentIndexed_(int levelIdx, Lines &lines, EndpointIndex &idx,
375 : const Point &start, const Point &end)
376 : {
377 : // attempt to merge the segment with an existing line whose endpoint
378 : // matches one of the segment's endpoints
379 1835300 : auto findLine = [&](const Point &p)
380 : {
381 606331 : auto f = idx.find(p);
382 1212660 : return f == idx.end() ? lines.end() : f->second.its[0];
383 : };
384 :
385 393300 : Point matched = end;
386 393300 : Point added = start;
387 393300 : typename Lines::iterator it = findLine(end);
388 393300 : if (it == lines.end())
389 : {
390 213031 : matched = start;
391 213031 : added = end;
392 213031 : it = findLine(start);
393 : }
394 :
395 393300 : if (it == lines.end())
396 : {
397 : // new line
398 16303 : lines.push_back(LineStringEx());
399 16303 : it = std::prev(lines.end());
400 16303 : it->ls.push_back(start);
401 16303 : it->ls.push_back(end);
402 16303 : it->isMerged = true;
403 16303 : registerEndpoint_(idx, it, start);
404 16303 : registerEndpoint_(idx, it, end);
405 16303 : return;
406 : }
407 :
408 : // extend the matched line with the segment's other endpoint
409 376997 : unregisterEndpoint_(idx, it, matched);
410 376997 : if (it->ls.back() == matched)
411 196728 : it->ls.push_back(added);
412 : else
413 180269 : it->ls.push_front(added);
414 376997 : it->isMerged = true;
415 376997 : registerEndpoint_(idx, it, added);
416 :
417 : // The extension may close a ring, or bring the line's new endpoint
418 : // onto another line's endpoint; merging two lines moves the free
419 : // endpoint again, so iterate until neither applies.
420 16378 : for (;;)
421 : {
422 393375 : if (polygonize && it->ls.front() == it->ls.back())
423 : {
424 : // ring closed
425 127 : emitLine_(levelIdx, it, /* closed */ true);
426 376997 : return;
427 : }
428 : // is there another line ending at `added`?
429 393248 : auto f = idx.find(added);
430 393248 : assert(f != idx.end());
431 393248 : typename Lines::iterator other = lines.end();
432 770118 : for (int i = 0; i < f->second.n; i++)
433 : {
434 393248 : if (f->second.its[i] != it)
435 : {
436 16378 : other = f->second.its[i];
437 16378 : break;
438 : }
439 : }
440 393248 : if (other == lines.end())
441 376870 : return;
442 :
443 : // merge `other` into `it` at the shared point `added`
444 16378 : unregisterEndpoint_(idx, it, added);
445 16378 : unregisterEndpoint_(idx, other, added);
446 16378 : const Point otherEnd = other->ls.front() == added
447 32756 : ? other->ls.back()
448 16378 : : other->ls.front();
449 16378 : unregisterEndpoint_(idx, other, otherEnd);
450 16378 : if (it->ls.back() == added)
451 : {
452 0 : it->ls.pop_back();
453 0 : if (other->ls.front() == added)
454 : {
455 0 : it->ls.splice(it->ls.end(), other->ls);
456 : }
457 : else
458 : {
459 : // opposite direction: append reversed
460 0 : for (auto rit = other->ls.rbegin(); rit != other->ls.rend();
461 0 : ++rit)
462 : {
463 0 : it->ls.push_back(*rit);
464 : }
465 : }
466 : }
467 : else
468 : {
469 16378 : it->ls.pop_front();
470 16378 : if (other->ls.back() == added)
471 : {
472 16378 : it->ls.splice(it->ls.begin(), other->ls);
473 : }
474 : else
475 : {
476 : // opposite direction: prepend reversed
477 0 : for (auto rit = other->ls.begin(); rit != other->ls.end();
478 0 : ++rit)
479 : {
480 0 : it->ls.push_front(*rit);
481 : }
482 : }
483 : }
484 16378 : lines.erase(other);
485 16378 : registerEndpoint_(idx, it, otherEnd);
486 16378 : added = otherEnd;
487 : }
488 : }
489 :
490 5512 : typename Lines::iterator emitLine_(int levelIdx,
491 : typename Lines::iterator it, bool closed)
492 : {
493 :
494 5512 : Lines &lines = lines_[levelIdx];
495 5512 : if (lines.empty())
496 0 : lines_.erase(levelIdx);
497 :
498 5512 : auto idxIt = endpointIndex_.find(levelIdx);
499 5512 : if (idxIt != endpointIndex_.end())
500 : {
501 127 : unregisterEndpoint_(idxIt->second, it, it->ls.front());
502 127 : unregisterEndpoint_(idxIt->second, it, it->ls.back());
503 : // An emptied index means no line is open at this level anymore:
504 : // drop it so the level returns to the linear path until it grows
505 : // past the threshold again.
506 127 : if (idxIt->second.empty())
507 2 : endpointIndex_.erase(idxIt);
508 : }
509 :
510 : // consume "it" and remove it from the list
511 : // but clear the line if the level should be skipped
512 5512 : if (std::find(m_anSkipLevels.begin(), m_anSkipLevels.end(), levelIdx) !=
513 11024 : m_anSkipLevels.end())
514 : {
515 32 : it->ls.clear();
516 : }
517 5512 : lineWriter_.addLine(levelGenerator_.level(levelIdx), it->ls, closed);
518 5512 : return lines.erase(it);
519 : }
520 : };
521 :
522 : } // namespace marching_squares
523 : #endif
|
