Drizzle

class drizzle.resample.Drizzle(kernel='square', fillval=None, fillval2=None, out_shape=None, out_img=None, out_wht=None, out_ctx=None, out_img2=None, out_dq=None, exptime=0.0, begin_ctx_id=0, max_ctx_id=None, disable_ctx=False)[source]

Bases: object

A class for managing resampling and co-adding of multiple images onto a common output grid. The main method of this class is add_image(). The main functionality of this class is to resample and co-add multiple images onto one output image using the “drizzle” algorithm described in Fruchter and Hook, PASP 2002. In the simplest terms, it redistributes flux from input pixels to one or more output pixels based on the chosen kernel, supplied weights, and input-to-output coordinate transformations as defined by the pixmap argument. For more details, see drizzle Documentation.

This class keeps track of the total exposure time of all co-added images and also of which input images have contributed to an output (resampled) pixel. This is accomplished via context image.

Main outputs of add_image() can be accessed as class properties out_img, out_img2, out_wht, out_ctx, and exptime.

Warning

Output arrays (out_img, out_img2, out_wht, and out_ctx) can be pre-allocated by the caller and be passed to the initializer or the class initializer can allocate these arrays based on other input parameters such as output_shape. If caller-supplied output arrays have the correct type (numpy.float32 for out_img, out_img2 and out_wht, numpy.int32 for the out_ctx array and numpy.uint32 for the out_dq array) and if out_ctx is large enough not to need to be resized, these arrays will be used as is and may be modified by the add_image() method. If not, a copy of these arrays will be made when converting to the expected type (or expanding the context array).

Scaling of input image data

It is important to highlight that the drizzle algorithm computes weighted mean of input pixel values – see equations (4) and (5) in Fruchter and Hook, PASP 2002. Therefore, it is important that all input pixel values that contribute to an output pixel are from the same distribution. In other words, input pixel values from different images must be on the same footing, i.e., they must be comparable and must be representative of the same physical quantity.

For example, for Hubble Space Telescope data, calibrated images (i.e., *_flt.fits, *_flc.fits) are in unit of counts, counts per second, electrons, or electrons per second. To convert them to flux units (e.g., erg/cm^2/s/Angstrom), one needs to multiply these images by the PHOTFLAM. Sometimes, images that are drizzle-combined have been observed at very different times (separated by many years) and the sensitivity of the instrument (represented by PHOTFLAM) may have changed significantly. Other times a source is observed in different chips, i.e., the two chips of the Wide Field Camera. In such cases detector’s sensitivity (PHOTFLAM) may be different for the images to be combined. Consequently, pixel values in these images may not be directly comparable and drizzle-combining such images would result in systematic errors.

In this case, it is important to rescale images to the same flux units either by multiplying by the appropriate PHOTFLAM values or some other appropriate scaling factor before combining them using drizzle. This can be accomplished by using the iscale parameter of add_image() which simply multiplies each input image by iscale.

Also, for the case of HST images that have flux units instead of surface brightness, if input images have different pixel scales, then the pixel values must be rescaled by the square of the pixel scale ratio (the linear dimension of a side of an output pixel as seen in the input image’s coordinate frame) in order to preserve flux. In this case iscale is equivalent to s**2 factor in equations (3) and (5) of Fruchter and Hook, PASP 2002 (s may be different for each input image).

Output Science Image

Output science image is obtained by computing weighted mean of input pixel values according to equations (4) and (5) in Fruchter and Hook, PASP 2002. The weights and coefficients in those equations will depend on the chosen kernel, input image weights, and pixel overlaps computed from pixmap.

Output Weight Image

Output weight image stores the total weight of output science pixels according to equation (4) in Fruchter and Hook, PASP 2002. It depends on the chosen kernel, input image weights, and pixel overlaps computed from pixmap.

Output Context Image

Each pixel in the context image is a bit field that encodes information about which input image has contributed to the corresponding pixel in the resampled data array. Context image uses 32 bit integers to encode this information and hence it can keep track of only 32 input images. The first bit corresponds to the first input image, the second bit corresponds to the second input image, and so on. We call this (0-indexed) order “context ID” which is represented by the ctx_id parameter/property. If the number of input images exceeds 32, then it is necessary to have multiple context images (“planes”) to hold information about all input images, with the first plane encoding which of the first 32 images contributed to the output data pixel, the second plane representing next 32 input images (number 33-64), etc. For this reason, context array is either a 2D array (if the total number of resampled images is less than 33) of the type numpy.int32 and shape (ny, nx) or a a 3D array of shape (np, ny, nx) where nx and ny are dimensions of the image data. np is the number of “planes” computed as (number of input images - 1) // 32 + 1. If a bit at position k in a pixel with coordinates (p, y, x) is 0, then input image number 32 * p + k (0-indexed) did not contribute to the output data pixel with array coordinates (y, x) and if that bit is 1, then input image number 32 * p + k did contribute to the pixel (y, x) in the resampled image.

