From 1f5d0be02c72ac438d311d87e39d3660a8f80273 Mon Sep 17 00:00:00 2001 From: Bradley Lowekamp Date: Fri, 7 Feb 2025 13:27:11 -0500 Subject: [PATCH] ENH: Add support for VectorImages to Projection filters Update the ProjectionImageFilter base class to propagate the number of components to the output image. Add multi-component support to max and min projection filters by performing operation on per commonents with initialization support for VariableLengthVectors. --- .../include/itkMaximumProjectionImageFilter.h | 34 ++++++++- .../include/itkMinimumProjectionImageFilter.h | 33 ++++++++- .../include/itkProjectionImageFilter.hxx | 7 ++ .../ImageStatistics/test/CMakeLists.txt | 19 ++++- .../itkMaximumProjectionImageFilterTest4.cxx | 69 +++++++++++++++++++ 5 files changed, 154 insertions(+), 8 deletions(-) create mode 100644 Modules/Filtering/ImageStatistics/test/itkMaximumProjectionImageFilterTest4.cxx diff --git a/Modules/Filtering/ImageStatistics/include/itkMaximumProjectionImageFilter.h b/Modules/Filtering/ImageStatistics/include/itkMaximumProjectionImageFilter.h index 807c9b18d77..d852db7d859 100644 --- a/Modules/Filtering/ImageStatistics/include/itkMaximumProjectionImageFilter.h +++ b/Modules/Filtering/ImageStatistics/include/itkMaximumProjectionImageFilter.h @@ -55,13 +55,40 @@ class MaximumAccumulator inline void Initialize() { - m_Maximum = NumericTraits::NonpositiveMin(); + + // check if scalar or fixed length array type + if constexpr (std::is_same::ValueType>::value) + { + m_Maximum = NumericTraits::NonpositiveMin(); + } + else + { + m_Maximum = TInputPixel(); + m_Maximum.Fill(NumericTraits::NonpositiveMin()); + } } inline void operator()(const TInputPixel & input) { - m_Maximum = std::max(m_Maximum, input); + if constexpr (std::is_same::ValueType>::value) + { + m_Maximum = std::max(m_Maximum, input); + } + else + { + if (itk::NumericTraits::GetLength(m_Maximum) == 0) + { + m_Maximum = input; + } + else + { + for (unsigned int i = 0; i < itk::NumericTraits::GetLength(m_Maximum); ++i) + { + m_Maximum[i] = std::max(m_Maximum[i], input[i]); + } + } + } } inline TInputPixel @@ -99,7 +126,8 @@ class MaximumProjectionImageFilter /** Method for creation through the object factory. */ itkNewMacro(Self); - itkConceptMacro(InputPixelTypeGreaterThanComparable, (Concept::GreaterThanComparable)); + itkConceptMacro(InputPixelTypeGreaterThanComparable, + (Concept::GreaterThanComparable::ValueType>)); itkConceptMacro(InputHasNumericTraitsCheck, (Concept::HasNumericTraits)); protected: diff --git a/Modules/Filtering/ImageStatistics/include/itkMinimumProjectionImageFilter.h b/Modules/Filtering/ImageStatistics/include/itkMinimumProjectionImageFilter.h index 1d374d1dcdd..3cb0425aea1 100644 --- a/Modules/Filtering/ImageStatistics/include/itkMinimumProjectionImageFilter.h +++ b/Modules/Filtering/ImageStatistics/include/itkMinimumProjectionImageFilter.h @@ -54,13 +54,39 @@ class MinimumAccumulator inline void Initialize() { - m_Minimum = NumericTraits::max(); + if constexpr (std::is_same::ValueType>::value) + { + m_Minimum = NumericTraits::max(); + } + else + { + m_Minimum = TInputPixel(); + m_Minimum.Fill(NumericTraits::max()); + } } inline void operator()(const TInputPixel & input) { - m_Minimum = std::min(m_Minimum, input); + + if constexpr (std::is_same::ValueType>::value) + { + m_Minimum = std::min(m_Minimum, input); + } + else + { + if (itk::NumericTraits::GetLength(m_Minimum) == 0) + { + m_Minimum = input; + } + else + { + for (unsigned int i = 0; i < itk::NumericTraits::GetLength(m_Minimum); ++i) + { + m_Minimum[i] = std::min(m_Minimum[i], input[i]); + } + } + } } inline TInputPixel @@ -98,7 +124,8 @@ class MinimumProjectionImageFilter /** Method for creation through the object factory. */ itkNewMacro(Self); - itkConceptMacro(InputPixelTypeGreaterThanComparable, (Concept::LessThanComparable)); + itkConceptMacro(InputPixelTypeGreaterThanComparable, + (Concept::LessThanComparable::ValueType>)); itkConceptMacro(InputHasNumericTraitsCheck, (Concept::HasNumericTraits)); protected: diff --git a/Modules/Filtering/ImageStatistics/include/itkProjectionImageFilter.hxx b/Modules/Filtering/ImageStatistics/include/itkProjectionImageFilter.hxx index 3ecfb0eb0cd..4185b8916dc 100644 --- a/Modules/Filtering/ImageStatistics/include/itkProjectionImageFilter.hxx +++ b/Modules/Filtering/ImageStatistics/include/itkProjectionImageFilter.hxx @@ -120,6 +120,13 @@ ProjectionImageFilter::GenerateOutputIn output->SetDirection(outDirection); output->SetLargestPossibleRegion(outputRegion); + // Support VectorImages by setting number of components on output. + const unsigned int numComponents = input->GetNumberOfComponentsPerPixel(); + if (numComponents != output->GetNumberOfComponentsPerPixel()) + { + output->SetNumberOfComponentsPerPixel(numComponents); + } + itkDebugMacro("GenerateOutputInformation End"); } diff --git a/Modules/Filtering/ImageStatistics/test/CMakeLists.txt b/Modules/Filtering/ImageStatistics/test/CMakeLists.txt index 37058a70a48..3aaa2b6bfe6 100644 --- a/Modules/Filtering/ImageStatistics/test/CMakeLists.txt +++ b/Modules/Filtering/ImageStatistics/test/CMakeLists.txt @@ -7,12 +7,13 @@ set(ITKImageStatisticsTests itkImageMomentsTest.cxx itkMinimumMaximumImageFilterTest.cxx itkImagePCAShapeModelEstimatorTest.cxx - itkMaximumProjectionImageFilterTest2.cxx - itkMaximumProjectionImageFilterTest3.cxx itkMinimumProjectionImageFilterTest.cxx itkMeanProjectionImageFilterTest.cxx itkImagePCADecompositionCalculatorTest.cxx itkMaximumProjectionImageFilterTest.cxx + itkMaximumProjectionImageFilterTest2.cxx + itkMaximumProjectionImageFilterTest3.cxx + itkMaximumProjectionImageFilterTest4.cxx itkMedianProjectionImageFilterTest.cxx itkHistogramToEntropyImageFilterTest1.cxx itkHistogramToEntropyImageFilterTest2.cxx @@ -177,6 +178,20 @@ itk_add_test( 1 DATA{${ITK_DATA_ROOT}/Input/HeadMRVolumeWithDirection.mhd,HeadMRVolume.raw} ${ITK_TEST_OUTPUT_DIR}/HeadMRVolumeMaximumProjection2D1.tif) +itk_add_test( + NAME + itkMaximumProjectionImageFilterTest4_1 + COMMAND + ITKImageStatisticsTestDriver + --compare-MD5 + ${ITK_TEST_OUTPUT_DIR}/VHFColorMaximumProjection2D1.tif + c0674d2b9a01cb2d6fae981c32c95877 + itkMaximumProjectionImageFilterTest4 + 2 + DATA{${ITK_DATA_ROOT}/Input/VHFColor.mhd,VHFColor.mhd.raw} + ${ITK_TEST_OUTPUT_DIR}/VHFColorMaximumProjection2D1.tif +) + itk_add_test( NAME itkMinimumProjectionImageFilterTest diff --git a/Modules/Filtering/ImageStatistics/test/itkMaximumProjectionImageFilterTest4.cxx b/Modules/Filtering/ImageStatistics/test/itkMaximumProjectionImageFilterTest4.cxx new file mode 100644 index 00000000000..c3dade34da3 --- /dev/null +++ b/Modules/Filtering/ImageStatistics/test/itkMaximumProjectionImageFilterTest4.cxx @@ -0,0 +1,69 @@ +/*========================================================================= + * + * 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 "itkImageFileReader.h" +#include "itkImageFileWriter.h" +#include "itkMetaImageIO.h" +#include "itkSimpleFilterWatcher.h" + +#include "itkVectorImage.h" +#include "itkMaximumProjectionImageFilter.h" +#include "itkTestingMacros.h" + +int +itkMaximumProjectionImageFilterTest4(int argc, char * argv[]) +{ + if (argc < 4) + { + std::cerr << "Missing parameters." << std::endl; + std::cerr << "Usage: " << itkNameOfTestExecutableMacro(argv); + std::cerr << "Dimension Inputimage Outputimage " << std::endl; + return EXIT_FAILURE; + } + + // Legacy compat with older MetaImages + itk::MetaImageIO::SetDefaultDoublePrecision(6); + const int dim = std::stoi(argv[1]); + + using PixelType = unsigned char; + + using ImageType = itk::VectorImage; + + using ReaderType = itk::ImageFileReader; + auto reader = ReaderType::New(); + reader->SetFileName(argv[2]); + + using FilterType = itk::MaximumProjectionImageFilter; + auto filter = FilterType::New(); + filter->SetInput(reader->GetOutput()); + filter->SetProjectionDimension(dim); + // to be sure that the result is ok with several threads, even on a single + // proc computer + filter->SetNumberOfWorkUnits(2); + + const itk::SimpleFilterWatcher watcher(filter, "filter"); + + using WriterType = itk::ImageFileWriter; + auto writer = WriterType::New(); + writer->SetInput(filter->GetOutput()); + writer->SetFileName(argv[3]); + + ITK_TRY_EXPECT_NO_EXCEPTION(writer->Update()); + + + return EXIT_SUCCESS; +}