Commit 355cba73 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

add simple (non w-stacking) high-level routines

parent a4356ac6
......@@ -991,6 +991,55 @@ template<typename T, typename T2> void apply_wcorr(const GridderConfig<T> &gconf
}
}
template<typename T, typename Serv> void x2dirty(
const GridderConfig<T> &gconf, const Serv &srv, const mav<T,2> &dirty, size_t verbosity=0)
{
auto nx_dirty=gconf.Nxdirty();
auto ny_dirty=gconf.Nydirty();
auto nu=gconf.Nu();
auto nv=gconf.Nv();
if (verbosity>0)
cout << "Gridding without w-stacking: " << srv.Nvis()
<< " visibilities" << endl;
tmpStorage<complex<T>,2> grid_({nu,nv});
auto grid=grid_.getMav();
grid_.fill(0.);
x2grid_c(gconf, srv, grid);
tmpStorage<complex<T>,2> cdirty_({nx_dirty,ny_dirty});
auto cdirty=cdirty_.getMav();
gconf.grid2dirty_c(grid, cdirty);
for (size_t i=0; i<nx_dirty; ++i)
for (size_t j=0; j<ny_dirty; ++j)
dirty(i,j) = cdirty(i,j).real();
}
template<typename T, typename Serv> void dirty2x(
const GridderConfig<T> &gconf, const const_mav<T,2> &dirty,
const Serv &srv, size_t verbosity=0)
{
auto nx_dirty=gconf.Nxdirty();
auto ny_dirty=gconf.Nydirty();
auto nu=gconf.Nu();
auto nv=gconf.Nv();
size_t nvis = srv.Nvis();
for (size_t i=0; i<nvis; ++i)
srv.setVis(i, 0.);
tmpStorage<complex<T>,2> cdirty_({nx_dirty,ny_dirty});
auto cdirty=cdirty_.getMav();
for (size_t i=0; i<nx_dirty; ++i)
for (size_t j=0; j<ny_dirty; ++j)
cdirty(i,j) = dirty(i,j);
tmpStorage<complex<T>,2> grid_({nu,nv});
auto grid=grid_.getMav();
gconf.dirty2grid_c(cdirty, grid);
grid2x_c(gconf, const_mav<complex<T>,2>(grid.data(),{nu,nv}), srv);
}
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, size_t verbosity)
......@@ -1317,6 +1366,35 @@ template<typename T> vector<uint32_t> getWgtIndices(const Baselines<T> &baseline
return res;
}
template<typename T> void gridding(const const_mav<T,2> &uvw,
const const_mav<T,1> &freq, const const_mav<complex<T>,2> &ms,
const const_mav<T,2> &wgt, T pixsize_x, T pixsize_y, double epsilon,
size_t nthreads, const mav<T,2> &dirty, size_t verbosity)
{
Baselines<T> baselines(uvw, freq);
GridderConfig<T> gconf(dirty.shape(0), dirty.shape(1), epsilon, pixsize_x, pixsize_y, nthreads);
auto idx = getWgtIndices(baselines, gconf, wgt);
auto idx2 = const_mav<uint32_t,1>(idx.data(),{idx.size()});
auto srv = MsServ<T,const_mav<complex<T>,2>,const_mav<T,2>>(baselines,idx2,ms,wgt);
x2dirty(gconf, srv, dirty, verbosity);
}
template<typename T> void degridding(const const_mav<T,2> &uvw,
const const_mav<T,1> &freq, const const_mav<T,2> &dirty,
const const_mav<T,2> &wgt, T pixsize_x, T pixsize_y, double epsilon,
size_t nthreads, const mav<complex<T>,2> &ms, size_t verbosity)
{
Baselines<T> baselines(uvw, freq);
GridderConfig<T> gconf(dirty.shape(0), dirty.shape(1), epsilon, pixsize_x, pixsize_y, nthreads);
auto idx = getWgtIndices(baselines, gconf, wgt);
auto idx2 = const_mav<uint32_t,1>(idx.data(),{idx.size()});
for (size_t i=0; i<ms.shape(0); ++i)
for (size_t j=0; j<ms.shape(1); ++j)
ms(i,j) = 0;
auto srv = MsServ<T,mav<complex<T>,2>,const_mav<T,2>>(baselines,idx2,ms,wgt);
dirty2x(gconf, dirty, srv, verbosity);
}
template<typename T> void full_gridding(const const_mav<T,2> &uvw,
const const_mav<T,1> &freq, const const_mav<complex<T>,2> &ms,
const const_mav<T,2> &wgt, T pixsize_x, T pixsize_y, double epsilon,
......
......@@ -643,6 +643,37 @@ template<typename T> pyarr<complex<T>> Pydirty2vis_wstack(const PyBaselines<T> &
return vis;
}
template<typename T> pyarr<T> Pygridding(const pyarr<T> &uvw,
const pyarr<T> &freq, const pyarr<complex<T>> &ms, const py::object &wgt_,
size_t npix_x, size_t npix_y, T pixsize_x, T pixsize_y, double epsilon,
size_t nthreads, size_t verbosity)
{
auto uvw2 = make_const_mav<2>(uvw);
auto freq2 = make_const_mav<1>(freq);
auto ms2 = make_const_mav<2>(ms);
auto wgt=providePotentialArray<T>(wgt_,{ms2.shape(0),ms2.shape(1)});
auto wgt2 = make_const_mav<2>(wgt);
auto dirty=makeArray<T>({npix_x,npix_y});
auto dirty2 =make_mav<2>(dirty);
gridding(uvw2,freq2,ms2,wgt2,pixsize_x,pixsize_y,epsilon,nthreads,dirty2,verbosity);
return dirty;
}
template<typename T> pyarr<complex<T>> Pydegridding(const pyarr<T> &uvw,
const pyarr<T> &freq, const pyarr<T> &dirty, const py::object &wgt_,
T pixsize_x, T pixsize_y, double epsilon, size_t nthreads, size_t verbosity)
{
auto uvw2 = make_const_mav<2>(uvw);
auto freq2 = make_const_mav<1>(freq);
auto dirty2 = make_const_mav<2>(dirty);
auto wgt=providePotentialArray<T>(wgt_,{uvw2.shape(0),freq2.shape(0)});
auto wgt2 = make_const_mav<2>(wgt);
auto ms=makeArray<complex<T>>({uvw2.shape(0),freq2.shape(0)});
auto ms2 =make_mav<2>(ms);
degridding(uvw2,freq2,dirty2,wgt2,pixsize_x,pixsize_y,epsilon,nthreads,ms2,verbosity);
return ms;
}
template<typename T> pyarr<T> Pyfull_gridding(const pyarr<T> &uvw,
const pyarr<T> &freq, const pyarr<complex<T>> &ms, const py::object &wgt_,
size_t npix_x, size_t npix_y, T pixsize_x, T pixsize_y, double epsilon,
......@@ -758,4 +789,16 @@ PYBIND11_MODULE(nifty_gridder, m)
m.def("full_degridding_f",&Pyfull_degridding<float>,"uvw"_a,"freq"_a,"dirty"_a,
"wgt"_a,"pixsize_x"_a,"pixsize_y"_a, "epsilon"_a, "nthreads"_a=1,
"verbosity"_a=0);
m.def("gridding",&Pygridding<double>,"uvw"_a,"freq"_a,"ms"_a,
"wgt"_a,"npix_x"_a,"npix_y"_a,"pixsize_x"_a,"pixsize_y"_a,"epsilon"_a,
"nthreads"_a=1, "verbosity"_a=0);
m.def("gridding_f",&Pygridding<float>,"uvw"_a,"freq"_a,"ms"_a,
"wgt"_a,"npix_x"_a,"npix_y"_a,"pixsize_x"_a,"pixsize_y"_a,"epsilon"_a,
"nthreads"_a=1, "verbosity"_a=0);
m.def("degridding",&Pydegridding<double>,"uvw"_a,"freq"_a,"dirty"_a,
"wgt"_a,"pixsize_x"_a,"pixsize_y"_a, "epsilon"_a,"nthreads"_a=1,
"verbosity"_a=0);
m.def("degridding_f",&Pydegridding<float>,"uvw"_a,"freq"_a,"dirty"_a,
"wgt"_a,"pixsize_x"_a,"pixsize_y"_a, "epsilon"_a, "nthreads"_a=1,
"verbosity"_a=0);
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment