/src/gdal/alg/marching_squares/square.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_SQUARE_H
#define MARCHING_SQUARE_SQUARE_H

#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <math.h>
#include "utility.h"
#include "point.h"

namespace marching_squares
{

struct Square
{
    // Bit flags to determine borders around pixel
    static const uint8_t NO_BORDER = 0;          // 0000 0000
    static const uint8_t LEFT_BORDER = 1 << 0;   // 0000 0001
    static const uint8_t LOWER_BORDER = 1 << 1;  // 0000 0010
    static const uint8_t RIGHT_BORDER = 1 << 2;  // 0000 0100
    static const uint8_t UPPER_BORDER = 1 << 3;  // 0000 1000

    // Bit flags for marching square case
    static const uint8_t ALL_LOW = 0;           // 0000 0000
    static const uint8_t UPPER_LEFT = 1 << 0;   // 0000 0001
    static const uint8_t LOWER_LEFT = 1 << 1;   // 0000 0010
    static const uint8_t LOWER_RIGHT = 1 << 2;  // 0000 0100
    static const uint8_t UPPER_RIGHT = 1 << 3;  // 0000 1000
    static const uint8_t ALL_HIGH =
        UPPER_LEFT | LOWER_LEFT | LOWER_RIGHT | UPPER_RIGHT;    // 0000 1111
    static const uint8_t SADDLE_NW = UPPER_LEFT | LOWER_RIGHT;  // 0000 0101
    static const uint8_t SADDLE_NE = UPPER_RIGHT | LOWER_LEFT;  // 0000 1010

    typedef std::pair<Point, Point> Segment;
    typedef std::pair<ValuedPoint, ValuedPoint> ValuedSegment;

    //
    // An array of segments, at most 3 segments
    struct Segments
    {
        Segments() : sz_(0)
        {
        }

        Segments(const Segment &first) : sz_(1), segs_()
        {
            segs_[0] = first;
        }

        Segments(const Segment &first, const Segment &second) : sz_(2), segs_()
        {
            segs_[0] = first;
            segs_[1] = second;
        }

        Segments(const Segment &first, const Segment &second,
                 const Segment &third)
            : sz_(3), segs_()
        {
            segs_[0] = first;
            segs_[1] = second;
            segs_[2] = third;
        }

        std::size_t size() const
        {
            return sz_;
        }

        const Segment &operator[](std::size_t idx) const
        {
            assert(idx < sz_);
            return segs_[idx];
        }

      private:
        const std::size_t sz_;
        /* const */ Segment segs_[3];
    };

    Square(const ValuedPoint &upperLeft_, const ValuedPoint &upperRight_,
           const ValuedPoint &lowerLeft_, const ValuedPoint &lowerRight_,
           uint8_t borders_ = NO_BORDER, bool split_ = false)
        : upperLeft(upperLeft_), lowerLeft(lowerLeft_), lowerRight(lowerRight_),
          upperRight(upperRight_),
          nanCount((std::isnan(upperLeft.value) ? 1 : 0) +
                   (std::isnan(upperRight.value) ? 1 : 0) +
                   (std::isnan(lowerLeft.value) ? 1 : 0) +
                   (std::isnan(lowerRight.value) ? 1 : 0)),
          borders(borders_), split(split_)
    {
        assert(upperLeft.y == upperRight.y);
        assert(lowerLeft.y == lowerRight.y);
        assert(lowerLeft.x == upperLeft.x);
        assert(lowerRight.x == upperRight.x);
        assert(!split || nanCount == 0);
    }

    Square upperLeftSquare() const
    {
        assert(!std::isnan(upperLeft.value));
        return Square(
            upperLeft, upperCenter(), leftCenter(), center(),
            (std::isnan(upperRight.value) ? RIGHT_BORDER : NO_BORDER) |
                (std::isnan(lowerLeft.value) ? LOWER_BORDER : NO_BORDER),
            true);
    }

    Square lowerLeftSquare() const
    {
        assert(!std::isnan(lowerLeft.value));
        return Square(
            leftCenter(), center(), lowerLeft, lowerCenter(),
            (std::isnan(lowerRight.value) ? RIGHT_BORDER : NO_BORDER) |
                (std::isnan(upperLeft.value) ? UPPER_BORDER : NO_BORDER),
            true);
    }

    Square lowerRightSquare() const
    {
        assert(!std::isnan(lowerRight.value));
        return Square(
            center(), rightCenter(), lowerCenter(), lowerRight,
            (std::isnan(lowerLeft.value) ? LEFT_BORDER : NO_BORDER) |
                (std::isnan(upperRight.value) ? UPPER_BORDER : NO_BORDER),
            true);
    }

    Square upperRightSquare() const
    {
        assert(!std::isnan(upperRight.value));
        return Square(
            upperCenter(), upperRight, center(), rightCenter(),
            (std::isnan(lowerRight.value) ? LOWER_BORDER : NO_BORDER) |
                (std::isnan(upperLeft.value) ? LEFT_BORDER : NO_BORDER),
            true);
    }

