CoordinateSequence.h

/**********************************************************************
 *
 * GEOS - Geometry Engine Open Source
 * http://geos.osgeo.org
 *
 * Copyright (C) 2022 ISciences LLC
 * Copyright (C) 2006 Refractions Research Inc.
 *
 * This is free software; you can redistribute and/or modify it under
 * the terms of the GNU Lesser General Public Licence as published
 * by the Free Software Foundation.
 * See the COPYING file for more information.
 *
 **********************************************************************/
 
#pragma once
 
#include <geos/export.h>
 
#include <geos/geom/Coordinate.h> // for applyCoordinateFilter
#include <geos/geom/CoordinateSequenceIterator.h>
 
#include <cassert>
#include <vector>
#include <iostream>
#include <iosfwd> // ostream
#include <memory> // for unique_ptr typedef
 
// Forward declarations
namespace geos {
namespace geom {
class Envelope;
class CoordinateFilter;
}
}
 
namespace geos {
namespace geom { // geos::geom
 
class GEOS_DLL CoordinateSequence {
 
public:
 
    enum { X, Y, Z, M };
 
    using iterator = CoordinateSequenceIterator<CoordinateSequence, Coordinate>;
    using const_iterator = CoordinateSequenceIterator<const CoordinateSequence, const Coordinate>;
 
    typedef std::unique_ptr<CoordinateSequence> Ptr;
 
 
    CoordinateSequence();
 
    CoordinateSequence(std::size_t size, std::size_t dim = 0);
 
    CoordinateSequence(std::size_t size, bool hasz, bool hasm, bool initialize = true);
 
    CoordinateSequence(const std::initializer_list<Coordinate>&);
 
    CoordinateSequence(const std::initializer_list<CoordinateXY>&);
 
    CoordinateSequence(const std::initializer_list<CoordinateXYM>&);
 
    CoordinateSequence(const std::initializer_list<CoordinateXYZM>&);
 
    static CoordinateSequence XY(std::size_t size)  {
        return CoordinateSequence(size, false, false);
    }
 
    static CoordinateSequence XYZ(std::size_t size)  {
        return CoordinateSequence(size, true, false);
    }
 
    static CoordinateSequence XYZM(std::size_t size)  {
        return CoordinateSequence(size, true, true);
    }
 
    static CoordinateSequence XYM(std::size_t size)  {
        return CoordinateSequence(size, false, true);
    }
 
    std::unique_ptr<CoordinateSequence> clone() const;
 
 
    Envelope getEnvelope() const;
 
    std::size_t getSize() const {
        return size();
    }
 
    size_t size() const
    {
        assert(stride() == 2 || stride() == 3 || stride() == 4);
        switch(stride()) {
            case 2: return m_vect.size() / 2;
            case 4: return m_vect.size() / 4;
            default : return m_vect.size() / 3;
        }
    }
 
    bool isEmpty() const {
        return m_vect.empty();
    }
 
    bool isRing() const;
 
    std::size_t getDimension() const;
 
    bool hasZ() const {
        return m_hasdim ? m_hasz : (m_vect.empty() || !std::isnan(m_vect[2]));
    }
 
    bool hasM() const {
        return m_hasm;
    }
 
    bool hasRepeatedPoints() const;
 
    bool hasRepeatedOrInvalidPoints() const;
 
    CoordinateType getCoordinateType() const {
        switch(stride()) {
            case 4: return CoordinateType::XYZM;
            case 2: return CoordinateType::XY;
            default: return hasM() ? CoordinateType::XYM : CoordinateType::XYZ;
        }
    }
 
 
    template<typename T=Coordinate>
    const T& getAt(std::size_t i) const {
        static_assert(std::is_base_of<CoordinateXY, T>::value, "Must be a Coordinate class");
        assert(sizeof(T) <= sizeof(double) * stride());
        assert(i*stride() < m_vect.size());
        const T* orig = reinterpret_cast<const T*>(&m_vect[i*stride()]);
        return *orig;
    }
 
    template<typename T=Coordinate>
    T& getAt(std::size_t i) {
        static_assert(std::is_base_of<CoordinateXY, T>::value, "Must be a Coordinate class");
        assert(sizeof(T) <= sizeof(double) * stride());
        assert(i*stride() < m_vect.size());
        T* orig = reinterpret_cast<T*>(&m_vect[i*stride()]);
        return *orig;
    }
 
    template<typename T>
    void getAt(std::size_t i, T& c) const {
        switch(getCoordinateType()) {
            case CoordinateType::XY: c = getAt<CoordinateXY>(i); break;
            case CoordinateType::XYZ: c = getAt<Coordinate>(i); break;
            case CoordinateType::XYZM: c = getAt<CoordinateXYZM>(i); break;
            case CoordinateType::XYM: c = getAt<CoordinateXYM>(i); break;
            default: getAt<Coordinate>(i);
        }
    }
 
