WIP: Add FFTHilbertImageFilter

include/rtkFFTHilbertImageFilter.h

@@ -0,0 +1,112 @@
+/*=========================================================================
+ *
+ *  Copyright RTK Consortium
+ *
+ *  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
+ *
+ *         http://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.
+ *
+ *=========================================================================*/
+
+#ifndef rtkFFTHilbertImageFilter_h
+#define rtkFFTHilbertImageFilter_h
+
+#include <itkConceptChecking.h>
+#include "rtkConfiguration.h"
+#include "rtkFFTProjectionsConvolutionImageFilter.h"
+#include "rtkMacro.h"
+
+namespace rtk
+{
+
+/** \class FFTHilbertImageFilter
+ * \brief Implements the Hilbert transform.
+ *
+ * \test rtkhilbertfiltertest.cxx
+ *
+ * \author Aurélien Coussat
+ *
+ * \ingroup RTK ImageToImageFilter
+ */
+
+template <class TInputImage, class TOutputImage = TInputImage, class TFFTPrecision = double>
+class ITK_EXPORT FFTHilbertImageFilter
+  : public rtk::FFTProjectionsConvolutionImageFilter<TInputImage, TOutputImage, TFFTPrecision>
+{
+public:
+  ITK_DISALLOW_COPY_AND_ASSIGN(FFTHilbertImageFilter);
+
+  /** Standard class type alias. */
+  using Self = FFTHilbertImageFilter;
+  using Superclass = rtk::FFTProjectionsConvolutionImageFilter<TInputImage, TOutputImage, TFFTPrecision>;
+  using Pointer = itk::SmartPointer<Self>;
+  using ConstPointer = itk::SmartPointer<const Self>;
+
+  /** Some convenient type alias. */
+  using InputImageType = TInputImage;
+  using OutputImageType = TOutputImage;
+  using FFTPrecisionType = TFFTPrecision;
+  using IndexType = typename InputImageType::IndexType;
+  using SizeType = typename InputImageType::SizeType;
+
+  using FFTInputImageType = typename Superclass::FFTInputImageType;
+  using FFTInputImagePointer = typename FFTInputImageType::Pointer;
+  using FFTOutputImageType = typename Superclass::FFTOutputImageType;
+  using FFTOutputImagePointer = typename FFTOutputImageType::Pointer;
+
+  /** Standard New method. */
+  itkNewMacro(Self);
+
+  itkGetMacro(PixelShift, double);
+  // The Set macro is redefined to clear the current FFT kernel when a parameter
+  // is modified.
+  virtual void
+  SetPixelShift(const double _arg)
+  {
+    itkDebugMacro("setting PixelShift to " << _arg);
+    CLANG_PRAGMA_PUSH
+    CLANG_SUPPRESS_Wfloat_equal if (this->m_PixelShift != _arg)
+    {
+      this->m_PixelShift = _arg;
+      this->Modified();
+      this->m_KernelFFT = nullptr;
+    }
+    CLANG_PRAGMA_POP
+  }
+
+  /** Runtime information support. */
+  itkTypeMacro(FFTHilbertImageFilter, FFTProjectionsConvolutionImageFilter);
+
+protected:
+  FFTHilbertImageFilter() = default;
+  ~FFTHilbertImageFilter() override = default;
+
+  void
+  GenerateOutputInformation() override;
+
+  /** Create and return a pointer to one line of the Hilbert kernel in Fourier space.
+   *  Used in generate data functions.  */
+  void
+  UpdateFFTProjectionsConvolutionKernel(const SizeType s) override;
+
+private:
+  SizeType m_PreviousKernelUpdateSize;
+  double   m_PixelShift;
+
+}; // end of class
+
+} // end namespace rtk
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#  include "rtkFFTHilbertImageFilter.hxx"
+#endif
+
+#endif

include/rtkFFTHilbertImageFilter.hxx

