MorphologicalTree.hpp

#pragma once
 
#include "../utils/RegularGridAdjacency2D.hpp"
#include "../utils/Assert.hpp"
#include "../utils/Altitude.hpp"
#include "../utils/Common.hpp"
#include "../utils/Contract.hpp"
#include "../dataStructure/FastQueue.hpp"
#include "MorphologicalTreeSemantics.hpp"
#include "ProperPartDomain.hpp"
#include "detail/MorphologicalTreeConstructionTag.hpp"
#include "detail/NativeHierarchyValidationDetail.hpp"
#include <algorithm>
#include <cassert>
#include <cstddef>
#include <cstdint>
#include <functional>
#include <limits>
#include <list>
#include <memory>
#include <optional>
#include <set>
#include <span>
#include <stdexcept>
#include <string>
#include <tuple>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>
 
namespace mmcfilters {
 
namespace detail {
class CommittedTreeAccess;
}
 
enum class NodeIdSpace { MorphologicalTree, Higra };
 
struct [[nodiscard]] TreeValidationResult {
    bool ok = false;
 
    std::string message;
 
    explicit operator bool() const noexcept { return ok; }
};
 
struct TreeEditValidationStatistics {
    std::size_t completeValidationCommits = 0;
    std::size_t incrementalValidationCommits = 0;
};
 
enum class TreeEditValidationMode { Complete, Incremental };
 
// Forward declaration for the edit-session wrapper.
class TreeEditor;
 
class MorphologicalTree {
  private:
    friend class TreeEditor;
    friend class MorphologicalTreeFactory;
    friend class detail::CommittedTreeAccess;
 
    class LCAEulerRMQ; // Forward declaration for the LCA cache implementation.
 
    // ========================= Private attributes ========================= //
    NodeId rootNodeId_ = InvalidNode;
    MorphologicalTreeSemantics semantics_;
    std::optional<GridDomain2D> gridDomain2D_;
    int numNodes_ = 0;
    std::optional<NodeId> preservedExternalNodeIdOffset_;
    bool editSessionOpen_ = false;
    TreeEditValidationStatistics editValidationStatistics_;
 
    // Smallest-node map, indexed by pixel id [0, numPixels()).
    std::vector<NodeId> smallestNodeMap_;
 
    // Parent links, indexed by local node-slot id [0, numInternalNodeSlots()).
    std::vector<NodeId> nodeParent_;
 
    // Internal hierarchy linked structure, indexed by local node-slot id [0, numInternalNodeSlots()).
    std::vector<NodeId> firstChild_;
    std::vector<NodeId> nextSibling_;
    std::vector<NodeId> prevSibling_;
    std::vector<NodeId> lastChild_;
    std::vector<int> numChildrenByNode_;
 
    // Free slot management and alive-node iteration, indexed by local node-slot id [0, numInternalNodeSlots()).
    std::vector<uint8_t> alive_;
    std::vector<NodeId> freeNodeIds_;
 
    // Per-node proper-part pixel lists, indexed by local node-slot id.
    std::vector<PixelId> properHead_;
    std::vector<PixelId> properTail_;
    std::vector<int> properPartCardinalityByNode_;
    std::vector<PixelId> nextProperPart_;
    std::vector<PixelId> prevProperPart_;
 
    // Structural caches for traversal-based queries, indexed by local node-slot id [0, numInternalNodeSlots()).
    struct DfsIntervalCache {
        std::vector<int> entryIndex;
        std::vector<int> exitIndex;
        bool valid = false;
        void invalidate() noexcept { valid = false; }
    };
    mutable DfsIntervalCache dfsIntervalCache_;
    mutable std::unique_ptr<LCAEulerRMQ> lcaCache_;
 
    struct NodeSupportMetadataCache {
        std::vector<std::int32_t> cardinalityByNode;
        std::vector<PixelId> smallestPixelByNode;
        std::size_t nodeStructureVersion = std::numeric_limits<std::size_t>::max();
        std::size_t topologyVersion = std::numeric_limits<std::size_t>::max();
        std::size_t properPartVersion = std::numeric_limits<std::size_t>::max();
 
        void invalidate() noexcept {
            nodeStructureVersion = std::numeric_limits<std::size_t>::max();
            topologyVersion = std::numeric_limits<std::size_t>::max();
            properPartVersion = std::numeric_limits<std::size_t>::max();
        }
    };
    mutable NodeSupportMetadataCache nodeSupportMetadataCache_;
 
    // Version counters for iterator invalidation.
    std::size_t nodeStructureVersion_ = 0;
    std::size_t topologyVersion_ = 0;
    std::size_t properPartVersion_ = 0;
    std::size_t mutationVersion_ = 0;
 
    enum class ChildSplicePolicy { AppendToTargetTail, ReplaceSourceSlotWhenDirectChild };
 
    // ========================= Private methods ========================= //
    inline NodeId allocateSlot() {
        if (!freeNodeIds_.empty()) {
            const NodeId slotId = freeNodeIds_.back();
            freeNodeIds_.pop_back();
 
            nodeParent_[slotId] = InvalidNode;
            firstChild_[slotId] = InvalidNode;
            nextSibling_[slotId] = InvalidNode;
            prevSibling_[slotId] = InvalidNode;
            lastChild_[slotId] = InvalidNode;
            numChildrenByNode_[slotId] = 0;
            alive_[slotId] = 1;
            properHead_[slotId] = InvalidPixel;
            properTail_[slotId] = InvalidPixel;
            properPartCardinalityByNode_[slotId] = 0;
            return slotId;
        }
 
        const NodeId slotId = static_cast<NodeId>(nodeParent_.size());
        nodeParent_.push_back(InvalidNode);
        firstChild_.push_back(InvalidNode);
        nextSibling_.push_back(InvalidNode);
        prevSibling_.push_back(InvalidNode);
        lastChild_.push_back(InvalidNode);
        numChildrenByNode_.push_back(0);
        alive_.push_back(1);
        properHead_.push_back(InvalidPixel);
        properTail_.push_back(InvalidPixel);
        properPartCardinalityByNode_.push_back(0);
        return slotId;
    }
 
    inline void initializeProperPartStorage(size_t numPixels) {
        nextProperPart_.assign(numPixels, InvalidPixel);
        prevProperPart_.assign(numPixels, InvalidPixel);
    }
 
    inline void invalidateHigraNodeIdSpace() noexcept { preservedExternalNodeIdOffset_.reset(); }
 
    inline void preserveExternalNodeIdOffset(NodeId internalNodeOffset) {
        if (internalNodeOffset < 0) {
            throw std::invalid_argument("An external internal-node offset must be non-negative.");
        }
        preservedExternalNodeIdOffset_ = internalNodeOffset;
    }
 
    inline void beginEditSession() {
        if (editSessionOpen_) {
            throw std::logic_error("A MorphologicalTree edit session is already open.");
        }
        editSessionOpen_ = true;
    }
 
    inline void endEditSession() noexcept { editSessionOpen_ = false; }
 
    inline void recordEditCommit(TreeEditValidationMode mode) noexcept {
        if (mode == TreeEditValidationMode::Incremental) {
            ++editValidationStatistics_.incrementalValidationCommits;
        } else {
            ++editValidationStatistics_.completeValidationCommits;
        }
    }
 
    inline void initializeEmptyStorage(size_t numPixels) {
        rootNodeId_ = InvalidNode;
        numNodes_ = 0;
        smallestNodeMap_.assign(numPixels, InvalidNode);
 
        nodeParent_.clear();
        firstChild_.clear();
        nextSibling_.clear();
        prevSibling_.clear();
        lastChild_.clear();
        numChildrenByNode_.clear();
        alive_.clear();
        freeNodeIds_.clear();
        properHead_.clear();
        properTail_.clear();
        properPartCardinalityByNode_.clear();
        initializeProperPartStorage(numPixels);
        invalidateHigraNodeIdSpace();
 
        dfsIntervalCache_.entryIndex.clear();
        dfsIntervalCache_.exitIndex.clear();
        invalidateDfsIntervalCache();
        nodeSupportMetadataCache_.cardinalityByNode.clear();
        nodeSupportMetadataCache_.smallestPixelByNode.clear();
        nodeSupportMetadataCache_.invalidate();
        invalidateAllIterators();
    }
 
    inline void releaseSlotStorage(NodeId slotId) {
        // Grow the free-list before changing the slot. If allocation fails, the
        // live topology is therefore left untouched instead of terminating from
        // inside a falsely-noexcept operation or publishing a half-released slot.
        freeNodeIds_.push_back(slotId);
        nodeParent_[slotId] = InvalidNode;
        firstChild_[slotId] = InvalidNode;
        nextSibling_[slotId] = InvalidNode;
        prevSibling_[slotId] = InvalidNode;
        lastChild_[slotId] = InvalidNode;
        numChildrenByNode_[slotId] = 0;
        alive_[slotId] = 0;
        properHead_[slotId] = InvalidPixel;
        properTail_[slotId] = InvalidPixel;
        properPartCardinalityByNode_[slotId] = 0;
    }
 
    inline bool isFreeSlot(NodeId slotId) const noexcept { return slotId >= 0 && slotId < static_cast<NodeId>(alive_.size()) && alive_[slotId] == 0; }
 
    inline void requireAliveNode(NodeId nodeId, const char* context) const {
        MMCFILTERS_CONTRACT_REQUIRE(isAlive(nodeId), throw std::invalid_argument(std::string(context) + " requires a live internal NodeId."));
    }
 
    inline void requireAliveNonRootNode(NodeId nodeId, const char* context) const {
        requireAliveNode(nodeId, context);
        MMCFILTERS_CONTRACT_REQUIRE(!isRoot(nodeId), throw std::invalid_argument(std::string(context) + " cannot target the root node."));
    }
 
    inline void rebuildProperPartLinksFromSmallestNodeMap() {
        properHead_.assign(nodeParent_.size(), InvalidPixel);
        properTail_.assign(nodeParent_.size(), InvalidPixel);
        properPartCardinalityByNode_.assign(nodeParent_.size(), 0);
        initializeProperPartStorage(smallestNodeMap_.size());
 
        for (PixelId pixel = 0; pixel < static_cast<PixelId>(smallestNodeMap_.size()); ++pixel) {
            const NodeId smallestNodeSlotId = smallestNodeMap_[static_cast<std::size_t>(pixel)];
            if (smallestNodeSlotId == InvalidNode || smallestNodeSlotId >= static_cast<NodeId>(nodeParent_.size()) || isFreeSlot(smallestNodeSlotId)) {
                continue;
            }
 
            const PixelId tailProperPartId = properTail_[smallestNodeSlotId];
            if (tailProperPartId == InvalidPixel) {
                properHead_[smallestNodeSlotId] = pixel;
                properTail_[smallestNodeSlotId] = pixel;
            } else {
                nextProperPart_[tailProperPartId] = pixel;
                prevProperPart_[pixel] = tailProperPartId;
                properTail_[smallestNodeSlotId] = pixel;
            }
            properPartCardinalityByNode_[smallestNodeSlotId]++;
        }
    }
 
    inline void unlinkPixelFromProperPart(NodeId smallestNodeSlotId, PixelId pixel) noexcept {
        const PixelId prev = prevProperPart_[static_cast<size_t>(pixel)];
        const PixelId next = nextProperPart_[static_cast<size_t>(pixel)];
 
        if (prev == InvalidPixel) {
            properHead_[static_cast<size_t>(smallestNodeSlotId)] = next;
        } else {
            nextProperPart_[static_cast<size_t>(prev)] = next;
        }
 
        if (next == InvalidPixel) {
            properTail_[static_cast<size_t>(smallestNodeSlotId)] = prev;
        } else {
            prevProperPart_[static_cast<size_t>(next)] = prev;
        }
 
        nextProperPart_[static_cast<size_t>(pixel)] = InvalidPixel;
        prevProperPart_[static_cast<size_t>(pixel)] = InvalidPixel;
        --properPartCardinalityByNode_[static_cast<size_t>(smallestNodeSlotId)];
    }
 
    inline void appendDetachedProperPart(NodeId targetSlotId, PixelId pixel) noexcept {
        smallestNodeMap_[static_cast<size_t>(pixel)] = targetSlotId;
        const PixelId tail = properTail_[static_cast<size_t>(targetSlotId)];
        if (tail == InvalidPixel) {
            properHead_[static_cast<size_t>(targetSlotId)] = pixel;
            properTail_[static_cast<size_t>(targetSlotId)] = pixel;
        } else {
            nextProperPart_[static_cast<size_t>(tail)] = pixel;
            prevProperPart_[static_cast<size_t>(pixel)] = tail;
            properTail_[static_cast<size_t>(targetSlotId)] = pixel;
        }
        ++properPartCardinalityByNode_[static_cast<size_t>(targetSlotId)];
    }
 
    inline void spliceProperPartsSlots(NodeId targetSlotId, NodeId sourceSlotId) noexcept {
        const PixelId sourceHead = properHead_[static_cast<size_t>(sourceSlotId)];
        if (sourceHead == InvalidPixel) {
            return;
        }
 
        for (PixelId pixel = sourceHead; pixel != InvalidPixel; pixel = nextProperPart_[static_cast<size_t>(pixel)]) {
            smallestNodeMap_[static_cast<size_t>(pixel)] = targetSlotId;
        }
 
        const PixelId sourceTail = properTail_[static_cast<size_t>(sourceSlotId)];
        const int sourceCount = properPartCardinalityByNode_[static_cast<size_t>(sourceSlotId)];
        const PixelId targetTail = properTail_[static_cast<size_t>(targetSlotId)];
 
        if (targetTail == InvalidPixel) {
            properHead_[static_cast<size_t>(targetSlotId)] = sourceHead;
            properTail_[static_cast<size_t>(targetSlotId)] = sourceTail;
        } else {
            nextProperPart_[static_cast<size_t>(targetTail)] = sourceHead;
            prevProperPart_[static_cast<size_t>(sourceHead)] = targetTail;
            properTail_[static_cast<size_t>(targetSlotId)] = sourceTail;
        }
 
        properPartCardinalityByNode_[static_cast<size_t>(targetSlotId)] += sourceCount;
        properHead_[static_cast<size_t>(sourceSlotId)] = InvalidPixel;
        properTail_[static_cast<size_t>(sourceSlotId)] = InvalidPixel;
        properPartCardinalityByNode_[static_cast<size_t>(sourceSlotId)] = 0;
    }
 
    inline void releaseSlotNode(NodeId slotNodeId) {
        releaseSlotStorage(slotNodeId);
        numNodes_--;
        invalidateDfsIntervalCache();
        bumpNodeStructureVersion();
    }
 
