RegularGridAdjacency2D.hpp

#pragma once
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <initializer_list>
#include <iterator>
#include <limits>
#include <numbers>
#include <span>
#include <stdexcept>
#include <utility>
#include <vector>
 
#include "Common.hpp"
#include "Contract.hpp"
 
namespace mmcfilters {
 
namespace detail {
class CommittedGridAccess;
}
 
struct GridOffset2D {
    int row = 0;
    int column = 0;
 
    bool operator==(const GridOffset2D&) const noexcept = default;
};
 
enum class RegularGridAdjacencyShape { EuclideanDisk, StructuringElement };
 
class RegularGridAdjacency2D {
  private:
    friend class detail::CommittedGridAccess;
    struct StructuringElementTag {};
 
    int numColumns;
    int numRows;
    double radius;
    double radius2;
    int n;
    RegularGridAdjacencyShape shape = RegularGridAdjacencyShape::EuclideanDisk;
 
    std::vector<int> offsetRow;
    std::vector<int> offsetColumn;
    std::vector<int> forwardOffsetIndices;
 
    static void requireDomainDimensions(int rows, int columns) {
        if (rows < 0 || columns < 0) {
            throw std::invalid_argument("RegularGridAdjacency2D grid dimensions must be non-negative.");
        }
        if (static_cast<std::int64_t>(rows) * columns > std::numeric_limits<PixelId>::max()) {
            throw std::overflow_error("RegularGridAdjacency2D grid dimensions exceed the supported pixel-indexed size.");
        }
    }
 
    static double checkedRadiusParameters(int rows, int columns, double radius) {
        requireDomainDimensions(rows, columns);
        const double maxSafeRadius = (std::sqrt(static_cast<double>(std::numeric_limits<int>::max())) - 1.0) / 2.0;
        if (!std::isfinite(radius) || radius < 1.0 || radius > maxSafeRadius) {
            throw std::invalid_argument("RegularGridAdjacency2D radius must be finite, at least 1.0, and representable by the integer stencil.");
        }
        return radius;
    }
 
    [[nodiscard]] bool containsOffset(int rowOffset, int columnOffset) const noexcept {
        for (int index = 0; index < n; ++index) {
            if (offsetRow[static_cast<std::size_t>(index)] == rowOffset && offsetColumn[static_cast<std::size_t>(index)] == columnOffset) {
                return true;
            }
        }
        return false;
    }
 
    void requireConnectedDomain() const {
        if (numRows == 0 || numColumns == 0) {
            throw std::invalid_argument("A structuring-element adjacency requires a nonempty connected grid domain.");
        }
        // Unit steps connect every row and column without a graph traversal.
        if ((numRows == 1 || containsOffset(1, 0)) && (numColumns == 1 || containsOffset(0, 1))) {
            return;
        }
 
        const auto numPixels = static_cast<std::size_t>(numRows) * static_cast<std::size_t>(numColumns);
        std::vector<std::uint8_t> visited(numPixels, 0);
        std::vector<PixelId> pending{0};
        visited[0] = 1;
        std::size_t reached = 1;
        while (!pending.empty()) {
            const PixelId pixel = pending.back();
            pending.pop_back();
            for (PixelId neighbor : NeighborIndexRange(*this, pixel / numColumns, pixel % numColumns)) {
                if (!visited[static_cast<std::size_t>(neighbor)]) {
                    visited[static_cast<std::size_t>(neighbor)] = 1;
                    ++reached;
                    pending.push_back(neighbor);
                }
            }
        }
        if (reached != numPixels) {
            throw std::invalid_argument("A structuring-element adjacency must induce one connected grid-domain graph.");
        }
    }
 
    [[nodiscard]] bool matchesOffsets(std::initializer_list<GridOffset2D> expected) const noexcept {
        if (static_cast<std::size_t>(n) != expected.size()) {
            return false;
        }
        for (const GridOffset2D offset : expected) {
            if (!containsOffset(offset.row, offset.column)) {
                return false;
            }
        }
        return true;
    }
 
