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 (dan.bramich@hotmail.co.uk)
;
; History:
;
; 24/04/2009 - Further module optimisation performed (dmb).
; 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