    void getAt(std::size_t i, CoordinateXY& c) const {
        c = getAt<CoordinateXY>(i);
    }
 
    // TODO: change to return CoordinateXY
    const Coordinate& operator[](std::size_t i) const
    {
        return getAt(i);
    }
 
    // TODO: change to return CoordinateXY
    Coordinate&
    operator[](std::size_t i)
    {
        return getAt(i);
    }
 
    double getOrdinate(std::size_t index, std::size_t ordinateIndex) const;
 
    double getX(std::size_t index) const
    {
        return m_vect[index * stride()];
    }
 
    double getY(std::size_t index) const
    {
        return m_vect[index * stride() + 1];
    }
 
    double getZ(std::size_t index) const
    {
        return getOrdinate(index, Z);
    }
 
    double getM(std::size_t index) const
    {
        return getOrdinate(index, M);
    }
 
    void setX(std::size_t index, double x)
    {
        m_vect[index * stride()] = x;
    }
 
    void setY(std::size_t index, double y)
    {
        m_vect[index * stride() + 1] = y;
    }
 
    void setZ(std::size_t index, double z)
    {
        if (hasZ())
            setOrdinate(index, Z, z);
    }
 
    void setM(std::size_t index, double m)
    {
        if (hasM())
            setOrdinate(index, M, m);
    }
 
    template<typename T=Coordinate>
    const T& back() const
    {
        return getAt<T>(size() - 1);
    }
 
    template<typename T=Coordinate>
    T& back()
    {
        return getAt<T>(size() - 1);
    }
 
    template<typename T=Coordinate>
    const T& front() const
    {
        return *(reinterpret_cast<const T*>(m_vect.data()));
    }
 
    template<typename T=Coordinate>
    T& front()
    {
        return *(reinterpret_cast<T*>(m_vect.data()));
    }
 
    void toVector(std::vector<Coordinate>& coords) const;
 
    void toVector(std::vector<CoordinateXY>& coords) const;
 
 
 
    template<typename T>
    void setAt(const T& c, std::size_t pos) {
        switch(getCoordinateType()) {
            case CoordinateType::XY: setAtImpl<CoordinateXY>(c, pos); break;
            case CoordinateType::XYZ: setAtImpl<Coordinate>(c, pos); break;
            case CoordinateType::XYZM: setAtImpl<CoordinateXYZM>(c, pos); break;
            case CoordinateType::XYM: setAtImpl<CoordinateXYM>(c, pos); break;
            default: setAtImpl<Coordinate>(c, pos);
        }
    }
 
    void setOrdinate(std::size_t index, std::size_t ordinateIndex, double value);
 
    void setPoints(const std::vector<Coordinate>& v);
 
    void setPoints(const std::vector<CoordinateXY>& v);
 
 
    template<typename T=Coordinate>
    void add(const T& c) {
        add(c, size());
    }
 
    template<typename T>
    void add(const T& c, bool allowRepeated)
    {
        if(!allowRepeated && !isEmpty()) {
            const CoordinateXY& last = back<CoordinateXY>();
            if(last.equals2D(c)) {
                return;
            }
        }
 
        add(c);
    }
 
    template<typename T>
    void add(const T& c, std::size_t pos)
    {
        static_assert(std::is_base_of<CoordinateXY, T>::value, "Must be a Coordinate class");
 
        // c may be a reference inside m_vect, so we make sure it will not
        // grow before adding it
        if (m_vect.size() + stride() <= m_vect.capacity()) {
            make_space(pos, 1);
            setAt(c, static_cast<std::size_t>(pos));
        } else {
            T tmp{c};
            make_space(pos, 1);
            setAt(tmp, static_cast<std::size_t>(pos));
        }
    }
 
    template<typename T>
    void add(std::size_t i, const T& coord, bool allowRepeated)
    {
        // don't add duplicate coordinates
        if(! allowRepeated) {
            std::size_t sz = size();
            if(sz > 0) {
                if(i > 0) {
                    const CoordinateXY& prev = getAt<CoordinateXY>(i - 1);
                    if(prev.equals2D(coord)) {
                        return;
                    }
                }
                if(i < sz) {
                    const CoordinateXY& next = getAt<CoordinateXY>(i);
                    if(next.equals2D(coord)) {
                        return;
                    }
                }
            }
        }
 
        add(coord, i);
    }
 
    void add(double x, double y) {
        CoordinateXY c(x, y);
        add(c);
    }
 
    void add(const CoordinateSequence& cs);
 
    void add(const CoordinateSequence& cs, bool allowRepeated);
 
