VolumeComputer.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>
 
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, RELATIVE_VOLUME)};
    }
 
private:
    [[nodiscard]] static bool containsVolumeAttribute(
        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 computeVolumeAttributeKernel(
    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 VolumeRequest request = VolumeRequest::from(requestedAttributes);
    if (!request.any()) {
        return;
    }
 
    auto indexOfVol = [&](NodeId idx) { return attrNames.linearIndex(idx, VOLUME); };
    auto indexOfRel = [&](NodeId idx) { return attrNames.linearIndex(idx, RELATIVE_VOLUME); };
    const DependencyResolver<Real> dependencies{dependencySources};
    const DependencySourceT<Real>* dependencyArea = request.needsAreaDependency()
        ? &dependencies.require(AREA)
        : nullptr;
    auto indexOfArea = [&](NodeId idx) { return dependencyArea->attrNames->linearIndex(idx, AREA); };
 
    ::mmcfilters::detail::traversePostOrder(
        tree,
        tree.getRoot(),
        [&](NodeId nodeId) {
            const NodeId node = detail::volumeSlotOf(tree, nodeId);
            const T nodeAltitude = TreeAltitudeAlgorithms::getAltitude(altitude, nodeId);
            if (request.volume)
                buffer[indexOfVol(node)] =
                    static_cast<Real>(tree.getNumProperParts(nodeId)) * static_cast<Real>(nodeAltitude);
            if (request.relative)
                buffer[indexOfRel(node)] = Real{0};
        },
        [&](NodeId parentNodeId, NodeId childNodeId) {
            const NodeId parent = detail::volumeSlotOf(tree, parentNodeId);
            const NodeId child = detail::volumeSlotOf(tree, childNodeId);
            if (request.volume)
                buffer[indexOfVol(parent)] += buffer[indexOfVol(child)];
            if (request.relative)
                buffer[indexOfRel(parent)] +=
                    buffer[indexOfRel(child)] +
                    static_cast<Real>(
                        dependencyArea->buffer[indexOfArea(child)] *
                        std::abs(
                            static_cast<Real>(TreeAltitudeAlgorithms::getAltitude(altitude, childNodeId)) -
                            static_cast<Real>(TreeAltitudeAlgorithms::getAltitude(altitude, parentNodeId))));
        },
        [&](NodeId nodeId) {
            const NodeId node = detail::volumeSlotOf(tree, nodeId);
            if (request.relative)
                buffer[indexOfRel(node)] += dependencyArea->buffer[indexOfArea(node)];
        });
}
 
}
 
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,
        RELATIVE_VOLUME};
 
    template <std::floating_point Real, AltitudeValue T>
    static void compute(const AltitudeAttributeComputeContext<Real, T>& context)
    {
        requireAttributeBufferShape(context.tree, context.buffer, context.attrNames);
        detail::computeVolumeAttributeKernel(
            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::VolumeRequest request = detail::VolumeRequest::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)];
            if (request.volume) {
                context.buffer[context.attrNames.linearIndex(leafIndex, VOLUME)] =
                    static_cast<Real>(unitAltitude(context.tree, context.altitude, properPart));
            }
            if (request.relative) {
                context.buffer[context.attrNames.linearIndex(leafIndex, RELATIVE_VOLUME)] = Real{1};
            }
        }
    }
 
};
 
} // namespace mmcfilters::attributes::computers