File: kernels.rst

package info (click to toggle)
astropy 5.2.1-2
  • links: PTS, VCS
  • area: main
  • in suites: bookworm
  • size: 41,972 kB
  • sloc: python: 219,331; ansic: 147,297; javascript: 13,556; lex: 8,496; sh: 3,319; xml: 1,622; makefile: 185
file content (387 lines) | stat: -rw-r--r-- 12,531 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
Convolution Kernels
*******************

Introduction and Concept
========================

The convolution module provides several built-in kernels to cover the most
common applications in astronomy. It is also possible to define custom kernels
from arrays or combine existing kernels to match specific applications.

Every filter kernel is characterized by its response function. For time series
we speak of an "impulse response function" or for images we call it "point
spread function." This response function is given for every kernel by a
`~astropy.modeling.FittableModel`, which is evaluated on a grid with
:func:`~astropy.convolution.discretize_model` to obtain a kernel
array, which can be used for discrete convolution with the binned data.


Examples
========

1D Kernels
----------

..
  EXAMPLE START
  Using 1D Kernels to Smooth Noisy Data

One application of filtering is to smooth noisy data. In this case we
consider a noisy Lorentz curve:

>>> import numpy as np
>>> from astropy.modeling.models import Lorentz1D
>>> from astropy.convolution import convolve, Gaussian1DKernel, Box1DKernel
>>> rng = np.random.default_rng()
>>> lorentz = Lorentz1D(1, 0, 1)
>>> x = np.linspace(-5, 5, 100)
>>> data_1D = lorentz(x) + 0.1 * (rng.random(100) - 0.5)

Smoothing the noisy data with a `~astropy.convolution.Gaussian1DKernel`
with a standard deviation of 2 pixels:

>>> gauss_kernel = Gaussian1DKernel(2)
>>> smoothed_data_gauss = convolve(data_1D, gauss_kernel)

Smoothing the same data with a `~astropy.convolution.Box1DKernel` of width 5
pixels:

>>> box_kernel = Box1DKernel(5)
>>> smoothed_data_box = convolve(data_1D, box_kernel)

The following plot illustrates the results:

.. plot::

    import numpy as np
    import matplotlib.pyplot as plt
    from astropy.modeling.models import Lorentz1D
    from astropy.convolution import convolve, Gaussian1DKernel, Box1DKernel

    # Fake Lorentz data including noise
    rng = np.random.default_rng()
    lorentz = Lorentz1D(1, 0, 1)
    x = np.linspace(-5, 5, 100)
    data_1D = lorentz(x) + 0.1 * (rng.random(100) - 0.5)

    # Smooth data
    gauss_kernel = Gaussian1DKernel(2)
    smoothed_data_gauss = convolve(data_1D, gauss_kernel)
    box_kernel = Box1DKernel(5)
    smoothed_data_box = convolve(data_1D, box_kernel)

    # Plot data and smoothed data
    plt.plot(x, data_1D, label='Original')
    plt.plot(x, smoothed_data_gauss, label='Smoothed with Gaussian1DKernel')
    plt.plot(x, smoothed_data_box, label='Smoothed with Box1DKernel')
    plt.xlabel('x [a.u.]')
    plt.ylabel('amplitude [a.u.]')
    plt.xlim(-5, 5)
    plt.ylim(-0.1, 1.5)
    plt.legend(prop={'size':12})
    plt.show()


Beside the ``astropy`` convolution functions
`~astropy.convolution.convolve` and
`~astropy.convolution.convolve_fft`, it is also possible to use
the kernels with ``numpy`` or ``scipy`` convolution by passing the ``array``
attribute. This will be faster in most cases than the ``astropy`` convolution,
but will not work properly if NaN values are present in the data.

>>> smoothed = np.convolve(data_1D, box_kernel.array)

..
  EXAMPLE END

2D Kernels
----------

..
  EXAMPLE START
  Using 2D Kernels to Smooth Noisy Data

As all 2D kernels are symmetric, it is sufficient to specify the width in one
direction. Therefore the use of 2D kernels is basically the same as for 1D
kernels. Here we consider a small Gaussian-shaped source of amplitude 1 in the
middle of the image and add 10% noise:

>>> import numpy as np
>>> from astropy.convolution import convolve, Gaussian2DKernel, Tophat2DKernel
>>> from astropy.modeling.models import Gaussian2D
>>> gauss = Gaussian2D(1, 0, 0, 3, 3)
>>> # Fake image data including noise
>>> rng = np.random.default_rng()
>>> x = np.arange(-100, 101)
>>> y = np.arange(-100, 101)
>>> x, y = np.meshgrid(x, y)
>>> data_2D = gauss(x, y) + 0.1 * (rng.random((201, 201)) - 0.5)

Smoothing the noisy data with a
:class:`~astropy.convolution.Gaussian2DKernel` with a standard
deviation of 2 pixels:

>>> gauss_kernel = Gaussian2DKernel(2)
>>> smoothed_data_gauss = convolve(data_2D, gauss_kernel)

Smoothing the noisy data with a
:class:`~astropy.convolution.Tophat2DKernel` of width 5 pixels:

>>> tophat_kernel = Tophat2DKernel(5)
>>> smoothed_data_tophat = convolve(data_2D, tophat_kernel)

This is what the original image looks like:

.. plot::

    import numpy as np
    import matplotlib.pyplot as plt
    from astropy.modeling.models import Gaussian2D
    gauss = Gaussian2D(1, 0, 0, 2, 2)
    # Fake image data including noise
    rng = np.random.default_rng()
    x = np.arange(-100, 101)
    y = np.arange(-100, 101)
    x, y = np.meshgrid(x, y)
    data_2D = gauss(x, y) + 0.1 * (rng.random((201, 201)) - 0.5)
    plt.imshow(data_2D, origin='lower')
    plt.xlabel('x [pixels]')
    plt.ylabel('y [pixels]')
    plt.colorbar()
    plt.show()

The following plot illustrates the differences between several 2D kernels
applied to the simulated data. Note that it has a slightly different color
scale compared to the original image.

.. plot::

    import numpy as np
    import matplotlib.pyplot as plt

    from astropy.convolution import *
    from astropy.modeling.models import Gaussian2D

    # Small Gaussian source in the middle of the image
    gauss = Gaussian2D(1, 0, 0, 2, 2)
    # Fake data including noise
    rng = np.random.default_rng()
    x = np.arange(-100, 101)
    y = np.arange(-100, 101)
    x, y = np.meshgrid(x, y)
    data_2D = gauss(x, y) + 0.1 * (rng.random((201, 201)) - 0.5)

    # Setup kernels, including unity kernel for original image
    # Choose normalization for linear scale space for RickerWavelet

    kernels = [TrapezoidDisk2DKernel(11, slope=0.2),
               Tophat2DKernel(11),
               Gaussian2DKernel(11),
               Box2DKernel(11),
               11 ** 2 * RickerWavelet2DKernel(11),
               AiryDisk2DKernel(11)]

    fig, axes = plt.subplots(nrows=2, ncols=3)

    # Plot kernels
    for kernel, ax in zip(kernels, axes.flat):
        smoothed = convolve(data_2D, kernel, normalize_kernel=False)
        im = ax.imshow(smoothed, vmin=-0.01, vmax=0.08, origin='lower',
                       interpolation='None')
        title = kernel.__class__.__name__
        ax.set_title(title, fontsize=12)
        ax.set_yticklabels([])
        ax.set_xticklabels([])

    cax = fig.add_axes([0.9, 0.1, 0.03, 0.8])
    fig.colorbar(im, cax=cax)
    plt.subplots_adjust(left=0.05, right=0.85, top=0.95, bottom=0.05)
    plt.show()


The Gaussian kernel has better smoothing properties compared to the Box and the
Top Hat. The Box filter is not isotropic and can produce artifacts (the source
appears rectangular). The Ricker Wavelet filter removes noise and slowly varying
structures (i.e., background), but produces a negative ring around the source.
The best choice for the filter strongly depends on the application.

..
  EXAMPLE END

Available Kernels
=================

.. currentmodule:: astropy.convolution

.. autosummary::

   AiryDisk2DKernel
   Box1DKernel
   Box2DKernel
   CustomKernel
   Gaussian1DKernel
   Gaussian2DKernel
   RickerWavelet1DKernel
   RickerWavelet2DKernel
   Model1DKernel
   Model2DKernel
   Moffat2DKernel
   Ring2DKernel
   Tophat2DKernel
   Trapezoid1DKernel
   TrapezoidDisk2DKernel

Kernel Arithmetics
==================

Addition and Subtraction
------------------------

As convolution is a linear operation, kernels can be added or subtracted from
each other. They can also be multiplied with some number.

Examples
^^^^^^^^

..
  EXAMPLE START
  Adding and Subtracting Kernels in astropy.convolution

One basic example of subtracting kernels would be the definition of a
Difference of Gaussian filter:

>>> from astropy.convolution import Gaussian1DKernel
>>> gauss_1 = Gaussian1DKernel(10)
>>> gauss_2 = Gaussian1DKernel(16)
>>> DoG = gauss_2 - gauss_1

Another application is to convolve faked data with an instrument response
function model. For example, if the response function can be described by
the weighted sum of two Gaussians:

>>> gauss_1 = Gaussian1DKernel(10)
>>> gauss_2 = Gaussian1DKernel(16)
>>> SoG = 4 * gauss_1 + gauss_2

Most times it will be necessary to normalize the resulting kernel by calling
explicitly:

>>> SoG.normalize()

..
  EXAMPLE END

Convolution
-----------

Furthermore, two kernels can be convolved with each other, which is useful when
data is filtered with two different kinds of kernels or to create a new,
special kernel.

Examples
^^^^^^^^

..
  EXAMPLE START
  Convolving Kernels in astropy.convolution

To convolve two kernels with each other:

>>> from astropy.convolution import Gaussian1DKernel, convolve
>>> gauss_1 = Gaussian1DKernel(10)
>>> gauss_2 = Gaussian1DKernel(16)
>>> broad_gaussian = convolve(gauss_2,  gauss_1)  # doctest: +IGNORE_WARNINGS

Or in case of multistage smoothing:

>>> import numpy as np
>>> from astropy.modeling.models import Lorentz1D
>>> from astropy.convolution import convolve, Gaussian1DKernel, Box1DKernel
>>> rng = np.random.default_rng()
>>> lorentz = Lorentz1D(1, 0, 1)
>>> x = np.linspace(-5, 5, 100)
>>> data_1D = lorentz(x) + 0.1 * (rng.random(100) - 0.5)

>>> gauss = Gaussian1DKernel(3)
>>> box = Box1DKernel(5)
>>> smoothed_gauss = convolve(data_1D, gauss)
>>> smoothed_gauss_box = convolve(smoothed_gauss, box)

You would rather do the following:

>>> gauss = Gaussian1DKernel(3)
>>> box = Box1DKernel(5)
>>> smoothed_gauss_box = convolve(data_1D, convolve(box, gauss))  # doctest: +IGNORE_WARNINGS

Which, in most cases, will also be faster than the first method because only
one convolution with the often times larger data array will be necessary.

..
  EXAMPLE END

Discretization
==============

To obtain the kernel array for discrete convolution, the kernel's response
function is evaluated on a grid with
:func:`~astropy.convolution.discretize_model`. For the
discretization step the following modes are available:

* Mode ``'center'`` (default) evaluates the response function on the grid by
  taking the value at the center of the bin.

   >>> from astropy.convolution import Gaussian1DKernel
   >>> gauss_center = Gaussian1DKernel(3, mode='center')

* Mode ``'linear_interp'`` takes the values at the corners of the bin and
  linearly interpolates the value at the center:

  >>> gauss_interp = Gaussian1DKernel(3, mode='linear_interp')

* Mode ``'oversample'`` evaluates the response function by taking the mean on an
  oversampled grid. The oversample factor can be specified with the ``factor``
  argument. If the oversample factor is too large, the evaluation becomes slow.

 >>> gauss_oversample = Gaussian1DKernel(3, mode='oversample', factor=10)

* Mode ``'integrate'`` integrates the function over the pixel using
  ``scipy.integrate.quad`` and ``scipy.integrate.dblquad``. This mode is very
  slow and is only recommended when the highest accuracy is required.

.. doctest-requires:: scipy

    >>> gauss_integrate = Gaussian1DKernel(3, mode='integrate')

Especially in the range where the kernel width is in order of only a few pixels,
it can be advantageous to use the mode ``oversample`` or ``integrate`` to
conserve the integral on a subpixel scale.

.. _kernel_normalization:

Normalization
=============

The kernel models are normalized per default (i.e.,
:math:`\int_{-\infty}^{\infty} f(x) dx = 1`). But because of the limited kernel
array size, the normalization for kernels with an infinite response can differ
from one. The value of this deviation is stored in the kernel's ``truncation``
attribute.

The normalization can also differ from one, especially for small kernels, due
to the discretization step. This can be partly controlled by the ``mode``
argument, when initializing the kernel. (See also
:func:`~astropy.convolution.discretize_model`.) Setting the
``mode`` to ``'oversample'`` allows us to conserve the normalization even on the
subpixel scale.

The kernel arrays can be renormalized explicitly by calling either the
``normalize()`` method or by setting the ``normalize_kernel`` argument in the
:func:`~astropy.convolution.convolve` and
:func:`~astropy.convolution.convolve_fft` functions. The latter
method leaves the kernel itself unchanged but works with an internal normalized
version of the kernel.

Note that for :class:`~astropy.convolution.RickerWavelet1DKernel`
and :class:`~astropy.convolution.RickerWavelet2DKernel` there is
:math:`\int_{-\infty}^{\infty} f(x) dx = 0`. To define a proper normalization,
both filters are derived from a normalized Gaussian function.