HierarchySaliencyMapProjection.hpp

#pragma once
 
#include "HierarchySaliencyMap.hpp"
#include "../../utils/Image.hpp"
#include "../MorphologicalTree.hpp"
 
#include <cstddef>
#include <cstdint>
#include <numeric>
#include <span>
#include <stdexcept>
#include <string>
#include <vector>
 
namespace mmcfilters {
 
struct EdgeContourMap {
    int numRows = 0;
    int numColumns = 0;
    double adjacencyRadius = 0.0;
    std::vector<NodeId> sources;
    std::vector<NodeId> targets;
 
    [[nodiscard]] std::size_t size() const noexcept { return sources.size(); }
 
    [[nodiscard]] bool empty() const noexcept { return sources.empty(); }
};
 
enum class EdgeToPixelReducer {
    Max,
    Mean
};
 
struct NodeContourEdgeMap {
    int numRows = 0;
    int numColumns = 0;
    double adjacencyRadius = 0.0;
    std::vector<NodeId> sources;
    std::vector<NodeId> targets;
    std::vector<NodeId> nodes;
 
    [[nodiscard]] std::size_t size() const noexcept { return nodes.size(); }
 
    [[nodiscard]] bool empty() const noexcept { return nodes.empty(); }
};
 
struct IncrementalNodeContourMap {
    int numRows = 0;
    int numColumns = 0;
    int numNodeSlots = 0;
    double adjacencyRadius = 0.0;
    std::vector<std::size_t> offsets;
    std::vector<NodeId> sources;
    std::vector<NodeId> targets;
 
    [[nodiscard]] std::size_t size() const noexcept { return sources.size(); }
 
    [[nodiscard]] bool empty() const noexcept { return sources.empty(); }
 
    [[nodiscard]] std::size_t nodeBegin(NodeId nodeId) const { return offsets.at(static_cast<std::size_t>(nodeId)); }
 
    [[nodiscard]] std::size_t nodeEnd(NodeId nodeId) const { return offsets.at(static_cast<std::size_t>(nodeId) + 1); }
 
    [[nodiscard]] std::size_t nodeSize(NodeId nodeId) const { return nodeEnd(nodeId) - nodeBegin(nodeId); }
};
 
class HierarchySaliencyMapProjection {
  private:
    template <class Value> static void validateEdgeMap(const EdgeSaliencyMap<Value>& edgeMap, const char* context) {
        if (edgeMap.sources.size() != edgeMap.targets.size() || edgeMap.sources.size() != edgeMap.values.size()) {
            throw std::invalid_argument(std::string(context) + " requires sources, targets, and values to have the same length.");
        }
        if (edgeMap.numRows <= 0 || edgeMap.numColumns <= 0) {
            throw std::invalid_argument(std::string(context) + " requires a non-empty image domain.");
        }
    }
 
    static void validateIncrementalContours(const IncrementalNodeContourMap& contours, const char* context) {
        if (contours.numRows <= 0 || contours.numColumns <= 0) {
            throw std::invalid_argument(std::string(context) + " requires a non-empty image domain.");
        }
        if (contours.numNodeSlots < 0) {
            throw std::invalid_argument(std::string(context) + " requires a non-negative node-slot count.");
        }
        if (contours.sources.size() != contours.targets.size()) {
            throw std::invalid_argument(std::string(context) + " requires sources and targets to have the same length.");
        }
        if (contours.offsets.size() != static_cast<std::size_t>(contours.numNodeSlots) + 1) {
            throw std::invalid_argument(std::string(context) + " requires offsets to have numNodeSlots + 1 entries.");
        }
        if (contours.offsets.empty() || contours.offsets.front() != 0 || contours.offsets.back() != contours.sources.size()) {
            throw std::invalid_argument(std::string(context) + " requires offsets to span exactly the contour edge arrays.");
        }
        for (std::size_t i = 1; i < contours.offsets.size(); ++i) {
            if (contours.offsets[i] < contours.offsets[i - 1] || contours.offsets[i] > contours.sources.size()) {
                throw std::invalid_argument(std::string(context) + " requires monotonically increasing offsets inside the edge-array bounds.");
            }
        }
    }
 
