ENH: Initial itkHalideDiscreteGaussianImageFilter (#10)

* WIP: ENH: Add trivial halide cast image filter

* WIP: FIX: Use GenerateData() rather than DynamicThreadedGenerateData()

* WIP: Simple itkHalideDiscreteGaussianImpl generator stub; integrated and working.

* WIP: ENH: Implement HalideDiscreteGaussian with naive schedule.

Implemented with separable convolution. Handles separate sigma and kernel size for each dimension to support anisotropic images.

Strip out the CastImageFilter and MinimalStandardRandomVariateGenerator placeholders.

* WIP: ENH: Add naive CUDA schedule

* ENH: Improved CUDA schedule.

About 200ms on my machine. Bandwidth-limited on reasonable kernel/image sizes.

* ENH: Generate kernels with itk::GaussianOperator

for parity with itk::DiscreteGaussianImageFilter.

* FIX: Disable multithreading.

It wasn't used anyway, but this should make the filter better-behaved.

Not sure how to disable/remove the MultiThreader object; setting it null is not valid.

* FIX: Verification Test

Change the test to blur CTChest. The reference output was generated with itk::DiscreteGaussianImageFilter directly, so the Halide output is verified to be numerically compatible.

---------

Co-authored-by: Anthony Lombardi <[email protected]>

Author: allemangD

CMakeLists.txt

@@ -65,7 +65,7 @@ FetchContent_MakeAvailable(halide)
 set(Halide_DIR "${halide_SOURCE_DIR}/lib/cmake/Halide")
 set(HalideHelpers_DIR "${halide_SOURCE_DIR}/lib/cmake/HalideHelpers")
 
-set(HalideFilters_LIBRARIES HalideFilters my_first_generator)
+set(HalideFilters_LIBRARIES HalideFilters)
 
 if(NOT ITK_SOURCE_DIR)
   find_package(ITK REQUIRED)

examples/SampleHalidePipeline.py

@@ -1,42 +1,34 @@
 #!/usr/bin/env python
 
-# 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.
-
-# Run with:
-# ./SampleHalidePipeline.py <input_image> <output_image>
-# <parameters>
-# e.g.
-# ./SampleHalidePipeline.py MyImage.mha Output.mha 2 0.2
-# (A rule of thumb is to set the Threshold to be about 1 / 100 of the Level.)
-#
-#  parameter_1: absolute minimum...
-#        The assumption is that...
-#  parameter_2: controls the..
-#        A tradeoff between...
-
 import argparse
+from time import perf_counter
 
 import itk
 
 
 parser = argparse.ArgumentParser(description="Example short description.")
 parser.add_argument("input_image")
+parser.add_argument("variance", type=float)
 parser.add_argument("output_image")
-parser.add_argument("parameter_1")
-args = parser.parse_args()
 
-# Please, write a complete, self-containted and useful example that
-# demonstrate a class when being used along with other ITK classes or in
-# the context of a wider or specific application.
+ImageType = itk.Image[itk.F, 3]
+FilterType = itk.HalideDiscreteGaussianImageFilter[ImageType, ImageType]
+
+if __name__ == '__main__':
+    args = parser.parse_args()
+
+    im = itk.imread(args.input_image, itk.F)
+
+    proc = FilterType.New()
+    proc.SetInput(im)
+    proc.SetVariance(args.variance)
+
+    start = perf_counter()
+    proc.Update()
+    end = perf_counter()
+
+    delta = end - start
+
+    print(f'took {delta * 1_000:.3f}ms')
+
+    itk.imwrite(proc.GetOutput(), args.output_image)

include/itkHalideDiscreteGaussianImageFilter.h

@@ -32,6 +32,10 @@ namespace itk
  *
  * \ingroup HalideFilters
  *
+ * Limitations compared te itkDiscreteGaussianImageFilter:
+ * - Only supports isotropic variance and maximum error (to simplify wrapper)
+ * - Only supports 3d images (to simplify wrapper)
+ *
  */
 template <typename TInputImage, typename TOutputImage>
 class HalideDiscreteGaussianImageFilter : public ImageToImageFilter<TInputImage, TOutputImage>
@@ -59,23 +63,42 @@ class HalideDiscreteGaussianImageFilter : public ImageToImageFilter<TInputImage,
   /** Standard New macro. */
   itkNewMacro(Self);
 
+  itkSetMacro(Variance, float);
+  itkGetMacro(Variance, float);
+
+  itkSetMacro(MaximumError, float);
+  itkGetMacro(MaximumError, float);
+
+  itkGetMacro(MaximumKernelWidth, unsigned int);
+  itkSetMacro(MaximumKernelWidth, unsigned int);
+
+  itkGetMacro(UseImageSpacing, bool);
+  itkSetMacro(UseImageSpacing, bool);
+  itkBooleanMacro(UseImageSpacing);
+
 protected:
   HalideDiscreteGaussianImageFilter();
-  ~HalideDiscreteGaussianImageFilter() override = default;
+  ~
+  HalideDiscreteGaussianImageFilter() override = default;
 
   void
   PrintSelf(std::ostream & os, Indent indent) const override;
 
   using OutputRegionType = typename OutputImageType::RegionType;
 
   void
-  DynamicThreadedGenerateData(const OutputRegionType & outputRegion) override;
+  GenerateData() override;
 
 private:
 #ifdef ITK_USE_CONCEPT_CHECKING
   // Add concept checking such as
-  // itkConceptMacro( FloatingPointPixel, ( itk::Concept::IsFloatingPoint< typename InputImageType::PixelType > ) );
+  itkConceptMacro(FloatingPointPixel, (itk::Concept::IsFloatingPoint<typename InputImageType::PixelType>));
 #endif
+
+  float        m_Variance = 0;
+  float        m_MaximumError = 0.01;
+  unsigned int m_MaximumKernelWidth = 32;
+  bool         m_UseImageSpacing = true;
 };
 } // namespace itk
 