    inline void invalidateAllIterators() noexcept {
        nodeStructureVersion_ = 0;
        topologyVersion_ = 0;
        properPartVersion_ = 0;
        ++mutationVersion_;
        lcaCache_.reset();
        nodeSupportMetadataCache_.invalidate();
    }
 
    inline void invalidateLcaCache() const noexcept { lcaCache_.reset(); }
 
    inline void bumpNodeStructureVersion() noexcept {
        ++nodeStructureVersion_;
        ++mutationVersion_;
        invalidateLcaCache();
        invalidateHigraNodeIdSpace();
    }
 
    inline void bumpTopologyVersion() noexcept {
        ++topologyVersion_;
        ++mutationVersion_;
        invalidateLcaCache();
        invalidateDfsIntervalCache();
        invalidateHigraNodeIdSpace();
    }
 
    inline void bumpProperPartVersion() noexcept {
        ++properPartVersion_;
        ++mutationVersion_;
        invalidateHigraNodeIdSpace();
    }
 
    inline void checkNodeIteratorVersion([[maybe_unused]] std::size_t expectedVersion) const {
        assert(expectedVersion == nodeStructureVersion_ && "Alive-node iterator invalidated by node-structure mutation.");
    }
 
    inline void checkTopologyIteratorVersion([[maybe_unused]] std::size_t expectedVersion) const {
        assert(expectedVersion == topologyVersion_ && "Topology iterator invalidated by tree-structure mutation.");
    }
 
    inline void checkProperPartIteratorVersion([[maybe_unused]] std::size_t expectedVersion) const {
        assert(expectedVersion == properPartVersion_ && "Proper-parts iterator invalidated by proper-part mutation.");
    }
 
    inline void invalidateDfsIntervalCache() const noexcept { dfsIntervalCache_.invalidate(); }
 
    inline void recomputeDfsIntervalCache() const {
        dfsIntervalCache_.entryIndex.assign(nodeParent_.size(), -1);
        dfsIntervalCache_.exitIndex.assign(nodeParent_.size(), -1);
 
        if (rootNodeId_ == InvalidNode) {
            dfsIntervalCache_.valid = true;
            return;
        }
 
        // Explicit frames preserve child order without using the call stack on deep trees.
        int nextDfsEventIndex = 0;
        std::vector<std::pair<NodeId, NodeId>> stack;
        stack.emplace_back(rootNodeId_, firstChild_[rootNodeId_]);
        dfsIntervalCache_.entryIndex[rootNodeId_] = nextDfsEventIndex++;
        while (!stack.empty()) {
            auto& [node, nextChild] = stack.back();
            if (nextChild == InvalidNode) {
                dfsIntervalCache_.exitIndex[node] = nextDfsEventIndex++;
                stack.pop_back();
            } else {
                const NodeId child = nextChild;
                nextChild = nextSibling_[child];
                dfsIntervalCache_.entryIndex[child] = nextDfsEventIndex++;
                stack.emplace_back(child, firstChild_[child]);
            }
        }
 
        dfsIntervalCache_.valid = true;
    }
 
    inline void ensureDfsIntervalCache() const {
        if (!dfsIntervalCache_.valid) {
            recomputeDfsIntervalCache();
        }
    }
 
    inline const LCAEulerRMQ& ensureLcaCache() const {
        if (!lcaCache_) {
            lcaCache_ = std::make_unique<LCAEulerRMQ>(this);
        }
        return *lcaCache_;
    }
 
    [[nodiscard]] inline std::optional<NodeId> lowestCommonAncestorFromDfsIntervalsEstablished(NodeId u, NodeId v) const {
        if (u == v) {
            return u;
        }
        ensureDfsIntervalCache();
        if (dfsIntervalCache_.entryIndex[static_cast<std::size_t>(u)] < 0 ||
            dfsIntervalCache_.entryIndex[static_cast<std::size_t>(v)] < 0) {
            return InvalidNode;
        }
        const auto isAncestorCached = [&](NodeId ancestor, NodeId node) {
            return dfsIntervalCache_.entryIndex[static_cast<std::size_t>(ancestor)] <=
                       dfsIntervalCache_.entryIndex[static_cast<std::size_t>(node)] &&
                   dfsIntervalCache_.exitIndex[static_cast<std::size_t>(ancestor)] >=
                       dfsIntervalCache_.exitIndex[static_cast<std::size_t>(node)];
        };
        if (isAncestorCached(u, v)) {
            return u;
        }
        if (isAncestorCached(v, u)) {
            return v;
        }
        return std::nullopt;
    }
 
    [[nodiscard]] inline NodeId lowestCommonAncestorEstablished(NodeId u, NodeId v) const {
        const std::optional<NodeId> intervalResult = lowestCommonAncestorFromDfsIntervalsEstablished(u, v);
        return intervalResult.has_value() ? *intervalResult : ensureLcaCache().findLowestCommonAncestor(u, v);
    }
 
    [[nodiscard]] inline long double estimatedLcaRmqStorageBytes() const noexcept {
        std::size_t rmqLevels = 1;
        const std::size_t eulerLength = numNodes_ > 0 ? static_cast<std::size_t>(numNodes_) * 2 - 1 : 0;
        for (std::size_t blockSize = 2; blockSize <= eulerLength && blockSize <= std::numeric_limits<std::size_t>::max() / 2;
             blockSize *= 2) {
            ++rmqLevels;
        }
        return static_cast<long double>(nodeParent_.size()) * static_cast<long double>(36 + 8 * rmqLevels);
    }
 
    template <typename QueryAccessor, typename ResultConsumer>
    inline void forEachLowestCommonAncestorEstablished(std::size_t numQueries, QueryAccessor&& queryForIndex, ResultConsumer&& consumeLca,
                                                       bool forceTarjan = false) const {
        if (numQueries == 0) {
            return;
        }
 
        if (lcaCache_ && !forceTarjan) {
            for (std::size_t queryIndex = 0; queryIndex < numQueries; ++queryIndex) {
                const auto [first, second] = queryForIndex(queryIndex);
                consumeLca(queryIndex, lowestCommonAncestorEstablished(first, second));
            }
            return;
        }
 
        const std::size_t numSlots = nodeParent_.size();
        std::size_t numUnresolvedQueries = numQueries;
        if (!forceTarjan) {
            numUnresolvedQueries = 0;
            for (std::size_t queryIndex = 0; queryIndex < numQueries; ++queryIndex) {
                const auto [first, second] = queryForIndex(queryIndex);
                const std::optional<NodeId> intervalResult = lowestCommonAncestorFromDfsIntervalsEstablished(first, second);
                if (intervalResult.has_value()) {
                    consumeLca(queryIndex, *intervalResult);
                    continue;
                }
                ++numUnresolvedQueries;
            }
        }
        if (numUnresolvedQueries == 0 || rootNodeId_ == InvalidNode) {
            return;
        }
 
        const long double estimatedRmqBytes = estimatedLcaRmqStorageBytes();
        const long double estimatedTarjanBytes =
            static_cast<long double>(numSlots) * 30.0L + static_cast<long double>(numUnresolvedQueries) * 16.0L;
        const bool batchFitsCompactIndices = numQueries <= static_cast<std::size_t>(std::numeric_limits<std::uint32_t>::max() / 2);
        if (!forceTarjan && (!batchFitsCompactIndices || estimatedRmqBytes <= estimatedTarjanBytes)) {
            for (std::size_t queryIndex = 0; queryIndex < numQueries; ++queryIndex) {
                const auto [first, second] = queryForIndex(queryIndex);
                if (lowestCommonAncestorFromDfsIntervalsEstablished(first, second).has_value()) {
                    continue;
                }
                consumeLca(queryIndex, ensureLcaCache().findLowestCommonAncestor(first, second));
            }
            return;
        }
 
        using QueryOffset = std::uint32_t;
        std::vector<QueryOffset> queryOffsets(numSlots + 1, QueryOffset{0});
        for (std::size_t queryIndex = 0; queryIndex < numQueries; ++queryIndex) {
            const auto [first, second] = queryForIndex(queryIndex);
            if (!forceTarjan && lowestCommonAncestorFromDfsIntervalsEstablished(first, second).has_value()) {
                continue;
            }
            ++queryOffsets[static_cast<std::size_t>(first) + 1];
            ++queryOffsets[static_cast<std::size_t>(second) + 1];
        }
        for (std::size_t index = 1; index < queryOffsets.size(); ++index) {
            queryOffsets[index] += queryOffsets[index - 1];
        }
 
        struct QueryReference {
            NodeId otherEndpoint = InvalidNode;
            std::uint32_t queryIndex = 0;
        };
        std::vector<QueryReference> queryReferences(queryOffsets.back());
        for (std::size_t queryIndex = 0; queryIndex < numQueries; ++queryIndex) {
            const auto [first, second] = queryForIndex(queryIndex);
            if (!forceTarjan && lowestCommonAncestorFromDfsIntervalsEstablished(first, second).has_value()) {
                continue;
            }
            const QueryOffset firstPosition = queryOffsets[static_cast<std::size_t>(first)]++;
            queryReferences[firstPosition] = {second, static_cast<std::uint32_t>(queryIndex)};
            const QueryOffset secondPosition = queryOffsets[static_cast<std::size_t>(second)]++;
            queryReferences[secondPosition] = {first, static_cast<std::uint32_t>(queryIndex)};
        }
        for (std::size_t index = numSlots; index > 0; --index) {
            queryOffsets[index] = queryOffsets[index - 1];
        }
        queryOffsets[0] = QueryOffset{0};
 
        const auto findSet = [](std::vector<NodeId>& setParents, NodeId node) {
            NodeId representative = node;
            while (setParents[static_cast<std::size_t>(representative)] != representative) {
                representative = setParents[static_cast<std::size_t>(representative)];
            }
            while (setParents[static_cast<std::size_t>(node)] != node) {
                const NodeId next = setParents[static_cast<std::size_t>(node)];
                setParents[static_cast<std::size_t>(node)] = representative;
                node = next;
            }
            return representative;
        };
        const auto uniteSets = [&](std::vector<NodeId>& setParents, std::vector<std::uint8_t>& setRanks, NodeId first, NodeId second) {
            NodeId firstRoot = findSet(setParents, first);
            NodeId secondRoot = findSet(setParents, second);
            if (firstRoot == secondRoot) {
                return firstRoot;
            }
            if (setRanks[static_cast<std::size_t>(firstRoot)] < setRanks[static_cast<std::size_t>(secondRoot)]) {
                std::swap(firstRoot, secondRoot);
            }
            setParents[static_cast<std::size_t>(secondRoot)] = firstRoot;
            if (setRanks[static_cast<std::size_t>(firstRoot)] == setRanks[static_cast<std::size_t>(secondRoot)]) {
                ++setRanks[static_cast<std::size_t>(firstRoot)];
            }
            return firstRoot;
        };
 
        struct TraversalFrame {
            NodeId node = InvalidNode;
            NodeId nextChild = InvalidNode;
        };
 
        std::vector<NodeId> setParents(numSlots, InvalidNode);
        std::vector<std::uint8_t> setRanks(numSlots, std::uint8_t{0});
        std::vector<NodeId> ancestors(numSlots, InvalidNode);
        std::vector<std::uint8_t> finishedNodes(numSlots, std::uint8_t{0});
        std::vector<TraversalFrame> traversalStack;
        traversalStack.reserve(static_cast<std::size_t>(numNodes_));
        setParents[static_cast<std::size_t>(rootNodeId_)] = rootNodeId_;
        ancestors[static_cast<std::size_t>(rootNodeId_)] = rootNodeId_;
        traversalStack.push_back({rootNodeId_, firstChild_[static_cast<std::size_t>(rootNodeId_)]});
 
        while (!traversalStack.empty()) {
            TraversalFrame& frame = traversalStack.back();
            if (frame.nextChild != InvalidNode) {
                const NodeId child = frame.nextChild;
                frame.nextChild = nextSibling_[static_cast<std::size_t>(child)];
                setParents[static_cast<std::size_t>(child)] = child;
                ancestors[static_cast<std::size_t>(child)] = child;
                traversalStack.push_back({child, firstChild_[static_cast<std::size_t>(child)]});
                continue;
            }
 
            const NodeId node = frame.node;
            const std::size_t nodeIndex = static_cast<std::size_t>(node);
            finishedNodes[nodeIndex] = std::uint8_t{1};
            for (QueryOffset referenceIndex = queryOffsets[nodeIndex]; referenceIndex < queryOffsets[nodeIndex + 1]; ++referenceIndex) {
                const QueryReference query = queryReferences[referenceIndex];
                if (finishedNodes[static_cast<std::size_t>(query.otherEndpoint)] != 0) {
                    const NodeId representative = findSet(setParents, query.otherEndpoint);
                    consumeLca(query.queryIndex, ancestors[static_cast<std::size_t>(representative)]);
                }
            }
 
            traversalStack.pop_back();
            if (!traversalStack.empty()) {
                const NodeId parent = traversalStack.back().node;
                const NodeId representative = uniteSets(setParents, setRanks, parent, node);
                ancestors[static_cast<std::size_t>(representative)] = parent;
            }
        }
    }
 
    [[nodiscard]] inline std::vector<NodeId>
    lowestCommonAncestorsEstablished(std::span<const std::pair<NodeId, NodeId>> queries) const {
        std::vector<NodeId> lcas(queries.size(), InvalidNode);
        forEachLowestCommonAncestorEstablished(
            queries.size(), [&](std::size_t queryIndex) { return queries[queryIndex]; },
            [&](std::size_t queryIndex, NodeId lca) { lcas[queryIndex] = lca; });
        return lcas;
    }
 
