HierarchySaliencyMapValidation.hpp

#pragma once
 
#include "../../utils/Altitude.hpp"
#include "../../utils/Common.hpp"
#include "../../utils/RegularGridAdjacency2D.hpp"
#include "../MorphologicalTree.hpp"
#include "../ValuedMorphologicalTree.hpp"
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <iterator>
#include <sstream>
#include <span>
#include <stdexcept>
#include <string>
#include <type_traits>
#include <utility>
#include <vector>
 
namespace mmcfilters {
 
namespace detail {
 
inline void requireCommittedRootedHierarchy(const MorphologicalTree& tree, const char* context) {
    if (tree.isEditing()) {
        throw std::invalid_argument(std::string(context) + " requires a committed tree; an edit session is still open.");
    }
    if (tree.root() == InvalidNode || !tree.isAlive(tree.root())) {
        throw std::invalid_argument(std::string(context) + " requires a non-empty connected rooted tree.");
    }
    if (tree.parent(tree.root()) != tree.root()) {
        throw std::invalid_argument(std::string(context) + " requires the root to point to itself.");
    }
 
    const std::size_t slotCount = static_cast<std::size_t>(tree.numInternalNodeSlots());
    std::vector<std::uint8_t> visited(slotCount, 0);
    std::vector<NodeId> stack{tree.root()};
    visited[static_cast<std::size_t>(tree.root())] = 1;
    std::size_t visitedCount = 0;
    std::size_t traversedEdges = 0;
 
    while (!stack.empty()) {
        const NodeId nodeId = stack.back();
        stack.pop_back();
        ++visitedCount;
 
        for (NodeId childId : tree.children(nodeId)) {
            ++traversedEdges;
            if (traversedEdges >= slotCount || !tree.isAlive(childId) || tree.parent(childId) != nodeId) {
                throw std::invalid_argument(std::string(context) + " requires consistent parent-child relations.");
            }
            const std::size_t childIndex = static_cast<std::size_t>(childId);
            if (visited[childIndex] != 0) {
                throw std::invalid_argument(std::string(context) + " requires an acyclic rooted tree.");
            }
            visited[childIndex] = 1;
            stack.push_back(childId);
        }
    }
 
    std::size_t aliveCount = 0;
    for (NodeId nodeId : tree.aliveNodeIds()) {
        ++aliveCount;
        if (visited[static_cast<std::size_t>(nodeId)] == 0) {
            throw std::invalid_argument(std::string(context) + " requires every live node to be connected to the root.");
        }
    }
    if (visitedCount != aliveCount) {
        throw std::invalid_argument(std::string(context) + " requires one connected rooted tree.");
    }
}
 
} // namespace detail
 
enum class HierarchyValuationPolicy {
    AllowLevelCollapse,
    RequireStrictHierarchy,
};
 
enum class HierarchyValuationRangePolicy {
    AllowAnyFinite,
    RequireNonNegative,
};
 
enum class HierarchyConnectivityPolicy {
    AssumeConnected,
    ValidateConnected,
};
 
class HierarchySaliencyMapValidation {
  private:
    template <class Value>
    static void validateValuationValue(const Value& value, NodeId nodeId, const char* context, HierarchyValuationRangePolicy rangePolicy) {
        if constexpr (std::is_floating_point_v<Value>) {
            if (!std::isfinite(value)) {
                std::ostringstream oss;
                oss << context << " requires finite valuation values; node " << nodeId << " has value " << value << ".";
                throw std::invalid_argument(oss.str());
            }
        }
        if (rangePolicy == HierarchyValuationRangePolicy::RequireNonNegative && value < Value{}) {
            std::ostringstream oss;
            oss << context << " requires non-negative valuation values; node " << nodeId << " has value " << value << ".";
            throw std::invalid_argument(oss.str());
        }
    }
 
