图像数据

在本章中,我们将讨论图像HDU中的数据组件。

图像数据作为数组

FITS主HDU或图像扩展HDU可能包含图像数据。以下讨论适用于这两个HDU类。在大多数情况下 astropy ,它是一个 numpy 数组,其形状由NAXIS关键字指定,数据类型由BITPIX关键字指定-除非数据被缩放,在这种情况下,请参阅下一节。下面是FITS图像中允许的BITPIX值与 numpy 数据类型:

BITPIX    Numpy Data Type
8         numpy.uint8 (note it is UNsigned integer)
16        numpy.int16
32        numpy.int32
64        numpy.int64
-32       numpy.float32
-64       numpy.float64

总而言之,在 numpy 数组为0索引,轴从慢到快排序。因此,如果FITS图像的NAXIS1=300和NAXIS2=400,则 numpy 其数据数组的形状为(400300)。

实例

以下是读取和更新图像数据值的摘要:

>>> from astropy.io import fits
>>> fits_image_filename = fits.util.get_testdata_filepath('test0.fits')

>>> with fits.open(fits_image_filename) as hdul:  # open a FITS file
...     data = hdul[1].data  # assume the first extension is an image
>>> print(data[1, 4])   # get the pixel value at x=5, y=2
313
>>> # get values of the subsection from x=11 to 20, y=31 to 40 (inclusive)
>>> data[30:40, 10:20]
array([[314, 314, 313, 312, 313, 313, 313, 313, 313, 312],
       [314, 314, 312, 313, 313, 311, 313, 312, 312, 314],
       [314, 315, 313, 313, 313, 313, 315, 312, 314, 312],
       [314, 313, 313, 314, 311, 313, 313, 313, 313, 313],
       [313, 314, 312, 314, 312, 314, 314, 315, 313, 313],
       [312, 311, 311, 312, 312, 312, 312, 313, 311, 312],
       [314, 314, 314, 314, 312, 313, 314, 314, 314, 311],
       [314, 313, 312, 313, 313, 314, 312, 312, 311, 314],
       [313, 313, 313, 314, 313, 313, 315, 313, 312, 313],
       [314, 313, 313, 314, 313, 312, 312, 314, 310, 314]], dtype=int16)
>>> data[1,4] = 999  # update a pixel value
>>> data[30:40, 10:20] = 0  # update values of a subsection
>>> data[3] = data[2]    # copy the 3rd row to the 4th row

这里有一些使用“掩模阵列”概念的更复杂的例子。第一个例子是改变 data 归零。第二种方法是取正像素值的对数:

>>> data[data < 0] = 0
>>> import numpy as np
>>> data[data > 0] = np.log(data[data > 0])

这些例子显示了 numpy 阵列操作。

比例数据

有时会缩放图像;也就是说,存储在文件中的数据不是图像的物理(真)值,而是根据以下等式进行线性变换:

physical value = BSCALE * (storage value) + BZERO

BSCALE和BZERO作为相同名称的关键字存储在同一HDU的头中。缩放图像最常见的用途是存储无符号16位整数数据,因为FITS标准不允许这样做。在这种情况下,存储的数据是有符号的16位整数(BITPIX=16),BZERO=32768 (\(2^{{15}}\) ),b刻度=1。

读取缩放图像数据

只有当标题中存在BSCALE/BZERO关键字之一且它们的值不是默认值(BSCALE=1,BZERO=0)时,才会缩放图像。

对于未标度数据,中HDU的数据属性 astropy 是一个 numpy 由BITPIX关键字指定的相同数据类型的数组。对于缩放图像 .data 属性将是物理数据(即已经从存储数据转换而来,并且可能与BITPIX中规定的数据类型不同)。这意味着需要额外的复制步骤,因此需要相应的内存需求。这也意味着,对于缩放数据,内存映射的优势将降低。

