convolution_operators.py 3.81 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# 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) 2013-2019 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.


19
20
from .. import utilities
from ..domain_tuple import DomainTuple
21
from ..domains.gl_space import GLSpace
22
23
24
from ..domains.hp_space import HPSpace
from ..domains.rg_space import RGSpace
from .diagonal_operator import DiagonalOperator
25
26
from .endomorphic_operator import EndomorphicOperator
from .harmonic_operators import HarmonicTransformOperator
27
from .simple_linear_operators import WeightApplier
28
29


30
def FuncConvolutionOperator(domain, func, space=None):
31
32
33
34
    """Convolves input with a radially symmetric kernel defined by `func`

    Parameters
    ----------
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
35
    domain: DomainTuple
36
        Domain of the operator.
Martin Reinecke's avatar
Martin Reinecke committed
37
38
39
40
    func: function
        This function needs to take exactly one argument, which is
        colatitude in radians, and return the kernel amplitude at that
        colatitude.
41
42
43
44
    space: int, optional
        The index of the subdomain on which the operator should act
        If None, it is set to 0 if `domain` contains exactly one space.
        `domain[space]` must be of type `RGSpace`, `HPSpace`, or `GLSpace`.
45
46
47

    Notes
    -----
Martin Reinecke's avatar
tweaks    
Martin Reinecke committed
48
49
50
51
    The operator assumes periodic boundaries in the input domain. This means
    for a sufficiently broad function a point source close to the boundary will
    blur into the opposite side of the image. Zero padding can be applied to
    avoid this behaviour.
52
    """
53
54
55
56
57
58
    domain = DomainTuple.make(domain)
    space = utilities.infer_space(domain, space)
    if not isinstance(domain[space], (RGSpace, HPSpace, GLSpace)):
        raise TypeError("unsupported domain")
    codomain = domain[space].get_default_codomain()
    kernel = codomain.get_conv_kernel_from_func(func)
59
    return _ConvolutionOperator(domain, kernel, space)
60
61


62
def _ConvolutionOperator(domain, kernel, space=None):
63
64
    domain = DomainTuple.make(domain)
    space = utilities.infer_space(domain, space)
Martin Reinecke's avatar
Martin Reinecke committed
65
    if len(kernel.target) != 1:
66
67
68
69
70
71
        raise ValueError("kernel needs exactly one domain")
    if not isinstance(domain[space], (HPSpace, GLSpace, RGSpace)):
        raise TypeError("need RGSpace, HPSpace, or GLSpace")
    lm = [d for d in domain]
    lm[space] = lm[space].get_default_codomain()
    lm = DomainTuple.make(lm)
Martin Reinecke's avatar
Martin Reinecke committed
72
    if lm[space] != kernel.target[0]:
73
74
75
76
        raise ValueError("Input domain and kernel are incompatible")
    HT = HarmonicTransformOperator(lm, domain[space], space)
    diag = DiagonalOperator(kernel*domain[space].total_volume, lm, (space,))
    wgt = WeightApplier(domain, space, 1)
77
    op = HT(diag(HT.adjoint(wgt)))
78
    return _ApplicationWithoutMeanOperator(op)
79
80
81


class _ApplicationWithoutMeanOperator(EndomorphicOperator):
Lukas Platz's avatar
cleanup    
Lukas Platz committed
82
    def __init__(self, op):
83
        self._capability = self.TIMES | self.ADJOINT_TIMES
Lukas Platz's avatar
cleanup    
Lukas Platz committed
84
85
86
        if op.domain != op.target:
            raise TypeError("Operator needs to be endomorphic")
        self._domain = op.domain
87
88
89
90
        self._op = op

    def apply(self, x, mode):
        self._check_input(x, mode)
Martin Reinecke's avatar
Martin Reinecke committed
91
        mean = x.s_mean()
92
93
94
95
96
97
98
        return mean + self._op.apply(x - mean, mode)

    def __repr__(self):
        from ..utilities import indent
        return "\n".join((
            "_ApplicationWithoutMeanOperator:",
            indent(self._op.__repr__())))