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