Index: /trunk/doc/modules/ModulesSDRS.tex =================================================================== --- /trunk/doc/modules/ModulesSDRS.tex (revision 3426) +++ /trunk/doc/modules/ModulesSDRS.tex (revision 3427) @@ -1,3 +1,3 @@ -%%% $Id: ModulesSDRS.tex,v 1.34 2005-03-15 23:18:29 price Exp $ +%%% $Id: ModulesSDRS.tex,v 1.35 2005-03-15 23:19:36 price Exp $ \documentclass[panstarrs]{panstarrs} @@ -1555,14 +1555,14 @@ \begin{verbatim} psImage *pmCombineImages(psImage *combined, // Combined image - psArray **questionablePixels, // Array of rejection masks - const psArray *images, // Array of input images - const psArray *errors, // Array of input error images - const psArray *masks,// Array of input masks - unsigned int maskVal, // Mask value - const psArray *pixels, // Pixels to combine - int numIter, // Number of rejection iterations - float sigmaClip, // Number of standard deviations at which to reject - const psStats *stats // Statistics to use in the combination - ); + psArray **questionablePixels, // Array of rejection masks + const psArray *images, // Array of input images + const psArray *errors, // Array of input error images + const psArray *masks,// Array of input masks + unsigned int maskVal, // Mask value + const psArray *pixels, // Pixels to combine + int numIter, // Number of rejection iterations + float sigmaClip, // Number of standard deviations at which to reject + const psStats *stats // Statistics to use in the combination + ); \end{verbatim} @@ -1601,10 +1601,10 @@ \begin{verbatim} psArray *pmRejectPixels(const psArray *images, // Array of input images - const psArray *pixels, // These are the pixels which were rejected in the combination - const psArray *inToOut, // Transformations from input to output system - const psArray *outToIn, // Transformations from output to input system - float rejThreshold, // Rejection threshold - float gradLimit // Gradient limit - ); + const psArray *pixels, // These are the pixels which were rejected in the combination + const psArray *inToOut, // Transformations from input to output system + const psArray *outToIn, // Transformations from output to input system + float rejThreshold, // Rejection threshold + float gradLimit // Gradient limit + ); \end{verbatim} @@ -1669,53 +1669,53 @@ for (int i = 0; i < nImages; i++) { - psPixels *blanks = NULL; // List of blank pixels - transformed->data[i] = psImageTransform(NULL, &blanks, inputs->data[i], - inputMask->data[i], inputMaskVal, NAN, skyToDetector, - combineRegion, NULL, PS_INTERPOLATE_BILINEAR); - transformedErr->data[i] = psImageTransform(NULL, NULL, inputsErr->data[i], inputMask->data[i], - inputMaskVal, NAN, skyToDetector, combineRegion, NULL, - PS_INTERPOLATE_BILINEAR_VARIANCE); - psImage *skyImage = transformed->data[i]; // Dereference the transformed image - psRegion *blankRegion = psRegionAlloc(0, 0, skyImage->numCols, skyImage->numRows); // Size of - // transformed - // image - transformedMask->data[i] = psPixelsToMask(NULL, blanks, blankRegion, PS_MASK_BLANK); - psFree(blankRegion); - psFree(blanks); + psPixels *blanks = NULL; // List of blank pixels + transformed->data[i] = psImageTransform(NULL, &blanks, inputs->data[i], + inputMask->data[i], inputMaskVal, NAN, skyToDetector, + combineRegion, NULL, PS_INTERPOLATE_BILINEAR); + transformedErr->data[i] = psImageTransform(NULL, NULL, inputsErr->data[i], inputMask->data[i], + inputMaskVal, NAN, skyToDetector, combineRegion, NULL, + PS_INTERPOLATE_BILINEAR_VARIANCE); + psImage *skyImage = transformed->data[i]; // Dereference the transformed image + psRegion *blankRegion = psRegionAlloc(0, 0, skyImage->numCols, skyImage->numRows); // Size of + // transformed + // image + transformedMask->data[i] = psPixelsToMask(NULL, blanks, blankRegion, PS_MASK_BLANK); + psFree(blankRegion); + psFree(blanks); } - psArray *rejected = NULL; // Array of rejected pixel lists + psArray *rejected = NULL; // Array of rejected pixel lists psStats *combineStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // Statistic to use in doing the combination psImage *combined = pmCombineImages(NULL, &rejected, transformed, transformedErr, transformedMask, 0, - NULL, numIter, sigmaClip, combineStats); // Combined image + NULL, numIter, sigmaClip, combineStats); // Combined image psArray *bad = pmRejectPixels(inputs, rejected, NULL, skyToDetector, rejThreshold, gradLimit); // Bad pix - psPixels *combinePixels = NULL; // Pixels to combine + psPixels *combinePixels = NULL; // Pixels to combine for (int i = 0; i < nImages; i++) { - psPixels *badSource = psPixelsTransform(NULL, bad->data[i], skyToDetector); // Bad pixels on the input - psImage *badMask = psPixelsToMask(NULL, badSource, PS_MASK_COSMICRAY); // Mask image for the input - (void)psBinaryOp(inputMask->data[i], inputMask->data[i], "|", badMask); // Put CRs into original mask - psFree(badSource); - psFree(badMask); - - combinePixels = psPixelsConcatenate(redo, bad->data[i]); - - // Update transformed image - psPixels *blanks = NULL; // List of blank pixels - transformed->data[i] = psImageTransform(transformed->data[i], &blanks, inputs->data[i], - inputMask->data[i], inputMaskVal | PS_MASK_COSMICRAY, NAN, - skyToDetector, combineRegion, bad->data[i], - PS_INTERPOLATE_BILINEAR); - transformedErr->data[i] = psImageTransform(transformedErr->data[i], NULL, inputsErr->data[i], - inputMask->data[i], inputMaskVal | PS_MASK_COSMICRAY, - NAN, skyToDetector, combineRegion, bad->data[i], - PS_INTERPOLATE_BILINEAR_VARIANCE); - psImage *skyImage = transformed->data[i]; // Dereference the transformed image - psRegion *blankRegion = psRegionAlloc(0, 0, skyImage->numCols, skyImage->numRows); // Size of - // transformed - // image - transformedMask->data[i] = psPixelsToMask(transformedMask->data[i], blanks, blankRegion, - PS_MASK_BLANK); - psFree(blankRegion); - psFree(blanks); + psPixels *badSource = psPixelsTransform(NULL, bad->data[i], skyToDetector); // Bad pixels on the input + psImage *badMask = psPixelsToMask(NULL, badSource, PS_MASK_COSMICRAY); // Mask image for the input + (void)psBinaryOp(inputMask->data[i], inputMask->data[i], "|", badMask); // Put CRs into original mask + psFree(badSource); + psFree(badMask); + + combinePixels = psPixelsConcatenate(redo, bad->data[i]); + + // Update transformed image + psPixels *blanks = NULL; // List of blank pixels + transformed->data[i] = psImageTransform(transformed->data[i], &blanks, inputs->data[i], + inputMask->data[i], inputMaskVal | PS_MASK_COSMICRAY, NAN, + skyToDetector, combineRegion, bad->data[i], + PS_INTERPOLATE_BILINEAR); + transformedErr->data[i] = psImageTransform(transformedErr->data[i], NULL, inputsErr->data[i], + inputMask->data[i], inputMaskVal | PS_MASK_COSMICRAY, + NAN, skyToDetector, combineRegion, bad->data[i], + PS_INTERPOLATE_BILINEAR_VARIANCE); + psImage *skyImage = transformed->data[i]; // Dereference the transformed image + psRegion *blankRegion = psRegionAlloc(0, 0, skyImage->numCols, skyImage->numRows); // Size of + // transformed + // image + transformedMask->data[i] = psPixelsToMask(transformedMask->data[i], blanks, blankRegion, + PS_MASK_BLANK); + psFree(blankRegion); + psFree(blanks); } psFree(bad); @@ -1723,5 +1723,5 @@ // Combine with no rejection combined = pmCombineImages(combined, NULL, transformed, transformedErr, transformedMask, - PS_MASK_BLANK, combinePixels, 0, 0.0, combineStats); + PS_MASK_BLANK, combinePixels, 0, 0.0, combineStats); psFree(combineStats); psFree(combinePixels); @@ -1763,6 +1763,6 @@ \begin{verbatim} typedef enum { - PM_SUBTRACTION_KERNEL_POIS, // POIS kernel --- delta functions - PM_SUBTRACTION_KERNEL_ISIS // ISIS kernel --- gaussians modified by polynomials + PM_SUBTRACTION_KERNEL_POIS, // POIS kernel --- delta functions + PM_SUBTRACTION_KERNEL_ISIS // ISIS kernel --- gaussians modified by polynomials } pmSubtractionKernelsType; \end{verbatim} @@ -1774,10 +1774,10 @@ \begin{verbatim} typedef struct { - pmSubtractionKernelType type; // Type of kernels --- allowing the use of multiple kernels - psVector *u, *v; // Offset (for POIS) or polynomial order (for ISIS) - psVector *sigma; // Width of Gaussian (for ISIS) - psVector *xOrder, *yOrder; // Spatial polynomial order (for all) + pmSubtractionKernelType type; // Type of kernels --- allowing the use of multiple kernels + psVector *u, *v; // Offset (for POIS) or polynomial order (for ISIS) + psVector *sigma; // Width of Gaussian (for ISIS) + psVector *xOrder, *yOrder; // Spatial polynomial order (for all) int subIndex; // Index of kernel to be subtracted to maintain flux conservation - psArray *preCalc; // Array of images containing pre-calculated kernel (to + psArray *preCalc; // Array of images containing pre-calculated kernel (to // accelerate ISIS; don't use for POIS) } pmSubtractionKernels; @@ -1846,8 +1846,8 @@ \begin{verbatim} typedef struct { - int x, y; // Position - psImage *matrix; // Associated matrix - psVector *vector; // Associated vector - psStampStatus status; // Status of stamp + int x, y; // Position + psImage *matrix; // Associated matrix + psVector *vector; // Associated vector + psStampStatus status; // Status of stamp } pmStamp; \end{verbatim} @@ -1867,8 +1867,8 @@ \begin{verbatim} typedef enum { - PM_STAMP_USED, // Use this stamp - PM_STAMP_REJECTED, // This stamp has been rejected - PM_STAMP_RECALC, // Having been reset, this stamp needs to be recalculated - PM_STAMP_NONE // No stamp in this region + PM_STAMP_USED, // Use this stamp + PM_STAMP_REJECTED, // This stamp has been rejected + PM_STAMP_RECALC, // Having been reset, this stamp needs to be recalculated + PM_STAMP_NONE // No stamp in this region } pmStampStatus; \end{verbatim} @@ -1876,11 +1876,11 @@ \begin{verbatim} psArary *pmSubtractionFindStamps(psArray *stamps, // Output stamps, or NULL - const psImage *image, // Image for which to find stamps - const psImage *mask, // Mask - unsigned int maskVal, // Value for mask - float threshold, // Threshold for stamps in the image - int xNum, int yNum, // Number of stamps in x and y - int border // Border around image to ignore (should be size of kernel) - ); + const psImage *image, // Image for which to find stamps + const psImage *mask, // Mask + unsigned int maskVal, // Value for mask + float threshold, // Threshold for stamps in the image + int xNum, int yNum, // Number of stamps in x and y + int border // Border around image to ignore (should be size of kernel) + ); \end{verbatim} @@ -1913,9 +1913,9 @@ \begin{verbatim} bool pmSubtractionCalculateEquation(psArray *stamps, // The stamps for which to calculate the equation - const psImage *reference, // Reference image - const psImage *input, // Input image - const psSubtractionKernels *kernels, // The kernel basis functions - int footprint // Half-size of region over which to calculate equation - ); + const psImage *reference, // Reference image + const psImage *input, // Input image + const psSubtractionKernels *kernels, // The kernel basis functions + int footprint // Half-size of region over which to calculate equation + ); \end{verbatim} @@ -1959,7 +1959,7 @@ \begin{verbatim} -psVector *pmSubtractionSolveEquation(psVector *solution, // Solution vector, or NULL - const psArray *stamps // Array of stamps - ); +psVector *pmSubtractionSolveEquation(psVector *solution, // Solution vector, or NULL + const psArray *stamps // Array of stamps + ); \end{verbatim} @@ -1981,13 +1981,13 @@ \begin{verbatim} bool pmSubtractionRejectStamps(psArray *stamps, // Array of stamps to check for rejection - psImage *mask, // Mask image - unsigned int badStampMaskVal, // Value to use in mask for bad stamp - int footprint, // Region to mask if stamp is bad - float sigmaRej, // Number of RMS deviations above zero at which to reject - const psImage *refImage, // Reference image - const psImage *inImage, // Input image - const psVector *solution, // Solution vector - const pmSubtractionKernels *kernels // Array of kernel parameters - ); + psImage *mask, // Mask image + unsigned int badStampMaskVal, // Value to use in mask for bad stamp + int footprint, // Region to mask if stamp is bad + float sigmaRej, // Number of RMS deviations above zero at which to reject + const psImage *refImage, // Reference image + const psImage *inImage, // Input image + const psVector *solution, // Solution vector + const pmSubtractionKernels *kernels // Array of kernel parameters + ); \end{verbatim} @@ -2059,35 +2059,35 @@ // Generate kernel basis functions - psArray *kernels = NULL; // Array of kernel basis functions + psArray *kernels = NULL; // Array of kernel basis functions switch (kernelType) { case PM_SUBTRACTION_KERNEL_POIS: - // Create the kernel basis functions - kernels = pmSubtractionKernelsGeneratePOIS(kernelHalfSize, spatialOrder); - break; + // Create the kernel basis functions + kernels = pmSubtractionKernelsGeneratePOIS(kernelHalfSize, spatialOrder); + break; case PM_SUBTRACTION_KERNEL_ISIS: - kernels = pmSubtractionKernelsGenerateISIS(sigmas, polyOrders, kernelHalfSize, spatialOrder); - break; + kernels = pmSubtractionKernelsGenerateISIS(sigmas, polyOrders, kernelHalfSize, spatialOrder); + break; default: - barf(); + barf(); } - psArray *stamps = NULL; // Array of stamps - psVector *kernelCoeffs = NULL; // Coefficients for the kernels - bool rejected = true; // Did we reject a stamp in the last iteration? + psArray *stamps = NULL; // Array of stamps + psVector *kernelCoeffs = NULL; // Coefficients for the kernels + bool rejected = true; // Did we reject a stamp in the last iteration? // Iterate for a solution for (int iter = 0; iter < numIter && rejected; iter++) { - // Find stamps - stamps = pmSubtractionFindStamps(stamps, reference, subMask, maskVal | PS_MASK_BAD_STAMP, - stampThreshold, nStampsX, nStampsY, stampSize, kernelHalfSize); - - // Generate and solve matrix equations - (void)pmSubtractionCalculateEquation(stamps, reference, input, kernels, stampSize); - kernelCoeffs = pmSubtractionSolveEquation(kernelCoeffs, stamps); - - // Reject bad stamps - rejected = pmSubtractionRejectStamps(stamps, subMask, PS_MASK_BAD_STAMP, stampSize, sigmaRej, - reference, input, kernelCoeffs, kernels); + // Find stamps + stamps = pmSubtractionFindStamps(stamps, reference, subMask, maskVal | PS_MASK_BAD_STAMP, + stampThreshold, nStampsX, nStampsY, stampSize, kernelHalfSize); + + // Generate and solve matrix equations + (void)pmSubtractionCalculateEquation(stamps, reference, input, kernels, stampSize); + kernelCoeffs = pmSubtractionSolveEquation(kernelCoeffs, stamps); + + // Reject bad stamps + rejected = pmSubtractionRejectStamps(stamps, subMask, PS_MASK_BAD_STAMP, stampSize, sigmaRej, + reference, input, kernelCoeffs, kernels); } @@ -2140,5 +2140,5 @@ keep = psRegionAlloc (XBORDER, image->nCol - YBORDER, - YBORDER, image->nRow - YBORDER); + YBORDER, image->nRow - YBORDER); Sky = stats->median;