    void buildForwardOffsetIndices() {
        forwardOffsetIndices.clear();
        forwardOffsetIndices.reserve(static_cast<std::size_t>(n / 2));
        for (int index = 1; index < n; ++index) {
            const int dx = offsetColumn[static_cast<std::size_t>(index)];
            const int dy = offsetRow[static_cast<std::size_t>(index)];
            if (dy > 0 || (dy == 0 && dx > 0)) {
                forwardOffsetIndices.push_back(index);
            }
        }
    }
 
    RegularGridAdjacency2D(int rows, int columns, std::vector<GridOffset2D> offsets, StructuringElementTag)
        : numColumns(columns), numRows(rows), radius(0.0), radius2(0.0), n(0), shape(RegularGridAdjacencyShape::StructuringElement) {
        requireDomainDimensions(rows, columns);
        if (offsets.empty()) {
            throw std::invalid_argument("A structuring-element adjacency requires at least the origin offset.");
        }
 
        for (const GridOffset2D offset : offsets) {
            if (offset.row == std::numeric_limits<int>::min() || offset.column == std::numeric_limits<int>::min()) {
                throw std::invalid_argument("Structuring-element offsets must be safely negatable.");
            }
        }
 
        std::sort(offsets.begin(), offsets.end(), [](const GridOffset2D& lhs, const GridOffset2D& rhs) {
            if (lhs.row != rhs.row) {
                return lhs.row < rhs.row;
            }
            return lhs.column < rhs.column;
        });
        const auto duplicate = std::adjacent_find(offsets.begin(), offsets.end());
        if (duplicate != offsets.end()) {
            throw std::invalid_argument("A structuring-element adjacency cannot contain duplicate offsets.");
        }
 
        const auto origin = std::find(offsets.begin(), offsets.end(), GridOffset2D{0, 0});
        if (origin == offsets.end()) {
            throw std::invalid_argument("A structuring-element adjacency must contain the origin offset.");
        }
        for (const GridOffset2D offset : offsets) {
            if (!std::binary_search(offsets.begin(), offsets.end(), GridOffset2D{-offset.row, -offset.column},
                                    [](const GridOffset2D& lhs, const GridOffset2D& rhs) {
                                        if (lhs.row != rhs.row) {
                                            return lhs.row < rhs.row;
                                        }
                                        return lhs.column < rhs.column;
                                    })) {
                throw std::invalid_argument("An adjacency-inducing structuring element must be centrally symmetric.");
            }
        }
 
        std::vector<GridOffset2D> ordered;
        ordered.reserve(offsets.size());
        ordered.push_back({0, 0});
        offsets.erase(origin);
        std::sort(offsets.begin(), offsets.end(), [](const GridOffset2D& lhs, const GridOffset2D& rhs) {
            auto angle = [](GridOffset2D offset) {
                double value = std::atan2(-static_cast<double>(offset.row), -static_cast<double>(offset.column));
                if (value < 0.0) {
                    value += 2.0 * std::numbers::pi;
                }
                return value;
            };
            const double lhsAngle = angle(lhs);
            const double rhsAngle = angle(rhs);
            if (lhsAngle != rhsAngle) {
                return lhsAngle < rhsAngle;
            }
            const std::int64_t lhsRadius = static_cast<std::int64_t>(lhs.row) * lhs.row + static_cast<std::int64_t>(lhs.column) * lhs.column;
            const std::int64_t rhsRadius = static_cast<std::int64_t>(rhs.row) * rhs.row + static_cast<std::int64_t>(rhs.column) * rhs.column;
            return lhsRadius < rhsRadius;
        });
        ordered.insert(ordered.end(), offsets.begin(), offsets.end());
 
        if (ordered.size() > static_cast<std::size_t>(std::numeric_limits<int>::max())) {
            throw std::length_error("Structuring-element stencil exceeds the supported offset count.");
        }
        n = static_cast<int>(ordered.size());
        offsetRow.reserve(ordered.size());
        offsetColumn.reserve(ordered.size());
        for (const GridOffset2D offset : ordered) {
            offsetRow.push_back(offset.row);
            offsetColumn.push_back(offset.column);
            const double squaredDistance = static_cast<double>(offset.row) * offset.row + static_cast<double>(offset.column) * offset.column;
            radius2 = std::max(radius2, squaredDistance);
        }
        radius = std::sqrt(radius2);
        buildForwardOffsetIndices();
        requireConnectedDomain();
    }
 
