MaxDistComputer.hpp

#pragma once
 
#include "AttributeComputerDomain.hpp"
#include "AttributeComputerFamily.hpp"
#include "../detail/AttributeKernelSupport.hpp"
#include "../../contours/ContoursComputedIncrementally.hpp"
#include "../../trees/MorphologicalTree.hpp"
#include "../../trees/TreeAltitudeAlgorithms.hpp"
#include "../../utils/Common.hpp"
 
#include "detail/maxdist/EdtDIFT.hpp"
 
#include <algorithm>
#include <array>
#include <cmath>
#include <concepts>
#include <cstddef>
#include <cstdint>
#include <span>
#include <stdexcept>
#include <string>
#include <string_view>
#include <type_traits>
#include <vector>
 
namespace mmcfilters::attributes::computers::detail {
 
class MaxDistAttributeKernel {
public:
    static void requireSupportedTreeKind(const MorphologicalTree& tree)
    {
        switch (tree.getTreeType()) {
            case MorphologicalTreeKind::MAX_TREE:
            case MorphologicalTreeKind::MIN_TREE:
                return;
            case MorphologicalTreeKind::TREE_OF_SHAPES:
                throw std::invalid_argument(
                    "MAX_DIST is currently defined only for MAX_TREE and MIN_TREE; TREE_OF_SHAPES is not supported.");
            case MorphologicalTreeKind::SELF_DUAL_RESIDUAL_TREE:
                throw std::invalid_argument(
                    "MAX_DIST is currently defined only for MAX_TREE and MIN_TREE; SELF_DUAL_RESIDUAL_TREE is not supported.");
        }
 
        throw std::invalid_argument("MAX_DIST received an unsupported morphological tree kind.");
    }
 
    template<std::floating_point Real, AltitudeValue T>
    [[nodiscard]] static std::vector<Real> compute(
        const MorphologicalTree& tree,
        std::span<const T> altitude)
    {
        requireSupportedTreeKind(tree);
        if (!tree.hasAdjacencyRelation()) {
            throw std::invalid_argument("MAX_DIST requires an adjacency relation.");
        }
        TreeAltitudeAlgorithms::validateAltitudeBufferShape(tree, altitude);
        validateFiniteAltitude(altitude);
 
        std::vector<Real> maxDist(static_cast<std::size_t>(tree.getNumInternalNodeSlots()), Real{0});
        maxdist::EdtDIFT edtDIFT(tree.getNumRowsOfImage(), tree.getNumColsOfImage());
 
        const ContourDeltaStore contourDeltas =
            ::mmcfilters::ContoursComputedIncrementally::extractContourDeltas(tree);
        std::vector<std::vector<int>> contours(static_cast<std::size_t>(tree.getNumInternalNodeSlots()));
        const std::size_t totalPixels = static_cast<std::size_t>(
            tree.getNumRowsOfImage() * tree.getNumColsOfImage());
        std::vector<std::uint8_t> removalMark(totalPixels, 0);
        std::vector<std::uint8_t> contourAdditionMark(totalPixels, 0);
 
        std::vector<NodeId> sortedNodes = sortedNodesByAltitude(tree, altitude);
        std::size_t groupBegin = 0;
        while (groupBegin < sortedNodes.size()) {
            std::size_t groupEnd = groupBegin + 1;
            while (groupEnd < sortedNodes.size() &&
                   sameAltitude(altitude, sortedNodes[groupBegin], sortedNodes[groupEnd])) {
                ++groupEnd;
            }
 
            processLevel(
                tree,
                std::span<const NodeId>(
                    sortedNodes.data() + static_cast<std::ptrdiff_t>(groupBegin),
                    groupEnd - groupBegin),
                contourDeltas,
                edtDIFT,
                contours,
                removalMark,
                contourAdditionMark,
                maxDist);
 
            groupBegin = groupEnd;
        }
 
        return maxDist;
    }
 
private:
    using ContourDeltaStore = ::mmcfilters::ContoursComputedIncrementally::LocalContourDeltas;
 
    template<AltitudeValue T>
    static bool isFiniteAltitude(T value) noexcept
    {
        if constexpr (std::is_floating_point_v<T>) {
            return std::isfinite(value);
        }
        return true;
    }
 
    template<AltitudeValue T>
    static void validateFiniteAltitude(std::span<const T> altitude)
    {
        for (std::size_t index = 0; index < altitude.size(); ++index) {
            if (!isFiniteAltitude(altitude[index])) {
                throw std::invalid_argument(
                    "MAX_DIST requires finite altitude values; node " +
                    std::to_string(index) + " has a non-finite altitude.");
            }
        }
    }
 
    static void markPixels(std::span<const int> pixels, std::vector<std::uint8_t>& marks)
    {
        for (int pixelId : pixels) {
            marks[static_cast<std::size_t>(pixelId)] = 1;
        }
    }
 
    static void clearPixelMarks(std::span<const int> pixels, std::vector<std::uint8_t>& marks)
    {
        for (int pixelId : pixels) {
            marks[static_cast<std::size_t>(pixelId)] = 0;
        }
    }
 
    template<AltitudeValue T>
    [[nodiscard]] static std::vector<NodeId> sortedNodesByAltitude(
        const MorphologicalTree& tree,
        std::span<const T> altitude)
    {
        std::vector<NodeId> nodes;
        nodes.reserve(static_cast<std::size_t>(tree.getNumNodes()));
        for (NodeId nodeId : tree.getPostOrderNodes()) {
            nodes.push_back(nodeId);
        }
 
        std::stable_sort(nodes.begin(), nodes.end(), [&](NodeId lhs, NodeId rhs) {
            const T lhsAltitude = TreeAltitudeAlgorithms::getAltitude(altitude, lhs);
            const T rhsAltitude = TreeAltitudeAlgorithms::getAltitude(altitude, rhs);
            return tree.getTreeType() == MorphologicalTreeKind::MAX_TREE
                ? lhsAltitude > rhsAltitude
                : lhsAltitude < rhsAltitude;
        });
 
        return nodes;
    }
 
