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.
|