    template <typename QueryAccessor, typename ResultConsumer>
    inline void forEachGeneratedLowestCommonAncestorEstablished(std::size_t numQueries, QueryAccessor&& queryForIndex,
                                                                 ResultConsumer&& consumeLca) const {
        if (numQueries == 0) {
            return;
        }
        if (lcaCache_) {
            for (std::size_t queryIndex = 0; queryIndex < numQueries; ++queryIndex) {
                const auto [first, second] = queryForIndex(queryIndex);
                consumeLca(queryIndex, lowestCommonAncestorEstablished(first, second));
            }
            return;
        }
 
        const bool fullBatchFitsCompactIndices = numQueries <= static_cast<std::size_t>(std::numeric_limits<std::uint32_t>::max() / 2);
        const long double estimatedFullBatchTarjanBytes =
            static_cast<long double>(nodeParent_.size()) * 30.0L + static_cast<long double>(numQueries) * 24.0L;
        if (fullBatchFitsCompactIndices && estimatedLcaRmqStorageBytes() > estimatedFullBatchTarjanBytes) {
            constexpr std::size_t MaxStrategyProbes = 1024;
            const std::size_t numProbes = std::min(numQueries, MaxStrategyProbes);
            bool foundIncomparablePair = false;
            for (std::size_t probe = 0; probe < numProbes; ++probe) {
                const std::size_t queryIndex = probe * numQueries / numProbes;
                const auto [first, second] = queryForIndex(queryIndex);
                if (!lowestCommonAncestorFromDfsIntervalsEstablished(first, second).has_value()) {
                    foundIncomparablePair = true;
                    break;
                }
            }
            if (foundIncomparablePair) {
                std::vector<std::pair<NodeId, NodeId>> allPairs(numQueries);
                for (std::size_t queryIndex = 0; queryIndex < numQueries; ++queryIndex) {
                    allPairs[queryIndex] = queryForIndex(queryIndex);
                }
                forEachLowestCommonAncestorEstablished(
                    numQueries, [&](std::size_t queryIndex) { return allPairs[queryIndex]; }, std::forward<ResultConsumer>(consumeLca), true);
                return;
            }
        }
 
        std::size_t numUnresolvedQueries = 0;
        for (std::size_t queryIndex = 0; queryIndex < numQueries; ++queryIndex) {
            const auto [first, second] = queryForIndex(queryIndex);
            const std::optional<NodeId> intervalResult = lowestCommonAncestorFromDfsIntervalsEstablished(first, second);
            if (intervalResult.has_value()) {
                consumeLca(queryIndex, *intervalResult);
            } else {
                ++numUnresolvedQueries;
            }
        }
        if (numUnresolvedQueries == 0) {
            return;
        }
 
        const bool filteredBatchFitsCompactIndices = numQueries <= static_cast<std::size_t>(std::numeric_limits<std::uint32_t>::max()) &&
                                                     numUnresolvedQueries <=
                                                         static_cast<std::size_t>(std::numeric_limits<std::uint32_t>::max() / 2);
        const long double estimatedFilteredTarjanBytes =
            static_cast<long double>(nodeParent_.size()) * 30.0L + static_cast<long double>(numUnresolvedQueries) * 28.0L;
        if (!filteredBatchFitsCompactIndices || estimatedLcaRmqStorageBytes() <= estimatedFilteredTarjanBytes) {
            for (std::size_t queryIndex = 0; queryIndex < numQueries; ++queryIndex) {
                const auto [first, second] = queryForIndex(queryIndex);
                if (!lowestCommonAncestorFromDfsIntervalsEstablished(first, second).has_value()) {
                    consumeLca(queryIndex, ensureLcaCache().findLowestCommonAncestor(first, second));
                }
            }
            return;
        }
 
        std::vector<std::pair<NodeId, NodeId>> unresolvedPairs;
        std::vector<std::uint32_t> unresolvedQueryIndices;
        unresolvedPairs.reserve(numUnresolvedQueries);
        unresolvedQueryIndices.reserve(numUnresolvedQueries);
        for (std::size_t queryIndex = 0; queryIndex < numQueries; ++queryIndex) {
            const auto query = queryForIndex(queryIndex);
            if (!lowestCommonAncestorFromDfsIntervalsEstablished(query.first, query.second).has_value()) {
                unresolvedPairs.push_back(query);
                unresolvedQueryIndices.push_back(static_cast<std::uint32_t>(queryIndex));
            }
        }
        forEachLowestCommonAncestorEstablished(
            unresolvedPairs.size(), [&](std::size_t queryIndex) { return unresolvedPairs[queryIndex]; },
            [&](std::size_t queryIndex, NodeId lca) { consumeLca(unresolvedQueryIndices[queryIndex], lca); }, true);
    }
 
    // ========================= Internal topology helpers ========================= //
    // These helpers mutate the dense node-indexed topology storage directly.
    void linkChildSlot(NodeId parentSlotId, NodeId childId) {
        if (parentSlotId < 0 || childId < 0)
            return;
 
        // If the child already has a parent, detach it before reading sibling links.
        if (nodeParent_[childId] != InvalidNode) {
            unlinkChildSlot(nodeParent_[childId], childId);
        }
 
        nodeParent_[childId] = parentSlotId;
        prevSibling_[childId] = lastChild_[parentSlotId];
        nextSibling_[childId] = InvalidNode;
 
        if (firstChild_[parentSlotId] == InvalidNode) {
            firstChild_[parentSlotId] = lastChild_[parentSlotId] = childId;
        } else {
            nextSibling_[lastChild_[parentSlotId]] = childId;
            lastChild_[parentSlotId] = childId;
        }
        ++numChildrenByNode_[parentSlotId];
        bumpTopologyVersion();
    }
 
    inline void unlinkChildSlot(NodeId parentSlotId, NodeId childId) {
        if (parentSlotId < 0 || childId < 0)
            return;
        if (nodeParent_[childId] != parentSlotId)
            return;
 
        const NodeId prev = prevSibling_[childId];
        const NodeId next = nextSibling_[childId];
 
        if (prev == InvalidNode)
            firstChild_[parentSlotId] = next;
        else
            nextSibling_[prev] = next;
 
        if (next == InvalidNode)
            lastChild_[parentSlotId] = prev;
        else
            prevSibling_[next] = prev;
 
        if (numChildrenByNode_[parentSlotId] > 0)
            --numChildrenByNode_[parentSlotId];
 
        nodeParent_[childId] = InvalidNode;
        prevSibling_[childId] = InvalidNode;
        nextSibling_[childId] = InvalidNode;
        bumpTopologyVersion();
    }
 
    inline void spliceChildrenSlots(NodeId toId, NodeId fromId, ChildSplicePolicy policy = ChildSplicePolicy::AppendToTargetTail) {
        if (toId < 0 || fromId < 0 || toId == fromId)
            return;
 
        const NodeId firstFrom = firstChild_[fromId];
        const NodeId lastFrom = lastChild_[fromId];
        const int movedCount = numChildrenByNode_[fromId];
        const bool replaceSourceSlot = policy == ChildSplicePolicy::ReplaceSourceSlotWhenDirectChild && nodeParent_[fromId] == toId;
 
        for (NodeId childId = firstFrom; childId != InvalidNode; childId = nextSibling_[childId]) {
            nodeParent_[childId] = toId;
        }
 
        if (replaceSourceSlot) {
            const NodeId prev = prevSibling_[fromId];
            const NodeId next = nextSibling_[fromId];
 
            if (firstFrom == InvalidNode) {
                if (prev == InvalidNode) {
                    firstChild_[toId] = next;
                } else {
                    nextSibling_[prev] = next;
                }
 
                if (next == InvalidNode) {
                    lastChild_[toId] = prev;
                } else {
                    prevSibling_[next] = prev;
                }
            } else {
                if (prev == InvalidNode) {
                    firstChild_[toId] = firstFrom;
                    prevSibling_[firstFrom] = InvalidNode;
                } else {
                    nextSibling_[prev] = firstFrom;
                    prevSibling_[firstFrom] = prev;
                }
 
                if (next == InvalidNode) {
                    lastChild_[toId] = lastFrom;
                    nextSibling_[lastFrom] = InvalidNode;
                } else {
                    prevSibling_[next] = lastFrom;
                    nextSibling_[lastFrom] = next;
                }
            }
 
            nodeParent_[fromId] = InvalidNode;
            prevSibling_[fromId] = InvalidNode;
            nextSibling_[fromId] = InvalidNode;
            numChildrenByNode_[toId] += movedCount - 1;
        } else {
            if (firstFrom == InvalidNode)
                return;
 
            if (firstChild_[toId] == InvalidNode) {
                firstChild_[toId] = firstFrom;
                lastChild_[toId] = lastFrom;
            } else {
                nextSibling_[lastChild_[toId]] = firstFrom;
                prevSibling_[firstFrom] = lastChild_[toId];
                lastChild_[toId] = lastFrom;
            }
 
            numChildrenByNode_[toId] += movedCount;
        }
 
        firstChild_[fromId] = InvalidNode;
        lastChild_[fromId] = InvalidNode;
        numChildrenByNode_[fromId] = 0;
        bumpTopologyVersion();
    }
 
    inline void releaseNode(NodeId nodeId) {
        if (!isNode(nodeId) || !isAlive(nodeId) || isRoot(nodeId)) {
            return;
        }
        const NodeId nodeSlot = nodeId;
        if (nodeParent_[nodeSlot] != nodeSlot) {
            return;
        }
        if (numChildrenByNode_[nodeSlot] != 0 || properPartCardinalityByNode_[nodeSlot] != 0) {
            return;
        }
        releaseSlotNode(nodeSlot);
    }
 
    inline NodeId allocateNode() {
        if (freeNodeIds_.empty()) {
            return InvalidNode;
        }
        const NodeId nodeSlot = allocateSlot();
        nodeParent_[nodeSlot] = nodeSlot;
        numNodes_++;
        invalidateDfsIntervalCache();
        bumpNodeStructureVersion();
        return nodeSlot;
    }
 
    inline NodeId createDetachedNode() {
        const NodeId nodeSlot = allocateSlot();
        nodeParent_[static_cast<size_t>(nodeSlot)] = nodeSlot;
        numNodes_++;
        invalidateDfsIntervalCache();
        bumpNodeStructureVersion();
        return nodeSlot;
    }
 
    inline void removeChild(NodeId parentNodeId, NodeId childId, bool releaseNodeFlag) {
        if (!isAlive(parentNodeId) || !isAlive(childId) || !hasChild(parentNodeId, childId)) {
            return;
        }
        const NodeId parentSlotId = parentNodeId;
        const NodeId childSlotId = childId;
        if (releaseNodeFlag && numChildrenByNode_[static_cast<std::size_t>(childSlotId)] == 0 &&
            properPartCardinalityByNode_[static_cast<std::size_t>(childSlotId)] == 0) {
            freeNodeIds_.reserve(freeNodeIds_.size() + 1);
        }
        unlinkChildSlot(parentSlotId, childSlotId);
        nodeParent_[childSlotId] = childSlotId;
        invalidateDfsIntervalCache();
        if (releaseNodeFlag) {
            releaseNode(childId);
        }
    }
 
    inline void attachNode(NodeId parentNodeId, NodeId nodeId) {
        if (!isAlive(parentNodeId) || !isAlive(nodeId) || isRoot(nodeId) || parentNodeId == nodeId) {
            return;
        }
        const NodeId parentSlotId = parentNodeId;
        const NodeId nodeSlotId = nodeId;
        const NodeId oldParentSlotId = nodeParent_[nodeSlotId];
        if (oldParentSlotId != InvalidNode && oldParentSlotId != nodeSlotId) {
            unlinkChildSlot(oldParentSlotId, nodeSlotId);
        }
        nodeParent_[nodeSlotId] = InvalidNode;
        linkChildSlot(parentSlotId, nodeSlotId);
    }
 
    inline void detachNode(NodeId nodeId) {
        if (!isAlive(nodeId) || isRoot(nodeId)) {
            return;
        }
        const NodeId nodeSlotId = nodeId;
        const NodeId parentSlotId = nodeParent_[nodeSlotId];
        if (parentSlotId == InvalidNode || parentSlotId == nodeSlotId) {
            return;
        }
        unlinkChildSlot(parentSlotId, nodeSlotId);
        nodeParent_[nodeSlotId] = nodeSlotId;
    }
 
    inline void moveNode(NodeId nodeId, NodeId newParentId) {
        if (!isAlive(nodeId) || !isAlive(newParentId) || isRoot(nodeId) || nodeId == newParentId) {
            return;
        }
        const NodeId nodeSlotId = nodeId;
        const NodeId newParentSlotId = newParentId;
        const NodeId oldParentSlotId = nodeParent_[nodeSlotId];
        if (oldParentSlotId == newParentSlotId) {
            return;
        }
        if (oldParentSlotId != InvalidNode && oldParentSlotId != nodeSlotId) {
            unlinkChildSlot(oldParentSlotId, nodeSlotId);
        }
        nodeParent_[nodeSlotId] = InvalidNode;
        linkChildSlot(newParentSlotId, nodeSlotId);
    }
 
    inline void moveChildren(NodeId parentNodeId, NodeId sourceId) {
        if (!isAlive(parentNodeId) || !isAlive(sourceId) || parentNodeId == sourceId) {
            return;
        }
        spliceChildrenSlots(parentNodeId, sourceId);
    }
 
    inline void movePixelToProperPart(NodeId targetNodeId, NodeId sourceNodeId, PixelId pixel) {
        if (!isAlive(targetNodeId) || !isAlive(sourceNodeId) || !isPixel(pixel) || targetNodeId == sourceNodeId) {
            return;
        }
        const NodeId sourceSlotId = sourceNodeId;
        if (smallestNodeMap_[pixel] != sourceSlotId) {
            return;
        }
        unlinkPixelFromProperPart(sourceSlotId, pixel);
        appendDetachedProperPart(targetNodeId, pixel);
        bumpProperPartVersion();
    }
 
    inline void mergeProperParts(NodeId targetNodeId, NodeId sourceNodeId) {
        if (!isAlive(targetNodeId) || !isAlive(sourceNodeId) || targetNodeId == sourceNodeId) {
            return;
        }
        if (properHead_[static_cast<size_t>(sourceNodeId)] == InvalidPixel) {
            return;
        }
        spliceProperPartsSlots(targetNodeId, sourceNodeId);
        bumpProperPartVersion();
    }
 
    inline void setRoot(NodeId nodeId) {
        if (!isAlive(nodeId)) {
            return;
        }
        const NodeId nodeSlot = nodeId;
        if (rootNodeId_ == nodeSlot) {
            return;
        }
        if (rootNodeId_ != InvalidNode && !isFreeSlot(rootNodeId_)) {
            nodeParent_[rootNodeId_] = rootNodeId_;
            prevSibling_[rootNodeId_] = InvalidNode;
            nextSibling_[rootNodeId_] = InvalidNode;
        }
        const NodeId oldParentSlot = nodeParent_[nodeSlot];
        if (oldParentSlot != InvalidNode && oldParentSlot != nodeSlot) {
            unlinkChildSlot(oldParentSlot, nodeSlot);
        }
        rootNodeId_ = nodeSlot;
        nodeParent_[nodeSlot] = nodeSlot;
        prevSibling_[nodeSlot] = InvalidNode;
        nextSibling_[nodeSlot] = InvalidNode;
        bumpTopologyVersion();
    }
 
