HierarchySaliencyMap.hpp

#pragma once
 
#include "HierarchySaliencyMapValidation.hpp"
#include "../../utils/RegularGridAdjacency2D.hpp"
#include "../../utils/Altitude.hpp"
#include "../../utils/Common.hpp"
#include "../MorphologicalTree.hpp"
#include "../detail/MorphologicalTreeConstructionContextQueries.hpp"
#include "../ValuedMorphologicalTree.hpp"
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <functional>
#include <iterator>
#include <sstream>
#include <span>
#include <stdexcept>
#include <string>
#include <type_traits>
#include <utility>
#include <vector>
 
namespace mmcfilters {
 
template <class Value> struct EdgeSaliencyMap {
    int numRows = 0;
    int numColumns = 0;
    double adjacencyRadius = 0.0;
    std::vector<NodeId> sources;
    std::vector<NodeId> targets;
    std::vector<Value> values;
 
    [[nodiscard]] std::size_t size() const noexcept { return values.size(); }
 
    [[nodiscard]] bool empty() const noexcept { return values.empty(); }
};
 
enum class HierarchyLevelConvention {
    EdgeSaliencyValue,
    PartitionAppearanceLevel,
};
 
class HierarchySaliencyMap {
  private:
    static void validateTreeAndAdjacency(const MorphologicalTree& tree, const RegularGridAdjacency2D& adjacency, const char* context) {
        detail::requireCommittedRootedHierarchy(tree, context);
        if (tree.numRows() <= 0 || tree.numColumns() <= 0 || tree.numPixels() <= 0) {
            throw std::invalid_argument(std::string(context) + " requires a non-empty image/pixel domain.");
        }
        if (adjacency.getNumRows() != tree.numRows() || adjacency.getNumColumns() != tree.numColumns()) {
            std::ostringstream oss;
            oss << context << " adjacency domain must match the tree image domain; got " << adjacency.getNumRows() << "x" << adjacency.getNumColumns()
                << " for a tree domain of " << tree.numRows() << "x" << tree.numColumns() << ".";
            throw std::invalid_argument(oss.str());
        }
    }
 
    static bool sameAdjacencyGraph(const RegularGridAdjacency2D& lhs, const RegularGridAdjacency2D& rhs) noexcept {
        if (lhs.getNumRows() != rhs.getNumRows() || lhs.getNumColumns() != rhs.getNumColumns() || lhs.getSize() != rhs.getSize()) {
            return false;
        }
        for (int index = 0; index < lhs.getSize(); ++index) {
            if (lhs.getOffsetRow(index) != rhs.getOffsetRow(index) || lhs.getOffsetColumn(index) != rhs.getOffsetColumn(index)) {
                return false;
            }
        }
        return true;
    }
 
    static RegularGridAdjacency2D requireStoredAdjacency(const MorphologicalTree& tree, const char* context) {
        if (const RegularGridAdjacency2D* adjacency = ::mmcfilters::detail::constructionAdjacency(tree)) {
            return *adjacency;
        }
        if (const ComplementaryAdjacencies* adjacencies = ::mmcfilters::detail::complementaryAdjacencies(tree)) {
            if (sameAdjacencyGraph(adjacencies->minAdjacency, adjacencies->maxAdjacency)) {
                return adjacencies->minAdjacency;
            }
            throw std::invalid_argument(
                std::string(context) +
                " cannot infer one image graph from distinct minimum/maximum adjacencies; pass the intended adjacency explicitly.");
        }
        if (const TopographicConvention* convention = tree.topographicConvention();
            convention != nullptr && std::holds_alternative<SelfDualSpanImmersion>(convention->immersion) && tree.hasGridDomain2D()) {
            const GridDomain2D& domain = *tree.gridDomain2D();
            return RegularGridAdjacency2D(domain.rows, domain.columns, 1.0);
        }
        throw std::invalid_argument(std::string(context) + " requires an attached adjacency relation; pass an explicit adjacency relation instead.");
    }
 