    void add(const CoordinateSequence& cl, bool allowRepeated, bool forwardDirection);
 
    void add(const CoordinateSequence& cs, std::size_t from, std::size_t to);
 
    void add(const CoordinateSequence& cs, std::size_t from, std::size_t to, bool allowRepeated);
 
    template<typename T, typename... Args>
    void add(T begin, T end, Args... args) {
        for (auto it = begin; it != end; ++it) {
            add(*it, args...);
        }
    }
 
    template<typename T>
    void add(std::size_t i, T from, T to) {
        auto npts = static_cast<std::size_t>(std::distance(from, to));
        make_space(i, npts);
 
        for (auto it = from; it != to; ++it) {
            setAt(*it, i);
            i++;
        }
    }
 
 
    void clear() {
        m_vect.clear();
    }
 
    void reserve(std::size_t capacity) {
        m_vect.reserve(capacity * stride());
    }
 
    void resize(std::size_t capacity) {
        m_vect.resize(capacity * stride());
    }
 
    void pop_back();
 
    std::string toString() const;
 
    const CoordinateXY* minCoordinate() const;
 
    static CoordinateSequence* atLeastNCoordinatesOrNothing(std::size_t n,
            CoordinateSequence* c);
 
    //
    static std::size_t indexOf(const CoordinateXY* coordinate,
                               const CoordinateSequence* cl);
 
    static bool equals(const CoordinateSequence* cl1,
                       const CoordinateSequence* cl2);
 
    bool equalsIdentical(const CoordinateSequence& other) const;
 
    static void scroll(CoordinateSequence* cl, const CoordinateXY* firstCoordinate);
 
    static int increasingDirection(const CoordinateSequence& pts);
 
    void reverse();
 
    void sort();
 
    void swap(std::size_t i, std::size_t j);
 
    void expandEnvelope(Envelope& env) const;
 
    void closeRing(bool allowRepeated = false);
 
 
    template<typename Filter>
    void apply_rw(const Filter* filter) {
        switch(getCoordinateType()) {
            case CoordinateType::XY:
                for (auto& c : items<CoordinateXY>()) {
                    if (filter->isDone()) break;
                    filter->filter_rw(&c);
                }
                break;
            case CoordinateType::XYZ:
                for (auto& c : items<Coordinate>()) {
                    if (filter->isDone()) break;
                    filter->filter_rw(&c);
                }
                break;
            case CoordinateType::XYM:
                for (auto& c : items<CoordinateXYM>()) {
                    if (filter->isDone()) break;
                    filter->filter_rw(&c);
                }
                break;
            case CoordinateType::XYZM:
                for (auto& c : items<CoordinateXYZM>()) {
                    if (filter->isDone()) break;
                    filter->filter_rw(&c);
                }
                break;
        }
        m_hasdim = m_hasz = false; // re-check (see http://trac.osgeo.org/geos/ticket/435)
    }
 
    template<typename Filter>
    void apply_ro(Filter* filter) const {
        switch(getCoordinateType()) {
            case CoordinateType::XY:
                for (const auto& c : items<CoordinateXY>()) {
                    if (filter->isDone()) break;
                    filter->filter_ro(&c);
                }
                break;
            case CoordinateType::XYZ:
                for (const auto& c : items<Coordinate>()) {
                    if (filter->isDone()) break;
                    filter->filter_ro(&c);
                }
                break;
            case CoordinateType::XYM:
                for (const auto& c : items<CoordinateXYM>()) {
                    if (filter->isDone()) break;
                    filter->filter_ro(&c);
                }
                break;
            case CoordinateType::XYZM:
                for (const auto& c : items<CoordinateXYZM>()) {
                    if (filter->isDone()) break;
                    filter->filter_ro(&c);
                }
                break;
        }
    }
 
    template<typename F>
    auto applyAt(size_t i, F&& fun) const {
        switch(getCoordinateType()) {
            case CoordinateType::XYZ:   return fun(getAt<Coordinate>(i));
            case CoordinateType::XYM:   return fun(getAt<CoordinateXYM>(i));
            case CoordinateType::XYZM:  return fun(getAt<CoordinateXYZM>(i));
            default:                    return fun(getAt<CoordinateXY>(i));
        }
    }
 
    template<typename F>
    void forEach(F&& fun) const {
        switch(getCoordinateType()) {
            case CoordinateType::XY:    for (const auto& c : items<CoordinateXY>())   { fun(c); } break;
            case CoordinateType::XYZ:   for (const auto& c : items<Coordinate>())     { fun(c); } break;
            case CoordinateType::XYM:   for (const auto& c : items<CoordinateXYM>())  { fun(c); } break;
            case CoordinateType::XYZM:  for (const auto& c : items<CoordinateXYZM>()) { fun(c); } break;
        }
    }
 