    void initializeNativeTopologyStorage(std::vector<NodeId> nodeParent, std::vector<NodeId> smallestNodeMap, NodeId root,
                                         std::optional<GridDomain2D> gridDomain2D, MorphologicalTreeSemantics semantics) {
        if (gridDomain2D && smallestNodeMap.size() != gridDomain2D->size("Native topology 2D pixel domain")) {
            throw std::invalid_argument("Native topology pixel domain must match the attached 2D grid.");
        }
        if (nodeParent.size() > static_cast<std::size_t>(std::numeric_limits<NodeId>::max())) {
            throw std::invalid_argument("Native topology internal-node domain exceeds NodeId range.");
        }
        if (smallestNodeMap.size() > static_cast<std::size_t>(std::numeric_limits<PixelId>::max())) {
            throw std::invalid_argument("Native topology pixel domain exceeds PixelId range.");
        }
        switch (semantics.kind) {
        case MorphologicalTreeKind::Generic:
        case MorphologicalTreeKind::MaxTree:
        case MorphologicalTreeKind::MinTree:
        case MorphologicalTreeKind::TreeOfShapes:
        case MorphologicalTreeKind::UnrestrictedResidualTree:
        case MorphologicalTreeKind::SaturatedResidualTree:
            break;
        default:
            throw std::invalid_argument("Native topology kind is not supported.");
        }
        const int numNodeSlots = static_cast<int>(nodeParent.size());
        const std::size_t numPixels = smallestNodeMap.size();
        if (numNodeSlots <= 0) {
            throw std::invalid_argument("Native topology import requires at least one internal node.");
        }
        if (numPixels == 0) {
            throw std::invalid_argument("Native topology import requires at least one pixel.");
        }
        if (root < 0 || root >= numNodeSlots) {
            throw std::invalid_argument("Native topology import requires a valid root node id.");
        }
        validateMorphologicalTreeSemantics(semantics);
        semantics_ = std::move(semantics);
        gridDomain2D_ = gridDomain2D;
        const auto requireMatchingGrid = [this](const RegularGridAdjacency2D& adjacency, const char* context) {
            if (!gridDomain2D_) {
                throw std::invalid_argument(std::string(context) + " requires an attached 2D grid domain.");
            }
            if (adjacency.getNumRows() != gridDomain2D_->rows || adjacency.getNumColumns() != gridDomain2D_->columns) {
                throw std::invalid_argument(std::string(context) + " must match the attached 2D grid.");
            }
        };
        if (const auto* context = std::get_if<SharedAdjacencyContext>(&semantics_.constructionContext)) {
            requireMatchingGrid(context->adjacency, "SharedAdjacencyContext adjacency");
        } else if (const auto* context = std::get_if<SaturatedResidualContext>(&semantics_.constructionContext)) {
            requireMatchingGrid(context->adjacency, "SaturatedResidualContext adjacency");
            if (context->infinityPixel < 0 || static_cast<std::size_t>(context->infinityPixel) >= gridDomain2D_->size("Saturated residual domain")) {
                throw std::invalid_argument("SaturatedResidualContext infinity pixel must belong to the attached 2D grid.");
            }
        } else if (const auto* convention = std::get_if<TopographicConvention>(&semantics_.constructionContext)) {
            if (convention->infinityPixel < 0) {
                throw std::invalid_argument("TopographicConvention infinity pixel must be non-negative.");
            }
            if (!gridDomain2D_) {
                throw std::invalid_argument("TopographicConvention requires an attached 2D grid domain.");
            }
            const std::int64_t extension = convention->domainExtension == TopographicDomainExtension::ExteriorRing ? 1 : -1;
            const std::int64_t activeRows = 2 * static_cast<std::int64_t>(gridDomain2D_->rows) + extension;
            const std::int64_t activeColumns = 2 * static_cast<std::int64_t>(gridDomain2D_->columns) + extension;
            if (activeRows <= 0 || activeColumns <= 0 || static_cast<std::int64_t>(convention->infinityPixel) >= activeRows * activeColumns) {
                throw std::invalid_argument("TopographicConvention infinity pixel must belong to the active topographic domain.");
            }
            if (const auto* immersion = std::get_if<ComplementaryGridImmersion>(&convention->immersion)) {
                requireMatchingGrid(immersion->complementaryAdjacencies.minAdjacency, "Topographic minimum adjacency");
                requireMatchingGrid(immersion->complementaryAdjacencies.maxAdjacency, "Topographic maximum adjacency");
            }
        }
 
        initializeEmptyStorage(numPixels);
        nodeParent_ = std::move(nodeParent);
        firstChild_.assign(static_cast<size_t>(numNodeSlots), InvalidNode);
        nextSibling_.assign(static_cast<size_t>(numNodeSlots), InvalidNode);
        prevSibling_.assign(static_cast<size_t>(numNodeSlots), InvalidNode);
        lastChild_.assign(static_cast<size_t>(numNodeSlots), InvalidNode);
        numChildrenByNode_.assign(static_cast<size_t>(numNodeSlots), 0);
        alive_.assign(static_cast<size_t>(numNodeSlots), 1);
        freeNodeIds_.clear();
        properHead_.assign(static_cast<size_t>(numNodeSlots), InvalidPixel);
        properTail_.assign(static_cast<size_t>(numNodeSlots), InvalidPixel);
        properPartCardinalityByNode_.assign(static_cast<size_t>(numNodeSlots), 0);
        smallestNodeMap_ = std::move(smallestNodeMap);
        initializeProperPartStorage(numPixels);
 
        rootNodeId_ = root;
        numNodes_ = static_cast<int>(numNodeSlots);
 
        int selfParentedRoots = 0;
        for (NodeId nodeId = 0; nodeId < numNodeSlots; ++nodeId) {
            const NodeId parentId = nodeParent_[static_cast<size_t>(nodeId)];
            if (parentId == nodeId) {
                if (nodeId != root) {
                    throw std::invalid_argument("Native topology import found a detached self-parented non-root node.");
                }
                ++selfParentedRoots;
                continue;
            }
            if (parentId < 0 || parentId >= numNodeSlots) {
                throw std::invalid_argument("Native topology import found a parent outside the internal-node domain.");
            }
 
            prevSibling_[static_cast<size_t>(nodeId)] = lastChild_[static_cast<size_t>(parentId)];
            if (lastChild_[static_cast<size_t>(parentId)] == InvalidNode) {
                firstChild_[static_cast<size_t>(parentId)] = nodeId;
            } else {
                nextSibling_[static_cast<size_t>(lastChild_[static_cast<size_t>(parentId)])] = nodeId;
            }
            lastChild_[static_cast<size_t>(parentId)] = nodeId;
            ++numChildrenByNode_[static_cast<size_t>(parentId)];
        }
        if (selfParentedRoots != 1) {
            throw std::invalid_argument("Native topology import must encode exactly one self-parented root.");
        }
 
        rebuildProperPartLinksFromSmallestNodeMap();
        invalidateDfsIntervalCache();
        invalidateAllIterators();
        preservedExternalNodeIdOffset_.reset();
    }
 
    void initializeNativeTopology(std::span<const NodeId> nodeParent, std::span<const NodeId> smallestNodeMap, NodeId root,
                                  std::optional<GridDomain2D> gridDomain2D, MorphologicalTreeSemantics semantics) {
        if (nodeParent.size() > static_cast<std::size_t>(std::numeric_limits<NodeId>::max())) {
            throw std::invalid_argument("Native topology internal-node domain exceeds NodeId range.");
        }
        const int numNodeSlots = static_cast<int>(nodeParent.size());
        for (NodeId smallestNodeId : smallestNodeMap) {
            if (smallestNodeId < 0 || smallestNodeId >= numNodeSlots) {
                throw std::invalid_argument("Native topology import found a smallest-node-map entry outside the internal-node domain.");
            }
        }
        initializeNativeTopologyStorage(std::vector<NodeId>(nodeParent.begin(), nodeParent.end()),
                                        std::vector<NodeId>(smallestNodeMap.begin(), smallestNodeMap.end()), root, gridDomain2D, std::move(semantics));
        validateConnectedRootedTree();
    }
 
    void initializeValidatedNativeTopology(std::vector<NodeId>&& nodeParent, std::vector<NodeId>&& smallestNodeMap, NodeId root,
                                           std::optional<GridDomain2D> gridDomain2D, MorphologicalTreeSemantics semantics,
                                           detail::NativeTopologyProof&& topologyProof) {
        topologyProof.requireMatches(nodeParent.size(), smallestNodeMap.size(), root);
        initializeNativeTopologyStorage(std::move(nodeParent), std::move(smallestNodeMap), root, gridDomain2D, std::move(semantics));
#ifndef NDEBUG
        validateConnectedRootedTree();
#endif
    }
 
    MorphologicalTree(detail::MorphologicalTreeConstructionTag, std::vector<NodeId>&& nodeParent, std::vector<NodeId>&& smallestNodeMap, NodeId root,
                      std::optional<GridDomain2D> gridDomain2D, MorphologicalTreeSemantics semantics, detail::NativeTopologyProof&& topologyProof) {
        initializeValidatedNativeTopology(std::move(nodeParent), std::move(smallestNodeMap), root, gridDomain2D, std::move(semantics),
                                          std::move(topologyProof));
    }
 
    MorphologicalTree() = default;
 
    void moveCommittedStateFrom(MorphologicalTree&& other) {
        rootNodeId_ = std::exchange(other.rootNodeId_, InvalidNode);
        semantics_ = std::move(other.semantics_);
        other.semantics_ = MorphologicalTreeSemantics{};
        gridDomain2D_ = std::move(other.gridDomain2D_);
        other.gridDomain2D_.reset();
        numNodes_ = std::exchange(other.numNodes_, 0);
        preservedExternalNodeIdOffset_ = std::move(other.preservedExternalNodeIdOffset_);
        other.preservedExternalNodeIdOffset_.reset();
        editSessionOpen_ = false;
        editValidationStatistics_ = other.editValidationStatistics_;
 
        smallestNodeMap_ = std::move(other.smallestNodeMap_);
        nodeParent_ = std::move(other.nodeParent_);
        firstChild_ = std::move(other.firstChild_);
        nextSibling_ = std::move(other.nextSibling_);
        prevSibling_ = std::move(other.prevSibling_);
        lastChild_ = std::move(other.lastChild_);
        numChildrenByNode_ = std::move(other.numChildrenByNode_);
        alive_ = std::move(other.alive_);
        freeNodeIds_ = std::move(other.freeNodeIds_);
        properHead_ = std::move(other.properHead_);
        properTail_ = std::move(other.properTail_);
        properPartCardinalityByNode_ = std::move(other.properPartCardinalityByNode_);
        nextProperPart_ = std::move(other.nextProperPart_);
        prevProperPart_ = std::move(other.prevProperPart_);
 
        dfsIntervalCache_ = std::move(other.dfsIntervalCache_);
        other.dfsIntervalCache_ = {};
        lcaCache_.reset();
        other.lcaCache_.reset();
        nodeSupportMetadataCache_ = std::move(other.nodeSupportMetadataCache_);
        other.nodeSupportMetadataCache_ = {};
 
        nodeStructureVersion_ = other.nodeStructureVersion_;
        topologyVersion_ = other.topologyVersion_;
        properPartVersion_ = other.properPartVersion_;
        mutationVersion_ = other.mutationVersion_;
 
        // Iterators into the moved-from owner must not accidentally validate
        // merely because its previous token happened to be zero.
        ++other.nodeStructureVersion_;
        ++other.topologyVersion_;
        ++other.properPartVersion_;
        ++other.mutationVersion_;
    }
 
  public:
    MorphologicalTree(const MorphologicalTree&) = delete;
 
    MorphologicalTree& operator=(const MorphologicalTree&) = delete;
 
    MorphologicalTree(MorphologicalTree&& other) {
        other.requireNotEditing("MorphologicalTree move construction");
        moveCommittedStateFrom(std::move(other));
    }
 
    MorphologicalTree& operator=(MorphologicalTree&& other) {
        requireNotEditing("MorphologicalTree move assignment destination");
        other.requireNotEditing("MorphologicalTree move assignment source");
        if (this != &other) {
            const std::size_t previousNodeVersion = nodeStructureVersion_;
            const std::size_t previousTopologyVersion = topologyVersion_;
            const std::size_t previousProperPartVersion = properPartVersion_;
            const std::size_t previousMutationVersion = mutationVersion_;
            moveCommittedStateFrom(std::move(other));
            nodeStructureVersion_ = std::max(nodeStructureVersion_, previousNodeVersion) + 1;
            topologyVersion_ = std::max(topologyVersion_, previousTopologyVersion) + 1;
            properPartVersion_ = std::max(properPartVersion_, previousProperPartVersion) + 1;
            mutationVersion_ = std::max(mutationVersion_, previousMutationVersion) + 1;
        }
        return *this;
    }
 
    virtual ~MorphologicalTree() = default;
 
    MorphologicalTree(detail::MorphologicalTreeConstructionTag, std::span<const NodeId> nodeParent, std::span<const NodeId> smallestNodeMap, NodeId root,
                      int rows, int columns, MorphologicalTreeSemantics semantics) {
        initializeNativeTopology(nodeParent, smallestNodeMap, root, GridDomain2D{rows, columns}, std::move(semantics));
    }
 
    MorphologicalTree(detail::MorphologicalTreeConstructionTag, std::span<const NodeId> nodeParent, std::span<const NodeId> smallestNodeMap, NodeId root,
                      MorphologicalTreeSemantics semantics) {
        initializeNativeTopology(nodeParent, smallestNodeMap, root, std::nullopt, std::move(semantics));
    }
 
