/astropy/io/fits/hdu/compressed.py
Python | 1133 lines | 804 code | 152 blank | 177 comment | 204 complexity | 6f517256587bde95146c9baf1fc92a38 MD5 | raw file
- # Licensed under a 3-clause BSD style license - see PYFITS.rst
- import ctypes
- import gc
- import math
- import re
- import time
- import warnings
- import numpy as np
- from .base import DELAYED, ExtensionHDU
- from .image import _ImageBaseHDU, ImageHDU
- from .table import BinTableHDU
- from ..card import Card
- from ..column import Column, ColDefs, TDEF_RE
- from ..column import KEYWORD_NAMES as TABLE_KEYWORD_NAMES
- from ..fitsrec import FITS_rec
- from ..header import Header
- from ..util import (_is_pseudo_unsigned, _unsigned_zero, _is_int,
- _get_array_mmap)
- from ....extern.six import string_types, iteritems
- from ....utils import lazyproperty, deprecated
- from ....utils.compat import ignored
- from ....utils.exceptions import (AstropyPendingDeprecationWarning,
- AstropyUserWarning)
- try:
- from .. import compression
- COMPRESSION_SUPPORTED = COMPRESSION_ENABLED = True
- except ImportError:
- COMPRESSION_SUPPORTED = COMPRESSION_ENABLED = False
- # Quantization dithering method constants; these are right out of fitsio.h
- NO_DITHER = -1
- SUBTRACTIVE_DITHER_1 = 1
- SUBTRACTIVE_DITHER_2 = 2
- QUANTIZE_METHOD_NAMES = {
- NO_DITHER: 'NO_DITHER',
- SUBTRACTIVE_DITHER_1: 'SUBTRACTIVE_DITHER_1',
- SUBTRACTIVE_DITHER_2: 'SUBTRACTIVE_DITHER_2'
- }
- DITHER_SEED_CLOCK = 0
- DITHER_SEED_CHECKSUM = -1
- # Default compression parameter values
- DEFAULT_COMPRESSION_TYPE = 'RICE_1'
- DEFAULT_QUANTIZE_LEVEL = 16.
- DEFAULT_QUANTIZE_METHOD = NO_DITHER
- DEFAULT_DITHER_SEED = DITHER_SEED_CLOCK
- DEFAULT_HCOMP_SCALE = 0
- DEFAULT_HCOMP_SMOOTH = 0
- DEFAULT_BLOCK_SIZE = 32
- DEFAULT_BYTE_PIX = 4
- CMTYPE_ALIASES = {}
- # CFITSIO version-specific features
- if COMPRESSION_SUPPORTED:
- try:
- CFITSIO_SUPPORTS_GZIPDATA = compression.CFITSIO_VERSION >= 3.28
- CFITSIO_SUPPORTS_Q_FORMAT = compression.CFITSIO_VERSION >= 3.35
- if compression.CFITSIO_VERSION >= 3.35:
- CMTYPE_ALIASES['RICE_ONE'] = 'RICE_1'
- except AttributeError:
- # This generally shouldn't happen unless running setup.py in an
- # environment where an old build of pyfits exists
- CFITSIO_SUPPORTS_GZIPDATA = True
- CFITSIO_SUPPORTS_Q_FORMAT = True
- COMPRESSION_KEYWORDS = set(['ZIMAGE', 'ZCMPTYPE', 'ZBITPIX', 'ZNAXIS',
- 'ZMASKCMP', 'ZSIMPLE', 'ZTENSION', 'ZEXTEND'])
- class CompImageHeader(Header):
- """
- Header object for compressed image HDUs designed to keep the compression
- header and the underlying image header properly synchronized.
- This essentially wraps the image header, so that all values are read from
- and written to the image header. However, updates to the image header will
- also update the table header where appropriate.
- """
- # TODO: The difficulty of implementing this screams a need to rewrite this
- # module
- _keyword_remaps = {
- 'SIMPLE': 'ZSIMPLE', 'XTENSION': 'ZTENSION', 'BITPIX': 'ZBITPIX',
- 'NAXIS': 'ZNAXIS', 'EXTEND': 'ZEXTEND', 'BLOCKED': 'ZBLOCKED',
- 'PCOUNT': 'ZPCOUNT', 'GCOUNT': 'ZGCOUNT', 'CHECKSUM': 'ZHECKSUM',
- 'DATASUM': 'ZDATASUM'
- }
- _zdef_re = re.compile(r'(?P<label>^[Zz][a-zA-Z]*)(?P<num>[1-9][0-9 ]*$)?')
- _compression_keywords = set(_keyword_remaps.values()).union(
- ['ZIMAGE', 'ZCMPTYPE', 'ZMASKCMP', 'ZQUANTIZ', 'ZDITHER0'])
- _indexed_compression_keywords = set(['ZNAXIS', 'ZTILE', 'ZNAME', 'ZVAL'])
- # TODO: Once it place it should be possible to manage some of this through
- # the schema system, but it's not quite ready for that yet. Also it still
- # makes more sense to change CompImageHDU to subclass ImageHDU :/
- def __init__(self, table_header, image_header=None):
- if image_header is None:
- image_header = Header()
- self._cards = image_header._cards
- self._keyword_indices = image_header._keyword_indices
- self._rvkc_indices = image_header._rvkc_indices
- self._modified = image_header._modified
- self._table_header = table_header
- # We need to override and Header methods that can modify the header, and
- # ensure that they sync with the underlying _table_header
- def __setitem__(self, key, value):
- # This isn't pretty, but if the `key` is either an int or a tuple we
- # need to figure out what keyword name that maps to before doing
- # anything else; these checks will be repeated later in the
- # super().__setitem__ call but I don't see another way around it
- # without some major refactoring
- if self._set_slice(key, value, self):
- return
- if isinstance(key, int):
- keyword, index = self._keyword_from_index(key)
- elif isinstance(key, tuple):
- keyword, index = key
- else:
- # We don't want to specify and index otherwise, because that will
- # break the behavior for new keywords and for commentary keywords
- keyword, index = key, None
- if self._is_reserved_keyword(keyword):
- return
- super(CompImageHeader, self).__setitem__(key, value)
- if index is not None:
- remapped_keyword = self._remap_keyword(keyword)
- self._table_header[remapped_keyword, index] = value
- # Else this will pass through to ._update
- def __delitem__(self, key):
- if isinstance(key, slice) or self._haswildcard(key):
- # If given a slice pass that on to the superclass and bail out
- # early; we only want to make updates to _table_header when given
- # a key specifying a single keyword
- return super(CompImageHeader, self).__delitem__(key)
- if isinstance(key, int):
- keyword, index = self._keyword_from_index(key)
- elif isinstance(key, tuple):
- keyword, index = key
- else:
- keyword, index = key, None
- if key not in self:
- raise KeyError("Keyword %r not found." % key)
- super(CompImageHeader, self).__delitem__(key)
- remapped_keyword = self._remap_keyword(keyword)
- if remapped_keyword in self._table_header:
- if index is not None:
- del self._table_header[(remapped_keyword, index)]
- else:
- del self._table_header[remapped_keyword]
- def append(self, card=None, useblanks=True, bottom=False, end=False):
- # This logic unfortunately needs to be duplicated from the base class
- # in order to determine the keyword
- if isinstance(card, string_types):
- card = Card(card)
- elif isinstance(card, tuple):
- card = Card(*card)
- elif card is None:
- card = Card()
- elif not isinstance(card, Card):
- raise ValueError(
- 'The value appended to a Header must be either a keyword or '
- '(keyword, value, [comment]) tuple; got: %r' % card)
- if self._is_reserved_keyword(card.keyword):
- return
- super(CompImageHeader, self).append(card=card, useblanks=useblanks,
- bottom=bottom, end=end)
- remapped_keyword = self._remap_keyword(card.keyword)
- card = Card(remapped_keyword, card.value, card.comment)
- self._table_header.append(card=card, useblanks=useblanks,
- bottom=bottom, end=end)
- def insert(self, key, card, useblanks=True, after=False):
- if isinstance(key, int):
- # Determine condition to pass through to append
- if after:
- if key == -1:
- key = len(self._cards)
- else:
- key += 1
- if key >= len(self._cards):
- self.append(card, end=True)
- return
- if isinstance(card, string_types):
- card = Card(card)
- elif isinstance(card, tuple):
- card = Card(*card)
- elif not isinstance(card, Card):
- raise ValueError(
- 'The value inserted into a Header must be either a keyword or '
- '(keyword, value, [comment]) tuple; got: %r' % card)
- if self._is_reserved_keyword(card.keyword):
- return
- # Now the tricky part is to determine where to insert in the table
- # header. If given a numerical index we need to map that to the
- # corresponding index in the table header. Although rare, there may be
- # cases where there is no mapping in which case we just try the same
- # index
- # NOTE: It is crucial that remapped_index in particular is figured out
- # before the image header is modified
- remapped_index = self._remap_index(key)
- remapped_keyword = self._remap_keyword(card.keyword)
- super(CompImageHeader, self).insert(key, card, useblanks=useblanks,
- after=after)
- card = Card(remapped_keyword, card.value, card.comment)
- self._table_header.insert(remapped_index, card, useblanks=useblanks,
- after=after)
- def _update(self, card):
- keyword = card[0]
- if self._is_reserved_keyword(keyword):
- return
- super(CompImageHeader, self)._update(card)
- if keyword in Card._commentary_keywords:
- # Otherwise this will result in a duplicate insertion
- return
- remapped_keyword = self._remap_keyword(keyword)
- self._table_header._update((remapped_keyword,) + card[1:])
- # Last piece needed (I think) for synchronizing with the real header
- # This one is tricky since _relativeinsert calls insert
- def _relativeinsert(self, card, before=None, after=None, replace=False):
- keyword = card[0]
- if self._is_reserved_keyword(keyword):
- return
- # Now we have to figure out how to remap 'before' and 'after'
- if before is None:
- if isinstance(after, int):
- remapped_after = self._remap_index(after)
- else:
- remapped_after = self._remap_keyword(after)
- remapped_before = None
- else:
- if isinstance(before, int):
- remapped_before = self._remap_index(before)
- else:
- remapped_before = self._remap_keyword(before)
- remapped_after = None
- super(CompImageHeader, self)._relativeinsert(card, before=before,
- after=after,
- replace=replace)
- remapped_keyword = self._remap_keyword(keyword)
- card = Card(remapped_keyword, card[1], card[2])
- self._table_header._relativeinsert(card, before=remapped_before,
- after=remapped_after,
- replace=replace)
- @classmethod
- def _is_reserved_keyword(cls, keyword, warn=True):
- msg = ('Keyword %r is reserved for use by the FITS Tiled Image '
- 'Convention and will not be stored in the header for the '
- 'image being compressed.' % keyword)
- if keyword == 'TFIELDS':
- if warn:
- warnings.warn(msg)
- return True
- m = TDEF_RE.match(keyword)
- if m and m.group('label').upper() in TABLE_KEYWORD_NAMES:
- if warn:
- warnings.warn(msg)
- return True
- m = cls._zdef_re.match(keyword)
- if m:
- label = m.group('label').upper()
- num = m.group('num')
- if num is not None and label in cls._indexed_compression_keywords:
- if warn:
- warnings.warn(msg)
- return True
- elif label in cls._compression_keywords:
- if warn:
- warnings.warn(msg)
- return True
- return False
- @classmethod
- def _remap_keyword(cls, keyword):
- # Given a keyword that one might set on an image, remap that keyword to
- # the name used for it in the COMPRESSED HDU header
- # This is mostly just a lookup in _keyword_remaps, but needs handling
- # for NAXISn keywords
- is_naxisn = False
- if keyword[:5] == 'NAXIS':
- with ignored(ValueError):
- index = int(keyword[5:])
- is_naxisn = index > 0
- if is_naxisn:
- return 'ZNAXIS%d' % index
- # If the keyword does not need to be remapped then just return the
- # original keyword
- return cls._keyword_remaps.get(keyword, keyword)
- def _remap_index(self, idx):
- # Given an integer index into this header, map that to the index in the
- # table header for the same card. If the card doesn't exist in the
- # table header (generally should *not* be the case) this will just
- # return the same index
- # This *does* also accept a keyword or (keyword, repeat) tuple and
- # obtains the associated numerical index with self._cardindex
- if not isinstance(idx, int):
- idx = self._cardindex(idx)
- keyword, repeat = self._keyword_from_index(idx)
- remapped_insert_keyword = self._remap_keyword(keyword)
- with ignored(IndexError, KeyError):
- idx = self._table_header._cardindex((remapped_insert_keyword,
- repeat))
- return idx
- # TODO: Fix this class so that it doesn't actually inherit from BinTableHDU,
- # but instead has an internal BinTableHDU reference
- class CompImageHDU(BinTableHDU):
- """
- Compressed Image HDU class.
- """
- # Maps deprecated keyword arguments to __init__ to their new names
- DEPRECATED_KWARGS = {
- 'compressionType': 'compression_type', 'tileSize': 'tile_size',
- 'hcompScale': 'hcomp_scale', 'hcompSmooth': 'hcomp_smooth',
- 'quantizeLevel': 'quantize_level'
- }
- _manages_own_heap = True
- """
- The calls to CFITSIO lay out the heap data in memory, and we write it out
- the same way CFITSIO organizes it. In principle this would break if a user
- manually changes the underlying compressed data by hand, but there is no
- reason they would want to do that (and if they do that's their
- responsibility).
- """
- def __init__(self, data=None, header=None, name=None,
- compression_type=DEFAULT_COMPRESSION_TYPE,
- tile_size=None,
- hcomp_scale=DEFAULT_HCOMP_SCALE,
- hcomp_smooth=DEFAULT_HCOMP_SMOOTH,
- quantize_level=DEFAULT_QUANTIZE_LEVEL,
- quantize_method=DEFAULT_QUANTIZE_METHOD,
- dither_seed=DEFAULT_DITHER_SEED,
- do_not_scale_image_data=False,
- uint=False, scale_back=False, **kwargs):
- """
- Parameters
- ----------
- data : array, optional
- Uncompressed image data
- header : Header instance, optional
- Header to be associated with the image; when reading the HDU from a
- file (data=DELAYED), the header read from the file
- name : str, optional
- The ``EXTNAME`` value; if this value is `None`, then the name from
- the input image header will be used; if there is no name in the
- input image header then the default name ``COMPRESSED_IMAGE`` is
- used.
- compression_type : str, optional
- Compression algorithm: one of
- ``'RICE_1'``, ``'RICE_ONE'``, ``'PLIO_1'``, ``'GZIP_1'``,
- ``'GZIP_2'``, ``'HCOMPRESS_1'``
- tile_size : int, optional
- Compression tile sizes. Default treats each row of image as a
- tile.
- hcomp_scale : float, optional
- HCOMPRESS scale parameter
- hcomp_smooth : float, optional
- HCOMPRESS smooth parameter
- quantize_level : float, optional
- Floating point quantization level; see note below
- quantize_method : int, optional
- Floating point quantization dithering method; can be either
- ``NO_DITHER`` (-1), ``SUBTRACTIVE_DITHER_1`` (1; default), or
- ``SUBTRACTIVE_DITHER_2`` (2); see note below
- dither_seed : int, optional
- Random seed to use for dithering; can be either an integer in the
- range 1 to 1000 (inclusive), ``DITHER_SEED_CLOCK`` (0; default), or
- ``DITHER_SEED_CHECKSUM`` (-1); see note below
- Notes
- -----
- The astropy.io.fits package supports 2 methods of image compression:
- 1) The entire FITS file may be externally compressed with the gzip
- or pkzip utility programs, producing a ``*.gz`` or ``*.zip``
- file, respectively. When reading compressed files of this type,
- Astropy first uncompresses the entire file into a temporary file
- before performing the requested read operations. The
- astropy.io.fits package does not support writing to these types
- of compressed files. This type of compression is supported in
- the ``_File`` class, not in the `CompImageHDU` class. The file
- compression type is recognized by the ``.gz`` or ``.zip`` file
- name extension.
- 2) The `CompImageHDU` class supports the FITS tiled image
- compression convention in which the image is subdivided into a
- grid of rectangular tiles, and each tile of pixels is
- individually compressed. The details of this FITS compression
- convention are described at the `FITS Support Office web site
- <http://fits.gsfc.nasa.gov/registry/tilecompression.html>`_.
- Basically, the compressed image tiles are stored in rows of a
- variable length array column in a FITS binary table. The
- astropy.io.fits recognizes that this binary table extension
- contains an image and treats it as if it were an image
- extension. Under this tile-compression format, FITS header
- keywords remain uncompressed. At this time, Astropy does not
- support the ability to extract and uncompress sections of the
- image without having to uncompress the entire image.
- The astropy.io.fits package supports 3 general-purpose compression
- algorithms plus one other special-purpose compression technique that is
- designed for data masks with positive integer pixel values. The 3
- general purpose algorithms are GZIP, Rice, and HCOMPRESS, and the
- special-purpose technique is the IRAF pixel list compression technique
- (PLIO). The ``compression_type`` parameter defines the compression
- algorithm to be used.
- The FITS image can be subdivided into any desired rectangular grid of
- compression tiles. With the GZIP, Rice, and PLIO algorithms, the
- default is to take each row of the image as a tile. The HCOMPRESS
- algorithm is inherently 2-dimensional in nature, so the default in this
- case is to take 16 rows of the image per tile. In most cases, it makes
- little difference what tiling pattern is used, so the default tiles are
- usually adequate. In the case of very small images, it could be more
- efficient to compress the whole image as a single tile. Note that the
- image dimensions are not required to be an integer multiple of the tile
- dimensions; if not, then the tiles at the edges of the image will be
- smaller than the other tiles. The ``tile_size`` parameter may be
- provided as a list of tile sizes, one for each dimension in the image.
- For example a ``tile_size`` value of ``[100,100]`` would divide a 300 X
- 300 image into 9 100 X 100 tiles.
- The 4 supported image compression algorithms are all 'lossless' when
- applied to integer FITS images; the pixel values are preserved exactly
- with no loss of information during the compression and uncompression
- process. In addition, the HCOMPRESS algorithm supports a 'lossy'
- compression mode that will produce larger amount of image compression.
- This is achieved by specifying a non-zero value for the ``hcomp_scale``
- parameter. Since the amount of compression that is achieved depends
- directly on the RMS noise in the image, it is usually more convenient
- to specify the ``hcomp_scale`` factor relative to the RMS noise.
- Setting ``hcomp_scale = 2.5`` means use a scale factor that is 2.5
- times the calculated RMS noise in the image tile. In some cases it may
- be desirable to specify the exact scaling to be used, instead of
- specifying it relative to the calculated noise value. This may be done
- by specifying the negative of the desired scale value (typically in the
- range -2 to -100).
- Very high compression factors (of 100 or more) can be achieved by using
- large ``hcomp_scale`` values, however, this can produce undesirable
- 'blocky' artifacts in the compressed image. A variation of the
- HCOMPRESS algorithm (called HSCOMPRESS) can be used in this case to
- apply a small amount of smoothing of the image when it is uncompressed
- to help cover up these artifacts. This smoothing is purely cosmetic
- and does not cause any significant change to the image pixel values.
- Setting the ``hcomp_smooth`` parameter to 1 will engage the smoothing
- algorithm.
- Floating point FITS images (which have ``BITPIX`` = -32 or -64) usually
- contain too much 'noise' in the least significant bits of the mantissa
- of the pixel values to be effectively compressed with any lossless
- algorithm. Consequently, floating point images are first quantized
- into scaled integer pixel values (and thus throwing away much of the
- noise) before being compressed with the specified algorithm (either
- GZIP, RICE, or HCOMPRESS). This technique produces much higher
- compression factors than simply using the GZIP utility to externally
- compress the whole FITS file, but it also means that the original
- floating point value pixel values are not exactly preserved. When done
- properly, this integer scaling technique will only discard the
- insignificant noise while still preserving all the real information in
- the image. The amount of precision that is retained in the pixel
- values is controlled by the ``quantize_level`` parameter. Larger
- values will result in compressed images whose pixels more closely match
- the floating point pixel values, but at the same time the amount of
- compression that is achieved will be reduced. Users should experiment
- with different values for this parameter to determine the optimal value
- that preserves all the useful information in the image, without
- needlessly preserving all the 'noise' which will hurt the compression
- efficiency.
- The default value for the ``quantize_level`` scale factor is 16, which
- means that scaled integer pixel values will be quantized such that the
- difference between adjacent integer values will be 1/16th of the noise
- level in the image background. An optimized algorithm is used to
- accurately estimate the noise in the image. As an example, if the RMS
- noise in the background pixels of an image = 32.0, then the spacing
- between adjacent scaled integer pixel values will equal 2.0 by default.
- Note that the RMS noise is independently calculated for each tile of
- the image, so the resulting integer scaling factor may fluctuate
- slightly for each tile. In some cases, it may be desirable to specify
- the exact quantization level to be used, instead of specifying it
- relative to the calculated noise value. This may be done by specifying
- the negative of desired quantization level for the value of
- ``quantize_level``. In the previous example, one could specify
- ``quantize_level = -2.0`` so that the quantized integer levels differ
- by 2.0. Larger negative values for ``quantize_level`` means that the
- levels are more coarsely-spaced, and will produce higher compression
- factors.
- The quantization algorithm can also apply one of two random dithering
- methods in order to reduce bias in the measured intensity of background
- regions. The default method, specified with the constant
- ``SUBTRACTIVE_DITHER_1`` adds dithering to the zero-point of the
- quantization array itself rather than adding noise to the actual image.
- The random noise is added on a pixel-by-pixel basis, so in order
- restore each pixel from its integer value to its floating point value
- it is necessary to replay the same sequence of random numbers for each
- pixel (see below). The other method, ``SUBTRACTIVE_DITHER_2``, is
- exactly like the first except that before dithering any pixel with a
- floating point value of ``0.0`` is replaced with the special integer
- value ``-2147483647``. When the image is uncompressed, pixels with
- this value are restored back to ``0.0`` exactly. Finally, a value of
- ``NO_DITHER`` disables dithering entirely.
- As mentioned above, when using the subtractive dithering algorithm it
- is necessary to be able to generate a (pseudo-)random sequence of noise
- for each pixel, and replay that same sequence upon decompressing. To
- facilitate this, a random seed between 1 and 10000 (inclusive) is used
- to seed a random number generator, and that seed is stored in the
- ``ZDITHER0`` keyword in the header of the compressed HDU. In order to
- use that seed to generate the same sequence of random numbers the same
- random number generator must be used at compression and decompression
- time; for that reason the tiled image convention provides an
- implementation of a very simple pseudo-random number generator. The
- seed itself can be provided in one of three ways, controllable by the
- ``dither_seed`` argument: It may be specified manually, or it may be
- generated arbitrarily based on the system's clock
- (``DITHER_SEED_CLOCK``) or based on a checksum of the pixels in the
- image's first tile (``DITHER_SEED_CHECKSUM``). The clock-based method
- is the default, and is sufficient to ensure that the value is
- reasonably "arbitrary" and that the same seed is unlikely to be
- generated sequentially. The checksum method, on the other hand,
- ensures that the same seed is used every time for a specific image.
- This is particularly useful for software testing as it ensures that the
- same image will always use the same seed.
- """
- if not COMPRESSION_SUPPORTED:
- # TODO: Raise a more specific Exception type
- raise Exception('The astropy.io.fits.compression module is not '
- 'available. Creation of compressed image HDUs is '
- 'disabled.')
- compression_type = CMTYPE_ALIASES.get(compression_type, compression_type)
- # Handle deprecated keyword arguments
- compression_opts = {}
- for oldarg, newarg in self.DEPRECATED_KWARGS.items():
- if oldarg in kwargs:
- warnings.warn('Keyword argument %s to %s is pending '
- 'deprecation; use %s instead' %
- (oldarg, self.__class__.__name__, newarg),
- AstropyPendingDeprecationWarning)
- compression_opts[newarg] = kwargs[oldarg]
- del kwargs[oldarg]
- else:
- compression_opts[newarg] = locals()[newarg]
- # Include newer compression options that don't required backwards
- # compatibility with deprecated spellings
- compression_opts['quantize_method'] = quantize_method
- compression_opts['dither_seed'] = dither_seed
- if data is DELAYED:
- # Reading the HDU from a file
- super(CompImageHDU, self).__init__(data=data, header=header)
- else:
- # Create at least a skeleton HDU that matches the input
- # header and data (if any were input)
- super(CompImageHDU, self).__init__(data=None, header=header)
- # Store the input image data
- self.data = data
- # Update the table header (_header) to the compressed
- # image format and to match the input data (if any);
- # Create the image header (_image_header) from the input
- # image header (if any) and ensure it matches the input
- # data; Create the initially empty table data array to
- # hold the compressed data.
- self._update_header_data(header, name, **compression_opts)
- # TODO: A lot of this should be passed on to an internal image HDU o
- # something like that, see ticket #88
- self._do_not_scale_image_data = do_not_scale_image_data
- self._uint = uint
- self._scale_back = scale_back
- self._axes = [self._header.get('ZNAXIS' + str(axis + 1), 0)
- for axis in range(self._header.get('ZNAXIS', 0))]
- # store any scale factors from the table header
- if do_not_scale_image_data:
- self._bzero = 0
- self._bscale = 1
- else:
- self._bzero = self._header.get('BZERO', 0)
- self._bscale = self._header.get('BSCALE', 1)
- self._bitpix = self._header['ZBITPIX']
- self._orig_bzero = self._bzero
- self._orig_bscale = self._bscale
- self._orig_bitpix = self._bitpix
- @classmethod
- def match_header(cls, header):
- card = header.cards[0]
- if card.keyword != 'XTENSION':
- return False
- xtension = card.value
- if isinstance(xtension, string_types):
- xtension = xtension.rstrip()
- if xtension not in ('BINTABLE', 'A3DTABLE'):
- return False
- if 'ZIMAGE' not in header or header['ZIMAGE'] != True:
- return False
- if COMPRESSION_SUPPORTED and COMPRESSION_ENABLED:
- return True
- elif not COMPRESSION_SUPPORTED:
- warnings.warn('Failure matching header to a compressed image '
- 'HDU: The compression module is not available.\n'
- 'The HDU will be treated as a Binary Table HDU.',
- AstropyUserWarning)
- return False
- else:
- # Compression is supported but disabled; just pass silently (#92)
- return False
- def _update_header_data(self, image_header,
- name=None,
- compression_type=None,
- tile_size=None,
- hcomp_scale=None,
- hcomp_smooth=None,
- quantize_level=None,
- quantize_method=None,
- dither_seed=None):
- """
- Update the table header (`_header`) to the compressed
- image format and to match the input data (if any). Create
- the image header (`_image_header`) from the input image
- header (if any) and ensure it matches the input
- data. Create the initially-empty table data array to hold
- the compressed data.
- This method is mainly called internally, but a user may wish to
- call this method after assigning new data to the `CompImageHDU`
- object that is of a different type.
- Parameters
- ----------
- image_header : Header instance
- header to be associated with the image
- name : str, optional
- the ``EXTNAME`` value; if this value is `None`, then the name from
- the input image header will be used; if there is no name in the
- input image header then the default name 'COMPRESSED_IMAGE' is used
- compression_type : str, optional
- compression algorithm 'RICE_1', 'PLIO_1', 'GZIP_1', 'HCOMPRESS_1';
- if this value is `None`, use value already in the header; if no
- value already in the header, use 'RICE_1'
- tile_size : sequence of int, optional
- compression tile sizes as a list; if this value is `None`, use
- value already in the header; if no value already in the header,
- treat each row of image as a tile
- hcomp_scale : float, optional
- HCOMPRESS scale parameter; if this value is `None`, use the value
- already in the header; if no value already in the header, use 1
- hcomp_smooth : float, optional
- HCOMPRESS smooth parameter; if this value is `None`, use the value
- already in the header; if no value already in the header, use 0
- quantize_level : float, optional
- floating point quantization level; if this value is `None`, use the
- value already in the header; if no value already in header, use 16
- quantize_method : int, optional
- floating point quantization dithering method; can be either
- NO_DITHER (-1), SUBTRACTIVE_DITHER_1 (1; default), or
- SUBTRACTIVE_DITHER_2 (2)
- dither_seed : int, optional
- random seed to use for dithering; can be either an integer in the
- range 1 to 1000 (inclusive), DITHER_SEED_CLOCK (0; default), or
- DITHER_SEED_CHECKSUM (-1)
- """
- image_hdu = ImageHDU(data=self.data, header=self._header)
- self._image_header = CompImageHeader(self._header, image_hdu.header)
- self._axes = image_hdu._axes
- del image_hdu
- # Determine based on the size of the input data whether to use the Q
- # column format to store compressed data or the P format.
- # The Q format is used only if the uncompressed data is larger than
- # 4 GB. This is not a perfect heuristic, as one can contrive an input
- # array which, when compressed, the entire binary table representing
- # the compressed data is larger than 4GB. That said, this is the same
- # heuristic used by CFITSIO, so this should give consistent results.
- # And the cases where this heuristic is insufficient are extreme and
- # almost entirely contrived corner cases, so it will do for now
- if self._has_data:
- huge_hdu = self.data.nbytes > 2 ** 32
- if huge_hdu and not CFITSIO_SUPPORTS_Q_FORMAT:
- raise IOError(
- "Astropy cannot compress images greater than 4 GB in size "
- "(%s is %s bytes) without CFITSIO >= 3.35" %
- ((self.name, self.ver), self.data.nbytes))
- else:
- huge_hdu = False
- # Update the extension name in the table header
- if not name and not 'EXTNAME' in self._header:
- name = 'COMPRESSED_IMAGE'
- if name:
- self._header.set('EXTNAME', name,
- 'name of this binary table extension',
- after='TFIELDS')
- self.name = name
- else:
- self.name = self._header['EXTNAME']
- # Set the compression type in the table header.
- if compression_type:
- if compression_type not in ['RICE_1', 'GZIP_1', 'PLIO_1',
- 'HCOMPRESS_1']:
- warnings.warn('Unknown compression type provided. Default '
- '(%s) compression used.' %
- DEFAULT_COMPRESSION_TYPE, AstropyUserWarning)
- compression_type = DEFAULT_COMPRESSION_TYPE
- self._header.set('ZCMPTYPE', compression_type,
- 'compression algorithm', after='TFIELDS')
- else:
- compression_type = self._header.get('ZCMPTYPE',
- DEFAULT_COMPRESSION_TYPE)
- compression_type = CMTYPE_ALIASES.get(compression_type,
- compression_type)
- # If the input image header had BSCALE/BZERO cards, then insert
- # them in the table header.
- if image_header:
- bzero = image_header.get('BZERO', 0.0)
- bscale = image_header.get('BSCALE', 1.0)
- after_keyword = 'EXTNAME'
- if bscale != 1.0:
- self._header.set('BSCALE', bscale, after=after_keyword)
- after_keyword = 'BSCALE'
- if bzero != 0.0:
- self._header.set('BZERO', bzero, after=after_keyword)
- bitpix_comment = image_header.comments['BITPIX']
- naxis_comment = image_header.comments['NAXIS']
- else:
- bitpix_comment = 'data type of original image'
- naxis_comment = 'dimension of original image'
- # Set the label for the first column in the table
- self._header.set('TTYPE1', 'COMPRESSED_DATA', 'label for field 1',
- after='TFIELDS')
- # Set the data format for the first column. It is dependent
- # on the requested compression type.
- if compression_type == 'PLIO_1':
- tform1 = '1QI' if huge_hdu else '1PI'
- else:
- tform1 = '1QB' if huge_hdu else '1PB'
- self._header.set('TFORM1', tform1,
- 'data format of field: variable length array',
- after='TTYPE1')
- # Create the first column for the table. This column holds the
- # compressed data.
- col1 = Column(name=self._header['TTYPE1'], format=tform1)
- # Create the additional columns required for floating point
- # data and calculate the width of the output table.
- zbitpix = self._image_header['BITPIX']
- if zbitpix < 0 and quantize_level != 0.0:
- # floating point image has 'COMPRESSED_DATA',
- # 'UNCOMPRESSED_DATA', 'ZSCALE', and 'ZZERO' columns (unless using
- # lossless compression, per CFITSIO)
- ncols = 4
- # CFITSIO 3.28 and up automatically use the GZIP_COMPRESSED_DATA
- # store floating point data that couldn't be quantized, instead
- # of the UNCOMPRESSED_DATA column. There's no way to control
- # this behavior so the only way to determine which behavior will
- # be employed is via the CFITSIO version
- if CFITSIO_SUPPORTS_GZIPDATA:
- ttype2 = 'GZIP_COMPRESSED_DATA'
- # The required format for the GZIP_COMPRESSED_DATA is actually
- # missing from the standard docs, but CFITSIO suggests it
- # should be 1PB, which is logical.
- tform2 = '1QB' if huge_hdu else '1PB'
- else:
- # Q format is not supported for UNCOMPRESSED_DATA columns.
- ttype2 = 'UNCOMPRESSED_DATA'
- if zbitpix == 8:
- tform2 = '1QB' if huge_hdu else '1PB'
- elif zbitpix == 16:
- tform2 = '1QI' if huge_hdu else '1PI'
- elif zbitpix == 32:
- tform2 = '1QJ' if huge_hdu else '1PJ'
- elif zbitpix == -32:
- tform2 = '1QE' if huge_hdu else '1PE'
- else:
- tform2 = '1QD' if huge_hdu else '1PD'
- # Set up the second column for the table that will hold any
- # uncompressable data.
- self._header.set('TTYPE2', ttype2, 'label for field 2',
- after='TFORM1')
- self._header.set('TFORM2', tform2,
- 'data format of field: variable length array',
- after='TTYPE2')
- col2 = Column(name=ttype2, format=tform2)
- # Set up the third column for the table that will hold
- # the scale values for quantized data.
- self._header.set('TTYPE3', 'ZSCALE', 'label for field 3',
- after='TFORM2')
- self._header.set('TFORM3', '1D',
- 'data format of field: 8-byte DOUBLE',
- after='TTYPE3')
- col3 = Column(name=self._header['TTYPE3'],
- format=self._header['TFORM3'])
- # Set up the fourth column for the table that will hold
- # the zero values for the quantized data.
- self._header.set('TTYPE4', 'ZZERO', 'label for field 4',
- after='TFORM3')
- self._header.set('TFORM4', '1D',
- 'data format of field: 8-byte DOUBLE',
- after='TTYPE4')
- after = 'TFORM4'
- col4 = Column(name=self._header['TTYPE4'],
- format=self._header['TFORM4'])
- # Create the ColDefs object for the table
- cols = ColDefs([col1, col2, col3, col4])
- else:
- # default table has just one 'COMPRESSED_DATA' column
- ncols = 1
- after = 'TFORM1'
- # remove any header cards for the additional columns that
- # may be left over from the previous data
- to_remove = ['TTYPE2', 'TFORM2', 'TTYPE3', 'TFORM3', 'TTYPE4',
- 'TFORM4']
- for k in to_remove:
- try:
- del self._header[k]
- except KeyError:
- pass
- # Create the ColDefs object for the table
- cols = ColDefs([col1])
- # Update the table header with the width of the table, the
- # number of fields in the table, the indicator for a compressed
- # image HDU, the data type of the image data and the number of
- # dimensions in the image data array.
- self._header.set('NAXIS1', cols.dtype.itemsize,
- 'width of table in bytes')
- self._header.set('TFIELDS', ncols, 'number of fields in each row')
- self._header.set('ZIMAGE', True, 'extension contains compressed image',
- after=after)
- self._header.set('ZBITPIX', zbitpix,
- bitpix_comment, after='ZIMAGE')
- self._header.set('ZNAXIS', self._image_header['NAXIS'], naxis_comment,
- after='ZBITPIX')
- # Strip the table header of all the ZNAZISn and ZTILEn keywords
- # that may be left over from the previous data
- idx = 1
- while True:
- try:
- del self._header['ZNAXIS' + str(idx)]
- del self._header['ZTILE' + str(idx)]
- idx += 1
- except KeyError:
- break
- # Verify that any input tile size parameter is the appropriate
- # size to match the HDU's data.
- naxis = self._image_header['NAXIS']
- if not tile_size:
- tile_size = []
- elif len(tile_size) != naxis:
- warnings.warn('Provided tile size not appropriate for the data. '
- 'Default tile size will be used.', AstropyUserWarning)
- tile_size = []
- # Set default tile dimensions for HCOMPRESS_1
- if compression_type == 'HCOMPRESS_1':
- if (self._image_header['NAXIS1'] < 4 or
- self._image_header['NAXIS2'] < 4):
- raise ValueError('Hcompress minimum image dimension is '
- '4 pixels')
- elif tile_size:
- if tile_size[0] < 4 or tile_size[1] < 4:
- # user specified tile size is too small
- raise ValueError('Hcompress minimum tile dimension is '
- '4 pixels')
- major_dims = len([ts for ts in tile_size if ts > 1])
- if major_dims > 2:
- raise ValueError(
- 'HCOMPRESS can only support 2-dimensional tile sizes.'
- 'All but two of the tile_size dimensions must be set '
- 'to 1.')
- if tile_size and (tile_size[0] == 0 and tile_size[1] == 0):
- # compress the whole image as a single tile
- tile_size[0] = self._image_header['NAXIS1']
- tile_size[1] = self._image_header['NAXIS2']
- for i in range(2, naxis):
- # set all higher tile dimensions = 1
- tile_size[i] = 1
- elif not tile_size:
- # The Hcompress algorithm is inherently 2D in nature, so the
- # row by row tiling that is used for other compression
- # algorithms is not appropriate. If the image has less than 30
- # rows, then the entire image will be compressed as a single
- # tile. Otherwise the tiles will consist of 16 rows of the
- # image. This keeps the tiles to a reasonable size, and it
- # also includes enough rows to allow good compression
- # efficiency. It the last tile of the image happens to contain
- # less than 4 rows, then find another tile size with between 14
- # and 30 rows (preferably even), so that the last tile has at
- # least 4 rows.
- # 1st tile dimension is the row length of the image
- tile_size.append(self._image_header['NAXIS1'])
- if self._image_header['NAXIS2'] <= 30:
- tile_size.append(self._image_header['NAXIS1'])
- else:
- # look for another good tile dimension
- naxis2 = self._image_header['NAXIS2']
- for dim in [16, 24, 20, 30, 28, 26, 22, 18, 14]:
- if naxis2 % dim == 0 or naxis2 % dim > 3:
- tile_size.append(dim)
- break
- else:
- tile_size.append(17)
- for i in range(2, naxis):
- # set all higher tile dimensions = 1
- tile_size.append(1)
- # check if requested tile size causes the last tile to have
- # less than 4 pixels
- remain = self._image_header['NAXIS1'] % tile_size[0] # 1st dimen
- if remain > 0 and remain < 4:
- tile_size[0] += 1 # try increasing tile size by 1
- remain = self._image_header['NAXIS1'] % tile_size[0]
- if remain > 0 and remain < 4:
- raise ValueError('Last tile along 1st dimension has '
- 'less than 4 pixels')
- remain = self._image_header['NAXIS2'] % tile_size[1] # 2nd dimen
- if remain > 0 and remain < 4:
- tile_size[1] += 1 # try increasing tile size by 1
- remain = self._image_header['NAXIS2'] % tile_size[1]
- if remain > 0 and remain < 4:
- raise ValueError('Last tile along 2nd dimension has '
- 'less than 4 pixels')
- # Set up locations for writing the next cards in the header.
- last_znaxis = 'ZNAXIS'
- if self._image_header['NAXIS'] > 0:
- after1 = 'ZNAXIS1'
- else:
- after1 = 'ZNAXIS'
- # Calculate the number of rows in the output table and
- # write the ZNAXISn and ZTILEn cards to the table header.
- nrows = 0
- for idx, axis in enumerate(self._axes):
- naxis = 'NAXIS' + str(idx + 1)
- znaxis = 'ZNAXIS' + str(idx + 1)
- ztile = 'ZTILE' + str(idx + 1)
- if tile_size and len(tile_size) >= idx + 1:
- ts = tile_size[idx]
- else:
- if not ztile in self._header:
- # Default tile size
- if not idx:
- ts = self._image_header['NAXIS1']
- else:
- ts = 1
- else:
- ts = self._header[ztile]
- tile_size.append(ts)
- if not nrows:
- nrows = (axis - 1) // ts + 1
- else:
- nrows *= ((axis - 1) // ts + 1)
- if image_header and naxis in image_header:
- self._header.set(znaxis, axis, image_header.comments[naxis],
- after=last_znaxis)
- else:
- self._header.set(znaxis, axis,
- 'length of original image axis',
- after=last_znaxis)
- self._header.set(ztile, ts, 'size of tiles to be compressed',
- after=after1)
- last_znaxis = znaxis
- after1 = ztile
- # Set the NAXIS2 header card in the table hdu to the number of
- # rows in the table.
- self._header.set('NAXIS2', nrows, 'number of rows in table')
- self.columns = cols
- # Set the compression parameters in the table header.
- # First, setup the values to be used for the compression parameters
- # in case none were passed in. This will be either the value
- # already in the table header for that parameter or the default
- # value.
- idx = 1
- while True:
- zname = 'ZNAME' + str(idx)
- if zname not in self._header:
- break
- zval = 'ZVAL' + str(idx)
- if