gridder.py 4.29 KB
Newer Older
Philipp Arras's avatar
Philipp Arras committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2019 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.

Philipp Arras's avatar
Philipp Arras committed
18
19
20
21
22
import numpy as np

from ..domain_tuple import DomainTuple
from ..domains.rg_space import RGSpace
from ..domains.unstructured_domain import UnstructuredDomain
23
from ..fft import hartley
Philipp Arras's avatar
Philipp Arras committed
24
25
26
27
28
from ..operators.linear_operator import LinearOperator
from ..sugar import from_global_data, makeDomain


class GridderMaker(object):
Martin Reinecke's avatar
Martin Reinecke committed
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
    def __init__(self, dirty_domain, uv, eps=2e-13):
        import nifty_gridder
        dirty_domain = makeDomain(dirty_domain)
        if (len(dirty_domain) != 1 or
                not isinstance(dirty_domain[0], RGSpace) or
                not len(dirty_domain.shape) == 2):
            raise ValueError("need dirty_domain with exactly one 2D RGSpace")
        bl = nifty_gridder.Baselines(uv, np.array([1.]));
        nxdirty, nydirty = dirty_domain.shape
        gconf = nifty_gridder.GridderConfig(nxdirty, nydirty, eps, 1., 1.)
        nu = gconf.Nu()
        nv = gconf.Nv()
        idx = bl.getIndices()
        idx = gconf.reorderIndices(idx, bl)

        grid_domain = RGSpace([nu, nv], distances=[1, 1], harmonic=False)

        self._rest = _RestOperator(dirty_domain, grid_domain, gconf)
        self._gridder = RadioGridder(grid_domain, bl, gconf, idx)

    def getGridder(self):
        return self._gridder
Philipp Arras's avatar
Philipp Arras committed
51
52
53
54

    def getRest(self):
        return self._rest

Martin Reinecke's avatar
Martin Reinecke committed
55
56
    def getFull(self):
        return self.getRest() @ self._gridder
Philipp Arras's avatar
Philipp Arras committed
57
58
59


class _RestOperator(LinearOperator):
Martin Reinecke's avatar
Martin Reinecke committed
60
61
62
63
64
65
66
67
    def __init__(self, dirty_domain, grid_domain, gconf):
        import nifty_gridder
        self._domain = makeDomain(grid_domain)
        self._target = makeDomain(dirty_domain)
        self._gconf = gconf
        fu = gconf.U_corrections()
        fv = gconf.V_corrections()
        nu, nv = dirty_domain.shape
Philipp Arras's avatar
Philipp Arras committed
68
        # compute deconvolution operator
Martin Reinecke's avatar
fixes    
Martin Reinecke committed
69
70
71
72
73
74
        rng = np.arange(nu)
        k = np.minimum(rng, nu-rng)
        self._deconv_u = np.roll(fu[k], -nu//2).reshape((-1, 1))
        rng = np.arange(nv)
        k = np.minimum(rng, nv-rng)
        self._deconv_v = np.roll(fv[k], -nv//2).reshape((1, -1))
Philipp Arras's avatar
Philipp Arras committed
75
76
77
        self._capability = self.TIMES | self.ADJOINT_TIMES

    def apply(self, x, mode):
Martin Reinecke's avatar
Martin Reinecke committed
78
        import nifty_gridder
Philipp Arras's avatar
Philipp Arras committed
79
80
        self._check_input(x, mode)
        nu, nv = self._target.shape
81
        res = x.to_global_data()
Philipp Arras's avatar
Philipp Arras committed
82
        if mode == self.TIMES:
83
            res = hartley(res)
Philipp Arras's avatar
Philipp Arras committed
84
85
            res = np.roll(res, (nu//2, nv//2), axis=(0, 1))
            res = res[:nu, :nv]
Martin Reinecke's avatar
fixes    
Martin Reinecke committed
86
87
            res *= self._deconv_u
            res *= self._deconv_v
Philipp Arras's avatar
Philipp Arras committed
88
        else:
Martin Reinecke's avatar
fixes    
Martin Reinecke committed
89
90
            res = res*self._deconv_u
            res *= self._deconv_v
Philipp Arras's avatar
Philipp Arras committed
91
            nu2, nv2 = self._domain.shape
92
            res = np.pad(res, ((0, nu2-nu), (0, nv2-nv)), mode='constant',
Philipp Arras's avatar
Philipp Arras committed
93
94
                         constant_values=0)
            res = np.roll(res, (-nu//2, -nv//2), axis=(0, 1))
95
            res = hartley(res)
Philipp Arras's avatar
Philipp Arras committed
96
97
98
99
        return from_global_data(self._tgt(mode), res)


class RadioGridder(LinearOperator):
Martin Reinecke's avatar
Martin Reinecke committed
100
101
102
103
104
105
    def __init__(self, grid_domain, bl, gconf, idx):
        self._domain = DomainTuple.make(UnstructuredDomain((idx.shape[0],)))
        self._target = DomainTuple.make(grid_domain)
        self._bl = bl
        self._gconf = gconf
        self._idx = idx
Philipp Arras's avatar
Philipp Arras committed
106
107
108
        self._capability = self.TIMES | self.ADJOINT_TIMES

    def apply(self, x, mode):
Martin Reinecke's avatar
Martin Reinecke committed
109
        import nifty_gridder
Philipp Arras's avatar
Philipp Arras committed
110
111
        self._check_input(x, mode)
        if mode == self.TIMES:
Martin Reinecke's avatar
Martin Reinecke committed
112
113
            res = nifty_gridder.ms2grid(
                self._bl, self._gconf, self._idx, x.to_global_data().reshape((-1,1)))
Philipp Arras's avatar
Philipp Arras committed
114
        else:
Martin Reinecke's avatar
Martin Reinecke committed
115
116
            res = nifty_gridder.grid2ms(
                self._bl, self._gconf, self._idx, x.to_global_data())
Philipp Arras's avatar
Philipp Arras committed
117
        return from_global_data(self._tgt(mode), res)