/src/gdal/alg/marching_squares/level_generator.h

Source
/******************************************************************************
 *
 * Project:  Marching square algorithm
 * Purpose:  Core algorithm implementation for contour line generation.
 * Author:   Oslandia <infos at oslandia dot com>
 *
 ******************************************************************************
 * Copyright (c) 2018, Oslandia <infos at oslandia dot com>
 *
 * SPDX-License-Identifier: MIT
 ****************************************************************************/
#ifndef MARCHING_SQUARE_LEVEL_GENERATOR_H
#define MARCHING_SQUARE_LEVEL_GENERATOR_H

#include <climits>
#include <vector>
#include <limits>
#include <cmath>
#include <cstdlib>
#include "utility.h"

namespace marching_squares
{

template <class Iterator> class Range
{
  public:
    Range(Iterator b, Iterator e) : begin_(b), end_(e)
    {
    }

    Iterator begin() const
    {
        return begin_;
    }

    Iterator end() const
    {
        return end_;
    }

  private:
    Iterator begin_;
    Iterator end_;
};

template <typename LevelIterator> class RangeIterator
{
  public:
    RangeIterator(const LevelIterator &parent, int idx)
        : parent_(parent), idx_(idx)
    {
    }

    // Warning: this is a "pseudo" iterator, since operator* returns a value,
    // not a reference. This means we cannot have operator->
    std::pair<int, double> operator*() const
    {
        return std::make_pair(idx_, parent_.level(idx_));
    }

    bool operator!=(const RangeIterator &other) const
    {
        return idx_ != other.idx_;
    }

    const RangeIterator &operator++()
    {
        idx_++;
        return *this;
    }

  private:
    const LevelIterator &parent_;
    int idx_;
};

class FixedLevelRangeIterator
{
  public:
    typedef RangeIterator<FixedLevelRangeIterator> Iterator;

    FixedLevelRangeIterator(const double *levels, size_t count, double minLevel,
                            double maxLevel)
        : levels_(levels), count_(count), minLevel_(minLevel),
          maxLevel_(maxLevel)
    {
    }

    Range<Iterator> range(double min, double max) const
    {
        if (min > max)
            std::swap(min, max);
        size_t b = 0;
        for (; b != count_ && levels_[b] < fudge(min, levels_[b]); b++)
            ;
        if (min == max)
            return Range<Iterator>(Iterator(*this, int(b)),
                                   Iterator(*this, int(b)));
        size_t e = b;
        for (; e != count_ && levels_[e] <= fudge(max, levels_[e]); e++)
            ;
        return Range<Iterator>(Iterator(*this, int(b)),
                               Iterator(*this, int(e)));
    }

    double level(int idx) const
    {
        if (idx >= int(count_))
            return maxLevel_;
        return levels_[size_t(idx)];
    }

    double minLevel() const
    {
        return minLevel_;
    }

    size_t levelsCount() const
    {
        return count_;
    }

  private:
    const double *levels_;
    size_t count_;
    const double minLevel_;
    const double maxLevel_;
};

#if defined(__clang__)
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wweak-vtables"
#endif

struct TooManyLevelsException : public std::exception
{
    const char *what() const throw() override
    {
        return "Input values and/or interval settings would lead to too many "
               "levels";
    }
};

#if defined(__clang__)
#pragma clang diagnostic pop
#endif

// Arbitrary threshold to avoid too much computation time and memory
// consumption
constexpr int knMAX_NUMBER_LEVELS = 100 * 1000;

struct IntervalLevelRangeIterator
{
    typedef RangeIterator<IntervalLevelRangeIterator> Iterator;

    // Construction by a offset and an interval
    IntervalLevelRangeIterator(double offset, double interval, double minLevel)
        : offset_(offset), interval_(interval), minLevel_(minLevel)
    {
    }

