UltimateAttributeOpening.hpp

#pragma once
 
#include "../utils/Image.hpp"
#include "../utils/Common.hpp"
#include "../trees/ValuedMorphologicalTree.hpp"
#include "../trees/ValuedMorphologicalTreeView.hpp"
#include "../trees/detail/CommittedTreeAccess.hpp"
#include "../utils/Contract.hpp"
#include "DepthStableRegionComputer.hpp"
#include "MSERComputer.hpp"
 
#include <cmath>
#include <concepts>
#include <memory>
#include <stack>
#include <stdexcept>
#include <string>
#include <vector>
 
namespace mmcfilters {
 
template <AltitudeValue T, std::floating_point Real = float> class UltimateAttributeOpening {
  public:
    using attribute_value_type = Real;
 
  protected:
 
    using AltitudeView = ValuedMorphologicalTreeView<T>;
 
    Real maximumAttributeThreshold = Real{0};
    const Real* const attrs_increasing;
    std::shared_ptr<Real[]> ownedAttrsIncreasing_;
    const ValuedMorphologicalTree<T>* valuedTree_ = nullptr;
    AltitudeView view_;
    const MorphologicalTree& tree;
    std::size_t treeMutationVersion_ = 0;
    std::vector<T> maxContrastLUT;
    std::vector<int> associatedIndexLUT;
    std::vector<uint8_t> selectedForFiltering;
 
    AltitudeView view() const { return valuedTree_ != nullptr ? valuedTree_->asView() : view_; }
 
    void requireStableTree(const char* context) const { tree.requireMutationVersion(treeMutationVersion_, context); }
 
    static T altitudeOf(const AltitudeView& view, NodeId nodeId) noexcept {
        return view.nodeAltitudes()[static_cast<std::size_t>(nodeId)];
    }
 
    static void requireAttributePointer(const Real* attr, const char* context) {
        MMCFILTERS_CONTRACT_REQUIRE(attr != nullptr,
                                    throw std::invalid_argument(std::string(context) + " requires a non-null attribute buffer."));
    }
 
    static const Real* requireAttributeBuffer(const MorphologicalTree& tree, const std::vector<Real>& attr, const char* context) {
        MMCFILTERS_CONTRACT_REQUIRE(attr.size() == static_cast<std::size_t>(tree.numInternalNodeSlots()),
                                    throw std::invalid_argument(std::string(context) + " attribute size must match the internal node slot count."));
        return attr.data();
    }
 
    static T absoluteAltitudeDifference(AltitudeDifference<T> lhs, AltitudeDifference<T> rhs) {
        const long double difference = std::abs(static_cast<long double>(lhs) - static_cast<long double>(rhs));
        // UAO stores the contrast in the same type as the altitude by API
        // decision. For signed integral altitudes, contrasts that do not fit in
        // T can lose information here; avoiding that requires a separate
        // contrast type, which this API intentionally does not introduce.
        return static_cast<T>(difference);
    }
 
    void computeUAO(const AltitudeView& view, NodeId currentNodeId, AltitudeDifference<T> altitudeNodeNotInNR, bool qPropag, bool isCalculateResidue) {
        const NodeId parentNodeId = detail::CommittedTreeAccess::nodeParent(tree, currentNodeId);
        const AltitudeDifference<T> altitudeNodeInNR = static_cast<AltitudeDifference<T>>(altitudeOf(view, currentNodeId));
        bool flagPropag = false;
        T contrast = T{};
 
        if (this->isSelectedForPruning(currentNodeId)) {
            altitudeNodeNotInNR = static_cast<AltitudeDifference<T>>(altitudeOf(view, parentNodeId));
            if (this->attrs_increasing[currentNodeId] <= this->maximumAttributeThreshold) {
                isCalculateResidue = hasNodeSelectedInPrimitive(currentNodeId);
            }
        }
 
        if (this->attrs_increasing[currentNodeId] <= this->maximumAttributeThreshold) {
            if (isCalculateResidue) {
                contrast = absoluteAltitudeDifference(altitudeNodeInNR, altitudeNodeNotInNR);
            }
 
            if (this->maxContrastLUT[parentNodeId] >= contrast) {
                this->maxContrastLUT[currentNodeId] = this->maxContrastLUT[parentNodeId];
                this->associatedIndexLUT[currentNodeId] = this->associatedIndexLUT[parentNodeId];
            } else {
                this->maxContrastLUT[currentNodeId] = contrast;
                this->associatedIndexLUT[currentNodeId] =
                    !qPropag ? static_cast<int>(this->attrs_increasing[currentNodeId] + 1) : this->associatedIndexLUT[parentNodeId];
                flagPropag = true;
            }
        }
 
        for (NodeId childNodeId : detail::CommittedTreeAccess::children(tree, currentNodeId)) {
            this->computeUAO(view, childNodeId, altitudeNodeNotInNR, flagPropag, isCalculateResidue);
        }
    }
 
    void executeImpl(Real maximumAttributeThreshold, const std::vector<uint8_t>& selectedForFiltering) {
        const AltitudeView altitudeView = view();
        this->maximumAttributeThreshold = maximumAttributeThreshold;
        this->selectedForFiltering = selectedForFiltering;
 
        for (NodeId id : tree.aliveNodeIds()) {
            maxContrastLUT[id] = T{};
            associatedIndexLUT[id] = 0;
        }
 
        const NodeId rootNodeId = tree.root();
        const AltitudeDifference<T> level = static_cast<AltitudeDifference<T>>(altitudeOf(altitudeView, rootNodeId));
        for (NodeId childNodeId : detail::CommittedTreeAccess::children(tree, rootNodeId)) {
            computeUAO(altitudeView, childNodeId, level, false, false);
        }
    }
 
    bool isSelectedForPruning(NodeId currentNodeId) const {
        const NodeId parentNodeId = detail::CommittedTreeAccess::nodeParent(tree, currentNodeId);
        if (parentNodeId == InvalidNode) {
            return false;
        }
        return this->attrs_increasing[currentNodeId] != this->attrs_increasing[parentNodeId];
    }
 
    bool hasNodeSelectedInPrimitive(NodeId currentNodeId) const {
        std::stack<NodeId> stack;
        stack.push(currentNodeId);
        while (!stack.empty()) {
            const NodeId nodeId = stack.top();
            stack.pop();
            if (selectedForFiltering[nodeId]) {
                return true;
            }
 
            for (NodeId childNodeId : detail::CommittedTreeAccess::children(tree, nodeId)) {
                if (this->attrs_increasing[childNodeId] == this->attrs_increasing[nodeId]) {
                    stack.push(childNodeId);
                }
            }
        }
        return false;
    }
 
  public:
    UltimateAttributeOpening(const AltitudeView& view, const std::shared_ptr<Real[]>& attrs_increasing)
        : UltimateAttributeOpening(view, attrs_increasing.get()) {
        this->ownedAttrsIncreasing_ = attrs_increasing;
    }
 
    UltimateAttributeOpening(const AltitudeView& view, const std::vector<Real>& attrs_increasing)
        : UltimateAttributeOpening(view, requireAttributeBuffer(view.topology(), attrs_increasing, "UltimateAttributeOpening")) {}
 
    UltimateAttributeOpening(const AltitudeView& view, const Real* attrs_increasing)
        : attrs_increasing(attrs_increasing), view_(view), tree(view_.topology()), treeMutationVersion_(tree.getMutationVersion()),
          maxContrastLUT(this->tree.numInternalNodeSlots()), associatedIndexLUT(this->tree.numInternalNodeSlots()) {
        view_.requireTopologyUnchanged("UltimateAttributeOpening");
        requireAttributePointer(attrs_increasing, "UltimateAttributeOpening");
        this->selectedForFiltering.assign(this->tree.numInternalNodeSlots(), true);
    }
 
    UltimateAttributeOpening(const ValuedMorphologicalTree<T>& valuedTree, const std::shared_ptr<Real[]>& attrs_increasing)
        : UltimateAttributeOpening(valuedTree, attrs_increasing.get()) {
        this->ownedAttrsIncreasing_ = attrs_increasing;
    }
 
    UltimateAttributeOpening(const ValuedMorphologicalTree<T>& valuedTree, const std::vector<Real>& attrs_increasing)
        : UltimateAttributeOpening(valuedTree.asView(), attrs_increasing) {
        valuedTree_ = &valuedTree;
    }
 
    UltimateAttributeOpening(const ValuedMorphologicalTree<T>& valuedTree, const Real* attrs_increasing)
        : UltimateAttributeOpening(valuedTree.asView(), attrs_increasing) {
        valuedTree_ = &valuedTree;
    }
 
    ~UltimateAttributeOpening() = default;
 
  public:
    void execute(Real maximumAttributeThreshold) {
        requireStableTree("UltimateAttributeOpening::execute");
        std::vector<uint8_t> tmp(this->tree.numInternalNodeSlots(), true);
        executeImpl(maximumAttributeThreshold, tmp);
    }
 
    void execute(Real maximumAttributeThreshold, const std::vector<uint8_t>& selectedForFiltering) {
        requireStableTree("UltimateAttributeOpening::execute");
        MMCFILTERS_CONTRACT_REQUIRE(
            selectedForFiltering.size() == static_cast<std::size_t>(this->tree.numInternalNodeSlots()),
            throw std::invalid_argument("UltimateAttributeOpening::execute selectedForFiltering size must match the internal node slot count."));
        executeImpl(maximumAttributeThreshold, selectedForFiltering);
    }
 
    void executeWithMSER(Real maximumAttributeThreshold, AltitudeDifference<T> altitudeWindowRadius) {
        requireStableTree("UltimateAttributeOpening::executeWithMSER");
        if (valuedTree_ == nullptr) {
            throw std::logic_error(
                "UltimateAttributeOpening::executeWithMSER requires a ValuedMorphologicalTree owner because MSER uses the tree-owned altitude.");
        }
        MSERComputer<T, Real> mser(*valuedTree_);
        executeImpl(maximumAttributeThreshold, mser.computeMSER(altitudeWindowRadius));
    }
 
    void executeWithDepthStability(Real maximumAttributeThreshold, int depthWindowRadius) {
        requireStableTree("UltimateAttributeOpening::executeWithDepthStability");
        DepthStableRegionComputer<Real> stabilityComputer(this->tree);
        executeImpl(maximumAttributeThreshold, stabilityComputer.computeByDepth(depthWindowRadius));
    }
 
    [[nodiscard]] ImagePtr<T> getMaxContrastImage() const {
        requireStableTree("UltimateAttributeOpening::getMaxContrastImage");
        const int size = this->tree.numColumns() * this->tree.numRows();
        ImagePtr<T> imgOut = Image<T>::create(this->tree.numRows(), this->tree.numColumns());
        auto out = imgOut->rawData();
 
        for (int pidx = 0; pidx < size; pidx++) {
            out[pidx] = this->maxContrastLUT[detail::CommittedTreeAccess::smallestNodeMap(tree, pidx)];
        }
        return imgOut;
    }
 
    [[nodiscard]] ImageInt32Ptr getAssociatedImage() const {
        requireStableTree("UltimateAttributeOpening::getAssociatedImage");
        const int size = this->tree.numColumns() * this->tree.numRows();
        ImageInt32Ptr imgOut = ImageInt32::create(this->tree.numRows(), this->tree.numColumns());
        auto out = imgOut->rawData();
 
        for (int pidx = 0; pidx < size; pidx++) {
            out[pidx] = this->associatedIndexLUT[detail::CommittedTreeAccess::smallestNodeMap(tree, pidx)];
        }
        return imgOut;
    }
 
    [[nodiscard]] ImageUInt8Ptr getAssociatedColorImage() const {
        return ImageUtils::createRandomColor(this->getAssociatedImage()->rawData(), this->tree.numRows(),
                                             this->tree.numColumns());
    }
};
 
} // namespace mmcfilters