BoundingBoxComputer.hpp

#pragma once
 
#include "AttributeComputerDomain.hpp"
#include "AttributeComputerFamily.hpp"
#include "../detail/AttributeKernelSupport.hpp"
#include "../../trees/detail/TreeKernelAccess.hpp"
#include "../../trees/MorphologicalTree.hpp"
#include "../../utils/Contract.hpp"
#include "../../utils/Image.hpp"
 
#include <algorithm>
#include <array>
#include <concepts>
#include <cmath>
#include <span>
#include <string_view>
#include <vector>
 
namespace mmcfilters::attributes::computers {
 
namespace detail {
inline NodeId boundingBoxSlotOf(const MorphologicalTree&, NodeId nodeId) noexcept { return nodeId; }
 
struct BoundingBoxRequest {
    bool width = false;
    bool height = false;
    bool rectangularity = false;
    bool ratioWH = false;
    bool columnMin = false;
    bool columnMax = false;
    bool rowMin = false;
    bool rowMax = false;
    bool diagonalLength = false;
 
    [[nodiscard]] bool any() const noexcept { return width || height || rectangularity || ratioWH || columnMin || columnMax || rowMin || rowMax || diagonalLength; }
 
    [[nodiscard]] bool needsAreaDependency() const noexcept { return rectangularity; }
 
    [[nodiscard]] static BoundingBoxRequest from(std::span<const Attribute> requestedAttributes) {
        return {.width = contains(requestedAttributes, BoundingBoxWidth),
                .height = contains(requestedAttributes, BoundingBoxHeight),
                .rectangularity = contains(requestedAttributes, Rectangularity),
                .ratioWH = contains(requestedAttributes, RatioWh),
                .columnMin = contains(requestedAttributes, BoxColumnMin),
                .columnMax = contains(requestedAttributes, BoxColumnMax),
                .rowMin = contains(requestedAttributes, BoxRowMin),
                .rowMax = contains(requestedAttributes, BoxRowMax),
                .diagonalLength = contains(requestedAttributes, DiagonalLength)};
    }
 
