bitmatrix.hpp

Go to the documentation of this file.

#pragma once
 
#include "algorithm/polyominoes.hpp"

 
#include <algorithm>
#include <array>
#include <bit>
#include <cassert>
#include <cmath>
#include <compare>
#include <concepts>
#include <cstddef>
#include <cstdint>
#include <functional>
#include <iterator>
#include <optional>
#include <ranges>
#include <span>
#include <stdexcept>
#include <type_traits>
#include <utility>
#include <vector>
 
namespace pgl {

enum class GridAdjacency { edge, vertex };
 
// Declared here so the point-range constructor can constrain on it: a matrix is
// itself a range of points, and it must not fall into that overload.
template <class TPointType>
    requires std::signed_integral<typename TPointType::NumberType>
class BitMatrix;
 
namespace detail {

template <class T>
inline constexpr bool is_bit_matrix_v = false;
 
template <class TPointType>
inline constexpr bool is_bit_matrix_v<BitMatrix<TPointType>> = true;

template <class Int, class Number>
[[nodiscard]] Int gridCoordinate(const Number& value) {
    if constexpr (std::same_as<Int, Number>) {
        return value;
    } else if constexpr (is_Rational_v<Number>) {
        if (!value.isInteger()) {
            throw std::logic_error("pgl::BitMatrix: a coordinate is not an integer");
        }
        return gridCoordinate<Int>(static_cast<rational_int_t<Number>>(value));
    } else if constexpr (std::is_floating_point_v<Number>) {
        if (!std::isfinite(value) || value != std::floor(value)) {
            throw std::logic_error("pgl::BitMatrix: a coordinate is not an integer");
        }
        // The bound is a power of two, so it and its negation are both exact in
        // any binary floating-point type: this range test never rounds, even
        // where Int carries more digits than Number does.
        const Number low = static_cast<Number>(numeric_limits<Int>::min());
        if (value < low || value >= -low) {
            throw std::logic_error("pgl::BitMatrix: a coordinate does not fit the grid");
        }
        return static_cast<Int>(value);
    } else {
        // An exact integer of some other width. Range-check it whenever Number
        // can spell Int's bounds at all; an unbounded BigInt reports no digits
        // and always can, a narrower fixed width cannot go out of range.
        constexpr bool comparable = !numeric_limits<Number>::is_specialized
                                    || numeric_limits<Number>::digits >= numeric_limits<Int>::digits;
        if constexpr (comparable) {
            if (value < static_cast<Number>(numeric_limits<Int>::min())
                || value > static_cast<Number>(numeric_limits<Int>::max())) {
                throw std::logic_error("pgl::BitMatrix: a coordinate does not fit the grid");
            }
        }
        return static_cast<Int>(value);
    }
}

template <class GridPointType, class OtherPointType>
[[nodiscard]] GridPointType gridPoint(const OtherPointType& point) {
    using Int = typename GridPointType::NumberType;
    return GridPointType(gridCoordinate<Int>(point.x()), gridCoordinate<Int>(point.y()));
}

template <class GridPointType, class OtherPointType, class TLabel>
[[nodiscard]] Polygon<GridPointType> gridRing(const Polygon<OtherPointType, TLabel>& ring) {
    std::vector<GridPointType> vertices;
    vertices.reserve(ring.size());
    for (const OtherPointType& vertex : ring) {
        vertices.push_back(gridPoint<GridPointType>(vertex));
    }
    return Polygon<GridPointType>(std::move(vertices), true);
}

template <class GridPointType, class OtherPointType, class TLabel>
[[nodiscard]] PolygonWithHoles<GridPointType> gridRegion(
    const PolygonWithHoles<OtherPointType, TLabel>& region) {
    std::vector<Polygon<GridPointType>> holes;
    holes.reserve(region.holeCount());
    for (const auto& hole : region.holes()) {
        holes.push_back(gridRing<GridPointType>(hole));
    }
    return PolygonWithHoles<GridPointType>(gridRing<GridPointType>(region.outer()), std::move(holes),
                                           true);
}

template <class GridPointType, class OtherPointType, class TLabel>
[[nodiscard]] PolygonSet<GridPointType> gridSet(const PolygonSet<OtherPointType, TLabel>& set) {
    std::vector<PolygonWithHoles<GridPointType>> components;
    components.reserve(set.componentCount());
    for (const auto& component : set.components()) {
        components.push_back(gridRegion<GridPointType>(component));
    }
    return PolygonSet<GridPointType>(std::move(components), true);
}

template <class PointType, class Int>
using grid_point_t = std::conditional_t<std::same_as<Int, typename PointType::NumberType>, PointType,
                                        Point<Int, typename PointType::LabelType>>;
 
}  // namespace detail

template <class TPointType = Point<int>>
    requires std::signed_integral<typename TPointType::NumberType>
class BitMatrix {
public:
    using PointType = TPointType;
    using NumberType = typename PointType::NumberType;
    using RectangleType = Rectangle<PointType>;
    using RegionType = PolygonWithHoles<PointType>;
    using PolygonSetType = PolygonSet<PointType>;
    using ConvexType = Convex<PointType>;
    using value_type = PointType;
 
    class Iterator;
    using iterator = Iterator;
    using const_iterator = Iterator;
 
    // -----------------------------------------------------------------------
    // Construction

    BitMatrix() = default;

    BitMatrix(PointType origin, int width, int height) {
        if (width > 0 && height > 0) {
            origin_ = std::move(origin);
            width_ = width;
            height_ = height;
        }
        words_ = wordsPerRow(width_);
        bits_.assign(words_ * static_cast<std::size_t>(height_), 0);
    }

    explicit BitMatrix(const RectangleType& box)
        : BitMatrix(box.empty() ? PointType() : box.min(),
                    box.empty() ? 0 : static_cast<int>(box.width()),
                    box.empty() ? 0 : static_cast<int>(box.height())) {}

    template <class TLabel>
    explicit BitMatrix(const PolygonWithHoles<PointType, TLabel>& region)
        : BitMatrix(RectangleType(region.bbox())) {
        Crossings crossings(static_cast<std::size_t>(height_));
        addRegionCrossings(region, crossings);
        fillCrossings(crossings);
    }

    // Unconstrained: partial ordering already prefers the same-type overload
    // above, and a constraint here would ride onto an implicit deduction guide,
    // which the CI clang mishandles.
    template <class OtherPointType, class TLabel>
    explicit BitMatrix(const PolygonWithHoles<OtherPointType, TLabel>& region)
        : BitMatrix(detail::gridRegion<PointType>(region)) {}

    template <class TLabel>
    explicit BitMatrix(const Polygon<PointType, TLabel>& polygon)
        : BitMatrix(RectangleType(polygon.bbox())) {
        Crossings crossings(static_cast<std::size_t>(height_));
        addRingCrossings(polygon, crossings);
        fillCrossings(crossings);
    }

    template <class OtherPointType, class TLabel>
    explicit BitMatrix(const Polygon<OtherPointType, TLabel>& polygon)
        : BitMatrix(detail::gridRing<PointType>(polygon)) {}

    template <class TLabel>
    explicit BitMatrix(const PolygonSet<PointType, TLabel>& set)
        : BitMatrix(RectangleType(set.bbox())) {
        Crossings crossings(static_cast<std::size_t>(height_));
        for (const PolygonWithHoles<PointType>& component : set.components()) {
            addRegionCrossings(component, crossings);
        }
        fillCrossings(crossings);
    }

    template <class OtherPointType, class TLabel>
    explicit BitMatrix(const PolygonSet<OtherPointType, TLabel>& set)
        : BitMatrix(detail::gridSet<PointType>(set)) {}

    template <std::ranges::input_range Range = std::initializer_list<PointType>>
        requires(detail::is_point_v<std::remove_cvref_t<std::ranges::range_value_t<Range>>>
                 && !AnyShapeConcept<std::remove_cvref_t<Range>>
                 && !detail::is_bit_matrix_v<std::remove_cvref_t<Range>>)
    explicit BitMatrix(Range&& points) {
        using SourcePoint = std::remove_cvref_t<std::ranges::range_value_t<Range>>;
        if constexpr (std::ranges::forward_range<Range>
                      && std::same_as<typename SourcePoint::NumberType, NumberType>) {
            *this = emptyOver(points);
            for (const auto& point : points) {
                set(point.x(), point.y());
            }
        } else {
            // A single-pass range cannot be walked twice, and a converting one
            // should not be: the window needs a pass of its own, so convert
            // once into a range that gives both.
            std::vector<PointType> cells;
            for (const auto& point : points) {
                cells.push_back(detail::gridPoint<PointType>(point));
            }
            *this = BitMatrix(cells);
        }
    }
 
    // -----------------------------------------------------------------------
    // The window

    [[nodiscard]] const PointType& origin() const { return origin_; }

    [[nodiscard]] int width() const { return width_; }

    [[nodiscard]] int height() const { return height_; }

    [[nodiscard]] RectangleType window() const {
        if (emptyWindow()) {
            return RectangleType();
        }
        return RectangleType(origin_, PointType(origin_.x() + static_cast<NumberType>(width_),
                                                origin_.y() + static_cast<NumberType>(height_)));
    }

    [[nodiscard]] bool emptyWindow() const { return width_ == 0; }

    [[nodiscard]] bool inWindow(NumberType x, NumberType y) const {
        const std::int64_t i = localX(x), j = localY(y);
        return i >= 0 && i < width_ && j >= 0 && j < height_;
    }

    [[nodiscard]] bool inWindow(const PointType& cell) const { return inWindow(cell.x(), cell.y()); }

