Index: branches/tap_branches/ppSub/src/Makefile.am =================================================================== --- branches/tap_branches/ppSub/src/Makefile.am (revision 25900) +++ branches/tap_branches/ppSub/src/Makefile.am (revision 27838) @@ -1,3 +1,3 @@ -bin_PROGRAMS = ppSub ppSubKernel +bin_PROGRAMS = ppSub ppSubKernel ppSubConvolve if HAVE_SVNVERSION @@ -34,4 +34,5 @@ ppSubData.c \ ppSubErrorCodes.c \ + ppSubExit.c \ ppSubFiles.c \ ppSubLoop.c \ @@ -53,4 +54,12 @@ ppSubKernel_SOURCES = \ ppSubKernel.c + +ppSubConvolve_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PPSTATS_CFLAGS) $(PSPHOT_CFLAGS) $(PPSUB_CFLAGS) +ppSubConvolve_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) $(PPSTATS_LIBS) $(PSPHOT_LIBS) $(PPSUB_LIBS) + +ppSubConvolve_SOURCES = \ + ppSubConvolve.c \ + ppSubExit.c \ + ppSubVersion.c noinst_HEADERS = \ Index: branches/tap_branches/ppSub/src/ppSub.c =================================================================== --- branches/tap_branches/ppSub/src/ppSub.c (revision 25900) +++ branches/tap_branches/ppSub/src/ppSub.c (revision 27838) @@ -24,13 +24,14 @@ int main(int argc, char *argv[]) { - psExit exitValue = PS_EXIT_SUCCESS; // Exit value + psLibInit(NULL); psTimerStart("ppSub"); - psLibInit(NULL); + + pmErrorRegister(); + ppSubErrorRegister(); + psphotErrorRegister(); ppSubData *data = NULL; // Processing data pmConfig *config = pmConfigRead(&argc, argv, PPSUB_RECIPE); // Configuration if (!config) { - psErrorStackPrint(stderr, "Error reading configuration."); - exitValue = PS_EXIT_CONFIG_ERROR; goto die; } @@ -39,6 +40,5 @@ if (!pmModelClassInit()) { - psErrorStackPrint(stderr, "Error initialising model classes.\n"); - exitValue = PS_EXIT_PROG_ERROR; + psError(PPSUB_ERR_PROG, false, "Unable to initialise model classes."); psFree(config); goto die; @@ -46,6 +46,5 @@ if (!psphotInit()) { - psErrorStackPrint(stderr, "Error initialising psphot.\n"); - exitValue = PS_EXIT_PROG_ERROR; + psError(PPSUB_ERR_PROG, false, "Error initialising psphot."); psFree(config); goto die; @@ -55,39 +54,63 @@ if (!ppSubArguments(argc, argv, data)) { - psErrorStackPrint(stderr, "Error reading arguments.\n"); - exitValue = PS_EXIT_CONFIG_ERROR; + psError(psErrorCodeLast(), false, "Error reading arguments."); goto die; } if (!ppSubCamera(data)) { - psErrorStackPrint(stderr, "Error setting up camera.\n"); - exitValue = PS_EXIT_CONFIG_ERROR; goto die; } if (!ppSubLoop(data)) { - psErrorStackPrint(stderr, "Error performing subtraction.\n"); - exitValue = PS_EXIT_SYS_ERROR; goto die; } die: - psTrace("ppSub", 1, "Finished at %f sec\n", psTimerMark("ppSub")); - psTimerStop(); + { + psExit exitValue = ppSubExitCode(PS_EXIT_SUCCESS); // Exit code - psString dump_file = psMetadataLookupStr(NULL, data->config->arguments, "-dumpconfig"); - if (dump_file) { - if (!pmConfigDump(data->config, dump_file)) { - psErrorStackPrint(stderr, "Unable to dump configuration.\n"); - exitValue = PS_EXIT_SYS_ERROR; + if (data && data->stats && data->statsFile) { + psString stats = psMetadataConfigFormat(data->stats); // Statistics to output + if (!stats || strlen(stats) == 0) { + psError(PPSUB_ERR_IO, false, "Unable to format statistics file"); + } else if (fprintf(data->statsFile, "%s", stats) != strlen(stats)) { + psError(PPSUB_ERR_IO, true, "Unable to write statistics file"); + } + psFree(stats); + if (fclose(data->statsFile) == EOF) { + psError(PPSUB_ERR_IO, true, "Unable to close statistics file"); + } + data->statsFile = NULL; + pmConfigRunFilenameAddWrite(data->config, "STATS", data->statsName); + exitValue = ppSubExitCode(exitValue); } + + if (config && !ppSubFilesIterateUp(config, PPSUB_FILES_ALL)) { + psError(psErrorCodeLast(), false, "Unable to close files."); + exitValue = ppSubExitCode(exitValue); + pmFPAfileFreeSetStrict(false); + } + + if (data) { + psString dump_file = psMetadataLookupStr(NULL, data->config->arguments, "-dumpconfig"); + if (dump_file) { + if (!pmConfigDump(data->config, dump_file)) { + psError(PPSUB_ERR_IO, false, "Unable to dump configuration."); + exitValue = ppSubExitCode(exitValue); + } + } + psFree(data); + } + + psTrace("ppSub", 1, "Finished at %f sec\n", psTimerMark("ppSub")); + psTimerStop(); + + pmVisualClose(); //close plot windows, if -visual is set + pmModelClassCleanup(); + pmConfigDone(); + psLibFinalize(); + + exitValue = ppSubExitCode(exitValue); + exit(exitValue); } - psFree(data); - - pmVisualClose(); //close plot windows, if -visual is set - pmModelClassCleanup(); - pmConfigDone(); - psLibFinalize(); - - exit(exitValue); } Index: branches/tap_branches/ppSub/src/ppSub.h =================================================================== --- branches/tap_branches/ppSub/src/ppSub.h (revision 25900) +++ branches/tap_branches/ppSub/src/ppSub.h (revision 27838) @@ -45,4 +45,6 @@ bool photometry; // Perform photometry? bool inverse; // Output inverse subtraction as well? + bool saveInConv; // Save convolved input? + bool saveRefConv; // Save convolved reference? psString stamps; // Stamps file pmPSF *psf; // Point Spread Function @@ -152,4 +154,7 @@ ); +/// Generate JPEG images +bool ppSubResidualSampleJpeg(pmConfig *config); + /// Generate inverse subtraction bool ppSubReadoutInverse(pmConfig *config // Configuration @@ -160,4 +165,10 @@ bool psMetadataCopySingle(psMetadata *target, psMetadata *source, const char *name); +bool ppSubCopyPSF (pmFPAfile *output, pmFPAfile *input, pmFPAview *view); + +/// Return appropriate exit code +psExit ppSubExitCode(psExit exitValue // Current exit value + ); + ///@} #endif Index: branches/tap_branches/ppSub/src/ppSubArguments.c =================================================================== --- branches/tap_branches/ppSub/src/ppSubArguments.c (revision 25900) +++ branches/tap_branches/ppSub/src/ppSubArguments.c (revision 27838) @@ -86,5 +86,8 @@ psMetadataAddS32(arguments, PS_LIST_TAIL, "-convolve", 0, "Image to convolve [1 or 2]", 0); psMetadataAddBool(arguments, PS_LIST_TAIL, "-photometry", 0, "Perform photometry?", NULL); - psMetadataAddBool(arguments, PS_LIST_TAIL, "-inverse", 0, "Generate inverse subtractions?", NULL); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp", 0, "Zero point for photometry", NAN); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-inverse", 0, "Generate inverse subtractions?", false); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-save-inconv", 0, "Save input convolved images?", false); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-save-refconv", 0, "Save reference convolved images?", false); psMetadataAddBool(arguments, PS_LIST_TAIL, "-visual", 0, "Show diagnostic plots", NULL); @@ -140,4 +143,8 @@ if (data->statsName && strlen(data->statsName) > 0) { psString resolved = pmConfigConvertFilename(data->statsName, config, true, true); // Resolved filename + if (!resolved) { + psError(psErrorCodeLast(), false, "Unable to resolve statistics file: %s", data->statsName); + return false; + } data->statsFile = fopen(resolved, "w"); if (!data->statsFile) { @@ -149,4 +156,7 @@ } + data->saveInConv = psMetadataLookupBool(NULL, arguments, "-save-inconv"); + data->saveRefConv = psMetadataLookupBool(NULL, arguments, "-save-refconv"); + if (psMetadataLookupBool(NULL, arguments, "-visual")) { pmVisualSetVisual(true); @@ -156,5 +166,5 @@ if (threads > 0) { if (!psThreadPoolInit(threads)) { - psError(PS_ERR_UNKNOWN, false, "Unable to setup %d threads", threads); + psError(psErrorCodeLast(), false, "Unable to setup %d threads", threads); return false; } Index: branches/tap_branches/ppSub/src/ppSubBackground.c =================================================================== --- branches/tap_branches/ppSub/src/ppSubBackground.c (revision 25900) +++ branches/tap_branches/ppSub/src/ppSubBackground.c (revision 27838) @@ -42,6 +42,6 @@ if (!modelRO) { // Create the background model - if (!psphotModelBackground(config, view, "PPSUB.OUTPUT")) { - psError(PS_ERR_UNKNOWN, false, "Unable to model background"); + if (!psphotModelBackgroundReadoutFileIndex(config, view, "PPSUB.OUTPUT", 0)) { + psError(psErrorCodeLast(), false, "Unable to model background"); psFree(view); return false; @@ -50,5 +50,5 @@ modelRO = pmFPAfileThisReadout(config->files, view, "PSPHOT.BACKMDL"); if (!modelRO) { - psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find background model"); + psError(psErrorCodeLast(), false, "Unable to find background model"); psFree(view); return false; Index: branches/tap_branches/ppSub/src/ppSubCamera.c =================================================================== --- branches/tap_branches/ppSub/src/ppSubCamera.c (revision 25900) +++ branches/tap_branches/ppSub/src/ppSubCamera.c (revision 27838) @@ -23,5 +23,6 @@ // Define an input file -static pmFPAfile *defineInputFile(pmConfig *config,// Configuration +static pmFPAfile *defineInputFile(bool *success, + pmConfig *config,// Configuration pmFPAfile *bind, // File to which to bind, or NULL char *filerule, // Name of file rule @@ -32,32 +33,40 @@ bool status; - // look for the file on the RUN metadata - pmFPAfile *file = pmFPAfileDefineFromRun(&status, bind, config, filerule); // File to return + *success = false; + + pmFPAfile *file = NULL; + + // look for the file on the argument list + if (bind) { + file = pmFPAfileBindFromArgs(&status, bind, config, filerule, argname); + } else { + file = pmFPAfileDefineFromArgs(&status, config, filerule, argname); + } + if (!status) { - psError(PS_ERR_UNKNOWN, false, "Failed to load file definition for %s", filerule); - return NULL; + psError(psErrorCodeLast(), false, "Failed to load file definition for %s", filerule); + return false; } if (!file) { - // look for the file on the argument list - if (bind) { - file = pmFPAfileBindFromArgs(&status, bind, config, filerule, argname); - } else { - file = pmFPAfileDefineFromArgs(&status, config, filerule, argname); - } + // look for the file on the RUN metadata + file = pmFPAfileDefineFromRun(&status, bind, config, filerule); // File to return if (!status) { - psError(PS_ERR_UNKNOWN, false, "Failed to load file definition for %s", filerule); - return false; + psError(psErrorCodeLast(), false, "Failed to load file definition for %s", filerule); + return NULL; } } if (!file) { + // no file defined + *success = true; return NULL; } if (file->type != fileType) { - psError(PS_ERR_UNKNOWN, true, "%s is not of type %s", filerule, pmFPAfileStringFromType(fileType)); - return NULL; - } - + psError(PPSUB_ERR_CONFIG, true, "%s is not of type %s", filerule, pmFPAfileStringFromType(fileType)); + return NULL; + } + + *success = true; return file; } @@ -75,9 +84,9 @@ pmFPAfileDefineOutput(config, template->fpa, filerule); if (!file) { - psError(PS_ERR_IO, false, _("Unable to generate output file from %s"), filerule); + psError(psErrorCodeLast(), false, _("Unable to generate output file from %s"), filerule); return NULL; } if (file->type != fileType) { - psError(PS_ERR_IO, true, "%s is not of type %s", filerule, pmFPAfileStringFromType(fileType)); + psError(PPSUB_ERR_CONFIG, true, "%s is not of type %s", filerule, pmFPAfileStringFromType(fileType)); return NULL; } @@ -92,4 +101,5 @@ pmFPAfile *bind, // File to which to bind, or NULL char *filerule, // Name of file rule + const char *argname, // Argument name for file pmFPAfileType fileType // Type of file ) @@ -98,7 +108,12 @@ // look for the file on the RUN metadata - pmFPAfile *file = pmFPAfileDefineFromRun(&status, bind, config, filerule); // File to return + pmFPAfile *file = pmFPAfileDefineFromRun(&status, NULL, config, filerule); // File to return if (!status) { - psError(PS_ERR_UNKNOWN, false, "Failed to load file definition for %s", filerule); + psError(psErrorCodeLast(), false, "Failed to load file definition for %s", filerule); + return NULL; + } + file = pmFPAfileBindFromArgs(&status, bind, config, filerule, argname); + if (!status) { + psError(psErrorCodeLast(), false, "Failed to load file definition for %s", filerule); return NULL; } @@ -112,5 +127,5 @@ file = pmFPAfileDefineOutput(config, bind ? bind->fpa : NULL, filerule); if (!status) { - psError(PS_ERR_UNKNOWN, false, "Failed to load file definition for %s", filerule); + psError(psErrorCodeLast(), false, "Failed to load file definition for %s", filerule); return false; } @@ -126,5 +141,5 @@ if (file->type != fileType) { - psError(PS_ERR_UNKNOWN, true, "%s is not of type %s", filerule, pmFPAfileStringFromType(fileType)); + psError(PPSUB_ERR_CONFIG, true, "%s is not of type %s", filerule, pmFPAfileStringFromType(fileType)); return NULL; } @@ -136,4 +151,6 @@ bool ppSubCamera(ppSubData *data) { + bool success = true; + psAssert(data, "Require processing data"); pmConfig *config = data->config; @@ -141,30 +158,57 @@ // Input image - pmFPAfile *input = defineInputFile(config, NULL, "PPSUB.INPUT", "INPUT", PM_FPA_FILE_IMAGE); - if (!input) { - psError(PS_ERR_IO, false, "Failed to build FPA from PPSUB.INPUT"); - return false; - } - defineInputFile(config, input, "PPSUB.INPUT.MASK", "INPUT.MASK", PM_FPA_FILE_MASK); - pmFPAfile *inVar = defineInputFile(config, input, "PPSUB.INPUT.VARIANCE", "INPUT.VARIANCE", - PM_FPA_FILE_VARIANCE); - defineInputFile(config, NULL, "PPSUB.INPUT.SOURCES", "INPUT.SOURCES", PM_FPA_FILE_CMF); + pmFPAfile *input = defineInputFile(&success, config, NULL, "PPSUB.INPUT", "INPUT", PM_FPA_FILE_IMAGE); + if (!success) { + psError(psErrorCodeLast(), false, "Failed to build FPA from PPSUB.INPUT"); + return false; + } + + defineInputFile(&success, config, input, "PPSUB.INPUT.MASK", "INPUT.MASK", PM_FPA_FILE_MASK); + if (!success) { + psError(psErrorCodeLast(), false, "Failed to build FPA from PPSUB.INPUT.MASK"); + return false; + } + + pmFPAfile *inVar = defineInputFile(&success, config, input, "PPSUB.INPUT.VARIANCE", "INPUT.VARIANCE", PM_FPA_FILE_VARIANCE); + if (!success) { + psError(psErrorCodeLast(), false, "Failed to build FPA from PPSUB.INPUT.VARIANCE"); + return false; + } + + defineInputFile(&success, config, NULL, "PPSUB.INPUT.SOURCES", "INPUT.SOURCES", PM_FPA_FILE_CMF); + if (!success) { + psError(psErrorCodeLast(), false, "Failed to build FPA from PPSUB.INPUT.SOURCES"); + return false; + } // Reference image - pmFPAfile *ref = defineInputFile(config, NULL, "PPSUB.REF", "REF", PM_FPA_FILE_IMAGE); - if (!ref) { - psError(PS_ERR_IO, false, "Failed to build FPA from PPSUB.REF"); - return false; - } - defineInputFile(config, ref, "PPSUB.REF.MASK", "REF.MASK", PM_FPA_FILE_MASK); - pmFPAfile *refVar = defineInputFile(config, ref, "PPSUB.REF.VARIANCE", "REF.VARIANCE", - PM_FPA_FILE_VARIANCE); - defineInputFile(config, NULL, "PPSUB.REF.SOURCES", "REF.SOURCES", PM_FPA_FILE_CMF); - + pmFPAfile *ref = defineInputFile(&success, config, NULL, "PPSUB.REF", "REF", PM_FPA_FILE_IMAGE); + if (!success) { + psError(psErrorCodeLast(), false, "Failed to build FPA from PPSUB.REF"); + return false; + } + + defineInputFile(&success, config, ref, "PPSUB.REF.MASK", "REF.MASK", PM_FPA_FILE_MASK); + if (!success) { + psError(psErrorCodeLast(), false, "Failed to build FPA from PPSUB.REF.MASK"); + return false; + } + + pmFPAfile *refVar = defineInputFile(&success, config, ref, "PPSUB.REF.VARIANCE", "REF.VARIANCE", PM_FPA_FILE_VARIANCE); + if (!success) { + psError(psErrorCodeLast(), false, "Failed to build FPA from PPSUB.REF.VARIANCE"); + return false; + } + + defineInputFile(&success, config, NULL, "PPSUB.REF.SOURCES", "REF.SOURCES", PM_FPA_FILE_CMF); + if (!success) { + psError(psErrorCodeLast(), false, "Failed to build FPA from PPSUB.REF.SOURCES"); + return false; + } // Now that the camera has been determined, we can read the recipe psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim if (!recipe) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find recipe %s", PPSUB_RECIPE); + psError(PPSUB_ERR_CONFIG, false, "Unable to find recipe %s", PPSUB_RECIPE); return false; } @@ -180,7 +224,4 @@ data->inverse = psMetadataLookupBool(NULL, recipe, "INVERSE"); data->photometry = psMetadataLookupBool(NULL, recipe, "PHOTOMETRY"); - - bool mdok; // Status of MD lookup - bool saveConv = psMetadataLookupBool(&mdok, recipe, "SAVE.CONVOLVED"); // Save convolved images? // Convolved input image @@ -189,21 +230,17 @@ PM_FPA_FILE_MASK); if (!inConvImage || !inConvMask) { - psError(PS_ERR_UNKNOWN, false, "Unable to define output files"); - return false; - } - if (saveConv) { - inConvImage->save = true; - inConvMask->save = true; - } + psError(psErrorCodeLast(), false, "Unable to define output files"); + return false; + } + inConvImage->save = data->saveInConv; + inConvMask->save = data->saveInConv; if (inVar) { pmFPAfile *inConvVar = defineOutputFile(config, inConvImage, false, "PPSUB.INPUT.CONV.VARIANCE", PM_FPA_FILE_VARIANCE); if (!inConvVar) { - psError(PS_ERR_UNKNOWN, false, "Unable to define output files"); - return false; - } - if (saveConv) { - inConvVar->save = true; - } + psError(psErrorCodeLast(), false, "Unable to define output files"); + return false; + } + inConvVar->save = data->saveInConv; } @@ -213,21 +250,17 @@ PM_FPA_FILE_MASK); if (!refConvImage || !refConvMask) { - psError(PS_ERR_UNKNOWN, false, "Unable to define output files"); - return false; - } - if (saveConv) { - refConvImage->save = true; - refConvMask->save = true; - } + psError(psErrorCodeLast(), false, "Unable to define output files"); + return false; + } + refConvImage->save = data->saveRefConv; + refConvMask->save = data->saveRefConv; if (refVar) { pmFPAfile *refConvVar = defineOutputFile(config, refConvImage, false, "PPSUB.REF.CONV.VARIANCE", PM_FPA_FILE_VARIANCE); if (!refConvVar) { - psError(PS_ERR_UNKNOWN, false, "Unable to define output files"); - return false; - } - if (saveConv) { - refConvVar->save = true; - } + psError(psErrorCodeLast(), false, "Unable to define output files"); + return false; + } + refConvVar->save = data->saveRefConv; } @@ -236,5 +269,5 @@ pmFPAfile *outMask = defineOutputFile(config, output, false, "PPSUB.OUTPUT.MASK", PM_FPA_FILE_MASK); if (!output || !outMask) { - psError(PS_ERR_UNKNOWN, false, "Unable to define output files"); + psError(psErrorCodeLast(), false, "Unable to define output files"); return false; } @@ -245,5 +278,5 @@ PM_FPA_FILE_VARIANCE); if (!outVar) { - psError(PS_ERR_UNKNOWN, false, "Unable to define output files"); + psError(psErrorCodeLast(), false, "Unable to define output files"); return false; } @@ -258,5 +291,5 @@ PM_FPA_FILE_MASK); if (!inverse || !invMask) { - psError(PS_ERR_UNKNOWN, false, "Unable to define output files"); + psError(psErrorCodeLast(), false, "Unable to define output files"); return false; } @@ -267,5 +300,5 @@ PM_FPA_FILE_VARIANCE); if (!invVar) { - psError(PS_ERR_UNKNOWN, false, "Unable to define output files"); + psError(psErrorCodeLast(), false, "Unable to define output files"); return false; } @@ -278,9 +311,9 @@ pmFPAfile *jpeg1 = pmFPAfileDefineOutput(config, NULL, "PPSUB.OUTPUT.JPEG1"); if (!jpeg1) { - psError(PS_ERR_IO, false, _("Unable to generate output file from PPSUB.OUTPUT.JPEG1")); + psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSUB.OUTPUT.JPEG1")); return false; } if (jpeg1->type != PM_FPA_FILE_JPEG) { - psError(PS_ERR_IO, true, "PPSUB.OUTPUT.JPEG1 is not of type JPEG"); + psError(psErrorCodeLast(), true, "PPSUB.OUTPUT.JPEG1 is not of type JPEG"); return false; } @@ -288,15 +321,28 @@ pmFPAfile *jpeg2 = pmFPAfileDefineOutput(config, NULL, "PPSUB.OUTPUT.JPEG2"); if (!jpeg2) { - psError(PS_ERR_IO, false, _("Unable to generate output file from PPSUB.OUTPUT.JPEG2")); + psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSUB.OUTPUT.JPEG2")); return false; } if (jpeg2->type != PM_FPA_FILE_JPEG) { - psError(PS_ERR_IO, true, "PPSUB.OUTPUT.JPEG2 is not of type JPEG"); + psError(psErrorCodeLast(), true, "PPSUB.OUTPUT.JPEG2 is not of type JPEG"); return false; } jpeg2->save = true; + // Output residual JPEG + pmFPAfile *jpeg3 = pmFPAfileDefineOutput(config, NULL, "PPSUB.OUTPUT.RESID.JPEG"); + if (!jpeg3) { + psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSUB.OUTPUT.RESID.JPEG")); + return false; + } + if (jpeg3->type != PM_FPA_FILE_JPEG) { + psError(psErrorCodeLast(), true, "PPSUB.OUTPUT.RESID.JPEG is not of type JPEG"); + return false; + } + jpeg3->save = true; + // Output subtraction kernel - pmFPAfile *kernel = defineCalcFile(config, output, "PPSUB.OUTPUT.KERNELS", PM_FPA_FILE_SUBKERNEL); + pmFPAfile *kernel = defineCalcFile(config, output, "PPSUB.OUTPUT.KERNELS", "KERNEL", + PM_FPA_FILE_SUBKERNEL); if (!kernel) { psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to define file PPSUB.OUTPUT.KERNELS"); @@ -305,26 +351,28 @@ // psPhot input - if (data->photometry) { + if (data->photometry || 1) { psphotModelClassInit(); // load implementation-specific models pmFPAfile *psphot = pmFPAfileDefineFromFPA(config, output->fpa, 1, 1, "PSPHOT.INPUT"); if (!psphot) { - psError(PS_ERR_IO, false, "Failed to build FPA from PSPHOT.INPUT"); + psError(psErrorCodeLast(), false, "Failed to build FPA from PSPHOT.INPUT"); return false; } if (psphot->type != PM_FPA_FILE_IMAGE) { - psError(PS_ERR_IO, true, "PSPHOT.INPUT is not of type IMAGE"); - return false; - } + psError(psErrorCodeLast(), true, "PSPHOT.INPUT is not of type IMAGE"); + return false; + } + // specify the number of psphot input images + psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "number of inputs", 1); pmFPAfileActivate(config->files, false, "PSPHOT.INPUT"); - // Internal-ish file for getting the PSF from the minuend - pmFPAfile *psf = pmFPAfileDefineOutputFromFile(config, psphot, "PSPHOT.PSF.LOAD"); + // Internal file for getting the PSF from the minuend + pmFPAfile *psf = pmFPAfileDefineFromFPA(config, output->fpa, 1, 1, "PSPHOT.PSF.LOAD"); if (!psf) { - psError(PS_ERR_IO, false, "Failed to build FPA from PSPHOT.PSF.LOAD"); + psError(psErrorCodeLast(), false, "Failed to build FPA from PSPHOT.PSF.LOAD"); return false; } if (psf->type != PM_FPA_FILE_PSF) { - psError(PS_ERR_IO, true, "PSPHOT.PSF.LOAD is not of type PSF"); + psError(psErrorCodeLast(), true, "PSPHOT.PSF.LOAD is not of type PSF"); return false; } @@ -332,5 +380,5 @@ if (!psphotDefineFiles(config, psphot)) { - psError(PS_ERR_UNKNOWN, false, "Unable to set up psphot files."); + psError(psErrorCodeLast(), false, "Unable to set up psphot files."); return false; } @@ -343,5 +391,5 @@ PM_FPA_FILE_CMF); if (!outSources) { - psError(PS_ERR_UNKNOWN, false, "Unable to set up output source file."); + psError(psErrorCodeLast(), false, "Unable to set up output source file."); return false; } @@ -352,5 +400,5 @@ PM_FPA_FILE_CMF); if (!invSources) { - psError(PS_ERR_UNKNOWN, false, "Unable to set up inverse source file."); + psError(psErrorCodeLast(), false, "Unable to set up inverse source file."); return false; } Index: branches/tap_branches/ppSub/src/ppSubConvolve.c =================================================================== --- branches/tap_branches/ppSub/src/ppSubConvolve.c (revision 27838) +++ branches/tap_branches/ppSub/src/ppSubConvolve.c (revision 27838) @@ -0,0 +1,285 @@ +#include +#include +#include + +#include "ppSub.h" + +// Print usage information and die +static void usage(const char *program, // Name of the program + psMetadata *arguments, // Command-line arguments + pmConfig *config // Configuration + ) +{ + fprintf(stderr, "\nPan-STARRS image convolution\n\n"); + fprintf(stderr, "Usage: %s OUTPUT_ROOT\n\n", + program); + fprintf(stderr, "\n"); + psArgumentHelp(arguments); + psFree(config); + pmConfigDone(); + psLibFinalize(); + exit(PS_EXIT_CONFIG_ERROR); +} + +// Put filename in it's own list +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; +} + + +int main(int argc, char *argv[]) +{ + pmFPAview *view = pmFPAviewAlloc(0); // Pointer into FPA hierarchy + + pmConfig *config = pmConfigRead(&argc, argv, PPSUB_RECIPE); // Configuration + if (!config) { + psError(psErrorCodeLast(), false, "Unable to read configuration"); + goto die; + } + + bool reference = false; // Input is actually the reference image? + int threads = 0; // Number of threads + { + psMetadata *arguments = config->arguments; // Command-line arguments + psMetadataAddStr(arguments, PS_LIST_TAIL, "-image", 0, "Input image", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-mask", 0, "Input mask", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-kernel", 0, "Convolution kernel", NULL); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-reference", 0, "Input is actually reference?", false); + psMetadataAddS32(arguments, PS_LIST_TAIL, "-threads", 0, "Threads to use", 0); + + if (argc == 1 || !psArgumentParse(arguments, &argc, argv) || argc != 2) { + usage(argv[0], arguments, config); + } + + const char *inImage = psMetadataLookupStr(NULL, arguments, "-image"); // Input image + const char *inMask = psMetadataLookupStr(NULL, arguments, "-mask"); // Input mask + const char *inKernel = psMetadataLookupStr(NULL, arguments, "-kernel"); // Input kernel + if (!inImage || !inKernel) { + usage(argv[0], arguments, config); + } + + reference = psMetadataLookupBool(NULL, arguments, "-reference"); + threads = psMetadataLookupS32(NULL, arguments, "-threads"); + fileList("PPSUB.INPUT", inImage, "Input image", config); + fileList("PPSUB.INPUT.MASK", inMask, "Input mask", config); + fileList("PPSUB.INPUT.KERNEL", inKernel, "Input kernel", config); + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "Name of the output image", argv[1]); + + bool status; // Status of file definition + pmFPAfile *image = pmFPAfileDefineFromArgs(&status, config, "PPSUB.INPUT", "PPSUB.INPUT"); + if (!image || !status || image->type != PM_FPA_FILE_IMAGE) { + psError(PPSUB_ERR_CONFIG, false, "Unable to define input image file."); + goto die; + } + pmFPAfile *mask = pmFPAfileBindFromArgs(&status, image, config, "PPSUB.INPUT.MASK", + "PPSUB.INPUT.MASK"); + if (!mask || !status || mask->type != PM_FPA_FILE_MASK) { + psError(PPSUB_ERR_CONFIG, false, "Unable to define input mask file."); + goto die; + } + + pmFPAfile *kernel = pmFPAfileBindFromArgs(&status, image, config, "PPSUB.INPUT.KERNEL", + "PPSUB.INPUT.KERNEL"); + if (!kernel || !status || kernel->type != PM_FPA_FILE_SUBKERNEL) { + psError(PPSUB_ERR_CONFIG, false, "Unable to define kernel file."); + goto die; + } + + pmFPAfile *output = pmFPAfileDefineFromFile(config, image, 1, 1, "PPSUB.INPUT.CONV"); + if (!output || output->type != PM_FPA_FILE_IMAGE) { + psError(PPSUB_ERR_CONFIG, false, "Unable to define output file."); + goto die; + } + output->save = true; + } + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim + if (!recipe) { + psError(PPSUB_ERR_CONFIG, false, "Unable to find recipe %s", PPSUB_RECIPE); + goto die; + } + + // Read everything + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "File checks failed."); + goto die; + } + view->chip = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "File checks failed."); + goto die; + } + view->cell = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "File checks failed."); + goto die; + } + view->readout = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "File checks failed."); + goto die; + } + + + pmReadout *inRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); // Input readout + pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, inRO->analysis, + PM_SUBTRACTION_ANALYSIS_KERNEL); // Convolution kernel + pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSUB.INPUT.CONV"); // Output cell + pmReadout *outRO = pmReadoutAlloc(outCell); // Output readout + pmConceptsCopyFPA(outRO->parent->parent->parent, inRO->parent->parent->parent, true, true); + pmHDU *inHDU = pmHDUFromCell(inRO->parent); // Input header + pmHDU *outHDU = pmHDUFromCell(outCell); // Output header + outHDU->header = psMetadataCopy(outHDU->header, inHDU->header); + ppSubVersionHeader(outHDU->header); + + pmReadout *input, *ref, *inConv, *refConv; // Input and output readouts + if (reference) { + ref = inRO; + refConv = outRO; + input = NULL; + inConv = NULL; + } else { + input = inRO; + inConv = outRO; + ref = NULL; + refConv = NULL; + } + + bool convolve = false; // Do we need to convolve? + float norm = 1.0; // Normalisation to apply + switch (kernels->mode) { + case PM_SUBTRACTION_MODE_1: + if (!reference) { + convolve = true; + } + norm = 1.0 / kernels->solution1->data.F64[PM_SUBTRACTION_INDEX_NORM(kernels)]; + break; + case PM_SUBTRACTION_MODE_2: + if (reference) { + convolve = true; + } + // No normalisation; should already be normalised appropriately + break; + case PM_SUBTRACTION_MODE_DUAL: { + // We only have a single input, so hack together a single-mode kernel + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index of normalisation term + norm = 1.0 / kernels->solution1->data.F64[normIndex]; + convolve = true; + if (reference) { + for (int i = 0; i < kernels->solution2->n; i++) { + kernels->solution1->data.F64[i] = -kernels->solution2->data.F64[i]; + } + kernels->solution1->data.F64[normIndex] = 1.0; + kernels->solution1->data.F64[PM_SUBTRACTION_INDEX_BG(kernels)] = 0.0; + kernels->mode = PM_SUBTRACTION_MODE_2; + } else { + kernels->mode = PM_SUBTRACTION_MODE_1; + } + psFree(kernels->solution2); + kernels->solution2 = NULL; + break; + } + default: + psAbort("Unrecognised kernel mode: %x", kernels->mode); + } + + if (convolve) { + int stride = psMetadataLookupS32(NULL, recipe, "STRIDE"); // Size of convolution patches + float kernelErr = psMetadataLookupF32(NULL, recipe, "KERNEL.ERR"); // Relative sys error in kernel + float covarFrac = psMetadataLookupF32(NULL, recipe, "COVAR.FRAC"); // Fraction for covar calculation + float badFrac = psMetadataLookupF32(NULL, recipe, "BADFRAC"); // Maximum bad fraction + float poorFrac = psMetadataLookupF32(NULL, recipe, "POOR.FRACTION"); // Fraction for "poor" + + psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config); // Bits to mask in inputs + psImageMaskType maskPoor = pmConfigMaskGet("CONV.POOR", config); // Bits to mask for poor pixels + psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask for bad pixels + + if (threads > 0) { + pmSubtractionThreadsInit(); + } + if (!pmSubtractionMatchPrecalc(inConv, refConv, input, ref, inRO->analysis, + stride, kernelErr, covarFrac, maskVal, maskBad, maskPoor, + poorFrac, badFrac)) { + psError(psErrorCodeLast(), false, "Unable to convolve images."); + if (threads > 0) { + pmSubtractionThreadsFinalize(); + } + goto die; + } + if (threads > 0) { + pmSubtractionThreadsFinalize(); + } + } else { + // Subtract background, since that's what ppSub will do with it + psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config); // Bits to mask in inputs + psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for background + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); + psStatsInit(bg); + if (!psImageBackground(bg, NULL, inRO->image, inRO->mask, maskVal, rng)) { + psError(PPSUB_ERR_DATA, false, "Unable to measure background statistics."); + psFree(bg); + psFree(rng); + goto die; + } + psFree(rng); + outRO->image = (psImage*)psBinaryOp(outRO->image, inRO->image, "-", + psScalarAlloc((float)bg->robustMedian, PS_TYPE_F32)); + psFree(bg); + outRO->data_exists = outRO->parent->data_exists = outRO->parent->parent->data_exists = true; + } + + if (norm != 1.0) { + psBinaryOp(outRO->image, outRO->image, "*", psScalarAlloc(norm, PS_TYPE_F32)); + } + + psFree(outRO); + + die: + { + psExit exitValue = ppSubExitCode(PS_EXIT_SUCCESS); // Exit value + + // Write everything + view->chip = view->cell = view->readout = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(psErrorCodeLast(), false, "File checks failed."); + exitValue = ppSubExitCode(exitValue); + pmFPAfileFreeSetStrict(false); + } + view->readout = -1; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(psErrorCodeLast(), false, "File checks failed."); + exitValue = ppSubExitCode(exitValue); + pmFPAfileFreeSetStrict(false); + } + view->cell = -1; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(psErrorCodeLast(), false, "File checks failed."); + exitValue = ppSubExitCode(exitValue); + pmFPAfileFreeSetStrict(false); + } + view->chip = -1; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(psErrorCodeLast(), false, "File checks failed."); + exitValue = ppSubExitCode(exitValue); + pmFPAfileFreeSetStrict(false); + } + psFree(view); + + psFree(config); + pmConceptsDone(); + pmConfigDone(); + psLibFinalize(); + exit(exitValue); + } + + psAbort("Should never reach here."); +} Index: branches/tap_branches/ppSub/src/ppSubData.c =================================================================== --- branches/tap_branches/ppSub/src/ppSubData.c (revision 25900) +++ branches/tap_branches/ppSub/src/ppSubData.c (revision 27838) @@ -13,15 +13,5 @@ static void subDataFree(ppSubData *data) { - if (data->statsFile) { - psString stats = psMetadataConfigFormat(data->stats); // Statistics to output - if (!stats || strlen(stats) == 0) { - psWarning("Unable to generate statistics file."); - } else { - fprintf(data->statsFile, "%s", stats); - } - psFree(stats); - fclose(data->statsFile); - pmConfigRunFilenameAddWrite(data->config, "STATS", data->statsName); - } + psAssert(!data->statsFile, "Statistics file still open."); psFree(data->statsName); psFree(data->stamps); @@ -43,4 +33,6 @@ data->photometry = false; data->inverse = false; + data->saveInConv = false; + data->saveRefConv = false; data->stamps = NULL; data->psf = NULL; Index: branches/tap_branches/ppSub/src/ppSubDefineOutput.c =================================================================== --- branches/tap_branches/ppSub/src/ppSubDefineOutput.c (revision 25900) +++ branches/tap_branches/ppSub/src/ppSubDefineOutput.c (revision 27838) @@ -51,17 +51,20 @@ bool mdok; // Status of MD lookup psMetadata *analysis = inConv->analysis; // Analysis metadata with kernel information - pmSubtractionKernels *kernels = psMetadataLookupPtr(&mdok, analysis, - PM_SUBTRACTION_ANALYSIS_KERNEL); // Subtraction kernel + pmHDU *hdu = pmHDUFromCell(inConv->parent); + pmSubtractionKernels *kernels = psMetadataLookupPtr(&mdok, analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); // Subtraction kernel + if (!kernels) { + hdu = pmHDUFromCell(refConv->parent); analysis = refConv->analysis; kernels = psMetadataLookupPtr(&mdok, analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); } if (!kernels) { - psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find SUBTRACTION.KERNEL"); + psError(PPSUB_ERR_UNKNOWN, true, "Unable to find SUBTRACTION.KERNEL"); psFree(outRO); return false; } - psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.KERNEL", 0, "Subtraction kernel", - kernels->description); + psAssert (hdu, "unable to find HDU"); + psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.KERNEL", 0, "Subtraction kernel", kernels->description); + outHDU->header = psMetadataCopy(outHDU->header, hdu->header); // Add additional data to the header Index: branches/tap_branches/ppSub/src/ppSubErrorCodes.dat =================================================================== --- branches/tap_branches/ppSub/src/ppSubErrorCodes.dat (revision 25900) +++ branches/tap_branches/ppSub/src/ppSubErrorCodes.dat (revision 27838) @@ -10,2 +10,3 @@ DATA Problem in data values NO_OVERLAP No overlap between input and skycell +PROG Programming error Index: branches/tap_branches/ppSub/src/ppSubExit.c =================================================================== --- branches/tap_branches/ppSub/src/ppSubExit.c (revision 27838) +++ branches/tap_branches/ppSub/src/ppSubExit.c (revision 27838) @@ -0,0 +1,64 @@ +#include +#include +#include + +#include "ppSub.h" + +psExit ppSubExitCode(psExit exitValue) +{ + if (exitValue != PS_EXIT_SUCCESS) { + return exitValue; + } + + psErrorCode errorCode = psErrorCodeLast(); // Error code + if (errorCode != PS_ERR_NONE) { + pmFPAfileFreeSetStrict(false); + psErrorStackPrint(stderr, "Error in subtraction:"); + switch (errorCode) { + case PPSUB_ERR_UNKNOWN: + case PS_ERR_UNKNOWN: + psLogMsg("ppSub", PS_LOG_WARN, "Unknown error code: %x", errorCode); + exitValue = PS_EXIT_UNKNOWN_ERROR; + break; + case PS_ERR_IO: + case PS_ERR_DB_CLIENT: + case PS_ERR_DB_SERVER: + case PS_ERR_BAD_FITS: + case PS_ERR_OS_CALL_FAILED: + case PM_ERR_SYS: + case PPSUB_ERR_IO: + psLogMsg("ppSub", PS_LOG_WARN, "I/O error code: %x", errorCode); + exitValue = PS_EXIT_SYS_ERROR; + break; + case PS_ERR_BAD_PARAMETER_VALUE: + case PS_ERR_BAD_PARAMETER_TYPE: + case PS_ERR_BAD_PARAMETER_NULL: + case PS_ERR_BAD_PARAMETER_SIZE: + case PPSUB_ERR_ARGUMENTS: + case PPSUB_ERR_CONFIG: + psLogMsg("ppSub", PS_LOG_WARN, "Configuration error code: %x", errorCode); + exitValue = PS_EXIT_CONFIG_ERROR; + break; + case PSPHOT_ERR_PSF: + case PPSUB_ERR_DATA: + case PPSUB_ERR_NO_OVERLAP: + psLogMsg("ppSub", PS_LOG_WARN, "Data error code: %x", errorCode); + exitValue = PS_EXIT_DATA_ERROR; + break; + case PS_ERR_UNEXPECTED_NULL: + case PS_ERR_PROGRAMMING: + case PPSUB_ERR_NOT_IMPLEMENTED: + case PPSUB_ERR_PROG: + psLogMsg("ppSub", PS_LOG_WARN, "Programming error code: %x", errorCode); + exitValue = PS_EXIT_PROG_ERROR; + break; + default: + // It's a programming error if we're not dealing with the error correctly + psLogMsg("ppSub", PS_LOG_WARN, "Unrecognised error code: %x", errorCode); + exitValue = PS_EXIT_PROG_ERROR; + break; + } + } + + return exitValue; +} Index: branches/tap_branches/ppSub/src/ppSubFiles.c =================================================================== --- branches/tap_branches/ppSub/src/ppSubFiles.c (revision 25900) +++ branches/tap_branches/ppSub/src/ppSubFiles.c (revision 27838) @@ -23,5 +23,5 @@ // Subtraction files static const char *subFiles[] = { "PPSUB.OUTPUT", "PPSUB.OUTPUT.MASK", "PPSUB.OUTPUT.VARIANCE", - "PPSUB.OUTPUT.JPEG1", "PPSUB.OUTPUT.JPEG2", + "PPSUB.OUTPUT.JPEG1", "PPSUB.OUTPUT.JPEG2", "PPSUB.OUTPUT.RESID.JPEG", NULL }; @@ -183,4 +183,41 @@ } + psFree(view); + + // Disable writing of data now that we've written things out + if (files & PPSUB_FILES_CONV) { + pmFPAview *view = ppSubViewReadout(); // View to readout + { + pmReadout *ro = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); + if (ro) { + ro->data_exists = ro->parent->data_exists = ro->parent->parent->data_exists = false; + } + } + { + pmReadout *ro = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); + if (ro) { + ro->data_exists = ro->parent->data_exists = ro->parent->parent->data_exists = false; + } + } + psFree(view); + } + if (files & PPSUB_FILES_SUB) { + pmFPAview *view = ppSubViewReadout(); // View to readout + pmReadout *ro = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT"); + if (ro) { + ro->data_exists = ro->parent->data_exists = ro->parent->parent->data_exists = false; + } + psFree(view); + } + // Disable writing of data now that we've written things out + if (files & PPSUB_FILES_INV) { + pmFPAview *view = ppSubViewReadout(); // View to readout + pmReadout *ro = pmFPAfileThisReadout(config->files, view, "PPSUB.INVERSE"); + if (ro) { + ro->data_exists = ro->parent->data_exists = ro->parent->parent->data_exists = false; + } + psFree(view); + } + ppSubFilesActivate(config, PPSUB_FILES_ALL, false); Index: branches/tap_branches/ppSub/src/ppSubLoop.c =================================================================== --- branches/tap_branches/ppSub/src/ppSubLoop.c (revision 25900) +++ branches/tap_branches/ppSub/src/ppSubLoop.c (revision 27838) @@ -27,5 +27,5 @@ pmConfigCamerasCull(config, NULL); - pmConfigRecipesCull(config, "PPSUB,PPSTATS,PSPHOT,MASKS,JPEG"); + pmConfigRecipesCull(config, "PPSUB,PPSTATS,PSPHOT,PSASTRO,MASKS,JPEG"); @@ -43,10 +43,10 @@ if (!ppSubSetMasks(config)) { - psError(PS_ERR_UNKNOWN, false, "Unable to set masks."); + psError(psErrorCodeLast(), false, "Unable to set masks."); return false; } if (!ppSubMatchPSFs(data)) { - psError(PS_ERR_UNKNOWN, false, "Unable to match PSFs."); + psError(psErrorCodeLast(), false, "Unable to match PSFs."); return false; } @@ -55,4 +55,9 @@ return true; } + // generate the residual stamp grid for visualization + if (!ppSubResidualSampleJpeg(config)) { + psError(psErrorCodeLast(), false, "Unable to update."); + return false; + } psMetadataAddF32(data->stats, PS_LIST_TAIL, "TIME_MATCH", 0, "Time to match PSFs", @@ -66,5 +71,5 @@ if (!ppSubLowThreshold(data)) { - psError(PS_ERR_UNKNOWN, false, "Unable to threshold images."); + psError(psErrorCodeLast(), false, "Unable to threshold images."); return false; } @@ -77,15 +82,25 @@ if (!ppSubDefineOutput("PPSUB.OUTPUT", config)) { - psError(PS_ERR_UNKNOWN, false, "Unable to define output."); + psError(psErrorCodeLast(), false, "Unable to define output."); return false; } if (!data->quality && !ppSubMakePSF(data)) { - psError(PS_ERR_UNKNOWN, false, "Unable to generate PSF."); - return false; + psError(psErrorCodeLast(), false, "Unable to generate PSF."); + return false; + } + + // Now we've got a PSF, blow away detections in case they're confused with real output detections + { + pmFPAview *view = ppSubViewReadout(); // View to readout + pmReadout *out = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT"); + if (psMetadataLookup(out->analysis, "PSPHOT.DETECTIONS")) { + psMetadataRemoveKey(out->analysis, "PSPHOT.DETECTIONS"); + } + psFree(view); } if (!ppSubReadoutSubtract(config)) { - psError(PS_ERR_UNKNOWN, false, "Unable to subtract images."); + psError(psErrorCodeLast(), false, "Unable to subtract images."); return false; } @@ -99,5 +114,5 @@ // Higher order background subtraction using psphot if (!ppSubBackground(config)) { - psError(PS_ERR_UNKNOWN, false, "Unable to subtract background."); + psError(psErrorCodeLast(), false, "Unable to subtract background."); return false; } @@ -105,5 +120,5 @@ // Perform Variance correction (rescale within a modest range) if (!ppSubVarianceRescale(config)) { - psError(PS_ERR_UNKNOWN, false, "Unable to rescale variance."); + psError(psErrorCodeLast(), false, "Unable to rescale variance."); return false; } @@ -115,5 +130,5 @@ if (!data->quality && !ppSubReadoutPhotometry("PPSUB.OUTPUT", data)) { - psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry."); + psError(psErrorCodeLast(), false, "Unable to perform photometry."); return false; } @@ -126,10 +141,11 @@ // Perform statistics on the cell if (!ppSubReadoutStats(data)) { - psError(PS_ERR_UNKNOWN, false, "Unable to collect statistics"); - return false; - } - + psError(psErrorCodeLast(), false, "Unable to collect statistics"); + return false; + } + + // generate the binned image used to write the jpeg if (!ppSubReadoutJpeg(config)) { - psError(PS_ERR_UNKNOWN, false, "Unable to update."); + psError(psErrorCodeLast(), false, "Unable to update."); return false; } @@ -143,10 +159,10 @@ if (data->inverse && !ppSubDefineOutput("PPSUB.INVERSE", config)) { - psError(PS_ERR_UNKNOWN, false, "Unable to define inverse."); + psError(psErrorCodeLast(), false, "Unable to define inverse."); return false; } if (!ppSubReadoutInverse(config)) { - psError(PS_ERR_UNKNOWN, false, "Unable to invert images."); + psError(psErrorCodeLast(), false, "Unable to invert images."); return false; } @@ -164,5 +180,5 @@ if (!data->quality && !ppSubReadoutPhotometry("PPSUB.INVERSE", data)) { - psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry."); + psError(psErrorCodeLast(), false, "Unable to perform photometry."); return false; } Index: branches/tap_branches/ppSub/src/ppSubMakePSF.c =================================================================== --- branches/tap_branches/ppSub/src/ppSubMakePSF.c (revision 25900) +++ branches/tap_branches/ppSub/src/ppSubMakePSF.c (revision 27838) @@ -21,4 +21,6 @@ #include "ppSub.h" + +psArray *ppSubSelectPSFSources(psArray *sources); bool ppSubMakePSF(ppSubData *data) @@ -57,17 +59,39 @@ pmFPAfile *photFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.INPUT"); // Photometry file if (!pmFPACopy(photFile->fpa, minuendFile->fpa)) { - psError(PS_ERR_UNKNOWN, false, "Unable to copy FPA for photometry"); + psError(PPSUB_ERR_CONFIG, false, "Unable to copy FPA for photometry"); psFree(view); return false; } + pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout to photometer - if (psMetadataLookup(photRO->analysis, "PSPHOT.SOURCES")) { - psMetadataRemoveKey(photRO->analysis, "PSPHOT.SOURCES"); + + // we need to remove any existing PSPHOT.DETECTIONS (why not do this in psphot?) + if (psMetadataLookup(photRO->analysis, "PSPHOT.DETECTIONS")) { + psMetadataRemoveKey(photRO->analysis, "PSPHOT.DETECTIONS"); } + if (psMetadataLookup(photRO->parent->parent->analysis, "PSPHOT.PSF")) { + psMetadataRemoveKey(photRO->parent->parent->analysis, "PSPHOT.PSF"); + } + +# ifdef TESTING + // XXX for testing, dump these images: + psphotSaveImage (NULL, photRO->image, "findpsf.im.fits"); + psphotSaveImage (NULL, photRO->variance, "findpsf.wt.fits"); + psphotSaveImage (NULL, photRO->mask, "findpsf.mk.fits"); +# endif // Extract the loaded sources from the associated readout, and generate PSF // Here, we assume the image is background-subtracted - psArray *sources = psMetadataLookupPtr(&mdok, minuend->analysis, "PSPHOT.SOURCES"); - if (!psphotReadoutFindPSF(config, view, sources)) { + pmDetections *detections = psMetadataLookupPtr(&mdok, minuend->analysis, "PSPHOT.DETECTIONS"); + if (!detections || !detections->allSources) { + psError(PPSUB_ERR_CONFIG, true, "No sources from which to determine PSF."); + return false; + } + psArray *sources = detections->allSources; + + // XXX filter sources? limit the total number and return only brighter objects? + // use flags to toss totally bogus entries? + psArray *goodSources = ppSubSelectPSFSources (sources); + if (!psphotReadoutFindPSF(config, view, goodSources)) { // This is likely a data quality issue // XXX Split into multiple cases using error codes? @@ -76,5 +100,14 @@ ppSubDataQuality(data, PSPHOT_ERR_PSF, PPSUB_FILES_PHOT_SUB | PPSUB_FILES_PHOT_INV); psFree(view); + psFree(goodSources); return true; + } + + // save the resulting PSF information on the pmFPAfile PSPHOT.PSF.LOAD + pmFPAfile *psfFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.PSF.LOAD"); // PSF file + if (!ppSubCopyPSF(psfFile, photFile, view)) { + psErrorStackPrint(stderr, "PSF was not generated"); + psWarning("PSF was not generated --- suspect bad data quality."); + ppSubDataQuality(data, psErrorCodeLast(), PPSUB_FILES_PHOT_SUB | PPSUB_FILES_PHOT_INV); } @@ -82,11 +115,70 @@ psMetadataRemoveKey(photRO->analysis, "PSPHOT.HEADER"); - if (!data->quality) { - data->psf = psMemIncrRefCounter(psMetadataLookupPtr(NULL, photRO->parent->parent->analysis, - "PSPHOT.PSF")); - } - + psFree(goodSources); psFree(view); return true; } + +bool ppSubCopyPSF(pmFPAfile *output, pmFPAfile *input, pmFPAview *view) +{ + pmChip *inputChip = pmFPAviewThisChip(view, input->fpa); // Chip with PSF info + pmChip *outputChip = pmFPAviewThisChip(view, output->fpa); // Chip to store PSF info + + pmReadout *inputRO = pmFPAviewThisReadout(view, input->fpa); // Readout with PSF info + pmReadout *outputRO = pmFPAviewThisReadout(view, output->fpa); // Readout to store PSF info + + if (!outputRO) { + pmCell *outputCell = pmFPAviewThisCell(view, output->fpa); + outputRO = pmReadoutAlloc(outputCell); + outputRO->image = psMemIncrRefCounter(inputRO->image); + } + + // Copy the PSF-related data + psMetadataIterator *iter = psMetadataIteratorAlloc(inputRO->analysis, PS_LIST_HEAD, "^PSF\\.CLUMP.*"); + psMetadataItem *item; // Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + psMetadataAddItem(outputRO->analysis, item, PS_LIST_TAIL, PS_META_REPLACE); + } + psFree(iter); + + // copy the PSF model data + pmPSF *psf = psMetadataLookupPtr(NULL, inputChip->analysis, "PSPHOT.PSF"); // PSF for photometry + if (!psf) { + psErrorStackPrint(stderr, "No PSF available"); + return false; + } + + psMetadataAddPtr(outputChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN | PS_META_REPLACE, + "PSF from ppSubMakePSF", psf); + + return true; +} + + +// XXX hardwired MAX for now +# define MAX_NPSF 500 + +psArray *ppSubSelectPSFSources(psArray *sources){ + + sources = psArraySort (sources, pmSourceSortBySN); + + psArray *subset = psArrayAllocEmpty(MAX_NPSF); + + int nPSF = 0; + for (int i = 0; (nPSF < MAX_NPSF) && (i < sources->n); i++) { + + pmSource *source = sources->data[i]; + if (!source) continue; + + // skip non-astronomical objects (very likely defects) + if (source->type == PM_SOURCE_TYPE_DEFECT) continue; + if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + if (source->mode & PM_SOURCE_MODE_SATSTAR) continue; + + psArrayAdd (subset, 100, source); + nPSF++; + } + + return subset; +} Index: branches/tap_branches/ppSub/src/ppSubMatchPSFs.c =================================================================== --- branches/tap_branches/ppSub/src/ppSubMatchPSFs.c (revision 25900) +++ branches/tap_branches/ppSub/src/ppSubMatchPSFs.c (revision 27838) @@ -18,6 +18,9 @@ #include #include +#include #include "ppSub.h" + +#define COVAR_FRAC 0.01 // Truncation fraction for covariance matrix // Normalise a region on an image @@ -37,4 +40,148 @@ } +// Measure the PSF for an image +static float subImagePSF(ppSubData *data, // Processing data + const pmReadout *ro, // Readout for which to measure PSF + psArray *sources // Sources with positions at which to measure PSF + ) +{ + psAssert(data, "Require processing data"); + pmConfig *config = data->config; // Configuration + psAssert(config, "Require configuration"); + + psAssert(ro, "Need readout"); + psAssert(sources, "Need sources."); + + pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT"); // Photometry file + psAssert(photFile, "Need photometry file."); + if (!pmFPACopy(photFile->fpa, ro->parent->parent->parent)) { + psError(PPSUB_ERR_CONFIG, false, "Unable to copy FPA for photometry"); + return false; + } + + pmFPAview *view = ppSubViewReadout(); // View to readout + pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout to photometer + + if (psMetadataLookup(photRO->analysis, "PSPHOT.DETECTIONS")) { + psMetadataRemoveKey(photRO->analysis, "PSPHOT.DETECTIONS"); + } + if (psMetadataLookup(photRO->parent->parent->analysis, "PSPHOT.PSF")) { + psMetadataRemoveKey(photRO->parent->parent->analysis, "PSPHOT.PSF"); + } + + // Extract the loaded sources from the associated readout, and generate PSF + // Here, we assume the image is background-subtracted + if (!psphotReadoutFindPSF(config, view, sources)) { + psErrorStackPrint(stderr, "Unable to determine PSF"); + psWarning("Unable to determine PSF."); + psFree(view); + return true; + } + + psFree(view); + + psMetadata *header = psMetadataLookupMetadata(NULL, photRO->analysis, "PSPHOT.HEADER"); + psAssert(header, "Require header."); + float fwhm = psMetadataLookupF32(NULL, header, "FWHM_MAJ"); + + if (psMetadataLookup(photRO->analysis, "PSPHOT.DETECTIONS")) { + psMetadataRemoveKey(photRO->analysis, "PSPHOT.DETECTIONS"); + } + if (psMetadataLookup(photRO->parent->parent->analysis, "PSPHOT.PSF")) { + psMetadataRemoveKey(photRO->parent->parent->analysis, "PSPHOT.PSF"); + } + + return fwhm; +} + +// Scale the kernel parameters according to the PSFs +static bool subScaleKernel(ppSubData *data, // Processing data + psVector *kernelWidths, // Widths for kernel + int *kernelSize, // Size of kernel + int *stampSize // Size of stamps (footprint) + ) +{ + psAssert(data, "Require processing data"); + pmConfig *config = data->config; // Configuration + psAssert(config, "Require configuration"); + + // Nothing to do if pre-calculated kernel exists + pmFPAview *view = ppSubViewReadout(); // View to readout + pmReadout *kernelRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT.KERNELS"); // RO with kernel + if (kernelRO) { + psFree(view); + return true; + } + + // Look up recipe values + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim + psAssert(recipe, "We checked this earlier, so it should be here."); + if (!psMetadataLookupBool(NULL, recipe, "SCALE")) { + // No scaling requested + psFree(view); + return true; + } + + // Input images + pmReadout *inRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); // Input readout + pmReadout *refRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF"); // Reference readout + + // Input sources + pmReadout *inSourceRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.SOURCES"); + pmReadout *refSourceRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.SOURCES"); + if (!inSourceRO || !refSourceRO) { + psWarning("Unable to scale kernel, since no sources were provided."); + return true; + } + + pmDetections *inDetections = psMetadataLookupPtr(NULL, inSourceRO->analysis, "PSPHOT.DETECTIONS"); // Input sources + pmDetections *refDetections = psMetadataLookupPtr(NULL, refSourceRO->analysis, "PSPHOT.DETECTIONS"); // Ref sources + + psFree(view); + + if (!inDetections || !refDetections) { + psWarning("Unable to scale kernel, since no sources were provided."); + return true; + } + + psArray *inSources = inDetections->allSources; + psAssert (inSources, "missing inSources?"); + + psArray *refSources = refDetections->allSources; + psAssert (refSources, "missing refSources?"); + + float inFWHM = psMetadataLookupF32(NULL, inRO->parent->parent->concepts, "CHIP.SEEING"); // FWHM for input + if (!isfinite(inFWHM) || inFWHM == 0.0) { + inFWHM = subImagePSF(data, inRO, inSources); + } + float refFWHM = psMetadataLookupF32(NULL, refRO->parent->parent->concepts, "CHIP.SEEING"); // FWHM for ref + if (!isfinite(refFWHM) || refFWHM == 0.0) { + refFWHM = subImagePSF(data, refRO, refSources); + } + psLogMsg("ppSub", PS_LOG_INFO, "Input FWHM: %f\nReference FWHM: %f\n", inFWHM, refFWHM); + if (!isfinite(inFWHM) || !isfinite(refFWHM)) { + psWarning("Unable to scale kernel, since unable to measure PSFs."); + return true; + } + + float scaleRef = psMetadataLookupF32(NULL, recipe, "SCALE.REF"); // Reference for scaling + float scaleMin = psMetadataLookupF32(NULL, recipe, "SCALE.MIN"); // Minimum for scaling + float scaleMax = psMetadataLookupF32(NULL, recipe, "SCALE.MAX"); // Maximum for scaling + if (!isfinite(scaleRef) || !isfinite(scaleMin) || !isfinite(scaleMax)) { + psError(PPSUB_ERR_ARGUMENTS, false, + "Scale parameters (SCALE.REF=%f, SCALE.MIN=%f, SCALE.MAX=%f) not set in recipe.", + scaleRef, scaleMin, scaleMax); + return false; + } + + if (!pmSubtractionParamsScale(kernelSize, stampSize, kernelWidths, inFWHM, refFWHM, + scaleRef, scaleMin, scaleMax)) { + psError(PPSUB_ERR_DATA, false, "Unable to scale parameters."); + return false; + } + + return true; +} + bool ppSubMatchPSFs(ppSubData *data) @@ -71,37 +218,44 @@ pmReadout *kernelRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT.KERNELS"); // RO with kernel - psFree(view); - // Sources in image, used for stamps: these must be loaded from previous analysis stages pmReadout *inSourceRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.SOURCES"); pmReadout *refSourceRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.SOURCES"); - psArray *inSources = inSourceRO ? psMetadataLookupPtr(&mdok, inSourceRO->analysis, "PSPHOT.SOURCES") : - NULL; // Source list from input image - psArray *refSources = refSourceRO ? psMetadataLookupPtr(&mdok, refSourceRO->analysis, "PSPHOT.SOURCES") : - NULL ; // Source list from reference image - - psArray *sources = NULL; // Merged list of sources - if (inSources && refSources) { + + pmDetections *inDetections = inSourceRO ? psMetadataLookupPtr(&mdok, inSourceRO->analysis, "PSPHOT.DETECTIONS") : NULL; // Input sources + pmDetections *refDetections = refSourceRO ? psMetadataLookupPtr(&mdok, refSourceRO->analysis, "PSPHOT.DETECTIONS") : NULL; // Ref sources + + psFree(view); + + pmDetections *detections = NULL; // Merged detection set + if (inDetections && refDetections) { + psArray *inSources = inDetections->allSources; + psArray *refSources = refDetections->allSources; + + psAssert (inSources, "missing in sources?"); + psAssert (refSources, "missing ref sources?"); + + detections = pmDetectionsAlloc(); float radius = psMetadataLookupF32(NULL, recipe, "SOURCE.RADIUS"); // Matching radius psArray *lists = psArrayAlloc(2); // Source lists lists->data[0] = psMemIncrRefCounter(inSources); lists->data[1] = psMemIncrRefCounter(refSources); - sources = pmSourceMatchMerge(lists, radius); + detections->allSources = pmSourceMatchMerge(lists, radius); psFree(lists); - if (!sources) { - psError(PS_ERR_UNKNOWN, false, "Unable to merge source lists"); + if (!detections->allSources) { + psError(PPSUB_ERR_DATA, false, "Unable to merge source lists"); + psFree(detections); return false; } - } else if (inSources) { - sources = psMemIncrRefCounter(inSources); - } else if (refSources) { - sources = psMemIncrRefCounter(refSources); - } - - psMetadataAddArray(inConv->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE, - "Merged source list", sources); - psMetadataAddArray(refConv->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE, - "Merged source list", sources); - psFree(sources); // Drop reference + } + if (!detections && inDetections) { + detections = psMemIncrRefCounter(inDetections); + } + if (!detections && refDetections) { + detections = psMemIncrRefCounter(refDetections); + } + + psMetadataAddPtr(inConv->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "Merged source list", detections); + psMetadataAddPtr(refConv->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "Merged source list", detections); + psFree(detections); // Drop reference int footprint = psMetadataLookupS32(NULL, recipe, "STAMP.FOOTPRINT"); // Stamp half-size @@ -115,5 +269,5 @@ pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Type of kernel if (type == PM_SUBTRACTION_KERNEL_NONE) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unrecognised kernel type: %s", typeStr); + psError(PPSUB_ERR_ARGUMENTS, true, "Unrecognised kernel type: %s", typeStr); return false; } @@ -130,5 +284,9 @@ int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations float rej = psMetadataLookupF32(NULL, recipe, "REJ"); // Rejection threshold - float sys = psMetadataLookupF32(NULL, recipe, "SYS"); // Relative systematic error in kernel + float kernelErr = psMetadataLookupF32(NULL, recipe, "KERNEL.ERR"); // Relative systematic error in kernel + float normFrac = psMetadataLookupF32(NULL, recipe, "NORM.FRAC"); // Fraction of window for normalisn windw + float sysErr = psMetadataLookupF32(NULL, recipe, "SYS.ERR"); // Relative systematic error in images + float skyErr = psMetadataLookupF32(NULL, recipe, "SKY.ERR"); // Additional error in sky + float covarFrac = psMetadataLookupF32(NULL, recipe, "COVAR.FRAC"); // Fraction for covariance calculation float badFrac = psMetadataLookupF32(NULL, recipe, "BADFRAC"); // Maximum bad fraction @@ -168,13 +326,29 @@ break; default: - psErrorStackPrint(stderr, "Invalid value for -convolve"); + psError(PPSUB_ERR_ARGUMENTS, false, "Invalid value for -convolve"); return false; } } + if (!subScaleKernel(data, widths, &size, &footprint)) { + psError(PPSUB_ERR_DATA, false, "Unable to scale kernel parameters"); + return false; + } + int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads if (threads > 0) { - pmSubtractionThreadsInit(inRO, refRO); - } + pmSubtractionThreadsInit(); + } + + if (inRO->covariance) { + psImageCovarianceTruncate(inRO->covariance, COVAR_FRAC); + } + if (refRO->covariance) { + psImageCovarianceTruncate(refRO->covariance, COVAR_FRAC); + } + + // Data doesn't exist until it's created in pmSubtraction... + inConv->data_exists = inConv->parent->data_exists = inConv->parent->parent->data_exists = false; + refConv->data_exists = refConv->parent->data_exists = refConv->parent->parent->data_exists = false; // Match the PSFs @@ -182,15 +356,36 @@ if (kernelRO) { success = pmSubtractionMatchPrecalc(inConv, refConv, inRO, refRO, kernelRO->analysis, - stride, sys, maskVal, maskBad, maskPoor, poorFrac, badFrac); + stride, kernelErr, covarFrac, maskVal, maskBad, maskPoor, + poorFrac, badFrac); } else { success = pmSubtractionMatch(inConv, refConv, inRO, refRO, footprint, stride, regionSize, - spacing, threshold, sources, data->stamps, type, size, order, - widths, orders, inner, ringsOrder, binning, penalty, optimum, - optWidths, optOrder, optThresh, iter, rej, sys, maskVal, - maskBad, maskPoor, poorFrac, badFrac, subMode); - } + spacing, threshold, detections ? detections->allSources : NULL, + data->stamps, type, size, order, widths, orders, inner, ringsOrder, + binning, penalty, optimum, optWidths, optOrder, optThresh, iter, rej, + normFrac, sysErr, skyErr, kernelErr, covarFrac, + maskVal, maskBad, maskPoor, poorFrac, badFrac, subMode); + } + +# ifdef TESTING + // XXX for testing + psphotSaveImage (NULL, refRO->image, "refRO.im.fits"); + psphotSaveImage (NULL, refRO->variance, "refRO.wt.fits"); + psphotSaveImage (NULL, refRO->mask, "refRO.mk.fits"); + + psphotSaveImage (NULL, inRO->image, "inRO.im.fits"); + psphotSaveImage (NULL, inRO->variance, "inRO.wt.fits"); + psphotSaveImage (NULL, inRO->mask, "inRO.mk.fits"); + + psphotSaveImage (NULL, inConv->image, "inConv.im.fits"); + psphotSaveImage (NULL, inConv->variance, "inConv.wt.fits"); + psphotSaveImage (NULL, inConv->mask, "inConv.mk.fits"); + + psphotSaveImage (NULL, refConv->image, "refConv.im.fits"); + psphotSaveImage (NULL, refConv->variance, "refConv.wt.fits"); + psphotSaveImage (NULL, refConv->mask, "refConv.mk.fits"); +# endif psFree(optWidths); - pmSubtractionThreadsFinalize(inRO, refRO); + pmSubtractionThreadsFinalize(); if (!success) { @@ -207,5 +402,5 @@ return true; } else { - psError(PS_ERR_UNKNOWN, false, "Unable to match images."); + psError(PPSUB_ERR_DATA, false, "Unable to match images."); return false; } @@ -260,4 +455,11 @@ pmConceptsCopyFPA(refConv->parent->parent->parent, refRO->parent->parent->parent, true, true); + if (inConv->covariance) { + psImageCovarianceTruncate(inConv->covariance, COVAR_FRAC); + } + if (refConv->covariance) { + psImageCovarianceTruncate(refConv->covariance, COVAR_FRAC); + } + if (inConv->variance) { psImageCovarianceTransfer(inConv->variance, inConv->covariance); Index: branches/tap_branches/ppSub/src/ppSubReadoutInverse.c =================================================================== --- branches/tap_branches/ppSub/src/ppSubReadoutInverse.c (revision 25900) +++ branches/tap_branches/ppSub/src/ppSubReadoutInverse.c (revision 27838) @@ -20,4 +20,5 @@ invRO->variance = psMemIncrRefCounter(outRO->variance); invRO->covariance = psMemIncrRefCounter(outRO->covariance); + invRO->analysis = psMetadataCopy(invRO->analysis, outRO->analysis); invRO->data_exists = invRO->parent->data_exists = invRO->parent->parent->data_exists = true; @@ -34,10 +35,10 @@ pmHDU *invHDU = invFPA->hdu; // Inverse HDU if (!pmAstromWriteWCS(invHDU->header, outFPA, outChip, WCS_TOLERANCE)) { - psError(PS_ERR_UNKNOWN, false, "Unable to write WCS astrometry to PPSUB.INVERSE."); + psError(psErrorCodeLast(), false, "Unable to write WCS astrometry to PPSUB.INVERSE."); return false; } // Read from newly written astrometry so that it exists in the "inverse" FPA (for sources) if (!pmAstromReadWCS(invFPA, invChip, invHDU->header, 1.0)) { - psError(PS_ERR_UNKNOWN, false, "Unable to read WCS astrometry."); + psError(psErrorCodeLast(), false, "Unable to read WCS astrometry."); return false; } Index: branches/tap_branches/ppSub/src/ppSubReadoutJpeg.c =================================================================== --- branches/tap_branches/ppSub/src/ppSubReadoutJpeg.c (revision 25900) +++ branches/tap_branches/ppSub/src/ppSubReadoutJpeg.c (revision 27838) @@ -33,5 +33,5 @@ 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)"); + psError(PPSUB_ERR_DATA, false, "Unable to bin output (1st binning)"); psFree(ro1); psFree(ro2); @@ -39,5 +39,5 @@ } if (!pmReadoutRebin(ro2, ro1, 0, bin2, bin2)) { - psError(PS_ERR_UNKNOWN, false, "Unable to bin output (2nd binning)"); + psError(PPSUB_ERR_DATA, false, "Unable to bin output (2nd binning)"); psFree(ro1); psFree(ro2); @@ -49,2 +49,94 @@ return true; } + +bool ppSubResidualSampleJpeg(pmConfig *config) +{ + return true; + + pmFPAview *view = ppSubViewReadout(); // View to readout + + // we save sample difference stamps on the convolved image readout->analysis metadata + pmReadout *inConv = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); // Input convolved + + psArray *samples = psArrayAllocEmpty(16); // Array of sample stamp images + + // we may have multiple entries with the same name; pull them off into a separate array: + psString regex = NULL; // Regular expression + psStringAppend(®ex, "^%s$", "SUBTRACTION.SAMPLE.STAMP.SET"); + psMetadataIterator *iter = psMetadataIteratorAlloc(inConv->analysis, PS_LIST_HEAD, regex); // Iterator + psFree(regex); + + psMetadataItem *item = NULL;// Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + assert(item->type == PS_DATA_ARRAY); + psArray *sampleStamps = item->data.V; + samples = psArrayAdd(samples, 16, sampleStamps); + } + psFree(iter); + psAssert (samples, "no sample stamps?"); + psAssert (samples->n, "no sample stamps?"); + + // get the kernel sizes + psArray *kernels = samples->data[0]; + psAssert (kernels, "no valid kernel?"); + psAssert (kernels->n, "no valid kernel?"); + + psImage *kernel = kernels->data[0]; + psAssert (kernel, "missing valid kernel?"); + + int DX = kernel->numCols; + int DY = kernel->numRows; + + // each array contains up to 9 sample stamps. generate an image mosaicking these together in 3x3 blocks. + int innerBorder = 1; + int outerBorder = 2; + int NXblock = sqrt(samples->n); + int NYblock = samples->n / NXblock; + if (samples->n % NXblock) NYblock ++; + + int NXpix = NXblock * (3 * (DX + innerBorder) + outerBorder); + int NYpix = NYblock * (3 * (DY + innerBorder) + outerBorder); + + // output cell + pmCell *cell = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT.RESID.JPEG"); + pmReadout *readout = pmReadoutAlloc(cell); + readout->image = psImageAlloc(NXpix, NYpix, PS_TYPE_F32); + + cell->data_exists = true; + cell->parent->data_exists = true; + + psImageInit (readout->image, 0.0); + + for (int i = 0; i < samples->n; i++) { + + int xBlock = i % NXblock; + int yBlock = i / NYblock; + + psArray *kernels = samples->data[i]; + + for (int j = 0; j < kernels->n; j++) { + + psImage *kernel = kernels->data[j]; + + int xSub = j % 3; + int ySub = j / 3; + + int xPix = xBlock * (3 * (DX + innerBorder) + outerBorder) + xSub * (DX + innerBorder); + int yPix = yBlock * (3 * (DX + innerBorder) + outerBorder) + ySub * (DY + innerBorder); + + for (int y = 0; y < kernel->numRows; y++) { + for (int x = 0; x < kernel->numCols; x++) { + readout->image->data.F32[y + yPix][x + xPix] = kernel->data.F32[y][x]; + } + } + } + } + + { + psFits *fits = psFitsOpen ("resid.stamps.fits", "w"); + psFitsWriteImage (fits, NULL, readout->image, 0, NULL); + psFitsClose (fits); + } + + return true; +} Index: branches/tap_branches/ppSub/src/ppSubReadoutPhotometry.c =================================================================== --- branches/tap_branches/ppSub/src/ppSubReadoutPhotometry.c (revision 25900) +++ branches/tap_branches/ppSub/src/ppSubReadoutPhotometry.c (revision 27838) @@ -24,4 +24,6 @@ bool ppSubReadoutPhotometry(const char *name, ppSubData *data) { + bool mdok = false; + psAssert(data, "Require processing data"); pmConfig *config = data->config; // Configuration @@ -42,48 +44,40 @@ } - // The PSF (measured in ppSubMakePSF) is stored on the chip->analysis of PSPHOT.INPUT - // In order to use an incoming PSF, it must be stored on the chip->analysis of PSPHOT.PSF.LOAD + // select the view of interest pmFPAview *view = ppSubViewReadout(); // View to readout - pmChip *psfInputChip = pmFPAfileThisChip(config->files, view, "PSPHOT.INPUT"); // Chip with PSF - psAssert (psfInputChip, "should have been generated for ppSubMakePSF"); - pmChip *psfLoadChip = pmFPAfileThisChip(config->files, view, "PSPHOT.PSF.LOAD"); // Chip to have PSF - psAssert (psfLoadChip, "PSPHOT.PSF.LOAD should have been defined in ppSubCamera"); - pmPSF *psf = psMetadataLookupPtr(NULL, psfInputChip->analysis, "PSPHOT.PSF"); // PSF for photometry - if (!psf) { - psErrorStackPrint(stderr, "No PSF available"); - psWarning("No PSF available --- suspect bad data quality."); - ppSubDataQuality(data, psErrorCodeLast(), PPSUB_FILES_PHOT_SUB | PPSUB_FILES_PHOT_INV); - return true; - } - psMetadataAddPtr(psfLoadChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN | PS_META_REPLACE, - "PSF from ppSubMakePSF", psf); - - bool mdok = false; // psphotReadoutMinimal performs the photometry analysis on PSPHOT.INPUT; we need to move - // around the pointers so PSPHOT.INPUT corresponds to the output image; previously, it was - // equivalent to the minuend image. - pmReadout *inRO = pmFPAfileThisReadout(config->files, view, name); // Readout with image and sources - if (psMetadataLookup(inRO->analysis, "PSPHOT.SOURCES")) { - psMetadataRemoveKey(inRO->analysis, "PSPHOT.SOURCES"); - } - + // around the pointers so PSPHOT.INPUT corresponds to the output image of interest (on one + // pass this is the subtraction image, in another it is negative of the subtraction pmFPAfile *photFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.INPUT"); // Photometry file pmFPAfile *inFile = psMetadataLookupPtr(&mdok, config->files, name); // Input file if (!pmFPACopy(photFile->fpa, inFile->fpa)) { - psError(PS_ERR_UNKNOWN, false, "Unable to copy FPA for photometry"); + psError(PPSUB_ERR_CONFIG, false, "Unable to copy FPA for photometry"); psFree(view); return false; } + + // drop references to PSPHOT.DETECTIONS on both of these (why is this needed for both??) pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout to photometer - if (psMetadataLookup(photRO->analysis, "PSPHOT.SOURCES")) { - psMetadataRemoveKey(photRO->analysis, "PSPHOT.SOURCES"); + if (psMetadataLookup(photRO->analysis, "PSPHOT.DETECTIONS")) { + psMetadataRemoveKey(photRO->analysis, "PSPHOT.DETECTIONS"); + } + pmReadout *inRO = pmFPAfileThisReadout(config->files, view, name); // Readout with image and sources + if (psMetadataLookup(inRO->analysis, "PSPHOT.DETECTIONS")) { + psMetadataRemoveKey(inRO->analysis, "PSPHOT.DETECTIONS"); } - psMetadataAddPtr(photRO->parent->parent->analysis, PS_LIST_TAIL, "PSPHOT.PSF", - PS_META_REPLACE | PS_DATA_UNKNOWN, "Point-spread function", data->psf); + // grab the PSF information from the pmFPAfile PSPHOT.PSF.LOAD + pmFPAfile *psfFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.PSF.LOAD"); // PSF file + ppSubCopyPSF (photFile, psfFile, view); - psFree(photRO->analysis); - photRO->analysis = psMetadataAlloc(); + // psphotSaveImage (photFile->fpa->hdu->header, photRO->image, "findsrc.im.fits"); + // psphotSaveImage (photFile->fpa->hdu->header, photRO->variance, "findsrc.wt.fits"); + // psphotSaveImage (photFile->fpa->hdu->header, photRO->mask, "findsrc.mk.fits"); + + // erase the overlays from a previous psphot-related step + if (pmVisualIsVisual()) { + // psphotVisualEraseOverlays (1, "all"); + } if (!psphotReadoutMinimal(config, view)) { @@ -96,8 +90,10 @@ // If no sources were found, there's no error, but we want to trigger 'bad quality' - psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources - if (!sources) { + pmDetections *detections = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.DETECTIONS"); // Sources + if (!detections) { ppSubDataQuality(data, PSPHOT_ERR_DATA, PPSUB_FILES_PHOT_SUB | PPSUB_FILES_PHOT_INV); } + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); if (data->stats) { @@ -114,11 +110,22 @@ if (!data->quality) { - if (!psMetadataCopySingle(inRO->analysis, photRO->analysis, "PSPHOT.SOURCES")) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to copy PSPHOT.SOURCES"); + if (!psMetadataCopySingle(inRO->analysis, photRO->analysis, "PSPHOT.DETECTIONS")) { + psError(PPSUB_ERR_PROG, false, "Unable to copy PSPHOT.DETECTIONS"); return false; } if (!psMetadataCopySingle(inRO->analysis, photRO->analysis, "PSPHOT.HEADER")) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to copy PSPHOT.HEADER"); + psError(PPSUB_ERR_PROG, false, "Unable to copy PSPHOT.HEADER"); return false; + } + if (!psMetadataCopySingle(inRO->analysis, photRO->analysis, PM_DETEFF_ANALYSIS)) { + psError(PPSUB_ERR_PROG, false, "Unable to copy Detection Efficiency"); + return false; + } + + // Ensure photometry information is put in the header + pmHDU *hdu = pmHDUFromReadout(inRO); // HDU for readout + if (hdu) { + psMetadata *photHeader = psMetadataLookupMetadata(NULL, inRO->analysis, "PSPHOT.HEADER"); // Header + hdu->header = psMetadataCopy(hdu->header, photHeader); } } @@ -130,13 +137,10 @@ } - - - - #ifdef TESTING // Record data about sources: not everything gets into the output CMF files { pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout with the sources - psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources + pmDetections *detections = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.DETECTIONS"); // Sources + psArray *sources = detections->allSources; FILE *sourceFile = fopen("sources.dat", "w"); // File for sources fprintf(sourceFile, Index: branches/tap_branches/ppSub/src/ppSubReadoutStats.c =================================================================== --- branches/tap_branches/ppSub/src/ppSubReadoutStats.c (revision 25900) +++ branches/tap_branches/ppSub/src/ppSubReadoutStats.c (revision 27838) @@ -23,5 +23,5 @@ 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"); + psError(PPSUB_ERR_PROG, true, "Unable to find file PPSUB.OUTPUT.\n"); return false; } Index: branches/tap_branches/ppSub/src/ppSubReadoutSubtract.c =================================================================== --- branches/tap_branches/ppSub/src/ppSubReadoutSubtract.c (revision 25900) +++ branches/tap_branches/ppSub/src/ppSubReadoutSubtract.c (revision 27838) @@ -71,5 +71,5 @@ pmFPA *outFPA = outFile->fpa; // Output FPA if (!pmConceptsCopyFPA(outFPA, inFPA, true, true)) { - psError(PS_ERR_UNKNOWN, false, "Unable to copy concepts from input to output."); + psError(PPSUB_ERR_CONFIG, false, "Unable to copy concepts from input to output."); psFree(outRO); psFree(view); @@ -84,13 +84,13 @@ psFree(view); if (!outHDU || !inHDU) { - psError(PS_ERR_UNKNOWN, false, "Unable to find HDU at FPA level to copy astrometry."); + psError(PPSUB_ERR_PROG, true, "Unable to find HDU at FPA level to copy astrometry."); return false; } if (!pmAstromReadWCS(outFPA, outChip, inHDU->header, 1.0)) { - psError(PS_ERR_UNKNOWN, false, "Unable to read WCS astrometry from input FPA."); + psError(psErrorCodeLast(), false, "Unable to read WCS astrometry from input FPA."); return false; } if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_TOLERANCE)) { - psError(PS_ERR_UNKNOWN, false, "Unable to write WCS astrometry to output FPA."); + psError(psErrorCodeLast(), false, "Unable to write WCS astrometry to output FPA."); return false; } Index: branches/tap_branches/ppSub/src/ppSubSetMasks.c =================================================================== --- branches/tap_branches/ppSub/src/ppSubSetMasks.c (revision 25900) +++ branches/tap_branches/ppSub/src/ppSubSetMasks.c (revision 27838) @@ -28,10 +28,10 @@ psImageMaskType maskValue, markValue; // Mask values if (!pmConfigMaskSetBits(&maskValue, &markValue, config)) { - psError(PS_ERR_UNKNOWN, false, "Unable to determine mask value."); + psError(PPSUB_ERR_CONFIG, false, "Unable to determine mask value."); return false; } // Set the mask bits needed by psphot (in psphot recipe) - psphotSetMaskRecipe (config, maskValue, markValue); + psphotSetMaskRecipe(config, maskValue, markValue); // Look up recipe values @@ -79,9 +79,9 @@ // Mask the NAN values if (!pmReadoutMaskNonfinite(inRO, satValue)) { - psError(PS_ERR_UNKNOWN, false, "Unable to mask non-finite pixels in input."); + psError(PPSUB_ERR_DATA, false, "Unable to mask non-finite pixels in input."); return false; } if (!pmReadoutMaskNonfinite(refRO, satValue)) { - psError(PS_ERR_UNKNOWN, false, "Unable to mask non-finite pixels in reference."); + psError(PPSUB_ERR_DATA, false, "Unable to mask non-finite pixels in reference."); return false; } @@ -94,5 +94,5 @@ psImageInterpolateMode interpMode = psImageInterpolateModeFromString(interpModeStr); // Interp if (interpMode == PS_INTERPOLATE_NONE) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unknown interpolation mode: %s", interpModeStr); + psError(PPSUB_ERR_CONFIG, false, "Unknown interpolation mode: %s", interpModeStr); return false; } @@ -104,9 +104,9 @@ // Interpolate over bad pixels, so the bad pixels don't explode if (!pmReadoutInterpolateBadPixels(inRO, maskVal, interpMode, poorFrac, maskPoor, maskBad)) { - psError(PS_ERR_UNKNOWN, false, "Unable to interpolate bad pixels for input image."); + psError(PPSUB_ERR_DATA, false, "Unable to interpolate bad pixels for input image."); return false; } if (!pmReadoutInterpolateBadPixels(refRO, maskVal, interpMode, poorFrac, maskPoor, maskBad)) { - psError(PS_ERR_UNKNOWN, false, "Unable to interpolate bad pixels for reference image."); + psError(PPSUB_ERR_DATA, false, "Unable to interpolate bad pixels for reference image."); return false; } Index: branches/tap_branches/ppSub/src/ppSubThreshold.c =================================================================== --- branches/tap_branches/ppSub/src/ppSubThreshold.c (revision 25900) +++ branches/tap_branches/ppSub/src/ppSubThreshold.c (revision 27838) @@ -27,4 +27,5 @@ psImageMaskType maskIgnore, // Ignore pixels with this mask psImageMaskType maskThresh, // Give pixels this mask if below threshold + psRegion *region, // Region of interest const char *description // Description of image ) @@ -37,10 +38,12 @@ psImage *image = ro->image; // Image psImage *mask = ro->mask; // Mask - int numCols = ro->image->numCols, numRows = ro->image->numRows; // Size of image + + psImage *subImage = psImageSubset(image, *region); // Image with region of interest + psImage *subMask = psImageSubset(mask, *region); // Maks with region of interest psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator - if (!psImageBackground(stats, NULL, image, mask, maskIgnore, rng)) { - psError(PS_ERR_UNKNOWN, false, "Unable to determine threshold."); + if (!psImageBackground(stats, NULL, subImage, subMask, maskIgnore, rng)) { + psError(PPSUB_ERR_DATA, false, "Unable to determine threshold."); psFree(rng); psFree(stats); @@ -49,10 +52,13 @@ psFree(rng); + psFree(subImage); + psFree(subMask); + float threshold = stats->robustMedian - thresh * stats->robustStdev; // Threshold below which to clip psFree(stats); psLogMsg("ppSub", PS_LOG_INFO, "Masking pixels below %f in %s", threshold, description); - for (int y = 0; y < numRows; y++) { - for (int x = 0; x < numCols; x++) { + for (int y = region->y0; y < region->y1; y++) { + for (int x = region->x0; x < region->x1; x++) { if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskIgnore) { continue; @@ -91,9 +97,5 @@ pmReadout *in = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); // Input image if (!in) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find readout."); - return false; - } - if (!lowThreshold(in, thresh, maskVal, maskThresh, "input convolved image")) { - psError(PS_ERR_UNKNOWN, false, "Unable to threshold input image."); + psError(PPSUB_ERR_UNKNOWN, false, "Unable to find readout."); return false; } @@ -101,11 +103,24 @@ pmReadout *ref = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); // Reference image if (!ref) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find readout."); + psError(PPSUB_ERR_UNKNOWN, false, "Unable to find readout."); return false; } - if (!lowThreshold(ref, thresh, maskVal, maskThresh, "reference convolved image")) { - psError(PS_ERR_UNKNOWN, false, "Unable to threshold input image."); - return false; + + psMetadataIterator *regIter = psMetadataIteratorAlloc(in->analysis, PS_LIST_HEAD, + "^" PM_SUBTRACTION_ANALYSIS_REGION "$"); + psMetadataItem *regItem; // Item with region + while ((regItem = psMetadataGetAndIncrement(regIter))) { + psAssert(regItem->type == PS_DATA_REGION && regItem->data.V, "Expect region type"); + psRegion *region = regItem->data.V; // Region of interest + if (!lowThreshold(in, thresh, maskVal, maskThresh, region, "input convolved image")) { + psError(psErrorCodeLast(), false, "Unable to threshold input image."); + return false; + } + if (!lowThreshold(ref, thresh, maskVal, maskThresh, region, "reference convolved image")) { + psError(psErrorCodeLast(), false, "Unable to threshold input image."); + return false; + } } + psFree(regIter); psFree(view); Index: branches/tap_branches/ppSub/src/ppSubVarianceRescale.c =================================================================== --- branches/tap_branches/ppSub/src/ppSubVarianceRescale.c (revision 25900) +++ branches/tap_branches/ppSub/src/ppSubVarianceRescale.c (revision 27838) @@ -35,15 +35,15 @@ int num = psMetadataLookupS32(&mdok, recipe, "RENORM.NUM"); if (!mdok) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "RENORM.NUM is not set in the recipe"); + psError(PPSUB_ERR_ARGUMENTS, true, "RENORM.NUM is not set in the recipe"); return false; } float minValid = psMetadataLookupF32(&mdok, recipe, "RENORM.MIN"); if (!mdok) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "RENORM.MIN is not set in the recipe"); + psError(PPSUB_ERR_ARGUMENTS, true, "RENORM.MIN is not set in the recipe"); return false; } float maxValid = psMetadataLookupF32(&mdok, recipe, "RENORM.MAX"); if (!mdok) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "RENORM.MAX is not set in the recipe"); + psError(PPSUB_ERR_ARGUMENTS, true, "RENORM.MAX is not set in the recipe"); return false; } @@ -54,4 +54,10 @@ pmReadout *readout = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT"); // Output image + if (!readout->variance) { + // Nothing to renormalise + psWarning("Renormalisation of the variance requested, but no variance provided."); + return true; + } + return pmReadoutVarianceRenormalise(readout, maskBad, num, minValid, maxValid); }