Skip to content
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

ComputeImageExtremaFilter: Remove UseMask property and let it only compute min and max #1046

Merged
merged 2 commits into from
Feb 19, 2024
Merged

ComputeImageExtremaFilter: Remove UseMask property and let it only compute min and max #1046

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
554d37d
STYLE: Remove `UseMask` property from ComputeImageExtremaFilter
N-Dekker Feb 16, 2024
476969a
PERF: Let ComputeImageExtremaFilter only compute min and max
N-Dekker Feb 16, 2024
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
13 changes: 2 additions & 11 deletions Common/CostFunctions/itkAdvancedImageToImageMetric.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -192,12 +192,7 @@ AdvancedImageToImageMetric<TFixedImage, TMovingImage>::InitializeLimiters()

const auto computeFixedImageExtrema = ComputeImageExtremaFilter<FixedImageType>::New();
computeFixedImageExtrema->SetInput(this->GetFixedImage());
if (const auto * const fMask = this->GetFixedImageMask())
{
computeFixedImageExtrema->SetUseMask(true);
computeFixedImageExtrema->SetImageSpatialMask(fMask);
}

computeFixedImageExtrema->SetImageSpatialMask(this->GetFixedImageMask());
computeFixedImageExtrema->Update();

this->m_FixedImageTrueMax = computeFixedImageExtrema->GetMaximum();
Expand Down Expand Up @@ -228,11 +223,7 @@ AdvancedImageToImageMetric<TFixedImage, TMovingImage>::InitializeLimiters()

const auto computeMovingImageExtrema = ComputeImageExtremaFilter<MovingImageType>::New();
computeMovingImageExtrema->SetInput(this->GetMovingImage());
if (const auto * const mask = this->GetMovingImageMask())
{
computeMovingImageExtrema->SetUseMask(true);
computeMovingImageExtrema->SetImageSpatialMask(mask);
}
computeMovingImageExtrema->SetImageSpatialMask(this->GetMovingImageMask());
computeMovingImageExtrema->Update();

