From e3c411403427ac5b485a87feb6bdebcb84db1b6c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Aur=C3=A9lien=20Coussat?=
 <aurelien.coussat@creatis.insa-lyon.fr>
Date: Tue, 24 Nov 2020 15:45:18 +0100
Subject: [PATCH] ENH: test for FFT Hilbert filter

---
 test/CMakeLists.txt           |  2 +
 test/rtkhilbertfiltertest.cxx | 87 +++++++++++++++++++++++++++++++++++
 2 files changed, 89 insertions(+)
 create mode 100644 test/rtkhilbertfiltertest.cxx

diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 197562c55..8a8d8f5e5 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -77,6 +77,8 @@ rtk_add_cuda_test(rtkRampFilterCudaTest rtkrampfiltertest.cxx)
 rtk_add_test(rtkRampFilterTest2 rtkrampfiltertest2.cxx)
 rtk_add_cuda_test(rtkRampFilterCudaTest2 rtkrampfiltertest2.cxx)
 
+rtk_add_test(rtkHilbertFilterTest rtkhilbertfiltertest.cxx)
+
 rtk_add_test(rtkScatterGlareFilterTest rtkscatterglarefiltertest.cxx)
 rtk_add_cuda_test(rtkScatterGlareFilterCudaTest rtkscatterglarefiltertest.cxx)
 
diff --git a/test/rtkhilbertfiltertest.cxx b/test/rtkhilbertfiltertest.cxx
new file mode 100644
index 000000000..888a1983f
--- /dev/null
+++ b/test/rtkhilbertfiltertest.cxx
@@ -0,0 +1,87 @@
+#include "rtkTest.h"
+
+#include "rtkFFTHilbertImageFilter.h"
+/**
+ * \file rtkhilbertfiltertest.cxx
+ *
+ * \brief Functional test for the FFT Hilbert filter.
+ *
+ * This test computes the FFT Hilbert transform of the signal t -> sin(20*pi*t).
+ * The computed Hilbert transform is compared against the known analytic Hilbert transform of the signal.
+ *
+ * \author Aurélien Coussat
+ */
+
+int
+main(int, char **)
+{
+  constexpr unsigned int Dimension = 3;
+  using OutputPixelType = float;
+#ifdef USE_CUDA
+  using OutputImageType = itk::CudaImage<OutputPixelType, Dimension>;
+#else
+  using OutputImageType = itk::Image<OutputPixelType, Dimension>;
+#endif
+
+  OutputImageType::SizeType size;
+  size[0] = 500;
+  size[1] = size[2] = 1;
+
+  OutputImageType::SpacingType spacing;
+  spacing[0] = spacing[1] = spacing[2] = .001;
+
+  OutputImageType::IndexType ix;
+
+  // Build the signal t -> sin(20*pi*t)
+  OutputImageType::Pointer signal = OutputImageType::New();
+  signal->SetSpacing(spacing);
+  signal->SetRegions(size);
+  signal->Allocate();
+
+  itk::ImageRegionIterator<OutputImageType> signalIt(signal, signal->GetLargestPossibleRegion());
+
+  while (!signalIt.IsAtEnd())
+  {
+    ix = signalIt.GetIndex();
+    double t = ix[0] * spacing[0];
+    double v = sin(20 * itk::Math::pi * t);
+    signalIt.Set(v);
+    ++signalIt;
+  }
+
+  // Build the analytic Hilbert transform of the signal
+  // It is t -> -cos(20*pi*t)
+  OutputImageType::Pointer analyticHilbertSignal = OutputImageType::New();
+  analyticHilbertSignal->SetSpacing(spacing);
+  analyticHilbertSignal->SetRegions(size);
+  analyticHilbertSignal->Allocate();
+
+  itk::ImageRegionIterator<OutputImageType> analyticHilbertSignalIt(analyticHilbertSignal,
+                                                                    analyticHilbertSignal->GetLargestPossibleRegion());
+
+  while (!analyticHilbertSignalIt.IsAtEnd())
+  {
+    ix = analyticHilbertSignalIt.GetIndex();
+    double t = ix[0] * spacing[0];
+    double v = -cos(20 * itk::Math::pi * t);
+    analyticHilbertSignalIt.Set(v);
+    ++analyticHilbertSignalIt;
+  }
+
+  // Compute the Hilbert transform of the signal using rtkFFTHilbertImageFilter
+  using HilbertType = rtk::FFTHilbertImageFilter<OutputImageType>;
+  HilbertType::Pointer hilbert = HilbertType::New();
+  hilbert->SetInput(signal);
+
+  using ZeroPadFactorsType = HilbertType::ZeroPadFactorsType;
+  ZeroPadFactorsType zeroPadFactors;
+  zeroPadFactors.Fill(1);
+  hilbert->SetZeroPadFactors(zeroPadFactors);
+
+  TRY_AND_EXIT_ON_ITK_EXCEPTION(hilbert->Update());
+
+  CheckImageQuality<OutputImageType>(hilbert->GetOutput(), analyticHilbertSignal, .1, 1., 1.);
+
+  std::cout << "\n\nTest PASSED! " << std::endl;
+  return EXIT_SUCCESS;
+}