FilledShapeAttributeComputer.hpp

#pragma once
 
#include "AttributeComputerDomain.hpp"
#include "AttributeComputerFamily.hpp"
#include "../detail/AttributeKernelSupport.hpp"
#include "../../contours/ContourTraceComputation.hpp"
 
#include <algorithm>
#include <array>
#include <cmath>
#include <concepts>
#include <limits>
#include <numbers>
#include <span>
#include <string_view>
 
namespace mmcfilters::attributes::computers {
 
namespace detail {
 
struct FilledShapeMoments {
    long double area = 0;       
    long double column = 0;     
    long double row = 0;        
    long double column2 = 0;    
    long double row2 = 0;       
    long double columnRow = 0;  
 
    void addBoundary(ContourTraceView trace, const ContourBoundary& boundary, int columns,
                     int originColumn, int originRow, bool secondOrder) {
        for (ContourEdge edge : trace.boundaryEdges(boundary)) {
            if (edge.side != ContourSide::East && edge.side != ContourSide::West) {
                continue;
            }
            const bool east = edge.side == ContourSide::East;
            const long double k = static_cast<long double>(edge.pixel % columns - originColumn) + (east ? 1 : 0);
            const long double y = edge.pixel / columns - originRow;
            const long double signedWidth = east ? k : -k;
            const long double columnSum = signedWidth * (k - 1) / 2;
            area += signedWidth;
            column += columnSum;
            row += signedWidth * y;
            if (secondOrder) {
                column2 += columnSum * (2 * k - 1) / 3;
                row2 += signedWidth * y * y;
                columnRow += columnSum * y;
            }
        }
    }
};
 
} // namespace detail
 
class FilledShapeAttributeComputer {
  public:
    static constexpr std::string_view familyName = "filled-shape"; 
    static constexpr AttributeComputerFamily family = AttributeComputerFamily::FilledShape; 
    static constexpr AttributeComputerDomain domain = AttributeComputerDomain::Topology; 
 
    inline static constexpr std::array<Attribute, 13> producedAttributes{
        FilledArea, FilledCentroidRow, FilledCentroidColumn, FilledLengthMajorAxis, FilledLengthMinorAxis,
        FilledAxisOrientation, FilledEccentricity, FilledInertia, HoleAreaFraction, FilledCentroidDisplacementNormalized,
        FilledCompactness, FilledCircularity, FilledAxisOrientationSigned};
 
    template <std::floating_point Real> static void compute(const AttributeComputeContext<Real>& context) {
        MMCFILTERS_CONTRACT_CHECKED_ONLY(validateContext(context));
        if (!context.requestedAttributes.empty()) {
            computeTraces(context, ContourTraceComputation(context.tree));
        }
    }
 
    template <std::floating_point Real, AltitudeValue T> static void compute(const AltitudeAttributeComputeContext<Real, T>& context) {
        const AttributeComputeContext<Real> topologyContext{context.tree, context.buffer, context.attrNames, context.requestedAttributes};
        MMCFILTERS_CONTRACT_CHECKED_ONLY(validateContext(topologyContext));
        if (!context.requestedAttributes.empty()) {
            const ValuedMorphologicalTreeView<T> view(context.tree, context.altitude);
            computeTraces(topologyContext, ContourTraceComputation(view));
        }
    }
 
    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 int columns = context.tree.requireGridDomain2D("FilledShapeAttributeComputer").columns;
        const auto offsets = outputOffsets(context.attrNames, context.requestedAttributes);
        for (NodeId row = 0; row < static_cast<NodeId>(context.unitPixels.size()); ++row) {
            const PixelId pixel = context.unitPixels[static_cast<std::size_t>(row)];
            const std::array<long double, 13> values{1, static_cast<long double>(pixel / columns),
                static_cast<long double>(pixel % columns), 0, 0, 0, 1, 0, 0, 0, 0, 1};
            writeRow(context.buffer, context.attrNames.NUM_ATTRIBUTES, row, offsets, values);
        }
    }
 
  private:
    template <std::floating_point Real> static void validateContext(const AttributeComputeContext<Real>& context) {
        requireAttributeBufferShape(context.tree, context.buffer, context.attrNames);
        requireRequestedAttributeColumns(context);
    }
 
