diff --git a/Modules/Filtering/Smoothing/include/itkSmoothingRecursiveGaussianImageFilter.h b/Modules/Filtering/Smoothing/include/itkSmoothingRecursiveGaussianImageFilter.h index e7834720cb3..7d287bcc258 100644 --- a/Modules/Filtering/Smoothing/include/itkSmoothingRecursiveGaussianImageFilter.h +++ b/Modules/Filtering/Smoothing/include/itkSmoothingRecursiveGaussianImageFilter.h @@ -37,6 +37,10 @@ namespace itk * filters. For multi-component images, the filter works on each * component independently. * + * Axes with fewer than four pixels are not smoothed. If all axes are shorter + * than four pixels, the input is only cast to the output pixel type. + * Image dimensionality and physical geometry are preserved. + * * For this filter to be able to run in-place the input and output * image types need to be the same and/or the same type as the * RealImageType. diff --git a/Modules/Filtering/Smoothing/include/itkSmoothingRecursiveGaussianImageFilter.hxx b/Modules/Filtering/Smoothing/include/itkSmoothingRecursiveGaussianImageFilter.hxx index 33a7eebe557..0afa12f222e 100644 --- a/Modules/Filtering/Smoothing/include/itkSmoothingRecursiveGaussianImageFilter.hxx +++ b/Modules/Filtering/Smoothing/include/itkSmoothingRecursiveGaussianImageFilter.hxx @@ -27,12 +27,7 @@ template SmoothingRecursiveGaussianImageFilter::SmoothingRecursiveGaussianImageFilter() : m_FirstSmoothingFilter(FirstGaussianFilterType::New()) { - // NB: The first filter is the last dimension because it does not - // always run in-place. As this dimension provides the least amount - // of cache coherency, it will provide the least amount of performance - // gains from running in-place. In fact, some performance tests - // indicate that the running in-place in the 3rd dimension the - // performance actually declines compared to not in-place methods. + // The highest eligible dimension benefits least from in-place execution, so process it first. m_FirstSmoothingFilter->SetOrder(GaussianOrderEnum::ZeroOrder); m_FirstSmoothingFilter->SetDirection(ImageDimension - 1); m_FirstSmoothingFilter->SetNormalizeAcrossScale(m_NormalizeAcrossScale); @@ -49,14 +44,7 @@ SmoothingRecursiveGaussianImageFilter::SmoothingRecur 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.back()->GetOutput()); m_CastingFilter->InPlaceOn(); this->InPlaceOff(); @@ -202,19 +190,68 @@ SmoothingRecursiveGaussianImageFilter::GenerateData() const typename TInputImage::ConstPointer inputImage(this->GetInput()); - const typename TInputImage::RegionType region = inputImage->GetRequestedRegion(); - const typename TInputImage::SizeType size = region.GetSize(); - + const auto size = inputImage->GetLargestPossibleRegion().GetSize(); + unsigned int numberOfActiveDimensions = 0; + unsigned int firstDirection = 0; for (unsigned int d = 0; d < ImageDimension; ++d) { - if (size[d] < 4) + if (m_Sigma[d] <= 0.0) { - 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."); + itkExceptionStringMacro("Sigma must be greater than zero."); + } + if (size[d] == 0) + { + itkExceptionMacro("The number of pixels along dimension " << d << " is zero."); + } + if (size[d] >= 4) + { + firstDirection = d; + ++numberOfActiveDimensions; } } + auto progress = ProgressAccumulator::New(); + progress->SetMiniPipelineFilter(this); + + m_FirstSmoothingFilter->SetInput(inputImage); + + if (numberOfActiveDimensions == 0) + { + m_CastingFilter->SetInput(nullptr); + m_CastingFilter->GetOutput()->ReleaseData(); + auto cast = CastImageFilter::New(); + cast->SetInput(inputImage); + cast->SetInPlace(this->GetInPlace()); + cast->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits()); + progress->RegisterInternalFilter(cast, 1.0); + cast->GraftOutput(this->GetOutput()); + cast->Update(); + this->GraftOutput(cast->GetOutput()); + return; + } + + m_FirstSmoothingFilter->SetDirection(firstDirection); + m_FirstSmoothingFilter->SetSigma(m_Sigma[firstDirection]); + progress->RegisterInternalFilter(m_FirstSmoothingFilter, 1.0 / numberOfActiveDimensions); + + typename RealImageType::Pointer previousOutput = m_FirstSmoothingFilter->GetOutput(); + for (unsigned int d = 0; d < ImageDimension - 1; ++d) + { + if (d != firstDirection && size[d] >= 4) + { + auto & filter = m_SmoothingFilters[d]; + filter->SetSigma(m_Sigma[d]); + filter->SetInput(previousOutput); + progress->RegisterInternalFilter(filter, 1.0 / numberOfActiveDimensions); + previousOutput = filter->GetOutput(); + } + else + { + m_SmoothingFilters[d]->SetInput(nullptr); + } + } + m_CastingFilter->SetInput(previousOutput); + // 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()) @@ -238,20 +275,6 @@ SmoothingRecursiveGaussianImageFilter::GenerateData() this->GetOutput()->ReleaseData(); } - // Create a process accumulator for tracking the progress of this minipipeline. - auto progress = ProgressAccumulator::New(); - progress->SetMiniPipelineFilter(this); - - // Register the filter with the with progress accumulator using - // equal weight proportion. - for (const auto & filter : m_SmoothingFilters) - { - progress->RegisterInternalFilter(filter, 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, and the bulk data which be be from the input due // to the in-place option. diff --git a/Modules/Filtering/Smoothing/test/CMakeLists.txt b/Modules/Filtering/Smoothing/test/CMakeLists.txt index 19f3cd69ba6..fbd1ffcad11 100644 --- a/Modules/Filtering/Smoothing/test/CMakeLists.txt +++ b/Modules/Filtering/Smoothing/test/CMakeLists.txt @@ -399,5 +399,6 @@ set( itkBoxSigmaImageFilterGTest.cxx itkMeanImageFilterGTest.cxx itkMedianImageFilterGTest.cxx + itkSmoothingRecursiveGaussianImageFilterGTest.cxx ) creategoogletestdriver(ITKSmoothing "${ITKSmoothing-Test_LIBRARIES}" "${ITKSmoothingGTests}") diff --git a/Modules/Filtering/Smoothing/test/itkSmoothingRecursiveGaussianImageFilterGTest.cxx b/Modules/Filtering/Smoothing/test/itkSmoothingRecursiveGaussianImageFilterGTest.cxx new file mode 100644 index 00000000000..81bd53f5739 --- /dev/null +++ b/Modules/Filtering/Smoothing/test/itkSmoothingRecursiveGaussianImageFilterGTest.cxx @@ -0,0 +1,247 @@ +/*========================================================================= + * + * 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 "itkSmoothingRecursiveGaussianImageFilter.h" +#include "itkImageRegionIteratorWithIndex.h" +#include "itkImageRegionIterator.h" +#include "itkImageRegionConstIterator.h" +#include "itkVectorImage.h" + +#include + +namespace +{ +using ImageType = itk::Image; +using FilterType = itk::SmoothingRecursiveGaussianImageFilter; + +ImageType::Pointer +MakeImage(const ImageType::SizeType & size) +{ + auto image = ImageType::New(); + image->SetRegions(ImageType::RegionType({ { -2, 3, 5 } }, size)); + image->SetOrigin(ImageType::PointType{ { 7.0, -11.0, 2.0 } }); + image->SetSpacing(ImageType::SpacingType{ { 0.7, 1.3, 2.1 } }); + ImageType::DirectionType direction; + direction.SetIdentity(); + direction[0][0] = 0; + direction[0][1] = -1; + direction[1][0] = 1; + direction[1][1] = 0; + image->SetDirection(direction); + image->Allocate(); + for (itk::ImageRegionIteratorWithIndex it(image, image->GetLargestPossibleRegion()); !it.IsAtEnd(); ++it) + { + const auto index = it.GetIndex(); + it.Set(static_cast((index[0] * index[0] + 3 * index[1] + index[2] * index[2]) % 17)); + } + return image; +} + +ImageType::Pointer +SmoothEligibleAxes(const ImageType * input, const FilterType::SigmaArrayType & sigma) +{ + ImageType::ConstPointer previous = input; + ImageType::Pointer output; + const auto size = input->GetLargestPossibleRegion().GetSize(); + unsigned int firstDirection = 0; + bool hasActiveDirection = false; + for (unsigned int d = 0; d < 3; ++d) + { + if (size[d] >= 4) + { + firstDirection = d; + hasActiveDirection = true; + } + } + if (!hasActiveDirection) + { + auto copy = itk::CastImageFilter::New(); + copy->SetInput(input); + copy->Update(); + output = copy->GetOutput(); + output->DisconnectPipeline(); + return output; + } + for (const unsigned int d : { firstDirection, 0u, 1u }) + { + if (size[d] < 4 || (d == firstDirection && output.IsNotNull())) + { + continue; + } + auto recursive = itk::RecursiveGaussianImageFilter::New(); + recursive->SetInput(previous); + recursive->SetDirection(d); + recursive->SetSigma(sigma[d]); + recursive->Update(); + output = recursive->GetOutput(); + output->DisconnectPipeline(); + previous = output; + } + return output; +} + +void +ExpectSameImage(const ImageType * actual, const ImageType * expected) +{ + ASSERT_EQ(actual->GetLargestPossibleRegion(), expected->GetLargestPossibleRegion()); + EXPECT_EQ(actual->GetOrigin(), expected->GetOrigin()); + EXPECT_EQ(actual->GetSpacing(), expected->GetSpacing()); + EXPECT_EQ(actual->GetDirection(), expected->GetDirection()); + itk::ImageRegionConstIterator actualIt(actual, actual->GetLargestPossibleRegion()); + itk::ImageRegionConstIterator expectedIt(expected, expected->GetLargestPossibleRegion()); + for (; !actualIt.IsAtEnd(); ++actualIt, ++expectedIt) + { + EXPECT_FLOAT_EQ(actualIt.Get(), expectedIt.Get()); + } +} +} // namespace + +TEST(SmoothingRecursiveGaussianImageFilter, ShortAxesAndRepeatedUpdates) +{ + const FilterType::SigmaArrayType sigma{ { 0.8, 1.7, 2.5 } }; + for (const bool inPlace : { false, true }) + { + auto filter = FilterType::New(); + filter->SetSigmaArray(sigma); + filter->SetInPlace(inPlace); + for (const unsigned int x : { 1u, 2u, 3u, 4u, 9u }) + { + for (const unsigned int y : { 1u, 2u, 3u, 4u, 9u }) + { + for (const unsigned int z : { 1u, 2u, 3u, 4u, 9u }) + { + const ImageType::SizeType size{ { x, y, z } }; + SCOPED_TRACE(::testing::Message() << "size=" << size << ", inPlace=" << inPlace); + const auto image = MakeImage(size); + const auto expected = SmoothEligibleAxes(image, sigma); + filter->SetInput(image); + ASSERT_NO_THROW(filter->Update()); + ExpectSameImage(filter->GetOutput(), expected); + } + } + } + } +} + +TEST(SmoothingRecursiveGaussianImageFilter, AllShortAxesPreserveDoublePrecisionAndCast) +{ + using DoubleImageType = itk::Image; + auto image = DoubleImageType::New(); + image->SetRegions(DoubleImageType::SizeType{ { 1, 2, 3 } }); + image->Allocate(); + constexpr double value = 1.0000000000000002; + image->FillBuffer(value); + for (const bool inPlace : { false, true }) + { + auto filter = itk::SmoothingRecursiveGaussianImageFilter::New(); + filter->SetInput(image); + filter->SetInPlace(inPlace); + ASSERT_NO_THROW(filter->Update()); + EXPECT_EQ(filter->GetOutput()->GetPixel({ { 0, 0, 0 } }), value); + } + + using OutputImageType = itk::Image; + const auto input = MakeImage({ { 2, 3, 1 } }); + auto filter = itk::SmoothingRecursiveGaussianImageFilter::New(); + filter->SetInput(input); + ASSERT_NO_THROW(filter->Update()); + itk::ImageRegionConstIterator inputIt(input, input->GetLargestPossibleRegion()); + itk::ImageRegionConstIterator outputIt(filter->GetOutput(), input->GetLargestPossibleRegion()); + for (; !inputIt.IsAtEnd(); ++inputIt, ++outputIt) + { + EXPECT_EQ(outputIt.Get(), static_cast(inputIt.Get())); + } +} + +TEST(SmoothingRecursiveGaussianImageFilter, VectorImageShortAxes) +{ + using VectorImageType = itk::VectorImage; + for (const ImageType::SizeType size : { ImageType::SizeType{ { 9, 2, 1 } }, ImageType::SizeType{ { 1, 2, 3 } } }) + { + const auto scalarImage = MakeImage(size); + auto vectorImage = VectorImageType::New(); + vectorImage->CopyInformation(scalarImage); + vectorImage->SetRegions(scalarImage->GetLargestPossibleRegion()); + vectorImage->SetNumberOfComponentsPerPixel(2); + vectorImage->Allocate(); + itk::ImageRegionConstIterator scalarIt(scalarImage, scalarImage->GetLargestPossibleRegion()); + itk::ImageRegionIterator vectorIt(vectorImage, vectorImage->GetLargestPossibleRegion()); + for (; !scalarIt.IsAtEnd(); ++scalarIt, ++vectorIt) + { + VectorImageType::PixelType pixel(2); + pixel[0] = scalarIt.Get(); + pixel[1] = 2 * scalarIt.Get(); + vectorIt.Set(pixel); + } + auto filter = itk::SmoothingRecursiveGaussianImageFilter::New(); + filter->SetInput(vectorImage); + ASSERT_NO_THROW(filter->Update()); + const auto expected = SmoothEligibleAxes(scalarImage, FilterType::SigmaArrayType(1.0)); + itk::ImageRegionConstIterator expectedIt(expected, expected->GetLargestPossibleRegion()); + itk::ImageRegionConstIterator outputIt(filter->GetOutput(), expected->GetLargestPossibleRegion()); + for (; !expectedIt.IsAtEnd(); ++expectedIt, ++outputIt) + { + ASSERT_EQ(outputIt.Get().GetSize(), 2u); + EXPECT_FLOAT_EQ(outputIt.Get()[0], expectedIt.Get()); + EXPECT_FLOAT_EQ(outputIt.Get()[1], 2 * expectedIt.Get()); + } + } +} + +TEST(SmoothingRecursiveGaussianImageFilter, IntegerPixelsWithShortFirstPipelineAxis) +{ + using IntegerImageType = itk::Image; + const auto image = MakeImage({ { 9, 9, 1 } }); + const FilterType::SigmaArrayType sigma{ { 0.8, 1.7, 2.5 } }; + const auto expected = SmoothEligibleAxes(image, sigma); + for (const bool inPlace : { false, true }) + { + auto cast = itk::CastImageFilter::New(); + cast->SetInput(image); + cast->Update(); + auto filter = itk::SmoothingRecursiveGaussianImageFilter::New(); + filter->SetInput(cast->GetOutput()); + filter->SetSigmaArray(sigma); + filter->SetInPlace(inPlace); + ASSERT_NO_THROW(filter->Update()); + itk::ImageRegionConstIterator expectedIt(expected, expected->GetLargestPossibleRegion()); + itk::ImageRegionConstIterator outputIt(filter->GetOutput(), expected->GetLargestPossibleRegion()); + for (; !expectedIt.IsAtEnd(); ++expectedIt, ++outputIt) + { + EXPECT_EQ(outputIt.Get(), static_cast(expectedIt.Get())); + } + } +} + +TEST(SmoothingRecursiveGaussianImageFilter, InvalidSigmaAndEmptyAxesStillThrow) +{ + for (const ImageType::SizeType size : { ImageType::SizeType{ { 9, 9, 1 } }, ImageType::SizeType{ { 1, 2, 3 } } }) + { + const auto image = MakeImage(size); + for (const double sigma : { 0.0, -1.0 }) + { + auto filter = FilterType::New(); + filter->SetInput(image); + filter->SetSigmaArray(FilterType::SigmaArrayType{ { 1.0, 1.0, sigma } }); + EXPECT_THROW(filter->Update(), itk::ExceptionObject); + } + } + auto filter = FilterType::New(); + filter->SetInput(MakeImage({ { 0, 9, 9 } })); + EXPECT_THROW(filter->Update(), itk::ExceptionObject); +} diff --git a/Modules/Registration/RegistrationMethodsv4/itk-module.cmake b/Modules/Registration/RegistrationMethodsv4/itk-module.cmake index be828f7905d..5ef2ef3fbf8 100644 --- a/Modules/Registration/RegistrationMethodsv4/itk-module.cmake +++ b/Modules/Registration/RegistrationMethodsv4/itk-module.cmake @@ -10,6 +10,7 @@ itk_module( ITKOptimizersv4 ITKMetricsv4 TEST_DEPENDS + ITKGoogleTest ITKTestKernel DESCRIPTION "${DOCUMENTATION}" ) diff --git a/Modules/Registration/RegistrationMethodsv4/test/CMakeLists.txt b/Modules/Registration/RegistrationMethodsv4/test/CMakeLists.txt index 1433166182f..cf199cc51d6 100644 --- a/Modules/Registration/RegistrationMethodsv4/test/CMakeLists.txt +++ b/Modules/Registration/RegistrationMethodsv4/test/CMakeLists.txt @@ -534,3 +534,7 @@ set_tests_properties( COST 30 ) + +creategoogletestdriver(ITKRegistrationMethodsv4 "${ITKRegistrationMethodsv4-Test_LIBRARIES}" + itkImageRegistrationMethodv4ShortDimensionsGTest.cxx +) diff --git a/Modules/Registration/RegistrationMethodsv4/test/itkImageRegistrationMethodv4ShortDimensionsGTest.cxx b/Modules/Registration/RegistrationMethodsv4/test/itkImageRegistrationMethodv4ShortDimensionsGTest.cxx new file mode 100644 index 00000000000..022035a1795 --- /dev/null +++ b/Modules/Registration/RegistrationMethodsv4/test/itkImageRegistrationMethodv4ShortDimensionsGTest.cxx @@ -0,0 +1,85 @@ +/*========================================================================= + * + * 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 "itkImageRegistrationMethodv4.h" +#include "itkImageRegionIteratorWithIndex.h" +#include "itkMeanSquaresImageToImageMetricv4.h" +#include "itkRegularStepGradientDescentOptimizerv4.h" +#include "itkTranslationTransform.h" + +#include + +#include + +TEST(ImageRegistrationMethodv4ShortDimensions, SmoothSingletonImagesWithoutRecursiveGradientFilters) +{ + using ImageType = itk::Image; + using TransformType = itk::TranslationTransform; + using RegistrationType = itk::ImageRegistrationMethodv4; + for (unsigned int shortAxis = 0; shortAxis < 3; ++shortAxis) + { + SCOPED_TRACE(shortAxis); + auto size = ImageType::SizeType::Filled(16); + size[shortAxis] = 1; + auto image = ImageType::New(); + image->SetRegions(size); + image->Allocate(); + for (itk::ImageRegionIteratorWithIndex it(image, image->GetLargestPossibleRegion()); !it.IsAtEnd(); ++it) + { + const auto index = it.GetIndex(); + double squaredDistance = 0; + for (unsigned int d = 0; d < 3; ++d) + { + if (d != shortAxis) + { + const double distance = index[d] - 7.0; + squaredDistance += distance * distance; + } + } + it.Set(static_cast(std::exp(-squaredDistance / 16.0))); + } + + auto metric = itk::MeanSquaresImageToImageMetricv4::New(); + metric->SetUseFixedImageGradientFilter(false); + metric->SetUseMovingImageGradientFilter(false); + auto optimizer = itk::RegularStepGradientDescentOptimizerv4::New(); + optimizer->SetLearningRate(1.0); + optimizer->SetMinimumStepLength(0.001); + optimizer->SetNumberOfIterations(2); + auto registration = RegistrationType::New(); + registration->SetFixedImage(image); + registration->SetMovingImage(image); + registration->SetMetric(metric); + registration->SetOptimizer(optimizer); + registration->SetInitialTransform(TransformType::New()); + registration->SetNumberOfLevels(2); + RegistrationType::ShrinkFactorsArrayType shrinkFactors(2); + shrinkFactors[0] = 2; + shrinkFactors[1] = 1; + registration->SetShrinkFactorsPerLevel(shrinkFactors); + RegistrationType::SmoothingSigmasArrayType sigmas(2); + sigmas[0] = 1.0; + sigmas[1] = 0.5; + registration->SetSmoothingSigmasPerLevel(sigmas); + ASSERT_NO_THROW(registration->Update()); + for (const double parameter : registration->GetTransform()->GetParameters()) + { + EXPECT_NEAR(parameter, 0.0, 1e-8); + } + } +}