ResidualEvolution.hpp

Go to the documentation of this file.

#pragma once
 
#include "../../utils/Altitude.hpp"
#include "../../utils/Common.hpp"
 
#include <cstddef>
#include <limits>
#include <set>
#include <span>
#include <stdexcept>
#include <utility>
#include <vector>
 
namespace mmcfilters::sdrt {
 
enum class Polarity { Maximum, Minimum };
 
class SpatialOrder {
  public:
    explicit SpatialOrder(std::vector<PixelId> pixelsInOrder)
        : rankByPixel_(pixelsInOrder.size(), invalidRank()) {
        for (std::size_t rank = 0; rank < pixelsInOrder.size(); ++rank) {
            const PixelId pixel = pixelsInOrder[rank];
            if (pixel < 0 || static_cast<std::size_t>(pixel) >= pixelsInOrder.size()) {
                throw std::invalid_argument("A spatial order must be a permutation of the pixel domain.");
            }
            auto& storedRank = rankByPixel_[static_cast<std::size_t>(pixel)];
            if (storedRank != invalidRank()) {
                throw std::invalid_argument("A spatial order cannot contain a pixel more than once.");
            }
            storedRank = rank;
        }
    }
 
    [[nodiscard]] bool isRowMajor() const noexcept { return rankByPixel_.empty(); }
 
    void validateDomainSize(std::size_t numPixels) const {
        if (!isRowMajor() && rankByPixel_.size() != numPixels) {
            throw std::invalid_argument("The explicit spatial order size differs from the image domain.");
        }
    }
 
    [[nodiscard]] bool precedes(PixelId lhs, PixelId rhs) const {
        return rank(lhs) < rank(rhs);
    }
 
    [[nodiscard]] PixelId spatialMinimum(std::span<const PixelId> support) const {
        if (support.empty()) {
            throw std::invalid_argument("A spatial minimum requires a non-empty support.");
        }
        PixelId minimum = support.front();
        for (PixelId pixel : support.subspan(1)) {
            if (precedes(pixel, minimum)) {
                minimum = pixel;
            }
        }
        return minimum;
    }
 
  protected:
    SpatialOrder() = default;
 
  private:
    [[nodiscard]] static constexpr std::size_t invalidRank() noexcept {
        return std::numeric_limits<std::size_t>::max();
    }
 
    [[nodiscard]] std::size_t rank(PixelId pixel) const {
        if (pixel < 0) {
            throw std::out_of_range("A spatial-order query received a negative pixel identifier.");
        }
        if (isRowMajor()) {
            return static_cast<std::size_t>(pixel);
        }
        if (static_cast<std::size_t>(pixel) >= rankByPixel_.size()) {
            throw std::out_of_range("A spatial-order query lies outside the configured domain.");
        }
        return rankByPixel_[static_cast<std::size_t>(pixel)];
    }
 
    std::vector<std::size_t> rankByPixel_; 
};
 
class RowMajorSpatialOrder final : public SpatialOrder {
  public:
    RowMajorSpatialOrder() = default;
};
 
struct SelfDualResidualKey {
    std::size_t supportCardinality = 0; 
    PixelId spatialMinimum = InvalidPixel; 
 
    friend bool operator==(const SelfDualResidualKey&, const SelfDualResidualKey&) = default;
};
 
class SelfDualResidualOrder {
  public:
    explicit SelfDualResidualOrder(SpatialOrder spatialOrder = RowMajorSpatialOrder{})
        : spatialOrder_(std::move(spatialOrder)) {}
 
    [[nodiscard]] bool compareResidualCandidates(const SelfDualResidualKey& lhs, const SelfDualResidualKey& rhs) const {
        if (lhs.supportCardinality != rhs.supportCardinality) {
            return lhs.supportCardinality < rhs.supportCardinality;
        }
        if (lhs.spatialMinimum == rhs.spatialMinimum) {
            return false;
        }
        return spatialOrder_.precedes(lhs.spatialMinimum, rhs.spatialMinimum);
    }
 
    [[nodiscard]] bool operator()(const SelfDualResidualKey& lhs, const SelfDualResidualKey& rhs) const {
        return compareResidualCandidates(lhs, rhs);
    }
 
    [[nodiscard]] const SpatialOrder& spatialOrder() const noexcept { return spatialOrder_; }
 
  private:
    SpatialOrder spatialOrder_; 
};
 
class SelfDualResidualSchedule {
  public:
    explicit SelfDualResidualSchedule(SpatialOrder spatialOrder = RowMajorSpatialOrder{})
        : residualOrder_(std::move(spatialOrder)) {}
 
    [[nodiscard]] std::size_t selectResidualCandidate(std::span<const SelfDualResidualKey> residualKeys) const {
        if (residualKeys.empty()) {
            throw std::invalid_argument("Residual-candidate selection requires a non-empty schedule.");
        }
        std::set<std::pair<std::size_t, PixelId>> seenKeys;
        std::size_t selected = 0;
        for (std::size_t index = 0; index < residualKeys.size(); ++index) {
            if (!seenKeys.emplace(residualKeys[index].supportCardinality, residualKeys[index].spatialMinimum).second) {
                throw std::invalid_argument("A self-dual residual schedule cannot contain duplicate candidate keys.");
            }
            if (index != 0 && residualOrder_.compareResidualCandidates(residualKeys[index], residualKeys[selected])) {
                selected = index;
            }
        }
        return selected;
    }
 
    [[nodiscard]] const SelfDualResidualOrder& residualOrder() const noexcept { return residualOrder_; }
 
  private:
    SelfDualResidualOrder residualOrder_; 
};
 
template <AltitudeValue T> struct ResidualCandidate {
    std::span<const PixelId> support; 
    Polarity polarity = Polarity::Maximum; 
    T candidateAltitude{}; 
    T firstMergingLevel{}; 
    SelfDualResidualKey residualKey; 
};
 
template <AltitudeValue T> struct ResidualEvent {
    std::size_t eventIndex = 0; 
    std::span<const PixelId> support; 
    Polarity polarity = Polarity::Maximum; 
    T nodeAltitude{}; 
    T firstMergingLevel{}; 
    AltitudeDifference<T> signedResidualValue{}; 
};
 
template <AltitudeValue T>
[[nodiscard]] ResidualEvent<T> recordResidualEvent(std::size_t eventIndex, const ResidualCandidate<T>& candidate) {
    if (candidate.support.empty()) {
        throw std::invalid_argument("A residual event requires a non-empty candidate support.");
    }
    if (candidate.residualKey.supportCardinality != candidate.support.size()) {
        throw std::invalid_argument("A residual candidate key has an inconsistent support cardinality.");
    }
    const auto signedResidualValue = static_cast<AltitudeDifference<T>>(candidate.candidateAltitude) -
                                     static_cast<AltitudeDifference<T>>(candidate.firstMergingLevel);
    if ((candidate.polarity == Polarity::Maximum && signedResidualValue <= AltitudeDifference<T>{}) ||
        (candidate.polarity == Polarity::Minimum && signedResidualValue >= AltitudeDifference<T>{})) {
        throw std::invalid_argument("A residual candidate polarity is inconsistent with its signed residual value.");
    }
    return ResidualEvent<T>{eventIndex, candidate.support, candidate.polarity, candidate.candidateAltitude,
                            candidate.firstMergingLevel, signedResidualValue};
}
 
} // namespace mmcfilters::sdrt