TreeAltitudeAlgorithms.hpp

#pragma once
 
#include "../utils/Altitude.hpp"
#include "../utils/Contract.hpp"
#include "../utils/Image.hpp"
#include "MorphologicalTree.hpp"
#include "detail/HigraExportLayoutDetail.hpp"
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <span>
#include <sstream>
#include <stdexcept>
#include <type_traits>
#include <utility>
#include <vector>
 
namespace mmcfilters {
 
namespace detail::tree_altitude {
 
template <class Contribution>
    requires(std::is_arithmetic_v<Contribution> && !std::is_same_v<std::remove_cv_t<Contribution>, bool>)
[[nodiscard]] std::vector<Contribution> reconstructNodeContributionValues(const MorphologicalTree& tree,
                                                                           std::span<const Contribution> nodeContributions,
                                                                           const char* context) {
    tree.requireNotEditing(context);
    MMCFILTERS_CONTRACT_REQUIRE(
        nodeContributions.size() == static_cast<std::size_t>(tree.numInternalNodeSlots()),
        throw std::invalid_argument(std::string(context) + " nodeContributions size must match the internal node slot count."));
    if constexpr (std::is_floating_point_v<Contribution> && contract::validationsEnabled) {
        for (Contribution contribution : nodeContributions) {
            if (!std::isfinite(contribution)) {
                throw std::invalid_argument(std::string(context) + " requires finite nodeContributions.");
            }
        }
    }
 
    std::vector<Contribution> accumulated(static_cast<std::size_t>(tree.numInternalNodeSlots()), Contribution{});
    const NodeId root = tree.root();
    accumulated[static_cast<std::size_t>(root)] = nodeContributions[static_cast<std::size_t>(root)];
 
    std::vector<NodeId> pending{root};
    while (!pending.empty()) {
        const NodeId nodeId = pending.back();
        pending.pop_back();
        for (NodeId childId : tree.children(nodeId)) {
            accumulated[static_cast<std::size_t>(childId)] =
                accumulated[static_cast<std::size_t>(nodeId)] + nodeContributions[static_cast<std::size_t>(childId)];
            pending.push_back(childId);
        }
    }
 
    std::vector<Contribution> pixels(static_cast<std::size_t>(tree.numPixels()), Contribution{});
    for (PixelId pixel = 0; pixel < tree.numPixels(); ++pixel) {
        pixels[static_cast<std::size_t>(pixel)] = accumulated[static_cast<std::size_t>(tree.smallestNode(pixel))];
    }
    return pixels;
}
 
} // namespace detail::tree_altitude
 
class TreeAltitudeAlgorithms {
  public:
    template <AltitudeValue T> static void validateNodeAltitudeBufferShape(const MorphologicalTree& tree, std::span<const T> altitude) {
        MMCFILTERS_CONTRACT_REQUIRE(altitude.size() == static_cast<std::size_t>(tree.numInternalNodeSlots()),
                                    throw std::runtime_error("Altitude buffer size must match the dense internal-node domain."));
    }
 
    template <AltitudeValue T> static void validateFiniteAltitudeValue(T altitude, std::size_t index, const char* context) {
        if constexpr (std::is_floating_point_v<T>) {
            const long double level = static_cast<long double>(altitude);
            MMCFILTERS_CONTRACT_REQUIRE(std::isfinite(level), {
                std::ostringstream oss;
                oss << context << " requires finite floating-point altitudes; value at index " << index << " is " << level << ".";
                throw std::invalid_argument(oss.str());
            });
        }
    }
 
    template <AltitudeValue T> static void validateFiniteAltitudeValues(std::span<const T> altitude, const char* context) {
        if constexpr (std::is_floating_point_v<T> && contract::validationsEnabled) {
            for (std::size_t index = 0; index < altitude.size(); ++index) {
                validateFiniteAltitudeValue(altitude[index], index, context);
            }
        }
    }
 
    template <AltitudeValue T> static void validateFiniteImageAltitudes(const ImagePtr<T>& image, const char* context) {
        if constexpr (std::is_floating_point_v<T> && contract::validationsEnabled) {
            if (!image) {
                throw std::invalid_argument("Image altitude validation requires a non-null image.");
            }
            validateFiniteAltitudeValues(std::span<const T>(image->rawData(), static_cast<std::size_t>(image->getSize())), context);
        }
    }
 