    void requireCoordinates(int row, int column) const {
        MMCFILTERS_CONTRACT_REQUIRE(row >= 0 && row < numRows && column >= 0 && column < numColumns, throw std::out_of_range("Index out of bounds."));
    }
 
    void requireLinearIndex(PixelId index) const {
        const std::int64_t domainSize = static_cast<std::int64_t>(numRows) * static_cast<std::int64_t>(numColumns);
        MMCFILTERS_CONTRACT_REQUIRE(index >= 0 && static_cast<std::int64_t>(index) < domainSize, throw std::out_of_range("Index out of bounds."));
    }
 
  public:
    RegularGridAdjacency2D(int numRows, int numColumns, double radius) {
        this->radius = checkedRadiusParameters(numRows, numColumns, radius);
        this->numRows = numRows;
        this->numColumns = numColumns;
        this->radius2 = radius * radius;
 
        int i, j, k, dx, dy, r0, r2, i0 = 0;
        this->n = 0;
        r0 = (int)radius;
        r2 = (int)radius2;
        for (dy = -r0; dy <= r0; dy++)
            for (dx = -r0; dx <= r0; dx++)
                if (((dx * dx) + (dy * dy)) <= r2)
                    this->n++;
 
        i = 0;
        this->offsetColumn.resize(this->n);
        this->offsetRow.resize(this->n);
 
        for (dy = -r0; dy <= r0; dy++) {
            for (dx = -r0; dx <= r0; dx++) {
                if (((dx * dx) + (dy * dy)) <= r2) {
                    this->offsetColumn[i] = dx;
                    this->offsetRow[i] = dy;
                    if ((dx == 0) && (dy == 0))
                        i0 = i;
                    i++;
                }
            }
        }
 
        float aux;
        std::vector<float> da(n);
        std::vector<float> dr(n);
 
        /* Set clockwise */
        for (i = 0; i < n; i++) {
            dx = this->offsetColumn[i];
            dy = this->offsetRow[i];
            dr[i] = std::sqrt((dx * dx) + (dy * dy));
            if (i != i0) {
                da[i] = (std::atan2(-dy, -dx) * 180.0 / std::numbers::pi);
                if (da[i] < 0.0)
                    da[i] += 360.0;
            }
        }
        da[i0] = 0.0;
        dr[i0] = 0.0;
 
        /* Place the central grid index first. */
        aux = da[i0];
        da[i0] = da[0];
        da[0] = aux;
 
        aux = dr[i0];
        dr[i0] = dr[0];
        dr[0] = aux;
 
        int auxX, auxY;
        auxX = this->offsetColumn[i0];
        auxY = this->offsetRow[i0];
        this->offsetColumn[i0] = this->offsetColumn[0];
        this->offsetRow[i0] = this->offsetRow[0];
 
        this->offsetColumn[0] = auxX;
        this->offsetRow[0] = auxY;
 
        /* sort by angle */
        for (i = 1; i < n - 1; i++) {
            k = i;
            for (j = i + 1; j < n; j++)
                if (da[j] < da[k]) {
                    k = j;
                }
            aux = da[i];
            da[i] = da[k];
            da[k] = aux;
            aux = dr[i];
            dr[i] = dr[k];
            dr[k] = aux;
 
            auxX = this->offsetColumn[i];
            auxY = this->offsetRow[i];
            this->offsetColumn[i] = this->offsetColumn[k];
            this->offsetRow[i] = this->offsetRow[k];
 
            this->offsetColumn[k] = auxX;
            this->offsetRow[k] = auxY;
        }
 
        /* sort by radius for each angle */
        for (i = 1; i < n - 1; i++) {
            k = i;
            for (j = i + 1; j < n; j++)
                if ((dr[j] < dr[k]) && (da[j] == da[k])) {
                    k = j;
                }
            aux = dr[i];
            dr[i] = dr[k];
            dr[k] = aux;
 
            auxX = this->offsetColumn[i];
            auxY = this->offsetRow[i];
            this->offsetColumn[i] = this->offsetColumn[k];
            this->offsetRow[i] = this->offsetRow[k];
 
            this->offsetColumn[k] = auxX;
            this->offsetRow[k] = auxY;
        }
 
        // Forward-only stencil keeps one orientation of every undirected edge
        // while preserving the complete clockwise stencil order.
        buildForwardOffsetIndices();
    }
 