    [[nodiscard]] bool sameWindow(const BitMatrix& other) const {
        return width_ == other.width_ && height_ == other.height_ && origin_ == other.origin_;
    }

    [[nodiscard]] BitMatrix resized(const RectangleType& box) const {
        BitMatrix result(box);
        result.combine(*this, [](std::uint64_t, std::uint64_t source) { return source; });
        return result;
    }

    [[nodiscard]] BitMatrix trimmed() const { return resized(bbox()); }
 
    // -----------------------------------------------------------------------
    // Individual cells

    [[nodiscard]] bool get(NumberType x, NumberType y) const {
        const std::int64_t i = localX(x), j = localY(y);
        if (i < 0 || i >= width_ || j < 0 || j >= height_) {
            return false;
        }
        return ((row(static_cast<int>(j))[static_cast<std::size_t>(i) / 64] >> (i % 64)) & 1) != 0;
    }

    [[nodiscard]] bool get(const PointType& cell) const { return get(cell.x(), cell.y()); }

    void set(NumberType x, NumberType y) {
        const std::int64_t i = localX(x), j = localY(y);
        if (i < 0 || i >= width_ || j < 0 || j >= height_) {
            return;
        }
        row(static_cast<int>(j))[static_cast<std::size_t>(i) / 64] |= std::uint64_t(1) << (i % 64);
    }

    void set(const PointType& cell) { set(cell.x(), cell.y()); }

    void set(NumberType x, NumberType y, bool value) {
        if (value) {
            set(x, y);
        } else {
            reset(x, y);
        }
    }

    template <std::ranges::input_range Range = std::initializer_list<PointType>>
        requires(detail::is_point_v<std::remove_cvref_t<std::ranges::range_value_t<Range>>>
                 && !AnyShapeConcept<std::remove_cvref_t<Range>>
                 && !detail::is_bit_matrix_v<std::remove_cvref_t<Range>>)
    void set(Range&& points) {
        for (const auto& point : points) {
            set(detail::gridPoint<PointType>(point));
        }
    }

    void reset(NumberType x, NumberType y) {
        const std::int64_t i = localX(x), j = localY(y);
        if (i < 0 || i >= width_ || j < 0 || j >= height_) {
            return;
        }
        row(static_cast<int>(j))[static_cast<std::size_t>(i) / 64] &= ~(std::uint64_t(1) << (i % 64));
    }

    void reset(const PointType& cell) { reset(cell.x(), cell.y()); }

    template <std::ranges::input_range Range = std::initializer_list<PointType>>
        requires(detail::is_point_v<std::remove_cvref_t<std::ranges::range_value_t<Range>>>
                 && !AnyShapeConcept<std::remove_cvref_t<Range>>
                 && !detail::is_bit_matrix_v<std::remove_cvref_t<Range>>)
    void reset(Range&& points) {
        for (const auto& point : points) {
            reset(detail::gridPoint<PointType>(point));
        }
    }

    void flip(NumberType x, NumberType y) {
        const std::int64_t i = localX(x), j = localY(y);
        if (i < 0 || i >= width_ || j < 0 || j >= height_) {
            return;
        }
        row(static_cast<int>(j))[static_cast<std::size_t>(i) / 64] ^= std::uint64_t(1) << (i % 64);
    }

    void flip(const PointType& cell) { flip(cell.x(), cell.y()); }

    template <std::ranges::input_range Range = std::initializer_list<PointType>>
        requires(detail::is_point_v<std::remove_cvref_t<std::ranges::range_value_t<Range>>>
                 && !AnyShapeConcept<std::remove_cvref_t<Range>>
                 && !detail::is_bit_matrix_v<std::remove_cvref_t<Range>>)
    void flip(Range&& points) {
        for (const auto& point : points) {
            flip(detail::gridPoint<PointType>(point));
        }
    }

    void setAll() {
        std::fill(bits_.begin(), bits_.end(), ~std::uint64_t(0));
        maskTails();
    }

    void clear() { std::fill(bits_.begin(), bits_.end(), 0); }
 
    // -----------------------------------------------------------------------
    // Cardinality

    [[nodiscard]] bool empty() const {
        for (std::uint64_t word : bits_) {
            if (word != 0) {
                return false;
            }
        }
        return true;
    }

    [[nodiscard]] std::size_t count() const {
        std::size_t total = 0;
        for (std::uint64_t word : bits_) {
            total += static_cast<std::size_t>(std::popcount(word));
        }
        return total;
    }

    template <class ResultNumber = NumberType>
    [[nodiscard]] ResultNumber area() const {
        return ResultNumber(static_cast<std::int64_t>(count()));
    }

    template <class ResultNumber = NumberType>
    [[nodiscard]] ResultNumber perimeter() const {
        std::size_t adjacent = 0;
        for (int j = 0; j < height_; ++j) {
            const std::uint64_t* here = row(j);
            const std::uint64_t* above = j + 1 < height_ ? row(j + 1) : nullptr;
            for (std::size_t w = 0; w < words_; ++w) {
                const std::uint64_t word = here[w];
                if (word == 0) {
                    continue;
                }
                adjacent += static_cast<std::size_t>(
                    std::popcount(word & shiftedWord(here, words_, static_cast<std::int64_t>(w) * 64 + 1)));
                if (above != nullptr) {
                    adjacent += static_cast<std::size_t>(std::popcount(word & above[w]));
                }
            }
        }
        return ResultNumber(static_cast<std::int64_t>(4 * count() - 2 * adjacent));
    }

    template <class ResultNumber = division_result_t<NumberType>>
    [[nodiscard]] Point<ResultNumber> centroid() const {
        // Summed a word at a time: the cells of a word share a row and a base
        // column, so only their bit offsets are walked.
        std::int64_t sumX = 0, sumY = 0, cells = 0;
        for (int j = 0; j < height_; ++j) {
            const std::uint64_t* here = row(j);
            const std::int64_t y = static_cast<std::int64_t>(origin_.y()) + j;
            for (std::size_t w = 0; w < words_; ++w) {
                std::uint64_t word = here[w];
                if (word == 0) {
                    continue;
                }
                const std::int64_t inWord = std::popcount(word);
                const std::int64_t base =
                    static_cast<std::int64_t>(origin_.x()) + static_cast<std::int64_t>(w) * 64;
                std::int64_t offsets = 0;
                for (; word != 0; word &= word - 1) {
                    offsets += std::countr_zero(word);
                }
                sumX += base * inWord + offsets;
                sumY += y * inWord;
                cells += inWord;
            }
        }
        if (cells == 0) {
            throw std::logic_error("pgl::BitMatrix::centroid: no cell is set");
        }
        // Twice the sums over twice the count, which is the cell-center offset
        // of one half folded into an exact integer ratio.
        const ResultNumber total = ResultNumber(2 * cells);
        return Point<ResultNumber>(ResultNumber(2 * sumX + cells) / total,
                                   ResultNumber(2 * sumY + cells) / total);
    }
 
    // -----------------------------------------------------------------------
    // Reading the cells out

    [[nodiscard]] Iterator begin() const { return Iterator(this); }

    [[nodiscard]] Iterator end() const { return Iterator(); }

    [[nodiscard]] auto latticeView() const { return std::ranges::subrange(begin(), end()); }

    [[nodiscard]] auto cellsView() const {
        return latticeView() |
               std::views::transform([](const PointType& cell) { return cellSquare(cell); });
    }

    [[nodiscard]] std::vector<PointType> lattice() const {
        std::vector<PointType> result;
        result.reserve(count());
        for (const PointType& cell : *this) {
            result.push_back(cell);
        }
        return result;
    }

    [[nodiscard]] std::vector<RectangleType> cells() const {
        std::vector<RectangleType> result;
        result.reserve(count());
        for (const PointType& cell : *this) {
            result.push_back(cellSquare(cell));
        }
        return result;
    }

    [[nodiscard]] std::vector<RectangleType> rectangles() const {
        std::vector<RectangleType> result;
        for (int j = 0; j < height_; ++j) {
            std::int64_t i = 0;
            while (i < width_) {
                if (!localGet(j, i)) {
                    ++i;
                    continue;
                }
                const std::int64_t start = i;
                while (i < width_ && localGet(j, i)) {
                    ++i;
                }
                result.emplace_back(
                    PointType(origin_.x() + static_cast<NumberType>(start),
                              origin_.y() + static_cast<NumberType>(j)),
                    PointType(origin_.x() + static_cast<NumberType>(i),
                              origin_.y() + static_cast<NumberType>(j + 1)));
            }
        }
        return result;
    }

    [[nodiscard]] RectangleType bbox() const {
        std::int64_t minX = width_, maxX = -1, minY = height_, maxY = -1;
        for (int j = 0; j < height_; ++j) {
            const std::uint64_t* here = row(j);
            std::size_t first = 0;
            while (first < words_ && here[first] == 0) {
                ++first;
            }
            if (first == words_) {
                continue;
            }
            std::size_t last = words_ - 1;
            while (here[last] == 0) {
                --last;
            }
            minX = std::min(minX, static_cast<std::int64_t>(first) * 64 + std::countr_zero(here[first]));
            maxX = std::max(maxX, static_cast<std::int64_t>(last) * 64 + 63 - std::countl_zero(here[last]));
            minY = std::min(minY, static_cast<std::int64_t>(j));
            maxY = std::max(maxY, static_cast<std::int64_t>(j));
        }
        if (maxX < 0) {
            return RectangleType();
        }
        return RectangleType(PointType(origin_.x() + static_cast<NumberType>(minX),
                                       origin_.y() + static_cast<NumberType>(minY)),
                             PointType(origin_.x() + static_cast<NumberType>(maxX + 1),
                                       origin_.y() + static_cast<NumberType>(maxY + 1)));
    }

