square.h

/******************************************************************************
 *
 * 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>
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 ****************************************************************************/
#ifndef MARCHING_SQUARE_SQUARE_H
#define MARCHING_SQUARE_SQUARE_H

#include <algorithm>
#include <cassert>
#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(level, upperLeft.value) ? UPPER_LEFT : ALL_LOW) |
               (level < fudge(level, lowerLeft.value) ? LOWER_LEFT : ALL_LOW) |
               (level < fudge(level, lowerRight.value) ? LOWER_RIGHT
                                                       : ALL_LOW) |
               (level < fudge(level, upperRight.value) ? 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(level, y1);
            const double fym = fudge(level, ym);
            if ((fy1 < level && level < fym) || (fy1 > level && level > fym))
            {
                x2 = xm;
                y2 = ym;
            }
            else
            {
                x1 = xm;
                y1 = ym;
            }
        }
        const double fy1 = fudge(level, y1);
        const double ratio = (level - fy1) / (fudge(level, y2) - 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