  public:
    static void validateHierarchyConnectivity(const MorphologicalTree& tree, const RegularGridAdjacency2D& adjacency,
                                              const char* context = "HierarchySaliencyMapValidation::validateHierarchyConnectivity") {
        detail::requireCommittedRootedHierarchy(tree, context);
        const int rows = tree.numRows();
        const int columns = tree.numColumns();
        const int numPixels = tree.numPixels();
        if (rows <= 0 || columns <= 0 || numPixels <= 0 || adjacency.getNumRows() != rows || adjacency.getNumColumns() != columns ||
            static_cast<std::size_t>(numPixels) != static_cast<std::size_t>(rows) * static_cast<std::size_t>(columns)) {
            throw std::invalid_argument(std::string(context) + " requires one graph vertex per 2D proper part and matching adjacency dimensions.");
        }
 
        struct DomainDisjointSet {
            std::vector<NodeId> parent;
            std::vector<int> size;
 
            explicit DomainDisjointSet(int count) : parent(static_cast<std::size_t>(count)), size(static_cast<std::size_t>(count), 1) {
                for (NodeId id = 0; id < count; ++id) {
                    parent[static_cast<std::size_t>(id)] = id;
                }
            }
 
            NodeId find(NodeId id) {
                NodeId root = id;
                while (parent[static_cast<std::size_t>(root)] != root) {
                    root = parent[static_cast<std::size_t>(root)];
                }
                while (parent[static_cast<std::size_t>(id)] != id) {
                    const NodeId next = parent[static_cast<std::size_t>(id)];
                    parent[static_cast<std::size_t>(id)] = root;
                    id = next;
                }
                return root;
            }
 
            void unite(NodeId lhs, NodeId rhs) {
                lhs = find(lhs);
                rhs = find(rhs);
                if (lhs == rhs) {
                    return;
                }
                if (size[static_cast<std::size_t>(lhs)] < size[static_cast<std::size_t>(rhs)]) {
                    std::swap(lhs, rhs);
                }
                parent[static_cast<std::size_t>(rhs)] = lhs;
                size[static_cast<std::size_t>(lhs)] += size[static_cast<std::size_t>(rhs)];
            }
        };
 
        using DomainEdge = std::pair<NodeId, NodeId>;
        std::vector<std::vector<DomainEdge>> edgesByLca(static_cast<std::size_t>(tree.numInternalNodeSlots()));
        for (NodeId source = 0; source < numPixels; ++source) {
            const NodeId sourceSmallestNode = tree.smallestNode(source);
            if (!tree.isAlive(sourceSmallestNode)) {
                throw std::invalid_argument(std::string(context) + " found a pixel without a live smallest node.");
            }
            for (int targetValue : adjacency.getForwardNeighborIndices(source)) {
                const NodeId target = static_cast<NodeId>(targetValue);
                const NodeId targetSmallestNode = tree.smallestNode(target);
                if (!tree.isAlive(targetSmallestNode)) {
                    throw std::invalid_argument(std::string(context) + " found a neighbour pixel without a live smallest node.");
                }
                const NodeId lca = sourceSmallestNode == targetSmallestNode ? sourceSmallestNode : tree.lowestCommonAncestor(sourceSmallestNode, targetSmallestNode);
                if (lca == InvalidNode || !tree.isAlive(lca)) {
                    throw std::invalid_argument(std::string(context) + " could not assign an adjacency edge to a live hierarchy node.");
                }
                edgesByLca[static_cast<std::size_t>(lca)].emplace_back(source, target);
            }
        }
 
        DomainDisjointSet components(numPixels);
        std::vector<NodeId> supportRepresentative(static_cast<std::size_t>(tree.numInternalNodeSlots()), InvalidNode);
        for (NodeId nodeId : tree.postOrder()) {
            for (const DomainEdge& edge : edgesByLca[static_cast<std::size_t>(nodeId)]) {
                components.unite(edge.first, edge.second);
            }
 
            NodeId representative = InvalidNode;
            auto acceptRepresentative = [&](NodeId candidate) {
                if (representative == InvalidNode) {
                    representative = candidate;
                    return;
                }
                if (components.find(representative) != components.find(candidate)) {
                    std::ostringstream oss;
                    oss << context << " requires every hierarchy region to be connected in the projection graph; node " << nodeId
                        << " has a disconnected support.";
                    throw std::invalid_argument(oss.str());
                }
            };
 
            for (PixelId pixel : tree.properPart(nodeId)) {
                acceptRepresentative(pixel);
            }
            for (NodeId childId : tree.children(nodeId)) {
                const NodeId childRepresentative = supportRepresentative[static_cast<std::size_t>(childId)];
                if (childRepresentative == InvalidNode) {
                    throw std::invalid_argument(std::string(context) + " found a live child without proper-part support.");
                }
                acceptRepresentative(childRepresentative);
            }
            if (representative == InvalidNode) {
                throw std::invalid_argument(std::string(context) + " found a live node without proper-part support.");
            }
            supportRepresentative[static_cast<std::size_t>(nodeId)] = components.find(representative);
        }
    }
 
