merge improvements from minimal branch

gridder_cxx.h

@@ -605,13 +605,8 @@ class GridderConfig
             img2(i,j) = img(i,j)*cfu[i]*cfv[j];
       }
 
-    template<typename T> void grid2dirty(const const_mav<T,2> &grid, const mav<T,2> &dirty) const
+    template<typename T> void grid2dirty_post(const mav<T,2> &tmav, const mav<T,2> &dirty) const
       {
-      checkShape(grid.shape(), {nu,nv});
-      checkShape(dirty.shape(), {nx_dirty,ny_dirty});
-      tmpStorage<T,2> tmpdat({nu,nv});
-      auto tmav = tmpdat.getMav();
-      hartley2_2D<T>(grid, tmav, nthreads);
       for (size_t i=0; i<nx_dirty; ++i)
         for (size_t j=0; j<ny_dirty; ++j)
           {
@@ -622,6 +617,15 @@ class GridderConfig
           dirty(i,j) = tmav(i2,j2)*T(cfu[i]*cfv[j]);
           }
       }
+    template<typename T> void grid2dirty(const const_mav<T,2> &grid, const mav<T,2> &dirty) const
+      {
+      checkShape(grid.shape(), {nu,nv});
+      checkShape(dirty.shape(), {nx_dirty,ny_dirty});
+      tmpStorage<T,2> tmpdat({nu,nv});
+      auto tmav = tmpdat.getMav();
+      hartley2_2D<T>(grid, tmav, nthreads);
+      grid2dirty_post(tmav, dirty);
+      }
 
     template<typename T> void grid2dirty_c(const mav<const complex<T>,2> &grid, mav<complex<T>,2> &dirty) const
       {
@@ -633,18 +637,10 @@ class GridderConfig
       pocketfft::c2c({nu,nv},{grid.stride(0)*sc,grid.stride(1)*sc},
         {tmp.stride(0)*sc, tmp.stride(1)*sc}, {0,1}, pocketfft::BACKWARD,
         grid.data(), tmp.data(), T(1), nthreads);
-      for (size_t i=0; i<nx_dirty; ++i)
-        for (size_t j=0; j<ny_dirty; ++j)
-          {
-          size_t i2 = nu-nx_dirty/2+i;
-          if (i2>=nu) i2-=nu;
-          size_t j2 = nv-ny_dirty/2+j;
-          if (j2>=nv) j2-=nv;
-          dirty(i,j) = tmp(i2,j2)*T(cfu[i]*cfv[j]);
-          }
+      grid2dirty_post(tmp, dirty);
       }
 
-    template<typename T> void dirty2grid(const const_mav<T,2> &dirty, mav<T,2> &grid) const
+    template<typename T> void dirty2grid_pre(const const_mav<T,2> &dirty, mav<T,2> &grid) const
       {
       checkShape(dirty.shape(), {nx_dirty, ny_dirty});
       checkShape(grid.shape(), {nu, nv});
@@ -658,24 +654,17 @@ class GridderConfig
           if (j2>=nv) j2-=nv;
           grid(i2,j2) = dirty(i,j)*T(cfu[i]*cfv[j]);
           }
+      }
+    template<typename T> void dirty2grid(const const_mav<T,2> &dirty, mav<T,2> &grid) const
+      {
+      dirty2grid_pre(dirty, grid);
       hartley2_2D<T>(cmav(grid), grid, nthreads);
       }
 
     template<typename T> void dirty2grid_c(const const_mav<complex<T>,2> &dirty,
       mav<complex<T>,2> &grid) const
       {
-      checkShape(dirty.shape(), {nx_dirty, ny_dirty});
-      checkShape(grid.shape(), {nu, nv});
-      grid.fill(0);
-      for (size_t i=0; i<nx_dirty; ++i)
-        for (size_t j=0; j<ny_dirty; ++j)
-          {
-          size_t i2 = nu-nx_dirty/2+i;
-          if (i2>=nu) i2-=nu;
-          size_t j2 = nv-ny_dirty/2+j;
-          if (j2>=nv) j2-=nv;
-          grid(i2,j2) = dirty(i,j)*T(cfu[i]*cfv[j]);
-          }
+      dirty2grid_pre(dirty, grid);
       constexpr auto sc = ptrdiff_t(sizeof(complex<T>));
       pocketfft::stride_t strides{grid.stride(0)*sc,grid.stride(1)*sc};
       pocketfft::c2c({nu,nv}, strides, strides, {0,1}, pocketfft::FORWARD,
@@ -1216,10 +1205,8 @@ template<typename Serv> void wstack_common(
   size_t nvis = srv.Nvis();
   size_t nthreads = gconf.Nthreads();
 
-  // determine w values for every visibility, and min/max w;
   wmin=1e38;
   double wmax=-1e38;
-
 // FIXME maybe this can be done more intelligently
 #pragma omp parallel for num_threads(nthreads) reduction(min:wmin) reduction(max:wmax)
   for (size_t ipart=0; ipart<nvis; ++ipart)
@@ -1283,7 +1270,8 @@ template<typename T, typename Serv> void x2dirty(
 
     dirty.fill(0);
     vector<complex<T>> vis_loc_;
-    vector<uint32_t> idx_loc_, newidx;
+    vector<uint32_t> idx_loc_;
+    vector<uint32_t>  newidx;
     tmpStorage<complex<T>,2> grid_({gconf.Nu(),gconf.Nv()});
     auto grid=grid_.getMav();
     tmpStorage<complex<T>,2> tdirty_(dirty.shape());
@@ -1304,7 +1292,8 @@ template<typename T, typename Serv> void x2dirty(
           if ((int(iw)>=minplane[ipart]) && (iw<minplane[ipart]+supp))
             newidx[cnt++] = srv.getIdx(ipart);
         grid.fill(0);
-        Serv newsrv(srv, const_mav<uint32_t, 1>(newidx.data(), {newidx.size()}));
+        auto newidx2 = const_mav<uint32_t, 1>(newidx.data(), {newidx.size()});
+        Serv newsrv(srv, newidx2);
         x2grid_c(gconf, newsrv, grid, wcur, dw);
         }
       else
@@ -1389,7 +1378,8 @@ template<typename T, typename Serv> void dirty2x(
     apply_wcorr(gconf, tdirty, ES_Kernel_with_correction(supp, nthreads), dw);
 
     vector<complex<T>> vis_loc_;
-    vector<uint32_t> idx_loc_, newidx;
+    vector<uint32_t> idx_loc_;
+    vector<uint32_t> newidx;
     tmpStorage<complex<T>,2> grid_({gconf.Nu(),gconf.Nv()});
     auto grid=grid_.getMav();
     tmpStorage<complex<T>,2> tdirty2_({nx_dirty,ny_dirty});