DiscreteFrechetDistance.h

/**********************************************************************
 *
 * GEOS - Geometry Engine Open Source
 * http://geos.osgeo.org
 *
 * Copyright (c) 2021 Felix Obermaier
 * Copyright (c) 2025 Paul Ramsey
 *
 * 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/algorithm/distance/PointPairDistance.h>
#include <geos/geom/CoordinateSequenceFilter.h>
#include <geos/util/math.h>
 
#include <algorithm> // for lower_bound
#include <cstdint>
#include <unordered_map>
#include <array>
#include <vector>
#include <memory>
 
// Forward declarations
namespace geos {
namespace geom {
class Geometry;
class Coordinate;
class CoordinateSequence;
}
}
 
namespace geos {
namespace algorithm { // geos::algorithm
namespace distance {  // geos::algorithm::distance
 
class GEOS_DLL DiscreteFrechetDistance {
 
    using DistanceToPairMap = std::unordered_map<double, std::array<std::size_t, 2>>;
 
public:
 
    DiscreteFrechetDistance(const geom::Geometry& geom0, const geom::Geometry& geom1)
        : g0(geom0)
        , g1(geom1) {};
 
 
    class GEOS_DLL MatrixStorage {
 
        public:
 
            std::size_t m_numRows;
            std::size_t m_numCols;
            double m_defaultValue;
 
            MatrixStorage(std::size_t numRows, std::size_t numCols, double defaultValue)
                : m_numRows(numRows)
                , m_numCols(numCols)
                , m_defaultValue(defaultValue) {};
 
            virtual ~MatrixStorage() = 0;
 
            virtual double get(std::size_t i, std::size_t j) const = 0;
 
            virtual void set(std::size_t i, std::size_t j, double value) = 0;
 
            virtual bool isValueSet(std::size_t i, std::size_t j) const = 0;
    };
 
 
    class RectMatrix : public MatrixStorage {
 
        private:
            std::vector<double> m_matrix;
 
        public:
 
            RectMatrix(std::size_t numRows, std::size_t numCols, double defaultValue)
                : MatrixStorage(numRows, numCols, defaultValue)
            {
                m_matrix.resize(numRows * numCols, defaultValue);
            };
 
            double get(std::size_t i, std::size_t j) const override {
                return m_matrix[i * m_numCols + j];
            };
 
            void set(std::size_t i, std::size_t j, double value) override {
                m_matrix[i * m_numCols + j] = value;
            };
 
            bool isValueSet(std::size_t i, std::size_t j) const override {
                return get(i, j) != m_defaultValue;
            };
    };
 
 
    class CsrMatrix : public MatrixStorage {
 
        private:
            std::vector<double> m_v;
            std::vector<std::size_t> m_ri;
            std::vector<std::size_t> m_ci;
 
        public:
 
            CsrMatrix(std::size_t numRows, std::size_t numCols, double defaultValue, std::size_t expectedValues)
                : MatrixStorage(numRows, numCols, defaultValue)
            {
                m_v.resize(expectedValues);
                m_ci.resize(expectedValues);
                m_ri.resize(numRows + 1);
            };
 
            CsrMatrix(std::size_t numRows, std::size_t numCols, double defaultValue)
                : CsrMatrix(numRows, numCols, defaultValue, expectedValuesHeuristic(numRows, numCols))
                {};
 
            static std::size_t expectedValuesHeuristic(std::size_t numRows, std::size_t numCols) {
                std::size_t max = std::max(numRows, numCols);
                return max * max / 10;
            };
 
            std::pair<bool, std::size_t>
            cppBinarySearch(const std::vector<std::size_t>& vec, std::size_t fromIndex, std::size_t toIndex, std::size_t key) const {
                // Return not found on invalid input
                if (fromIndex > toIndex || toIndex > vec.size())
                    return {false, 0};
 
                // Define the iterators for the sub-range
                auto first_it = vec.begin() + static_cast<std::ptrdiff_t>(fromIndex);
                auto last_it = vec.begin() + static_cast<std::ptrdiff_t>(toIndex);
 
                // Perform the binary search using std::lower_bound
                auto it = std::lower_bound(first_it, last_it, key);
 
                // Calculate the index relative to the *beginning of the original vector*
                auto result_index = std::distance(vec.begin(), it);
 
                // Determine if the element was found and return the appropriate value
                if (it != last_it && *it == key) {
                    // Element found, return its index
                    return {true, result_index};
                } else {
                    // Element not found, return the insertion point
                    return {false, result_index};
                }
            };
 
            std::pair<bool, std::size_t> indexOf(std::size_t i, std::size_t j) const {
                std::size_t cLow = m_ri[i];
                std::size_t cHigh = m_ri[i+1];
                if (cHigh <= cLow) return {false, cLow};
                return cppBinarySearch(m_ci, cLow, cHigh, j);
            };
 
            double get(std::size_t i, std::size_t j) const override {
                // get the index in the vector
                std::pair<bool, std::size_t> vi = indexOf(i, j);
                // if the vector index is negative, return default value
                if (!vi.first)
                    return m_defaultValue;
 
                return m_v[vi.second];
            };
 
            void set(std::size_t i, std::size_t j, double value) override {
 
                // get the index in the vector
                std::pair<bool, std::size_t> vi = indexOf(i, j);
 
                // do we already have a value?
                if (!vi.first)
                {
                    // no, we don't, we need to ensure space!
                    ensureCapacity(m_ri[m_numRows] + 1);
 
                    // update row indices
                    for (std::size_t ii = i + 1; ii <= m_numRows; ii++)
                        m_ri[ii] += 1;
 
                    // move and update column indices, move values
                    std::size_t viv = vi.second;
                    for (std::size_t ii = m_ri[m_numRows]; ii > viv; ii--) {
                        m_ci[ii] = m_ci[ii - 1];
                        m_v[ii] = m_v[ii - 1];
                    }
 
                    // insert column index
                    m_ci[viv] = j;
                }
 
                // set the new value
                m_v[vi.second] = value;
            };
 
            bool isValueSet(std::size_t i, std::size_t j) const override {
                auto r = indexOf(i, j);
                return r.first;
            };
 
            void ensureCapacity(std::size_t required) {
                if (required < m_v.size())
                    return;
 
                std::size_t increment = std::max(m_numRows, m_numCols);
                m_v.resize(m_v.size() + increment);
                m_ci.resize(m_v.size() + increment);
            };
    };
 
 
    class HashMapMatrix : public MatrixStorage {
 
        private:
            std::unordered_map<std::size_t, double> m_matrix;
 
        public:
 
        HashMapMatrix(std::size_t numRows, std::size_t numCols, double defaultValue)
            : MatrixStorage(numRows, numCols, defaultValue)
            {};
 
        double get(std::size_t i, std::size_t j) const override {
            static constexpr std::size_t halfsz = sizeof(std::size_t) * 4;
            std::size_t key = i << halfsz | j;
            auto it_found = m_matrix.find(key);
            if (it_found != m_matrix.end()) {
                return (*it_found).second;
            }
            else {
                return m_defaultValue;
            }
        };
 
        void set(std::size_t i, std::size_t j, double value) override {
            static constexpr std::size_t halfsz = sizeof(std::size_t) * 4;
            std::size_t key = i << halfsz | j;
            m_matrix[key] = value;
        };
 
        bool isValueSet(std::size_t i, std::size_t j) const override {
            static constexpr std::size_t halfsz = sizeof(std::size_t) * 4;
            std::size_t key = i << halfsz | j;
            auto it_found = m_matrix.find(key);
            return it_found != m_matrix.end();
        };
    };
 
 
    static double distance(const geom::Geometry& geom0, const geom::Geometry& geom1);
 
    static double distance(const geom::Geometry& geom0, const geom::Geometry& geom1, double densityFrac);
 
    std::array<geom::CoordinateXY, 2> getCoordinates();
 
    void setDensifyFraction(double dFrac);
 
 
private:
 
    // Members
    const geom::Geometry& g0;
    const geom::Geometry& g1;
    std::unique_ptr<PointPairDistance> ptDist;
    double densifyFraction = -1.0;
 
    /* private */
    double distance();
 
    /*
     * Utility method to ape Java behaviour
     */
    void putIfAbsent(
        DistanceToPairMap& distanceToPair,
        double key, std::array<std::size_t, 2> val);
 
    static std::unique_ptr<MatrixStorage> createMatrixStorage(
        std::size_t rows, std::size_t cols);
 
    static std::unique_ptr<PointPairDistance> computeFrechet(
        const geom::CoordinateSequence& coords0,
        const geom::CoordinateSequence& coords1,
        std::vector<std::size_t>& diagonal,
        MatrixStorage& distances,
        DistanceToPairMap& distanceToPair);
 
    /* private */
    static double getMinDistanceAtCorner(
        MatrixStorage& matrix, std::size_t i, std::size_t j);
 
    void computeCoordinateDistances(
        const geom::CoordinateSequence& coords0, const geom::CoordinateSequence& coords1,
        std::vector<std::size_t>& diagonal,
        MatrixStorage& distances, DistanceToPairMap& distanceToPair);
 
    static std::vector<std::size_t> bresenhamDiagonal(
        std::size_t numCols, std::size_t numRows);
 
    static std::unique_ptr<geom::CoordinateSequence> getDensifiedCoordinates(
        const geom::Geometry& geom, double densifyFrac);
 
 
    class DensifiedCoordinatesFilter : public geom::CoordinateSequenceFilter {
 
        public:
 
            DensifiedCoordinatesFilter(double fraction, geom::CoordinateSequence& coords)
                : m_numSubSegs(static_cast<uint32_t>(util::round(1.0 / fraction)))
                , m_coords(coords)
            {};
 
            void filter_ro(
                const geom::CoordinateSequence& seq,
                std::size_t index) override;
 
            bool isGeometryChanged() const override {
                return false;
            }
 
            bool isDone() const override {
                return false;
            }
 
        private:
 
            uint32_t m_numSubSegs;
            geom::CoordinateSequence& m_coords;
 
    };
 
 
 
}; // DiscreteFrechetDistance
 
} // geos::algorithm::distance
} // geos::algorithm
} // geos