mmcfilters
Public API documentation
Loading...
Searching...
No Matches
MomentBasedAttributeComputer.hpp
1#pragma once
2
3#include "AttributeComputerDomain.hpp"
4#include "AttributeComputerFamily.hpp"
5#include "../detail/AttributeKernelSupport.hpp"
6#include "../../trees/detail/TreeKernelAccess.hpp"
7#include "../../trees/MorphologicalTree.hpp"
8#include "../../utils/Contract.hpp"
9
10#include <array>
11#include <concepts>
12#include <cmath>
13#include <limits>
14#include <numbers>
15#include <string_view>
16#include <vector>
17
18namespace mmcfilters::attributes::computers {
19
20namespace detail {
21
23struct CentralMomentsRequest {
24 bool mu20 = false;
25 bool mu02 = false;
26 bool mu11 = false;
27 bool mu30 = false;
28 bool mu03 = false;
29 bool mu21 = false;
30 bool mu12 = false;
31 bool rawInertia = false;
32
34 [[nodiscard]] bool any() const noexcept { return mu20 || mu02 || mu11 || mu30 || mu03 || mu21 || mu12 || rawInertia; }
35
41 [[nodiscard]] static CentralMomentsRequest from(std::span<const Attribute> requestedAttributes) {
42 return {.mu20 = requestsAttribute(requestedAttributes, CentralMoment20),
43 .mu02 = requestsAttribute(requestedAttributes, CentralMoment02),
44 .mu11 = requestsAttribute(requestedAttributes, CentralMoment11),
45 .mu30 = requestsAttribute(requestedAttributes, CentralMoment30),
46 .mu03 = requestsAttribute(requestedAttributes, CentralMoment03),
47 .mu21 = requestsAttribute(requestedAttributes, CentralMoment21),
48 .mu12 = requestsAttribute(requestedAttributes, CentralMoment12),
49 .rawInertia = requestsAttribute(requestedAttributes, RawInertia)};
50 }
51};
52
54struct HuMomentsRequest {
55 bool hu1 = false;
56 bool hu2 = false;
57 bool hu3 = false;
58 bool hu4 = false;
59 bool hu5 = false;
60 bool hu6 = false;
61 bool hu7 = false;
62
64 [[nodiscard]] bool any() const noexcept { return hu1 || hu2 || hu3 || hu4 || hu5 || hu6 || hu7; }
65
71 [[nodiscard]] static HuMomentsRequest from(std::span<const Attribute> requestedAttributes) {
72 return {.hu1 = requestsAttribute(requestedAttributes, HuMoment1),
73 .hu2 = requestsAttribute(requestedAttributes, HuMoment2),
74 .hu3 = requestsAttribute(requestedAttributes, HuMoment3),
75 .hu4 = requestsAttribute(requestedAttributes, HuMoment4),
76 .hu5 = requestsAttribute(requestedAttributes, HuMoment5),
77 .hu6 = requestsAttribute(requestedAttributes, HuMoment6),
78 .hu7 = requestsAttribute(requestedAttributes, HuMoment7)};
79 }
80};
81
83struct MomentDerivedRequest {
84 bool inertia = false;
85 bool compactness = false;
86 bool eccentricity = false;
87 bool majorAxis = false;
88 bool minorAxis = false;
89 bool axisOrientation = false;
90 bool axisOrientationSigned = false;
91 bool circularity = false;
92
94 [[nodiscard]] bool any() const noexcept {
95 return inertia || compactness || eccentricity || majorAxis || minorAxis || axisOrientation || axisOrientationSigned || circularity;
96 }
97
103 [[nodiscard]] static MomentDerivedRequest from(std::span<const Attribute> requestedAttributes) {
104 return {.inertia = requestsAttribute(requestedAttributes, Inertia),
105 .compactness = requestsAttribute(requestedAttributes, Compactness),
106 .eccentricity = requestsAttribute(requestedAttributes, Eccentricity),
107 .majorAxis = requestsAttribute(requestedAttributes, LengthMajorAxis),
108 .minorAxis = requestsAttribute(requestedAttributes, LengthMinorAxis),
109 .axisOrientation = requestsAttribute(requestedAttributes, AxisOrientation),
110 .axisOrientationSigned = requestsAttribute(requestedAttributes, AxisOrientationSigned),
111 .circularity = requestsAttribute(requestedAttributes, Circularity)};
112 }
113};
114
115template <std::floating_point Real> inline constexpr Real maximumFiniteEccentricity() noexcept { return static_cast<Real>(1.0e6); }
116
117namespace kernel {
118
124template <std::floating_point Real>
125inline void computeCentralMoments(const AttributeComputeContext<Real>& context, const CentralMomentsRequest& request) {
126 if (!request.any()) {
127 return;
128 }
129
130 struct RawMoments {
131 Real m00 = Real{0};
132 Real m10 = Real{0};
133 Real m01 = Real{0};
134 Real m20 = Real{0};
135 Real m02 = Real{0};
136 Real m11 = Real{0};
137 Real m30 = Real{0};
138 Real m03 = Real{0};
139 Real m21 = Real{0};
140 Real m12 = Real{0};
141
142 void add(const RawMoments& other) {
143 m00 += other.m00;
144 m10 += other.m10;
145 m01 += other.m01;
146 m20 += other.m20;
147 m02 += other.m02;
148 m11 += other.m11;
149 m30 += other.m30;
150 m03 += other.m03;
151 m21 += other.m21;
152 m12 += other.m12;
153 }
154 };
155
156 const int stride = context.attrNames.NUM_ATTRIBUTES;
157 const auto offsetOf = [&](Attribute attribute) { return context.attrNames.indexMap.find(attribute)->second; };
158 const int mu20Offset = request.mu20 ? offsetOf(CentralMoment20) : 0;
159 const int mu02Offset = request.mu02 ? offsetOf(CentralMoment02) : 0;
160 const int mu11Offset = request.mu11 ? offsetOf(CentralMoment11) : 0;
161 const int mu30Offset = request.mu30 ? offsetOf(CentralMoment30) : 0;
162 const int mu03Offset = request.mu03 ? offsetOf(CentralMoment03) : 0;
163 const int mu21Offset = request.mu21 ? offsetOf(CentralMoment21) : 0;
164 const int mu12Offset = request.mu12 ? offsetOf(CentralMoment12) : 0;
165 const int rawInertiaOffset = request.rawInertia ? offsetOf(RawInertia) : 0;
166 const auto outputIndex = [&](NodeId node, int offset) { return static_cast<std::size_t>(node * stride + offset); };
167
168 const GridDomain2D& domain = ::mmcfilters::detail::TreeKernelAccess::gridDomain2D(context.tree);
169 std::vector<RawMoments> raw(static_cast<std::size_t>(context.tree.numInternalNodeSlots()));
170 ::mmcfilters::detail::TreeKernelAccess::traversePostOrder(
171 context.tree, context.tree.root(),
172 [&](NodeId node) {
173 RawMoments& moments = raw[static_cast<std::size_t>(node)];
174 moments = RawMoments{};
175 for (PixelId pixel : ::mmcfilters::detail::TreeKernelAccess::properParts(context.tree, node)) {
176 const Real x = static_cast<Real>(pixel % domain.columns);
177 const Real y = static_cast<Real>(pixel / domain.columns);
178 const Real x2 = x * x;
179 const Real y2 = y * y;
180 moments.m00 += Real{1};
181 moments.m10 += x;
182 moments.m01 += y;
183 moments.m20 += x2;
184 moments.m02 += y2;
185 moments.m11 += x * y;
186 moments.m30 += x2 * x;
187 moments.m03 += y2 * y;
188 moments.m21 += x2 * y;
189 moments.m12 += x * y2;
190 }
191 },
192 [&](NodeId parent, NodeId child) { raw[static_cast<std::size_t>(parent)].add(raw[static_cast<std::size_t>(child)]); },
193 [&](NodeId node) {
194 const RawMoments& moments = raw[static_cast<std::size_t>(node)];
195 if (moments.m00 <= Real{0}) {
196 return;
197 }
198 const Real cx = moments.m10 / moments.m00;
199 const Real cy = moments.m01 / moments.m00;
200 const bool needsMu20 = request.mu20 || request.rawInertia;
201 const bool needsMu02 = request.mu02 || request.rawInertia;
202 const Real mu20 = needsMu20 ? moments.m20 - Real{2} * cx * moments.m10 + cx * cx * moments.m00 : Real{0};
203 const Real mu02 = needsMu02 ? moments.m02 - Real{2} * cy * moments.m01 + cy * cy * moments.m00 : Real{0};
204 if (request.mu20)
205 context.buffer[outputIndex(node, mu20Offset)] = mu20;
206 if (request.mu02)
207 context.buffer[outputIndex(node, mu02Offset)] = mu02;
208 if (request.mu11)
209 context.buffer[outputIndex(node, mu11Offset)] = moments.m11 - cx * moments.m01 - cy * moments.m10 + cx * cy * moments.m00;
210 if (request.mu30)
211 context.buffer[outputIndex(node, mu30Offset)] =
212 moments.m30 - Real{3} * cx * moments.m20 + Real{3} * cx * cx * moments.m10 - cx * cx * cx * moments.m00;
213 if (request.mu03)
214 context.buffer[outputIndex(node, mu03Offset)] =
215 moments.m03 - Real{3} * cy * moments.m02 + Real{3} * cy * cy * moments.m01 - cy * cy * cy * moments.m00;
216 if (request.mu21)
217 context.buffer[outputIndex(node, mu21Offset)] = moments.m21 - Real{2} * cx * moments.m11 - cy * moments.m20 +
218 Real{2} * cx * cy * moments.m10 + cx * cx * moments.m01 -
219 cx * cx * cy * moments.m00;
220 if (request.mu12)
221 context.buffer[outputIndex(node, mu12Offset)] = moments.m12 - Real{2} * cy * moments.m11 - cx * moments.m02 +
222 Real{2} * cx * cy * moments.m01 + cy * cy * moments.m10 -
223 cx * cy * cy * moments.m00;
224 if (request.rawInertia)
225 context.buffer[outputIndex(node, rawInertiaOffset)] = mu20 + mu02;
226 });
227}
228
236template <std::floating_point Real>
237inline void computeHuMoments(const AttributeComputeContext<Real>& context, const HuMomentsRequest& request,
238 const DependencySourceT<Real>& momentDependency, const DependencySourceT<Real>& areaDependency) {
239 if (!request.any()) {
240 return;
241 }
242
243 const int stride = context.attrNames.NUM_ATTRIBUTES;
244 const auto offsetOf = [&](Attribute attribute) { return context.attrNames.indexMap.find(attribute)->second; };
245 const int hu1Offset = request.hu1 ? offsetOf(HuMoment1) : 0;
246 const int hu2Offset = request.hu2 ? offsetOf(HuMoment2) : 0;
247 const int hu3Offset = request.hu3 ? offsetOf(HuMoment3) : 0;
248 const int hu4Offset = request.hu4 ? offsetOf(HuMoment4) : 0;
249 const int hu5Offset = request.hu5 ? offsetOf(HuMoment5) : 0;
250 const int hu6Offset = request.hu6 ? offsetOf(HuMoment6) : 0;
251 const int hu7Offset = request.hu7 ? offsetOf(HuMoment7) : 0;
252 const auto outputIndex = [&](NodeId node, int offset) { return static_cast<std::size_t>(node * stride + offset); };
253
254 const int momentStride = momentDependency.attrNames->NUM_ATTRIBUTES;
255 const auto momentOffset = [&](Attribute attribute) { return momentDependency.attrNames->indexMap.find(attribute)->second; };
256 const auto momentIndex = [&](NodeId node, Attribute attribute) {
257 return static_cast<std::size_t>(node * momentStride + momentOffset(attribute));
258 };
259 const int areaStride = areaDependency.attrNames->NUM_ATTRIBUTES;
260 const int areaOffset = areaDependency.attrNames->indexMap.find(Area)->second;
261 const auto areaIndex = [&](NodeId node) { return static_cast<std::size_t>(node * areaStride + areaOffset); };
262 const auto normMoment = [](Real area, Real moment, int p, int q) {
263 return ::mmcfilters::attributes::numeric::safeDivide(moment, std::pow(area, static_cast<Real>(p + q + 2) / Real{2}));
264 };
265
266 ::mmcfilters::detail::TreeKernelAccess::traversePostOrder(
267 context.tree, context.tree.root(), [](NodeId) {}, [](NodeId, NodeId) {},
268 [&](NodeId node) {
269 const Real area = areaDependency.buffer[areaIndex(node)];
270 if (area <= Real{0}) {
271 return;
272 }
273 const Real eta20 = normMoment(area, momentDependency.buffer[momentIndex(node, CentralMoment20)], 2, 0);
274 const Real eta02 = normMoment(area, momentDependency.buffer[momentIndex(node, CentralMoment02)], 0, 2);
275 const Real eta11 = normMoment(area, momentDependency.buffer[momentIndex(node, CentralMoment11)], 1, 1);
276 const Real eta30 = normMoment(area, momentDependency.buffer[momentIndex(node, CentralMoment30)], 3, 0);
277 const Real eta03 = normMoment(area, momentDependency.buffer[momentIndex(node, CentralMoment03)], 0, 3);
278 const Real eta21 = normMoment(area, momentDependency.buffer[momentIndex(node, CentralMoment21)], 2, 1);
279 const Real eta12 = normMoment(area, momentDependency.buffer[momentIndex(node, CentralMoment12)], 1, 2);
280
281 if (request.hu1)
282 context.buffer[outputIndex(node, hu1Offset)] = eta20 + eta02;
283 if (request.hu2)
284 context.buffer[outputIndex(node, hu2Offset)] = std::pow(eta20 - eta02, 2) + Real{4} * std::pow(eta11, 2);
285 if (request.hu3)
286 context.buffer[outputIndex(node, hu3Offset)] = std::pow(eta30 - Real{3} * eta12, 2) + std::pow(Real{3} * eta21 - eta03, 2);
287 if (request.hu4)
288 context.buffer[outputIndex(node, hu4Offset)] = std::pow(eta30 + eta12, 2) + std::pow(eta21 + eta03, 2);
289 if (request.hu5)
290 context.buffer[outputIndex(node, hu5Offset)] =
291 (eta30 - Real{3} * eta12) * (eta30 + eta12) * (std::pow(eta30 + eta12, 2) - Real{3} * std::pow(eta21 + eta03, 2)) +
292 (Real{3} * eta21 - eta03) * (eta21 + eta03) * (Real{3} * std::pow(eta30 + eta12, 2) - std::pow(eta21 + eta03, 2));
293 if (request.hu6)
294 context.buffer[outputIndex(node, hu6Offset)] =
295 (eta20 - eta02) * (std::pow(eta30 + eta12, 2) - std::pow(eta21 + eta03, 2)) +
296 Real{4} * eta11 * (eta30 + eta12) * (eta21 + eta03);
297 if (request.hu7)
298 context.buffer[outputIndex(node, hu7Offset)] =
299 (Real{3} * eta21 - eta03) * (eta30 + eta12) * (std::pow(eta30 + eta12, 2) - Real{3} * std::pow(eta21 + eta03, 2)) -
300 (eta30 - Real{3} * eta12) * (eta21 + eta03) * (Real{3} * std::pow(eta30 + eta12, 2) - std::pow(eta21 + eta03, 2));
301 });
302}
303
311template <std::floating_point Real>
312inline void computeMomentDerived(const AttributeComputeContext<Real>& context, const MomentDerivedRequest& request,
313 const DependencySourceT<Real>& momentDependency, const DependencySourceT<Real>& areaDependency) {
314 if (!request.any()) {
315 return;
316 }
317
318 const int stride = context.attrNames.NUM_ATTRIBUTES;
319 const auto offsetOf = [&](Attribute attribute) { return context.attrNames.indexMap.find(attribute)->second; };
320 const int inertiaOffset = request.inertia ? offsetOf(Inertia) : 0;
321 const int compactnessOffset = request.compactness ? offsetOf(Compactness) : 0;
322 const int eccentricityOffset = request.eccentricity ? offsetOf(Eccentricity) : 0;
323 const int majorAxisOffset = request.majorAxis ? offsetOf(LengthMajorAxis) : 0;
324 const int minorAxisOffset = request.minorAxis ? offsetOf(LengthMinorAxis) : 0;
325 const int orientationOffset = request.axisOrientation ? offsetOf(AxisOrientation) : 0;
326 const int signedOrientationOffset = request.axisOrientationSigned ? offsetOf(AxisOrientationSigned) : 0;
327 const int circularityOffset = request.circularity ? offsetOf(Circularity) : 0;
328 const auto outputIndex = [&](NodeId node, int offset) { return static_cast<std::size_t>(node * stride + offset); };
329
330 const int momentStride = momentDependency.attrNames->NUM_ATTRIBUTES;
331 const int mu20Offset = momentDependency.attrNames->indexMap.find(CentralMoment20)->second;
332 const int mu02Offset = momentDependency.attrNames->indexMap.find(CentralMoment02)->second;
333 const int mu11Offset = momentDependency.attrNames->indexMap.find(CentralMoment11)->second;
334 const auto momentIndex = [&](NodeId node, int offset) { return static_cast<std::size_t>(node * momentStride + offset); };
335 const int areaStride = areaDependency.attrNames->NUM_ATTRIBUTES;
336 const int areaOffset = areaDependency.attrNames->indexMap.find(Area)->second;
337 const auto areaIndex = [&](NodeId node) { return static_cast<std::size_t>(node * areaStride + areaOffset); };
338
339 ::mmcfilters::detail::TreeKernelAccess::traversePostOrder(
340 context.tree, context.tree.root(), [](NodeId) {}, [](NodeId, NodeId) {},
341 [&](NodeId node) {
342 const Real mu20 = momentDependency.buffer[momentIndex(node, mu20Offset)];
343 const Real mu02 = momentDependency.buffer[momentIndex(node, mu02Offset)];
344 const Real mu11 = momentDependency.buffer[momentIndex(node, mu11Offset)];
345 const Real area = areaDependency.buffer[areaIndex(node)];
346 const Real discriminant = std::pow(mu20 - mu02, 2) + Real{4} * std::pow(mu11, 2);
347 const Real sqrtDiscriminant = ::mmcfilters::attributes::numeric::safeSqrt(discriminant);
348 const Real lambda1 = mu20 + mu02 + sqrtDiscriminant;
349 const Real lambda2 = mu20 + mu02 - sqrtDiscriminant;
350
351 if (request.majorAxis)
352 context.buffer[outputIndex(node, majorAxisOffset)] =
353 ::mmcfilters::attributes::numeric::safeSqrt(::mmcfilters::attributes::numeric::safeDivide(Real{2} * lambda1, area));
354 if (request.minorAxis)
355 context.buffer[outputIndex(node, minorAxisOffset)] =
356 ::mmcfilters::attributes::numeric::safeSqrt(::mmcfilters::attributes::numeric::safeDivide(Real{2} * lambda2, area));
357 if (request.eccentricity) {
358 const Real eps = std::numeric_limits<Real>::epsilon();
359 if (lambda1 <= eps && std::abs(lambda2) <= eps) {
360 context.buffer[outputIndex(node, eccentricityOffset)] = Real{1};
361 } else if (lambda2 <= eps) {
362 context.buffer[outputIndex(node, eccentricityOffset)] = maximumFiniteEccentricity<Real>();
363 } else {
364 context.buffer[outputIndex(node, eccentricityOffset)] = ::mmcfilters::attributes::numeric::clampUpper(
365 ::mmcfilters::attributes::numeric::safeDivide(lambda1, lambda2, maximumFiniteEccentricity<Real>()),
366 maximumFiniteEccentricity<Real>());
367 }
368 }
369 if (request.compactness)
370 context.buffer[outputIndex(node, compactnessOffset)] = (Real{1} / (Real{2} * std::numbers::pi_v<Real>)) *
371 ::mmcfilters::attributes::numeric::safeDivide(area, mu20 + mu02);
372 if (request.axisOrientation || request.axisOrientationSigned) {
373 const Real degrees = (mu20 != mu02 || mu11 != Real{0})
374 ? Real{0.5} * std::atan2(Real{2} * mu11, mu20 - mu02) * (Real{180} / std::numbers::pi_v<Real>)
375 : Real{0};
376 if (request.axisOrientation)
377 context.buffer[outputIndex(node, orientationOffset)] = std::fmod(std::abs(degrees), Real{360});
378 if (request.axisOrientationSigned)
379 context.buffer[outputIndex(node, signedOrientationOffset)] = degrees;
380 }
381 if (request.inertia) {
382 const Real areaSquared = area * area;
383 context.buffer[outputIndex(node, inertiaOffset)] = ::mmcfilters::attributes::numeric::safeDivide(mu20, areaSquared) +
384 ::mmcfilters::attributes::numeric::safeDivide(mu02, areaSquared);
385 }
386 if (request.circularity) {
387 const Real eps = std::numeric_limits<Real>::epsilon();
388 if (lambda1 <= eps && std::abs(lambda2) <= eps) {
389 context.buffer[outputIndex(node, circularityOffset)] = Real{1};
390 } else if (lambda1 <= eps || lambda2 <= eps) {
391 context.buffer[outputIndex(node, circularityOffset)] = Real{0};
392 } else {
393 context.buffer[outputIndex(node, circularityOffset)] = ::mmcfilters::attributes::numeric::safeDivide(lambda2, lambda1);
394 }
395 }
396 });
397}
398
399} // namespace kernel
400
401template <std::floating_point Real> inline void validateCentralMomentsContext(const AttributeComputeContext<Real>& context) {
402 requireAttributeBufferShape(context.tree, context.buffer, context.attrNames);
403 requireRequestedAttributeColumns(context);
404 static_cast<void>(context.tree.requireGridDomain2D("CentralMomentsComputer"));
405}
406
407template <std::floating_point Real> inline void validateMomentDependencyContext(const AttributeComputeContext<Real>& context) {
408 requireAttributeBufferShape(context.tree, context.buffer, context.attrNames);
409 requireRequestedAttributeColumns(context);
410}
411} // namespace detail
412
434class CentralMomentsComputer {
435 public:
437 static constexpr std::string_view familyName = "central-moments";
438
440 static constexpr AttributeComputerFamily family = AttributeComputerFamily::CentralMoments;
441
443 static constexpr AttributeComputerDomain domain = AttributeComputerDomain::Topology;
444
448 inline static constexpr std::array<Attribute, 8> producedAttributes{CentralMoment20, CentralMoment02, CentralMoment11, CentralMoment30,
449 CentralMoment03, CentralMoment21, CentralMoment12, RawInertia};
450
463 template <std::floating_point Real> static void compute(const AttributeComputeContext<Real>& context) {
464 const detail::CentralMomentsRequest request = detail::CentralMomentsRequest::from(context.requestedAttributes);
465 MMCFILTERS_CONTRACT_CHECKED_ONLY(detail::validateCentralMomentsContext(context));
466 detail::kernel::computeCentralMoments(context, request);
467 }
468
469 public:
477 template <std::floating_point Real> static void computeUnitRows(const UnitAttributeComputeContext<Real>& context) {
478 requireUnitAttributeBufferShape(context.tree, context.unitPixels, context.buffer, context.attrNames);
479 computeUnitRowsEstablished(context);
480 }
481
482 private:
489 template <class Computer, class Context> friend void mmcfilters::detail::kernel::computeUnitRows(const Context& context);
490
492 template <std::floating_point Real> static void computeUnitRowsEstablished(const UnitAttributeComputeContext<Real>& context) {
493 constexpr std::array<Attribute, 8> zeroAttributes{CentralMoment20, CentralMoment02, CentralMoment11, CentralMoment30,
494 CentralMoment03, CentralMoment21, CentralMoment12, RawInertia};
495 for (const Attribute attribute : zeroAttributes) {
496 if (!requestsAttribute(context.requestedAttributes, attribute)) {
497 continue;
498 }
499 for (NodeId leafIndex = 0; leafIndex < static_cast<NodeId>(context.unitPixels.size()); ++leafIndex) {
500 context.buffer[context.attrNames.linearIndex(leafIndex, attribute)] = Real{0};
501 }
502 }
503 }
504};
505
519class HuMomentsComputer {
520 public:
522 static constexpr std::string_view familyName = "hu-moments";
523
525 static constexpr AttributeComputerFamily family = AttributeComputerFamily::HuMoments;
526
528 static constexpr AttributeComputerDomain domain = AttributeComputerDomain::Topology;
529
533 inline static constexpr std::array<Attribute, 7> producedAttributes{HuMoment1, HuMoment2, HuMoment3, HuMoment4,
534 HuMoment5, HuMoment6, HuMoment7};
535
547 template <std::floating_point Real> static void compute(const AttributeComputeContext<Real>& context) {
548 const detail::HuMomentsRequest request = detail::HuMomentsRequest::from(context.requestedAttributes);
549 MMCFILTERS_CONTRACT_CHECKED_ONLY(detail::validateMomentDependencyContext(context));
550 if (!request.any()) {
551 return;
552 }
553
554 const DependencySourceT<Real>* momentDependency = nullptr;
555 const DependencySourceT<Real>* areaDependency = nullptr;
556 if constexpr (contract::validationsEnabled) {
557 momentDependency = &context.dependencies.requireAll({CentralMoment20, CentralMoment02, CentralMoment11, CentralMoment30,
558 CentralMoment03, CentralMoment21, CentralMoment12});
559 areaDependency = &context.dependencies.require(Area);
560 } else {
561 momentDependency = ::mmcfilters::findDependencySource(context.dependencySources, CentralMoment20);
562 areaDependency = ::mmcfilters::findDependencySource(context.dependencySources, Area);
563 }
564 detail::kernel::computeHuMoments(context, request, *momentDependency, *areaDependency);
565 }
566
567 public:
576 template <std::floating_point Real> static void computeUnitRows(const UnitAttributeComputeContext<Real>& context) {
577 requireUnitAttributeBufferShape(context.tree, context.unitPixels, context.buffer, context.attrNames);
578 computeUnitRowsEstablished(context);
579 }
580
581 private:
588 template <class Computer, class Context> friend void mmcfilters::detail::kernel::computeUnitRows(const Context& context);
589
591 template <std::floating_point Real> static void computeUnitRowsEstablished(const UnitAttributeComputeContext<Real>& context) {
592 constexpr std::array<Attribute, 7> zeroAttributes{HuMoment1, HuMoment2, HuMoment3, HuMoment4, HuMoment5, HuMoment6, HuMoment7};
593 for (const Attribute attribute : zeroAttributes) {
594 if (!requestsAttribute(context.requestedAttributes, attribute)) {
595 continue;
596 }
597 for (NodeId leafIndex = 0; leafIndex < static_cast<NodeId>(context.unitPixels.size()); ++leafIndex) {
598 context.buffer[context.attrNames.linearIndex(leafIndex, attribute)] = Real{0};
599 }
600 }
601 }
602};
603
625class MomentBasedAttributeComputer {
626 public:
628 static constexpr std::string_view familyName = "moment-derived";
629
631 static constexpr AttributeComputerFamily family = AttributeComputerFamily::MomentDerived;
632
634 static constexpr AttributeComputerDomain domain = AttributeComputerDomain::Topology;
635
641 template <std::floating_point Real> static constexpr Real maxFiniteEccentricity() noexcept {
642 return detail::maximumFiniteEccentricity<Real>();
643 }
644
648 inline static constexpr std::array<Attribute, 8> producedAttributes{Inertia, Compactness, Eccentricity, LengthMajorAxis,
649 LengthMinorAxis, AxisOrientation, Circularity, AxisOrientationSigned};
650
662 template <std::floating_point Real> static void compute(const AttributeComputeContext<Real>& context) {
663 const detail::MomentDerivedRequest request = detail::MomentDerivedRequest::from(context.requestedAttributes);
664 MMCFILTERS_CONTRACT_CHECKED_ONLY(detail::validateMomentDependencyContext(context));
665 if (!request.any()) {
666 return;
667 }
668
669 const DependencySourceT<Real>* momentDependency = nullptr;
670 const DependencySourceT<Real>* areaDependency = nullptr;
671 if constexpr (contract::validationsEnabled) {
672 momentDependency = &context.dependencies.requireAll({CentralMoment20, CentralMoment02, CentralMoment11});
673 areaDependency = &context.dependencies.require(Area);
674 } else {
675 momentDependency = ::mmcfilters::findDependencySource(context.dependencySources, CentralMoment20);
676 areaDependency = ::mmcfilters::findDependencySource(context.dependencySources, Area);
677 }
678 detail::kernel::computeMomentDerived(context, request, *momentDependency, *areaDependency);
679 }
680
681 public:
690 template <std::floating_point Real> static void computeUnitRows(const UnitAttributeComputeContext<Real>& context) {
691 requireUnitAttributeBufferShape(context.tree, context.unitPixels, context.buffer, context.attrNames);
692 computeUnitRowsEstablished(context);
693 }
694
695 private:
702 template <class Computer, class Context> friend void mmcfilters::detail::kernel::computeUnitRows(const Context& context);
703
705 template <std::floating_point Real> static void computeUnitRowsEstablished(const UnitAttributeComputeContext<Real>& context) {
706 constexpr std::array<Attribute, 8> zeroAttributes{Compactness, Eccentricity, LengthMajorAxis, LengthMinorAxis,
707 AxisOrientation, Inertia, Circularity, AxisOrientationSigned};
708 for (const Attribute attribute : zeroAttributes) {
709 if (!requestsAttribute(context.requestedAttributes, attribute)) {
710 continue;
711 }
712 for (NodeId leafIndex = 0; leafIndex < static_cast<NodeId>(context.unitPixels.size()); ++leafIndex) {
713 context.buffer[context.attrNames.linearIndex(leafIndex, attribute)] =
714 (attribute == Eccentricity || attribute == Circularity) ? Real{1} : Real{0};
715 }
716 }
717 }
718};
719
720} // namespace mmcfilters::attributes::computers
mmcfilters::NodeId
int NodeId
Node identifier type used throughout the project.
Definition Common.hpp:17
MMCFILTERS_CONTRACT_CHECKED_ONLY
#define MMCFILTERS_CONTRACT_CHECKED_ONLY(...)
Executes validation statements only when defensive checks are enabled.
Definition Contract.hpp:67
mmcfilters::attributes::computers::CentralMomentsComputer
Computes geometric central moments up to third order.
Definition MomentBasedAttributeComputer.hpp:434
mmcfilters::attributes::computers::CentralMomentsComputer::compute
static void compute(const AttributeComputeContext< Real > &context)
Computes the requested central moments.
Definition MomentBasedAttributeComputer.hpp:463
mmcfilters::attributes::computers::CentralMomentsComputer::computeUnitRows
static void computeUnitRows(const UnitAttributeComputeContext< Real > &context)
Materializes central moments for one-pixel unit supports.
Definition MomentBasedAttributeComputer.hpp:477
mmcfilters::attributes::computers::HuMomentsComputer::compute
static void compute(const AttributeComputeContext< Real > &context)
Computes the requested Hu invariant moments.
Definition MomentBasedAttributeComputer.hpp:547
mmcfilters::attributes::computers::HuMomentsComputer::computeUnitRows
static void computeUnitRows(const UnitAttributeComputeContext< Real > &context)
Materializes Hu moments for one-pixel unit supports.
Definition MomentBasedAttributeComputer.hpp:576
mmcfilters::attributes::computers::MomentBasedAttributeComputer
Computes higher-level shape descriptors derived from second-order central moments.
Definition MomentBasedAttributeComputer.hpp:625
mmcfilters::attributes::computers::MomentBasedAttributeComputer::computeUnitRows
static void computeUnitRows(const UnitAttributeComputeContext< Real > &context)
Materializes moment-derived descriptors for one-pixel unit supports.
Definition MomentBasedAttributeComputer.hpp:690
mmcfilters::attributes::computers::MomentBasedAttributeComputer::maxFiniteEccentricity
static constexpr Real maxFiniteEccentricity() noexcept
Largest finite eccentricity emitted for degenerate one-dimensional supports.
Definition MomentBasedAttributeComputer.hpp:641
mmcfilters::attributes::computers::MomentBasedAttributeComputer::compute
static void compute(const AttributeComputeContext< Real > &context)
Computes the requested moment-derived descriptors.
Definition MomentBasedAttributeComputer.hpp:662
mmcfilters::ComputedAttributeData
Owning result for one computed scalar attribute layout and buffer.
Definition AttributeResultTypes.hpp:32
mmcfilters::ComputedAttributeData::second
std::vector< Real > second
Flat per-node attribute buffer indexed through first.
Definition AttributeResultTypes.hpp:37
Generated by doxygen 1.9.8