  private:
    [[nodiscard]] static bool contains(std::span<const Attribute> requestedAttributes, Attribute attribute) {
        return std::find(requestedAttributes.begin(), requestedAttributes.end(), attribute) != requestedAttributes.end();
    }
};
 
namespace kernel {
 
template <std::floating_point Real>
inline void computeBoundingBox(const AttributeComputeContext<Real>& context, const BoundingBoxRequest& 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 widthOffset = request.width ? offsetOf(BoundingBoxWidth) : 0;
    const int heightOffset = request.height ? offsetOf(BoundingBoxHeight) : 0;
    const int rectangularityOffset = request.rectangularity ? offsetOf(Rectangularity) : 0;
    const int ratioOffset = request.ratioWH ? offsetOf(RatioWh) : 0;
    const int columnMinOffset = request.columnMin ? offsetOf(BoxColumnMin) : 0;
    const int columnMaxOffset = request.columnMax ? offsetOf(BoxColumnMax) : 0;
    const int rowMinOffset = request.rowMin ? offsetOf(BoxRowMin) : 0;
    const int rowMaxOffset = request.rowMax ? offsetOf(BoxRowMax) : 0;
    const int diagonalOffset = request.diagonalLength ? offsetOf(DiagonalLength) : 0;
    auto outputIndex = [&](NodeId node, int offset) { return static_cast<std::size_t>(node * stride + offset); };
 
    const int areaStride = areaDependency != nullptr ? areaDependency->attrNames->NUM_ATTRIBUTES : 0;
    const int areaOffset = areaDependency != nullptr ? areaDependency->attrNames->indexMap.find(Area)->second : 0;
    auto areaIndex = [&](NodeId node) { return static_cast<std::size_t>(node * areaStride + areaOffset); };
 
    const int numNodes = context.tree.numInternalNodeSlots();
    const GridDomain2D& domain = ::mmcfilters::detail::TreeKernelAccess::gridDomain2D(context.tree);
    std::vector<int> columnMin(static_cast<std::size_t>(numNodes), domain.columns);
    std::vector<int> columnMax(static_cast<std::size_t>(numNodes), 0);
    std::vector<int> rowMin(static_cast<std::size_t>(numNodes), domain.rows);
    std::vector<int> rowMax(static_cast<std::size_t>(numNodes), 0);
 
    ::mmcfilters::detail::TreeKernelAccess::traversePostOrder(
        context.tree, context.tree.root(),
        [&](NodeId node) {
            for (PixelId pixel : ::mmcfilters::detail::TreeKernelAccess::properParts(context.tree, node)) {
                const int row = pixel / domain.columns;
                const int column = pixel % domain.columns;
                columnMin[static_cast<std::size_t>(node)] = std::min(columnMin[static_cast<std::size_t>(node)], column);
                columnMax[static_cast<std::size_t>(node)] = std::max(columnMax[static_cast<std::size_t>(node)], column);
                rowMin[static_cast<std::size_t>(node)] = std::min(rowMin[static_cast<std::size_t>(node)], row);
                rowMax[static_cast<std::size_t>(node)] = std::max(rowMax[static_cast<std::size_t>(node)], row);
            }
        },
        [&](NodeId parent, NodeId child) {
            columnMin[static_cast<std::size_t>(parent)] =
                std::min(columnMin[static_cast<std::size_t>(parent)], columnMin[static_cast<std::size_t>(child)]);
            columnMax[static_cast<std::size_t>(parent)] =
                std::max(columnMax[static_cast<std::size_t>(parent)], columnMax[static_cast<std::size_t>(child)]);
            rowMin[static_cast<std::size_t>(parent)] = std::min(rowMin[static_cast<std::size_t>(parent)], rowMin[static_cast<std::size_t>(child)]);
            rowMax[static_cast<std::size_t>(parent)] = std::max(rowMax[static_cast<std::size_t>(parent)], rowMax[static_cast<std::size_t>(child)]);
        },
        [&](NodeId node) {
            const Real width = static_cast<Real>(columnMax[static_cast<std::size_t>(node)] - columnMin[static_cast<std::size_t>(node)] + 1);
            const Real height = static_cast<Real>(rowMax[static_cast<std::size_t>(node)] - rowMin[static_cast<std::size_t>(node)] + 1);
            if (request.width)
                context.buffer[outputIndex(node, widthOffset)] = width;
            if (request.height)
                context.buffer[outputIndex(node, heightOffset)] = height;
            if (request.rectangularity)
                context.buffer[outputIndex(node, rectangularityOffset)] =
                    ::mmcfilters::attributes::numeric::safeDivide(areaDependency->buffer[areaIndex(node)], width * height);
            if (request.ratioWH)
                context.buffer[outputIndex(node, ratioOffset)] =
                    ::mmcfilters::attributes::numeric::safeDivide(std::max(width, height), std::min(width, height));
            if (request.columnMin)
                context.buffer[outputIndex(node, columnMinOffset)] = static_cast<Real>(columnMin[static_cast<std::size_t>(node)]);
            if (request.columnMax)
                context.buffer[outputIndex(node, columnMaxOffset)] = static_cast<Real>(columnMax[static_cast<std::size_t>(node)]);
            if (request.rowMin)
                context.buffer[outputIndex(node, rowMinOffset)] = static_cast<Real>(rowMin[static_cast<std::size_t>(node)]);
            if (request.rowMax)
                context.buffer[outputIndex(node, rowMaxOffset)] = static_cast<Real>(rowMax[static_cast<std::size_t>(node)]);
            if (request.diagonalLength)
                context.buffer[outputIndex(node, diagonalOffset)] = ::mmcfilters::attributes::numeric::safeSqrt(width * width + height * height);
        });
}
 
} // namespace kernel
 
template <std::floating_point Real>
inline void validateBoundingBoxContext(const AttributeComputeContext<Real>& context) {
    requireAttributeBufferShape(context.tree, context.buffer, context.attrNames);
    requireRequestedAttributeColumns(context);
    static_cast<void>(context.tree.requireGridDomain2D("BoundingBoxComputer"));
}
} // namespace detail
 
class BoundingBoxComputer {
  public:
    static constexpr std::string_view familyName = "bounding-box";
 