As an example, let’s assume we have 8 input images. Then, when out_ctx pixel values are displayed using binary representation (and decimal in parenthesis), one could see values like this:

00000001 (1) - only first input image contributed to this output pixel;
00000010 (2) - 2nd input image contributed;
00000100 (4) - 3rd input image contributed;
10000000 (128) - 8th input image contributed;
10000100 (132=128+4) - 3rd and 8th input images contributed;
11001101 (205=1+4+8+64+128) - input images 1, 3, 4, 7, 8 have contributed
to this output pixel.

In order to test if a specific input image contributed to an output pixel, one needs to use bitwise operations. Using the example above, to test whether input images number 4 and 5 have contributed to the output pixel whose corresponding out_ctx value is 205 (11001101 in binary form) we can do the following:

>>> bool(205 & (1 << (5 - 1)))  # (205 & 16) = 0 (== 0 => False): did NOT contribute
False
>>> bool(205 & (1 << (4 - 1)))  # (205 & 8) = 8 (!= 0 => True): did contribute
True

In general, to get a list of all input images that have contributed to an output resampled pixel with image coordinates (x, y), and given a context array ctx, one can do something like this:

>>> import numpy as np
>>> np.flatnonzero([v & (1 << k) for v in ctx[:, y, x] for k in range(32)])

For convenience, this functionality was implemented in the decode_context() function.

Output DQ Image

If DQ array of input image pixels is provided via dq parameter of add_image(), then an output DQ array will be computed by combining (using bitwise-OR) DQ bitfields of all input pixels that contribute to a given output pixel.

References

A full description of the drizzling algorithm can be found in Fruchter and Hook, PASP 2002.

Examples

# wcs1 - WCS of the input image usually with distortions (to be resampled)
# wcs2 - WCS of the output image without distortions

import numpy as np
from drizzle.resample import Drizzle
from drizzle.utils import calc_pixmap

# simulate some data and a pixel map:
data = np.ones((240, 570))
pixmap = calc_pixmap(wcs1, wcs2)
# or simulate a mapping from input image to output image frame:
# y, x = np.indices((240, 570), dtype=np.float64)
# pixmap = np.dstack([x, y])

# initialize Drizzle object
d = Drizzle(out_shape=(240, 570))
d.add_image(data, exptime=15, pixmap=pixmap)

# access outputs:
d.out_img
d.out_ctx
d.out_wht

Attributes Summary

ctx_id

Context image "ID" (0-based ) of the next image to be resampled.

fillval

Fill value for output pixels without contributions from input images.

fillval2

Fill value for output pixels in out_img2 without contributions from input images.

kernel

Resampling kernel.

out_ctx

Output "context" image.

out_dq

Output DQ image computed by OR-combining DQ bitfields of input images' pixels that have contributed to a given output pixel.

out_img

Output resampled image.

out_img2

Output resampled image(s) obtained with squared weights.

out_wht

Output weight image.

total_exptime

Total exposure time of all resampled images.

Methods Summary

add_image(data, exptime, pixmap[, data2, ...])

Resample and add an image to the cumulative output image.

Attributes Documentation

ctx_id

Context image “ID” (0-based ) of the next image to be resampled.

fillval

Fill value for output pixels without contributions from input images.

fillval2

Fill value for output pixels in out_img2 without contributions from input images.

kernel

Resampling kernel.

out_ctx

Output “context” image.

out_dq

Output DQ image computed by OR-combining DQ bitfields of input images’ pixels that have contributed to a given output pixel.

out_img

Output resampled image.

out_img2

Output resampled image(s) obtained with squared weights. It is always a list of one or more 2D arrays.

out_wht

Output weight image.

total_exptime

Total exposure time of all resampled images.

Methods Documentation

add_image(data, exptime, pixmap, data2=None, dq=None, scale=<object object>, iscale=1.0, pixel_scale_ratio=1.0, weight_map=None, wht_scale=1.0, pixfrac=1.0, in_units='cps', xmin=None, xmax=None, ymin=None, ymax=None)[source]

Resample and add an image to the cumulative output image. Also, update output total weight image and context images.