    template <class ResultNumber = double>
    [[nodiscard]] Rectangle<Point<ResultNumber>> fbox() const {
        return bbox().template fbox<ResultNumber>();
    }

    template <class ResultNumber = division_result_t<NumberType>>
    [[nodiscard]] Point<ResultNumber> pointInside() const {
        const std::optional<PointType> cell = firstSetCell();
        if (!cell) {
            throw std::logic_error("pgl::BitMatrix::pointInside: no cell is set");
        }
        const ResultNumber half = ResultNumber(1) / ResultNumber(2);
        return Point<ResultNumber>(ResultNumber(cell->x()) + half, ResultNumber(cell->y()) + half);
    }

    [[nodiscard]] RegionType asPolygonWithHoles() const {
        if (componentCount() > 1) {
            throw std::logic_error("pgl::BitMatrix::asPolygonWithHoles: the cells are not edge-connected");
        }
        return regionFromLoops(boundaryLoops(
            [](const detail::PolyCell&, const detail::PolyCell&) { return true; }));
    }

    [[nodiscard]] PolygonSetType asPolygonSet() const {
        if (emptyWindow()) {
            return PolygonSetType();
        }
 
        // Label the components first, so that each traced loop can be handed to
        // the one it bounds: a loop's first edge names a set cell on its left,
        // and the sorted runs of that cell's row say which component holds it.
        struct LabeledRun {
            std::int64_t x0, x1;
            std::size_t component;
        };
        std::vector<std::vector<LabeledRun>> rowRuns(static_cast<std::size_t>(height_));
        std::size_t total = 0;
        visitComponentRuns(GridAdjacency::edge, [&](const std::vector<CellRun>& runs) {
            for (const CellRun& run : runs) {
                rowRuns[static_cast<std::size_t>(run.y)].push_back({run.x0, run.x1, total});
            }
            ++total;
        });
        if (total == 0) {
            return PolygonSetType();
        }
        for (std::vector<LabeledRun>& runs : rowRuns) {
            std::sort(runs.begin(), runs.end(),
                      [](const LabeledRun& a, const LabeledRun& b) { return a.x0 < b.x0; });
        }
 
        auto componentAt = [&](const detail::PolyCell& cell) {
            const std::vector<LabeledRun>& runs = rowRuns[static_cast<std::size_t>(cell.second)];
            const auto after = std::upper_bound(
                runs.begin(), runs.end(), static_cast<std::int64_t>(cell.first),
                [](std::int64_t at, const LabeledRun& run) { return at < run.x0; });
            assert(after != runs.begin() && cell.first < std::prev(after)->x1 &&
                   "pgl::BitMatrix::asPolygonSet: a boundary edge borders no cell");
            return std::prev(after)->component;
        };
 
        // Two groups touching only at a corner are two components, so a pinch
        // there stays on its own side; a pinch within one group crosses, which
        // is what keeps a hole's loop apart from the loop around it.
        std::vector<std::vector<std::vector<detail::PolyCell>>> grouped(total);
        for (std::vector<detail::PolyCell>& loop :
             boundaryLoops([&](const detail::PolyCell& here, const detail::PolyCell& across) {
                 return componentAt(here) == componentAt(across);
             })) {
            grouped[componentAt(loopSeedCell(loop))].push_back(std::move(loop));
        }
 
        std::vector<RegionType> components;
        components.reserve(total);
        for (const std::vector<std::vector<detail::PolyCell>>& loops : grouped) {
            components.push_back(regionFromLoops(loops));
        }
        return PolygonSetType(components);
    }

    [[nodiscard]] ConvexType convexHull() const {
        std::vector<PointType> corners;
        corners.reserve(static_cast<std::size_t>(height_) * 4);
        for (int j = 0; j < height_; ++j) {
            const std::optional<std::pair<std::int64_t, std::int64_t>> extent = rowExtent(j);
            if (!extent) {
                continue;
            }
            const NumberType low = origin_.x() + static_cast<NumberType>(extent->first);
            const NumberType high = origin_.x() + static_cast<NumberType>(extent->second + 1);
            const NumberType bottom = origin_.y() + static_cast<NumberType>(j);
            const NumberType top = bottom + NumberType(1);
            corners.emplace_back(low, bottom);
            corners.emplace_back(low, top);
            corners.emplace_back(high, bottom);
            corners.emplace_back(high, top);
        }
        return ConvexType(corners);
    }

    friend Canvas& operator<<(Canvas& canvas, const BitMatrix& matrix) {
        return canvas << matrix.asPolygonSet();
    }
 
    // -----------------------------------------------------------------------
    // Set algebra
    //
    // Every binary operator returns the smallest window that provably loses no
    // cell: the overlap of the two windows for an intersection, their hull for a
    // union or a symmetric difference, the left window for a difference. The
    // compound assignments instead never move their window, and drop whatever
    // falls outside it, exactly as `set` does.

    BitMatrix operator~() const {
        BitMatrix result(origin_, width_, height_);
        for (std::size_t w = 0; w < bits_.size(); ++w) {
            result.bits_[w] = ~bits_[w];
        }
        result.maskTails();
        return result;
    }

    BitMatrix& operator&=(const BitMatrix& other) {
        combine(other, [](std::uint64_t left, std::uint64_t right) { return left & right; });
        return *this;
    }

    BitMatrix& operator|=(const BitMatrix& other) {
        combine(other, [](std::uint64_t left, std::uint64_t right) { return left | right; });
        return *this;
    }

    BitMatrix& operator^=(const BitMatrix& other) {
        combine(other, [](std::uint64_t left, std::uint64_t right) { return left ^ right; });
        return *this;
    }

    BitMatrix operator&(const BitMatrix& other) const {
        BitMatrix result = resized(overlapWindow(*this, other));
        result &= other;
        return result;
    }

    BitMatrix operator|(const BitMatrix& other) const {
        BitMatrix result = resized(hullWindow(*this, other));
        result |= other;
        return result;
    }

    BitMatrix operator^(const BitMatrix& other) const {
        BitMatrix result = resized(hullWindow(*this, other));
        result ^= other;
        return result;
    }

    [[nodiscard]] BitMatrix difference(const BitMatrix& other) const {
        BitMatrix result = *this;
        result.combine(other, [](std::uint64_t left, std::uint64_t right) { return left & ~right; });
        return result;
    }

    [[nodiscard]] BitMatrix symmetricDifference(const BitMatrix& other) const { return *this ^ other; }

    bool operator==(const BitMatrix& other) const {
        return width_ == other.width_ && height_ == other.height_ && origin_ == other.origin_ &&
               bits_ == other.bits_;
    }

    std::strong_ordering operator<=>(const BitMatrix& other) const {
        if (const std::strong_ordering order = origin_ <=> other.origin_; order != 0) {
            return order;
        }
        if (const std::strong_ordering order = width_ <=> other.width_; order != 0) {
            return order;
        }
        if (const std::strong_ordering order = height_ <=> other.height_; order != 0) {
            return order;
        }
        return bits_ <=> other.bits_;
    }

    [[nodiscard]] bool samePointSet(const BitMatrix& other) const {
        const std::size_t here = count();
        return here == other.count() && andCount(other) == here;
    }

    [[nodiscard]] bool contains(const BitMatrix& other) const { return other.andCount(*this) == other.count(); }

    [[nodiscard]] bool interiorContains(const BitMatrix& other) const {
        return interior(GridAdjacency::vertex).contains(other);
    }

    [[nodiscard]] bool boundaryContains(const BitMatrix& other) const { return other.empty(); }

    [[nodiscard]] bool intersects(const BitMatrix& other) const {
        for (std::int64_t dy = -1; dy <= 1; ++dy) {
            for (std::int64_t dx = -1; dx <= 1; ++dx) {
                if (anyCommonCell(other, dx, dy)) {
                    return true;
                }
            }
        }
        return false;
    }

    [[nodiscard]] bool interiorsIntersect(const BitMatrix& other) const { return anyCommonCell(other, 0, 0); }

    [[nodiscard]] std::size_t andCount(const BitMatrix& other) const {
        std::size_t total = 0;
        forEachAlignedWord(other, [&](std::uint64_t left, std::uint64_t right) {
            total += static_cast<std::size_t>(std::popcount(left & right));
        });
        return total;
    }

    [[nodiscard]] std::size_t orCount(const BitMatrix& other) const {
        return count() + other.count() - andCount(other);
    }

    [[nodiscard]] std::size_t xorCount(const BitMatrix& other) const {
        return count() + other.count() - 2 * andCount(other);
    }
 
    // -----------------------------------------------------------------------
    // Translations, reflections and morphology
    //
    // These read a cell as the lattice point at its lower-left corner, so the
    // sum of two cells is one cell and a reflection maps cell `c` to cell `-c`.

    [[nodiscard]] BitMatrix translated(const PointType& vector) const {
        BitMatrix result = *this;
        result.origin_ = PointType(origin_.x() + vector.x(), origin_.y() + vector.y());
        return result;
    }

    BitMatrix operator+(const PointType& vector) const { return translated(vector); }

    BitMatrix operator-(const PointType& vector) const {
        return translated(PointType(-vector.x(), -vector.y()));
    }

    BitMatrix& operator+=(const PointType& vector) {
        origin_ = PointType(origin_.x() + vector.x(), origin_.y() + vector.y());
        return *this;
    }

