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_LEVEL_GENERATOR_H
13 : #define MARCHING_SQUARE_LEVEL_GENERATOR_H
14 :
15 : #include <climits>
16 : #include <vector>
17 : #include <limits>
18 : #include <cmath>
19 : #include <cstdlib>
20 : #include "utility.h"
21 :
22 : namespace marching_squares
23 : {
24 :
25 : template <class Iterator> class Range
26 : {
27 : public:
28 830735 : Range(Iterator b, Iterator e) : begin_(b), end_(e)
29 : {
30 830735 : }
31 :
32 1635416 : Iterator begin() const
33 : {
34 1635416 : return begin_;
35 : }
36 :
37 831616 : Iterator end() const
38 : {
39 831616 : return end_;
40 : }
41 :
42 : private:
43 : Iterator begin_;
44 : Iterator end_;
45 : };
46 :
47 : template <typename LevelIterator> class RangeIterator
48 : {
49 : public:
50 1661261 : RangeIterator(const LevelIterator &parent, int idx)
51 1661261 : : parent_(parent), idx_(idx)
52 : {
53 1661261 : }
54 :
55 : // Warning: this is a "pseudo" iterator, since operator* returns a value,
56 : // not a reference. This means we cannot have operator->
57 224316 : std::pair<int, double> operator*() const
58 : {
59 224316 : return std::make_pair(idx_, parent_.level(idx_));
60 : }
61 :
62 922171 : bool operator!=(const RangeIterator &other) const
63 : {
64 922171 : return idx_ != other.idx_;
65 : }
66 :
67 986678 : const RangeIterator &operator++()
68 : {
69 986678 : idx_++;
70 986678 : return *this;
71 : }
72 :
73 : private:
74 : const LevelIterator &parent_;
75 : int idx_;
76 : };
77 :
78 : class FixedLevelRangeIterator
79 : {
80 : public:
81 : typedef RangeIterator<FixedLevelRangeIterator> Iterator;
82 :
83 74 : FixedLevelRangeIterator(const double *levels, size_t count, double minLevel,
84 : double maxLevel)
85 74 : : levels_(levels), count_(count), minLevel_(minLevel),
86 74 : maxLevel_(maxLevel)
87 : {
88 74 : }
89 :
90 830329 : Range<Iterator> range(double min, double max) const
91 : {
92 830329 : if (min > max)
93 1597 : std::swap(min, max);
94 830329 : size_t b = 0;
95 2205874 : for (; b != count_ && levels_[b] < fudge(min, minLevel_, levels_[b]);
96 : b++)
97 : ;
98 830329 : if (min == max)
99 : return Range<Iterator>(Iterator(*this, int(b)),
100 714742 : Iterator(*this, int(b)));
101 115587 : size_t e = b;
102 206232 : for (; e != count_ && levels_[e] <= fudge(max, minLevel_, levels_[e]);
103 : e++)
104 : ;
105 : return Range<Iterator>(Iterator(*this, int(b)),
106 115587 : Iterator(*this, int(e)));
107 : }
108 :
109 228560 : double level(int idx) const
110 : {
111 228560 : if (idx >= int(count_))
112 1456 : return maxLevel_;
113 227104 : return levels_[size_t(idx)];
114 : }
115 :
116 90644 : double minLevel() const
117 : {
118 90644 : return minLevel_;
119 : }
120 :
121 29 : size_t levelsCount() const
122 : {
123 29 : return count_;
124 : }
125 :
126 : private:
127 : const double *levels_;
128 : size_t count_;
129 : const double minLevel_;
130 : const double maxLevel_;
131 : };
132 :
133 : #if defined(__clang__)
134 : #pragma clang diagnostic push
135 : #pragma clang diagnostic ignored "-Wweak-vtables"
136 : #endif
137 :
138 : struct TooManyLevelsException : public std::exception
139 : {
140 4 : const char *what() const throw() override
141 : {
142 : return "Input values and/or interval settings would lead to too many "
143 4 : "levels";
144 : }
145 : };
146 :
147 : #if defined(__clang__)
148 : #pragma clang diagnostic pop
149 : #endif
150 :
151 : // Arbitrary threshold to avoid too much computation time and memory
152 : // consumption
153 : constexpr int knMAX_NUMBER_LEVELS = 100 * 1000;
154 :
155 : struct IntervalLevelRangeIterator
156 : {
157 : typedef RangeIterator<IntervalLevelRangeIterator> Iterator;
158 :
159 : // Construction by a offset and an interval
160 60 : IntervalLevelRangeIterator(double offset, double interval, double minLevel)
161 60 : : offset_(offset), interval_(interval), minLevel_(minLevel)
162 : {
163 60 : }
164 :
165 401 : Range<Iterator> range(double min, double max) const
166 : {
167 401 : if (min > max)
168 25 : std::swap(min, max);
169 :
170 : // compute the min index, adjusted to the fudged value if needed
171 401 : double df_i1 = ceil((min - offset_) / interval_);
172 401 : if (!(df_i1 >= INT_MIN && df_i1 < INT_MAX))
173 0 : throw TooManyLevelsException();
174 401 : int i1 = static_cast<int>(df_i1);
175 401 : double l1 = fudge(min, minLevel_, level(i1));
176 401 : if (l1 > min)
177 : {
178 28 : df_i1 = ceil((l1 - offset_) / interval_);
179 28 : if (!(df_i1 >= INT_MIN && df_i1 < INT_MAX))
180 0 : throw TooManyLevelsException();
181 28 : i1 = static_cast<int>(df_i1);
182 : }
183 401 : Iterator b(*this, i1);
184 :
185 401 : if (min == max)
186 220 : return Range<Iterator>(b, b);
187 :
188 : // compute the max index, adjusted to the fudged value if needed
189 181 : double df_i2 = floor((max - offset_) / interval_) + 1;
190 181 : if (!(df_i2 >= INT_MIN && df_i2 < INT_MAX))
191 1 : throw TooManyLevelsException();
192 180 : int i2 = static_cast<int>(df_i2);
193 180 : double l2 = fudge(max, minLevel_, level(i2));
194 180 : if (l2 > max)
195 : {
196 0 : df_i2 = floor((l2 - offset_) / interval_) + 1;
197 0 : if (!(df_i2 >= INT_MIN && df_i2 < INT_MAX))
198 0 : throw TooManyLevelsException();
199 0 : i2 = static_cast<int>(df_i2);
200 : }
201 180 : Iterator e(*this, i2);
202 :
203 : // Arbitrary threshold to avoid too much computation time and memory
204 : // consumption
205 180 : if (i2 > i1 + static_cast<double>(knMAX_NUMBER_LEVELS))
206 1 : throw TooManyLevelsException();
207 :
208 179 : return Range<Iterator>(b, e);
209 : }
210 :
211 1723 : double level(int idx) const
212 : {
213 1723 : return idx * interval_ + offset_;
214 : }
215 :
216 79 : double minLevel() const
217 : {
218 79 : return minLevel_;
219 : }
220 :
221 : private:
222 : const double offset_;
223 : const double interval_;
224 : const double minLevel_;
225 : };
226 :
227 : class ExponentialLevelRangeIterator
228 : {
229 : public:
230 : typedef RangeIterator<ExponentialLevelRangeIterator> Iterator;
231 :
232 9 : ExponentialLevelRangeIterator(double base, double minLevel)
233 9 : : base_(base), base_ln_(std::log(base_)), minLevel_(minLevel)
234 : {
235 9 : }
236 :
237 37 : double level(int idx) const
238 : {
239 37 : if (idx <= 0)
240 0 : return 0.0;
241 37 : return std::pow(base_, idx - 1);
242 : }
243 :
244 9 : Range<Iterator> range(double min, double max) const
245 : {
246 9 : if (min > max)
247 0 : std::swap(min, max);
248 :
249 9 : int i1 = index1(min);
250 9 : double l1 = fudge(min, minLevel_, level(i1));
251 9 : if (l1 > min)
252 0 : i1 = index1(l1);
253 9 : Iterator b(*this, i1);
254 :
255 9 : if (min == max)
256 0 : return Range<Iterator>(b, b);
257 :
258 9 : int i2 = index2(max);
259 9 : double l2 = fudge(max, minLevel_, level(i2));
260 9 : if (l2 > max)
261 0 : i2 = index2(l2);
262 9 : Iterator e(*this, i2);
263 :
264 : // Arbitrary threshold to avoid too much computation time and memory
265 : // consumption
266 9 : if (i2 > i1 + static_cast<double>(knMAX_NUMBER_LEVELS))
267 2 : throw TooManyLevelsException();
268 :
269 7 : return Range<Iterator>(b, e);
270 : }
271 :
272 : double minLevel() const
273 : {
274 : return minLevel_;
275 : }
276 :
277 : private:
278 9 : int index1(double plevel) const
279 : {
280 9 : if (plevel < 1.0)
281 3 : return 1;
282 6 : const double dfVal = ceil(std::log(plevel) / base_ln_) + 1;
283 6 : if (!(dfVal >= INT_MIN && dfVal < INT_MAX))
284 0 : throw TooManyLevelsException();
285 6 : return static_cast<int>(dfVal);
286 : }
287 :
288 9 : int index2(double plevel) const
289 : {
290 9 : if (plevel < 1.0)
291 0 : return 0;
292 9 : const double dfVal = floor(std::log(plevel) / base_ln_) + 1 + 1;
293 9 : if (!(dfVal >= INT_MIN && dfVal < INT_MAX))
294 0 : throw TooManyLevelsException();
295 9 : return static_cast<int>(dfVal);
296 : }
297 :
298 : // exponentiation base
299 : const double base_;
300 : const double base_ln_;
301 : const double minLevel_;
302 : };
303 :
304 : } // namespace marching_squares
305 : #endif
|