@@ -0,0 +1,101 @@
+/*=========================================================================
+ *
+ *  Copyright RTK Consortium
+ *
+ *  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
+ *
+ *         http://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.
+ *
+ *=========================================================================*/
+
+#ifndef rtkFFTHilbertImageFilter_hxx
+#define rtkFFTHilbertImageFilter_hxx
+
+#include "rtkFFTHilbertImageFilter.h"
+#include "itkForwardFFTImageFilter.h"
+
+namespace rtk
+{
+
+template <class TInputImage, class TOutputImage, class TFFTPrecision>
+void
+FFTHilbertImageFilter<TInputImage, TOutputImage, TFFTPrecision>::GenerateOutputInformation()
+{
+  FFTProjectionsConvolutionImageFilter<TInputImage, TOutputImage, TFFTPrecision>::GenerateOutputInformation();
+
+  auto   origin = this->GetOutput()->GetOrigin();
+  double spacing = this->GetOutput()->GetSpacing()[0];
+  origin[0] += m_PixelShift * spacing;
+  this->GetOutput()->SetOrigin(origin);
+}
+
+template <class TInputImage, class TOutputImage, class TFFTPrecision>
+void
+FFTHilbertImageFilter<TInputImage, TOutputImage, TFFTPrecision>::UpdateFFTProjectionsConvolutionKernel(const SizeType s)
+{
+  if (this->m_KernelFFT.GetPointer() != nullptr && s == this->m_PreviousKernelUpdateSize)
+  {
+    return;
+  }
+  m_PreviousKernelUpdateSize = s;
+
+  const int width = s[0];
+  const int height = s[1];
+
+  // Allocate kernel
+  SizeType size;
+  size.Fill(1);
+  size[0] = width;
+  FFTInputImagePointer kernel = FFTInputImageType::New();
+  kernel->SetRegions(size);
+  kernel->Allocate();
+  kernel->FillBuffer(0.);
+
+  itk::ImageRegionIterator<FFTInputImageType> it(kernel, kernel->GetLargestPossibleRegion());
+
+  // Compute band-limited kernel in space domain
+  double    spacing = this->GetInput()->GetSpacing()[0];
+  IndexType ix;
+
+  while (!it.IsAtEnd())
+  {
+    ix = it.GetIndex();
+
+    double x = (ix[0] + m_PixelShift) * spacing;
+    if (x > (size[0] / 2) * spacing)
+      x -= size[0] * spacing;
+
+    if (x == 0.)
+    {
+      it.Set(0.);
+    }
+    else
+    {
+      double v = spacing * (1. - std::cos(itk::Math::pi * x / spacing)) / (itk::Math::pi * x);
+      it.Set(v);
+    }
+
+    ++it;
+  }
+
+  // FFT kernel
+  using FFTType = itk::RealToHalfHermitianForwardFFTImageFilter<FFTInputImageType, FFTOutputImageType>;
+  typename FFTType::Pointer fftK = FFTType::New();
+  fftK->SetInput(kernel);
+  fftK->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
+  fftK->Update();
+
+  this->m_KernelFFT = fftK->GetOutput();
+  this->m_KernelFFT->DisconnectPipeline();
+}
+
+} // end namespace rtk
+#endif

wrapping/rtkFFTHilbertImageFilter.wrap

@@ -0,0 +1,12 @@
+if(RTK_USE_CUDA)
+  itk_wrap_include(itkCudaImage.h)
+endif()
+
+itk_wrap_class("rtk::FFTHilbertImageFilter"  POINTER)
+  foreach(t ${WRAP_ITK_REAL})
+    itk_wrap_template("I${ITKM_${t}}3I${ITKM_${t}}3${ITKM_D}" "itk::Image<${ITKT_${t}}, 3>, itk::Image<${ITKT_${t}}, 3>, ${ITKT_D}")
+  endforeach()
+  if(RTK_USE_CUDA)
+    itk_wrap_template("CI${ITKM_F}3CI${ITKM_F}3${ITKM_F}" "itk::CudaImage<${ITKT_F}, 3>, itk::CudaImage<${ITKT_F}, 3>, ${ITKT_F}")
+  endif()
+itk_end_wrap_class()