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.
InsightSoftwareConsortium · Feb 8, 2025 · 1f5d0be · 1f5d0be
1 parent f28ff80
commit 1f5d0be
Show file tree

Hide file tree

Showing 5 changed files with 154 additions and 8 deletions.
diff --git 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<TInputPixel>::NonpositiveMin();
+
+    // check if scalar or fixed length array type
+    if constexpr (std::is_same<TInputPixel, typename NumericTraits<TInputPixel>::ValueType>::value)
+    {
+      m_Maximum = NumericTraits<TInputPixel>::NonpositiveMin();
+    }
+    else
+    {
+      m_Maximum = TInputPixel();
+      m_Maximum.Fill(NumericTraits<typename TInputPixel::ValueType>::NonpositiveMin());
+    }
   }
 
   inline void
   operator()(const TInputPixel & input)
   {
-    m_Maximum = std::max(m_Maximum, input);
+    if constexpr (std::is_same<TInputPixel, typename NumericTraits<TInputPixel>::ValueType>::value)
+    {
+      m_Maximum = std::max(m_Maximum, input);
+    }
+    else
+    {
+      if (itk::NumericTraits<TInputPixel>::GetLength(m_Maximum) == 0)
+      {
+        m_Maximum = input;
+      }
+      else
+      {
+        for (unsigned int i = 0; i < itk::NumericTraits<TInputPixel>::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<InputPixelType>));
+  itkConceptMacro(InputPixelTypeGreaterThanComparable,
+                  (Concept::GreaterThanComparable<typename itk::NumericTraits<InputPixelType>::ValueType>));
   itkConceptMacro(InputHasNumericTraitsCheck, (Concept::HasNumericTraits<InputPixelType>));
 
 protected:

diff --git 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<TInputPixel>::max();
+    if constexpr (std::is_same<TInputPixel, typename NumericTraits<TInputPixel>::ValueType>::value)
+    {
+      m_Minimum = NumericTraits<TInputPixel>::max();
+    }
+    else
+    {
+      m_Minimum = TInputPixel();
+      m_Minimum.Fill(NumericTraits<typename TInputPixel::ValueType>::max());
+    }
   }
 
   inline void
   operator()(const TInputPixel & input)
   {
-    m_Minimum = std::min(m_Minimum, input);
+
+    if constexpr (std::is_same<TInputPixel, typename NumericTraits<TInputPixel>::ValueType>::value)
+    {
+      m_Minimum = std::min(m_Minimum, input);
+    }
+    else
+    {
+      if (itk::NumericTraits<TInputPixel>::GetLength(m_Minimum) == 0)
+      {
+        m_Minimum = input;
+      }
+      else
+      {
+        for (unsigned int i = 0; i < itk::NumericTraits<TInputPixel>::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<InputPixelType>));
+  itkConceptMacro(InputPixelTypeGreaterThanComparable,
+                  (Concept::LessThanComparable<typename itk::NumericTraits<InputPixelType>::ValueType>));
   itkConceptMacro(InputHasNumericTraitsCheck, (Concept::HasNumericTraits<InputPixelType>));
 
 protected:

diff --git a/Modules/Filtering/ImageStatistics/include/itkProjectionImageFilter.hxx b/Modules/Filtering/ImageStatistics/include/itkProjectionImageFilter.hxx
@@ -120,6 +120,13 @@ ProjectionImageFilter<TInputImage, TOutputImage, TAccumulator>::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
@@ -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
@@ -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<PixelType, 3>;
+
+  using ReaderType = itk::ImageFileReader<ImageType>;
+  auto reader = ReaderType::New();
+  reader->SetFileName(argv[2]);
+
+  using FilterType = itk::MaximumProjectionImageFilter<ImageType, ImageType>;
+  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<ImageType>;
+  auto writer = WriterType::New();
+  writer->SetInput(filter->GetOutput());
+  writer->SetFileName(argv[3]);
+
+  ITK_TRY_EXPECT_NO_EXCEPTION(writer->Update());
+
+
+  return EXIT_SUCCESS;
+}