ENH: Add VectorImage support to MeanImageFilter

Add support for the runtime length vector images of
VectorImageFilter. The GenerateDataInSubregion method is overloaded to
specially handle VariableLengthVector Images. The specialization
avoids memory allocation multiple times per pixel  from temporary
variables, by pulling variable outside the loop.

When the change of variable scope was applied to scalars a 20%
performance decreases was observed with MSVC, so multiple methods are
used.

Modules/Filtering/Smoothing/include/itkMeanImageFilter.h

@@ -21,6 +21,7 @@
 #include "itkBoxImageFilter.h"
 #include "itkImage.h"
 #include "itkNumericTraits.h"
+#include "itkVariableLengthVector.h"
 
 #include <vector>
 
@@ -106,12 +107,21 @@ class ITK_TEMPLATE_EXPORT MeanImageFilter : public BoxImageFilter<TInputImage, T
   DynamicThreadedGenerateData(const OutputImageRegionType & outputRegionForThread) override;
 
 private:
-  template <typename TPixelAccessPolicy>
+  template <typename TPixelAccessPolicy, typename TPixelType>
   static void
   GenerateDataInSubregion(const TInputImage &                              inputImage,
                           TOutputImage &                                   outputImage,
                           const ImageRegion<InputImageDimension> &         imageRegion,
-                          const std::vector<Offset<InputImageDimension>> & neighborhoodOffsets);
+                          const std::vector<Offset<InputImageDimension>> & neighborhoodOffsets,
+                          const TPixelType *);
+
+  template <typename TPixelAccessPolicy, typename TValue>
+  static void
+  GenerateDataInSubregion(const TInputImage &                              inputImage,
+                          TOutputImage &                                   outputImage,
+                          const ImageRegion<InputImageDimension> &         imageRegion,
+                          const std::vector<Offset<InputImageDimension>> & neighborhoodOffsets,
+                          const VariableLengthVector<TValue> *);
 };
 } // end namespace itk
 

Modules/Filtering/Smoothing/include/itkMeanImageFilter.hxx

@@ -26,6 +26,7 @@
 #include "itkNeighborhoodAlgorithm.h"
 #include "itkOffset.h"
 #include "itkShapedImageNeighborhoodRange.h"
+#include "itkDefaultConvertPixelTraits.h"
 
 namespace itk
 {
@@ -54,26 +55,31 @@ MeanImageFilter<TInputImage, TOutputImage>::DynamicThreadedGenerateData(
 
   // Process the non-boundary subregion, using a faster pixel access policy without boundary extrapolation.
   GenerateDataInSubregion<BufferedImageNeighborhoodPixelAccessPolicy<InputImageType>>(
-    *input, *output, calculatorResult.GetNonBoundaryRegion(), neighborhoodOffsets);
+    *input,
+    *output,
+    calculatorResult.GetNonBoundaryRegion(),
+    neighborhoodOffsets,
+    static_cast<InputPixelType *>(nullptr));
 
   // Process each of the boundary faces. These are N-d regions which border
   // the edge of the buffer.
   for (const auto & boundaryFace : calculatorResult.GetBoundaryFaces())
   {
     GenerateDataInSubregion<ZeroFluxNeumannImageNeighborhoodPixelAccessPolicy<InputImageType>>(
-      *input, *output, boundaryFace, neighborhoodOffsets);
+      *input, *output, boundaryFace, neighborhoodOffsets, static_cast<InputPixelType *>(nullptr));
   }
 }
 
 
 template <typename TInputImage, typename TOutputImage>
-template <typename TPixelAccessPolicy>
+template <typename TPixelAccessPolicy, typename TPixelType>
 void
 MeanImageFilter<TInputImage, TOutputImage>::GenerateDataInSubregion(
   const TInputImage &                              inputImage,
   TOutputImage &                                   outputImage,
   const ImageRegion<InputImageDimension> &         imageRegion,
-  const std::vector<Offset<InputImageDimension>> & neighborhoodOffsets)
+  const std::vector<Offset<InputImageDimension>> & neighborhoodOffsets,
+  const TPixelType *)
 {
   const auto neighborhoodSize = static_cast<double>(neighborhoodOffsets.size());
 
@@ -98,6 +104,57 @@ MeanImageFilter<TInputImage, TOutputImage>::GenerateDataInSubregion(
     ++outputIterator;
   }
 }