    double maxValue() const
    {
        assert(nanCount == 0);
        return std::max(std::max(upperLeft.value, upperRight.value),
                        std::max(lowerLeft.value, lowerRight.value));
    }

    double minValue() const
    {
        assert(nanCount == 0);
        return std::min(std::min(upperLeft.value, upperRight.value),
                        std::min(lowerLeft.value, lowerRight.value));
    }

    ValuedSegment segment(uint8_t border) const
    {
        switch (border)
        {
            case LEFT_BORDER:
                return ValuedSegment(upperLeft, lowerLeft);
            case LOWER_BORDER:
                return ValuedSegment(lowerLeft, lowerRight);
            case RIGHT_BORDER:
                return ValuedSegment(lowerRight, upperRight);
            case UPPER_BORDER:
                return ValuedSegment(upperRight, upperLeft);
        }
        assert(false);
        return ValuedSegment(upperLeft, upperLeft);
    }

    // returns segments of contour
    //
    // segments are oriented:
    //   - they form a vector from their first point to their second point.
    //   - when looking at the vector upward, values greater than the level are
    //   on the right
    //
    //     ^
    //  -  |  +
    Segments segments(double level) const
    {
        switch (marchingCase(level))
        {
            case (ALL_LOW):
                // debug("ALL_LOW");
                return Segments();
            case (ALL_HIGH):
                // debug("ALL_HIGH");
                return Segments();
            case (UPPER_LEFT):
                // debug("UPPER_LEFT");
                return Segments(Segment(interpolate(UPPER_BORDER, level),
                                        interpolate(LEFT_BORDER, level)));
            case (LOWER_LEFT):
                // debug("LOWER_LEFT");
                return Segments(Segment(interpolate(LEFT_BORDER, level),
                                        interpolate(LOWER_BORDER, level)));
            case (LOWER_RIGHT):
                // debug("LOWER_RIGHT");
                return Segments(Segment(interpolate(LOWER_BORDER, level),
                                        interpolate(RIGHT_BORDER, level)));
            case (UPPER_RIGHT):
                // debug("UPPER_RIGHT");
                return Segments(Segment(interpolate(RIGHT_BORDER, level),
                                        interpolate(UPPER_BORDER, level)));
            case (UPPER_LEFT | LOWER_LEFT):
                // debug("UPPER_LEFT | LOWER_LEFT");
                return Segments(Segment(interpolate(UPPER_BORDER, level),
                                        interpolate(LOWER_BORDER, level)));
            case (LOWER_LEFT | LOWER_RIGHT):
                // debug("LOWER_LEFT | LOWER_RIGHT");
                return Segments(Segment(interpolate(LEFT_BORDER, level),
                                        interpolate(RIGHT_BORDER, level)));
            case (LOWER_RIGHT | UPPER_RIGHT):
                // debug("LOWER_RIGHT | UPPER_RIGHT");
                return Segments(Segment(interpolate(LOWER_BORDER, level),
                                        interpolate(UPPER_BORDER, level)));
            case (UPPER_RIGHT | UPPER_LEFT):
                // debug("UPPER_RIGHT | UPPER_LEFT");
                return Segments(Segment(interpolate(RIGHT_BORDER, level),
                                        interpolate(LEFT_BORDER, level)));
            case (ALL_HIGH & ~UPPER_LEFT):
                // debug("ALL_HIGH & ~UPPER_LEFT");
                return Segments(Segment(interpolate(LEFT_BORDER, level),
                                        interpolate(UPPER_BORDER, level)));
            case (ALL_HIGH & ~LOWER_LEFT):
                // debug("ALL_HIGH & ~LOWER_LEFT");
                return Segments(Segment(interpolate(LOWER_BORDER, level),
                                        interpolate(LEFT_BORDER, level)));
            case (ALL_HIGH & ~LOWER_RIGHT):
                // debug("ALL_HIGH & ~LOWER_RIGHT");
                return Segments(Segment(interpolate(RIGHT_BORDER, level),
                                        interpolate(LOWER_BORDER, level)));
            case (ALL_HIGH & ~UPPER_RIGHT):
                // debug("ALL_HIGH & ~UPPER_RIGHT");
                return Segments(Segment(interpolate(UPPER_BORDER, level),
                                        interpolate(RIGHT_BORDER, level)));
            case (SADDLE_NE):
            case (SADDLE_NW):
                // From the two possible saddle configurations, we always return
                // the same one.
                //
                // The classical marching square algorithm says the ambiguity
                // should be resolved between the two possible configurations by
                // looking at the value of the center point. But in certain
                // cases, this may lead to line contours from different levels
                // that cross each other and then gives invalid polygons.
                //
                // Arbitrarily choosing one of the two possible configurations
                // is not really that worse than deciding based on the center
                // point.
                return Segments(Segment(interpolate(LEFT_BORDER, level),
                                        interpolate(LOWER_BORDER, level)),
                                Segment(interpolate(RIGHT_BORDER, level),
                                        interpolate(UPPER_BORDER, level)));
        }
        assert(false);
        return Segments();
    }

