gridder.py 4.52 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
29
30
31
32
33
34
35
36
from ..operators.linear_operator import LinearOperator
from ..sugar import from_global_data, makeDomain


class GridderMaker(object):
    def __init__(self, domain, eps=1e-15):
        domain = makeDomain(domain)
        if (len(domain) != 1 or not isinstance(domain[0], RGSpace) or
                not len(domain.shape) == 2):
            raise ValueError("need domain with exactly one 2D RGSpace")
        nu, nv = domain.shape
        if nu % 2 != 0 or nv % 2 != 0:
            raise ValueError("dimensions must be even")
37
38
39
40
41
        nu2, nv2 = 2*nu, 2*nv
        w = int(-np.log10(eps)+1.9999)
        nsafe = (w+1)//2
        nu2 = max([nu2, 2*nsafe])
        nv2 = max([nv2, 2*nsafe])
Philipp Arras's avatar
Philipp Arras committed
42
43
44
45

        oversampled_domain = RGSpace(
            [nu2, nv2], distances=[1, 1], harmonic=False)

46
47
        self._eps = eps
        self._rest = _RestOperator(domain, oversampled_domain, eps)
Philipp Arras's avatar
Philipp Arras committed
48
49

    def getReordering(self, uv):
50
        from nifty_gridder import peanoindex
Philipp Arras's avatar
Philipp Arras committed
51
52
53
54
        nu2, nv2 = self._rest._domain.shape
        return peanoindex(uv, nu2, nv2)

    def getGridder(self, uv):
55
        return RadioGridder(self._rest.domain, self._eps, uv)
Philipp Arras's avatar
Philipp Arras committed
56
57
58
59
60
61
62
63
64

    def getRest(self):
        return self._rest

    def getFull(self, uv):
        return self.getRest() @ self.getGridder(uv)


class _RestOperator(LinearOperator):
65
    def __init__(self, domain, oversampled_domain, eps):
Martin Reinecke's avatar
fixes    
Martin Reinecke committed
66
        from nifty_gridder import correction_factors
Philipp Arras's avatar
Philipp Arras committed
67
68
69
70
71
        self._domain = makeDomain(oversampled_domain)
        self._target = domain
        nu, nv = domain.shape
        nu2, nv2 = oversampled_domain.shape

Martin Reinecke's avatar
fixes    
Martin Reinecke committed
72
73
        fu = correction_factors(nu2, nu//2+1, eps)
        fv = correction_factors(nv2, nv//2+1, eps)
Philipp Arras's avatar
Philipp Arras committed
74
        # compute deconvolution operator
Martin Reinecke's avatar
fixes    
Martin Reinecke committed
75
76
77
78
79
80
        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
81
82
83
84
85
        self._capability = self.TIMES | self.ADJOINT_TIMES

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


class RadioGridder(LinearOperator):
105
    def __init__(self, target, eps, uv):
Philipp Arras's avatar
Philipp Arras committed
106
107
108
109
        self._domain = DomainTuple.make(
            UnstructuredDomain((uv.shape[0],)))
        self._target = DomainTuple.make(target)
        self._capability = self.TIMES | self.ADJOINT_TIMES
110
        self._eps = float(eps)
Philipp Arras's avatar
Philipp Arras committed
111
112
113
        self._uv = uv  # FIXME: should we write-protect this?

    def apply(self, x, mode):
Martin Reinecke's avatar
Martin Reinecke committed
114
115
        from nifty_gridder import (to_grid, to_grid_post,
                                   from_grid, from_grid_pre)
Philipp Arras's avatar
Philipp Arras committed
116
117
118
119
        self._check_input(x, mode)
        nu2, nv2 = self._target.shape
        x = x.to_global_data()
        if mode == self.TIMES:
120
            res = to_grid(self._uv, x, nu2, nv2, self._eps)
Martin Reinecke's avatar
Martin Reinecke committed
121
            res = to_grid_post(res)
Philipp Arras's avatar
Philipp Arras committed
122
        else:
Martin Reinecke's avatar
Martin Reinecke committed
123
            x = from_grid_pre(x)
124
            res = from_grid(self._uv, x, nu2, nv2, self._eps)
Philipp Arras's avatar
Philipp Arras committed
125
        return from_global_data(self._tgt(mode), res)