    template <class Value> static void validateEdgeEndpoints(const EdgeSaliencyMap<Value>& edgeMap, const char* context) {
        const auto pixelCount = static_cast<std::size_t>(edgeMap.numRows) * static_cast<std::size_t>(edgeMap.numColumns);
        for (std::size_t i = 0; i < edgeMap.size(); ++i) {
            const NodeId source = edgeMap.sources[i];
            const NodeId target = edgeMap.targets[i];
            if (source < 0 || static_cast<std::size_t>(source) >= pixelCount || target < 0 || static_cast<std::size_t>(target) >= pixelCount) {
                throw std::invalid_argument(std::string(context) + " requires sources and targets inside the image domain.");
            }
        }
    }
 
    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()) {
            throw std::invalid_argument(std::string(context) + " adjacency domain must match the tree image domain.");
        }
    }
 
    static IncrementalNodeContourMap groupNodeContourEdges(NodeContourEdgeMap&& flat, int numNodeSlots) {
        IncrementalNodeContourMap contours;
        contours.numRows = flat.numRows;
        contours.numColumns = flat.numColumns;
        contours.numNodeSlots = numNodeSlots;
        contours.adjacencyRadius = flat.adjacencyRadius;
        contours.offsets.assign(static_cast<std::size_t>(contours.numNodeSlots) + 1, 0);
 
        for (NodeId nodeId : flat.nodes) {
            ++contours.offsets[static_cast<std::size_t>(nodeId) + 1];
        }
        std::partial_sum(contours.offsets.begin(), contours.offsets.end(), contours.offsets.begin());
 
        contours.sources.resize(flat.size());
        contours.targets.resize(flat.size());
        std::vector<std::size_t> cursor = contours.offsets;
        for (std::size_t i = 0; i < flat.size(); ++i) {
            const NodeId nodeId = flat.nodes[i];
            const std::size_t dst = cursor[static_cast<std::size_t>(nodeId)]++;
            contours.sources[dst] = flat.sources[i];
            contours.targets[dst] = flat.targets[i];
        }
 
        return contours;
    }
 
  public:
    template <class Value>
    [[nodiscard]] static ImagePtr<double> edgeMapToPixelImage(const EdgeSaliencyMap<Value>& edgeMap, EdgeToPixelReducer reducer = EdgeToPixelReducer::Max) {
        constexpr const char* context = "HierarchySaliencyMapProjection::edgeMapToPixelImage";
        validateEdgeMap(edgeMap, context);
        validateEdgeEndpoints(edgeMap, context);
 
        ImagePtr<double> image = Image<double>::create(edgeMap.numRows, edgeMap.numColumns, 0.0);
        double* pixels = image->rawData();
        const auto pixelCount = static_cast<std::size_t>(edgeMap.numRows) * static_cast<std::size_t>(edgeMap.numColumns);
 
        switch (reducer) {
        case EdgeToPixelReducer::Max: {
            std::vector<uint8_t> initialized(pixelCount, false);
            for (std::size_t i = 0; i < edgeMap.size(); ++i) {
                const double value = static_cast<double>(edgeMap.values[i]);
                const auto source = static_cast<std::size_t>(edgeMap.sources[i]);
                const auto target = static_cast<std::size_t>(edgeMap.targets[i]);
                if (!initialized[source] || pixels[source] < value) {
                    pixels[source] = value;
                    initialized[source] = true;
                }
                if (!initialized[target] || pixels[target] < value) {
                    pixels[target] = value;
                    initialized[target] = true;
                }
            }
            return image;
        }
        case EdgeToPixelReducer::Mean: {
            std::vector<std::size_t> count(pixelCount, 0);
            for (std::size_t i = 0; i < edgeMap.size(); ++i) {
                const double value = static_cast<double>(edgeMap.values[i]);
                const auto source = static_cast<std::size_t>(edgeMap.sources[i]);
                const auto target = static_cast<std::size_t>(edgeMap.targets[i]);
                pixels[source] += value;
                pixels[target] += value;
                ++count[source];
                ++count[target];
            }
            for (std::size_t i = 0; i < pixelCount; ++i) {
                if (count[i] != 0) {
                    pixels[i] /= static_cast<double>(count[i]);
                }
            }
            return image;
        }
        default:
            throw std::invalid_argument(std::string(context) + " received an unknown edge-to-pixel reducer.");
        }
    }
 
    template <class Value, class Threshold> [[nodiscard]] static EdgeContourMap thresholdCut(const EdgeSaliencyMap<Value>& edgeMap, Threshold threshold) {
        constexpr const char* context = "HierarchySaliencyMapProjection::thresholdCut";
        validateEdgeMap(edgeMap, context);
 
        EdgeContourMap contours;
        contours.numRows = edgeMap.numRows;
        contours.numColumns = edgeMap.numColumns;
        contours.adjacencyRadius = edgeMap.adjacencyRadius;
        contours.sources.reserve(edgeMap.size());
        contours.targets.reserve(edgeMap.size());
 
        for (std::size_t i = 0; i < edgeMap.size(); ++i) {
            if (edgeMap.values[i] >= threshold) {
                contours.sources.push_back(edgeMap.sources[i]);
                contours.targets.push_back(edgeMap.targets[i]);
            }
        }
        return contours;
    }
 
    [[nodiscard]] static NodeContourEdgeMap nodeContourEdges(const MorphologicalTree& tree, const RegularGridAdjacency2D& adjacency) {
        constexpr const char* context = "HierarchySaliencyMapProjection::nodeContourEdges";
        validateTreeAndAdjacency(tree, adjacency, context);
 
        NodeContourEdgeMap contours;
        contours.numRows = tree.numRows();
        contours.numColumns = tree.numColumns();
        contours.adjacencyRadius = adjacency.getRadius();
 
        const int numPixels = tree.numPixels();
        for (NodeId source = 0; source < numPixels; ++source) {
            const NodeId sourceSmallestNode = detail::TreeKernelAccess::smallestNodeMap(tree, source);
 
            for (int target : detail::CommittedGridAccess::forwardNeighbors(adjacency, source)) {
                const NodeId targetSmallestNode = detail::TreeKernelAccess::smallestNodeMap(tree, target);
                if (sourceSmallestNode == targetSmallestNode) {
                    continue;
                }
 
                const NodeId lca = detail::TreeKernelAccess::lowestCommonAncestor(tree, sourceSmallestNode, targetSmallestNode);
 
                contours.sources.push_back(source);
                contours.targets.push_back(target);
                contours.nodes.push_back(lca);
            }
        }
        return contours;
    }
 
    [[nodiscard]] static NodeContourEdgeMap nodeContourEdges(const MorphologicalTree& tree) {
        const RegularGridAdjacency2D* adjacency = tree.constructionAdjacency();
        if (adjacency == nullptr) {
            throw std::invalid_argument(
                "HierarchySaliencyMapProjection::nodeContourEdges requires an attached adjacency relation; pass an explicit adjacency relation instead.");
        }
        return nodeContourEdges(tree, *adjacency);
    }
 
    template <AltitudeValue T>
    [[nodiscard]] static NodeContourEdgeMap nodeContourEdges(const ValuedMorphologicalTree<T>& tree, const RegularGridAdjacency2D& adjacency) {
        return nodeContourEdges(tree.topology(), adjacency);
    }
 
    template <AltitudeValue T> [[nodiscard]] static NodeContourEdgeMap nodeContourEdges(const ValuedMorphologicalTree<T>& tree) {
        return nodeContourEdges(tree.topology());
    }
 
    [[nodiscard]] static IncrementalNodeContourMap computeIncrementalNodeContours(const MorphologicalTree& tree, const RegularGridAdjacency2D& adjacency) {
        return groupNodeContourEdges(nodeContourEdges(tree, adjacency), tree.numInternalNodeSlots());
    }
 
    [[nodiscard]] static IncrementalNodeContourMap computeIncrementalNodeContours(const MorphologicalTree& tree) {
        return groupNodeContourEdges(nodeContourEdges(tree), tree.numInternalNodeSlots());
    }
 
    template <AltitudeValue T>
    [[nodiscard]] static IncrementalNodeContourMap computeIncrementalNodeContours(const ValuedMorphologicalTree<T>& tree,
                                                                                  const RegularGridAdjacency2D& adjacency) {
        return computeIncrementalNodeContours(tree.topology(), adjacency);
    }
 
    template <AltitudeValue T> [[nodiscard]] static IncrementalNodeContourMap computeIncrementalNodeContours(const ValuedMorphologicalTree<T>& tree) {
        return computeIncrementalNodeContours(tree.topology());
    }
 
    template <class Value>
    [[nodiscard]] static EdgeSaliencyMap<Value> projectNodeValuation(const IncrementalNodeContourMap& contours, std::span<const Value> nodeValuation) {
        constexpr const char* context = "HierarchySaliencyMapProjection::projectNodeValuation";
        validateIncrementalContours(contours, context);
        if (nodeValuation.size() != static_cast<std::size_t>(contours.numNodeSlots)) {
            throw std::invalid_argument(std::string(context) + " requires one valuation value per dense node slot.");
        }
 
        EdgeSaliencyMap<Value> edgeMap;
        edgeMap.numRows = contours.numRows;
        edgeMap.numColumns = contours.numColumns;
        edgeMap.adjacencyRadius = contours.adjacencyRadius;
        edgeMap.sources.reserve(contours.size());
        edgeMap.targets.reserve(contours.size());
        edgeMap.values.reserve(contours.size());
 
        for (NodeId nodeId = 0; nodeId < contours.numNodeSlots; ++nodeId) {
            const Value value = nodeValuation[static_cast<std::size_t>(nodeId)];
            for (std::size_t i = contours.offsets[static_cast<std::size_t>(nodeId)]; i < contours.offsets[static_cast<std::size_t>(nodeId) + 1]; ++i) {
                edgeMap.sources.push_back(contours.sources[i]);
                edgeMap.targets.push_back(contours.targets[i]);
                edgeMap.values.push_back(value);
            }
        }
        return edgeMap;
    }
 
    template <class Value, class Threshold>
    [[nodiscard]] static EdgeContourMap thresholdByNodeValuation(const IncrementalNodeContourMap& contours, std::span<const Value> nodeValuation,
                                                                 Threshold threshold) {
        constexpr const char* context = "HierarchySaliencyMapProjection::thresholdByNodeValuation";
        validateIncrementalContours(contours, context);
        if (nodeValuation.size() != static_cast<std::size_t>(contours.numNodeSlots)) {
            throw std::invalid_argument(std::string(context) + " requires one valuation value per dense node slot.");
        }
 
        EdgeContourMap cut;
        cut.numRows = contours.numRows;
        cut.numColumns = contours.numColumns;
        cut.adjacencyRadius = contours.adjacencyRadius;
        cut.sources.reserve(contours.size());
        cut.targets.reserve(contours.size());
 
        for (NodeId nodeId = 0; nodeId < contours.numNodeSlots; ++nodeId) {
            if (nodeValuation[static_cast<std::size_t>(nodeId)] >= threshold) {
                for (std::size_t i = contours.offsets[static_cast<std::size_t>(nodeId)]; i < contours.offsets[static_cast<std::size_t>(nodeId) + 1]; ++i) {
                    cut.sources.push_back(contours.sources[i]);
                    cut.targets.push_back(contours.targets[i]);
                }
            }
        }
        return cut;
    }
};
 
} // namespace mmcfilters