    template <AltitudeValue T> [[nodiscard]] static T nodeAltitude(std::span<const T> altitude, NodeId nodeId) {
        MMCFILTERS_CONTRACT_REQUIRE(nodeId >= 0 && static_cast<std::size_t>(nodeId) < altitude.size(),
                                    throw std::invalid_argument("Altitude access requires a valid internal NodeId."));
        return altitude[static_cast<std::size_t>(nodeId)];
    }
 
    template <AltitudeValue T> [[nodiscard]] static AltitudeDifference<T> nodeResidue(const MorphologicalTree& tree, std::span<const T> altitude, NodeId nodeId) {
        validateNodeAltitudeBufferShape(tree, altitude);
        MMCFILTERS_CONTRACT_REQUIRE(tree.isAlive(nodeId), throw std::invalid_argument("Node residue requires a live internal NodeId."));
        const NodeId parentNodeId = tree.parent(nodeId);
        if (parentNodeId == InvalidNode || parentNodeId == nodeId) {
            return static_cast<AltitudeDifference<T>>(nodeAltitude(altitude, nodeId));
        }
        return static_cast<AltitudeDifference<T>>(nodeAltitude(altitude, nodeId)) - static_cast<AltitudeDifference<T>>(nodeAltitude(altitude, parentNodeId));
    }
 
    template <AltitudeValue T> [[nodiscard]] static std::uint8_t requireUInt8AltitudeValue(T altitude, NodeId nodeId, const char* context) {
        const long double level = static_cast<long double>(altitude);
        if constexpr (std::is_floating_point_v<T>) {
            if (!std::isfinite(level)) {
                std::ostringstream oss;
                oss << context << " requires finite node altitudes in the uint8 domain [0, 255]; node " << nodeId << " has altitude " << level << ".";
                throw std::invalid_argument(oss.str());
            }
        }
        if (level < 0.0L || level > 255.0L) {
            std::ostringstream oss;
            oss << context << " requires node altitudes in the uint8 domain [0, 255]; node " << nodeId << " has altitude " << level << ".";
            throw std::invalid_argument(oss.str());
        }
        return static_cast<std::uint8_t>(altitude);
    }
 
    template <AltitudeValue T> static void validateUInt8AltitudeDomain(const MorphologicalTree& tree, std::span<const T> altitude, const char* context) {
        validateNodeAltitudeBufferShape(tree, altitude);
        for (NodeId nodeId : tree.aliveNodeIds()) {
            (void)requireUInt8AltitudeValue(nodeAltitude(altitude, nodeId), nodeId, context);
        }
    }
 
    template <AltitudeValue T>
    [[nodiscard]] static ImagePtr<T> reconstructFromNodeAltitudes(const MorphologicalTree& tree, std::span<const T> altitude,
                                                      const char* context = "TreeAltitudeAlgorithms::reconstructFromNodeAltitudes") {
        (void)context;
        tree.requireNotEditing(context);
        validateNodeAltitudeBufferShape(tree, altitude);
        ImagePtr<T> image = Image<T>::create(tree.numRows(), tree.numColumns());
        auto imgBuffer = image->rawData();
        for (PixelId pixel = 0; pixel < tree.numPixels(); ++pixel) {
            const NodeId nodeId = tree.smallestNode(pixel);
            imgBuffer[static_cast<std::size_t>(pixel)] = nodeAltitude(altitude, nodeId);
        }
        return image;
    }
 
    template <class Contribution>
        requires(std::is_arithmetic_v<Contribution> && !std::is_same_v<std::remove_cv_t<Contribution>, bool>)
    [[nodiscard]] static ImagePtr<Contribution> reconstructFromNodeContributions(
        const MorphologicalTree& tree, std::span<const Contribution> nodeContributions,
        const char* context = "TreeAltitudeAlgorithms::reconstructFromNodeContributions") {
        std::vector<Contribution> pixelValues = detail::tree_altitude::reconstructNodeContributionValues(tree, nodeContributions, context);
        ImagePtr<Contribution> image = Image<Contribution>::create(tree.numRows(), tree.numColumns());
        Contribution* pixels = image->rawData();
        std::copy(pixelValues.begin(), pixelValues.end(), pixels);
        return image;
    }
 