    BitMatrix& operator-=(const PointType& vector) {
        origin_ = PointType(origin_.x() - vector.x(), origin_.y() - vector.y());
        return *this;
    }

    [[nodiscard]] BitMatrix reflected() const { return mapped([](std::int64_t x, std::int64_t y) {
        return std::pair<std::int64_t, std::int64_t>(-x - 1, -y - 1);
    }); }

    [[nodiscard]] BitMatrix operator-() const { return reflected(); }

    [[nodiscard]] BitMatrix reflectedX() const { return mapped([](std::int64_t x, std::int64_t y) {
        return std::pair<std::int64_t, std::int64_t>(x, -y - 1);
    }); }

    [[nodiscard]] BitMatrix reflectedY() const { return mapped([](std::int64_t x, std::int64_t y) {
        return std::pair<std::int64_t, std::int64_t>(-x - 1, y);
    }); }

    [[nodiscard]] BitMatrix transposed() const { return mapped([](std::int64_t x, std::int64_t y) {
        return std::pair<std::int64_t, std::int64_t>(y, x);
    }); }

    [[nodiscard]] BitMatrix rotated90(int k = 1) const {
        const int turns = ((k % 4) + 4) % 4;
        return mapped([turns](std::int64_t x, std::int64_t y) {
            switch (turns) {
                case 1: return std::pair<std::int64_t, std::int64_t>(-y - 1, x);
                case 2: return std::pair<std::int64_t, std::int64_t>(-x - 1, -y - 1);
                case 3: return std::pair<std::int64_t, std::int64_t>(y, -x - 1);
                default: return std::pair<std::int64_t, std::int64_t>(x, y);
            }
        });
    }

    BitMatrix& rotate90(int k = 1) {
        *this = rotated90(k);
        return *this;
    }

    [[nodiscard]] BitMatrix latticeReflected() const { return mapped([](std::int64_t x, std::int64_t y) {
        return std::pair<std::int64_t, std::int64_t>(-x, -y);
    }); }

    [[nodiscard]] BitMatrix latticeReflectedX() const { return mapped([](std::int64_t x, std::int64_t y) {
        return std::pair<std::int64_t, std::int64_t>(x, -y);
    }); }

    [[nodiscard]] BitMatrix latticeReflectedY() const { return mapped([](std::int64_t x, std::int64_t y) {
        return std::pair<std::int64_t, std::int64_t>(-x, y);
    }); }

    [[nodiscard]] BitMatrix latticeTransposed() const { return transposed(); }

    [[nodiscard]] BitMatrix latticeRotated90(int k = 1) const {
        const int turns = ((k % 4) + 4) % 4;
        return mapped([turns](std::int64_t x, std::int64_t y) {
            switch (turns) {
                case 1: return std::pair<std::int64_t, std::int64_t>(-y, x);
                case 2: return std::pair<std::int64_t, std::int64_t>(-x, -y);
                case 3: return std::pair<std::int64_t, std::int64_t>(y, -x);
                default: return std::pair<std::int64_t, std::int64_t>(x, y);
            }
        });
    }

    BitMatrix& latticeRotate90(int k = 1) {
        *this = latticeRotated90(k);
        return *this;
    }

    [[nodiscard]] BitMatrix latticeMinkowskiSum(const BitMatrix& other) const {
        const BitMatrix left = trimmed(), right = other.trimmed();
        if (left.emptyWindow() || right.emptyWindow()) {
            return BitMatrix();
        }
 
        const bool stampIsRight = right.count() < left.count();
        const BitMatrix& stamp = stampIsRight ? right : left;
        const BitMatrix& canvas = stampIsRight ? left : right;
 
        BitMatrix result(PointType(left.origin_.x() + right.origin_.x(),
                                   left.origin_.y() + right.origin_.y()),
                         left.width_ + right.width_ - 1, left.height_ + right.height_ - 1);
        for (const PointType& cell : stamp) {
            result.orShifted(canvas, static_cast<std::int64_t>(cell.x() - stamp.origin_.x()),
                             static_cast<std::int64_t>(cell.y() - stamp.origin_.y()));
        }
        return result;
    }

    [[nodiscard]] BitMatrix minkowskiSum(const BitMatrix& other) const {
        if (empty() || other.empty()) {
            return BitMatrix();
        }
        return latticeMinkowskiSum(other).latticeMinkowskiSum(unitSquareSum());
    }

    [[nodiscard]] BitMatrix operator+(const BitMatrix& other) const { return minkowskiSum(other); }

    [[nodiscard]] BitMatrix minkowskiErosion(const BitMatrix& other) const {
        return latticeMinkowskiErosion(other.latticeMinkowskiSum(unitSquareSum()));
    }

    [[nodiscard]] BitMatrix latticeMinkowskiErosion(const BitMatrix& other) const {
        const BitMatrix stamp = other.trimmed();
        if (stamp.emptyWindow()) {
            BitMatrix result(origin_, width_, height_);
            result.setAll();
            return result;
        }
        const int resultWidth = width_ - stamp.width_ + 1;
        const int resultHeight = height_ - stamp.height_ + 1;
        BitMatrix result(PointType(origin_.x() - stamp.origin_.x(), origin_.y() - stamp.origin_.y()),
                         resultWidth, resultHeight);
        if (result.emptyWindow()) {
            return result;
        }
        result.setAll();
        for (const PointType& cell : stamp) {
            result.andShifted(*this, static_cast<std::int64_t>(cell.x() - stamp.origin_.x()),
                              static_cast<std::int64_t>(cell.y() - stamp.origin_.y()));
        }
        return result;
    }

    [[nodiscard]] BitMatrix latticeOpening(const BitMatrix& other) const {
        return latticeMinkowskiErosion(other).latticeMinkowskiSum(other);
    }

    [[nodiscard]] BitMatrix latticeClosing(const BitMatrix& other) const {
        return latticeMinkowskiSum(other).latticeMinkowskiErosion(other);
    }

    [[nodiscard]] BitMatrix interior(GridAdjacency adjacency = GridAdjacency::edge) const {
        BitMatrix result = *this;
        for (const auto& [dx, dy] : neighborOffsets(adjacency)) {
            result.andShifted(*this, dx, dy);
        }
        return result;
    }

    [[nodiscard]] BitMatrix boundary(GridAdjacency adjacency = GridAdjacency::edge) const {
        return difference(interior(adjacency));
    }
 
    // -----------------------------------------------------------------------
    // Connectivity

    [[nodiscard]] std::vector<BitMatrix> connectedComponents(GridAdjacency adjacency = GridAdjacency::edge) const {
        std::vector<BitMatrix> result;
        visitComponentRuns(adjacency, [&](const std::vector<CellRun>& runs) {
            std::int64_t minX = runs.front().x0, maxX = runs.front().x1 - 1;
            int minY = runs.front().y, maxY = minY;
            for (const CellRun& run : runs) {
                minX = std::min(minX, run.x0);
                maxX = std::max(maxX, run.x1 - 1);
                minY = std::min(minY, run.y);
                maxY = std::max(maxY, run.y);
            }
            BitMatrix component(PointType(origin_.x() + static_cast<NumberType>(minX),
                                          origin_.y() + static_cast<NumberType>(minY)),
                                static_cast<int>(maxX - minX) + 1, maxY - minY + 1);
            for (const CellRun& run : runs) {
                component.setLocalRange(run.y - minY, run.x0 - minX, run.x1 - minX);
            }
            result.push_back(std::move(component));
        });
        return result;
    }

    [[nodiscard]] std::size_t componentCount(GridAdjacency adjacency = GridAdjacency::edge) const {
        std::size_t total = 0;
        visitComponentRuns(adjacency, [&](const std::vector<CellRun>&) { ++total; });
        return total;
    }

    [[nodiscard]] bool isConnected(GridAdjacency adjacency = GridAdjacency::edge) const {
        return !empty() && componentCount(adjacency) == 1;
    }

    [[nodiscard]] BitMatrix fillHoles(GridAdjacency adjacency = GridAdjacency::edge) const {
        if (empty()) {
            return *this;
        }
        const RectangleType box = bbox();
        const RectangleType padded(PointType(box.min().x() - NumberType(1), box.min().y() - NumberType(1)),
                                   PointType(box.max().x() + NumberType(1), box.max().y() + NumberType(1)));
        const BitMatrix outside =
            (~resized(padded)).floodFrom(padded.min(), complementary(adjacency));
        return (~outside).resized(window());
    }

    [[nodiscard]] std::size_t holeCount(GridAdjacency adjacency = GridAdjacency::edge) const {
        return static_cast<std::size_t>(static_cast<std::int64_t>(componentCount(adjacency)) -
                                        eulerNumber(adjacency));
    }