    [[nodiscard]] static RegularGridAdjacency2D fromStructuringElement(int numRows, int numColumns, std::span<const GridOffset2D> offsets) {
        return RegularGridAdjacency2D(numRows, numColumns, std::vector<GridOffset2D>(offsets.begin(), offsets.end()), StructuringElementTag{});
    }
 
    [[nodiscard]] static RegularGridAdjacency2D rectangular(int numRows, int numColumns, int rowRadius, int columnRadius) {
        if (rowRadius < 0 || columnRadius < 0) {
            throw std::invalid_argument("Rectangular adjacency radii must be non-negative.");
        }
        const std::int64_t height = 2 * static_cast<std::int64_t>(rowRadius) + 1;
        const std::int64_t width = 2 * static_cast<std::int64_t>(columnRadius) + 1;
        const std::int64_t count = height * width;
        if (count > static_cast<std::int64_t>(std::numeric_limits<int>::max())) {
            throw std::length_error("Rectangular adjacency exceeds the supported offset count.");
        }
 
        std::vector<GridOffset2D> offsets;
        offsets.reserve(static_cast<std::size_t>(count));
        for (int rowOffset = -rowRadius; rowOffset <= rowRadius; ++rowOffset) {
            for (int columnOffset = -columnRadius; columnOffset <= columnRadius; ++columnOffset) {
                offsets.push_back({rowOffset, columnOffset});
            }
        }
        return fromStructuringElement(numRows, numColumns, offsets);
    }
 
    [[nodiscard]] static RegularGridAdjacency2D line(int numRows, int numColumns, int rowExtent, int columnExtent) {
        if (rowExtent == std::numeric_limits<int>::min() || columnExtent == std::numeric_limits<int>::min()) {
            throw std::invalid_argument("Line adjacency extents must be safely negatable.");
        }
        const std::int64_t absoluteRow = std::abs(static_cast<std::int64_t>(rowExtent));
        const std::int64_t absoluteColumn = std::abs(static_cast<std::int64_t>(columnExtent));
        const std::int64_t steps = std::max(absoluteRow, absoluteColumn);
        if (2 * steps + 1 > static_cast<std::int64_t>(std::numeric_limits<int>::max())) {
            throw std::length_error("Line adjacency exceeds the supported offset count.");
        }
        if (steps == 0) {
            const GridOffset2D origin{0, 0};
            return fromStructuringElement(numRows, numColumns, std::span<const GridOffset2D>(&origin, 1));
        }
 
        auto roundedRatio = [](std::int64_t numerator, std::int64_t denominator) {
            if (numerator >= 0) {
                return (numerator + denominator / 2) / denominator;
            }
            return -((-numerator + denominator / 2) / denominator);
        };
 
        std::vector<GridOffset2D> offsets;
        offsets.reserve(static_cast<std::size_t>(2 * steps + 1));
        for (std::int64_t sample = -steps; sample <= steps; ++sample) {
            offsets.push_back({static_cast<int>(roundedRatio(sample * rowExtent, steps)), static_cast<int>(roundedRatio(sample * columnExtent, steps))});
        }
        return fromStructuringElement(numRows, numColumns, offsets);
    }
 
    [[nodiscard]] static RegularGridAdjacency2D horizontalLine(int numRows, int numColumns, int halfLength) {
        if (halfLength < 0) {
            throw std::invalid_argument("Horizontal-line half-length must be non-negative.");
        }
        return line(numRows, numColumns, 0, halfLength);
    }
 
    [[nodiscard]] static RegularGridAdjacency2D verticalLine(int numRows, int numColumns, int halfLength) {
        if (halfLength < 0) {
            throw std::invalid_argument("Vertical-line half-length must be non-negative.");
        }
        return line(numRows, numColumns, halfLength, 0);
    }
 
