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/CommittedTreeAccess.hpp"
7#include "../../trees/detail/TreeTraversalDetail.hpp"
8#include "../../trees/MorphologicalTree.hpp"
9#include "../../utils/Contract.hpp"
10
11#include <array>
12#include <concepts>
13#include <cmath>
14#include <limits>
15#include <numbers>
16#include <string_view>
17#include <vector>
18
19namespace mmcfilters::attributes::computers {
20
21namespace detail {
22
24struct CentralMomentsRequest {
25 bool mu20 = false;
26 bool mu02 = false;
27 bool mu11 = false;
28 bool mu30 = false;
29 bool mu03 = false;
30 bool mu21 = false;
31 bool mu12 = false;
32
34 [[nodiscard]] bool any() const noexcept { return mu20 || mu02 || mu11 || mu30 || mu03 || mu21 || mu12; }
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 }
50};
51
53struct HuMomentsRequest {
54 bool hu1 = false;
55 bool hu2 = false;
56 bool hu3 = false;
57 bool hu4 = false;
58 bool hu5 = false;
59 bool hu6 = false;
60 bool hu7 = false;
61
63 [[nodiscard]] bool any() const noexcept { return hu1 || hu2 || hu3 || hu4 || hu5 || hu6 || hu7; }
64
70 [[nodiscard]] static HuMomentsRequest from(std::span<const Attribute> requestedAttributes) {
71 return {.hu1 = requestsAttribute(requestedAttributes, HuMoment1),
72 .hu2 = requestsAttribute(requestedAttributes, HuMoment2),
73 .hu3 = requestsAttribute(requestedAttributes, HuMoment3),
74 .hu4 = requestsAttribute(requestedAttributes, HuMoment4),
75 .hu5 = requestsAttribute(requestedAttributes, HuMoment5),
76 .hu6 = requestsAttribute(requestedAttributes, HuMoment6),
77 .hu7 = requestsAttribute(requestedAttributes, HuMoment7)};
78 }
79};
80
82struct MomentDerivedRequest {
83 bool inertia = false;
84 bool compactness = false;
85 bool eccentricity = false;
86 bool majorAxis = false;
87 bool minorAxis = false;
88 bool axisOrientation = false;
89 bool circularity = false;
90
92 [[nodiscard]] bool any() const noexcept {
93 return inertia || compactness || eccentricity || majorAxis || minorAxis || axisOrientation || circularity;
94 }
95
101 [[nodiscard]] static MomentDerivedRequest from(std::span<const Attribute> requestedAttributes) {
102 return {.inertia = requestsAttribute(requestedAttributes, Inertia),
103 .compactness = requestsAttribute(requestedAttributes, Compactness),
104 .eccentricity = requestsAttribute(requestedAttributes, Eccentricity),
105 .majorAxis = requestsAttribute(requestedAttributes, LengthMajorAxis),
106 .minorAxis = requestsAttribute(requestedAttributes, LengthMinorAxis),
107 .axisOrientation = requestsAttribute(requestedAttributes, AxisOrientation),
108 .circularity = requestsAttribute(requestedAttributes, Circularity)};
109 }
110};
111
112template <std::floating_point Real> inline constexpr Real maximumFiniteEccentricity() noexcept { return static_cast<Real>(1.0e6); }
113
114namespace kernel {
115
121template <std::floating_point Real>
122inline void computeCentralMoments(const AttributeComputeContext<Real>& context, const CentralMomentsRequest& request) {
123 if (!request.any()) {
124 return;
125 }
126
127 struct RawMoments {
128 Real m00 = Real{0};
129 Real m10 = Real{0};
130 Real m01 = Real{0};
131 Real m20 = Real{0};
132 Real m02 = Real{0};
133 Real m11 = Real{0};
134 Real m30 = Real{0};
135 Real m03 = Real{0};
136 Real m21 = Real{0};
137 Real m12 = Real{0};
138
139 void add(const RawMoments& other) {
140 m00 += other.m00;
141 m10 += other.m10;
142 m01 += other.m01;
143 m20 += other.m20;
144 m02 += other.m02;
145 m11 += other.m11;
146 m30 += other.m30;
147 m03 += other.m03;
148 m21 += other.m21;
149 m12 += other.m12;
150 }
151 };
152
153 const int stride = context.attrNames.NUM_ATTRIBUTES;
154 const auto offsetOf = [&](Attribute attribute) { return context.attrNames.indexMap.find(attribute)->second; };
155 const int mu20Offset = request.mu20 ? offsetOf(CentralMoment20) : 0;
156 const int mu02Offset = request.mu02 ? offsetOf(CentralMoment02) : 0;
157 const int mu11Offset = request.mu11 ? offsetOf(CentralMoment11) : 0;
158 const int mu30Offset = request.mu30 ? offsetOf(CentralMoment30) : 0;
159 const int mu03Offset = request.mu03 ? offsetOf(CentralMoment03) : 0;
160 const int mu21Offset = request.mu21 ? offsetOf(CentralMoment21) : 0;
161 const int mu12Offset = request.mu12 ? offsetOf(CentralMoment12) : 0;
162 const auto outputIndex = [&](NodeId node, int offset) { return static_cast<std::size_t>(node * stride + offset); };
163
164 const GridDomain2D& domain = ::mmcfilters::detail::CommittedTreeAccess::gridDomain2D(context.tree);
165 std::vector<RawMoments> raw(static_cast<std::size_t>(context.tree.numInternalNodeSlots()));
166 ::mmcfilters::detail::kernel::traversePostOrder(
167 context.tree, context.tree.root(),
168 [&](NodeId node) {
169 RawMoments& moments = raw[static_cast<std::size_t>(node)];
170 moments = RawMoments{};
171 for (PixelId pixel : ::mmcfilters::detail::CommittedTreeAccess::properParts(context.tree, node)) {
172 const Real x = static_cast<Real>(pixel % domain.columns);
173 const Real y = static_cast<Real>(pixel / domain.columns);
174 const Real x2 = x * x;
175 const Real y2 = y * y;
176 moments.m00 += Real{1};
177 moments.m10 += x;
178 moments.m01 += y;
179 moments.m20 += x2;
180 moments.m02 += y2;
181 moments.m11 += x * y;
182 moments.m30 += x2 * x;
183 moments.m03 += y2 * y;
184 moments.m21 += x2 * y;
185 moments.m12 += x * y2;
186 }
187 },
188 [&](NodeId parent, NodeId child) { raw[static_cast<std::size_t>(parent)].add(raw[static_cast<std::size_t>(child)]); },
189 [&](NodeId node) {
190 const RawMoments& moments = raw[static_cast<std::size_t>(node)];
191 if (moments.m00 <= Real{0}) {
192 return;
193 }
194 const Real cx = moments.m10 / moments.m00;
195 const Real cy = moments.m01 / moments.m00;
196 if (request.mu20)
197 context.buffer[outputIndex(node, mu20Offset)] = moments.m20 - Real{2} * cx * moments.m10 + cx * cx * moments.m00;
198 if (request.mu02)
199 context.buffer[outputIndex(node, mu02Offset)] = moments.m02 - Real{2} * cy * moments.m01 + cy * cy * moments.m00;
200 if (request.mu11)
201 context.buffer[outputIndex(node, mu11Offset)] = moments.m11 - cx * moments.m01 - cy * moments.m10 + cx * cy * moments.m00;
202 if (request.mu30)
203 context.buffer[outputIndex(node, mu30Offset)] =
204 moments.m30 - Real{3} * cx * moments.m20 + Real{3} * cx * cx * moments.m10 - cx * cx * cx * moments.m00;
205 if (request.mu03)
206 context.buffer[outputIndex(node, mu03Offset)] =
207 moments.m03 - Real{3} * cy * moments.m02 + Real{3} * cy * cy * moments.m01 - cy * cy * cy * moments.m00;
208 if (request.mu21)
209 context.buffer[outputIndex(node, mu21Offset)] = moments.m21 - Real{2} * cx * moments.m11 - cy * moments.m20 +
210 Real{2} * cx * cy * moments.m10 + cx * cx * moments.m01 -
211 cx * cx * cy * moments.m00;
212 if (request.mu12)
213 context.buffer[outputIndex(node, mu12Offset)] = moments.m12 - Real{2} * cy * moments.m11 - cx * moments.m02 +
214 Real{2} * cx * cy * moments.m01 + cy * cy * moments.m10 -
215 cx * cy * cy * moments.m00;
216 });
217}
218
226template <std::floating_point Real>
227inline void computeHuMoments(const AttributeComputeContext<Real>& context, const HuMomentsRequest& request,
228 const DependencySourceT<Real>& momentDependency, const DependencySourceT<Real>& areaDependency) {
229 if (!request.any()) {
230 return;
231 }
232
233 const int stride = context.attrNames.NUM_ATTRIBUTES;
234 const auto offsetOf = [&](Attribute attribute) { return context.attrNames.indexMap.find(attribute)->second; };
235 const int hu1Offset = request.hu1 ? offsetOf(HuMoment1) : 0;
236 const int hu2Offset = request.hu2 ? offsetOf(HuMoment2) : 0;
237 const int hu3Offset = request.hu3 ? offsetOf(HuMoment3) : 0;
238 const int hu4Offset = request.hu4 ? offsetOf(HuMoment4) : 0;
239 const int hu5Offset = request.hu5 ? offsetOf(HuMoment5) : 0;
240 const int hu6Offset = request.hu6 ? offsetOf(HuMoment6) : 0;
241 const int hu7Offset = request.hu7 ? offsetOf(HuMoment7) : 0;
242 const auto outputIndex = [&](NodeId node, int offset) { return static_cast<std::size_t>(node * stride + offset); };
243
244 const int momentStride = momentDependency.attrNames->NUM_ATTRIBUTES;
245 const auto momentOffset = [&](Attribute attribute) { return momentDependency.attrNames->indexMap.find(attribute)->second; };
246 const auto momentIndex = [&](NodeId node, Attribute attribute) {
247 return static_cast<std::size_t>(node * momentStride + momentOffset(attribute));
248 };
249 const int areaStride = areaDependency.attrNames->NUM_ATTRIBUTES;
250 const int areaOffset = areaDependency.attrNames->indexMap.find(Area)->second;
251 const auto areaIndex = [&](NodeId node) { return static_cast<std::size_t>(node * areaStride + areaOffset); };
252 const auto normMoment = [](Real area, Real moment, int p, int q) {
253 return ::mmcfilters::attributes::numeric::safeDivide(moment, std::pow(area, static_cast<Real>(p + q + 2) / Real{2}));
254 };
255
256 ::mmcfilters::detail::kernel::traversePostOrder(
257 context.tree, context.tree.root(), [](NodeId) {}, [](NodeId, NodeId) {},
258 [&](NodeId node) {
259 const Real area = areaDependency.buffer[areaIndex(node)];
260 if (area <= Real{0}) {
261 return;
262 }
263 const Real eta20 = normMoment(area, momentDependency.buffer[momentIndex(node, CentralMoment20)], 2, 0);
264 const Real eta02 = normMoment(area, momentDependency.buffer[momentIndex(node, CentralMoment02)], 0, 2);
265 const Real eta11 = normMoment(area, momentDependency.buffer[momentIndex(node, CentralMoment11)], 1, 1);
266 const Real eta30 = normMoment(area, momentDependency.buffer[momentIndex(node, CentralMoment30)], 3, 0);
267 const Real eta03 = normMoment(area, momentDependency.buffer[momentIndex(node, CentralMoment03)], 0, 3);
268 const Real eta21 = normMoment(area, momentDependency.buffer[momentIndex(node, CentralMoment21)], 2, 1);
269 const Real eta12 = normMoment(area, momentDependency.buffer[momentIndex(node, CentralMoment12)], 1, 2);
270
271 if (request.hu1)
272 context.buffer[outputIndex(node, hu1Offset)] = eta20 + eta02;
273 if (request.hu2)
274 context.buffer[outputIndex(node, hu2Offset)] = std::pow(eta20 - eta02, 2) + Real{4} * std::pow(eta11, 2);
275 if (request.hu3)
276 context.buffer[outputIndex(node, hu3Offset)] = std::pow(eta30 - Real{3} * eta12, 2) + std::pow(Real{3} * eta21 - eta03, 2);
277 if (request.hu4)
278 context.buffer[outputIndex(node, hu4Offset)] = std::pow(eta30 + eta12, 2) + std::pow(eta21 + eta03, 2);
279 if (request.hu5)
280 context.buffer[outputIndex(node, hu5Offset)] =
281 (eta30 - Real{3} * eta12) * (eta30 + eta12) * (std::pow(eta30 + eta12, 2) - Real{3} * std::pow(eta21 + eta03, 2)) +
282 (Real{3} * eta21 - eta03) * (eta21 + eta03) * (Real{3} * std::pow(eta30 + eta12, 2) - std::pow(eta21 + eta03, 2));
283 if (request.hu6)
284 context.buffer[outputIndex(node, hu6Offset)] =
285 (eta20 - eta02) * (std::pow(eta30 + eta12, 2) - std::pow(eta21 + eta03, 2)) +
286 Real{4} * eta11 * (eta30 + eta12) * (eta21 + eta03);
287 if (request.hu7)
288 context.buffer[outputIndex(node, hu7Offset)] =
289 (Real{3} * eta21 - eta03) * (eta30 + eta12) * (std::pow(eta30 + eta12, 2) - Real{3} * std::pow(eta21 + eta03, 2)) -
290 (eta30 - Real{3} * eta12) * (eta21 + eta03) * (Real{3} * std::pow(eta30 + eta12, 2) - std::pow(eta21 + eta03, 2));
291 });
292}
293
301template <std::floating_point Real>
302inline void computeMomentDerived(const AttributeComputeContext<Real>& context, const MomentDerivedRequest& request,
303 const DependencySourceT<Real>& momentDependency, const DependencySourceT<Real>& areaDependency) {
304 if (!request.any()) {
305 return;
306 }
307
308 const int stride = context.attrNames.NUM_ATTRIBUTES;
309 const auto offsetOf = [&](Attribute attribute) { return context.attrNames.indexMap.find(attribute)->second; };
310 const int inertiaOffset = request.inertia ? offsetOf(Inertia) : 0;
311 const int compactnessOffset = request.compactness ? offsetOf(Compactness) : 0;
312 const int eccentricityOffset = request.eccentricity ? offsetOf(Eccentricity) : 0;
313 const int majorAxisOffset = request.majorAxis ? offsetOf(LengthMajorAxis) : 0;
314 const int minorAxisOffset = request.minorAxis ? offsetOf(LengthMinorAxis) : 0;
315 const int orientationOffset = request.axisOrientation ? offsetOf(AxisOrientation) : 0;
316 const int circularityOffset = request.circularity ? offsetOf(Circularity) : 0;
317 const auto outputIndex = [&](NodeId node, int offset) { return static_cast<std::size_t>(node * stride + offset); };
318
319 const int momentStride = momentDependency.attrNames->NUM_ATTRIBUTES;
320 const int mu20Offset = momentDependency.attrNames->indexMap.find(CentralMoment20)->second;
321 const int mu02Offset = momentDependency.attrNames->indexMap.find(CentralMoment02)->second;
322 const int mu11Offset = momentDependency.attrNames->indexMap.find(CentralMoment11)->second;
323 const auto momentIndex = [&](NodeId node, int offset) { return static_cast<std::size_t>(node * momentStride + offset); };
324 const int areaStride = areaDependency.attrNames->NUM_ATTRIBUTES;
325 const int areaOffset = areaDependency.attrNames->indexMap.find(Area)->second;
326 const auto areaIndex = [&](NodeId node) { return static_cast<std::size_t>(node * areaStride + areaOffset); };
327
328 ::mmcfilters::detail::kernel::traversePostOrder(
329 context.tree, context.tree.root(), [](NodeId) {}, [](NodeId, NodeId) {},
330 [&](NodeId node) {
331 const Real mu20 = momentDependency.buffer[momentIndex(node, mu20Offset)];
332 const Real mu02 = momentDependency.buffer[momentIndex(node, mu02Offset)];
333 const Real mu11 = momentDependency.buffer[momentIndex(node, mu11Offset)];
334 const Real area = areaDependency.buffer[areaIndex(node)];
335 const Real discriminant = std::pow(mu20 - mu02, 2) + Real{4} * std::pow(mu11, 2);
336 const Real sqrtDiscriminant = ::mmcfilters::attributes::numeric::safeSqrt(discriminant);
337 const Real lambda1 = mu20 + mu02 + sqrtDiscriminant;
338 const Real lambda2 = mu20 + mu02 - sqrtDiscriminant;
339
340 if (request.majorAxis)
341 context.buffer[outputIndex(node, majorAxisOffset)] =
342 ::mmcfilters::attributes::numeric::safeSqrt(::mmcfilters::attributes::numeric::safeDivide(Real{2} * lambda1, area));
343 if (request.minorAxis)
344 context.buffer[outputIndex(node, minorAxisOffset)] =
345 ::mmcfilters::attributes::numeric::safeSqrt(::mmcfilters::attributes::numeric::safeDivide(Real{2} * lambda2, area));
346 if (request.eccentricity) {
347 const Real eps = std::numeric_limits<Real>::epsilon();
348 if (lambda1 <= eps && std::abs(lambda2) <= eps) {
349 context.buffer[outputIndex(node, eccentricityOffset)] = Real{1};
350 } else if (lambda2 <= eps) {
351 context.buffer[outputIndex(node, eccentricityOffset)] = maximumFiniteEccentricity<Real>();
352 } else {
353 context.buffer[outputIndex(node, eccentricityOffset)] = ::mmcfilters::attributes::numeric::clampUpper(
354 ::mmcfilters::attributes::numeric::safeDivide(lambda1, lambda2, maximumFiniteEccentricity<Real>()),
355 maximumFiniteEccentricity<Real>());
356 }
357 }
358 if (request.compactness)
359 context.buffer[outputIndex(node, compactnessOffset)] = (Real{1} / (Real{2} * std::numbers::pi_v<Real>)) *
360 ::mmcfilters::attributes::numeric::safeDivide(area, mu20 + mu02);
361 if (request.axisOrientation) {
362 if (mu20 != mu02 || mu11 != Real{0}) {
363 const Real radians = Real{0.5} * std::atan2(Real{2} * mu11, mu20 - mu02);
364 const Real degrees = radians * (Real{180} / std::numbers::pi_v<Real>);
365 context.buffer[outputIndex(node, orientationOffset)] = std::fmod(std::abs(degrees), Real{360});
366 } else {
367 context.buffer[outputIndex(node, orientationOffset)] = Real{0};
368 }
369 }
370 if (request.inertia) {
371 const Real areaSquared = area * area;
372 context.buffer[outputIndex(node, inertiaOffset)] = ::mmcfilters::attributes::numeric::safeDivide(mu20, areaSquared) +
373 ::mmcfilters::attributes::numeric::safeDivide(mu02, areaSquared);
374 }
375 if (request.circularity) {
376 const Real eps = std::numeric_limits<Real>::epsilon();
377 if (lambda1 <= eps && std::abs(lambda2) <= eps) {
378 context.buffer[outputIndex(node, circularityOffset)] = Real{1};
379 } else if (lambda1 <= eps || lambda2 <= eps) {
380 context.buffer[outputIndex(node, circularityOffset)] = Real{0};
381 } else {
382 context.buffer[outputIndex(node, circularityOffset)] = ::mmcfilters::attributes::numeric::safeDivide(lambda2, lambda1);
383 }
384 }
385 });
386}
387
388} // namespace kernel
389
390template <std::floating_point Real> inline void validateCentralMomentsContext(const AttributeComputeContext<Real>& context) {
391 requireAttributeBufferShape(context.tree, context.buffer, context.attrNames);
392 requireRequestedAttributeColumns(context);
393 static_cast<void>(context.tree.requireGridDomain2D("CentralMomentsComputer"));
394}
395
396template <std::floating_point Real> inline void validateMomentDependencyContext(const AttributeComputeContext<Real>& context) {
397 requireAttributeBufferShape(context.tree, context.buffer, context.attrNames);
398 requireRequestedAttributeColumns(context);
399}
400} // namespace detail
401
419class CentralMomentsComputer {
420 public:
422 static constexpr std::string_view familyName = "central-moments";
423
425 static constexpr AttributeComputerFamily family = AttributeComputerFamily::CentralMoments;
426
428 static constexpr AttributeComputerDomain domain = AttributeComputerDomain::Topology;
429
433 inline static constexpr std::array<Attribute, 7> producedAttributes{CentralMoment20, CentralMoment02, CentralMoment11, CentralMoment30,
434 CentralMoment03, CentralMoment21, CentralMoment12};
435
448 template <std::floating_point Real> static void compute(const AttributeComputeContext<Real>& context) {
449 const detail::CentralMomentsRequest request = detail::CentralMomentsRequest::from(context.requestedAttributes);
450 MMCFILTERS_CONTRACT_CHECKED_ONLY(detail::validateCentralMomentsContext(context));
451 detail::kernel::computeCentralMoments(context, request);
452 }
453
454 public:
462 template <std::floating_point Real> static void computeUnitRows(const UnitAttributeComputeContext<Real>& context) {
463 requireUnitAttributeBufferShape(context.tree, context.unitPixels, context.buffer, context.attrNames);
464 constexpr std::array<Attribute, 7> zeroAttributes{CentralMoment20, CentralMoment02, CentralMoment11, CentralMoment30,
465 CentralMoment03, CentralMoment21, CentralMoment12};
466 for (const Attribute attribute : zeroAttributes) {
467 if (!requestsAttribute(context.requestedAttributes, attribute)) {
468 continue;
469 }
470 for (NodeId leafIndex = 0; leafIndex < static_cast<NodeId>(context.unitPixels.size()); ++leafIndex) {
471 context.buffer[context.attrNames.linearIndex(leafIndex, attribute)] = Real{0};
472 }
473 }
474 }
475};
476
490class HuMomentsComputer {
491 public:
493 static constexpr std::string_view familyName = "hu-moments";
494
496 static constexpr AttributeComputerFamily family = AttributeComputerFamily::HuMoments;
497
499 static constexpr AttributeComputerDomain domain = AttributeComputerDomain::Topology;
500
504 inline static constexpr std::array<Attribute, 7> producedAttributes{HuMoment1, HuMoment2, HuMoment3, HuMoment4,
505 HuMoment5, HuMoment6, HuMoment7};
506
518 template <std::floating_point Real> static void compute(const AttributeComputeContext<Real>& context) {
519 const detail::HuMomentsRequest request = detail::HuMomentsRequest::from(context.requestedAttributes);
520 MMCFILTERS_CONTRACT_CHECKED_ONLY(detail::validateMomentDependencyContext(context));
521 if (!request.any()) {
522 return;
523 }
524
525 const DependencySourceT<Real>* momentDependency = nullptr;
526 const DependencySourceT<Real>* areaDependency = nullptr;
527 if constexpr (contract::validationsEnabled) {
528 momentDependency = &context.dependencies.requireAll({CentralMoment20, CentralMoment02, CentralMoment11, CentralMoment30,
529 CentralMoment03, CentralMoment21, CentralMoment12});
530 areaDependency = &context.dependencies.require(Area);
531 } else {
532 momentDependency = ::mmcfilters::findDependencySource(context.dependencySources, CentralMoment20);
533 areaDependency = ::mmcfilters::findDependencySource(context.dependencySources, Area);
534 }
535 detail::kernel::computeHuMoments(context, request, *momentDependency, *areaDependency);
536 }
537
538 public:
547 template <std::floating_point Real> static void computeUnitRows(const UnitAttributeComputeContext<Real>& context) {
548 requireUnitAttributeBufferShape(context.tree, context.unitPixels, context.buffer, context.attrNames);
549 constexpr std::array<Attribute, 7> zeroAttributes{HuMoment1, HuMoment2, HuMoment3, HuMoment4, HuMoment5, HuMoment6, HuMoment7};
550 for (const Attribute attribute : zeroAttributes) {
551 if (!requestsAttribute(context.requestedAttributes, attribute)) {
552 continue;
553 }
554 for (NodeId leafIndex = 0; leafIndex < static_cast<NodeId>(context.unitPixels.size()); ++leafIndex) {
555 context.buffer[context.attrNames.linearIndex(leafIndex, attribute)] = Real{0};
556 }
557 }
558 }
559};
560
582class MomentBasedAttributeComputer {
583 public:
585 static constexpr std::string_view familyName = "moment-derived";
586
588 static constexpr AttributeComputerFamily family = AttributeComputerFamily::MomentDerived;
589
591 static constexpr AttributeComputerDomain domain = AttributeComputerDomain::Topology;
592
598 template <std::floating_point Real> static constexpr Real maxFiniteEccentricity() noexcept {
599 return detail::maximumFiniteEccentricity<Real>();
600 }
601
605 inline static constexpr std::array<Attribute, 7> producedAttributes{Inertia, Compactness, Eccentricity, LengthMajorAxis,
606 LengthMinorAxis, AxisOrientation, Circularity};
607
619 template <std::floating_point Real> static void compute(const AttributeComputeContext<Real>& context) {
620 const detail::MomentDerivedRequest request = detail::MomentDerivedRequest::from(context.requestedAttributes);
621 MMCFILTERS_CONTRACT_CHECKED_ONLY(detail::validateMomentDependencyContext(context));
622 if (!request.any()) {
623 return;
624 }
625
626 const DependencySourceT<Real>* momentDependency = nullptr;
627 const DependencySourceT<Real>* areaDependency = nullptr;
628 if constexpr (contract::validationsEnabled) {
629 momentDependency = &context.dependencies.requireAll({CentralMoment20, CentralMoment02, CentralMoment11});
630 areaDependency = &context.dependencies.require(Area);
631 } else {
632 momentDependency = ::mmcfilters::findDependencySource(context.dependencySources, CentralMoment20);
633 areaDependency = ::mmcfilters::findDependencySource(context.dependencySources, Area);
634 }
635 detail::kernel::computeMomentDerived(context, request, *momentDependency, *areaDependency);
636 }
637
638 public:
647 template <std::floating_point Real> static void computeUnitRows(const UnitAttributeComputeContext<Real>& context) {
648 requireUnitAttributeBufferShape(context.tree, context.unitPixels, context.buffer, context.attrNames);
649 constexpr std::array<Attribute, 7> zeroAttributes{Compactness, Eccentricity, LengthMajorAxis, LengthMinorAxis,
650 AxisOrientation, Inertia, Circularity};
651 for (const Attribute attribute : zeroAttributes) {
652 if (!requestsAttribute(context.requestedAttributes, attribute)) {
653 continue;
654 }
655 for (NodeId leafIndex = 0; leafIndex < static_cast<NodeId>(context.unitPixels.size()); ++leafIndex) {
656 context.buffer[context.attrNames.linearIndex(leafIndex, attribute)] =
657 (attribute == Eccentricity || attribute == Circularity) ? Real{1} : Real{0};
658 }
659 }
660 }
661};
662
663} // 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:419
mmcfilters::attributes::computers::CentralMomentsComputer::compute
static void compute(const AttributeComputeContext< Real > &context)
Computes the requested central moments.
Definition MomentBasedAttributeComputer.hpp:448
mmcfilters::attributes::computers::CentralMomentsComputer::computeUnitRows
static void computeUnitRows(const UnitAttributeComputeContext< Real > &context)
Materializes central moments for one-pixel unit supports.
Definition MomentBasedAttributeComputer.hpp:462
mmcfilters::attributes::computers::HuMomentsComputer::compute
static void compute(const AttributeComputeContext< Real > &context)
Computes the requested Hu invariant moments.
Definition MomentBasedAttributeComputer.hpp:518
mmcfilters::attributes::computers::HuMomentsComputer::computeUnitRows
static void computeUnitRows(const UnitAttributeComputeContext< Real > &context)
Materializes Hu moments for one-pixel unit supports.
Definition MomentBasedAttributeComputer.hpp:547
mmcfilters::attributes::computers::MomentBasedAttributeComputer
Computes higher-level shape descriptors derived from second-order central moments.
Definition MomentBasedAttributeComputer.hpp:582
mmcfilters::attributes::computers::MomentBasedAttributeComputer::computeUnitRows
static void computeUnitRows(const UnitAttributeComputeContext< Real > &context)
Materializes moment-derived descriptors for one-pixel unit supports.
Definition MomentBasedAttributeComputer.hpp:647
mmcfilters::attributes::computers::MomentBasedAttributeComputer::maxFiniteEccentricity
static constexpr Real maxFiniteEccentricity() noexcept
Largest finite eccentricity emitted for degenerate one-dimensional supports.
Definition MomentBasedAttributeComputer.hpp:598
mmcfilters::attributes::computers::MomentBasedAttributeComputer::compute
static void compute(const AttributeComputeContext< Real > &context)
Computes the requested moment-derived descriptors.
Definition MomentBasedAttributeComputer.hpp:619
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