Skip to content

Add GetRangeOfSamples(ThreadIdType) to AdvancedImageToImageMetric, remove duplicate code #1441

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
N-Dekker merged 3 commits into from
May 8, 2026
Merged

Add GetRangeOfSamples(ThreadIdType) to AdvancedImageToImageMetric, remove duplicate code #1441

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
6e437a5
STYLE: Add GetRangeOfSamples(ThreadIdType) to AdvancedImageToImageMetric
N-Dekker May 7, 2026
d1a3f1f
STYLE: Use `rangeOfSamples = GetRangeOfSamples(threadId)` in PCAMetric
N-Dekker May 7, 2026
1f3cd70
STYLE: Use GetRangeOfSamples(ThreadIdType) in Metrics and CostFunctions
N-Dekker May 7, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
38 changes: 38 additions & 0 deletions Common/CostFunctions/itkAdvancedImageToImageMetric.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -386,6 +386,44 @@ class ITK_TEMPLATE_EXPORT AdvancedImageToImageMetric : public ImageToImageMetric
MovingImageLimiterOutputType m_MovingImageMinLimit{ 0 };
MovingImageLimiterOutputType m_MovingImageMaxLimit{ 1 };


/** Retrieves a subrange of the samples from its ImageSampler, for the specified work unit. */
auto
GetRangeOfSamples(const ThreadIdType workUnitID) const
{
const auto & sampleContainer = m_ImageSampler->GetOutput()->CastToSTLConstContainer();

const size_t sampleContainerSize{ sampleContainer.size() };
const size_t numberOfWorkUnits{ this->Superclass::GetNumberOfWorkUnits() };
const auto numberOfSamplesPerWorkUnit =
static_cast<size_t>(std::ceil(static_cast<double>(sampleContainerSize) / static_cast<double>(numberOfWorkUnits)));

using IteratorType = decltype(sampleContainer.cbegin());

struct Range
{
IteratorType Begin{};
IteratorType End{};

auto
begin() const
{
return Begin;
}

auto
end() const
{
return End;
}
};

const auto beginOfSampleContainer = sampleContainer.cbegin();
return Range{ beginOfSampleContainer + std::min(numberOfSamplesPerWorkUnit * workUnitID, sampleContainerSize),
beginOfSampleContainer +
std::min(numberOfSamplesPerWorkUnit * (workUnitID + 1), sampleContainerSize) };
}

/** Multi-threaded metric computation. */

/** Multi-threaded version of GetValue(). */
Expand Down
22 changes: 3 additions & 19 deletions Common/CostFunctions/itkParzenWindowHistogramImageToImageMetric.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1005,30 +1005,14 @@ ParzenWindowHistogramImageToImageMetric<TFixedImage, TMovingImage>::ThreadedComp
JointPDFPointer & jointPDF = m_ParzenWindowHistogramGetValueAndDerivativePerThreadVariables[threadId].st_JointPDF;
jointPDF->FillBuffer(PDFValueType{});

/** Get a handle to the sample container. */
ImageSampleContainerPointer sampleContainer = this->GetImageSampler()->GetOutput();
const size_t sampleContainerSize{ sampleContainer->size() };

/** Get the samples for this thread. */
const auto nrOfSamplesPerThreads = static_cast<unsigned long>(
std::ceil(static_cast<double>(sampleContainerSize) / static_cast<double>(Self::GetNumberOfWorkUnits())));

const auto pos_begin = std::min<size_t>(nrOfSamplesPerThreads * threadId, sampleContainerSize);
const auto pos_end = std::min<size_t>(nrOfSamplesPerThreads * (threadId + 1), sampleContainerSize);

/** Create iterator over the sample container. */
const auto beginOfSampleContainer = sampleContainer->cbegin();
const auto fbegin = beginOfSampleContainer + pos_begin;
const auto fend = beginOfSampleContainer + pos_end;

/** Create variables to store intermediate results. circumvent false sharing */
unsigned long numberOfPixelsCounted = 0;