    template <AltitudeValue T>
    [[nodiscard]] static std::pair<std::vector<NodeId>, std::vector<T>> exportHigraHierarchy(const MorphologicalTree& tree, std::span<const T> altitude) {
        tree.requireNotEditing("TreeAltitudeAlgorithms::exportHigraHierarchy");
        const detail::ExportedHigraLayout layout = detail::computeExportedHigraLayout(tree, altitude);
        const int numLeaves = layout.numLeaves;
        const int numVertices = layout.numVertices;
 
        std::vector<NodeId> parent(static_cast<std::size_t>(numVertices), InvalidNode);
        std::vector<T> exportedAltitude(static_cast<std::size_t>(numVertices), T{});
 
        for (NodeId oldNodeId : layout.sortedNodes) {
            const NodeId newNodeId = layout.nodeToHigra[static_cast<std::size_t>(oldNodeId)];
            exportedAltitude[static_cast<std::size_t>(newNodeId)] = nodeAltitude(altitude, oldNodeId);
        }
 
        for (NodeId leafIndex = 0; leafIndex < numLeaves; ++leafIndex) {
            const PixelId pixel = layout.properParts[static_cast<std::size_t>(leafIndex)];
            const NodeId smallestNodeId = tree.smallestNode(pixel);
            if (smallestNodeId == InvalidNode || !tree.isAlive(smallestNodeId)) {
                throw std::runtime_error("Each proper part must belong to one alive node when exporting a compact Higra hierarchy.");
            }
            parent[static_cast<std::size_t>(leafIndex)] = layout.nodeToHigra[static_cast<std::size_t>(smallestNodeId)];
            exportedAltitude[static_cast<std::size_t>(leafIndex)] = nodeAltitude(altitude, smallestNodeId);
        }
 
        for (NodeId oldNodeId : layout.sortedNodes) {
            const NodeId newNodeId = layout.nodeToHigra[static_cast<std::size_t>(oldNodeId)];
            const NodeId oldParentNodeId = tree.parent(oldNodeId);
            parent[static_cast<std::size_t>(newNodeId)] =
                oldParentNodeId == oldNodeId ? newNodeId : layout.nodeToHigra[static_cast<std::size_t>(oldParentNodeId)];
        }
 
        return {std::move(parent), std::move(exportedAltitude)};
    }
 
    template <AltitudeValue T> static void validateMonotoneNodeAltitudes(const MorphologicalTree& tree, std::span<const T> altitude) {
        validateNodeAltitudeBufferShape(tree, altitude);
        const NodeAltitudeOrder nodeAltitudeOrder = tree.nodeAltitudeOrder();
        bool increasingFromRoot = false;
        switch (nodeAltitudeOrder) {
        case NodeAltitudeOrder::Increasing:
            increasingFromRoot = true;
            break;
        case NodeAltitudeOrder::Decreasing:
            increasingFromRoot = false;
            break;
        case NodeAltitudeOrder::Unconstrained:
            return;
        }
 
        for (NodeId nodeId : tree.aliveNodeIds()) {
            if (tree.isRoot(nodeId)) {
                continue;
            }
 
            const NodeId parentNodeId = tree.parent(nodeId);
            if (parentNodeId == InvalidNode || !tree.isAlive(parentNodeId)) {
                throw std::runtime_error("Monotonic validation requires every alive non-root node to have an alive parent.");
            }
 
            if (increasingFromRoot) {
                if (nodeAltitude(altitude, parentNodeId) >= nodeAltitude(altitude, nodeId)) {
                    throw std::runtime_error("Hierarchy altitude buffer must be strictly increasing from parent to child.");
                }
            } else if (nodeAltitude(altitude, parentNodeId) <= nodeAltitude(altitude, nodeId)) {
                throw std::runtime_error("Hierarchy altitude buffer must be strictly decreasing from parent to child.");
            }
        }
    }
};
 
} // namespace mmcfilters