对于浮点存储数据,缩放数据将具有相同的数据类型。对于整数数据类型,缩放数据将始终是单精度浮点 (numpy.float32

例子

下面是一个缩放数据在接触数据之前和之后发生的情况的示例:

>>> fits_scaledimage_filename = fits.util.get_testdata_filepath('scale.fits')

>>> hdul = fits.open(fits_scaledimage_filename)
>>> hdu = hdul[0]
>>> hdu.header['bitpix']
16
>>> hdu.header['bzero']
1500.0
>>> hdu.data[0, 0]  # once data is touched, it is scaled  
557.7563
>>> hdu.data.dtype.name
'float32'
>>> hdu.header['bitpix']  # BITPIX is also updated
-32
>>> # BZERO and BSCALE are removed after the scaling
>>> hdu.header['bzero']
Traceback (most recent call last):
    ...
KeyError: "Keyword 'BZERO' not found."

警告

在中处理缩放数据时需要注意的一个重要注意事项 astropy ,是在通过 .data 属性,则使用BZERO和BSCALE参数自动缩放数据。如果文件是在“更新”模式下打开的,它将与重新缩放的数据一起保存。这种令人惊讶的行为是在不丢失数据的情况下出错的一种折衷:如果对数据进行了一些浮点计算,则在保存时重新调整其大小可能会导致信息丢失。

为了防止文件被自动缩放 do_not_scale_image_data=True 参数 fits.open() . 这对于更新某些头值特别有用,同时确保数据不被修改。

也可以手动重新应用比例参数,方法是使用 hdu.scale() (见下文)。或者,您可以使用 scale_back=True 争论。即使在保存物理值时,也能确保在保存时保留原始缩放比例。换句话说,它在保存时将缩放重新应用到新的物理值。

写入缩放图像数据

由于需要额外的处理和内存,我们尽量不使用缩放数据。然而, astropy 提供了用 scale 方法。

实例

使用 scale 方法:

>>> # scale the data to Int16 with user specified bscale/bzero
>>> hdu.scale('int16', bzero=32768)
>>> # scale the data to Int32 with the min/max of the data range, emits
>>> # RuntimeWarning: overflow encountered in short_scalars
>>> hdu.scale('int32', 'minmax')  
>>> # scale the data, using the original BSCALE/BZERO, emits
>>> # RuntimeWarning: invalid value encountered in add
>>> hdu.scale('int32', 'old')  
>>> hdul.close()

上面的第一个示例演示了如何存储无符号短整数数组。

在使用 scale() 方法。这个 data 图像HDU的属性,位于 scale() 存储不会变成物理值。所以,只打电话 scale() 在写入FITS文件之前(即 writeto()flush()close() ). 不应进一步使用数据。下面是一个示例,说明 data 属性之后 scale() 呼叫:

>>> hdu = fits.PrimaryHDU(np.array([0., 1, 2, 3]))
>>> print(hdu.data)  
[0. 1. 2. 3.]
>>> hdu.scale('int16', bzero=32768)
>>> print(hdu.data)  # now the data has storage values
[-32768 -32767 -32766 -32765]
>>> hdu.writeto('new.fits')

数据段

当a适合图像HDU data 在数据被复制到内存的情况下(或者在数据被复制到内存的情况下)。如果同时访问多个非常大的图像hdu,系统可能会耗尽内存。

如果用户不需要同时处理整个图像(例如,一次处理十行图像),则 section HDU的属性可以用来缓解这种内存问题。

astropy 改进了对内存映射的支持,sections功能不再像以前那样需要处理非常大的图像。但是,如果图像的数据是用非平凡的b缩放/b零值来缩放的,则在当前的实现中可能仍然需要分段访问数据。Memmap也不足以在32位系统上加载大于2到4 GB的映像—在这种情况下,可能需要使用节。

例子

下面是一个从大小为5000x5000的三个输入图像中获取中值图像的示例。

hdul1 = fits.open('file1.fits')
hdul2 = fits.open('file2.fits')
hdul3 = fits.open('file3.fits')
output = np.zeros((5000, 5000))
for i in range(50):
    j = i * 100
    k = j + 100
    x1 = hdul1[0].section[j:k,:]
    x2 = hdul2[0].section[j:k,:]
    x3 = hdul3[0].section[j:k,:]
    output[j:k, :] = np.median([x1, x2, x3], axis=0)

每个 section 不需要是连续的才能节省内存。 astropy 将尽最大努力将数组中不连续的部分连接在一起,同时尽可能少地读入主内存。

当前无法分配分区。对数据节所做的任何修改都不会保存回原始文件。