/** Loop over sample container and compute contribution of each sample to pdfs. */
for (auto fiter = fbegin; fiter != fend; ++fiter)
for (const auto & sample : this->Superclass::GetRangeOfSamples(threadId))
{
/** Read fixed coordinates and initialize some variables. */
const FixedImagePointType & fixedPoint = fiter->m_ImageCoordinates;
const FixedImagePointType & fixedPoint = sample.m_ImageCoordinates;
RealType movingImageValue;

/** Transform point. */
Expand All @@ -1050,7 +1034,7 @@ ParzenWindowHistogramImageToImageMetric<TFixedImage, TMovingImage>::ThreadedComp
++numberOfPixelsCounted;

/** Get the fixed image value. */
auto fixedImageValue = static_cast<RealType>(fiter->m_ImageValue);
auto fixedImageValue = static_cast<RealType>(sample.m_ImageValue);

/** Make sure the values fall within the histogram range. */
fixedImageValue = this->GetFixedImageLimiter()->Evaluate(fixedImageValue);
Expand Down
22 changes: 3 additions & 19 deletions Components/Metrics/AdvancedKappaStatistic/itkAdvancedKappaStatisticImageToImageMetric.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -449,34 +449,18 @@ AdvancedKappaStatisticImageToImageMetric<TFixedImage, TMovingImage>::ThreadedGet
DerivativeType & vecSum1 = this->m_KappaGetValueAndDerivativePerThreadVariables[threadId].st_DerivativeSum1;
DerivativeType & vecSum2 = this->m_KappaGetValueAndDerivativePerThreadVariables[threadId].st_DerivativeSum2;

/** Get a handle to the sample container. */
ImageSampleContainerPointer sampleContainer = this->GetImageSampler()->GetOutput();
const size_t sampleContainerSize{ sampleContainer->size() };

/** Get the samples for this thread. */
const auto nrOfSamplesPerThreads = static_cast<unsigned long>(
std::ceil(static_cast<double>(sampleContainerSize) / static_cast<double>(Self::GetNumberOfWorkUnits())));

const auto pos_begin = std::min<size_t>(nrOfSamplesPerThreads * threadId, sampleContainerSize);
const auto pos_end = std::min<size_t>(nrOfSamplesPerThreads * (threadId + 1), sampleContainerSize);

/** Some variables. */
RealType movingImageValue;
std::size_t fixedForegroundArea = 0; // or unsigned long
std::size_t movingForegroundArea = 0;
std::size_t intersection = 0;
unsigned long numberOfPixelsCounted = 0;

/** Create iterator over the sample container. */
const auto beginOfSampleContainer = sampleContainer->cbegin();
const auto fbegin = beginOfSampleContainer + pos_begin;
const auto fend = beginOfSampleContainer + pos_end;

/** Loop over the fixed image to calculate the kappa statistic. */
for (auto fiter = fbegin; fiter != fend; ++fiter)
for (const auto & sample : this->Superclass::GetRangeOfSamples(threadId))
{
/** Read fixed coordinates. */
const FixedImagePointType & fixedPoint = fiter->m_ImageCoordinates;
const FixedImagePointType & fixedPoint = sample.m_ImageCoordinates;

/** Transform point. */
const MovingImagePointType mappedPoint = this->TransformPoint(fixedPoint);
Expand All @@ -500,7 +484,7 @@ AdvancedKappaStatisticImageToImageMetric<TFixedImage, TMovingImage>::ThreadedGet
++numberOfPixelsCounted;

/** Get the fixed image value. */
const auto fixedImageValue = static_cast<RealType>(fiter->m_ImageValue);
const auto fixedImageValue = static_cast<RealType>(sample.m_ImageValue);

#if 0
/** Get the TransformJacobian dT/dmu. */
Expand Down
22 changes: 3 additions & 19 deletions ...cs/AdvancedMattesMutualInformation/itkParzenWindowMutualInformationImageToImageMetric.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -433,27 +433,11 @@ ParzenWindowMutualInformationImageToImageMetric<TFixedImage, TMovingImage>::Thre
preconditioningDivisor = DerivativeType(this->GetNumberOfParameters(), 0.0);
}

/** Get a handle to the sample container. */
ImageSampleContainerPointer sampleContainer = this->GetImageSampler()->GetOutput();
const size_t sampleContainerSize{ sampleContainer->size() };

/** Get the samples for this thread. */
const auto nrOfSamplesPerThreads = static_cast<unsigned long>(
std::ceil(static_cast<double>(sampleContainerSize) / static_cast<double>(Self::GetNumberOfWorkUnits())));

const auto pos_begin = std::min<size_t>(nrOfSamplesPerThreads * threadId, sampleContainerSize);
const auto pos_end = std::min<size_t>(nrOfSamplesPerThreads * (threadId + 1), sampleContainerSize);

/** Create iterator over the sample container. */
const auto beginOfSampleContainer = sampleContainer->cbegin();
const auto fbegin = beginOfSampleContainer + pos_begin;
const auto fend = beginOfSampleContainer + pos_end;

/** Loop over sample container and compute contribution of each sample to pdfs. */
for (auto fiter = fbegin; fiter != fend; ++fiter)
for (const auto & sample : this->Superclass::GetRangeOfSamples(threadId))
{
/** Read fixed coordinates and create some variables. */
const FixedImagePointType & fixedPoint = fiter->m_ImageCoordinates;
const FixedImagePointType & fixedPoint = sample.m_ImageCoordinates;
RealType movingImageValue;
MovingImageDerivativeType movingImageDerivative;

Expand All @@ -475,7 +459,7 @@ ParzenWindowMutualInformationImageToImageMetric<TFixedImage, TMovingImage>::Thre
if (sampleOk)
{
/** Get the fixed image value. */
auto fixedImageValue = static_cast<RealType>(fiter->m_ImageValue);
auto fixedImageValue = static_cast<RealType>(sample.m_ImageValue);

/** Make sure the values fall within the histogram range. */
fixedImageValue = this->GetFixedImageLimiter()->Evaluate(fixedImageValue);
Expand Down
44 changes: 6 additions & 38 deletions Components/Metrics/AdvancedMeanSquares/itkAdvancedMeanSquaresImageToImageMetric.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -251,31 +251,15 @@ template <typename TFixedImage, typename TMovingImage>
void
AdvancedMeanSquaresImageToImageMetric<TFixedImage, TMovingImage>::ThreadedGetValue(ThreadIdType threadId) const
{
/** Get a handle to the sample container. */
ImageSampleContainerPointer sampleContainer = this->GetImageSampler()->GetOutput();
const size_t sampleContainerSize{ sampleContainer->size() };

/** Get the samples for this thread. */
const auto nrOfSamplesPerThreads = static_cast<unsigned long>(
std::ceil(static_cast<double>(sampleContainerSize) / static_cast<double>(Self::GetNumberOfWorkUnits())));

const auto pos_begin = std::min<size_t>(nrOfSamplesPerThreads * threadId, sampleContainerSize);
const auto pos_end = std::min<size_t>(nrOfSamplesPerThreads * (threadId + 1), sampleContainerSize);

/** Create iterator over the sample container. */
const auto beginOfSampleContainer = sampleContainer->cbegin();
const auto threader_fbegin = beginOfSampleContainer + pos_begin;
const auto threader_fend = beginOfSampleContainer + pos_end;

/** Create variables to store intermediate results. circumvent false sharing */
unsigned long numberOfPixelsCounted = 0;
MeasureType measure{};

/** Loop over the fixed image to calculate the mean squares. */
for (auto threader_fiter = threader_fbegin; threader_fiter != threader_fend; ++threader_fiter)
for (const auto & sample : this->Superclass::GetRangeOfSamples(threadId))
{
/** Read fixed coordinates and initialize some variables. */
const FixedImagePointType & fixedPoint = threader_fiter->m_ImageCoordinates;
const FixedImagePointType & fixedPoint = sample.m_ImageCoordinates;
RealType movingImageValue;

/** Transform point. */
Expand All @@ -297,7 +281,7 @@ AdvancedMeanSquaresImageToImageMetric<TFixedImage, TMovingImage>::ThreadedGetVal
++numberOfPixelsCounted;

/** Get the fixed image value. */
const auto fixedImageValue = static_cast<RealType>(threader_fiter->m_ImageValue);
const auto fixedImageValue = static_cast<RealType>(sample.m_ImageValue);

/** The difference squared. */
const RealType diff = movingImageValue - fixedImageValue;
Expand Down Expand Up @@ -547,31 +531,15 @@ AdvancedMeanSquaresImageToImageMetric<TFixedImage, TMovingImage>::ThreadedGetVal
*/
DerivativeType & derivative = Superclass::m_GetValueAndDerivativePerThreadVariables[threadId].st_Derivative;

/** Get a handle to the sample container. */
ImageSampleContainerPointer sampleContainer = this->GetImageSampler()->GetOutput();
const size_t sampleContainerSize{ sampleContainer->size() };

/** Get the samples for this thread. */
const auto nrOfSamplesPerThreads = static_cast<unsigned long>(
std::ceil(static_cast<double>(sampleContainerSize) / static_cast<double>(Self::GetNumberOfWorkUnits())));

const auto pos_begin = std::min<size_t>(nrOfSamplesPerThreads * threadId, sampleContainerSize);
const auto pos_end = std::min<size_t>(nrOfSamplesPerThreads * (threadId + 1), sampleContainerSize);

/** Create iterator over the sample container. */
const auto beginOfSampleContainer = sampleContainer->cbegin();
const auto threader_fbegin = beginOfSampleContainer + pos_begin;
const auto threader_fend = beginOfSampleContainer + pos_end;

/** Create variables to store intermediate results. circumvent false sharing */
unsigned long numberOfPixelsCounted = 0;
MeasureType measure{};

/** Loop over the fixed image to calculate the mean squares. */
for (auto threader_fiter = threader_fbegin; threader_fiter != threader_fend; ++threader_fiter)
for (const auto & sample : this->Superclass::GetRangeOfSamples(threadId))
{
/** Read fixed coordinates and initialize some variables. */
const FixedImagePointType & fixedPoint = threader_fiter->m_ImageCoordinates;
const FixedImagePointType & fixedPoint = sample.m_ImageCoordinates;
RealType movingImageValue;
MovingImageDerivativeType movingImageDerivative;

Expand All @@ -595,7 +563,7 @@ AdvancedMeanSquaresImageToImageMetric<TFixedImage, TMovingImage>::ThreadedGetVal
++numberOfPixelsCounted;

/** Get the fixed image value. */
const auto fixedImageValue = static_cast<RealType>(threader_fiter->m_ImageValue);
const auto fixedImageValue = static_cast<RealType>(sample.m_ImageValue);

#if 0
/** Get the TransformJacobian dT/dmu. */
Expand Down
22 changes: 3 additions & 19 deletions ...rics/AdvancedNormalizedCorrelation/itkAdvancedNormalizedCorrelationImageToImageMetric.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -474,22 +474,6 @@ AdvancedNormalizedCorrelationImageToImageMetric<TFixedImage, TMovingImage>::Thre
DerivativeType & derivativeM = this->m_CorrelationGetValueAndDerivativePerThreadVariables[threadId].st_DerivativeM;
DerivativeType & differential = this->m_CorrelationGetValueAndDerivativePerThreadVariables[threadId].st_Differential;

/** Get a handle to the sample container. */
ImageSampleContainerPointer sampleContainer = this->GetImageSampler()->GetOutput();
const size_t sampleContainerSize{ sampleContainer->size() };

/** Get the samples for this thread. */
const auto nrOfSamplesPerThreads = static_cast<unsigned long>(
std::ceil(static_cast<double>(sampleContainerSize) / static_cast<double>(Self::GetNumberOfWorkUnits())));

const auto pos_begin = std::min<size_t>(nrOfSamplesPerThreads * threadId, sampleContainerSize);
const auto pos_end = std::min<size_t>(nrOfSamplesPerThreads * (threadId + 1), sampleContainerSize);

/** Create iterator over the sample container. */
const auto beginOfSampleContainer = sampleContainer->cbegin();
const auto threader_fbegin = beginOfSampleContainer + pos_begin;
const auto threader_fend = beginOfSampleContainer + pos_end;

/** Create variables to store intermediate results. */
AccumulateType sff{};
AccumulateType smm{};
Expand All @@ -499,10 +483,10 @@ AdvancedNormalizedCorrelationImageToImageMetric<TFixedImage, TMovingImage>::Thre
unsigned long numberOfPixelsCounted = 0;

/** Loop over the fixed image to calculate the mean squares. */
for (auto threader_fiter = threader_fbegin; threader_fiter != threader_fend; ++threader_fiter)
for (const auto & sample : this->Superclass::GetRangeOfSamples(threadId))
{
/** Read fixed coordinates and initialize some variables. */
const FixedImagePointType & fixedPoint = threader_fiter->m_ImageCoordinates;
const FixedImagePointType & fixedPoint = sample.m_ImageCoordinates;
RealType movingImageValue;
MovingImageDerivativeType movingImageDerivative;

Expand All @@ -526,7 +510,7 @@ AdvancedNormalizedCorrelationImageToImageMetric<TFixedImage, TMovingImage>::Thre
++numberOfPixelsCounted;

/** Get the fixed image value. */
const auto fixedImageValue = static_cast<RealType>(threader_fiter->m_ImageValue);
const auto fixedImageValue = static_cast<RealType>(sample.m_ImageValue);

#if 0
/** Get the TransformJacobian dT/dmu. */
Expand Down
20 changes: 2 additions & 18 deletions Components/Metrics/BendingEnergyPenalty/itkTransformBendingEnergyPenaltyTerm.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -397,31 +397,15 @@ TransformBendingEnergyPenaltyTerm<TFixedImage, TScalarType>::ThreadedGetValueAnd
*/
DerivativeType & derivative = Superclass::m_GetValueAndDerivativePerThreadVariables[threadId].st_Derivative;

/** Get a handle to the sample container. */
ImageSampleContainerPointer sampleContainer = this->GetImageSampler()->GetOutput();
const size_t sampleContainerSize{ sampleContainer->size() };

/** Get the samples for this thread. */
const auto nrOfSamplesPerThreads = static_cast<unsigned long>(
std::ceil(static_cast<double>(sampleContainerSize) / static_cast<double>(Self::GetNumberOfWorkUnits())));

const auto pos_begin = std::min<size_t>(nrOfSamplesPerThreads * threadId, sampleContainerSize);
const auto pos_end = std::min<size_t>(nrOfSamplesPerThreads * (threadId + 1), sampleContainerSize);

/** Create iterator over the sample container. */
const auto beginOfSampleContainer = sampleContainer->cbegin();
const auto fbegin = beginOfSampleContainer + pos_begin;
const auto fend = beginOfSampleContainer + pos_end;

/** Create variables to store intermediate results. circumvent false sharing */
unsigned long numberOfPixelsCounted = 0;
MeasureType measure{};

/** Loop over the fixed image to calculate the penalty term and its derivative. */
for (auto fiter = fbegin; fiter != fend; ++fiter)
for (const auto & sample : this->Superclass::GetRangeOfSamples(threadId))
{
/** Read fixed coordinates and initialize some variables. */
const FixedImagePointType & fixedPoint = fiter->m_ImageCoordinates;
const FixedImagePointType & fixedPoint = sample.m_ImageCoordinates;

/** Although the mapped point is not needed to compute the penalty term,
* we compute in order to check if it maps inside the support region of
Expand Down
21 changes: 4 additions & 17 deletions Components/Metrics/PCAMetric/itkPCAMetric.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -644,30 +644,17 @@ template <typename TFixedImage, typename TMovingImage>
void
PCAMetric<TFixedImage, TMovingImage>::ThreadedGetSamples(ThreadIdType threadId)
{
/** Get a handle to the sample container. */
ImageSampleContainerPointer sampleContainer = this->GetImageSampler()->GetOutput();
const size_t sampleContainerSize{ sampleContainer->size() };

/** Get the samples for this thread. */
const auto nrOfSamplesPerThreads = static_cast<unsigned long>(
std::ceil(static_cast<double>(sampleContainerSize) / static_cast<double>(Self::GetNumberOfWorkUnits())));

const auto pos_begin = std::min<size_t>(nrOfSamplesPerThreads * threadId, sampleContainerSize);
const auto pos_end = std::min<size_t>(nrOfSamplesPerThreads * (threadId + 1), sampleContainerSize);

/** Create iterator over the sample container. */
const auto beginOfSampleContainer = sampleContainer->cbegin();
const auto threader_fbegin = beginOfSampleContainer + pos_begin;
const auto threader_fend = beginOfSampleContainer + pos_end;
const auto & rangeOfSamples = this->Superclass::GetRangeOfSamples(threadId);

std::vector<FixedImagePointType> SamplesOK;
MatrixType datablock(nrOfSamplesPerThreads, m_G);
MatrixType datablock(rangeOfSamples.end() - rangeOfSamples.begin(), m_G);

unsigned int pixelIndex = 0;
for (auto threader_fiter = threader_fbegin; threader_fiter != threader_fend; ++threader_fiter)
for (const auto & sample : rangeOfSamples)
{
/** Read fixed coordinates. */
FixedImagePointType fixedPoint = threader_fiter->m_ImageCoordinates;
FixedImagePointType fixedPoint = sample.m_ImageCoordinates;

/** Transform sampled point to voxel coordinates. */
auto voxelCoord =
Expand Down
Loading
Loading