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