diff --git a/Modules/Filtering/Smoothing/include/itkDiscreteGaussianImageFilter.h b/Modules/Filtering/Smoothing/include/itkDiscreteGaussianImageFilter.h index 7d856781dd2..d9da5834660 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 @@ -293,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 @@ -326,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 @@ -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/include/itkFFTDiscreteGaussianImageFilter.h b/Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilter.h new file mode 100644 index 00000000000..5b282797a26 --- /dev/null +++ b/Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilter.h @@ -0,0 +1,193 @@ +/*========================================================================= + * + * 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 "itkFFTConvolutionImageFilter.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 + { + 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 + * 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 typename Superclass::RealOutputPixelType; + using typename Superclass::RealOutputImageType; + using typename Superclass::RealOutputPixelValueType; + using RealPixelType = RealOutputPixelType; + using RealImageType = RealOutputImageType; + + /** Typedef to describe the boundary condition. */ + 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 typename Superclass::ArrayType; + using typename Superclass::SigmaArrayType; + using typename Superclass::ScalarRealType; + + /** Typedef to describe kernel parameters */ + using typename Superclass::KernelType; + using typename Superclass::RadiusType; + + /** Typedef for convolution */ + using ConvolutionImageFilterType = FFTConvolutionImageFilter; + + /** Overridden accessors for unused parameters */ + + void + SetInputBoundaryCondition(const InputBoundaryConditionPointerType) override; + + 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 */ + 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; + + /** Create a kernel image matching user input specifications. + * Returns a reference to the member m_KernelImage. */ + auto + GenerateKernelImage() -> RealImageType *; + +private: + /* 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 + +#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..e3fbe1b5333 --- /dev/null +++ b/Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilter.hxx @@ -0,0 +1,227 @@ +/*========================================================================= + * + * 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 "itkImageRegionIteratorWithIndex.h" +#include "itkMacro.h" +#include "itkVariableLengthVector.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 +auto +FFTDiscreteGaussianImageFilter::GenerateKernelImage() -> RealImageType * +{ + m_KernelImage = RealImageType::New(); + if (m_KernelSource == FFTDiscreteGaussianImageFilterEnums::KernelSource::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 + typename RealImageType::IndexType index; + index.Fill(0); + typename RealImageType::RegionType region; + region.SetSize(kernelSize); + region.SetIndex(index); + m_KernelImage->SetRegions(region); + m_KernelImage->Allocate(); + m_KernelImage->CopyInformation(this->GetInput()); + + // Compute kernel image as product of vectors + itk::ImageRegionIteratorWithIndex kernelIt(m_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(); + m_KernelImage = kernelSource->GetOutput(); + m_KernelImage->DisconnectPipeline(); + } + else + { + itkExceptionMacro("Unknown kernel source enum"); + } + + return m_KernelImage.GetPointer(); +} + +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); + + RealImageType * kernelImage = GenerateKernelImage(); + + // Perform image convolution by FFT + + m_ConvolutionImageFilter->SetInput(this->GetInput()); + m_ConvolutionImageFilter->SetKernelImage(kernelImage); + m_ConvolutionImageFilter->SetBoundaryCondition(this->GetRealBoundaryCondition()); + m_ConvolutionImageFilter->SetNormalize(false); // Kernel is already normalized + + 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. + m_ConvolutionImageFilter->GraftOutput(output); + + // Update the last filter in the mini-pipeline + 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 + // bulk data + this->GraftOutput(output); +} +} // namespace itk + +#endif 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/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/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..fc518ce2fa6 --- /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::OPERATORS: + return "itk::FFTDiscreteGaussianImageFilterEnums::KernelSource::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/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/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/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/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..ea7c3809688 --- /dev/null +++ b/Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterTestKernelDimensionalityOutput.mha.sha512 @@ -0,0 +1 @@ +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 db9c444cff2..8d0b4e3623f 100644 --- a/Modules/Filtering/Smoothing/test/CMakeLists.txt +++ b/Modules/Filtering/Smoothing/test/CMakeLists.txt @@ -3,6 +3,8 @@ set(ITKSmoothingTests itkBoxMeanImageFilterTest.cxx itkBoxSigmaImageFilterTest.cxx itkDiscreteGaussianImageFilterTest2.cxx +itkFFTDiscreteGaussianImageFilterTest.cxx +itkFFTDiscreteGaussianImageFilterFactoryTest.cxx itkSmoothingRecursiveGaussianImageFilterTest.cxx itkSmoothingRecursiveGaussianImageFilterOnVectorImageTest.cxx itkSmoothingRecursiveGaussianImageFilterOnImageOfVectorTest.cxx @@ -41,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 @@ -59,6 +64,101 @@ 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 itkDiscreteGaussianImageFilter2DComparisonTest + 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 itkFFTDiscreteGaussianImageFilter2DComparisonTest + COMMAND ITKSmoothingTestDriver + --compare DATA{Baseline/DiscreteGaussianImageFilterComparisonTestOutput.mha} + ${ITK_TEST_OUTPUT_DIR}/FFTDiscreteGaussianImageFilterComparisonTestOutput.mha + 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 2 + 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 2 + DATA{${ITK_DATA_ROOT}/Input/2th_cthead1.png} + ${ITK_TEST_OUTPUT_DIR}/FFTDiscreteGaussianImageFilterTestKernelDimensionalityOutput.mha + 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 COMMAND ITKSmoothingTestDriver itkMedianImageFilterTest) itk_add_test(NAME itkRecursiveGaussianImageFiltersOnTensorsTest 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/itkDiscreteGaussianImageFilterTest.cxx b/Modules/Filtering/Smoothing/test/itkDiscreteGaussianImageFilterTest.cxx index 09f731c4d65..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[]) { @@ -92,7 +94,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,7 +126,6 @@ itkDiscreteGaussianImageFilterTest(int argc, char * argv[]) ITK_TRY_EXPECT_NO_EXCEPTION(test1.Execute()); - 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..333d406e051 100644 --- a/Modules/Filtering/Smoothing/test/itkDiscreteGaussianImageFilterTest2.cxx +++ b/Modules/Filtering/Smoothing/test/itkDiscreteGaussianImageFilterTest2.cxx @@ -22,9 +22,16 @@ #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 unsigned int kernelWidth, + const unsigned int filterDimensionality) { using ImageType = TIMAGE; @@ -36,7 +43,10 @@ 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->SetFilterDimensionality(filterDimensionality); filter->Update(); using WriterType = itk::ImageFileWriter; @@ -53,43 +63,56 @@ 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] [filterDimensionality]" << 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) + 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; } + 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 > 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) + { + 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) { - case 1: - itkDiscreteGaussianImageFilterTestA(argv[3], std::stod(argv[4]), argv[5]); - break; - case 3: - itkDiscreteGaussianImageFilterTestA(argv[3], std::stod(argv[4]), argv[5]); - break; + itkDiscreteGaussianImageFilterTestA>( + argv[3], argv[4], sigma, kernelError, kernelWidth, filterDimensionality); } return EXIT_SUCCESS; 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/test/itkFFTDiscreteGaussianImageFilterTest.cxx b/Modules/Filtering/Smoothing/test/itkFFTDiscreteGaussianImageFilterTest.cxx new file mode 100644 index 00000000000..7ae6fa4a5f2 --- /dev/null +++ b/Modules/Filtering/Smoothing/test/itkFFTDiscreteGaussianImageFilterTest.cxx @@ -0,0 +1,131 @@ +/*========================================================================= + * + * 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" + +template +int +itkFFTDiscreteGaussianImageFilterTestProcedure(int argc, char ** argv) +{ + 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; + + typename ImageType::Pointer inputImage = itk::ReadImage(argv[2]); + + using FilterType = itk::FFTDiscreteGaussianImageFilter; + auto filter = FilterType::New(); + ITK_EXERCISE_BASIC_OBJECT_METHODS(filter, FFTDiscreteGaussianImageFilter, DiscreteGaussianImageFilter); + + // Test setting inputs + + filter->SetInput(inputImage); + + filter->SetSigma(sigma); + for (auto & param : filter->GetSigmaArray()) + { + ITK_TEST_EXPECT_EQUAL(sigma, param); + } + for (auto & param : filter->GetVariance()) + { + double tolerance = 1e-6; + ITK_TEST_EXPECT_TRUE(std::fabs((sigma * sigma) - param) < tolerance); + } + + filter->SetMaximumError(kernelError); + for (size_t dim = 0; dim < ImageType::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()); + + itk::FFTDiscreteGaussianImageFilterEnums::KernelSource source; + if (kernelSource == 0) + { + source = itk::FFTDiscreteGaussianImageFilterEnums::KernelSource::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 + + itk::ConstantBoundaryCondition constantBoundaryCondition; + filter->SetInputBoundaryCondition(&constantBoundaryCondition); + + // Run convolution + + ITK_TRY_EXPECT_NO_EXCEPTION(filter->Update()); + + 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; +} diff --git a/Modules/Filtering/Smoothing/wrapping/itkFFTDiscreteGaussianImageFilter.wrap b/Modules/Filtering/Smoothing/wrapping/itkFFTDiscreteGaussianImageFilter.wrap new file mode 100644 index 00000000000..7992fd06f48 --- /dev/null +++ b/Modules/Filtering/Smoothing/wrapping/itkFFTDiscreteGaussianImageFilter.wrap @@ -0,0 +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() 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)