    MorphologicalTree clone() const {
        requireNotEditing("MorphologicalTree::clone");
        MorphologicalTree cloned;
        cloned.rootNodeId_ = rootNodeId_;
        cloned.semantics_ = semantics_;
        cloned.gridDomain2D_ = gridDomain2D_;
        cloned.numNodes_ = numNodes_;
        cloned.preservedExternalNodeIdOffset_ = preservedExternalNodeIdOffset_;
        cloned.editSessionOpen_ = false;
        cloned.editValidationStatistics_ = editValidationStatistics_;
        cloned.smallestNodeMap_ = smallestNodeMap_;
        cloned.nodeParent_ = nodeParent_;
        cloned.firstChild_ = firstChild_;
        cloned.nextSibling_ = nextSibling_;
        cloned.prevSibling_ = prevSibling_;
        cloned.lastChild_ = lastChild_;
        cloned.numChildrenByNode_ = numChildrenByNode_;
        cloned.alive_ = alive_;
        cloned.freeNodeIds_ = freeNodeIds_;
        cloned.properHead_ = properHead_;
        cloned.properTail_ = properTail_;
        cloned.properPartCardinalityByNode_ = properPartCardinalityByNode_;
        cloned.nextProperPart_ = nextProperPart_;
        cloned.prevProperPart_ = prevProperPart_;
        cloned.dfsIntervalCache_ = dfsIntervalCache_;
        cloned.lcaCache_.reset();
        cloned.nodeSupportMetadataCache_ = nodeSupportMetadataCache_;
        cloned.nodeStructureVersion_ = nodeStructureVersion_;
        cloned.topologyVersion_ = topologyVersion_;
        cloned.properPartVersion_ = properPartVersion_;
        cloned.mutationVersion_ = mutationVersion_;
        return cloned;
    }
 
    // Forward declarations for nested iterator/range types whose definitions stay at the end of the class.
    class AliveNodeIterator;
    class AliveNodeRange;
    class ChildrenIterator;
    class ChildrenRange;
    class ProperPartIterator;
    class ProperPartRange;
    class NodeSupportIterator;
    class NodeSupportRange;
    class PostOrderNodeIterator;
    class PostOrderNodeRange;
    class BreadthFirstNodeIterator;
    class BreadthFirstNodeRange;
    class AncestorNodeIterator;
    class AncestorNodeRange;
    class PathBetweenNodesIterator;
    class PathBetweenNodesRange;
    class SubtreeNodeIterator;
    class SubtreeNodeRange;
    class DescendantNodeRange;
 
  public:
    // ========================= Public methods ========================= //
 
    std::size_t getMutationVersion() const noexcept { return mutationVersion_; }
 
    void requireMutationVersion(std::size_t expectedVersion, const char* context) const {
        MMCFILTERS_CONTRACT_REQUIRE(mutationVersion_ == expectedVersion,
                                    throw std::logic_error(std::string(context) + " cannot be used after the referenced tree topology has changed."));
    }
 
    inline int numInternalNodeSlots() const { return static_cast<int>(nodeParent_.size()); }
 
    inline int numPixels() const { return static_cast<int>(smallestNodeMap_.size()); }
 
    inline int getNumHigraNodes() const {
        if (!preservedExternalNodeIdOffset_) {
            throw std::runtime_error("This tree does not preserve an imported Higra node-id space.");
        }
        return *preservedExternalNodeIdOffset_ + numInternalNodeSlots();
    }
 
    inline int getNodeIdSpaceSize(NodeIdSpace outputSpace) const {
        switch (outputSpace) {
        case NodeIdSpace::MorphologicalTree:
            return numInternalNodeSlots();
        case NodeIdSpace::Higra:
            return getNumHigraNodes();
        }
        throw std::runtime_error("Unknown NodeIdSpace.");
    }
 
    inline NodeId getHigraNodeId(NodeId nodeId) const noexcept {
        if (!preservedExternalNodeIdOffset_ || !isNode(nodeId) || !isAlive(nodeId)) {
            return InvalidNode;
        }
        return *preservedExternalNodeIdOffset_ + nodeId;
    }
 
    inline NodeId root() const { return rootNodeId_; }
 
    inline bool isNode(NodeId nodeId) const noexcept { return nodeId >= 0 && nodeId < static_cast<int>(nodeParent_.size()); }
 
    inline bool isPixel(PixelId pixel) const noexcept { return pixel >= 0 && pixel < static_cast<int>(smallestNodeMap_.size()); }
 
    inline bool isAlive(NodeId nodeId) const {
        if (!isNode(nodeId)) {
            return false;
        }
        const NodeId localId = nodeId;
        return localId >= 0 && localId < static_cast<NodeId>(nodeParent_.size()) && !isFreeSlot(localId) &&
               (nodeParent_[localId] != InvalidNode || localId == rootNodeId_);
    }
 
    inline bool isRoot(NodeId nodeId) const { return nodeId == root(); }
 
    inline int getNumFreeNodeSlots() const { return static_cast<int>(freeNodeIds_.size()); }
 
    inline int numLeafNodes() const {
        int count = 0;
        for (NodeId nodeId : aliveNodeIds()) {
            if (isLeaf(nodeId)) {
                ++count;
            }
        }
        return count;
    }
 
    inline int numChildren(NodeId nodeId) const {
        requireAliveNode(nodeId, "MorphologicalTree::numChildren");
        return numChildrenByNode_[nodeId];
    }
 
    inline int numDescendants(NodeId nodeId) const {
        requireAliveNode(nodeId, "MorphologicalTree::numDescendants");
        ensureDfsIntervalCache();
        return (dfsIntervalCache_.exitIndex[nodeId] - dfsIntervalCache_.entryIndex[nodeId] - 1) / 2;
    }
 
    inline int numSiblings(NodeId nodeId) const {
        requireAliveNode(nodeId, "MorphologicalTree::numSiblings");
        if (isRoot(nodeId)) {
            return 0;
        }
        const NodeId parentNodeId = parent(nodeId);
        return std::max(0, numChildren(parentNodeId) - 1);
    }
 
    inline int dfsEntryIndex(NodeId nodeId) const {
        requireAliveNode(nodeId, "MorphologicalTree::dfsEntryIndex");
        ensureDfsIntervalCache();
        return dfsIntervalCache_.entryIndex[nodeId];
    }
 
    inline int dfsExitIndex(NodeId nodeId) const {
        requireAliveNode(nodeId, "MorphologicalTree::dfsExitIndex");
        ensureDfsIntervalCache();
        return dfsIntervalCache_.exitIndex[nodeId];
    }
 
    inline NodeId getFirstChild(NodeId nodeId) const {
        requireAliveNode(nodeId, "MorphologicalTree::getFirstChild");
        return firstChild_[nodeId];
    }
 
    inline NodeId getNextSibling(NodeId nodeId) const {
        requireAliveNode(nodeId, "MorphologicalTree::getNextSibling");
        return nextSibling_[nodeId];
    }
 
    inline bool isLeaf(NodeId nodeId) const { return getFirstChild(nodeId) == InvalidNode; }
 
    inline int properPartCardinality(NodeId nodeId) const {
        requireAliveNode(nodeId, "MorphologicalTree::properPartCardinality");
        return properPartCardinalityByNode_[nodeId];
    }
 
    inline bool hasEmptyProperPart(NodeId nodeId) const { return properPartCardinality(nodeId) == 0; }
 
    inline bool hasChild(NodeId parentNodeId, NodeId childId) const {
        requireAliveNode(parentNodeId, "MorphologicalTree::hasChild");
        requireAliveNode(childId, "MorphologicalTree::hasChild");
        return parent(childId) == parentNodeId;
    }
 
    inline NodeId parent(NodeId nodeId) const {
        requireAliveNode(nodeId, "MorphologicalTree::parent");
        if (isRoot(nodeId)) {
            return nodeId;
        }
        const NodeId parentSlot = nodeParent_[nodeId];
        if (parentSlot == InvalidNode) {
            return InvalidNode;
        }
        if (parentSlot == nodeId) {
            return nodeId;
        }
        return parentSlot;
    }
 
    inline NodeId smallestNode(PixelId pixel) const {
        if constexpr (contract::validationsEnabled) {
            return isPixel(pixel) ? smallestNodeMap_[static_cast<size_t>(pixel)] : InvalidNode;
        }
        return smallestNodeMap_[static_cast<size_t>(pixel)];
    }
 
    [[nodiscard]] inline std::span<const NodeId> smallestNodeMap() const noexcept { return smallestNodeMap_; }
 
    std::vector<NodeId> leaves() const {
        std::vector<NodeId> leaves;
        if (rootNodeId_ == InvalidNode) {
            return leaves;
        }
        FastQueue<NodeId> s;
        s.push(this->rootNodeId_);
 
        while (!s.empty()) {
            const NodeId id = s.pop();
            if (numChildrenByNode_[id] == 0) {
                leaves.push_back(id);
            } else {
                for (NodeId c = firstChild_[id]; c != InvalidNode; c = nextSibling_[c]) {
                    s.push(c);
                }
            }
        }
        return leaves;
    }
 
    inline MorphologicalTreeKind kind() const noexcept { return semantics_.kind; }
 
    inline const MorphologicalTreeSemantics& semantics() const noexcept { return semantics_; }
 
    inline NodeAltitudeOrder nodeAltitudeOrder() const noexcept { return semantics_.nodeAltitudeOrder; }
 
    inline const MorphologicalTreeConstructionContext& constructionContext() const noexcept { return semantics_.constructionContext; }
 
    inline const SharedAdjacencyContext* sharedAdjacencyContext() const noexcept {
        return std::get_if<SharedAdjacencyContext>(&semantics_.constructionContext);
    }
 
    inline const SaturatedResidualContext* saturatedResidualContext() const noexcept {
        return std::get_if<SaturatedResidualContext>(&semantics_.constructionContext);
    }
 
    inline const TopographicConvention* topographicConvention() const noexcept {
        return std::get_if<TopographicConvention>(&semantics_.constructionContext);
    }
 
    inline int numNodes() const noexcept { return numNodes_; }
 
    inline bool isEditing() const noexcept { return editSessionOpen_; }
 
    [[nodiscard]] const TreeEditValidationStatistics& getEditValidationStatistics() const noexcept { return editValidationStatistics_; }
 
    inline void requireNotEditing(const char* context) const {
        MMCFILTERS_CONTRACT_REQUIRE(!isEditing(),
                                    throw std::logic_error(std::string(context) + " requires a committed MorphologicalTree; an edit session is still open."));
    }
 
    inline bool hasDetachedAliveNodes() const noexcept {
        if (numNodes_ == 0) {
            return false;
        }
        if (rootNodeId_ == InvalidNode || !isAlive(rootNodeId_)) {
            return true;
        }
 
        for (NodeId nodeId = 0; nodeId < static_cast<NodeId>(nodeParent_.size()); ++nodeId) {
            if (!isAlive(nodeId) || nodeId == rootNodeId_) {
                continue;
            }
            if (nodeParent_[nodeId] == nodeId) {
                return true;
            }
        }
        return false;
    }
 
    [[nodiscard]] inline bool hasGridDomain2D() const noexcept { return gridDomain2D_.has_value(); }
 
    [[nodiscard]] inline const std::optional<GridDomain2D>& gridDomain2D() const noexcept { return gridDomain2D_; }
 
    inline const GridDomain2D& requireGridDomain2D(const char* context) const {
        if (!gridDomain2D_) {
            throw std::invalid_argument(std::string(context) + " requires a regular 2D pixel domain.");
        }
        return *gridDomain2D_;
    }
 
    inline int numRows() const { return requireGridDomain2D("MorphologicalTree::numRows").rows; }
 
    inline int numColumns() const { return requireGridDomain2D("MorphologicalTree::numColumns").columns; }
 
