Grid.h
/**********************************************************************
*
* GEOS - Geometry Engine Open Source
* http://geos.osgeo.org
*
* Copyright (C) 2018-2025 ISciences, LLC
*
* 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 <numeric>
#include <stdexcept>
#include <vector>
#include <geos/geom/Envelope.h>
constexpr double DEFAULT_GRID_COMPAT_TOL = 1e-6;
namespace geos::operation::grid {
struct infinite_extent
{
static constexpr size_t padding = 1;
};
struct bounded_extent
{
static constexpr size_t padding = 0;
};
template<typename extent_tag>
class Grid
{
public:
Grid(const geom::Envelope& extent, double dx, double dy)
: m_extent{ extent }
, m_domain{}
, m_dx{ dx }
, m_dy{ dy }
, m_num_rows{ 2 * extent_tag::padding + (extent.getMaxY() > extent.getMinY() ? static_cast<size_t>(std::round(extent.getHeight() / dy)) : 0) }
, m_num_cols{ 2 * extent_tag::padding + (extent.getMaxX() > extent.getMinX() ? static_cast<size_t>(std::round(extent.getWidth() / dx)) : 0) } {
static_assert(std::is_same_v<extent_tag, bounded_extent>);
}
Grid(const geom::Envelope& extent, double dx, double dy, const geom::Envelope& domain)
: m_extent{ extent }
, m_domain{ domain }
, m_dx{ dx }
, m_dy{ dy }
, m_num_rows{ 2 * extent_tag::padding + (extent.getMaxY() > extent.getMinY() ? static_cast<size_t>(std::round(extent.getHeight() / dy)) : 0) }
, m_num_cols{ 2 * extent_tag::padding + (extent.getMaxX() > extent.getMinX() ? static_cast<size_t>(std::round(extent.getWidth() / dx)) : 0) }
{
static_assert(std::is_same_v<extent_tag, infinite_extent>);
}
static Grid makeEmpty()
{
if constexpr (std::is_same_v<extent_tag, bounded_extent>) {
return Grid({ 0, 0, 0, 0 }, 0, 0);
} else {
return Grid({ 0, 0, 0, 0 }, 0, 0, geom::Envelope());
}
}
size_t getColumn(double x) const
{
if (extent_tag::padding) {
if (x < m_extent.getMinX())
return 0;
if (x > m_extent.getMaxX())
return m_num_cols - 1;
if (x == m_extent.getMaxX()) // special case, returning the cell for which xmax is the right
return m_num_cols - 2;
} else {
if (x < m_extent.getMinX() || x > m_extent.getMaxX())
throw std::out_of_range("x");
if (x == m_extent.getMaxX())
return m_num_cols - 1;
}
// Since we've already range-checked our x value, make sure that
// the computed column index is no greater than the column index
// associated with xmax.
return std::min(
extent_tag::padding + static_cast<size_t>(std::floor((x - m_extent.getMinX()) / m_dx)),
getColumn(m_extent.getMaxX()));
}
size_t getRow(double y) const
{
if (extent_tag::padding) {
if (y > m_extent.getMaxY())
return 0;
if (y < m_extent.getMinY())
return m_num_rows - 1;
if (y == m_extent.getMinY()) // special case, returning the cell for which ymin is the bottom
return m_num_rows - 2;
} else {
if (y < m_extent.getMinY() || y > m_extent.getMaxY())
throw std::out_of_range("y");
if (y == m_extent.getMinY())
return m_num_rows - 1;
}
// Since we've already range-checked our y value, make sure that
// the computed row index is no greater than the column index
// associated with ymin.
return std::min(
extent_tag::padding + static_cast<size_t>(std::floor((m_extent.getMaxY() - y) / m_dy)),
getRow(m_extent.getMinY()));
}
std::size_t getCell(double x, double y) const
{
return getRow(y) * getNumCols() + getColumn(x);
}
bool isEmpty() const { return m_num_rows <= 2 * extent_tag::padding && m_num_cols <= 2 * extent_tag::padding; }
size_t getNumRows() const { return m_num_rows; }
size_t getNumCols() const { return m_num_cols; }
size_t getSize() const { return getNumRows() * getNumCols(); }
double xmin() const { return m_extent.getMinX(); }
double xmax() const { return m_extent.getMaxX(); }
double ymin() const { return m_extent.getMinY(); }
double ymax() const { return m_extent.getMaxY(); }
double dx() const { return m_dx; }
double dy() const { return m_dy; }
const geom::Envelope& getExtent() const { return m_extent; }
size_t getRowOffset(const Grid& other) const { return static_cast<size_t>(std::round(std::abs(other.m_extent.getMaxY() - m_extent.getMaxY()) / m_dy)); }
size_t getColOffset(const Grid& other) const { return static_cast<size_t>(std::round(std::abs(m_extent.getMinX() - other.m_extent.getMinX()) / m_dx)); }
double getColX(size_t col) const { return m_extent.getMinX() + (static_cast<double>(col - extent_tag::padding) + 0.5) * m_dx; }
double getRowY(size_t row) const { return m_extent.getMaxY() - (static_cast<double>(row - extent_tag::padding) + 0.5) * m_dy; }
Grid<extent_tag> shrinkToFit(const geom::Envelope& b) const
{
if (b.getArea() == 0) {
return makeEmpty();
}
if (b.getMinX() < m_extent.getMinX() || b.getMinY() < m_extent.getMinY() || b.getMaxX() > m_extent.getMaxX() || b.getMaxY() > m_extent.getMaxY()) {
throw std::range_error("Cannot shrink extent to bounds larger than original.");
}
size_t col0 = getColumn(b.getMinX());
size_t row1 = getRow(b.getMaxY());
// Shrink xmin and ymax to fit the upper-left corner of the supplied extent
double snapped_xmin = m_extent.getMinX() + static_cast<double>(col0 - extent_tag::padding) * m_dx;
double snapped_ymax = m_extent.getMaxY() - static_cast<double>(row1 - extent_tag::padding) * m_dy;
// Make sure x0 and y1 are within the reduced extent. Because of
// floating point round-off errors, this is not always the case.
if (b.getMinX() < snapped_xmin) {
snapped_xmin -= m_dx;
col0--;
}
if (b.getMaxY() > snapped_ymax) {
snapped_ymax += m_dy;
row1--;
}
size_t col1 = getColumn(b.getMaxX());
size_t row0 = getRow(b.getMinY());
size_t num_rows = 1 + (row0 - row1);
size_t num_cols = 1 + (col1 - col0);
// If xmax or ymin falls cleanly on a cell boundary, we don't
// need as many rows or columns as we otherwise would, because
// we assume that the rightmost cell of the grid is a closed
// interval in x, and the lowermost cell of the grid is a
// closed interval in y.
if (num_rows > 2 && (snapped_ymax - static_cast<double>(num_rows - 1) * m_dy <= b.getMinY())) {
num_rows--;
}
if (num_cols > 2 && (snapped_xmin + static_cast<double>(num_cols - 1) * m_dx >= b.getMaxX())) {
num_cols--;
}
// Perform offsets relative to the new xmin, ymax origin
// points. If this is not done, then floating point roundoff
// error can cause progressive shrink() calls with the same
// inputs to produce different results.
geom::Envelope reduced_box = {
snapped_xmin,
std::max(snapped_xmin + static_cast<double>(num_cols) * m_dx, b.getMaxX()),
std::min(snapped_ymax - static_cast<double>(num_rows) * m_dy, b.getMinY()),
snapped_ymax
};
// Fudge computed xmax and ymin, if needed, to prevent extent
// from growing during a shrink operation.
if (reduced_box.getMaxX() > m_extent.getMaxX()) {
if (std::round(reduced_box.getWidth() / m_dx) ==
std::round((m_extent.getMaxX() - reduced_box.getMinX()) / m_dx)) {
reduced_box = geom::Envelope(reduced_box.getMinX(), m_extent.getMaxX(), reduced_box.getMinY(), reduced_box.getMaxY());
} else {
throw std::runtime_error("Shrink operation failed.");
}
}
if (reduced_box.getMinY() < m_extent.getMinY()) {
if (std::round(reduced_box.getHeight() / m_dy) ==
std::round((reduced_box.getMaxY() - m_extent.getMinY()) / m_dy)) {
reduced_box = geom::Envelope(reduced_box.getMinX(), reduced_box.getMaxX(), m_extent.getMinY(), reduced_box.getMaxY());
} else {
throw std::runtime_error("Shrink operation failed.");
}
}
if constexpr (std::is_same_v<extent_tag, bounded_extent>) {
Grid<extent_tag> reduced{ reduced_box, m_dx, m_dy };
if (b.getMinX() < reduced.xmin() || b.getMinY() < reduced.ymin() || b.getMaxX() > reduced.xmax() || b.getMaxY() > reduced.ymax()) {
throw std::runtime_error("Shrink operation failed.");
}
return reduced;
} else {
Grid<extent_tag> reduced{ reduced_box, m_dx, m_dy , m_domain };
if (b.getMinX() < reduced.xmin() || b.getMinY() < reduced.ymin() || b.getMaxX() > reduced.xmax() || b.getMaxY() > reduced.ymax()) {
throw std::runtime_error("Shrink operation failed.");
}
return reduced;
}
}
bool operator==(const Grid<extent_tag>& b) const
{
return m_extent == b.m_extent &&
m_dx == b.m_dx &&
m_dy == b.m_dy;
}
bool operator!=(const Grid<extent_tag>& b) const
{
return !(*this == b);
}
geom::Envelope getCellEnvelope(size_t row, size_t col) const;
private:
geom::Envelope m_extent;
geom::Envelope m_domain;
double m_dx;
double m_dy;
size_t m_num_rows;
size_t m_num_cols;
};
Grid<infinite_extent>
make_infinite(const Grid<bounded_extent>& grid, const geom::Envelope& domain);
Grid<bounded_extent>
make_finite(const Grid<infinite_extent>& grid);
}
geos::util::round
double round(double val)
Definition math.h:36