GrayLevelStatsComputer.hpp

#pragma once
 
#include "AttributeComputerDomain.hpp"
#include "AttributeComputerFamily.hpp"
#include "../detail/AttributeKernelSupport.hpp"
#include "../../trees/detail/TreeTraversalDetail.hpp"
#include "../../trees/TreeAltitudeAlgorithms.hpp"
#include "../../utils/Altitude.hpp"
 
#include <algorithm>
#include <array>
#include <cmath>
#include <span>
#include <stdexcept>
#include <string_view>
#include <vector>
 
namespace mmcfilters::attributes::computers {
 
namespace detail {
inline NodeId grayStatsSlotOf(const MorphologicalTree&, NodeId nodeId) noexcept {
    return nodeId;
}
 
struct GrayLevelStatsRequest {
    bool level = false;
    bool meanLevel = false;
    bool varianceLevel = false;
    bool grayHeight = false;
 
    [[nodiscard]] bool any() const noexcept {
        return level || meanLevel || varianceLevel || grayHeight;
    }
 
    [[nodiscard]] bool needsAggregateDependencies() const noexcept {
        return meanLevel || varianceLevel;
    }
 
    [[nodiscard]] static GrayLevelStatsRequest from(std::span<const Attribute> requestedAttributes) {
        return {
            .level = containsGrayStatsAttribute(requestedAttributes, LEVEL),
            .meanLevel = containsGrayStatsAttribute(requestedAttributes, MEAN_LEVEL),
            .varianceLevel = containsGrayStatsAttribute(requestedAttributes, VARIANCE_LEVEL),
            .grayHeight = containsGrayStatsAttribute(requestedAttributes, GRAY_HEIGHT)};
    }
 
private:
    [[nodiscard]] static bool containsGrayStatsAttribute(
        std::span<const Attribute> requestedAttributes,
        Attribute attribute)
    {
        return std::find(requestedAttributes.begin(), requestedAttributes.end(), attribute) != requestedAttributes.end();
    }
};
 
template<std::floating_point Real, AltitudeValue T>
void computeGrayLevelStatsAttributeKernel(
    const MorphologicalTree& tree,
    std::span<const T> altitude,
    std::span<Real> buffer,
    const AttributeNames& attrNames,
    std::span<const Attribute> requestedAttributes,
    std::span<const DependencySourceT<Real>> dependencySources)
{
    TreeAltitudeAlgorithms::validateAltitudeBufferShape(tree, altitude);
 
    const GrayLevelStatsRequest request = GrayLevelStatsRequest::from(requestedAttributes);
    if (!request.any()) {
        return;
    }
 
    const bool needsAggregateDependencies = request.needsAggregateDependencies();
 
    auto indexOfMean = [&](NodeId idx) { return attrNames.linearIndex(idx, MEAN_LEVEL); };
    auto indexOfLevel = [&](NodeId idx) { return attrNames.linearIndex(idx, LEVEL); };
    auto indexOfVariance = [&](NodeId idx) { return attrNames.linearIndex(idx, VARIANCE_LEVEL); };
    auto indexOfGrayHeight = [&](NodeId idx) { return attrNames.linearIndex(idx, GRAY_HEIGHT); };
 
    const DependencyResolver<Real> dependencies{dependencySources};
    const DependencySourceT<Real>* dependencyVol = needsAggregateDependencies
        ? &dependencies.require(VOLUME)
        : nullptr;
    const DependencySourceT<Real>* dependencyArea = needsAggregateDependencies
        ? &dependencies.require(AREA)
        : nullptr;
    auto indexOfVol = [&](NodeId idx) { return dependencyVol->attrNames->linearIndex(idx, VOLUME); };
    auto indexOfArea = [&](NodeId idx) { return dependencyArea->attrNames->linearIndex(idx, AREA); };
 
    std::vector<double> sumGrayLevelSquare;
    if (request.varianceLevel) {
        sumGrayLevelSquare.assign(tree.getNumInternalNodeSlots(), 0.0);
    }
 
    ::mmcfilters::detail::traversePostOrder(
        tree,
        tree.getRoot(),
        [&](NodeId nodeId) {
            const NodeId node = detail::grayStatsSlotOf(tree, nodeId);
            const T nodeAltitude = TreeAltitudeAlgorithms::getAltitude(altitude, nodeId);
            const Real nodeAltitudeAsReal = static_cast<Real>(nodeAltitude);
            if (request.varianceLevel) {
                const double nodeAltitudeAsDouble = static_cast<double>(nodeAltitude);
                sumGrayLevelSquare[node] =
                    static_cast<double>(tree.getNumProperParts(nodeId)) *
                    nodeAltitudeAsDouble *
                    nodeAltitudeAsDouble;
            }
            if (request.level)
                buffer[indexOfLevel(node)] = nodeAltitudeAsReal;
            if (request.grayHeight)
                buffer[indexOfGrayHeight(node)] = nodeAltitudeAsReal;
        },
        [&](NodeId parentNodeId, NodeId childNodeId) {
            const NodeId parent = detail::grayStatsSlotOf(tree, parentNodeId);
            const NodeId child = detail::grayStatsSlotOf(tree, childNodeId);
            if (request.varianceLevel)
                sumGrayLevelSquare[parent] += sumGrayLevelSquare[child];
            if (request.grayHeight) {
                Real childValue = buffer[indexOfGrayHeight(child)];
                Real& parentValue = buffer[indexOfGrayHeight(parent)];
                if (tree.getTreeType() == MorphologicalTreeKind::MAX_TREE)
                    parentValue = std::max(parentValue, childValue);
                else
                    parentValue = std::min(parentValue, childValue);
            }
        },
        [&](NodeId nodeId) {
            const NodeId node = detail::grayStatsSlotOf(tree, nodeId);
            Real area = needsAggregateDependencies ? dependencyArea->buffer[indexOfArea(node)] : Real{0};
            if (request.meanLevel)
                buffer[indexOfMean(node)] =
                    ::mmcfilters::attributes::numeric::safeDivide(dependencyVol->buffer[indexOfVol(node)], area);
            if (request.varianceLevel) {
                Real meanGrayLevel =
                    ::mmcfilters::attributes::numeric::safeDivide(dependencyVol->buffer[indexOfVol(node)], area);
                double meanGrayLevelSquare =
                    ::mmcfilters::attributes::numeric::safeDivide(sumGrayLevelSquare[node], static_cast<double>(area));
                Real var = static_cast<Real>(meanGrayLevelSquare - (static_cast<double>(meanGrayLevel) * static_cast<double>(meanGrayLevel)));
                buffer[indexOfVariance(node)] = ::mmcfilters::attributes::numeric::clampNonNegative(var);
            }
        });
 
    if (request.grayHeight) {
        for (NodeId nodeId : tree.getPostOrderNodes()) {
            const NodeId node = detail::grayStatsSlotOf(tree, nodeId);
            if (tree.isLeaf(nodeId)) {
                buffer[indexOfGrayHeight(node)] = 0.0f;
            } else {
                buffer[indexOfGrayHeight(node)] =
                    std::abs(static_cast<Real>(TreeAltitudeAlgorithms::getAltitude(altitude, nodeId)) -
                             buffer[indexOfGrayHeight(node)]);
            }
        }
    }
}
 
}
 
class GrayLevelStatsComputer {
public:
    static constexpr std::string_view familyName = "gray-level-stats";
 
