用椭圆圈出天文图像中的物体

发表于 | 分类于 |
字数统计: 374

内容如题，结果如下：

所用工具：astropy，photutils包

image-20211029174406755

代码：

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
import numpy as np
from astropy.io import fits
import cv2
import matplotlib.pyplot as plt
from photutils.segmentation import detect_threshold
from astropy.convolution import Gaussian2DKernel
from astropy.stats import gaussian_fwhm_to_sigma
from photutils.segmentation import detect_sources, deblend_sources, SourceCatalog,source_properties
import warnings
from astropy.visualization import simple_norm
warnings.filterwarnings("ignore")


def plot_circle(img_rgb, img_w):
    """
    将图像中的物体用椭圆圈出来
    :param img_rgb: 伪彩图
    :param img_w: 单通道图
    :return: 圈出来的
    """
    img_rgb = cv2.imread(img_rgb)[:,:,::-1]
    # cv2.imshow('w',img_rgb)
    # cv2.waitKey(0)
    hdu = fits.open(img_w)
    data = hdu[0].data
    data = data[::-1, :]
    data = np.array(data)

    # nsigma越大，越能忽略掉小物体
    threshold = detect_threshold(data, nsigma=5.)
    sigma = 3.0 * gaussian_fwhm_to_sigma
    # print(sigma)
    kernal = Gaussian2DKernel(sigma, x_size=5, y_size=5)
    kernal.normalize()
    npixels = 5  # 最小物体的像素数
    # 检测物体
    segm = detect_sources(data, threshold, npixels=5, filter_kernel=kernal, connectivity=4)
    # 将两个连一块被判成一个的物体分割开
    segm_deblend = deblend_sources(data, segm, npixels=npixels, filter_kernel=kernal, nlevels=1,
                                   contrast=0.001)

    # cmap = segm_deblend.make_cmap(seed=123)

    # 用椭圆圈
    cat = source_properties(data, segm)
    a = cat.kron_aperture  # 圈出的椭圆
    # tbl = cat.to_table()

    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 12.5))
    norm = simple_norm(data, 'sqrt')
    ax1.imshow(img_rgb)
    ax1.set_title('Data')

    for aperture in a:
        # 里面会有一些None
        try:
            # 画出椭圆
            aperture.plot(axes=ax2,color='white', lw=0.6)
        except:
            pass
    ax2.imshow(img_rgb, )
    ax2.set_title('Segmentation Image')
    plt.show()

# 要用r波段的，其他波段会有偏差，因为伪彩图十一r波段为基准的
plot_circle(r'img\1237648705133805706.jpg',
            r'img\1237648705133805706-r.fits.bz2', )

这里是在r波段圈的椭圆，直接在伪彩图上也能圈，需要将其转成灰度图。

伪彩图转成灰度图圈出来的没有在r波段上圈的效果好。