    [[nodiscard]] std::int64_t eulerNumber(GridAdjacency adjacency = GridAdjacency::edge) const {
        if (emptyWindow()) {
            return 0;
        }
        const std::size_t blockWords = static_cast<std::size_t>((width_ + 1 + 63) / 64);
        const int tail = (width_ + 1) % 64;
        const std::uint64_t lastMask = tail == 0 ? ~std::uint64_t(0) : (std::uint64_t(1) << tail) - 1;
        const std::vector<std::uint64_t> zeros(words_, 0);
 
        std::int64_t ones = 0, threes = 0, diagonals = 0;
        for (int j = -1; j < height_; ++j) {
            const std::uint64_t* lower = j >= 0 ? row(j) : zeros.data();
            const std::uint64_t* upper = j + 1 < height_ ? row(j + 1) : zeros.data();
            for (std::size_t w = 0; w < blockWords; ++w) {
                const std::int64_t base = static_cast<std::int64_t>(w) * 64 - 1;
                const std::uint64_t a = shiftedWord(lower, words_, base);
                const std::uint64_t b = shiftedWord(lower, words_, base + 1);
                const std::uint64_t c = shiftedWord(upper, words_, base);
                const std::uint64_t d = shiftedWord(upper, words_, base + 1);
                const std::uint64_t valid = w + 1 == blockWords ? lastMask : ~std::uint64_t(0);
                const std::uint64_t one = (a & ~b & ~c & ~d) | (~a & b & ~c & ~d) |
                                          (~a & ~b & c & ~d) | (~a & ~b & ~c & d);
                const std::uint64_t three = (~a & b & c & d) | (a & ~b & c & d) |
                                            (a & b & ~c & d) | (a & b & c & ~d);
                const std::uint64_t diagonal = (a & ~b & ~c & d) | (~a & b & c & ~d);
                ones += std::popcount(one & valid);
                threes += std::popcount(three & valid);
                diagonals += std::popcount(diagonal & valid);
            }
        }
        const std::int64_t sign = adjacency == GridAdjacency::edge ? 2 : -2;
        return (ones - threes + sign * diagonals) / 4;
    }
 
    // -----------------------------------------------------------------------
    // Convexity

    bool fillRows() {
        bool changed = false;
        for (int j = 0; j < height_; ++j) {
            if (const std::optional<std::pair<std::int64_t, std::int64_t>> extent = rowExtent(j)) {
                changed |= setLocalRange(j, extent->first, extent->second + 1);
            }
        }
        return changed;
    }

    bool fillColumns() {
        if (emptyWindow()) {
            return false;
        }
 
        // below[j] holds the columns with a set cell in some row under j.
        std::vector<std::uint64_t> below(bits_.size(), 0);
        for (int j = 1; j < height_; ++j) {
            for (std::size_t w = 0; w < words_; ++w) {
                below[static_cast<std::size_t>(j) * words_ + w] =
                    below[static_cast<std::size_t>(j - 1) * words_ + w] | row(j - 1)[w];
            }
        }
 
        bool changed = false;
        std::vector<std::uint64_t> above(words_, 0);  // Columns set in some row over j.
        for (int j = height_ - 1; j >= 0; --j) {
            std::uint64_t* here = row(j);
            for (std::size_t w = 0; w < words_; ++w) {
                const std::uint64_t original = here[w];
                const std::uint64_t fill =
                    above[w] & below[static_cast<std::size_t>(j) * words_ + w] & ~original;
                here[w] |= fill;
                above[w] |= original;
                changed |= fill != 0;
            }
        }
        return changed;
    }

    std::size_t makeHvConvex() {
        const std::size_t before = count();
        // Both fills run on every pass: filling the columns can open a row gap
        // and the other way round, so neither may be short-circuited away.
        for (bool changed = true; changed;) {
            const bool rows = fillRows();
            const bool columns = fillColumns();
            changed = rows || columns;
        }
        return count() - before;
    }

    [[nodiscard]] bool isRowConvex() const {
        for (int j = 0; j < height_; ++j) {
            const std::optional<std::pair<std::int64_t, std::int64_t>> extent = rowExtent(j);
            if (extent && static_cast<std::int64_t>(countRow(j)) !=
                              extent->second - extent->first + 1) {
                return false;
            }
        }
        return true;
    }

    [[nodiscard]] bool isColumnConvex() const {
        std::vector<std::uint64_t> seen(words_, 0), closed(words_, 0);
        for (int j = 0; j < height_; ++j) {
            const std::uint64_t* here = row(j);
            for (std::size_t w = 0; w < words_; ++w) {
                if ((here[w] & closed[w]) != 0) {
                    return false;
                }
                closed[w] |= seen[w] & ~here[w];
                seen[w] |= here[w];
            }
        }
        return true;
    }

    [[nodiscard]] bool isHvConvex() const { return isRowConvex() && isColumnConvex(); }
 
    // -----------------------------------------------------------------------

    class Iterator {
    public:
        using iterator_category = std::forward_iterator_tag;
        using value_type = PointType;
        using difference_type = std::ptrdiff_t;
        using pointer = const PointType*;
        using reference = const PointType&;
 
        Iterator() = default;

        reference operator*() const { return current_; }

        pointer operator->() const { return &current_; }

        Iterator& operator++() {
            advance();
            return *this;
        }

        Iterator operator++(int) {
            Iterator previous = *this;
            advance();
            return previous;
        }

        bool operator==(const Iterator& other) const {
            return owner_ == other.owner_ && (owner_ == nullptr || (rowIndex_ == other.rowIndex_ &&
                                                                    wordIndex_ == other.wordIndex_ &&
                                                                    rest_ == other.rest_));
        }

        bool operator!=(const Iterator& other) const { return !(*this == other); }
 
    private:
        friend class BitMatrix;
 
        explicit Iterator(const BitMatrix* owner) : owner_(owner) {
            if (owner_->words_ == 0) {
                owner_ = nullptr;
                return;
            }
            rest_ = owner_->row(0)[0];
            advance();
        }
 
        void advance() {
            while (true) {
                if (rest_ != 0) {
                    const int bit = std::countr_zero(rest_);
                    rest_ &= rest_ - 1;
                    current_ = PointType(
                        owner_->origin_.x() +
                            static_cast<NumberType>(static_cast<std::int64_t>(wordIndex_) * 64 + bit),
                        owner_->origin_.y() + static_cast<NumberType>(rowIndex_));
                    return;
                }
                ++wordIndex_;
                if (wordIndex_ >= owner_->words_) {
                    wordIndex_ = 0;
                    ++rowIndex_;
                }
                if (rowIndex_ >= owner_->height_) {
                    owner_ = nullptr;
                    return;
                }
                rest_ = owner_->row(rowIndex_)[wordIndex_];
            }
        }
 
        const BitMatrix* owner_ = nullptr;
        int rowIndex_ = 0;
        std::size_t wordIndex_ = 0;
        std::uint64_t rest_ = 0;
        PointType current_{};
    };
 
private:
    friend struct std::hash<BitMatrix>;

    struct CellRun {
        int y;            
        std::int64_t x0;  
        std::int64_t x1;  
    };

    static constexpr std::array<std::pair<int, int>, 4> boundarySteps{
        {{1, 0}, {0, 1}, {-1, 0}, {0, -1}}};

    [[nodiscard]] static RectangleType cellSquare(const PointType& cell) {
        return RectangleType(cell, PointType(cell.x() + NumberType(1), cell.y() + NumberType(1)),
                             true);
    }
 
    static std::size_t wordsPerRow(int width) {
        return static_cast<std::size_t>((width + 63) / 64);
    }

    static int cellSpan(NumberType low, NumberType high) {
        // Unsigned subtraction is exact here even where the difference overflows
        // the coordinate type, since high is at least low.
        const std::uint64_t span =
            static_cast<std::uint64_t>(high) - static_cast<std::uint64_t>(low);
        if (span >= static_cast<std::uint64_t>(detail::numeric_limits<int>::max())) {
            throw std::logic_error("pgl::BitMatrix: the points do not fit a window");
        }
        return static_cast<int>(span) + 1;
    }

    template <class Range>
    static BitMatrix emptyOver(Range&& points) {
        auto it = std::ranges::begin(points);
        const auto last = std::ranges::end(points);
        if (it == last) {
            return BitMatrix();
        }
        NumberType minX = (*it).x(), maxX = minX;
        NumberType minY = (*it).y(), maxY = minY;
        for (++it; it != last; ++it) {
            minX = std::min<NumberType>(minX, (*it).x());
            maxX = std::max<NumberType>(maxX, (*it).x());
            minY = std::min<NumberType>(minY, (*it).y());
            maxY = std::max<NumberType>(maxY, (*it).y());
        }
        return BitMatrix(PointType(minX, minY), cellSpan(minX, maxX), cellSpan(minY, maxY));
    }
 
    std::uint64_t* row(int j) { return bits_.data() + static_cast<std::size_t>(j) * words_; }
    const std::uint64_t* row(int j) const {
        return bits_.data() + static_cast<std::size_t>(j) * words_;
    }
 
    [[nodiscard]] std::int64_t localX(NumberType x) const {
        return static_cast<std::int64_t>(x) - static_cast<std::int64_t>(origin_.x());
    }
    [[nodiscard]] std::int64_t localY(NumberType y) const {
        return static_cast<std::int64_t>(y) - static_cast<std::int64_t>(origin_.y());
    }
 
    [[nodiscard]] bool localGet(int j, std::int64_t i) const {
        return ((row(j)[static_cast<std::size_t>(i) / 64] >> (i % 64)) & 1) != 0;
    }
 
    [[nodiscard]] std::size_t countRow(int j) const {
        std::size_t total = 0;
        for (std::size_t w = 0; w < words_; ++w) {
            total += static_cast<std::size_t>(std::popcount(row(j)[w]));
        }
        return total;
    }

    [[nodiscard]] std::optional<std::pair<std::int64_t, std::int64_t>> rowExtent(int j) const {
        const std::uint64_t* here = row(j);
        std::size_t first = 0;
        while (first < words_ && here[first] == 0) {
            ++first;
        }
        if (first == words_) {
            return std::nullopt;
        }
        std::size_t last = words_ - 1;
        while (here[last] == 0) {
            --last;
        }
        return std::pair{static_cast<std::int64_t>(first) * 64 + std::countr_zero(here[first]),
                         static_cast<std::int64_t>(last) * 64 + 63 - std::countl_zero(here[last])};
    }