    void validateConnectedRootedTree() const {
        if (numNodes_ <= 0) {
            throw std::runtime_error("Connected-tree validation requires at least one live node.");
        }
        if (!isAlive(rootNodeId_)) {
            throw std::runtime_error("Connected-tree validation requires a live root.");
        }
        if (nodeParent_[rootNodeId_] != rootNodeId_) {
            throw std::runtime_error("Connected-tree validation requires the root to point to itself.");
        }
 
        const int numSlots = numInternalNodeSlots();
        const int domainPixelCount = numPixels();
        std::vector<int> expectedChildrenByNode(static_cast<size_t>(numSlots), 0);
        std::vector<int> expectedProperPartCardinalityByNode(static_cast<size_t>(numSlots), 0);
        int aliveCount = 0;
 
        for (NodeId nodeId = 0; nodeId < numSlots; ++nodeId) {
            if (!isAlive(nodeId)) {
                continue;
            }
            ++aliveCount;
            const NodeId parentNodeId = nodeParent_[static_cast<size_t>(nodeId)];
            if (nodeId == rootNodeId_) {
                continue;
            }
            if (parentNodeId == InvalidNode) {
                throw std::runtime_error("Connected-tree validation found an alive node with no parent.");
            }
            if (parentNodeId == nodeId) {
                throw std::runtime_error("Connected-tree validation found a detached alive node.");
            }
            if (!isNode(parentNodeId) || !isAlive(parentNodeId)) {
                throw std::runtime_error("Connected-tree validation found an alive node whose parent is outside the alive node domain.");
            }
            expectedChildrenByNode[static_cast<size_t>(parentNodeId)] += 1;
        }
 
        if (aliveCount != numNodes_) {
            throw std::runtime_error("Connected-tree validation found numNodes() out of sync with the alive node slots.");
        }
 
        std::vector<uint8_t> seenAsChild(static_cast<size_t>(numSlots), 0);
        for (NodeId nodeId = 0; nodeId < numSlots; ++nodeId) {
            if (!isAlive(nodeId)) {
                continue;
            }
 
            int actualChildCount = 0;
            NodeId previousChildId = InvalidNode;
            for (NodeId childNodeId = firstChild_[static_cast<size_t>(nodeId)]; childNodeId != InvalidNode;
                 childNodeId = nextSibling_[static_cast<size_t>(childNodeId)]) {
                if (!isNode(childNodeId) || !isAlive(childNodeId)) {
                    throw std::runtime_error("Connected-tree validation found a child list that references a non-alive node.");
                }
                if (nodeParent_[static_cast<size_t>(childNodeId)] != nodeId) {
                    throw std::runtime_error("Connected-tree validation found a child whose parent pointer disagrees with the child list.");
                }
                if (prevSibling_[static_cast<size_t>(childNodeId)] != previousChildId) {
                    throw std::runtime_error("Connected-tree validation found broken previous-sibling links.");
                }
                if (seenAsChild[static_cast<size_t>(childNodeId)] != 0) {
                    throw std::runtime_error("Connected-tree validation found a node referenced by multiple child lists.");
                }
                seenAsChild[static_cast<size_t>(childNodeId)] = 1;
                previousChildId = childNodeId;
                ++actualChildCount;
                if (actualChildCount > aliveCount) {
                    throw std::runtime_error("Connected-tree validation found a cycle in a child list.");
                }
            }
 
            if (firstChild_[static_cast<size_t>(nodeId)] == InvalidNode) {
                if (lastChild_[static_cast<size_t>(nodeId)] != InvalidNode) {
                    throw std::runtime_error("Connected-tree validation found an empty child list with a non-empty tail pointer.");
                }
            } else {
                if (lastChild_[static_cast<size_t>(nodeId)] != previousChildId) {
                    throw std::runtime_error("Connected-tree validation found an incorrect last-child pointer.");
                }
                if (nextSibling_[static_cast<size_t>(previousChildId)] != InvalidNode) {
                    throw std::runtime_error("Connected-tree validation found a child list whose tail still points to a next sibling.");
                }
            }
 
            if (actualChildCount != numChildrenByNode_[static_cast<size_t>(nodeId)]) {
                throw std::runtime_error("Connected-tree validation found an incorrect child count cache.");
            }
            if (actualChildCount != expectedChildrenByNode[static_cast<size_t>(nodeId)]) {
                throw std::runtime_error("Connected-tree validation found child lists out of sync with parent pointers.");
            }
        }
 
        if (seenAsChild[static_cast<size_t>(rootNodeId_)] != 0) {
            throw std::runtime_error("Connected-tree validation found the root inside a child list.");
        }
        for (NodeId nodeId = 0; nodeId < numSlots; ++nodeId) {
            if (!isAlive(nodeId) || nodeId == rootNodeId_) {
                continue;
            }
            if (seenAsChild[static_cast<size_t>(nodeId)] == 0) {
                throw std::runtime_error("Connected-tree validation found a non-root alive node missing from the child lists.");
            }
        }
 
        std::vector<NodeId> traversal;
        traversal.reserve(static_cast<size_t>(aliveCount));
        traversal.push_back(rootNodeId_);
        for (size_t index = 0; index < traversal.size(); ++index) {
            const NodeId nodeId = traversal[index];
            for (NodeId childId = firstChild_[static_cast<size_t>(nodeId)]; childId != InvalidNode; childId = nextSibling_[static_cast<size_t>(childId)]) {
                traversal.push_back(childId);
            }
        }
        if (traversal.size() != static_cast<size_t>(aliveCount)) {
            throw std::runtime_error("Connected-tree validation found nodes outside the rooted component or a cycle.");
        }
 
        for (PixelId pixel = 0; pixel < domainPixelCount; ++pixel) {
            const NodeId smallestNodeId = smallestNodeMap_[static_cast<size_t>(pixel)];
            if (!isAlive(smallestNodeId)) {
                throw std::runtime_error("Connected-tree validation found a pixel without a live smallest node.");
            }
            expectedProperPartCardinalityByNode[static_cast<size_t>(smallestNodeId)] += 1;
        }
 
        std::vector<uint8_t> seenProperPart(static_cast<size_t>(domainPixelCount), 0);
        for (NodeId nodeId = 0; nodeId < numSlots; ++nodeId) {
            if (!isAlive(nodeId)) {
                continue;
            }
 
            int actualProperPartCardinality = 0;
            PixelId previousProperPartId = InvalidPixel;
            for (PixelId pixel = properHead_[static_cast<size_t>(nodeId)]; pixel != InvalidPixel;
                 pixel = nextProperPart_[static_cast<size_t>(pixel)]) {
                if (!isPixel(pixel)) {
                    throw std::runtime_error("Connected-tree validation found an invalid pixel id in a proper-part list.");
                }
                if (smallestNodeMap_[static_cast<size_t>(pixel)] != nodeId) {
                    throw std::runtime_error("Connected-tree validation found a proper-part list that disagrees with the smallest-node map.");
                }
                if (prevProperPart_[static_cast<size_t>(pixel)] != previousProperPartId) {
                    throw std::runtime_error("Connected-tree validation found broken previous-proper-part links.");
                }
                if (seenProperPart[static_cast<size_t>(pixel)] != 0) {
                    throw std::runtime_error("Connected-tree validation found a proper part referenced multiple times.");
                }
                seenProperPart[static_cast<size_t>(pixel)] = 1;
                previousProperPartId = pixel;
                ++actualProperPartCardinality;
                if (actualProperPartCardinality > domainPixelCount) {
                    throw std::runtime_error("Connected-tree validation found a cycle in a proper-part list.");
                }
            }
 
            if (properHead_[static_cast<size_t>(nodeId)] == InvalidPixel) {
                if (properTail_[static_cast<size_t>(nodeId)] != InvalidPixel) {
                    throw std::runtime_error("Connected-tree validation found an empty proper-part list with a non-empty tail pointer.");
                }
            } else {
                if (properTail_[static_cast<size_t>(nodeId)] != previousProperPartId) {
                    throw std::runtime_error("Connected-tree validation found an incorrect proper-part tail pointer.");
                }
                if (nextProperPart_[static_cast<size_t>(previousProperPartId)] != InvalidPixel) {
                    throw std::runtime_error("Connected-tree validation found a proper-part list whose tail still points forward.");
                }
            }
 
            if (actualProperPartCardinality != properPartCardinalityByNode_[static_cast<size_t>(nodeId)]) {
                throw std::runtime_error("Connected-tree validation found an incorrect direct proper-part count cache.");
            }
            if (actualProperPartCardinality != expectedProperPartCardinalityByNode[static_cast<size_t>(nodeId)]) {
                throw std::runtime_error("Connected-tree validation found proper-part lists out of sync with the smallest-node map.");
            }
        }
 
        for (PixelId pixel = 0; pixel < domainPixelCount; ++pixel) {
            if (seenProperPart[static_cast<size_t>(pixel)] == 0) {
                throw std::runtime_error("Connected-tree validation found a pixel missing from the proper-part lists.");
            }
        }
 
        std::vector<int> subtreeSupport = expectedProperPartCardinalityByNode;
        for (auto it = traversal.rbegin(); it != traversal.rend(); ++it) {
            const NodeId nodeId = *it;
            if (subtreeSupport[static_cast<size_t>(nodeId)] == 0) {
                throw std::runtime_error("Connected-tree validation found a live node whose subtree support is empty.");
            }
            if (nodeId != rootNodeId_) {
                const NodeId parentId = nodeParent_[static_cast<size_t>(nodeId)];
                subtreeSupport[static_cast<size_t>(parentId)] += subtreeSupport[static_cast<size_t>(nodeId)];
            }
        }
    }
 
    [[nodiscard]] TreeValidationResult validateConnectedRootedTreeResult() const noexcept {
        try {
            validateConnectedRootedTree();
            return {true, ""};
        } catch (const std::exception& ex) {
            return {false, ex.what()};
        } catch (...) {
            return {false, "Connected-tree validation failed with an unknown error."};
        }
    }
 
    [[nodiscard]] bool isTreeOfPartialPartitions() const {
        if (numNodes_ <= 0) {
            return false;
        }
        for (NodeId nodeId = 0; nodeId < numInternalNodeSlots(); ++nodeId) {
            if (isAlive(nodeId) && properPartCardinalityByNode_[static_cast<std::size_t>(nodeId)] == 0) {
                return false;
            }
        }
        return true;
    }
 
    void validateTreeOfPartialPartitions() const {
        validateConnectedRootedTree();
        for (NodeId nodeId = 0; nodeId < numInternalNodeSlots(); ++nodeId) {
            if (isAlive(nodeId) && properPartCardinalityByNode_[static_cast<std::size_t>(nodeId)] == 0) {
                throw std::runtime_error("Tree-of-partial-partitions validation found a node with an empty proper part: " + std::to_string(nodeId) + ".");
            }
        }
    }
 
    inline bool isAncestor(NodeId u, NodeId v) const {
        requireAliveNode(u, "MorphologicalTree::isAncestor");
        requireAliveNode(v, "MorphologicalTree::isAncestor");
        ensureDfsIntervalCache();
        const NodeId slotU = u;
        const NodeId slotV = v;
        return dfsIntervalCache_.entryIndex[slotU] <= dfsIntervalCache_.entryIndex[slotV] &&
               dfsIntervalCache_.exitIndex[slotU] >= dfsIntervalCache_.exitIndex[slotV];
    }
 
    inline bool isDescendant(NodeId u, NodeId v) const {
        requireAliveNode(u, "MorphologicalTree::isDescendant");
        requireAliveNode(v, "MorphologicalTree::isDescendant");
        ensureDfsIntervalCache();
        const NodeId slotU = u;
        const NodeId slotV = v;
        return dfsIntervalCache_.entryIndex[slotV] <= dfsIntervalCache_.entryIndex[slotU] &&
               dfsIntervalCache_.exitIndex[slotV] >= dfsIntervalCache_.exitIndex[slotU];
    }
    inline bool isComparable(NodeId u, NodeId v) const { return isAncestor(u, v) || isAncestor(v, u); }
 
    inline bool isStrictAncestor(NodeId u, NodeId v) const { return u != v && isAncestor(u, v); }
 
    inline bool isStrictDescendant(NodeId u, NodeId v) const { return u != v && isDescendant(u, v); }
 
    inline bool isStrictComparable(NodeId u, NodeId v) const { return isStrictAncestor(u, v) || isStrictAncestor(v, u); }
 
    inline NodeId lowestCommonAncestor(NodeId u, NodeId v) const {
        if (!isAlive(u) || !isAlive(v)) {
            return InvalidNode;
        }
        return lowestCommonAncestorEstablished(u, v);
    }
 
    [[nodiscard]] inline std::vector<NodeId>
    lowestCommonAncestors(std::span<const std::pair<NodeId, NodeId>> queries) const {
        std::vector<std::pair<NodeId, NodeId>> liveQueries;
        std::vector<std::size_t> liveQueryIndices;
        std::vector<NodeId> lcas(queries.size(), InvalidNode);
        liveQueries.reserve(queries.size());
        liveQueryIndices.reserve(queries.size());
        for (std::size_t queryIndex = 0; queryIndex < queries.size(); ++queryIndex) {
            const auto [first, second] = queries[queryIndex];
            if (!isAlive(first) || !isAlive(second)) {
                continue;
            }
            liveQueries.emplace_back(first, second);
            liveQueryIndices.push_back(queryIndex);
        }
 
        const std::vector<NodeId> liveLcas = lowestCommonAncestorsEstablished(liveQueries);
        for (std::size_t liveIndex = 0; liveIndex < liveLcas.size(); ++liveIndex) {
            lcas[liveQueryIndices[liveIndex]] = liveLcas[liveIndex];
        }
        return lcas;
    }
 
    inline AliveNodeRange aliveNodeIds() const { return AliveNodeRange(this, 0, numInternalNodeSlots(), nodeStructureVersion_); }
 
    inline ChildrenRange children(NodeId nodeId) const {
        requireAliveNode(nodeId, "MorphologicalTree::children");
        return ChildrenRange(this, firstChild_[nodeId], topologyVersion_);
    }
 
    inline ProperPartRange properPart(NodeId nodeId) const {
        requireAliveNode(nodeId, "MorphologicalTree::properPart");
        return ProperPartRange(this, properHead_[nodeId], properPartVersion_);
    }
 
    inline NodeSupportRange nodeSupport(NodeId nodeId) const {
        requireAliveNode(nodeId, "MorphologicalTree::nodeSupport");
        return NodeSupportRange(this, nodeId, topologyVersion_, properPartVersion_);
    }
 
    inline PostOrderNodeRange postOrder() const { return PostOrderNodeRange(this, root(), topologyVersion_); }
 
    inline PostOrderNodeRange postOrder(NodeId rootNodeId) const {
        requireAliveNode(rootNodeId, "MorphologicalTree::postOrder");
        return PostOrderNodeRange(this, rootNodeId, topologyVersion_);
    }
 
    inline BreadthFirstNodeRange breadthFirstTraversal() const { return BreadthFirstNodeRange(this, root(), topologyVersion_); }
 
    inline BreadthFirstNodeRange breadthFirstTraversal(NodeId rootNodeId) const {
        requireAliveNode(rootNodeId, "MorphologicalTree::breadthFirstTraversal");
        return BreadthFirstNodeRange(this, rootNodeId, topologyVersion_);
    }
 
    inline AncestorNodeRange ancestors(NodeId nodeId) const {
        requireAliveNode(nodeId, "MorphologicalTree::ancestors");
        return AncestorNodeRange(this, nodeId, topologyVersion_);
    }
 
    inline PathBetweenNodesRange getPathBetweenNodes(NodeId sourceNodeId, NodeId targetNodeId) const {
        requireAliveNode(sourceNodeId, "MorphologicalTree::getPathBetweenNodes");
        requireAliveNode(targetNodeId, "MorphologicalTree::getPathBetweenNodes");
        return PathBetweenNodesRange(this, sourceNodeId, targetNodeId, topologyVersion_);
    }
 
    inline SubtreeNodeRange subtreeNodes(NodeId nodeId) const {
        requireAliveNode(nodeId, "MorphologicalTree::subtreeNodes");
        return SubtreeNodeRange(this, nodeId, topologyVersion_);
    }
 
    inline DescendantNodeRange descendants(NodeId nodeId) const {
        requireAliveNode(nodeId, "MorphologicalTree::descendants");
        return DescendantNodeRange(this, nodeId, topologyVersion_);
    }
 
    [[nodiscard]] TreeEditor edit();
 
    inline void pruneNode(NodeId nodeId) {
        requireNotEditing("MorphologicalTree::pruneNode");
        requireAliveNonRootNode(nodeId, "MorphologicalTree::pruneNode");
        const NodeId parentNodeId = parent(nodeId);
        if (parentNodeId == InvalidNode || parentNodeId == nodeId) {
            throw std::invalid_argument("MorphologicalTree::pruneNode requires an attached non-root node.");
        }
        const NodeId parentSlotId = parentNodeId;
        std::size_t releaseCount = 0;
        for ([[maybe_unused]] NodeId subtreeNode : subtreeNodes(nodeId)) {
            ++releaseCount;
        }
        freeNodeIds_.reserve(freeNodeIds_.size() + releaseCount);
 
        const auto descendToDeepestLastChild = [this](NodeId startId) {
            NodeId currentId = startId;
            while (firstChild_[currentId] != InvalidNode) {
                currentId = lastChild_[currentId];
            }
            return currentId;
        };
 
        NodeId currentId = descendToDeepestLastChild(nodeId);
        while (true) {
            const NodeId currentParentSlotId = nodeParent_[currentId];
            mergeProperParts(parentSlotId, currentId);
            if (currentParentSlotId != InvalidNode && currentParentSlotId != currentId) {
                unlinkChildSlot(currentParentSlotId, currentId);
            }
            nodeParent_[currentId] = currentId;
            releaseNode(currentId);
 
            if (currentId == nodeId) {
                break;
            }
            currentId = firstChild_[currentParentSlotId] == InvalidNode ? currentParentSlotId : descendToDeepestLastChild(currentParentSlotId);
        }
    }
 
