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 <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
geos::algorithm::distance::DiscreteFrechetDistance::CsrMatrix::expectedValuesHeuristic
static std::size_t expectedValuesHeuristic(std::size_t numRows, std::size_t numCols)
Definition DiscreteFrechetDistance.h:213
geos::algorithm::distance::DiscreteFrechetDistance::CsrMatrix::ensureCapacity
void ensureCapacity(std::size_t required)
Definition DiscreteFrechetDistance.h:313
geos::algorithm::distance::DiscreteFrechetDistance::CsrMatrix::set
void set(std::size_t i, std::size_t j, double value) override
Definition DiscreteFrechetDistance.h:274
geos::algorithm::distance::DiscreteFrechetDistance::CsrMatrix::isValueSet
bool isValueSet(std::size_t i, std::size_t j) const override
Definition DiscreteFrechetDistance.h:304
geos::algorithm::distance::DiscreteFrechetDistance::CsrMatrix::get
double get(std::size_t i, std::size_t j) const override
Definition DiscreteFrechetDistance.h:264
geos::algorithm::distance::DiscreteFrechetDistance::CsrMatrix::cppBinarySearch
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
Definition DiscreteFrechetDistance.h:232
geos::algorithm::distance::DiscreteFrechetDistance::HashMapMatrix::get
double get(std::size_t i, std::size_t j) const override
Definition DiscreteFrechetDistance.h:344
geos::algorithm::distance::DiscreteFrechetDistance::HashMapMatrix::set
void set(std::size_t i, std::size_t j, double value) override
Definition DiscreteFrechetDistance.h:356
geos::algorithm::distance::DiscreteFrechetDistance::HashMapMatrix::HashMapMatrix
HashMapMatrix(std::size_t numRows, std::size_t numCols, double defaultValue)
Definition DiscreteFrechetDistance.h:340
geos::algorithm::distance::DiscreteFrechetDistance::HashMapMatrix::isValueSet
bool isValueSet(std::size_t i, std::size_t j) const override
Definition DiscreteFrechetDistance.h:362
geos::algorithm::distance::DiscreteFrechetDistance::MatrixStorage
Definition DiscreteFrechetDistance.h:94
geos::algorithm::distance::DiscreteFrechetDistance::MatrixStorage::set
virtual void set(std::size_t i, std::size_t j, double value)=0
geos::algorithm::distance::DiscreteFrechetDistance::MatrixStorage::get
virtual double get(std::size_t i, std::size_t j) const =0
geos::algorithm::distance::DiscreteFrechetDistance::MatrixStorage::MatrixStorage
MatrixStorage(std::size_t numRows, std::size_t numCols, double defaultValue)
Definition DiscreteFrechetDistance.h:108
geos::algorithm::distance::DiscreteFrechetDistance::MatrixStorage::isValueSet
virtual bool isValueSet(std::size_t i, std::size_t j) const =0
geos::algorithm::distance::DiscreteFrechetDistance::RectMatrix::set
void set(std::size_t i, std::size_t j, double value) override
Definition DiscreteFrechetDistance.h:170
geos::algorithm::distance::DiscreteFrechetDistance::RectMatrix::get
double get(std::size_t i, std::size_t j) const override
Definition DiscreteFrechetDistance.h:166
geos::algorithm::distance::DiscreteFrechetDistance::RectMatrix::isValueSet
bool isValueSet(std::size_t i, std::size_t j) const override
Definition DiscreteFrechetDistance.h:174
geos::algorithm::distance::DiscreteFrechetDistance::RectMatrix::RectMatrix
RectMatrix(std::size_t numRows, std::size_t numCols, double defaultValue)
Definition DiscreteFrechetDistance.h:160
geos::algorithm::distance::DiscreteFrechetDistance::distance
static double distance(const geom::Geometry &geom0, const geom::Geometry &geom1, double densityFrac)
geos::algorithm::distance::DiscreteFrechetDistance::distance
static double distance(const geom::Geometry &geom0, const geom::Geometry &geom1)
geos::algorithm::distance::DiscreteFrechetDistance::DiscreteFrechetDistance
DiscreteFrechetDistance(const geom::Geometry &geom0, const geom::Geometry &geom1)
Definition DiscreteFrechetDistance.h:86
geos::algorithm::distance::DiscreteFrechetDistance::getCoordinates
std::array< geom::CoordinateXY, 2 > getCoordinates()
geos::geom::CoordinateSequenceFilter
Interface for classes which provide operations that can be applied to the coordinates in a Coordinate...
Definition CoordinateSequenceFilter.h:56
geos::geom::CoordinateSequence
The internal representation of a list of coordinates inside a Geometry.
Definition CoordinateSequence.h:56
geos::geom::Coordinate
Coordinate is the lightweight class used to store coordinates.
Definition Coordinate.h:217
geos::geom::Geometry
Basic implementation of Geometry, constructed and destructed by GeometryFactory.
Definition Geometry.h:196
geos::algorithm::distance
Classes to compute distance metrics between geometries.
Definition DiscreteFrechetDistance.h:40
geos::algorithm
Contains classes and interfaces implementing fundamental computational geometry algorithms.
Definition Angle.h:32
geos::geom
Definition Angle.h:26
geos::index
Provides classes for various kinds of spatial indexes.
Definition CleanCoverage.h:39
geos::util
Utility classes for GEOS.
Definition namespaces.h:314
geos
Basic namespace for all GEOS functionalities.
Definition geos.h:38