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