    inline void mergeNodeIntoParent(NodeId nodeId) {
        requireNotEditing("MorphologicalTree::mergeNodeIntoParent");
        requireAliveNonRootNode(nodeId, "MorphologicalTree::mergeNodeIntoParent");
        const NodeId parentNodeId = parent(nodeId);
        if (parentNodeId == InvalidNode || parentNodeId == nodeId) {
            throw std::invalid_argument("MorphologicalTree::mergeNodeIntoParent requires an attached non-root node.");
        }
        const NodeId parentSlotId = parentNodeId;
        const NodeId nodeSlotId = nodeId;
        freeNodeIds_.reserve(freeNodeIds_.size() + 1);
 
        spliceProperPartsSlots(parentSlotId, nodeSlotId);
        spliceChildrenSlots(parentSlotId, nodeSlotId, ChildSplicePolicy::ReplaceSourceSlotWhenDirectChild);
        nodeParent_[nodeSlotId] = nodeSlotId;
        bumpProperPartVersion();
        releaseNode(nodeId);
    }
 
  private:
    // ========================= Internal classes ========================= //
    class LCAEulerRMQ {
      private:
        std::vector<NodeId> euler_;
        std::vector<int> depth_;
        std::vector<int> firstOccurrence_;
        std::vector<int> log2_;
        std::vector<int> sparseTable_;
        int sparseTableStride_ = 0;
        const MorphologicalTree* tree_ = nullptr;
 
        void depthFirstTraversal(NodeId nodeId, int currentDepth) {
            std::vector<std::pair<NodeId, NodeId>> stack;
            stack.emplace_back(nodeId, tree_->firstChild_[nodeId]);
            firstOccurrence_[static_cast<size_t>(nodeId)] = static_cast<int>(euler_.size());
            euler_.push_back(nodeId);
            depth_.push_back(currentDepth);
            while (!stack.empty()) {
                const NodeId child = stack.back().second;
                if (child == InvalidNode) {
                    stack.pop_back();
                    if (!stack.empty()) {
                        euler_.push_back(stack.back().first);
                        depth_.push_back(currentDepth + static_cast<int>(stack.size()) - 1);
                    }
                } else {
                    stack.back().second = tree_->nextSibling_[child];
                    firstOccurrence_[static_cast<size_t>(child)] = static_cast<int>(euler_.size());
                    euler_.push_back(child);
                    depth_.push_back(currentDepth + static_cast<int>(stack.size()));
                    stack.emplace_back(child, tree_->firstChild_[child]);
                }
            }
        }
 
        void buildSparseTable() {
            const int n = static_cast<int>(depth_.size());
            if (n == 0) {
                log2_.clear();
                sparseTable_.clear();
                sparseTableStride_ = 0;
                return;
            }
 
            log2_.assign(static_cast<size_t>(n + 1), 0);
            for (int i = 2; i <= n; ++i) {
                log2_[static_cast<size_t>(i)] = log2_[static_cast<size_t>(i / 2)] + 1;
            }
 
            sparseTableStride_ = 1;
            while ((1 << sparseTableStride_) <= n) {
                ++sparseTableStride_;
            }
 
            sparseTable_.assign(static_cast<size_t>(n) * static_cast<size_t>(sparseTableStride_), 0);
 
            for (int i = 0; i < n; ++i) {
                sparseTable_[sparseTableIndex(i, 0)] = i;
            }
 
            for (int j = 1; j < sparseTableStride_; ++j) {
                const int blockSize = 1 << j;
                const int halfBlock = blockSize >> 1;
                for (int i = 0; i + blockSize <= n; ++i) {
                    const int leftIndex = sparseTable_[sparseTableIndex(i, j - 1)];
                    const int rightIndex = sparseTable_[sparseTableIndex(i + halfBlock, j - 1)];
                    sparseTable_[sparseTableIndex(i, j)] =
                        depth_[static_cast<size_t>(leftIndex)] <= depth_[static_cast<size_t>(rightIndex)] ? leftIndex : rightIndex;
                }
            }
        }
 
        std::size_t sparseTableIndex(int row, int column) const {
            return static_cast<std::size_t>(row) * static_cast<std::size_t>(sparseTableStride_) + static_cast<std::size_t>(column);
        }
 
        int rmq(int left, int right) const {
            const int length = right - left + 1;
            const int logLength = log2_[static_cast<size_t>(length)];
            const int leftIndex = sparseTable_[sparseTableIndex(left, logLength)];
            const int rightIndex = sparseTable_[sparseTableIndex(right - (1 << logLength) + 1, logLength)];
            return depth_[static_cast<size_t>(leftIndex)] <= depth_[static_cast<size_t>(rightIndex)] ? leftIndex : rightIndex;
        }
 
      public:
        explicit LCAEulerRMQ(const MorphologicalTree* tree) : tree_(tree) {
            if (tree_ == nullptr || tree_->root() == InvalidNode) {
                return;
            }
 
            euler_.reserve(static_cast<size_t>(std::max(0, 2 * tree_->numNodes() - 1)));
            depth_.reserve(static_cast<size_t>(std::max(0, 2 * tree_->numNodes() - 1)));
            firstOccurrence_.assign(static_cast<size_t>(tree_->numInternalNodeSlots()), -1);
 
            depthFirstTraversal(tree_->root(), 0);
            buildSparseTable();
        }
 
        NodeId findLowestCommonAncestor(NodeId u, NodeId v) const {
            const int firstU = firstOccurrence_[static_cast<size_t>(u)];
            const int firstV = firstOccurrence_[static_cast<size_t>(v)];
            if (firstU < 0 || firstV < 0) {
                return InvalidNode;
            }
 
            int left = firstU;
            int right = firstV;
            if (left > right) {
                std::swap(left, right);
            }
            return euler_[static_cast<size_t>(rmq(left, right))];
        }
    };
 
  public:
    // ========================= Internal classes and iterators ========================= //
 
    class AliveNodeIterator {
      private:
        const MorphologicalTree* T_ = nullptr;
        NodeId current_ = InvalidNode;
        NodeId end_ = InvalidNode;
        std::size_t expectedVersion_ = 0;
 
        void settle_() {
            if (!T_) {
                current_ = InvalidNode;
                return;
            }
            T_->checkNodeIteratorVersion(expectedVersion_);
            while (current_ != InvalidNode && current_ < end_ && !T_->isAlive(current_)) {
                ++current_;
            }
            if (current_ >= end_) {
                current_ = InvalidNode;
            }
        }
 
      public:
        using iterator_category = std::input_iterator_tag;
        using value_type = NodeId;
        using difference_type = std::ptrdiff_t;
        using pointer = const NodeId*;
        using reference = const NodeId&;
 
        AliveNodeIterator() = default;
 
        AliveNodeIterator(const MorphologicalTree* tree, NodeId current, NodeId end, std::size_t expectedVersion)
            : T_(tree), current_(current), end_(end), expectedVersion_(expectedVersion) {
            settle_();
        }
 
        AliveNodeIterator& operator++() {
            T_->checkNodeIteratorVersion(expectedVersion_);
            if (current_ != InvalidNode) {
                ++current_;
                settle_();
            }
            return *this;
        }
 
        NodeId operator*() const {
            T_->checkNodeIteratorVersion(expectedVersion_);
            return current_;
        }
 
        bool operator==(const AliveNodeIterator& other) const { return current_ == other.current_; }
 
        bool operator!=(const AliveNodeIterator& other) const { return !(*this == other); }
    };
 
    class AliveNodeRange {
      private:
        const MorphologicalTree* T_ = nullptr;
        NodeId begin_ = InvalidNode;
        NodeId end_ = InvalidNode;
        std::size_t expectedVersion_ = 0;
 
      public:
        AliveNodeRange() = default;
 
        AliveNodeRange(const MorphologicalTree* tree, NodeId begin, NodeId end, std::size_t expectedVersion)
            : T_(tree), begin_(begin), end_(end), expectedVersion_(expectedVersion) {}
 
        AliveNodeIterator begin() const { return AliveNodeIterator(T_, begin_, end_, expectedVersion_); }
 
        AliveNodeIterator end() const { return AliveNodeIterator(T_, InvalidNode, end_, expectedVersion_); }
    };
 
    class ChildrenIterator {
      private:
        const MorphologicalTree* T_ = nullptr;
        NodeId currentLocal_ = InvalidNode;
        std::size_t expectedVersion_ = 0;
 
      public:
        using iterator_category = std::forward_iterator_tag;
        using value_type = NodeId;
        using difference_type = std::ptrdiff_t;
        using pointer = const NodeId*;
        using reference = const NodeId&;
 
        ChildrenIterator() = default;
 
        ChildrenIterator(const MorphologicalTree* tree, NodeId currentLocal, std::size_t expectedVersion)
            : T_(tree), currentLocal_(currentLocal), expectedVersion_(expectedVersion) {}
 
        ChildrenIterator& operator++() {
            T_->checkTopologyIteratorVersion(expectedVersion_);
            if (T_ && currentLocal_ != InvalidNode) {
                currentLocal_ = T_->nextSibling_[currentLocal_];
            }
            return *this;
        }
 
        NodeId operator*() const {
            T_->checkTopologyIteratorVersion(expectedVersion_);
            return currentLocal_;
        }
 
        bool operator==(const ChildrenIterator& other) const { return currentLocal_ == other.currentLocal_; }
 
        bool operator!=(const ChildrenIterator& other) const { return !(*this == other); }
    };
 
    class ChildrenRange {
      private:
        const MorphologicalTree* T_ = nullptr;
        NodeId firstLocal_ = InvalidNode;
        std::size_t expectedVersion_ = 0;
 
      public:
        ChildrenRange() = default;
 
        ChildrenRange(const MorphologicalTree* tree, NodeId firstLocal, std::size_t expectedVersion)
            : T_(tree), firstLocal_(firstLocal), expectedVersion_(expectedVersion) {}
 
        ChildrenIterator begin() const { return ChildrenIterator(T_, firstLocal_, expectedVersion_); }
 
        ChildrenIterator end() const { return ChildrenIterator(T_, InvalidNode, expectedVersion_); }
    };
 
    class ProperPartIterator {
      private:
        const MorphologicalTree* T_ = nullptr;
        PixelId currentProperPart_ = InvalidPixel;
        std::size_t expectedVersion_ = 0;
 
      public:
        using iterator_category = std::forward_iterator_tag;
        using value_type = PixelId;
        using difference_type = std::ptrdiff_t;
        using pointer = const PixelId*;
        using reference = const PixelId&;
 
        ProperPartIterator() = default;
 
        ProperPartIterator(const MorphologicalTree* tree, PixelId currentProperPart, std::size_t expectedVersion)
            : T_(tree), currentProperPart_(currentProperPart), expectedVersion_(expectedVersion) {}
 
        ProperPartIterator& operator++() {
            T_->checkProperPartIteratorVersion(expectedVersion_);
            if (T_ && currentProperPart_ != InvalidPixel) {
                currentProperPart_ = T_->nextProperPart_[currentProperPart_];
            }
            return *this;
        }
 
        PixelId operator*() const {
            T_->checkProperPartIteratorVersion(expectedVersion_);
            return currentProperPart_;
        }
 
        bool operator==(const ProperPartIterator& other) const { return currentProperPart_ == other.currentProperPart_; }
 
        bool operator!=(const ProperPartIterator& other) const { return !(*this == other); }
    };
 
    class ProperPartRange {
      private:
        const MorphologicalTree* T_ = nullptr;
        PixelId firstProperPart_ = InvalidPixel;
        std::size_t expectedVersion_ = 0;
 
      public:
        ProperPartRange() = default;
 
        ProperPartRange(const MorphologicalTree* tree, PixelId firstProperPart, std::size_t expectedVersion)
            : T_(tree), firstProperPart_(firstProperPart), expectedVersion_(expectedVersion) {}
 
        ProperPartIterator begin() const { return ProperPartIterator(T_, firstProperPart_, expectedVersion_); }
 
        ProperPartIterator end() const { return ProperPartIterator(T_, InvalidPixel, expectedVersion_); }
    };
 
    class NodeSupportIterator {
      private:
        const MorphologicalTree* T_ = nullptr;
        std::vector<NodeId> nodeStack_;
        PixelId currentProperPart_ = InvalidPixel;
        std::size_t expectedTopologyVersion_ = 0;
        std::size_t expectedProperPartVersion_ = 0;
 
        void checkVersions() const {
            T_->checkTopologyIteratorVersion(expectedTopologyVersion_);
            T_->checkProperPartIteratorVersion(expectedProperPartVersion_);
        }
 
        void pushChildren(NodeId nodeId) {
            std::vector<NodeId> children;
            for (NodeId childId : T_->children(nodeId)) {
                children.push_back(childId);
            }
            for (auto it = children.rbegin(); it != children.rend(); ++it) {
                nodeStack_.push_back(*it);
            }
        }
 
        void settle() {
            checkVersions();
            while (currentProperPart_ == InvalidPixel && !nodeStack_.empty()) {
                const NodeId nodeId = nodeStack_.back();
                nodeStack_.pop_back();
                pushChildren(nodeId);
                currentProperPart_ = T_->properHead_[static_cast<size_t>(nodeId)];
            }
        }
 
      public:
        using iterator_category = std::forward_iterator_tag;
        using value_type = PixelId;
        using difference_type = std::ptrdiff_t;
        using pointer = const PixelId*;
        using reference = const PixelId&;
 
        NodeSupportIterator() = default;
 
        NodeSupportIterator(const MorphologicalTree* tree, NodeId rootNodeId, std::size_t expectedTopologyVersion, std::size_t expectedProperPartVersion)
            : T_(tree), expectedTopologyVersion_(expectedTopologyVersion), expectedProperPartVersion_(expectedProperPartVersion) {
            if (T_ && rootNodeId != InvalidNode) {
                nodeStack_.push_back(rootNodeId);
                settle();
            }
        }
 
        NodeSupportIterator& operator++() {
            checkVersions();
            if (currentProperPart_ != InvalidPixel) {
                currentProperPart_ = T_->nextProperPart_[static_cast<size_t>(currentProperPart_)];
            }
            settle();
            return *this;
        }
 
        PixelId operator*() const {
            checkVersions();
            return currentProperPart_;
        }
 
        bool operator==(const NodeSupportIterator& other) const { return currentProperPart_ == other.currentProperPart_; }
 
        bool operator!=(const NodeSupportIterator& other) const { return !(*this == other); }
    };
 