    using Crossings = std::vector<std::vector<NumberType>>;

    template <class Ring>
    void addRingCrossings(const Ring& ring, Crossings& crossings) const {
        for (std::size_t i = 0; i < ring.size(); ++i) {
            const PointType p = ring[i];
            const PointType q = ring[(i + 1) % ring.size()];
            if (p.x() == q.x()) {
                for (NumberType y = std::min(p.y(), q.y()); y < std::max(p.y(), q.y()); ++y) {
                    crossings[static_cast<std::size_t>(y - origin_.y())].push_back(p.x());
                }
            } else if (p.y() != q.y()) {
                throw std::logic_error("pgl::BitMatrix: the shape is not rectilinear");
            }
        }
    }

    template <class Region>
    void addRegionCrossings(const Region& region, Crossings& crossings) const {
        addRingCrossings(region.outer(), crossings);
        for (const Polygon<PointType>& hole : region.holes()) {
            addRingCrossings(hole, crossings);
        }
    }

    void fillCrossings(Crossings& crossings) {
        for (int j = 0; j < height_; ++j) {
            std::vector<NumberType>& crossingsInRow = crossings[static_cast<std::size_t>(j)];
            std::sort(crossingsInRow.begin(), crossingsInRow.end());
            for (std::size_t i = 0; i + 1 < crossingsInRow.size(); i += 2) {
                setLocalRange(j, static_cast<std::int64_t>(crossingsInRow[i] - origin_.x()),
                              static_cast<std::int64_t>(crossingsInRow[i + 1] - origin_.x()));
            }
        }
    }

    bool setLocalRange(int j, std::int64_t low, std::int64_t high) {
        low = std::max<std::int64_t>(low, 0);
        high = std::min<std::int64_t>(high, width_);
        bool changed = false;
        for (std::int64_t i = low; i < high;) {
            const std::size_t w = static_cast<std::size_t>(i) / 64;
            const int offset = static_cast<int>(i % 64);
            const int bits = static_cast<int>(std::min<std::int64_t>(64 - offset, high - i));
            const std::uint64_t mask =
                bits == 64 ? ~std::uint64_t(0) : ((std::uint64_t(1) << bits) - 1) << offset;
            changed |= (row(j)[w] & mask) != mask;
            row(j)[w] |= mask;
            i += bits;
        }
        return changed;
    }

    static void clearWordRange(std::uint64_t* here, std::int64_t low, std::int64_t high) {
        for (std::int64_t i = low; i < high;) {
            const std::size_t w = static_cast<std::size_t>(i) / 64;
            const int offset = static_cast<int>(i % 64);
            const int bits = static_cast<int>(std::min<std::int64_t>(64 - offset, high - i));
            const std::uint64_t mask =
                bits == 64 ? ~std::uint64_t(0) : ((std::uint64_t(1) << bits) - 1) << offset;
            here[w] &= ~mask;
            i += bits;
        }
    }
 
    // The three searches below read a row of `words_` words holding `width_`
    // cells, the bits past the width being clear, and skip whole words. They
    // take the row rather than its index so a flood fill can run them over its
    // own copy of the bits.

    [[nodiscard]] std::int64_t nextSetInRow(const std::uint64_t* here, std::int64_t from) const {
        assert(from >= 0 && from <= width_);
        if (from >= width_) {
            return width_;
        }
        std::size_t w = static_cast<std::size_t>(from) / 64;
        std::uint64_t rest = here[w] & (~std::uint64_t(0) << (from % 64));
        while (rest == 0) {
            if (++w >= words_) {
                return width_;
            }
            rest = here[w];
        }
        return static_cast<std::int64_t>(w) * 64 + std::countr_zero(rest);
    }

    [[nodiscard]] std::int64_t nextClearInRow(const std::uint64_t* here, std::int64_t from) const {
        assert(from >= 0 && from < width_);
        std::size_t w = static_cast<std::size_t>(from) / 64;
        std::uint64_t rest = ~here[w] & (~std::uint64_t(0) << (from % 64));
        while (rest == 0) {
            if (++w >= words_) {
                return width_;
            }
            rest = ~here[w];
        }
        return std::min<std::int64_t>(
            static_cast<std::int64_t>(w) * 64 + std::countr_zero(rest), width_);
    }

    [[nodiscard]] std::int64_t runStartInRow(const std::uint64_t* here, std::int64_t at) const {
        assert(at >= 0 && at < width_ && ((here[static_cast<std::size_t>(at) / 64] >>
                                           (at % 64)) & 1) != 0);
        std::size_t w = static_cast<std::size_t>(at) / 64;
        const int bit = static_cast<int>(at % 64);
        std::uint64_t rest = ~here[w] & (~std::uint64_t(0) >> (63 - bit));
        while (rest == 0) {
            if (w == 0) {
                return 0;
            }
            rest = ~here[--w];
        }
        return static_cast<std::int64_t>(w) * 64 + 64 - std::countl_zero(rest);
    }

    void maskTails() {
        const int tail = width_ % 64;
        if (emptyWindow() || tail == 0) {
            return;
        }
        const std::uint64_t mask = (std::uint64_t(1) << tail) - 1;
        for (int j = 0; j < height_; ++j) {
            row(j)[words_ - 1] &= mask;
        }
    }
 
    static std::uint64_t wordAt(const std::uint64_t* words, std::size_t count, std::int64_t index) {
        return index < 0 || index >= static_cast<std::int64_t>(count) ? 0
                                                                     : words[static_cast<std::size_t>(index)];
    }

    static std::uint64_t shiftedWord(const std::uint64_t* words, std::size_t count,
                                     std::int64_t offset) {
        const std::int64_t index = offset >= 0 ? offset / 64 : -((-offset + 63) / 64);
        const int rest = static_cast<int>(offset - index * 64);
        const std::uint64_t low = wordAt(words, count, index) >> rest;
        if (rest == 0) {
            return low;
        }
        return low | (wordAt(words, count, index + 1) << (64 - rest));
    }

    template <class WordOp>
    void combine(const BitMatrix& other, WordOp op) {
        const std::int64_t shift =
            static_cast<std::int64_t>(other.origin_.x()) - static_cast<std::int64_t>(origin_.x());
        for (int j = 0; j < height_; ++j) {
            const std::int64_t sourceRow = static_cast<std::int64_t>(origin_.y()) + j -
                                           static_cast<std::int64_t>(other.origin_.y());
            const bool inside = sourceRow >= 0 && sourceRow < other.height_;
            const std::uint64_t* source = inside ? other.row(static_cast<int>(sourceRow)) : nullptr;
            std::uint64_t* destination = row(j);
            for (std::size_t w = 0; w < words_; ++w) {
                const std::uint64_t word =
                    source == nullptr
                        ? 0
                        : shiftedWord(source, other.words_, static_cast<std::int64_t>(w) * 64 - shift);
                destination[w] = op(destination[w], word);
            }
        }
        maskTails();
    }

    template <class Fn>
    void forEachAlignedWord(const BitMatrix& other, Fn fn) const {
        const std::int64_t shift =
            static_cast<std::int64_t>(other.origin_.x()) - static_cast<std::int64_t>(origin_.x());
        for (int j = 0; j < height_; ++j) {
            const std::int64_t sourceRow = static_cast<std::int64_t>(origin_.y()) + j -
                                           static_cast<std::int64_t>(other.origin_.y());
            if (sourceRow < 0 || sourceRow >= other.height_) {
                continue;
            }
            const std::uint64_t* source = other.row(static_cast<int>(sourceRow));
            const std::uint64_t* here = row(j);
            for (std::size_t w = 0; w < words_; ++w) {
                if (here[w] == 0) {
                    continue;
                }
                fn(here[w], shiftedWord(source, other.words_, static_cast<std::int64_t>(w) * 64 - shift));
            }
        }
    }

    [[nodiscard]] bool anyCommonCell(const BitMatrix& other, std::int64_t dx, std::int64_t dy) const {
        const std::int64_t shift = static_cast<std::int64_t>(other.origin_.x()) + dx -
                                   static_cast<std::int64_t>(origin_.x());
        for (int j = 0; j < height_; ++j) {
            const std::int64_t sourceRow = static_cast<std::int64_t>(origin_.y()) + j -
                                           static_cast<std::int64_t>(other.origin_.y()) - dy;
            if (sourceRow < 0 || sourceRow >= other.height_) {
                continue;
            }
            const std::uint64_t* source = other.row(static_cast<int>(sourceRow));
            const std::uint64_t* here = row(j);
            for (std::size_t w = 0; w < words_; ++w) {
                if (here[w] == 0) {
                    continue;
                }
                if ((here[w] & shiftedWord(source, other.words_,
                                           static_cast<std::int64_t>(w) * 64 - shift)) != 0) {
                    return true;
                }
            }
        }
        return false;
    }

    void orShifted(const BitMatrix& source, std::int64_t dx, std::int64_t dy) {
        assert(dx >= 0 && dy >= 0);
        assert(dy + source.height_ <= height_);
        const std::size_t wordShift = static_cast<std::size_t>(dx / 64);
        const int bitShift = static_cast<int>(dx % 64);
        for (int j = 0; j < source.height_; ++j) {
            const std::uint64_t* from = source.row(j);
            std::uint64_t* to = row(j + static_cast<int>(dy));
            for (std::size_t w = 0; w < source.words_; ++w) {
                const std::uint64_t word = from[w];
                if (word == 0) {
                    continue;
                }
                if (w + wordShift < words_) {
                    to[w + wordShift] |= word << bitShift;
                }
                if (bitShift != 0 && w + wordShift + 1 < words_) {
                    to[w + wordShift + 1] |= word >> (64 - bitShift);
                }
            }
        }
    }

