Commit 4475d484 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

getting uglier before it gets better

parent 4e95fdc1
Pipeline #81386 passed with stages
in 12 minutes and 20 seconds
...@@ -230,27 +230,48 @@ class Baselines ...@@ -230,27 +230,48 @@ class Baselines
}; };
template<typename T> class GridderConfig constexpr int logsquare=4;
template<typename T> class Params
{ {
protected: private:
size_t nx_dirty, ny_dirty, nu, nv; bool gridding;
double epsilon, ofactor; TimerHierarchy timers;
const mav<complex<T>,2> &ms_in;
mav<complex<T>,2> &ms_out;
const mav<T,2> &dirty_in;
mav<T,2> &dirty_out;
const mav<T,2> &wgt;
const mav<uint8_t,2> &mask;
double pixsize_x, pixsize_y;
size_t nxdirty, nydirty;
double epsilon;
bool do_wgridding;
size_t nthreads;
size_t verbosity;
bool negate_v, divide_by_n;
Baselines bl;
VVR ranges;
double wmin_d, wmax_d;
size_t nvis;
double wmin, dw;
size_t nplanes;
double nm1min;
vector<uint8_t> active;
size_t nu, nv;
double ofactor;
// FIXME: this should probably be done more cleanly
public:
TimerHierarchy &timers;
shared_ptr<HornerKernel<T>> krn; shared_ptr<HornerKernel<T>> krn;
protected:
double psx, psy;
size_t supp, nsafe; size_t supp, nsafe;
size_t nthreads;
double ushift, vshift; double ushift, vshift;
int maxiu0, maxiv0; int maxiu0, maxiv0;
size_t vlim; size_t vlim;
bool uv_side_fast; bool uv_side_fast;
T phase (T x, T y, T w, bool adjoint) const static T phase (T x, T y, T w, bool adjoint)
{ {
constexpr T pi = T(3.141592653589793238462643383279502884197); constexpr T pi = T(3.141592653589793238462643383279502884197);
T tmp = 1-x-y; T tmp = 1-x-y;
...@@ -261,68 +282,23 @@ template<typename T> class GridderConfig ...@@ -261,68 +282,23 @@ template<typename T> class GridderConfig
return phs; return phs;
} }
public:
GridderConfig(size_t nxdirty, size_t nydirty, size_t nu_, size_t nv_,
size_t kidx, double epsilon_, double pixsize_x, double pixsize_y,
const Baselines &baselines, size_t nthreads_, TimerHierarchy &timers_)
: nx_dirty(nxdirty), ny_dirty(nydirty), nu(nu_), nv(nv_),
epsilon(epsilon_),
ofactor(min(double(nu)/nxdirty, double(nv)/nydirty)),
timers(timers_),
krn(selectKernel<T>(ofactor, epsilon,kidx)),
psx(pixsize_x), psy(pixsize_y),
supp(krn->support()), nsafe((supp+1)/2),
nthreads(nthreads_),
ushift(supp*(-0.5)+1+nu), vshift(supp*(-0.5)+1+nv),
maxiu0((nu+nsafe)-supp), maxiv0((nv+nsafe)-supp),
vlim(min(nv/2, size_t(nv*baselines.Vmax()*psy+0.5*supp+1))),
uv_side_fast(true)
{
size_t vlim2 = (nydirty+1)/2+(supp+1)/2;
if (vlim2<vlim)
{
vlim = vlim2;
uv_side_fast = false;
}
MR_assert(nu>=2*nsafe, "nu too small");
MR_assert(nv>=2*nsafe, "nv too small");
MR_assert((nx_dirty&1)==0, "nx_dirty must be even");
MR_assert((ny_dirty&1)==0, "ny_dirty must be even");
MR_assert((nu&1)==0, "nu must be even");
MR_assert((nv&1)==0, "nv must be even");
MR_assert(epsilon>0, "epsilon must be positive");
MR_assert(pixsize_x>0, "pixsize_x must be positive");
MR_assert(pixsize_y>0, "pixsize_y must be positive");
}
size_t Nxdirty() const { return nx_dirty; }
size_t Nydirty() const { return ny_dirty; }
double Pixsize_x() const { return psx; }
double Pixsize_y() const { return psy; }
size_t Nu() const { return nu; }
size_t Nv() const { return nv; }
size_t Supp() const { return supp; }
size_t Nsafe() const { return nsafe; }
size_t Nthreads() const { return nthreads; }
double Epsilon() const { return epsilon; }
double Ofactor() const { return ofactor; }
void grid2dirty_post(mav<T,2> &tmav, void grid2dirty_post(mav<T,2> &tmav,
mav<T,2> &dirty) const mav<T,2> &dirty) const
{ {
checkShape(dirty.shape(), {nx_dirty, ny_dirty}); checkShape(dirty.shape(), {nxdirty, nydirty});
auto cfu = krn->corfunc(nx_dirty/2+1, 1./nu, nthreads); auto cfu = krn->corfunc(nxdirty/2+1, 1./nu, nthreads);
auto cfv = krn->corfunc(ny_dirty/2+1, 1./nv, nthreads); auto cfv = krn->corfunc(nydirty/2+1, 1./nv, nthreads);
execStatic(nx_dirty, nthreads, 0, [&](Scheduler &sched) execStatic(nxdirty, nthreads, 0, [&](Scheduler &sched)
{ {
while (auto rng=sched.getNext()) for(auto i=rng.lo; i<rng.hi; ++i) while (auto rng=sched.getNext()) for(auto i=rng.lo; i<rng.hi; ++i)
{ {
int icfu = abs(int(nx_dirty/2)-int(i)); int icfu = abs(int(nxdirty/2)-int(i));
for (size_t j=0; j<ny_dirty; ++j) for (size_t j=0; j<nydirty; ++j)
{ {
int icfv = abs(int(ny_dirty/2)-int(j)); int icfv = abs(int(nydirty/2)-int(j));
size_t i2 = nu-nx_dirty/2+i; size_t i2 = nu-nxdirty/2+i;
if (i2>=nu) i2-=nu; if (i2>=nu) i2-=nu;
size_t j2 = nv-ny_dirty/2+j; size_t j2 = nv-nydirty/2+j;
if (j2>=nv) j2-=nv; if (j2>=nv) j2-=nv;
dirty.v(i,j) = tmav(i2,j2)*T(cfu[icfu]*cfv[icfv]); dirty.v(i,j) = tmav(i2,j2)*T(cfu[icfu]*cfv[icfv]);
} }
...@@ -332,27 +308,27 @@ template<typename T> class GridderConfig ...@@ -332,27 +308,27 @@ template<typename T> class GridderConfig
void grid2dirty_post2( void grid2dirty_post2(
mav<complex<T>,2> &tmav, mav<T,2> &dirty, T w) const mav<complex<T>,2> &tmav, mav<T,2> &dirty, T w) const
{ {
checkShape(dirty.shape(), {nx_dirty,ny_dirty}); checkShape(dirty.shape(), {nxdirty,nydirty});
double x0 = -0.5*nx_dirty*psx, double x0 = -0.5*nxdirty*pixsize_x,
y0 = -0.5*ny_dirty*psy; y0 = -0.5*nydirty*pixsize_y;
execStatic(nx_dirty/2+1, nthreads, 0, [&](Scheduler &sched) execStatic(nxdirty/2+1, nthreads, 0, [&](Scheduler &sched)
{ {
using vtype = native_simd<T>; using vtype = native_simd<T>;
constexpr size_t vlen=vtype::size(); constexpr size_t vlen=vtype::size();
size_t nvec = (ny_dirty/2+1+(vlen-1))/vlen; size_t nvec = (nydirty/2+1+(vlen-1))/vlen;
vector<vtype> ph(nvec), sp(nvec), cp(nvec); vector<vtype> ph(nvec), sp(nvec), cp(nvec);
while (auto rng=sched.getNext()) for(auto i=rng.lo; i<rng.hi; ++i) while (auto rng=sched.getNext()) for(auto i=rng.lo; i<rng.hi; ++i)
{ {
T fx = T(x0+i*psx); T fx = T(x0+i*pixsize_x);
fx *= fx; fx *= fx;
size_t ix = nu-nx_dirty/2+i; size_t ix = nu-nxdirty/2+i;
if (ix>=nu) ix-=nu; if (ix>=nu) ix-=nu;
size_t i2 = nx_dirty-i; size_t i2 = nxdirty-i;
size_t ix2 = nu-nx_dirty/2+i2; size_t ix2 = nu-nxdirty/2+i2;
if (ix2>=nu) ix2-=nu; if (ix2>=nu) ix2-=nu;
for (size_t j=0; j<=ny_dirty/2; ++j) for (size_t j=0; j<=nydirty/2; ++j)
{ {
T fy = T(y0+j*psy); T fy = T(y0+j*pixsize_y);
ph[j/vlen][j%vlen] = phase(fx, fy*fy, w, true); ph[j/vlen][j%vlen] = phase(fx, fy*fy, w, true);
} }
for (size_t j=0; j<nvec; ++j) for (size_t j=0; j<nvec; ++j)
...@@ -362,17 +338,17 @@ template<typename T> class GridderConfig ...@@ -362,17 +338,17 @@ template<typename T> class GridderConfig
for (size_t k=0; k<vlen; ++k) for (size_t k=0; k<vlen; ++k)
cp[j][k]=cos(ph[j][k]); cp[j][k]=cos(ph[j][k]);
if ((i>0)&&(i<i2)) if ((i>0)&&(i<i2))
for (size_t j=0, jx=nv-ny_dirty/2; j<ny_dirty; ++j, jx=(jx+1>=nv)? jx+1-nv : jx+1) for (size_t j=0, jx=nv-nydirty/2; j<nydirty; ++j, jx=(jx+1>=nv)? jx+1-nv : jx+1)
{ {
size_t j2 = min(j, ny_dirty-j); size_t j2 = min(j, nydirty-j);
T re = cp[j2/vlen][j2%vlen], im = sp[j2/vlen][j2%vlen]; T re = cp[j2/vlen][j2%vlen], im = sp[j2/vlen][j2%vlen];
dirty.v(i,j) += tmav(ix,jx).real()*re - tmav(ix,jx).imag()*im; dirty.v(i,j) += tmav(ix,jx).real()*re - tmav(ix,jx).imag()*im;
dirty.v(i2,j) += tmav(ix2,jx).real()*re - tmav(ix2,jx).imag()*im; dirty.v(i2,j) += tmav(ix2,jx).real()*re - tmav(ix2,jx).imag()*im;
} }
else else
for (size_t j=0, jx=nv-ny_dirty/2; j<ny_dirty; ++j, jx=(jx+1>=nv)? jx+1-nv : jx+1) for (size_t j=0, jx=nv-nydirty/2; j<nydirty; ++j, jx=(jx+1>=nv)? jx+1-nv : jx+1)
{ {
size_t j2 = min(j, ny_dirty-j); size_t j2 = min(j, nydirty-j);
T re = cp[j2/vlen][j2%vlen], im = sp[j2/vlen][j2%vlen]; T re = cp[j2/vlen][j2%vlen], im = sp[j2/vlen][j2%vlen];
dirty.v(i,j) += tmav(ix,jx).real()*re - tmav(ix,jx).imag()*im; // lower left dirty.v(i,j) += tmav(ix,jx).real()*re - tmav(ix,jx).imag()*im; // lower left
} }
...@@ -380,7 +356,7 @@ template<typename T> class GridderConfig ...@@ -380,7 +356,7 @@ template<typename T> class GridderConfig
}); });
} }
void grid2dirty_overwrite(mav<T,2> &grid, mav<T,2> &dirty) const void grid2dirty_overwrite(mav<T,2> &grid, mav<T,2> &dirty)
{ {
timers.push("FFT"); timers.push("FFT");
checkShape(grid.shape(), {nu,nv}); checkShape(grid.shape(), {nu,nv});
...@@ -391,7 +367,7 @@ template<typename T> class GridderConfig ...@@ -391,7 +367,7 @@ template<typename T> class GridderConfig
} }
void grid2dirty_c_overwrite_wscreen_add void grid2dirty_c_overwrite_wscreen_add
(mav<complex<T>,2> &grid, mav<T,2> &dirty, T w) const (mav<complex<T>,2> &grid, mav<T,2> &dirty, T w)
{ {
timers.push("FFT"); timers.push("FFT");
checkShape(grid.shape(), {nu,nv}); checkShape(grid.shape(), {nu,nv});
...@@ -417,23 +393,23 @@ template<typename T> class GridderConfig ...@@ -417,23 +393,23 @@ template<typename T> class GridderConfig
void dirty2grid_pre(const mav<T,2> &dirty, void dirty2grid_pre(const mav<T,2> &dirty,
mav<T,2> &grid) const mav<T,2> &grid) const
{ {
checkShape(dirty.shape(), {nx_dirty, ny_dirty}); checkShape(dirty.shape(), {nxdirty, nydirty});
checkShape(grid.shape(), {nu, nv}); checkShape(grid.shape(), {nu, nv});
auto cfu = krn->corfunc(nx_dirty/2+1, 1./nu, nthreads); auto cfu = krn->corfunc(nxdirty/2+1, 1./nu, nthreads);
auto cfv = krn->corfunc(ny_dirty/2+1, 1./nv, nthreads); auto cfv = krn->corfunc(nydirty/2+1, 1./nv, nthreads);
// FIXME: maybe we don't have to fill everything and can save some time // FIXME: maybe we don't have to fill everything and can save some time
grid.fill(0); grid.fill(0);
execStatic(nx_dirty, nthreads, 0, [&](Scheduler &sched) execStatic(nxdirty, nthreads, 0, [&](Scheduler &sched)
{ {
while (auto rng=sched.getNext()) for(auto i=rng.lo; i<rng.hi; ++i) while (auto rng=sched.getNext()) for(auto i=rng.lo; i<rng.hi; ++i)
{ {
int icfu = abs(int(nx_dirty/2)-int(i)); int icfu = abs(int(nxdirty/2)-int(i));
for (size_t j=0; j<ny_dirty; ++j) for (size_t j=0; j<nydirty; ++j)
{ {
int icfv = abs(int(ny_dirty/2)-int(j)); int icfv = abs(int(nydirty/2)-int(j));
size_t i2 = nu-nx_dirty/2+i; size_t i2 = nu-nxdirty/2+i;
if (i2>=nu) i2-=nu; if (i2>=nu) i2-=nu;
size_t j2 = nv-ny_dirty/2+j; size_t j2 = nv-nydirty/2+j;
if (j2>=nv) j2-=nv; if (j2>=nv) j2-=nv;
grid.v(i2,j2) = dirty(i,j)*T(cfu[icfu]*cfv[icfv]); grid.v(i2,j2) = dirty(i,j)*T(cfu[icfu]*cfv[icfv]);
} }
...@@ -443,31 +419,31 @@ template<typename T> class GridderConfig ...@@ -443,31 +419,31 @@ template<typename T> class GridderConfig
void dirty2grid_pre2(const mav<T,2> &dirty, void dirty2grid_pre2(const mav<T,2> &dirty,
mav<complex<T>,2> &grid, T w) const mav<complex<T>,2> &grid, T w) const
{ {
checkShape(dirty.shape(), {nx_dirty, ny_dirty}); checkShape(dirty.shape(), {nxdirty, nydirty});
checkShape(grid.shape(), {nu, nv}); checkShape(grid.shape(), {nu, nv});
// FIXME: maybe we don't have to fill everything and can save some time // FIXME: maybe we don't have to fill everything and can save some time
grid.fill(0); grid.fill(0);
double x0 = -0.5*nx_dirty*psx, double x0 = -0.5*nxdirty*pixsize_x,
y0 = -0.5*ny_dirty*psy; y0 = -0.5*nydirty*pixsize_y;
execStatic(nx_dirty/2+1, nthreads, 0, [&](Scheduler &sched) execStatic(nxdirty/2+1, nthreads, 0, [&](Scheduler &sched)
{ {
using vtype = native_simd<T>; using vtype = native_simd<T>;
constexpr size_t vlen=vtype::size(); constexpr size_t vlen=vtype::size();
size_t nvec = (ny_dirty/2+1+(vlen-1))/vlen; size_t nvec = (nydirty/2+1+(vlen-1))/vlen;
vector<vtype> ph(nvec), sp(nvec), cp(nvec); vector<vtype> ph(nvec), sp(nvec), cp(nvec);
while (auto rng=sched.getNext()) for(auto i=rng.lo; i<rng.hi; ++i) while (auto rng=sched.getNext()) for(auto i=rng.lo; i<rng.hi; ++i)
{ {
T fx = T(x0+i*psx); T fx = T(x0+i*pixsize_x);
fx *= fx; fx *= fx;
size_t ix = nu-nx_dirty/2+i; size_t ix = nu-nxdirty/2+i;
if (ix>=nu) ix-=nu; if (ix>=nu) ix-=nu;
size_t i2 = nx_dirty-i; size_t i2 = nxdirty-i;
size_t ix2 = nu-nx_dirty/2+i2; size_t ix2 = nu-nxdirty/2+i2;
if (ix2>=nu) ix2-=nu; if (ix2>=nu) ix2-=nu;
for (size_t j=0; j<=ny_dirty/2; ++j) for (size_t j=0; j<=nydirty/2; ++j)
{ {
T fy = T(y0+j*psy); T fy = T(y0+j*pixsize_y);
ph[j/vlen][j%vlen] = phase(fx, fy*fy, w, false); ph[j/vlen][j%vlen] = phase(fx, fy*fy, w, false);
} }
for (size_t j=0; j<nvec; ++j) for (size_t j=0; j<nvec; ++j)
...@@ -477,17 +453,17 @@ template<typename T> class GridderConfig ...@@ -477,17 +453,17 @@ template<typename T> class GridderConfig
for (size_t k=0; k<vlen; ++k) for (size_t k=0; k<vlen; ++k)
cp[j][k]=cos(ph[j][k]); cp[j][k]=cos(ph[j][k]);
if ((i>0)&&(i<i2)) if ((i>0)&&(i<i2))
for (size_t j=0, jx=nv-ny_dirty/2; j<ny_dirty; ++j, jx=(jx+1>=nv)? jx+1-nv : jx+1) for (size_t j=0, jx=nv-nydirty/2; j<nydirty; ++j, jx=(jx+1>=nv)? jx+1-nv : jx+1)
{ {
size_t j2 = min(j, ny_dirty-j); size_t j2 = min(j, nydirty-j);
auto ws = complex<T>(cp[j2/vlen][j2%vlen],sp[j2/vlen][j2%vlen]); auto ws = complex<T>(cp[j2/vlen][j2%vlen],sp[j2/vlen][j2%vlen]);
grid.v(ix,jx) = dirty(i,j)*ws; // lower left grid.v(ix,jx) = dirty(i,j)*ws; // lower left
grid.v(ix2,jx) = dirty(i2,j)*ws; // lower right grid.v(ix2,jx) = dirty(i2,j)*ws; // lower right
} }
else else
for (size_t j=0, jx=nv-ny_dirty/2; j<ny_dirty; ++j, jx=(jx+1>=nv)? jx+1-nv : jx+1) for (size_t j=0, jx=nv-nydirty/2; j<nydirty; ++j, jx=(jx+1>=nv)? jx+1-nv : jx+1)
{ {
size_t j2 = min(j, ny_dirty-j); size_t j2 = min(j, nydirty-j);
auto ws = complex<T>(cp[j2/vlen][j2%vlen],sp[j2/vlen][j2%vlen]); auto ws = complex<T>(cp[j2/vlen][j2%vlen],sp[j2/vlen][j2%vlen]);
grid.v(ix,jx) = dirty(i,j)*ws; // lower left grid.v(ix,jx) = dirty(i,j)*ws; // lower left
} }
...@@ -495,8 +471,7 @@ template<typename T> class GridderConfig ...@@ -495,8 +471,7 @@ template<typename T> class GridderConfig
}); });
} }
void dirty2grid(const mav<T,2> &dirty, void dirty2grid(const mav<T,2> &dirty, mav<T,2> &grid)
mav<T,2> &grid) const
{ {
timers.push("grid correction"); timers.push("grid correction");
dirty2grid_pre(dirty, grid); dirty2grid_pre(dirty, grid);
...@@ -506,7 +481,7 @@ template<typename T> class GridderConfig ...@@ -506,7 +481,7 @@ template<typename T> class GridderConfig
} }
void dirty2grid_c_wscreen(const mav<T,2> &dirty, void dirty2grid_c_wscreen(const mav<T,2> &dirty,
mav<complex<T>,2> &grid, T w) const mav<complex<T>,2> &grid, T w)
{ {
timers.push("wscreen+grid correction"); timers.push("wscreen+grid correction");
dirty2grid_pre2(dirty, grid, w); dirty2grid_pre2(dirty, grid, w);
...@@ -530,54 +505,23 @@ template<typename T> class GridderConfig ...@@ -530,54 +505,23 @@ template<typename T> class GridderConfig
[[gnu::always_inline]] void getpix(double u_in, double v_in, double &u, double &v, int &iu0, int &iv0) const [[gnu::always_inline]] void getpix(double u_in, double v_in, double &u, double &v, int &iu0, int &iv0) const
{ {
u=fmod1(u_in*psx)*nu; u=fmod1(u_in*pixsize_x)*nu;
iu0 = min(int(u+ushift)-int(nu), maxiu0); iu0 = min(int(u+ushift)-int(nu), maxiu0);
v=fmod1(v_in*psy)*nv; v=fmod1(v_in*pixsize_y)*nv;
iv0 = min(int(v+vshift)-int(nv), maxiv0); iv0 = min(int(v+vshift)-int(nv), maxiv0);
} }
};
constexpr int logsquare=4;
template<typename T> class Params void report()
{
private:
bool gridding;
TimerHierarchy timers;
const mav<complex<T>,2> &ms_in;
mav<complex<T>,2> &ms_out;
const mav<T,2> &dirty_in;
mav<T,2> &dirty_out;
const mav<T,2> &wgt;
const mav<uint8_t,2> &mask;
double pixsize_x, pixsize_y;
size_t nxdirty, nydirty;
double epsilon;
bool do_wgridding;
size_t nthreads;
size_t verbosity;
bool negate_v, divide_by_n;
Baselines bl;
VVR ranges;
double wmin_d, wmax_d;
size_t nvis;
double wmin, dw;
size_t nplanes;
double nm1min;
vector<uint8_t> active;
void report(const GridderConfig<T> &gconf)
{ {
if (verbosity==0) return; if (verbosity==0) return;
cout << (gridding ? "Gridding" : "Degridding") cout << (gridding ? "Gridding" : "Degridding")
<< ": nthreads=" << nthreads << ", " << ": nthreads=" << nthreads << ", "
<< "dirty=(" << nxdirty << "x" << nydirty << "), " << "dirty=(" << nxdirty << "x" << nydirty << "), "
<< "grid=(" << gconf.Nu() << "x" << gconf.Nv(); << "grid=(" << nu << "x" << nv;
if (nplanes>0) cout << "x" << nplanes; if (nplanes>0) cout << "x" << nplanes;
cout << "), nvis=" << nvis cout << "), nvis=" << nvis
<< ", supp=" << gconf.Supp() << ", supp=" << supp
<< ", eps=" << (gconf.Epsilon() * ((nplanes==0) ? 2 : 3)) << ", eps=" << (epsilon * ((nplanes==0) ? 2 : 3))
<< endl; << endl;
cout << " w=[" << wmin_d << "; " << wmax_d << "], min(n-1)=" << nm1min << ", dw=" << dw cout << " w=[" << wmin_d << "; " << wmax_d << "], min(n-1)=" << nm1min << ", dw=" << dw
<< ", wmax/dw=" << wmax_d/dw << ", nranges=" << ranges.size() << endl; << ", wmax/dw=" << wmax_d/dw << ", nranges=" << ranges.size() << endl;
...@@ -671,19 +615,18 @@ auto getNuNv() ...@@ -671,19 +615,18 @@ auto getNuNv()
} }
} }
timers.pop(); timers.pop();
return make_tuple(minnu, minnv, minidx); nu = minnu;
nv = minnv;
return minidx;
} }
void countRanges(const GridderConfig<T> &gconf) void countRanges()
{ {
timers.push("range count"); timers.push("range count");
size_t nrow=bl.Nrows(), size_t nrow=bl.Nrows(),
nchan=bl.Nchannels(), nchan=bl.Nchannels();
nsafe=gconf.Nsafe(),
nthreads = gconf.Nthreads(),
supp = gconf.Supp();
dw = 0.5/gconf.Ofactor()/abs(nm1min); dw = 0.5/ofactor/abs(nm1min);
nplanes = size_t((wmax_d-wmin_d)/dw+supp); nplanes = size_t((wmax_d-wmin_d)/dw+supp);
wmin = (wmin_d+wmax_d)*0.5 - 0.5*(nplanes-1)*dw; wmin = (wmin_d+wmax_d)*0.5 - 0.5*(nplanes-1)*dw;
...@@ -712,7 +655,7 @@ void countRanges(const GridderConfig<T> &gconf) ...@@ -712,7 +655,7 @@ void countRanges(const GridderConfig<T> &gconf)
if (uvw.w<0) uvw.Flip(); if (uvw.w<0) uvw.Flip();
double u, v; double u, v;
int iu0, iv0, iw; int iu0, iv0, iw;
gconf.getpix(uvw.u, uvw.v, u, v, iu0, iv0); getpix(uvw.u, uvw.v, u, v, iu0, iv0);
iu0 = (iu0+nsafe)>>logsquare; iu0 = (iu0+nsafe)>>logsquare;
iv0 = (iv0+nsafe)>>logsquare; iv0 = (iv0+nsafe)>>logsquare;
iw = max(0,int(1+(abs(uvw.w)-(0.5*supp*dw)-wmin)/dw)); iw = max(0,int(1+(abs(uvw.w)-(0.5*supp*dw)-wmin)/dw));
...@@ -764,13 +707,13 @@ void countRanges(const GridderConfig<T> &gconf) ...@@ -764,13 +707,13 @@ void countRanges(const GridderConfig<T> &gconf)
timers.pop(); timers.pop();
} }
void apply_global_corrections(const GridderConfig<T> &gconf, mav<T,2> &dirty) void apply_global_corrections(mav<T,2> &dirty)
{ {
timers.push("global corrections"); timers.push("global corrections");
double x0 = -0.5*nxdirty*pixsize_x, double x0 = -0.5*nxdirty*pixsize_x,
y0 = -0.5*nydirty*pixsize_y; y0 = -0.5*nydirty*pixsize_y;
auto cfu = gconf.krn->corfunc(nxdirty/2+1, 1./gconf.Nu(), nthreads); auto cfu = krn->corfunc(nxdirty/2+1, 1./nu, nthreads);
auto cfv = gconf.krn->corfunc(nydirty/2+1, 1./gconf.Nv(), nthreads); auto cfv = krn->corfunc(nydirty/2+1, 1./nv, nthreads);
execStatic(nxdirty/2+1, nthreads, 0, [&](Scheduler &sched) execStatic(nxdirty/2+1, nthreads, 0, [&](Scheduler &sched)
{ {
while (auto rng=sched.getNext()) for(auto i=rng.lo; i<rng.hi; ++i) while (auto rng=sched.getNext()) for(auto i=rng.lo; i<rng.hi; ++i)
...@@ -786,7 +729,7 @@ void apply_global_corrections(const GridderConfig<T> &gconf, mav<T,2> &dirty) ...@@ -786,7 +729,7 @@ void apply_global_corrections(const GridderConfig<T> &gconf, mav<T,2> &dirty)
if (tmp>=0) if (tmp>=0)
{ {
auto nm1 = (-fx-fy)/(sqrt(tmp)+1); // accurate form of sqrt(1-x-y)-1 auto nm1 = (-fx-fy)/(sqrt(tmp)+1); // accurate form of sqrt(1-x-y)-1
fct = T(gconf.krn->corfunc(nm1*dw));