convolution.py 1.76 KB
Newer Older
1
import numpy as np
Martin Reinecke's avatar
5->6  
Martin Reinecke committed
2
import nifty6 as ift
3
4
5
6
7
8
9
10
11


def convtest(test_signal, delta, func):
    domain = test_signal.domain

    # Create Convolution Operator
    conv_op = ift.FuncConvolutionOperator(domain, func)

    # Convolve, Adjoint-Convolve
12
    conv_signal = conv_op(test_signal)
13
14
    cac_signal = conv_op.adjoint_times(conv_signal)

15
16
    print(test_signal.integrate(), conv_signal.integrate(),
          cac_signal.integrate())
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

    # generate kernel image
    conv_delta = conv_op(delta)

    # Plot results
    plot = ift.Plot()
    plot.add(signal, title='Signal')
    plot.add(conv_signal, title='Signal Convolved')
    plot.add(cac_signal, title='Signal, Conv, Adj-Conv')
    plot.add(conv_delta, title='Kernel')
    plot.output()


# Healpix test
nside = 64
npix = 12 * nside * nside

domain = ift.HPSpace(nside)

# Define test signal (some point sources)
signal_vals = np.zeros(npix, dtype=np.float64)
for i in range(0, npix, npix//12 + 27):
39
    signal_vals[i] = 500.
40

Martin Reinecke's avatar
Martin Reinecke committed
41
signal = ift.makeField(domain, signal_vals)
42
43
44

delta_vals = np.zeros(npix, dtype=np.float64)
delta_vals[0] = 1.0
Martin Reinecke's avatar
Martin Reinecke committed
45
delta = ift.makeField(domain, delta_vals)
46
47
48
49
50
51
52
53
54
55
56
57
58


# Define kernel function
def func(theta):
    ct = np.cos(theta)
    return 1. * np.logical_and(ct > 0.7, ct <= 0.8)


convtest(signal, delta, func)

domain = ift.RGSpace((100, 100))
# Define test signal (some point sources)
signal_vals = np.zeros(domain.shape, dtype=np.float64)
59
60
signal_vals[35, 70] = 5000.
signal_vals[45, 8] = 5000.
Martin Reinecke's avatar
Martin Reinecke committed
61
signal = ift.makeField(domain, signal_vals)
62
63
64
65

# Define delta signal, generate kernel image
delta_vals = np.zeros(domain.shape, dtype=np.float64)
delta_vals[0, 0] = 1.0
Martin Reinecke's avatar
Martin Reinecke committed
66
delta = ift.makeField(domain, delta_vals)
67
68
69
70
71
72
73
74


# Define kernel function
def func(dist):
    return 1. * np.logical_and(dist > 0.1, dist <= 0.2)


convtest(signal, delta, func)