alternative approach

nifty_gridder.cc

@@ -12,7 +12,7 @@
  *  GNU General Public License for more details.
  *
  *  You should have received a copy of the GNU General Public License
- *  along with nifty_fridder; if not, write to the Free Software
+ *  along with nifty_gridder; if not, write to the Free Software
  *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  */
 
@@ -310,31 +310,6 @@ vector<double> correction_factors (size_t n, size_t nval, size_t w)
   return res;
   }
 
-constexpr auto get_taper_DS = R"""(
-Returns the correction coefficients.
-
-Parameters
-==========
-ndirty : integer
-    the number of pixels on the dirty image in the desired direction
-epsilon : np.float64
-    the gridder epsilon
-
-Returns
-=======
-np.array((ndirty,), dtype=np.float64)
-    the correction factors.
-)""";
-pyarr_c<double> get_taper(size_t ndirty, double epsilon)
-  {
-  auto tmp = correction_factors(2*ndirty, ndirty, get_w(epsilon));
-  auto res = makeArray<double>({ndirty});
-  auto vres = res.mutable_data();
-  for (size_t i=0; i<ndirty; ++i)
-    vres[i] = tmp[i];
-  return res;
-  }
-
 template<typename T> struct UVW
   {
   T u, v, w;
@@ -507,6 +482,22 @@ np.array((nu, nv), dtype=np.float64)
     gridded UV data
 )""";
 
+constexpr auto apply_taper_DS = R"""(
+Applies the taper (or its inverse) to an image
+
+Parameters
+==========
+img: nd.array((nxdirty, nydirty), dtype=np.float64)
+    the image
+divide: bool
+    if True, the routine dividex by the taper, otherwise it multiplies by it
+
+Returns
+=======
+np.array((nxdirty, nydirty), dtype=np.float64)
+    the image with the taper applied
+)""";
+
 constexpr auto GridderConfig_DS = R"""(
 Class storing information related to the gridding/degridding process.
 
@@ -596,6 +587,25 @@ template<typename T> class GridderConfig
       }
       return res;
       }
+    pyarr_c<T> apply_taper(const pyarr_c<T> &img, bool divide) const
+      {
+      checkArray(img, "img", {nx_dirty, ny_dirty});
+      auto pin = img.data();
+      auto res = makeArray<T>({nx_dirty, ny_dirty});
+      auto pout = res.mutable_data();
+      {
+      py::gil_scoped_release release;
+      if (divide)
+        for (size_t i=0; i<nx_dirty; ++i)
+          for (size_t j=0; j<ny_dirty; ++j)
+            pout[ny_dirty*i + j] = pin[ny_dirty*i + j]/(cfu[i]*cfv[j]);
+      else
+        for (size_t i=0; i<nx_dirty; ++i)
+          for (size_t j=0; j<ny_dirty; ++j)
+            pout[ny_dirty*i + j] = pin[ny_dirty*i + j]*cfu[i]*cfv[j];
+      }
+      return res;
+      }
     pyarr_c<complex<T>> grid2dirty_c(const pyarr_c<complex<T>> &grid) const
       {
       checkArray(grid, "grid", {nu, nv});
@@ -1382,6 +1392,8 @@ PYBIND11_MODULE(nifty_gridder, m)
     .def("Pixsize_y", &GridderConfig<double>::Pixsize_y)
     .def("Nu", &GridderConfig<double>::Nu)
     .def("Nv", &GridderConfig<double>::Nv)
+    .def("apply_taper", &GridderConfig<double>::apply_taper, apply_taper_DS,
+      "img"_a, "divide"_a=false)
     .def("grid2dirty", &GridderConfig<double>::grid2dirty,
       grid2dirty_DS, "grid"_a)
     .def("grid2dirty_c", &GridderConfig<double>::grid2dirty_c, "grid"_a)
@@ -1437,5 +1449,4 @@ PYBIND11_MODULE(nifty_gridder, m)
     "idx"_a, "grid"_a, "wgt"_a, "ms_in"_a=None);
   m.def("apply_holo",&apply_holo<double>, "baselines"_a, "gconf"_a, "idx"_a,
     "grid"_a);
-  m.def("get_taper",&get_taper, get_taper_DS, "n"_a, "epsilon"_a);
   }

setup.py

@@ -32,6 +32,7 @@ else:
 def get_extension_modules():
     return [Extension('nifty_gridder',
                       sources=['nifty_gridder.cc'],
+                      depends=['pocketfft_hdronly.h', 'setup.py'],
                       include_dirs=include_dirs,
                       extra_compile_args=extra_compile_args,
                       extra_link_args=python_module_link_args)]