  public:
    [[nodiscard]] static RegularGridAdjacency2D requireProjectionAdjacency(
        const MorphologicalTree& tree, const char* context = "HierarchySaliencyMap::requireProjectionAdjacency") {
        return requireStoredAdjacency(tree, context);
    }
 
    [[nodiscard]] static std::vector<int> computeTopologicalLevels(const MorphologicalTree& tree) {
        detail::requireCommittedRootedHierarchy(tree, "HierarchySaliencyMap::computeTopologicalLevels");
 
        std::vector<int> levels(static_cast<std::size_t>(tree.numInternalNodeSlots()), 0);
        for (NodeId nodeId : tree.postOrder()) {
            int level = 0;
            for (NodeId childId : tree.children(nodeId)) {
                level = std::max(level, levels[static_cast<std::size_t>(childId)] + 1);
            }
            levels[static_cast<std::size_t>(nodeId)] = level;
        }
        return levels;
    }
 
    [[nodiscard]] static std::vector<int> computePartitionAppearanceLevels(const MorphologicalTree& tree) {
        std::vector<int> levels = computeTopologicalLevels(tree);
        for (NodeId nodeId : tree.aliveNodeIds()) {
            ++levels[static_cast<std::size_t>(nodeId)];
        }
        return levels;
    }
 
    template <class NodeValue>
    [[nodiscard]] static auto computeEdgeMap(const MorphologicalTree& tree, const RegularGridAdjacency2D& adjacency, NodeValue&& nodeValue)
        -> EdgeSaliencyMap<std::decay_t<std::invoke_result_t<NodeValue, NodeId>>> {
        using Value = std::decay_t<std::invoke_result_t<NodeValue, NodeId>>;
        static_assert(!std::is_void_v<Value>, "HierarchySaliencyMap nodeValue must return a value.");
 
        constexpr const char* context = "HierarchySaliencyMap::computeEdgeMap";
        validateTreeAndAdjacency(tree, adjacency, context);
 
        EdgeSaliencyMap<Value> saliency;
        saliency.numRows = tree.numRows();
        saliency.numColumns = tree.numColumns();
        saliency.adjacencyRadius = adjacency.getRadius();
 
        auto&& scorer = nodeValue;
        const int numPixels = tree.numPixels();
        for (NodeId source = 0; source < numPixels; ++source) {
            const NodeId sourceSmallestNode = tree.smallestNode(source);
            if (!tree.isAlive(sourceSmallestNode)) {
                throw std::runtime_error("HierarchySaliencyMap::computeEdgeMap found a pixel without a live smallest node.");
            }
 
            for (int target : adjacency.getForwardNeighborIndices(source)) {
                const NodeId targetSmallestNode = tree.smallestNode(target);
                if (!tree.isAlive(targetSmallestNode)) {
                    throw std::runtime_error("HierarchySaliencyMap::computeEdgeMap found a neighbour pixel without a live smallest node.");
                }
 
                const NodeId lca = tree.lowestCommonAncestor(sourceSmallestNode, targetSmallestNode);
                if (lca == InvalidNode || !tree.isAlive(lca)) {
                    throw std::runtime_error("HierarchySaliencyMap::computeEdgeMap could not find a live LCA for an adjacency edge.");
                }
 
                saliency.sources.push_back(source);
                saliency.targets.push_back(target);
                saliency.values.push_back(std::invoke(scorer, lca));
            }
        }
 
        return saliency;
    }
 
