复杂的二维背景#
天空背景估计#
使用案例： 在复杂场景中估计天空背景并评估天空估计的质量。
数据： 在笔记本中创建的具有病态背景模式的图像。
工具： photutils。
跨仪器： 所有仪器。
文档： 本笔记本是STScI更大后处理数据分析工具生态系统的一部分。
介绍#
天空估计是图像处理中最棘手的方面之一。这部分是因为“天空”是天文场景的一部分。某些被视为背景的内容可能位于前景（探测器中的热背景、散射光或黄道光）。有时感兴趣的物体会重叠（例如，前面有星系的其他星系）。本笔记本不解决重叠星系的“去混叠”问题，但概述了一些处理天空中大规模模式的技术。
Photutils手册对背景估计进行了广泛讨论，这也是本笔记本中许多内容的基础。
导入模块#
Numpy 用于一般数组计算
Scipy 用于一些统计操作和插值
Photutils 用于光度计算
Astropy 表格用于查看天空背景瑕疵的参数
Astropy 卷积用于平滑图像
Photutils 中的例程用于背景减除和掩膜处理
Astropy 表格用于操作注入到图像中的源（星系）列表
Astropy 卷积用于扩展天空图像和源掩膜
Astropy 建模用于拟合背景减除残差的直线
Astropy block_reduce 用于查看残差的均方根（RMS）作为尺度（分块因子）的函数
import numpy as np # 导入NumPy库，用于数值计算
from astropy.convolution import (Box2DKernel, Gaussian2DKernel, Ring2DKernel,
Tophat2DKernel, convolve) # 导入Astropy的卷积模块，用于图像卷积
from astropy.modeling import fitting, models # 导入Astropy的建模模块，用于模型拟合
from astropy.nddata.blocks import block_reduce # 导入块减少函数，用于降采样
from astropy.stats import SigmaClip # 导入SigmaClip，用于统计分析中的剪切
from astropy.table import Table # 导入表格模块，用于数据表的创建和操作
from IPython.display import Image # 导入Image模块，用于在IPython中显示图像
from jdaviz import Imviz # 导入Jdaviz的Imviz模块，用于图像可视化
from photutils.background import (Background2D,
BkgIDWInterpolator, BkgZoomInterpolator,
MedianBackground) # 导入Photutils的背景处理模块
from photutils.datasets import make_model_image # 导入数据集模块，用于生成模型图像
from photutils.segmentation import detect_sources, make_2dgaussian_kernel # 导入分割模块，用于源检测和生成高斯核
from scipy.ndimage import median_filter # 导入SciPy的中值滤波函数
Matplotlib 绘图设置#
有两个版本：
notebook – 提供交互式图形，但会使整个笔记本的滚动变得稍微困难
inline – 提供非交互式图形，以便更好地滚动整体内容
import matplotlib.pyplot as plt # 导入绘图库，用于绘制图形
import matplotlib as mpl # 导入matplotlib库的其他功能
# 使用此版本进行非交互式绘图（便于在笔记本中滚动）
%matplotlib inline # 设置为内联模式，图形将在笔记本中直接显示
# 如果你想要交互式图形，可以使用此版本
# %matplotlib notebook # 设置为交互模式，允许与图形进行交互
# 这些设置是为了确保在内联和笔记本版本中图形的大小相同
%config InlineBackend.print_figure_kwargs = {'bbox_inches': None} # 配置图形输出时的边界框设置
mpl.rcParams['savefig.dpi'] = 80 # 设置保存图形时的分辨率为80 dpi
mpl.rcParams['figure.dpi'] = 80 # 设置图形显示时的分辨率为80 dpi
np.random.seed(seed=4) # 设置随机种子为4，以确保结果的可重复性
创建具有某种病态背景模式的图像#
该图像的像素间均方根（RMS）为0.1，我们添加的背景非均匀性与像素间噪声的水平相当。因此，无法在不进行相当程度的平滑处理的情况下去除这些非均匀性。我们将添加：
一个背景渐变
几个散布在图像中的sinc函数瑕疵。使sinc函数的周期足够大，以至于这些瑕疵看起来不太像源（假设它们的大小仅为几个像素到几十个像素）。
这个特定的测试案例被巧妙地选择为一种情况，仅仅通过添加更高阶的（例如）切比雪夫多项式表面拟合可能效果不佳。同样，增加双三次样条的阶数也不太可能令人满意。
设置网格和随机数种子以进行模拟#
该测试最好通过一个2000x2000的图像来说明，但可以使用shrink_factor来加快笔记本的运行速度以便于测试。随机数种子被设置为允许结果的可重复性。
default_size = 2000 # 默认图像大小设为2000
shrink_factor = 2 # 缩小因子，用于加快执行速度
nrow = ncol = default_size // shrink_factor # 计算缩小后的图像行列数
row, col = np.mgrid[0:nrow, 0:ncol] # 创建一个网格，用于生成行和列的坐标
创建恶劣的天空背景#
瑕疵作为带有随机中心、宽度和幅度的 sinc 函数插入。
nblemish = 16 // shrink_factor # 计算要添加的sinc函数瑕疵数量
row_centers = ncol * np.random.random_sample(size=nblemish) # 随机生成瑕疵的行中心位置
col_centers = nrow * np.random.random_sample(size=nblemish) # 随机生成瑕疵的列中心位置
# 使sinc波纹的宽度相对较大
width = 30 + 100. * np.random.random_sample(size=nblemish) # 随机生成瑕疵的宽度
amplitude = 1. * np.random.random_sample(size=nblemish) # 随机生成瑕疵的幅度
noiseless_sky = np.zeros((nrow, ncol), dtype=np.float32) # 创建一个全零的图像数组，用于存储无噪声的天空图像
for rr, cc, w, a in zip(row_centers, col_centers, width, amplitude): # 遍历每个瑕疵的参数
r = np.sqrt((row - rr) ** 2 + (col - cc) ** 2) # 计算每个像素到瑕疵中心的距离
noiseless_sky = noiseless_sky + a * np.sinc(r / (w * np.pi)) # 将sinc函数的值添加到无噪声天空图像中
blemishes = Table([col_centers, row_centers, width, amplitude], # 创建一个表格，存储瑕疵的参数
names=['x', 'y', 'width', 'amplitude']) # 指定表格的列名
blemishes # 输出瑕疵参数表
在瑕疵上添加渐变
# 计算梯度，row和col分别表示行和列的索引
gradient = (1.*row)/nrow + (0.25*col)/ncol
# 将梯度添加到无噪声的天空背景中
noiseless_sky += gradient
# 生成与无噪声天空背景相同形状的随机噪声，标准差为0.1
noise = np.random.normal(scale=0.1, size=np.shape(noiseless_sky))
# 将噪声添加到无噪声天空背景中，得到有噪声的天空背景
sky_bkgd = noiseless_sky + noise
# 计算有噪声天空背景的平均值
mean_bkgd = sky_bkgd.mean()
import matplotlib.pyplot as plt # 导入matplotlib库中的pyplot模块
plt.figure(figsize=(6, 6)) # 创建一个6x6英寸的图形
plt.imshow(noiseless_sky, origin='lower') # 显示无噪声的天空图像，原点位于图像的下方
plt.title('The noiseless sky background') # 设置图形的标题为“无噪声的天空背景”
或者，尝试在 Imviz 中打开图像，这是一个新的交互式图像查看工具。有关 Imviz 的更多信息，请访问此链接：Imviz 文档
imviz = Imviz() # 创建Imviz对象，用于数据可视化
viewer = imviz.default_viewer # 获取Imviz的默认查看器
imviz.load_data(noiseless_sky, data_label='noiseless_sky_background') # 加载无噪声的天空数据，并为其指定标签
# 显示JWST望远镜数据的可视化界面
imviz.show()
考虑一些可以从笔记本运行的示例 API 命令，尽管许多操作可以在 Imviz 中完成。
# 获取当前查看器的拉伸选项
viewer.stretch_options
# 设置查看器的拉伸方式为平方根拉伸
viewer.stretch = 'sqrt'
# 获取当前查看器的颜色映射选项
viewer.colormap_options
# 设置查看器的颜色映射为'Viridis'
viewer.set_colormap('Viridis')
加入噪声后观察，并稍作平滑#
使用10x10的箱形核进行平滑
# 对背景天空数据进行卷积处理，使用Box2DKernel进行平滑
sky_smoothed = convolve(sky_bkgd-mean_bkgd, Box2DKernel(10))
# 计算平滑后的天空数据的最小和最大值，用于后续显示的色阶
zmin, zmax = np.percentile(sky_smoothed, (0.1, 99.9))
# 创建一个图形，设置图形大小
plt.figure(figsize=(10, 5))
# 创建第一个子图，显示原始天空数据
ax1 = plt.subplot(121)
# 显示原始天空数据，设置颜色范围
ax1.imshow(sky_bkgd-mean_bkgd, vmin=zmin, vmax=zmax, origin='lower')
# 设置第一个子图的标题
ax1.set_title('The sky with noise')
# 创建第二个子图，分享x和y轴与第一个子图
ax2 = plt.subplot(122, sharex=ax1, sharey=ax1)
# 设置第二个子图的标题
ax2.set_title('The sky smoothed')
# 显示平滑后的天空数据，设置颜色范围
ax2.imshow(sky_smoothed, vmin=zmin, vmax=zmax, origin='lower')
在Imviz中执行相同操作，并并排比较图像，按像素锁定。
imviz2 = Imviz() # 创建Imviz对象，用于处理JWST数据
viewer = imviz2.default_viewer # 获取Imviz的默认查看器
# 调用imviz2对象的show方法以显示图像
imviz2.show()
# 加载数据，计算天空背景减去平均背景，并为数据设置标签
imviz2.load_data(sky_bkgd - mean_bkgd, data_label='sky_with_noise')
# 获取当前的拉伸选项
viewer.stretch_options
# 设置视图的拉伸方式为平方根拉伸
viewer.stretch = 'sqrt'
# 获取当前的颜色映射选项
viewer.colormap_options
# 设置视图的颜色映射为'Viridis'
viewer.set_colormap('Viridis')
# 加载平滑后的天空数据，并指定数据标签
imviz2.load_data(sky_smoothed, data_label='sky_smoothed')
# 定义第二个查看器的名称
viewer2_name = 'viewer2'
# 创建一个新的图像查看器，并将其命名为 viewer2
viewer2 = imviz2.create_image_viewer(viewer_name=viewer2_name)
# 将数据添加到指定的查看器中
imviz2.app.add_data_to_viewer(viewer2_name, 'sky_smoothed')
# 获取查看器的拉伸选项
viewer2.stretch_options
# 设置查看器的拉伸方式为平方根
viewer2.stretch = 'sqrt'
# 获取查看器的颜色映射选项
viewer2.colormap_options
# 设置查看器的颜色映射为'Viridis'
viewer2.set_colormap('Viridis')
这两幅图像默认是通过像素链接的。点击缩放/平移工具可以更仔细地查看图像。这两幅图像将同步缩放和平移。
# 导入Image类，用于显示图像
Image('./imviz_2images.png', alt='Imviz with two viewers. Highlighted tool for linked zoom and pan.')
添加一些随机源#
我们将创建一个源目录，这些源具有椭圆高斯（elliptical Gaussian）轮廓，具有多种通量（flux）、大小（size）和位置角（position angle）。这些源的大小相对于背景结构保持相对较小。为了方便起见，将这些源存储在Astropy表中。
sources = Table() # 创建一个空的表格用于存储源数据
nsources = 5000//(shrink_factor**2) # 计算源的数量，考虑缩放因子
rand_sample = np.random.random_sample # 保存随机数生成函数以减少输入
random_sample = rand_sample(nsources) # 生成 nsources 个随机数
sources['x_mean'] = 2.+(ncol-4.)*rand_sample(nsources) # 计算 x 坐标的均值，避免边缘
sources['y_mean'] = 2.+(nrow-4.)*rand_sample(nsources) # 计算 y 坐标的均值，避免边缘
sources['flux'] = 50*rand_sample(nsources) # 计算源的流量，范围在 0 到 50 之间
sources['x_stddev'] = 5.*rand_sample(nsources) # 计算 x 坐标的标准差
sources['y_stddev'] = 5.*rand_sample(nsources) # 计算 y 坐标的标准差
sources['re'] = np.sqrt(sources['x_stddev']*sources['y_stddev']) # 计算有效半径
sources['theta'] = 180.*rand_sample(nsources)*np.pi / 180. # 计算角度，转换为弧度
将这些来源添加到图像中以创建一个场景。
model = models.Gaussian2D() # 创建一个二维高斯模型
source_img = make_model_image(sky_bkgd.shape, model, sources, # 生成源图像，使用背景图像的形状和源数据
x_name='x_mean', y_name='y_mean') # 指定高斯模型的x和y均值参数名称
scene = source_img + sky_bkgd # 将源图像与背景图像相加，形成最终场景
plt.figure(figsize=(6, 6)) # 创建一个6x6英寸的图形
plt.imshow(scene-mean_bkgd, vmin=zmin, vmax=zmax, origin='lower') # 显示场景图像，减去背景，设置颜色范围和原点位置
plt.title('The scene, with the sources and nasty sky background') # 设置图形标题
并排比较两幅图像的绘图例程#
def plot_two(img1, img2, vmin, vmax, figsize=(10, 6), titles=['', ''],
cmap=plt.cm.viridis):
# 创建一个新的图形，设置图形大小
ax1 = plt.figure(figsize=figsize)
# 在图形中添加第一个子图，位置为121
ax1 = plt.subplot(121)
# 显示第一个图像，设置颜色范围和颜色映射
ax1.imshow(img1, vmin=vmin, vmax=vmax, cmap=cmap, origin='lower')
# 设置第一个子图的标题
ax1.set_title(titles[0])
# 在图形中添加第二个子图，位置为122，并共享x轴和y轴
ax2 = plt.subplot(122, sharex=ax1, sharey=ax1)
# 显示第二个图像，设置颜色范围和颜色映射
ax2.imshow(img2, vmin=vmin, vmax=vmax, cmap=cmap, origin='lower')
# 设置第二个子图的标题
ax2.set_title(titles[1])
尝试使用 Photutils 的 Background2D 网格作为第一次尝试#
现在我们需要对背景进行初步估计。Photutils 提供了一个灵活的接口，用于在一个网格单元中估计背景。我们在这里使用 Background2D 例程，并设置了一组相当典型的参数。有关更多信息，请参见 文档。
此任务创建一个覆盖图像的矩形单元网格。每个单元内未被遮罩的像素使用可配置的稳健统计量进行“平均”。然后，这些平均值可以在更大范围内使用可配置的稳健统计量进行“平均”。这个经过滤波的平均值集随后用于输入插值例程，以生成平滑的背景。默认的插值方法是双三次样条插值，但我们将在后面的笔记本中展示反距离加权插值。
尝试调整 Background2D 的参数以查看效果。
第二个参数（下面的 50）是网格大小。如果需要，这也可以是一个矩形，例如 (50,40)。
sigma_clip是在网格单元内进行稳健平均时使用的方法。filter_size是在进行插值之前要“平均”的单元数量。如果需要，这也可以是一个矩形，例如 (3,4)。bkg_estimator是用于平均单元中值的方法。exclude_percentile如果一个网格的像素遮罩超过此百分比，则它将被排除在低分辨率图中。interpolator是用于插值背景的方法（在这种情况下为双三次样条插值）。
# 创建一个背景估计对象，使用2D背景估计
bkg = Background2D(scene, box_size=50, # 设置背景估计的场景和盒子大小为50
sigma_clip=SigmaClip(sigma=3.), # 使用SigmaClip进行噪声剪切，设置sigma为3
filter_size=3, # 设置滤波器大小为3
exclude_percentile=50, # 排除前50%的像素值
bkg_estimator=MedianBackground(), # 使用中位数背景估计
interpolator=BkgZoomInterpolator(order=3)) # 使用三阶插值器进行背景插值
plt.figure(figsize=(6, 6)) # 创建一个6x6英寸的图形
bkgimg = bkg.background # 获取估计的背景图像
# 显示背景图像与无噪声天空图像的差异，设置颜色范围
plt.imshow(bkgimg - noiseless_sky, vmin=0.1 * zmin, vmax=0.1 * zmax, origin='lower') # 显示图像，设置最小值和最大值
制作几个掩膜。make_source_mask 例程有几个选项。尝试更改它们以查看效果。这里的目标是找到并掩膜源，而不在背景较高的区域看到明显的源增强。可以调整几个参数：
nsigma– 我们将尝试使用 2 和 3 sigma 截断npixels– 在掩膜之前，要求超过阈值的连接像素数量dilate_size– 用于扩展掩膜的方形箱形滤波器的大小filter_fwhm– 在阈值处理之前，图像与此 FWHM 的高斯函数进行卷积。
def make_source_mask(data, nsigma, npixels, filter_fwhm=3.0, dilate_size=11):
# 创建一个二维高斯核，使用给定的滤波器全宽半最大值（FWHM）和大小
kernel = make_2dgaussian_kernel(filter_fwhm, size=3)
# 使用高斯核对数据进行卷积
convolved_data = convolve(data, kernel)
# 创建一个SigmaClip对象，用于进行sigma剪切，设置sigma和最大迭代次数
sigma_clip = SigmaClip(sigma=3.0, maxiters=5)
# 对数据进行sigma剪切，返回剪切后的数据
clipped_data = sigma_clip(data, masked=False)
# 计算剪切后数据的均值
mean_ = np.nanmean(clipped_data)
# 计算剪切后数据的标准差
std_ = np.nanstd(clipped_data)
# 计算阈值，使用均值和指定的nsigma倍数的标准差
threshold_2sigma = mean_ + nsigma * std_
# 检测源，返回分割图像
segm = detect_sources(convolved_data, threshold_2sigma, npixels)
# 返回源掩膜，使用指定的膨胀大小
return segm.make_source_mask(size=dilate_size)
# 创建一个2σ的源掩膜
mask_2sigma = make_source_mask(scene-bkgimg, nsigma=2, npixels=3,
dilate_size=5, filter_fwhm=3)
# 创建一个3σ的源掩膜
mask_3sigma = make_source_mask(scene-bkgimg, nsigma=3, npixels=3,
dilate_size=5, filter_fwhm=3)
# 绘制2σ和3σ掩膜的比较图
plot_two(mask_2sigma, mask_3sigma, 0, 1)
尝试在 Background2D 任务中减小背景估计的网格大小。较小的网格大小确实有帮助，但在右侧的残差图像中有一个迹象表明，源对象正在影响背景估计。
sigma_clip = SigmaClip(sigma=3.) # 创建一个SigmaClip对象，设置sigma为3，用于背景噪声的剪切
bkg_estimator = MedianBackground() # 创建一个MedianBackground对象，用于估计背景
# 创建第一个背景估计对象bkg1，使用30x30的区域和3x3的滤波器
bkg1 = Background2D(scene, (30, 30), filter_size=(3, 3), mask=mask_3sigma,
sigma_clip=sigma_clip, bkg_estimator=bkg_estimator)
# 创建第二个背景估计对象bkg2，使用15x15的区域和3x3的滤波器
bkg2 = Background2D(scene, (15, 15), filter_size=(3, 3), mask=mask_3sigma,
sigma_clip=sigma_clip, bkg_estimator=bkg_estimator)
bkg1img = bkg1.background # 获取第一个背景估计的图像
bkg2img = bkg2.background # 获取第二个背景估计的图像
# 绘制两个背景图像与无噪声天空图像的差异，设置z轴范围和标题
plot_two(bkg1img-noiseless_sky, bkg2img-noiseless_sky,
0.1*zmin, 0.1*zmax, titles=['bkg1', 'bkg2'])
环形中值滤波图像#
显然，我们需要进行更多的迭代来去除源。在此之前，让我们尝试另一种去除源的方法：环形中值滤波。为此，创建一个比感兴趣的源更大的滤波核。使用 scipy 的中值滤波例程，将其滑动到图像上，计算环内所有像素的中值。这基本上为每个像素提供了局部背景估计。
（虽然我们在这里没有说明，但可以使用 scipy.ndimage.generic_filter 并传递一个函数来计算其他内容，而不仅仅是中值。）
下面的单元格设置了滤波核。
ring = Ring2DKernel(15, 5) # 创建一个半径为15，宽度为5的环形2D核
plt.imshow(ring.array) # 显示环形核的数组图像
比较场景（去除均值背景）与过滤后的场景。
# 对场景数据应用中值滤波，使用环形的足迹
filtered = median_filter(scene, footprint=ring.array)
# 绘制两个图像：一个是去除背景后的源图像，另一个是去除背景后的中值滤波图像
plot_two(scene-mean_bkgd, filtered-mean_bkgd, zmin, zmax,
titles=['Sources', 'ring-median filtered']) # 设置图像标题
从场景中减去环形中值滤波图像，作为背景减除的第一步。用一个卷积核对其进行卷积，以平滑处理以便进行源检测，然后对其进行阈值处理，以识别出属于源的像素。
# 将场景过滤后的数据赋值给变量difference
difference = scene-filtered
# 使用高斯2D卷积核对difference进行平滑处理，卷积核大小为3
smoothed = convolve(difference, Gaussian2DKernel(3))
# 创建一个掩膜，掩膜的条件是平滑后的数据大于其标准差的1倍
mask = smoothed > 1.*smoothed.std()
# 绘制两个图像：原始差异图和掩膜图
plot_two(difference, mask, zmin, zmax,
titles=['scene minus ring-median filtered scene', 'mask'])
通常，扩展掩模是很有用的。这可以通过与滤波器进行卷积来实现。在这里，我们采用圆形顶帽滤波器。
def dilate_mask(mask, tophat_size):
''' 取一个掩膜并使被掩膜区域变大。'''
area = np.pi * tophat_size**2. # 计算圆形区域的面积
kernel = Tophat2DKernel(tophat_size) # 创建一个二维的顶帽核
dilated_mask = convolve(mask, kernel) >= 1. / area # 使用卷积扩展掩膜区域
return dilated_mask # 返回扩展后的掩膜
# 使用dilate_mask函数对mask进行膨胀处理，膨胀程度为2
dilated_mask = dilate_mask(mask, 2)
# 调用plot_two函数绘制原始mask和膨胀后的mask，设置子图的索引和标题
plot_two(mask, dilated_mask, 0, 1, titles=['mask', 'dilated mask'])
对图像进行遮罩处理，查看仍然剩余多少源。
# 绘制两个场景图像
plot_two(scene-noiseless_sky, # 第一个图像：无噪声的天空场景
(scene-noiseless_sky)*(1.-mask), # 第二个图像：应用掩膜后的无噪声天空场景
zmin, # 图像的最小值
zmax, # 图像的最大值
titles=['scene without background', # 图像标题：无背景场景
'masked scene without background']) # 图像标题：无背景的掩膜场景
让我们制作一个更精美的图表，展示背景残差、包含源的残差以及叠加源掩模的图像，以及背景减去后的场景。
def plot_bkgd(scene, mask, bkgd, noiseless_bkgd, sources,
zmin, zmax, factor=0.1, # 控制颜色映射的拉伸
figsize=(20, 10)):
'''制作一个三面板的图像：
* 残余天空背景
* 残余天空背景乘以掩膜，掩膜源以'+'符号叠加，
未掩膜源以圆圈叠加，以及
* 背景减去的场景。
'''
residual = bkgd - noiseless_bkgd # 计算残余背景
plt.figure(figsize=figsize) # 创建一个指定大小的图形
# 绘制背景残余
ax1 = plt.subplot(131) # 创建第一个子图
ax1.imshow(residual, vmin=factor*zmin, vmax=factor*zmax, origin='lower') # 显示残余图像
ax1.set_title('residual sky background') # 设置标题
ax2 = plt.subplot(132, sharex=ax1, sharey=ax1) # 创建第二个子图，分享x和y轴
ax2.imshow(residual * (1 - mask), vmin=factor*zmin, vmax=factor*zmax, origin='lower') # 显示掩膜后的残余图像
xs = np.rint(sources['x_mean']).astype(np.int32) # 获取源的x坐标并四舍五入为整数
ys = np.rint(sources['y_mean']).astype(np.int32) # 获取源的y坐标并四舍五入为整数
masked = mask[ys, xs] # 获取被掩膜的源
unmasked = np.invert(masked) # 获取未被掩膜的源
ax2.scatter(xs[masked], ys[masked], s=sources['flux'][masked], marker='+', # 绘制被掩膜源
c='red', alpha=0.2) # 设置颜色和透明度
ax2.scatter(xs[unmasked], ys[unmasked], s=sources['flux'][unmasked], # 绘制未被掩膜源
c='black', alpha=0.5) # 设置颜色和透明度
ax2.set_xlim(0, scene.shape[1]) # 设置x轴范围
ax2.set_ylim(0, scene.shape[0]) # 设置y轴范围
ax2.set_title('with masked + unmasked O sources') # 设置标题
ax3 = plt.subplot(133, sharex=ax1, sharey=ax1) # 创建第三个子图，分享x和y轴
ax3.imshow(scene - bkgd, vmin=zmin, vmax=zmax, origin='lower') # 显示背景减去的场景图像
ax3.set_title('Scene minus estimated background') # 设置标题
print(f"Residual RMS, peak = {residual.std():.4f}, {residual.max():.4f}") # 打印残余的标准差和最大值
试试这个更复杂的图。前两个图的默认缩放在颜色映射上有更大的拉伸，因此我们现在可以看到所有明亮源周围都有环状残差。这在第三个图中并不明显，因为它使用的是我们之前图形的缩放方式。
def plot_bkgd(scene, mask, filtered, noiseless_sky, sources, zmin, zmax):
"""
绘制背景图像和源图像的函数
:param scene: 场景图像数据
:param mask: 掩膜数据
:param filtered: 过滤后的图像数据
:param noiseless_sky: 无噪声天空图像数据
:param sources: 源图像数据
:param zmin: z轴最小值
:param zmax: z轴最大值
"""
import matplotlib.pyplot as plt # 导入绘图库
import numpy as np # 导入数值计算库
# 创建一个新的图形
plt.figure(figsize=(15, 10)) # 设置图形大小
# 绘制无噪声天空图像
plt.subplot(2, 2, 1) # 创建2x2的子图，选择第1个子图
plt.imshow(noiseless_sky, vmin=zmin, vmax=zmax, cmap='gray') # 显示无噪声天空图像
plt.title('Noiseless Sky') # 设置标题
plt.colorbar() # 添加颜色条
# 绘制场景图像
plt.subplot(2, 2, 2) # 选择第2个子图
plt.imshow(scene, vmin=zmin, vmax=zmax, cmap='gray') # 显示场景图像
plt.title('Scene') # 设置标题
plt.colorbar() # 添加颜色条
# 绘制过滤后的图像
plt.subplot(2, 2, 3) # 选择第3个子图
plt.imshow(filtered, vmin=zmin, vmax=zmax, cmap='gray') # 显示过滤后的图像
plt.title('Filtered') # 设置标题
plt.colorbar() # 添加颜色条
# 绘制源图像
plt.subplot(2, 2, 4) # 选择第4个子图
plt.imshow(sources, vmin=zmin, vmax=zmax, cmap='gray') # 显示源图像
plt.title('Sources') # 设置标题
plt.colorbar() # 添加颜色条
plt.tight_layout() # 调整子图间距
plt.show() # 显示图形
设置一些便于迭代掩膜的例程#
接下来，我们将尝试在拟合背景时进行迭代，并逐步去除信噪比（S/N）较低的源。我们希望在每一步都进行检查。以下是一些减少输入的函数：
SourceMask– 这是一个类，用于设置掩膜的一些参数，并提供几个方法：single– 应用现有掩膜，然后使用 photutils 的 make_source_mask 创建一个新掩膜；用圆形顶帽（circular tophat）卷积掩膜，并进行阈值处理以扩展它。multiple– 重复应用single方法以生成新掩膜。或者将其与之前的掩膜进行逻辑或（OR）操作。
find_worst_residual_near_center– 在绘制放大图像以进行检查时，我们希望选择远离边缘的区域，该区域具有最差的残差。plot_mask– 绘制整个图像的掩膜，并将其作为轮廓叠加在小子区域上进行绘制。
def my_background(img, box_size, mask, interp=None, filter_size=1,
exclude_percentile=90):
''' 运行 photutils 背景处理，使用 SigmaClip 和 MedianBackground '''
if interp is None:
interp = BkgZoomInterpolator() # 如果没有提供插值器，则使用默认的 BkgZoomInterpolator
return Background2D(img, box_size, # 创建一个 2D 背景对象
sigma_clip=SigmaClip(sigma=3.), # 使用 SigmaClip 进行背景估计，设置 sigma 为 3
filter_size=filter_size, # 设置滤波器的大小
bkg_estimator=MedianBackground(), # 使用中值背景估计
exclude_percentile=exclude_percentile, # 设置排除的百分位数
mask=mask, # 应用掩膜
interpolator=interp) # 使用指定的插值器
class SourceMask:
def __init__(self, img, nsigma=3., npixels=3):
'''
初始化源掩膜的帮助类。
参数：
img: 输入图像
nsigma: 用于掩膜的标准差阈值
npixels: 用于掩膜的像素数量
'''
self.img = img # 保存输入图像
self.nsigma = nsigma # 保存标准差阈值
self.npixels = npixels # 保存像素数量
def single(self, filter_fwhm=3., tophat_size=5., mask=None):
'''
在单一尺度上创建掩膜
参数：
filter_fwhm: 高斯滤波器的全宽半最大值
tophat_size: 顶帽滤波器的大小
mask: 可选的掩膜
'''
if mask is None:
image = self.img # 如果没有提供掩膜，则使用原始图像
else:
image = self.img * (1 - mask) # 使用掩膜来调整图像
# 创建源掩膜
mask = make_source_mask(image, nsigma=self.nsigma,
npixels=self.npixels,
dilate_size=1, filter_fwhm=filter_fwhm)
return dilate_mask(mask, tophat_size) # 返回扩展后的掩膜
def multiple(self, filter_fwhm=[3.], tophat_size=[3.], mask=None):
'''
在不同尺度上重复创建掩膜
参数：
filter_fwhm: 高斯滤波器的全宽半最大值列表
tophat_size: 顶帽滤波器的大小列表
mask: 可选的掩膜
'''
if mask is None:
self.mask = np.zeros(self.img.shape, dtype=bool) # 初始化掩膜为全False
for fwhm, tophat in zip(filter_fwhm, tophat_size):
smask = self.single(filter_fwhm=fwhm, tophat_size=tophat) # 创建单尺度掩膜
self.mask = self.mask | smask # 将当前掩膜与新掩膜进行逻辑或操作
return self.mask # 返回最终的综合掩膜
def find_worst_residual_near_center(resid):
'''查找最差残差的像素位置，避免边缘区域'''
yc, xc = resid.shape[0]/2., resid.shape[1]/2. # 计算图像中心的y和x坐标
radius = resid.shape[0]/3. # 定义一个半径，取图像高度的三分之一
y, x = np.mgrid[0:resid.shape[0], 0:resid.shape[1]] # 创建一个网格，包含所有像素的y和x坐标
mask = np.sqrt((y-yc)**2+(x-xc)**2) < radius # 创建一个掩膜，选择距离中心在半径内的像素
rmasked = resid*mask # 将残差图像与掩膜相乘，保留中心区域的残差值，其余区域为0
return np.unravel_index(np.argmax(rmasked), resid.shape) # 返回最差残差的像素位置
def plot_mask(scene, bkgd, mask, zmin, zmax, worst=None, smooth=0):
'''制作一个三面板图：
* 整个图像的掩膜,
* 场景乘以掩膜,
* 放大区域，掩膜以轮廓显示
'''
if worst is None: # 如果没有提供最差的残差位置
y, x = find_worst_residual_near_center(bkgd-noiseless_sky) # 查找中心附近的最差残差位置
else: # 如果提供了最差的残差位置
x, y = worst # 使用提供的位置
plt.figure(figsize=(20, 10)) # 创建一个20x10英寸的图形
plt.subplot(131) # 创建第一个子图
plt.imshow(mask, vmin=0, vmax=1, cmap=plt.cm.gray, origin='lower') # 显示掩膜图像
plt.subplot(132) # 创建第二个子图
if smooth == 0: # 如果不需要平滑
plt.imshow((scene-bkgd)*(1-mask), vmin=zmin, vmax=zmax, origin='lower') # 显示场景与背景的差异
else: # 如果需要平滑
smoothed = convolve((scene-bkgd)*(1-mask), Gaussian2DKernel(smooth)) # 对差异进行高斯平滑
plt.imshow(smoothed*(1-mask), vmin=zmin/smooth, vmax=zmax/smooth, origin='lower') # 显示平滑后的图像
plt.subplot(133) # 创建第三个子图
plt.imshow(scene-bkgd, vmin=zmin, vmax=zmax) # 显示场景与背景的差异
plt.contour(mask, colors='red', alpha=0.2) # 在图像上绘制掩膜的轮廓
plt.ylim(y-100, y+100) # 设置y轴的范围
plt.xlim(x-100, x+100) # 设置x轴的范围
return x, y # 返回最差残差的位置
在掩蔽源后更精细地估计背景#
# 计算背景，使用指定的场景、盒子大小、滤波器大小和掩膜
bkg3 = my_background(scene, box_size=10, filter_size=5, mask=mask)
# 获取计算出的背景图像
bkg3img = bkg3.background
# 绘制背景图像，传入场景、掩膜、背景图像、无噪声天空、源、最小值和最大值
plot_bkgd(scene, mask, bkg3img, noiseless_sky, sources, zmin, zmax)
尝试不同的插值算法。#
Shephard反距离加权法（Shephard inverse-distance weighting）搜索离感兴趣像素最近的N个网格点，并将它们的权重设置为 \(1/(D^p + \lambda)\)，其中 \(D\) 是到邻居的距离，\(p\) 是一个幂次，\(\lambda\) 是一个参数，用于使邻居的权重在靠近感兴趣像素时更加平滑。
# 创建一个反距离加权插值器，设置邻居数量为20，权重为1，正则化参数为30
interpolator = BkgIDWInterpolator(n_neighbors=20, power=1, reg=30)
# 计算背景，使用指定的场景、盒子大小、滤波器大小和掩膜
bkg4 = my_background(scene, box_size=10, filter_size=5, mask=mask,
interp=interpolator, exclude_percentile=90)
# 提取计算得到的背景图像
bkg4img = bkg4.background
# 绘制背景图像，包含场景、掩膜、背景图像、无噪声天空、源、最小和最大z值
plot_bkgd(scene, mask, bkg4img, noiseless_sky, sources, zmin, zmax)
制作更深的源掩膜#
此过程运行三次迭代，使用不同的核宽度和不同的顶帽尺寸来扩展掩膜。尝试调整 sigma、filter_fwhm 和 tophat_size，观察它们如何影响源的检测。目标是掩盖那些明显可见的源，同时保留足够的像素以便追踪背景。
# 创建源掩模，使用场景图像和背景图像，设置标准差为1.5
sm = SourceMask(scene-bkg4img, nsigma=1.5)
# 生成多个掩模，使用不同的滤波器全宽半高（FWHM）和方形滤波器大小
mask = sm.multiple(filter_fwhm=[1, 3, 5],
tophat_size=[4, 2, 1])
# 绘制掩模，显示场景图像、背景图像和生成的掩模，设置最小值和最大值
plot_mask(scene, bkg4img, mask, zmin, zmax)
# 计算掩模中所有元素的总和
mask.sum()
使用这个新掩模重新进行背景估计。#
调整 n_neighbors、box_size、filter_size 和 exclude_percentile，观察它们如何影响残差。
# 创建一个反距离加权插值器，设置邻居数量为20，权重为0，正则化参数为30
interpolator = BkgIDWInterpolator(n_neighbors=20, power=0, reg=30)
# 计算背景，使用指定的场景、盒子大小、滤波器大小和掩膜
bkg5 = my_background(scene, box_size=6, filter_size=3, mask=mask,
interp=interpolator, exclude_percentile=90)
# 获取计算得到的背景图像
bkg5img = bkg5.background
# 绘制背景图像，包含场景、掩膜、背景图像、无噪声天空、源、最小值和最大值
plot_bkgd(scene, mask, bkg5img, noiseless_sky, sources, zmin, zmax)
检查由于星系引起的偏差#
由于我们处理的是一个已知真实情况的模拟，我们可以直接评估偏差。（在真实图像的情况下，我们无法这样做；我们将在后面建议另一种测试。）
将每个星系下方中央3x3像素的残余背景与星系的通量进行制表，分别针对被遮罩和未被遮罩的情况。首先，我们需要获取这些值。对每个源中心的3x3像素取平均值。记录在拟合背景时被遮罩的像素和未被遮罩的像素，以便我们可以分别检查任何偏差。
bkgd = bkg5img # 将背景图像赋值给变量bkgd
residual_image = bkgd - noiseless_sky # 计算残差图像，等于背景减去无噪声天空
# 创建表格中的列，用于存储背景估计、残差和是否被遮罩的标志
sources['bkgd'] = 0. * sources['flux'] # 初始化背景列为0
sources['resid'] = 0. * sources['flux'] # 初始化残差列为0
sources['masked'] = np.zeros(len(sources), dtype=bool) # 初始化遮罩列为False
# 计算每个源中心的3x3像素值
for i in range(len(sources)): # 遍历每个源
s = sources[i] # 获取当前源
x = np.rint(s['x_mean']).astype(np.int32) # 获取源的x坐标并四舍五入为整数
y = np.rint(s['y_mean']).astype(np.int32) # 获取源的y坐标并四舍五入为整数
xmin, xmax = max(0, x - 1), min(ncol, x + 1) # 计算x坐标的最小和最大值
ymin, ymax = max(0, y - 1), min(nrow, y + 1) # 计算y坐标的最小和最大值
sources['resid'][i] = residual_image[ymin:ymax, xmin:xmax].mean() # 计算当前源的残差平均值
sources['bkgd'][i] = bkgd[ymin:ymax, xmin:xmax].mean() # 计算当前源的背景平均值
sources['masked'][i] = mask[y, x] # 如果中心像素被遮罩，则标记为True
制作一份随机位置的列表，并进行相同的测量，作为对照。
# 设置随机位置的数组
N_random = 5 * len(sources) # 随机位置的数量为源数量的5倍
resid_under_random = np.zeros(N_random, dtype=np.float64) # 初始化残差数组
bkgd_under_random = np.zeros(N_random, dtype=np.float64) # 初始化背景数组
masked_random = np.zeros(N_random, dtype=bool) # 初始化掩膜数组
# 从目录中随机获取一个通量和半径，以便我们可以绘图
# 随机源的结果与真实源在同一坐标轴上
random_flux = sources['flux'][np.random.randint(0, len(sources), N_random)] # 随机通量
random_re = sources['re'][np.random.randint(0, len(sources), N_random)] # 随机半径
rnd_x = np.random.randint(2, ncol - 2, size=N_random) # 随机生成x坐标
rnd_y = np.random.randint(2, nrow - 2, size=N_random) # 随机生成y坐标
# 计算以每个源为中心的3x3像素的值
for i, x, y in zip(range(len(rnd_x)), rnd_x, rnd_y): # 遍历随机坐标
resid_under_random[i] = residual_image[y - 1:y + 1, x - 1:x + 1].mean() # 计算残差的平均值
bkgd_under_random[i] = bkgd[y - 1:y + 1, x - 1:x + 1].mean() # 计算背景的平均值
masked_random[i] = mask[y - 1:y + 1, x - 1:x + 1].min().astype(bool) # 检查掩膜的最小值
# 仅保留未被掩膜的位置的结果
resid_under_random = resid_under_random[~masked_random] # 过滤未掩膜的残差
random_flux = random_flux[~masked_random] # 过滤未掩膜的随机通量
random_re = random_re[~masked_random] # 过滤未掩膜的随机半径
print(f"Nsources, Nrandom: {len(sources)} {len(resid_under_random)}") # 打印源数量和随机数量
列出被遮蔽和未被遮蔽的源#
被遮蔽的源 (Masked Sources)#
源1
源2
源3
未被遮蔽的源 (Unmasked Sources)#
源A
源B
源C
flux = sources['flux'] # 获取源的光通量数据
resid = sources['resid'] # 获取源的残差数据
re = sources['re'] # 获取源的有效半径数据
masked = sources['masked'] # 获取源的掩蔽状态数据
unmasked = np.invert(sources['masked']) # 反转掩蔽状态，获取未掩蔽的源
import numpy as np import matplotlib.pyplot as plt
def fit_line(x, y): “”” 拟合一条直线到数据点 :param x: 自变量数据点 (numpy array) :param y: 因变量数据点 (numpy array) :return: 斜率和截距 (slope, intercept) “”” # 计算斜率和截距 A = np.vstack([x, np.ones(len(x))]).T slope, intercept = np.linalg.lstsq(A, y, rcond=None)[0] return slope, intercept
def fit_and_plot(x, y): “”” 拟合数据点并绘制数据点和最佳拟合直线 :param x: 自变量数据点 (numpy array) :param y: 因变量数据点 (numpy array) “”” # 拟合直线 slope, intercept = fit_line(x, y)
# 生成拟合直线的y值
y_fit = slope * x + intercept
# 绘制数据点
plt.scatter(x, y, color='blue', label='数据点 (Data Points)')
# 绘制拟合直线
plt.plot(x, y_fit, color='red', label='最佳拟合直线 (Best Fit Line)')
# 添加标签和图例
plt.xlabel('自变量 (X)')
plt.ylabel('因变量 (Y)')
plt.title('数据点与最佳拟合直线 (Data Points and Best Fit Line)')
plt.legend()
# 显示图形
plt.show()
说明：#
fit_line函数使用最小二乘法来计算给定数据点的最佳拟合直线的斜率和截距。fit_and_plot函数首先调用fit_line来获取拟合参数，然后绘制数据点和拟合直线，并添加适当的标签和图例。
def fit_line(x, y):
# 创建一个线性模型，初始斜率和截距均为1
line = models.Linear1D(1., 1.)
# 创建一个线性最小二乘拟合器
fit = fitting.LinearLSQFitter()
# 使用拟合器对给定的x和y数据进行拟合，返回最佳拟合结果
best_fit = fit(line, x, y)
return best_fit # 返回最佳拟合模型
def fit_and_plot(x, y, alpha=0.2, color='red', s=10,
marker='o', label=''):
# 生成用于拟合的x值，从0到x的最大值，分成10个点
xfit = np.linspace(0, x.max(), 10)
# 调用fit_line函数拟合x和y数据，返回拟合的线性函数
line = fit_line(x, y)
# 绘制散点图，显示原始数据点
plt.scatter(x, y, alpha=alpha, color=color, s=s, marker=marker,
label=label)
# 绘制拟合线，使用拟合的线性函数和生成的xfit值
plt.plot(xfit, line(xfit), color=color, alpha=0.7)
绘制残差与源通量的关系图，分别针对被遮罩（masked）、未遮罩（unmasked）和随机化（randomized）源位置。将标记的大小与源的通量成正比。理想情况下，这三个情况的残差应以0.0为中心，并且与源通量没有相关性。当然，避免在未遮罩的源位置出现偏差几乎是不可能的，因为它们会影响背景估计。在这种特定情况下，当源未被遮罩时，残差存在小的偏移，并且有证据表明与源通量存在趋势。
plt.figure(figsize=(10, 5)) # 创建一个大小为10x5的图形
# 绘制未掩蔽数据的拟合和散点图，使用黑色方形标记
fit_and_plot(flux[unmasked], resid[unmasked], s=10.*re[unmasked],
label='unmasked', color='black', alpha=0.5, marker='s')
# 绘制已掩蔽数据的拟合和散点图，使用红色圆形标记
fit_and_plot(flux[masked], resid[masked], s=10.*re[masked],
label='masked', color='red', alpha=0.3, marker='o')
# 绘制随机数据的拟合和散点图，使用蓝色加号标记
fit_and_plot(random_flux, resid_under_random, color='blue', alpha=0.1,
marker='+', label='random')
plt.legend() # 显示图例
plt.ylim(-0.05, 0.05) # 设置y轴的范围
plt.title('residuals vs. flux') # 设置图表标题
绘制残差与源半径的关系图，分别针对被遮罩、未被遮罩和随机化的源位置。标记的大小与源的通量成正比。
plt.figure(figsize=(10, 6)) # 创建一个10x6英寸的图形
# 绘制未掩蔽数据的拟合和残差，设置点的大小为2倍的flux值
fit_and_plot(re[unmasked], resid[unmasked], s=2*flux[unmasked],
label='unmasked', color='black', alpha=0.5, marker='s') # 未掩蔽数据，黑色方形标记
# 绘制掩蔽数据的拟合和残差，设置点的大小为flux值
fit_and_plot(re[masked], resid[masked], s=flux[masked],
label='masked', color='red', alpha=0.2, marker='o') # 掩蔽数据，红色圆形标记
# 绘制随机数据的拟合和残差，设置透明度和标记样式
fit_and_plot(random_re, resid_under_random, color='blue', alpha=0.1,
marker='+', label='random') # 随机数据，蓝色加号标记
plt.legend() # 显示图例
plt.ylim(-0.05, 0.05) # 设置y轴的范围
plt.title('residuals vs. radius') # 设置图形标题
在真实数据上测试偏差#
对于真实数据，我们无法相对于无噪声天空来计算残差。然而，我们可以检查被遮罩区域下的背景是否在统计上高于随机区域的背景。考虑到天文源具有翼状结构，而这些结构在没有建模的情况下无法被减去，因此在这里实现完美是非常困难的。这个测试将告诉你潜在偏差的大小（以每像素的通量为单位）以及它是否看起来显著。
mean_masked = bkgd[mask].mean() # 计算被掩膜像素的均值
std_masked = bkgd[mask].std() # 计算被掩膜像素的标准差
stderr_masked = mean_masked/(np.sqrt(len(bkgd[mask]))*std_masked) # 计算被掩膜像素的标准误差
mean_unmasked = bkgd[~mask].mean() # 计算未被掩膜像素的均值
std_unmasked = bkgd[~mask].std() # 计算未被掩膜像素的标准差
stderr_unmasked = mean_unmasked/(np.sqrt(len(bkgd[~mask]))*std_unmasked) # 计算未被掩膜像素的标准误差
diff = mean_masked - mean_unmasked # 计算被掩膜和未被掩膜像素均值的差异
significance = diff/np.sqrt(stderr_masked**2 + stderr_unmasked**2) # 计算显著性
print(f"Mean under masked pixels = {mean_masked:.4f} +- {stderr_masked:.4f}") # 打印被掩膜像素的均值和标准误差
print(f"Mean under unmasked pixels = "
f"{mean_unmasked:.4f} +- {stderr_unmasked:.4f}") # 打印未被掩膜像素的均值和标准误差
print(f"Difference = {diff:.4f} at {significance:.2f} sigma significance") # 打印均值差异及其显著性
便于观察残余背景的RMS与尺度的关系的常规#
评估天空减除是否合理的一种方法是查看RMS是否随着尺度的变化而合理下降。它应该随着采样图像的区域大小线性下降。
为了进行此测试，我们希望查看未被遮罩的区域，但这些区域的信噪比（S/N）应相同。因此，我们需要一种方法来设置一个网格，以确保每个网格区域内没有被遮罩的像素。
def smoothsample_masked(img, mask, width=9):
'''对没有被掩膜的像素进行平方区域的均值计算'''
nrow, ncol = img.shape # 获取图像的行数和列数
w = width / 2 # 计算宽度的一半
# 掩膜掉边界区域
row, col = np.mgrid[0:nrow, 0:ncol] # 创建行列网格
# 更新掩膜，标记出边界区域
mask = mask | (row < w) | (row > nrow - w) | (col < w) | (col > ncol - w)
# 创建一个不包含被掩膜像素的平方区域列表
mm = block_reduce(mask, width) # 对掩膜进行块减少
means = block_reduce(img, width, func=np.mean) # 对图像进行块减少并计算均值
good = np.where(mm == 0) # 找到没有被掩膜的区域
return means[good] # 返回均值数组中未被掩膜区域的均值
def calculate_stats(residual, mask, widths):
stats = [] # 初始化统计数据列表
for w in widths: # 遍历每个宽度
val = smoothsample_masked(residual, mask, w) # 计算平滑样本
stats += [val.std()] # 计算标准差并添加到统计数据列表
return np.array(stats) # 返回标准差的数组
计算不同背景估计的统计数据。这些是经过掩膜的“场景”的残差——即包含源的场景，但经过掩膜处理，就像在实际场景中所做的那样。这些残差将与“完美”情况进行比较，即有噪声的天空背景减去无噪声的天空背景，在相同的未掩膜单元中计算统计数据。
rms = Table() # 创建一个新的表格对象
widths = rms['widths'] = np.arange(3, 30, 1) # 生成一个从3到29的数组，并将其赋值给'rms'表中的'widths'列
rms['bkg1'] = calculate_stats(scene - bkg1img, mask, widths) # 计算背景1的统计数据，并存储在'rms'表中的'bkg1'列
rms['bkg2'] = calculate_stats(scene - bkg2img, mask, widths) # 计算背景2的统计数据，并存储在'rms'表中的'bkg2'列
rms['bkg3'] = calculate_stats(scene - bkg3img, mask, widths) # 计算背景3的统计数据，并存储在'rms'表中的'bkg3'列
rms['bkg4'] = calculate_stats(scene - bkg4img, mask, widths) # 计算背景4的统计数据，并存储在'rms'表中的'bkg4'列
rms['bkg5'] = calculate_stats(scene - bkg5img, mask, widths) # 计算背景5的统计数据，并存储在'rms'表中的'bkg5'列
rms['perfect'] = calculate_stats(sky_bkgd - noiseless_sky, mask, widths) # 计算完美背景的统计数据，并存储在'rms'表中的'perfect'列
绘制相对于完美案例的结果（我们已将初始的 bkg1 估计注释掉，以便更清晰地显示更好估计的比例）。值得注意的是，bkg1 是从 SExtractor 获得的典型结果，SExtractor 仅提供一次背景估计。
plt.figure(figsize=(10, 6)) # 创建一个10x6英寸的图形
perfect = rms['perfect'] # 获取完美背景数据
plt.plot(rms['widths'], rms['bkg1']/perfect, alpha=0.5, label='bkg1') # 绘制bkg1与完美背景的比值
plt.plot(rms['widths'], rms['bkg2']/perfect, alpha=0.5, label='bkg2') # 绘制bkg2与完美背景的比值
plt.plot(rms['widths'], rms['bkg3']/perfect, alpha=0.5, label='bkg3') # 绘制bkg3与完美背景的比值
plt.plot(rms['widths'], rms['bkg4']/perfect, alpha=0.5, label='bkg4') # 绘制bkg4与完美背景的比值
plt.plot(rms['widths'], rms['bkg5']/perfect, alpha=0.5, label='bkg5') # 绘制bkg5与完美背景的比值
plt.plot(rms['widths'], rms['perfect']/perfect, 'k-', alpha=1, label='perfect') # 绘制完美背景的比值，使用黑色实线
plt.legend() # 显示图例