Parameters:
  • data (2D numpy.ndarray) – A 2D numpy array containing the input image to be drizzled.

  • exptime (float) – The exposure time of the input image, a positive number. The exposure time is used to scale the image if the units are counts.

  • pixmap (3D array) – A mapping from input image (data) coordinates to resampled (out_img) coordinates. pixmap must be an array of shape (Ny, Nx, 2) where (Ny, Nx) is the shape of the input image. pixmap[..., 0] forms a 2D array of X-coordinates of input pixels in the output frame and pixmap[..., 1] forms a 2D array of Y-coordinates of input pixels in the output coordinate frame.

  • data2 (2D array of float32, list of 2D arrays of float32 or None, None, optional) –

    A 2D numpy array (or a list of 2D arrays) with image data to be resampled and co-added using squared weights. The resampled output image can be accessed via out_img2 property of the Drizzle object. This is useful for performing standard error propagation using variance arrays.

    Multiple data arrays (of the same shape as that of data) can be provided as a list of 2D arrays. The number of arrays must match the number of output data arrays provided during initialization via argument out_img2. If an item in the list is None, that item will not be resampled to the corresponding out_img2 element.

    Note

    It is assumed that data in data2 have squared units of data. Therefore, when in_units are “counts”, data2 arrays will be rescaled by exptime**2 to convert to rate units before resampling.

  • dq (2D array, None, optional) –

    A 2D numpy array of type numpy.uint32 (unsigned 32-bit integer type) containing DQ bitfields of input pixels. It must have the same shape as data. If provided, output DQ array (accessible via out_dq property) will be computed by combining (using bitwise-OR) DQ bitfields of input pixels that contributed to the output pixel. If None, DQ array of the output image will not be computed.

    Warning

    64-bit integer type is not supported and will raise an exception. Contact the authors to add support for 64-bit DQ if you need it.

  • scale (float, optional) – Deprecated: use iscale and pixel_scale_ratio instead. It is a factor used both to rescale input image data by scale**2 AND to compute the correct kernel size for some kernels (“turbo”, “gaussian”, and “lanczos”). It is recommended scale be set to pixel scale ratio: the linear dimension of a side of an output pixel relative to the size of an input pixel (or size of an output pixel in the input image’s coordinate system).

  • iscale (float, optional) – It is a multiplicative factor used to rescale input image data by iscale value. data2 images will be rescaled by iscale**2. It may make sense to rescale input image (data) by squared pixel scale ratio (the linear dimension of a side of an output pixel as seen in the input image’s coordinate frame) depending on the units of the input image, i.e., counts vs brightness. For more details see section “Scaling of input image data” in Drizzle.

  • pixel_scale_ratio (float, None, optional) – It is a factor used to compute the correct kernel size in output image’s coordinate system for some of the kernels (“turbo”, “gaussian”, and “lanczos”) from their nominal sizes in input image pixels. For example, for the “lanczos3” kernel, the nominal size is 3 input pixels. It is recommended that pixel_scale_ratio be set to pixel scale ratio: the linear dimension of output pixel relative to the size of an input pixel. When pixel_scale_ratio is None, it will be estimated from pixmap but this can impose a performance penalty.

  • weight_map (2D array, None, optional) –

    A 2D numpy array containing the pixel by pixel weighting. Must have the same dimensions as data.

    When weight_map is None, the weight of input data pixels will be assumed to be 1.

  • wht_scale (float) – A scaling factor applied to the pixel by pixel weighting.

  • pixfrac (float, optional) – The fraction of a pixel that the pixel flux is confined to. The default value of 1 has the pixel flux evenly spread across the image. A value of 0.5 confines it to half a pixel in the linear dimension, so the flux is confined to a quarter of the pixel area when the square kernel is used.

  • in_units (str) – The units of the input image. The units can either be “counts” or “cps” (counts per second.)

  • xmin (float, optional) – This and the following three parameters set a bounding rectangle on the input image. Only pixels on the input image inside this rectangle will have their flux added to the output image. Xmin sets the minimum value of the x dimension. The x dimension is the dimension that varies quickest on the image. If the value is zero, no minimum will be set in the x dimension. All four parameters are zero based, counting starts at zero.

  • xmax (float, optional) – Sets the maximum value of the x dimension on the bounding box of the input image. If the value is zero, no maximum will be set in the x dimension, the full x dimension of the output image is the bounding box.

  • ymin (float, optional) – Sets the minimum value in the y dimension on the bounding box. The y dimension varies less rapidly than the x and represents the line index on the input image. If the value is zero, no minimum will be set in the y dimension.

  • ymax (float, optional) – Sets the maximum value in the y dimension. If the value is zero, no maximum will be set in the y dimension, the full x dimension of the output image is the bounding box.

Returns:

  • nskip (float) – The number of lines from the box defined by ((xmin, xmax), (ymin, ymax)) in the input image that were ignored and did not contribute to the output image.

  • nmiss (float) – The number of pixels from the box defined by ((xmin, xmax), (ymin, ymax)) in the input image that were ignored and did not contribute to the output image.