VolumeComputer.hpp

#pragma once
 
#include "AttributeComputerDomain.hpp"
#include "AttributeComputerFamily.hpp"
#include "../detail/AttributeKernelSupport.hpp"
#include "../../trees/detail/TreeTraversalDetail.hpp"
#include "../../trees/detail/CommittedTreeAccess.hpp"
#include "../../trees/TreeAltitudeAlgorithms.hpp"
#include "../../utils/Altitude.hpp"
#include "../../utils/Contract.hpp"
 
#include <algorithm>
#include <array>
#include <cmath>
#include <span>
#include <stdexcept>
#include <string_view>
 
namespace mmcfilters::attributes::computers {
 
namespace detail {
inline NodeId volumeSlotOf(const MorphologicalTree&, NodeId nodeId) noexcept { return nodeId; }
 
struct VolumeRequest {
    bool volume = false;
    bool relative = false;
 
    [[nodiscard]] bool any() const noexcept { return volume || relative; }
 
    [[nodiscard]] bool needsAreaDependency() const noexcept { return relative; }
 
    [[nodiscard]] static VolumeRequest from(std::span<const Attribute> requestedAttributes) {
        return {.volume = containsVolumeAttribute(requestedAttributes, Volume), .relative = containsVolumeAttribute(requestedAttributes, RelativeVolume)};
    }
 
  private:
    [[nodiscard]] static bool containsVolumeAttribute(std::span<const Attribute> requestedAttributes, Attribute attribute) {
        return std::find(requestedAttributes.begin(), requestedAttributes.end(), attribute) != requestedAttributes.end();
    }
};
 
namespace kernel {
 
template <std::floating_point Real, AltitudeValue T>
void computeVolume(const AltitudeAttributeComputeContext<Real, T>& context, const VolumeRequest& request,
                   const DependencySourceT<Real>* areaDependency) {
    if (!request.any()) {
        return;
    }
 
    const int stride = context.attrNames.NUM_ATTRIBUTES;
    const auto offsetOf = [&](Attribute attribute) { return context.attrNames.indexMap.find(attribute)->second; };
    const int volumeOffset = request.volume ? offsetOf(Volume) : 0;
    const int relativeOffset = request.relative ? offsetOf(RelativeVolume) : 0;
    auto indexOfVolume = [&](NodeId node) { return static_cast<std::size_t>(node * stride + volumeOffset); };
    auto indexOfRelative = [&](NodeId node) { return static_cast<std::size_t>(node * stride + relativeOffset); };
    const int areaStride = areaDependency != nullptr ? areaDependency->attrNames->NUM_ATTRIBUTES : 0;
    const int areaOffset = areaDependency != nullptr ? areaDependency->attrNames->indexMap.find(Area)->second : 0;
    auto indexOfArea = [&](NodeId node) { return static_cast<std::size_t>(node * areaStride + areaOffset); };
 
    ::mmcfilters::detail::kernel::traversePostOrder(
        context.tree, context.tree.root(),
        [&](NodeId nodeId) {
            const NodeId node = detail::volumeSlotOf(context.tree, nodeId);
            const T nodeAltitude = context.altitude[static_cast<std::size_t>(nodeId)];
            if (request.volume)
                context.buffer[indexOfVolume(node)] =
                    static_cast<Real>(::mmcfilters::detail::CommittedTreeAccess::properPartCardinality(context.tree, nodeId)) * static_cast<Real>(nodeAltitude);
            if (request.relative)
                context.buffer[indexOfRelative(node)] = Real{0};
        },
        [&](NodeId parentNodeId, NodeId childNodeId) {
            const NodeId parent = detail::volumeSlotOf(context.tree, parentNodeId);
            const NodeId child = detail::volumeSlotOf(context.tree, childNodeId);
            if (request.volume)
                context.buffer[indexOfVolume(parent)] += context.buffer[indexOfVolume(child)];
            if (request.relative)
                context.buffer[indexOfRelative(parent)] +=
                    context.buffer[indexOfRelative(child)] +
                    static_cast<Real>(areaDependency->buffer[indexOfArea(child)] *
                                      std::abs(static_cast<Real>(context.altitude[static_cast<std::size_t>(childNodeId)]) -
                                               static_cast<Real>(context.altitude[static_cast<std::size_t>(parentNodeId)])));
        },
        [&](NodeId nodeId) {
            const NodeId node = detail::volumeSlotOf(context.tree, nodeId);
            if (request.relative)
                context.buffer[indexOfRelative(node)] += areaDependency->buffer[indexOfArea(node)];
        });
}
 
} // namespace kernel
 
template <std::floating_point Real>
inline const DependencySourceT<Real>* findDependency(std::span<const DependencySourceT<Real>> sources, Attribute attribute) {
    for (const DependencySourceT<Real>& source : sources) {
        if (source.attrNames->contains(attribute)) {
            return &source;
        }
    }
    return nullptr;
}
 
template <std::floating_point Real, AltitudeValue T>
inline void validateVolumeContext(const AltitudeAttributeComputeContext<Real, T>& context, const VolumeRequest& request) {
    requireAttributeBufferShape(context.tree, context.buffer, context.attrNames);
    TreeAltitudeAlgorithms::validateNodeAltitudeBufferShape(context.tree, context.altitude);
    if (request.volume && !context.attrNames.contains(Volume)) {
        throw std::invalid_argument("VOLUME computation requires a VOLUME column in the output layout.");
    }
    if (request.relative && !context.attrNames.contains(RelativeVolume)) {
        throw std::invalid_argument("RELATIVE_VOLUME computation requires a RELATIVE_VOLUME column in the output layout.");
    }
}
 
} // namespace detail
 
class VolumeComputer {
  public:
    static constexpr std::string_view familyName = "volume";
 
