Merge branch 'feature-more-filters' into 'master'

Adds Gaussian type filters with general power exponent.

See merge request !110

cpp/kspace.cpp

@@ -345,7 +345,62 @@ void kspace<be, dt>::Gauss_filter(
         typename fftw_interface<rnumber>::complex *__restrict__ a,
         const double sigma)
 {
-    const double prefactor = - sigma*sigma/2;
+    this->template Gauss_n_filter<rnumber, fc, 2>(a, sigma);
+}
+
+template<int n>
+double hack_second_power(const double base_squared);
+
+template<>
+double hack_second_power<2>(const double base_squared)
+{
+    return base_squared;
+}
+template<>
+double hack_second_power<4>(const double base_squared)
+{
+    return base_squared*base_squared;
+}
+template<>
+double hack_second_power<6>(const double base_squared)
+{
+    return base_squared*base_squared*base_squared;
+}
+
+template<>
+double hack_second_power<8>(const double base_squared)
+{
+    const double base_to_fourth = base_squared*base_squared;
+    return base_to_fourth*base_to_fourth;
+}
+
+template<int n>
+double hack_second_power(const double base_squared)
+{
+     return std::pow(base_squared, n / 2.0);
+}
+
+/** \brief Filter a field using a Gaussian kernel with exponent to the n-th power.
+ *
+ * \tparam rnumber type of real number, float or double.
+ * \tparam fc field components, ONE, THREE or THREExTHREE.
+ * \return nothing
+ *
+ *  Filter's mathematical expression in Fourier space is as follows:
+ *  \f[
+ *      \hat{g}_\ell(\mathbf{k}) = \exp(-k^n \sigma^n / 2)
+ *  \f]
+ */
+template <field_backend be,
+          kspace_dealias_type dt>
+template <typename rnumber,
+          field_components fc,
+          int n>
+void kspace<be, dt>::Gauss_n_filter(
+        typename fftw_interface<rnumber>::complex *__restrict__ a,
+        const double sigma)
+{
+    const double prefactor = -pow(sigma, n)/2;
     this->CLOOP_K2(
             [&](const ptrdiff_t cindex,
                 const ptrdiff_t xindex,
@@ -354,11 +409,13 @@ void kspace<be, dt>::Gauss_filter(
                 const double k2){
                 {
                     for (unsigned int tcounter=0; tcounter<2*ncomp(fc); tcounter++)
-                        ((rnumber*)a)[2*ncomp(fc)*cindex + tcounter] *= std::exp(prefactor*k2);
+                        ((rnumber*)a)[2*ncomp(fc)*cindex + tcounter] *= std::exp(
+                            prefactor*hack_second_power<n>(k2));
                 }
                 });
 }
 
+
 /** \brief Filter a field.
  *
  *  This is a wrapper that can choose between a sharp Fourier spherical filter,
@@ -425,12 +482,34 @@ int kspace<be, dt>::filter(
                 a,
                 2*acos(0.)/wavenumber);
     }
+    else if (filter_type == std::string("Gauss4"))
+    {
+        this->template Gauss_n_filter<rnumber, fc, 4>(
+                a,
+                2*acos(0.)/wavenumber);
+    }
+    else if (filter_type == std::string("Gauss6"))
+    {
+        this->template Gauss_n_filter<rnumber, fc, 6>(
+                a,
+                2*acos(0.)/wavenumber);
+    }
+    else if (filter_type == std::string("Gauss8"))
+    {
+        this->template Gauss_n_filter<rnumber, fc, 8>(
+                a,
+                2*acos(0.)/wavenumber);
+    }
     else if (filter_type == std::string("ball"))
     {
         this->template ball_filter<rnumber, fc>(
                 a,
                 2*acos(0.)/wavenumber);
     }
+    else{
+        throw std::runtime_error(
+            "Filter type not available.");
+    }
     return EXIT_SUCCESS;
 }
 
@@ -513,6 +592,9 @@ int kspace<be, dt>::filter_calibrated_ell(
         this->template ball_filter<rnumber, fc>(
                 a,
                 ell);
+    }else{
+        throw std::runtime_error(
+            "Filter type not available.");
     }
     return EXIT_SUCCESS;
 }
@@ -958,23 +1040,103 @@ template void kspace<FFTW, SMOOTH>::low_pass<double, THREExTHREE>(
 
 template void kspace<FFTW, SMOOTH>::Gauss_filter<float, ONE>(
         typename fftw_interface<float>::complex *__restrict__ a,
-        const double kmax);
+        const double sigma);
 template void kspace<FFTW, SMOOTH>::Gauss_filter<float, THREE>(
         typename fftw_interface<float>::complex *__restrict__ a,
-        const double kmax);
+        const double sigma);
 template void kspace<FFTW, SMOOTH>::Gauss_filter<float, THREExTHREE>(
         typename fftw_interface<float>::complex *__restrict__ a,
-        const double kmax);
+        const double sigma);
 
 template void kspace<FFTW, SMOOTH>::Gauss_filter<double, ONE>(
         typename fftw_interface<double>::complex *__restrict__ a,
-        const double kmax);
+        const double sigma);
 template void kspace<FFTW, SMOOTH>::Gauss_filter<double, THREE>(
         typename fftw_interface<double>::complex *__restrict__ a,
-        const double kmax);
+        const double sigma);
 template void kspace<FFTW, SMOOTH>::Gauss_filter<double, THREExTHREE>(
         typename fftw_interface<double>::complex *__restrict__ a,
-        const double kmax);
+        const double sigma);
+
+template void kspace<FFTW, SMOOTH>::Gauss_n_filter<float, ONE, 4>(
+        typename fftw_interface<float>::complex *__restrict__ a,
+        const double sigma);
+template void kspace<FFTW, SMOOTH>::Gauss_n_filter<float, THREE, 4>(
+        typename fftw_interface<float>::complex *__restrict__ a,
+        const double sigma);
+template void kspace<FFTW, SMOOTH>::Gauss_n_filter<float, THREExTHREE, 4>(
+        typename fftw_interface<float>::complex *__restrict__ a,
+        const double sigma);
+
+template void kspace<FFTW, SMOOTH>::Gauss_n_filter<double, ONE, 4>(
+        typename fftw_interface<double>::complex *__restrict__ a,
+        const double sigma);
+template void kspace<FFTW, SMOOTH>::Gauss_n_filter<double, THREE, 4>(
+        typename fftw_interface<double>::complex *__restrict__ a,
+        const double sigma);
+template void kspace<FFTW, SMOOTH>::Gauss_n_filter<double, THREExTHREE, 4>(
+        typename fftw_interface<double>::complex *__restrict__ a,
+        const double sigma);
+
+template void kspace<FFTW, SMOOTH>::Gauss_n_filter<float, ONE, 6>(
+        typename fftw_interface<float>::complex *__restrict__ a,
+        const double sigma);
+template void kspace<FFTW, SMOOTH>::Gauss_n_filter<float, THREE, 6>(
+        typename fftw_interface<float>::complex *__restrict__ a,
+        const double sigma);
+template void kspace<FFTW, SMOOTH>::Gauss_n_filter<float, THREExTHREE, 6>(
+        typename fftw_interface<float>::complex *__restrict__ a,
+        const double sigma);
+
+template void kspace<FFTW, SMOOTH>::Gauss_n_filter<double, ONE, 6>(
+        typename fftw_interface<double>::complex *__restrict__ a,
+        const double sigma);
+template void kspace<FFTW, SMOOTH>::Gauss_n_filter<double, THREE, 6>(
+        typename fftw_interface<double>::complex *__restrict__ a,
+        const double sigma);
+template void kspace<FFTW, SMOOTH>::Gauss_n_filter<double, THREExTHREE, 6>(
+        typename fftw_interface<double>::complex *__restrict__ a,
+        const double sigma);
+
+template void kspace<FFTW, SMOOTH>::Gauss_n_filter<float, ONE, 8>(
+        typename fftw_interface<float>::complex *__restrict__ a,
+        const double sigma);
+template void kspace<FFTW, SMOOTH>::Gauss_n_filter<float, THREE, 8>(
+        typename fftw_interface<float>::complex *__restrict__ a,
+        const double sigma);
+template void kspace<FFTW, SMOOTH>::Gauss_n_filter<float, THREExTHREE, 8>(
+        typename fftw_interface<float>::complex *__restrict__ a,
+        const double sigma);
+
+template void kspace<FFTW, SMOOTH>::Gauss_n_filter<double, ONE, 8>(
+        typename fftw_interface<double>::complex *__restrict__ a,
+        const double sigma);
+template void kspace<FFTW, SMOOTH>::Gauss_n_filter<double, THREE, 8>(
+        typename fftw_interface<double>::complex *__restrict__ a,
+        const double sigma);
+template void kspace<FFTW, SMOOTH>::Gauss_n_filter<double, THREExTHREE, 8>(
+        typename fftw_interface<double>::complex *__restrict__ a,
+        const double sigma);
+
+template void kspace<FFTW, SMOOTH>::ball_filter<float, ONE>(
+        typename fftw_interface<float>::complex *__restrict__ a,
+        const double ell);
+template void kspace<FFTW, SMOOTH>::ball_filter<float, THREE>(
+        typename fftw_interface<float>::complex *__restrict__ a,
+        const double ell);
+template void kspace<FFTW, SMOOTH>::ball_filter<float, THREExTHREE>(
+        typename fftw_interface<float>::complex *__restrict__ a,
+        const double ell);
+
+template void kspace<FFTW, SMOOTH>::ball_filter<double, ONE>(
+        typename fftw_interface<double>::complex *__restrict__ a,
+        const double ell);
+template void kspace<FFTW, SMOOTH>::ball_filter<double, THREE>(
+        typename fftw_interface<double>::complex *__restrict__ a,
+        const double ell);
+template void kspace<FFTW, SMOOTH>::ball_filter<double, THREExTHREE>(
+        typename fftw_interface<double>::complex *__restrict__ a,
+        const double ell);
 
 template int kspace<FFTW, SMOOTH>::filter<float, ONE>(
         typename fftw_interface<float>::complex *__restrict__ a,

cpp/kspace.hpp

@@ -91,6 +91,13 @@ class kspace
                 typename fftw_interface<rnumber>::complex *__restrict__ a,
                 const double sigma);
 
+        template <typename rnumber,
+                  field_components fc,
+                  int n>
+        void Gauss_n_filter(
+                typename fftw_interface<rnumber>::complex *__restrict__ a,
+                const double sigma);
+
         template <typename rnumber,
                   field_components fc>
         void ball_filter(