    template<typename T, typename F>
    void forEach(F&& fun) const
    {
        for (std::size_t i = 0; i < size(); i++) {
            fun(getAt<T>(i));
        }
    }
 
    template<typename F>
    void forEachSegment(F&& fun) const {
        switch(getCoordinateType()) {
            case CoordinateType::XY:    forEachSegmentImpl<CoordinateXY, F>(std::move(fun)); break;
            case CoordinateType::XYZ:   forEachSegmentImpl<Coordinate, F>(std::move(fun)); break;
            case CoordinateType::XYM:   forEachSegmentImpl<CoordinateXYM, F>(std::move(fun)); break;
            case CoordinateType::XYZM:  forEachSegmentImpl<CoordinateXYZM, F>(std::move(fun)); break;
        }
    }
 
    template<typename T, typename F>
    void forEach(std::size_t from, std::size_t to, F&& fun) const
    {
        for (std::size_t i = from; i <= to; i++) {
            fun(getAt<T>(i));
        }
    }
 
    template<typename T>
    class Coordinates {
    public:
        using SequenceType = typename std::conditional<std::is_const<T>::value, const CoordinateSequence, CoordinateSequence>::type;
 
        explicit Coordinates(SequenceType* seq) : m_seq(seq) {}
 
        CoordinateSequenceIterator<SequenceType, T> begin() {
            return {m_seq};
        }
 
        CoordinateSequenceIterator<SequenceType, T> end() {
            return {m_seq, m_seq->getSize()};
        }
 
        CoordinateSequenceIterator<const SequenceType, typename std::add_const<T>::type>
        begin() const {
            return CoordinateSequenceIterator<const SequenceType, typename std::add_const<T>::type>{m_seq};
        }
 
        CoordinateSequenceIterator<const SequenceType, typename std::add_const<T>::type>
        end() const {
            return CoordinateSequenceIterator<const SequenceType, typename std::add_const<T>::type>{m_seq, m_seq->getSize()};
        }
 
        CoordinateSequenceIterator<const SequenceType, typename std::add_const<T>::type>
        cbegin() const {
            return CoordinateSequenceIterator<const SequenceType, typename std::add_const<T>::type>{m_seq};
        }
 
        CoordinateSequenceIterator<const SequenceType, typename std::add_const<T>::type>
        cend() const {
            return CoordinateSequenceIterator<const SequenceType, typename std::add_const<T>::type>{m_seq, m_seq->getSize()};
        }
 
    private:
        SequenceType* m_seq;
    };
 
    template<typename T>
    Coordinates<typename std::add_const<T>::type> items() const {
        return Coordinates<typename std::add_const<T>::type>(this);
    }
 
    template<typename T>
    Coordinates<T> items() {
        return Coordinates<T>(this);
    }
 
 
 
    double* data() {
        return m_vect.data();
    }
 
    const double* data() const {
        return m_vect.data();
    }
 
private:
    std::vector<double> m_vect; // Vector to store values
 
    uint8_t m_stride;           // Stride of stored values, corresponding to underlying type
 
    mutable bool m_hasdim;      // Has the dimension of this sequence been determined? Or was it created with no
                                // explicit dimensionality, and we're waiting for getDimension() to be called
                                // after some coordinates have been added?
    mutable bool m_hasz;
    bool m_hasm;
 
    void initialize();
 
    template<typename T1, typename T2>
    void setAtImpl(const T2& c, std::size_t pos) {
        auto& orig = getAt<T1>(pos);
        orig = c;
    }
 
    template<typename CoordType, typename F>
    void forEachSegmentImpl(F&& fun) const {
        auto p0it = items<CoordType>().cbegin();
        const auto end = items<CoordType>().end();
 
        if (p0it == end) {
            return;
        }
 
        auto p1it = std::next(p0it);
 
        while (p1it != end) {
            fun(*p0it, *p1it);
            ++p0it;
            ++p1it;
        }
    }
 
    void make_space(std::size_t pos, std::size_t n) {
        m_vect.insert(std::next(m_vect.begin(), static_cast<std::ptrdiff_t>(pos * stride())),
                      m_stride * n,
                      DoubleNotANumber);
    }
 
    std::uint8_t stride() const {
        return m_stride;
    }
 
};
 
GEOS_DLL std::ostream& operator<< (std::ostream& os, const CoordinateSequence& cs);
 
GEOS_DLL bool operator== (const CoordinateSequence& s1, const CoordinateSequence& s2);
 
GEOS_DLL bool operator!= (const CoordinateSequence& s1, const CoordinateSequence& s2);
 
} // namespace geos::geom
} // namespace geos