include/itkHalideDiscreteGaussianImageFilter.hxx

@@ -20,44 +20,79 @@
 
 #include "itkHalideDiscreteGaussianImageFilter.h"
 
-#include "itkImageRegionIterator.h"
-#include "itkImageRegionConstIterator.h"
+#include "itkHalideSeparableConvolutionImpl.h"
+
+#include "itkGaussianOperator.h"
+
+#include <Halide.h>
+#include <HalideBuffer.h>
+#include <iomanip>
 
 namespace itk
 {
 
 template <typename TInputImage, typename TOutputImage>
-HalideDiscreteGaussianImageFilter<TInputImage, TOutputImage>
-::HalideDiscreteGaussianImageFilter()
-{}
+HalideDiscreteGaussianImageFilter<TInputImage, TOutputImage>::HalideDiscreteGaussianImageFilter()
+{
+  this->DynamicMultiThreadingOff();
+  this->ThreaderUpdateProgressOff();
+}
 
 
 template <typename TInputImage, typename TOutputImage>
 void
-HalideDiscreteGaussianImageFilter<TInputImage, TOutputImage>
-::PrintSelf(std::ostream & os, Indent indent) const
+HalideDiscreteGaussianImageFilter<TInputImage, TOutputImage>::PrintSelf(std::ostream & os, Indent indent) const
 {
   Superclass::PrintSelf(os, indent);
 }
 
 
 template <typename TInputImage, typename TOutputImage>
 void
-HalideDiscreteGaussianImageFilter<TInputImage, TOutputImage>
-::DynamicThreadedGenerateData(const OutputRegionType & outputRegion)
+HalideDiscreteGaussianImageFilter<TInputImage, TOutputImage>::GenerateData()
 {
-  OutputImageType *      output = this->GetOutput();
-  const InputImageType * input = this->GetInput();
-  using InputRegionType = typename InputImageType::RegionType;
-  InputRegionType inputRegion = InputRegionType(outputRegion.GetSize());
+  const InputImageType *               input = this->GetInput();
+  typename InputImageType::RegionType  inputRegion = input->GetBufferedRegion();
+  typename InputImageType::SizeType    inputSize = inputRegion.GetSize();
+  typename InputImageType::SpacingType inputSpacing = input->GetSpacing();
 
-  itk::ImageRegionConstIterator<InputImageType> in(input, inputRegion);
-  itk::ImageRegionIterator<OutputImageType>     out(output, outputRegion);
+  std::vector<Halide::Runtime::Buffer<float, 1>> kernel_buffers{};
 
-  for (in.GoToBegin(), out.GoToBegin(); !in.IsAtEnd() && !out.IsAtEnd(); ++in, ++out)
+  // compute kernel coefficients with itk::GaussianOperator to match behavior with itk::DiscreteGaussianImageFilter
+  for (int dim = 0; dim < InputImageDimension; ++dim)
   {
-    out.Set(in.Get());
+    GaussianOperator<float, 1> oper{};
+    oper.SetMaximumError(m_MaximumError);
+    oper.SetMaximumKernelWidth(m_MaximumKernelWidth);
+
+    float variance = m_Variance;
+    if (m_UseImageSpacing)
+    {
+      variance /= inputSpacing[dim];
+    }
+    oper.SetVariance(variance);
+
+    oper.CreateDirectional();
+
+    Halide::Runtime::Buffer<float, 1> & buf = kernel_buffers.emplace_back(static_cast<int>(oper.GetSize(0)));
+    buf.set_min(-static_cast<int>(oper.GetRadius(0)));
+    std::copy(oper.Begin(), oper.End(), buf.begin());
+    buf.set_host_dirty();
   }
+
+  OutputImageType * output = this->GetOutput();
+  output->SetRegions(inputRegion);
+  output->Allocate();
+
+  std::vector<int> sizes(3, 1);
+  std::copy(inputSize.begin(), inputSize.end(), sizes.begin());
+
+  Halide::Runtime::Buffer<const InputPixelType> inputBuffer(input->GetBufferPointer(), sizes);
+  Halide::Runtime::Buffer<OutputPixelType>      outputBuffer(output->GetBufferPointer(), sizes);
+
+  inputBuffer.set_host_dirty();
+  itkHalideSeparableConvolutionImpl(inputBuffer, kernel_buffers[0], kernel_buffers[1], kernel_buffers[2], outputBuffer);
+  outputBuffer.copy_to_host();
 }
 
 } // end namespace itk

include/itkMinimalStandardRandomVariateGenerator.h

@@ -1,13 +1,18 @@
 find_package(Halide REQUIRED shared)
 
-add_executable(my_generators generators.cpp)
-target_link_libraries(my_generators PRIVATE Halide::Generator)
-add_halide_library(my_first_generator FROM my_generators HEADER my_first_generator_h)
+add_executable(itkHalideGenerators generators.cpp)
+target_link_libraries(itkHalideGenerators PRIVATE Halide::Generator)
+
+add_halide_library(itkHalideSeparableConvolutionImpl
+  FROM itkHalideGenerators
+  HEADER itkHalideSeparableConvolutionImpl_h
+  FEATURES cuda
+  )
 
 set(HalideFilters_SRCS
-  itkMinimalStandardRandomVariateGenerator.cxx
-  ${my_first_generator_h}
+  ${itkHalideSeparableConvolutionImpl_h}
   )
 
 itk_module_add_library(HalideFilters ${HalideFilters_SRCS})
 target_include_directories(HalideFilters PRIVATE ${CMAKE_CURRENT_BINARY_DIR})
+target_link_libraries(HalideFilters PUBLIC itkHalideSeparableConvolutionImpl)