factor out common code

gridder_cxx.h

@@ -983,21 +983,20 @@ template<typename T, typename T2> void apply_wcorr(const GridderConfig<T> &gconf
     }
   }
 
-template<typename T, typename Serv> void x2dirty_wstack(const Baselines<T> &baselines,
-  const GridderConfig<T> &gconf, const Serv &srv, const mav<complex<T>,2> &dirty)
+template<typename T, typename Serv> void wstack_common(
+  const GridderConfig<T> &gconf, const Serv &srv, T &wmin, vector<T> &wval,
+  double &dw, size_t &nplanes, vector<size_t> &nvis_plane, vector<int> &minplane)
   {
   auto nx_dirty=gconf.Nxdirty();
   auto ny_dirty=gconf.Nydirty();
-  auto nu=gconf.Nu();
-  auto nv=gconf.Nv();
   auto psx=gconf.Pixsize_x();
   auto psy=gconf.Pixsize_y();
   size_t nvis = srv.Nvis();
-  checkShape(dirty.shape(),{nx_dirty,ny_dirty});
 
   // determine w values for every visibility, and min/max w;
-  T wmin=T(1e38), wmax=T(-1e38);
-  vector<T> wval(nvis);
+  wmin=T(1e38);
+  T wmax=T(-1e38);
+  wval.resize(nvis);
 
 #pragma omp parallel for num_threads(nthreads) reduction(min:wmin) reduction(max:wmax)
   for (size_t ipart=0; ipart<nvis; ++ipart)
@@ -1010,17 +1009,16 @@ template<typename T, typename Serv> void x2dirty_wstack(const Baselines<T> &base
   double x0 = -0.5*nx_dirty*psx,
          y0 = -0.5*ny_dirty*psy;
   double nmin = sqrt(max(1.-x0*x0-y0*y0,0.))-1.;
-  double dw = 0.25/abs(nmin);
-  size_t nplanes = max<size_t>(2,((wmax-wmin)/dw+2));
+  dw = 0.25/abs(nmin);
+  nplanes = max<size_t>(2,((wmax-wmin)/dw+2));
   dw=(wmax-wmin)/(nplanes-1);
 
   auto w_supp = gconf.Supp();
-  EC_Kernel_with_correction<T> kernel(w_supp);
   wmin -= (0.5*w_supp-1)*dw;
   wmax += (0.5*w_supp-1)*dw;
   nplanes += w_supp-2;
-  vector<size_t> nvis_plane(nplanes,0);
-  vector<int> minplane(nvis);
+  nvis_plane.resize(nplanes,0);
+  minplane.resize(nvis);
 #pragma omp parallel num_threads(nthreads)
 {
   vector<size_t> nvploc(nplanes,0);
@@ -1036,6 +1034,26 @@ template<typename T, typename Serv> void x2dirty_wstack(const Baselines<T> &base
   for (size_t i=0; i<nplanes; ++i)
     nvis_plane[i]+=nvploc[i];
 }
+  }
+
+template<typename T, typename Serv> void x2dirty_wstack(const Baselines<T> &baselines,
+  const GridderConfig<T> &gconf, const Serv &srv, const mav<complex<T>,2> &dirty)
+  {
+  auto nx_dirty=gconf.Nxdirty();
+  auto ny_dirty=gconf.Nydirty();
+  auto nu=gconf.Nu();
+  auto nv=gconf.Nv();
+  size_t nvis = srv.Nvis();
+  auto w_supp = gconf.Supp();
+  EC_Kernel_with_correction<T> kernel(w_supp);
+  T wmin;
+  vector<T> wval;
+  double dw;
+  size_t nplanes;
+  vector<size_t> nvis_plane;
+  vector<int> minplane;
+  wstack_common(gconf, srv, wmin, wval, dw, nplanes, nvis_plane, minplane);
+
   for (size_t i=0; i<nx_dirty; ++i)
     for (size_t j=0; j<ny_dirty; ++j)
       dirty(i,j) = 0.;
@@ -1098,52 +1116,16 @@ template<typename T, typename Serv> void dirty2x_wstack(const Baselines<T> &base
   auto ny_dirty=gconf.Nydirty();
   auto nu=gconf.Nu();
   auto nv=gconf.Nv();
-  auto psx=gconf.Pixsize_x();
-  auto psy=gconf.Pixsize_y();
   size_t nvis = srv.Nvis();
-  checkShape(dirty.shape(),{nx_dirty,ny_dirty});
-
-  // determine w values for every visibility, and min/max w;
-  T wmin=T(1e38), wmax=T(-1e38);
-  vector<T> wval(nvis);
-
-#pragma omp parallel for num_threads(nthreads) reduction(min:wmin) reduction(max:wmax)
-  for (size_t ipart=0; ipart<nvis; ++ipart)
-    {
-    wval[ipart] = srv.getCoord(ipart).w;
-    wmin = min(wmin,wval[ipart]);
-    wmax = max(wmax,wval[ipart]);
-    }
-
-  double x0 = -0.5*nx_dirty*psx,
-         y0 = -0.5*ny_dirty*psy;
-  double nmin = sqrt(max(1.-x0*x0-y0*y0,0.))-1.;
-  double dw = 0.25/abs(nmin);
-  size_t nplanes = max<size_t>(2,((wmax-wmin)/dw+2));
-  dw=(wmax-wmin)/(nplanes-1);
-
   auto w_supp = gconf.Supp();
   EC_Kernel_with_correction<T> kernel(w_supp);
-  wmin -= (0.5*w_supp-1)*dw;
-  wmax += (0.5*w_supp-1)*dw;
-  nplanes += w_supp-2;
-  vector<size_t> nvis_plane(nplanes,0);
-  vector<int> minplane(nvis);
-#pragma omp parallel num_threads(nthreads)
-{
-  vector<size_t> nvploc(nplanes,0);
-#pragma omp for
-  for (size_t ipart=0; ipart<nvis; ++ipart)
-    {
-    int plane0 = int((wval[ipart]-wmin)/dw+0.5*w_supp)-int(w_supp-1);
-    minplane[ipart]=plane0;
-    for (int i=max<int>(0,plane0); i<min<int>(nplanes,plane0+w_supp); ++i)
-      ++nvploc[i];
-    }
-#pragma omp critical
-  for (size_t i=0; i<nplanes; ++i)
-    nvis_plane[i]+=nvploc[i];
-}
+  T wmin;
+  vector<T> wval;
+  double dw;
+  size_t nplanes;
+  vector<size_t> nvis_plane;
+  vector<int> minplane;
+  wstack_common(gconf, srv, wmin, wval, dw, nplanes, nvis_plane, minplane);
 
   for (size_t i=0; i<nvis; ++i)
     srv.setVis(i, 0.);