diff --git a/Modules/Core/Common/include/itkGreedyReduceAlgorithm.h b/Modules/Core/Common/include/itkGreedyReduceAlgorithm.h new file mode 100644 index 00000000000..5b042eb9276 --- /dev/null +++ b/Modules/Core/Common/include/itkGreedyReduceAlgorithm.h @@ -0,0 +1,123 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkGreedyReduceAlgorithm_h +#define itkGreedyReduceAlgorithm_h + +#include "itkReduceAlgorithm.h" + +namespace itk +{ +/** \class GreedyReduceAlgorithm + * \brief Thread-safe parallel reduce using a greedy swap-and-merge strategy. + * + * Implements the parallel reduction pattern used in several ITK filters + * (e.g. LabelStatisticsImageFilter, LabelOverlapMeasuresImageFilter, + * ImageToHistogramFilter). Each thread calls Merge() with its local + * partial result; the implementation combines them without holding the + * global mutex while performing the (potentially expensive) merge step. + * + * \par Non-deterministic merge order + * The order in which partial results are combined depends on thread + * scheduling and is not deterministic across runs. Operations whose + * result is independent of evaluation order (e.g. integer counts, + * minimum, maximum) produce bit-identical outputs. Floating-point + * reductions (e.g. running sums, means) may produce results that differ + * in the last few ULPs between runs due to non-associativity of + * floating-point arithmetic. + * + * Algorithm (per Merge() call): + * -# Acquire the mutex. + * -# If no result is accumulated yet, store the local result and return. + * -# Otherwise, atomically take ownership of the current accumulated result, + * clearing the shared state to allow other threads to proceed immediately. + * -# Release the mutex. + * -# Merge the taken result into the local result (off the critical path). + * -# Repeat from step 1 until the local result is successfully deposited. + * + * A merge function must be supplied via SetMergeFunction() before calling + * Merge(). The function signature is: + * \code + * void mergeFunction(T & target, T & source); + * \endcode + * It must merge \c source into \c target. + * + * \tparam T The type of the object being reduced. + * + * \ingroup ITKCommon + */ +template +class ITK_TEMPLATE_EXPORT GreedyReduceAlgorithm : public ReduceAlgorithm +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(GreedyReduceAlgorithm); + + /** Standard class type aliases. */ + using Self = GreedyReduceAlgorithm; + using Superclass = ReduceAlgorithm; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** \see LightObject::GetNameOfClass() */ + itkOverrideGetNameOfClassMacro(GreedyReduceAlgorithm); + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** The value type being reduced. */ + using typename Superclass::ValueType; + + /** Bring the chunk-ID overload into scope; Greedy ignores the ID and + * forwards to Merge(T&&) (see ReduceAlgorithm::Merge(SizeValueType, T&&)). */ + using Superclass::Merge; + + /** Merge a per-thread \p localResult into the accumulated output. + * This method is thread-safe; multiple threads may call it concurrently. */ + void + Merge(T && localResult) override; + + /** Return the accumulated result. + * Should only be called after all concurrent Merge() calls have finished. */ + const T & + GetResult() const override; + + /** Reset accumulated state so this object can be reused. */ + void + Clear() override; + +protected: + GreedyReduceAlgorithm() = default; + ~GreedyReduceAlgorithm() override = default; + + void + PrintSelf(std::ostream & os, Indent indent) const override; + +private: + /** Accumulated result. Default-constructed value serves as initial state. */ + T m_Result{}; + + /** True once the first Merge() has deposited a value. */ + bool m_HasResult{ false }; +}; + +} // namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkGreedyReduceAlgorithm.hxx" +#endif + +#endif // itkGreedyReduceAlgorithm_h diff --git a/Modules/Core/Common/include/itkGreedyReduceAlgorithm.hxx b/Modules/Core/Common/include/itkGreedyReduceAlgorithm.hxx new file mode 100644 index 00000000000..6d2229c39fe --- /dev/null +++ b/Modules/Core/Common/include/itkGreedyReduceAlgorithm.hxx @@ -0,0 +1,88 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkGreedyReduceAlgorithm_hxx +#define itkGreedyReduceAlgorithm_hxx + +#include "itkGreedyReduceAlgorithm.h" + +#include + +namespace itk +{ + +template +void +GreedyReduceAlgorithm::Merge(T && localResult) +{ + // Greedy swap-and-merge strategy: + // Acquire the mutex only long enough to swap ownership of the shared + // accumulator. The actual merge work is done outside the critical section + // so other threads can proceed concurrently. + while (true) + { + T tomerge{}; + { + const std::lock_guard lockGuard(this->m_Mutex); + + if (!m_HasResult) + { + // No accumulated result yet: store the local result and return. + m_Result = std::move(localResult); + m_HasResult = true; + return; + } + + // Take ownership of the current accumulated result so other threads + // can deposit their own results immediately after we release the lock. + std::swap(m_Result, tomerge); + m_HasResult = false; + } // release lock + + // Merge the taken result into localResult outside the critical section. + this->m_MergeFunction(localResult, tomerge); + } +} + +template +const T & +GreedyReduceAlgorithm::GetResult() const +{ + return m_Result; +} + +template +void +GreedyReduceAlgorithm::Clear() +{ + const std::lock_guard lockGuard(this->m_Mutex); + m_Result = T{}; + m_HasResult = false; + Superclass::Clear(); +} + +template +void +GreedyReduceAlgorithm::PrintSelf(std::ostream & os, Indent indent) const +{ + Superclass::PrintSelf(os, indent); + os << indent << "HasResult: " << (m_HasResult ? "true" : "false") << std::endl; +} + +} // namespace itk + +#endif // itkGreedyReduceAlgorithm_hxx diff --git a/Modules/Core/Common/include/itkLinearReduceAlgorithm.h b/Modules/Core/Common/include/itkLinearReduceAlgorithm.h new file mode 100644 index 00000000000..c0d724231cf --- /dev/null +++ b/Modules/Core/Common/include/itkLinearReduceAlgorithm.h @@ -0,0 +1,120 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkLinearReduceAlgorithm_h +#define itkLinearReduceAlgorithm_h + +#include "itkReduceAlgorithm.h" + +#include +#include + +namespace itk +{ +/** \class LinearReduceAlgorithm + * \brief Deterministic reduce using an ordered linear merge strategy. + * + * Each work unit deposits its partial result into a chunk-indexed array. + * The actual merge is deferred until GetResult() is called: it walks the + * array once, in ascending chunk-ID order, and combines the values under a + * single lock. Because the merge order is fixed by chunk ID rather than by + * thread scheduling, the final result is bit-identical across runs. + * + * \par Usage + * Call SetNumberOfWorkUnits() \em before any Merge() calls; this allocates + * the internal array. Each work unit then calls Merge(chunkId, value) + * exactly once with its 0-based chunk identifier, in any order. Call + * GetResult() only after every chunk has called Merge(). + * + * \tparam T The type of the object being reduced. + * + * \ingroup ITKCommon + */ +template +class ITK_TEMPLATE_EXPORT LinearReduceAlgorithm : public ReduceAlgorithm +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(LinearReduceAlgorithm); + + /** Standard class type aliases. */ + using Self = LinearReduceAlgorithm; + using Superclass = ReduceAlgorithm; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** \see LightObject::GetNameOfClass() */ + itkOverrideGetNameOfClassMacro(LinearReduceAlgorithm); + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** The value type being reduced. */ + using typename Superclass::ValueType; + + /** Set the total number of work units. Allocates the internal array. + * Must be called before the first Merge(). */ + void + SetNumberOfWorkUnits(SizeValueType numberOfWorkUnits) override; + + /** Not supported — use Merge(SizeValueType chunkId, T &&) instead. + * Always throws itk::ExceptionObject. */ + void + Merge(T && localResult) override; + + /** Deposit \p localResult for chunk \p chunkId. Thread-safe; merely + * stores the value, the actual merge happens lazily in GetResult(). */ + void + Merge(SizeValueType chunkId, T && localResult) override; + + /** Merge every deposited chunk, in ascending chunk-ID order, and return + * the result. The merge is performed once, under a lock, and cached; + * later calls return the cached result. Behaviour is defined only after + * every chunk has called Merge(). */ + const T & + GetResult() const override; + + /** Reset the array so this object can be reused with the same N. + * Does not change the merge function or work-unit count. */ + void + Clear() override; + +protected: + LinearReduceAlgorithm() = default; + ~LinearReduceAlgorithm() override = default; + + void + PrintSelf(std::ostream & os, Indent indent) const override; + +private: + /** Chunk-indexed array of partial results awaiting the merge in + * GetResult(). Mutable because GetResult() consumes it lazily. */ + mutable std::vector> m_Values{}; + + /** Cached result of merging m_Values, computed lazily by GetResult(). */ + mutable T m_Result{}; + + /** True once GetResult() has performed the merge. */ + mutable bool m_Merged{ false }; +}; + +} // namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkLinearReduceAlgorithm.hxx" +#endif + +#endif // itkLinearReduceAlgorithm_h diff --git a/Modules/Core/Common/include/itkLinearReduceAlgorithm.hxx b/Modules/Core/Common/include/itkLinearReduceAlgorithm.hxx new file mode 100644 index 00000000000..ce732c7e91a --- /dev/null +++ b/Modules/Core/Common/include/itkLinearReduceAlgorithm.hxx @@ -0,0 +1,121 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkLinearReduceAlgorithm_hxx +#define itkLinearReduceAlgorithm_hxx + +#include "itkLinearReduceAlgorithm.h" + +#include + +namespace itk +{ + +template +void +LinearReduceAlgorithm::SetNumberOfWorkUnits(SizeValueType numberOfWorkUnits) +{ + Superclass::SetNumberOfWorkUnits(numberOfWorkUnits); + + const std::lock_guard lockGuard(this->m_Mutex); + m_Values.assign(numberOfWorkUnits, std::nullopt); + m_Result = T{}; + m_Merged = false; +} + +template +void +LinearReduceAlgorithm::Merge(T && /*localResult*/) +{ + itkExceptionMacro("LinearReduceAlgorithm::Merge(T&&) requires a chunk ID. " + "Call Merge(SizeValueType chunkId, T&&) instead."); +} + +template +void +LinearReduceAlgorithm::Merge(SizeValueType chunkId, T && localResult) +{ + const std::lock_guard lockGuard(this->m_Mutex); + + if (m_Values.empty()) + { + itkExceptionMacro("LinearReduceAlgorithm::Merge() called before SetNumberOfWorkUnits()."); + } + if (chunkId >= this->m_NumberOfWorkUnits) + { + itkExceptionMacro("LinearReduceAlgorithm::Merge(): chunkId " << chunkId << " >= NumberOfWorkUnits " + << this->m_NumberOfWorkUnits); + } + + // Just deposit the value; GetResult() performs the actual ordered merge. + m_Values[chunkId] = std::move(localResult); +} + +template +const T & +LinearReduceAlgorithm::GetResult() const +{ + const std::lock_guard lockGuard(this->m_Mutex); + + if (!m_Merged) + { + for (auto & value : m_Values) + { + if (!value.has_value()) + { + continue; + } + if (m_Merged) + { + this->m_MergeFunction(m_Result, *value); + } + else + { + m_Result = std::move(*value); + m_Merged = true; + } + value.reset(); + } + } + return m_Result; +} + +template +void +LinearReduceAlgorithm::Clear() +{ + const std::lock_guard lockGuard(this->m_Mutex); + for (auto & v : m_Values) + { + v.reset(); + } + m_Result = T{}; + m_Merged = false; + this->Modified(); +} + +template +void +LinearReduceAlgorithm::PrintSelf(std::ostream & os, Indent indent) const +{ + Superclass::PrintSelf(os, indent); + os << indent << "Merged: " << (m_Merged ? "true" : "false") << std::endl; +} + +} // namespace itk + +#endif // itkLinearReduceAlgorithm_hxx diff --git a/Modules/Core/Common/include/itkReduceAlgorithm.h b/Modules/Core/Common/include/itkReduceAlgorithm.h new file mode 100644 index 00000000000..b9c7b7172e4 --- /dev/null +++ b/Modules/Core/Common/include/itkReduceAlgorithm.h @@ -0,0 +1,144 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkReduceAlgorithm_h +#define itkReduceAlgorithm_h + +#include "itkObject.h" +#include "itkObjectFactory.h" + +#include +#include + +namespace itk +{ +/** \class ReduceAlgorithm + * \brief Abstract base class for thread-safe parallel reduce (merge) algorithms. + * + * Defines the interface for combining per-thread local results into a single + * accumulated output. Concrete subclasses implement `Merge()` and + * `GetResult()` and are free to choose any reduction strategy (greedy, + * tree-based, etc.). + * + * The merge function has the signature: + * \code + * void mergeFunction(T & target, T & source); + * \endcode + * and must merge \c source \em into \c target. + * + * \tparam T The type of the object being reduced. + * + * \ingroup ITKCommon + */ +template +class ITK_TEMPLATE_EXPORT ReduceAlgorithm : public Object +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(ReduceAlgorithm); + + /** Standard class type aliases. */ + using Self = ReduceAlgorithm; + using Superclass = Object; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** \see LightObject::GetNameOfClass() */ + itkOverrideGetNameOfClassMacro(ReduceAlgorithm); + + /** The value type being reduced. */ + using ValueType = T; + + /** Merge function signature: merges \p source into \p target. */ + using MergeFunctionType = std::function; + + /** Set the function used to merge two objects. */ + virtual void + SetMergeFunction(MergeFunctionType mergeFunction) + { + m_MergeFunction = std::move(mergeFunction); + this->Modified(); + } + + /** Get the merge function. */ + const MergeFunctionType & + GetMergeFunction() const + { + return m_MergeFunction; + } + + /** Set the expected number of work units (chunks) that will call Merge(). + * Stored for use by concrete subclasses or future algorithms; not enforced + * by all implementations. */ + itkSetMacro(NumberOfWorkUnits, SizeValueType); + + /** Get the number of work units. */ + itkGetConstMacro(NumberOfWorkUnits, SizeValueType); + + /** Merge a local result into the accumulated output. + * Implementations must be thread-safe. */ + virtual void + Merge(T && localResult) = 0; + + /** Merge a local result identified by \p chunkId (0-based) into the + * accumulated output. The default implementation ignores the ID and + * delegates to Merge(T&&). Subclasses that need deterministic ordering + * (e.g. TreeReduceAlgorithm) override this to exploit the chunk ID. */ + virtual void + Merge(SizeValueType /*chunkId*/, T && localResult) + { + this->Merge(std::move(localResult)); + } + + /** Return the accumulated result after all Merge() calls have completed. + * The returned reference is only stable while no concurrent Merge() or + * Clear() is in progress. */ + virtual const T & + GetResult() const = 0; + + /** Reset internal state so the object can be reused for a new reduction. */ + virtual void + Clear() + { + m_NumberOfWorkUnits = 0; + this->Modified(); + } + +protected: + ReduceAlgorithm() = default; + ~ReduceAlgorithm() override = default; + + void + PrintSelf(std::ostream & os, Indent indent) const override + { + Superclass::PrintSelf(os, indent); + os << indent << "NumberOfWorkUnits: " << m_NumberOfWorkUnits << std::endl; + os << indent << "MergeFunction: " << (m_MergeFunction ? "set" : "not set") << std::endl; + } + + /** Merge function provided by the caller. */ + MergeFunctionType m_MergeFunction{}; + + /** Expected number of work units. */ + SizeValueType m_NumberOfWorkUnits{ 0 }; + + /** Mutex protecting the accumulated result. */ + mutable std::mutex m_Mutex{}; +}; + +} // namespace itk + +#endif // itkReduceAlgorithm_h diff --git a/Modules/Core/Common/include/itkTreeReduceAlgorithm.h b/Modules/Core/Common/include/itkTreeReduceAlgorithm.h new file mode 100644 index 00000000000..1cf8f7d4aeb --- /dev/null +++ b/Modules/Core/Common/include/itkTreeReduceAlgorithm.h @@ -0,0 +1,168 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkTreeReduceAlgorithm_h +#define itkTreeReduceAlgorithm_h + +#include "itkReduceAlgorithm.h" + +#include +#include +#include +#include + +namespace itk +{ +/** \class TreeReduceAlgorithm + * \brief Deterministic parallel reduce using a binary-tree merge strategy. + * + * Merges per-work-unit partial results in a fixed binary-tree order, so the + * final value is bit-identical across runs regardless of thread scheduling. + * This makes it suitable for floating-point reductions where reproducibility + * matters. + * + * \par Usage + * Call SetNumberOfWorkUnits() \em before any Merge() calls; this allocates + * and initialises the internal heap-like tree. Each work unit then calls + * Merge(chunkId, value) exactly once with its 0-based chunk identifier. + * GetResult() returns the root value once all chunks have merged. + * + * \par Tree layout + * The tree is stored 1-indexed in a flat array of size \c 2*paddedN, where + * \c paddedN is the smallest power of two ≥ N. Chunk \c k occupies leaf + * index \c paddedN+k; the root is at index 1. Internal nodes are populated + * only by merge operations; leaf data is moved up the tree eagerly after + * each chunk arrives. + * + * \par Merge direction + * At every internal node the \em left child (lower chunk-ID subtree) is the + * target and the \em right child is the source: + * \code + * mergeFunction(leftValue, rightValue); + * \endcode + * This order is independent of which thread executes the merge. + * + * \par Non-blocking behaviour + * Merge(chunkId, value) returns as soon as it cannot continue up the tree + * because the sibling has not yet arrived. There is no spinning or blocking. + * GetResult() is undefined until every chunk has called Merge(). + * + * \par Arbitrary N + * When N is not a power of two the tree is padded with phantom leaves that + * carry no value. Phantom subtrees are pre-processed at build time so that + * real chunks still walk up the tree correctly without performing spurious + * merges. + * + * \par Memory + * Each node's value is a \c std::optional; after both children of a node + * have been merged the children's storage is released (\c reset()), so peak + * memory is O(N) rather than O(2N). + * + * \tparam T The type of the object being reduced. + * + * \ingroup ITKCommon + */ +template +class ITK_TEMPLATE_EXPORT TreeReduceAlgorithm : public ReduceAlgorithm +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(TreeReduceAlgorithm); + + /** Standard class type aliases. */ + using Self = TreeReduceAlgorithm; + using Superclass = ReduceAlgorithm; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** \see LightObject::GetNameOfClass() */ + itkOverrideGetNameOfClassMacro(TreeReduceAlgorithm); + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** The value type being reduced. */ + using typename Superclass::ValueType; + + /** Set the total number of work units. Rebuilds the internal binary tree. + * Must be called before the first Merge(). */ + void + SetNumberOfWorkUnits(SizeValueType numberOfWorkUnits) override; + + /** Not supported — use Merge(SizeValueType chunkId, T &&) instead. + * Always throws itk::ExceptionObject. */ + void + Merge(T && localResult) override; + + /** Deposit \p localResult for chunk \p chunkId and walk up the tree, + * merging completed sibling pairs until either the root is reached or + * the sibling has not yet arrived. Thread-safe; non-blocking. */ + void + Merge(SizeValueType chunkId, T && localResult) override; + + /** Return the fully reduced result at the tree root. + * Behaviour is defined only after every chunk has called Merge(). */ + const T & + GetResult() const override; + + /** Reset the tree so this object can be reused with the same N. + * Does not change the merge function or work-unit count. */ + void + Clear() override; + +protected: + TreeReduceAlgorithm() = default; + ~TreeReduceAlgorithm() override = default; + + void + PrintSelf(std::ostream & os, Indent indent) const override; + +private: + /** Build (or rebuild) the internal heap-like tree based on + * m_NumberOfWorkUnits. Pre-processes phantom leaves so that the atomic + * counters are initialised correctly. */ + void + BuildTree(); + + /** Smallest power of two ≥ m_NumberOfWorkUnits. Leaves live at indices + * [m_PaddedSize, 2*m_PaddedSize). */ + SizeValueType m_PaddedSize{ 0 }; + + /** 1-indexed heap array. Index 0 is unused; index 1 is the root; + * leaves are at [m_PaddedSize, 2*m_PaddedSize). */ + std::vector> m_Values{}; + + /** Per-node atomic counter: incremented by each child that completes. + * Reaches 2 when both children are ready; the second thread to arrive + * performs the merge. Stored as a unique_ptr because std::atomic is + * not copyable/movable. */ + std::unique_ptr[]> m_ChildrenReady{}; + + /** Initial values of m_ChildrenReady (includes phantom pre-increments). + * Used by Clear() to reset the counters without re-analysing phantoms. */ + std::vector m_InitialCounts{}; + + /** Returned by GetResult() when no chunks have merged yet (tree is empty). */ + T m_DefaultResult{}; +}; + +} // namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkTreeReduceAlgorithm.hxx" +#endif + +#endif // itkTreeReduceAlgorithm_h diff --git a/Modules/Core/Common/include/itkTreeReduceAlgorithm.hxx b/Modules/Core/Common/include/itkTreeReduceAlgorithm.hxx new file mode 100644 index 00000000000..6c745aada4a --- /dev/null +++ b/Modules/Core/Common/include/itkTreeReduceAlgorithm.hxx @@ -0,0 +1,251 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ +#ifndef itkTreeReduceAlgorithm_hxx +#define itkTreeReduceAlgorithm_hxx + +#include "itkTreeReduceAlgorithm.h" + +#include + +namespace itk +{ + +// --------------------------------------------------------------------------- +// SetNumberOfWorkUnits +// --------------------------------------------------------------------------- + +template +void +TreeReduceAlgorithm::SetNumberOfWorkUnits(SizeValueType numberOfWorkUnits) +{ + // Store and propagate the Modified() flag via the base-class macro. + Superclass::SetNumberOfWorkUnits(numberOfWorkUnits); + BuildTree(); +} + +// --------------------------------------------------------------------------- +// BuildTree +// --------------------------------------------------------------------------- + +template +void +TreeReduceAlgorithm::BuildTree() +{ + const SizeValueType n = this->m_NumberOfWorkUnits; + + if (n == 0) + { + m_PaddedSize = 0; + m_Values.clear(); + m_ChildrenReady.reset(); + m_InitialCounts.clear(); + return; + } + + // Compute the smallest power of two >= n. + m_PaddedSize = 1; + while (m_PaddedSize < n) + { + m_PaddedSize <<= 1; + } + + // Heap array is 1-indexed; total size = 2 * paddedSize (index 0 unused). + const SizeValueType totalNodes = 2 * m_PaddedSize; + + // Reset all values to empty. + m_Values.assign(totalNodes, std::nullopt); + + // Allocate atomic counters and zero-initialise. + m_ChildrenReady = std::make_unique[]>(totalNodes); + for (SizeValueType i = 0; i < totalNodes; ++i) + { + m_ChildrenReady[i].store(0, std::memory_order_relaxed); + } + + // --- Phantom leaf pre-processing --- + // + // For each phantom leaf (chunk k where n <= k < paddedSize) we increment + // the parent's counter. If a parent's counter reaches 2 (both children + // are phantom) the parent is also fully phantom: walk further up the tree. + // + // We use a plain vector to accumulate the increments first, then apply + // them atomically so the loop logic stays simple. + std::vector preCount(totalNodes, 0); + for (SizeValueType k = n; k < m_PaddedSize; ++k) + { + SizeValueType nodeIdx = m_PaddedSize + k; // phantom leaf index + while (nodeIdx > 1) + { + const SizeValueType parentIdx = nodeIdx / 2; + ++preCount[parentIdx]; + if (preCount[parentIdx] < 2) + { + break; // sibling is real; stop propagation + } + // Both children of this parent are phantom: continue propagation up. + nodeIdx = parentIdx; + } + } + + // Apply pre-counts and store for use by Clear(). + m_InitialCounts.assign(totalNodes, 0); + for (SizeValueType i = 0; i < totalNodes; ++i) + { + m_InitialCounts[i] = preCount[i]; + m_ChildrenReady[i].store(preCount[i], std::memory_order_relaxed); + } +} + +// --------------------------------------------------------------------------- +// Merge(T&&) — not supported, always throws +// --------------------------------------------------------------------------- + +template +void +TreeReduceAlgorithm::Merge(T && /*localResult*/) +{ + itkExceptionMacro("TreeReduceAlgorithm::Merge(T&&) requires a chunk ID. " + "Call Merge(SizeValueType chunkId, T&&) instead."); +} + +// --------------------------------------------------------------------------- +// Merge(SizeValueType, T&&) — tree algorithm +// --------------------------------------------------------------------------- + +template +void +TreeReduceAlgorithm::Merge(SizeValueType chunkId, T && localResult) +{ + if (m_PaddedSize == 0) + { + itkExceptionMacro("TreeReduceAlgorithm::Merge() called before SetNumberOfWorkUnits()."); + } + if (chunkId >= this->m_NumberOfWorkUnits) + { + itkExceptionMacro("TreeReduceAlgorithm::Merge(): chunkId " << chunkId << " >= NumberOfWorkUnits " + << this->m_NumberOfWorkUnits); + } + + // Deposit value at the leaf. + SizeValueType nodeIdx = m_PaddedSize + chunkId; + m_Values[nodeIdx] = std::move(localResult); + + // Walk up the tree. At each internal node: + // - Atomically increment the parent's ready counter. + // - If we are the first child (counter was 0 before increment): return. + // The sibling's thread will complete the merge when it arrives. + // - If we are the second child (counter was 1): perform the merge and + // continue up to the grandparent. + while (nodeIdx > 1) + { + const SizeValueType parentIdx = nodeIdx / 2; + + const int prev = m_ChildrenReady[parentIdx].fetch_add(1, std::memory_order_acq_rel); + if (prev < 1) + { + // First child to arrive; the sibling will do the merge. + return; + } + + // Second child: perform the merge at this parent node. + // Merge direction: left child is target (lower chunk IDs), right is source. + const SizeValueType leftIdx = parentIdx * 2; + const SizeValueType rightIdx = parentIdx * 2 + 1; + + if (m_Values[leftIdx].has_value() && m_Values[rightIdx].has_value()) + { + // Both children have real values — merge right into left. + this->m_MergeFunction(*m_Values[leftIdx], *m_Values[rightIdx]); + m_Values[parentIdx] = std::move(m_Values[leftIdx]); + } + else if (m_Values[leftIdx].has_value()) + { + // Only left has a value (right subtree was phantom). + m_Values[parentIdx] = std::move(m_Values[leftIdx]); + } + else if (m_Values[rightIdx].has_value()) + { + // Only right has a value (left subtree was phantom). + m_Values[parentIdx] = std::move(m_Values[rightIdx]); + } + // else: both phantom; parent stays nullopt. + + // Release children's memory eagerly. + m_Values[leftIdx].reset(); + m_Values[rightIdx].reset(); + + nodeIdx = parentIdx; + } + // Reached the root (nodeIdx == 1); result is in m_Values[1]. +} + +// --------------------------------------------------------------------------- +// GetResult +// --------------------------------------------------------------------------- + +template +const T & +TreeReduceAlgorithm::GetResult() const +{ + if (m_PaddedSize > 0 && m_Values[1].has_value()) + { + return *m_Values[1]; + } + return m_DefaultResult; +} + +// --------------------------------------------------------------------------- +// Clear +// --------------------------------------------------------------------------- + +template +void +TreeReduceAlgorithm::Clear() +{ + // Reset all values to empty. + for (auto & v : m_Values) + { + v.reset(); + } + + // Restore counters to their initial (post-phantom-analysis) state. + const SizeValueType totalNodes = static_cast(m_InitialCounts.size()); + for (SizeValueType i = 0; i < totalNodes; ++i) + { + m_ChildrenReady[i].store(m_InitialCounts[i], std::memory_order_relaxed); + } + + this->Modified(); +} + +// --------------------------------------------------------------------------- +// PrintSelf +// --------------------------------------------------------------------------- + +template +void +TreeReduceAlgorithm::PrintSelf(std::ostream & os, Indent indent) const +{ + Superclass::PrintSelf(os, indent); + os << indent << "PaddedSize: " << m_PaddedSize << std::endl; + os << indent << "HasResult: " << (m_PaddedSize > 0 && m_Values[1].has_value() ? "true" : "false") << std::endl; +} + +} // namespace itk + +#endif // itkTreeReduceAlgorithm_hxx diff --git a/Modules/Core/Common/test/CMakeLists.txt b/Modules/Core/Common/test/CMakeLists.txt index 18e63994a63..99209f13c1d 100644 --- a/Modules/Core/Common/test/CMakeLists.txt +++ b/Modules/Core/Common/test/CMakeLists.txt @@ -1479,6 +1479,7 @@ set( itkIndexRangeGTest.cxx itkIntTypesGTest.cxx itkLightObjectGTest.cxx + itkLinearReduceAlgorithmGTest.cxx itkMakeUniqueForOverwriteGTest.cxx itkMathCastWithRangeCheckGTest.cxx itkMathGTest.cxx diff --git a/Modules/Core/Common/test/itkGreedyReduceAlgorithmGTest.cxx b/Modules/Core/Common/test/itkGreedyReduceAlgorithmGTest.cxx new file mode 100644 index 00000000000..472b2a916a6 --- /dev/null +++ b/Modules/Core/Common/test/itkGreedyReduceAlgorithmGTest.cxx @@ -0,0 +1,368 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +#include "itkGreedyReduceAlgorithm.h" +#include "itkGTest.h" + +#include +#include +#include +#include + +namespace +{ + +// Merge function for int: add source into target. +void +IntMerge(int & target, int & source) +{ + target += source; +} + +// Merge function for std::map: accumulate values per key. +void +MapMerge(std::map & target, std::map & source) +{ + for (auto & [key, value] : source) + { + target[key] += value; + } +} + +// Merge function for std::vector: concatenate source into target. +void +VectorMerge(std::vector & target, std::vector & source) +{ + target.insert(target.end(), source.begin(), source.end()); +} + +} // namespace + +// --------------------------------------------------------------------------- +// Basic object methods (Print, GetNameOfClass) +// --------------------------------------------------------------------------- + +TEST(GreedyReduceAlgorithm, BasicObjectMethods) +{ + using ReducerType = itk::GreedyReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + reducer->SetNumberOfWorkUnits(4); + + ITK_GTEST_EXERCISE_BASIC_OBJECT_METHODS(reducer, GreedyReduceAlgorithm, ReduceAlgorithm); +} + +// --------------------------------------------------------------------------- +// GetResult before any Merge returns default-constructed value +// --------------------------------------------------------------------------- + +TEST(GreedyReduceAlgorithm, GetResultBeforeMerge) +{ + using ReducerType = itk::GreedyReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + + EXPECT_EQ(reducer->GetResult(), int{}); +} + +TEST(GreedyReduceAlgorithm, GetResultBeforeMergeMap) +{ + using MapType = std::map; + using ReducerType = itk::GreedyReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(MapMerge); + + EXPECT_TRUE(reducer->GetResult().empty()); +} + +// --------------------------------------------------------------------------- +// Single Merge stores the value as-is +// --------------------------------------------------------------------------- + +TEST(GreedyReduceAlgorithm, SingleMerge) +{ + using ReducerType = itk::GreedyReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + + int value = 42; + reducer->Merge(std::move(value)); + + EXPECT_EQ(reducer->GetResult(), 42); +} + +// --------------------------------------------------------------------------- +// Sequential merges accumulate correctly +// --------------------------------------------------------------------------- + +TEST(GreedyReduceAlgorithm, SequentialMerges) +{ + using ReducerType = itk::GreedyReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + + constexpr int N = 10; + for (int i = 1; i <= N; ++i) + { + int v = i; + reducer->Merge(std::move(v)); + } + + // sum 1..10 = 55 + EXPECT_EQ(reducer->GetResult(), 55); +} + +// --------------------------------------------------------------------------- +// Clear resets state; subsequent merges start fresh +// --------------------------------------------------------------------------- + +TEST(GreedyReduceAlgorithm, ClearResetsState) +{ + using ReducerType = itk::GreedyReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + + int v1 = 100; + reducer->Merge(std::move(v1)); + EXPECT_EQ(reducer->GetResult(), 100); + + reducer->Clear(); + EXPECT_EQ(reducer->GetResult(), int{}); + + int v2 = 7; + reducer->Merge(std::move(v2)); + EXPECT_EQ(reducer->GetResult(), 7); +} + +// --------------------------------------------------------------------------- +// Custom merge function: take the maximum instead of summing +// --------------------------------------------------------------------------- + +TEST(GreedyReduceAlgorithm, CustomMergeFunction) +{ + using ReducerType = itk::GreedyReduceAlgorithm; + auto reducer = ReducerType::New(); + + // Merge function that keeps the maximum + reducer->SetMergeFunction([](int & target, int & source) { target = std::max(target, source); }); + + for (int v : { 5, 3, 9, 1, 7 }) + { + reducer->Merge(std::move(v)); + } + + EXPECT_EQ(reducer->GetResult(), 9); +} + +// --------------------------------------------------------------------------- +// Concurrent merges (integer): sum N threads * value == expected total +// --------------------------------------------------------------------------- + +TEST(GreedyReduceAlgorithm, ConcurrentMergesInt) +{ + using ReducerType = itk::GreedyReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + + constexpr int numThreads = 100; + constexpr int valuePerThread = 1; + reducer->SetNumberOfWorkUnits(numThreads); + + std::vector threads; + threads.reserve(numThreads); + + for (int i = 0; i < numThreads; ++i) + { + threads.emplace_back([&]() { + int v = valuePerThread; + reducer->Merge(std::move(v)); + }); + } + + for (auto & t : threads) + { + t.join(); + } + + EXPECT_EQ(reducer->GetResult(), numThreads * valuePerThread); +} + +// --------------------------------------------------------------------------- +// Concurrent merges (map): each thread contributes unique keys +// --------------------------------------------------------------------------- + +TEST(GreedyReduceAlgorithm, ConcurrentMergesMap) +{ + using MapType = std::map; + using ReducerType = itk::GreedyReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(MapMerge); + + constexpr int numThreads = 50; + // Thread i inserts key=i with value=i + reducer->SetNumberOfWorkUnits(numThreads); + + std::vector threads; + threads.reserve(numThreads); + + for (int i = 0; i < numThreads; ++i) + { + threads.emplace_back([&, i]() { + MapType local; + local[i] = i; + reducer->Merge(std::move(local)); + }); + } + + for (auto & t : threads) + { + t.join(); + } + + const MapType & result = reducer->GetResult(); + ASSERT_EQ(static_cast(result.size()), numThreads); + for (int i = 0; i < numThreads; ++i) + { + auto it = result.find(i); + ASSERT_NE(it, result.end()) << "Key " << i << " missing from result"; + EXPECT_EQ(it->second, i) << "Wrong value for key " << i; + } +} + +// --------------------------------------------------------------------------- +// Concurrent merges (map, shared keys): values are correctly summed +// --------------------------------------------------------------------------- + +TEST(GreedyReduceAlgorithm, ConcurrentMergesMapSharedKeys) +{ + using MapType = std::map; + using ReducerType = itk::GreedyReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(MapMerge); + + constexpr int numThreads = 50; + constexpr int numKeys = 5; + // Every thread contributes 1 to each of the numKeys keys. + reducer->SetNumberOfWorkUnits(numThreads); + + std::vector threads; + threads.reserve(numThreads); + + for (int i = 0; i < numThreads; ++i) + { + threads.emplace_back([&]() { + MapType local; + for (int k = 0; k < numKeys; ++k) + { + local[k] = 1; + } + reducer->Merge(std::move(local)); + }); + } + + for (auto & t : threads) + { + t.join(); + } + + const MapType & result = reducer->GetResult(); + ASSERT_EQ(static_cast(result.size()), numKeys); + for (int k = 0; k < numKeys; ++k) + { + auto it = result.find(k); + ASSERT_NE(it, result.end()) << "Key " << k << " missing"; + EXPECT_EQ(it->second, numThreads) << "Wrong count for key " << k; + } +} + +// --------------------------------------------------------------------------- +// Concurrent merges (vector): all elements are present after reduction +// --------------------------------------------------------------------------- + +TEST(GreedyReduceAlgorithm, ConcurrentMergesVector) +{ + using VecType = std::vector; + using ReducerType = itk::GreedyReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(VectorMerge); + + constexpr int numThreads = 40; + reducer->SetNumberOfWorkUnits(numThreads); + + std::vector threads; + threads.reserve(numThreads); + + for (int i = 0; i < numThreads; ++i) + { + threads.emplace_back([&, i]() { + VecType local = { i }; + reducer->Merge(std::move(local)); + }); + } + + for (auto & t : threads) + { + t.join(); + } + + const VecType & result = reducer->GetResult(); + ASSERT_EQ(static_cast(result.size()), numThreads); + + // All thread indices 0..numThreads-1 must appear exactly once + VecType sorted = result; + std::sort(sorted.begin(), sorted.end()); + for (int i = 0; i < numThreads; ++i) + { + EXPECT_EQ(sorted[i], i); + } +} + +// --------------------------------------------------------------------------- +// NumberOfWorkUnits is stored and retrieved correctly +// --------------------------------------------------------------------------- + +TEST(GreedyReduceAlgorithm, NumberOfWorkUnits) +{ + using ReducerType = itk::GreedyReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + + EXPECT_EQ(reducer->GetNumberOfWorkUnits(), itk::SizeValueType{ 0 }); + + reducer->SetNumberOfWorkUnits(8); + EXPECT_EQ(reducer->GetNumberOfWorkUnits(), itk::SizeValueType{ 8 }); + + // Clear resets work unit count to zero + reducer->Clear(); + EXPECT_EQ(reducer->GetNumberOfWorkUnits(), itk::SizeValueType{ 0 }); +} + +// --------------------------------------------------------------------------- +// MergeFunction getter returns the stored function +// --------------------------------------------------------------------------- + +TEST(GreedyReduceAlgorithm, MergeFunctionRoundTrip) +{ + using ReducerType = itk::GreedyReduceAlgorithm; + auto reducer = ReducerType::New(); + + EXPECT_FALSE(reducer->GetMergeFunction()); + + reducer->SetMergeFunction(IntMerge); + EXPECT_TRUE(reducer->GetMergeFunction()); +} diff --git a/Modules/Core/Common/test/itkLinearReduceAlgorithmGTest.cxx b/Modules/Core/Common/test/itkLinearReduceAlgorithmGTest.cxx new file mode 100644 index 00000000000..b33ae442236 --- /dev/null +++ b/Modules/Core/Common/test/itkLinearReduceAlgorithmGTest.cxx @@ -0,0 +1,366 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +#include "itkLinearReduceAlgorithm.h" +#include "itkGTest.h" + +#include +#include +#include + +namespace +{ + +void +IntMerge(int & target, int & source) +{ + target += source; +} + +void +FloatMerge(float & target, float & source) +{ + target += source; +} + +void +MapMerge(std::map & target, std::map & source) +{ + for (auto & [key, val] : source) + { + target[key] += val; + } +} + +} // namespace + +// --------------------------------------------------------------------------- +// Basic object methods +// --------------------------------------------------------------------------- + +TEST(LinearReduceAlgorithm, BasicObjectMethods) +{ + using ReducerType = itk::LinearReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + reducer->SetNumberOfWorkUnits(4); + + ITK_GTEST_EXERCISE_BASIC_OBJECT_METHODS(reducer, LinearReduceAlgorithm, ReduceAlgorithm); +} + +// --------------------------------------------------------------------------- +// GetResult before any Merge returns default-constructed value +// --------------------------------------------------------------------------- + +TEST(LinearReduceAlgorithm, GetResultBeforeMerge) +{ + using ReducerType = itk::LinearReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + reducer->SetNumberOfWorkUnits(4); + + EXPECT_EQ(reducer->GetResult(), int{}); +} + +// --------------------------------------------------------------------------- +// Single chunk (N=1): no merge needed, value is the result +// --------------------------------------------------------------------------- + +TEST(LinearReduceAlgorithm, SingleChunk) +{ + using ReducerType = itk::LinearReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + reducer->SetNumberOfWorkUnits(1); + + int v = 42; + reducer->Merge(0, std::move(v)); + + EXPECT_EQ(reducer->GetResult(), 42); +} + +// --------------------------------------------------------------------------- +// Chunks submitted in ascending order +// --------------------------------------------------------------------------- + +TEST(LinearReduceAlgorithm, ChunksInOrder) +{ + using ReducerType = itk::LinearReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + reducer->SetNumberOfWorkUnits(4); + + for (int i = 0; i < 4; ++i) + { + reducer->Merge(static_cast(i), std::move(i)); + } + // 0+1+2+3 = 6 + EXPECT_EQ(reducer->GetResult(), 6); +} + +// --------------------------------------------------------------------------- +// Chunks submitted in reverse order — same result (deterministic) +// --------------------------------------------------------------------------- + +TEST(LinearReduceAlgorithm, ChunksReversed) +{ + using ReducerType = itk::LinearReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + reducer->SetNumberOfWorkUnits(4); + + for (int i = 3; i >= 0; --i) + { + reducer->Merge(static_cast(i), std::move(i)); + } + EXPECT_EQ(reducer->GetResult(), 6); +} + +// --------------------------------------------------------------------------- +// Chunks submitted out of order +// --------------------------------------------------------------------------- + +TEST(LinearReduceAlgorithm, ChunksOutOfOrder) +{ + using ReducerType = itk::LinearReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + reducer->SetNumberOfWorkUnits(5); + + for (int i : { 2, 0, 4, 1, 3 }) + { + reducer->Merge(static_cast(i), std::move(i)); + } + // 0+1+2+3+4 = 10 + EXPECT_EQ(reducer->GetResult(), 10); +} + +// --------------------------------------------------------------------------- +// Arbitrary N (no power-of-two restriction) +// --------------------------------------------------------------------------- + +TEST(LinearReduceAlgorithm, ArbitrarySizes) +{ + for (int n : { 1, 2, 3, 5, 6, 7, 9, 10, 13, 15 }) + { + using ReducerType = itk::LinearReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + reducer->SetNumberOfWorkUnits(static_cast(n)); + + const int expected = n * (n - 1) / 2; + for (int i = 0; i < n; ++i) + { + reducer->Merge(static_cast(i), std::move(i)); + } + EXPECT_EQ(reducer->GetResult(), expected) << "N=" << n; + } +} + +// --------------------------------------------------------------------------- +// Clear and reuse with the same N +// --------------------------------------------------------------------------- + +TEST(LinearReduceAlgorithm, ClearAndReuse) +{ + using ReducerType = itk::LinearReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + reducer->SetNumberOfWorkUnits(4); + + for (int i = 0; i < 4; ++i) + { + reducer->Merge(static_cast(i), std::move(i)); + } + EXPECT_EQ(reducer->GetResult(), 6); + + reducer->Clear(); + EXPECT_EQ(reducer->GetResult(), int{}); + + for (int i = 0; i < 4; ++i) + { + int v = i * 10; + reducer->Merge(static_cast(i), std::move(v)); + } + // 0+10+20+30 = 60 + EXPECT_EQ(reducer->GetResult(), 60); +} + +// --------------------------------------------------------------------------- +// Calling Merge(T&&) without chunk ID must throw +// --------------------------------------------------------------------------- + +TEST(LinearReduceAlgorithm, MergeWithoutIdThrows) +{ + using ReducerType = itk::LinearReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + reducer->SetNumberOfWorkUnits(2); + + int v = 5; + EXPECT_THROW(reducer->Merge(std::move(v)), itk::ExceptionObject); +} + +// --------------------------------------------------------------------------- +// Merge called before SetNumberOfWorkUnits must throw +// --------------------------------------------------------------------------- + +TEST(LinearReduceAlgorithm, MergeBeforeSetNumberOfWorkUnitsThrows) +{ + using ReducerType = itk::LinearReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + + int v = 5; + EXPECT_THROW(reducer->Merge(0, std::move(v)), itk::ExceptionObject); +} + +// --------------------------------------------------------------------------- +// chunkId out of range must throw +// --------------------------------------------------------------------------- + +TEST(LinearReduceAlgorithm, ChunkIdOutOfRangeThrows) +{ + using ReducerType = itk::LinearReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + reducer->SetNumberOfWorkUnits(2); + + int v = 5; + EXPECT_THROW(reducer->Merge(2, std::move(v)), itk::ExceptionObject); +} + +// --------------------------------------------------------------------------- +// Concurrent merges: N threads each contribute a known integer. +// Result must equal the expected sum, regardless of arrival order. +// --------------------------------------------------------------------------- + +TEST(LinearReduceAlgorithm, ConcurrentMergesInt) +{ + constexpr int numChunks = 16; + + using ReducerType = itk::LinearReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + reducer->SetNumberOfWorkUnits(numChunks); + + std::vector threads; + threads.reserve(numChunks); + for (int i = 0; i < numChunks; ++i) + { + threads.emplace_back([&reducer, i]() { + int v = i; + reducer->Merge(static_cast(i), std::move(v)); + }); + } + for (auto & t : threads) + { + t.join(); + } + + const int expected = numChunks * (numChunks - 1) / 2; + EXPECT_EQ(reducer->GetResult(), expected); +} + +// --------------------------------------------------------------------------- +// Determinism: floating-point sum is bit-identical across many repeated runs +// with concurrent threads, since the merge order is fixed by chunk ID. +// --------------------------------------------------------------------------- + +TEST(LinearReduceAlgorithm, DeterministicFloat) +{ + constexpr int numChunks = 8; + const std::vector values = { 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f }; + + using ReducerType = itk::LinearReduceAlgorithm; + + // Compute the reference result from a single-threaded run. + float referenceResult{}; + { + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(FloatMerge); + reducer->SetNumberOfWorkUnits(numChunks); + for (int i = 0; i < numChunks; ++i) + { + float v = values[static_cast(i)]; + reducer->Merge(static_cast(i), std::move(v)); + } + referenceResult = reducer->GetResult(); + } + + // Repeat with concurrent threads many times and verify exact equality. + constexpr int numRuns = 20; + for (int run = 0; run < numRuns; ++run) + { + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(FloatMerge); + reducer->SetNumberOfWorkUnits(numChunks); + + std::vector threads; + threads.reserve(numChunks); + for (int i = 0; i < numChunks; ++i) + { + threads.emplace_back([&reducer, &values, i]() { + float v = values[static_cast(i)]; + reducer->Merge(static_cast(i), std::move(v)); + }); + } + for (auto & t : threads) + { + t.join(); + } + + EXPECT_EQ(reducer->GetResult(), referenceResult) << "Run " << run << " produced a different result"; + } +} + +// --------------------------------------------------------------------------- +// Map reduction: each chunk contributes unique keys; verify all keys present +// --------------------------------------------------------------------------- + +TEST(LinearReduceAlgorithm, ConcurrentMergesMap) +{ + constexpr int numChunks = 12; + + using MapType = std::map; + using ReducerType = itk::LinearReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(MapMerge); + reducer->SetNumberOfWorkUnits(numChunks); + + std::vector threads; + threads.reserve(numChunks); + for (int i = 0; i < numChunks; ++i) + { + threads.emplace_back([&reducer, i]() { + MapType m{ { i, i * i } }; + reducer->Merge(static_cast(i), std::move(m)); + }); + } + for (auto & t : threads) + { + t.join(); + } + + const MapType & result = reducer->GetResult(); + ASSERT_EQ(result.size(), static_cast(numChunks)); + for (int i = 0; i < numChunks; ++i) + { + EXPECT_EQ(result.at(i), i * i); + } +} diff --git a/Modules/Core/Common/test/itkTreeReduceAlgorithmGTest.cxx b/Modules/Core/Common/test/itkTreeReduceAlgorithmGTest.cxx new file mode 100644 index 00000000000..d862d252746 --- /dev/null +++ b/Modules/Core/Common/test/itkTreeReduceAlgorithmGTest.cxx @@ -0,0 +1,381 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +#include "itkTreeReduceAlgorithm.h" +#include "itkGTest.h" + +#include +#include +#include +#include +#include + +namespace +{ + +void +IntMerge(int & target, int & source) +{ + target += source; +} + +void +FloatMerge(float & target, float & source) +{ + target += source; +} + +void +MapMerge(std::map & target, std::map & source) +{ + for (auto & [key, val] : source) + { + target[key] += val; + } +} + +} // namespace + +// --------------------------------------------------------------------------- +// Basic object methods +// --------------------------------------------------------------------------- + +TEST(TreeReduceAlgorithm, BasicObjectMethods) +{ + using ReducerType = itk::TreeReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + reducer->SetNumberOfWorkUnits(4); + + ITK_GTEST_EXERCISE_BASIC_OBJECT_METHODS(reducer, TreeReduceAlgorithm, ReduceAlgorithm); +} + +// --------------------------------------------------------------------------- +// GetResult before any Merge returns default-constructed value +// --------------------------------------------------------------------------- + +TEST(TreeReduceAlgorithm, GetResultBeforeMerge) +{ + using ReducerType = itk::TreeReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + reducer->SetNumberOfWorkUnits(4); + + EXPECT_EQ(reducer->GetResult(), int{}); +} + +// --------------------------------------------------------------------------- +// Single chunk (N=1): no merge needed, value is the result +// --------------------------------------------------------------------------- + +TEST(TreeReduceAlgorithm, SingleChunk) +{ + using ReducerType = itk::TreeReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + reducer->SetNumberOfWorkUnits(1); + + int v = 42; + reducer->Merge(0, std::move(v)); + + EXPECT_EQ(reducer->GetResult(), 42); +} + +// --------------------------------------------------------------------------- +// Two chunks: N=2 (already a power of two) +// --------------------------------------------------------------------------- + +TEST(TreeReduceAlgorithm, TwoChunks) +{ + using ReducerType = itk::TreeReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + reducer->SetNumberOfWorkUnits(2); + + int v0 = 10, v1 = 32; + reducer->Merge(0, std::move(v0)); + reducer->Merge(1, std::move(v1)); + + EXPECT_EQ(reducer->GetResult(), 42); +} + +// --------------------------------------------------------------------------- +// Four chunks (power of two): chunks submitted in ascending order +// --------------------------------------------------------------------------- + +TEST(TreeReduceAlgorithm, FourChunksInOrder) +{ + using ReducerType = itk::TreeReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + reducer->SetNumberOfWorkUnits(4); + + for (int i = 0; i < 4; ++i) + { + reducer->Merge(static_cast(i), std::move(i)); + } + // 0+1+2+3 = 6 + EXPECT_EQ(reducer->GetResult(), 6); +} + +// --------------------------------------------------------------------------- +// Four chunks: submitted in reverse order — same result (deterministic) +// --------------------------------------------------------------------------- + +TEST(TreeReduceAlgorithm, FourChunksReversed) +{ + using ReducerType = itk::TreeReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + reducer->SetNumberOfWorkUnits(4); + + for (int i = 3; i >= 0; --i) + { + reducer->Merge(static_cast(i), std::move(i)); + } + EXPECT_EQ(reducer->GetResult(), 6); +} + +// --------------------------------------------------------------------------- +// Non-power-of-two: N=5 (padded to 8) +// --------------------------------------------------------------------------- + +TEST(TreeReduceAlgorithm, NonPowerOfTwo_N5) +{ + using ReducerType = itk::TreeReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + reducer->SetNumberOfWorkUnits(5); + + for (int i = 0; i < 5; ++i) + { + reducer->Merge(static_cast(i), std::move(i)); + } + // 0+1+2+3+4 = 10 + EXPECT_EQ(reducer->GetResult(), 10); +} + +// --------------------------------------------------------------------------- +// Various non-power-of-two sizes +// --------------------------------------------------------------------------- + +TEST(TreeReduceAlgorithm, NonPowerOfTwo_Various) +{ + for (int n : { 3, 5, 6, 7, 9, 10, 13, 15 }) + { + using ReducerType = itk::TreeReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + reducer->SetNumberOfWorkUnits(static_cast(n)); + + const int expected = n * (n - 1) / 2; + for (int i = 0; i < n; ++i) + { + reducer->Merge(static_cast(i), std::move(i)); + } + EXPECT_EQ(reducer->GetResult(), expected) << "N=" << n; + } +} + +// --------------------------------------------------------------------------- +// Clear and reuse with the same N +// --------------------------------------------------------------------------- + +TEST(TreeReduceAlgorithm, ClearAndReuse) +{ + using ReducerType = itk::TreeReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + reducer->SetNumberOfWorkUnits(4); + + for (int i = 0; i < 4; ++i) + { + reducer->Merge(static_cast(i), std::move(i)); + } + EXPECT_EQ(reducer->GetResult(), 6); + + reducer->Clear(); + EXPECT_EQ(reducer->GetResult(), int{}); + + for (int i = 0; i < 4; ++i) + { + int v = i * 10; + reducer->Merge(static_cast(i), std::move(v)); + } + // 0+10+20+30 = 60 + EXPECT_EQ(reducer->GetResult(), 60); +} + +// --------------------------------------------------------------------------- +// Calling Merge(T&&) without chunk ID must throw +// --------------------------------------------------------------------------- + +TEST(TreeReduceAlgorithm, MergeWithoutIdThrows) +{ + using ReducerType = itk::TreeReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + reducer->SetNumberOfWorkUnits(2); + + int v = 5; + EXPECT_THROW(reducer->Merge(std::move(v)), itk::ExceptionObject); +} + +// --------------------------------------------------------------------------- +// Concurrent merges: N threads each contribute a known integer. +// Result must equal the expected sum. +// --------------------------------------------------------------------------- + +TEST(TreeReduceAlgorithm, ConcurrentMergesInt) +{ + constexpr int numChunks = 16; + + using ReducerType = itk::TreeReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + reducer->SetNumberOfWorkUnits(numChunks); + + std::vector threads; + threads.reserve(numChunks); + for (int i = 0; i < numChunks; ++i) + { + threads.emplace_back([&reducer, i]() { + int v = i; + reducer->Merge(static_cast(i), std::move(v)); + }); + } + for (auto & t : threads) + { + t.join(); + } + + const int expected = numChunks * (numChunks - 1) / 2; + EXPECT_EQ(reducer->GetResult(), expected); +} + +// --------------------------------------------------------------------------- +// Determinism: floating-point sum is bit-identical across many repeated runs +// with concurrent threads, regardless of thread scheduling. +// --------------------------------------------------------------------------- + +TEST(TreeReduceAlgorithm, DeterministicFloat) +{ + constexpr int numChunks = 8; + const std::vector values = { 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f }; + + using ReducerType = itk::TreeReduceAlgorithm; + + // Compute the reference result from the first run. + float referenceResult{}; + { + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(FloatMerge); + reducer->SetNumberOfWorkUnits(numChunks); + for (int i = 0; i < numChunks; ++i) + { + float v = values[static_cast(i)]; + reducer->Merge(static_cast(i), std::move(v)); + } + referenceResult = reducer->GetResult(); + } + + // Repeat with concurrent threads many times and verify exact equality. + constexpr int numRuns = 20; + for (int run = 0; run < numRuns; ++run) + { + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(FloatMerge); + reducer->SetNumberOfWorkUnits(numChunks); + + std::vector threads; + threads.reserve(numChunks); + for (int i = 0; i < numChunks; ++i) + { + threads.emplace_back([&reducer, &values, i]() { + float v = values[static_cast(i)]; + reducer->Merge(static_cast(i), std::move(v)); + }); + } + for (auto & t : threads) + { + t.join(); + } + + EXPECT_EQ(reducer->GetResult(), referenceResult) << "Run " << run << " produced a different result"; + } +} + +// --------------------------------------------------------------------------- +// Map reduction: each chunk contributes unique keys; verify all keys present +// --------------------------------------------------------------------------- + +TEST(TreeReduceAlgorithm, ConcurrentMergesMap) +{ + constexpr int numChunks = 12; + + using MapType = std::map; + using ReducerType = itk::TreeReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(MapMerge); + reducer->SetNumberOfWorkUnits(numChunks); + + std::vector threads; + threads.reserve(numChunks); + for (int i = 0; i < numChunks; ++i) + { + threads.emplace_back([&reducer, i]() { + MapType local; + local[i] = i * 10; + reducer->Merge(static_cast(i), std::move(local)); + }); + } + for (auto & t : threads) + { + t.join(); + } + + const MapType & result = reducer->GetResult(); + ASSERT_EQ(static_cast(result.size()), numChunks); + for (int i = 0; i < numChunks; ++i) + { + auto it = result.find(i); + ASSERT_NE(it, result.end()) << "Key " << i << " missing"; + EXPECT_EQ(it->second, i * 10) << "Wrong value for key " << i; + } +} + +// --------------------------------------------------------------------------- +// NumberOfWorkUnits getter reflects the set value +// --------------------------------------------------------------------------- + +TEST(TreeReduceAlgorithm, NumberOfWorkUnits) +{ + using ReducerType = itk::TreeReduceAlgorithm; + auto reducer = ReducerType::New(); + reducer->SetMergeFunction(IntMerge); + + EXPECT_EQ(reducer->GetNumberOfWorkUnits(), itk::SizeValueType{ 0 }); + + reducer->SetNumberOfWorkUnits(8); + EXPECT_EQ(reducer->GetNumberOfWorkUnits(), itk::SizeValueType{ 8 }); + + // Clear does NOT reset the work unit count. + reducer->Merge(0, std::move(int{ 1 })); + reducer->Clear(); + EXPECT_EQ(reducer->GetNumberOfWorkUnits(), itk::SizeValueType{ 8 }); +} diff --git a/Modules/Filtering/ImageStatistics/include/itkLabelOverlapMeasuresImageFilter.h b/Modules/Filtering/ImageStatistics/include/itkLabelOverlapMeasuresImageFilter.h index 1e893152b11..cb49636961e 100644 --- a/Modules/Filtering/ImageStatistics/include/itkLabelOverlapMeasuresImageFilter.h +++ b/Modules/Filtering/ImageStatistics/include/itkLabelOverlapMeasuresImageFilter.h @@ -21,7 +21,7 @@ #include "itkImageSink.h" #include "itkNumericTraits.h" #include "itkLabelOverlapLabelSetMeasures.h" -#include +#include "itkGreedyReduceAlgorithm.h" #include namespace itk @@ -235,6 +235,9 @@ class ITK_TEMPLATE_EXPORT LabelOverlapMeasuresImageFilter : public ImageSink::Pointer m_Reducer{}; }; // end of class } // end namespace itk diff --git a/Modules/Filtering/ImageStatistics/include/itkLabelOverlapMeasuresImageFilter.hxx b/Modules/Filtering/ImageStatistics/include/itkLabelOverlapMeasuresImageFilter.hxx index a565f4d2dc7..4ab3adec56c 100644 --- a/Modules/Filtering/ImageStatistics/include/itkLabelOverlapMeasuresImageFilter.hxx +++ b/Modules/Filtering/ImageStatistics/include/itkLabelOverlapMeasuresImageFilter.hxx @@ -33,6 +33,9 @@ LabelOverlapMeasuresImageFilter::LabelOverlapMeasuresImageFilter() // This filter requires two input images this->SetNumberOfRequiredInputs(2); + + m_Reducer = GreedyReduceAlgorithm::New(); + m_Reducer->SetMergeFunction([this](MapType & target, MapType & source) { this->MergeMap(target, source); }); } template @@ -43,6 +46,17 @@ LabelOverlapMeasuresImageFilter::BeforeStreamedGenerateData() // Initialize the final map this->m_LabelSetMeasures.clear(); + m_Reducer->Clear(); +} + +template +void +LabelOverlapMeasuresImageFilter::AfterStreamedGenerateData() +{ + Superclass::AfterStreamedGenerateData(); + + // Retrieve the merged per-label measures accumulated across all threads. + this->m_LabelSetMeasures = m_Reducer->GetResult(); } template @@ -116,28 +130,8 @@ LabelOverlapMeasuresImageFilter::ThreadedStreamedGenerateData(const } - // Merge localStatistics and m_LabelSetMeasures concurrently safe in a - // local copy, this thread may do multiple merges. - while (true) - { - MapType tomerge{}; - { - const std::lock_guard lockGuard(m_Mutex); - - if (m_LabelSetMeasures.empty()) - { - swap(m_LabelSetMeasures, localStatistics); - break; - } - - // Move the data of the output map to the local `tomerge` and clear the output map. - swap(m_LabelSetMeasures, tomerge); - - } // release lock, allow other threads to merge data - - // Merge tomerge into localStatistics, locally - MergeMap(localStatistics, tomerge); - } + // Merge localStatistics into the accumulated per-label measures. + m_Reducer->Merge(std::move(localStatistics)); } template diff --git a/Modules/Filtering/ImageStatistics/include/itkLabelStatisticsImageFilter.h b/Modules/Filtering/ImageStatistics/include/itkLabelStatisticsImageFilter.h index e0b213a4502..3816456e479 100644 --- a/Modules/Filtering/ImageStatistics/include/itkLabelStatisticsImageFilter.h +++ b/Modules/Filtering/ImageStatistics/include/itkLabelStatisticsImageFilter.h @@ -23,7 +23,7 @@ #include "itkSimpleDataObjectDecorator.h" #include "itkHistogram.h" #include "itkPrintHelper.h" -#include +#include "itkGreedyReduceAlgorithm.h" #include #include @@ -380,6 +380,7 @@ class ITK_TEMPLATE_EXPORT LabelStatisticsImageFilter : public ImageSinkAllocateOutputs(); m_LabelStatistics.clear(); + m_Reducer->Clear(); } /** Do final mean and variance computation from data accumulated in threads. @@ -404,7 +405,7 @@ class ITK_TEMPLATE_EXPORT LabelStatisticsImageFilter : public ImageSink::Pointer m_Reducer{}; }; // end of class } // end namespace itk diff --git a/Modules/Filtering/ImageStatistics/include/itkLabelStatisticsImageFilter.hxx b/Modules/Filtering/ImageStatistics/include/itkLabelStatisticsImageFilter.hxx index dfdef2285df..654f4c80653 100644 --- a/Modules/Filtering/ImageStatistics/include/itkLabelStatisticsImageFilter.hxx +++ b/Modules/Filtering/ImageStatistics/include/itkLabelStatisticsImageFilter.hxx @@ -34,6 +34,9 @@ LabelStatisticsImageFilter::LabelStatisticsImageFilter Self::AddRequiredInputName("LabelInput"); m_NumBins.SetSize(1); m_NumBins[0] = 256; + + m_Reducer = GreedyReduceAlgorithm::New(); + m_Reducer->SetMergeFunction([this](MapType & target, MapType & source) { this->MergeMap(target, source); }); } template @@ -112,6 +115,9 @@ LabelStatisticsImageFilter::AfterStreamedGenerateData( { Superclass::AfterStreamedGenerateData(); + // Retrieve the merged per-label statistics accumulated across all threads. + m_LabelStatistics = m_Reducer->GetResult(); + // compute the remainder of the statistics for (auto & mapValue : m_LabelStatistics) { @@ -241,28 +247,9 @@ LabelStatisticsImageFilter::ThreadedStreamedGenerateDa } - // Merge localStatistics and m_LabelStatistics concurrently safe in a - // local copy, this thread may do multiple merges. - while (true) - { - MapType tomerge{}; - { - const std::lock_guard lockGuard(m_Mutex); - - if (m_LabelStatistics.empty()) - { - swap(m_LabelStatistics, localStatistics); - break; - } - - // Move the data of the output map to the local `tomerge` and clear the output map. - swap(m_LabelStatistics, tomerge); - - } // release lock, allow other threads to merge data - - // Merge tomerge into localStatistics, locally - MergeMap(localStatistics, tomerge); - } + // Merge localStatistics into the accumulated per-label statistics, + // deterministically ordered is not required here. + m_Reducer->Merge(std::move(localStatistics)); } template diff --git a/Modules/Numerics/Statistics/include/itkImageToHistogramFilter.h b/Modules/Numerics/Statistics/include/itkImageToHistogramFilter.h index 535e7e3f995..125b12d272e 100644 --- a/Modules/Numerics/Statistics/include/itkImageToHistogramFilter.h +++ b/Modules/Numerics/Statistics/include/itkImageToHistogramFilter.h @@ -24,6 +24,7 @@ #include "itkImageSink.h" #include "itkSimpleDataObjectDecorator.h" #include "itkProgressReporter.h" +#include "itkGreedyReduceAlgorithm.h" namespace itk::Statistics { @@ -176,7 +177,7 @@ class ITK_TEMPLATE_EXPORT ImageToHistogramFilter : public ImageSink std::mutex m_Mutex{}; - HistogramPointer m_MergeHistogram{}; + typename GreedyReduceAlgorithm::Pointer m_Reducer{}; HistogramMeasurementVectorType m_Minimum{}; HistogramMeasurementVectorType m_Maximum{}; diff --git a/Modules/Numerics/Statistics/include/itkImageToHistogramFilter.hxx b/Modules/Numerics/Statistics/include/itkImageToHistogramFilter.hxx index ee936927ee4..7c810c39288 100644 --- a/Modules/Numerics/Statistics/include/itkImageToHistogramFilter.hxx +++ b/Modules/Numerics/Statistics/include/itkImageToHistogramFilter.hxx @@ -42,6 +42,22 @@ ImageToHistogramFilter::ImageToHistogramFilter() { this->Self::SetAutoMinimumMaximum(true); } + + m_Reducer = GreedyReduceAlgorithm::New(); + m_Reducer->SetMergeFunction([](HistogramPointer & target, HistogramPointer & source) { + using HistogramIterator = typename HistogramType::ConstIterator; + + HistogramIterator hit = source->Begin(); + const HistogramIterator end = source->End(); + + typename HistogramType::IndexType index; + while (hit != end) + { + target->GetIndex(hit.GetMeasurementVector(), index); + target->IncreaseFrequencyOfIndex(index, hit.GetFrequency()); + ++hit; + } + }); } template @@ -119,7 +135,7 @@ ImageToHistogramFilter::InitializeOutputHistogram() m_Minimum.Fill(NumericTraits::max()); m_Maximum.Fill(NumericTraits::NonpositiveMin()); - m_MergeHistogram = nullptr; + m_Reducer->Clear(); HistogramType * outputHistogram = this->GetOutput(); outputHistogram->SetClipBinsAtEnds(true); @@ -186,8 +202,8 @@ ImageToHistogramFilter::AfterStreamedGenerateData() Superclass::AfterStreamedGenerateData(); HistogramType * outputHistogram = this->GetOutput(); - outputHistogram->Graft(m_MergeHistogram); - m_MergeHistogram = nullptr; + outputHistogram->Graft(m_Reducer->GetResult()); + m_Reducer->Clear(); } @@ -257,39 +273,7 @@ template void ImageToHistogramFilter::ThreadedMergeHistogram(HistogramPointer && histogram) { - while (true) - { - HistogramPointer tomergeHistogram{}; - { - const std::lock_guard lockGuard(m_Mutex); - - if (m_MergeHistogram.IsNull()) - { - m_MergeHistogram = std::move(histogram); - return; - } - - // merge/reduce the local results with current values in m_MergeHistogram - - // take ownership locally - swap(m_MergeHistogram, tomergeHistogram); - - } // release lock, allow other threads to merge data - - using HistogramIterator = typename HistogramType::ConstIterator; - - HistogramIterator hit = tomergeHistogram->Begin(); - const HistogramIterator end = tomergeHistogram->End(); - - typename HistogramType::IndexType index; - - while (hit != end) - { - histogram->GetIndex(hit.GetMeasurementVector(), index); - histogram->IncreaseFrequencyOfIndex(index, hit.GetFrequency()); - ++hit; - } - } + m_Reducer->Merge(std::move(histogram)); } template