    static constexpr AttributeComputerFamily family = AttributeComputerFamily::Volume;
 
    static constexpr AttributeComputerDomain domain = AttributeComputerDomain::Altitude;
 
    inline static constexpr std::array<Attribute, 2> producedAttributes{Volume, RelativeVolume};
 
    template <std::floating_point Real, AltitudeValue T> static void compute(const AltitudeAttributeComputeContext<Real, T>& context) {
        const detail::VolumeRequest request = detail::VolumeRequest::from(context.requestedAttributes);
        MMCFILTERS_CONTRACT_CHECKED_ONLY(detail::validateVolumeContext(context, request));
        const DependencySourceT<Real>* areaDependency = nullptr;
        if (request.needsAreaDependency()) {
            if constexpr (contract::validationsEnabled) {
                areaDependency = &context.dependencies.require(Area);
            } else {
                areaDependency = detail::findDependency(context.dependencySources, Area);
            }
        }
        detail::kernel::computeVolume(context, request, areaDependency);
    }
 
    template <std::floating_point Real, AltitudeValue T> static void computeUnitRows(const AltitudeUnitAttributeComputeContext<Real, T>& context) {
        requireUnitAttributeBufferShape(context.tree, context.unitPixels, context.buffer, context.attrNames);
        TreeAltitudeAlgorithms::validateNodeAltitudeBufferShape(context.tree, context.altitude);
 
        const detail::VolumeRequest request = detail::VolumeRequest::from(context.requestedAttributes);
        if (!request.any()) {
            return;
        }
 
        for (NodeId leafIndex = 0; leafIndex < static_cast<NodeId>(context.unitPixels.size()); ++leafIndex) {
            const PixelId pixel = context.unitPixels[static_cast<size_t>(leafIndex)];
            if (request.volume) {
                context.buffer[context.attrNames.linearIndex(leafIndex, Volume)] =
                    static_cast<Real>(::mmcfilters::detail::kernel::unitAltitude(context.tree, context.altitude, pixel));
            }
            if (request.relative) {
                context.buffer[context.attrNames.linearIndex(leafIndex, RelativeVolume)] = Real{1};
            }
        }
    }
};
 
} // namespace mmcfilters::attributes::computers