    int getSize() const noexcept { return this->n; }
 
    int getNumRows() const noexcept { return numRows; }
 
    int getNumColumns() const noexcept { return numColumns; }
 
    RegularGridAdjacencyShape getShape() const noexcept { return shape; }
 
    inline bool isAdjacent(PixelId p, PixelId q) const noexcept {
        if (numColumns <= 0) {
            return false;
        }
        int py = p / numColumns, px = p % numColumns;
        int qy = q / numColumns, qx = q % numColumns;
 
        return isAdjacent(px, py, qx, qy);
    }
 
    inline bool isAdjacent(int px, int py, int qx, int qy) const noexcept {
        const std::int64_t dx = static_cast<std::int64_t>(px) - qx;
        const std::int64_t dy = static_cast<std::int64_t>(py) - qy;
        if (shape == RegularGridAdjacencyShape::EuclideanDisk) {
            return static_cast<double>(dx) * dx + static_cast<double>(dy) * dy <= radius2;
        }
        if (dx < std::numeric_limits<int>::min() || dx > std::numeric_limits<int>::max() || dy < std::numeric_limits<int>::min() ||
            dy > std::numeric_limits<int>::max()) {
            return false;
        }
        return containsOffset(static_cast<int>(dy), static_cast<int>(dx));
    }
 
    double getRadius() const noexcept { return this->radius; }
 
    bool is4connectivity() const noexcept {
        if (shape == RegularGridAdjacencyShape::EuclideanDisk) {
            return radius == 1.0;
        }
        return matchesOffsets({{0, 0}, {-1, 0}, {0, -1}, {1, 0}, {0, 1}});
    }
 
    bool is8connectivity() const noexcept {
        if (shape == RegularGridAdjacencyShape::EuclideanDisk) {
            return radius == 1.5;
        }
        return matchesOffsets({{0, 0}, {-1, -1}, {-1, 0}, {-1, 1}, {0, -1}, {0, 1}, {1, -1}, {1, 0}, {1, 1}});
    }
 
    bool isCanonical4Or8Connectivity() const noexcept { return is4connectivity() || is8connectivity(); }
 
    bool isGridBoundary(PixelId index) const {
        requireLinearIndex(index);
        return isGridBoundary(index / numColumns, index % numColumns);
    }
 
    bool isGridBoundary(int row, int column) const noexcept { return row == 0 || column == 0 || row == this->numRows - 1 || column == this->numColumns - 1; }
 
    int getOffsetRow(int index) const noexcept { return offsetRow[index]; }
 
    int getOffsetColumn(int index) const noexcept { return offsetColumn[index]; }
 
