diff --git a/drizzle/__init__.py b/drizzle/__init__.py index 9418e47b..66fcd87b 100644 --- a/drizzle/__init__.py +++ b/drizzle/__init__.py @@ -8,4 +8,4 @@ __version__ = version(__name__) except PackageNotFoundError: # package is not installed - __version__ = 'unknown' + __version__ = '' diff --git a/drizzle/doblot.py b/drizzle/doblot.py deleted file mode 100644 index 818ebedd..00000000 --- a/drizzle/doblot.py +++ /dev/null @@ -1,80 +0,0 @@ -""" -STScI Python compatable blot module -""" -import numpy as np - -from drizzle import calc_pixmap -from drizzle import cdrizzle - - -def doblot(source, source_wcs, blot_wcs, exptime, coeffs=True, - interp='poly5', sinscl=1.0, stepsize=10, wcsmap=None): - """ - Low level routine for performing the 'blot' operation. - - Create a single blotted image from a single source image. The - interface is compatible with STScI code. All distortion information - is assumed to be included in the WCS specification of the 'output' - blotted image given in 'blot_wcs'. - - Parameters - ---------- - - source : 2d array - Input numpy array of the source image in units of 'cps'. - - source_wcs : wcs - The source image WCS. - - blot_wcs : wcs - The blotted image WCS. The WCS that the source image will be - resampled to. - - exptime : float - The exposure time of the input image. - - interp : str, optional - The type of interpolation used in the blotting. The - possible values are "nearest" (nearest neighbor interpolation), - "linear" (bilinear interpolation), "poly3" (cubic polynomial - interpolation), "poly5" (quintic polynomial interpolation), - "sinc" (sinc interpolation), "lan3" (3rd order Lanczos - interpolation), and "lan5" (5th order Lanczos interpolation). - - sincscl : float, optional - The scaling factor for sinc interpolation. - - Returns - ------- - - A 2d numpy array with the blotted image - - Other Parameters - ---------------- - - coeffs : bool, optional - Not used. Only kept for backwards compatibility. - - stepsize : float, optional - Was used when input to output mapping was computed - internally. Is no longer used and only here for backwards compatibility. - - wcsmap : function, optional - Was used when input to output mapping was computed - internally. Is no longer used and only here for backwards compatibility. - """ - _outsci = np.zeros(blot_wcs.pixel_shape[::-1], dtype=np.float32) - - # compute the undistorted 'natural' plate scale - blot_wcs.sip = None - blot_wcs.cpdis1 = None - blot_wcs.cpdis2 = None - blot_wcs.det2im = None - - pixmap = calc_pixmap.calc_pixmap(blot_wcs, source_wcs) - pix_ratio = source_wcs.pscale / blot_wcs.pscale - - cdrizzle.tblot(source, pixmap, _outsci, scale=pix_ratio, kscale=1.0, - interp=interp, exptime=exptime, misval=0.0, sinscl=sinscl) - - return _outsci diff --git a/drizzle/dodrizzle.py b/drizzle/dodrizzle.py deleted file mode 100644 index ced4d044..00000000 --- a/drizzle/dodrizzle.py +++ /dev/null @@ -1,186 +0,0 @@ -""" -STScI Python compatable drizzle module -""" -import numpy as np - -from . import util -from . import calc_pixmap -from . import cdrizzle - - -def dodrizzle(insci, input_wcs, inwht, - output_wcs, outsci, outwht, outcon, - expin, in_units, wt_scl, - wcslin_pscale=1.0, uniqid=1, - xmin=0, xmax=0, ymin=0, ymax=0, - pixfrac=1.0, kernel='square', fillval="INDEF"): - """ - Low level routine for performing 'drizzle' operation.on one image. - - The interface is compatible with STScI code. All images are Python - ndarrays, instead of filenames. File handling (input and output) is - performed by the calling routine. - - Parameters - ---------- - - insci : 2d array - A 2d numpy array containing the input image to be drizzled. - it is an error to not supply an image. - - input_wcs : 2d array - The world coordinate system of the input image. - - inwht : 2d array - A 2d numpy array containing the pixel by pixel weighting. - Must have the same dimensions as insci. If none is supplied, - the weghting is set to one. - - output_wcs : wcs - The world coordinate system of the output image. - - outsci : 2d array - A 2d numpy array containing the output image produced by - drizzling. On the first call it should be set to zero. - Subsequent calls it will hold the intermediate results - - outwht : 2d array - A 2d numpy array containing the output counts. On the first - call it should be set to zero. On subsequent calls it will - hold the intermediate results. - - outcon : 2d or 3d array, optional - A 2d or 3d numpy array holding a bitmap of which image was an input - for each output pixel. Should be integer zero on first call. - Subsequent calls hold intermediate results. - - expin : 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. - - in_units : str - The units of the input image. The units can either be "counts" - or "cps" (counts per second.) - - wt_scl : float - A scaling factor applied to the pixel by pixel weighting. - - wcslin_pscale : float, optional - The pixel scale of the input image. Conceptually, this is the - linear dimension of a side of a pixel in the input image, but it - is not limited to this and can be set to change how the drizzling - algorithm operates. - - uniqid : int, optional - The id number of the input image. Should be one the first time - this function is called and incremented by one on each subsequent - call. - - 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. - - 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. - - kernel: str, optional - The name of the kernel used to combine the input. The choice of - kernel controls the distribution of flux over the kernel. The kernel - names are: "square", "gaussian", "point", "tophat", "turbo", "lanczos2", - and "lanczos3". The square kernel is the default. - - fillval: str, optional - The value a pixel is set to in the output if the input image does - not overlap it. The default value of INDEF does not set a value. - - Returns - ------- - A tuple with three values: a version string, the number of pixels - on the input image that do not overlap the output image, and the - number of complete lines on the input image that do not overlap the - output input image. - - """ - - # Ensure that the fillval parameter gets properly interpreted - # for use with tdriz - if util.is_blank(fillval): - fillval = 'INDEF' - else: - fillval = str(fillval) - - if in_units == 'cps': - expscale = 1.0 - else: - expscale = expin - - # Add input weight image if it was not passed in - - if (insci.dtype > np.float32): - insci = insci.astype(np.float32) - - if inwht is None: - inwht = np.ones_like(insci) - - # Compute what plane of the context image this input would - # correspond to: - planeid = int((uniqid - 1) / 32) - - # Check if the context image has this many planes - if outcon.ndim == 3: - nplanes = outcon.shape[0] - elif outcon.ndim == 2: - nplanes = 1 - else: - nplanes = 0 - - if nplanes <= planeid: - raise IndexError("Not enough planes in drizzle context image") - - # Alias context image to the requested plane if 3d - if outcon.ndim == 3: - outcon = outcon[planeid] - - pix_ratio = output_wcs.pscale / wcslin_pscale - - # Compute the mapping between the input and output pixel coordinates - pixmap = calc_pixmap.calc_pixmap(input_wcs, output_wcs) - - # - # Call 'drizzle' to perform image combination - # This call to 'cdriz.tdriz' uses the new C syntax - # - _vers, nmiss, nskip = cdrizzle.tdriz( - insci, inwht, pixmap, outsci, outwht, outcon, - uniqid=uniqid, xmin=xmin, xmax=xmax, - ymin=ymin, ymax=ymax, scale=pix_ratio, pixfrac=pixfrac, - kernel=kernel, in_units=in_units, expscale=expscale, - wtscale=wt_scl, fillstr=fillval) - - return _vers, nmiss, nskip diff --git a/drizzle/drizzle.py b/drizzle/drizzle.py deleted file mode 100644 index 2cbfb9e4..00000000 --- a/drizzle/drizzle.py +++ /dev/null @@ -1,566 +0,0 @@ -""" -The `drizzle` module defines the `Drizzle` class, for combining input -images into a single output image using the drizzle algorithm. -""" -import os -import os.path - -import numpy as np -from astropy import wcs -from astropy.io import fits - -from . import util -from . import doblot -from . import dodrizzle - - -class Drizzle(object): - """ - Combine images using the drizzle algorithm - """ - def __init__(self, infile="", outwcs=None, - wt_scl="exptime", pixfrac=1.0, kernel="square", - fillval="INDEF"): - """ - Create a new Drizzle output object and set the drizzle parameters. - - All parameters are optional, but either infile or outwcs must be supplied. - If infile initializes the object from a file written after a - previous run of drizzle. Results from the previous run will be combined - with new results. The value passed in outwcs will be ignored. If infile is - not set, outwcs will be used to initilize a new run of drizzle. - - Parameters - ---------- - - infile : str, optional - A fits file containing results from a previous run. The three - extensions SCI, WHT, and CTX contain the combined image, total counts - and image id bitmap, repectively. The WCS of the combined image is - also read from the SCI extension. - - outwcs : wcs, optional - The world coordinate system (WCS) of the combined image. This - parameter must be present if no input file is given and is ignored if - one is. - - wt_scl : str, optional - How each input image should be scaled. The choices are `exptime` - which scales each image by its exposure time, `expsq` which scales - each image by the exposure time squared, or an empty string, which - allows each input image to be scaled individually. - - 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. - - kernel : str, optional - The name of the kernel used to combine the inputs. The choice of - kernel controls the distribution of flux over the kernel. The kernel - names are: "square", "gaussian", "point", "tophat", "turbo", "lanczos2", - and "lanczos3". The square kernel is the default. - - fillval : str, otional - The value a pixel is set to in the output if the input image does - not overlap it. The default value of INDEF does not set a value. - """ - - # Initialize the object fields - - self.outsci = None - self.outwht = None - self.outcon = None - - self.outexptime = 0.0 - self.uniqid = 0 - - self.outwcs = outwcs - self.wt_scl = wt_scl - self.kernel = kernel - self.fillval = fillval - self.pixfrac = float(pixfrac) - - self.sciext = "SCI" - self.whtext = "WHT" - self.ctxext = "CTX" - - out_units = "cps" - - if not util.is_blank(infile): - if os.path.exists(infile): - handle = fits.open(infile) - - # Read parameters from image header - self.outexptime = util.get_keyword(handle, "DRIZEXPT", default=0.0) - self.uniqid = util.get_keyword(handle, "NDRIZIM", default=0) - - self.sciext = util.get_keyword(handle, "DRIZOUDA", default="SCI") - self.whtext = util.get_keyword(handle, "DRIZOUWE", default="WHT") - self.ctxext = util.get_keyword(handle, "DRIZOUCO", default="CTX") - - self.wt_scl = util.get_keyword(handle, "DRIZWTSC", default=wt_scl) - self.kernel = util.get_keyword(handle, "DRIZKERN", default=kernel) - self.fillval = util.get_keyword(handle, "DRIZFVAL", default=fillval) - self.pixfrac = float(util.get_keyword(handle, - "DRIZPIXF", default=pixfrac)) - - out_units = util.get_keyword(handle, "DRIZOUUN", default="cps") - - try: - hdu = handle[self.sciext] - self.outsci = hdu.data.copy().astype(np.float32) - self.outwcs = wcs.WCS(hdu.header, fobj=handle) - except KeyError: - pass - - try: - hdu = handle[self.whtext] - self.outwht = hdu.data.copy().astype(np.float32) - except KeyError: - pass - - try: - hdu = handle[self.ctxext] - self.outcon = hdu.data.copy().astype(np.int32) - if self.outcon.ndim == 2: - self.outcon = np.reshape(self.outcon, (1, - self.outcon.shape[0], - self.outcon.shape[1])) - - elif self.outcon.ndim == 3: - pass - - else: - msg = ("Drizzle context image has wrong dimensions: " + - infile) - raise ValueError(msg) - - except KeyError: - pass - - handle.close() - - # Check field values - - if self.outwcs: - util.set_pscale(self.outwcs) - else: - raise ValueError("Either an existing file or wcs must be supplied to Drizzle") - - if util.is_blank(self.wt_scl): - self.wt_scl = '' - elif self.wt_scl != "exptime" and self.wt_scl != "expsq": - raise ValueError("Illegal value for wt_scl: %s" % out_units) - - if out_units == "counts": - np.divide(self.outsci, self.outexptime, self.outsci) - elif out_units != "cps": - raise ValueError("Illegal value for wt_scl: %s" % out_units) - - # Initialize images if not read from a file - outwcs_naxis1, outwcs_naxis2 = self.outwcs.pixel_shape - if self.outsci is None: - self.outsci = np.zeros(self.outwcs.pixel_shape[::-1], - dtype=np.float32) - - if self.outwht is None: - self.outwht = np.zeros(self.outwcs.pixel_shape[::-1], - dtype=np.float32) - if self.outcon is None: - self.outcon = np.zeros((1, outwcs_naxis2, outwcs_naxis1), - dtype=np.int32) - - def add_fits_file(self, infile, inweight="", - xmin=0, xmax=0, ymin=0, ymax=0, - unitkey="", expkey="", wt_scl=1.0): - """ - Combine a fits file with the output drizzled image. - - Parameters - ---------- - - infile : str - The name of the fits file, possibly including an extension. - - inweight : str, otional - The name of a file containing a pixel by pixel weighting - of the input data. If it is not set, an array will be generated - where all values are set to one. - - xmin : float, otional - This and the following three parameters set a bounding rectangle - on the output image. Only pixels on the output image inside this - rectangle will have their flux updated. 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 or less, no minimum will - be set in the x dimension. All four parameters are zero based, - counting starts at zero. - - xmax : float, otional - Sets the maximum value of the x dimension on the bounding box - of the ouput image. If the value is zero or less, no maximum will - be set in the x dimension. - - 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 output image. If the value is zero or less, 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 or - less, no maximum will be set in the y dimension. - - unitkey : string, optional - The name of the header keyword containing the image units. The - units can either be "counts" or "cps" (counts per second.) If it is - left blank, the value is assumed to be "cps." If the value is counts, - before using the input image it is scaled by dividing it by the - exposure time. - - expkey : string, optional - The name of the header keyword containing the exposure time. The - exposure time is used to scale the image if the units are counts and - to scale the image weighting if the drizzle was initialized with - wt_scl equal to "exptime" or "expsq." If the value of this parameter - is blank, the exposure time is set to one, implying no scaling. - - wt_scl : float, optional - If drizzle was initialized with wt_scl left blank, this value will - set a scaling factor for the pixel weighting. If drizzle was - initialized with wt_scl set to "exptime" or "expsq", the exposure time - will be used to set the weight scaling and the value of this parameter - will be ignored. - """ - - insci = None - inwht = None - - if not util.is_blank(infile): - fileroot, extn = util.parse_filename(infile) - - if os.path.exists(fileroot): - handle = fits.open(fileroot) - hdu = util.get_extn(handle, extn=extn) - - if hdu is not None: - insci = hdu.data - inwcs = wcs.WCS(header=hdu.header) - insci = hdu.data.copy() - handle.close() - - if insci is None: - raise ValueError("Drizzle cannot find input file: %s" % infile) - - if not util.is_blank(inweight): - fileroot, extn = util.parse_filename(inweight) - - if os.path.exists(fileroot): - handle = fits.open(fileroot) - hdu = util.get_extn(handle, extn=extn) - - if hdu is not None: - inwht = hdu.data.copy() - handle.close() - - in_units = util.get_keyword(fileroot, unitkey, "cps") - expin = util.get_keyword(fileroot, expkey, 1.0) - - self.add_image(insci, inwcs, inwht=inwht, - xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, - expin=expin, in_units=in_units, wt_scl=wt_scl) - - def add_image(self, insci, inwcs, inwht=None, - xmin=0, xmax=0, ymin=0, ymax=0, - expin=1.0, in_units="cps", wt_scl=1.0): - """ - Combine an input image with the output drizzled image. - - Instead of reading the parameters from a fits file, you can set - them by calling this lower level method. `Add_fits_file` calls - this method after doing its setup. - - Parameters - ---------- - - insci : array - A 2d numpy array containing the input image to be drizzled. - it is an error to not supply an image. - - inwcs : wcs - The world coordinate system of the input image. This is - used to convert the pixels to the output coordinate system. - - inwht : array, optional - A 2d numpy array containing the pixel by pixel weighting. - Must have the same dimenstions as insci. If none is supplied, - the weghting is set to one. - - xmin : float, optional - This and the following three parameters set a bounding rectangle - on the output image. Only pixels on the output image inside this - rectangle will have their flux updated. 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 or less, 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 ouput image. If the value is zero or less, no maximum will - be set in the x dimension. - - 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 output image. If the value is zero or less, 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 or - less, no maximum will be set in the y dimension. - - expin : float, optional - The exposure time of the input image, a positive number. The - exposure time is used to scale the image if the units are counts and - to scale the image weighting if the drizzle was initialized with - wt_scl equal to "exptime" or "expsq." - - in_units : str, optional - The units of the input image. The units can either be "counts" - or "cps" (counts per second.) If the value is counts, before using - the input image it is scaled by dividing it by the exposure time. - - wt_scl : float, optional - If drizzle was initialized with wt_scl left blank, this value will - set a scaling factor for the pixel weighting. If drizzle was - initialized with wt_scl set to "exptime" or "expsq", the exposure time - will be used to set the weight scaling and the value of this parameter - will be ignored. - """ - - insci = insci.astype(np.float32) - util.set_pscale(inwcs) - - if inwht is None: - inwht = np.ones(insci.shape, dtype=insci.dtype) - else: - inwht = inwht.astype(np.float32) - - if self.wt_scl == "exptime": - wt_scl = expin - elif self.wt_scl == "expsq": - wt_scl = expin * expin - - self.increment_id() - self.outexptime += expin - - dodrizzle.dodrizzle(insci, inwcs, inwht, self.outwcs, - self.outsci, self.outwht, self.outcon, - expin, in_units, wt_scl, - wcslin_pscale=inwcs.pscale, uniqid=self.uniqid, - xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, - pixfrac=self.pixfrac, kernel=self.kernel, - fillval=self.fillval) - - def blot_fits_file(self, infile, interp='poly5', sinscl=1.0): - """ - Resample the output using another image's world coordinate system. - - Parameters - ---------- - - infile : str - The name of the fits file containing the world coordinate - system that the output file will be resampled to. The name may - possibly include an extension. - - interp : str, optional - The type of interpolation used in the resampling. The - possible values are "nearest" (nearest neighbor interpolation), - "linear" (bilinear interpolation), "poly3" (cubic polynomial - interpolation), "poly5" (quintic polynomial interpolation), - "sinc" (sinc interpolation), "lan3" (3rd order Lanczos - interpolation), and "lan5" (5th order Lanczos interpolation). - - sincscl : float, optional - The scaling factor for sinc interpolation. - """ - blotwcs = None - - fileroot, extn = util.parse_filename(infile) - - if os.path.exists(fileroot): - handle = fits.open(fileroot) - hdu = util.get_extn(handle, extn=extn) - - if hdu is not None: - blotwcs = wcs.WCS(header=hdu.header) - handle.close() - - if not blotwcs: - raise ValueError("Drizzle did not get a blot reference image") - - self.blot_image(blotwcs, interp=interp, sinscl=sinscl) - - def blot_image(self, blotwcs, interp='poly5', sinscl=1.0): - """ - Resample the output image using an input world coordinate system. - - Parameters - ---------- - - blotwcs : wcs - The world coordinate system to resample on. - - interp : str, optional - The type of interpolation used in the resampling. The - possible values are "nearest" (nearest neighbor interpolation), - "linear" (bilinear interpolation), "poly3" (cubic polynomial - interpolation), "poly5" (quintic polynomial interpolation), - "sinc" (sinc interpolation), "lan3" (3rd order Lanczos - interpolation), and "lan5" (5th order Lanczos interpolation). - - sincscl : float, optional - The scaling factor for sinc interpolation. - """ - - util.set_pscale(blotwcs) - self.outsci = doblot.doblot(self.outsci, self.outwcs, blotwcs, - 1.0, interp=interp, sinscl=sinscl) - - self.outwcs = blotwcs - - def increment_id(self): - """ - Increment the id count and add a plane to the context image if needed - - Drizzle tracks which input images contribute to the output image - by setting a bit in the corresponding pixel in the context image. - The uniqid indicates which bit. So it must be incremented each time - a new image is added. Each plane in the context image can hold 32 bits, - so after each 32 images, a new plane is added to the context. - """ - - # Compute what plane of the context image this input would - # correspond to: - planeid = int(self.uniqid / 32) - - # Add a new plane to the context image if planeid overflows - - if self.outcon.shape[0] == planeid: - plane = np.zeros_like(self.outcon[0]) - self.outcon = np.append(self.outcon, [plane], axis=0) - - # Increment the id - self.uniqid += 1 - - def write(self, outfile, out_units="cps", outheader=None): - """ - Write the output from a set of drizzled images to a file. - - The output file will contain three extensions. The "SCI" extension - contains the resulting image. The "WHT" extension contains the - combined weights. The "CTX" extension is a bit map. The nth bit - is set to one if the nth input image contributed non-zero flux - to the output image. The "CTX" image is three dimensionsional - to account for the possibility that there are more than 32 input - images. - - Parameters - ---------- - - outfile : str - The name of the output file. If the file already exists, - the old file is deleted after writing the new file. - - out_units : str, optional - The units of the output image, either `counts` or `cps` - (counts per second.) If the units are counts, the resulting - image will be multiplied by the computed exposure time. - - outheader : header, optional - A fits header containing cards to be added to the primary - header of the output image. - """ - - if out_units != "counts" and out_units != "cps": - raise ValueError("Illegal value for out_units: %s" % str(out_units)) - - # Write the WCS to the output image - - handle = self.outwcs.to_fits() - phdu = handle[0] - - # Write the class fields to the primary header - phdu.header['DRIZOUDA'] = \ - (self.sciext, 'Drizzle, output data image') - phdu.header['DRIZOUWE'] = \ - (self.whtext, 'Drizzle, output weighting image') - phdu.header['DRIZOUCO'] = \ - (self.ctxext, 'Drizzle, output context image') - phdu.header['DRIZWTSC'] = \ - (self.wt_scl, 'Drizzle, weighting factor for input image') - phdu.header['DRIZKERN'] = \ - (self.kernel, 'Drizzle, form of weight distribution kernel') - phdu.header['DRIZPIXF'] = \ - (self.pixfrac, 'Drizzle, linear size of drop') - phdu.header['DRIZFVAL'] = \ - (self.fillval, 'Drizzle, fill value for zero weight output pix') - phdu.header['DRIZOUUN'] = \ - (out_units, 'Drizzle, units of output image - counts or cps') - - # Update header keyword NDRIZIM to keep track of how many images have - # been combined in this product so far - phdu.header['NDRIZIM'] = (self.uniqid, 'Drizzle, number of images') - - # Update header of output image with exptime used to scale the output data - # if out_units is not counts, this will simply be a value of 1.0 - # the keyword 'exptime' will always contain the total exposure time - # of all input image regardless of the output units - - phdu.header['EXPTIME'] = \ - (self.outexptime, 'Drizzle, total exposure time') - - outexptime = 1.0 - if out_units == 'counts': - np.multiply(self.outsci, self.outexptime, self.outsci) - outexptime = self.outexptime - phdu.header['DRIZEXPT'] = (outexptime, 'Drizzle, exposure time scaling factor') - - # Copy the optional header to the primary header - - if outheader: - phdu.header.extend(outheader, unique=True) - - # Add three extensions containing, the drizzled output image, - # the total counts, and the context bitmap, in that order - - extheader = self.outwcs.to_header() - - ehdu = fits.ImageHDU() - ehdu.data = self.outsci - ehdu.header['EXTNAME'] = (self.sciext, 'Extension name') - ehdu.header['EXTVER'] = (1, 'Extension version') - ehdu.header.extend(extheader, unique=True) - handle.append(ehdu) - - whdu = fits.ImageHDU() - whdu.data = self.outwht - whdu.header['EXTNAME'] = (self.whtext, 'Extension name') - whdu.header['EXTVER'] = (1, 'Extension version') - whdu.header.extend(extheader, unique=True) - handle.append(whdu) - - xhdu = fits.ImageHDU() - xhdu.data = self.outcon - xhdu.header['EXTNAME'] = (self.ctxext, 'Extension name') - xhdu.header['EXTVER'] = (1, 'Extension version') - xhdu.header.extend(extheader, unique=True) - handle.append(xhdu) - - handle.writeto(outfile, overwrite=True) - handle.close() diff --git a/drizzle/resample.py b/drizzle/resample.py new file mode 100644 index 00000000..db276b3c --- /dev/null +++ b/drizzle/resample.py @@ -0,0 +1,384 @@ +""" +The `drizzle` module defines the `Drizzle` class, for combining input +images into a single output image using the drizzle algorithm. +""" +import numpy as np + +from . import cdrizzle + + +SUPPORTED_DRIZZLE_KERNELS = [ + "square", + "gaussian", + "point", + "tophat", + "turbo", + "lanczos2", + "lanczos3" +] + +CTX_PLANE_BITS = 32 + + +class Drizzle(): + def __init__(self, n_images=1, kernel="square", + fillval=None, out_shape=None, out_sci=None, out_wht=None, + out_ctx=None): + """ + n_images : int, optional + The number of images expected to be added to the resampled + output. When it is a positive number and ``out_ctx`` is `None`, + it allows to pre-allocate the necessary array for the output context + image. If the actual number of input images that will be resampled + will exceed initial allocation for the context image, additional + context planes will be added as needed (context array will "grow" + in the third dimention as new input images are added.) + This parameter is ignored when ``out_ctx`` is provided. + + kernel: str, optional + The name of the kernel used to combine the input. The choice of + kernel controls the distribution of flux over the kernel. The kernel + names are: "square", "gaussian", "point", "tophat", "turbo", + "lanczos2", and "lanczos3". The square kernel is the default. + + out_shape : tuple, None, optional + Shape of the output images (context image will have a third + dimension of size proportional to the number of input images). + This parameter is helpful when neither ``out_sci``, ``out_wht``, + nor ``out_ctx`` images are provided. + + fillval: float, None, str, optional + The value of output pixels that did not have contributions from + input images' pixels. When ``fillval`` is either `None` or + ``"INDEF"`` and ``out_sci`` is provided, the values of ``out_sci`` + will not be modified. When ``fillval`` is either `None` or + ``"INDEF"`` and ``out_sci`` is **not provided**, the values of + ``out_sci`` will be initialized to 0. If ``fillval`` is a number + of a string that can be converted to a number, then the output + pixels with no contributions from input images will be set to this + ``fillval`` value. + + out_sci : 2d array of float32, optional + A 2D numpy array containing the output image produced by + drizzling. On the first call it should be set to zero. + Subsequent calls it will hold the intermediate results + + out_wht : 2D array of float32, optional + A 2D numpy array containing the output counts. On the first + call it should be set to zero. On subsequent calls it will + hold the intermediate results. + + out_ctx : 2D or 3D array of int32, optional + A 2D or 3D numpy array holding a bitmap of which image was an input + for each output pixel. Should be integer zero on first call. + Subsequent calls hold intermediate results. + + """ + if not isinstance(ctx_id, numbers.Integral) or ctx_id < 0: + raise TypeError("'ctx_id' must be a non-negative integer number") + self.ctx_id = -1 # the ID of the *last* image to be resampled + self.exptime = 0.0 + + if kernel.lower() not in SUPPORTED_DRIZZLE_KERNELS: + raise ValueError(f"Kernel '{kernel}' is not supported.") + self.kernel = kernel + + if fillval is None: + fillval = "INDEF" + elif isinstance(str, fillval): + fillval = fillval.strip() + if fillval == "": + fillval = "INDEF" + elif fillval.upper() == "INDEF": + pass + else: + float(fillval) + fillval = str(fillval) + self.fillval = fillval + + if (out_shape is None and out_sci is None and out_wht is None and + out_ctx is None): + raise ValueError( + "'out_shape' cannot be None when all output arrays are None." + ) + + shapes = set() + + if out_sci is None: + if fillval == "INDEF": + raise ValueError( + "Fill value 'INDEF' is not allowed when 'out_sci' is None." + ) + else: + shapes.add(out_sci.shape) + + if out_wht is not None: + shapes.add(out_wht.shape) + + if out_ctx is not None: + if len(out_ctx.shape) == 3: + shapes.add(out_ctx.shape[1:]) + else: + shapes.add(out_ctx.shape) + + if out_shape is None: + if not shapes: + raise ValueError( + "'out_shape' cannot be None when all output arrays are " + "None." + ) + out_shape = shapes[0] + else: + shapes.add(out_shape) + + if len(shapes) == 1: + self.out_shape = shapes.pop() + elif len(shapes) > 1: + raise ValueError( + "Inconsistent data shapes specified: 'out_shape' and/or " + "out_sci, out_wht, out_ctx have different shapes." + ) + else: + raise ValueError( + "Either 'out_shape' and/or out_sci, out_wht, out_ctx must be " + "provided." + ) + + self._alloc_output_arrays(n_images=n_images) + + def _alloc_output_arrays(self, n_images): + # allocate arrays as needed: + if out_sci is None: + self.out_sci = np.empty(self.out_shape, dtype=np.float32) + self.out_sci.fill(float(self.fillval)) + + if out_wht is None: + self.out_wht = np.ones(self.out_shape, dtype=np.float32) + + if out_ctx is None: + n_ctx_planes = (n_images - 1) // CTX_PLANE_BITS + 1 + if n_ctx_planes == 1: + ctx_shape = self.out_shape + else: + ctx_shape = (n_ctx_planes, ) + self.out_shape + self.out_ctx = np.zeros(ctx_shape, dtype=np.int32) + + def _increment_ctx_id(self): + self.ctx_id += 1 + if self.ctx_id < 0: + ValueError("Invalid context image ID") + + self._plane_no = self.ctx_id // CTX_PLANE_BITS + depth = 1 if len(self.out_ctx.shape) == 2 else self.out_ctx.shape[0] + + if self._plane_no >= depth: + # Add a new plane to the context image if planeid overflows + plane = np.zeros_like(self.out_shape, np.int32) + self.outcon = np.append(self.out_ctx, [plane], axis=0) + + return (self._plane_no, self.ctx_id % CTX_PLANE_BITS) + + def add_image(self, data, exptime, pixmap, scale=1.0, + weight_map=None, wht_scale=1.0, pixfrac=1.0, in_units='cps', + xmin=None, xmax=None, ymin=None, ymax=None): + """ + + Resample and add an image to the cumulative output image. Also, update + output total weight image and context images. + + TODO: significantly expand this with examples, especially for + exptime, expsq, and ivm weightings. + + Parameters + ---------- + + data : 2d array + A 2d numpy array containing the input image to be drizzled. + it is an error to not supply an image. + + 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 : array + TODO: add description. + + scale : float, optional + The pixel scale of the input image. Conceptually, this is the + linear dimension of a side of a pixel in the input image, but it + is not limited to this and can be set to change how the drizzling + algorithm operates. + + weight_map : 2d array + A 2d numpy array containing the pixel by pixel weighting. + Must have the same dimensions as ``data``. + + TODO: I think this is wrong: If none is supplied, the weghting is set to one. + + 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. + + """ + plane_no, id_in_plane = self._increment_ctx_id() + + # Ensure that the fillval parameter gets properly interpreted + # for use with tdriz + if in_units == 'cps': + expscale = 1.0 + else: + expscale = exptime + + self.exptime += exptime + + data = np.asarray(data, dtype=np.float32) + + # TODO: this is code that would enable initializer to not need output + # image shape at all and set output image shape based on + # output coordinates from the pixmap. If this is enabled, then + # some checks in __init__ may be removed: + # if self.out_shape is None: + # pmap_xmin = int(np.floor(np.min(pixmap[:, :, 0]))) + # pmap_xmax = int(np.ceil(np.max(pixmap[:, :, 0]))) + # pmap_ymin = int(np.floor(np.min(pixmap[:, :, 1]))) + # pmap_ymax = int(np.ceil(np.max(pixmap[:, :, 1]))) + # pixmap = pixmap.copy() + # pixmap[:, :, 0] -= pmap_xmin + # pixmap[:, :, 1] -= pmap_ymin + # self.out_shape = ( + # pmap_xmax - pmap_xmin + 1, + # pmap_ymax - pmap_ymin + 1 + # ) + + if xmin is None or xmin < 0: + xmin = 0 + + if ymin is None or ymin < 0: + ymin = 0 + + if xmax > self.out_shape[1] - 1: + xmax = self.out_shape[1] - 1 + + if ymax > self.out_shape[0] - 1: + ymax = self.out_shape[0] - 1 + + if weight_map is not None: + weight_map = np.asarray(weight_map, dtype=np.float32) + + pixmap = np.asarray(pixmap, dtype=np.float64) + + # TODO: probably tdriz should be modified to not return version. + # we should not have git, Python, C, ... versions + + # TODO: While drizzle code in cdrizzlebox.c supports weight_map=None, + # cdrizzleapi.c does not. It should be modified to support this + # for performance reasons. + _vers, nmiss, nskip = cdrizzle.tdriz( + data, + weight_map, + pixmap, + self.out_sci, + self.out_wht, + self.out_ctx[plane_no], + uniqid=id_in_plane, + xmin=xmin, + xmax=xmax, + ymin=ymin, + ymax=ymax, + scale=scale, # scales image intensity. usually equal to pixel scale + pixfrac=pixfrac, + kernel=self.kernel, + in_units=in_units, + expscale=expscale, + wtscale=wht_scale, + fillstr=self.fillval + ) + self._cversion = _vers # TODO: probably not needed + + return nmiss, nskip + + +# pixmap = calc_pixmap.calc_pixmap(blot_wcs, source_wcs) +# pix_ratio = source_wcs.pscale / blot_wcs.pscale +def blot_image(data, pixmap, pix_ratio, exptime, output_pixel_shape, + interp='poly5', sinscl=1.0): + """ + Resample input image using interpolation to an output grid. + Typically, this is used to resample the resampled image that is the + result of multiple applications of ``Drizzle.add_image`` **if it is well + sampled** back to the pixel coordinate system of input (with distorted + WCS) images of ``Drizzle.add_image``. + + Parameters + ---------- + + data : 2D array + Input numpy array of the source image in units of 'cps'. + + pixmap : 3D array + The world coordinate system to resample on. + + output_pixel_shape : tuple of int + A tuple of two integer numbers indicating the dimensions of the output + image ``(Nx, Ny)``. + + pix_ratio : float + Ratio of the input image pixel scale to the ouput image pixel scale. + + exptime : float + The exposure time of the input image. + + interp : str, optional + The type of interpolation used in the resampling. The + possible values are "nearest" (nearest neighbor interpolation), + "linear" (bilinear interpolation), "poly3" (cubic polynomial + interpolation), "poly5" (quintic polynomial interpolation), + "sinc" (sinc interpolation), "lan3" (3rd order Lanczos + interpolation), and "lan5" (5th order Lanczos interpolation). + + sincscl : float, optional + The scaling factor for sinc interpolation. + """ + + out_sci = np.zeros(output_pixel_shape[::-1], dtype=np.float32) + + cdrizzle.tblot(data, pixmap, out_sci, scale=pix_ratio, kscale=1.0, + interp=interp, exptime=exptime, misval=0.0, sinscl=sinscl) + + return out_sci diff --git a/drizzle/util.py b/drizzle/util.py deleted file mode 100644 index 090de2ba..00000000 --- a/drizzle/util.py +++ /dev/null @@ -1,256 +0,0 @@ -import numpy as np - - -def find_keyword_extn(fimg, keyword, value=None): - """ - This function will return the index of the extension in a multi-extension - FITS file which contains the desired keyword with the given value. - - Parameters - ---------- - - fimg : hdulist - A list of header data units - - keyword : str - The keyword to search for - - value : str or number, optional - If set, the value the keyword must have to match - - Returns - ------- - - The index of the extension - """ - - i = 0 - extnum = -1 - # Search through all the extensions in the FITS object - for chip in fimg: - hdr = chip.header - # Check to make sure the extension has the given keyword - if keyword in hdr: - if value is not None: - # If it does, then does the value match the desired value - # MUST use 'str.strip' to match against any input string! - if hdr[keyword].strip() == value: - extnum = i - break - else: - extnum = i - break - i += 1 - # Return the index of the extension which contained the - # desired EXTNAME value. - return extnum - - -def get_extn(fimg, extn=''): - """ - Returns the FITS extension corresponding to extension specified in - filename. Defaults to returning the first extension with data or the - primary extension, if none have data. - - Parameters - ---------- - - fimg : hdulist - A list of header data units - - extn : str - The extension name and version to match - - Returns - ------- - - The matching header data unit - """ - - if extn: - try: - _extn = parse_extn(extn) - if _extn[0]: - _e = fimg.index_of(_extn) - else: - _e = _extn[1] - - except KeyError: - _e = None - - if _e is None: - _extn = None - else: - _extn = fimg[_e] - - else: - # If no extension is provided, search for first extension - # in FITS file with data associated with it. - - # Set up default to point to PRIMARY extension. - _extn = fimg[0] - # then look for first extension with data. - for _e in fimg: - if _e.data is not None: - _extn = _e - break - - return _extn - - -def get_keyword(fimg, keyword, default=None): - """ - Return a keyword value from the header of an image, - or the default if the keyword is not found. - - Parameters - ---------- - - fimg : hdulist - A list of header data units - - keyword : hdulist - The keyword value to search for - - default : str or number, optional - The default value if not found - - Returns - ------- - - The value if found or default if not - """ - - value = None - if keyword: - _nextn = find_keyword_extn(fimg, keyword) - try: - value = fimg[_nextn].header[keyword] - except KeyError: - value = None - - if value is None and default is not None: - value = default - - return value - - -def is_blank(value): - """ - Determines whether or not a value is considered 'blank'. - - Parameters - ---------- - - value : str - The value to check - - Returns - ------- - - True or False - """ - return value.strip() == "" - - -def parse_extn(extn=''): - """ - Parse a string representing a qualified fits extension name as in the - output of parse_filename and return a tuple (str(extname), int(extver)), - which can be passed to pyfits functions using the 'ext' kw. - Default return is the first extension in a fits file. - - Examples - -------- - >>> parse_extn('sci,2') - ('sci', 2) - >>> parse_extn('2') - ('', 2) - >>> parse_extn('sci') - ('sci', 1) - - Parameters - ---------- - extn : str - The extension name - - Returns - ------- - A tuple of the extension name and value - """ - if not extn: - return ('', 0) - - try: - lext = extn.split(',') - except Exception: - return ('', 1) - - if len(lext) == 1 and lext[0].isdigit(): - return ("", int(lext[0])) - elif len(lext) == 2: - return (lext[0], int(lext[1])) - else: - return (lext[0], 1) - - -def parse_filename(filename): - """ - Parse out filename from any specified extensions. - Returns rootname and string version of extension name. - - Parameters - ---------- - - filename : str - The filename to be parsed - - Returns - ------- - - A tuple with the filename root and extension - """ - # Parse out any extension specified in filename - _indx = filename.find('[') - if _indx > 0: - # Read extension name provided - _fname = filename[:_indx] - _extn = filename[_indx + 1:-1] - else: - _fname = filename - _extn = '' - - return _fname, _extn - - -def set_pscale(the_wcs): - """ - Calculates the plate scale from cdelt and the pc matrix and adds - it to the WCS. Plate scale is not part of the WCS standard, but is - required by the drizzle code - - Parameters - ---------- - - the_wcs : wcs - A WCS object - """ - try: - cdelt = the_wcs.wcs.get_cdelt() - pc = the_wcs.wcs.get_pc() - - except Exception: - try: - # for non-celestial axes, get_cdelt doesnt work - cdelt = the_wcs.wcs.cd * the_wcs.wcs.cdelt - except AttributeError: - cdelt = the_wcs.wcs.cdelt - - try: - pc = the_wcs.wcs.pc - except AttributeError: - pc = 1 - - pccd = np.array(cdelt * pc) - scales = np.sqrt((pccd ** 2).sum(axis=0, dtype=float)) - the_wcs.pscale = scales[0] diff --git a/drizzle/calc_pixmap.py b/drizzle/utils.py similarity index 100% rename from drizzle/calc_pixmap.py rename to drizzle/utils.py