    template <class Value>
    static void validateHierarchyValuation(const MorphologicalTree& tree, std::span<const Value> valuation,
                                           HierarchyValuationPolicy policy = HierarchyValuationPolicy::AllowLevelCollapse,
                                           HierarchyValuationRangePolicy rangePolicy = HierarchyValuationRangePolicy::AllowAnyFinite,
                                           const char* context = "HierarchySaliencyMapValidation::validateHierarchyValuation") {
        detail::requireCommittedRootedHierarchy(tree, context);
        if (valuation.size() != static_cast<std::size_t>(tree.numInternalNodeSlots())) {
            std::ostringstream oss;
            oss << context << " requires one valuation value per dense internal node slot; expected " << tree.numInternalNodeSlots() << " values but got "
                << valuation.size() << ".";
            throw std::invalid_argument(oss.str());
        }
 
        for (NodeId nodeId : tree.aliveNodeIds()) {
            validateValuationValue(valuation[static_cast<std::size_t>(nodeId)], nodeId, context, rangePolicy);
        }
 
        for (NodeId parentId : tree.aliveNodeIds()) {
            const Value& parentValue = valuation[static_cast<std::size_t>(parentId)];
            for (NodeId childId : tree.children(parentId)) {
                const Value& childValue = valuation[static_cast<std::size_t>(childId)];
                const bool validOrder = policy == HierarchyValuationPolicy::RequireStrictHierarchy ? childValue < parentValue : !(parentValue < childValue);
                if (!validOrder) {
                    std::ostringstream oss;
                    oss << context << " requires "
                        << (policy == HierarchyValuationPolicy::RequireStrictHierarchy ? "valuation(parent) > valuation(child)"
                                                                                       : "valuation(parent) >= valuation(child)")
                        << "; parent node " << parentId << " has value " << parentValue << " and child node " << childId << " has value " << childValue << ".";
                    throw std::invalid_argument(oss.str());
                }
            }
        }
    }
 
    template <class Value>
    [[nodiscard]] static std::vector<int> rankHierarchyValuation(const MorphologicalTree& tree, std::span<const Value> valuation,
                                                                 HierarchyValuationPolicy policy = HierarchyValuationPolicy::AllowLevelCollapse) {
        validateHierarchyValuation(tree, valuation, policy, HierarchyValuationRangePolicy::AllowAnyFinite,
                                   "HierarchySaliencyMapValidation::rankHierarchyValuation");
 
        std::vector<Value> uniqueValues;
        uniqueValues.reserve(static_cast<std::size_t>(tree.numInternalNodeSlots()));
        for (NodeId nodeId : tree.aliveNodeIds()) {
            uniqueValues.push_back(valuation[static_cast<std::size_t>(nodeId)]);
        }
 
        std::sort(uniqueValues.begin(), uniqueValues.end());
        uniqueValues.erase(std::unique(uniqueValues.begin(), uniqueValues.end()), uniqueValues.end());
 
        std::vector<int> ranks(static_cast<std::size_t>(tree.numInternalNodeSlots()), 0);
        for (NodeId nodeId : tree.aliveNodeIds()) {
            const Value& value = valuation[static_cast<std::size_t>(nodeId)];
            const auto it = std::lower_bound(uniqueValues.begin(), uniqueValues.end(), value);
            ranks[static_cast<std::size_t>(nodeId)] = static_cast<int>(std::distance(uniqueValues.begin(), it));
        }
        return ranks;
    }
 