    template <class Value>
    [[nodiscard]] static EdgeSaliencyMap<Value>
    computeSaliencyEdgeMap(const MorphologicalTree& tree, const RegularGridAdjacency2D& adjacency, std::span<const Value> valuation,
                           HierarchyValuationPolicy policy = HierarchyValuationPolicy::AllowLevelCollapse,
                           HierarchyLevelConvention levelConvention = HierarchyLevelConvention::EdgeSaliencyValue,
                           HierarchyConnectivityPolicy connectivityPolicy = HierarchyConnectivityPolicy::ValidateConnected) {
        constexpr const char* context = "HierarchySaliencyMap::computeSaliencyEdgeMap";
        HierarchySaliencyMapValidation::validateHierarchyValuation(tree, valuation, policy, HierarchyValuationRangePolicy::RequireNonNegative, context);
        validateTreeAndAdjacency(tree, adjacency, context);
        if (connectivityPolicy == HierarchyConnectivityPolicy::ValidateConnected) {
            HierarchySaliencyMapValidation::validateHierarchyConnectivity(tree, adjacency, context);
        }
        if (levelConvention == HierarchyLevelConvention::PartitionAppearanceLevel) {
            for (NodeId nodeId : tree.aliveNodeIds()) {
                const Value& level = valuation[static_cast<std::size_t>(nodeId)];
                if (!(Value{} < level)) {
                    throw std::invalid_argument(std::string(context) + " requires positive partition-appearance levels.");
                }
                if constexpr (std::is_floating_point_v<Value>) {
                    if (std::floor(level) != level) {
                        throw std::invalid_argument(std::string(context) + " requires integer-valued partition-appearance levels.");
                    }
                }
            }
        }
 
        EdgeSaliencyMap<Value> saliency;
        saliency.numRows = tree.numRows();
        saliency.numColumns = tree.numColumns();
        saliency.adjacencyRadius = adjacency.getRadius();
 
        const int numPixels = tree.numPixels();
        for (NodeId source = 0; source < numPixels; ++source) {
            const NodeId sourceSmallestNode = tree.smallestNode(source);
            if (!tree.isAlive(sourceSmallestNode)) {
                throw std::runtime_error(std::string(context) + " found a pixel without a live smallest node.");
            }
 
            for (int target : adjacency.getForwardNeighborIndices(source)) {
                const NodeId targetSmallestNode = tree.smallestNode(target);
                if (!tree.isAlive(targetSmallestNode)) {
                    throw std::runtime_error(std::string(context) + " found a neighbour pixel without a live smallest node.");
                }
 
                Value edgeValue{};
                if (sourceSmallestNode != targetSmallestNode) {
                    const NodeId lca = tree.lowestCommonAncestor(sourceSmallestNode, targetSmallestNode);
                    if (lca == InvalidNode || !tree.isAlive(lca)) {
                        throw std::runtime_error(std::string(context) + " could not find a live LCA for an adjacency edge.");
                    }
                    edgeValue = valuation[static_cast<std::size_t>(lca)];
                    if (levelConvention == HierarchyLevelConvention::PartitionAppearanceLevel) {
                        edgeValue -= Value{1};
                    }
                }
 
                saliency.sources.push_back(source);
                saliency.targets.push_back(target);
                saliency.values.push_back(edgeValue);
            }
        }
 
        return saliency;
    }
 
    template <class Value>
    [[nodiscard]] static EdgeSaliencyMap<Value>
    computeSaliencyEdgeMap(const MorphologicalTree& tree, std::span<const Value> valuation,
                           HierarchyValuationPolicy policy = HierarchyValuationPolicy::AllowLevelCollapse,
                           HierarchyLevelConvention levelConvention = HierarchyLevelConvention::EdgeSaliencyValue,
                           HierarchyConnectivityPolicy connectivityPolicy = HierarchyConnectivityPolicy::ValidateConnected) {
        return computeSaliencyEdgeMap(tree, requireStoredAdjacency(tree, "HierarchySaliencyMap::computeSaliencyEdgeMap"), valuation, policy, levelConvention,
                                      connectivityPolicy);
    }
 