    void andShifted(const BitMatrix& source, std::int64_t dx, std::int64_t dy) {
        for (int j = 0; j < height_; ++j) {
            const std::int64_t sourceRow = static_cast<std::int64_t>(j) + dy;
            std::uint64_t* to = row(j);
            if (sourceRow < 0 || sourceRow >= source.height_) {
                std::fill(to, to + words_, std::uint64_t(0));
                continue;
            }
            const std::uint64_t* from = source.row(static_cast<int>(sourceRow));
            for (std::size_t w = 0; w < words_; ++w) {
                to[w] &= shiftedWord(from, source.words_, static_cast<std::int64_t>(w) * 64 + dx);
            }
        }
    }

    template <class Map>
    [[nodiscard]] BitMatrix mapped(Map map) const {
        if (emptyWindow()) {
            return BitMatrix();
        }
        const std::int64_t x0 = static_cast<std::int64_t>(origin_.x());
        const std::int64_t y0 = static_cast<std::int64_t>(origin_.y());
        const std::int64_t x1 = x0 + width_ - 1, y1 = y0 + height_ - 1;
        std::int64_t minX = 0, maxX = 0, minY = 0, maxY = 0;
        bool first = true;
        for (const std::int64_t x : {x0, x1}) {
            for (const std::int64_t y : {y0, y1}) {
                const auto [imageX, imageY] = map(x, y);
                minX = first ? imageX : std::min(minX, imageX);
                maxX = first ? imageX : std::max(maxX, imageX);
                minY = first ? imageY : std::min(minY, imageY);
                maxY = first ? imageY : std::max(maxY, imageY);
                first = false;
            }
        }
        BitMatrix result(PointType(static_cast<NumberType>(minX), static_cast<NumberType>(minY)),
                         static_cast<int>(maxX - minX) + 1, static_cast<int>(maxY - minY) + 1);
        for (const PointType& cell : *this) {
            const auto [imageX, imageY] =
                map(static_cast<std::int64_t>(cell.x()), static_cast<std::int64_t>(cell.y()));
            result.set(static_cast<NumberType>(imageX), static_cast<NumberType>(imageY));
        }
        return result;
    }

    static constexpr std::array<std::pair<std::int64_t, std::int64_t>, 8> allNeighborOffsets{
        {{1, 0}, {-1, 0}, {0, 1}, {0, -1}, {1, 1}, {1, -1}, {-1, 1}, {-1, -1}}};
 
    static std::span<const std::pair<std::int64_t, std::int64_t>> neighborOffsets(
        GridAdjacency adjacency) {
        return std::span(allNeighborOffsets)
            .first(adjacency == GridAdjacency::vertex ? 8u : 4u);
    }

    [[nodiscard]] static BitMatrix unitSquareSum() {
        BitMatrix stamp(PointType(), 2, 2);
        stamp.setAll();
        return stamp;
    }
 
    static GridAdjacency complementary(GridAdjacency adjacency) {
        return adjacency == GridAdjacency::edge ? GridAdjacency::vertex : GridAdjacency::edge;
    }

    [[nodiscard]] std::optional<PointType> firstSetCell() const {
        for (int j = 0; j < height_; ++j) {
            const std::uint64_t* here = row(j);
            for (std::size_t w = 0; w < words_; ++w) {
                if (here[w] != 0) {
                    return PointType(
                        origin_.x() + static_cast<NumberType>(static_cast<std::int64_t>(w) * 64 +
                                                              std::countr_zero(here[w])),
                        origin_.y() + static_cast<NumberType>(j));
                }
            }
        }
        return std::nullopt;
    }

    void floodRuns(std::uint64_t* remaining, int j, std::int64_t i, std::int64_t reach,
                   std::vector<CellRun>& runs) const {
        const std::size_t first = runs.size();
        auto rowOf = [&](int at) { return remaining + static_cast<std::size_t>(at) * words_; };
        auto take = [&](int at, std::int64_t cell) {
            std::uint64_t* here = rowOf(at);
            const CellRun run{at, runStartInRow(here, cell), nextClearInRow(here, cell)};
            clearWordRange(here, run.x0, run.x1);
            runs.push_back(run);
            return run;
        };
        auto scan = [&](int at, std::int64_t low, std::int64_t high) {
            if (at < 0 || at >= height_) {
                return;
            }
            const std::uint64_t* here = rowOf(at);
            high = std::min<std::int64_t>(high, width_);
            for (std::int64_t cell = nextSetInRow(here, std::max<std::int64_t>(low, 0)); cell < high;
                 cell = nextSetInRow(here, cell)) {
                cell = take(at, cell).x1;
            }
        };
 
        take(j, i);
        for (std::size_t next = first; next < runs.size(); ++next) {
            const CellRun run = runs[next];  // Copied: taking a run may reallocate.
            scan(run.y - 1, run.x0 - reach, run.x1 + reach);
            scan(run.y + 1, run.x0 - reach, run.x1 + reach);
        }
    }

    [[nodiscard]] BitMatrix floodFrom(const PointType& seed, GridAdjacency adjacency) const {
        BitMatrix result(origin_, width_, height_);
        if (!get(seed)) {
            return result;
        }
        std::vector<std::uint64_t> remaining = bits_;
        std::vector<CellRun> runs;
        floodRuns(remaining.data(), static_cast<int>(localY(seed.y())), localX(seed.x()),
                  adjacency == GridAdjacency::vertex ? 1 : 0, runs);
        for (const CellRun& run : runs) {
            result.setLocalRange(run.y, run.x0, run.x1);
        }
        return result;
    }

    template <class Fn>
    void visitComponentRuns(GridAdjacency adjacency, Fn fn) const {
        if (emptyWindow()) {
            return;
        }
        std::vector<std::uint64_t> remaining = bits_;
        // A diagonal neighbor of a run sits one cell past either of its ends.
        const std::int64_t reach = adjacency == GridAdjacency::vertex ? 1 : 0;
        std::vector<CellRun> runs;
        for (int j = 0; j < height_; ++j) {
            const std::uint64_t* here = remaining.data() + static_cast<std::size_t>(j) * words_;
            for (std::int64_t i = nextSetInRow(here, 0); i < width_; i = nextSetInRow(here, i)) {
                runs.clear();
                floodRuns(remaining.data(), j, i, reach, runs);
                fn(std::as_const(runs));
            }
        }
    }

    template <class PinchCrosses>
    [[nodiscard]] std::vector<std::vector<detail::PolyCell>> boundaryLoops(
        PinchCrosses crosses) const {
        std::vector<std::vector<detail::PolyCell>> loops;
        if (emptyWindow()) {
            return loops;
        }
        const std::size_t vertexWords = wordsPerRow(width_ + 1);
        const std::size_t vertexRows = static_cast<std::size_t>(height_) + 1;
        const std::size_t plane = vertexRows * vertexWords;
        std::vector<std::uint64_t> present(4 * plane, 0);
        std::vector<std::uint64_t> used(4 * plane, 0);
 
        auto cellWord = [&](int j, std::size_t w) -> std::uint64_t {
            return j < 0 || j >= height_ || w >= words_ ? 0 : row(j)[w];
        };
        // The north and west edges of a vertex row come from the cells one to
        // the left, so they are built raw and then shifted across the words.
        std::vector<std::uint64_t> rawNorth(vertexWords), rawWest(vertexWords);
        for (std::size_t vertexRow = 0; vertexRow < vertexRows; ++vertexRow) {
            const int above = static_cast<int>(vertexRow);  // Cells resting on the row.
            const int below = above - 1;                    // Cells hanging under it.
            for (std::size_t w = 0; w < vertexWords; ++w) {
                const std::uint64_t upper = cellWord(above, w);
                const std::uint64_t lower = cellWord(below, w);
                const std::uint64_t upperRight = (upper >> 1) | (cellWord(above, w + 1) << 63);
                const std::uint64_t lowerLeft =
                    (lower << 1) | (w == 0 ? 0 : cellWord(below, w - 1) >> 63);
                present[0 * plane + vertexRow * vertexWords + w] = upper & ~lower;      // East.
                present[3 * plane + vertexRow * vertexWords + w] = lower & ~lowerLeft;  // South.
                rawNorth[w] = upper & ~upperRight;
                rawWest[w] = lower & ~upper;
            }
            std::uint64_t carryNorth = 0, carryWest = 0;
            for (std::size_t w = 0; w < vertexWords; ++w) {
                present[1 * plane + vertexRow * vertexWords + w] = (rawNorth[w] << 1) | carryNorth;
                present[2 * plane + vertexRow * vertexWords + w] = (rawWest[w] << 1) | carryWest;
                carryNorth = rawNorth[w] >> 63;
                carryWest = rawWest[w] >> 63;
            }
        }
 
        auto index = [&](int direction, int y, int x) {
            return static_cast<std::size_t>(direction) * plane +
                   static_cast<std::size_t>(y) * vertexWords + static_cast<std::size_t>(x) / 64;
        };
        auto hasEdge = [&](int direction, int y, int x) {
            return ((present[index(direction, y, x)] >> (x % 64)) & 1) != 0;
        };
        auto markUsed = [&](int direction, int y, int x) {
            used[index(direction, y, x)] |= std::uint64_t(1) << (x % 64);
        };
 
        auto trace = [&](int direction, int startX, int startY) {
            std::vector<detail::PolyCell> loop;
            int x = startX, y = startY;
            for (;;) {
                markUsed(direction, y, x);
                loop.emplace_back(x, y);
                x += boundarySteps[static_cast<std::size_t>(direction)].first;
                y += boundarySteps[static_cast<std::size_t>(direction)].second;
                if (x == startX && y == startY) {
                    break;
                }
                int next = (direction + 3) % 4;        // Sharpest right.
                const int left = (direction + 1) % 4;  // Sharpest left.
                if (hasEdge(next, y, x) && hasEdge(left, y, x)) {
                    // The cell just walked along is on the left of the edge that
                    // arrived here, which left the vertex one step back.
                    const auto& step = boundarySteps[static_cast<std::size_t>(direction)];
                    if (!crosses(leftCell(direction, x - step.first, y - step.second),
                                 leftCell(next, x, y))) {
                        next = left;
                    }
                } else {
                    // Preference order: sharpest right, straight, left, reverse.
                    for (int turn = 0; turn < 3 && !hasEdge(next, y, x); ++turn) {
                        next = (next + 1) % 4;
                    }
                }
                assert(hasEdge(next, y, x) && "pgl::BitMatrix: the boundary walk ran off an edge");
                direction = next;
            }
            return loop;
        };
 
        for (int direction = 0; direction < 4; ++direction) {
            for (int y = 0; y <= height_; ++y) {
                for (std::size_t w = 0; w < vertexWords; ++w) {
                    const std::size_t at =
                        static_cast<std::size_t>(direction) * plane +
                        static_cast<std::size_t>(y) * vertexWords + w;
                    while (const std::uint64_t rest = present[at] & ~used[at]) {
                        loops.push_back(trace(
                            direction, static_cast<int>(w) * 64 + std::countr_zero(rest), y));
                    }
                }
            }
        }
        return loops;
    }

