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 (dmb@ing.iac.es) ; ; History: ; ; 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