PRO dandia_subtract, imdata, im_bpm, refdata, ref_bpm, masterflat, mflat_bpm, subframes_x, subframes_y, grow, hksize, hires_hksize, ps_var, back_var, diffpro, $
; Description: This module determines the convolution kernel solution (spatially variant or invariant) matching
; a pair of images of the same field. The module also creates a difference image using the derived
; solution. The derived convolution kernel matches a reference image "refdata" with another image
; "imdata", where both images are registered such that, in the worst case, only a translation
; is required to fully align the images. The kernel is modelled as a discrete pixel array of side
; "(2*hksize) + 1" pixels, and the kernel pixel values are solved for directly using linear least
; squares. The pixel kernel model is sufficiently flexible to correct for image misalignments
; consisting of simple translations, hence the relaxed image registration requirement of this
; module. Spatial variance of the kernel solution is achieved by solving for individual pixel
; kernels in a grid of image regions and interpolating the solutions to obtain the kernel at any
; one pixel. For more information about the method please see the paper by Daniel Bramich
; (Bramich D.M., 2008, MNRAS, 386, 77).
; If the parameter "niter" is greater than "1", then the module will iterate each kernel solution,
; using the image model to determine the image pixel variances for the second iteration and beyond.
; For the third iteration and beyond, a 4-sigma clip algorithm will be applied to the residuals
; to identify and reject outlier image pixel values. The module will stop iterating the current
; kernel solution when no more image pixels are rejected or when the maximum number of iterations
; "niter" has been reached.
; The module also supplies the option of solving for a kernel containing lower resolution kernel
; elements. A lower resolution kernel element consists of a 3x3 array of "single-pixel" kernel
; pixels. See the parameters "hksize" and "hires_hksize" for more details. Note that the module
; solves for a circular kernel shape in all cases.
; The module returns the information about the kernel solution for each image region via the IDL
; structure "reg_ksol", including the photometric scale factor and differential background for each
; region. The module also returns the image of kernel solutions "ker_im" representing the kernel
; solution for each image region in a grid layout that has exactly the same form as the image region
; layout. Finally the module also returns the model for the image data "im_model", the difference
; image "diffim", and the bad pixel mask "diffim_bpm" corresponding to both "im_model" and "diffim".
; This module will use C++ code if possible to achieve a faster execution than the default IDL
; code. The C++ code can be up to a factor of ~5 times faster than the IDL code when solving for the
; kernel and differential background, and it can be up to a factor of ~??? times faster than the IDL
; code when constructing the image model.
;
; Input Parameters:
;
; imdata - FLOAT/DOUBLE ARRAY - An image array for which the difference image is required.
; im_bpm - BYTE/INTEGER/LONG ARRAY - A bad pixel mask array of the same size as "imdata" which flags bad pixels
; in "imdata" with a value of "0" and good pixels with a value of "1". If
; this array does not have the same dimensions as "imdata", then all image
; pixels will be considered good, and "im_bpm" will be set to an array of
; the same dimensions as "imdata" with all values equal to "1".
; refdata - FLOAT/DOUBLE ARRAY - An image array of the same dimensions as "imdata" that represents the
; reference image.
; ref_bpm - BYTE/INTEGER/LONG ARRAY - A bad pixel mask array of the same size as "imdata" which flags bad pixels
; in "refdata" with a value of "0" and good pixels with a value of "1". If
; this array does not have the same dimensions as "imdata", then all reference
; image pixels will be considered good, and "ref_bpm" will be set to an array
; of the same dimensions as "imdata" with all values equal to "1". Note that
; bad pixels flagged in this bad pixel mask will be grown by the kernel size
; and shape during the kernel solution process.
; masterflat - FLOAT/DOUBLE ARRAY - An image array of the same dimensions as "imdata" that represents the
; normalised master flat field used to calibrate the image data "imdata".
; Note that if a geometric transformation was applied to "imdata" to
; register the image array with the reference image "refdata", then the
; master flat field image "masterflat" supplied here must also be the
; result of the same transformation applied to the original master flat
; field image. If this array does not have the same dimensions as "imdata",
; then all image pixels will be considered to be of equal sensitivity, and
; "masterflat" will be set to an array of the same dimensions as "imdata"
; with all values equal to "1.0". Any zero or negative values will be
; flagged as bad pixels and set to the value "-1.0" to avoid division by
; zero.
; mflat_bpm - BYTE/INTEGER/LONG ARRAY - A bad pixel mask array of the same size as "imdata" which flags bad
; pixels in "masterflat" with a value of "0" and good pixels with a value
; of "1". If this array does not have the same dimensions as "imdata", then
; all master flat field pixels will be considered good, and "mflat_bpm" will
; be set to an array of the same dimensions as "imdata" with all values equal
; to "1".
; subframes_x - INTEGER/LONG - The number of subdivisions of the image in the x direction used in defining the
; grid of image regions and kernel solutions (Default value = 1). If this input
; parameter is zero or negative, then it will be reset to the default value.
; subframes_y - INTEGER/LONG - The number of subdivisions of the image in the y direction used in defining the
; grid of image regions and kernel solutions (Default value = 1). If this input
; parameter is zero or negative, then it will be reset to the default value.
; grow - FLOAT/DOUBLE - This parameter controls the amount of overlap between the image regions used for
; determining the kernel solutions and it is specified as a fraction of the size of an
; image region before expansion (Default value = 0.0 : Maximum value = 1.0). After
; expansion, the final size of an image region in any one direction is given by
; "(2*grow + 1)*floor(S/M)" where S is the image size along the relevant axis (pix) and
; M is the number of grid subdivisions along the same axis. If this input parameter is
; negative, then it will be reset to the default value of "0.0", and if this parameter
; is greater than "1.0", then it will be reset to the value of "1.0".
; hksize - INTEGER/LONG - The half kernel size (pix) such that the side of the kernel pixel array is of length
; "(2*hksize) + 1" pixels.
; hires_hksize - INTEGER/LONG - The kernel radius (pix) for the single-pixel kernel elements such that the kernel
; diameter for the single-pixel kernel elements is given by "(2*hires_hksize) + 1"
; pixels. If this parameter has a negative or a zero value, then the kernel array
; will consist exclusively of nine-pixel kernel elements. If this parameter has a
; value that is greater than or equal to "hksize", then the kernel array will
; consist exclusively of single-pixel kernel elements. If this parameter has a
; value that is greater than zero but less than "hksize", then the kernel array
; will consist of a mixture of single-pixel and nine-pixel kernel elements.
; ps_var - INTEGER/LONG - If this parameter is set to "0", then the module assumes that the photometric scale
; factor is spatially invariant and that it is fully specified by a single value. If this
; parameter is set to "1", then the module assumes that the photometric scale factor is
; spatially variant, and the module will determine the photometric scale factor for each
; pixel in the image model "im_model". The photometric scale factor is important when
; producing the difference image, and unless there is variable extinction across the image
; area (or the flat fielding was not done correctly), then it should be spatially invariant.
; Hence the default value is "0".
; back_var - INTEGER/LONG - If this parameter is set to "0", then the module assumes that the differential sky
; background is spatially invariant and that it is fully specified by a single value. If
; this parameter is set to "1", then the module assumes that the differential sky
; background is spatially variant, and the module will determine the differential sky
; background for each pixel in the image model "im_model". The differential sky background
; is important when producing the difference image, and generally there is a variable sky
; background across the image area. Hence the default value is "1".
; diffpro - INTEGER/LONG - This parameter is a switch controlling the method used to produce the image model
; "im_model". The following parameter values are accepted (Default value = 0):
;
; (i) A value of "0" means that for each image pixel, the corresponding kernel is formed as
; a bilinear interpolation of the kernel solutions for the nearest image regions to the
; image pixel. If "ps_var" is set to "1", then the corresponding photometric scale
; factor is determined in the same way. If "back_var" is set to "1", then the
; corresponding differential sky background is determined in the same way.
;
; (ii) A value of "1" means that for each group of image pixels of size equal to a tenth of
; the size of an image region used for determining the kernel solution, the
; corresponding kernel is formed as a bilinear interpolation of the kernel solutions
; for the nearest image regions to the current group of image pixels. If "ps_var" is
; set to "1", then the corresponding photometric scale factor is determined in the same
; way. If "back_var" is set to "1", then the corresponding differential sky background
; is determined in the same way. This method is much faster than the one outlined in
; (ii), and it is similar in accuracy if the assumption that the kernel solution is
; approximately constant over the group of image pixels is valid.
;
; ron - FLOAT/DOUBLE - The CCD readout noise (ADU). A negative value of "ron" will be reset to "0.0".
; gain - FLOAT/DOUBLE - The CCD gain (e-/ADU). A zero or negative value of "gain" will be reset to "1.0".
; coeff2 - FLOAT/DOUBLE - The quadratic coefficient in the CCD linearisation equation:
;
; Xnew = X + (coeff2 * X^2) + (coeff3 * X^3)
;
; where X represents the image counts (ADU) after overscan and bias correction, and
; Xnew represents the non-linearity corrected image counts (ADU).
; coeff3 - FLOAT/DOUBLE - The cubic coefficient in the CCD linearisation equation:
;
; Xnew = X + (coeff2 * X^2) + (coeff3 * X^3)
;
; where X represents the image counts (ADU) after overscan and bias correction, and
; Xnew represents the non-linearity corrected image counts (ADU).
; niter - INTEGER/LONG - The maximum number of iterations to be performed when solving for the kernel and
; differential background (Default value = 3). A negative value of "niter" will be reset to
; the default value "3". Note that for the first iteration, the image pixel variances
; are calculated using the image pixel values in the CCD noise model, which introduces a bias
; into the kernel solution for this iteration. For the second iteration, the image pixel
; variances are calculated using the image model created from the first iteration in the CCD
; noise model, and hence the kernel solution for this iteration is less biased. For the third
; iteration and beyond, the image pixel variances are calculated using the image model created
; from the previous iteration in the CCD noise model, and a 4-sigma clip algorithm employing
; these variances is applied to the residuals from the previous iteration, rejecting outlier
; image pixel values from the current iteration. By setting "niter" to at least "3", one is
; guaranteeing that the image pixel variances have been calculated from the image model and
; not the image data, and that most outlier image pixel values have been excluded from the
; kernel solution. The iterative kernel solution process will stop when no more image pixels
; are rejected or when the maximum number of iterations "niter" has been reached.
; DanIDL_C_code - STRING - The full directory path indicating where the DanIDL C++ code is installed. The shared
; library of DanIDL C++ routines should exist as "/dist/libDanIDL.so"
; within the installation.
;
; Output Parameters:
;
; reg_ksol - IDL STRUCTURE - An IDL structure containing tags that record the details about the kernel solution
; for each image region. This structure contains the following tags, with one record
; for each image region:
;
; rlx - LONG - The x-axis pixel index of the lower left-hand corner of the rectangular
; image region.
; rux - LONG - The x-axis pixel index of the upper right-hand corner of the rectangular
; image region.
; rly - LONG - The y-axis pixel index of the lower left-hand corner of the rectangular
; image region.
; ruy - LONG - The y-axis pixel index of the upper right-hand corner of the rectangular
; image region.
; xc - DOUBLE - The mean x coordinate of the good pixels in the image region.
; yc - DOUBLE - The mean y coordinate of the good pixels in the image region.
; ngoodpix - LONG - The number of good pixels in the image region.
; npix_out - LONG - The final number of image pixels used in the kernel and differential
; background solution for the image region.
; npix_rej - LONG - The number of image pixels that were rejected while determining the
; kernel and differential background solution for the image region.
; ps - DOUBLE - The photometric scale factor which matches the flux in the corresponding
; reference image region with the flux in the input image region.
; back - DOUBLE - The differential background (ADU) which matches the sky background in
; the corresponding reference image region with the sky background in
; the input image region.
; chisqpp - DOUBLE - The reduced chi squared (or chi squared per degree of freedom) of
; the kernel and differential background solution for the image
; region. The number of degrees of freedom is given by
; "npix_out - knpix - 1", where "knpix" is the number of kernel
; elements.
; klx - LONG - The x-axis pixel index in the image of kernel solutions "ker_im" of the
; lower left-hand corner of the kernel pixel array that represents the
; kernel solution for the image region.
; kux - LONG - The x-axis pixel index in the image of kernel solutions "ker_im" of the
; upper right-hand corner of the kernel pixel array that represents the
; kernel solution for the image region.
; kly - LONG - The y-axis pixel index in the image of kernel solutions "ker_im" of the
; lower left-hand corner of the kernel pixel array that represents the
; kernel solution for the image region.
; kuy - LONG - The y-axis pixel index in the image of kernel solutions "ker_im" of the
; upper right-hand corner of the kernel pixel array that represents the
; kernel solution for the image region.
; tag - INTEGER - If this flag is set to a value of "1", then the kernel solution for
; the image region was successful, and if it is set to a value of "0",
; then the module failed in determining a kernel solution for the image
; region.
;
; nregions - LONG - The number of image regions used when determining the kernel solution and creating the
; difference image. This parameter is also the number of entries for each tag in the IDL
; structure "reg_ksol".
; ker_im - DOUBLE ARRAY - An image array of dimensions "subframes_x*ksize" by "subframes_y*ksize" pixels, where
; "ksize" is given by "(2*hksize) + 1", that stores the grid of normalised kernel
; solutions. The layout of the grid of kernel solutions corresponds directly with the
; layout of the image regions used to solve for the kernels. To extract the kernel
; solution for the image region "i", where the image region "i" is defined by
; "[reg_ksol[i].rlx:reg_ksol[i].rux, reg_ksol[i].rly:reg_ksol[i].ruy]", simply extract
; the region "[reg_ksol[i].klx:reg_ksol[i].kux, reg_ksol[i].kly:reg_ksol[i].kuy]" from
; "ker_im". Note that each kernel solution is normalised and has a sum of exactly "1.0".
; im_model - DOUBLE ARRAY - An image array of the the same dimensions as "imdata" that represents the model for
; the image data "imdata". The image model "im_model" is constructed by convolving
; the reference image with the kernel solution(s), scaling by the photometric scale
; factor, and adding in the differential background.
; diffim - DOUBLE ARRAY - An image array of the the same dimensions as "imdata" that represents the difference
; image, formed by simply subtracting the model for the image data "im_model" from the
; image data "imdata" itself.
; diffim_bpm - INTEGER ARRAY - A bad pixel mask array of the same size as "imdata" which flags bad pixels in
; both the image model "im_model" and the difference image "diffim" with a value of
; "0", and good pixels with a value of "1".
; status - INTEGER - If the module successfully determined the kernel solution and created the difference image,
; then "status" is returned with a value of "1", otherwise it is returned with a value of "0".
;
; Keywords:
;
; If the keyword QUIET is set (as "/QUIET"), then the program will not print any progress messages on screen.
; This functionality is useful if the module is to be called many times within an image processing routine for
; instance.
;
; Author: Dan Bramich (dan.bramich@hotmail.co.uk)
;
; History:
;
; 17/04/2009 - The module was failing to accept a supplied master flat field image of the correct dimensions and
; number type. The bug is now fixed (dmb).
; 16/09/2008 - Module created (dmb).
;
; This module does the following, treating bad pixels correctly at all stages:
;
; 1) Checks that the input image array "imdata" is an image of the correct number type, and determines its
; dimensions. If "imdata" is not an image of the correct number type, then the module returns without doing
; anything.
;
; 2) Checks that the bad pixel mask "im_bpm" has the same dimensions as "imdata" and is of the correct number type.
; If this is not the case, then the module creates a bad pixel mask that flags every pixel in "imdata" as good,
; using this new bad pixel mask to replace the values stored in "im_bpm".
;
; 3) Checks that the reference image array "refdata" is an image of the correct number type and has the same
; dimensions as "imdata". If this is not the case, then the module returns without doing anything.
;
; 4) Checks that the reference image bad pixel mask "ref_bpm" has the same dimensions as "imdata" and is of the
; correct number type. If this is not the case, then the module creates a reference image bad pixel mask that
; flags every pixel in "refdata" as good, using this new bad pixel mask to replace the values stored in "ref_bpm".
;
; 5) Checks that the normalised master flat field image "masterflat" is an image of the correct number type and has
; the same dimensions as "imdata". If this is not the case, then the module creates a master flat field image
; under the assumption that every image pixel has the same sensitivity, using this new master flat field image
; to replace the values stored in "masterflat".
;
; 6) Checks that the master flat field bad pixel mask "mflat_bpm" has the same dimensions as "imdata" and is of the
; correct number type. If this is not the case, then the module creates a master flat field bad pixel mask that
; flags every pixel in "masterflat" as good, using this new bad pixel mask to replace the values stored in
; "mflat_bpm".
;
; 7) Replaces any negative or zero values in the master flat field image "masterflat" with the value "-1.0", and
; flags these pixels as bad pixels in the master flat field bad pixel mask "mflat_bpm".
;
; 8) Checks the values of some of the input parameters, resetting them to default values where necessary, and
; returning without doing anything when parameter values are out of range.
;
; 9) Splits the input image into a grid of "subframes_x" by "subframes_y" regions with the required overlap as
; specified by the input parameter "grow". A kernel and differential background solution will be determined
; independently for each image region.
;
; 10) Defines the size and shape of the kernel pixel array, and sets up an image to hold the kernel solutions.
;
; 11) Checks that theoretically each image region contains enough pixels in order to solve for the kernel.
;
; 12) Grows the bad pixels in the reference image bad pixel mask "ref_bpm" by the kernel size and shape.
;
; 13) Creates a master bad pixel mask that includes the bad pixels from "im_bpm", the grown bad pixels from "ref_bpm"
; and the bad pixels from "mflat_bpm".
;
; 14) For each image region for which a kernel solution is required, the module carries out steps (15)-(19), and then
; moves on to step (20).
;
; 15) Extracts the current region from the input image, and calculates the mean x and y coordinates of the good
; pixels in the current image region.
;
; 16) Checks that the current image region contains enough good pixels in order to solve for the kernel. If this is
; not the case, then the module moves on to the next image region in step (14), and records a failed kernel
; solution for the current image region.
;
; 17) Extracts the current region from both the reference image and the master flat field image.
;
; 18) Solves for the kernel pixel array, photometric scale factor and differential sky background for the current
; image region by calling the DanIDL module "dandia_ksol.pro". For more details, please see the module itself.
;
; 19) If the kernel solution process for the current image region was successful, then the module stores the kernel
; pixel array solution, photometric scale factor and differential sky background for the current image region.
; If the kernel solution process for the current image region was not successful, then the module moves on to the
; next image region in step (14), and records a failed kernel solution for the current image region.
;
; 20) Checks that the kernel solution process was successful for at least one image region. If this is not the case,
; then the module reports an error and exits.
;
; 21) Constructs a model for the input image using the kernel and differential background solutions. This is done by
; convolving the reference image with the normalised kernel solution (spatially invariant for an image split into
; only one image region, and spatially variant for an image split into multiple image regions), multiplying by the
; photometric scale factor and adding the differential background solution (again spatially invariant or variant).
; The DanIDL module "dandia_model_image.pro" is used for this step.
;
; 22) Constructs the difference image, which is simply the input image minus the model for the input image that was
; constructed in step (21).
;If the keyword QUIET is set, then do not print progress messages
;Check that "imdata" is an image
;If "im_bpm" does not have the same dimensions as "imdata", then assume that all image pixels are good
;Check that "refdata" is an image of the same dimensions as "imdata"
;If "ref_bpm" does not have the same dimensions as "imdata", then assume that all reference image pixels are good
;If "masterflat" is not an image of the same dimensions as "imdata", then assume that all image pixels are of equal
;sensitivity
;If "mflat_bpm" does not have the same dimensions as "imdata", then assume that all master flat field image
;pixels are good
;Ensure that the master flat field contains no zero values, and that negative and zero values are considered
;as bad pixels
;Check some of the input parameters, and reset them if necessary
;Define the image regions to be used in determining the kernel solution
;Define the size and shape of the kernel pixel array, and set up an image to hold the kernel solutions
;Check that theoretically the image regions contain enough pixels in order to solve for the kernel
;Grow the bad pixels in the reference image bad pixel mask by the kernel size and shape
;Create a master bad pixel mask that includes the bad pixels from "im_bpm", the grown bad pixels from
;"ref_bpm" and the bad pixels from "mflat_bpm"
;For each image region for which a kernel solution is required
;Extract the current region limits
;Extract the current region from the input image
;Calculate the mean x and y coordinates of the good pixels in the current region
;Check that the current region contains enough good pixels in order to solve for the kernel
;Extract the current region from the master flat field image
;Extract the current region from the reference image
;Solve for the kernel pixel array, photometric scale factor and differential sky background
;Save the kernel pixel array, photometric scale factor and differential sky background
;Check that the kernel solution process succeeded for at least one image region
;Construct a model for the current image
;Create the difference image
;Finish