From 3c1b13d7a30336675c95b0fc60487a356f645844 Mon Sep 17 00:00:00 2001 From: Tom Birdsong Date: Fri, 21 Jan 2022 13:16:45 -0500 Subject: [PATCH 1/5] ENH: Expose kernel parameters in `DiscreteGaussianImageFilter` interface `itk::DiscreteGaussianImageFilter` exposes parameters for generating a Gaussian kernel to error and size specifications, but previously did not make information on the kernel that was actually generated available to other classes. These changes expose protected methods to allow subclasses to create kernels to the same specifications and also publicly exposes information on the kernel sigma when image spacing is being accounted for. --- .../include/itkDiscreteGaussianImageFilter.h | 23 ++- .../itkDiscreteGaussianImageFilter.hxx | 135 ++++++++---------- .../itkDiscreteGaussianImageFilterTest.cxx | 28 +++- 3 files changed, 110 insertions(+), 76 deletions(-) diff --git a/Modules/Filtering/Smoothing/include/itkDiscreteGaussianImageFilter.h b/Modules/Filtering/Smoothing/include/itkDiscreteGaussianImageFilter.h index 7d856781dd2..d574c9c7e75 100644 --- a/Modules/Filtering/Smoothing/include/itkDiscreteGaussianImageFilter.h +++ b/Modules/Filtering/Smoothing/include/itkDiscreteGaussianImageFilter.h @@ -18,6 +18,7 @@ #ifndef itkDiscreteGaussianImageFilter_h #define itkDiscreteGaussianImageFilter_h +#include "itkGaussianOperator.h" #include "itkImageToImageFilter.h" #include "itkImage.h" #include "itkZeroFluxNeumannBoundaryCondition.h" @@ -113,6 +114,10 @@ class ITK_TEMPLATE_EXPORT DiscreteGaussianImageFilter : public ImageToImageFilte using SigmaArrayType = ArrayType; using ScalarRealType = double; + /** Type of kernel to be used in blurring */ + using KernelType = GaussianOperator; + using RadiusType = typename KernelType::RadiusType; + /** The variance for the discrete Gaussian kernel. Sets the variance * independently for each dimension, but * see also SetVariance(const double v). The default is 0.0 in each @@ -130,8 +135,8 @@ class ITK_TEMPLATE_EXPORT DiscreteGaussianImageFilter : public ImageToImageFilte /** Set the kernel to be no wider than MaximumKernelWidth pixels, * even if MaximumError demands it. The default is 32 pixels. */ - itkGetConstMacro(MaximumKernelWidth, int); - itkSetMacro(MaximumKernelWidth, int); + itkGetConstMacro(MaximumKernelWidth, unsigned int); + itkSetMacro(MaximumKernelWidth, unsigned int); /** Set the number of dimensions to smooth. Defaults to the image * dimension. Can be set to less than ImageDimension, smoothing all @@ -334,6 +339,18 @@ class ITK_TEMPLATE_EXPORT DiscreteGaussianImageFilter : public ImageToImageFilte void GenerateData() override; + /** Build a directional kernel to match user specifications */ + void + GenerateKernel(const unsigned int dimension, KernelType & oper) const; + + /** Get the radius of the generated directional kernel */ + unsigned int + GetKernelRadius(const unsigned int dimension) const; + + /** Get the variance, optionally adjusted for pixel spacing */ + ArrayType + GetKernelVarianceArray() const; + private: /** The variance of the gaussian blurring kernel in each dimensional direction. */ @@ -346,7 +363,7 @@ class ITK_TEMPLATE_EXPORT DiscreteGaussianImageFilter : public ImageToImageFilte /** Maximum allowed kernel width for any dimension of the discrete Gaussian approximation */ - int m_MaximumKernelWidth; + unsigned int m_MaximumKernelWidth; /** Number of dimensions to process. Default is all dimensions */ unsigned int m_FilterDimensionality; diff --git a/Modules/Filtering/Smoothing/include/itkDiscreteGaussianImageFilter.hxx b/Modules/Filtering/Smoothing/include/itkDiscreteGaussianImageFilter.hxx index 47d6d575c3c..6be73303430 100644 --- a/Modules/Filtering/Smoothing/include/itkDiscreteGaussianImageFilter.hxx +++ b/Modules/Filtering/Smoothing/include/itkDiscreteGaussianImageFilter.hxx @@ -26,6 +26,59 @@ namespace itk { +template +void +DiscreteGaussianImageFilter::GenerateKernel(const unsigned int dimension, + KernelType & oper) const +{ + // Determine the size of the operator in this dimension. Note that the + // Gaussian is built as a 1D operator in each of the specified directions. + oper.SetDirection(dimension); + oper.SetMaximumError(m_MaximumError[dimension]); + oper.SetMaximumKernelWidth(m_MaximumKernelWidth); + + oper.SetVariance(this->GetKernelVarianceArray()[dimension]); + + oper.CreateDirectional(); +} + +template +unsigned int +DiscreteGaussianImageFilter::GetKernelRadius(const unsigned int dimension) const +{ + KernelType oper; + this->GenerateKernel(dimension, oper); + return oper.GetRadius(dimension); +} + +template +typename DiscreteGaussianImageFilter::ArrayType +DiscreteGaussianImageFilter::GetKernelVarianceArray() const +{ + if (m_UseImageSpacing) + { + if (this->GetInput() == nullptr) + { + itkExceptionMacro("UseImageSpacing is ON but no input image was provided"); + } + + ArrayType adjustedVariance; + // Adjusted variance = var / (spacing ^ 2) + for (unsigned int dim = 0; dim < ImageDimension; ++dim) + { + // convert the variance from physical units to pixels + double s = this->GetInput()->GetSpacing()[dim]; + s = s * s; + adjustedVariance[dim] = m_Variance[dim] / s; + } + return adjustedVariance; + } + else + { + return this->GetVariance(); + } +} + template void DiscreteGaussianImageFilter::GenerateInputRequestedRegion() @@ -42,39 +95,18 @@ DiscreteGaussianImageFilter::GenerateInputRequestedRe return; } - // Build an operator so that we can determine the kernel size - GaussianOperator oper; - - typename TInputImage::SizeType radius; - + // Determine the kernel size in each direction + RadiusType radius; for (unsigned int i = 0; i < TInputImage::ImageDimension; ++i) { - // Determine the size of the operator in this dimension. Note that the - // Gaussian is built as a 1D operator in each of the specified directions. - oper.SetDirection(i); - if (m_UseImageSpacing == true) + if (i < m_FilterDimensionality) { - if (this->GetInput()->GetSpacing()[i] == 0.0) - { - itkExceptionMacro(<< "Pixel spacing cannot be zero"); - } - else - { - // convert the variance from physical units to pixels - double s = this->GetInput()->GetSpacing()[i]; - s = s * s; - oper.SetVariance(m_Variance[i] / s); - } + radius[i] = GetKernelRadius(i); } else { - oper.SetVariance(m_Variance[i]); + radius[i] = 0; } - oper.SetMaximumError(m_MaximumError[i]); - oper.SetMaximumKernelWidth(m_MaximumKernelWidth); - oper.CreateDirectional(); - - radius[i] = oper.GetRadius(i); } // get a copy of the input requested region (should equal the output @@ -86,26 +118,9 @@ DiscreteGaussianImageFilter::GenerateInputRequestedRe inputRequestedRegion.PadByRadius(radius); // crop the input requested region at the input's largest possible region - if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion())) - { - inputPtr->SetRequestedRegion(inputRequestedRegion); - return; - } - else - { - // Couldn't crop the region (requested region is outside the largest - // possible region). Throw an exception. - - // store what we tried to request (prior to trying to crop) - inputPtr->SetRequestedRegion(inputRequestedRegion); - - // build an exception - InvalidRequestedRegionError e(__FILE__, __LINE__); - e.SetLocation(ITK_LOCATION); - e.SetDescription("Requested region is (at least partially) outside the largest possible region."); - e.SetDataObject(inputPtr); - throw e; - } + inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()); + + inputPtr->SetRequestedRegion(inputRequestedRegion); } template @@ -160,9 +175,8 @@ DiscreteGaussianImageFilter::GenerateData() using SingleFilterPointer = typename SingleFilterType::Pointer; // Create a series of operators - using OperatorType = GaussianOperator; - std::vector oper; + std::vector oper; oper.resize(filterDimensionality); // Create a process accumulator for tracking the progress of minipipeline @@ -177,30 +191,7 @@ DiscreteGaussianImageFilter::GenerateData() // the largest dimension will be split slice wise for streaming unsigned int reverse_i = filterDimensionality - i - 1; - // Set up the operator for this dimension - oper[reverse_i].SetDirection(i); - if (m_UseImageSpacing == true) - { - if (localInput->GetSpacing()[i] == 0.0) - { - itkExceptionMacro(<< "Pixel spacing cannot be zero"); - } - else - { - // convert the variance from physical units to pixels - double s = localInput->GetSpacing()[i]; - s = s * s; - oper[reverse_i].SetVariance(m_Variance[i] / s); - } - } - else - { - oper[reverse_i].SetVariance(m_Variance[i]); - } - - oper[reverse_i].SetMaximumKernelWidth(m_MaximumKernelWidth); - oper[reverse_i].SetMaximumError(m_MaximumError[i]); - oper[reverse_i].CreateDirectional(); + this->GenerateKernel(i, oper[reverse_i]); } // Create a chain of filters diff --git a/Modules/Filtering/Smoothing/test/itkDiscreteGaussianImageFilterTest.cxx b/Modules/Filtering/Smoothing/test/itkDiscreteGaussianImageFilterTest.cxx index 09f731c4d65..6caf8d8526f 100644 --- a/Modules/Filtering/Smoothing/test/itkDiscreteGaussianImageFilterTest.cxx +++ b/Modules/Filtering/Smoothing/test/itkDiscreteGaussianImageFilterTest.cxx @@ -71,6 +71,10 @@ itkDiscreteGaussianImageFilterTest(int argc, char * argv[]) filter->SetSigmaArray(sigma); ITK_TEST_SET_GET_VALUE(sigma, filter->GetSigmaArray()); + // Verify spacing-adjusted variance is not available without input image + filter->SetUseImageSpacing(true); + ITK_TRY_EXPECT_EXCEPTION(filter->GetImageSpacingVarianceArray()); + // Set the boundary condition used by the filter itk::ConstantBoundaryCondition constantBoundaryCondition; filter->SetInputBoundaryCondition(&constantBoundaryCondition); @@ -92,7 +96,7 @@ itkDiscreteGaussianImageFilterTest(int argc, char * argv[]) filter->SetMaximumError(maximumError); ITK_TEST_SET_GET_VALUE(maximumError, filter->GetMaximumError()); - int maximumKernelWidth = 32; + unsigned int maximumKernelWidth = 32; filter->SetMaximumKernelWidth(maximumKernelWidth); ITK_TEST_SET_GET_VALUE(maximumKernelWidth, filter->GetMaximumKernelWidth()); @@ -124,6 +128,28 @@ itkDiscreteGaussianImageFilterTest(int argc, char * argv[]) ITK_TRY_EXPECT_NO_EXCEPTION(test1.Execute()); + // Test variance adjustments for image spacing + ImageType::Pointer inputImage = ImageType::New(); + typename ImageType::SpacingType spacing; + spacing[0] = 1.5; + spacing[1] = 1.0; + spacing[2] = 0.5; + inputImage->SetSpacing(spacing); + filter->SetInput(inputImage); + + filter->SetUseImageSpacing(false); + for (itk::IndexValueType dim = 0; dim < Dimension; ++dim) + { + ITK_TEST_EXPECT_EQUAL(sigmaValue, filter->GetVariance()[dim]) + ITK_TEST_EXPECT_EQUAL(sigmaValue, filter->GetImageSpacingVarianceArray()[dim]); + } + + filter->SetUseImageSpacing(true); + for (itk::IndexValueType dim = 0; dim < Dimension; ++dim) + { + ITK_TEST_EXPECT_EQUAL(sigmaValue, filter->GetVariance()[dim]); + ITK_TEST_EXPECT_EQUAL((sigmaValue / (spacing[dim] * spacing[dim])), filter->GetImageSpacingVarianceArray()[dim]); + } std::cout << "Test finished" << std::endl; return EXIT_SUCCESS; From a18d2e488725ed252558121ab66849580ad5ee25 Mon Sep 17 00:00:00 2001 From: Tom Birdsong Date: Fri, 14 Jan 2022 17:48:10 -0500 Subject: [PATCH 2/5] ENH: Add `FFTDiscreteGaussianImageFilter` Implements `FFTDiscreteGaussianImageFilter` as a subclass of `DiscreteGaussianImageFilter` performing convolution of an ND single-channel input image with a Gaussian kernel of fixed width in the Fourier domain. Includes wrappings, test and baseline image. --- .../include/itkDiscreteGaussianImageFilter.h | 18 +- .../itkFFTDiscreteGaussianImageFilter.h | 134 ++++++++++++++ .../itkFFTDiscreteGaussianImageFilter.hxx | 173 ++++++++++++++++++ Modules/Filtering/Smoothing/itk-module.cmake | 4 + ...teGaussianImageFilterTestOutput.mha.sha512 | 1 + .../Filtering/Smoothing/test/CMakeLists.txt | 10 + .../itkFFTDiscreteGaussianImageFilterTest.cxx | 115 ++++++++++++ .../itkFFTDiscreteGaussianImageFilter.wrap | 3 + 8 files changed, 449 insertions(+), 9 deletions(-) create mode 100644 Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilter.h create mode 100644 Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilter.hxx create mode 100644 Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterTestOutput.mha.sha512 create mode 100644 Modules/Filtering/Smoothing/test/itkFFTDiscreteGaussianImageFilterTest.cxx create mode 100644 Modules/Filtering/Smoothing/wrapping/itkFFTDiscreteGaussianImageFilter.wrap diff --git a/Modules/Filtering/Smoothing/include/itkDiscreteGaussianImageFilter.h b/Modules/Filtering/Smoothing/include/itkDiscreteGaussianImageFilter.h index d574c9c7e75..d9da5834660 100644 --- a/Modules/Filtering/Smoothing/include/itkDiscreteGaussianImageFilter.h +++ b/Modules/Filtering/Smoothing/include/itkDiscreteGaussianImageFilter.h @@ -298,15 +298,6 @@ class ITK_TEMPLATE_EXPORT DiscreteGaussianImageFilter : public ImageToImageFilte itkLegacyMacro(unsigned int GetInternalNumberOfStreamDivisions() const); itkLegacyMacro(void SetInternalNumberOfStreamDivisions(unsigned int)); - /** DiscreteGaussianImageFilter needs a larger input requested region - * than the output requested region (larger by the size of the - * Gaussian kernel). As such, DiscreteGaussianImageFilter needs to - * provide an implementation for GenerateInputRequestedRegion() in - * order to inform the pipeline execution model. - * \sa ImageToImageFilter::GenerateInputRequestedRegion() */ - void - GenerateInputRequestedRegion() override; - #ifdef ITK_USE_CONCEPT_CHECKING // Begin concept checking @@ -331,6 +322,15 @@ class ITK_TEMPLATE_EXPORT DiscreteGaussianImageFilter : public ImageToImageFilte void PrintSelf(std::ostream & os, Indent indent) const override; + /** DiscreteGaussianImageFilter needs a larger input requested region + * than the output requested region (larger by the size of the + * Gaussian kernel). As such, DiscreteGaussianImageFilter needs to + * provide an implementation for GenerateInputRequestedRegion() in + * order to inform the pipeline execution model. + * \sa ImageToImageFilter::GenerateInputRequestedRegion() */ + void + GenerateInputRequestedRegion() override; + /** Standard pipeline method. While this class does not implement a * ThreadedGenerateData(), its GenerateData() delegates all * calculations to an NeighborhoodOperatorImageFilter. Since the diff --git a/Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilter.h b/Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilter.h new file mode 100644 index 00000000000..cfc371fa896 --- /dev/null +++ b/Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilter.h @@ -0,0 +1,134 @@ +/*========================================================================= + * + * 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 + * + * http://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 itkFFTDiscreteGaussianImageFilter_h +#define itkFFTDiscreteGaussianImageFilter_h + +#include "itkDiscreteGaussianImageFilter.h" +#include "itkImage.h" +#include "itkMacro.h" + +namespace itk +{ +/** + * \class FFTDiscreteGaussianImageFilter + * \brief Blurs an image by convolution with a discrete gaussian kernel + * in the frequency domain. + * + * This filter performs Gaussian blurring by convolution of an image + * and a discrete Gaussian operator (kernel) in the frequency domain + * by way of Fast Fourier Transform forward and inverse operations. + * + * \sa DiscreteGaussianImageFilter + * \sa GaussianImageSource + * \sa FFTConvolutionImageFilter + * \sa RecursiveGaussianImageFilter + * \sa ZeroFluxNeumannBoundaryCondition + * + * \ingroup ImageEnhancement + * \ingroup ImageFeatureExtraction + * \ingroup ITKSmoothing + */ + +template +class ITK_TEMPLATE_EXPORT FFTDiscreteGaussianImageFilter : public DiscreteGaussianImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(FFTDiscreteGaussianImageFilter); + + /** Standard class type aliases. */ + using Self = FFTDiscreteGaussianImageFilter; + using Superclass = DiscreteGaussianImageFilter; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(FFTDiscreteGaussianImageFilter, DiscreteGaussianImageFilter); + + /** Image type information. */ + using InputImageType = typename Superclass::InputImageType; + using OutputImageType = typename Superclass::OutputImageType; + + /** Extract some information from the image types. Dimensionality + * of the two images is assumed to be the same. */ + using OutputPixelType = typename Superclass::OutputPixelType; + using OutputInternalPixelType = typename Superclass::OutputInternalPixelType; + using InputPixelType = typename Superclass::InputPixelType; + using InputInternalPixelType = typename Superclass::InputInternalPixelType; + + /** Pixel value type for Vector pixel types **/ + using InputPixelValueType = typename Superclass::InputPixelValueType; + using OutputPixelValueType = typename Superclass::OutputPixelValueType; + + /** Extract some information from the image types. Dimensionality + * of the two images is assumed to be the same. */ + static constexpr unsigned int ImageDimension = Superclass::ImageDimension; + + /** Type of the pixel to use for intermediate results */ + using RealOutputPixelType = typename Superclass::RealOutputPixelType; + using RealOutputImageType = typename Superclass::RealOutputImageType; + using RealOutputPixelValueType = typename Superclass::RealOutputPixelValueType; + using RealPixelType = RealOutputPixelType; + using RealImageType = RealOutputImageType; + + /** Typedef to describe the boundary condition. */ + using BoundaryConditionType = typename Superclass::BoundaryConditionType; + using InputBoundaryConditionPointerType = typename Superclass::InputBoundaryConditionPointerType; + using InputDefaultBoundaryConditionType = typename Superclass::InputDefaultBoundaryConditionType; + using RealBoundaryConditionPointerType = typename Superclass::RealBoundaryConditionPointerType; + using RealDefaultBoundaryConditionType = typename Superclass::RealDefaultBoundaryConditionType; + + /** Typedef of double containers */ + using ArrayType = typename Superclass::ArrayType; + using SigmaArrayType = typename Superclass::SigmaArrayType; + using ScalarRealType = typename Superclass::ScalarRealType; + + /** Typedef to describe kernel parameters */ + using typename Superclass::KernelType; + using typename Superclass::RadiusType; + + /** Overridden accessors for unused parameters */ + + void + SetInputBoundaryCondition(const InputBoundaryConditionPointerType) override; + +protected: + FFTDiscreteGaussianImageFilter() = default; + ~FFTDiscreteGaussianImageFilter() override = default; + + // Pad input region to kernel image size + void + GenerateInputRequestedRegion() override; + + /** Standard pipeline method. While this class does not implement a + * ThreadedGenerateData(), its GenerateData() delegates all + * calculations to an NeighborhoodOperatorImageFilter. Since the + * NeighborhoodOperatorImageFilter is multithreaded, this filter is + * multithreaded by default. */ + void + GenerateData() override; +}; +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkFFTDiscreteGaussianImageFilter.hxx" +#endif + +#endif diff --git a/Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilter.hxx b/Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilter.hxx new file mode 100644 index 00000000000..5a352bd3817 --- /dev/null +++ b/Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilter.hxx @@ -0,0 +1,173 @@ +/*========================================================================= + * + * 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 + * + * http://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 itkFFTDiscreteGaussianImageFilter_hxx +#define itkFFTDiscreteGaussianImageFilter_hxx + +#include "itkGaussianOperator.h" +#include "itkGaussianImageSource.h" +#include "itkFFTConvolutionImageFilter.h" + +#include "itkProgressAccumulator.h" +#include "itkImageAlgorithm.h" +#include "itkMacro.h" + +namespace itk +{ + +template +void +FFTDiscreteGaussianImageFilter::GenerateInputRequestedRegion() +{ + // call the superclass' implementation of this method. this should + // copy the output requested region to the input requested region + ImageToImageFilter::GenerateInputRequestedRegion(); + + // Get pointer to input + typename Superclass::InputImagePointer inputPtr = const_cast(this->GetInput()); + if (inputPtr.IsNull()) + { + return; + } + + // get a copy of the input requested region (should equal the output + // requested region) + typename TInputImage::RegionType inputRequestedRegion; + inputRequestedRegion = inputPtr->GetRequestedRegion(); + + // pad the input requested region by the operator radius + RadiusType radius; + radius.Fill(0); + for (size_t dim = 0; dim < ImageDimension; ++dim) + { + radius[dim] = this->GetKernelRadius(dim); + } + inputRequestedRegion.PadByRadius(radius); + + // crop the input requested region at the input's largest possible region + inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()); + + inputPtr->SetRequestedRegion(inputRequestedRegion); +} + +template +void +FFTDiscreteGaussianImageFilter::SetInputBoundaryCondition( + const InputBoundaryConditionPointerType) +{ + itkWarningMacro("FFTDiscreteGaussianImageFilter ignores InputBoundaryCondition, use RealBoundaryCondition instead"); +} + +template +void +FFTDiscreteGaussianImageFilter::GenerateData() +{ + TOutputImage * output = this->GetOutput(); + + output->SetBufferedRegion(output->GetRequestedRegion()); + output->Allocate(); + + // Create an internal image to protect the input image's metdata + // (e.g. RequestedRegion). The StreamingImageFilter changes the + // requested region as part of its normal processing. + auto localInput = TInputImage::New(); + localInput->Graft(this->GetInput()); + + // Determine the dimensionality to filter + unsigned int filterDimensionality = this->GetFilterDimensionality(); + if (filterDimensionality > ImageDimension) + { + filterDimensionality = ImageDimension; + } + if (filterDimensionality == 0) + { + // no smoothing, copy input to output + ImageAlgorithm::Copy(localInput.GetPointer(), + output, + this->GetOutput()->GetRequestedRegion(), + this->GetOutput()->GetRequestedRegion()); + return; + } + + // Create a process accumulator for tracking the progress of minipipeline + auto progress = ProgressAccumulator::New(); + progress->SetMiniPipelineFilter(this); + + // Create kernel image for blurring in requested dimensions + using GaussianImageSourceType = GaussianImageSource; + using KernelSizeType = typename GaussianImageSourceType::SizeType; + using KernelMeanType = typename GaussianImageSourceType::ArrayType; + + typename GaussianImageSourceType::Pointer kernelSource = GaussianImageSourceType::New(); + + kernelSource->SetScale(1.0); + kernelSource->SetNormalized(true); + + kernelSource->SetSpacing(localInput->GetSpacing()); + kernelSource->SetOrigin(localInput->GetOrigin()); + kernelSource->SetDirection(localInput->GetDirection()); + + KernelSizeType kernelSize; + kernelSize.Fill(1); + for (size_t dim = 0; dim < filterDimensionality; ++dim) + { + kernelSize[dim] = static_cast(this->GetKernelRadius(dim)) * 2 + 1; + } + kernelSource->SetSize(kernelSize); + + KernelMeanType mean; + auto inputSpacing = localInput->GetSpacing(); + auto inputOrigin = localInput->GetOrigin(); + for (size_t dim = 0; dim < ImageDimension; ++dim) + { + double radius = (kernelSize[dim] - 1) / 2; + mean[dim] = inputSpacing[dim] * radius + inputOrigin[dim]; // center pixel pos + } + kernelSource->SetMean(mean); + kernelSource->SetSigma(this->GetSigmaArray()); + + // Estimate kernel generation to be roughly 5% of effort + progress->RegisterInternalFilter(kernelSource, 0.05f); + + // Perform image convolution by FFT + using ConvolutionImageFilterType = FFTConvolutionImageFilter; + + typename ConvolutionImageFilterType::Pointer convolutionFilter = ConvolutionImageFilterType::New(); + convolutionFilter->SetInput(this->GetInput()); + convolutionFilter->SetKernelImage(kernelSource->GetOutput()); + convolutionFilter->SetBoundaryCondition(this->GetRealBoundaryCondition()); + convolutionFilter->SetNormalize(false); // Kernel is already normalized + + // Estimate kernel generation to be roughly 95% of effort + progress->RegisterInternalFilter(convolutionFilter, 0.95f); + + // Graft this filters output onto the mini-pipeline so the mini-pipeline + // has the correct region ivars and will write to this filters bulk data + // output. + convolutionFilter->GraftOutput(output); + + // Update the last filter in the mini-pipeline + convolutionFilter->Update(); + + // Graft the last output of the mini-pipeline onto this filters output so + // the final output has the correct region ivars and a handle to the final + // bulk data + this->GraftOutput(output); +} +} // namespace itk + +#endif diff --git a/Modules/Filtering/Smoothing/itk-module.cmake b/Modules/Filtering/Smoothing/itk-module.cmake index 567ecc3b526..2909d21e5c5 100644 --- a/Modules/Filtering/Smoothing/itk-module.cmake +++ b/Modules/Filtering/Smoothing/itk-module.cmake @@ -5,8 +5,12 @@ interesting to look at the ITKAnisotropicSmoothing group of filters.") itk_module(ITKSmoothing ENABLE_SHARED COMPILE_DEPENDS + ITKConvolution + ITKFFT ITKImageFunction + ITKImageSources TEST_DEPENDS + ITKConvolution ITKTestKernel DESCRIPTION "${DOCUMENTATION}" diff --git a/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterTestOutput.mha.sha512 b/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterTestOutput.mha.sha512 new file mode 100644 index 00000000000..8f1ea54c767 --- /dev/null +++ b/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterTestOutput.mha.sha512 @@ -0,0 +1 @@ +0eeb11ed3467d7af854c70506ea3e476fb73f57aaa1779ba3eb8abbfc30035bf4b0c11a7e47d12c65f01f1682dcac11a01e8ac20f1a0398ddb5f47f668a24a75 diff --git a/Modules/Filtering/Smoothing/test/CMakeLists.txt b/Modules/Filtering/Smoothing/test/CMakeLists.txt index db9c444cff2..e44cfcaf0ff 100644 --- a/Modules/Filtering/Smoothing/test/CMakeLists.txt +++ b/Modules/Filtering/Smoothing/test/CMakeLists.txt @@ -3,6 +3,7 @@ set(ITKSmoothingTests itkBoxMeanImageFilterTest.cxx itkBoxSigmaImageFilterTest.cxx itkDiscreteGaussianImageFilterTest2.cxx +itkFFTDiscreteGaussianImageFilterTest.cxx itkSmoothingRecursiveGaussianImageFilterTest.cxx itkSmoothingRecursiveGaussianImageFilterOnVectorImageTest.cxx itkSmoothingRecursiveGaussianImageFilterOnImageOfVectorTest.cxx @@ -59,6 +60,15 @@ itk_add_test(NAME itkDiscreteGaussianImageFilterTest1a COMMAND ITKSmoothingTestDriver itkDiscreteGaussianImageFilterTest 1) itk_add_test(NAME itkDiscreteGaussianImageFilterTest1b COMMAND ITKSmoothingTestDriver itkDiscreteGaussianImageFilterTest 0) +itk_add_test(NAME itkFFTDiscreteGaussianImageFilterTest + COMMAND ITKSmoothingTestDriver + --compare DATA{Baseline/FFTDiscreteGaussianImageFilterTestOutput.mha} + ${ITK_TEST_OUTPUT_DIR}/FFTDiscreteGaussianImageFilterTestOutput.mha + itkFFTDiscreteGaussianImageFilterTest + DATA{${ITK_DATA_ROOT}/Input/2th_cthead1.png} + 3.5 + 35 + ${ITK_TEST_OUTPUT_DIR}/FFTDiscreteGaussianImageFilterTestOutput.mha) itk_add_test(NAME itkMedianImageFilterTest COMMAND ITKSmoothingTestDriver itkMedianImageFilterTest) itk_add_test(NAME itkRecursiveGaussianImageFiltersOnTensorsTest diff --git a/Modules/Filtering/Smoothing/test/itkFFTDiscreteGaussianImageFilterTest.cxx b/Modules/Filtering/Smoothing/test/itkFFTDiscreteGaussianImageFilterTest.cxx new file mode 100644 index 00000000000..2ded25decb7 --- /dev/null +++ b/Modules/Filtering/Smoothing/test/itkFFTDiscreteGaussianImageFilterTest.cxx @@ -0,0 +1,115 @@ +/*========================================================================= + * + * 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 + * + * http://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 + +#include "itkConstantBoundaryCondition.h" +#include "itkFFTDiscreteGaussianImageFilter.h" +#include "itkImageFileReader.h" +#include "itkImageFileWriter.h" +#include "itkTestingMacros.h" +#include "itkZeroFluxNeumannBoundaryCondition.h" + +int +itkFFTDiscreteGaussianImageFilterTest(int argc, char * argv[]) +{ + if (argc < 5) + { + std::cerr << "Missing parameters." << std::endl; + std::cerr << "Usage:" << std::endl; + std::cerr << itkNameOfTestExecutableMacro(argv) + << " inputFilename sigma kernelWidth outputFilename filterDimensionality" << std::endl; + return EXIT_FAILURE; + } + + using ScalarPixelType = float; + const size_t ImageDimension = 2; + using ImageType = itk::Image; + + double sigma = std::stod(argv[2]); + unsigned int kernelWidth = std::stoi(argv[3]); + unsigned int filterDimensionality = (argc < 6 ? ImageDimension : std::stoi(argv[5])); + + typename ImageType::Pointer inputImage = itk::ReadImage(argv[1]); + + using FilterType = itk::FFTDiscreteGaussianImageFilter; + auto filter = FilterType::New(); + ITK_EXERCISE_BASIC_OBJECT_METHODS(filter, FFTDiscreteGaussianImageFilter, DiscreteGaussianImageFilter); + + // Test setting inputs + + filter->SetInput(inputImage); + + filter->SetVariance(sigma); + for (auto & param : filter->GetVariance()) + { + ITK_TEST_EXPECT_EQUAL(sigma, param); + } + + filter->SetUseImageSpacingOff(); + for (auto & param : filter->GetImageSpacingVarianceArray()) + { + ITK_TEST_EXPECT_EQUAL(sigma, param); + } + + using SpacingType = typename ImageType::SpacingType; + SpacingType spacing; + spacing[0] = 0.5; + spacing[1] = 0.25; + inputImage->SetSpacing(spacing); + filter->SetUseImageSpacingOn(); + for (size_t dim = 0; dim < ImageDimension; ++dim) + { + ITK_TEST_EXPECT_EQUAL(sigma, filter->GetVariance()[dim]) + ITK_TEST_EXPECT_EQUAL(sigma / (spacing[dim] * spacing[dim]), filter->GetImageSpacingVarianceArray()[dim]); + } + + spacing[0] = 1.0; + spacing[1] = 1.0; + inputImage->SetSpacing(spacing); + + const float kernelError = 0.05; + filter->SetMaximumError(kernelError); + for (size_t dim = 0; dim < ImageDimension; ++dim) + { + ITK_TEST_SET_GET_VALUE(kernelError, filter->GetMaximumError()[dim]); + } + + filter->SetMaximumKernelWidth(kernelWidth); + ITK_TEST_SET_GET_VALUE(kernelWidth, filter->GetMaximumKernelWidth()); + + filter->SetFilterDimensionality(filterDimensionality); + ITK_TEST_SET_GET_VALUE(filterDimensionality, filter->GetFilterDimensionality()); + + itk::ZeroFluxNeumannBoundaryCondition zfnBoundaryCondition; + filter->SetRealBoundaryCondition(&zfnBoundaryCondition); + ITK_TEST_SET_GET_VALUE(&zfnBoundaryCondition, filter->GetRealBoundaryCondition()); + + // Test that attempting to set unused interface parameters throws warnings + // Ignore values from "Get" methods as these are unused + + itk::ConstantBoundaryCondition constantBoundaryCondition; + filter->SetInputBoundaryCondition(&constantBoundaryCondition); + + // Run convolution + + ITK_TRY_EXPECT_NO_EXCEPTION(filter->Update()); + + itk::WriteImage(filter->GetOutput(), argv[4], true); + + return EXIT_SUCCESS; +} diff --git a/Modules/Filtering/Smoothing/wrapping/itkFFTDiscreteGaussianImageFilter.wrap b/Modules/Filtering/Smoothing/wrapping/itkFFTDiscreteGaussianImageFilter.wrap new file mode 100644 index 00000000000..6c638893481 --- /dev/null +++ b/Modules/Filtering/Smoothing/wrapping/itkFFTDiscreteGaussianImageFilter.wrap @@ -0,0 +1,3 @@ +itk_wrap_class("itk::FFTDiscreteGaussianImageFilter" POINTER) + itk_wrap_image_filter("${WRAP_ITK_SCALAR}" 2) +itk_end_wrap_class() From 50e5e9c12462a0ac14545b0053d3f989e956b00d Mon Sep 17 00:00:00 2001 From: Tom Birdsong Date: Fri, 21 Jan 2022 16:07:00 -0500 Subject: [PATCH 3/5] ENH: Add `FFTDiscreteGaussianImageFilter` factory override Adds a factory class for manual override of `DiscreteGaussianImageFilter` with `FFTDiscreteGaussianImageFilter` through object factory instantiation. This provides an optional method for speeding up image blurring over very large images or using a very large kernel via GPU-accelerated FFTs inside other classes relying on `DiscreteGaussianImageFilter` for internal blurring (see `itkMultiResolutionPyramidImageFilter.h`). --- ...itkFFTDiscreteGaussianImageFilterFactory.h | 113 ++++++++++++++++++ .../Filtering/Smoothing/test/CMakeLists.txt | 3 + ...DiscreteGaussianImageFilterFactoryTest.cxx | 79 ++++++++++++ ...FFTDiscreteGaussianImageFilterFactory.wrap | 1 + 4 files changed, 196 insertions(+) create mode 100644 Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilterFactory.h create mode 100644 Modules/Filtering/Smoothing/test/itkFFTDiscreteGaussianImageFilterFactoryTest.cxx create mode 100644 Modules/Filtering/Smoothing/wrapping/itkFFTDiscreteGaussianImageFilterFactory.wrap diff --git a/Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilterFactory.h b/Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilterFactory.h new file mode 100644 index 00000000000..fc78ea68a9a --- /dev/null +++ b/Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilterFactory.h @@ -0,0 +1,113 @@ +/*========================================================================= + * + * 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 + * + * http://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 itkFFTDiscreteGaussianImageFilterFactory_h +#define itkFFTDiscreteGaussianImageFilterFactory_h +#include "ITKSmoothingExport.h" + +#include "itkFFTDiscreteGaussianImageFilter.h" +#include "itkImage.h" +#include "itkObjectFactoryBase.h" +#include "itkVersion.h" + +namespace itk +{ +/** \class FFTDiscreteGaussianImageFilterFactory + * + * \brief Object Factory implementation for overriding + * DiscreteGaussianImageFilter with FFTDiscreteGaussianImageFilter + * + * \sa ObjectFactoryBase + * \sa FFTDiscreteGaussianImageFilter + * \sa DiscreteGaussianImageFilter + * + * \ingroup ITKSmoothing + * \ingroup FourierTransform + * \ingroup ITKFFT + */ +class FFTDiscreteGaussianImageFilterFactory : public itk::ObjectFactoryBase +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(FFTDiscreteGaussianImageFilterFactory); + + using Self = FFTDiscreteGaussianImageFilterFactory; + using Superclass = ObjectFactoryBase; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Class methods used to interface with the registered factories. */ + const char * + GetITKSourceVersion() const override + { + return ITK_SOURCE_VERSION; + } + const char * + GetDescription() const override + { + return "An FFTDiscreteGaussianImageFilterFactory factory"; + } + + /** Method for class instantiation. */ + itkFactorylessNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(FFTDiscreteGaussianImageFilterFactory, itk::ObjectFactoryBase); + + /** Register one factory of this type */ + static void + RegisterOneFactory() + { + FFTDiscreteGaussianImageFilterFactory::Pointer factory = FFTDiscreteGaussianImageFilterFactory::New(); + + ObjectFactoryBase::RegisterFactoryInternal(factory); + } + +protected: + /** Override base DiscreteGaussianImageFilter constructor at runtime to return + * an upcast FFTDiscreteGaussianImageFilter instance through the object factory + */ + template + void + OverrideSuperclassType(const std::integer_sequence &) + { + using InputImageType = Image; + using OutputImageType = Image; + this->RegisterOverride( + typeid(typename FFTDiscreteGaussianImageFilter::Superclass).name(), + typeid(FFTDiscreteGaussianImageFilter).name(), + "FFTDiscreteGaussianImageFilter Override", + true, + CreateObjectFunction>::New()); + OverrideSuperclassType(std::integer_sequence{}); + } + template + void + OverrideSuperclassType(const std::integer_sequence &) + {} + + FFTDiscreteGaussianImageFilterFactory() + { + OverrideSuperclassType(std::integer_sequence{}); + + OverrideSuperclassType(std::integer_sequence{}); + } +}; + +} // namespace itk + +#endif // itkFFTDiscreteGaussianImageFilterFactory_h diff --git a/Modules/Filtering/Smoothing/test/CMakeLists.txt b/Modules/Filtering/Smoothing/test/CMakeLists.txt index e44cfcaf0ff..3b4f995dc1b 100644 --- a/Modules/Filtering/Smoothing/test/CMakeLists.txt +++ b/Modules/Filtering/Smoothing/test/CMakeLists.txt @@ -4,6 +4,7 @@ itkBoxMeanImageFilterTest.cxx itkBoxSigmaImageFilterTest.cxx itkDiscreteGaussianImageFilterTest2.cxx itkFFTDiscreteGaussianImageFilterTest.cxx +itkFFTDiscreteGaussianImageFilterFactoryTest.cxx itkSmoothingRecursiveGaussianImageFilterTest.cxx itkSmoothingRecursiveGaussianImageFilterOnVectorImageTest.cxx itkSmoothingRecursiveGaussianImageFilterOnImageOfVectorTest.cxx @@ -69,6 +70,8 @@ itk_add_test(NAME itkFFTDiscreteGaussianImageFilterTest 3.5 35 ${ITK_TEST_OUTPUT_DIR}/FFTDiscreteGaussianImageFilterTestOutput.mha) +itk_add_test(NAME itkFFTDiscreteGaussianImageFilterFactoryTest + COMMAND ITKSmoothingTestDriver itkFFTDiscreteGaussianImageFilterFactoryTest) itk_add_test(NAME itkMedianImageFilterTest COMMAND ITKSmoothingTestDriver itkMedianImageFilterTest) itk_add_test(NAME itkRecursiveGaussianImageFiltersOnTensorsTest diff --git a/Modules/Filtering/Smoothing/test/itkFFTDiscreteGaussianImageFilterFactoryTest.cxx b/Modules/Filtering/Smoothing/test/itkFFTDiscreteGaussianImageFilterFactoryTest.cxx new file mode 100644 index 00000000000..106e219ccb6 --- /dev/null +++ b/Modules/Filtering/Smoothing/test/itkFFTDiscreteGaussianImageFilterFactoryTest.cxx @@ -0,0 +1,79 @@ +/*========================================================================= + * + * 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 + * + * http://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 + +#include "itkFFTDiscreteGaussianImageFilterFactory.h" +#include "itkFFTDiscreteGaussianImageFilter.h" +#include "itkDiscreteGaussianImageFilter.h" +#include "itkImageFileReader.h" +#include "itkImageFileWriter.h" +#include "itkTestingMacros.h" + +int +itkFFTDiscreteGaussianImageFilterFactoryTest(int, char *[]) +{ + const unsigned int ImageDimension = 3; + using PixelType = float; + using ImageType = itk::Image; + using BaseFilterType = itk::DiscreteGaussianImageFilter; + using OverrideFilterType = typename itk::FFTDiscreteGaussianImageFilter; + + BaseFilterType::Pointer baseFilter; + + // Verify expected instance is returned before factory registration + ITK_TRY_EXPECT_NO_EXCEPTION(baseFilter = BaseFilterType::New()); + if (dynamic_cast(baseFilter.GetPointer()) != nullptr) + { + std::cout << "Object factory instantiation succeeded before a factory was registered" << std::endl; + return EXIT_FAILURE; + } + + // Register factory and verify instantiation override + using FactoryType = itk::FFTDiscreteGaussianImageFilterFactory; + FactoryType::Pointer overrideFactory = FactoryType::New(); + itk::ObjectFactoryBase::RegisterFactory(overrideFactory); + + ITK_TRY_EXPECT_NO_EXCEPTION(baseFilter = BaseFilterType::New()); + if (dynamic_cast(baseFilter.GetPointer()) == nullptr) + { + std::cout << "Object factory instantiation failed" << std::endl; + return EXIT_FAILURE; + } + + // Unregister factory and verify behavior returns to default + itk::ObjectFactoryBase::UnRegisterFactory(overrideFactory); + + ITK_TRY_EXPECT_NO_EXCEPTION(baseFilter = BaseFilterType::New()); + if (dynamic_cast(baseFilter.GetPointer()) != nullptr) + { + std::cout << "Object factory instantiation succeeded after factory was removed" << std::endl; + return EXIT_FAILURE; + } + + // Register factory through static method + itk::FFTDiscreteGaussianImageFilterFactory::RegisterOneFactory(); + + ITK_TRY_EXPECT_NO_EXCEPTION(baseFilter = BaseFilterType::New()); + if (dynamic_cast(baseFilter.GetPointer()) == nullptr) + { + std::cout << "Object factory instantiation failed" << std::endl; + return EXIT_FAILURE; + } + + return EXIT_SUCCESS; +} diff --git a/Modules/Filtering/Smoothing/wrapping/itkFFTDiscreteGaussianImageFilterFactory.wrap b/Modules/Filtering/Smoothing/wrapping/itkFFTDiscreteGaussianImageFilterFactory.wrap new file mode 100644 index 00000000000..f3531ff04b2 --- /dev/null +++ b/Modules/Filtering/Smoothing/wrapping/itkFFTDiscreteGaussianImageFilterFactory.wrap @@ -0,0 +1 @@ +itk_wrap_simple_class("itk::FFTDiscreteGaussianImageFilterFactory" POINTER) From 93ad716eebdedf1d73b6b54968ab9bbaf7c219c1 Mon Sep 17 00:00:00 2001 From: Tom Birdsong Date: Sat, 22 Jan 2022 15:09:53 -0500 Subject: [PATCH 4/5] ENH: Expand FFT Gaussian kernel options + improve test coverage --- .../itkFFTDiscreteGaussianImageFilter.h | 66 +++++++-- .../itkFFTDiscreteGaussianImageFilter.hxx | 134 +++++++++++++----- .../Filtering/Smoothing/src/CMakeLists.txt | 1 + .../src/itkFFTDiscreteGaussianImageFilter.cxx | 38 +++++ ...ImageFilterComparisonTestOutput.mha.sha512 | 1 + ...mageFilterTestImageSourceOutput.mha.sha512 | 1 + ...rTestKernelDimensionalityOutput.mha.sha512 | 1 + ...teGaussianImageFilterTestOutput.mha.sha512 | 1 - .../Filtering/Smoothing/test/CMakeLists.txt | 47 ++++-- .../itkDiscreteGaussianImageFilterTest.cxx | 29 +--- .../itkDiscreteGaussianImageFilterTest2.cxx | 29 ++-- .../itkFFTDiscreteGaussianImageFilterTest.cxx | 53 ++++--- .../itkFFTDiscreteGaussianImageFilter.wrap | 5 + 13 files changed, 284 insertions(+), 122 deletions(-) create mode 100644 Modules/Filtering/Smoothing/src/itkFFTDiscreteGaussianImageFilter.cxx create mode 100644 Modules/Filtering/Smoothing/test/Baseline/DiscreteGaussianImageFilterComparisonTestOutput.mha.sha512 create mode 100644 Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterTestImageSourceOutput.mha.sha512 create mode 100644 Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterTestKernelDimensionalityOutput.mha.sha512 delete mode 100644 Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterTestOutput.mha.sha512 diff --git a/Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilter.h b/Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilter.h index cfc371fa896..3332b9adc10 100644 --- a/Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilter.h +++ b/Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilter.h @@ -21,9 +21,42 @@ #include "itkDiscreteGaussianImageFilter.h" #include "itkImage.h" #include "itkMacro.h" +#include "ITKSmoothingExport.h" namespace itk { +/**\class FFTDiscreteGaussianImageFilterEnums + * \brief Contains all enum classes used by FFTDiscreteGaussianImageFilter class. + * \ingroup ITKSmoothing + */ +class FFTDiscreteGaussianImageFilterEnums +{ +public: + /** + * \class KernelSource + * \ingroup ITKSmoothing + * ITK defines multiple possible sources for generating a + * Gaussian kernel for smoothing. + * + * Generating an ND kernel by multiplying + * N 1D directional operators produces nearly identical + * results to the DiscreteGaussianImageFilter class. + * + * Generating an ND kernel image directly with + * GaussianImageSource relies on a slightly different + * spatial function and will produce different smoothing, + * but kernel generation will be faster. + */ + enum class KernelSource : uint8_t + { + COMBINED_OPERATORS = 0, + IMAGE_SOURCE, + }; +}; +// Define how to print enumeration +extern ITKSmoothing_EXPORT std::ostream & + operator<<(std::ostream & out, const FFTDiscreteGaussianImageFilterEnums::KernelSource value); + /** * \class FFTDiscreteGaussianImageFilter * \brief Blurs an image by convolution with a discrete gaussian kernel @@ -82,23 +115,24 @@ class ITK_TEMPLATE_EXPORT FFTDiscreteGaussianImageFilter : public DiscreteGaussi static constexpr unsigned int ImageDimension = Superclass::ImageDimension; /** Type of the pixel to use for intermediate results */ - using RealOutputPixelType = typename Superclass::RealOutputPixelType; - using RealOutputImageType = typename Superclass::RealOutputImageType; - using RealOutputPixelValueType = typename Superclass::RealOutputPixelValueType; + using typename Superclass::RealOutputPixelType; + using typename Superclass::RealOutputImageType; + using typename Superclass::RealOutputPixelValueType; using RealPixelType = RealOutputPixelType; using RealImageType = RealOutputImageType; + using RealImagePointerType = typename RealImageType::Pointer; /** Typedef to describe the boundary condition. */ - using BoundaryConditionType = typename Superclass::BoundaryConditionType; - using InputBoundaryConditionPointerType = typename Superclass::InputBoundaryConditionPointerType; - using InputDefaultBoundaryConditionType = typename Superclass::InputDefaultBoundaryConditionType; - using RealBoundaryConditionPointerType = typename Superclass::RealBoundaryConditionPointerType; - using RealDefaultBoundaryConditionType = typename Superclass::RealDefaultBoundaryConditionType; + using typename Superclass::BoundaryConditionType; + using typename Superclass::InputBoundaryConditionPointerType; + using typename Superclass::InputDefaultBoundaryConditionType; + using typename Superclass::RealBoundaryConditionPointerType; + using typename Superclass::RealDefaultBoundaryConditionType; /** Typedef of double containers */ - using ArrayType = typename Superclass::ArrayType; - using SigmaArrayType = typename Superclass::SigmaArrayType; - using ScalarRealType = typename Superclass::ScalarRealType; + using typename Superclass::ArrayType; + using typename Superclass::SigmaArrayType; + using typename Superclass::ScalarRealType; /** Typedef to describe kernel parameters */ using typename Superclass::KernelType; @@ -109,6 +143,13 @@ class ITK_TEMPLATE_EXPORT FFTDiscreteGaussianImageFilter : public DiscreteGaussi void SetInputBoundaryCondition(const InputBoundaryConditionPointerType) override; + auto + GenerateKernelImage() -> RealImagePointerType; + + + itkSetMacro(KernelSource, FFTDiscreteGaussianImageFilterEnums::KernelSource); + itkGetConstMacro(KernelSource, FFTDiscreteGaussianImageFilterEnums::KernelSource); + protected: FFTDiscreteGaussianImageFilter() = default; ~FFTDiscreteGaussianImageFilter() override = default; @@ -124,6 +165,9 @@ class ITK_TEMPLATE_EXPORT FFTDiscreteGaussianImageFilter : public DiscreteGaussi * multithreaded by default. */ void GenerateData() override; + +private: + FFTDiscreteGaussianImageFilterEnums::KernelSource m_KernelSource; }; } // end namespace itk diff --git a/Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilter.hxx b/Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilter.hxx index 5a352bd3817..b06d6629c0a 100644 --- a/Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilter.hxx +++ b/Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilter.hxx @@ -24,7 +24,9 @@ #include "itkProgressAccumulator.h" #include "itkImageAlgorithm.h" +#include "itkImageRegionIteratorWithIndex.h" #include "itkMacro.h" +#include "itkVariableLengthVector.h" namespace itk { @@ -72,6 +74,97 @@ FFTDiscreteGaussianImageFilter::SetInputBoundaryCondi itkWarningMacro("FFTDiscreteGaussianImageFilter ignores InputBoundaryCondition, use RealBoundaryCondition instead"); } +template +auto +FFTDiscreteGaussianImageFilter::GenerateKernelImage() -> RealImagePointerType +{ + RealImagePointerType kernelImage; + if (m_KernelSource == FFTDiscreteGaussianImageFilterEnums::KernelSource::COMBINED_OPERATORS) + { + // Get directional 1D Gaussian kernels to compose image + itk::VariableLengthVector directionalOperators; + directionalOperators.SetSize(this->GetFilterDimensionality()); + RadiusType kernelSize; + kernelSize.Fill(1); + for (size_t dim = 0; dim < this->GetFilterDimensionality(); ++dim) + { + this->GenerateKernel(dim, directionalOperators[dim]); + kernelSize[dim] = directionalOperators[dim].GetRadius(dim) * 2 + 1; + } + + // Set up kernel image + kernelImage = RealImageType::New(); + typename RealImageType::IndexType index; + index.Fill(0); + typename RealImageType::RegionType region; + region.SetSize(kernelSize); + region.SetIndex(index); + kernelImage->SetRegions(region); + kernelImage->Allocate(); + kernelImage->CopyInformation(this->GetInput()); + + // Compute kernel image as product of vectors + itk::ImageRegionIteratorWithIndex kernelIt(kernelImage, region); + while (!kernelIt.IsAtEnd()) + { + auto imageIndex = kernelIt.GetIndex(); + double val = 1; + for (size_t dim = 0; dim < directionalOperators.GetSize(); ++dim) + { + val *= directionalOperators[dim].GetElement(imageIndex[dim]); + } + kernelIt.Set(val); + + ++kernelIt; + } + } + else if (m_KernelSource == FFTDiscreteGaussianImageFilterEnums::KernelSource::IMAGE_SOURCE) + { + // Create kernel image for blurring in requested dimensions + using GaussianImageSourceType = GaussianImageSource; + using KernelSizeType = typename GaussianImageSourceType::SizeType; + using KernelMeanType = typename GaussianImageSourceType::ArrayType; + + typename GaussianImageSourceType::Pointer kernelSource = GaussianImageSourceType::New(); + + auto inputSpacing = this->GetInput()->GetSpacing(); + auto inputOrigin = this->GetInput()->GetOrigin(); + + kernelSource->SetScale(1.0); + kernelSource->SetNormalized(true); + + kernelSource->SetSpacing(inputSpacing); + kernelSource->SetOrigin(inputOrigin); + kernelSource->SetDirection(this->GetInput()->GetDirection()); + + KernelSizeType kernelSize; + kernelSize.Fill(1); + for (size_t dim = 0; dim < this->GetFilterDimensionality(); ++dim) + { + kernelSize[dim] = static_cast(this->GetKernelRadius(dim)) * 2 + 1; + } + kernelSource->SetSize(kernelSize); + + KernelMeanType mean; + for (size_t dim = 0; dim < ImageDimension; ++dim) + { + double radius = (kernelSize[dim] - 1) / 2; + mean[dim] = inputSpacing[dim] * radius + inputOrigin[dim]; // center pixel pos + } + kernelSource->SetMean(mean); + kernelSource->SetSigma(this->GetSigmaArray()); + kernelSource->Update(); + kernelImage = kernelSource->GetOutput(); + kernelImage->DisconnectPipeline(); + } + else + { + itkExceptionMacro("Unknown kernel source enum"); + } + + return kernelImage; +} + template void FFTDiscreteGaussianImageFilter::GenerateData() @@ -107,53 +200,18 @@ FFTDiscreteGaussianImageFilter::GenerateData() auto progress = ProgressAccumulator::New(); progress->SetMiniPipelineFilter(this); - // Create kernel image for blurring in requested dimensions - using GaussianImageSourceType = GaussianImageSource; - using KernelSizeType = typename GaussianImageSourceType::SizeType; - using KernelMeanType = typename GaussianImageSourceType::ArrayType; - - typename GaussianImageSourceType::Pointer kernelSource = GaussianImageSourceType::New(); - - kernelSource->SetScale(1.0); - kernelSource->SetNormalized(true); - - kernelSource->SetSpacing(localInput->GetSpacing()); - kernelSource->SetOrigin(localInput->GetOrigin()); - kernelSource->SetDirection(localInput->GetDirection()); - - KernelSizeType kernelSize; - kernelSize.Fill(1); - for (size_t dim = 0; dim < filterDimensionality; ++dim) - { - kernelSize[dim] = static_cast(this->GetKernelRadius(dim)) * 2 + 1; - } - kernelSource->SetSize(kernelSize); - - KernelMeanType mean; - auto inputSpacing = localInput->GetSpacing(); - auto inputOrigin = localInput->GetOrigin(); - for (size_t dim = 0; dim < ImageDimension; ++dim) - { - double radius = (kernelSize[dim] - 1) / 2; - mean[dim] = inputSpacing[dim] * radius + inputOrigin[dim]; // center pixel pos - } - kernelSource->SetMean(mean); - kernelSource->SetSigma(this->GetSigmaArray()); - - // Estimate kernel generation to be roughly 5% of effort - progress->RegisterInternalFilter(kernelSource, 0.05f); + RealImagePointerType kernelImage = GenerateKernelImage(); // Perform image convolution by FFT using ConvolutionImageFilterType = FFTConvolutionImageFilter; typename ConvolutionImageFilterType::Pointer convolutionFilter = ConvolutionImageFilterType::New(); convolutionFilter->SetInput(this->GetInput()); - convolutionFilter->SetKernelImage(kernelSource->GetOutput()); + convolutionFilter->SetKernelImage(kernelImage); convolutionFilter->SetBoundaryCondition(this->GetRealBoundaryCondition()); convolutionFilter->SetNormalize(false); // Kernel is already normalized - // Estimate kernel generation to be roughly 95% of effort - progress->RegisterInternalFilter(convolutionFilter, 0.95f); + progress->RegisterInternalFilter(convolutionFilter, 1.0f); // Graft this filters output onto the mini-pipeline so the mini-pipeline // has the correct region ivars and will write to this filters bulk data diff --git a/Modules/Filtering/Smoothing/src/CMakeLists.txt b/Modules/Filtering/Smoothing/src/CMakeLists.txt index 5f672a06b1b..07cda34113d 100644 --- a/Modules/Filtering/Smoothing/src/CMakeLists.txt +++ b/Modules/Filtering/Smoothing/src/CMakeLists.txt @@ -1,4 +1,5 @@ set(ITKSmoothing_SRCS + itkFFTDiscreteGaussianImageFilter.cxx itkRecursiveGaussianImageFilter.cxx ) itk_module_add_library(ITKSmoothing ${ITKSmoothing_SRCS}) diff --git a/Modules/Filtering/Smoothing/src/itkFFTDiscreteGaussianImageFilter.cxx b/Modules/Filtering/Smoothing/src/itkFFTDiscreteGaussianImageFilter.cxx new file mode 100644 index 00000000000..b34c22ee591 --- /dev/null +++ b/Modules/Filtering/Smoothing/src/itkFFTDiscreteGaussianImageFilter.cxx @@ -0,0 +1,38 @@ +/*========================================================================= + * + * 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 + * + * http://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 "itkFFTDiscreteGaussianImageFilter.h" + +namespace itk +{ +/** Print enum values */ +std::ostream & +operator<<(std::ostream & out, const FFTDiscreteGaussianImageFilterEnums::KernelSource value) +{ + return out << [value] { + switch (value) + { + case FFTDiscreteGaussianImageFilterEnums::KernelSource::COMBINED_OPERATORS: + return "itk::FFTDiscreteGaussianImageFilterEnums::KernelSource::COMBINED_OPERATORS"; + case FFTDiscreteGaussianImageFilterEnums::KernelSource::IMAGE_SOURCE: + return "itk::FFTDiscreteGaussianImageFilterEnums::KernelSource::IMAGE_SOURCE"; + default: + return "INVALID VALUE FOR itk::FFTDiscreteGaussianImageFilterEnums::KernelSource"; + } + }(); +} +} // namespace itk diff --git a/Modules/Filtering/Smoothing/test/Baseline/DiscreteGaussianImageFilterComparisonTestOutput.mha.sha512 b/Modules/Filtering/Smoothing/test/Baseline/DiscreteGaussianImageFilterComparisonTestOutput.mha.sha512 new file mode 100644 index 00000000000..d07c0909f41 --- /dev/null +++ b/Modules/Filtering/Smoothing/test/Baseline/DiscreteGaussianImageFilterComparisonTestOutput.mha.sha512 @@ -0,0 +1 @@ +1c08c1470588ceb69a58a2e7612f44b6a7d147d61383b37299a73c45a37e1cc58726194ae344a30c47eb8d4b50ee3dea97805da15150412038e084a768723978 diff --git a/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterTestImageSourceOutput.mha.sha512 b/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterTestImageSourceOutput.mha.sha512 new file mode 100644 index 00000000000..a6c997629b7 --- /dev/null +++ b/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterTestImageSourceOutput.mha.sha512 @@ -0,0 +1 @@ +240f7a6dcd5c2a59867200314ff4b0ec7c78d9ca4c4ec27e2949d3829c430066092c1c4c1fd6f395ab818d7fb8e9575a3487ca73e1e3f48a783aa6cc5b6e28a8 diff --git a/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterTestKernelDimensionalityOutput.mha.sha512 b/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterTestKernelDimensionalityOutput.mha.sha512 new file mode 100644 index 00000000000..ea8cedd2252 --- /dev/null +++ b/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterTestKernelDimensionalityOutput.mha.sha512 @@ -0,0 +1 @@ +15c62fb6743991414dee0cb1340b23f42c14bde50a7976fa65c0d47be694c1c1793e9fbac86adc3e29b62f61bd266feed5643b37d3204abe75bd5d7f5a93dffc diff --git a/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterTestOutput.mha.sha512 b/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterTestOutput.mha.sha512 deleted file mode 100644 index 8f1ea54c767..00000000000 --- a/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterTestOutput.mha.sha512 +++ /dev/null @@ -1 +0,0 @@ -0eeb11ed3467d7af854c70506ea3e476fb73f57aaa1779ba3eb8abbfc30035bf4b0c11a7e47d12c65f01f1682dcac11a01e8ac20f1a0398ddb5f47f668a24a75 diff --git a/Modules/Filtering/Smoothing/test/CMakeLists.txt b/Modules/Filtering/Smoothing/test/CMakeLists.txt index 3b4f995dc1b..7d6a43790fd 100644 --- a/Modules/Filtering/Smoothing/test/CMakeLists.txt +++ b/Modules/Filtering/Smoothing/test/CMakeLists.txt @@ -43,7 +43,10 @@ itk_add_test(NAME itkDiscreteGaussianImageFilterTest2 COMMAND ITKSmoothingTestDriver --compare DATA{${ITK_DATA_ROOT}/Baseline/BasicFilters/DiscreteGaussianImageFilterTest2_OutputA.mha} ${ITK_TEST_OUTPUT_DIR}/DiscreteGaussianImageFilterTest2_OutputA.mha - itkDiscreteGaussianImageFilterTest2 2 3 DATA{${ITK_DATA_ROOT}/Input/RGBTestImage.tif} 3.5 ${ITK_TEST_OUTPUT_DIR}/DiscreteGaussianImageFilterTest2_OutputA.mha) + itkDiscreteGaussianImageFilterTest2 2 3 + DATA{${ITK_DATA_ROOT}/Input/RGBTestImage.tif} + ${ITK_TEST_OUTPUT_DIR}/DiscreteGaussianImageFilterTest2_OutputA.mha + 1.870829) itk_add_test(NAME itkSmoothingRecursiveGaussianImageFilterTestNormalizeScalesOn COMMAND ITKSmoothingTestDriver --compare-MD5 ${ITK_TEST_OUTPUT_DIR}/itkSmoothingRecursiveGaussianImageFilterTestNormalizeScalesOn.png @@ -61,15 +64,43 @@ itk_add_test(NAME itkDiscreteGaussianImageFilterTest1a COMMAND ITKSmoothingTestDriver itkDiscreteGaussianImageFilterTest 1) itk_add_test(NAME itkDiscreteGaussianImageFilterTest1b COMMAND ITKSmoothingTestDriver itkDiscreteGaussianImageFilterTest 0) -itk_add_test(NAME itkFFTDiscreteGaussianImageFilterTest + +# Use equivalent input parameters to compare standard and FFT +# procedures to a common baseline for equivalent output +itk_add_test(NAME itkDiscreteGaussianImageFilterComparisonTest + COMMAND ITKSmoothingTestDriver + --compare DATA{Baseline/DiscreteGaussianImageFilterComparisonTestOutput.mha} + ${ITK_TEST_OUTPUT_DIR}/DiscreteGaussianImageFilterComparisonTestOutput.mha + itkDiscreteGaussianImageFilterTest2 2 1 + DATA{${ITK_DATA_ROOT}/Input/2th_cthead1.png} + ${ITK_TEST_OUTPUT_DIR}/DiscreteGaussianImageFilterComparisonTestOutput.mha + 3.5 0.05 35) +itk_add_test(NAME itkFFTDiscreteGaussianImageFilterComparisonTest COMMAND ITKSmoothingTestDriver - --compare DATA{Baseline/FFTDiscreteGaussianImageFilterTestOutput.mha} - ${ITK_TEST_OUTPUT_DIR}/FFTDiscreteGaussianImageFilterTestOutput.mha - itkFFTDiscreteGaussianImageFilterTest + --compare DATA{Baseline/DiscreteGaussianImageFilterComparisonTestOutput.mha} + ${ITK_TEST_OUTPUT_DIR}/FFTDiscreteGaussianImageFilterComparisonTestOutput.mha + itkFFTDiscreteGaussianImageFilterTest DATA{${ITK_DATA_ROOT}/Input/2th_cthead1.png} - 3.5 - 35 - ${ITK_TEST_OUTPUT_DIR}/FFTDiscreteGaussianImageFilterTestOutput.mha) + ${ITK_TEST_OUTPUT_DIR}/FFTDiscreteGaussianImageFilterComparisonTestOutput.mha + 3.5 0.05 35 2 0) + +itk_add_test(NAME itkFFTDiscreteGaussianImageFilterImageSourceTest + COMMAND ITKSmoothingTestDriver + --compare DATA{Baseline/FFTDiscreteGaussianImageFilterTestImageSourceOutput.mha} + ${ITK_TEST_OUTPUT_DIR}/FFTDiscreteGaussianImageFilterTestImageSourceOutput.mha + itkFFTDiscreteGaussianImageFilterTest + DATA{${ITK_DATA_ROOT}/Input/2th_cthead1.png} + ${ITK_TEST_OUTPUT_DIR}/FFTDiscreteGaussianImageFilterTestImageSourceOutput.mha + 3.5 0.05 35 2 1) +itk_add_test(NAME itkFFTDiscreteGaussianImageFilterKernelDimensionalityTest + COMMAND ITKSmoothingTestDriver + --compare DATA{Baseline/FFTDiscreteGaussianImageFilterTestKernelDimensionalityOutput.mha} + ${ITK_TEST_OUTPUT_DIR}/FFTDiscreteGaussianImageFilterTestKernelDimensionalityOutput.mha + itkFFTDiscreteGaussianImageFilterTest + DATA{${ITK_DATA_ROOT}/Input/2th_cthead1.png} + ${ITK_TEST_OUTPUT_DIR}/FFTDiscreteGaussianImageFilterTestKernelDimensionalityOutput.mha + 1.2 0.25 11 1 0) + itk_add_test(NAME itkFFTDiscreteGaussianImageFilterFactoryTest COMMAND ITKSmoothingTestDriver itkFFTDiscreteGaussianImageFilterFactoryTest) itk_add_test(NAME itkMedianImageFilterTest diff --git a/Modules/Filtering/Smoothing/test/itkDiscreteGaussianImageFilterTest.cxx b/Modules/Filtering/Smoothing/test/itkDiscreteGaussianImageFilterTest.cxx index 6caf8d8526f..c0971ba9d79 100644 --- a/Modules/Filtering/Smoothing/test/itkDiscreteGaussianImageFilterTest.cxx +++ b/Modules/Filtering/Smoothing/test/itkDiscreteGaussianImageFilterTest.cxx @@ -23,6 +23,8 @@ #include "itkTestingMacros.h" #include "itkConstantBoundaryCondition.h" +/** Check basic image filter parameters and operations */ + int itkDiscreteGaussianImageFilterTest(int argc, char * argv[]) { @@ -71,10 +73,6 @@ itkDiscreteGaussianImageFilterTest(int argc, char * argv[]) filter->SetSigmaArray(sigma); ITK_TEST_SET_GET_VALUE(sigma, filter->GetSigmaArray()); - // Verify spacing-adjusted variance is not available without input image - filter->SetUseImageSpacing(true); - ITK_TRY_EXPECT_EXCEPTION(filter->GetImageSpacingVarianceArray()); - // Set the boundary condition used by the filter itk::ConstantBoundaryCondition constantBoundaryCondition; filter->SetInputBoundaryCondition(&constantBoundaryCondition); @@ -128,29 +126,6 @@ itkDiscreteGaussianImageFilterTest(int argc, char * argv[]) ITK_TRY_EXPECT_NO_EXCEPTION(test1.Execute()); - // Test variance adjustments for image spacing - ImageType::Pointer inputImage = ImageType::New(); - typename ImageType::SpacingType spacing; - spacing[0] = 1.5; - spacing[1] = 1.0; - spacing[2] = 0.5; - inputImage->SetSpacing(spacing); - filter->SetInput(inputImage); - - filter->SetUseImageSpacing(false); - for (itk::IndexValueType dim = 0; dim < Dimension; ++dim) - { - ITK_TEST_EXPECT_EQUAL(sigmaValue, filter->GetVariance()[dim]) - ITK_TEST_EXPECT_EQUAL(sigmaValue, filter->GetImageSpacingVarianceArray()[dim]); - } - - filter->SetUseImageSpacing(true); - for (itk::IndexValueType dim = 0; dim < Dimension; ++dim) - { - ITK_TEST_EXPECT_EQUAL(sigmaValue, filter->GetVariance()[dim]); - ITK_TEST_EXPECT_EQUAL((sigmaValue / (spacing[dim] * spacing[dim])), filter->GetImageSpacingVarianceArray()[dim]); - } - std::cout << "Test finished" << std::endl; return EXIT_SUCCESS; } diff --git a/Modules/Filtering/Smoothing/test/itkDiscreteGaussianImageFilterTest2.cxx b/Modules/Filtering/Smoothing/test/itkDiscreteGaussianImageFilterTest2.cxx index 7d8f5fd04a5..d4544f76b64 100644 --- a/Modules/Filtering/Smoothing/test/itkDiscreteGaussianImageFilterTest2.cxx +++ b/Modules/Filtering/Smoothing/test/itkDiscreteGaussianImageFilterTest2.cxx @@ -22,9 +22,15 @@ #include "itkImageFileWriter.h" #include "itkTestingMacros.h" +/** Execute the filter on user-provided input */ + template int -itkDiscreteGaussianImageFilterTestA(const char * inputFilename, float sigma, const char * outputFilename) +itkDiscreteGaussianImageFilterTestA(const char * inputFilename, + const char * outputFilename, + const float sigma, + const float kernelError, + const float kernelWidth) { using ImageType = TIMAGE; @@ -36,7 +42,9 @@ itkDiscreteGaussianImageFilterTestA(const char * inputFilename, float sigma, con auto filter = FilterType::New(); filter->SetInput(reader->GetOutput()); - filter->SetVariance(sigma); + filter->SetSigma(sigma); + filter->SetMaximumError(kernelError); + filter->SetMaximumKernelWidth(kernelWidth); filter->Update(); using WriterType = itk::ImageFileWriter; @@ -53,30 +61,33 @@ itkDiscreteGaussianImageFilterTestA(const char * inputFilename, float sigma, con int itkDiscreteGaussianImageFilterTest2(int argc, char * argv[]) { - if (argc != 6) + if (argc < 5) { std::cerr << "Missing parameters." << std::endl; std::cerr << "Usage:" << std::endl; - std::cerr << itkNameOfTestExecutableMacro(argv) - << " imageDimension vectorDimension inputFilename sigma outputFilename" << std::endl; + std::cerr << itkNameOfTestExecutableMacro(argv) << " imageDimension vectorDimension inputFilename outputFilename" + << " [sigma] [kernelError] [kernelWidth] " << std::endl; return EXIT_FAILURE; } unsigned int img_dim = std::stoi(argv[1]); - unsigned int vec_dim = std::stoi(argv[2]); - if (img_dim != 2) { std::cerr << "This test only supports 2D image for demo! exiting ..." << std::endl; return EXIT_FAILURE; } + unsigned int vec_dim = std::stoi(argv[2]); if (vec_dim != 1 && vec_dim != 3) { std::cerr << "This test only supports 3-channel image or 1-channel image for demo! Exiting ... " << std::endl; return EXIT_FAILURE; } + float sigma = (argc >= 6) ? std::stof(argv[5]) : 0.0; + float kernelError = (argc >= 7) ? std::stof(argv[6]) : 0.01; + unsigned int kernelWidth = (argc >= 8) ? static_cast(std::stoi(argv[7])) : 32; + using ScalarPixelType = float; using ScalarImageType = itk::Image; using VectorPixelType = itk::Vector; @@ -85,10 +96,10 @@ itkDiscreteGaussianImageFilterTest2(int argc, char * argv[]) switch (vec_dim) { case 1: - itkDiscreteGaussianImageFilterTestA(argv[3], std::stod(argv[4]), argv[5]); + itkDiscreteGaussianImageFilterTestA(argv[3], argv[4], sigma, kernelError, kernelWidth); break; case 3: - itkDiscreteGaussianImageFilterTestA(argv[3], std::stod(argv[4]), argv[5]); + itkDiscreteGaussianImageFilterTestA(argv[3], argv[4], sigma, kernelError, kernelWidth); break; } diff --git a/Modules/Filtering/Smoothing/test/itkFFTDiscreteGaussianImageFilterTest.cxx b/Modules/Filtering/Smoothing/test/itkFFTDiscreteGaussianImageFilterTest.cxx index 2ded25decb7..7ab74326abf 100644 --- a/Modules/Filtering/Smoothing/test/itkFFTDiscreteGaussianImageFilterTest.cxx +++ b/Modules/Filtering/Smoothing/test/itkFFTDiscreteGaussianImageFilterTest.cxx @@ -32,7 +32,8 @@ itkFFTDiscreteGaussianImageFilterTest(int argc, char * argv[]) std::cerr << "Missing parameters." << std::endl; std::cerr << "Usage:" << std::endl; std::cerr << itkNameOfTestExecutableMacro(argv) - << " inputFilename sigma kernelWidth outputFilename filterDimensionality" << std::endl; + << " inputFilename outputFilename sigma kernelError kernelWidth filterDimensionality kernelSource" + << std::endl; return EXIT_FAILURE; } @@ -40,9 +41,11 @@ itkFFTDiscreteGaussianImageFilterTest(int argc, char * argv[]) const size_t ImageDimension = 2; using ImageType = itk::Image; - double sigma = std::stod(argv[2]); - unsigned int kernelWidth = std::stoi(argv[3]); - unsigned int filterDimensionality = (argc < 6 ? ImageDimension : std::stoi(argv[5])); + float sigma = (argc >= 4) ? std::stof(argv[3]) : 0.0; + float kernelError = (argc >= 5) ? std::stof(argv[4]) : 0.01; + unsigned int kernelWidth = (argc >= 6) ? std::stoi(argv[5]) : 32; + unsigned int filterDimensionality = (argc >= 7) ? std::stoi(argv[6]) : ImageDimension; + unsigned int kernelSource = (argc >= 8) ? std::stoi(argv[7]) : ImageDimension; typename ImageType::Pointer inputImage = itk::ReadImage(argv[1]); @@ -54,35 +57,17 @@ itkFFTDiscreteGaussianImageFilterTest(int argc, char * argv[]) filter->SetInput(inputImage); - filter->SetVariance(sigma); - for (auto & param : filter->GetVariance()) + filter->SetSigma(sigma); + for (auto & param : filter->GetSigmaArray()) { ITK_TEST_EXPECT_EQUAL(sigma, param); } - - filter->SetUseImageSpacingOff(); - for (auto & param : filter->GetImageSpacingVarianceArray()) - { - ITK_TEST_EXPECT_EQUAL(sigma, param); - } - - using SpacingType = typename ImageType::SpacingType; - SpacingType spacing; - spacing[0] = 0.5; - spacing[1] = 0.25; - inputImage->SetSpacing(spacing); - filter->SetUseImageSpacingOn(); - for (size_t dim = 0; dim < ImageDimension; ++dim) + for (auto & param : filter->GetVariance()) { - ITK_TEST_EXPECT_EQUAL(sigma, filter->GetVariance()[dim]) - ITK_TEST_EXPECT_EQUAL(sigma / (spacing[dim] * spacing[dim]), filter->GetImageSpacingVarianceArray()[dim]); + double tolerance = 1e-6; + ITK_TEST_EXPECT_TRUE(std::fabs((sigma * sigma) - param) < tolerance); } - spacing[0] = 1.0; - spacing[1] = 1.0; - inputImage->SetSpacing(spacing); - - const float kernelError = 0.05; filter->SetMaximumError(kernelError); for (size_t dim = 0; dim < ImageDimension; ++dim) { @@ -99,6 +84,18 @@ itkFFTDiscreteGaussianImageFilterTest(int argc, char * argv[]) filter->SetRealBoundaryCondition(&zfnBoundaryCondition); ITK_TEST_SET_GET_VALUE(&zfnBoundaryCondition, filter->GetRealBoundaryCondition()); + itk::FFTDiscreteGaussianImageFilterEnums::KernelSource source; + if (kernelSource == 0) + { + source = itk::FFTDiscreteGaussianImageFilterEnums::KernelSource::COMBINED_OPERATORS; + } + else + { + source = itk::FFTDiscreteGaussianImageFilterEnums::KernelSource::IMAGE_SOURCE; + } + filter->SetKernelSource(source); + ITK_TEST_SET_GET_VALUE(filter->GetKernelSource(), source); + // Test that attempting to set unused interface parameters throws warnings // Ignore values from "Get" methods as these are unused @@ -109,7 +106,7 @@ itkFFTDiscreteGaussianImageFilterTest(int argc, char * argv[]) ITK_TRY_EXPECT_NO_EXCEPTION(filter->Update()); - itk::WriteImage(filter->GetOutput(), argv[4], true); + itk::WriteImage(filter->GetOutput(), argv[2], true); return EXIT_SUCCESS; } diff --git a/Modules/Filtering/Smoothing/wrapping/itkFFTDiscreteGaussianImageFilter.wrap b/Modules/Filtering/Smoothing/wrapping/itkFFTDiscreteGaussianImageFilter.wrap index 6c638893481..7992fd06f48 100644 --- a/Modules/Filtering/Smoothing/wrapping/itkFFTDiscreteGaussianImageFilter.wrap +++ b/Modules/Filtering/Smoothing/wrapping/itkFFTDiscreteGaussianImageFilter.wrap @@ -1,3 +1,8 @@ +set(WRAPPER_AUTO_INCLUDE_HEADERS OFF) +itk_wrap_include("itkFFTDiscreteGaussianImageFilter.h") + +itk_wrap_simple_class("itk::FFTDiscreteGaussianImageFilterEnums") + itk_wrap_class("itk::FFTDiscreteGaussianImageFilter" POINTER) itk_wrap_image_filter("${WRAP_ITK_SCALAR}" 2) itk_end_wrap_class() From d1e7c9bac5ce2e9286f62690906f585bd208b7a4 Mon Sep 17 00:00:00 2001 From: Tom Birdsong Date: Wed, 26 Jan 2022 11:23:03 -0500 Subject: [PATCH 5/5] ENH: Increase test coverage for `FFTDiscreteGaussianImageFilter` --- .../itkFFTDiscreteGaussianImageFilter.h | 31 ++++++-- .../itkFFTDiscreteGaussianImageFilter.hxx | 40 +++++----- .../src/itkFFTDiscreteGaussianImageFilter.cxx | 4 +- ...ageFilter3DComparisonTestOutput.mha.sha512 | 1 + ...geFilter3DComparisonTestOutput2.mha.sha512 | 1 + ...onunitImageComparisonTestOutput.mha.sha512 | 1 + ...ageFilter3DComparisonTestOutput.mha.sha512 | 1 + ...geFilter3DComparisonTestOutput2.mha.sha512 | 1 + ...geFilter3DImageSourceTestOutput.mha.sha512 | 1 + ...anImageFilter3DVolumeTestOutput.mha.sha512 | 1 + ...ImageFilterComparisonTestOutput.mha.sha512 | 1 + ...onunitImageComparisonTestOutput.mha.sha512 | 1 + ...rTestKernelDimensionalityOutput.mha.sha512 | 2 +- ...teGaussianImageFilterTestOutput.mha.sha512 | 1 + .../Filtering/Smoothing/test/CMakeLists.txt | 74 ++++++++++++++++--- ...ageFilterTestNonNullOriginInput.mha.sha512 | 1 + .../itkDiscreteGaussianImageFilterTest2.cxx | 52 ++++++++----- .../itkFFTDiscreteGaussianImageFilterTest.cxx | 67 +++++++++++------ 18 files changed, 195 insertions(+), 86 deletions(-) create mode 100644 Modules/Filtering/Smoothing/test/Baseline/DiscreteGaussianImageFilter3DComparisonTestOutput.mha.sha512 create mode 100644 Modules/Filtering/Smoothing/test/Baseline/DiscreteGaussianImageFilter3DComparisonTestOutput2.mha.sha512 create mode 100644 Modules/Filtering/Smoothing/test/Baseline/DiscreteGaussianImageFilterNonunitImageComparisonTestOutput.mha.sha512 create mode 100644 Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilter3DComparisonTestOutput.mha.sha512 create mode 100644 Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilter3DComparisonTestOutput2.mha.sha512 create mode 100644 Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilter3DImageSourceTestOutput.mha.sha512 create mode 100644 Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilter3DVolumeTestOutput.mha.sha512 create mode 100644 Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterComparisonTestOutput.mha.sha512 create mode 100644 Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterNonunitImageComparisonTestOutput.mha.sha512 create mode 100644 Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterTestOutput.mha.sha512 create mode 100644 Modules/Filtering/Smoothing/test/Input/itkDiscreteGaussianImageFilterTestNonNullOriginInput.mha.sha512 diff --git a/Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilter.h b/Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilter.h index 3332b9adc10..5b282797a26 100644 --- a/Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilter.h +++ b/Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilter.h @@ -19,6 +19,7 @@ #define itkFFTDiscreteGaussianImageFilter_h #include "itkDiscreteGaussianImageFilter.h" +#include "itkFFTConvolutionImageFilter.h" #include "itkImage.h" #include "itkMacro.h" #include "ITKSmoothingExport.h" @@ -49,7 +50,7 @@ class FFTDiscreteGaussianImageFilterEnums */ enum class KernelSource : uint8_t { - COMBINED_OPERATORS = 0, + OPERATORS = 0, IMAGE_SOURCE, }; }; @@ -120,7 +121,6 @@ class ITK_TEMPLATE_EXPORT FFTDiscreteGaussianImageFilter : public DiscreteGaussi using typename Superclass::RealOutputPixelValueType; using RealPixelType = RealOutputPixelType; using RealImageType = RealOutputImageType; - using RealImagePointerType = typename RealImageType::Pointer; /** Typedef to describe the boundary condition. */ using typename Superclass::BoundaryConditionType; @@ -138,23 +138,24 @@ class ITK_TEMPLATE_EXPORT FFTDiscreteGaussianImageFilter : public DiscreteGaussi using typename Superclass::KernelType; using typename Superclass::RadiusType; + /** Typedef for convolution */ + using ConvolutionImageFilterType = FFTConvolutionImageFilter; + /** Overridden accessors for unused parameters */ void SetInputBoundaryCondition(const InputBoundaryConditionPointerType) override; - auto - GenerateKernelImage() -> RealImagePointerType; - - itkSetMacro(KernelSource, FFTDiscreteGaussianImageFilterEnums::KernelSource); itkGetConstMacro(KernelSource, FFTDiscreteGaussianImageFilterEnums::KernelSource); + itkGetConstMacro(KernelImage, typename RealImageType::Pointer); + protected: FFTDiscreteGaussianImageFilter() = default; ~FFTDiscreteGaussianImageFilter() override = default; - // Pad input region to kernel image size + /** Pad input region to kernel image size */ void GenerateInputRequestedRegion() override; @@ -166,8 +167,22 @@ class ITK_TEMPLATE_EXPORT FFTDiscreteGaussianImageFilter : public DiscreteGaussi void GenerateData() override; + /** Create a kernel image matching user input specifications. + * Returns a reference to the member m_KernelImage. */ + auto + GenerateKernelImage() -> RealImageType *; + private: - FFTDiscreteGaussianImageFilterEnums::KernelSource m_KernelSource; + /* Enum for choosing among available kernel generation methods */ + FFTDiscreteGaussianImageFilterEnums::KernelSource m_KernelSource = + FFTDiscreteGaussianImageFilterEnums::KernelSource::OPERATORS; + + /* Kernel image is allocated with GenerateKernelImage() */ + typename RealImageType::Pointer m_KernelImage; + + /* Persist mini-pipeline filter to minimize construction costs + * on repeated calls to GenerateData() */ + typename ConvolutionImageFilterType::Pointer m_ConvolutionImageFilter = ConvolutionImageFilterType::New(); }; } // end namespace itk diff --git a/Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilter.hxx b/Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilter.hxx index b06d6629c0a..e3fbe1b5333 100644 --- a/Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilter.hxx +++ b/Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilter.hxx @@ -30,7 +30,6 @@ namespace itk { - template void FFTDiscreteGaussianImageFilter::GenerateInputRequestedRegion() @@ -76,10 +75,10 @@ FFTDiscreteGaussianImageFilter::SetInputBoundaryCondi template auto -FFTDiscreteGaussianImageFilter::GenerateKernelImage() -> RealImagePointerType +FFTDiscreteGaussianImageFilter::GenerateKernelImage() -> RealImageType * { - RealImagePointerType kernelImage; - if (m_KernelSource == FFTDiscreteGaussianImageFilterEnums::KernelSource::COMBINED_OPERATORS) + m_KernelImage = RealImageType::New(); + if (m_KernelSource == FFTDiscreteGaussianImageFilterEnums::KernelSource::OPERATORS) { // Get directional 1D Gaussian kernels to compose image itk::VariableLengthVector directionalOperators; @@ -93,18 +92,17 @@ FFTDiscreteGaussianImageFilter::GenerateKernelImage() } // Set up kernel image - kernelImage = RealImageType::New(); typename RealImageType::IndexType index; index.Fill(0); typename RealImageType::RegionType region; region.SetSize(kernelSize); region.SetIndex(index); - kernelImage->SetRegions(region); - kernelImage->Allocate(); - kernelImage->CopyInformation(this->GetInput()); + m_KernelImage->SetRegions(region); + m_KernelImage->Allocate(); + m_KernelImage->CopyInformation(this->GetInput()); // Compute kernel image as product of vectors - itk::ImageRegionIteratorWithIndex kernelIt(kernelImage, region); + itk::ImageRegionIteratorWithIndex kernelIt(m_KernelImage, region); while (!kernelIt.IsAtEnd()) { auto imageIndex = kernelIt.GetIndex(); @@ -154,15 +152,15 @@ FFTDiscreteGaussianImageFilter::GenerateKernelImage() kernelSource->SetMean(mean); kernelSource->SetSigma(this->GetSigmaArray()); kernelSource->Update(); - kernelImage = kernelSource->GetOutput(); - kernelImage->DisconnectPipeline(); + m_KernelImage = kernelSource->GetOutput(); + m_KernelImage->DisconnectPipeline(); } else { itkExceptionMacro("Unknown kernel source enum"); } - return kernelImage; + return m_KernelImage.GetPointer(); } template @@ -200,26 +198,24 @@ FFTDiscreteGaussianImageFilter::GenerateData() auto progress = ProgressAccumulator::New(); progress->SetMiniPipelineFilter(this); - RealImagePointerType kernelImage = GenerateKernelImage(); + RealImageType * kernelImage = GenerateKernelImage(); // Perform image convolution by FFT - using ConvolutionImageFilterType = FFTConvolutionImageFilter; - typename ConvolutionImageFilterType::Pointer convolutionFilter = ConvolutionImageFilterType::New(); - convolutionFilter->SetInput(this->GetInput()); - convolutionFilter->SetKernelImage(kernelImage); - convolutionFilter->SetBoundaryCondition(this->GetRealBoundaryCondition()); - convolutionFilter->SetNormalize(false); // Kernel is already normalized + m_ConvolutionImageFilter->SetInput(this->GetInput()); + m_ConvolutionImageFilter->SetKernelImage(kernelImage); + m_ConvolutionImageFilter->SetBoundaryCondition(this->GetRealBoundaryCondition()); + m_ConvolutionImageFilter->SetNormalize(false); // Kernel is already normalized - progress->RegisterInternalFilter(convolutionFilter, 1.0f); + progress->RegisterInternalFilter(m_ConvolutionImageFilter, 1.0f); // Graft this filters output onto the mini-pipeline so the mini-pipeline // has the correct region ivars and will write to this filters bulk data // output. - convolutionFilter->GraftOutput(output); + m_ConvolutionImageFilter->GraftOutput(output); // Update the last filter in the mini-pipeline - convolutionFilter->Update(); + m_ConvolutionImageFilter->Update(); // Graft the last output of the mini-pipeline onto this filters output so // the final output has the correct region ivars and a handle to the final diff --git a/Modules/Filtering/Smoothing/src/itkFFTDiscreteGaussianImageFilter.cxx b/Modules/Filtering/Smoothing/src/itkFFTDiscreteGaussianImageFilter.cxx index b34c22ee591..fc518ce2fa6 100644 --- a/Modules/Filtering/Smoothing/src/itkFFTDiscreteGaussianImageFilter.cxx +++ b/Modules/Filtering/Smoothing/src/itkFFTDiscreteGaussianImageFilter.cxx @@ -26,8 +26,8 @@ operator<<(std::ostream & out, const FFTDiscreteGaussianImageFilterEnums::Kernel return out << [value] { switch (value) { - case FFTDiscreteGaussianImageFilterEnums::KernelSource::COMBINED_OPERATORS: - return "itk::FFTDiscreteGaussianImageFilterEnums::KernelSource::COMBINED_OPERATORS"; + case FFTDiscreteGaussianImageFilterEnums::KernelSource::OPERATORS: + return "itk::FFTDiscreteGaussianImageFilterEnums::KernelSource::OPERATORS"; case FFTDiscreteGaussianImageFilterEnums::KernelSource::IMAGE_SOURCE: return "itk::FFTDiscreteGaussianImageFilterEnums::KernelSource::IMAGE_SOURCE"; default: diff --git a/Modules/Filtering/Smoothing/test/Baseline/DiscreteGaussianImageFilter3DComparisonTestOutput.mha.sha512 b/Modules/Filtering/Smoothing/test/Baseline/DiscreteGaussianImageFilter3DComparisonTestOutput.mha.sha512 new file mode 100644 index 00000000000..8e5ccfa470f --- /dev/null +++ b/Modules/Filtering/Smoothing/test/Baseline/DiscreteGaussianImageFilter3DComparisonTestOutput.mha.sha512 @@ -0,0 +1 @@ +b0a82461048deba245eff84a0d174028fca230393fe87393b376bc2b59ea25bce540c6f036d5b62a063939eb1e6170d94aafd5e06b783dc111cf900d905ac165 diff --git a/Modules/Filtering/Smoothing/test/Baseline/DiscreteGaussianImageFilter3DComparisonTestOutput2.mha.sha512 b/Modules/Filtering/Smoothing/test/Baseline/DiscreteGaussianImageFilter3DComparisonTestOutput2.mha.sha512 new file mode 100644 index 00000000000..d7219996290 --- /dev/null +++ b/Modules/Filtering/Smoothing/test/Baseline/DiscreteGaussianImageFilter3DComparisonTestOutput2.mha.sha512 @@ -0,0 +1 @@ +fa76505faae56b05bbee62ceeac1b959bb5cbdda555795a6b2f7763301a449077d9348a7642f647e6b9a17ceabd136b1658dbe632acec989deafb9b1bd2ea6ad diff --git a/Modules/Filtering/Smoothing/test/Baseline/DiscreteGaussianImageFilterNonunitImageComparisonTestOutput.mha.sha512 b/Modules/Filtering/Smoothing/test/Baseline/DiscreteGaussianImageFilterNonunitImageComparisonTestOutput.mha.sha512 new file mode 100644 index 00000000000..57f84bcde8e --- /dev/null +++ b/Modules/Filtering/Smoothing/test/Baseline/DiscreteGaussianImageFilterNonunitImageComparisonTestOutput.mha.sha512 @@ -0,0 +1 @@ +685f2ad5b2f79b50d64e496a1fb631e2e6c66ef95625ba1fe007f64f4708d6092a2a72c9adc7813273e283120ff6e9676b65c607ca7897a77b2a180f9840596b diff --git a/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilter3DComparisonTestOutput.mha.sha512 b/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilter3DComparisonTestOutput.mha.sha512 new file mode 100644 index 00000000000..609a6ad4694 --- /dev/null +++ b/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilter3DComparisonTestOutput.mha.sha512 @@ -0,0 +1 @@ +6ddf34b759488a8f4cfb75b4f7f49a0f09d234bd3c0c033167274d3d2b58043d1981112e7020d79ef28c57a91f5817b897845071ae69c715a21d816fd162d699 diff --git a/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilter3DComparisonTestOutput2.mha.sha512 b/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilter3DComparisonTestOutput2.mha.sha512 new file mode 100644 index 00000000000..d0e09eb0842 --- /dev/null +++ b/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilter3DComparisonTestOutput2.mha.sha512 @@ -0,0 +1 @@ +33074175fb567145ff5714606c8bddecf5f45d5dd11b51a63bd83f406f79a785c41afc1d8864a2b5d82780c62eeef073aa8243676516daf1d300ba6dfd796385 diff --git a/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilter3DImageSourceTestOutput.mha.sha512 b/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilter3DImageSourceTestOutput.mha.sha512 new file mode 100644 index 00000000000..a7406a36093 --- /dev/null +++ b/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilter3DImageSourceTestOutput.mha.sha512 @@ -0,0 +1 @@ +afd93c5868e8d240550bac56b5d121a9fdc4e4d9802a5de45854c73cc6991da96a5a4d0f9e2e415fd621762ed31754adcc91cd1d6b94f437b7fdbb9530146600 diff --git a/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilter3DVolumeTestOutput.mha.sha512 b/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilter3DVolumeTestOutput.mha.sha512 new file mode 100644 index 00000000000..a7406a36093 --- /dev/null +++ b/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilter3DVolumeTestOutput.mha.sha512 @@ -0,0 +1 @@ +afd93c5868e8d240550bac56b5d121a9fdc4e4d9802a5de45854c73cc6991da96a5a4d0f9e2e415fd621762ed31754adcc91cd1d6b94f437b7fdbb9530146600 diff --git a/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterComparisonTestOutput.mha.sha512 b/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterComparisonTestOutput.mha.sha512 new file mode 100644 index 00000000000..a138fdc17a1 --- /dev/null +++ b/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterComparisonTestOutput.mha.sha512 @@ -0,0 +1 @@ +51ae3a28f4bd6a922b6a15f138bf610fef68cefc3f1ee0da28428249a151f35a993238f49ee3da8825a12c5637210f56194dbbf8b2df29a07b7aa024a7163372 diff --git a/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterNonunitImageComparisonTestOutput.mha.sha512 b/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterNonunitImageComparisonTestOutput.mha.sha512 new file mode 100644 index 00000000000..d40dea1e625 --- /dev/null +++ b/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterNonunitImageComparisonTestOutput.mha.sha512 @@ -0,0 +1 @@ +8bce8d04ab5b9c44371bd254d5bdba38443fb8e9ee975eab30b440bcc7dbc69dc9acc41deba2e94df2f2a45bb3b1b0e259313c897bcefd203d2f4b2dcbf20b5c diff --git a/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterTestKernelDimensionalityOutput.mha.sha512 b/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterTestKernelDimensionalityOutput.mha.sha512 index ea8cedd2252..ea7c3809688 100644 --- a/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterTestKernelDimensionalityOutput.mha.sha512 +++ b/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterTestKernelDimensionalityOutput.mha.sha512 @@ -1 +1 @@ -15c62fb6743991414dee0cb1340b23f42c14bde50a7976fa65c0d47be694c1c1793e9fbac86adc3e29b62f61bd266feed5643b37d3204abe75bd5d7f5a93dffc +c458ca85b1a891bae65d2e5d6e684a0b60593bbf18a7f6bf756ab47d42464d3fa73908e9fbf197e08eea32872972b84676973bc9b178aa9c47225b761c25b96f diff --git a/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterTestOutput.mha.sha512 b/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterTestOutput.mha.sha512 new file mode 100644 index 00000000000..a6c997629b7 --- /dev/null +++ b/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterTestOutput.mha.sha512 @@ -0,0 +1 @@ +240f7a6dcd5c2a59867200314ff4b0ec7c78d9ca4c4ec27e2949d3829c430066092c1c4c1fd6f395ab818d7fb8e9575a3487ca73e1e3f48a783aa6cc5b6e28a8 diff --git a/Modules/Filtering/Smoothing/test/CMakeLists.txt b/Modules/Filtering/Smoothing/test/CMakeLists.txt index 7d6a43790fd..8d0b4e3623f 100644 --- a/Modules/Filtering/Smoothing/test/CMakeLists.txt +++ b/Modules/Filtering/Smoothing/test/CMakeLists.txt @@ -64,10 +64,9 @@ itk_add_test(NAME itkDiscreteGaussianImageFilterTest1a COMMAND ITKSmoothingTestDriver itkDiscreteGaussianImageFilterTest 1) itk_add_test(NAME itkDiscreteGaussianImageFilterTest1b COMMAND ITKSmoothingTestDriver itkDiscreteGaussianImageFilterTest 0) - # Use equivalent input parameters to compare standard and FFT # procedures to a common baseline for equivalent output -itk_add_test(NAME itkDiscreteGaussianImageFilterComparisonTest +itk_add_test(NAME itkDiscreteGaussianImageFilter2DComparisonTest COMMAND ITKSmoothingTestDriver --compare DATA{Baseline/DiscreteGaussianImageFilterComparisonTestOutput.mha} ${ITK_TEST_OUTPUT_DIR}/DiscreteGaussianImageFilterComparisonTestOutput.mha @@ -75,20 +74,70 @@ itk_add_test(NAME itkDiscreteGaussianImageFilterComparisonTest DATA{${ITK_DATA_ROOT}/Input/2th_cthead1.png} ${ITK_TEST_OUTPUT_DIR}/DiscreteGaussianImageFilterComparisonTestOutput.mha 3.5 0.05 35) -itk_add_test(NAME itkFFTDiscreteGaussianImageFilterComparisonTest +itk_add_test(NAME itkFFTDiscreteGaussianImageFilter2DComparisonTest COMMAND ITKSmoothingTestDriver --compare DATA{Baseline/DiscreteGaussianImageFilterComparisonTestOutput.mha} ${ITK_TEST_OUTPUT_DIR}/FFTDiscreteGaussianImageFilterComparisonTestOutput.mha - itkFFTDiscreteGaussianImageFilterTest + itkFFTDiscreteGaussianImageFilterTest 2 DATA{${ITK_DATA_ROOT}/Input/2th_cthead1.png} ${ITK_TEST_OUTPUT_DIR}/FFTDiscreteGaussianImageFilterComparisonTestOutput.mha 3.5 0.05 35 2 0) - +itk_add_test(NAME itkDiscreteGaussianImageFilter3DComparisonTest + COMMAND ITKSmoothingTestDriver + --compare DATA{Baseline/DiscreteGaussianImageFilter3DComparisonTestOutput.mha} + ${ITK_TEST_OUTPUT_DIR}/DiscreteGaussianImageFilter3DComparisonTestOutput.mha + itkDiscreteGaussianImageFilterTest2 3 1 + DATA{${ITK_DATA_ROOT}/Input/HeadMRVolume.mha} + ${ITK_TEST_OUTPUT_DIR}/DiscreteGaussianImageFilter3DComparisonTestOutput.mha + 5.0 0.03 25) +itk_add_test(NAME itkFFTDiscreteGaussianImageFilter3DComparisonTest + COMMAND ITKSmoothingTestDriver + --compare DATA{Baseline/DiscreteGaussianImageFilter3DComparisonTestOutput.mha} + ${ITK_TEST_OUTPUT_DIR}/FFTDiscreteGaussianImageFilter3DComparisonTestOutput.mha + itkFFTDiscreteGaussianImageFilterTest 3 + DATA{${ITK_DATA_ROOT}/Input/HeadMRVolume.mha} + ${ITK_TEST_OUTPUT_DIR}/FFTDiscreteGaussianImageFilter3DComparisonTestOutput.mha + 5.0 0.03 25 3 0) +# Compare results on images with adjusted origin and spacing +itk_add_test(NAME itkDiscreteGaussianImageFilterNonunitImageComparisonTest + COMMAND ITKSmoothingTestDriver + --compare DATA{Baseline/DiscreteGaussianImageFilterNonunitImageComparisonTestOutput.mha} + ${ITK_TEST_OUTPUT_DIR}/DiscreteGaussianImageFilterNonunitImageComparisonTestOutput.mha + itkDiscreteGaussianImageFilterTest2 2 1 + DATA{Input/itkDiscreteGaussianImageFilterTestNonNullOriginInput.mha} + ${ITK_TEST_OUTPUT_DIR}/DiscreteGaussianImageFilterNonunitImageComparisonTestOutput.mha + 10.0 0.03 31 2) +itk_add_test(NAME itkFFTDiscreteGaussianImageFilterNonunitImageComparisonTest + COMMAND ITKSmoothingTestDriver + --compare DATA{Baseline/DiscreteGaussianImageFilterNonunitImageComparisonTestOutput.mha} + ${ITK_TEST_OUTPUT_DIR}/FFTDiscreteGaussianImageFilterNonunitImageComparisonTestOutput.mha + itkFFTDiscreteGaussianImageFilterTest 2 + DATA{Input/itkDiscreteGaussianImageFilterTestNonNullOriginInput.mha} + ${ITK_TEST_OUTPUT_DIR}/FFTDiscreteGaussianImageFilterNonunitImageComparisonTestOutput.mha + 10.0 0.03 31 2 0) +# Comparing results of lower-dimensionality smoothing +itk_add_test(NAME itkDiscreteGaussianImageFilter3DComparisonTest2 + COMMAND ITKSmoothingTestDriver + --compare DATA{Baseline/DiscreteGaussianImageFilter3DComparisonTestOutput2.mha} + ${ITK_TEST_OUTPUT_DIR}/DiscreteGaussianImageFilter3DComparisonTestOutput2.mha + itkDiscreteGaussianImageFilterTest2 3 1 + DATA{${ITK_DATA_ROOT}/Input/HeadMRVolume.mha} + ${ITK_TEST_OUTPUT_DIR}/DiscreteGaussianImageFilter3DComparisonTestOutput2.mha + 4.0 0.02 31 2) +itk_add_test(NAME itkFFTDiscreteGaussianImageFilter3DComparisonTest2 + COMMAND ITKSmoothingTestDriver + --compare DATA{Baseline/DiscreteGaussianImageFilter3DComparisonTestOutput2.mha} + ${ITK_TEST_OUTPUT_DIR}/FFTDiscreteGaussianImageFilter3DComparisonTestOutput2.mha + itkFFTDiscreteGaussianImageFilterTest 3 + DATA{${ITK_DATA_ROOT}/Input/HeadMRVolume.mha} + ${ITK_TEST_OUTPUT_DIR}/FFTDiscreteGaussianImageFilter3DComparisonTestOutput2.mha + 4.0 0.03 31 2 0) +# Test smoothing with GaussianImageSource itk_add_test(NAME itkFFTDiscreteGaussianImageFilterImageSourceTest COMMAND ITKSmoothingTestDriver --compare DATA{Baseline/FFTDiscreteGaussianImageFilterTestImageSourceOutput.mha} ${ITK_TEST_OUTPUT_DIR}/FFTDiscreteGaussianImageFilterTestImageSourceOutput.mha - itkFFTDiscreteGaussianImageFilterTest + itkFFTDiscreteGaussianImageFilterTest 2 DATA{${ITK_DATA_ROOT}/Input/2th_cthead1.png} ${ITK_TEST_OUTPUT_DIR}/FFTDiscreteGaussianImageFilterTestImageSourceOutput.mha 3.5 0.05 35 2 1) @@ -96,11 +145,18 @@ itk_add_test(NAME itkFFTDiscreteGaussianImageFilterKernelDimensionalityTest COMMAND ITKSmoothingTestDriver --compare DATA{Baseline/FFTDiscreteGaussianImageFilterTestKernelDimensionalityOutput.mha} ${ITK_TEST_OUTPUT_DIR}/FFTDiscreteGaussianImageFilterTestKernelDimensionalityOutput.mha - itkFFTDiscreteGaussianImageFilterTest + itkFFTDiscreteGaussianImageFilterTest 2 DATA{${ITK_DATA_ROOT}/Input/2th_cthead1.png} ${ITK_TEST_OUTPUT_DIR}/FFTDiscreteGaussianImageFilterTestKernelDimensionalityOutput.mha - 1.2 0.25 11 1 0) - + 1.2 0.25 11 1 1) +itk_add_test(NAME itkFFTDiscreteGaussianImageFilter3DImageSourceTest + COMMAND ITKSmoothingTestDriver + --compare DATA{Baseline/FFTDiscreteGaussianImageFilter3DImageSourceTestOutput.mha} + ${ITK_TEST_OUTPUT_DIR}/FFTDiscreteGaussianImageFilter3DImageSourceTestOutput.mha + itkFFTDiscreteGaussianImageFilterTest 3 + DATA{${ITK_DATA_ROOT}/Input/HeadMRVolume.mha} + ${ITK_TEST_OUTPUT_DIR}/FFTDiscreteGaussianImageFilter3DImageSourceTestOutput.mha + 6.0 0.10 25 3 1) itk_add_test(NAME itkFFTDiscreteGaussianImageFilterFactoryTest COMMAND ITKSmoothingTestDriver itkFFTDiscreteGaussianImageFilterFactoryTest) itk_add_test(NAME itkMedianImageFilterTest diff --git a/Modules/Filtering/Smoothing/test/Input/itkDiscreteGaussianImageFilterTestNonNullOriginInput.mha.sha512 b/Modules/Filtering/Smoothing/test/Input/itkDiscreteGaussianImageFilterTestNonNullOriginInput.mha.sha512 new file mode 100644 index 00000000000..9a96601c48d --- /dev/null +++ b/Modules/Filtering/Smoothing/test/Input/itkDiscreteGaussianImageFilterTestNonNullOriginInput.mha.sha512 @@ -0,0 +1 @@ +a67a85418bde40c5dfb25966b9a9dbe14ee5f2e2a43388c4404d7edcb5b7239de121c7740fb751da158937624cac6a0ffb3bbfcac8563d8a0a8f1381c02c857c diff --git a/Modules/Filtering/Smoothing/test/itkDiscreteGaussianImageFilterTest2.cxx b/Modules/Filtering/Smoothing/test/itkDiscreteGaussianImageFilterTest2.cxx index d4544f76b64..333d406e051 100644 --- a/Modules/Filtering/Smoothing/test/itkDiscreteGaussianImageFilterTest2.cxx +++ b/Modules/Filtering/Smoothing/test/itkDiscreteGaussianImageFilterTest2.cxx @@ -26,11 +26,12 @@ template int -itkDiscreteGaussianImageFilterTestA(const char * inputFilename, - const char * outputFilename, - const float sigma, - const float kernelError, - const float kernelWidth) +itkDiscreteGaussianImageFilterTestA(const char * inputFilename, + const char * outputFilename, + const float sigma, + const float kernelError, + const unsigned int kernelWidth, + const unsigned int filterDimensionality) { using ImageType = TIMAGE; @@ -45,6 +46,7 @@ itkDiscreteGaussianImageFilterTestA(const char * inputFilename, filter->SetSigma(sigma); filter->SetMaximumError(kernelError); filter->SetMaximumKernelWidth(kernelWidth); + filter->SetFilterDimensionality(filterDimensionality); filter->Update(); using WriterType = itk::ImageFileWriter; @@ -66,14 +68,14 @@ itkDiscreteGaussianImageFilterTest2(int argc, char * argv[]) std::cerr << "Missing parameters." << std::endl; std::cerr << "Usage:" << std::endl; std::cerr << itkNameOfTestExecutableMacro(argv) << " imageDimension vectorDimension inputFilename outputFilename" - << " [sigma] [kernelError] [kernelWidth] " << std::endl; + << " [sigma] [kernelError] [kernelWidth] [filterDimensionality]" << std::endl; return EXIT_FAILURE; } unsigned int img_dim = std::stoi(argv[1]); - if (img_dim != 2) + if (img_dim < 2 || img_dim > 3) { - std::cerr << "This test only supports 2D image for demo! exiting ..." << std::endl; + std::cerr << "This test only supports 2D or 3D images for demo! exiting ..." << std::endl; return EXIT_FAILURE; } @@ -84,23 +86,33 @@ itkDiscreteGaussianImageFilterTest2(int argc, char * argv[]) return EXIT_FAILURE; } - float sigma = (argc >= 6) ? std::stof(argv[5]) : 0.0; - float kernelError = (argc >= 7) ? std::stof(argv[6]) : 0.01; - unsigned int kernelWidth = (argc >= 8) ? static_cast(std::stoi(argv[7])) : 32; + float sigma = (argc > 5) ? std::stof(argv[5]) : 0.0; + float kernelError = (argc > 6) ? std::stof(argv[6]) : 0.01; + unsigned int kernelWidth = (argc > 7) ? static_cast(std::stoi(argv[7])) : 32; + unsigned int filterDimensionality = (argc > 8) ? static_cast(std::stoi(argv[8])) : img_dim; using ScalarPixelType = float; - using ScalarImageType = itk::Image; using VectorPixelType = itk::Vector; - using VectorImageType = itk::Image; - switch (vec_dim) + if (img_dim == 2 && vec_dim == 1) { - case 1: - itkDiscreteGaussianImageFilterTestA(argv[3], argv[4], sigma, kernelError, kernelWidth); - break; - case 3: - itkDiscreteGaussianImageFilterTestA(argv[3], argv[4], sigma, kernelError, kernelWidth); - break; + itkDiscreteGaussianImageFilterTestA>( + argv[3], argv[4], sigma, kernelError, kernelWidth, filterDimensionality); + } + else if (img_dim == 2 && vec_dim == 3) + { + itkDiscreteGaussianImageFilterTestA>( + argv[3], argv[4], sigma, kernelError, kernelWidth, filterDimensionality); + } + else if (img_dim == 3 && vec_dim == 1) + { + itkDiscreteGaussianImageFilterTestA>( + argv[3], argv[4], sigma, kernelError, kernelWidth, filterDimensionality); + } + else if (img_dim == 3 && vec_dim == 1) + { + itkDiscreteGaussianImageFilterTestA>( + argv[3], argv[4], sigma, kernelError, kernelWidth, filterDimensionality); } return EXIT_SUCCESS; diff --git a/Modules/Filtering/Smoothing/test/itkFFTDiscreteGaussianImageFilterTest.cxx b/Modules/Filtering/Smoothing/test/itkFFTDiscreteGaussianImageFilterTest.cxx index 7ab74326abf..7ae6fa4a5f2 100644 --- a/Modules/Filtering/Smoothing/test/itkFFTDiscreteGaussianImageFilterTest.cxx +++ b/Modules/Filtering/Smoothing/test/itkFFTDiscreteGaussianImageFilterTest.cxx @@ -24,30 +24,17 @@ #include "itkTestingMacros.h" #include "itkZeroFluxNeumannBoundaryCondition.h" +template int -itkFFTDiscreteGaussianImageFilterTest(int argc, char * argv[]) +itkFFTDiscreteGaussianImageFilterTestProcedure(int argc, char ** argv) { - if (argc < 5) - { - std::cerr << "Missing parameters." << std::endl; - std::cerr << "Usage:" << std::endl; - std::cerr << itkNameOfTestExecutableMacro(argv) - << " inputFilename outputFilename sigma kernelError kernelWidth filterDimensionality kernelSource" - << std::endl; - return EXIT_FAILURE; - } - - using ScalarPixelType = float; - const size_t ImageDimension = 2; - using ImageType = itk::Image; + float sigma = (argc > 4) ? std::stof(argv[4]) : 0.0; + float kernelError = (argc > 5) ? std::stof(argv[5]) : 0.01; + unsigned int kernelWidth = (argc > 6) ? std::stoi(argv[6]) : 32; + unsigned int filterDimensionality = (argc > 7) ? std::stoi(argv[7]) : ImageType::ImageDimension; + unsigned int kernelSource = (argc > 8) ? std::stoi(argv[8]) : 0; - float sigma = (argc >= 4) ? std::stof(argv[3]) : 0.0; - float kernelError = (argc >= 5) ? std::stof(argv[4]) : 0.01; - unsigned int kernelWidth = (argc >= 6) ? std::stoi(argv[5]) : 32; - unsigned int filterDimensionality = (argc >= 7) ? std::stoi(argv[6]) : ImageDimension; - unsigned int kernelSource = (argc >= 8) ? std::stoi(argv[7]) : ImageDimension; - - typename ImageType::Pointer inputImage = itk::ReadImage(argv[1]); + typename ImageType::Pointer inputImage = itk::ReadImage(argv[2]); using FilterType = itk::FFTDiscreteGaussianImageFilter; auto filter = FilterType::New(); @@ -69,7 +56,7 @@ itkFFTDiscreteGaussianImageFilterTest(int argc, char * argv[]) } filter->SetMaximumError(kernelError); - for (size_t dim = 0; dim < ImageDimension; ++dim) + for (size_t dim = 0; dim < ImageType::ImageDimension; ++dim) { ITK_TEST_SET_GET_VALUE(kernelError, filter->GetMaximumError()[dim]); } @@ -87,7 +74,7 @@ itkFFTDiscreteGaussianImageFilterTest(int argc, char * argv[]) itk::FFTDiscreteGaussianImageFilterEnums::KernelSource source; if (kernelSource == 0) { - source = itk::FFTDiscreteGaussianImageFilterEnums::KernelSource::COMBINED_OPERATORS; + source = itk::FFTDiscreteGaussianImageFilterEnums::KernelSource::OPERATORS; } else { @@ -106,7 +93,39 @@ itkFFTDiscreteGaussianImageFilterTest(int argc, char * argv[]) ITK_TRY_EXPECT_NO_EXCEPTION(filter->Update()); - itk::WriteImage(filter->GetOutput(), argv[2], true); + itk::WriteImage(filter->GetOutput(), argv[3], true); + + return EXIT_SUCCESS; +} +int +itkFFTDiscreteGaussianImageFilterTest(int argc, char * argv[]) +{ + if (argc < 5) + { + std::cerr << "Missing parameters." << std::endl; + std::cerr << "Usage:" << std::endl; + std::cerr + << itkNameOfTestExecutableMacro(argv) + << " imageDimension inputFilename outputFilename sigma kernelError kernelWidth filterDimensionality kernelSource" + << std::endl; + return EXIT_FAILURE; + } + + const unsigned int ImageDimension = static_cast(std::stoi(argv[1])); + + if (ImageDimension == 2) + { + itkFFTDiscreteGaussianImageFilterTestProcedure>(argc, &argv[0]); + } + else if (ImageDimension == 3) + { + itkFFTDiscreteGaussianImageFilterTestProcedure>(argc, &argv[0]); + } + else + { + std::cout << "Did not recognize image dimension argument!" << std::endl; + return EXIT_FAILURE; + } return EXIT_SUCCESS; }