卷积和滤波 (astropy.convolution

介绍

astropy.convolution 提供卷积函数和内核,与SciPy相比有改进 scipy.ndimage 卷积程序,包括:

  • 正确处理NaN值(在卷积期间忽略它们,并用插值值替换NaN像素)

  • 用于一维、二维和三维卷积的单个函数

  • 改进了边缘处理的选项

  • 直接和快速傅里叶变换(FFT)版本

  • 天文学中常用的内置内核

下面的缩略图显示了 scipyastropy 包含NaN值的天文图像上的卷积函数。 scipy 的函数本质上是为任何NaN值的内核内的所有像素返回NaN,这通常不是期望的结果。

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 Gaussian2DKernel
from scipy.signal import convolve as scipy_convolve
from astropy.convolution import convolve


# 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 do not have too many digits)
# Also, we crop it so you can see individual pixels
img = hdu.data[50:90, 60:100] * 1e5

# This example is intended to demonstrate how astropy.convolve and
# scipy.convolve handle missing data, so we start by setting the
# brightest pixels to NaN to simulate a "saturated" data set
img[img > 2e1] = np.nan

# We also create a copy of the data and set those NaNs to zero.  We will
# use this for the scipy convolution
img_zerod = img.copy()
img_zerod[np.isnan(img)] = 0

# We smooth with a Gaussian kernel with x_stddev=1 (and y_stddev=1)
# It is a 9x9 array
kernel = Gaussian2DKernel(x_stddev=1)

# Convolution: scipy's direct convolution mode spreads out NaNs (see
# panel 2 below)
scipy_conv = scipy_convolve(img, kernel, mode='same', method='direct')

# scipy's direct convolution mode run on the 'zero'd' image will not
# have NaNs, but will have some very low value zones where the NaNs were
# (see panel 3 below)
scipy_conv_zerod = scipy_convolve(img_zerod, kernel, mode='same',
                                  method='direct')

# astropy's convolution replaces the NaN pixels with a kernel-weighted
# interpolation from their neighbors
astropy_conv = convolve(img, kernel)


# Now we do a bunch of plots.  In the first two plots, the originally masked
# values are marked with red X's
plt.figure(1, figsize=(12, 12)).clf()
ax1 = plt.subplot(2, 2, 1)
im = ax1.imshow(img, vmin=-2., vmax=2.e1, origin='lower',
                interpolation='nearest', cmap='viridis')
y, x = np.where(np.isnan(img))
ax1.set_autoscale_on(False)
ax1.plot(x, y, 'rx', markersize=4)
ax1.set_title("Original")
ax1.set_xticklabels([])
ax1.set_yticklabels([])

ax2 = plt.subplot(2, 2, 2)
im = ax2.imshow(scipy_conv, vmin=-2., vmax=2.e1, origin='lower',
                interpolation='nearest', cmap='viridis')
ax2.set_autoscale_on(False)
ax2.plot(x, y, 'rx', markersize=4)
ax2.set_title("Scipy")
ax2.set_xticklabels([])
ax2.set_yticklabels([])

ax3 = plt.subplot(2, 2, 3)
im = ax3.imshow(scipy_conv_zerod, vmin=-2., vmax=2.e1, origin='lower',
                interpolation='nearest', cmap='viridis')
ax3.set_title("Scipy nan->zero")
ax3.set_xticklabels([])
ax3.set_yticklabels([])

ax4 = plt.subplot(2, 2, 4)
im = ax4.imshow(astropy_conv, vmin=-2., vmax=2.e1, origin='lower',
                interpolation='nearest', cmap='viridis')
ax4.set_title("Default astropy")
ax4.set_xticklabels([])
ax4.set_yticklabels([])

# we make a second plot of the amplitudes vs offset position to more
# clearly illustrate the value differences
plt.figure(2).clf()
plt.plot(img[:, 25], label='input', drawstyle='steps-mid', linewidth=2,
         alpha=0.5)
plt.plot(scipy_conv[:, 25], label='scipy', drawstyle='steps-mid',
         linewidth=2, alpha=0.5, marker='s')
plt.plot(scipy_conv_zerod[:, 25], label='scipy nan->zero',
         drawstyle='steps-mid', linewidth=2, alpha=0.5, marker='s')
plt.plot(astropy_conv[:, 25], label='astropy', drawstyle='steps-mid',
         linewidth=2, alpha=0.5)
plt.ylabel("Amplitude")
plt.ylabel("Position Offset")
plt.legend(loc='best')
plt.show()
../_images/index-1_00.png

(png _, svgpdf

../_images/index-1_01.png

(png _, svgpdf

以下各节介绍如何使用卷积函数,以及如何使用内置卷积内核:

入门

提供两个卷积函数。它们被导入为:

from astropy.convolution import convolve, convolve_fft

它们都用作:

result = convolve(image, kernel)
result = convolve_fft(image, kernel)

convolve() 实现为直接卷积算法,而 convolve_fft() 使用快速傅立叶变换(FFT)。因此,前者更适合于小内核,而后者对于较大的内核更有效。

例子

要使用用户指定的内核对1D数据集进行卷积,可以执行以下操作:

>>> from astropy.convolution import convolve
>>> convolve([1, 4, 5, 6, 5, 7, 8], [0.2, 0.6, 0.2])  
array([1.4, 3.6, 5. , 5.6, 5.6, 6.8, 6.2])

请注意,端点设置为零-默认情况下,过于接近边界而无法计算卷积值的点将设置为零。但是 convolve() 函数允许 boundary 可用于指定替代行为的参数。例如,设置 boundary='extend' 使计算边缘附近的值,假设原始数据仅使用边界外的常量外推进行扩展:

>>> from astropy.convolution import convolve
>>> convolve([1, 4, 5, 6, 5, 7, 8], [0.2, 0.6, 0.2], boundary='extend')  
array([1.6, 3.6, 5. , 5.6, 5.6, 6.8, 7.8])

在计算末尾的值时,假设第一个点以下的任何值都是 1 ,任何高于最后一点的值都是 8 . 有关边界处理的详细讨论,请参见 使用卷积函数 .

例子

卷积模块还包括可以导入的内置内核,例如:

>>> from astropy.convolution import Gaussian1DKernel

要使用内核,首先创建一个特定的内核实例:

>>> gauss = Gaussian1DKernel(stddev=2)

gauss 不是数组,而是内核对象。可以使用以下命令检索基础数组:

>>> gauss.array  
array([6.69151129e-05, 4.36341348e-04, 2.21592421e-03,
       8.76415025e-03, 2.69954833e-02, 6.47587978e-02,
       1.20985362e-01, 1.76032663e-01, 1.99471140e-01,
       1.76032663e-01, 1.20985362e-01, 6.47587978e-02,
       2.69954833e-02, 8.76415025e-03, 2.21592421e-03,
       4.36341348e-04, 6.69151129e-05])

然后可以在调用时直接使用内核 convolve()

import numpy as np
import matplotlib.pyplot as plt

from astropy.convolution import Gaussian1DKernel, convolve

plt.figure(3).clf()

# Generate fake data
x = np.arange(1000).astype(float)
y = np.sin(x / 100.) + np.random.normal(0., 1., x.shape)
y[::3] = np.nan

# Create kernel
g = Gaussian1DKernel(stddev=50)

# Convolve data
z = convolve(y, g)

# Plot data before and after convolution
plt.plot(x, y, 'k-', label='Before')
plt.plot(x, z, 'b-', label='After', alpha=0.5, linewidth=2)
plt.legend(loc='best')
plt.show()

(png _, svgpdf

../_images/index-2.png

使用 astropy 的卷积来替换坏数据

astropy 的卷积方法可以用来用从相邻数据中插值的值替换坏数据。基于核的插值对于处理含有少量坏像素的图像或对稀疏采样的图像进行插值非常有用。

插值工具的实现和用途如下:

from astropy.convolution import interpolate_replace_nans
result = interpolate_replace_nans(image, kernel)

您可能希望使用基于内核的插值的一些上下文包括:

  • 像素饱和的图像。一般来说,这些是成像区域中强度最高的区域,插值值不可靠,但这对于显示目的是有用的。

  • 带有标记像素的图像(例如,受宇宙射线或其他伪信号影响的一些小区域需要标记这些像素)。如果受影响的区域足够小,所得到的插值将对源统计数据产生较小的影响,并允许在结果数据上运行健壮的源查找算法。

  • 稀疏采样的图像,例如用单像素探测器构建的图像。这样的图像只有几个离散点在整个成像区域内采样,但是仍然可以构造出扩展天空发射的近似值。

注解

必须注意确保内核足够大,足以完全覆盖NaN值的潜在相邻区域。安 AstropyUserWarning 如果在后卷积检测到NaN值,则引发,在这种情况下,应增大内核大小。

例子

下面的脚本显示了填充标记的像素的内核插值示例:

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 Gaussian2DKernel, interpolate_replace_nans

# Load the data from data.astropy.org
filename = get_pkg_data_filename('galactic_center/gc_msx_e.fits')

hdu = fits.open(filename)[0]
img = hdu.data[50:90, 60:100] * 1e5

# This example is intended to demonstrate how astropy.convolve and
# scipy.convolve handle missing data, so we start by setting the brightest
# pixels to NaN to simulate a "saturated" data set
img[img > 2e1] = 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)

# create a "fixed" image with NaNs replaced by interpolated values
fixed_image = interpolate_replace_nans(img, kernel)

# Now we do a bunch of plots.  In the first two plots, the originally masked
# values are marked with red X's
plt.figure(1, figsize=(12, 6)).clf()
plt.close(2) # close the second plot from above

ax1 = plt.subplot(1, 2, 1)
im = ax1.imshow(img, vmin=-2., vmax=2.e1, origin='lower',
                interpolation='nearest', cmap='viridis')
y, x = np.where(np.isnan(img))
ax1.set_autoscale_on(False)
ax1.plot(x, y, 'rx', markersize=4)
ax1.set_title("Original")
ax1.set_xticklabels([])
ax1.set_yticklabels([])

ax2 = plt.subplot(1, 2, 2)
im = ax2.imshow(fixed_image, vmin=-2., vmax=2.e1, origin='lower',
                interpolation='nearest', cmap='viridis')
ax2.set_title("Fixed")
ax2.set_xticklabels([])
ax2.set_yticklabels([])

(png _, svgpdf

../_images/index-3.png

例子

这个脚本展示了这种技术从稀疏采样重建图像的能力。请注意,图像并不完美:点状源有时会丢失,但扩展的结构可以通过眼睛很好地恢复。

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 Gaussian2DKernel, interpolate_replace_nans

# Load the data from data.astropy.org
filename = get_pkg_data_filename('galactic_center/gc_msx_e.fits')

hdu = fits.open(filename)[0]
img = hdu.data[50:90, 60:100] * 1e5

indices = np.random.randint(low=0, high=img.size, size=300)

sampled_data = img.flat[indices]

# Build a new, sparsely sampled version of the original image
new_img = np.tile(np.nan, img.shape)
new_img.flat[indices] = sampled_data

# We smooth with a Gaussian kernel with x_stddev=1 (and y_stddev=1)
# It is a 9x9 array
kernel = Gaussian2DKernel(x_stddev=1)

# create a "reconstructed" image with NaNs replaced by interpolated values
reconstructed_image = interpolate_replace_nans(new_img, kernel)

# Now we do a bunch of plots.  In the first two plots, the originally masked
# values are marked with red X's
plt.figure(1, figsize=(12, 6)).clf()
ax1 = plt.subplot(1, 3, 1)
im = ax1.imshow(img, vmin=-2., vmax=2.e1, origin='lower',
                interpolation='nearest', cmap='viridis')
y, x = np.where(np.isnan(img))
ax1.set_autoscale_on(False)
ax1.set_title("Original")
ax1.set_xticklabels([])
ax1.set_yticklabels([])

ax2 = plt.subplot(1, 3, 2)
im = ax2.imshow(new_img, vmin=-2., vmax=2.e1, origin='lower',
                interpolation='nearest', cmap='viridis')
ax2.set_title("Sparsely Sampled")
ax2.set_xticklabels([])
ax2.set_yticklabels([])

ax2 = plt.subplot(1, 3, 3)
im = ax2.imshow(reconstructed_image, vmin=-2., vmax=2.e1, origin='lower',
                interpolation='nearest', cmap='viridis')
ax2.set_title("Reconstructed")
ax2.set_xticklabels([])
ax2.set_yticklabels([])

(png _, svgpdf

../_images/index-4.png

关于向后兼容性的说明(v2.0之前)

行为 astropy 的直接卷积 (convolve() )在版本2.0中更改。一般来说,旧版本是不可取的。但是,要恢复旧的行为 (astropy 版本<2.0)直接卷积函数,您可以插值然后卷积,例如:

from astropy.convolution import interpolate_replace_nans, convolve
interped_result = interpolate_replace_nans(image, kernel)
result = convolve(interped_image, kernel)

请注意,这两个 convolveconvolve_fft 是为了表演 归一化卷积 在这个过程中插入nan。本文中给出的示例,以及以前在旧版本中仅在直接卷积中完成的操作 astropy 现在执行两个步骤:首先,它用它们的插值值替换NaN,同时保持所有非NaN值不变,然后用指定的内核卷积得到的图像。

性能提示

这个 convolve() 函数最适合于小内核,而对于较大的内核则可能变得非常慢。在这种情况下,考虑使用 convolve_fft() (请注意,此函数使用更多内存)。

参考/API

星形卷积包裹

功能

convolve \(数组,内核)[, boundary, ...] )

用内核卷积数组。

convolve_fft \(数组,内核)[, boundary, ...] )

用nd核卷积ndarray。

convolve_models \(型号,内核)[, mode] )

使用卷积两个模型 convolve_fft .

discretize_model \(型号,xU系列[, y_range, ...] )

函数在网格上计算分析模型函数。

interpolate_replace_nans \(数组,内核)[, ...] )

给定一个包含nan的数据集,用给定的核从相邻的数据点插值来替换nan。

kernel_arithmetics \(内核,值,操作)

加、减或乘两个核。

Classes

AiryDisk2DKernel \(半径, *  * 克瓦格斯)

2D Airy磁盘内核。

Box1DKernel \(宽度, *  * 克瓦格斯)

1D盒过滤器内核。

Box2DKernel \(宽度, *  * 克瓦格斯)

二维盒过滤器内核。

CustomKernel \(数组)

从列表或数组创建筛选器内核。

Gaussian1DKernel \(标准偏差, *  * 克瓦格斯)

一维高斯滤波器核。

Gaussian2DKernel \(x u标准偏差[, y_stddev, theta] )

二维高斯滤波核。

Kernel \(数组)

卷积内核基类。

Kernel1D \ [model, x_size, array] )

1D筛选器内核的基类。

Kernel2D \ [model, x_size, y_size, array] )

二维筛选器内核的基类。

Model1DKernel \(型号, *  * 克瓦格斯)

从1D模型创建内核。

Model2DKernel \(型号, *  * 克瓦格斯)

从二维模型创建内核。

Moffat2DKernel \(伽马,阿尔法, *  * 克瓦格斯)

二维Moffat内核。

RickerWavelet1DKernel \(宽度, *  * 克瓦格斯)

一维Ricker小波滤波器核(有时称为“墨西哥帽”核)。

RickerWavelet2DKernel \(宽度, *  * 克瓦格斯)

二维Ricker小波滤波器核(有时称为“墨西哥帽”核)。

Ring2DKernel \(半径u in,宽度, *  * 克瓦格斯)

二维环形滤波器内核。

Tophat2DKernel \(半径, *  * 克瓦格斯)

2D Tophat过滤器内核。

Trapezoid1DKernel (宽度) [, slope] )

一维梯形核。

TrapezoidDisk2DKernel \(半径[, slope] )

二维梯形核。