    template <typename Writer, typename LevelGenerator>
    void process(const LevelGenerator &levelGenerator, Writer &writer) const
    {
        if (nanCount == 4)  // nothing to do
            return;

        if (nanCount)  // split in 4
        {
            if (!std::isnan(upperLeft.value))
                upperLeftSquare().process(levelGenerator, writer);
            if (!std::isnan(upperRight.value))
                upperRightSquare().process(levelGenerator, writer);
            if (!std::isnan(lowerLeft.value))
                lowerLeftSquare().process(levelGenerator, writer);
            if (!std::isnan(lowerRight.value))
                lowerRightSquare().process(levelGenerator, writer);
            return;
        }

        if (writer.polygonize && borders)
        {
            for (uint8_t border :
                 {UPPER_BORDER, LEFT_BORDER, RIGHT_BORDER, LOWER_BORDER})
            {
                // bitwise AND to test which borders we have on the square
                if ((border & borders) == 0)
                    continue;

                // convention: for a level = L, store borders for the previous
                // level up to (and including) L in the border of level "L". For
                // fixed sets of level, this means there is an "Inf" slot for
                // borders of the highest level
                const ValuedSegment s(segment(border));

                Point lastPoint(s.first.x, s.first.y);
                Point endPoint(s.second.x, s.second.y);
                if (s.first.value > s.second.value)
                    std::swap(lastPoint, endPoint);
                bool reverse =
                    (s.first.value > s.second.value) &&
                    ((border == UPPER_BORDER) || (border == LEFT_BORDER));

                auto levelIt =
                    levelGenerator.range(s.first.value, s.second.value);

                auto it = levelIt.begin();  // reused after the for
                for (; it != levelIt.end(); ++it)
                {
                    const int levelIdx = (*it).first;
                    const double level = (*it).second;

                    const Point nextPoint(interpolate(border, level));
                    if (reverse)
                        writer.addBorderSegment(levelIdx, nextPoint, lastPoint);
                    else
                        writer.addBorderSegment(levelIdx, lastPoint, nextPoint);
                    lastPoint = nextPoint;
                }
                // last level (past the end)
                if (reverse)
                    writer.addBorderSegment((*it).first, endPoint, lastPoint);
                else
                    writer.addBorderSegment((*it).first, lastPoint, endPoint);
            }
        }

        auto range = levelGenerator.range(minValue(), maxValue());
        auto it = range.begin();
        auto itEnd = range.end();
        auto next = range.begin();
        ++next;

        for (; it != itEnd; ++it, ++next)
        {
            const int levelIdx = (*it).first;
            const double level = (*it).second;

            const Segments segments_ = segments(level);

            for (std::size_t i = 0; i < segments_.size(); i++)
            {
                const Segment &s = segments_[i];
                writer.addSegment(levelIdx, s.first, s.second);

                if (writer.polygonize)
                {
                    // the contour is used in the polygon of higher level as
                    // well
                    //
                    // TODO: copying the segment to the higher level is easy,
                    // but it involves too much memory. We should reuse segment
                    // contours when constructing polygon rings.
                    writer.addSegment((*next).first, s.first, s.second);
                }
            }
        }
    }

    const ValuedPoint upperLeft;
    const ValuedPoint lowerLeft;
    const ValuedPoint lowerRight;
    const ValuedPoint upperRight;
    const int nanCount;
    const uint8_t borders;
    const bool split;

  private:
    ValuedPoint center() const
    {
        return ValuedPoint(
            .5 * (upperLeft.x + lowerRight.x),
            .5 * (upperLeft.y + lowerRight.y),
            ((std::isnan(lowerLeft.value) ? 0 : lowerLeft.value) +
             (std::isnan(upperLeft.value) ? 0 : upperLeft.value) +
             (std::isnan(lowerRight.value) ? 0 : lowerRight.value) +
             (std::isnan(upperRight.value) ? 0 : upperRight.value)) /
                (4 - nanCount));
    }

    ValuedPoint leftCenter() const
    {
        return ValuedPoint(
            upperLeft.x, .5 * (upperLeft.y + lowerLeft.y),
            std::isnan(upperLeft.value)
                ? lowerLeft.value
                : (std::isnan(lowerLeft.value)
                       ? upperLeft.value
                       : .5 * (upperLeft.value + lowerLeft.value)));
    }

    ValuedPoint lowerCenter() const
    {
        return ValuedPoint(
            .5 * (lowerLeft.x + lowerRight.x), lowerLeft.y,
            std::isnan(lowerRight.value)
                ? lowerLeft.value
                : (std::isnan(lowerLeft.value)
                       ? lowerRight.value
                       : .5 * (lowerRight.value + lowerLeft.value)));
    }

    ValuedPoint rightCenter() const
    {
        return ValuedPoint(
            upperRight.x, .5 * (upperRight.y + lowerRight.y),
            std::isnan(lowerRight.value)
                ? upperRight.value
                : (std::isnan(upperRight.value)
                       ? lowerRight.value
                       : .5 * (lowerRight.value + upperRight.value)));
    }

