diff --git a/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/CMakeLists.txt b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/CMakeLists.txt new file mode 100644 index 00000000000..a403dc68615 --- /dev/null +++ b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/CMakeLists.txt @@ -0,0 +1,3 @@ +project(SmoothingRecursiveYvvGaussianFilter) + +itk_module_impl() diff --git a/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/README.md b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/README.md new file mode 100644 index 00000000000..a46431a06b8 --- /dev/null +++ b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/README.md @@ -0,0 +1,29 @@ +# SmoothingRecursiveYvvGaussianFilter + +In-tree ITK module providing the Vliet/Young/Verbeek recursive +Gaussian smoothing approximation as a drop-in alternative to the +canonical `itk::SmoothingRecursiveGaussianImageFilter`. The flagship +classes are +`itk::SmoothingRecursiveYvvGaussianImageFilter`, +`itk::RecursiveLineYvvGaussianImageFilter`, and the +`ITK_USE_GPU`-gated +`itk::GPUSmoothingRecursiveYvvGaussianImageFilter`. + +## Origin + +Ingested from the standalone remote module +[**InsightSoftwareConsortium/ITKSmoothingRecursiveYvvGaussianFilter**](https://github.com/InsightSoftwareConsortium/ITKSmoothingRecursiveYvvGaussianFilter) +on 2026-05-08, following the v4 ingestion guidelines defined in +InsightSoftwareConsortium/ITK#6204. The upstream repository will be +archived read-only after this PR merges; it remains reachable at the +URL above for historical reference (notably the `doc/` directory and +the upstream's bespoke top-level OpenCL detection, which were +intentionally not ingested — ITK's in-tree GPU support is provided +by `ITKGPUCommon`). + +## References + +- Bouilhol G., Pop C., Sarrut D. + *Recursive Implementation of Vliet/Young/Verbeek Gaussian Smoothing.* + The Insight Journal. January-December. 2013. + diff --git a/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/include/itkGPUSmoothingRecursiveYvvGaussianImageFilter.h b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/include/itkGPUSmoothingRecursiveYvvGaussianImageFilter.h new file mode 100644 index 00000000000..998fcebcd47 --- /dev/null +++ b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/include/itkGPUSmoothingRecursiveYvvGaussianImageFilter.h @@ -0,0 +1,211 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +#ifndef itkGPUSmoothingRecursiveYvvGaussianImageFilter_h +#define itkGPUSmoothingRecursiveYvvGaussianImageFilter_h + +#ifdef GPU + +# include "itkImage.h" +# include "itkPixelTraits.h" +# include "itkCommand.h" +# include "itkFixedArray.h" +# include "itkSmoothingRecursiveYvvGaussianImageFilter.h" +# include "itkOpenCLUtil.h" +# include "itkGPUImageToImageFilter.h" +# include + +namespace itk +{ +/** + * \class GPUSmoothingRecursiveYvvGaussianImageFilter + * \brief Recursive Gaussian blur based on Young-Van Vliet's algorithm and + * implemented for GPU. + * + * The GPU implementation is more efficient than the CPU + * the bigger the image is; use the benchmark tests to establish the size for which + * the GPU performs better for your particular hardware configuration. + * + * More information in the Insight Journal publication: + * https://hdl.handle.net/10380/3425 + * + * \ingroup SmoothingRecursiveYvvGaussianFilter + */ + +itkGPUKernelClassMacro(GPUSmoothingRecursiveYvvGaussianImageFilterKernel); + +template +class ITK_EXPORT GPUSmoothingRecursiveYvvGaussianImageFilter + : public GPUImageToImageFilter> +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(GPUSmoothingRecursiveYvvGaussianImageFilter); + + /** Standard class type alias. */ + using Self = GPUSmoothingRecursiveYvvGaussianImageFilter; + // typedef SmoothingRecursiveYvvGaussianImageFilter + // CPUSuperclass; + using Superclass = GPUImageToImageFilter>; + using GPUSuperclass = GPUImageToImageFilter>; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Pixel Type of the input image */ + using InputImageType = TInputImage; + using OutputImageType = TOutputImage; + using PixelType = typename TInputImage::PixelType; +# ifdef WITH_DOUBLE + using RealType = typename NumericTraits::RealType; + using ScalarRealType = typename NumericTraits::ScalarRealType; +# else + using RealType = typename NumericTraits::FloatType; + using ScalarRealType = typename NumericTraits::FloatType; +# endif + + using GPUInputImage = typename itk::GPUTraits::Type; + using GPUOutputImage = typename itk::GPUTraits::Type; + + /** Runtime information support. */ + itkTypeMacro(GPUSmoothingRecursiveYvvGaussianImageFilter, GPUImageToImageFilter); + + /** Image dimension. */ + static constexpr unsigned int ImageDimension = TInputImage::ImageDimension; + + /** Define the type for the sigma array */ + using SigmaArrayType = FixedArray; + + /** Define the image type for internal computations + RealType is usually 'double' in NumericTraits. + Here we prefer float in order to save memory. */ + + using InternalRealType = typename NumericTraits::FloatType; + using RealImageType = GPUImage; + + /** Pointer to the Output Image */ + using OutputImagePointer = typename OutputImageType::Pointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Set Sigma value. Sigma is measured in the units of image spacing. You + may use the method SetSigma to set the same value across each axis or + use the method SetSigmaArray if you need different values along each + axis. */ + void + SetSigmaArray(const SigmaArrayType & sigmas); + + void + SetSigma(ScalarRealType sigma); + + SigmaArrayType + GetSigmaArray() const; + + ScalarRealType + GetSigma() const; + + /** Define which normalization factor will be used for the Gaussian */ + void + SetNormalizeAcrossScale(bool normalizeInScaleSpace); + + itkGetConstMacro(NormalizeAcrossScale, bool); + + virtual void + SetUp(ScalarRealType spacing); + +# ifdef ITK_USE_CONCEPT_CHECKING + /** Begin concept checking */ + itkConceptMacro(InputHasNumericTraitsCheck, (Concept::HasNumericTraits)); + /** End concept checking */ +# endif + /** Get OpenCL Kernel source as a string, creates a GetOpenCLSource method */ + itkGetOpenCLSourceFromKernelMacro(GPUSmoothingRecursiveYvvGaussianImageFilterKernel); + void + SetInput(const TInputImage * input) override; + using Superclass::SetInput; + +protected: + GPUSmoothingRecursiveYvvGaussianImageFilter(); + ~GPUSmoothingRecursiveYvvGaussianImageFilter() override = default; + void + PrintSelf(std::ostream & os, Indent indent) const override; + + /** Generate Data */ + void + GPUGenerateData() override; + + /** GPUSmoothingRecursiveYvvGaussianImageFilter needs all of the input to produce an + * output. Therefore, GPUSmoothingRecursiveYvvGaussianImageFilter needs to provide + * an implementation for GenerateInputRequestedRegion in order to inform + * the pipeline execution model. + * \sa ImageToImageFilter::GenerateInputRequestedRegion() */ + void + GenerateInputRequestedRegion() ITK_NOEXCEPT override; + + // Override since the filter produces the entire dataset + void + EnlargeOutputRequestedRegion(DataObject * output) override; + + void + AllocateGPUCoefficients(); + + std::ostringstream defines; + + ScalarRealType m_B1; + ScalarRealType m_B2; + ScalarRealType m_B3; + ScalarRealType m_B; + std::vector m_Bvalues; + + // Initialization matrix for anti-causal pass + std::vector m_CPUMatrix; + + using GPUDataPointer = GPUDataManager::Pointer; + + GPUDataPointer m_GPUMMatrixDataManager; + GPUDataPointer m_GPUBCoefficientsDataManager; + GPUDataPointer m_GPULocalDataManager; + +private: + void + BuildKernel(); + + /** Normalize the image across scale space */ + bool m_NormalizeAcrossScale; + + int m_FilterGPUKernelHandle; + typename GPUInputImage::Pointer inPtr; + typename GPUOutputImage::Pointer otPtr; + typename GPUOutputImage::SizeType m_requestedSize; + /** Standard deviation of the gaussian used for smoothing */ + SigmaArrayType m_Sigma; +}; +} // end namespace itk + +# ifndef ITK_MANUAL_INSTANTIATION +# include "itkGPUSmoothingRecursiveYvvGaussianImageFilter.hxx" +# endif + +#endif // GPU + +#endif // itkGPUSmoothingRecursiveYvvGaussianImageFilter_h diff --git a/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/include/itkGPUSmoothingRecursiveYvvGaussianImageFilter.hxx b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/include/itkGPUSmoothingRecursiveYvvGaussianImageFilter.hxx new file mode 100644 index 00000000000..da4178dc663 --- /dev/null +++ b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/include/itkGPUSmoothingRecursiveYvvGaussianImageFilter.hxx @@ -0,0 +1,359 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +#ifndef itkGPUSmoothingRecursiveYvvGaussianImageFilter_hxx +#define itkGPUSmoothingRecursiveYvvGaussianImageFilter_hxx + +#ifdef GPU + +# include "itkGPUSmoothingRecursiveYvvGaussianImageFilter.h" + +namespace itk +{ +template +GPUSmoothingRecursiveYvvGaussianImageFilter::GPUSmoothingRecursiveYvvGaussianImageFilter() +{ + m_NormalizeAcrossScale = false; + otPtr = dynamic_cast(this->ProcessObject::GetOutput(0)); + + // NB: We must call SetSigma in order to initialize the smoothing + // filters with the default scale. However, m_Sigma must first be + // initialized (it is used inside SetSigma) and it must be >= 0.5 + // or the call will be ignored. + + this->m_Sigma.Fill(0.0); + this->SetSigma(0.1); + m_FilterGPUKernelHandle = -1; +} + +template +void +GPUSmoothingRecursiveYvvGaussianImageFilter::BuildKernel() +{ + const char * GPUSource = GPUSmoothingRecursiveYvvGaussianImageFilterKernel::GetOpenCLSource(); + + typename GPUOutputImage::SizeType lineSizes = inPtr->GetRequestedRegion().GetSize(); + unsigned int maxLineSize = lineSizes[0]; + for (unsigned int d = 1; d < ImageDimension; ++d) + { + maxLineSize = maxLineSize < lineSizes[d] ? lineSizes[d] : maxLineSize; + } + +# ifdef WITH_DOUBLE + defines << "#pragma OPENCL EXTENSION cl_khr_fp64 : enable \n"; +# endif + +# ifdef NVIDIA + defines << "#pragma OPENCL EXTENSION cl_nv_pragma_unroll : enable \n"; +# endif + + defines << "#define DIM_" << TInputImage::ImageDimension << "\n"; + defines << "#define INPIXELTYPE "; + GetTypenameInString(typeid(typename TInputImage::PixelType), defines); + defines << "#define OUTPIXELTYPE "; + GetTypenameInString(typeid(typename TOutputImage::PixelType), defines); + defines << "#define REALTYPE "; + GetTypenameInString(typeid(ScalarRealType), defines); + defines << "#define MAX_LINE_LENGTH " << maxLineSize << "\n"; + + // load, build, create kernel + this->m_GPUKernelManager->LoadProgramFromString(GPUSource, defines.str().c_str()); + m_FilterGPUKernelHandle = this->m_GPUKernelManager->CreateKernel("YvvFilter"); +} + +template +void +GPUSmoothingRecursiveYvvGaussianImageFilter::SetUp(ScalarRealType spacing) +{ + const ScalarRealType sigmad = this->GetSigma() / spacing; + + if (this->GetSigma() >= 0.5) + { + // Compute q according to 16 in Young et al on Gabor filering + ScalarRealType q = 0; + if (sigmad >= 3.556) + { + q = 0.9804 * (sigmad - 3.556) + 2.5091; + } + else + { + if (sigmad < 0.5) + { + std::cerr << "Too low sigma value (< 0.5), computation will not be precise." << std::endl; + } + + q = 0.0561 * sigmad * sigmad + 0.5784 * sigmad - 0.2568; + } + + // Compute B and B0 to B3 according to Young et al 1995 + + ScalarRealType m0 = 1.16680; + ScalarRealType m1 = 1.10783; + ScalarRealType m2 = 1.40586; + ScalarRealType scale = (m0 + q) * (m1 * m1 + m2 * m2 + 2 * m1 * q + q * q); + m_Bvalues.assign(4, ScalarRealType{}); + + m_B1 = m_Bvalues[1] = q * (2 * m0 * m1 + m1 * m1 + m2 * m2 + (2 * m0 + 4 * m1) * q + 3 * q * q) / scale; + m_B2 = m_Bvalues[2] = -q * q * (m0 + 2 * m1 + 3 * q) / scale; + m_B3 = m_Bvalues[3] = q * q * q / scale; + + ScalarRealType baseB = (m0 * (m1 * m1 + m2 * m2)) / scale; + m_B = m_Bvalues[0] = baseB * baseB; + + // M Matrix for initialization on backward pass, from Triggs and Sdika, IEEE + // TSP + m_CPUMatrix.assign(9, ScalarRealType{}); + const ScalarRealType factor = + (1 + m_B1 - m_B2 + m_B3) * (1 - m_B1 - m_B2 - m_B3) * (1 + m_B2 + (m_B1 - m_B3) * m_B3); + + m_CPUMatrix[0] = (-m_B3 * m_B1 + 1 - m_B3 * m_B3 - m_B2) / factor; + m_CPUMatrix[1] = ((m_B3 + m_B1) * (m_B2 + m_B3 * m_B1)) / factor; + m_CPUMatrix[2] = (m_B3 * (m_B1 + m_B3 * m_B2)) / factor; + + m_CPUMatrix[3] = (m_B1 + m_B3 * m_B2) / factor; + m_CPUMatrix[4] = ((1 - m_B2) * (m_B2 + m_B3 * m_B1)) / factor; + m_CPUMatrix[5] = (-m_B3 * (m_B3 * m_B1 + m_B3 * m_B3 + m_B2 - 1)) / factor; + + m_CPUMatrix[6] = (m_B3 * m_B1 + m_B2 + m_B1 * m_B1 - m_B2 * m_B2) / factor; + m_CPUMatrix[7] = + (m_B1 * m_B2 + m_B3 * m_B2 * m_B2 - m_B1 * m_B3 * m_B3 - m_B3 * m_B3 * m_B3 - m_B3 * m_B2 + m_B3) / factor; + m_CPUMatrix[8] = (m_B3 * (m_B1 + m_B3 * m_B2)) / factor; + + if (this->GetDebug()) + { + for (int i = 0; i < 4; ++i) + { + std::cout << "B" << i << " " << m_Bvalues[i] << std::endl; + } + + for (int i = 0; i < 9; ++i) + { + std::cout << "M" << i << " " << m_CPUMatrix[i] << std::endl; + } + } + this->AllocateGPUCoefficients(); + this->Modified(); + } +} + +template +void +GPUSmoothingRecursiveYvvGaussianImageFilter::SetInput(const TInputImage * input) +{ + // ProcessObject is not const_correct so this const_cast is required + ProcessObject::SetNthInput(0, const_cast(input)); + + inPtr = dynamic_cast(this->ProcessObject::GetInput(0)); + inPtr->GetGPUDataManager()->SetBufferFlag(CL_MEM_READ_ONLY); + + m_requestedSize = inPtr->GetRequestedRegion().GetSize(); + + itkDebugMacro(<< "GPUSmoothingRecursiveYvvGaussianImageFilter generating data "); + + for (unsigned int d = 0; d < ImageDimension; d++) + { + if (m_requestedSize[d] < 4) + { + itkExceptionMacro( + "The number of pixels along dimension " + << d << " is less than 4. This filter requires a minimum of four pixels along the dimension to be processed."); + } + } +} + +template +void +GPUSmoothingRecursiveYvvGaussianImageFilter::SetSigma(ScalarRealType sigma) +{ + SigmaArrayType sigmas(sigma); + this->SetSigmaArray(sigmas); +} + +template +void +GPUSmoothingRecursiveYvvGaussianImageFilter::SetSigmaArray(const SigmaArrayType & sigma) +{ + if (this->m_Sigma != sigma) + { + this->m_Sigma = sigma; + + const typename InputImageType::SpacingType & pixelSize = this->GetOutput()->GetSpacing(); + if (pixelSize[0] != 0) + { + this->SetUp(pixelSize[0]); + } + this->Modified(); + } +} + +template +typename GPUSmoothingRecursiveYvvGaussianImageFilter::SigmaArrayType +GPUSmoothingRecursiveYvvGaussianImageFilter::GetSigmaArray() const +{ + return m_Sigma; +} + +// Get the sigma scalar. If the sigma is anisotropic, we will just +// return the sigma along the first dimension. +template +typename GPUSmoothingRecursiveYvvGaussianImageFilter::ScalarRealType +GPUSmoothingRecursiveYvvGaussianImageFilter::GetSigma() const +{ + return m_Sigma[0]; +} + +template +void +GPUSmoothingRecursiveYvvGaussianImageFilter::SetNormalizeAcrossScale(bool normalize) +{ + m_NormalizeAcrossScale = normalize; + this->Modified(); +} + +template +void +GPUSmoothingRecursiveYvvGaussianImageFilter::GenerateInputRequestedRegion() ITK_NOEXCEPT +{ + // call the superclass' implementation of this method. this should + // copy the output requested region to the input requested region + Superclass::GenerateInputRequestedRegion(); + + // This filter needs all of the input + typename GPUSmoothingRecursiveYvvGaussianImageFilter::InputImagePointer image = + const_cast(this->GetInput()); + if (image) + { + image->SetRequestedRegion(this->GetInput()->GetLargestPossibleRegion()); + } +} + +template +void +GPUSmoothingRecursiveYvvGaussianImageFilter::EnlargeOutputRequestedRegion( + DataObject * output) +{ + auto * out = dynamic_cast(output); + + if (out) + { + out->SetRequestedRegion(out->GetLargestPossibleRegion()); + } +} + +template +void +GPUSmoothingRecursiveYvvGaussianImageFilter::AllocateGPUCoefficients() +{ + m_GPUBCoefficientsDataManager = GPUDataManager::New(); + m_GPUBCoefficientsDataManager->SetBufferSize(4 * sizeof(ScalarRealType)); + m_GPUBCoefficientsDataManager->SetCPUBufferPointer(m_Bvalues.data()); + m_GPUBCoefficientsDataManager->SetBufferFlag(CL_MEM_READ_ONLY); + m_GPUBCoefficientsDataManager->Allocate(); + if (!m_Bvalues.empty()) + { + m_GPUBCoefficientsDataManager->SetGPUDirtyFlag(true); + } + + m_GPUMMatrixDataManager = GPUDataManager::New(); + m_GPUMMatrixDataManager->SetBufferSize(9 * sizeof(ScalarRealType)); + m_GPUMMatrixDataManager->SetCPUBufferPointer(m_CPUMatrix.data()); + m_GPUMMatrixDataManager->SetBufferFlag(CL_MEM_READ_ONLY); + m_GPUMMatrixDataManager->Allocate(); + if (!m_CPUMatrix.empty()) + { + m_GPUMMatrixDataManager->SetGPUDirtyFlag(true); + } +} + +/** + * Compute filter for Gaussian kernel + */ +template +void +GPUSmoothingRecursiveYvvGaussianImageFilter::GPUGenerateData() +{ + if (m_FilterGPUKernelHandle == -1) + { + this->BuildKernel(); + } + + constexpr unsigned int X = 0; + constexpr unsigned int Y = 1; + constexpr unsigned int Z = 2; + // arguments set up + int argidx = 0; + const int ndim = (int)TInputImage::ImageDimension; + + this->m_GPUKernelManager->SetKernelArgWithImage(m_FilterGPUKernelHandle, argidx++, inPtr->GetGPUDataManager()); + this->m_GPUKernelManager->SetKernelArgWithImage(m_FilterGPUKernelHandle, argidx++, otPtr->GetGPUDataManager()); + this->m_GPUKernelManager->SetKernelArgWithImage(m_FilterGPUKernelHandle, argidx++, m_GPUBCoefficientsDataManager); + this->m_GPUKernelManager->SetKernelArgWithImage(m_FilterGPUKernelHandle, argidx++, m_GPUMMatrixDataManager); + + for (int i = 0; i < ndim; ++i) + { + this->m_GPUKernelManager->SetKernelArg(m_FilterGPUKernelHandle, argidx++, sizeof(int), &(m_requestedSize[i])); + } + + const unsigned int dimArg = argidx; + + if (ndim < 3) + { + int globalSize1D = (m_requestedSize[1] > m_requestedSize[0] ? m_requestedSize[1] : m_requestedSize[0]); + + this->m_GPUKernelManager->SetKernelArg(m_FilterGPUKernelHandle, dimArg, sizeof(unsigned int), &(X)); + this->m_GPUKernelManager->LaunchKernel1D(m_FilterGPUKernelHandle, globalSize1D, 16); + + // change ONLY input and direction of filter + this->m_GPUKernelManager->SetKernelArgWithImage(m_FilterGPUKernelHandle, 0, otPtr->GetGPUDataManager()); + this->m_GPUKernelManager->SetKernelArg(m_FilterGPUKernelHandle, dimArg, sizeof(unsigned int), &(Y)); + this->m_GPUKernelManager->LaunchKernel1D(m_FilterGPUKernelHandle, globalSize1D, 16); + } + else + { + // We must optimize our 2D workgroup sizes to go over our 3D volume in each + // dimension. + + this->m_GPUKernelManager->SetKernelArg(m_FilterGPUKernelHandle, argidx++, sizeof(unsigned int), &(X)); + this->m_GPUKernelManager->LaunchKernel2D(m_FilterGPUKernelHandle, m_requestedSize[2], m_requestedSize[1], 16, 16); + + // change ONLY input and direction of filter + this->m_GPUKernelManager->SetKernelArgWithImage(m_FilterGPUKernelHandle, 0, otPtr->GetGPUDataManager()); + this->m_GPUKernelManager->SetKernelArg(m_FilterGPUKernelHandle, dimArg, sizeof(unsigned int), &(Y)); + this->m_GPUKernelManager->LaunchKernel2D(m_FilterGPUKernelHandle, m_requestedSize[0], m_requestedSize[2], 16, 16); + + // input is already pointing to previous output; change ONLY direction of + // filter + this->m_GPUKernelManager->SetKernelArg(m_FilterGPUKernelHandle, dimArg, sizeof(unsigned int), &(Z)); + this->m_GPUKernelManager->LaunchKernel2D(m_FilterGPUKernelHandle, m_requestedSize[0], m_requestedSize[1], 16, 16); + } +} + +template +void +GPUSmoothingRecursiveYvvGaussianImageFilter::PrintSelf(std::ostream & os, + Indent indent) const +{ + Superclass::PrintSelf(os, indent); + + os << "NormalizeAcrossScale: " << m_NormalizeAcrossScale << std::endl; + os << "Sigma: " << m_Sigma << std::endl; +} +} // end namespace itk + +#endif // GPU +#endif // itkGPUSmoothingRecursiveYvvGaussianImageFilter_hxx diff --git a/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/include/itkRecursiveLineYvvGaussianImageFilter.h b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/include/itkRecursiveLineYvvGaussianImageFilter.h new file mode 100644 index 00000000000..4298c5ed2a5 --- /dev/null +++ b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/include/itkRecursiveLineYvvGaussianImageFilter.h @@ -0,0 +1,197 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +#ifndef itkRecursiveLineYvvGaussianImageFilter_h +#define itkRecursiveLineYvvGaussianImageFilter_h + +#include "itkInPlaceImageFilter.h" +#include "itkNumericTraits.h" +#include "itkImageRegionSplitterDirection.h" + +namespace itk +{ +/** + * \class RecursiveLineYvvGaussianImageFilter + * \brief 1D recursive Gaussian blur based on Young-Van Vliet's algorithm, + * implemented for CPU. + * + * This CPU implementation is more efficient than the GPU implamentation for + * smaller images (e.g. 512 and smaller for quadcores at over 3GHz); use + * the benchmark tests to establish the size for which this implementation + * performs better for your particular hardware configuration. + * + * More information in the Insight Journal publication: + * https://hdl.handle.net/10380/3425 + * + * \ingroup SmoothingRecursiveYvvGaussianFilter + */ + +template +class ITK_EXPORT RecursiveLineYvvGaussianImageFilter : public InPlaceImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(RecursiveLineYvvGaussianImageFilter); + + /** Standard class type alias. */ + using Self = RecursiveLineYvvGaussianImageFilter; + using Superclass = InPlaceImageFilter; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Type macro that defines a name for this class. */ + itkTypeMacro(RecursiveLineYvvGaussianImageFilter, InPlaceImageFilter); + + /** Smart pointer type alias support. */ + using InputImagePointer = typename TInputImage::Pointer; + using InputImageConstPointer = typename TInputImage::ConstPointer; + + /** Real type to be used in internal computations. RealType in general is + * templated over the pixel type. (For example for vector or tensor pixels, + * RealType is a vector or a tensor of doubles.) ScalarRealType is a type + * meant for scalars. + */ + using InputPixelType = typename TInputImage::PixelType; +#ifdef WITH_DOUBLE + using RealType = typename NumericTraits::RealType; + using ScalarRealType = typename NumericTraits::ScalarRealType; +#else + using RealType = typename NumericTraits::FloatType; + using ScalarRealType = typename NumericTraits::FloatType; +#endif + + using OutputImageRegionType = typename TOutputImage::RegionType; + + /** Type of the input image */ + using InputImageType = TInputImage; + + /** Type of the output image */ + using OutputImageType = TOutputImage; + + /** Get the direction in which the filter is to be applied. */ + itkGetConstMacro(Direction, unsigned int); + + /** Set the direction in which the filter is to be applied. */ + itkSetMacro(Direction, unsigned int); + + /** Set Input Image. */ + void + SetInputImage(const TInputImage *); + + /** Get Input Image. */ + const TInputImage * + GetInputImage(); + + /** Set/Get the flag for normalizing the gaussian over scale space. + When this flag is ON the filter will be normalized in such a way + that larger sigmas will not result in the image fading away. + + \f[ + \frac{ 1 }{ \sqrt{ 2 \pi } }; + \f] + + When the flag is OFF the normalization will conserve contant the + integral of the image intensity. + \f[ + \frac{ 1 }{ \sigma \sqrt{ 2 \pi } }; + \f] + For analyzing an image across Scale Space you want to enable + this flag. It is disabled by default. */ + itkSetMacro(NormalizeAcrossScale, bool); + itkGetConstMacro(NormalizeAcrossScale, bool); + + /** Set/Get the Sigma, measured in world coordinates, of the Gaussian + * kernel. The default is 1.0. */ + itkGetConstMacro(Sigma, ScalarRealType); + itkSetMacro(Sigma, ScalarRealType); + +protected: + RecursiveLineYvvGaussianImageFilter(); + ~RecursiveLineYvvGaussianImageFilter() override = default; + void + PrintSelf(std::ostream & os, Indent indent) const override; + + /** GenerateData (apply) the filter. */ + void + BeforeThreadedGenerateData() override; + + void + ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, ThreadIdType threadId) override; + + const ImageRegionSplitterBase * + GetImageRegionSplitter() const override; + + /** RecursiveLineYvvGaussianImageFilter needs all of the input only in the + * "Direction" dimension. Therefore we enlarge the output's + * RequestedRegion to this. Then the superclass's + * GenerateInputRequestedRegion method will copy the output region + * to the input. + * + * \sa ImageToImageFilter::GenerateInputRequestedRegion() + */ + void + EnlargeOutputRequestedRegion(DataObject * output) override; + + /** Set up the coefficients of the filter to approximate a specific kernel. + * Typically it can be used to approximate a Gaussian or one of its + * derivatives. Parameter is the spacing along the dimension to + * filter. */ + virtual void + SetUp(ScalarRealType spacing); + + /** Apply the Recursive Filter to an array of data. This method is called + * for each line of the volume. Parameter "scratch" is a scratch + * area used for internal computations that is the same size as the + * parameters "outs" and "data". The scratch area must be allocated + * outside of this routine (this avoids memory allocation and + * deallocation in the inner loop of the overall algorithm. */ + void + FilterDataArray(RealType * outs, const RealType * data, RealType * scratch, unsigned int ln); + +protected: + /** Causal and anti-causal coefficients that multiply the input data. These + are already divided by B0 */ + ScalarRealType m_B1; + ScalarRealType m_B2; + ScalarRealType m_B3; + ScalarRealType m_B; + + // Initialization matrix for anti-causal pass + vnl_matrix m_MMatrix; + +private: + /** Direction in which the filter is to be applied + * this should be in the range [0,ImageDimension-1]. */ + unsigned int m_Direction; + + /** Sigma of the gaussian kernel. */ + ScalarRealType m_Sigma; + + /** Normalize the image across scale space */ + bool m_NormalizeAcrossScale; + ImageRegionSplitterDirection::Pointer m_ImageRegionSplitter; +}; +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkRecursiveLineYvvGaussianImageFilter.hxx" +#endif + +#endif // itkRecursiveLineYvvGaussianImageFilter_h diff --git a/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/include/itkRecursiveLineYvvGaussianImageFilter.hxx b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/include/itkRecursiveLineYvvGaussianImageFilter.hxx new file mode 100644 index 00000000000..e535166cbcd --- /dev/null +++ b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/include/itkRecursiveLineYvvGaussianImageFilter.hxx @@ -0,0 +1,441 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +#ifndef itkRecursiveLineYvvGaussianImageFilter_hxx +#define itkRecursiveLineYvvGaussianImageFilter_hxx + +#include "itkObjectFactory.h" +#include "itkImageLinearIteratorWithIndex.h" +#include "itkImageLinearConstIteratorWithIndex.h" +#include "itkProgressReporter.h" + +namespace itk +{ +template +RecursiveLineYvvGaussianImageFilter::RecursiveLineYvvGaussianImageFilter() +{ + m_Direction = 0; + this->SetNumberOfRequiredOutputs(1); + this->SetNumberOfRequiredInputs(1); + + this->InPlaceOff(); + this->DynamicMultiThreadingOff(); + + m_ImageRegionSplitter = ImageRegionSplitterDirection::New(); + + if (this->GetDebug()) + { + std::cout << "-----------Line filter TYPES\n"; + + if (typeid(typename TInputImage::PixelType) == typeid(double)) + { + std::cout << "InputPixelType double\n"; + } + if (typeid(typename TOutputImage::PixelType) == typeid(double)) + { + std::cout << "OutputPixelType double\n"; + } + + if (typeid(ScalarRealType) == typeid(double)) + { + std::cout << "ScalarRealType double\n"; + } + + if (typeid(RealType) == typeid(double)) + { + std::cout << "RealType double\n"; + } + } +} + +/** + * Set Input Image + */ +template +void +RecursiveLineYvvGaussianImageFilter::SetInputImage(const TInputImage * input) +{ + // ProcessObject is not const_correct so this const_cast is required + ProcessObject::SetNthInput(0, const_cast(input)); +} + + +template +const TInputImage * +RecursiveLineYvvGaussianImageFilter::GetInputImage() +{ + return dynamic_cast((ProcessObject::GetInput(0))); +} + +template +void +RecursiveLineYvvGaussianImageFilter::SetUp(ScalarRealType spacing) +{ + const ScalarRealType sigmad = m_Sigma / spacing; + + // Compute q according to 16 in Young et al on Gabor filering + ScalarRealType q = 0; + + if (sigmad >= 3.556) + { + q = 0.9804 * (sigmad - 3.556) + 2.5091; + } + else + { + if (sigmad < 0.5) + { + std::cerr << "Too low sigma value (< 0.5), computation will not be precise." << std::endl; + } + + q = 0.0561 * sigmad * sigmad + 0.5784 * sigmad - 0.2568; + } + + // Compute B and B1 to B3 according to Young et al 2003 + ScalarRealType m0 = 1.16680; + ScalarRealType m1 = 1.10783; + ScalarRealType m2 = 1.40586; + ScalarRealType scale = (m0 + q) * (m1 * m1 + m2 * m2 + 2 * m1 * q + q * q); + + m_B1 = q * (2 * m0 * m1 + m1 * m1 + m2 * m2 + (2 * m0 + 4 * m1) * q + 3 * q * q) / scale; + m_B2 = -q * q * (m0 + 2 * m1 + 3 * q) / scale; + m_B3 = q * q * q / scale; + + ScalarRealType baseB = (m0 * (m1 * m1 + m2 * m2)) / scale; + m_B = baseB * baseB; + + // M Matrix for initialization on backward pass, from Triggs and Sdika, IEEE + // TSP + m_MMatrix = vnl_matrix(3, 3); + + m_MMatrix(0, 0) = -m_B3 * m_B1 + 1 - m_B3 * m_B3 - m_B2; + m_MMatrix(0, 1) = (m_B3 + m_B1) * (m_B2 + m_B3 * m_B1); + m_MMatrix(0, 2) = m_B3 * (m_B1 + m_B3 * m_B2); + + m_MMatrix(1, 0) = m_B1 + m_B3 * m_B2; + m_MMatrix(1, 1) = (1 - m_B2) * (m_B2 + m_B3 * m_B1); + m_MMatrix(1, 2) = -m_B3 * (m_B3 * m_B1 + m_B3 * m_B3 + m_B2 - 1); + + m_MMatrix(2, 0) = m_B3 * m_B1 + m_B2 + m_B1 * m_B1 - m_B2 * m_B2; + m_MMatrix(2, 1) = m_B1 * m_B2 + m_B3 * m_B2 * m_B2 - m_B1 * m_B3 * m_B3 - m_B3 * m_B3 * m_B3 - m_B3 * m_B2 + m_B3; + m_MMatrix(2, 2) = m_B3 * (m_B1 + m_B3 * m_B2); + + m_MMatrix /= (1 + m_B1 - m_B2 + m_B3) * (1 - m_B1 - m_B2 - m_B3) * (1 + m_B2 + (m_B1 - m_B3) * m_B3); + + if (this->GetDebug()) + { + std::cout << "cB " << m_B << std::endl; + std::cout << "cB1 " << m_B1 << std::endl; + std::cout << "cB2 " << m_B2 << std::endl; + std::cout << "cB3 " << m_B3 << std::endl; + + for (int i = 0; i < 3; ++i) + { + for (int j = 0; j < 3; ++j) + { + std::cout << "cM(" << i << "," << j << ") " << m_MMatrix(i, j) << std::endl; + } + } + } +} + + +// Apply Recursive Filter +template +void +RecursiveLineYvvGaussianImageFilter::FilterDataArray(RealType * outs, + const RealType * data, + RealType * scratch, + unsigned int ln) +{ + /** + * Causal direction pass + */ + + // this value is assumed to exist from the border to infinity. + const RealType outV1 = data[0] / (1.0 - m_B1 - m_B2 - m_B3); + RealType sV0 = outV1, sV1 = outV1, sV2 = outV1; + + /** + * Recursively filter the rest + */ + + for (unsigned int i = 0; i < ln; i++) + { + scratch[i] = RealType(data[i] + sV0 * m_B1 + sV1 * m_B2 + sV2 * m_B3); + sV2 = sV1; + sV1 = sV0; + sV0 = scratch[i]; + } + /** + * Store the causal result + */ + + for (unsigned int i = 0; i < ln; ++i) + { + outs[i] = scratch[i]; + } + + // AntiCausal direction pass + + // Handle outside values according to Triggs and Sdika + const RealType u_p = data[ln - 1] / (1.0 - m_B1 - m_B2 - m_B3); + const RealType v_p = u_p / (1.0 - m_B1 - m_B2 - m_B3); + + RealType Vn0 = v_p; + RealType Vn1 = v_p; + RealType Vn2 = v_p; + for (unsigned int i = 0; i < 3; ++i) + { + Vn0 += (outs[ln - 1 - i] - u_p) * m_MMatrix(0, i); + Vn1 += (outs[ln - 1 - i] - u_p) * m_MMatrix(1, i); + Vn2 += (outs[ln - 1 - i] - u_p) * m_MMatrix(2, i); + } + + // This was not in the 2006 Triggs paper but sounds quite logical since m_B is + // not one + Vn0 *= m_B; + Vn1 *= m_B; + Vn2 *= m_B; + + scratch[ln - 1] = Vn0; + + // Recursively filter the rest + + for (int i = ln - 2; i >= 0; i--) + { + scratch[i] = RealType(outs[i] * m_B + Vn0 * m_B1 + Vn1 * m_B2 + Vn2 * m_B3); + Vn2 = Vn1; + Vn1 = Vn0; + Vn0 = scratch[i]; + } + + /** + * Roll the antiCausal part into the output + */ + for (unsigned int i = 0; i < ln; ++i) + { + outs[i] = scratch[i]; + } +} + +// we need all of the image in just the "Direction" we are separated into +template +void +RecursiveLineYvvGaussianImageFilter::EnlargeOutputRequestedRegion(DataObject * output) +{ + auto * out = dynamic_cast(output); + + if (out) + { + OutputImageRegionType outputRegion = out->GetRequestedRegion(); + const OutputImageRegionType & largestOutputRegion = out->GetLargestPossibleRegion(); + + // verify sane parameter + if (this->m_Direction >= outputRegion.GetImageDimension()) + { + itkExceptionMacro("Direction selected for filtering is greater than ImageDimension"); + } + + // expand output region to match largest in the "Direction" dimension + outputRegion.SetIndex(m_Direction, largestOutputRegion.GetIndex(m_Direction)); + outputRegion.SetSize(m_Direction, largestOutputRegion.GetSize(m_Direction)); + + out->SetRequestedRegion(outputRegion); + } +} + +template +const ImageRegionSplitterBase * +RecursiveLineYvvGaussianImageFilter::GetImageRegionSplitter() const +{ + return this->m_ImageRegionSplitter; +} + +template +void +RecursiveLineYvvGaussianImageFilter::BeforeThreadedGenerateData() +{ + using RegionType = ImageRegion; + + typename TInputImage::ConstPointer inputImage(this->GetInputImage()); + typename TOutputImage::Pointer outputImage(this->GetOutput()); + + const unsigned int imageDimension = inputImage->GetImageDimension(); + + if (this->m_Direction >= imageDimension) + { + itkExceptionMacro("Direction selected for filtering is greater than ImageDimension"); + } + + const typename InputImageType::SpacingType & pixelSize = inputImage->GetSpacing(); + + this->m_ImageRegionSplitter->SetDirection(m_Direction); + this->SetUp(pixelSize[m_Direction]); + + RegionType region = outputImage->GetRequestedRegion(); + + const unsigned int ln = region.GetSize()[this->m_Direction]; + + if (ln < 4) + { + itkExceptionMacro( + "The number of pixels along direction " + << this->m_Direction + << " is less than 4. This filter requires a minimum of four pixels along the dimension to be processed."); + } +} + +/** + * Compute Recursive filter + * line by line in one of the dimensions + */ +template +void +RecursiveLineYvvGaussianImageFilter::ThreadedGenerateData( + const OutputImageRegionType & outputRegionForThread, + ThreadIdType threadId) +{ + using OutputPixelType = typename TOutputImage::PixelType; + + using InputConstIteratorType = ImageLinearConstIteratorWithIndex; + using OutputIteratorType = ImageLinearIteratorWithIndex; + + using RegionType = ImageRegion; + + typename TInputImage::ConstPointer inputImage(this->GetInputImage()); + typename TOutputImage::Pointer outputImage(this->GetOutput()); + + RegionType region = outputRegionForThread; + + InputConstIteratorType inputIterator(inputImage, region); + OutputIteratorType outputIterator(outputImage, region); + + inputIterator.SetDirection(this->m_Direction); + outputIterator.SetDirection(this->m_Direction); + + const unsigned int ln = region.GetSize()[this->m_Direction]; + + if (ln < 4) + { + itkExceptionMacro("The number of pixels along direction " + << this->m_Direction << " in this work-unit image region is less than 4." + << " Note: TBBMultiThreader can divide small regions into really small pieces."); + } + + RealType * inps = nullptr; + RealType * outs = nullptr; + RealType * scratch = nullptr; + + try + { + inps = new RealType[ln]; + } + catch (std::bad_alloc &) + { + itkExceptionMacro("Problem allocating memory for internal computations"); + } + + try + { + outs = new RealType[ln]; + } + catch (std::bad_alloc &) + { + delete[] inps; + itkExceptionMacro("Problem allocating memory for internal computations"); + } + + try + { + scratch = new RealType[ln]; + } + catch (std::bad_alloc &) + { + delete[] inps; + delete[] outs; + itkExceptionMacro("Problem allocating memory for internal computations"); + } + + inputIterator.GoToBegin(); + outputIterator.GoToBegin(); + + const unsigned int numberOfLinesToProcess = + outputRegionForThread.GetNumberOfPixels() / outputRegionForThread.GetSize(this->m_Direction); + ProgressReporter progress(this, threadId, numberOfLinesToProcess, 10); + + try // this try is intended to catch an eventual AbortException. + { + while (!inputIterator.IsAtEnd() && !outputIterator.IsAtEnd()) + { + unsigned int i = 0; + while (!inputIterator.IsAtEndOfLine()) + { + inps[i++] = inputIterator.Get(); + ++inputIterator; + } + + this->FilterDataArray(outs, inps, scratch, ln); + + unsigned int j = 0; + while (!outputIterator.IsAtEndOfLine()) + { + outputIterator.Set(static_cast(outs[j++])); + ++outputIterator; + } + + inputIterator.NextLine(); + outputIterator.NextLine(); + + // Although the method name is CompletedPixel(), + // this is being called after each line is processed + progress.CompletedPixel(); + } + } + catch (ProcessAborted &) + { + // User aborted filter excecution Here we catch an exception thrown by the + // progress reporter and rethrow it with the correct line number and file + // name. We also invoke AbortEvent in case some observer was interested on + // it. + // release locally allocated memory + delete[] outs; + delete[] inps; + delete[] scratch; + // Throw the final exception. + ProcessAborted e(__FILE__, __LINE__); + e.SetDescription("Process aborted."); + e.SetLocation(ITK_LOCATION); + throw e; + } + + delete[] outs; + delete[] inps; + delete[] scratch; +} + +template +void +RecursiveLineYvvGaussianImageFilter::PrintSelf(std::ostream & os, Indent indent) const +{ + Superclass::PrintSelf(os, indent); + + os << indent << "Direction: " << m_Direction << std::endl; +} +} // end namespace itk + +#endif diff --git a/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/include/itkSmoothingRecursiveYvvGaussianImageFilter.h b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/include/itkSmoothingRecursiveYvvGaussianImageFilter.h new file mode 100644 index 00000000000..12ad0b66cd0 --- /dev/null +++ b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/include/itkSmoothingRecursiveYvvGaussianImageFilter.h @@ -0,0 +1,178 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +#ifndef itkSmoothingRecursiveYvvGaussianImageFilter_h +#define itkSmoothingRecursiveYvvGaussianImageFilter_h + +#include "itkRecursiveLineYvvGaussianImageFilter.h" +#include "itkCastImageFilter.h" +#include "itkImage.h" +#include "itkPixelTraits.h" +#include "itkCommand.h" +#include "itkFixedArray.h" + +namespace itk +{ +/** + * \class SmoothingRecursiveYvvGaussianImageFilter + * \brief Recursive Gaussian blurring filter based on Young-Van Vliet's + * algorithm, implemented for CPU. + * + * This CPU implementation is more efficient than the GPU implamentation for + * smaller images (e.g. 512x512 and smaller for quadcores at over 3GHz); use + * the benchmark tests to establish the size for which this implementation + * performs better for your particular hardware configuration. + * + * More information in the Insight Journal publication: + * https://hdl.handle.net/10380/3425 + * + * \ingroup SmoothingRecursiveYvvGaussianFilter + */ +template +class ITK_EXPORT SmoothingRecursiveYvvGaussianImageFilter : public InPlaceImageFilter +{ +public: + ITK_DISALLOW_COPY_AND_MOVE(SmoothingRecursiveYvvGaussianImageFilter); + + /** Standard class type alias. */ + using Self = SmoothingRecursiveYvvGaussianImageFilter; + using Superclass = InPlaceImageFilter; + using Pointer = SmartPointer; + using ConstPointer = SmartPointer; + + /** Pixel Type of the input image */ + using InputImageType = TInputImage; + using OutputImageType = TOutputImage; + using PixelType = typename TInputImage::PixelType; +#ifdef WITH_DOUBLE + using RealType = typename NumericTraits::RealType; + using ScalarRealType = typename NumericTraits::ScalarRealType; +#else + using RealType = typename NumericTraits::FloatType; + using ScalarRealType = typename NumericTraits::FloatType; +#endif + + /** Runtime information support. */ + itkTypeMacro(SmoothingRecursiveYvvGaussianImageFilter, InPlaceImageFilter); + + /** Image dimension. */ + static constexpr unsigned int ImageDimension = TInputImage::ImageDimension; + + /** Define the type for the sigma array */ + using SigmaArrayType = FixedArray; + + /** Define the image type for internal computations + RealType is usually 'double' in NumericTraits. + Here we prefer float in order to save memory. */ + + using InternalRealType = typename NumericTraits::FloatType; + using RealImageType = typename InputImageType::template Rebind::Type; + + /** The first in the pipeline */ + typedef RecursiveLineYvvGaussianImageFilter FirstGaussianFilterType; + + /** Smoothing filter type */ + typedef RecursiveLineYvvGaussianImageFilter InternalGaussianFilterType; + + /** The last in the pipeline */ + typedef CastImageFilter CastingFilterType; + + /** Pointer to a gaussian filter. */ + using InternalGaussianFilterPointer = typename InternalGaussianFilterType::Pointer; + + /** Pointer to the first gaussian filter. */ + using FirstGaussianFilterPointer = typename FirstGaussianFilterType::Pointer; + + /** Pointer to the last filter, casting */ + using CastingFilterPointer = typename CastingFilterType::Pointer; + + /** Pointer to the Output Image */ + using OutputImagePointer = typename OutputImageType::Pointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Set Sigma value. Sigma is measured in the units of image spacing. You + may use the method SetSigma to set the same value across each axis or + use the method SetSigmaArray if you need different values along each + axis. */ + void + SetSigmaArray(const SigmaArrayType & sigmas); + void + SetSigma(ScalarRealType sigma); + SigmaArrayType + GetSigmaArray() const; + ScalarRealType + GetSigma() const; + + /** Define which normalization factor will be used for the Gaussian */ + void + SetNormalizeAcrossScale(bool normalizeInScaleSpace); + itkGetConstMacro(NormalizeAcrossScale, bool); + + void + SetNumberOfWorkUnits(ThreadIdType nb) override; + + bool + CanRunInPlace() const override; + +#ifdef ITK_USE_CONCEPT_CHECKING + /** Begin concept checking */ + itkConceptMacro(InputHasNumericTraitsCheck, (Concept::HasNumericTraits)); + /** End concept checking */ +#endif +protected: + SmoothingRecursiveYvvGaussianImageFilter(); + ~SmoothingRecursiveYvvGaussianImageFilter() override = default; + void + PrintSelf(std::ostream & os, Indent indent) const override; + + /** Generate Data */ + void + GenerateData() override; + + /** SmoothingRecursiveYvvGaussianImageFilter needs all of the input to produce an + * output. Therefore, SmoothingRecursiveYvvGaussianImageFilter needs to provide + * an implementation for GenerateInputRequestedRegion in order to inform + * the pipeline execution model. + * \sa ImageToImageFilter::GenerateInputRequestedRegion() */ + void + GenerateInputRequestedRegion() ITK_NOEXCEPT override; + + // Override since the filter produces the entire dataset + void + EnlargeOutputRequestedRegion(DataObject * output) override; + +private: + InternalGaussianFilterPointer m_SmoothingFilters[ImageDimension - 1]; + FirstGaussianFilterPointer m_FirstSmoothingFilter; + CastingFilterPointer m_CastingFilter; + + /** Normalize the image across scale space */ + bool m_NormalizeAcrossScale; + + /** Standard deviation of the gaussian used for smoothing */ + SigmaArrayType m_Sigma; +}; +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +# include "itkSmoothingRecursiveYvvGaussianImageFilter.hxx" +#endif + +#endif diff --git a/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/include/itkSmoothingRecursiveYvvGaussianImageFilter.hxx b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/include/itkSmoothingRecursiveYvvGaussianImageFilter.hxx new file mode 100644 index 00000000000..261674fd98e --- /dev/null +++ b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/include/itkSmoothingRecursiveYvvGaussianImageFilter.hxx @@ -0,0 +1,291 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +#ifndef itkSmoothingRecursiveYvvGaussianImageFilter_hxx +#define itkSmoothingRecursiveYvvGaussianImageFilter_hxx + +#include "itkImageRegionIteratorWithIndex.h" +#include "itkProgressAccumulator.h" + +namespace itk +{ +template +SmoothingRecursiveYvvGaussianImageFilter::SmoothingRecursiveYvvGaussianImageFilter() +{ + m_NormalizeAcrossScale = false; + + m_FirstSmoothingFilter = FirstGaussianFilterType::New(); + m_FirstSmoothingFilter->SetDirection(ImageDimension - 1); + m_FirstSmoothingFilter->SetNormalizeAcrossScale(m_NormalizeAcrossScale); + m_FirstSmoothingFilter->ReleaseDataFlagOn(); + + for (unsigned int i = 0; i < ImageDimension - 1; i++) + { + m_SmoothingFilters[i] = InternalGaussianFilterType::New(); + m_SmoothingFilters[i]->SetNormalizeAcrossScale(m_NormalizeAcrossScale); + m_SmoothingFilters[i]->SetDirection(i); + m_SmoothingFilters[i]->ReleaseDataFlagOn(); + m_SmoothingFilters[i]->InPlaceOn(); + } + + m_SmoothingFilters[0]->SetInput(m_FirstSmoothingFilter->GetOutput()); + for (unsigned int i = 1; i < ImageDimension - 1; i++) + { + m_SmoothingFilters[i]->SetInput(m_SmoothingFilters[i - 1]->GetOutput()); + } + + m_CastingFilter = CastingFilterType::New(); + m_CastingFilter->SetInput(m_SmoothingFilters[ImageDimension - 2]->GetOutput()); + m_CastingFilter->InPlaceOn(); + + this->InPlaceOff(); + + // + // NB: We must call SetSigma in order to initialize the smoothing + // filters with the default scale. However, m_Sigma must first be + // initialized (it is used inside SetSigma) and it must be different + // from 1.0 or the call will be ignored. + this->m_Sigma.Fill(0.0); + this->SetSigma(1.0); + + if (this->GetDebug()) + { + std::cout << "-----------Smoothing filter TYPES\n"; + + if (typeid(typename TInputImage::PixelType) == typeid(double)) + { + std::cout << "PixelType double\n"; + } + if (typeid(typename TOutputImage::PixelType) == typeid(double)) + { + std::cout << "Output PixelType double\n"; + } + + if (typeid(ScalarRealType) == typeid(double)) + { + std::cout << "ScalarRealType double\n"; + } + + if (typeid(RealType) == typeid(double)) + { + std::cout << "RealType double\n"; + } + + if (typeid(InternalRealType) == typeid(double)) + { + std::cout << "InternalRealType double\n"; + } + } +} + +template +void +SmoothingRecursiveYvvGaussianImageFilter::SetNumberOfWorkUnits(ThreadIdType nb) +{ + Superclass::SetNumberOfWorkUnits(nb); + + for (unsigned int i = 0; i < ImageDimension - 1; i++) + { + m_SmoothingFilters[i]->SetNumberOfWorkUnits(nb); + } + m_FirstSmoothingFilter->SetNumberOfWorkUnits(nb); +} + +template +bool +SmoothingRecursiveYvvGaussianImageFilter::CanRunInPlace() const +{ + // Note: There are two different ways this filter may try to run + // in-place: + // 1) Similar to the standard way, when the input and output image + // are of the same type, they can share the bulk data. The output + // will be grafted onto the last filter. In this fashion the input + // and output will be the same bulk data, but the intermediate + // mini-pipeline will use different data. + // 2) If the input image is the same type as the RealImage used for + // the mini-pipeline, then all the filters may re-use the same + // bulk data, stealing it from the input then moving it down the + // pipeline filter by filter. Additionally, if the output is also + // RealType then the last filter will run in-place making the entire + // pipeline in-place and only utilizing on copy of the bulk data. + + return m_FirstSmoothingFilter->CanRunInPlace() || this->Superclass::CanRunInPlace(); +} + +template +void +SmoothingRecursiveYvvGaussianImageFilter::SetSigma(ScalarRealType sigma) +{ + SigmaArrayType sigmas(sigma); + this->SetSigmaArray(sigmas); +} + +template +void +SmoothingRecursiveYvvGaussianImageFilter::SetSigmaArray(const SigmaArrayType & sigma) +{ + if (this->m_Sigma != sigma) + { + this->m_Sigma = sigma; + for (unsigned int i = 0; i < ImageDimension - 1; i++) + { + m_SmoothingFilters[i]->SetSigma(m_Sigma[i]); + m_SmoothingFilters[i]->Modified(); + } + m_FirstSmoothingFilter->SetSigma(m_Sigma[ImageDimension - 1]); + m_FirstSmoothingFilter->Modified(); + + this->Modified(); + } +} + +template +typename SmoothingRecursiveYvvGaussianImageFilter::SigmaArrayType +SmoothingRecursiveYvvGaussianImageFilter::GetSigmaArray() const +{ + return m_Sigma; +} + +template +typename SmoothingRecursiveYvvGaussianImageFilter::ScalarRealType +SmoothingRecursiveYvvGaussianImageFilter::GetSigma() const +{ + return m_Sigma[0]; +} + +template +void +SmoothingRecursiveYvvGaussianImageFilter::SetNormalizeAcrossScale(bool normalize) +{ + m_NormalizeAcrossScale = normalize; + + for (unsigned int i = 0; i < ImageDimension - 1; i++) + { + m_SmoothingFilters[i]->SetNormalizeAcrossScale(normalize); + } + m_FirstSmoothingFilter->SetNormalizeAcrossScale(normalize); + + this->Modified(); +} + +template +void +SmoothingRecursiveYvvGaussianImageFilter::GenerateInputRequestedRegion() ITK_NOEXCEPT +{ + // call the superclass' implementation of this method. this should + // copy the output requested region to the input requested region + Superclass::GenerateInputRequestedRegion(); + + // This filter needs all of the input + typename SmoothingRecursiveYvvGaussianImageFilter::InputImagePointer image = + const_cast(this->GetInput()); + if (image) + { + image->SetRequestedRegion(this->GetInput()->GetLargestPossibleRegion()); + } +} + +template +void +SmoothingRecursiveYvvGaussianImageFilter::EnlargeOutputRequestedRegion(DataObject * output) +{ + auto * out = dynamic_cast(output); + + if (out) + { + out->SetRequestedRegion(out->GetLargestPossibleRegion()); + } +} + +// Compute filter for Gaussian kernel +template +void +SmoothingRecursiveYvvGaussianImageFilter::GenerateData() +{ + itkDebugMacro(<< "SmoothingRecursiveYvvGaussianImageFilter generating data "); + + const typename TInputImage::ConstPointer inputImage(this->GetInput()); + + const typename TInputImage::RegionType region = inputImage->GetRequestedRegion(); + const typename TInputImage::SizeType size = region.GetSize(); + + for (unsigned int d = 0; d < ImageDimension; d++) + { + if (size[d] < 4) + { + itkExceptionMacro( + "The number of pixels along dimension " + << d << " is less than 4. This filter requires a minimum of four pixels along the dimension to be processed."); + } + } + + // If this filter is running in-place, then set the first smoothing + // filter to steal the bulk data, by running in-place. + if (this->CanRunInPlace() && this->GetInPlace()) + { + m_FirstSmoothingFilter->InPlaceOn(); + + // To make this filter's input and out share the same data, call + // the InPlace's/Superclass's allocate methods, which takes care + // of the needed bulk data sharing. + this->AllocateOutputs(); + } + else + { + m_FirstSmoothingFilter->InPlaceOff(); + } + + // If the last filter is running in-place then this bulk data is not + // needed, release it to save memory. + if (m_CastingFilter->CanRunInPlace()) + { + this->GetOutput()->ReleaseData(); + } + + // Create a process accumulator for tracking the progress of this minipipeline + ProgressAccumulator::Pointer progress = ProgressAccumulator::New(); + progress->SetMiniPipelineFilter(this); + + // Register the filter with the with progress accumulator using + // equal weight proportion + for (unsigned int i = 0; i < ImageDimension - 1; ++i) + { + progress->RegisterInternalFilter(m_SmoothingFilters[i], 1.0 / (ImageDimension)); + } + + progress->RegisterInternalFilter(m_FirstSmoothingFilter, 1.0 / (ImageDimension)); + m_FirstSmoothingFilter->SetInput(inputImage); + // graft our output to the internal filter to force the proper regions + // to be generated + m_CastingFilter->GraftOutput(this->GetOutput()); + m_CastingFilter->Update(); + this->GraftOutput(m_CastingFilter->GetOutput()); +} + +template +void +SmoothingRecursiveYvvGaussianImageFilter::PrintSelf(std::ostream & os, Indent indent) const +{ + Superclass::PrintSelf(os, indent); + + os << "NormalizeAcrossScale: " << m_NormalizeAcrossScale << std::endl; + os << "Sigma: " << m_Sigma << std::endl; +} +} // end namespace itk + +#endif diff --git a/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/itk-module.cmake b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/itk-module.cmake new file mode 100644 index 00000000000..52be8dbfcbb --- /dev/null +++ b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/itk-module.cmake @@ -0,0 +1,32 @@ +# the top-level README is used for describing this module, just +# re-used it for documentation here +# itk_module() defines the module dependencies in +# SmoothingRecursiveYvvGaussianFilter +# The testing module in SmoothingRecursiveYvvGaussianFilter depends on +# ITKTestKernel +# By convention those modules outside of ITK are not prefixed with +# ITK. + +# define the dependencies of the include module and the tests +set(ModuleName "SmoothingRecursiveYvvGaussianFilter") +if(ITK_USE_GPU) + set(${ModuleName}_GPU_DEPENDENCIES "ITKGPUSmoothing") + set(${ModuleName}_GPU_COMMON "ITKGPUCommon") +endif() +itk_module( + ${ModuleName} + DEPENDS + ITKCommon + ITKIOImageBase + ITKImageFilterBase + ITKSmoothing + ${${ModuleName}_GPU_DEPENDENCIES} + ${${ModuleName}_GPU_COMMON} + TEST_DEPENDS + ITKTestKernel #to handle IO in src + ${${ModuleName}_GPU_COMMON} + ITKSmoothing + DESCRIPTION "Module ingested from upstream." + EXCLUDE_FROM_DEFAULT + # Header only library does not use ENABLE_SHARED +) diff --git a/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/src/CMakeLists.txt b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/src/CMakeLists.txt new file mode 100644 index 00000000000..356ebde7aee --- /dev/null +++ b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/src/CMakeLists.txt @@ -0,0 +1,11 @@ +if(ITK_USE_GPU) + set(GPU_SRC) + set(GPU_Kernels GPUSmoothingRecursiveYvvGaussianImageFilter.cl) + + #essentially a #define GPU. + add_definitions(-DGPU) + + write_gpu_kernels("${GPU_Kernels}" GPU_SRC) + + itk_module_add_library(${itk-module} ${GPU_SRC}) +endif() diff --git a/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/src/GPUSmoothingRecursiveYvvGaussianImageFilter.cl b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/src/GPUSmoothingRecursiveYvvGaussianImageFilter.cl new file mode 100644 index 00000000000..12ffa074aba --- /dev/null +++ b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/src/GPUSmoothingRecursiveYvvGaussianImageFilter.cl @@ -0,0 +1,235 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +#define BLOCK_DIM 16 // set by ITK for 2D + +// Basic filtering brick, used both for +// 1D filtering of 2D images +// and for 2D filtering of 3D volumes + +__kernel void +partialFilter(const __global INPIXELTYPE * data, + __global OUTPIXELTYPE * outs, + __constant REALTYPE * Bvalues, + __constant REALTYPE * MMatrix, + unsigned int start, + unsigned int step, + unsigned int length) +{ + REALTYPE scratch[MAX_LINE_LENGTH]; + const REALTYPE outV1 = data[start] / (1.0 - Bvalues[1] - Bvalues[2] - Bvalues[3]); + REALTYPE sV0 = outV1, sV1 = outV1, sV2 = outV1; + + // Unrolling does not work well here (truly recursive). + for (unsigned int i = 0; i < length; i++) + { + scratch[i] = (REALTYPE)(data[start + i * step] + sV0 * Bvalues[1] + sV1 * Bvalues[2] + sV2 * Bvalues[3]); + + sV2 = sV1; + sV1 = sV0; + sV0 = scratch[i]; + } + + const REALTYPE remB = (1.0 - Bvalues[1] - Bvalues[2] - Bvalues[3]); + const REALTYPE u_p = data[start + (length - 1) * step] / remB; + const REALTYPE v_p = u_p / remB; + + // Unrolling the loops explicitly helps on certain compilers (e.g. ATI/AMD architecture) and does not + // worsen performance on NVidia GPUs, for which the compiler has been reported to unroll the loops by default). + REALTYPE Vn0 = (REALTYPE)(v_p //; + + (scratch[length - 1] - u_p) * MMatrix[0] + (scratch[length - 2] - u_p) * MMatrix[1] + + (scratch[length - 3] - u_p) * MMatrix[2]); + REALTYPE Vn1 = (REALTYPE)(v_p //; + + (scratch[length - 1] - u_p) * MMatrix[3] + (scratch[length - 2] - u_p) * MMatrix[4] + + (scratch[length - 3] - u_p) * MMatrix[5]); + + REALTYPE Vn2 = (REALTYPE)(v_p //; + + (scratch[length - 1] - u_p) * MMatrix[6] + (scratch[length - 2] - u_p) * MMatrix[7] + + (scratch[length - 3] - u_p) * MMatrix[8]); + + /*REALTYPE Vn0 = v_p; + REALTYPE Vn1 = v_p; + REALTYPE Vn2 = v_p; + + //#pragma unroll 3 + for (unsigned int i = 0;i < 3;++i) + { + Vn0 += (scratch[length - 1 - i] - u_p) * MMatrix[i]; + Vn1 += (scratch[length - 1 - i] - u_p) * MMatrix[3+i]; + Vn2 += (scratch[length - 1 - i] - u_p) * MMatrix[6+i]; + }*/ + + // This was not in the 2006 Triggs paper but sounds quite logical since m_B is not one + Vn0 *= Bvalues[0]; + Vn1 *= Bvalues[0]; + Vn2 *= Bvalues[0]; + + scratch[length - 1] = Vn0; + + // Very bad idea to unroll here as well (see causal pass). + for (int i = length - 2; i >= 0; i--) + { + scratch[i] = (REALTYPE)(scratch[i] * Bvalues[0] + Vn0 * Bvalues[1] + Vn1 * Bvalues[2] + Vn2 * Bvalues[3]); + Vn2 = Vn1; + Vn1 = Vn0; + Vn0 = scratch[i]; + } + +// Roll the anticausal part into the output +#pragma unroll 16 + for (unsigned int i = 0; i < length; ++i) + { + outs[start + i * step] = (OUTPIXELTYPE)(scratch[i]); + } +} + +/* +#ifdef DIM_1 +//This means (recursively) filtering values in one single vector. +//Code present only for completeness (1D->3D). +__kernel void YvvFilter(__global const INPIXELTYPE* data, + __global OUTPIXELTYPE* outs, + __constant REALTYPE* m_Bvalues, + __constant REALTYPE* m_MMatrix, + int sizeX ) +{ + if (sizeX>MAX_LINE_LENGTH) {return;} + + int gix = get_global_id(0); + + if(gix == 0) + { + partialFilter(data, outs, m_Bvalues, m_MMatrix, 0, 1, sizeX); + } + + return; +} + +#endif +*/ + +#ifdef DIM_2 + +__kernel void +YvvFilter(const __global INPIXELTYPE * data, + __global OUTPIXELTYPE * outs, + __constant REALTYPE * m_Bvalues, + __constant REALTYPE * m_MMatrix, + int sizeX, + int sizeY, + unsigned int dim) +{ + if (sizeX > MAX_LINE_LENGTH || sizeY > MAX_LINE_LENGTH) + { + return; + } + + int WIDTH = sizeX; + int HEIGHT = sizeY; + + int gix = get_global_id(0); + + switch (dim) + { + case 0: // X + if (gix < HEIGHT) + { + // start= (id(0)*height + id(1))*width + 0; + partialFilter(data, outs, m_Bvalues, m_MMatrix, gix * WIDTH, 1, WIDTH); + } + break; + case 1: // Y + if (gix < WIDTH) + { + // start= (id(1)*height + 0)*width + id(0); + partialFilter(data, outs, m_Bvalues, m_MMatrix, gix, WIDTH, HEIGHT); + } + break; + } + return; +} + +#endif + + +#ifdef DIM_3 + +__kernel void +YvvFilter(const __global INPIXELTYPE * data, + __global OUTPIXELTYPE * outs, + __constant REALTYPE * m_Bvalues, + __constant REALTYPE * m_MMatrix, + int sizeX, + int sizeY, + int sizeZ, + unsigned int dim // to know which dimension to do +) +{ + if (sizeX > MAX_LINE_LENGTH || sizeY > MAX_LINE_LENGTH) + { + return; + } + + int WIDTH = sizeX; + int HEIGHT = sizeY; + int DEPTH = sizeZ; + + int gix = get_global_id(0); + int giy = get_global_id(1); + /* + __local REALTYPE Bvalues[4]; + __local REALTYPE MMatrix[9]; + + if(gix<9) + { + if(gix<4) + { + Bvalues[gix] = m_Bvalues[gix]; + } + MMatrix[gix] = m_MMatrix[gix]; + } + barrier(CLK_LOCAL_MEM_FENCE); + */ + + switch (dim) + { // generic start= (giz*height + giy)*width + gix; + case 0: // X + if (gix < DEPTH && giy < HEIGHT) + { + // start= (id(0)*height + id(1))*width + 0; + partialFilter(data, outs, m_Bvalues, m_MMatrix, (gix * HEIGHT + giy) * WIDTH, 1, WIDTH); + } + break; + case 1: // Y + if (giy < DEPTH && gix < WIDTH) + { + // start= (id(1)*height + 0)*width + id(0); + partialFilter(data, outs, m_Bvalues, m_MMatrix, giy * WIDTH * HEIGHT + gix, WIDTH, HEIGHT); + } + break; + case 2: // Z + if (gix < WIDTH && giy < HEIGHT) + { + // start= (0*height + id(1))*width + id(0); + partialFilter(data, outs, m_Bvalues, m_MMatrix, giy * WIDTH + gix, HEIGHT * WIDTH, DEPTH); + } + break; + } + return; +} +#endif diff --git a/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/test/CMakeLists.txt b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/test/CMakeLists.txt new file mode 100644 index 00000000000..69fbbf62797 --- /dev/null +++ b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/test/CMakeLists.txt @@ -0,0 +1,123 @@ +itk_module_test() + +if(ITK_USE_GPU) + add_definitions(-DGPU) + + include_directories("${${itk-module}_SOURCE_DIR}/include") + include_directories("${${itk-module}KernelDir}") + + set( + ${itk-module}TestSRC + itkCPURecursiveYvvGaussianImageFilterTest.cxx + itkGPURecursiveYvvGaussianImageFilterTest.cxx + itkYvvGpuCpuSimilarityTest.cxx + itkYvvBenchmark.cxx + itkYvvWhiteImageTest.cxx + ) +else() + #No GPU, but CPU user will probably want double-precision. + add_definitions(-DWITH_DOUBLE) + set( + ${itk-module}TestSRC + itkCPURecursiveYvvGaussianImageFilterTest.cxx + itkYvvBenchmark.cxx + itkYvvWhiteImageTest.cxx + ) +endif() + +#set(ITK_TEST_DRIVER itkTestDriver) + +createtestdriver(${itk-module} "${${itk-module}-Test_LIBRARIES}" "${${itk-module}TestSRC}") + +#Common tests for CPU and/or GPU +itk_add_test( + NAME itkCPURecursiveYvvGaussianImageFilterTest2D + COMMAND + ${itk-module}TestDriver + itkCPURecursiveYvvGaussianImageFilterTest + DATA{Input/512ex.jpg} + 2 + 12.0 + 4 +) + +itk_add_test( + NAME itkYvvBenchmark2D + COMMAND + ${itk-module}TestDriver + itkYvvBenchmark + DATA{Input/512ex.jpg} + 2 + 12.0 + 4 +) + +itk_add_test( + NAME itkYvvBenchmark3D + COMMAND + ${itk-module}TestDriver + itkYvvBenchmark + DATA{Input/256x256x64.tif} + 3 + 12.0 + 2 +) + +itk_add_test( + NAME itkYvvWhiteImageTest2D + COMMAND + ${itk-module}TestDriver + itkYvvWhiteImageTest + 2 + 12.0 + 6 + 512 + 512 +) + +itk_add_test( + NAME itkYvvWhiteImageTest3D + COMMAND + ${itk-module}TestDriver + itkYvvWhiteImageTest + 3 + 12.0 + 2 + 256 + 256 + 16 +) + +if(ITK_USE_GPU) + itk_add_test( + NAME itkGPURecursiveYvvGaussianImageFilterTest2DLow + COMMAND + ${itk-module}TestDriver + itkGPURecursiveYvvGaussianImageFilterTest + DATA{Input/512ex.jpg} + 2 + 0.5 + 1 + ) + + itk_add_test( + NAME itkGPURecursiveYvvGaussianImageFilterTest2D + COMMAND + ${itk-module}TestDriver + itkGPURecursiveYvvGaussianImageFilterTest + DATA{Input/512ex.jpg} + 2 + 20 + 1 + ) + + itk_add_test( + NAME itkYvvGpuCpuSimilarityTest2D + COMMAND + ${itk-module}TestDriver + itkYvvGpuCpuSimilarityTest + DATA{Input/512ex.jpg} + 2 + 14 + ) +endif() diff --git a/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/test/Input/256x256x64.tif b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/test/Input/256x256x64.tif new file mode 100644 index 00000000000..f19e9c33486 Binary files /dev/null and b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/test/Input/256x256x64.tif differ diff --git a/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/test/Input/512ex.jpg b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/test/Input/512ex.jpg new file mode 100644 index 00000000000..b41adabfcd6 Binary files /dev/null and b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/test/Input/512ex.jpg differ diff --git a/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/test/itkCPURecursiveYvvGaussianImageFilterTest.cxx b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/test/itkCPURecursiveYvvGaussianImageFilterTest.cxx new file mode 100644 index 00000000000..95d66a43990 --- /dev/null +++ b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/test/itkCPURecursiveYvvGaussianImageFilterTest.cxx @@ -0,0 +1,64 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +#include "yvvFilter.hxx" +#include +#include + +#ifdef WITH_DOUBLE +using PixelType = double; +#else +using PixelType = float; +#endif + +int +itkCPURecursiveYvvGaussianImageFilterTest(int argc, char * argv[]) +{ + if (argc < 5) + { + std::cerr << "Error: missing arguments" << std::endl; + std::cerr << "[Usage] " << argv[0] << " image_file ndimension sigma num_runs" << std::endl; + return EXIT_FAILURE; + } + + std::string inputFilename(argv[1]); + unsigned int dim = std::stoi(argv[2]); + float sigma = std::stod(argv[3]); + unsigned int ntests = std::stoi(argv[4]); + + int res = EXIT_SUCCESS; + itk::TimeProbesCollectorBase timeCollector; + + if (dim == 2) + { + res = testImage>(inputFilename, sigma, &timeCollector, ntests); + } + else if (dim == 3) + { + res = testImage>(inputFilename, sigma, &timeCollector, ntests); + } + else + { + std::cerr << "Error: only 2 or 3 dimensions allowed, " << dim << " selected." << std::endl; + res = EXIT_FAILURE; + } + + timeCollector.Report(); + std::cout << "\n(!) GPU results will only be shown if GPU support has been detected and activated by the user.\n"; + return res; +} diff --git a/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/test/itkGPURecursiveYvvGaussianImageFilterTest.cxx b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/test/itkGPURecursiveYvvGaussianImageFilterTest.cxx new file mode 100644 index 00000000000..4cedc54216b --- /dev/null +++ b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/test/itkGPURecursiveYvvGaussianImageFilterTest.cxx @@ -0,0 +1,81 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +#include "yvvFilter.hxx" +#include +#include + +#ifdef GPU +# include "itkGPUContextManager.h" +#endif + +#ifdef WITH_DOUBLE +using PixelType = double; +#else +using PixelType = float; +#endif + +int +itkGPURecursiveYvvGaussianImageFilterTest(int argc, char * argv[]) +{ +#ifdef GPU + if (!itk::IsGPUAvailable()) + { + std::cerr << "OpenCL-enabled GPU is not present." << std::endl; + return EXIT_FAILURE; + } +#endif + + if (argc < 5) + { + std::cerr << "Error: missing arguments" << std::endl; + std::cerr << "[Usage] " << argv[0] << " image_file ndimension sigma num_runs" << std::endl; + return EXIT_FAILURE; + } + + std::string inputFilename(argv[1]); + unsigned int dim = std::stoi(argv[2]); + float sigma = std::stod(argv[3]); + unsigned int ntests = std::stoi(argv[4]); + + + int res = EXIT_SUCCESS; + itk::TimeProbesCollectorBase timeCollector; + + if (dim == 2) + { + res = testImage>(inputFilename, sigma, &timeCollector, ntests); + } + else if (dim == 3) + { + res = testImage>(inputFilename, sigma, &timeCollector, ntests); + } + else + { + std::cerr << "Error: only 2 or 3 dimensions allowed, " << dim << " selected." << std::endl; + res = EXIT_FAILURE; + } + + timeCollector.Report(); +#ifndef GPU + std::cout + << "-- ITK GPU support was not detected and/or not configured, so no GPU filters were tested. --\n"; +#endif + + return res; +} diff --git a/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/test/itkYvvBenchmark.cxx b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/test/itkYvvBenchmark.cxx new file mode 100644 index 00000000000..9fbe4166fad --- /dev/null +++ b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/test/itkYvvBenchmark.cxx @@ -0,0 +1,81 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +#include "yvvFilter.hxx" +#include +#include + +#ifdef GPU +# include "itkGPUContextManager.h" +#endif + +#ifdef WITH_DOUBLE +using PixelType = double; +#else +using PixelType = float; +#endif + +int +itkYvvBenchmark(int argc, char * argv[]) +{ +#ifdef GPU + if (!itk::IsGPUAvailable()) + { + std::cerr << "OpenCL-enabled GPU is not present." << std::endl; + return EXIT_FAILURE; + } +#endif + + if (argc < 5) + { + std::cerr << "Error: missing arguments" << std::endl; + std::cerr << "[Usage] " << argv[0] << " image_file ndimension sigma num_runs" << std::endl; + return EXIT_FAILURE; + } + + std::string inputFilename(argv[1]); + unsigned int dim = std::stoi(argv[2]); + float sigma = std::stod(argv[3]); + unsigned int ntests = std::stoi(argv[4]); + + + int res = EXIT_SUCCESS; + itk::TimeProbesCollectorBase timeCollector; + + if (dim == 2) + { + res = testImage>(inputFilename, sigma, &timeCollector, ntests); + } + else if (dim == 3) + { + res = testImage>(inputFilename, sigma, &timeCollector, ntests); + } + else + { + std::cerr << "Error: only 2 or 3 dimensions allowed, " << dim << " selected." << std::endl; + res = EXIT_FAILURE; + } + + timeCollector.Report(); +#ifndef GPU + std::cout + << "-- ITK GPU support was not detected and/or not configured, so no GPU filters were tested. --\n"; +#endif + + return res; +} diff --git a/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/test/itkYvvGpuCpuSimilarityTest.cxx b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/test/itkYvvGpuCpuSimilarityTest.cxx new file mode 100644 index 00000000000..298f68f2fea --- /dev/null +++ b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/test/itkYvvGpuCpuSimilarityTest.cxx @@ -0,0 +1,237 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +#include +#include + +#include "itkCastImageFilter.h" +#include "itkImageFileReader.h" +#include "itkImageFileWriter.h" +#include "itkMinimumMaximumImageCalculator.h" +#include "itkSmoothingRecursiveYvvGaussianImageFilter.h" +#include "itkTimeProbe.h" + +#include "itkGPUContextManager.h" +#include "itkGPUImage.h" +#include "itkGPUKernelManager.h" +#include "itkGPUSmoothingRecursiveYvvGaussianImageFilter.h" + + +// The double precision is not well-supported on most GPUs +// and by many drivers at this time. +// Relax the RMS threshold to allow for errors due to +// differences in precision when the GPU does not support +// double-precision. +// TODO: Test on Mac/Windows. + +#ifdef WITH_DOUBLE +using PixelType = double; +# define NRMSTH 5e-07 +#else +using PixelType = float; +# define NRMSTH 5e-04 +#endif + +template +int +runYvvGpuCpuSimilarityTest(const std::string & inFile, float mySigma) +{ + using InputPixelType = PixelType; + using OutputPixelType = PixelType; + + using InputImageType = itk::GPUImage; + using OutputImageType = itk::GPUImage; + using UnsignedCharImageType = itk::Image; + using CastFilterType = itk::CastImageFilter; + + using CPUYvvFilterType = itk::SmoothingRecursiveYvvGaussianImageFilter; + using GPUYvvFilterType = itk::GPUSmoothingRecursiveYvvGaussianImageFilter; + + using ReaderType = itk::ImageFileReader; + using WriterType = itk::ImageFileWriter; + + typename ReaderType::Pointer reader = ReaderType::New(); + typename WriterType::Pointer writer = WriterType::New(); + typename CastFilterType::Pointer castFilter = CastFilterType::New(); + + typename WriterType::Pointer writerCPU = WriterType::New(); + typename CastFilterType::Pointer castFilterCPU = CastFilterType::New(); + + std::string outFile = "gpuYvvCompOutput"; + std::string outFileCPU = "cpuYvvCompOutput"; + + if (ImageType::ImageDimension == 2) + { + outFile += ".jpg"; + outFileCPU += ".jpg"; + } + else + { + outFile += ".tif"; + outFileCPU += ".tif"; + } + + reader->SetFileName(inFile); + writer->SetFileName(outFile); + writerCPU->SetFileName(outFileCPU); + + typename InputImageType::Pointer src = InputImageType::New(); + src = reader->GetOutput(); + reader->Update(); + + bool someDiffFlag = false; + + std::cout << "\nUsing " << inFile << std::endl; + + for (float sigma = 0.5; sigma <= mySigma; sigma *= 2) + { + typename CPUYvvFilterType::Pointer CPUFilter = CPUYvvFilterType::New(); + + itk::TimeProbe cputimer; + + cputimer.Start(); + CPUFilter->SetInput(src); + CPUFilter->SetSigma(sigma); + CPUFilter->Update(); + cputimer.Stop(); + + + std::cout << "\nSize: " << src->GetLargestPossibleRegion().GetSize() << std::endl; + std::cout << "CPU Recursive YVV Gaussian Filter took " << cputimer.GetMean() << " seconds with " + << CPUFilter->GetNumberOfWorkUnits() << " threads." << std::endl; + + castFilterCPU->SetInput(CPUFilter->GetOutput()); + writerCPU->SetInput(castFilterCPU->GetOutput()); + writerCPU->Update(); + + typename GPUYvvFilterType::Pointer GPUFilter = GPUYvvFilterType::New(); + + itk::TimeProbe gputimer; + + gputimer.Start(); + GPUFilter->SetInput(reader->GetOutput()); + GPUFilter->SetSigma(sigma); + GPUFilter->Update(); + GPUFilter->GetOutput()->UpdateBuffers(); // synchronization point (GPU->CPU memcpy) + gputimer.Stop(); + + using ImageCalculatorFilterType = itk::MinimumMaximumImageCalculator; + + castFilter->SetInput(GPUFilter->GetOutput()); + writer->SetInput(castFilter->GetOutput()); + writer->Update(); + + // --------------- + // RMS Error check + // --------------- + + using ImageCalculatorFilterType = itk::MinimumMaximumImageCalculator; + + typename ImageCalculatorFilterType::Pointer imageCalculatorFilter = ImageCalculatorFilterType::New(); + imageCalculatorFilter->SetImage(CPUFilter->GetOutput()); + imageCalculatorFilter->Compute(); + + double maxPx = imageCalculatorFilter->GetMaximum(); + double minPx = imageCalculatorFilter->GetMinimum(); + double diff = 0; + unsigned int nPix = 0; + itk::ImageRegionIterator cit(CPUFilter->GetOutput(), + CPUFilter->GetOutput()->GetLargestPossibleRegion()); + itk::ImageRegionIterator git(GPUFilter->GetOutput(), + GPUFilter->GetOutput()->GetLargestPossibleRegion()); + + for (cit.GoToBegin(), git.GoToBegin(); !cit.IsAtEnd(); ++cit, ++git) + { + double err = (double)(cit.Get()) - (double)(git.Get()); + // if(err > 0.1 || (double)cit.Get() < 0.1) std::cout << "CPU : " << (double)(cit.Get()) << ", GPU : " + // << (double)(git.Get()) << std::endl; + diff += err * err; + nPix++; + } + + + if (nPix > 0) + { + double NormRMSError = sqrt(diff / (double)nPix) / (maxPx - minPx); // + std::cout << "Normalised RMS Error with sigma = " << sigma << " : " << NormRMSError << std::endl; + + if (std::isnan(NormRMSError)) + { + std::cout << "Normalised RMS Error with sigma = " << sigma << " is NaN! nPix: " << nPix << std::endl; + return EXIT_FAILURE; + } + if (NormRMSError > NRMSTH) + { + std::cout << "Normalised RMS Error with sigma = " << sigma << " exceeds threshold (" << NRMSTH << ")" + << std::endl; + someDiffFlag = true; + } + } + else + { + std::cout << "No pixels in output!" << std::endl; + return EXIT_FAILURE; + } + } + + if (someDiffFlag) + return EXIT_FAILURE; + + return EXIT_SUCCESS; +} + + +int +itkYvvGpuCpuSimilarityTest(int argc, char * argv[]) +{ + if (!itk::IsGPUAvailable()) + { + std::cerr << "OpenCL-enabled GPU is not present." << std::endl; + return EXIT_FAILURE; + } + + if (argc < 3) + { + std::cerr << argc << ": Error: missing arguments: " << argv[1] << std::endl; + std::cerr << "[Usage] " << argv[0] << " image_file ndimension [sigma]" << std::endl; + return EXIT_FAILURE; + } + + std::string inputFilename(argv[1]); + unsigned int dim = std::stoi(argv[2]); + + float sigma = 0.5; + if (argc > 3 && std::stoi(argv[3]) >= 0.5) + { + sigma = std::stod(argv[3]); + } + + if (dim == 2) + { + return runYvvGpuCpuSimilarityTest>(inputFilename, sigma); + } + else if (dim == 3) + { + return runYvvGpuCpuSimilarityTest>(inputFilename, sigma); + } + else + { + std::cerr << "Error: only 2 or 3 dimensions allowed, " << dim << " selected." << std::endl; + } + return EXIT_FAILURE; +} diff --git a/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/test/itkYvvWhiteImageTest.cxx b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/test/itkYvvWhiteImageTest.cxx new file mode 100644 index 00000000000..e0261f005f2 --- /dev/null +++ b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/test/itkYvvWhiteImageTest.cxx @@ -0,0 +1,119 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +#include "yvvFilter.hxx" +#include +#include + +#ifdef GPU +# include "itkGPUContextManager.h" +#endif + +#ifdef WITH_DOUBLE +using PixelType = double; +#else +using PixelType = float; +#endif + +#define die(error_msg) \ + std::cerr << "Error: " << error_msg << std::endl; \ + std::cerr << "Usage: " << argv[0] << " ndimension sigma num_runs width [height] [depth]" << std::endl; \ + return EXIT_FAILURE; + +int +itkYvvWhiteImageTest(int argc, char * argv[]) +{ +#ifdef GPU + if (!itk::IsGPUAvailable()) + { + std::cerr << "OpenCL-enabled GPU is not present." << std::endl; + return EXIT_FAILURE; + } +#endif + + if (argc < 4) + { + die("missing arguments."); + } + + int dim = std::stoi(argv[1]); + if (argc < 4 + dim - 1) + { + die("missing arguments."); + } + if (dim == 3 && argc != 7) + { + die("missing arguments for a 3D image."); + } + + unsigned int ntests; + auto * size = new unsigned int[dim]; + float sigma; + try + { + sigma = std::stod(argv[2]); + ntests = std::stoi(argv[3]); + for (int i = 0; i < dim; ++i) + { + size[i] = std::stoi(argv[4 + i]); + } + } + catch (...) + { + delete[] size; + die("invalid size arguments."); + } + + std::cout << "::: Received: dim =" << dim << ", sigma= " << sigma << ", ntests= " << ntests << ", size= "; + for (int i = 0; i < dim; ++i) + { + std::cout << size[i] << " "; + } + std::cout << ":::" << std::endl << std::endl; + + int result; + + itk::TimeProbesCollectorBase timeCollector; + if (dim == 2) + { + using ImageType = itk::Image; + ImageType::SizeType size2D = { { size[0], size[1] } }; + result = testWhite(size2D, sigma, &timeCollector, ntests); + } + else if (dim == 3) + { + using ImageType = itk::Image; + ImageType::SizeType size3D = { { size[0], size[1], size[2] } }; + result = testWhite(size3D, sigma, &timeCollector, ntests); + } + else + { + std::cerr << "Error: only 2 or 3 dimensions allowed, " << dim << " selected." << std::endl; + result = EXIT_FAILURE; + } + timeCollector.Report(); + + delete[] size; + +#ifndef GPU + std::cout << "-- ITK GPU support was not detected and/or not configured, so no GPU filters were tested. --" + << std::endl; +#endif + + return result; +} diff --git a/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/test/yvvFilter.hxx b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/test/yvvFilter.hxx new file mode 100644 index 00000000000..d3498e494bc --- /dev/null +++ b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/test/yvvFilter.hxx @@ -0,0 +1,356 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +#include "itkCastImageFilter.h" +#include "itkImageFileReader.h" +#include "itkImageFileWriter.h" +#include "itkSmoothingRecursiveGaussianImageFilter.h" +#include "itkSmoothingRecursiveYvvGaussianImageFilter.h" +#include "itkTimeProbesCollectorBase.h" + +#ifdef GPU +# include "itkGPUImage.h" +# include "itkGPUSmoothingRecursiveYvvGaussianImageFilter.h" +#endif + + +template +void +writeImage(std::string filterLabel, ImageType * result) +{ + using UnsignedCharImageType = itk::Image; + using CastFilterType = itk::CastImageFilter; + using WriterType = itk::ImageFileWriter; + + std::string outputFilename = filterLabel; + + typename CastFilterType::Pointer castFilter = CastFilterType::New(); + castFilter->SetInput(result); + castFilter->Update(); + + typename WriterType::Pointer writer = WriterType::New(); + + if (ImageType::ImageDimension == 2) + { + outputFilename += ".jpg"; + } + else + { + outputFilename += ".tif"; + } + writer->SetFileName(outputFilename); + writer->SetInput(castFilter->GetOutput()); + try + { + writer->Update(); + } + catch (itk::ExceptionObject & e) + { + std::cerr << e << std::endl; + } +} + + +template +int +getSourceImage(void * sourceImagePtr, std::string inputFilename) +{ + using ReaderType = itk::ImageFileReader; + typename ReaderType::Pointer readerGPU = ReaderType::New(); + readerGPU->SetFileName(inputFilename); + try + { + readerGPU->Update(); + } + catch (const itk::ImageFileReaderException & e) + { + std::cout << "Error : " << e.GetDescription() << std::endl; + return -1; + } + + auto * kimg = (typename ImageType::Pointer *)sourceImagePtr; + *kimg = readerGPU->GetOutput(); + + (*kimg)->DisconnectPipeline(); + return EXIT_SUCCESS; +} + + +template +int +createWhiteImage(void * sourceImagePtr, typename ImageType::SizeType size) +{ + auto * kimg = (typename ImageType::Pointer *)sourceImagePtr; + typename ImageType::RegionType region; + region.SetSize(size); + + // create + *kimg = ImageType::New(); + (*kimg)->SetRegions(region); + (*kimg)->Allocate(); + (*kimg)->FillBuffer(255.0); + + return EXIT_SUCCESS; +} + + +template +int +testCpuFilter(std::string & filterLabel, + std::string & inputFilename, + typename FilterType::InputImageType::SizeType size, + float sigma, + std::string parameters, + itk::TimeProbesCollectorBase * timeCollector) +{ + using InputImage = typename FilterType::InputImageType; + typename InputImage::Pointer src; + void * imgPtr = &src; + + if (inputFilename.empty()) + { + createWhiteImage(imgPtr, size); + // createStepImage< InputImage >( imgPtr, size ); + + std::ostringstream sizeStream; + + sizeStream << size[0]; + for (unsigned int i = 1; i < InputImage::ImageDimension; ++i) + { + sizeStream << "x" << size[i]; + } + parameters += sizeStream.str(); + } + else + { + getSourceImage(imgPtr, inputFilename); + } + + typename FilterType::Pointer filter = FilterType::New(); + filter->SetSigma(sigma); + filter->SetInput(src); + filter->Update(); + + { + src->Modified(); + filter->Modified(); + + timeCollector->Start(filterLabel.c_str()); + filter->Update(); + filter->GetOutput(); + timeCollector->Stop(filterLabel.c_str()); + + writeImage(filterLabel + parameters, filter->GetOutput()); + } + + src->DisconnectPipeline(); + src = nullptr; + filter = nullptr; + imgPtr = nullptr; + return EXIT_SUCCESS; +} + + +#ifdef GPU + +template +int +testGpuFilter(std::string & filterLabel, + std::string & inputFilename, + typename FilterType::InputImageType::SizeType size, + float sigma, + std::string parameters, + itk::TimeProbesCollectorBase * timeCollector, + bool measureWithSync) +{ + using InputImage = typename FilterType::InputImageType; + + typename InputImage::Pointer src; + void * imgPtr = &src; + + if (inputFilename.empty()) + { + createWhiteImage(imgPtr, size); + // createStepImage(imgPtr, size); + + std::ostringstream sizeStream; + + sizeStream << size[0]; + for (int i = 1; i < InputImage::ImageDimension; ++i) + { + sizeStream << "x" << size[i]; + } + parameters += sizeStream.str(); + } + else + { + getSourceImage(imgPtr, inputFilename); + } + + typename FilterType::Pointer filter = FilterType::New(); + filter->SetSigma(sigma); + filter->SetInput(src); + filter->Update(); + filter->GetOutput()->UpdateBuffers(); // force first CPU-GPU synchronisation + + { + // src->Modified(); //No need, filter->Update forces call to GPUGenerateData() + filter->Modified(); + if (measureWithSync) + { + timeCollector->Start(filterLabel.c_str()); + filter->Update(); + filter->GetOutput()->UpdateBuffers(); + timeCollector->Stop(filterLabel.c_str()); + } + else + { + timeCollector->Start(filterLabel.c_str()); + filter->Update(); + timeCollector->Stop(filterLabel.c_str()); + filter->GetOutput()->UpdateBuffers(); + } + writeImage(filterLabel + parameters, filter->GetOutput()); + } + + src->DisconnectPipeline(); + src = nullptr; + filter = nullptr; + imgPtr = nullptr; + return EXIT_SUCCESS; +} + +#endif + + +template +int +testImage(std::string inputFilename, float sigma, itk::TimeProbesCollectorBase * timeCollector, unsigned int ntests) +{ +#ifdef GPU + using GPUImageType = itk::GPUImage; + using GPUrecursiveYVVFilterType = itk::GPUSmoothingRecursiveYvvGaussianImageFilter; +#endif + using CPUImageType = ImageType; + using RecursiveYVVFilterType = itk::SmoothingRecursiveYvvGaussianImageFilter; + using DericheFilterType = itk::SmoothingRecursiveGaussianImageFilter; + + std::cout << ":::: Testing on " << inputFilename << ", using sigma = " << sigma << " ::::" << std::endl; + typename ImageType::SizeType size; + size.Fill(0); // empty, since we'll be using an actual file + std::ostringstream parameterStream; + parameterStream << "_s" << sigma << "_"; + + //////////////////////////////// + std::string dericheLabel = "cpu_Deriche"; + for (unsigned int i = 0; i < ntests; ++i) + { + testCpuFilter(dericheLabel, inputFilename, size, sigma, parameterStream.str(), timeCollector); + } + //////////////////////////////// + std::string cpuYvvLabel = "cpu_Yvv"; + for (unsigned int i = 0; i < ntests; ++i) + { + testCpuFilter( + cpuYvvLabel, inputFilename, size, sigma, parameterStream.str(), timeCollector); + } + +////////////////////////////// +#ifdef GPU + std::string gpuYvvLabel = "gpu_Yvv_wSync"; + for (unsigned int i = 0; i < ntests; ++i) + { + testGpuFilter( + gpuYvvLabel, inputFilename, size, sigma, parameterStream.str(), timeCollector, true); + /* Stress tests (e.g. high ntests) may cause CL_MEM_OBJECT_ALLOCATION_FAILURE, + * even when re-instantiating input _and_ filter. Not all GPUs support resetting, so adding a sleep helps.*/ + sleep(1); + } + + std::string gpuYvvLabelNoSync = "gpu_Yvv_NoSyn"; + for (unsigned int i = 0; i < ntests; ++i) + { + testGpuFilter( + gpuYvvLabelNoSync, inputFilename, size, sigma, parameterStream.str(), timeCollector, false); + sleep(1); // helps for stress tests + } + +#endif + return EXIT_SUCCESS; +} + + +template +int +testWhite(typename ImageType::SizeType size, + float sigma, + itk::TimeProbesCollectorBase * timeCollector, + unsigned int ntests) +{ + std::cout << "Testing: " << size << " with sigma = " << sigma << ". Average over " << ntests << " runs." << std::endl; + +#ifdef GPU + using GPUImageType = itk::GPUImage; + using GPUrecursiveYVVFilterType = itk::GPUSmoothingRecursiveYvvGaussianImageFilter; +#endif + using CPUImageType = ImageType; + using RecursiveYVVFilterType = itk::SmoothingRecursiveYvvGaussianImageFilter; + using DericheFilterType = itk::SmoothingRecursiveGaussianImageFilter; + + std::string emptyFilename; + std::ostringstream parameterStream; + parameterStream << "_s" << sigma << "_"; + + std::cout << ":::: Testing on white image: " << size << ", using sigma = " << sigma << " ::::" << std::endl; + + //////////////////////////////// + std::string dericheLabel = "cpu_Deriche"; + for (unsigned int i = 0; i < ntests; ++i) + { + testCpuFilter(dericheLabel, emptyFilename, size, sigma, parameterStream.str(), timeCollector); + } + //////////////////////////////// + std::string cpuYvvLabel = "cpu_Yvv"; + for (unsigned int i = 0; i < ntests; ++i) + { + testCpuFilter( + cpuYvvLabel, emptyFilename, size, sigma, parameterStream.str(), timeCollector); + } + +////////////////////////////// +#ifdef GPU + std::string gpuYvvLabelSync = "gpu_Yvv_wSync"; + for (unsigned int i = 0; i < ntests; ++i) + { + testGpuFilter( + gpuYvvLabelSync, emptyFilename, size, sigma, parameterStream.str(), timeCollector, true); + /* Stress tests (e.g. high ntests) may cause CL_MEM_OBJECT_ALLOCATION_FAILURE, + * even when re-instantiating input _and_ filter. Not all GPUs support resetting, so adding a sleep helps.*/ + sleep(1); + } + + std::string gpuYvvLabelNoSync = "gpu_Yvv_NoSyn"; + for (unsigned int i = 0; i < ntests; ++i) + { + testGpuFilter( + gpuYvvLabelNoSync, emptyFilename, size, sigma, parameterStream.str(), timeCollector, false); + sleep(1); // helps for stress tests + } +#endif + return EXIT_SUCCESS; +} diff --git a/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/wrapping/CMakeLists.txt b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/wrapping/CMakeLists.txt new file mode 100644 index 00000000000..cd6f4293c3e --- /dev/null +++ b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/wrapping/CMakeLists.txt @@ -0,0 +1,3 @@ +itk_wrap_module(SmoothingRecursiveYvvGaussianFilter) +itk_auto_load_submodules() +itk_end_wrap_module() diff --git a/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/wrapping/itkSmoothingRecursiveYvvGaussianImageFilter.wrap b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/wrapping/itkSmoothingRecursiveYvvGaussianImageFilter.wrap new file mode 100644 index 00000000000..a6f80d241a6 --- /dev/null +++ b/Modules/Filtering/SmoothingRecursiveYvvGaussianFilter/wrapping/itkSmoothingRecursiveYvvGaussianImageFilter.wrap @@ -0,0 +1,3 @@ +itk_wrap_class("itk::SmoothingRecursiveYvvGaussianImageFilter" POINTER) +itk_wrap_image_filter("${WRAP_ITK_SCALAR}" 2 2+) +itk_end_wrap_class() diff --git a/Modules/Remote/SmoothingRecursiveYvvGaussianFilter.remote.cmake b/Modules/Remote/SmoothingRecursiveYvvGaussianFilter.remote.cmake deleted file mode 100644 index cc4ccccc0d6..00000000000 --- a/Modules/Remote/SmoothingRecursiveYvvGaussianFilter.remote.cmake +++ /dev/null @@ -1,52 +0,0 @@ -#-- # Grading Level Criteria Report -#-- EVALUATION DATE: 2020-03-01 -#-- EVALUATORS: [<>,<>] -#-- -#-- ## Compliance level 5 star (AKA ITK main modules, or remote modules that could become core modules) -#-- - [ ] Widespread community dependance -#-- - [ ] Above 90% code coverage -#-- - [ ] CI dashboards and testing monitored rigorously -#-- - [ ] Key API features are exposed in wrapping interface -#-- - [ ] All requirements of Levels 4,3,2,1 -#-- -#-- ## Compliance Level 4 star (Very high-quality code, perhaps small community dependance) -#-- - [ ] Meets all ITK code style standards -#-- - [ ] No external requirements beyond those needed by ITK proper -#-- - [ ] Builds and passes tests on all supported platforms within 1 month of each core tagged release -#-- - [ ] Windows Shared Library Build with Visual Studio -#-- - [ ] Mac with clang compiller -#-- - [ ] Linux with gcc compiler -#-- - [ ] Active developer community dedicated to maintaining code-base -#-- - [ ] 75% code coverage demonstrated for testing suite -#-- - [ ] Continuous integration testing performed -#-- - [ ] All requirements of Levels 3,2,1 -#-- -#-- ## Compliance Level 3 star (Quality beta code) -#-- - [ ] API | executable interface is considered mostly stable and feature complete -#-- - [ ] 10% C0-code coverage demonstrated for testing suite -#-- - [ ] Some tests exist and pass on at least some platform -#-- - [X] All requirements of Levels 2,1 -#-- -#-- ## Compliance Level 2 star (Alpha code feature API development or niche community/execution environment dependance ) -#-- - [X] Compiles for at least 1 niche set of execution envirionments, and perhaps others -#-- (may depend on specific external tools like a java environment, or specific external libraries to work ) -#-- - [X] All requirements of Levels 1 -#-- -#-- ## Compliance Level 1 star (Pre-alpha features under development and code of unknown quality) -#-- - [X] Code complies on at least 1 platform -#-- -#-- ## Compliance Level 0 star ( Code/Feature of known poor-quality or deprecated status ) -#-- - [ ] Code reviewed and explicitly identified as not recommended for use -#-- -#-- ### Please document here any justification for the criteria above -# Code style enforced by clang-format on 2020-02-19, and clang-tidy modernizations completed - -#Contact: I. Vidal-Migallon -itk_fetch_module( - SmoothingRecursiveYvvGaussianFilter - "GPU and CPU Young & Van Vliet Recursive Gaussian Smoothing Filter: https://doi.org/10.54294/cpyaig" - MODULE_COMPLIANCE_LEVEL 2 - #UPSTREAM_REPO GIT_REPOSITORY https://github.com/Inria-Asclepios/SmoothingRecursiveYvvGaussianFilter - GIT_REPOSITORY https://github.com/InsightSoftwareConsortium/ITKSmoothingRecursiveYvvGaussianFilter.git - GIT_TAG 8cfb58e73158d6ab779d2171b624c73929c0b35d - ) diff --git a/pyproject.toml b/pyproject.toml index 876691e7964..a7a073a93a1 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -64,6 +64,7 @@ cmd = '''cmake \ -DModule_MGHIO:BOOL=ON \ -DModule_PolarTransform:BOOL=ON \ -DModule_PrincipalComponentsAnalysis:BOOL=ON \ + -DModule_SmoothingRecursiveYvvGaussianFilter:BOOL=ON \ -DModule_SplitComponents:BOOL=ON \ -DModule_IOMeshMZ3:BOOL=ON \ -DModule_IOFDF:BOOL=ON \