this->m_MovingImageTrueMax = computeMovingImageExtrema->GetMaximum();
Expand Down
4 changes: 0 additions & 4 deletions Common/GTesting/itkComputeImageExtremaFilterGTest.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -196,7 +196,6 @@ Expect_one_non_zero_pixel_value_masked_in(const typename TImage::SizeType & imag
const auto filter = CheckNew<FilterType>();
filter->SetInput(image);
filter->SetImageSpatialMask(maskSpatialObject);
filter->SetUseMask(true);
filter->Update();

EXPECT_EQ(filter->GetMinimum(), pixelValue);
Expand Down Expand Up @@ -247,7 +246,6 @@ Expect_one_positive_pixel_value_all_pixels_masked_in(const typename TImage::Size
const auto filter = CheckNew<FilterType>();
filter->SetInput(image);
filter->SetImageSpatialMask(maskSpatialObject);
filter->SetUseMask(true);
filter->Update();

EXPECT_EQ(filter->GetMinimum(), 0);
Expand Down Expand Up @@ -356,8 +354,6 @@ GTEST_TEST(ComputeImageExtremaFilter, MaskSpacingAndDirectionAffectResults)
elastix::DefaultConstruct<FilterType> filter{};
filter.SetInput(image);
filter.SetImageSpatialMask(maskSpatialObject);
filter.SetUseMask(true);

filter.Update();
return filter.GetMaximum();
};
Expand Down
75 changes: 32 additions & 43 deletions Common/itkComputeImageExtremaFilter.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,40 +18,26 @@
#ifndef itkComputeImageExtremaFilter_h
#define itkComputeImageExtremaFilter_h

#include "itkStatisticsImageFilter.h"
#include "itkImageSink.h"
#include "itkImageMaskSpatialObject.h"

namespace itk
{
/** \class ComputeImageExtremaFilter
* \brief Compute min. max, variance and mean of an Image.
*
* StatisticsImageFilter computes the minimum, maximum, sum, mean, variance
* sigma of an image. The filter needs all of its input image. It
* behaves as a filter with an input and output. Thus it can be inserted
* in a pipline with other filters and the statistics will only be
* recomputed if a downstream filter changes.
*
* The filter passes its input through unmodified. The filter is
* threaded. It computes statistics in each thread then combines them in
* its AfterThreadedGenerate method.
* \brief Compute minimum and maximum pixel value of an Image.
*
* \ingroup MathematicalStatisticsImageFilters
* \ingroup ITKImageStatistics
*
* \wiki
* \wikiexample{Statistics/StatisticsImageFilter,Compute min\, max\, variance and mean of an Image.}
* \endwiki
*/
template <typename TInputImage>
class ITK_TEMPLATE_EXPORT ComputeImageExtremaFilter : public StatisticsImageFilter<TInputImage>
class ITK_TEMPLATE_EXPORT ComputeImageExtremaFilter : public ImageSink<TInputImage>
{
public:
ITK_DISALLOW_COPY_AND_MOVE(ComputeImageExtremaFilter);

/** Standard Self typedef */
using Self = ComputeImageExtremaFilter;
using Superclass = StatisticsImageFilter<TInputImage>;
using Superclass = ImageSink<TInputImage>;
using Pointer = SmartPointer<Self>;
using ConstPointer = SmartPointer<const Self>;

Expand All @@ -64,55 +50,58 @@ class ITK_TEMPLATE_EXPORT ComputeImageExtremaFilter : public StatisticsImageFilt
/** Image related typedefs. */
using InputImagePointer = typename TInputImage::Pointer;

using typename Superclass::RegionType;
using typename Superclass::SizeType;
using typename Superclass::IndexType;
using typename Superclass::PixelType;
using PointType = typename TInputImage::PointType;
using Superclass::InputImageDimension;
using typename Superclass::InputImageRegionType;
using PixelType = typename Superclass::InputImagePixelType;

/** Image related typedefs. */
itkStaticConstMacro(ImageDimension, unsigned int, TInputImage::ImageDimension);

/** Type to use for computations. */
using typename Superclass::RealType;

itkSetMacro(UseMask, bool);

using ImageSpatialMaskType = ImageMaskSpatialObject<Self::ImageDimension>;
using ImageSpatialMaskPointer = typename ImageSpatialMaskType::Pointer;
using ImageSpatialMaskConstPointer = typename ImageSpatialMaskType::ConstPointer;
itkSetConstObjectMacro(ImageSpatialMask, ImageSpatialMaskType);
itkGetConstObjectMacro(ImageSpatialMask, ImageSpatialMaskType);

PixelType
GetMinimum() const
{
return m_ThreadMin;
}

PixelType
GetMaximum() const
{
return m_ThreadMax;
}

protected:
ComputeImageExtremaFilter() = default;
~ComputeImageExtremaFilter() override = default;

/** Initialize some accumulators before the threads run. */
/** Initialize minimum and maximum before the threads run. */
void
BeforeStreamedGenerateData() override;

/** Do final mean and variance computation from data accumulated in threads.
*/
void
AfterStreamedGenerateData() override;

/** Multi-thread version GenerateData. */
void
ThreadedStreamedGenerateData(const RegionType &) override;
virtual void
ThreadedGenerateDataImageSpatialMask(const RegionType &);
ThreadedStreamedGenerateData(const InputImageRegionType &) override;

private:
struct MinMaxResult
{
PixelType Min;
PixelType Max;
};

static MinMaxResult
RetrieveMinMax(const TInputImage &, const InputImageRegionType &, const ImageSpatialMaskType *, bool);

ImageSpatialMaskConstPointer m_ImageSpatialMask{};
bool m_UseMask{ false };
bool m_SameGeometry{ false };

CompensatedSummation<RealType> m_ThreadSum{ 1 };
CompensatedSummation<RealType> m_SumOfSquares{ 1 };
SizeValueType m_Count{ 1 };
PixelType m_ThreadMin{ 1 };
PixelType m_ThreadMax{ 1 };
PixelType m_ThreadMin{ 1 };
PixelType m_ThreadMax{ 1 };

std::mutex m_Mutex{};
}; // end of class
Expand Down
177 changes: 61 additions & 116 deletions Common/itkComputeImageExtremaFilter.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,9 +17,12 @@
*=========================================================================*/
#ifndef itkComputeImageExtremaFilter_hxx
#define itkComputeImageExtremaFilter_hxx

#include "itkComputeImageExtremaFilter.h"

#include <itkImageRegionConstIterator.h>
#include <itkImageScanlineIterator.h>

#include "elxMaskHasSameImageDomain.h"

namespace itk
Expand All @@ -29,144 +32,86 @@ template <typename TInputImage>
void
ComputeImageExtremaFilter<TInputImage>::BeforeStreamedGenerateData()
{
if (!this->m_UseMask)
{
Superclass::BeforeStreamedGenerateData();
}
else
{
// Resize the thread temporaries
m_Count = SizeValueType{};
m_SumOfSquares = RealType{};
m_ThreadSum = RealType{};
m_ThreadMin = NumericTraits<PixelType>::max();
m_ThreadMax = NumericTraits<PixelType>::NonpositiveMin();

if (this->GetImageSpatialMask())
{
this->m_SameGeometry = elastix::MaskHasSameImageDomain(*m_ImageSpatialMask, *(this->GetInput()));
}
else
{
this->m_SameGeometry = false;
}
}
m_ThreadMin = NumericTraits<PixelType>::max();
m_ThreadMax = NumericTraits<PixelType>::NonpositiveMin();
m_SameGeometry =
(m_ImageSpatialMask != nullptr) && elastix::MaskHasSameImageDomain(*m_ImageSpatialMask, *(this->GetInput()));
}

template <typename TInputImage>
void
ComputeImageExtremaFilter<TInputImage>::AfterStreamedGenerateData()
{
if (!this->m_UseMask)
{
Superclass::AfterStreamedGenerateData();
}
else
{
const SizeValueType count = m_Count;
const RealType sumOfSquares(m_SumOfSquares);
const PixelType minimum = m_ThreadMin;
const PixelType maximum = m_ThreadMax;
const RealType sum(m_ThreadSum);

const RealType mean = sum / static_cast<RealType>(count);
const RealType variance =
(sumOfSquares - (sum * sum / static_cast<RealType>(count))) / (static_cast<RealType>(count) - 1);
const RealType sigma = std::sqrt(variance);

// Set the outputs
this->SetMinimum(minimum);
this->SetMaximum(maximum);
this->SetMean(mean);
this->SetSigma(sigma);
this->SetVariance(variance);
this->SetSum(sum);
this->SetSumOfSquares(sumOfSquares);
}
}

template <typename TInputImage>
void
ComputeImageExtremaFilter<TInputImage>::ThreadedStreamedGenerateData(const RegionType & regionForThread)
{
if (!this->m_UseMask)
{
Superclass::ThreadedStreamedGenerateData(regionForThread);
}
else
{
if (this->GetImageSpatialMask())
{
this->ThreadedGenerateDataImageSpatialMask(regionForThread);
}
}
} // end ThreadedGenerateData()

template <typename TInputImage>
void
ComputeImageExtremaFilter<TInputImage>::ThreadedGenerateDataImageSpatialMask(const RegionType & regionForThread)
auto
ComputeImageExtremaFilter<TInputImage>::RetrieveMinMax(const TInputImage & inputImage,
const InputImageRegionType & regionForThread,
const ImageSpatialMaskType * const imageSpatialMask,
const bool sameGeometry) -> MinMaxResult
{
if (regionForThread.GetSize(0) == 0)
{
return;
}
RealType sum{};
RealType sumOfSquares{};
SizeValueType count{};
PixelType min = NumericTraits<PixelType>::max();
PixelType max = NumericTraits<PixelType>::NonpositiveMin();

const auto & inputImage = *(this->GetInput());
PixelType min = NumericTraits<PixelType>::max();
PixelType max = NumericTraits<PixelType>::NonpositiveMin();

if (this->m_SameGeometry)
if (imageSpatialMask)
{
const auto & maskImage = *(this->m_ImageSpatialMask->GetImage());
if (sameGeometry)
{
const auto & maskImage = *(imageSpatialMask->GetImage());

for (ImageRegionConstIterator<TInputImage> it(&inputImage, regionForThread); !it.IsAtEnd(); ++it)
for (ImageRegionConstIterator<TInputImage> it(&inputImage, regionForThread); !it.IsAtEnd(); ++it)
{
if (maskImage.GetPixel(it.GetIndex()) != PixelType{})
{
const PixelType value = it.Get();
min = std::min(min, value);
max = std::max(max, value);
}
}
}
else
{
if (maskImage.GetPixel(it.GetIndex()) != PixelType{})
for (ImageRegionConstIterator<TInputImage> it(&inputImage, regionForThread); !it.IsAtEnd(); ++it)
{
const PixelType value = it.Get();
const auto realValue = static_cast<RealType>(value);

min = std::min(min, value);
max = std::max(max, value);

sum += realValue;
sumOfSquares += (realValue * realValue);
++count;
typename ImageSpatialMaskType::PointType point;
inputImage.TransformIndexToPhysicalPoint(it.GetIndex(), point);
if (imageSpatialMask->IsInsideInWorldSpace(point))
{
const PixelType value = it.Get();
min = std::min(min, value);
max = std::max(max, value);
}
}
} // end for
}
}
else
{
for (ImageRegionConstIterator<TInputImage> it(&inputImage, regionForThread); !it.IsAtEnd(); ++it)
for (ImageScanlineConstIterator<TInputImage> it(&inputImage, regionForThread); !it.IsAtEnd(); it.NextLine())
{
PointType point;
inputImage.TransformIndexToPhysicalPoint(it.GetIndex(), point);
if (this->m_ImageSpatialMask->IsInsideInWorldSpace(point))
while (!it.IsAtEndOfLine())
{
const PixelType value = it.Get();
const auto realValue = static_cast<RealType>(value);

min = std::min(min, value);
max = std::max(max, value);

sum += realValue;
sumOfSquares += (realValue * realValue);
++count;
++it;
}
} // end for
} // end if
}
}
return { min, max };
}

const std::lock_guard<std::mutex> lockGuard(m_Mutex);
m_ThreadSum += sum;
m_SumOfSquares += sumOfSquares;
m_Count += count;
m_ThreadMin = std::min(min, m_ThreadMin);
m_ThreadMax = std::max(max, m_ThreadMax);

} // end ThreadedGenerateDataImageSpatialMask()
template <typename TInputImage>
void
ComputeImageExtremaFilter<TInputImage>::ThreadedStreamedGenerateData(const InputImageRegionType & regionForThread)
{
if (regionForThread.GetSize(0) > 0)
{
const MinMaxResult minMaxResult =
RetrieveMinMax(*(this->GetInput()), regionForThread, m_ImageSpatialMask, m_SameGeometry);

// Lock after calling RetrieveMinMax.
const std::lock_guard<std::mutex> lockGuard(m_Mutex);
m_ThreadMin = std::min(minMaxResult.Min, m_ThreadMin);
m_ThreadMax = std::max(minMaxResult.Max, m_ThreadMax);
}
}

} // end namespace itk
#endif
Loading
Loading