    static constexpr AttributeComputerFamily family = AttributeComputerFamily::BoundingBox;
 
    static constexpr AttributeComputerDomain domain = AttributeComputerDomain::Topology;
 
    inline static constexpr std::array<Attribute, 9> producedAttributes{BoundingBoxWidth,   BoundingBoxHeight,  DiagonalLength, Rectangularity, RatioWh,
                                                                        BoxColumnMin, BoxColumnMax, BoxRowMin,     BoxRowMax};
 
    template <std::floating_point Real> static void compute(const AttributeComputeContext<Real>& context) {
        const detail::BoundingBoxRequest request = detail::BoundingBoxRequest::from(context.requestedAttributes);
        MMCFILTERS_CONTRACT_CHECKED_ONLY(detail::validateBoundingBoxContext(context));
 
        const DependencySourceT<Real>* areaDependency = nullptr;
        if (request.needsAreaDependency()) {
            if constexpr (contract::validationsEnabled) {
                areaDependency = &context.dependencies.require(Area);
            } else {
                areaDependency = ::mmcfilters::findDependencySource(context.dependencySources, Area);
            }
        }
        detail::kernel::computeBoundingBox(context, request, areaDependency);
    }
 
    template <std::floating_point Real> static void computeUnitRows(const UnitAttributeComputeContext<Real>& context) {
        requireUnitAttributeBufferShape(context.tree, context.unitPixels, context.buffer, context.attrNames);
        computeUnitRowsEstablished(context);
    }
 
  private:
    template <class Computer, class Context> friend void mmcfilters::detail::kernel::computeUnitRows(const Context& context);
 
    template <std::floating_point Real> static void computeUnitRowsEstablished(const UnitAttributeComputeContext<Real>& context) {
 
        const detail::BoundingBoxRequest request = detail::BoundingBoxRequest::from(context.requestedAttributes);
        if (!request.any()) {
            return;
        }
 
        const int numColumns = context.tree.numColumns();
        for (NodeId leafIndex = 0; leafIndex < static_cast<NodeId>(context.unitPixels.size()); ++leafIndex) {
            const PixelId pixel = context.unitPixels[static_cast<size_t>(leafIndex)];
            const auto [row, column] = ImageUtils::to2D(pixel, numColumns);
            if (request.width) {
                context.buffer[context.attrNames.linearIndex(leafIndex, BoundingBoxWidth)] = Real{1};
            }
            if (request.height) {
                context.buffer[context.attrNames.linearIndex(leafIndex, BoundingBoxHeight)] = Real{1};
            }
            if (request.rectangularity) {
                context.buffer[context.attrNames.linearIndex(leafIndex, Rectangularity)] = Real{1};
            }
            if (request.ratioWH) {
                context.buffer[context.attrNames.linearIndex(leafIndex, RatioWh)] = Real{1};
            }
            if (request.columnMin) {
                context.buffer[context.attrNames.linearIndex(leafIndex, BoxColumnMin)] = static_cast<Real>(column);
            }
            if (request.columnMax) {
                context.buffer[context.attrNames.linearIndex(leafIndex, BoxColumnMax)] = static_cast<Real>(column);
            }
            if (request.rowMin) {
                context.buffer[context.attrNames.linearIndex(leafIndex, BoxRowMin)] = static_cast<Real>(row);
            }
            if (request.rowMax) {
                context.buffer[context.attrNames.linearIndex(leafIndex, BoxRowMax)] = static_cast<Real>(row);
            }
            if (request.diagonalLength) {
                context.buffer[context.attrNames.linearIndex(leafIndex, DiagonalLength)] = std::sqrt(Real{2});
            }
        }
    }
};
 
} // namespace mmcfilters::attributes::computers