Index: /trunk/ppSub/src/ppSubArguments.c =================================================================== --- /trunk/ppSub/src/ppSubArguments.c (revision 14106) +++ /trunk/ppSub/src/ppSubArguments.c (revision 14107) @@ -17,7 +17,5 @@ { fprintf(stderr, "\nPan-STARRS PSF-matched image subtraction\n\n"); - fprintf(stderr, "Usage: %s INPUT.fits REFERENCE.fits OUTPUT_ROOT \n" - "\t[-psf REFERENCE.psf.fits|-psflist REFERENCE.list]\n", - program); + fprintf(stderr, "Usage: %s INPUT.fits REFERENCE.fits OUTPUT_ROOT\n", program); fprintf(stderr, "\n"); psArgumentHelp(arguments); @@ -57,20 +55,4 @@ } \ psMetadataAdd##TYPE(config->arguments, PS_LIST_TAIL, RECIPENAME, 0, NULL, value); \ -} - -// Get a mask value from the command-line or recipe, and add it to the arguments -#define VALUE_ARG_RECIPE_MASK(ARGNAME, RECIPENAME) { \ - bool mdok; \ - const char *name = psMetadataLookupStr(&mdok, arguments, ARGNAME); \ - if (!mdok || !name || strlen(name) == 0) { \ - name = psMetadataLookupStr(NULL, recipe, RECIPENAME); \ - if (!name) { \ - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to find %s in recipe %s", \ - RECIPENAME, PPSUB_RECIPE); \ - goto ERROR; \ - } \ - } \ - psMaskType value = pmConfigMask(name, config); \ - psMetadataAddU8(config->arguments, PS_LIST_TAIL, RECIPENAME, 0, NULL, value); \ } @@ -152,28 +134,8 @@ } -// Add a single filename to the arguments as an array, so that it can be used with pmFPAfileBindFromArgs, etc -static void fileList(const char *file, // The symbolic name for the file - const char *name, // The name of the file - const char *comment, // Description of the file - pmConfig *config // Configuration - ) -{ - psArray *files = psArrayAlloc(1); // Array with file names - files->data[0] = psStringCopy(name); - psMetadataAddArray(config->arguments, PS_LIST_TAIL, file, 0, comment, files); - psFree(files); - return; -} - -////////////////////////////////////////////////////////////////////////////////////////////////////////////// bool ppSubArguments(int argc, char *argv[], pmConfig *config) { assert(config); - - psMetadataAddBool(config->arguments, PS_LIST_TAIL, "PHOTOMETRY", 0, "Do photometry?", - (psArgumentGet(argc, argv, "-psf") || psArgumentGet(argc, argv, "-psflist")) ? - true : false); - pmConfigFileSetsMD(config->arguments, &argc, argv, "PSPHOT.PSF", "-psf", "-psflist"); psMetadata *arguments = psMetadataAlloc(); // Command-line arguments @@ -182,11 +144,13 @@ psMetadataAddStr(arguments, PS_LIST_TAIL, "-refmask", 0, "Referencemask image", NULL); psMetadataAddStr(arguments, PS_LIST_TAIL, "-refweight", 0, "Referenceweight image", NULL); - psMetadataAddStr(arguments, PS_LIST_TAIL, "-stats", 0, "Statistics file", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-stat", 0, "Statistics file", NULL); psMetadataAddS32(arguments, PS_LIST_TAIL, "-size", 0, "Kernel half-size (pixels)", 0); psMetadataAddS32(arguments, PS_LIST_TAIL, "-order", 0, "Spatial polynomial order", 0); - psMetadataAddStr(arguments, PS_LIST_TAIL, "-type", 0, "Kernel type (POIS or ISIS)", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-type", 0, "Kernel type (ISIS|POIS|SPAM|FRIES)", NULL); psMetadataAddStr(arguments, PS_LIST_TAIL, "-isis-widths", 0, "ISIS Gaussian widths (comma-separated)", NULL); psMetadataAddStr(arguments, PS_LIST_TAIL, "-isis-orders", 0, "ISIS polynomial orders (comma-separated)", NULL); + psMetadataAddS32(arguments, PS_LIST_TAIL, "-inner", 0, "SPAM and FRIES inner radius", 0); + psMetadataAddS32(arguments, PS_LIST_TAIL, "-spam-binning", 0, "SPAM kernel binning", 2); psMetadataAddF32(arguments, PS_LIST_TAIL, "-spacing", 0, "Typical stamp spacing (pixels)", NAN); psMetadataAddS32(arguments, PS_LIST_TAIL, "-footprint", 0, "Stamp footprint half-size (pixels)", 0); @@ -202,27 +166,19 @@ } - fileList("INPUT", argv[1], "Name of the input image", config); - fileList("REF", argv[2], "Name of the reference image", config); + psArray *files = psArrayAlloc(1); // Array with file names + files->data[0] = psStringCopy(argv[1]); + psMetadataAddArray(config->arguments, PS_LIST_TAIL, "INPUT", 0, "Name of the input image", files); + psFree(files); + files = psArrayAlloc(1); + files->data[0] = psStringCopy(argv[2]); + psMetadataAddArray(config->arguments, PS_LIST_TAIL, "REF", 0, "Name of the reference image", files); + psFree(files); psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "Name of the output image", argv[3]); - const char *inMask = psMetadataLookupStr(NULL, arguments, "-inmask"); // Name of input mask - if (inMask && strlen(inMask) > 0) { - fileList("INPUT.MASK", inMask, "Name of the input mask image", config); - } - const char *inWeight = psMetadataLookupStr(NULL, arguments, "-inweight"); // Name of input weight - if (inWeight && strlen(inWeight) > 0) { - fileList("INPUT.WEIGHT", inWeight, "Name of the input weight image", config); - } - - const char *refMask = psMetadataLookupStr(NULL, arguments, "-refmask"); // Name of reference mask - if (refMask && strlen(refMask) > 0) { - fileList("REF.MASK", refMask, "Name of the reference mask image", config); - } - const char *refWeight = psMetadataLookupStr(NULL, arguments, "-refweight"); // Name of reference weight - if (refWeight && strlen(refWeight) > 0) { - fileList("REF.WEIGHT", refWeight, "Name of the reference weight image", config); - } - - valueArgStr(config, arguments, "-stats", "STATS", config->arguments); + valueArgStr(config, arguments, "-inmask", "INPUT.MASK", config->arguments); + valueArgStr(config, arguments, "-inweight", "INPUT.WEIGHT", config->arguments); + valueArgStr(config, arguments, "-refmask", "REF.MASK", config->arguments); + valueArgStr(config, arguments, "-refweight", "REF.WEIGHT", config->arguments); + valueArgStr(config, arguments, "-stat", "STATS", config->arguments); psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim @@ -232,13 +188,15 @@ } - VALUE_ARG_RECIPE_INT("-size", "KERNEL.SIZE", S32, 0); - VALUE_ARG_RECIPE_INT("-order", "SPATIAL.ORDER", S32, 0); - VALUE_ARG_RECIPE_FLOAT("-spacing", "STAMP.SPACING", F32); - VALUE_ARG_RECIPE_INT("-footprint", "STAMP.FOOTPRINT", S32, 0); - VALUE_ARG_RECIPE_FLOAT("-threshold", "STAMP.THRESHOLD", F32); - VALUE_ARG_RECIPE_INT("-iter", "ITER", S32, 0); - VALUE_ARG_RECIPE_FLOAT("-rej", "REJ", F32); - VALUE_ARG_RECIPE_MASK("-mask-bad", "MASK.BAD"); - VALUE_ARG_RECIPE_MASK("-mask-blank", "MASK.BLANK"); + VALUE_ARG_RECIPE_INT("-size", "KERNEL.SIZE", S32, 0); + VALUE_ARG_RECIPE_INT("-order", "SPATIAL.ORDER", S32, 0); + VALUE_ARG_RECIPE_FLOAT("-spacing", "STAMP.SPACING", F32); + VALUE_ARG_RECIPE_INT("-inner", "INNER", S32, 0); + VALUE_ARG_RECIPE_INT("-spam-binning", "SPAM.BINNING", S32, 0); + VALUE_ARG_RECIPE_INT("-footprint", "STAMP.FOOTPRINT", S32, 0); + VALUE_ARG_RECIPE_FLOAT("-threshold", "STAMP.THRESHOLD", F32); + VALUE_ARG_RECIPE_INT("-iter", "ITER", S32, 0); + VALUE_ARG_RECIPE_FLOAT("-rej", "REJ", F32); + VALUE_ARG_RECIPE_INT("-mask-bad", "MASK.BAD", U8, 0); + VALUE_ARG_RECIPE_INT("-mask-blank", "MASK.BLANK", U8, 0); vectorArgRecipe(config, arguments, "-isis-widths", recipe, "ISIS.WIDTHS", config->arguments, PS_TYPE_F32); @@ -269,4 +227,8 @@ } else if (strcasecmp(type, "ISIS") == 0) { kernelType = PM_SUBTRACTION_KERNEL_ISIS; + } else if (strcasecmp(type, "SPAM") == 0) { + kernelType = PM_SUBTRACTION_KERNEL_SPAM; + } else if (strcasecmp(type, "FRIES") == 0) { + kernelType = PM_SUBTRACTION_KERNEL_FRIES; } else { psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unrecognised kernel type: %s", type); Index: /trunk/ppSub/src/ppSubReadout.c =================================================================== --- /trunk/ppSub/src/ppSubReadout.c (revision 14106) +++ /trunk/ppSub/src/ppSubReadout.c (revision 14107) @@ -6,7 +6,9 @@ #include #include -#include #include "ppSub.h" + +#define MASK_BAD 0x01 // Mask value for bad pixel +#define MASK_STAMP 0x02 // Mask value for bad stamp (and bad stamp region) bool ppSubReadout(pmConfig *config, const pmFPAview *view) @@ -14,6 +16,8 @@ pmReadout *inRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); // Input readout pmReadout *refRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF"); // Reference readout +#if 0 pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT"); // Output cell pmReadout *outRO = pmReadoutAlloc(outCell); // Output readout +#endif psImage *input = inRO->image; // Input image @@ -34,19 +38,35 @@ psVector *widths = psMetadataLookupPtr(NULL, config->arguments, "ISIS.WIDTHS"); // ISIS Gaussian widths psVector *orders = psMetadataLookupPtr(NULL, config->arguments, "ISIS.ORDERS"); // ISIS Polynomial orders + int inner = psMetadataLookupS32(NULL, config->arguments, "INNER"); // Inner radius + int binning = psMetadataLookupS32(NULL, config->arguments, "SPAM.BINNING"); // Binning for SPAM kernel psMaskType maskBad = psMetadataLookupU8(NULL, config->arguments, "MASK.BAD"); // Value to mask psMaskType maskBlank = psMetadataLookupU8(NULL, config->arguments, "MASK.BLANK"); // Mask for blank reg. - int numCols = input->numCols, numRows = input->numRows; // Image dimensions - if (!inRO->mask) { - inRO->mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); - psImageInit(inRO->mask, 0); + if (!inRO->mask && !pmReadoutGenerateMask(inRO)) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate mask for input image"); + return false; } - if (!refRO->mask) { - refRO->mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); - psImageInit(refRO->mask, 0); + if (!inRO->weight && !pmReadoutGenerateWeight(inRO, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate weight map for input image"); + return false; + } + if (!refRO->mask && !pmReadoutGenerateMask(refRO)) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate mask for reference image"); + return false; + } + if (!refRO->weight && !pmReadoutGenerateWeight(refRO, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate weight map for reference image"); + return false; } - // Mask for subtraction - psImage *subMask = pmSubtractionMask(inRO->mask, refRO->mask, maskBad, size, footprint); + // Worried about the masks for bad pixels and bad stamps colliding, so make our own mask + int numCols = input->numCols, numRows = input->numRows; // Image dimensions + psImage *stampMask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); // Mask to use for stamps + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + stampMask->data.PS_TYPE_MASK_DATA[y][x] = + (refRO->mask->data.PS_TYPE_MASK_DATA[y][x] & maskBad) ? MASK_BAD : 0; + } + } #if 0 @@ -66,6 +86,6 @@ #endif - pmSubtractionKernels *kernels = pmSubtractionKernelsGenerate(type, size, order, - widths, orders); // Kernel basis functions + pmSubtractionKernels *kernels = pmSubtractionKernelsGenerate(type, size, order, widths, orders, + inner, binning); // Kernel basis functions psArray *stamps = NULL; // Stamps for matching PSF psVector *solution = NULL; // Solution to match PSF @@ -73,5 +93,6 @@ int numRejected = -1; // Number of rejected stamps in each iteration for (int i = 0; i < iter && numRejected != 0; i++) { - stamps = pmSubtractionFindStamps(stamps, refRO->image, subMask, threshold, spacing); + stamps = pmSubtractionFindStamps(stamps, refRO->image, stampMask, MASK_BAD, MASK_STAMP, + threshold, spacing, size + footprint); if (!stamps) { psError(PS_ERR_UNKNOWN, false, "Unable to find stamps on reference image."); @@ -91,5 +112,5 @@ } - numRejected = pmSubtractionRejectStamps(stamps, refRO->image, inRO->image, subMask, + numRejected = pmSubtractionRejectStamps(stamps, refRO->image, inRO->image, stampMask, MASK_STAMP, solution, footprint, rej, kernels); if (numRejected < 0) { @@ -99,44 +120,26 @@ psLogMsg("ppSub", PS_LOG_INFO, "%d stamps rejected on iteration %d.", numRejected, i); } - - if (numRejected > 0) { - solution = pmSubtractionSolveEquation(solution, stamps); - if (!solution) { - psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); - goto ERROR; - } - } + psFree(stampMask); psImage *convImage = NULL, *convWeight = NULL, *convMask = NULL; // Convolved images - if (!pmSubtractionConvolve(&convImage, &convWeight, &convMask, refRO->image, refRO->weight, subMask, - maskBlank, solution, kernels)) { + if (!pmSubtractionConvolve(&convImage, &convWeight, &convMask, + refRO->image, refRO->weight, refRO->mask, + MASK_BAD, maskBlank, solution, kernels)) { psError(PS_ERR_UNKNOWN, false, "Unable to convolve reference image."); goto ERROR; } - psFree(subMask); // Do the subtraction if (reverse) { - outRO->image = (psImage*)psBinaryOp(NULL, convImage, "-", inRO->image); + (void)psBinaryOp(inRO->image, convImage, "-", inRO->image); } else { - outRO->image = (psImage*)psBinaryOp(NULL, inRO->image, "-", convImage); + (void)psBinaryOp(inRO->image, inRO->image, "-", convImage); } - outRO->mask = (psImage*)psBinaryOp(NULL, convMask, "|", inRO->mask); - if (convWeight) { - outRO->weight = (psImage*)psBinaryOp(NULL, convWeight, "+", inRO->weight); - } + (void)psBinaryOp(inRO->mask, convMask, "|", inRO->mask); + (void)psBinaryOp(inRO->weight, convWeight, "+", inRO->weight); psFree(convImage); psFree(convMask); psFree(convWeight); - - outRO->data_exists = true; - outCell->data_exists = true; - outCell->parent->data_exists = true; - - if (!pmFPACopyConcepts(outCell->parent->parent, inRO->parent->parent->parent)) { - psError(PS_ERR_UNKNOWN, false, "Unable to copy concepts from input to output."); - return false; - } #if 0 @@ -173,23 +176,14 @@ #endif - psFree(kernels); - psFree(stamps); - psFree(solution); - if (psMetadataLookupBool(NULL, config->arguments, "PHOTOMETRY")) { - pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT"); - pmFPACopy(photFile->fpa, inRO->parent->parent->parent); + psFree(kernels); + psFree(stamps); + psFree(solution); + return true; - psphotReadout(config, view); - - pmFPAfileActivate(config->files, false, "PSPHOT.INPUT"); - } - - return true; - -ERROR: - psFree(kernels); - psFree(stamps); - psFree(solution); - return false; + ERROR: + psFree(kernels); + psFree(stamps); + psFree(solution); + return false; }