PRO dandia_ksol, imreg, refreg, imreg_bpm, mflatreg, hksize, hires_knpix, hires_kxind, hires_kyind, lores_knpix, lores_kxind, lores_kyind, $ ; Description: This module determines the spatially invariant convolution kernel matching a pair of images of ; the same field. The derived convolution kernel matches a reference image "refreg" with another ; image "imreg", 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. 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 the 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 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. The normal "single-pixel" kernel pixels and lower resolution "nine-pixel" kernel pixels ; to be solved for are specified via vectors of pixel indices in the x and y directions. The use ; of pixel indices to specify the kernel elements allows the flexibility to solve for any kernel ; shape within the kernel pixel array of side "(2*hksize) + 1" pixels. ; The module returns the kernel solution "ker_out", the photometric scale factor "ps_out", the ; differential background "back_out", the number of image pixels used "npix_out", the number of ; image pixels rejected "npix_rej", and the chi squared per degree of freedom "chisqpp_out". ; 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 in this case. ; ; Input Parameters: ; ; imreg - FLOAT/DOUBLE ARRAY - An image array of dimensions Nx by Ny pixels. ; refreg - FLOAT/DOUBLE ARRAY - An image array of dimensions "Nx + (2*hksize)" by "Ny + (2*hksize)" pixels, ; where "hksize" is the kernel radius in pixels (see below). This image array ; is the reference image, and the central region of this image array of size ; Nx by Ny pixels must already have been registered, at least approximately, ; with the image array "imreg". The extra border of width "hksize" is necessary ; due to the nature of the kernel convolution process. ; imreg_bpm - BYTE/INTEGER/LONG ARRAY - A bad pixel mask array of the same size as "imreg" which flags bad ; pixels with a value of "0" and good pixels with a value of "1". If ; this array does not have the same dimensions as "imreg", then all ; image pixels will be considered good, and "imreg_bpm" will be set to ; an array of the same dimensions as "imreg" with all values equal to ; "1". Note that this bad pixel mask must include bad pixels from ; "imreg", bad pixels from "refreg" grown by the kernel size and shape, ; and bad pixels from "mflatreg" (see next parameter). ; mflatreg - FLOAT/DOUBLE ARRAY - An image array of the same size as "imreg" representing the normalised ; master flat field used to calibrate the image data "imreg". Note that if ; a geometric transformation was applied to "imreg" to register the image ; array with the reference image "refreg", then the master flat field image ; "mflatreg" 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 "imreg", then all image pixels will be considered ; to be of equal sensitivity, and "mflatreg" will be set to an array of the ; same dimensions as "imreg" with all values equal to "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_knpix - INTEGER/LONG - The number of "single-pixel" kernel pixels in the kernel pixel array. ; hires_kxind - INTEGER/LONG VECTOR - A one-dimensional vector of pixel indices along the x-axis where each ; element corresponds to a "single-pixel" kernel pixel. Hence there are ; "hires_knpix" elements in this vector. ; hires_kyind - INTEGER/LONG VECTOR - A one-dimensional vector of pixel indices along the y-axis where each ; element corresponds to a "single-pixel" kernel pixel. Hence there are ; "hires_knpix" elements in this vector. ; lores_knpix - INTEGER/LONG - The number of "nine-pixel" kernel pixels in the kernel pixel array. A "nine-pixel" ; kernel pixel consists of a 3x3 array of "single-pixel" kernel pixels, and acts as ; a low resolution kernel element (it can also be referred to as a basis function). ; lores_kxind - INTEGER/LONG VECTOR - A one-dimensional vector of pixel indices along the x-axis where each ; element corresponds to the central pixel of a "nine-pixel" kernel pixel. ; Hence there are "lores_knpix" elements in this vector. ; lores_kyind - INTEGER/LONG VECTOR - A one-dimensional vector of pixel indices along the y-axis where each ; element corresponds to the central pixel of a "nine-pixel" kernel pixel. ; Hence there are "lores_knpix" elements in this vector. ; 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 zero or 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: ; ; ker_out - DOUBLE ARRAY - The output normalised kernel pixel array of side "(2*hksize) + 1" pixels. The sum of ; this array is "1.0" exactly. ; ps_out - DOUBLE - The photometric scale factor which matches the flux in the reference image "refreg" with the ; flux in the image "imreg". ; back_out - DOUBLE - The differential background (ADU) which matches the sky background in the reference image ; "refreg" with the sky background in the image "imreg". ; npix_out - LONG - The final number of image pixels used in the kernel and differential background solution. ; npix_rej - LONG - The number of image pixels that were rejected while solving for the kernel and differential ; background. ; chisqpp_out - DOUBLE - The reduced chi squared (or chi squared per degree of freedom) of the kernel and ; differential background solution. The number of degrees of freedom is given by ; "npix_out - hires_knpix - lores_knpix - 1". ; status - INTEGER - If the module successfully solved for the kernel pixel array, the photometric scale factor ; and the differential background, 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: ; ; 05/09/2008 - Module created (dmb) ; ; This module does the following, treating bad pixels correctly at all stages: ; ; 1) Checks that the input image array "imreg" is an image of the correct number type, and determines its dimensions. ; If "imreg" is not an image of the correct number type, then the module returns without doing anything. ; ; 2) 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. ; ; 3) Checks that the reference image array "refreg" is an image of the correct dimensions and number type. If this ; is not the case, then the module returns without doing anything. ; ; 4) Checks that the bad pixel mask "imreg_bpm" has the same dimensions as "imreg" and is of the correct number type. ; If this is not the case, then the program creates a bad pixel mask that flags every pixel in "imreg" as good, ; using this new bad pixel mask to replace the values stored in "imreg_bpm". ; ; 5) Checks that the master flat field image "mflatreg" has the same dimensions as "imreg" and is of the correct ; number type. If this is not the case, then the program assumes that all image pixels are of the same sensitivity, ; and "mflatreg" will be set to an array of the same dimensions as "imreg" with all values equal to "1.0". ; ; 6) Checks the values and number types (and dimensions where necessary) of the parameters "hires_knpix", "hires_kxind", ; "hires_kyind", "lores_knpix", "lores_kxind" and "lores_kyind" specifying the size, shape and composition of the ; kernel to be solved for. ; ; 7) Calculates an initial set of inverse pixel variances based on the image data "imreg". ; ; 8) If "hires_knpix" is greater than zero and "lores_knpix" is equal to zero, then the module carries out steps ; (9)-(16) and then finishes. If "hires_knpix" is equal to zero and "lores_knpix" is greater than zero, then the ; module carries out steps (17)-(25) and then finishes. If "hires_knpix" is greater than zero and "lores_knpix" is ; greater than zero, then the module carries out steps (26)-(34) and then finishes. ; ; 9) Repeats steps (10)-(15) until the number of iterations reaches "niter", or the current number of image pixels ; being used to solve for the kernel becomes less than or equal to the number of kernel elements plus one ; (the number of free parameters), or the current number of rejected pixels is equal to zero. The module then ; moves on to step (16). ; ; 10) Calculates the least square matrix and vector on the right hand side. ; ; 11) Calculates the LU decomposition of the matrix and solves the system of linear equations. ; ; 12) Constructs the current kernel and differential background solution. ; ; 13) Convolves the reference image "refreg" with the current kernel solution and adds the differential background, ; to create a current model for the image data. ; ; 14) Uses the current model for the image data to calculate a new set of inverse pixel variances. ; ; 15) If this is not the first iteration, then the module will perform 4-sigma clipping on the residuals, ; rejecting a maximum of 1% of image pixels with absolute residuals greater than or equal to 4 sigma. ; ; 16) Reports the success or failure of the kernel solution process, and in the case of success, the module returns ; the normalised kernel pixel array, the photometric scale factor, the differential background, the final number ; of pixels used in the kernel solution, the number of pixels rejected during the kernel solution process, and ; the reduced chi squared. ; ; 17) Creates a boxcar smoothed version of the reference image "refreg" using a 3x3 pixel kernel with all values ; equal to "1.0". ; ; 18) Repeats steps (19)-(24) until the number of iterations reaches "niter", or the current number of image pixels ; being used to solve for the kernel becomes less than or equal to the number of kernel elements plus one ; (the number of free parameters), or the current number of rejected pixels is equal to zero. The module then ; moves on to step (25). ; ; 19) Calculates the least square matrix and vector on the right hand side. ; ; 20) Calculates the LU decomposition of the matrix and solves the system of linear equations. ; ; 21) Constructs the current kernel and differential background solution. ; ; 22) Convolves the reference image "refreg" with the current kernel solution and adds the differential background, ; to create a current model for the image data. ; ; 23) Uses the current model for the image data to calculate a new set of inverse pixel variances. ; ; 24) If this is not the first iteration, then the module will perform 4-sigma clipping on the residuals, ; rejecting a maximum of 1% of image pixels with absolute residuals greater than or equal to 4 sigma. ; ; 25) Reports the success or failure of the kernel solution process, and in the case of success, the module returns ; the normalised kernel pixel array, the photometric scale factor, the differential background, the final number ; of pixels used in the kernel solution, the number of pixels rejected during the kernel solution process, and ; the reduced chi squared. ; ; 26) Creates a boxcar smoothed version of the reference image "refreg" using a 3x3 pixel kernel with all values ; equal to "1.0". ; ; 27) Repeats steps (28)-(33) until the number of iterations reaches "niter", or the current number of image pixels ; being used to solve for the kernel becomes less than or equal to the number of kernel elements plus one ; (the number of free parameters), or the current number of rejected pixels is equal to zero. The module then ; moves on to step (34). ; ; 28) Calculates the least square matrix and vector on the right hand side. ; ; 29) Calculates the LU decomposition of the matrix and solves the system of linear equations. ; ; 30) Constructs the current kernel and differential background solution. ; ; 31) Convolves the reference image "refreg" with the current kernel solution and adds the differential background, ; to create a current model for the image data. ; ; 32) Uses the current model for the image data to calculate a new set of inverse pixel variances. ; ; 33) If this is not the first iteration, then the module will perform 4-sigma clipping on the residuals, ; rejecting a maximum of 1% of image pixels with absolute residuals greater than or equal to 4 sigma. ; ; 34) Reports the success or failure of the kernel solution process, and in the case of success, the module returns ; the normalised kernel pixel array, the photometric scale factor, the differential background, the final number ; of pixels used in the kernel solution, the number of pixels rejected during the kernel solution process, and ; the reduced chi squared. ;If the keyword QUIET is set, then do not print progress messages ;Check that "imreg" is an image ;Check some of the input parameters, and reset them if necessary ;Determine if the shared library of DanIDL C++ routines is installed ;Check that "refreg" is an image of the correct size ;If "imreg_bpm" does not have the same dimensions as "imreg" or the correct number type, then assume ;that all image pixels are good ;If "mflatreg" does not have the same dimensions as "imreg", then assume that all image pixels are of ;equal sensitivity ;Check the parameters specifying the type of kernel solution that is required ;Calculate an initial set of inverse pixel variances based on the image data ;Set up some necessary variables ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;Solve for a kernel array consisting exclusively of "single-pixel" kernel pixels ;Set up some upper indices ;Iterate the kernel solution process, using a 4-sigma clip algorithm on the residuals and only rejecting ;a maximum of 1% of pixels during any one iteration ;If the shared library of DanIDL C++ routines is installed, then use the C++ code to speed up ;the calculation of the kernel and differential background solution ;Solve for the kernel and differential background solution ;If the shared library of DanIDL C++ routines is not installed, then use the default IDL code ;Set up the relevant least squares matrix and vector on the right hand side ;Pre-calculate a temporary region ;Calculate various matrix elements ;Copy across symmetrical elements ;Calculate a vector element ;Calculate the remaining four matrix elements and two vector elements ;Calculate the LU decomposition of the matrix ;Solve the system of linear equations ;Construct the current kernel and differential background solution ;Convolve the reference image with the current kernel solution and add the differential background ;Calculate the reduced chi-squared ;Use the current model image to calculate a new set of inverse pixel variances ;If this is not the first iteration, then perform sigma clipping on the residuals ;Calculate the normalised residuals for the current kernel and differential background solution ;Reject a maximum of 1% of pixels with absolute residuals greater than or equal to 4 sigma ;If too many outlier image pixels have been found during the kernel solution process, then report ;that the kernel solution process failed ;If the kernel and differential background solutions are acceptable, then determine the photometric ;scale factor and normalise the kernel ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;Solve for a kernel array consisting exclusively of "nine-pixel" kernel pixels ;Set up some upper indices ;Create a boxcar smoothed version of the reference image ;Iterate the kernel solution process, using a 4-sigma clip algorithm on the residuals and only rejecting ;a maximum of 1% of pixels during any one iteration ;If the shared library of DanIDL C++ routines is installed, then use the C++ code to speed up ;the calculation of the kernel and differential background solution ;Solve for the kernel and differential background solution ;If the shared library of DanIDL C++ routines is not installed, then use the default IDL code ;Set up the relevant least squares matrix and vector on the right hand side ;Pre-calculate a temporary region ;Calculate various matrix elements ;Copy across symmetrical elements ;Calculate a vector element ;Calculate the remaining four matrix elements and two vector elements ;Calculate the LU decomposition of the matrix ;Solve the system of linear equations ;Construct the current kernel and differential background solution ;Convolve the reference image with the current kernel solution and add the differential background ;Calculate the reduced chi-squared ;Use the current model image to calculate a new set of inverse pixel variances ;If this is not the first iteration, then perform sigma clipping on the residuals ;Calculate the normalised residuals for the current kernel and differential background solution ;Reject a maximum of 1% of pixels with absolute residuals greater than or equal to 4 sigma ;If too many outlier image pixels have been found during the kernel solution process, then report ;that the kernel solution process failed ;If the kernel and differential background solutions are acceptable, then determine the photometric ;scale factor and normalise the kernel ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;Solve for a kernel array consisting of a mixture of "single-pixel" and "nine-pixel" kernel pixels ;Set up some upper indices ;Create a boxcar smoothed version of the reference image ;Iterate the kernel solution process, using a 4-sigma clip algorithm on the residuals and only rejecting ;a maximum of 1% of pixels during any one iteration ;If the shared library of DanIDL C++ routines is installed, then use the C++ code to speed up ;the calculation of the kernel and differential background solution ;Solve for the kernel and differential background solution ;If the shared library of DanIDL C++ routines is not installed, then use the default IDL code ;Set up the relevant least squares matrix and vector on the right hand side, starting with ;those elements relating to the "single-pixel" kernel pixels ;Pre-calculate a temporary region ;Calculate the matrix elements relating to the "single-pixel" kernel pixels ;Calculate a vector element relating to the "single-pixel" kernel pixels ;Calculate some more matrix and vector elements relating to the "single-pixel" kernel pixels ;Calculate the matrix and vector elements relating to the "nine-pixel" kernel pixels ;Pre-calculate a temporary region ;Calculate the matrix elements relating to the "nine-pixel" kernel pixels ;Calculate a vector element relating to the "nine-pixel" kernel pixels ;Calculate some more matrix and vector elements relating to the "nine-pixel" kernel pixels ;Calculate the final matrix element and the final vector element ;Copy across symmetrical elements ;Calculate the LU decomposition of the matrix ;Solve the system of linear equations ;Construct the current kernel and differential background solution ;Convolve the reference image with the current kernel solution and add the differential background ;Calculate the reduced chi-squared ;Use the current model image to calculate a new set of inverse pixel variances ;If this is not the first iteration, then perform sigma clipping on the residuals ;Calculate the normalised residuals for the current kernel and differential background solution ;Reject a maximum of 1% of pixels with absolute residuals greater than or equal to 4 sigma ;If too many outlier image pixels have been found during the kernel solution process, then report ;that the kernel solution process failed ;If the kernel and differential background solutions are acceptable, then determine the photometric ;scale factor and normalise the kernel ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;Finish ;Finish