    static constexpr AttributeComputerFamily family = AttributeComputerFamily::GrayLevelStats;
 
    static constexpr AttributeComputerDomain domain = AttributeComputerDomain::Altitude;
 
    inline static constexpr std::array<Attribute, 4> producedAttributes{
        LEVEL,
        MEAN_LEVEL,
        VARIANCE_LEVEL,
        GRAY_HEIGHT};
 
    template <std::floating_point Real, AltitudeValue T>
    static void compute(const AltitudeAttributeComputeContext<Real, T>& context)
    {
        requireAttributeBufferShape(context.tree, context.buffer, context.attrNames);
        detail::computeGrayLevelStatsAttributeKernel(
            context.tree,
            context.altitude,
            context.buffer,
            context.attrNames,
            context.requestedAttributes,
            context.dependencySources);
    }
 
    template <std::floating_point Real, AltitudeValue T>
    static void computeUnitRows(const AltitudeUnitAttributeComputeContext<Real, T>& context)
    {
        requireUnitAttributeBufferShape(context.tree, context.unitProperParts, context.buffer, context.attrNames);
        TreeAltitudeAlgorithms::validateAltitudeBufferShape(context.tree, context.altitude);
 
        const detail::GrayLevelStatsRequest request = detail::GrayLevelStatsRequest::from(context.requestedAttributes);
        if (!request.any()) {
            return;
        }
 
        for (NodeId leafIndex = 0; leafIndex < static_cast<NodeId>(context.unitProperParts.size()); ++leafIndex) {
            const NodeId properPart = context.unitProperParts[static_cast<size_t>(leafIndex)];
            const Real altitudeValue =
                (request.level || request.meanLevel)
                    ? static_cast<Real>(unitAltitude(context.tree, context.altitude, properPart))
                    : Real{0};
 
            if (request.level) {
                context.buffer[context.attrNames.linearIndex(leafIndex, LEVEL)] = altitudeValue;
            }
            if (request.meanLevel) {
                context.buffer[context.attrNames.linearIndex(leafIndex, MEAN_LEVEL)] = altitudeValue;
            }
            if (request.varianceLevel) {
                context.buffer[context.attrNames.linearIndex(leafIndex, VARIANCE_LEVEL)] = Real{0};
            }
            if (request.grayHeight) {
                context.buffer[context.attrNames.linearIndex(leafIndex, GRAY_HEIGHT)] = Real{0};
            }
        }
    }
 
};
 
} // namespace mmcfilters::attributes::computers