    Range<Iterator> range(double min, double max) const
    {
        if (min > max)
            std::swap(min, max);

        // compute the min index, adjusted to the fudged value if needed
        double df_i1 = ceil((min - offset_) / interval_);
        if (!(df_i1 >= INT_MIN && df_i1 < INT_MAX))
            throw TooManyLevelsException();
        int i1 = static_cast<int>(df_i1);
        double l1 = fudge(min, level(i1));
        if (l1 > min)
        {
            df_i1 = ceil((l1 - offset_) / interval_);
            if (!(df_i1 >= INT_MIN && df_i1 < INT_MAX))
                throw TooManyLevelsException();
            i1 = static_cast<int>(df_i1);
        }
        Iterator b(*this, i1);

        if (min == max)
            return Range<Iterator>(b, b);

        // compute the max index, adjusted to the fudged value if needed
        double df_i2 = floor((max - offset_) / interval_) + 1;
        if (!(df_i2 >= INT_MIN && df_i2 < INT_MAX))
            throw TooManyLevelsException();
        int i2 = static_cast<int>(df_i2);
        const double l2 = fudge(max, level(i2));
        if (l2 > max)
        {
            df_i2 = floor((l2 - offset_) / interval_) + 1;
            if (!(df_i2 >= INT_MIN && df_i2 < INT_MAX))
                throw TooManyLevelsException();
            i2 = static_cast<int>(df_i2);
        }
        Iterator e(*this, i2);

        // Arbitrary threshold to avoid too much computation time and memory
        // consumption
        if (i2 > i1 + static_cast<double>(knMAX_NUMBER_LEVELS))
            throw TooManyLevelsException();

        return Range<Iterator>(b, e);
    }

    double level(int idx) const
    {
        return idx * interval_ + offset_;
    }

    double minLevel() const
    {
        return minLevel_;
    }

  private:
    const double offset_;
    const double interval_;
    const double minLevel_;
};

class ExponentialLevelRangeIterator
{
  public:
    typedef RangeIterator<ExponentialLevelRangeIterator> Iterator;

    ExponentialLevelRangeIterator(double base, double minLevel)
        : base_(base), base_ln_(std::log(base_)), minLevel_(minLevel)
    {
    }

    double level(int idx) const
    {
        if (idx <= 0)
            return 0.0;
        return std::pow(base_, idx - 1);
    }

    Range<Iterator> range(double min, double max) const
    {
        if (min > max)
            std::swap(min, max);

        int i1 = index1(min);
        const double l1 = fudge(min, level(i1));
        if (l1 > min)
            i1 = index1(l1);
        Iterator b(*this, i1);

        if (min == max)
            return Range<Iterator>(b, b);

        int i2 = index2(max);
        const double l2 = fudge(max, level(i2));
        if (l2 > max)
            i2 = index2(l2);
        Iterator e(*this, i2);

        // Arbitrary threshold to avoid too much computation time and memory
        // consumption
        if (i2 > i1 + static_cast<double>(knMAX_NUMBER_LEVELS))
            throw TooManyLevelsException();

        return Range<Iterator>(b, e);
    }

    double minLevel() const
    {
        return minLevel_;
    }

  private:
    int index1(double plevel) const
    {
        if (plevel < 1.0)
            return 1;
        const double dfVal = ceil(std::log(plevel) / base_ln_) + 1;
        if (!(dfVal >= INT_MIN && dfVal < INT_MAX))
            throw TooManyLevelsException();
        return static_cast<int>(dfVal);
    }

    int index2(double plevel) const
    {
        if (plevel < 1.0)
            return 0;
        const double dfVal = floor(std::log(plevel) / base_ln_) + 1 + 1;
        if (!(dfVal >= INT_MIN && dfVal < INT_MAX))
            throw TooManyLevelsException();
        return static_cast<int>(dfVal);
    }

    // exponentiation base
    const double base_;
    const double base_ln_;
    const double minLevel_;
};

}  // namespace marching_squares
#endif

Coverage Report

Created: 2025-12-31 06:48

