/src/gdal/alg/marching_squares/square.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_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, double minLevel) const
    {
        switch (marchingCase(level, minLevel))
        {
            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, minLevel),
                            interpolate(LEFT_BORDER, level, minLevel)));
            case (LOWER_LEFT):
                // debug("LOWER_LEFT");
                return Segments(
                    Segment(interpolate(LEFT_BORDER, level, minLevel),
                            interpolate(LOWER_BORDER, level, minLevel)));
            case (LOWER_RIGHT):
                // debug("LOWER_RIGHT");
                return Segments(
                    Segment(interpolate(LOWER_BORDER, level, minLevel),
                            interpolate(RIGHT_BORDER, level, minLevel)));
            case (UPPER_RIGHT):
                // debug("UPPER_RIGHT");
                return Segments(
                    Segment(interpolate(RIGHT_BORDER, level, minLevel),
                            interpolate(UPPER_BORDER, level, minLevel)));
            case (UPPER_LEFT | LOWER_LEFT):
                // debug("UPPER_LEFT | LOWER_LEFT");
                return Segments(
                    Segment(interpolate(UPPER_BORDER, level, minLevel),
                            interpolate(LOWER_BORDER, level, minLevel)));
            case (LOWER_LEFT | LOWER_RIGHT):
                // debug("LOWER_LEFT | LOWER_RIGHT");
                return Segments(
                    Segment(interpolate(LEFT_BORDER, level, minLevel),
                            interpolate(RIGHT_BORDER, level, minLevel)));
            case (LOWER_RIGHT | UPPER_RIGHT):
                // debug("LOWER_RIGHT | UPPER_RIGHT");
                return Segments(
                    Segment(interpolate(LOWER_BORDER, level, minLevel),
                            interpolate(UPPER_BORDER, level, minLevel)));
            case (UPPER_RIGHT | UPPER_LEFT):
                // debug("UPPER_RIGHT | UPPER_LEFT");
                return Segments(
                    Segment(interpolate(RIGHT_BORDER, level, minLevel),
                            interpolate(LEFT_BORDER, level, minLevel)));
            case (ALL_HIGH & ~UPPER_LEFT):
                // debug("ALL_HIGH & ~UPPER_LEFT");
                return Segments(
                    Segment(interpolate(LEFT_BORDER, level, minLevel),
                            interpolate(UPPER_BORDER, level, minLevel)));
            case (ALL_HIGH & ~LOWER_LEFT):
                // debug("ALL_HIGH & ~LOWER_LEFT");
                return Segments(
                    Segment(interpolate(LOWER_BORDER, level, minLevel),
                            interpolate(LEFT_BORDER, level, minLevel)));
            case (ALL_HIGH & ~LOWER_RIGHT):
                // debug("ALL_HIGH & ~LOWER_RIGHT");
                return Segments(
                    Segment(interpolate(RIGHT_BORDER, level, minLevel),
                            interpolate(LOWER_BORDER, level, minLevel)));
            case (ALL_HIGH & ~UPPER_RIGHT):
                // debug("ALL_HIGH & ~UPPER_RIGHT");
                return Segments(
                    Segment(interpolate(UPPER_BORDER, level, minLevel),
                            interpolate(RIGHT_BORDER, level, minLevel)));
            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, minLevel),
                            interpolate(LOWER_BORDER, level, minLevel)),
                    Segment(interpolate(RIGHT_BORDER, level, minLevel),
                            interpolate(UPPER_BORDER, level, minLevel)));
        }
        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, levelGenerator.minLevel()));
                    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, levelGenerator.minLevel());

            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, double minLevel) const
    {
        return (level < fudge(upperLeft.value, minLevel, level) ? UPPER_LEFT
                                                                : ALL_LOW) |
               (level < fudge(lowerLeft.value, minLevel, level) ? LOWER_LEFT
                                                                : ALL_LOW) |
               (level < fudge(lowerRight.value, minLevel, level) ? LOWER_RIGHT
                                                                 : ALL_LOW) |
               (level < fudge(upperRight.value, minLevel, level) ? UPPER_RIGHT
                                                                 : ALL_LOW);
    }

    static double interpolate_(double level, double x1, double x2, double y1,
                               double y2, bool need_split, double minLevel)
    {
        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, minLevel, level);
            const double fym = fudge(ym, minLevel, level);
            if ((fy1 < level && level < fym) || (fy1 > level && level > fym))
            {
                x2 = xm;
                y2 = ym;
            }
            else
            {
                x1 = xm;
                y1 = ym;
            }
        }
        const double fy1 = fudge(y1, minLevel, level);
        const double ratio = (level - fy1) / (fudge(y2, minLevel, level) - fy1);
        return x1 * (1. - ratio) + x2 * ratio;
    }

    Point interpolate(uint8_t border, double level, double minLevel) const
    {
        switch (border)
        {
            case LEFT_BORDER:
                return Point(upperLeft.x,
                             interpolate_(level, lowerLeft.y, upperLeft.y,
                                          lowerLeft.value, upperLeft.value,
                                          !split, minLevel));
            case LOWER_BORDER:
                return Point(interpolate_(level, lowerLeft.x, lowerRight.x,
                                          lowerLeft.value, lowerRight.value,
                                          !split, minLevel),
                             lowerLeft.y);
            case RIGHT_BORDER:
                return Point(upperRight.x,
                             interpolate_(level, lowerRight.y, upperRight.y,
                                          lowerRight.value, upperRight.value,
                                          !split, minLevel));
            case UPPER_BORDER:
                return Point(interpolate_(level, upperLeft.x, upperRight.x,
                                          upperLeft.value, upperRight.value,
                                          !split, minLevel),
                             upperLeft.y);
        }
        assert(false);
        return Point();
    }
};
}  // namespace marching_squares
#endif