+
+template <typename TInputImage, typename TOutputImage>
+template <typename TPixelAccessPolicy, typename TValueType>
+void
+MeanImageFilter<TInputImage, TOutputImage>::GenerateDataInSubregion(
+  const TInputImage &                              inputImage,
+  TOutputImage &                                   outputImage,
+  const ImageRegion<InputImageDimension> &         imageRegion,
+  const std::vector<Offset<InputImageDimension>> & neighborhoodOffsets,
+  const VariableLengthVector<TValueType> *)
+{
+  const auto neighborhoodSize = static_cast<double>(neighborhoodOffsets.size());
+
+  using namespace Experimental;
+  auto neighborhoodRange = ShapedImageNeighborhoodRange<const InputImageType, TPixelAccessPolicy>(
+    inputImage, Index<InputImageDimension>(), neighborhoodOffsets);
+  auto outputIterator = ImageRegionRange<OutputImageType>(outputImage, imageRegion).begin();
+
+  // These temp variable are needed outside the loop for
+  // VariableLengthVectors to avoid memory allocations on a per-pixel
+  // basis.
+  InputRealType   sum(inputImage.GetNumberOfComponentsPerPixel());
+  OutputPixelType out(inputImage.GetNumberOfComponentsPerPixel());
+
+  for (const auto & index : ImageRegionIndexRange<InputImageDimension>(imageRegion))
+  {
+    neighborhoodRange.SetLocation(index);
+
+    using PixelComponentType = typename NumericTraits<InputRealType>::ValueType;
+    sum.Fill(NumericTraits<PixelComponentType>::Zero);
+
+    for (const InputPixelType pixelValue : neighborhoodRange)
+    {
+      sum += pixelValue;
+    }
+
+    sum /= neighborhoodSize;
+
+    // The following line is an implicit conversion from the
+    // InputRealType to the output pixel. If an explicit construction
+    // or static_cast was used a new VariableLengthVector temporary would be
+    // constructed requiring the array to be dynamic allocated. This
+    // implicit assignment reuses the array allocated in the variable out.
+    // *DO NOT USE static_cast*
+    out = sum;
+    *outputIterator = out;
+    ++outputIterator;
+  }
+}
+
+
 } // end namespace itk
 
 #endif

Modules/Filtering/Smoothing/test/itkMeanImageFilterGTest.cxx

@@ -20,6 +20,7 @@
 #include "itkMeanImageFilter.h"
 
 #include "itkImage.h"
+#include "itkVectorImage.h"
 #include "itkImageBufferRange.h"
 
 #include <numeric> // For iota.
@@ -39,6 +40,7 @@ Expect_output_pixels_have_same_value_as_input_when_input_image_is_uniform(
 
   const auto image = TImage::New();
   image->SetRegions(imageRegion);
+  image->SetNumberOfComponentsPerPixel(itk::NumericTraits<PixelType>::GetLength(inputPixelValue));
   image->Allocate();
   image->FillBuffer(inputPixelValue);
 
@@ -102,6 +104,11 @@ TEST(MeanImageFilter, OutputSameAsInputForUniformImage)
                                                                                              64);
   Expect_output_pixels_have_same_value_as_input_when_input_image_is_uniform<itk::Image<float, 3>>(
     itk::Size<3>{ { 3, 4, 5 } }, 0.5f);
+
+  float                            array[] = { 3.14f, 2.71f, 1.41f };
+  itk::VariableLengthVector<float> v{ array, 3 };
+  Expect_output_pixels_have_same_value_as_input_when_input_image_is_uniform<itk::VectorImage<float, 3>>(
+    itk::Size<3>{ { 3, 4, 5 } }, v);
 }
 
 

Modules/Filtering/Smoothing/test/itkMeanImageFilterTest.cxx

@@ -20,10 +20,14 @@
 #include "itkMeanImageFilter.h"
 #include "itkTextOutput.h"
 #include "itkTestingMacros.h"
+#include "itkVectorImage.h"
 
 int
 itkMeanImageFilterTest(int, char *[])
 {
+  using VectorImageType = itk::VectorImage<float, 3>;
+  itk::MeanImageFilter<VectorImageType, VectorImageType>::New();
+
   // Comment the following if you want to use the itk text output window
   itk::OutputWindow::SetInstance(itk::TextOutput::New());