Index: branches/pap/ppSub/src/ppSub.h =================================================================== --- branches/pap/ppSub/src/ppSub.h (revision 23704) +++ branches/pap/ppSub/src/ppSub.h (revision 23711) @@ -24,11 +24,12 @@ // Output files, for activation/deactivation typedef enum { - PPSUB_FILES_INPUT = 0x01, // Input files - PPSUB_FILES_CONV = 0x02, // Convolved files (output) - PPSUB_FILES_SUB = 0x04, // Subtracted files (output) - PPSUB_FILES_INV = 0x08, // Inverse subtracted files (output) - PPSUB_FILES_PHOT_SUB = 0x10, // Photometry files (output) - PPSUB_FILES_PHOT_INV = 0x20, // Photometry files ( - PPSUB_FILES_ALL = 0xFF, // All files + PPSUB_FILES_INPUT = 0x01, // Input files + PPSUB_FILES_CONV = 0x02, // Convolved files (output) + PPSUB_FILES_SUB = 0x04, // Subtracted files (output) + PPSUB_FILES_INV = 0x08, // Inverse subtracted files (output) + PPSUB_FILES_PHOT = 0x10, // psphot photometry files + PPSUB_FILES_PHOT_SUB = 0x20, // Subtraction photometry files (output) + PPSUB_FILES_PHOT_INV = 0x40, // Inverse subtraction photometry files (output) + PPSUB_FILES_ALL = 0xFF, // All files } ppSubFiles; Index: branches/pap/ppSub/src/ppSubCamera.c =================================================================== --- branches/pap/ppSub/src/ppSubCamera.c (revision 23704) +++ branches/pap/ppSub/src/ppSubCamera.c (revision 23711) @@ -167,5 +167,5 @@ return false; } - data->invert = psMetadataLookupBool(NULL, recipe, "INVERSE"); + data->inverse = psMetadataLookupBool(NULL, recipe, "INVERSE"); data->photometry = psMetadataLookupBool(NULL, recipe, "PHOTOMETRY"); @@ -192,21 +192,21 @@ if (data->inverse) { // Inverse output image - pmFPAfile *outinv = defineOutputFile(config, input, true, "PPSUB.OUTINV", PM_FPA_FILE_IMAGE); - pmFPAfile *outinvMask = defineOutputFile(config, outinv, false, "PPSUB.OUTINV.MASK", - PM_FPA_FILE_MASK); - if (!outinv || !outinvMask) { + pmFPAfile *inverse = defineOutputFile(config, input, true, "PPSUB.INV", PM_FPA_FILE_IMAGE); + pmFPAfile *invMask = defineOutputFile(config, inverse, false, "PPSUB.INVERSE.MASK", + PM_FPA_FILE_MASK); + if (!inv || !invMask) { psError(PS_ERR_UNKNOWN, false, "Unable to define output files"); return false; } - outinv->save = true; - outinvMask->save = true; + inverse->save = true; + invMask->save = true; if (inVar && refVar) { - pmFPAfile *outinvVar = defineOutputFile(config, outinv, false, "PPSUB.OUTINV.VARIANCE", - PM_FPA_FILE_VARIANCE); - if (!outinvVar) { + pmFPAfile *invVar = defineOutputFile(config, inverse, false, "PPSUB.INVERSE.VARIANCE", + PM_FPA_FILE_VARIANCE); + if (!invVar) { psError(PS_ERR_UNKNOWN, false, "Unable to define output files"); return false; } - outinvVar->save = true; + invVar->save = true; } } @@ -312,4 +312,26 @@ return false; } + + // Deactivate psphot output sources --- we want to define output source files of our own + pmFPAfile *psphotOutput = pmFPAfileSelectSingle(config->files, "PSPHOT.OUTPUT", 0); + psphotOutput->save = false; + + pmFPAfile *outSources = defineOutputFile(config, output, false, "PPSUB.OUTPUT.SOURCES", + PM_FPA_FILE_CMF); + if (!outSources) { + psError(PS_ERR_UNKNOWN, false, "Unable to set up output source file."); + return false; + } + outSources->save = true; + + if (data->inverse) { + pmFPAfile *invSources = defineOutputFile(config, inverse, false, "PPSUB.OUTPUT.SOURCES", + PM_FPA_FILE_CMF); + if (!invSources) { + psError(PS_ERR_UNKNOWN, false, "Unable to set up inverse source file."); + return false; + } + } + invSources->save = true; } Index: branches/pap/ppSub/src/ppSubDefineOutput.c =================================================================== --- branches/pap/ppSub/src/ppSubDefineOutput.c (revision 23704) +++ branches/pap/ppSub/src/ppSubDefineOutput.c (revision 23711) @@ -29,9 +29,10 @@ pmCell *outCell = pmFPAfileThisCell(config->files, view, name); // Output cell pmFPA *outFPA = outCell->parent->parent; // Output FPA - pmHDU *outHDU = outFPA->hdu; // Output HDU + pmHDU *outHDU = outFPA->hdu; // Output HDU if (!outHDU->header) { outHDU->header = psMetadataAlloc(); } + // The output readout may already be present if we read in the convolution kernel pmReadout *outRO = NULL; // Output readout if (outCell->readouts && outCell->readouts->n > 0 && outCell->readouts->data[0]) { @@ -41,5 +42,5 @@ } - // convolved input images + // Convolved input images pmReadout *inConv = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); // Input readout pmReadout *refConv = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); // Reference readout @@ -63,4 +64,14 @@ kernels->description); + // Add additional data to the header + pmFPAfile *refFile = psMetadataLookupPtr(NULL, config->files, "PPSUB.REF"); // Reference file + pmFPAfile *inFile = psMetadataLookupPtr(NULL, config->files, "PPSUB.INPUT"); // Input file + psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.REFERENCE", 0, + "Subtraction reference", refFile->filename); + psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.INPUT", 0, + "Subtraction input", inFile->filename); + ppSubVersionHeader(outHDU->header); + + outRO->analysis = psMetadataCopy(outRO->analysis, analysis); Index: branches/pap/ppSub/src/ppSubFiles.c =================================================================== --- branches/pap/ppSub/src/ppSubFiles.c (revision 23704) +++ branches/pap/ppSub/src/ppSubFiles.c (revision 23711) @@ -26,8 +26,14 @@ NULL }; +// Subtraction photometry +static const char *subPhotFiles[] = { "PPSUB.OUTPUT.SOURCES", NULL }; + // Inverse subtraction files static const char *invFiles[] = { "PPSUB.OUTINV", "PPSUB.OUTINV.MASK", "PPSUB.OUTINV.VARIANCE", "PPSUB.OUTINV.KERNELS", "PPSUB.OUTINV.JPEG1", "PPSUB.OUTINV.JPEG2", NULL }; + +// Inverse subtraction photometry +static const char *subPhotFiles[] = { "PPSUB.OUTINV.SOURCES", NULL }; // Calculation (may be either input or output) files @@ -36,5 +42,6 @@ // Activate/deactivate a list of files -static void filesActivate(const char **files, // List of files +static void filesActivate(pmConfig *config, // Configuration + const char **files, // List of files bool state // Activation status to set ) @@ -49,5 +56,6 @@ // This is so we can activate/deactivate the 'calculation' files, which may be either input or output, which // is indicated by their 'save' boolean. -static void filesActivateSave(const char **files, // List of files +static void filesActivateSave(pmConfig *config, // Configuration + const char **files, // List of files bool save, // Activate when this save state is set bool state // Activation status to set @@ -69,18 +77,24 @@ if (files & PPSUB_FILES_INPUT) { filesActivate(inputFiles, state); - filesActivateSave(calcFiles, false, state); + filesActivateSave(config, calcFiles, false, state); } if (files & PPSUB_FILES_CONV) { - filesActivate(convFiles, state); + filesActivate(config, convFiles, state); } if (files & PPSUB_FILES_SUB) { - filesActivate(subFiles, state); - filesActivateSave(calcFiles, true, state); + filesActivate(config, subFiles, state); + filesActivateSave(config, calcFiles, true, state); } if (files & PPSUB_FILES_INV) { - filesActivate(invFiles, state); + filesActivate(config, invFiles, state); + } + if (files & PPSUB_FILES_PHOT_SUB) { + filesActivate(config, subPhotFiles, state); + } + if (files & PPSUB_FILES_PHOT_INV) { + filesActivate(config, invPhotFiles, state); } if (files & PPSUB_FILES_PHOT) { - psphotFilesActivate(config, state); + psphotFilesActivate(config, false); } Index: branches/pap/ppSub/src/ppSubLoop.c =================================================================== --- branches/pap/ppSub/src/ppSubLoop.c (revision 23704) +++ branches/pap/ppSub/src/ppSubLoop.c (revision 23711) @@ -55,4 +55,7 @@ } + psMetadataAddF32(data->stats, PS_LIST_TAIL, "TIME_MATCH", 0, "Time to match PSFs", + psTimerClear("PPSUB_MATCH")); + // Close input files if (!ppSubFilesIterateUp(config, PPSUB_FILES_INPUT)) { @@ -72,11 +75,11 @@ } - if (data->inverse && !ppSubDefineOutput("PPSUB.OUTINV", config, data, view)) { - psError(PS_ERR_UNKNOWN, false, "Unable to define inverse."); + if (!data->quality && !ppSubMakePSF("PPSUB.OUTPUT", "PPSUB.INVERSE", config, data, view)) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate PSF."); return false; } - if (!data->quality && !ppSubMakePSF("PPSUB.OUTPUT", "PPSUB.OUTINV", config, data, view)) { - psError(PS_ERR_UNKNOWN, false, "Unable to generate PSF."); + if (!ppSubReadoutSubtract("PPSUB.OUTPUT", config, view)) { + psError(PS_ERR_UNKNOWN, false, "Unable to subtract images."); return false; } @@ -88,34 +91,24 @@ } - if (!ppSubReadoutSubtract("PPSUB.OUTPUT", config, view)) { - psError(PS_ERR_UNKNOWN, false, "Unable to subtract images."); - return false; - } - // Higher order background subtraction using psphot - if (!ppSubBackground("PPSUB.OUTPUT", config, view)) { + if (!ppSubBackground(config, view)) { psError(PS_ERR_UNKNOWN, false, "Unable to subtract background."); return false; } - if (!data->quality && data->!ppSubReadoutPhotometry("PPSUB.OUTPUT", config, data, view)) { + if (!data->quality && ppSubReadoutPhotometry("PPSUB.OUTPUT", PPSUB_FILES_PHOT_SUB, config, data, view)) { psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry."); return false; } - if (!ppSubReadoutUpdate("PPSUB.OUTPUT", config, data, view)) { - psError(PS_ERR_UNKNOWN, false, "Unable to update."); + // Perform statistics on the cell + if (!ppSubReadoutStats(config, data, view)) { + psError(PS_ERR_UNKNOWN, false, "Unable to collect statistics"); return false; } - // Perform statistics on the cell - if (statsFile) { - pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PPSUB.OUTPUT"); // Output file - if (!output) { - psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find file PPSUB.OUTPUT.\n"); - goto ERROR; - } - psImageMaskType maskValue = pmConfigMaskGet("MASK.VALUE", config); - ppStatsFPA(data->stats, output->fpa, view, maskValue, config); + if (!ppSubReadoutJpeg("PPSUB.OUTPUT", config, data, view)) { + psError(PS_ERR_UNKNOWN, false, "Unable to update."); + return false; } @@ -127,5 +120,10 @@ } - if (!ppSubReadoutInverse("PPSUB.OUTINV", "PPSUB.OUTPUT", config, view)) { + if (data->inverse && !ppSubDefineOutput("PPSUB.INVERSE", config, data, view)) { + psError(PS_ERR_UNKNOWN, false, "Unable to define inverse."); + return false; + } + + if (!ppSubReadoutInverse("PPSUB.INVERSE", "PPSUB.OUTPUT", config, view)) { psError(PS_ERR_UNKNOWN, false, "Unable to invert images."); return false; @@ -138,7 +136,14 @@ } - if (!data->quality && data->!ppSubReadoutPhotometry("PPSUB.OUTINV", config, data, view)) { - psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry."); - return false; + if (!data->quality && data->!ppSubReadoutPhotometry("PPSUB.INVERSE", PPSUB_FILES_PHOT_INV, + config, data, view)) { + psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry."); + return false; + } + + // Close inverse subtraction files + if (!ppSubFilesIterateUp(config, PPSUB_FILES_INV)) { + psError(PPSUB_ERR_IO, false, "Unable to close subtraction files."); + return false; } } else { @@ -149,11 +154,4 @@ } } - - - if (!ppSubReadoutUpdate("PPSUB.OUTPUT", config, data, view)) { - psError(PS_ERR_UNKNOWN, false, "Unable to update."); - return false; - } - psFree(view); Index: branches/pap/ppSub/src/ppSubReadoutJpeg.c =================================================================== --- branches/pap/ppSub/src/ppSubReadoutJpeg.c (revision 23711) +++ branches/pap/ppSub/src/ppSubReadoutJpeg.c (revision 23711) @@ -0,0 +1,51 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "ppSub.h" + +bool ppSubReadoutJpeg(pmConfig *config, ppSubData *data, const pmFPAview *view) +{ + psAssert(config, "Require configuration"); + psAssert(view, "Require view"); + + bool mdok = false; // Status of MD lookup + + // Look up recipe values + psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSUB_RECIPE); // Recipe for ppSim + psAssert(recipe, "We checked this earlier, so it should be here."); + + // Generate binned JPEGs + { + psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask for bad pixels + + int bin1 = psMetadataLookupS32(NULL, recipe, "BIN1"); // First binning level + int bin2 = psMetadataLookupS32(NULL, recipe, "BIN2"); // Second binning level + + // Target cells + pmCell *cell1 = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT.JPEG1"); // Rebinned cell once + pmCell *cell2 = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT.JPEG2"); // Rebinned cell twice + + pmReadout *ro1 = pmReadoutAlloc(cell1), *ro2 = pmReadoutAlloc(cell2); // Binned readouts + if (!pmReadoutRebin(ro1, outRO, maskBad, bin1, bin1)) { + psError(PS_ERR_UNKNOWN, false, "Unable to bin output (1st binning)"); + psFree(ro1); + psFree(ro2); + return false; + } + if (!pmReadoutRebin(ro2, ro1, 0, bin2, bin2)) { + psError(PS_ERR_UNKNOWN, false, "Unable to bin output (2nd binning)"); + psFree(ro1); + psFree(ro2); + return false; + } + psFree(ro1); + psFree(ro2); + } + + return true; +} Index: branches/pap/ppSub/src/ppSubReadoutStats.c =================================================================== --- branches/pap/ppSub/src/ppSubReadoutStats.c (revision 23711) +++ branches/pap/ppSub/src/ppSubReadoutStats.c (revision 23711) @@ -0,0 +1,35 @@ + +bool ppSubReadoutStats(pmConfig *config, ppSubData *data, const pmFPAview *view) +{ + psAssert(config, "Require configuration"); + psAssert(data, "Require data"); + psAssert(view, "Require view"); + + if (!data->statsFile) { + // Nothing to do + return true; + } + + pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PPSUB.OUTPUT"); // Output file + if (!output) { + psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find file PPSUB.OUTPUT.\n"); + goto ERROR; + } + psImageMaskType maskValue = pmConfigMaskGet("MASK.VALUE", config); + ppStatsFPA(data->stats, output->fpa, view, maskValue, config); + + // Statistics on the matching + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MODE); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_STAMPS); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_DEV_MEAN); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_DEV_RMS); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_NORM); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_BGDIFF); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MX); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MY); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MXX); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MXY); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MYY); + + return true; +} Index: branches/pap/ppSub/src/ppSubReadoutUpdate.c =================================================================== --- branches/pap/ppSub/src/ppSubReadoutUpdate.c (revision 23704) +++ branches/pap/ppSub/src/ppSubReadoutUpdate.c (revision 23711) @@ -32,36 +32,4 @@ psAssert(recipe, "We checked this earlier, so it should be here."); - pmFPAfile *outFile = psMetadataLookupPtr(&mdok, config->files, "PPSUB.OUTPUT"); // Output file - pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT"); // Output image - pmFPA *outFPA = outFile->fpa; // Output FPA - pmHDU *outHDU = outFPA->hdu; // Output HDU - - // Add additional data to the header - pmFPAfile *refFile = psMetadataLookupPtr(NULL, config->files, "PPSUB.REF"); // Reference file - pmFPAfile *inFile = psMetadataLookupPtr(NULL, config->files, "PPSUB.INPUT"); // Input file - psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.REFERENCE", 0, - "Subtraction reference", refFile->filename); - psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.INPUT", 0, - "Subtraction input", inFile->filename); - ppSubVersionHeader(outHDU->header); - - // Statistics on the matching - if (data->stats) { - psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MODE); - psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_STAMPS); - psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_DEV_MEAN); - psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_DEV_RMS); - psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_NORM); - psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_BGDIFF); - psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MX); - psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MY); - psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MXX); - psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MXY); - psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MYY); - - psMetadataAddF32(data->stats, PS_LIST_TAIL, "TIME_MATCH", 0, "Time to match PSFs", - psTimerClear("PPSUB_MATCH")); - } - // Generate binned JPEGs { @@ -92,16 +60,4 @@ } -#ifdef TESTING - // Significance image - { - psImage *sig = (psImage*)psBinaryOp(NULL, outRO->image, "*", outRO->image); - psBinaryOp(sig, sig, "/", outRO->variance); - psFits *fits = psFitsOpen("significance.fits", "w"); - psFitsWriteImage(fits, NULL, sig, 0, NULL); - psFitsClose(fits); - psFree(sig); - } -#endif - return true; }