    template <bool ForwardOnly> class IteratorAdjacencyT {
      private:
        const RegularGridAdjacency2D* relation_ = nullptr;
        int row_ = 0;
        int column_ = 0;
        int index_ = 0;
 
        [[nodiscard]] int stencilSize() const noexcept {
            if constexpr (ForwardOnly) {
                return static_cast<int>(relation_->forwardOffsetIndices.size());
            }
            return relation_->n;
        }
 
        [[nodiscard]] int stencilIndex() const noexcept {
            if constexpr (ForwardOnly) {
                return relation_->forwardOffsetIndices[static_cast<std::size_t>(index_)];
            }
            return index_;
        }
 
        void seekValid() noexcept {
            const int size = stencilSize();
            while (index_ < size) {
                const int offsetIndex = stencilIndex();
                const std::int64_t neighborRow = static_cast<std::int64_t>(row_) + relation_->offsetRow[static_cast<std::size_t>(offsetIndex)];
                const std::int64_t neighborColumn = static_cast<std::int64_t>(column_) + relation_->offsetColumn[static_cast<std::size_t>(offsetIndex)];
                if (neighborRow >= 0 && neighborRow < relation_->numRows && neighborColumn >= 0 && neighborColumn < relation_->numColumns) {
                    return;
                }
                ++index_;
            }
        }
 
      public:
        using iterator_category = std::forward_iterator_tag;
        using iterator_concept = std::forward_iterator_tag;
 
        using value_type = PixelId;
        using difference_type = std::ptrdiff_t;
        using reference = PixelId;
        using pointer = void;
 
        IteratorAdjacencyT() noexcept = default;
 
        IteratorAdjacencyT(const RegularGridAdjacency2D* relation, int row, int column, int index) noexcept
            : relation_(relation), row_(row), column_(column), index_(index) {
            seekValid();
        }
 
        IteratorAdjacencyT& operator++() noexcept {
            ++index_;
            seekValid();
            return *this;
        }
 
        IteratorAdjacencyT operator++(int) noexcept {
            IteratorAdjacencyT previous = *this;
            ++(*this);
            return previous;
        }
 
        bool operator==(const IteratorAdjacencyT& other) const noexcept {
            return relation_ == other.relation_ && row_ == other.row_ && column_ == other.column_ && index_ == other.index_;
        }
 
        bool operator!=(const IteratorAdjacencyT& other) const noexcept { return !(*this == other); }
 
        PixelId operator*() const noexcept {
            const int offsetIndex = stencilIndex();
            const std::int64_t neighborRow = static_cast<std::int64_t>(row_) + relation_->offsetRow[static_cast<std::size_t>(offsetIndex)];
            const std::int64_t neighborColumn = static_cast<std::int64_t>(column_) + relation_->offsetColumn[static_cast<std::size_t>(offsetIndex)];
            return static_cast<PixelId>(neighborRow * relation_->numColumns + neighborColumn);
        }
    };
 
    using IteratorAdjacency = IteratorAdjacencyT<false>;
    using ForwardIteratorAdjacency = IteratorAdjacencyT<true>;
 
    template <bool ForwardOnly, int FirstOffset> class GridIndexRangeT {
      private:
        const RegularGridAdjacency2D* relation_;
        int row_;
        int column_;
 
      public:
        GridIndexRangeT(const RegularGridAdjacency2D& relation, int row, int column) noexcept : relation_(&relation), row_(row), column_(column) {}
 
        using iterator = IteratorAdjacencyT<ForwardOnly>;
 
        [[nodiscard]] iterator begin() const noexcept { return iterator(relation_, row_, column_, FirstOffset); }
 
        [[nodiscard]] iterator end() const noexcept {
            const int endIndex = [&] {
                if constexpr (ForwardOnly) {
                    return static_cast<int>(relation_->forwardOffsetIndices.size());
                }
                return relation_->n;
            }();
            return iterator(relation_, row_, column_, endIndex);
        }
    };
 
    using AdjacentIndexRange = GridIndexRangeT<false, 0>;
    using NeighborIndexRange = GridIndexRangeT<false, 1>;
    using ForwardNeighborIndexRange = GridIndexRangeT<true, 0>;
 
    [[nodiscard]] AdjacentIndexRange getAdjacentIndices(int row, int column) const {
        requireCoordinates(row, column);
        return AdjacentIndexRange(*this, row, column);
    }
 
    [[nodiscard]] AdjacentIndexRange getAdjacentIndices(PixelId gridIndex) const {
        requireLinearIndex(gridIndex);
        return getAdjacentIndices(gridIndex / numColumns, gridIndex % numColumns);
    }
 
    [[nodiscard]] NeighborIndexRange getNeighborIndices(int row, int column) const {
        requireCoordinates(row, column);
        return NeighborIndexRange(*this, row, column);
    }
 
    [[nodiscard]] NeighborIndexRange getNeighborIndices(PixelId gridIndex) const {
        requireLinearIndex(gridIndex);
        return getNeighborIndices(gridIndex / numColumns, gridIndex % numColumns);
    }
 
    [[nodiscard]] ForwardNeighborIndexRange getForwardNeighborIndices(int row, int column) const {
        requireCoordinates(row, column);
        return ForwardNeighborIndexRange(*this, row, column);
    }
 
    [[nodiscard]] ForwardNeighborIndexRange getForwardNeighborIndices(PixelId gridIndex) const {
        requireLinearIndex(gridIndex);
        return getForwardNeighborIndices(gridIndex / numColumns, gridIndex % numColumns);
    }
};
 
} // namespace mmcfilters