TreeAltitudeAlgorithms.hpp

#pragma once
 
#include "../utils/Altitude.hpp"
#include "MorphologicalTree.hpp"
#include "detail/HigraExportLayoutDetail.hpp"
 
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <span>
#include <sstream>
#include <stdexcept>
#include <type_traits>
#include <utility>
#include <vector>
 
namespace mmcfilters {
 
class TreeAltitudeAlgorithms {
public:
    template<AltitudeValue T>
    static void validateAltitudeBufferShape(const MorphologicalTree& tree, std::span<const T> altitude) {
        if (altitude.size() != static_cast<std::size_t>(tree.getNumInternalNodeSlots())) {
            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);
            if (!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>) {
            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>) {
            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 getAltitude(std::span<const T> altitude, NodeId nodeId) {
        if (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 AltitudeDiff<T> getNodeResidue(const MorphologicalTree& tree, std::span<const T> altitude, NodeId nodeId) {
        validateAltitudeBufferShape(tree, altitude);
        if (!tree.isAlive(nodeId)) {
            throw std::invalid_argument("Node residue requires a live internal NodeId.");
        }
        const NodeId parentNodeId = tree.getNodeParent(nodeId);
        if (parentNodeId == InvalidNode || parentNodeId == nodeId) {
            return static_cast<AltitudeDiff<T>>(getAltitude(altitude, nodeId));
        }
        return static_cast<AltitudeDiff<T>>(getAltitude(altitude, nodeId)) -
               static_cast<AltitudeDiff<T>>(getAltitude(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) {
        validateAltitudeBufferShape(tree, altitude);
        for (NodeId nodeId : tree.getAliveNodeIds()) {
            (void)requireUInt8AltitudeValue(getAltitude(altitude, nodeId), nodeId, context);
        }
    }
 
    template<AltitudeValue T>
    [[nodiscard]] static ImagePtr<T> reconstructImage(const MorphologicalTree& tree, std::span<const T> altitude, const char* context = "TreeAltitudeAlgorithms::reconstructImage") {
        (void)context;
        tree.requireNotEditing(context);
        validateAltitudeBufferShape(tree, altitude);
        ImagePtr<T> image = Image<T>::create(tree.getNumRowsOfImage(), tree.getNumColsOfImage());
        auto imgBuffer = image->rawData();
        for (int pixelId = 0; pixelId < tree.getNumTotalProperParts(); ++pixelId) {
            const NodeId nodeId = tree.getProperPartOwner(pixelId);
            imgBuffer[static_cast<std::size_t>(pixelId)] = getAltitude(altitude, nodeId);
        }
        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 NodeId numLeaves = layout.numLeaves;
        const NodeId 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)] = getAltitude(altitude, oldNodeId);
        }
 
        for (NodeId leafIndex = 0; leafIndex < numLeaves; ++leafIndex) {
            const NodeId properPart = layout.properParts[static_cast<std::size_t>(leafIndex)];
            const NodeId ownerNodeId = tree.getProperPartOwner(properPart);
            if (ownerNodeId == InvalidNode || !tree.isAlive(ownerNodeId)) {
                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>(ownerNodeId)];
            exportedAltitude[static_cast<std::size_t>(leafIndex)] = getAltitude(altitude, ownerNodeId);
        }
 
        for (NodeId oldNodeId : layout.sortedNodes) {
            const NodeId newNodeId = layout.nodeToHigra[static_cast<std::size_t>(oldNodeId)];
            const NodeId oldParentNodeId = tree.getNodeParent(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 validateMonotoneAltitude(const MorphologicalTree& tree, std::span<const T> altitude) {
        validateAltitudeBufferShape(tree, altitude);
        const MorphologicalTreeKind treeType = tree.getTreeType();
        bool increasingTowardLeaves = false;
        switch (treeType) {
            case MorphologicalTreeKind::MAX_TREE:
                increasingTowardLeaves = true;
                break;
            case MorphologicalTreeKind::MIN_TREE:
                increasingTowardLeaves = false;
                break;
            case MorphologicalTreeKind::TREE_OF_SHAPES:
            case MorphologicalTreeKind::SELF_DUAL_RESIDUAL_TREE:
                return;
            default:
                throw std::invalid_argument("Unsupported tree type for monotone altitude validation.");
        }
 
        for (NodeId nodeId : tree.getAliveNodeIds()) {
            if (tree.isRoot(nodeId)) {
                continue;
            }
 
            const NodeId parentNodeId = tree.getNodeParent(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 (increasingTowardLeaves) {
                if (getAltitude(altitude, parentNodeId) > getAltitude(altitude, nodeId)) {
                    throw std::runtime_error("Max-tree altitude buffer must be non-decreasing from parent to child.");
                }
            } else if (getAltitude(altitude, parentNodeId) < getAltitude(altitude, nodeId)) {
                throw std::runtime_error("Min-tree altitude buffer must be non-increasing from parent to child.");
            }
        }
    }
};
 
} // namespace mmcfilters