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
|
************************************
Convolving with Unnormalized Kernels
************************************
There are some tasks, such as source finding, where you want to apply a filter
with a kernel that is not normalized.
For data that are well-behaved (contain no missing or infinite values), this
can be done in one step::
convolve(image, kernel)
Examples
--------
..
EXAMPLE START
Convolving with Unnormalized Kernels
For an example of applying a filter with a kernel that is not normalized, we
can try to run a commonly used peak enhancing kernel:
.. plot::
:context: reset
:include-source:
:align: center
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.utils.data import get_pkg_data_filename
from astropy.convolution import CustomKernel
from scipy.signal import convolve as scipy_convolve
from astropy.convolution import convolve, convolve_fft
# Load the data from data.astropy.org
filename = get_pkg_data_filename('galactic_center/gc_msx_e.fits')
hdu = fits.open(filename)[0]
# Scale the file to have reasonable numbers
# (this is mostly so that colorbars don't have too many digits)
# Also, we crop it so you can see individual pixels
img = hdu.data[50:90, 60:100] * 1e5
kernel = CustomKernel([[-1,-1,-1], [-1, 8, -1], [-1,-1,-1]])
astropy_conv = convolve(img, kernel, normalize_kernel=False, nan_treatment='fill')
#astropy_conv_fft = convolve_fft(img, kernel, normalize_kernel=False, nan_treatment='fill')
plt.figure(1, figsize=(12, 12)).clf()
ax1 = plt.subplot(1, 2, 1)
im = ax1.imshow(img, vmin=-6., vmax=5.e1, origin='lower',
interpolation='nearest', cmap='viridis')
ax2 = plt.subplot(1, 2, 2)
im = ax2.imshow(astropy_conv, vmin=-6., vmax=5.e1, origin='lower',
interpolation='nearest', cmap='viridis')
..
EXAMPLE END
..
EXAMPLE START
Replacing NaN Values with Interpolated Values Using Kernels
If you have an image with missing values (NaNs), you have to replace them with
real values first. Often, the best way to do this is to replace the NaN values
with interpolated values. In the example below, we use a Gaussian kernel
with a size similar to that of our peak-finding kernel to replace the bad data
before applying the peak-finding kernel.
.. plot::
:context:
:include-source:
:align: center
from astropy.convolution import Gaussian2DKernel, interpolate_replace_nans
# Select a random set of pixels that were affected by some sort of artifact
# and replaced with NaNs (e.g., cosmic-ray-affected pixels)
rng = np.random.default_rng(42)
yinds, xinds = np.indices(img.shape)
img[rng.choice(yinds.flat, 50), rng.choice(xinds.flat, 50)] = np.nan
# We smooth with a Gaussian kernel with x_stddev=1 (and y_stddev=1)
# It is a 9x9 array
kernel = Gaussian2DKernel(x_stddev=1)
# interpolate away the NaNs
reconstructed_image = interpolate_replace_nans(img, kernel)
# apply peak-finding
kernel = CustomKernel([[-1,-1,-1], [-1, 8, -1], [-1,-1,-1]])
# Use the peak-finding kernel
# We have to turn off kernel normalization and set nan_treatment to "fill"
# here because `nan_treatment='interpolate'` is incompatible with non-
# normalized kernels
peaked_image = convolve(reconstructed_image, kernel,
normalize_kernel=False,
nan_treatment='fill')
plt.figure(1, figsize=(12, 12)).clf()
ax1 = plt.subplot(1, 3, 1)
ax1.set_title("Image with missing data")
im = ax1.imshow(img, vmin=-6., vmax=5.e1, origin='lower',
interpolation='nearest', cmap='viridis')
ax2 = plt.subplot(1, 3, 2)
ax2.set_title("Interpolated")
im = ax2.imshow(reconstructed_image, vmin=-6., vmax=5.e1, origin='lower',
interpolation='nearest', cmap='viridis')
ax3 = plt.subplot(1, 3, 3)
ax3.set_title("Peak-Finding")
im = ax3.imshow(peaked_image, vmin=-6., vmax=5.e1, origin='lower',
interpolation='nearest', cmap='viridis')
..
EXAMPLE END
|