    template<AltitudeValue T>
    static bool sameAltitude(
        std::span<const T> altitude,
        NodeId lhs,
        NodeId rhs)
    {
        return TreeAltitudeAlgorithms::getAltitude(altitude, lhs) ==
               TreeAltitudeAlgorithms::getAltitude(altitude, rhs);
    }
 
    template <std::floating_point Real>
    static void processLevel(
        const MorphologicalTree& tree,
        std::span<const NodeId> nodes,
        const ContourDeltaStore& contourDeltas,
        maxdist::EdtDIFT& edtDIFT,
        std::vector<std::vector<int>>& contours,
        std::vector<std::uint8_t>& removalMark,
        std::vector<std::uint8_t>& contourAdditionMark,
        std::vector<Real>& maxDist)
    {
        if (nodes.empty()) {
            return;
        }
 
        std::vector<int> toRemove;
        toRemove.reserve(64);
        for (NodeId nodeId : nodes) {
            std::vector<int>& nodeContour = contours[static_cast<std::size_t>(nodeId)];
            nodeContour.clear();
 
            const auto removals = contourDeltas.removals(nodeId);
            toRemove.clear();
            toRemove.insert(toRemove.end(), removals.begin(), removals.end());
            markPixels(removals, removalMark);
 
            const auto additions = contourDeltas.additions(nodeId);
            std::size_t reserveSize = additions.size();
            for (NodeId childNodeId : tree.getChildren(nodeId)) {
                reserveSize += contours[static_cast<std::size_t>(childNodeId)].size();
            }
            nodeContour.reserve(reserveSize);
 
            for (NodeId childNodeId : tree.getChildren(nodeId)) {
                std::vector<int>& childContour = contours[static_cast<std::size_t>(childNodeId)];
                for (int pixelId : childContour) {
                    if (!removalMark[static_cast<std::size_t>(pixelId)]) {
                        nodeContour.push_back(pixelId);
                    }
                }
                std::vector<int>().swap(childContour);
            }
            clearPixelMarks(removals, removalMark);
 
            if (!toRemove.empty()) {
                edtDIFT.treeRemoval(toRemove);
            }
 
            markPixels(additions, contourAdditionMark);
            for (int pixelId : tree.getProperParts(nodeId)) {
                edtDIFT.addPixelToBinaryImage(pixelId);
 
                if (contourAdditionMark[static_cast<std::size_t>(pixelId)]) {
                    nodeContour.push_back(pixelId);
                    edtDIFT.seed(pixelId);
                }
                else {
                    edtDIFT.open(pixelId);
                    edtDIFT.insertNeighborsPQueue(pixelId);
                }
            }
            clearPixelMarks(additions, contourAdditionMark);
        }
 
        // All nodes at the same altitude must enter the binary image before the
        // distance transform propagates. This preserves the mathematical
        // per-level schedule without requiring fixed 0..255 buckets.
        edtDIFT.run();
 
        for (NodeId nodeId : nodes) {
            maxDist[static_cast<std::size_t>(nodeId)] =
                static_cast<Real>(edtDIFT.maxBedt(contours[static_cast<std::size_t>(nodeId)]));
        }
    }
};
 
} // namespace mmcfilters::attributes::computers::detail
 
namespace mmcfilters::attributes::computers {
 
class MaxDistComputer {
public:
    static constexpr std::string_view familyName = "max-dist";
 
    static constexpr AttributeComputerFamily family = AttributeComputerFamily::MaxDist;
 
    static constexpr AttributeComputerDomain domain = AttributeComputerDomain::Altitude;
 
    inline static constexpr std::array<Attribute, 1> producedAttributes{MAX_DIST};
 
    static void requireSupportedTreeKind(const MorphologicalTree& tree)
    {
        detail::MaxDistAttributeKernel::requireSupportedTreeKind(tree);
    }
 
    template<std::floating_point Real, AltitudeValue T>
    static void compute(const AltitudeAttributeComputeContext<Real, T>& context)
    {
        requireAttributeBufferShape(context.tree, context.buffer, context.attrNames);
        if (!requestsAttribute(context.requestedAttributes, MAX_DIST)) {
            return;
        }
 
        const std::vector<Real> maxDist =
            detail::MaxDistAttributeKernel::compute<Real>(context.tree, context.altitude);
        for (NodeId nodeId : context.tree.getAliveNodeIds()) {
            context.buffer[context.attrNames.linearIndex(nodeId, MAX_DIST)] =
                maxDist[static_cast<std::size_t>(nodeId)];
        }
    }
 
    template <std::floating_point Real, AltitudeValue T>
    static void computeUnitRows(const AltitudeUnitAttributeComputeContext<Real, T>& context)
    {
        requireUnitAttributeBufferShape(context.tree, context.unitProperParts, context.buffer, context.attrNames);
        TreeAltitudeAlgorithms::validateAltitudeBufferShape(context.tree, context.altitude);
        if (!requestsAttribute(context.requestedAttributes, MAX_DIST)) {
            return;
        }
 
        for (NodeId leafIndex = 0; leafIndex < static_cast<NodeId>(context.unitProperParts.size()); ++leafIndex) {
            context.buffer[context.attrNames.linearIndex(leafIndex, MAX_DIST)] = Real{0};
        }
    }
};
 
} // namespace mmcfilters::attributes::computers