    static std::array<int, 13> outputOffsets(const AttributeNames& names, std::span<const Attribute> requested) {
        std::array<int, 13> offsets;
        offsets.fill(-1);
        for (Attribute attribute : requested) {
            const auto it = std::find(producedAttributes.begin(), producedAttributes.end(), attribute);
            if (it == producedAttributes.end()) {
                throw std::invalid_argument("FilledShapeAttributeComputer received an attribute from another family.");
            }
            offsets[static_cast<std::size_t>(it - producedAttributes.begin())] = names.linearIndex(0, attribute);
        }
        return offsets;
    }
 
    template <std::floating_point Real>
    static void writeRow(std::span<Real> buffer, int stride, NodeId node, const std::array<int, 13>& offsets,
                         const std::array<long double, 13>& values) {
        const std::size_t row = static_cast<std::size_t>(node) * static_cast<std::size_t>(stride);
        for (std::size_t i = 0; i < offsets.size(); ++i) {
            if (offsets[i] >= 0) {
                buffer[row + static_cast<std::size_t>(offsets[i])] = static_cast<Real>(values[i]);
            }
        }
    }
 
    template <std::floating_point Real>
    static void computeTraces(const AttributeComputeContext<Real>& context, const ContourTraceComputation& traces) {
        const auto offsets = outputOffsets(context.attrNames, context.requestedAttributes);
        const bool secondOrder = offsets[10] >= 0 || offsets[11] >= 0 || offsets[12] >= 0 ||
            std::any_of(offsets.begin() + 3, offsets.begin() + 8, [](int offset) { return offset >= 0; });
        const bool compareSupport = offsets[8] >= 0 || offsets[9] >= 0;
        const int columns = context.tree.numColumns();
        traces.forEachTrace([&](NodeId node, ContourTraceView trace) {
            const ContourBoundary external = trace.externalBoundary();
            const PixelId origin = (*trace.boundaryEdges(external).begin()).pixel;
            const int originColumn = origin % columns;
            const int originRow = origin / columns;
            detail::FilledShapeMoments filled;
            filled.addBoundary(trace, external, columns, originColumn, originRow, secondOrder);
            const long double area = filled.area;
            const long double centroidColumn = filled.column / area;
            const long double centroidRow = filled.row / area;
            std::array<long double, 13> values{area, centroidRow + originRow, centroidColumn + originColumn, 0, 0, 0, 1, 0, 0, 0, 0, 1};
            if (secondOrder) {
                const long double mu20 = std::max(0.L, filled.column2 - filled.column * centroidColumn);
                const long double mu02 = std::max(0.L, filled.row2 - filled.row * centroidRow);
                const long double mu11 = filled.columnRow - filled.column * centroidRow;
                const long double discriminant = std::hypot(mu20 - mu02, 2 * mu11);
                const long double lambda1 = mu20 + mu02 + discriminant;
                const long double lambda2 = std::max(0.L, mu20 + mu02 - discriminant);
                values[3] = std::sqrt(2 * lambda1 / area);
                values[4] = std::sqrt(2 * lambda2 / area);
                if (mu20 != mu02 || mu11 != 0) {
                    values[12] = std::atan2(2 * mu11, mu20 - mu02) * 90 / std::numbers::pi_v<long double>;
                    values[5] = std::abs(values[12]);
                }
                const long double epsilon = std::numeric_limits<Real>::epsilon();
                values[6] = lambda1 <= epsilon ? 1 : (lambda2 <= epsilon ? 1e6L : std::min(lambda1 / lambda2, 1e6L));
                values[7] = (mu20 + mu02) / (area * area);
                values[10] = mu20 + mu02 > epsilon ? area / (2 * std::numbers::pi_v<long double> * (mu20 + mu02)) : 0;
                values[11] = lambda1 <= epsilon ? 1 : (lambda2 <= epsilon ? 0 : lambda2 / lambda1);
            }
            if (compareSupport) {
                auto support = filled;
                for (const ContourBoundary& boundary : trace.boundaries()) {
                    if (boundary.kind == ContourBoundaryKind::Internal) {
                        support.addBoundary(trace, boundary, columns, originColumn, originRow, false);
                    }
                }
                values[8] = std::clamp(1 - support.area / area, 0.L, 1.L);
                values[9] = std::hypot(centroidColumn - support.column / support.area,
                                       centroidRow - support.row / support.area) / std::sqrt(area);
            }
            writeRow(context.buffer, context.attrNames.NUM_ATTRIBUTES, node, offsets, values);
        });
    }
};
 
} // namespace mmcfilters::attributes::computers