Coverage Report

Created: 2025-06-13 06:18

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_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, double minLevel) const
192	0	{
193	0	switch (marchingCase(level, minLevel))
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(
204	0	Segment(interpolate(UPPER_BORDER, level, minLevel),
205	0	interpolate(LEFT_BORDER, level, minLevel)));
206	0	case (LOWER_LEFT):
207		// debug("LOWER_LEFT");
208	0	return Segments(
209	0	Segment(interpolate(LEFT_BORDER, level, minLevel),
210	0	interpolate(LOWER_BORDER, level, minLevel)));
211	0	case (LOWER_RIGHT):
212		// debug("LOWER_RIGHT");
213	0	return Segments(
214	0	Segment(interpolate(LOWER_BORDER, level, minLevel),
215	0	interpolate(RIGHT_BORDER, level, minLevel)));
216	0	case (UPPER_RIGHT):
217		// debug("UPPER_RIGHT");
218	0	return Segments(
219	0	Segment(interpolate(RIGHT_BORDER, level, minLevel),
220	0	interpolate(UPPER_BORDER, level, minLevel)));
221	0	case (UPPER_LEFT \| LOWER_LEFT):
222		// debug("UPPER_LEFT \| LOWER_LEFT");
223	0	return Segments(
224	0	Segment(interpolate(UPPER_BORDER, level, minLevel),
225	0	interpolate(LOWER_BORDER, level, minLevel)));
226	0	case (LOWER_LEFT \| LOWER_RIGHT):
227		// debug("LOWER_LEFT \| LOWER_RIGHT");
228	0	return Segments(
229	0	Segment(interpolate(LEFT_BORDER, level, minLevel),
230	0	interpolate(RIGHT_BORDER, level, minLevel)));
231	0	case (LOWER_RIGHT \| UPPER_RIGHT):
232		// debug("LOWER_RIGHT \| UPPER_RIGHT");
233	0	return Segments(
234	0	Segment(interpolate(LOWER_BORDER, level, minLevel),
235	0	interpolate(UPPER_BORDER, level, minLevel)));
236	0	case (UPPER_RIGHT \| UPPER_LEFT):
237		// debug("UPPER_RIGHT \| UPPER_LEFT");
238	0	return Segments(
239	0	Segment(interpolate(RIGHT_BORDER, level, minLevel),
240	0	interpolate(LEFT_BORDER, level, minLevel)));
241	0	case (ALL_HIGH & ~UPPER_LEFT):
242		// debug("ALL_HIGH & ~UPPER_LEFT");
243	0	return Segments(
244	0	Segment(interpolate(LEFT_BORDER, level, minLevel),
245	0	interpolate(UPPER_BORDER, level, minLevel)));
246	0	case (ALL_HIGH & ~LOWER_LEFT):
247		// debug("ALL_HIGH & ~LOWER_LEFT");
248	0	return Segments(
249	0	Segment(interpolate(LOWER_BORDER, level, minLevel),
250	0	interpolate(LEFT_BORDER, level, minLevel)));
251	0	case (ALL_HIGH & ~LOWER_RIGHT):
252		// debug("ALL_HIGH & ~LOWER_RIGHT");
253	0	return Segments(
254	0	Segment(interpolate(RIGHT_BORDER, level, minLevel),
255	0	interpolate(LOWER_BORDER, level, minLevel)));
256	0	case (ALL_HIGH & ~UPPER_RIGHT):
257		// debug("ALL_HIGH & ~UPPER_RIGHT");
258	0	return Segments(
259	0	Segment(interpolate(UPPER_BORDER, level, minLevel),
260	0	interpolate(RIGHT_BORDER, level, minLevel)));
261	0	case (SADDLE_NE):
262	0	case (SADDLE_NW):
263		// From the two possible saddle configurations, we always return
264		// the same one.
265		//
266		// The classical marching square algorithm says the ambiguity
267		// should be resolved between the two possible configurations by
268		// looking at the value of the center point. But in certain
269		// cases, this may lead to line contours from different levels
270		// that cross each other and then gives invalid polygons.
271		//
272		// Arbitrarily choosing one of the two possible configurations
273		// is not really that worse than deciding based on the center
274		// point.
275	0	return Segments(
276	0	Segment(interpolate(LEFT_BORDER, level, minLevel),
277	0	interpolate(LOWER_BORDER, level, minLevel)),
278	0	Segment(interpolate(RIGHT_BORDER, level, minLevel),
279	0	interpolate(UPPER_BORDER, level, minLevel)));
280	0	}
281	0	assert(false);
282	0	return Segments();
283	0	}
284
285		template <typename Writer, typename LevelGenerator>
286		void process(const LevelGenerator &levelGenerator, Writer &writer) const
287	0	{
288	0	if (nanCount == 4) // nothing to do
289	0	return;
290
291	0	if (nanCount) // split in 4
292	0	{
293	0	if (!std::isnan(upperLeft.value))
294	0	upperLeftSquare().process(levelGenerator, writer);
295	0	if (!std::isnan(upperRight.value))
296	0	upperRightSquare().process(levelGenerator, writer);
297	0	if (!std::isnan(lowerLeft.value))
298	0	lowerLeftSquare().process(levelGenerator, writer);
299	0	if (!std::isnan(lowerRight.value))
300	0	lowerRightSquare().process(levelGenerator, writer);
301	0	return;
302	0	}
303
304	0	if (writer.polygonize && borders)
305	0	{
306	0	for (uint8_t border :
307	0	{UPPER_BORDER, LEFT_BORDER, RIGHT_BORDER, LOWER_BORDER})
308	0	{
309		// bitwise AND to test which borders we have on the square
310	0	if ((border & borders) == 0)
311	0	continue;
312
313		// convention: for a level = L, store borders for the previous
314		// level up to (and including) L in the border of level "L". For
315		// fixed sets of level, this means there is an "Inf" slot for
316		// borders of the highest level
317	0	const ValuedSegment s(segment(border));
318
319	0	Point lastPoint(s.first.x, s.first.y);
320	0	Point endPoint(s.second.x, s.second.y);
321	0	if (s.first.value > s.second.value)
322	0	std::swap(lastPoint, endPoint);
323	0	bool reverse =
324	0	(s.first.value > s.second.value) &&
325	0	((border == UPPER_BORDER) \|\| (border == LEFT_BORDER));
326
327	0	auto levelIt =
328	0	levelGenerator.range(s.first.value, s.second.value);
329
330	0	auto it = levelIt.begin(); // reused after the for
331	0	for (; it != levelIt.end(); ++it)
332	0	{
333	0	const int levelIdx = (*it).first;
334	0	const double level = (*it).second;
335
336	0	const Point nextPoint(
337	0	interpolate(border, level, levelGenerator.minLevel()));
338	0	if (reverse)
339	0	writer.addBorderSegment(levelIdx, nextPoint, lastPoint);
340	0	else
341	0	writer.addBorderSegment(levelIdx, lastPoint, nextPoint);
342	0	lastPoint = nextPoint;
343	0	}
344		// last level (past the end)
345	0	if (reverse)
346	0	writer.addBorderSegment((*it).first, endPoint, lastPoint);
347	0	else
348	0	writer.addBorderSegment((*it).first, lastPoint, endPoint);
349	0	}
350	0	}
351
352	0	auto range = levelGenerator.range(minValue(), maxValue());
353	0	auto it = range.begin();
354	0	auto itEnd = range.end();
355	0	auto next = range.begin();
356	0	++next;
357
358	0	for (; it != itEnd; ++it, ++next)
359	0	{
360	0	const int levelIdx = (*it).first;
361	0	const double level = (*it).second;
362
363	0	const Segments segments_ =
364	0	segments(level, levelGenerator.minLevel());
365
366	0	for (std::size_t i = 0; i < segments_.size(); i++)
367	0	{
368	0	const Segment &s = segments_[i];
369	0	writer.addSegment(levelIdx, s.first, s.second);
370
371	0	if (writer.polygonize)
372	0	{
373		// the contour is used in the polygon of higher level as
374		// well
375		//
376		// TODO: copying the segment to the higher level is easy,
377		// but it involves too much memory. We should reuse segment
378		// contours when constructing polygon rings.
379	0	writer.addSegment((*next).first, s.first, s.second);
380	0	}
381	0	}
382	0	}
383	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
384
385		const ValuedPoint upperLeft;
386		const ValuedPoint lowerLeft;
387		const ValuedPoint lowerRight;
388		const ValuedPoint upperRight;
389		const int nanCount;
390		const uint8_t borders;
391		const bool split;
392
393		private:
394		ValuedPoint center() const
395	0	{
396	0	return ValuedPoint(
397	0	.5 * (upperLeft.x + lowerRight.x),
398	0	.5 * (upperLeft.y + lowerRight.y),
399	0	((std::isnan(lowerLeft.value) ? 0 : lowerLeft.value) +
400	0	(std::isnan(upperLeft.value) ? 0 : upperLeft.value) +
401	0	(std::isnan(lowerRight.value) ? 0 : lowerRight.value) +
402	0	(std::isnan(upperRight.value) ? 0 : upperRight.value)) /
403	0	(4 - nanCount));
404	0	}
405
406		ValuedPoint leftCenter() const
407	0	{
408	0	return ValuedPoint(
409	0	upperLeft.x, .5 * (upperLeft.y + lowerLeft.y),
410	0	std::isnan(upperLeft.value)
411	0	? lowerLeft.value
412	0	: (std::isnan(lowerLeft.value)
413	0	? upperLeft.value
414	0	: .5 * (upperLeft.value + lowerLeft.value)));
415	0	}
416
417		ValuedPoint lowerCenter() const
418	0	{
419	0	return ValuedPoint(
420	0	.5 * (lowerLeft.x + lowerRight.x), lowerLeft.y,
421	0	std::isnan(lowerRight.value)
422	0	? lowerLeft.value
423	0	: (std::isnan(lowerLeft.value)
424	0	? lowerRight.value
425	0	: .5 * (lowerRight.value + lowerLeft.value)));
426	0	}
427
428		ValuedPoint rightCenter() const
429	0	{
430	0	return ValuedPoint(
431	0	upperRight.x, .5 * (upperRight.y + lowerRight.y),
432	0	std::isnan(lowerRight.value)
433	0	? upperRight.value
434	0	: (std::isnan(upperRight.value)
435	0	? lowerRight.value
436	0	: .5 * (lowerRight.value + upperRight.value)));
437	0	}
438
439		ValuedPoint upperCenter() const
440	0	{
441	0	return ValuedPoint(
442	0	.5 * (upperLeft.x + upperRight.x), upperLeft.y,
443	0	std::isnan(upperLeft.value)
444	0	? upperRight.value
445	0	: (std::isnan(upperRight.value)
446	0	? upperLeft.value
447	0	: .5 * (upperLeft.value + upperRight.value)));
448	0	}
449
450		uint8_t marchingCase(double level, double minLevel) const
451	0	{
452	0	return (level < fudge(upperLeft.value, minLevel, level) ? UPPER_LEFT
453	0	: ALL_LOW) \|
454	0	(level < fudge(lowerLeft.value, minLevel, level) ? LOWER_LEFT
455	0	: ALL_LOW) \|
456	0	(level < fudge(lowerRight.value, minLevel, level) ? LOWER_RIGHT
457	0	: ALL_LOW) \|
458	0	(level < fudge(upperRight.value, minLevel, level) ? UPPER_RIGHT
459	0	: ALL_LOW);
460	0	}
461
462		static double interpolate_(double level, double x1, double x2, double y1,
463		double y2, bool need_split, double minLevel)
464	0	{
465	0	if (need_split)
466	0	{
467		// The two cases are here to avoid numerical roundup errors, for two
468		// points, we always compute the same interpolation. This condition
469		// is ensured by the order left->right bottom->top in interpole
470		// calls
471		//
472		// To obtain the same value for border (split) and non-border
473		// element, we take the middle value and interpolate from this to
474		// the end
475	0	const double xm = .5 * (x1 + x2);
476	0	const double ym = .5 * (y1 + y2);
477	0	const double fy1 = fudge(y1, minLevel, level);
478	0	const double fym = fudge(ym, minLevel, level);
479	0	if ((fy1 < level && level < fym) \|\| (fy1 > level && level > fym))
480	0	{
481	0	x2 = xm;
482	0	y2 = ym;
483	0	}
484	0	else
485	0	{
486	0	x1 = xm;
487	0	y1 = ym;
488	0	}
489	0	}
490	0	const double fy1 = fudge(y1, minLevel, level);
491	0	const double ratio = (level - fy1) / (fudge(y2, minLevel, level) - fy1);
492	0	return x1 * (1. - ratio) + x2 * ratio;
493	0	}
494
495		Point interpolate(uint8_t border, double level, double minLevel) const
496	0	{
497	0	switch (border)
498	0	{
499	0	case LEFT_BORDER:
500	0	return Point(upperLeft.x,
501	0	interpolate_(level, lowerLeft.y, upperLeft.y,
502	0	lowerLeft.value, upperLeft.value,
503	0	!split, minLevel));
504	0	case LOWER_BORDER:
505	0	return Point(interpolate_(level, lowerLeft.x, lowerRight.x,
506	0	lowerLeft.value, lowerRight.value,
507	0	!split, minLevel),
508	0	lowerLeft.y);
509	0	case RIGHT_BORDER:
510	0	return Point(upperRight.x,
511	0	interpolate_(level, lowerRight.y, upperRight.y,
512	0	lowerRight.value, upperRight.value,
513	0	!split, minLevel));
514	0	case UPPER_BORDER:
515	0	return Point(interpolate_(level, upperLeft.x, upperRight.x,
516	0	upperLeft.value, upperRight.value,
517	0	!split, minLevel),
518	0	upperLeft.y);
519	0	}
520	0	assert(false);
521	0	return Point();
522	0	}
523		};
524		} // namespace marching_squares
525		#endif