mmcfilters
Public API documentation
Loading...
Searching...
No Matches
MaxDistComputer.hpp
1#pragma once
2
3#include "AttributeComputerDomain.hpp"
4#include "AttributeComputerFamily.hpp"
5#include "../detail/AttributeKernelSupport.hpp"
6#include "detail/distance_transform_approx/MorphologicalTreeApproximateDistanceTransform.hpp"
7#include "../../trees/MorphologicalTree.hpp"
8
9#include <algorithm>
10#include <array>
11#include <cmath>
12#include <concepts>
13#include <span>
14#include <stdexcept>
15#include <string_view>
16
18namespace mmcfilters::attributes::computers::detail {
19
20inline void requireMaxDistCapabilities(const MorphologicalTree& tree) {
21 if (!tree.hasGridDomain2D()) {
22 throw std::invalid_argument("Approximate distance-field attributes require a regular 2D pixel domain.");
23 }
24}
25
26struct ApproximateDistanceTransformAttributeRequest {
27 bool maxDist = false;
28 bool maxSquaredDist = false;
29 bool squaredSum = false;
30 bool squaredMean = false;
31 bool rms = false;
32 bool squaredVariance = false;
33 bool centerRow = false;
34 bool centerColumn = false;
35 bool plateauArea = false;
36 bool plateauCentroidRow = false;
37 bool plateauCentroidColumn = false;
38 bool distanceSum = false;
39 bool distanceMean = false;
40 bool distanceVariance = false;
41 bool median = false;
42 bool mode = false;
43 bool q25 = false;
44 bool q75 = false;
45 bool q90 = false;
46 bool entropy = false;
47 bool positiveArea = false;
48 bool levelCount = false;
49 bool weightedCentroidRow = false;
50 bool weightedCentroidColumn = false;
51 bool weightedMu20 = false;
52 bool weightedMu02 = false;
53 bool weightedMu11 = false;
54 bool weightedAxisOrientation = false;
55 bool weightedEccentricity = false;
56
57 [[nodiscard]] bool any() const noexcept { return maximum() || moments() || profile() || spatial() || geometry(); }
58 [[nodiscard]] bool maximum() const noexcept { return maxDist || maxSquaredDist; }
59 [[nodiscard]] bool moments() const noexcept {
60 return squaredSum || squaredMean || rms || squaredVariance || distanceSum || distanceMean || distanceVariance;
61 }
62 [[nodiscard]] bool profile() const noexcept { return median || mode || q25 || q75 || q90 || entropy || positiveArea || levelCount; }
63 [[nodiscard]] bool spatial() const noexcept {
64 return weightedCentroidRow || weightedCentroidColumn || weightedMu20 || weightedMu02 || weightedMu11 || weightedAxisOrientation || weightedEccentricity;
65 }
66 [[nodiscard]] bool localization() const noexcept { return centerRow || centerColumn; }
67 [[nodiscard]] bool plateau() const noexcept { return plateauArea || plateauCentroidRow || plateauCentroidColumn; }
68 [[nodiscard]] bool geometry() const noexcept { return localization() || plateau(); }
69 [[nodiscard]] bool needsCoordinates() const noexcept { return geometry() || spatial(); }
70
71 [[nodiscard]] static ApproximateDistanceTransformAttributeRequest from(std::span<const Attribute> requestedAttributes) {
72 const auto contains = [requestedAttributes](Attribute attribute) {
73 return std::find(requestedAttributes.begin(), requestedAttributes.end(), attribute) != requestedAttributes.end();
74 };
75 return {.maxDist = contains(MaxDist),
76 .maxSquaredDist = contains(MaxSquaredDist),
77 .squaredSum = contains(DistSquaredSum),
78 .squaredMean = contains(DistSquaredMean),
79 .rms = contains(DistRms),
80 .squaredVariance = contains(DistSquaredVariance),
81 .centerRow = contains(MaxDistCenterRow),
82 .centerColumn = contains(MaxDistCenterColumn),
83 .plateauArea = contains(MaxDistPlateauArea),
84 .plateauCentroidRow = contains(MaxDistPlateauCentroidRow),
85 .plateauCentroidColumn = contains(MaxDistPlateauCentroidColumn),
86 .distanceSum = contains(DistSum),
87 .distanceMean = contains(DistMean),
88 .distanceVariance = contains(DistVariance),
89 .median = contains(DistMedian),
90 .mode = contains(DistMode),
91 .q25 = contains(DistQ25),
92 .q75 = contains(DistQ75),
93 .q90 = contains(DistQ90),
94 .entropy = contains(DistEntropy),
95 .positiveArea = contains(DistPositiveArea),
96 .levelCount = contains(DistLevelCount),
97 .weightedCentroidRow = contains(DistWeightedCentroidRow),
98 .weightedCentroidColumn = contains(DistWeightedCentroidColumn),
99 .weightedMu20 = contains(DistWeightedCentralMoment20),
100 .weightedMu02 = contains(DistWeightedCentralMoment02),
101 .weightedMu11 = contains(DistWeightedCentralMoment11),
102 .weightedAxisOrientation = contains(DistWeightedAxisOrientation),
103 .weightedEccentricity = contains(DistWeightedEccentricity)};
104 }
105};
106
107template <std::floating_point Real, class SquaredDistance>
108inline void materializeApproximateMaximum(const AttributeComputeContext<Real>& context, const ApproximateDistanceTransformAttributeRequest& request,
109 NodeId node, SquaredDistance squaredDistance) {
110 if (request.maxDist) {
111 context.buffer[context.attrNames.linearIndex(node, MaxDist)] = static_cast<Real>(std::sqrt(static_cast<long double>(squaredDistance)));
112 }
113 if (request.maxSquaredDist) {
114 context.buffer[context.attrNames.linearIndex(node, MaxSquaredDist)] = static_cast<Real>(squaredDistance);
115 }
116}
117
118template <std::floating_point Real>
119inline void materializeApproximateDistanceFieldExtremum(const AttributeComputeContext<Real>& context,
120 const ApproximateDistanceTransformAttributeRequest& request, NodeId node,
121 const distance_transform::DistanceFieldExtremum& extremum, int numColumns) {
122 if (extremum.pixel == InvalidPixel) {
123 throw std::logic_error("Approximate distance-field localization produced an empty live-node support.");
124 }
125 const auto materialize = [&context, node](Attribute attribute, auto value) {
126 context.buffer[context.attrNames.linearIndex(node, attribute)] = static_cast<Real>(value);
127 };
128 materializeApproximateMaximum(context, request, node, extremum.squaredDistance);
129 if (request.centerRow) {
130 materialize(MaxDistCenterRow, extremum.pixel / numColumns);
131 }
132 if (request.centerColumn) {
133 materialize(MaxDistCenterColumn, extremum.pixel % numColumns);
134 }
135}
136
137template <std::floating_point Real>
138inline void materializeApproximateDistanceFieldPlateau(const AttributeComputeContext<Real>& context,
139 const ApproximateDistanceTransformAttributeRequest& request, NodeId node,
140 const distance_transform::DistanceFieldMaximumPlateau& plateau, int numColumns) {
141 if (plateau.pixel == InvalidPixel || plateau.count == 0) {
142 throw std::logic_error("Approximate distance-field plateau reduction produced an empty live-node support.");
143 }
144 const auto materialize = [&context, node](Attribute attribute, auto value) {
145 context.buffer[context.attrNames.linearIndex(node, attribute)] = static_cast<Real>(value);
146 };
147 materializeApproximateMaximum(context, request, node, plateau.squaredDistance);
148 if (request.centerRow) {
149 materialize(MaxDistCenterRow, plateau.pixel / numColumns);
150 }
151 if (request.centerColumn) {
152 materialize(MaxDistCenterColumn, plateau.pixel % numColumns);
153 }
154 if (request.plateauArea) {
155 materialize(MaxDistPlateauArea, plateau.count);
156 }
157 if (request.plateauCentroidRow) {
158 materialize(MaxDistPlateauCentroidRow, plateau.centroidRow());
159 }
160 if (request.plateauCentroidColumn) {
161 materialize(MaxDistPlateauCentroidColumn, plateau.centroidColumn());
162 }
163}
164
165template <std::floating_point Real>
166inline void materializeApproximateDistanceFieldMoments(const AttributeComputeContext<Real>& context,
167 const ApproximateDistanceTransformAttributeRequest& request, NodeId node,
168 const distance_transform::DistanceFieldMoments& moments) {
169 const auto materialize = [&context, node](Attribute attribute, long double value) {
170 context.buffer[context.attrNames.linearIndex(node, attribute)] = static_cast<Real>(value);
171 };
172 if (request.squaredSum) {
173 materialize(DistSquaredSum, moments.sum());
174 }
175 if (request.squaredMean) {
176 materialize(DistSquaredMean, moments.mean());
177 }
178 if (request.rms) {
179 materialize(DistRms, moments.rms());
180 }
181 if (request.squaredVariance) {
182 materialize(DistSquaredVariance, moments.populationVariance());
183 }
184 if (request.distanceSum) {
185 materialize(DistSum, moments.distanceSum());
186 }
187 if (request.distanceMean) {
188 materialize(DistMean, moments.distanceMean());
189 }
190 if (request.distanceVariance) {
191 materialize(DistVariance, moments.distancePopulationVariance());
192 }
193}
194
195template <std::floating_point Real>
196inline void materializeApproximateDistanceFieldHistogram(const AttributeComputeContext<Real>& context,
197 const ApproximateDistanceTransformAttributeRequest& request, NodeId node,
198 const distance_transform::DistanceFieldHistogram& histogram) {
199 const auto materialize = [&context, node](Attribute attribute, auto value) {
200 context.buffer[context.attrNames.linearIndex(node, attribute)] = static_cast<Real>(value);
201 };
202 if (request.median) {
203 materialize(DistMedian, histogram.quantile(0.5L));
204 }
205 if (request.mode) {
206 materialize(DistMode, histogram.mode());
207 }
208 if (request.q25) {
209 materialize(DistQ25, histogram.quantile(0.25L));
210 }
211 if (request.q75) {
212 materialize(DistQ75, histogram.quantile(0.75L));
213 }
214 if (request.q90) {
215 materialize(DistQ90, histogram.quantile(0.9L));
216 }
217 if (request.entropy) {
218 materialize(DistEntropy, histogram.entropyBits());
219 }
220 if (request.positiveArea) {
221 materialize(DistPositiveArea, histogram.positiveArea());
222 }
223 if (request.levelCount) {
224 materialize(DistLevelCount, histogram.levelCount());
225 }
226}
227
228template <std::floating_point Real>
229inline void materializeApproximateDistanceFieldSpatialMoments(const AttributeComputeContext<Real>& context,
230 const ApproximateDistanceTransformAttributeRequest& request, NodeId node,
231 const distance_transform::DistanceWeightedSpatialMoments& spatialMoments) {
232 const auto materialize = [&context, node](Attribute attribute, auto value) {
233 context.buffer[context.attrNames.linearIndex(node, attribute)] = static_cast<Real>(value);
234 };
235 if (request.weightedCentroidRow) {
236 materialize(DistWeightedCentroidRow, spatialMoments.centroidRow());
237 }
238 if (request.weightedCentroidColumn) {
239 materialize(DistWeightedCentroidColumn, spatialMoments.centroidColumn());
240 }
241 if (request.weightedMu20) {
242 materialize(DistWeightedCentralMoment20, spatialMoments.centralMoment20());
243 }
244 if (request.weightedMu02) {
245 materialize(DistWeightedCentralMoment02, spatialMoments.centralMoment02());
246 }
247 if (request.weightedMu11) {
248 materialize(DistWeightedCentralMoment11, spatialMoments.centralMoment11());
249 }
250 if (request.weightedAxisOrientation) {
251 materialize(DistWeightedAxisOrientation, spatialMoments.axisOrientationDegrees());
252 }
253 if (request.weightedEccentricity) {
254 materialize(DistWeightedEccentricity, spatialMoments.eccentricity());
255 }
256}
257
258template <bool TrackMoments, bool TrackHistogram, bool TrackSpatialMoments, std::floating_point Real>
259inline void computeSelectedApproximateDistanceTransformAttributes(const AttributeComputeContext<Real>& context,
260 const ApproximateDistanceTransformAttributeRequest& request, int numColumns) {
261 using Transform = distance_transform_approx::MorphologicalTreeApproximateDistanceTransform;
262 const auto consume = [&context, request, numColumns](NodeId node, const distance_transform::DistanceFieldExtremum& extremum,
263 const distance_transform::DistanceFieldMaximumPlateau& plateau, const auto& observer,
264 PixelId representative) {
265 if (request.plateau()) {
266 materializeApproximateDistanceFieldPlateau(context, request, node, plateau, numColumns);
267 } else if (request.localization()) {
268 materializeApproximateDistanceFieldExtremum(context, request, node, extremum, numColumns);
269 } else {
270 materializeApproximateMaximum(context, request, node, extremum.squaredDistance);
271 }
272 if constexpr (TrackMoments) {
273 materializeApproximateDistanceFieldMoments(context, request, node, observer.momentsFor(representative));
274 }
275 if constexpr (TrackHistogram) {
276 materializeApproximateDistanceFieldHistogram(context, request, node, observer.histogramFor(representative));
277 }
278 if constexpr (TrackSpatialMoments) {
279 materializeApproximateDistanceFieldSpatialMoments(context, request, node, observer.spatialMomentsFor(representative));
280 }
281 };
282
283 if (request.plateau()) {
284 Transform::template forEachNodeSelectedStatistics<TrackMoments, TrackHistogram, TrackSpatialMoments, false, true>(context.tree, consume);
285 return;
286 }
287 if (request.localization()) {
288 Transform::template forEachNodeSelectedStatistics<TrackMoments, TrackHistogram, TrackSpatialMoments, true, false>(context.tree, consume);
289 return;
290 }
291 Transform::template forEachNodeSelectedStatistics<TrackMoments, TrackHistogram, TrackSpatialMoments>(context.tree, consume);
292}
293
294template <std::floating_point Real>
295inline void dispatchSelectedApproximateDistanceTransformAttributes(const AttributeComputeContext<Real>& context,
296 const ApproximateDistanceTransformAttributeRequest& request, int numColumns) {
297 if (request.profile()) {
298 if (request.spatial()) {
299 if (request.moments()) {
300 computeSelectedApproximateDistanceTransformAttributes<true, true, true>(context, request, numColumns);
301 return;
302 }
303 computeSelectedApproximateDistanceTransformAttributes<false, true, true>(context, request, numColumns);
304 return;
305 }
306 if (request.moments()) {
307 computeSelectedApproximateDistanceTransformAttributes<true, true, false>(context, request, numColumns);
308 return;
309 }
310 computeSelectedApproximateDistanceTransformAttributes<false, true, false>(context, request, numColumns);
311 return;
312 }
313 if (request.moments()) {
314 computeSelectedApproximateDistanceTransformAttributes<true, false, true>(context, request, numColumns);
315 return;
316 }
317 computeSelectedApproximateDistanceTransformAttributes<false, false, true>(context, request, numColumns);
318}
319
320template <std::floating_point Real>
321inline void computeApproximateDistanceTransformAttributes(const AttributeComputeContext<Real>& context,
322 const ApproximateDistanceTransformAttributeRequest& request) {
323 using Transform = distance_transform_approx::MorphologicalTreeApproximateDistanceTransform;
324 const int numColumns = request.needsCoordinates() ? context.tree.numColumns() : 0;
325 if (request.profile() || request.spatial()) {
326 dispatchSelectedApproximateDistanceTransformAttributes(context, request, numColumns);
327 return;
328 }
329 if (request.moments()) {
330 if (request.plateau()) {
331 Transform::forEachNodeSummaryAndPlateau(context.tree,
332 [&context, request, numColumns](NodeId node, const distance_transform::DistanceFieldMaximumPlateau& plateau,
333 const distance_transform::DistanceFieldMoments& moments) {
334 materializeApproximateDistanceFieldPlateau(context, request, node, plateau, numColumns);
335 materializeApproximateDistanceFieldMoments(context, request, node, moments);
336 });
337 return;
338 }
339 if (request.localization()) {
340 Transform::forEachNodeSummaryAndExtremum(context.tree,
341 [&context, request, numColumns](NodeId node, const distance_transform::DistanceFieldExtremum& extremum,
342 const distance_transform::DistanceFieldMoments& moments) {
343 materializeApproximateDistanceFieldExtremum(context, request, node, extremum, numColumns);
344 materializeApproximateDistanceFieldMoments(context, request, node, moments);
345 });
346 return;
347 }
348 Transform::forEachNodeSummary(context.tree, [&context, request](NodeId node, distance_transform_approx::ApproxSquaredDistance maximum,
349 const distance_transform::DistanceFieldMoments& moments) {
350 materializeApproximateMaximum(context, request, node, maximum);
351 materializeApproximateDistanceFieldMoments(context, request, node, moments);
352 });
353 return;
354 }
355 if (request.plateau()) {
356 Transform::forEachNodePlateau(context.tree,
357 [&context, request, numColumns](NodeId node, const distance_transform::DistanceFieldMaximumPlateau& plateau) {
358 materializeApproximateDistanceFieldPlateau(context, request, node, plateau, numColumns);
359 });
360 return;
361 }
362 if (request.localization()) {
363 Transform::forEachNodeExtremum(context.tree, [&context, request, numColumns](NodeId node, const distance_transform::DistanceFieldExtremum& extremum) {
364 materializeApproximateDistanceFieldExtremum(context, request, node, extremum, numColumns);
365 });
366 return;
367 }
368 Transform::forEachNodeMaximum(context.tree, [&context, request](NodeId node, distance_transform_approx::ApproxSquaredDistance value) {
369 materializeApproximateMaximum(context, request, node, value);
370 });
371}
372
373} // namespace mmcfilters::attributes::computers::detail
376namespace mmcfilters::attributes::computers {
377
379class MaxDistComputer {
380 public:
382 static constexpr std::string_view familyName = "max-dist";
383
385 static constexpr AttributeComputerFamily family = AttributeComputerFamily::MaxDist;
386
388 static constexpr AttributeComputerDomain domain = AttributeComputerDomain::Topology;
389
391 inline static constexpr std::array<Attribute, 29> producedAttributes{MaxDist,
392 MaxSquaredDist,
393 DistSquaredSum,
394 DistSquaredMean,
395 DistRms,
396 DistSquaredVariance,
397 MaxDistCenterRow,
398 MaxDistCenterColumn,
399 MaxDistPlateauArea,
400 MaxDistPlateauCentroidRow,
401 MaxDistPlateauCentroidColumn,
402 DistSum,
403 DistMean,
404 DistVariance,
405 DistMedian,
406 DistMode,
407 DistQ25,
408 DistQ75,
409 DistQ90,
410 DistEntropy,
411 DistPositiveArea,
412 DistLevelCount,
413 DistWeightedCentroidRow,
414 DistWeightedCentroidColumn,
415 DistWeightedCentralMoment20,
416 DistWeightedCentralMoment02,
417 DistWeightedCentralMoment11,
418 DistWeightedAxisOrientation,
419 DistWeightedEccentricity};
420
422 static void requireSupportedTreeKind(const MorphologicalTree& tree) { detail::requireMaxDistCapabilities(tree); }
423
425 template <std::floating_point Real> static void compute(const AttributeComputeContext<Real>& context) {
426 const detail::ApproximateDistanceTransformAttributeRequest request =
427 detail::ApproximateDistanceTransformAttributeRequest::from(context.requestedAttributes);
428 MMCFILTERS_CONTRACT_CHECKED_ONLY(requireAttributeBufferShape(context.tree, context.buffer, context.attrNames);
429 requireRequestedAttributeColumns(context); detail::requireMaxDistCapabilities(context.tree));
430 if (request.any()) {
431 static_cast<void>(detail::computeApproximateDistanceTransformAttributes(context, request));
432 }
433 }
434
436 template <std::floating_point Real> static void computeUnitRows(const UnitAttributeComputeContext<Real>& context) {
437 requireUnitAttributeBufferShape(context.tree, context.unitPixels, context.buffer, context.attrNames);
438 const int numColumns = context.tree.numColumns();
439 for (NodeId leafIndex = 0; leafIndex < static_cast<NodeId>(context.unitPixels.size()); ++leafIndex) {
440 const PixelId pixel = context.unitPixels[static_cast<std::size_t>(leafIndex)];
441 for (Attribute attribute : producedAttributes) {
442 if (!requestsAttribute(context.requestedAttributes, attribute)) {
443 continue;
444 }
445 Real value = Real{0};
446 if (attribute == MaxDistCenterRow) {
447 value = static_cast<Real>(pixel / numColumns);
448 } else if (attribute == MaxDistCenterColumn) {
449 value = static_cast<Real>(pixel % numColumns);
450 } else if (attribute == MaxDistPlateauArea) {
451 value = Real{1};
452 } else if (attribute == MaxDistPlateauCentroidRow) {
453 value = static_cast<Real>(pixel / numColumns);
454 } else if (attribute == MaxDistPlateauCentroidColumn) {
455 value = static_cast<Real>(pixel % numColumns);
456 } else if (attribute == DistLevelCount) {
457 value = Real{1};
458 } else if (attribute == DistWeightedCentroidRow) {
459 value = static_cast<Real>(pixel / numColumns);
460 } else if (attribute == DistWeightedCentroidColumn) {
461 value = static_cast<Real>(pixel % numColumns);
462 } else if (attribute == DistWeightedEccentricity) {
463 value = Real{1};
464 }
465 context.buffer[context.attrNames.linearIndex(leafIndex, attribute)] = value;
466 }
467 }
468 }
469};
470
471} // namespace mmcfilters::attributes::computers
mmcfilters::PixelId
int PixelId
Pixel identifier type used by source and active construction domains.
Definition Common.hpp:26
mmcfilters::NodeId
int NodeId
Node identifier type used throughout the project.
Definition Common.hpp:17
mmcfilters::contract::mode
constexpr Mode mode
Contract policy selected when this translation unit was compiled.
Definition Contract.hpp:40
MMCFILTERS_CONTRACT_CHECKED_ONLY
#define MMCFILTERS_CONTRACT_CHECKED_ONLY(...)
Executes validation statements only when defensive checks are enabled.
Definition Contract.hpp:67
mmcfilters::MorphologicalTree
Mutable connected-subset tree on a finite pixel domain.
Definition MorphologicalTree.hpp:114
mmcfilters::attributes::computers::MaxDistComputer
Topology-only approximate distance-transform attribute computer.
Definition MaxDistComputer.hpp:379
mmcfilters::attributes::computers::MaxDistComputer::familyName
static constexpr std::string_view familyName
Family name used in dependency-plan diagnostics.
Definition MaxDistComputer.hpp:382
mmcfilters::attributes::computers::MaxDistComputer::producedAttributes
static constexpr std::array< Attribute, 29 > producedAttributes
Canonical list of scalar descriptors materialized by this computer.
Definition MaxDistComputer.hpp:391
mmcfilters::attributes::computers::MaxDistComputer::compute
static void compute(const AttributeComputeContext< Real > &context)
Computes requested approximate distance-transform projections.
Definition MaxDistComputer.hpp:425
mmcfilters::attributes::computers::MaxDistComputer::computeUnitRows
static void computeUnitRows(const UnitAttributeComputeContext< Real > &context)
Materializes approximate unit-support values.
Definition MaxDistComputer.hpp:436
mmcfilters::attributes::computers::MaxDistComputer::domain
static constexpr AttributeComputerDomain domain
Approximate distance-transform attributes require topology/support, not altitudes.
Definition MaxDistComputer.hpp:388
mmcfilters::attributes::computers::MaxDistComputer::family
static constexpr AttributeComputerFamily family
Stable family id used by the scheduler.
Definition MaxDistComputer.hpp:385
mmcfilters::attributes::computers::MaxDistComputer::requireSupportedTreeKind
static void requireSupportedTreeKind(const MorphologicalTree &tree)
Validates the topology-only geometric capability contract.
Definition MaxDistComputer.hpp:422
mmcfilters::ComputedAttributeData
Owning result for one computed scalar attribute layout and buffer.
Definition AttributeResultTypes.hpp:32
Generated by doxygen 1.9.8