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

Source (jump to first uncovered line)
/******************************************************************************
 *
 * 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, minLevel_, 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, minLevel_, 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, minLevel_, 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);
        double l2 = fudge(max, minLevel_, 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);
        double l1 = fudge(min, minLevel_, level(i1));
        if (l1 > min)
            i1 = index1(l1);
        Iterator b(*this, i1);

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

        int i2 = index2(max);
        double l2 = fudge(max, minLevel_, 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-06-13 06:29

Line	Count	Source (jump to first uncovered line)
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, minLevel_, levels_[b]);
96	0	b++)
97	0	;
98	0	if (min == max)
99	0	return Range<Iterator>(Iterator(*this, int(b)),
100	0	Iterator(*this, int(b)));
101	0	size_t e = b;
102	0	for (; e != count_ && levels_[e] <= fudge(max, minLevel_, levels_[e]);
103	0	e++)
104	0	;
105	0	return Range<Iterator>(Iterator(*this, int(b)),
106	0	Iterator(*this, int(e)));
107	0	}
108
109		double level(int idx) const
110	0	{
111	0	if (idx >= int(count_))
112	0	return maxLevel_;
113	0	return levels_[size_t(idx)];
114	0	}
115
116		double minLevel() const
117	0	{
118	0	return minLevel_;
119	0	}
120
121		size_t levelsCount() const
122	0	{
123	0	return count_;
124	0	}
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		const char *what() const throw() override
141	0	{
142	0	return "Input values and/or interval settings would lead to too many "
143	0	"levels";
144	0	}
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		IntervalLevelRangeIterator(double offset, double interval, double minLevel)
161	0	: offset_(offset), interval_(interval), minLevel_(minLevel)
162	0	{
163	0	}
164
165		Range<Iterator> range(double min, double max) const
166	0	{
167	0	if (min > max)
168	0	std::swap(min, max);
169
170		// compute the min index, adjusted to the fudged value if needed
171	0	double df_i1 = ceil((min - offset_) / interval_);
172	0	if (!(df_i1 >= INT_MIN && df_i1 < INT_MAX))
173	0	throw TooManyLevelsException();
174	0	int i1 = static_cast<int>(df_i1);
175	0	double l1 = fudge(min, minLevel_, level(i1));
176	0	if (l1 > min)
177	0	{
178	0	df_i1 = ceil((l1 - offset_) / interval_);
179	0	if (!(df_i1 >= INT_MIN && df_i1 < INT_MAX))
180	0	throw TooManyLevelsException();
181	0	i1 = static_cast<int>(df_i1);
182	0	}
183	0	Iterator b(*this, i1);
184
185	0	if (min == max)
186	0	return Range<Iterator>(b, b);
187
188		// compute the max index, adjusted to the fudged value if needed
189	0	double df_i2 = floor((max - offset_) / interval_) + 1;
190	0	if (!(df_i2 >= INT_MIN && df_i2 < INT_MAX))
191	0	throw TooManyLevelsException();
192	0	int i2 = static_cast<int>(df_i2);
193	0	double l2 = fudge(max, minLevel_, level(i2));
194	0	if (l2 > max)
195	0	{
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	0	}
201	0	Iterator e(*this, i2);
202
203		// Arbitrary threshold to avoid too much computation time and memory
204		// consumption
205	0	if (i2 > i1 + static_cast<double>(knMAX_NUMBER_LEVELS))
206	0	throw TooManyLevelsException();
207
208	0	return Range<Iterator>(b, e);
209	0	}
210
211		double level(int idx) const
212	0	{
213	0	return idx * interval_ + offset_;
214	0	}
215
216		double minLevel() const
217	0	{
218	0	return minLevel_;
219	0	}
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		ExponentialLevelRangeIterator(double base, double minLevel)
233	0	: base_(base), base_ln_(std::log(base_)), minLevel_(minLevel)
234	0	{
235	0	}
236
237		double level(int idx) const
238	0	{
239	0	if (idx <= 0)
240	0	return 0.0;
241	0	return std::pow(base_, idx - 1);
242	0	}
243
244		Range<Iterator> range(double min, double max) const
245	0	{
246	0	if (min > max)
247	0	std::swap(min, max);
248
249	0	int i1 = index1(min);
250	0	double l1 = fudge(min, minLevel_, level(i1));
251	0	if (l1 > min)
252	0	i1 = index1(l1);
253	0	Iterator b(*this, i1);
254
255	0	if (min == max)
256	0	return Range<Iterator>(b, b);
257
258	0	int i2 = index2(max);
259	0	double l2 = fudge(max, minLevel_, level(i2));
260	0	if (l2 > max)
261	0	i2 = index2(l2);
262	0	Iterator e(*this, i2);
263
264		// Arbitrary threshold to avoid too much computation time and memory
265		// consumption
266	0	if (i2 > i1 + static_cast<double>(knMAX_NUMBER_LEVELS))
267	0	throw TooManyLevelsException();
268
269	0	return Range<Iterator>(b, e);
270	0	}
271
272		double minLevel() const
273	0	{
274	0	return minLevel_;
275	0	}
276
277		private:
278		int index1(double plevel) const
279	0	{
280	0	if (plevel < 1.0)
281	0	return 1;
282	0	const double dfVal = ceil(std::log(plevel) / base_ln_) + 1;
283	0	if (!(dfVal >= INT_MIN && dfVal < INT_MAX))
284	0	throw TooManyLevelsException();
285	0	return static_cast<int>(dfVal);
286	0	}
287
288		int index2(double plevel) const
289	0	{
290	0	if (plevel < 1.0)
291	0	return 0;
292	0	const double dfVal = floor(std::log(plevel) / base_ln_) + 1 + 1;
293	0	if (!(dfVal >= INT_MIN && dfVal < INT_MAX))
294	0	throw TooManyLevelsException();
295	0	return static_cast<int>(dfVal);
296	0	}
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