    template <class Value>
    [[nodiscard]] static EdgeSaliencyMap<int>
    computeCanonicalRankedSaliencyEdgeMap(const MorphologicalTree& tree, const RegularGridAdjacency2D& adjacency, std::span<const Value> valuation,
                                          HierarchyValuationPolicy policy = HierarchyValuationPolicy::AllowLevelCollapse,
                                          HierarchyConnectivityPolicy connectivityPolicy = HierarchyConnectivityPolicy::ValidateConnected) {
        constexpr const char* context = "HierarchySaliencyMap::computeCanonicalRankedSaliencyEdgeMap";
        HierarchySaliencyMapValidation::validateHierarchyValuation(tree, valuation, policy, HierarchyValuationRangePolicy::AllowAnyFinite, context);
        validateTreeAndAdjacency(tree, adjacency, context);
        if (connectivityPolicy == HierarchyConnectivityPolicy::ValidateConnected) {
            HierarchySaliencyMapValidation::validateHierarchyConnectivity(tree, adjacency, context);
        }
 
        EdgeSaliencyMap<int> saliency;
        saliency.numRows = tree.numRows();
        saliency.numColumns = tree.numColumns();
        saliency.adjacencyRadius = adjacency.getRadius();
 
        struct ProjectedLevel {
            NodeId source;
            NodeId target;
            NodeId lca;
        };
        std::vector<ProjectedLevel> projected;
        std::vector<Value> effectiveValues;
        bool hasBaseEdge = false;
        const int numPixels = tree.numPixels();
        for (NodeId source = 0; source < numPixels; ++source) {
            const NodeId sourceSmallestNode = tree.smallestNode(source);
            for (int targetValue : adjacency.getForwardNeighborIndices(source)) {
                const NodeId target = static_cast<NodeId>(targetValue);
                const NodeId targetSmallestNode = tree.smallestNode(target);
                NodeId lca = InvalidNode;
                if (sourceSmallestNode == targetSmallestNode) {
                    hasBaseEdge = true;
                } else {
                    lca = tree.lowestCommonAncestor(sourceSmallestNode, targetSmallestNode);
                    if (lca == InvalidNode || !tree.isAlive(lca)) {
                        throw std::runtime_error(std::string(context) + " could not find a live LCA for an adjacency edge.");
                    }
                    effectiveValues.push_back(valuation[static_cast<std::size_t>(lca)]);
                }
                projected.push_back(ProjectedLevel{source, target, lca});
            }
        }
 
        std::sort(effectiveValues.begin(), effectiveValues.end());
        effectiveValues.erase(std::unique(effectiveValues.begin(), effectiveValues.end()), effectiveValues.end());
        saliency.sources.reserve(projected.size());
        saliency.targets.reserve(projected.size());
        saliency.values.reserve(projected.size());
        const int transitionOffset = hasBaseEdge ? 1 : 0;
        for (const ProjectedLevel& edge : projected) {
            int rank = 0;
            if (edge.lca != InvalidNode) {
                const Value& value = valuation[static_cast<std::size_t>(edge.lca)];
                rank = transitionOffset +
                       static_cast<int>(std::distance(effectiveValues.begin(), std::lower_bound(effectiveValues.begin(), effectiveValues.end(), value)));
            }
            saliency.sources.push_back(edge.source);
            saliency.targets.push_back(edge.target);
            saliency.values.push_back(rank);
        }
        return saliency;
    }
 
    template <class Value>
    [[nodiscard]] static EdgeSaliencyMap<int>
    computeCanonicalRankedSaliencyEdgeMap(const MorphologicalTree& tree, std::span<const Value> valuation,
                                          HierarchyValuationPolicy policy = HierarchyValuationPolicy::AllowLevelCollapse,
                                          HierarchyConnectivityPolicy connectivityPolicy = HierarchyConnectivityPolicy::ValidateConnected) {
        return computeCanonicalRankedSaliencyEdgeMap(tree, requireStoredAdjacency(tree, "HierarchySaliencyMap::computeCanonicalRankedSaliencyEdgeMap"),
                                                     valuation, policy, connectivityPolicy);
    }
 