    ValuedPoint upperCenter() const
    {
        return ValuedPoint(
            .5 * (upperLeft.x + upperRight.x), upperLeft.y,
            std::isnan(upperLeft.value)
                ? upperRight.value
                : (std::isnan(upperRight.value)
                       ? upperLeft.value
                       : .5 * (upperLeft.value + upperRight.value)));
    }

    uint8_t marchingCase(double level) const
    {
        return (level < fudge(upperLeft.value, level) ? UPPER_LEFT : ALL_LOW) |
               (level < fudge(lowerLeft.value, level) ? LOWER_LEFT : ALL_LOW) |
               (level < fudge(lowerRight.value, level) ? LOWER_RIGHT
                                                       : ALL_LOW) |
               (level < fudge(upperRight.value, level) ? UPPER_RIGHT : ALL_LOW);
    }

    static double interpolate_(double level, double x1, double x2, double y1,
                               double y2, bool need_split)
    {
        if (need_split)
        {
            // The two cases are here to avoid numerical roundup errors, for two
            // points, we always compute the same interpolation. This condition
            // is ensured by the order left->right bottom->top in interpole
            // calls
            //
            // To obtain the same value for border (split) and non-border
            // element, we take the middle value and interpolate from this to
            // the end
            const double xm = .5 * (x1 + x2);
            const double ym = .5 * (y1 + y2);
            const double fy1 = fudge(y1, level);
            const double fym = fudge(ym, level);
            if ((fy1 < level && level < fym) || (fy1 > level && level > fym))
            {
                x2 = xm;
                y2 = ym;
            }
            else
            {
                x1 = xm;
                y1 = ym;
            }
        }
        const double fy1 = fudge(y1, level);
        const double ratio = (level - fy1) / (fudge(y2, level) - fy1);
        return x1 * (1. - ratio) + x2 * ratio;
    }