    template <class Value>
    [[nodiscard]] static std::vector<double>
    computeNormalizedScores(const MorphologicalTree& tree, std::span<const Value> valuation,
                            HierarchyValuationPolicy policy = HierarchyValuationPolicy::AllowLevelCollapse,
                            HierarchyValuationRangePolicy rangePolicy = HierarchyValuationRangePolicy::AllowAnyFinite) {
        using BareValue = std::remove_cv_t<Value>;
        static_assert(std::is_arithmetic_v<BareValue> && !std::is_same_v<BareValue, bool>,
                      "HierarchySaliencyMapValidation::computeNormalizedScores requires a numeric non-bool valuation type.");
        validateHierarchyValuation(tree, valuation, policy, rangePolicy, "HierarchySaliencyMapValidation::computeNormalizedScores");
 
        bool initialized = false;
        BareValue minValue{};
        BareValue maxValue{};
        for (NodeId nodeId : tree.aliveNodeIds()) {
            const BareValue value = valuation[static_cast<std::size_t>(nodeId)];
            if (!initialized) {
                minValue = value;
                maxValue = value;
                initialized = true;
            } else {
                minValue = std::min(minValue, value);
                maxValue = std::max(maxValue, value);
            }
        }
 
        std::vector<double> scores(static_cast<std::size_t>(tree.numInternalNodeSlots()), 0.0);
        if (maxValue == minValue) {
            return scores;
        }
 
        for (NodeId nodeId : tree.aliveNodeIds()) {
            const BareValue value = valuation[static_cast<std::size_t>(nodeId)];
            long double normalized = 0.0L;
 
            if constexpr (std::is_integral_v<BareValue>) {
                using UnsignedValue = std::make_unsigned_t<BareValue>;
                const UnsignedValue offset = static_cast<UnsignedValue>(value) - static_cast<UnsignedValue>(minValue);
                const UnsignedValue range = static_cast<UnsignedValue>(maxValue) - static_cast<UnsignedValue>(minValue);
                normalized = static_cast<long double>(offset) / static_cast<long double>(range);
            } else {
                const long double low = static_cast<long double>(minValue);
                const long double high = static_cast<long double>(maxValue);
                const long double current = static_cast<long double>(value);
 
                if (low < 0.0L && high > 0.0L) {
                    const long double scale = std::max(-low, high);
                    const long double scaledLow = low / scale;
                    const long double scaledHigh = high / scale;
                    normalized = (current / scale - scaledLow) / (scaledHigh - scaledLow);
                } else {
                    normalized = (current - low) / (high - low);
                }
            }
 
            scores[static_cast<std::size_t>(nodeId)] = static_cast<double>(std::clamp(normalized, 0.0L, 1.0L));
        }
        return scores;
    }
 
    template <AltitudeValue T> [[nodiscard]] static std::vector<double> computeNormalizedScores(const ValuedMorphologicalTree<T>& tree) {
        const MorphologicalTree& topology = tree.topology();
        detail::requireCommittedRootedHierarchy(topology, "HierarchySaliencyMapValidation::computeNormalizedScores");
        const NodeAltitudeOrder nodeAltitudeOrder = topology.nodeAltitudeOrder();
        if (nodeAltitudeOrder == NodeAltitudeOrder::Unconstrained) {
            throw std::invalid_argument("HierarchySaliencyMapValidation::computeNormalizedScores requires a globally monotone altitude order.");
        }
        // Use at least double precision while retaining long double when it is
        // the input type. An unconditional cast to double can collapse distinct
        // long-double hierarchy levels and change the induced hierarchy.
        using OrientedAltitude = std::common_type_t<T, double>;
        std::vector<OrientedAltitude> orientedAltitude(static_cast<std::size_t>(topology.numInternalNodeSlots()), OrientedAltitude{});
        for (NodeId nodeId : topology.aliveNodeIds()) {
            const OrientedAltitude altitude = static_cast<OrientedAltitude>(tree.nodeAltitude(nodeId));
            orientedAltitude[static_cast<std::size_t>(nodeId)] = nodeAltitudeOrder == NodeAltitudeOrder::Increasing ? -altitude : altitude;
        }
        return computeNormalizedScores(topology, std::span<const OrientedAltitude>(orientedAltitude), HierarchyValuationPolicy::AllowLevelCollapse);
    }
};
 
} // namespace mmcfilters