    template <class Value> [[nodiscard]] static EdgeSaliencyMap<int> rankEdgeSaliencyMap(const EdgeSaliencyMap<Value>& edgeMap) {
        EdgeSaliencyMap<int> ranked;
        ranked.numRows = edgeMap.numRows;
        ranked.numColumns = edgeMap.numColumns;
        ranked.adjacencyRadius = edgeMap.adjacencyRadius;
        ranked.sources = edgeMap.sources;
        ranked.targets = edgeMap.targets;
        std::vector<Value> uniqueValues = edgeMap.values;
        std::sort(uniqueValues.begin(), uniqueValues.end());
        uniqueValues.erase(std::unique(uniqueValues.begin(), uniqueValues.end()), uniqueValues.end());
        ranked.values.reserve(edgeMap.values.size());
        for (const Value& value : edgeMap.values) {
            ranked.values.push_back(static_cast<int>(std::distance(uniqueValues.begin(), std::lower_bound(uniqueValues.begin(), uniqueValues.end(), value))));
        }
        return ranked;
    }
 
    template <class NodeValue>
    [[nodiscard]] static auto computeEdgeMap(const MorphologicalTree& tree, NodeValue&& nodeValue)
        -> EdgeSaliencyMap<std::decay_t<std::invoke_result_t<NodeValue, NodeId>>> {
        return computeEdgeMap(tree, requireStoredAdjacency(tree, "HierarchySaliencyMap::computeEdgeMap"), std::forward<NodeValue>(nodeValue));
    }
 
    [[nodiscard]] static EdgeSaliencyMap<int> computeTopologicalLevelEdgeMap(const MorphologicalTree& tree, const RegularGridAdjacency2D& adjacency) {
        std::vector<int> levels = computeTopologicalLevels(tree);
        return computeCanonicalRankedSaliencyEdgeMap(tree, adjacency, std::span<const int>(levels), HierarchyValuationPolicy::RequireStrictHierarchy);
    }
 
    [[nodiscard]] static EdgeSaliencyMap<int> computeTopologicalLevelEdgeMap(const MorphologicalTree& tree) {
        return computeTopologicalLevelEdgeMap(tree, requireStoredAdjacency(tree, "HierarchySaliencyMap::computeTopologicalLevelEdgeMap"));
    }
 
    template <AltitudeValue T>
    [[nodiscard]] static EdgeSaliencyMap<double> computeNormalizedAltitudeEdgeMap(const ValuedMorphologicalTree<T>& tree,
                                                                                  const RegularGridAdjacency2D& adjacency) {
        std::vector<double> scores = HierarchySaliencyMapValidation::computeNormalizedScores(tree);
        return computeSaliencyEdgeMap(tree.topology(), adjacency, std::span<const double>(scores), HierarchyValuationPolicy::AllowLevelCollapse);
    }
 
    template <AltitudeValue T> [[nodiscard]] static EdgeSaliencyMap<double> computeNormalizedAltitudeEdgeMap(const ValuedMorphologicalTree<T>& tree) {
        return computeNormalizedAltitudeEdgeMap(tree, requireStoredAdjacency(tree.topology(), "HierarchySaliencyMap::computeNormalizedAltitudeEdgeMap"));
    }
};
 
class ComponentTreePartitionHierarchyAdapter {
  public:
    static void validate(const MorphologicalTree& tree, const RegularGridAdjacency2D& adjacency) {
        HierarchySaliencyMapValidation::validateHierarchyConnectivity(tree, adjacency, "ComponentTreePartitionHierarchyAdapter::validate");
    }
 
    [[nodiscard]] static std::vector<int> computePartitionAppearanceLevels(const MorphologicalTree& tree) {
        return HierarchySaliencyMap::computePartitionAppearanceLevels(tree);
    }
 
    template <class Value>
    [[nodiscard]] static EdgeSaliencyMap<Value> computeSaliencyEdgeMap(const MorphologicalTree& tree, const RegularGridAdjacency2D& adjacency,
                                                                       std::span<const Value> partitionAppearanceLevels,
                                                                       HierarchyValuationPolicy policy = HierarchyValuationPolicy::AllowLevelCollapse) {
        return HierarchySaliencyMap::computeSaliencyEdgeMap(tree, adjacency, partitionAppearanceLevels, policy,
                                                            HierarchyLevelConvention::PartitionAppearanceLevel, HierarchyConnectivityPolicy::ValidateConnected);
    }
};
 
} // namespace mmcfilters