    class NodeSupportRange {
      private:
        const MorphologicalTree* T_ = nullptr;
        NodeId rootNodeId_ = InvalidNode;
        std::size_t expectedTopologyVersion_ = 0;
        std::size_t expectedProperPartVersion_ = 0;
 
      public:
        NodeSupportRange() = default;
 
        NodeSupportRange(const MorphologicalTree* tree, NodeId rootNodeId, std::size_t expectedTopologyVersion, std::size_t expectedProperPartVersion)
            : T_(tree), rootNodeId_(rootNodeId), expectedTopologyVersion_(expectedTopologyVersion), expectedProperPartVersion_(expectedProperPartVersion) {}
 
        NodeSupportIterator begin() const { return NodeSupportIterator(T_, rootNodeId_, expectedTopologyVersion_, expectedProperPartVersion_); }
 
        NodeSupportIterator end() const { return NodeSupportIterator(); }
    };
 
    class PostOrderNodeIterator {
      private:
        struct Item {
            NodeId id;
            bool expanded;
        };
        const MorphologicalTree* T_ = nullptr;
        std::vector<Item> stack_;
        NodeId current_ = InvalidNode;
        std::size_t expectedVersion_ = 0;
 
        void settle_() {
            T_->checkTopologyIteratorVersion(expectedVersion_);
            while (!stack_.empty()) {
                Item& top = stack_.back();
                if (!top.expanded) {
                    top.expanded = true;
                    // Reverse sibling links preserve child order without a temporary vector.
                    for (NodeId child = T_->lastChild_[top.id]; child != InvalidNode; child = T_->prevSibling_[child]) {
                        stack_.push_back({child, false});
                    }
                } else {
                    current_ = top.id;
                    return;
                }
            }
            current_ = InvalidNode;
        }
 
      public:
        using iterator_category = std::input_iterator_tag;
        using value_type = NodeId;
        using difference_type = std::ptrdiff_t;
        using pointer = const NodeId*;
        using reference = const NodeId&;
 
        PostOrderNodeIterator() = default;
 
        PostOrderNodeIterator(const MorphologicalTree* tree, NodeId rootNodeId, std::size_t expectedVersion) : T_(tree), expectedVersion_(expectedVersion) {
            if (T_ && rootNodeId != InvalidNode) {
                stack_.push_back({rootNodeId, false});
                settle_();
            }
        }
 
        PostOrderNodeIterator& operator++() {
            T_->checkTopologyIteratorVersion(expectedVersion_);
            if (!stack_.empty()) {
                stack_.pop_back();
                settle_();
            }
            return *this;
        }
 
        NodeId operator*() const {
            T_->checkTopologyIteratorVersion(expectedVersion_);
            return current_;
        }
 
        bool operator==(const PostOrderNodeIterator& other) const { return current_ == other.current_; }
 
        bool operator!=(const PostOrderNodeIterator& other) const { return !(*this == other); }
    };
 
    class PostOrderNodeRange {
      private:
        const MorphologicalTree* T_ = nullptr;
        NodeId rootNodeId_ = InvalidNode;
        std::size_t expectedVersion_ = 0;
 
      public:
        PostOrderNodeRange() = default;
 
        PostOrderNodeRange(const MorphologicalTree* tree, NodeId rootNodeId, std::size_t expectedVersion)
            : T_(tree), rootNodeId_(rootNodeId), expectedVersion_(expectedVersion) {}
 
        PostOrderNodeIterator begin() const { return PostOrderNodeIterator(T_, rootNodeId_, expectedVersion_); }
 
        PostOrderNodeIterator end() const { return PostOrderNodeIterator(); }
    };
 
    class BreadthFirstNodeIterator {
      private:
        const MorphologicalTree* T_ = nullptr;
        FastQueue<NodeId> queue_;
        std::size_t expectedVersion_ = 0;
 
      public:
        using iterator_category = std::input_iterator_tag;
        using value_type = NodeId;
        using difference_type = std::ptrdiff_t;
        using pointer = const NodeId*;
        using reference = const NodeId&;
 
        BreadthFirstNodeIterator() = default;
 
        BreadthFirstNodeIterator(const MorphologicalTree* tree, NodeId rootNodeId, std::size_t expectedVersion) : T_(tree), expectedVersion_(expectedVersion) {
            if (T_ && rootNodeId != InvalidNode) {
                queue_.push(rootNodeId);
            }
        }
 
        BreadthFirstNodeIterator& operator++() {
            T_->checkTopologyIteratorVersion(expectedVersion_);
            if (!queue_.empty()) {
                NodeId current = queue_.pop();
                for (NodeId childId : T_->children(current)) {
                    queue_.push(childId);
                }
            }
            return *this;
        }
 
        NodeId operator*() const {
            T_->checkTopologyIteratorVersion(expectedVersion_);
            return queue_.front();
        }
 
        bool operator==(const BreadthFirstNodeIterator& other) const { return queue_.empty() == other.queue_.empty(); }
 
        bool operator!=(const BreadthFirstNodeIterator& other) const { return !(*this == other); }
    };
 
    class BreadthFirstNodeRange {
      private:
        const MorphologicalTree* T_ = nullptr;
        NodeId rootNodeId_ = InvalidNode;
        std::size_t expectedVersion_ = 0;
 
      public:
        BreadthFirstNodeRange() = default;
 
        BreadthFirstNodeRange(const MorphologicalTree* tree, NodeId rootNodeId, std::size_t expectedVersion)
            : T_(tree), rootNodeId_(rootNodeId), expectedVersion_(expectedVersion) {}
 
        BreadthFirstNodeIterator begin() const { return BreadthFirstNodeIterator(T_, rootNodeId_, expectedVersion_); }
 
        BreadthFirstNodeIterator end() const { return BreadthFirstNodeIterator(); }
    };
 
    class AncestorNodeIterator {
      private:
        const MorphologicalTree* T_ = nullptr;
        NodeId current_ = InvalidNode;
        std::size_t expectedVersion_ = 0;
 
      public:
        using iterator_category = std::input_iterator_tag;
        using value_type = NodeId;
        using difference_type = std::ptrdiff_t;
        using pointer = const NodeId*;
        using reference = const NodeId&;
 
        AncestorNodeIterator() = default;
 
        AncestorNodeIterator(const MorphologicalTree* tree, NodeId current, std::size_t expectedVersion)
            : T_(tree), current_(current), expectedVersion_(expectedVersion) {}
 
        AncestorNodeIterator& operator++() {
            T_->checkTopologyIteratorVersion(expectedVersion_);
            if (T_ && current_ != InvalidNode) {
                if (T_->isRoot(current_)) {
                    current_ = InvalidNode;
                } else {
                    current_ = T_->parent(current_);
                }
            }
            return *this;
        }
 
        NodeId operator*() const {
            T_->checkTopologyIteratorVersion(expectedVersion_);
            return current_;
        }
 
        bool operator==(const AncestorNodeIterator& other) const { return current_ == other.current_; }
 
        bool operator!=(const AncestorNodeIterator& other) const { return !(*this == other); }
    };
 
    class AncestorNodeRange {
      private:
        const MorphologicalTree* T_ = nullptr;
        NodeId start_ = InvalidNode;
        std::size_t expectedVersion_ = 0;
 
      public:
        AncestorNodeRange() = default;
 
        AncestorNodeRange(const MorphologicalTree* tree, NodeId start, std::size_t expectedVersion)
            : T_(tree), start_(start), expectedVersion_(expectedVersion) {}
 
        AncestorNodeIterator begin() const { return AncestorNodeIterator(T_, start_, expectedVersion_); }
 
        AncestorNodeIterator end() const { return AncestorNodeIterator(); }
    };
 
    class PathBetweenNodesIterator {
      private:
        const MorphologicalTree* T_ = nullptr;
        const std::vector<NodeId>* path_ = nullptr;
        std::size_t index_ = 0;
        std::size_t expectedVersion_ = 0;
 
      public:
        using iterator_category = std::input_iterator_tag;
        using value_type = NodeId;
        using difference_type = std::ptrdiff_t;
        using pointer = const NodeId*;
        using reference = const NodeId&;
 
        PathBetweenNodesIterator() = default;
 
        PathBetweenNodesIterator(const MorphologicalTree* tree, const std::vector<NodeId>* path, std::size_t index, std::size_t expectedVersion)
            : T_(tree), path_(path), index_(index), expectedVersion_(expectedVersion) {}
 
        PathBetweenNodesIterator& operator++() {
            T_->checkTopologyIteratorVersion(expectedVersion_);
            if (path_ && index_ < path_->size()) {
                ++index_;
            }
            return *this;
        }
 
        NodeId operator*() const {
            T_->checkTopologyIteratorVersion(expectedVersion_);
            return (*path_)[index_];
        }
 
        bool operator==(const PathBetweenNodesIterator& other) const { return path_ == other.path_ && index_ == other.index_; }
 
        bool operator!=(const PathBetweenNodesIterator& other) const { return !(*this == other); }
    };
 
    class PathBetweenNodesRange {
      private:
        const MorphologicalTree* T_ = nullptr;
        std::vector<NodeId> path_;
        std::size_t expectedVersion_ = 0;
 
        static std::vector<NodeId> buildPath(const MorphologicalTree* tree, NodeId sourceNodeId, NodeId targetNodeId) {
            if (tree == nullptr || !tree->isAlive(sourceNodeId) || !tree->isAlive(targetNodeId)) {
                return {};
            }
 
            auto componentAnchor = [tree](NodeId nodeId) {
                NodeId currentNodeId = nodeId;
                while (true) {
                    const NodeId parentNodeId = tree->parent(currentNodeId);
                    if (parentNodeId == InvalidNode || parentNodeId == currentNodeId) {
                        return currentNodeId;
                    }
                    currentNodeId = parentNodeId;
                }
            };
 
            if (componentAnchor(sourceNodeId) != componentAnchor(targetNodeId)) {
                return {};
            }
 
            const NodeId lcaNodeId = tree->lowestCommonAncestor(sourceNodeId, targetNodeId);
            if (lcaNodeId == InvalidNode) {
                return {};
            }
 
            std::vector<NodeId> path;
            for (NodeId currentNodeId = sourceNodeId;; currentNodeId = tree->parent(currentNodeId)) {
                path.push_back(currentNodeId);
                if (currentNodeId == lcaNodeId) {
                    break;
                }
 
                const NodeId parentNodeId = tree->parent(currentNodeId);
                if (parentNodeId == InvalidNode || parentNodeId == currentNodeId) {
                    return {};
                }
            }
 
            std::vector<NodeId> descendingTail;
            for (NodeId currentNodeId = targetNodeId; currentNodeId != lcaNodeId; currentNodeId = tree->parent(currentNodeId)) {
                descendingTail.push_back(currentNodeId);
 
                const NodeId parentNodeId = tree->parent(currentNodeId);
                if (parentNodeId == InvalidNode || parentNodeId == currentNodeId) {
                    return {};
                }
            }
 
            path.insert(path.end(), descendingTail.rbegin(), descendingTail.rend());
            return path;
        }
 
      public:
        PathBetweenNodesRange() = default;
 
        PathBetweenNodesRange(const MorphologicalTree* tree, NodeId sourceNodeId, NodeId targetNodeId, std::size_t expectedVersion)
            : T_(tree), path_(buildPath(tree, sourceNodeId, targetNodeId)), expectedVersion_(expectedVersion) {}
 
        PathBetweenNodesIterator begin() const { return PathBetweenNodesIterator(T_, &path_, 0, expectedVersion_); }
 
        PathBetweenNodesIterator end() const { return PathBetweenNodesIterator(T_, &path_, path_.size(), expectedVersion_); }
    };
 
    class SubtreeNodeIterator {
      private:
        const MorphologicalTree* T_ = nullptr;
        std::vector<NodeId> stack_;
        std::size_t expectedVersion_ = 0;
 
      public:
        using iterator_category = std::input_iterator_tag;
        using value_type = NodeId;
        using difference_type = std::ptrdiff_t;
        using pointer = const NodeId*;
        using reference = const NodeId&;
 
        SubtreeNodeIterator() = default;
 
        SubtreeNodeIterator(const MorphologicalTree* tree, NodeId rootNodeId, std::size_t expectedVersion) : T_(tree), expectedVersion_(expectedVersion) {
            if (T_ && rootNodeId != InvalidNode) {
                stack_.push_back(rootNodeId);
            }
        }
 
        SubtreeNodeIterator& operator++() {
            T_->checkTopologyIteratorVersion(expectedVersion_);
            if (!stack_.empty()) {
                NodeId current = stack_.back();
                stack_.pop_back();
                // Push in reverse sibling order so the first child is visited next.
                for (NodeId child = T_->lastChild_[current]; child != InvalidNode; child = T_->prevSibling_[child]) {
                    stack_.push_back(child);
                }
            }
            return *this;
        }
 
        NodeId operator*() const {
            T_->checkTopologyIteratorVersion(expectedVersion_);
            return stack_.back();
        }
 
        bool operator==(const SubtreeNodeIterator& other) const { return stack_.empty() == other.stack_.empty(); }
 
        bool operator!=(const SubtreeNodeIterator& other) const { return !(*this == other); }
    };
 
    class SubtreeNodeRange {
      private:
        const MorphologicalTree* T_ = nullptr;
        NodeId rootNodeId_ = InvalidNode;
        std::size_t expectedVersion_ = 0;
 
      public:
        SubtreeNodeRange() = default;
 
        SubtreeNodeRange(const MorphologicalTree* tree, NodeId rootNodeId, std::size_t expectedVersion)
            : T_(tree), rootNodeId_(rootNodeId), expectedVersion_(expectedVersion) {}
 
        SubtreeNodeIterator begin() const { return SubtreeNodeIterator(T_, rootNodeId_, expectedVersion_); }
 
        SubtreeNodeIterator end() const { return SubtreeNodeIterator(); }
    };
 
    class DescendantNodeRange {
      private:
        const MorphologicalTree* T_ = nullptr;
        NodeId rootNodeId_ = InvalidNode;
        std::size_t expectedVersion_ = 0;
 
      public:
        DescendantNodeRange() = default;
 
        DescendantNodeRange(const MorphologicalTree* tree, NodeId rootNodeId, std::size_t expectedVersion)
            : T_(tree), rootNodeId_(rootNodeId), expectedVersion_(expectedVersion) {}
 
        SubtreeNodeIterator begin() const {
            auto it = SubtreeNodeIterator(T_, rootNodeId_, expectedVersion_);
            ++it;
            return it;
        }
 
        SubtreeNodeIterator end() const { return SubtreeNodeIterator(); }
    };
};
 
} // namespace mmcfilters