    Point interpolate(uint8_t border, double level) const
    {
        switch (border)
        {
            case LEFT_BORDER:
                return Point(upperLeft.x,
                             interpolate_(level, lowerLeft.y, upperLeft.y,
                                          lowerLeft.value, upperLeft.value,
                                          !split));
            case LOWER_BORDER:
                return Point(interpolate_(level, lowerLeft.x, lowerRight.x,
                                          lowerLeft.value, lowerRight.value,
                                          !split),
                             lowerLeft.y);
            case RIGHT_BORDER:
                return Point(upperRight.x,
                             interpolate_(level, lowerRight.y, upperRight.y,
                                          lowerRight.value, upperRight.value,
                                          !split));
            case UPPER_BORDER:
                return Point(interpolate_(level, upperLeft.x, upperRight.x,
                                          upperLeft.value, upperRight.value,
                                          !split),
                             upperLeft.y);
        }
        assert(false);
        return Point();
    }
};
}  // 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_SQUARE_H
13		#define MARCHING_SQUARE_SQUARE_H
14
15		#include <algorithm>
16		#include <cassert>
17		#include <cmath>
18		#include <cstdint>
19		#include <math.h>
20		#include "utility.h"
21		#include "point.h"
22
23		namespace marching_squares
24		{
25
26		struct Square
27		{
28		// Bit flags to determine borders around pixel
29		static const uint8_t NO_BORDER = 0; // 0000 0000
30		static const uint8_t LEFT_BORDER = 1 << 0; // 0000 0001
31		static const uint8_t LOWER_BORDER = 1 << 1; // 0000 0010
32		static const uint8_t RIGHT_BORDER = 1 << 2; // 0000 0100
33		static const uint8_t UPPER_BORDER = 1 << 3; // 0000 1000
34
35		// Bit flags for marching square case
36		static const uint8_t ALL_LOW = 0; // 0000 0000
37		static const uint8_t UPPER_LEFT = 1 << 0; // 0000 0001
38		static const uint8_t LOWER_LEFT = 1 << 1; // 0000 0010
39		static const uint8_t LOWER_RIGHT = 1 << 2; // 0000 0100
40		static const uint8_t UPPER_RIGHT = 1 << 3; // 0000 1000
41		static const uint8_t ALL_HIGH =
42		UPPER_LEFT \| LOWER_LEFT \| LOWER_RIGHT \| UPPER_RIGHT; // 0000 1111
43		static const uint8_t SADDLE_NW = UPPER_LEFT \| LOWER_RIGHT; // 0000 0101
44		static const uint8_t SADDLE_NE = UPPER_RIGHT \| LOWER_LEFT; // 0000 1010
45
46		typedef std::pair<Point, Point> Segment;
47		typedef std::pair<ValuedPoint, ValuedPoint> ValuedSegment;
48
49		//
50		// An array of segments, at most 3 segments
51		struct Segments
52		{
53	0	Segments() : sz_(0)
54	0	{
55	0	}
56
57	0	Segments(const Segment &first) : sz_(1), segs_()
58	0	{
59	0	segs_[0] = first;
60	0	}
61
62	0	Segments(const Segment &first, const Segment &second) : sz_(2), segs_()
63	0	{
64	0	segs_[0] = first;
65	0	segs_[1] = second;
66	0	}
67
68		Segments(const Segment &first, const Segment &second,
69		const Segment &third)
70		: sz_(3), segs_()
71	0	{
72	0	segs_[0] = first;
73	0	segs_[1] = second;
74	0	segs_[2] = third;
75	0	}
76
77		std::size_t size() const
78	0	{
79	0	return sz_;
80	0	}
81
82		const Segment &operator[](std::size_t idx) const
83	0	{
84	0	assert(idx < sz_);
85	0	return segs_[idx];
86	0	}
87
88		private:
89		const std::size_t sz_;
90		/* const */ Segment segs_[3];
91		};
92
93		Square(const ValuedPoint &upperLeft_, const ValuedPoint &upperRight_,
94		const ValuedPoint &lowerLeft_, const ValuedPoint &lowerRight_,
95		uint8_t borders_ = NO_BORDER, bool split_ = false)
96	0	: upperLeft(upperLeft_), lowerLeft(lowerLeft_), lowerRight(lowerRight_),
97	0	upperRight(upperRight_),
98	0	nanCount((std::isnan(upperLeft.value) ? 1 : 0) +
99	0	(std::isnan(upperRight.value) ? 1 : 0) +
100	0	(std::isnan(lowerLeft.value) ? 1 : 0) +
101	0	(std::isnan(lowerRight.value) ? 1 : 0)),
102	0	borders(borders_), split(split_)
103	0	{
104	0	assert(upperLeft.y == upperRight.y);
105	0	assert(lowerLeft.y == lowerRight.y);
106	0	assert(lowerLeft.x == upperLeft.x);
107	0	assert(lowerRight.x == upperRight.x);
108	0	assert(!split \|\| nanCount == 0);
109	0	}
110
111		Square upperLeftSquare() const
112	0	{
113	0	assert(!std::isnan(upperLeft.value));
114	0	return Square(
115	0	upperLeft, upperCenter(), leftCenter(), center(),
116	0	(std::isnan(upperRight.value) ? RIGHT_BORDER : NO_BORDER) \|
117	0	(std::isnan(lowerLeft.value) ? LOWER_BORDER : NO_BORDER),
118	0	true);
119	0	}
120
121		Square lowerLeftSquare() const
122	0	{
123	0	assert(!std::isnan(lowerLeft.value));
124	0	return Square(
125	0	leftCenter(), center(), lowerLeft, lowerCenter(),
126	0	(std::isnan(lowerRight.value) ? RIGHT_BORDER : NO_BORDER) \|
127	0	(std::isnan(upperLeft.value) ? UPPER_BORDER : NO_BORDER),
128	0	true);
129	0	}
130
131		Square lowerRightSquare() const
132	0	{
133	0	assert(!std::isnan(lowerRight.value));
134	0	return Square(
135	0	center(), rightCenter(), lowerCenter(), lowerRight,
136	0	(std::isnan(lowerLeft.value) ? LEFT_BORDER : NO_BORDER) \|
137	0	(std::isnan(upperRight.value) ? UPPER_BORDER : NO_BORDER),
138	0	true);
139	0	}
140
141		Square upperRightSquare() const
142	0	{
143	0	assert(!std::isnan(upperRight.value));
144	0	return Square(
145	0	upperCenter(), upperRight, center(), rightCenter(),
146	0	(std::isnan(lowerRight.value) ? LOWER_BORDER : NO_BORDER) \|
147	0	(std::isnan(upperLeft.value) ? LEFT_BORDER : NO_BORDER),
148	0	true);
149	0	}
150
151		double maxValue() const
152	0	{
153	0	assert(nanCount == 0);
154	0	return std::max(std::max(upperLeft.value, upperRight.value),
155	0	std::max(lowerLeft.value, lowerRight.value));
156	0	}
157
158		double minValue() const
159	0	{
160	0	assert(nanCount == 0);
161	0	return std::min(std::min(upperLeft.value, upperRight.value),
162	0	std::min(lowerLeft.value, lowerRight.value));
163	0	}
164
165		ValuedSegment segment(uint8_t border) const
166	0	{
167	0	switch (border)
168	0	{
169	0	case LEFT_BORDER:
170	0	return ValuedSegment(upperLeft, lowerLeft);
171	0	case LOWER_BORDER:
172	0	return ValuedSegment(lowerLeft, lowerRight);
173	0	case RIGHT_BORDER:
174	0	return ValuedSegment(lowerRight, upperRight);
175	0	case UPPER_BORDER:
176	0	return ValuedSegment(upperRight, upperLeft);
177	0	}
178	0	assert(false);
179	0	return ValuedSegment(upperLeft, upperLeft);
180	0	}
181
182		// returns segments of contour
183		//
184		// segments are oriented:
185		// - they form a vector from their first point to their second point.
186		// - when looking at the vector upward, values greater than the level are
187		// on the right
188		//
189		// ^
190		// - \| +
191		Segments segments(double level) const
192	0	{
193	0	switch (marchingCase(level))
194	0	{
195	0	case (ALL_LOW):
196		// debug("ALL_LOW");
197	0	return Segments();
198	0	case (ALL_HIGH):
199		// debug("ALL_HIGH");
200	0	return Segments();
201	0	case (UPPER_LEFT):
202		// debug("UPPER_LEFT");
203	0	return Segments(Segment(interpolate(UPPER_BORDER, level),
204	0	interpolate(LEFT_BORDER, level)));
205	0	case (LOWER_LEFT):
206		// debug("LOWER_LEFT");
207	0	return Segments(Segment(interpolate(LEFT_BORDER, level),
208	0	interpolate(LOWER_BORDER, level)));
209	0	case (LOWER_RIGHT):
210		// debug("LOWER_RIGHT");
211	0	return Segments(Segment(interpolate(LOWER_BORDER, level),
212	0	interpolate(RIGHT_BORDER, level)));
213	0	case (UPPER_RIGHT):
214		// debug("UPPER_RIGHT");
215	0	return Segments(Segment(interpolate(RIGHT_BORDER, level),
216	0	interpolate(UPPER_BORDER, level)));
217	0	case (UPPER_LEFT \| LOWER_LEFT):
218		// debug("UPPER_LEFT \| LOWER_LEFT");
219	0	return Segments(Segment(interpolate(UPPER_BORDER, level),
220	0	interpolate(LOWER_BORDER, level)));
221	0	case (LOWER_LEFT \| LOWER_RIGHT):
222		// debug("LOWER_LEFT \| LOWER_RIGHT");
223	0	return Segments(Segment(interpolate(LEFT_BORDER, level),
224	0	interpolate(RIGHT_BORDER, level)));
225	0	case (LOWER_RIGHT \| UPPER_RIGHT):
226		// debug("LOWER_RIGHT \| UPPER_RIGHT");
227	0	return Segments(Segment(interpolate(LOWER_BORDER, level),
228	0	interpolate(UPPER_BORDER, level)));
229	0	case (UPPER_RIGHT \| UPPER_LEFT):
230		// debug("UPPER_RIGHT \| UPPER_LEFT");
231	0	return Segments(Segment(interpolate(RIGHT_BORDER, level),
232	0	interpolate(LEFT_BORDER, level)));
233	0	case (ALL_HIGH & ~UPPER_LEFT):
234		// debug("ALL_HIGH & ~UPPER_LEFT");
235	0	return Segments(Segment(interpolate(LEFT_BORDER, level),
236	0	interpolate(UPPER_BORDER, level)));
237	0	case (ALL_HIGH & ~LOWER_LEFT):
238		// debug("ALL_HIGH & ~LOWER_LEFT");
239	0	return Segments(Segment(interpolate(LOWER_BORDER, level),
240	0	interpolate(LEFT_BORDER, level)));
241	0	case (ALL_HIGH & ~LOWER_RIGHT):
242		// debug("ALL_HIGH & ~LOWER_RIGHT");
243	0	return Segments(Segment(interpolate(RIGHT_BORDER, level),
244	0	interpolate(LOWER_BORDER, level)));
245	0	case (ALL_HIGH & ~UPPER_RIGHT):
246		// debug("ALL_HIGH & ~UPPER_RIGHT");
247	0	return Segments(Segment(interpolate(UPPER_BORDER, level),
248	0	interpolate(RIGHT_BORDER, level)));
249	0	case (SADDLE_NE):
250	0	case (SADDLE_NW):
251		// From the two possible saddle configurations, we always return
252		// the same one.
253		//
254		// The classical marching square algorithm says the ambiguity
255		// should be resolved between the two possible configurations by
256		// looking at the value of the center point. But in certain
257		// cases, this may lead to line contours from different levels
258		// that cross each other and then gives invalid polygons.
259		//
260		// Arbitrarily choosing one of the two possible configurations
261		// is not really that worse than deciding based on the center
262		// point.
263	0	return Segments(Segment(interpolate(LEFT_BORDER, level),
264	0	interpolate(LOWER_BORDER, level)),
265	0	Segment(interpolate(RIGHT_BORDER, level),
266	0	interpolate(UPPER_BORDER, level)));
267	0	}
268	0	assert(false);
269	0	return Segments();
270	0	}
271
272		template <typename Writer, typename LevelGenerator>
273		void process(const LevelGenerator &levelGenerator, Writer &writer) const
274	0	{
275	0	if (nanCount == 4) // nothing to do
276	0	return;
277
278	0	if (nanCount) // split in 4
279	0	{
280	0	if (!std::isnan(upperLeft.value))
281	0	upperLeftSquare().process(levelGenerator, writer);
282	0	if (!std::isnan(upperRight.value))
283	0	upperRightSquare().process(levelGenerator, writer);
284	0	if (!std::isnan(lowerLeft.value))
285	0	lowerLeftSquare().process(levelGenerator, writer);
286	0	if (!std::isnan(lowerRight.value))
287	0	lowerRightSquare().process(levelGenerator, writer);
288	0	return;
289	0	}
290
291	0	if (writer.polygonize && borders)
292	0	{
293	0	for (uint8_t border :
294	0	{UPPER_BORDER, LEFT_BORDER, RIGHT_BORDER, LOWER_BORDER})
295	0	{
296		// bitwise AND to test which borders we have on the square
297	0	if ((border & borders) == 0)
298	0	continue;
299
300		// convention: for a level = L, store borders for the previous
301		// level up to (and including) L in the border of level "L". For
302		// fixed sets of level, this means there is an "Inf" slot for
303		// borders of the highest level
304	0	const ValuedSegment s(segment(border));
305
306	0	Point lastPoint(s.first.x, s.first.y);
307	0	Point endPoint(s.second.x, s.second.y);
308	0	if (s.first.value > s.second.value)
309	0	std::swap(lastPoint, endPoint);
310	0	bool reverse =
311	0	(s.first.value > s.second.value) &&
312	0	((border == UPPER_BORDER) \|\| (border == LEFT_BORDER));
313
314	0	auto levelIt =
315	0	levelGenerator.range(s.first.value, s.second.value);
316
317	0	auto it = levelIt.begin(); // reused after the for
318	0	for (; it != levelIt.end(); ++it)
319	0	{
320	0	const int levelIdx = (*it).first;
321	0	const double level = (*it).second;
322
323	0	const Point nextPoint(interpolate(border, level));
324	0	if (reverse)
325	0	writer.addBorderSegment(levelIdx, nextPoint, lastPoint);
326	0	else
327	0	writer.addBorderSegment(levelIdx, lastPoint, nextPoint);
328	0	lastPoint = nextPoint;
329	0	}
330		// last level (past the end)
331	0	if (reverse)
332	0	writer.addBorderSegment((*it).first, endPoint, lastPoint);
333	0	else
334	0	writer.addBorderSegment((*it).first, lastPoint, endPoint);
335	0	}
336	0	}
337
338	0	auto range = levelGenerator.range(minValue(), maxValue());
339	0	auto it = range.begin();
340	0	auto itEnd = range.end();
341	0	auto next = range.begin();
342	0	++next;
343
344	0	for (; it != itEnd; ++it, ++next)
345	0	{
346	0	const int levelIdx = (*it).first;
347	0	const double level = (*it).second;
348
349	0	const Segments segments_ = segments(level);
350
351	0	for (std::size_t i = 0; i < segments_.size(); i++)
352	0	{
353	0	const Segment &s = segments_[i];
354	0	writer.addSegment(levelIdx, s.first, s.second);
355
356	0	if (writer.polygonize)
357	0	{
358		// the contour is used in the polygon of higher level as
359		// well
360		//
361		// TODO: copying the segment to the higher level is easy,
362		// but it involves too much memory. We should reuse segment
363		// contours when constructing polygon rings.
364	0	writer.addSegment((*next).first, s.first, s.second);
365	0	}
366	0	}
367	0	}
368	0	} Unexecuted instantiation: void marching_squares::Square::process<marching_squares::SegmentMerger<marching_squares::PolygonRingAppender<PolygonContourWriter>, marching_squares::FixedLevelRangeIterator>, marching_squares::FixedLevelRangeIterator>(marching_squares::FixedLevelRangeIterator const&, marching_squares::SegmentMerger<marching_squares::PolygonRingAppender<PolygonContourWriter>, marching_squares::FixedLevelRangeIterator>&) const Unexecuted instantiation: void marching_squares::Square::process<marching_squares::SegmentMerger<GDALRingAppender, marching_squares::FixedLevelRangeIterator>, marching_squares::FixedLevelRangeIterator>(marching_squares::FixedLevelRangeIterator const&, marching_squares::SegmentMerger<GDALRingAppender, marching_squares::FixedLevelRangeIterator>&) const Unexecuted instantiation: void marching_squares::Square::process<marching_squares::SegmentMerger<GDALRingAppender, marching_squares::IntervalLevelRangeIterator>, marching_squares::IntervalLevelRangeIterator>(marching_squares::IntervalLevelRangeIterator const&, marching_squares::SegmentMerger<GDALRingAppender, marching_squares::IntervalLevelRangeIterator>&) const
369
370		const ValuedPoint upperLeft;
371		const ValuedPoint lowerLeft;
372		const ValuedPoint lowerRight;
373		const ValuedPoint upperRight;
374		const int nanCount;
375		const uint8_t borders;
376		const bool split;
377
378		private:
379		ValuedPoint center() const
380	0	{
381	0	return ValuedPoint(
382	0	.5 * (upperLeft.x + lowerRight.x),
383	0	.5 * (upperLeft.y + lowerRight.y),
384	0	((std::isnan(lowerLeft.value) ? 0 : lowerLeft.value) +
385	0	(std::isnan(upperLeft.value) ? 0 : upperLeft.value) +
386	0	(std::isnan(lowerRight.value) ? 0 : lowerRight.value) +
387	0	(std::isnan(upperRight.value) ? 0 : upperRight.value)) /
388	0	(4 - nanCount));
389	0	}
390
391		ValuedPoint leftCenter() const
392	0	{
393	0	return ValuedPoint(
394	0	upperLeft.x, .5 * (upperLeft.y + lowerLeft.y),
395	0	std::isnan(upperLeft.value)
396	0	? lowerLeft.value
397	0	: (std::isnan(lowerLeft.value)
398	0	? upperLeft.value
399	0	: .5 * (upperLeft.value + lowerLeft.value)));
400	0	}
401
402		ValuedPoint lowerCenter() const
403	0	{
404	0	return ValuedPoint(
405	0	.5 * (lowerLeft.x + lowerRight.x), lowerLeft.y,
406	0	std::isnan(lowerRight.value)
407	0	? lowerLeft.value
408	0	: (std::isnan(lowerLeft.value)
409	0	? lowerRight.value
410	0	: .5 * (lowerRight.value + lowerLeft.value)));
411	0	}
412
413		ValuedPoint rightCenter() const
414	0	{
415	0	return ValuedPoint(
416	0	upperRight.x, .5 * (upperRight.y + lowerRight.y),
417	0	std::isnan(lowerRight.value)
418	0	? upperRight.value
419	0	: (std::isnan(upperRight.value)
420	0	? lowerRight.value
421	0	: .5 * (lowerRight.value + upperRight.value)));
422	0	}
423
424		ValuedPoint upperCenter() const
425	0	{
426	0	return ValuedPoint(
427	0	.5 * (upperLeft.x + upperRight.x), upperLeft.y,
428	0	std::isnan(upperLeft.value)
429	0	? upperRight.value
430	0	: (std::isnan(upperRight.value)
431	0	? upperLeft.value
432	0	: .5 * (upperLeft.value + upperRight.value)));
433	0	}
434
435		uint8_t marchingCase(double level) const
436	0	{
437	0	return (level < fudge(upperLeft.value, level) ? UPPER_LEFT : ALL_LOW) \|
438	0	(level < fudge(lowerLeft.value, level) ? LOWER_LEFT : ALL_LOW) \|
439	0	(level < fudge(lowerRight.value, level) ? LOWER_RIGHT
440	0	: ALL_LOW) \|
441	0	(level < fudge(upperRight.value, level) ? UPPER_RIGHT : ALL_LOW);
442	0	}
443
444		static double interpolate_(double level, double x1, double x2, double y1,
445		double y2, bool need_split)
446	0	{
447	0	if (need_split)
448	0	{
449		// The two cases are here to avoid numerical roundup errors, for two
450		// points, we always compute the same interpolation. This condition
451		// is ensured by the order left->right bottom->top in interpole
452		// calls
453		//
454		// To obtain the same value for border (split) and non-border
455		// element, we take the middle value and interpolate from this to
456		// the end
457	0	const double xm = .5 * (x1 + x2);
458	0	const double ym = .5 * (y1 + y2);
459	0	const double fy1 = fudge(y1, level);
460	0	const double fym = fudge(ym, level);
461	0	if ((fy1 < level && level < fym) \|\| (fy1 > level && level > fym))
462	0	{
463	0	x2 = xm;
464	0	y2 = ym;
465	0	}
466	0	else
467	0	{
468	0	x1 = xm;
469	0	y1 = ym;
470	0	}
471	0	}
472	0	const double fy1 = fudge(y1, level);
473	0	const double ratio = (level - fy1) / (fudge(y2, level) - fy1);
474	0	return x1 * (1. - ratio) + x2 * ratio;
475	0	}
476
477		Point interpolate(uint8_t border, double level) const
478	0	{
479	0	switch (border)
480	0	{
481	0	case LEFT_BORDER:
482	0	return Point(upperLeft.x,
483	0	interpolate_(level, lowerLeft.y, upperLeft.y,
484	0	lowerLeft.value, upperLeft.value,
485	0	!split));
486	0	case LOWER_BORDER:
487	0	return Point(interpolate_(level, lowerLeft.x, lowerRight.x,
488	0	lowerLeft.value, lowerRight.value,
489	0	!split),
490	0	lowerLeft.y);
491	0	case RIGHT_BORDER:
492	0	return Point(upperRight.x,
493	0	interpolate_(level, lowerRight.y, upperRight.y,
494	0	lowerRight.value, upperRight.value,
495	0	!split));
496	0	case UPPER_BORDER:
497	0	return Point(interpolate_(level, upperLeft.x, upperRight.x,
498	0	upperLeft.value, upperRight.value,
499	0	!split),
500	0	upperLeft.y);
501	0	}
502	0	assert(false);
503	0	return Point();
504	0	}
505		};
506		} // namespace marching_squares
507		#endif