Index: /trunk/ppStack/src/Makefile.am =================================================================== --- /trunk/ppStack/src/Makefile.am (revision 14625) +++ /trunk/ppStack/src/Makefile.am (revision 14626) @@ -10,5 +10,6 @@ ppStackLoop.c \ ppStackReadout.c \ - ppStackVersion.c + ppStackVersion.c \ + ppStackMatch.c noinst_HEADERS = \ Index: /trunk/ppStack/src/ppStack.h =================================================================== --- /trunk/ppStack/src/ppStack.h (revision 14625) +++ /trunk/ppStack/src/ppStack.h (revision 14626) @@ -34,3 +34,9 @@ ); +/// Convolve image to match specified seeing +bool ppStackMatch(pmReadout *output, ///< Convolved readout + const pmReadout *input, ///< Readout to be convolved + const pmConfig *config ///< Configuration + ); + #endif Index: /trunk/ppStack/src/ppStackArguments.c =================================================================== --- /trunk/ppStack/src/ppStackArguments.c (revision 14625) +++ /trunk/ppStack/src/ppStackArguments.c (revision 14626) @@ -19,10 +19,9 @@ { fprintf(stderr, "\nPan-STARRS Image combination\n\n"); - fprintf(stderr, "Usage: %s INPUTS.mdc OUTPUT_ROOT\n" + fprintf(stderr, "Usage: %s INPUTS.mdc OUTPUT_ROOT -target FWHM -stamps FILENAME\n" "where INPUTS.mdc contains various METADATAs, each with:\n" "\tIMAGE(STR): Image filename\n" "\tMASK(STR): Mask filename\n" "\tWEIGHT(STR) Weight map filename\n" - "\tSEEING(F32): Seeing FWHM (pixels)\n" "\tWEIGHTING(F32): Relative weighting to be applied\n" "\tSCALE(F32): Relative scaling to be applied\n", @@ -103,4 +102,6 @@ psMetadata *arguments = psMetadataAlloc(); // Command-line arguments + psMetadataAddStr(arguments, PS_LIST_TAIL, "-stamps", 0, "Stamps file with x,y,flux per line", NULL); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-target", 0, "Target PSF FWHM", NAN); psMetadataAddStr(arguments, PS_LIST_TAIL, "-stats", 0, "Statistics file", NULL); psMetadataAddS32(arguments, PS_LIST_TAIL, "-iter", 0, "Number of rejection iterations", 0); @@ -110,4 +111,5 @@ psMetadataAddU8(arguments, PS_LIST_TAIL, "-mask-bad", 0, "Mask value for bad pixels", 0); psMetadataAddU8(arguments, PS_LIST_TAIL, "-mask-blank", 0, "Mask value for blank region", 0); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-threshold-mask", 0, "Threshold for mask deconvolution", NAN); // XXX I want to get this from the recipe as well. (same for everything else...) @@ -117,4 +119,19 @@ usage(argv[0], arguments, config); } + + const char *stampsName = psMetadataLookupStr(NULL, arguments, "-stamps"); // Name of stamps file + if (!stampsName || strlen(stampsName) == 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Stamps file not specified with -stamps."); + goto ERROR; + } + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "STAMPS", 0, "Stamps file", stampsName); + + float target = psMetadataLookupF32(NULL, arguments, "-target"); // Target PSF width + if (!isfinite(target)) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Target PSF width not specified with -target."); + goto ERROR; + } + psMetadataAddF32(config->arguments, PS_LIST_TAIL, "TARGET", 0, "Target PSF width", target); + unsigned int numBad = 0; // Number of bad lines @@ -138,10 +155,11 @@ } - VALUE_ARG_RECIPE_INT("-iter", "ITER", S32, 0); - VALUE_ARG_RECIPE_FLOAT("-combine-rej", "COMBINE.REJ", F32); - VALUE_ARG_RECIPE_FLOAT("-convolve-rej", "CONVOLVE.REJ", F32); - VALUE_ARG_RECIPE_FLOAT("-extent", "EXTENT", F32); - VALUE_ARG_RECIPE_MASK("-mask-bad", "MASK.BAD"); - VALUE_ARG_RECIPE_MASK("-mask-blank", "MASK.BLANK"); + VALUE_ARG_RECIPE_INT("-iter", "ITER", S32, 0); + VALUE_ARG_RECIPE_FLOAT("-combine-rej", "COMBINE.REJ", F32); + VALUE_ARG_RECIPE_FLOAT("-convolve-rej", "CONVOLVE.REJ", F32); + VALUE_ARG_RECIPE_FLOAT("-extent", "EXTENT", F32); + VALUE_ARG_RECIPE_MASK("-mask-bad", "MASK.BAD"); + VALUE_ARG_RECIPE_MASK("-mask-blank", "MASK.BLANK"); + VALUE_ARG_RECIPE_FLOAT("-threshold-mask", "THRESHOLD.MASK", F32); psTrace("ppStack", 1, "Done reading command-line arguments\n"); Index: /trunk/ppStack/src/ppStackCamera.c =================================================================== --- /trunk/ppStack/src/ppStackCamera.c (revision 14625) +++ /trunk/ppStack/src/ppStackCamera.c (revision 14626) @@ -40,11 +40,4 @@ psString weight = psMetadataLookupStr(&mdok, input, "WEIGHT"); // Name of weight map - float seeing = psMetadataLookupF32(&mdok, input, "SEEING"); // Seeing FWHM - if (!mdok || isnan(seeing) || seeing == 0.0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Component %s lacks SEEING of type F32", item->name); - psFree(iter); - return false; - } - float weighting = psMetadataLookupF32(&mdok, input, "WEIGHTING"); // Relative weighting if (!mdok || isnan(weighting) || weighting == 0.0) { @@ -122,6 +115,4 @@ } - psMetadataAddF32(imageFile->fpa->analysis, PS_LIST_TAIL, "PPSTACK.SEEING", 0, - "Seeing for image", seeing); psMetadataAddF32(imageFile->fpa->analysis, PS_LIST_TAIL, "PPSTACK.WEIGHTING", 0, "Relative weighting for image", weighting); Index: /trunk/ppStack/src/ppStackMatch.c =================================================================== --- /trunk/ppStack/src/ppStackMatch.c (revision 14626) +++ /trunk/ppStack/src/ppStackMatch.c (revision 14626) @@ -0,0 +1,50 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "ppStack.h" + +//#define TESTING + +bool ppStackMatch(pmReadout *output, const pmReadout *input, const pmConfig *config) +{ + // Look up appropriate values from the ppSub recipe + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe + int size = psMetadataLookupS32(NULL, recipe, "KERNEL.SIZE"); // Kernel half-size + int order = psMetadataLookupS32(NULL, recipe, "SPATIAL.ORDER"); // Spatial polynomial order + float regionSize = psMetadataLookupF32(NULL, recipe, "REGION.SIZE"); // Size of iso-kernel regs + float spacing = psMetadataLookupF32(NULL, recipe, "STAMP.SPACING"); // Typical stamp spacing + int footprint = psMetadataLookupS32(NULL, recipe, "STAMP.FOOTPRINT"); // Stamp half-size + float threshold = psMetadataLookupF32(NULL, recipe, "STAMP.THRESHOLD"); // Threshold for stmps + int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations + float rej = psMetadataLookupF32(NULL, recipe, "REJ"); // Rejection threshold + pmSubtractionKernelsType type = + pmSubtractionKernelsTypeFromString(psMetadataLookupStr(NULL, recipe, "KERNEL.TYPE")); // Kernel type + psVector *widths = psMetadataLookupPtr(NULL, recipe, "ISIS.WIDTHS"); // ISIS Gaussian widths + psVector *orders = psMetadataLookupPtr(NULL, recipe, "ISIS.ORDERS"); // ISIS Polynomial orders + int inner = psMetadataLookupS32(NULL, recipe, "INNER"); // Inner radius + int ringsOrder = psMetadataLookupS32(NULL, recipe, "RINGS.ORDER"); // RINGS polynomial order + int binning = psMetadataLookupS32(NULL, recipe, "SPAM.BINNING"); // Binning for SPAM kernel + psMaskType maskBad = pmConfigMask(psMetadataLookupStr(NULL, recipe, "MASK.BAD"), + config); // Value to mask + psMaskType maskBlank = pmConfigMask(psMetadataLookupStr(NULL, recipe, "MASK.BLANK"), + config); // Mask for blank reg. + + // These values are specified specifically for stacking + const char *stampsName = psMetadataLookupStr(NULL, config->arguments, "STAMPS"); // Filename for stamps + float target = psMetadataLookupF32(NULL, config->arguments, "TARGET"); // Target PSF width + + // Do the image matching + if (!pmSubtractionMatch(output, input, NULL, footprint, regionSize, spacing, threshold, stampsName, + target, type, size, order, widths, orders, inner, ringsOrder, binning, iter, + rej, maskBad, maskBlank)) { + psError(PS_ERR_UNKNOWN, false, "Unable to match images."); + return false; + } + + return true; +} Index: /trunk/ppStack/src/ppStackReadout.c =================================================================== --- /trunk/ppStack/src/ppStackReadout.c (revision 14625) +++ /trunk/ppStack/src/ppStackReadout.c (revision 14626) @@ -15,12 +15,9 @@ { // Get the recipe values - bool mdok; // Status of MD lookup int iter = psMetadataLookupS32(NULL, config->arguments, "ITER"); // Rejection iterations float combineRej = psMetadataLookupF32(NULL, config->arguments, "COMBINE.REJ"); // Combination threshold - float convolveRej = psMetadataLookupF32(NULL, config->arguments, "CONVOLVE.REJ"); // Convolution threshold - float extent = psMetadataLookupF32(NULL, config->arguments, "EXTENT"); // Extent of convolution psMaskType maskBad = psMetadataLookupU8(NULL, config->arguments, "MASK.BAD"); // Value to mask psMaskType maskBlank = psMetadataLookupU8(NULL, config->arguments, "MASK.BLANK"); // Mask for blank reg. - psVector *seeing = psMetadataLookupPtr(&mdok, config->arguments, "SEEING"); // Seeing in each image + float threshold = psMetadataLookupF32(NULL, config->arguments, "THRESHOLD.MASK"); // Threshold for mask deconvolution // Get the output target @@ -41,5 +38,4 @@ int fileNum = 0; // Number of file float totExposure = 0.0; // Total exposure time - pmReadout *templateRO = NULL; // Template readout, for copy WCS psList *fpaList = psListAlloc(NULL); // List of input FPAs, for concept averaging psList *cellList = psListAlloc(NULL); // List of input cells, for concept averaging @@ -53,9 +49,4 @@ bool mdok; // Status of MD lookup - float seeing = psMetadataLookupF32(&mdok, inputFile->fpa->analysis, "PPSTACK.SEEING"); // Seeing FWHM - if (!mdok || !isfinite(seeing)) { - psWarning("No SEEING supplied for image %d --- set to unity.", fileNum); - seeing = 1.0; - } float weighting = psMetadataLookupF32(&mdok, inputFile->fpa->analysis, "PPSTACK.WEIGHTING"); // Relative weighting @@ -69,4 +60,50 @@ scale = 1.0; } + + float exposure = psMetadataLookupF32(NULL, ro->parent->concepts, "CELL.EXPOSURE"); // Exposure time +#if 0 + if (!isfinite(exposure)) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "CELL.EXPOSURE is not set for input file %ld", stack->n); + psFree(stats); + psFree(rng); + psFree(fileIter); + psFree(fpaList); + psFree(cellList); + psFree(outRO); + psFree(stack); + return false; + } +#endif + totExposure += exposure; // Total exposure time + // Generate convolved version of input + pmReadout *convolved = pmReadoutAlloc(NULL); // Convolved version of input readout + if (!ppStackMatch(convolved, ro, config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to match image %d.", fileNum); + psFree(stats); + psFree(rng); + psFree(fileIter); + psFree(fpaList); + psFree(cellList); + psFree(stack); + psFree(outRO); + psFree(convolved); + return false; + } + +#ifdef REJECTION_FILES + { + psString name = NULL; // Name of image + psStringAppend(&name, "convolved%03d.fits", fileNum); + psFits *fits = psFitsOpen(name, "w"); + psFree(name); + psFitsWriteImage(fits, NULL, convolved->image, 0, NULL); + psFitsClose(fits); + } +#endif + + +#if 0 + // Background subtraction, scaling and normalisation is performed automatically by the image matching // Brain-dead background subtraction @@ -80,4 +117,5 @@ psFree(stack); psFree(outRO); + psFree(convolved); return false; } @@ -99,25 +137,55 @@ psFree(cellList); psFree(outRO); + psFree(convolved); psFree(stack); return false; } - totExposure += exposure; // Total exposure time + (void)psBinaryOp(ro->image, ro->image, "/", psScalarAlloc(exposure, PS_TYPE_F32)); if (ro->weight) { (void)psBinaryOp(ro->weight, ro->weight, "/", psScalarAlloc(exposure, PS_TYPE_F32)); } - pmStackData *data = pmStackDataAlloc(ro, seeing, weighting); // Data to stack +#endif + + if (fileNum == 0) { + // Copy astrometry over + pmFPA *fpa = ro->parent->parent->parent; // Template FPA + pmHDU *hdu = fpa->hdu; // Template HDU + pmHDU *outHDU = outFPA->hdu; // Output HDU + if (!outHDU || !hdu) { + psWarning("Unable to find HDU at FPA level to copy astrometry."); + } else { + if (!pmAstromReadWCS(outFPA, outCell->parent, hdu->header, 1.0)) { + psError(PS_ERR_UNKNOWN, false, "Unable to read WCS astrometry from input FPA."); + psFree(stack); + psFree(outRO); + return false; + } + if (!outHDU->header) { + outHDU->header = psMetadataAlloc(); + } + if (!pmAstromWriteWCS(outHDU->header, outFPA, outCell->parent, WCS_TOLERANCE)) { + psError(PS_ERR_UNKNOWN, false, "Unable to write WCS astrometry to output FPA."); + psFree(stack); + psFree(outRO); + return false; + } + } + } + + // Don't need the original any more! + psFree(ro); + + // Ensure there is a mask, or pmStackCombine will complain + if (!convolved->mask) { + convolved->mask = psImageAlloc(convolved->image->numCols, convolved->image->numRows, + PS_TYPE_MASK); + psImageInit(convolved->mask, 0); + } + + pmStackData *data = pmStackDataAlloc(convolved, weighting); // Data to stack + psFree(convolved); psArrayAdd(stack, ARRAY_BUFFER, data); psFree(data); // Drop reference - - // Ensure there is a mask, or pmStackCombine will complain - if (!ro->mask) { - ro->mask = psImageAlloc(ro->image->numCols, ro->image->numRows, PS_TYPE_MASK); - psImageInit(ro->mask, 0); - } - - if (fileNum == 0) { - templateRO = ro; - } fileNum++; @@ -159,12 +227,40 @@ #endif - // Only perform the additional rejection if we have seeing information - if (seeing && !pmStackReject(stack, maskBad, extent, convolveRej)) { - psError(PS_ERR_UNKNOWN, false, "Unable to reject input pixels."); - psFree(fpaList); - psFree(cellList); - psFree(stack); - psFree(outRO); - return false; + for (int i = 0; i < num; i++) { + pmStackData *data = stack->data[i]; // Data for this image + pmReadout *readout = data->readout; // Readout for this image + + // Extract the regions and solutions used in the image matching + psArray *regions = psArrayAllocEmpty(ARRAY_BUFFER); // Array of regions + { + psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD, + "^SUBTRACTION.REGION$"); // Iterator + psMetadataItem *item = NULL;// Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + assert(item->type == PS_DATA_REGION); + regions = psArrayAdd(regions, ARRAY_BUFFER, item->data.V); + } + psFree(iter); + } + psArray *solutions = psArrayAllocEmpty(ARRAY_BUFFER); // Array of solutions + { + psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD, + "^SUBTRACTION.SOLUTION$"); // Iterator + psMetadataItem *item = NULL;// Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + assert(item->type == PS_DATA_VECTOR); + solutions = psArrayAdd(solutions, ARRAY_BUFFER, item->data.V); + } + psFree(iter); + } + assert(regions->n == solutions->n); + + pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, readout->analysis, + "SUBTRACTION.KERNEL"); // Kernels + + psPixels *reject = pmSubtractionDeconvolveMask(data->pixels, threshold, regions, + solutions, kernels); // List of pixels to reject + psFree(data->pixels); + data->pixels = reject; } @@ -195,4 +291,5 @@ } +#if 0 // Restore image to counts using the total exposure time (void)psBinaryOp(outRO->image, outRO->image, "*", psScalarAlloc(totExposure, PS_TYPE_F32)); @@ -200,4 +297,5 @@ (void)psBinaryOp(outRO->weight, outRO->weight, "*", psScalarAlloc(totExposure, PS_TYPE_F32)); } +#endif psMetadataAddF32(outCell->concepts, PS_LIST_TAIL, "CELL.EXPOSURE", PS_META_REPLACE, "Summed exposure time (sec)", totExposure); @@ -215,28 +313,4 @@ psFree(cellList); - // Copy astrometry over - pmFPA *templateFPA = templateRO->parent->parent->parent; // Template FPA - pmHDU *templateHDU = templateFPA->hdu; // Template HDU - pmHDU *outHDU = outFPA->hdu; // Output HDU - if (!outHDU || !templateHDU) { - psWarning("Unable to find HDU at FPA level to copy astrometry."); - } else { - if (!pmAstromReadWCS(outFPA, outCell->parent, templateHDU->header, 1.0)) { - psError(PS_ERR_UNKNOWN, false, "Unable to read WCS astrometry from input FPA."); - psFree(stack); - psFree(outRO); - return false; - } - if (!outHDU->header) { - outHDU->header = psMetadataAlloc(); - } - if (!pmAstromWriteWCS(outHDU->header, outFPA, outCell->parent, WCS_TOLERANCE)) { - psError(PS_ERR_UNKNOWN, false, "Unable to write WCS astrometry to output FPA."); - psFree(stack); - psFree(outRO); - return false; - } - } - psFree(stack); psFree(outRO); // Drop reference