Line	Count	Source
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	0	Range(Iterator b, Iterator e) : begin_(b), end_(e)
29	0	{
30	0	} Unexecuted instantiation: marching_squares::Range<marching_squares::RangeIterator<marching_squares::ExponentialLevelRangeIterator> >::Range(marching_squares::RangeIterator<marching_squares::ExponentialLevelRangeIterator>, marching_squares::RangeIterator<marching_squares::ExponentialLevelRangeIterator>) Unexecuted instantiation: marching_squares::Range<marching_squares::RangeIterator<marching_squares::IntervalLevelRangeIterator> >::Range(marching_squares::RangeIterator<marching_squares::IntervalLevelRangeIterator>, marching_squares::RangeIterator<marching_squares::IntervalLevelRangeIterator>) Unexecuted instantiation: marching_squares::Range<marching_squares::RangeIterator<marching_squares::FixedLevelRangeIterator> >::Range(marching_squares::RangeIterator<marching_squares::FixedLevelRangeIterator>, marching_squares::RangeIterator<marching_squares::FixedLevelRangeIterator>)
31
32		Iterator begin() const
33	0	{
34	0	return begin_;
35	0	} Unexecuted instantiation: marching_squares::Range<marching_squares::RangeIterator<marching_squares::ExponentialLevelRangeIterator> >::begin() const Unexecuted instantiation: marching_squares::Range<marching_squares::RangeIterator<marching_squares::IntervalLevelRangeIterator> >::begin() const Unexecuted instantiation: marching_squares::Range<marching_squares::RangeIterator<marching_squares::FixedLevelRangeIterator> >::begin() const
36
37		Iterator end() const
38	0	{
39	0	return end_;
40	0	} Unexecuted instantiation: marching_squares::Range<marching_squares::RangeIterator<marching_squares::ExponentialLevelRangeIterator> >::end() const Unexecuted instantiation: marching_squares::Range<marching_squares::RangeIterator<marching_squares::IntervalLevelRangeIterator> >::end() const Unexecuted instantiation: marching_squares::Range<marching_squares::RangeIterator<marching_squares::FixedLevelRangeIterator> >::end() const
41
42		private:
43		Iterator begin_;
44		Iterator end_;
45		};
46
47		template <typename LevelIterator> class RangeIterator
48		{
49		public:
50		RangeIterator(const LevelIterator &parent, int idx)
51	0	: parent_(parent), idx_(idx)
52	0	{
53	0	} Unexecuted instantiation: marching_squares::RangeIterator<marching_squares::ExponentialLevelRangeIterator>::RangeIterator(marching_squares::ExponentialLevelRangeIterator const&, int) Unexecuted instantiation: marching_squares::RangeIterator<marching_squares::IntervalLevelRangeIterator>::RangeIterator(marching_squares::IntervalLevelRangeIterator const&, int) Unexecuted instantiation: marching_squares::RangeIterator<marching_squares::FixedLevelRangeIterator>::RangeIterator(marching_squares::FixedLevelRangeIterator const&, int)
54
55		// Warning: this is a "pseudo" iterator, since operator* returns a value,
56		// not a reference. This means we cannot have operator->
57		std::pair<int, double> operator*() const
58	0	{
59	0	return std::make_pair(idx_, parent_.level(idx_));
60	0	} Unexecuted instantiation: marching_squares::RangeIterator<marching_squares::ExponentialLevelRangeIterator>::operator() const Unexecuted instantiation: marching_squares::RangeIterator<marching_squares::IntervalLevelRangeIterator>::operator() const Unexecuted instantiation: marching_squares::RangeIterator<marching_squares::FixedLevelRangeIterator>::operator*() const
61
62		bool operator!=(const RangeIterator &other) const
63	0	{
64	0	return idx_ != other.idx_;
65	0	} Unexecuted instantiation: marching_squares::RangeIterator<marching_squares::ExponentialLevelRangeIterator>::operator!=(marching_squares::RangeIterator<marching_squares::ExponentialLevelRangeIterator> const&) const Unexecuted instantiation: marching_squares::RangeIterator<marching_squares::IntervalLevelRangeIterator>::operator!=(marching_squares::RangeIterator<marching_squares::IntervalLevelRangeIterator> const&) const Unexecuted instantiation: marching_squares::RangeIterator<marching_squares::FixedLevelRangeIterator>::operator!=(marching_squares::RangeIterator<marching_squares::FixedLevelRangeIterator> const&) const
66
67		const RangeIterator &operator++()
68	0	{
69	0	idx_++;
70	0	return *this;
71	0	} Unexecuted instantiation: marching_squares::RangeIterator<marching_squares::ExponentialLevelRangeIterator>::operator++() Unexecuted instantiation: marching_squares::RangeIterator<marching_squares::IntervalLevelRangeIterator>::operator++() Unexecuted instantiation: marching_squares::RangeIterator<marching_squares::FixedLevelRangeIterator>::operator++()
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		FixedLevelRangeIterator(const double *levels, size_t count, double minLevel,
84		double maxLevel)
85	0	: levels_(levels), count_(count), minLevel_(minLevel),
86	0	maxLevel_(maxLevel)
87	0	{
88	0	}
89
90		Range<Iterator> range(double min, double max) const
91	0	{
92	0	if (min > max)
93	0	std::swap(min, max);
94	0	size_t b = 0;
95	0	for (; b != count_ && levels_[b] < fudge(min, levels_[b]); b++)
96	0	;
97	0	if (min == max)
98	0	return Range<Iterator>(Iterator(*this, int(b)),
99	0	Iterator(*this, int(b)));
100	0	size_t e = b;
101	0	for (; e != count_ && levels_[e] <= fudge(max, levels_[e]); e++)
102	0	;
103	0	return Range<Iterator>(Iterator(*this, int(b)),
104	0	Iterator(*this, int(e)));
105	0	}
106
107		double level(int idx) const
108	0	{
109	0	if (idx >= int(count_))
110	0	return maxLevel_;
111	0	return levels_[size_t(idx)];
112	0	}
113
114		double minLevel() const
115	0	{
116	0	return minLevel_;
117	0	}
118
119		size_t levelsCount() const
120	0	{
121	0	return count_;
122	0	}
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		const char *what() const throw() override
139	0	{
140	0	return "Input values and/or interval settings would lead to too many "
141	0	"levels";
142	0	}
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		IntervalLevelRangeIterator(double offset, double interval, double minLevel)
159	0	: offset_(offset), interval_(interval), minLevel_(minLevel)
160	0	{
161	0	}
162
163		Range<Iterator> range(double min, double max) const
164	0	{
165	0	if (min > max)
166	0	std::swap(min, max);
167
168		// compute the min index, adjusted to the fudged value if needed
169	0	double df_i1 = ceil((min - offset_) / interval_);
170	0	if (!(df_i1 >= INT_MIN && df_i1 < INT_MAX))
171	0	throw TooManyLevelsException();
172	0	int i1 = static_cast<int>(df_i1);
173	0	double l1 = fudge(min, level(i1));
174	0	if (l1 > min)
175	0	{
176	0	df_i1 = ceil((l1 - offset_) / interval_);
177	0	if (!(df_i1 >= INT_MIN && df_i1 < INT_MAX))
178	0	throw TooManyLevelsException();
179	0	i1 = static_cast<int>(df_i1);
180	0	}
181	0	Iterator b(*this, i1);
182
183	0	if (min == max)
184	0	return Range<Iterator>(b, b);
185
186		// compute the max index, adjusted to the fudged value if needed
187	0	double df_i2 = floor((max - offset_) / interval_) + 1;
188	0	if (!(df_i2 >= INT_MIN && df_i2 < INT_MAX))
189	0	throw TooManyLevelsException();
190	0	int i2 = static_cast<int>(df_i2);
191	0	const double l2 = fudge(max, level(i2));
192	0	if (l2 > max)
193	0	{
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	0	}
199	0	Iterator e(*this, i2);
200
201		// Arbitrary threshold to avoid too much computation time and memory
202		// consumption
203	0	if (i2 > i1 + static_cast<double>(knMAX_NUMBER_LEVELS))
204	0	throw TooManyLevelsException();
205
206	0	return Range<Iterator>(b, e);
207	0	}
208
209		double level(int idx) const
210	0	{
211	0	return idx * interval_ + offset_;
212	0	}
213
214		double minLevel() const
215	0	{
216	0	return minLevel_;
217	0	}
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		ExponentialLevelRangeIterator(double base, double minLevel)
231	0	: base_(base), base_ln_(std::log(base_)), minLevel_(minLevel)
232	0	{
233	0	}
234
235		double level(int idx) const
236	0	{
237	0	if (idx <= 0)
238	0	return 0.0;
239	0	return std::pow(base_, idx - 1);
240	0	}
241
242		Range<Iterator> range(double min, double max) const
243	0	{
244	0	if (min > max)
245	0	std::swap(min, max);
246
247	0	int i1 = index1(min);
248	0	const double l1 = fudge(min, level(i1));
249	0	if (l1 > min)
250	0	i1 = index1(l1);
251	0	Iterator b(*this, i1);
252
253	0	if (min == max)
254	0	return Range<Iterator>(b, b);
255
256	0	int i2 = index2(max);
257	0	const double l2 = fudge(max, level(i2));
258	0	if (l2 > max)
259	0	i2 = index2(l2);
260	0	Iterator e(*this, i2);
261
262		// Arbitrary threshold to avoid too much computation time and memory
263		// consumption
264	0	if (i2 > i1 + static_cast<double>(knMAX_NUMBER_LEVELS))
265	0	throw TooManyLevelsException();
266
267	0	return Range<Iterator>(b, e);
268	0	}
269
270		double minLevel() const
271	0	{
272	0	return minLevel_;
273	0	}
274
275		private:
276		int index1(double plevel) const
277	0	{
278	0	if (plevel < 1.0)
279	0	return 1;
280	0	const double dfVal = ceil(std::log(plevel) / base_ln_) + 1;
281	0	if (!(dfVal >= INT_MIN && dfVal < INT_MAX))
282	0	throw TooManyLevelsException();
283	0	return static_cast<int>(dfVal);
284	0	}
285
286		int index2(double plevel) const
287	0	{
288	0	if (plevel < 1.0)
289	0	return 0;
290	0	const double dfVal = floor(std::log(plevel) / base_ln_) + 1 + 1;
291	0	if (!(dfVal >= INT_MIN && dfVal < INT_MAX))
292	0	throw TooManyLevelsException();
293	0	return static_cast<int>(dfVal);
294	0	}
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