    static detail::PolyCell leftCell(int direction, int x, int y) {
        switch (direction) {
            case 0:
                return {x, y};  // East, along the bottom edge of the cell.
            case 1:
                return {x - 1, y};  // North, along the right edge.
            case 2:
                return {x - 1, y - 1};  // West, along the top edge.
            default:
                return {x, y - 1};  // South, along the left edge.
        }
    }

    static detail::PolyCell loopSeedCell(const std::vector<detail::PolyCell>& loop) {
        const auto [x, y] = loop.front();
        const auto [nextX, nextY] = loop[1];
        return leftCell(nextX > x ? 0 : nextY > y ? 1 : nextX < x ? 2 : 3, x, y);
    }

    [[nodiscard]] RegionType regionFromLoops(
        const std::vector<std::vector<detail::PolyCell>>& loops) const {
        Polygon<Point<NumberType>> outer;
        std::vector<Polygon<Point<NumberType>>> holes;
        for (const std::vector<detail::PolyCell>& loop : loops) {
            if (detail::loopTwiceArea(loop) > 0) {
                outer = Polygon<Point<NumberType>>(detail::loopCorners<NumberType>(loop));
            } else {
                holes.emplace_back(detail::loopCorners<NumberType>(loop));
            }
        }
        return RegionType(
            Transformation<NumberType>::translation(origin_.x(), origin_.y()) *
            PolygonWithHoles<Point<NumberType>>(std::move(outer), std::move(holes)));
    }
 
    static RectangleType hullWindow(const BitMatrix& left, const BitMatrix& right) {
        if (left.emptyWindow()) {
            return right.window();
        }
        if (right.emptyWindow()) {
            return left.window();
        }
        const RectangleType a = left.window(), b = right.window();
        return RectangleType(PointType(std::min(a.min().x(), b.min().x()), std::min(a.min().y(), b.min().y())),
                             PointType(std::max(a.max().x(), b.max().x()), std::max(a.max().y(), b.max().y())),
                             true);
    }
 
    static RectangleType overlapWindow(const BitMatrix& left, const BitMatrix& right) {
        if (left.emptyWindow() || right.emptyWindow()) {
            return RectangleType();
        }
        const RectangleType a = left.window(), b = right.window();
        return RectangleType(PointType(std::max(a.min().x(), b.min().x()), std::max(a.min().y(), b.min().y())),
                             PointType(std::min(a.max().x(), b.max().x()), std::min(a.max().y(), b.max().y())),
                             true);
    }
 
    PointType origin_{};
    int width_ = 0;
    int height_ = 0;
    std::size_t words_ = 0;
    std::vector<std::uint64_t> bits_;
};

template <class PointType>
BitMatrix<PointType> operator+(const PointType& vector, const BitMatrix<PointType>& matrix) {
    return matrix.translated(vector);
}
 
// Written out rather than left to the implicit guides: the class template is
// constrained, and clang 18 mishandles a constraint carried onto an implicit
// guide.
template <class PointType>
BitMatrix(PointType, int, int) -> BitMatrix<PointType>;
 
template <class PointType, class LabelType>
BitMatrix(const Rectangle<PointType, LabelType>&) -> BitMatrix<PointType>;
 
template <class PointType, class LabelType>
BitMatrix(const PolygonWithHoles<PointType, LabelType>&) -> BitMatrix<PointType>;
 
template <class PointType, class LabelType>
BitMatrix(const Polygon<PointType, LabelType>&) -> BitMatrix<PointType>;
 
template <class PointType, class LabelType>
BitMatrix(const PolygonSet<PointType, LabelType>&) -> BitMatrix<PointType>;
 
// The point-range constructor deduces nothing on its own -- PointType does not
// appear in its signature -- so without this guide a range of points would take
// the default point type whatever it holds.
template <std::ranges::input_range Range>
    requires(detail::is_point_v<std::remove_cvref_t<std::ranges::range_value_t<Range>>>
             && !AnyShapeConcept<std::remove_cvref_t<Range>>
             && !detail::is_bit_matrix_v<std::remove_cvref_t<Range>>)
BitMatrix(Range&&) -> BitMatrix<std::remove_cvref_t<std::ranges::range_value_t<Range>>>;
 
// Out-of-line: asBitMatrix is declared in the shape headers (which precede this
// one in the layering) but can only be defined once BitMatrix is visible. Each
// is the rasterizing constructor for that shape, so it is the constructor that
// documents the window, the fill rule and the rectilinear requirement.
template <class PointType_, class TLabel>
template <class ResultNumber>
    requires(std::signed_integral<ResultNumber>)
auto Polygon<PointType_, TLabel>::asBitMatrix() const {
    return BitMatrix<detail::grid_point_t<PointType_, ResultNumber>>(*this);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber>
    requires(std::signed_integral<ResultNumber>)
auto PolygonWithHoles<PointType_, TLabel>::asBitMatrix() const {
    return BitMatrix<detail::grid_point_t<PointType_, ResultNumber>>(*this);
}
 
template <class PointType_, class TLabel>
template <class ResultNumber>
    requires(std::signed_integral<ResultNumber>)
auto PolygonSet<PointType_, TLabel>::asBitMatrix() const {
    return BitMatrix<detail::grid_point_t<PointType_, ResultNumber>>(*this);
}
 
namespace detail {

template <class PointType, class Predicate>
BitMatrix<PointType> rasterize(const Rectangle<PointType>& window, Predicate keep) {
    using Number = typename PointType::NumberType;
    BitMatrix<PointType> result(window);
    for (int j = 0; j < result.height(); ++j) {
        for (int i = 0; i < result.width(); ++i) {
            const Number x = result.origin().x() + static_cast<Number>(i);
            const Number y = result.origin().y() + static_cast<Number>(j);
            if (keep(Rectangle<PointType>(PointType(x, y),
                                          PointType(x + Number(1), y + Number(1)), true))) {
                result.set(x, y);
            }
        }
    }
    return result;
}
 
}  // namespace detail

template <class PointType, class ShapeType>
BitMatrix<PointType> outerRaster(const ShapeType& shape, const Rectangle<PointType>& window) {
    return detail::rasterize<PointType>(
        window, [&shape](const Rectangle<PointType>& cell) { return shape.intersects(cell); });
}

template <class PointType, class ShapeType>
BitMatrix<PointType> innerRaster(const ShapeType& shape, const Rectangle<PointType>& window) {
    return detail::rasterize<PointType>(
        window, [&shape](const Rectangle<PointType>& cell) { return shape.contains(cell); });
}

template <class PointType = Point<int>, class ShapeType>
    requires std::signed_integral<
        std::remove_cvref_t<decltype(std::declval<const ShapeType&>().bbox().min().x())>>
BitMatrix<PointType> outerRaster(const ShapeType& shape) {
    return outerRaster<PointType>(shape, Rectangle<PointType>(shape.bbox()));
}

template <class PointType = Point<int>, class ShapeType>
    requires std::signed_integral<
        std::remove_cvref_t<decltype(std::declval<const ShapeType&>().bbox().min().x())>>
BitMatrix<PointType> innerRaster(const ShapeType& shape) {
    return innerRaster<PointType>(shape, Rectangle<PointType>(shape.bbox()));
}
 
}  // namespace pgl
 
namespace std {

template <class PointType>
struct hash<pgl::BitMatrix<PointType>> {
    std::size_t operator()(const pgl::BitMatrix<PointType>& matrix) const {
        std::size_t seed = 0;
        pgl::detail::hashCombine(seed, matrix.origin());
        pgl::detail::hashCombine(seed, matrix.width());
        pgl::detail::hashCombine(seed, matrix.height());
        for (const std::uint64_t word : matrix.bits_) {
            pgl::detail::hashCombine(seed, word);
        }
        return seed;
    }
};
 
}  // namespace std