Index: branches/tap_branches/ppStack/src/Makefile.am =================================================================== --- branches/tap_branches/ppStack/src/Makefile.am (revision 25900) +++ branches/tap_branches/ppStack/src/Makefile.am (revision 27838) @@ -47,5 +47,6 @@ ppStackCleanup.c \ ppStackPhotometry.c \ - ppStackFinish.c + ppStackFinish.c \ + ppStackErrorCodes.c noinst_HEADERS = \ @@ -53,7 +54,17 @@ ppStackLoop.h \ ppStackOptions.h \ - ppStackThread.h + ppStackThread.h \ + ppStackErrorCodes.h -CLEANFILES = *~ +### Error codes. +BUILT_SOURCES = ppStackErrorCodes.h ppStackErrorCodes.c +CLEANFILES = ppStackErrorCodes.h ppStackErrorCodes.c + +ppStackErrorCodes.h : ppStackErrorCodes.dat ppStackErrorCodes.h.in + $(ERRORCODES) --data=ppStackErrorCodes.dat --outdir=. ppStackErrorCodes.h + +ppStackErrorCodes.c : ppStackErrorCodes.dat ppStackErrorCodes.c.in ppStackErrorCodes.h + $(ERRORCODES) --data=ppStackErrorCodes.dat --outdir=. ppStackErrorCodes.c + clean-local: Index: branches/tap_branches/ppStack/src/ppStack.c =================================================================== --- branches/tap_branches/ppStack/src/ppStack.c (revision 25900) +++ branches/tap_branches/ppStack/src/ppStack.c (revision 27838) @@ -15,21 +15,15 @@ int main(int argc, char *argv[]) { - psExit exitValue = PS_EXIT_SUCCESS; // Exit value + psLibInit(NULL); psTimerStart(TIMER_NAME); psTimerStart("PPSTACK_STEPS"); - psLibInit(NULL); + + pmErrorRegister(); + ppStackErrorRegister(); + psphotErrorRegister(); pmConfig *config = pmConfigRead(&argc, argv, PPSTACK_RECIPE); // Configuration + ppStackOptions *options = NULL; // Options for stacking if (!config) { - psErrorStackPrint(stderr, "Error reading configuration."); - exitValue = PS_EXIT_CONFIG_ERROR; - goto die; - } - - (void) psTraceSetLevel("ppStack", 5); - - if (!ppStackArgumentsSetup(argc, argv, config)) { - psErrorStackPrint(stderr, "Error reading arguments.\n"); - exitValue = PS_EXIT_CONFIG_ERROR; goto die; } @@ -38,45 +32,100 @@ if (!pmModelClassInit()) { - psErrorStackPrint(stderr, "Error initialising model classes.\n"); - exitValue = PS_EXIT_PROG_ERROR; + psError(PPSTACK_ERR_PROG, false, "Unable to initialise model classes."); goto die; } if (!psphotInit()) { - psErrorStackPrint(stderr, "Error initialising psphot.\n"); - exitValue = PS_EXIT_PROG_ERROR; + psError(PPSTACK_ERR_PROG, false, "Error initialising psphot."); + goto die; + } + + (void)psTraceSetLevel("ppStack", 5); + + if (!ppStackArgumentsSetup(argc, argv, config)) { goto die; } if (!ppStackCamera(config)) { - psErrorStackPrint(stderr, "Error setting up input files.\n"); - exitValue = PS_EXIT_CONFIG_ERROR; goto die; } if (!ppStackArgumentsParse(config)) { - psErrorStackPrint(stderr, "Error reading arguments.\n"); - exitValue = PS_EXIT_CONFIG_ERROR; goto die; } - if (!ppStackLoop(config)) { - psErrorStackPrint(stderr, "Error performing combination.\n"); - exitValue = PS_EXIT_DATA_ERROR; + options = ppStackOptionsAlloc(); + if (!ppStackLoop(config, options)) { goto die; } - // Common code for the death. -die: - psTrace("ppStack", 1, "Finished at %f sec\n", psTimerMark(TIMER_NAME)); - psTimerStop(); + die: + // Common code for the death. + { + psExit exitValue = ppStackExitCode(PS_EXIT_SUCCESS); // Exit code - psFree(config); - pmModelClassCleanup(); - pmConfigDone(); - psLibFinalize(); - pmVisualClose(); - exit(exitValue); + // Ensure everything closes + if (config) { + ppStackFileActivation(config, PPSTACK_FILES_PREPARE, true); + ppStackFileActivation(config, PPSTACK_FILES_CONVOLVE, true); + ppStackFileActivation(config, PPSTACK_FILES_STACK, true); + ppStackFileActivation(config, PPSTACK_FILES_UNCONV, true); + ppStackFileActivation(config, PPSTACK_FILES_PHOT, true); + if (!ppStackFilesIterateUp(config)) { + psError(psErrorCodeLast(), false, "Unable to close files."); + exitValue = ppStackExitCode(exitValue); + pmFPAfileFreeSetStrict(false); + } + } + + // Write out summary statistics + if (options && options->stats) { + + psMetadataAddS32(options->stats, PS_LIST_TAIL, "QUALITY", PS_META_REPLACE, + "Bad data quality flag", options->quality); + psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_STACK", 0, + "Total time", psTimerClear("PPSTACK_TOTAL")); + + const char *statsMDC = psMetadataConfigFormat(options->stats); + if (!statsMDC || strlen(statsMDC) == 0) { + psError(PS_ERR_IO, false, "Unable to get statistics MDC file."); + return false; + } + if (fprintf(options->statsFile, "%s", statsMDC) != strlen(statsMDC)) { + psError(PS_ERR_IO, false, "Unable to write statistics MDC file."); + return false; + } + psFree(statsMDC); + if (fclose(options->statsFile) == EOF) { + psError(PS_ERR_IO, false, "Unable to close statistics MDC file."); + return false; + } + options->statsFile = NULL; + pmConfigRunFilenameAddWrite(config, "STATS", + psMetadataLookupStr(NULL, config->arguments, "STATS")); + } + psFree(options); + + // Dump configuration + bool mdok; // Status of MD lookup + psString dump = psMetadataLookupStr(&mdok, config->arguments, "DUMP_CONFIG"); // File for config + if (dump && !pmConfigDump(config, dump)) { + psError(psErrorCodeLast(), false, "Unable to dump configuration."); + exitValue = ppStackExitCode(exitValue); + } + + psTrace("ppStack", 1, "Finished at %f sec\n", psTimerMark(TIMER_NAME)); + psTimerStop(); + + psFree(config); + pmModelClassCleanup(); + pmConfigDone(); + psLibFinalize(); + pmVisualClose(); + + exitValue = ppStackExitCode(exitValue); + exit(exitValue); + } } Index: branches/tap_branches/ppStack/src/ppStack.h =================================================================== --- branches/tap_branches/ppStack/src/ppStack.h (revision 25900) +++ branches/tap_branches/ppStack/src/ppStack.h (revision 27838) @@ -9,12 +9,15 @@ #include "ppStackOptions.h" +#include "ppStackErrorCodes.h" // Mask values for inputs typedef enum { + PPSTACK_MASK_NONE = 0x00, // Nothing wrong PPSTACK_MASK_CAL = 0x01, // Photometric calibration failed - PPSTACK_MASK_MATCH = 0x02, // PSF-matching failed - PPSTACK_MASK_CHI2 = 0x04, // Chi^2 too deviant - PPSTACK_MASK_REJECT = 0x08, // Rejection failed - PPSTACK_MASK_BAD = 0x10, // Bad image (too many pixels rejected) + PPSTACK_MASK_PSF = 0x02, // PSF measurement failed + PPSTACK_MASK_MATCH = 0x04, // PSF-matching failed + PPSTACK_MASK_CHI2 = 0x08, // Chi^2 too deviant + PPSTACK_MASK_REJECT = 0x10, // Rejection failed + PPSTACK_MASK_BAD = 0x20, // Bad image (too many pixels rejected) PPSTACK_MASK_ALL = 0xff // All errors } ppStackMask; @@ -24,5 +27,6 @@ PPSTACK_FILES_PREPARE, // Files for preparation PPSTACK_FILES_CONVOLVE, // Files for convolution - PPSTACK_FILES_COMBINE, // Files for combination + PPSTACK_FILES_STACK, // Stack files + PPSTACK_FILES_UNCONV, // Unconvolved stack files PPSTACK_FILES_PHOT // Files for photometry } ppStackFileList; @@ -59,5 +63,5 @@ // Perform stacking on a readout // -// Returns an array of pixels to inspect for each input image +// Returns two arrays: pixels to inspect for each input image, and pixels to reject for each input image. psArray *ppStackReadoutInitial(const pmConfig *config, // Configuration pmReadout *outRO, // Output readout @@ -65,4 +69,5 @@ const psVector *mask, // Mask for input readouts const psVector *weightings, // Weighting factors for each image + const psVector *exposures, // Exposure time for each image const psVector *addVariance // Additional variance for rejection ); @@ -79,9 +84,13 @@ bool ppStackReadoutFinal(const pmConfig *config, // Configuration pmReadout *outRO, // Output readout + pmReadout *expRO, // Exposure readout const psArray *readouts, // Input readouts const psVector *mask, // Mask for input readouts const psArray *rejected, // Array with pixels rejected in each image const psVector *weightings, // Weighting factors for each image - const psVector *addVariance // Additional variance for rejection + const psVector *exposures, // Exposure times for each image + const psVector *addVariance, // Additional variance for rejection + bool safety, // Enable safety switch? + const psVector *norm // Normalisations to apply ); @@ -164,4 +173,6 @@ ); +/// Return an appropriate exit code based on the error code +psExit ppStackExitCode(psExit exitValue); #endif Index: branches/tap_branches/ppStack/src/ppStackArguments.c =================================================================== --- branches/tap_branches/ppStack/src/ppStackArguments.c (revision 25900) +++ branches/tap_branches/ppStack/src/ppStackArguments.c (revision 27838) @@ -43,5 +43,5 @@ value = psMetadataLookup##TYPE(&mdok, recipe, RECIPENAME); \ if (!mdok) { \ - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to find %s in recipe %s", \ + psError(PPSTACK_ERR_CONFIG, true, "Unable to find %s in recipe %s", \ RECIPENAME, PPSTACK_RECIPE); \ goto ERROR; \ @@ -58,5 +58,5 @@ value = psMetadataLookup##TYPE(&mdok, recipe, RECIPENAME); \ if (!mdok) { \ - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to find %s in recipe %s", \ + psError(PPSTACK_ERR_CONFIG, true, "Unable to find %s in recipe %s", \ RECIPENAME, PPSTACK_RECIPE); \ goto ERROR; \ @@ -73,6 +73,6 @@ name = psMetadataLookupStr(NULL, recipe, RECIPENAME); \ if (!name) { \ - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to find %s in recipe %s", \ - RECIPENAME, PPSTACK_RECIPE); \ + psError(PPSTACK_ERR_CONFIG, true, "Unable to find %s in recipe %s", \ + RECIPENAME, PPSTACK_RECIPE); \ goto ERROR; \ } \ @@ -108,5 +108,5 @@ value = psMetadataLookupStr(NULL, recipe, mdName); if (!value) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to find %s in recipe %s", + psError(PPSTACK_ERR_CONFIG, true, "Unable to find %s in recipe %s", mdName, PPSTACK_RECIPE); return false; @@ -148,5 +148,5 @@ psMetadataAddStr(arguments, PS_LIST_TAIL, "-stamps", 0, "Stamps file with x,y,flux per line", NULL); psMetadataAddStr(arguments, PS_LIST_TAIL, "-stats", 0, "Statistics file", NULL); - psMetadataAddS32(arguments, PS_LIST_TAIL, "-iter", 0, "Number of rejection iterations", 0); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-combine-iter", 0, "Number of rejection iterations per input", NAN); psMetadataAddF32(arguments, PS_LIST_TAIL, "-combine-rej", 0, "Combination rejection thresold (sigma)", NAN); @@ -185,4 +185,5 @@ psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-star-sys-1", 0, "Estimated systematic error; pass 1", NAN); psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp-star-sys-2", 0, "Estimated systematic error; pass 2", NAN); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-temp-dir", 0, "Directory for temporary images", NULL); psMetadataAddStr(arguments, PS_LIST_TAIL, "-temp-image", 0, "Suffix for temporary images", NULL); psMetadataAddStr(arguments, PS_LIST_TAIL, "-temp-mask", 0, "Suffix for temporary masks", NULL); @@ -206,5 +207,5 @@ psMetadata *inputs = psMetadataConfigRead(NULL, &numBad, argv[argNum], false); // Input file info if (!inputs || numBad > 0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to cleanly read MDC file with inputs."); + psError(PPSTACK_ERR_ARGUMENTS, false, "Unable to cleanly read MDC file with inputs."); return false; } @@ -228,5 +229,5 @@ int numThreads = psMetadataLookupS32(NULL, arguments, "-threads"); // Number of threads if (numThreads > 0 && !psThreadPoolInit(numThreads)) { - psError(PS_ERR_UNKNOWN, false, "Unable to setup %d threads", numThreads); + psError(PPSTACK_ERR_ARGUMENTS, false, "Unable to setup %d threads", numThreads); return false; } @@ -246,9 +247,9 @@ psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // Recipe if (!recipe) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find recipe %s", PPSTACK_RECIPE); + psError(PPSTACK_ERR_CONFIG, false, "Unable to find recipe %s", PPSTACK_RECIPE); goto ERROR; } - VALUE_ARG_RECIPE_INT("-iter", "ITER", S32, 0); + VALUE_ARG_RECIPE_FLOAT("-combine-iter", "COMBINE.ITER", F32); VALUE_ARG_RECIPE_FLOAT("-combine-rej", "COMBINE.REJ", F32); VALUE_ARG_RECIPE_FLOAT("-combine-sys", "COMBINE.SYS", F32); @@ -260,5 +261,5 @@ VALUE_ARG_RECIPE_FLOAT("-poor-frac", "POOR.FRACTION", F32); - valueArgRecipeStr(arguments, recipe, "-mask-val", "MASK.VAL", recipe); + valueArgRecipeStr(arguments, recipe, "-mask-val", "MASK.IN", recipe); valueArgRecipeStr(arguments, recipe, "-mask-bad", "MASK.BAD", recipe); valueArgRecipeStr(arguments, recipe, "-mask-poor", "MASK.POOR", recipe); @@ -316,5 +317,5 @@ pmConfigCamerasCull(config, NULL); - pmConfigRecipesCull(config, "PPSTACK,PPSUB,PPSTATS,PSPHOT,MASKS,JPEG"); + pmConfigRecipesCull(config, "PPSTACK,PPSUB,PPSTATS,PSPHOT,PSASTRO,MASKS,JPEG"); return true; Index: branches/tap_branches/ppStack/src/ppStackCamera.c =================================================================== --- branches/tap_branches/ppStack/src/ppStackCamera.c (revision 25900) +++ branches/tap_branches/ppStack/src/ppStackCamera.c (revision 27838) @@ -30,9 +30,9 @@ pmFPAfileDefineFromArgs(&found, config, name, "FILENAMES"); if (!file || !found) { - psError(PS_ERR_UNKNOWN, false, "Unable to define file %s from %s", name, filename); + psError(psErrorCodeLast(), false, "Unable to define file %s from %s", name, filename); return NULL; } if (file->type != type) { - psError(PS_ERR_IO, true, "%s is not of type %s", name, pmFPAfileStringFromType(type)); + psError(PS_ERR_IO, PPSTACK_ERR_CONFIG, "%s is not of type %s", name, pmFPAfileStringFromType(type)); return NULL; } @@ -53,5 +53,5 @@ psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // Recipe for ppSim if (!recipe) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find recipe %s", PPSTACK_RECIPE); + psError(PPSTACK_ERR_CONFIG, false, "Unable to find recipe %s", PPSTACK_RECIPE); return false; } @@ -67,5 +67,5 @@ "PPSTACK.INPUT.MASK"); // Input masks if (!status) { - psError(PS_ERR_UNKNOWN, false, "Unable to define input masks from RUN metadata."); + psError(psErrorCodeLast(), false, "Unable to define input masks from RUN metadata."); psFree(runImages); return false; @@ -76,5 +76,5 @@ "PPSTACK.INPUT.VARIANCE"); // Input variances if (!status) { - psError(PS_ERR_UNKNOWN, false, "Unable to define input variances from RUN metadata."); + psError(psErrorCodeLast(), false, "Unable to define input variances from RUN metadata."); psFree(runImages); return false; @@ -88,10 +88,10 @@ "PPSTACK.INPUT.SOURCES"); // Input sources if (!status) { - psError(PS_ERR_UNKNOWN, false, "Unable to define input sources from RUN metadata."); + psError(psErrorCodeLast(), false, "Unable to define input sources from RUN metadata."); psFree(runImages); return false; } if (!runSrc) { - psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to define input sources from RUN metadata."); + psError(PPSTACK_ERR_CONFIG, true, "Unable to define input sources from RUN metadata."); psFree(runImages); return false; @@ -104,5 +104,5 @@ "PPSTACK.INPUT.PSF"); // Input PSFs if (!status) { - psError(PS_ERR_UNKNOWN, false, "Unable to define input PSFs from RUN metadata."); + psError(psErrorCodeLast(), false, "Unable to define input PSFs from RUN metadata."); psFree(runImages); return false; @@ -118,10 +118,11 @@ "PPSTACK.CONV.KERNEL"); // Conv'n kernels if (!status) { - psError(PS_ERR_UNKNOWN, false, "Unable to define convolution kernels from RUN metadata."); + psError(psErrorCodeLast(), false, + "Unable to define convolution kernels from RUN metadata."); psFree(runImages); return false; } if (!runKernel) { - psError(PS_ERR_UNEXPECTED_NULL, true, + psError(PPSTACK_ERR_CONFIG, true, "Unable to define convolution kernels from RUN metadata."); psFree(runImages); @@ -137,5 +138,5 @@ psMetadata *inputs = psMetadataLookupMetadata(NULL, config->arguments, "INPUTS"); // The inputs info if (!inputs) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find inputs."); + psError(PPSTACK_ERR_ARGUMENTS, false, "Unable to find inputs."); return false; } @@ -145,5 +146,5 @@ while ((item = psMetadataGetAndIncrement(iter))) { if (item->type != PS_DATA_METADATA) { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, + psError(PPSTACK_ERR_ARGUMENTS, true, "Component %s of the input metadata is not of type METADATA", item->name); psFree(iter); @@ -155,5 +156,5 @@ psString image = psMetadataLookupStr(NULL, input, "IMAGE"); // Name of image if (!image || strlen(image) == 0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Component %s lacks IMAGE of type STR", item->name); + psError(PPSTACK_ERR_ARGUMENTS, true, "Component %s lacks IMAGE of type STR", item->name); psFree(iter); return false; @@ -169,5 +170,5 @@ image, PM_FPA_FILE_IMAGE); // File for image if (!imageFile) { - psError(PS_ERR_UNKNOWN, false, "Unable to define file from image %d (%s)", i, image); + psError(psErrorCodeLast(), false, "Unable to define file from image %d (%s)", i, image); return false; } @@ -175,5 +176,5 @@ if (mask && strlen(mask) > 0 && !defineFile(config, imageFile, "PPSTACK.INPUT.MASK", mask, PM_FPA_FILE_MASK)) { - psError(PS_ERR_UNKNOWN, false, "Unable to define file from mask %d (%s)", i, mask); + psError(psErrorCodeLast(), false, "Unable to define file from mask %d (%s)", i, mask); return false; } @@ -183,5 +184,5 @@ if (!defineFile(config, imageFile, "PPSTACK.INPUT.VARIANCE", variance, PM_FPA_FILE_VARIANCE)) { - psError(PS_ERR_UNKNOWN, false, + psError(psErrorCodeLast(), false, "Unable to define file from variance %d (%s)", i, variance); return false; @@ -195,19 +196,19 @@ havePSFs = true; if (!defineFile(config, imageFile, "PPSTACK.INPUT.PSF", psf, PM_FPA_FILE_PSF)) { - psError(PS_ERR_UNKNOWN, false, "Unable to define file from psf %d (%s)", i, psf); + psError(psErrorCodeLast(), false, "Unable to define file from psf %d (%s)", i, psf); return false; } } } else if (havePSFs) { - psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find PSF %d", i); + psError(PPSTACK_ERR_CONFIG, true, "Unable to find PSF %d", i); return false; } if (!sources || strlen(sources) == 0) { - psError(PS_ERR_UNEXPECTED_NULL, true, "SOURCES not provided for file %d", i); + psError(PPSTACK_ERR_CONFIG, true, "SOURCES not provided for file %d", i); return false; } if (!defineFile(config, imageFile, "PPSTACK.INPUT.SOURCES", sources, PM_FPA_FILE_CMF)) { - psError(PS_ERR_UNKNOWN, false, "Unable to define file from sources %d (%s)", + psError(psErrorCodeLast(), false, "Unable to define file from sources %d (%s)", i, sources); return false; @@ -217,5 +218,6 @@ pmFPAfile *kernel = pmFPAfileDefineOutput(config, imageFile->fpa, "PPSTACK.CONV.KERNEL"); if (!kernel) { - psError(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.CONV.KERNEL")); + psError(psErrorCodeLast(), false, + "Unable to generate output file from PPSTACK.CONV.KERNEL"); return false; } @@ -233,24 +235,23 @@ // Output image - pmFPA *fpa = pmFPAConstruct(config->camera, config->cameraName); // FPA to contain the output - if (!fpa) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to construct an FPA from camera configuration."); - return false; - } - pmFPAfile *output = pmFPAfileDefineOutput(config, fpa, "PPSTACK.OUTPUT"); - psFree(fpa); // Drop reference + pmFPA *outFPA = pmFPAConstruct(config->camera, config->cameraName); // FPA to contain the output + if (!outFPA) { + psError(psErrorCodeLast(), false, "Unable to construct an FPA from camera configuration."); + return false; + } + pmFPAfile *output = pmFPAfileDefineOutput(config, outFPA, "PPSTACK.OUTPUT"); + psFree(outFPA); // Drop reference if (!output) { - psError(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.OUTPUT")); + psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.OUTPUT")); return false; } if (output->type != PM_FPA_FILE_IMAGE) { - psError(PS_ERR_IO, true, "PPSTACK.OUTPUT is not of type IMAGE"); + psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.OUTPUT is not of type IMAGE"); return false; } output->save = true; - if (!pmFPAAddSourceFromFormat(fpa, "Stack", output->format)) { - psError(PS_ERR_UNKNOWN, false, "Unable to generate output FPA."); - psFree(fpa); + if (!pmFPAAddSourceFromFormat(outFPA, "Stack", output->format)) { + psError(psErrorCodeLast(), false, "Unable to generate output FPA."); return false; } @@ -259,9 +260,9 @@ pmFPAfile *outMask = pmFPAfileDefineOutput(config, output->fpa, "PPSTACK.OUTPUT.MASK"); if (!outMask) { - psError(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.OUTPUT.MASK")); + psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.OUTPUT.MASK")); return false; } if (outMask->type != PM_FPA_FILE_MASK) { - psError(PS_ERR_IO, true, "PPSTACK.OUTPUT.MASK is not of type MASK"); + psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.OUTPUT.MASK is not of type MASK"); return false; } @@ -272,9 +273,9 @@ pmFPAfile *outVariance = pmFPAfileDefineOutput(config, output->fpa, "PPSTACK.OUTPUT.VARIANCE"); if (!outVariance) { - psError(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.OUTPUT.VARIANCE")); + psError(psErrorCodeLast(), false, "Unable to generate output file from PPSTACK.OUTPUT.VARIANCE"); return false; } if (outVariance->type != PM_FPA_FILE_VARIANCE) { - psError(PS_ERR_IO, true, "PPSTACK.OUTPUT.VARIANCE is not of type VARIANCE"); + psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.OUTPUT.VARIANCE is not of type VARIANCE"); return false; } @@ -282,25 +283,177 @@ } + + // Exposure image + pmFPA *expFPA = pmFPAConstruct(config->camera, config->cameraName); // FPA to contain the output + if (!expFPA) { + psError(psErrorCodeLast(), false, "Unable to construct an FPA from camera configuration."); + return false; + } + pmFPAfile *exp = pmFPAfileDefineOutput(config, expFPA, "PPSTACK.OUTPUT.EXP"); + psFree(expFPA); // Drop reference + if (!exp) { + psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.OUTPUT.EXP")); + return false; + } + if (exp->type != PM_FPA_FILE_IMAGE) { + psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.OUTPUT.EXP is not of type IMAGE"); + return false; + } + exp->save = true; + + if (!pmFPAAddSourceFromFormat(expFPA, "Stack", exp->format)) { + psError(psErrorCodeLast(), false, "Unable to generate output FPA."); + return false; + } + + // Exposure numbers + pmFPAfile *expNum = pmFPAfileDefineOutput(config, exp->fpa, "PPSTACK.OUTPUT.EXPNUM"); + if (!expNum) { + psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.OUTPUT.EXPNUM")); + return false; + } + if (expNum->type != PM_FPA_FILE_MASK) { + psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.OUTPUT.EXPNUM is not of type MASK"); + return false; + } + expNum->save = true; + + // Weighted exposure + pmFPAfile *expWt = pmFPAfileDefineOutput(config, exp->fpa, "PPSTACK.OUTPUT.EXPWT"); + if (!expWt) { + psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.OUTPUT.EXPWT")); + return false; + } + if (expWt->type != PM_FPA_FILE_VARIANCE) { + psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.OUTPUT.EXPWT is not of type VARIANCE"); + return false; + } + expWt->save = true; + + if (havePSFs) { - pmFPAfile *targetPSF = pmFPAfileDefineOutput(config, output->fpa, "PPSTACK.TARGET.PSF"); + pmFPA *psfFPA = pmFPAConstruct(config->camera, config->cameraName); // FPA to contain PSF + if (!psfFPA) { + psError(psErrorCodeLast(), false, "Unable to construct an FPA from camera configuration."); + return false; + } + pmFPAfile *targetPSF = pmFPAfileDefineOutput(config, psfFPA, "PPSTACK.TARGET.PSF"); if (!targetPSF) { - psError(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.TARGET.PSF")); - return false; - } + psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.TARGET.PSF")); + return false; + } + psFree(psfFPA); if (targetPSF->type != PM_FPA_FILE_PSF) { - psError(PS_ERR_IO, true, "PPSTACK.TARGET.PSF is not of type PSF"); + psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.TARGET.PSF is not of type PSF"); return false; } targetPSF->save = true; } + + // Unconvolved stack + pmFPA *unconvFPA = pmFPAConstruct(config->camera, config->cameraName); // FPA to contain unconvolved output + if (!unconvFPA) { + psError(psErrorCodeLast(), false, "Unable to construct an FPA from camera configuration."); + return false; + } + pmFPAfile *unConv = pmFPAfileDefineOutput(config, unconvFPA, "PPSTACK.UNCONV"); + psFree(unconvFPA); // Drop reference + if (!unConv) { + psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.UNCONV")); + return false; + } + if (unConv->type != PM_FPA_FILE_IMAGE) { + psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.UNCONV is not of type IMAGE"); + return false; + } + unConv->save = true; + + if (!pmFPAAddSourceFromFormat(unconvFPA, "Stack", unConv->format)) { + psError(psErrorCodeLast(), false, "Unable to generate output FPA."); + return false; + } + + // Unconvolved mask + pmFPAfile *unconvMask = pmFPAfileDefineOutput(config, unconvFPA, "PPSTACK.UNCONV.MASK"); + if (!unconvMask) { + psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.UNCONV.MASK")); + return false; + } + if (unconvMask->type != PM_FPA_FILE_MASK) { + psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.UNCONV.MASK is not of type MASK"); + return false; + } + unconvMask->save = true; + + // Unconvolved variance + if (haveVariances) { + pmFPAfile *unconvVariance = pmFPAfileDefineOutput(config, unconvFPA, "PPSTACK.UNCONV.VARIANCE"); + if (!unconvVariance) { + psError(psErrorCodeLast(), false, "Unable to generate output file from PPSTACK.UNCONV.VARIANCE"); + return false; + } + if (unconvVariance->type != PM_FPA_FILE_VARIANCE) { + psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.UNCONV.VARIANCE is not of type VARIANCE"); + return false; + } + unconvVariance->save = true; + } + + + // Exposure image + pmFPA *unconvExpFPA = pmFPAConstruct(config->camera, config->cameraName); // FPA to contain the output + if (!unconvExpFPA) { + psError(psErrorCodeLast(), false, "Unable to construct an FPA from camera configuration."); + return false; + } + pmFPAfile *unconvExp = pmFPAfileDefineOutput(config, unconvExpFPA, "PPSTACK.UNCONV.EXP"); + psFree(unconvExpFPA); // Drop reference + if (!unconvExp) { + psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.UNCONV.EXP")); + return false; + } + if (unconvExp->type != PM_FPA_FILE_IMAGE) { + psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.UNCONV.EXP is not of type IMAGE"); + return false; + } + unconvExp->save = true; + + if (!pmFPAAddSourceFromFormat(unconvExpFPA, "Stack", unconvExp->format)) { + psError(psErrorCodeLast(), false, "Unable to generate output FPA."); + return false; + } + + // Exposure numbers + pmFPAfile *unconvExpNum = pmFPAfileDefineOutput(config, unconvExp->fpa, "PPSTACK.UNCONV.EXPNUM"); + if (!unconvExpNum) { + psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.UNCONV.MASK")); + return false; + } + if (unconvExpNum->type != PM_FPA_FILE_MASK) { + psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.UNCONV.EXPNUM is not of type MASK"); + return false; + } + unconvExpNum->save = true; + + // Weighted exposure + pmFPAfile *unconvExpWt = pmFPAfileDefineOutput(config, unconvExp->fpa, "PPSTACK.UNCONV.EXPWT"); + if (!unconvExpWt) { + psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.UNCONV.EXPWT")); + return false; + } + if (unconvExpWt->type != PM_FPA_FILE_VARIANCE) { + psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.UNCONV.EXPWT is not of type VARIANCE"); + return false; + } + unconvExpWt->save = true; // Output JPEGs pmFPAfile *jpeg1 = pmFPAfileDefineOutput(config, NULL, "PPSTACK.OUTPUT.JPEG1"); if (!jpeg1) { - psError(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.OUTPUT.JPEG1")); + psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.OUTPUT.JPEG1")); return false; } if (jpeg1->type != PM_FPA_FILE_JPEG) { - psError(PS_ERR_IO, true, "PPSTACK.OUTPUT.JPEG1 is not of type JPEG"); + psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.OUTPUT.JPEG1 is not of type JPEG"); return false; } @@ -308,9 +461,9 @@ pmFPAfile *jpeg2 = pmFPAfileDefineOutput(config, NULL, "PPSTACK.OUTPUT.JPEG2"); if (!jpeg2) { - psError(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.OUTPUT.JPEG2")); + psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.OUTPUT.JPEG2")); return false; } if (jpeg2->type != PM_FPA_FILE_JPEG) { - psError(PS_ERR_IO, true, "PPSTACK.OUTPUT.JPEG2 is not of type JPEG"); + psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.OUTPUT.JPEG2 is not of type JPEG"); return false; } @@ -323,15 +476,17 @@ psMetadataLookupBool(&mdok, config->arguments, "-photometry"); // perform photometry if (doPhotom) { - // This file, PSPHOT.INPUT, is just used as a carrier; output files (eg, PSPHOT.RESID) are defined by - // psphotDefineFiles + // This pmFPAfile, PSPHOT.INPUT, is just used as a carrier; output files (eg, + // PSPHOT.RESID) are defined by psphotDefineFiles pmFPAfile *psphotInput = pmFPAfileDefineFromFPA(config, output->fpa, 1, 1, "PSPHOT.INPUT"); if (!psphotInput) { - psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.INPUT")); - return false; - } + psError(psErrorCodeLast(), false, _("Unable to generate output file from PSPHOT.INPUT")); + 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); // Define associated psphot input/output files if (!psphotDefineFiles(config, psphotInput)) { - psError(PSPHOT_ERR_CONFIG, false, + psError(psErrorCodeLast(), false, "Trouble defining the additional input/output files for psphot"); return false; @@ -341,9 +496,9 @@ pmFPAfile *outPSF = pmFPAfileDefineOutputFromFile(config, output, "PSPHOT.PSF.SAVE"); if (!outPSF) { - psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.PSF.SAVE")); + psError(psErrorCodeLast(), false, _("Unable to generate output file from PSPHOT.PSF.SAVE")); return false; } if (outPSF->type != PM_FPA_FILE_PSF) { - psError(PS_ERR_IO, true, "PSPHOT.PSF.SAVE is not of type PSF"); + psError(PPSTACK_ERR_CONFIG, true, "PSPHOT.PSF.SAVE is not of type PSF"); return false; } Index: branches/tap_branches/ppStack/src/ppStackCleanup.c =================================================================== --- branches/tap_branches/ppStack/src/ppStackCleanup.c (revision 25900) +++ branches/tap_branches/ppStack/src/ppStackCleanup.c (revision 27838) @@ -4,12 +4,10 @@ #include -#include #include #include +#include #include "ppStack.h" #include "ppStackLoop.h" - -#define WCS_TOLERANCE 0.001 // Tolerance for WCS @@ -22,7 +20,4 @@ psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe psAssert(recipe, "We've thrown an error on this before."); - - bool mdok; // Status of MD lookup - bool tempDelete = psMetadataLookupBool(&mdok, recipe, "TEMP.DELETE"); // Delete temporary files? #if 0 @@ -37,59 +32,4 @@ #endif - // Close up - bool wcsDone = false; // Have we done the WCS? - for (int i = 0; i < options->num; i++) { - if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) { - continue; - } - - ppStackThread *thread = stack->threads->data[0]; // Representative stack - pmReadout *inRO = thread->readouts->data[i]; // Template readout - if (inRO && !wcsDone) { - // Copy astrometry over - wcsDone = true; - pmHDU *inHDU = pmHDUFromCell(inRO->parent); // Template HDU - pmHDU *outHDU = pmHDUFromCell(options->outRO->parent); // Output HDU - pmChip *outChip = options->outRO->parent->parent; // Output chip - pmFPA *outFPA = outChip->parent; // Output FPA - if (!outHDU || !inHDU) { - psWarning("Unable to find HDU at FPA level to copy astrometry."); - } else { - if (!pmAstromReadWCS(outFPA, outChip, inHDU->header, 1.0)) { - psErrorClear(); - psWarning("Unable to read WCS astrometry from input FPA."); - wcsDone = false; - } else { - if (!outHDU->header) { - outHDU->header = psMetadataAlloc(); - } - if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_TOLERANCE)) { - psErrorClear(); - psWarning("Unable to write WCS astrometry to output FPA."); - wcsDone = false; - } - } - } - } - - if (tempDelete) { - psString imageResolved = pmConfigConvertFilename(options->imageNames->data[i], - config, false, false); - psString maskResolved = pmConfigConvertFilename(options->maskNames->data[i], - config, false, false); - psString varianceResolved = pmConfigConvertFilename(options->varianceNames->data[i], - config, false, false); - if (unlink(imageResolved) == -1 || unlink(maskResolved) == -1 || - unlink(varianceResolved) == -1) { - psWarning("Unable to delete temporary files for image %d", i); - } - psFree(imageResolved); - psFree(maskResolved); - psFree(varianceResolved); - } - } - - ppStackMemDump("cleanup"); - // Generate binned JPEGs { @@ -103,9 +43,10 @@ pmCell *cell2 = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT.JPEG2"); psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask + psFree(view); pmReadout *ro1 = pmReadoutAlloc(cell1), *ro2 = pmReadoutAlloc(cell2); // Binned readouts if (!pmReadoutRebin(ro1, options->outRO, maskValue, bin1, bin1) || !pmReadoutRebin(ro2, ro1, 0, bin2, bin2)) { - psError(PS_ERR_UNKNOWN, false, "Unable to bin output."); + psError(PPSTACK_ERR_DATA, false, "Unable to bin output."); psFree(ro1); psFree(ro2); @@ -116,14 +57,62 @@ } - // Put version information into the header - pmHDU *hdu = pmHDUFromCell(options->outRO->parent); - if (!hdu) { - psError(PS_ERR_UNKNOWN, false, "Unable to find HDU for output."); + // Statistics on output + if (options->stats) { + psTrace("ppStack", 1, "Gathering statistics on stacked image....\n"); + psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits for bad + psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels + + pmFPAview *view = pmFPAviewAlloc(0); // View to readout + view->chip = view->cell = view->readout = 0; + + ppStatsFPA(options->stats, options->outRO->parent->parent->parent, view, maskBad, config); + + psFree(view); + } + + if (!ppStackFilesIterateUp(config)) { + psError(psErrorCodeLast(), false, "Unable to close files."); return false; } - if (!hdu->header) { - hdu->header = psMetadataAlloc(); + ppStackFileActivation(config, PPSTACK_FILES_STACK, false); + ppStackFileActivation(config, PPSTACK_FILES_PHOT, false); + + // Ensure files are freed + { + options->outRO->data_exists = false; + options->outRO->parent->data_exists = false; + options->outRO->parent->parent->data_exists = false; + psFree(options->outRO); + options->outRO = NULL; + + options->expRO->data_exists = false; + options->expRO->parent->data_exists = false; + options->expRO->parent->parent->data_exists = false; + psFree(options->expRO); + options->expRO = NULL; + + pmFPAview *view = pmFPAviewAlloc(0);// Pointer into FPA hierarchy + view->chip = view->cell = 0; // pmFPAviewFreeData doesn't want to deal with readouts + pmFPAfile *jpeg1 = pmFPAfileSelectSingle(config->files, "PPSTACK.OUTPUT.JPEG1", 0); // JPEG file + pmFPAviewFreeData(view, jpeg1); + pmFPAfile *jpeg2 = pmFPAfileSelectSingle(config->files, "PPSTACK.OUTPUT.JPEG2", 0); // JPEG file + pmFPAviewFreeData(view, jpeg2); + pmFPAfile *phot = NULL; // Photometry file + if (options->photometry) { + phot = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", 0); // Photometry file + pmFPAviewFreeData(view, phot); + } + + view->readout = 0; + pmReadout *ro1 = pmFPAviewThisReadout(view, jpeg1->fpa); // JPEG readout + ro1->data_exists = ro1->parent->data_exists = ro1->parent->parent->data_exists = false; + pmReadout *ro2 = pmFPAviewThisReadout(view, jpeg2->fpa); // JPEG readout + ro2->data_exists = ro2->parent->data_exists = ro2->parent->parent->data_exists = false; + if (options->photometry) { + pmReadout *ro = pmFPAviewThisReadout(view, phot->fpa); // Photometry readout + ro->data_exists = ro->parent->data_exists = ro->parent->parent->data_exists = false; + } + psFree(view); } - ppStackVersionHeader(hdu->header); return true; Index: branches/tap_branches/ppStack/src/ppStackCombineFinal.c =================================================================== --- branches/tap_branches/ppStack/src/ppStackCombineFinal.c (revision 25900) +++ branches/tap_branches/ppStack/src/ppStackCombineFinal.c (revision 27838) @@ -10,9 +10,48 @@ #include "ppStackLoop.h" -bool ppStackCombineFinal(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config) +#define WCS_TOLERANCE 0.001 // Tolerance for WCS + +//#define TESTING // Enable test output + +bool ppStackCombineFinal(ppStackThreadData *stack, psArray *covariances, ppStackOptions *options, + pmConfig *config, bool safe, bool normalise, bool grow) { psAssert(stack, "Require stack"); psAssert(options, "Require options"); psAssert(config, "Require configuration"); + + pmReadout *outRO = options->outRO; // Output readout + pmReadout *expRO = options->expRO; // Exposure readout + int numCols = outRO->image->numCols, numRows = outRO->image->numRows; // Size of image + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe + psAssert(recipe, "We've thrown an error on this before."); + float poorFrac = psMetadataLookupF32(NULL, recipe, "POOR.FRACTION"); // Fraction for "poor" + + // Grow the list of rejected pixels, if desired + psArray *reject = psArrayAlloc(options->num); // Pixels rejected for each image + for (int i = 0; i < options->num; i++) { + if (options->inputMask->data.U8[i]) { + continue; + } + if (grow) { + reject->data[i] = pmStackRejectGrow(options->rejected->data[i], numCols, numRows, poorFrac, + options->regions->data[i], options->kernels->data[i]); + if (!reject->data[i]) { + psError(psErrorCodeLast(), false, "Unable to grow rejected pixels for image %d", i); + psFree(reject); + return false; + } + } else { + reject->data[i] = psMemIncrRefCounter(options->rejected->data[i]); + } + } + + if (!outRO->mask) { + outRO->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); + } + if (!expRO->mask) { + expRO->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); + } stack->lastScan = 0; // Reset read @@ -22,5 +61,6 @@ if (!status) { // Something went wrong - psError(PS_ERR_UNKNOWN, false, "Unable to read chunk %d", numChunk); + psError(psErrorCodeLast(), false, "Unable to read chunk %d", numChunk); + psFree(reject); return false; } @@ -33,8 +73,12 @@ psThreadJob *job = psThreadJobAlloc("PPSTACK_FINAL_COMBINE"); // Job to start psArrayAdd(job->args, 1, thread); + psArrayAdd(job->args, 1, reject); psArrayAdd(job->args, 1, options); psArrayAdd(job->args, 1, config); + PS_ARRAY_ADD_SCALAR(job->args, safe, PS_TYPE_U8); + PS_ARRAY_ADD_SCALAR(job->args, normalise, PS_TYPE_U8); if (!psThreadJobAddPending(job)) { psFree(job); + psFree(reject); return false; } @@ -43,36 +87,100 @@ if (!psThreadPoolWait(true)) { - psError(PS_ERR_UNKNOWN, false, "Unable to do final combination."); + psError(psErrorCodeLast(), false, "Unable to do final combination."); + psFree(reject); return false; } + psFree(reject); + // Sum covariance matrices - double sumWeights = 0.0; // Sum of weights - for (int i = 0; i < options->num; i++) { - if (options->inputMask->data.U8[i]) { - psFree(options->covariances->data[i]); - options->covariances->data[i] = NULL; + if (covariances) { + outRO->covariance = psImageCovarianceAverageWeighted(covariances, options->weightings); + } else { + outRO->covariance = psImageCovarianceNone(); + } + + // Propagate WCS + bool wcsDone = false; // Have we done the WCS? + for (int i = 0; i < options->num && !wcsDone; i++) { + if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) { continue; } - psKernel *covar = options->covariances->data[i]; // Covariance matrix - if (!covar) { - continue; + + ppStackThread *thread = stack->threads->data[0]; // Representative stack + pmReadout *inRO = thread->readouts->data[i]; // Template readout + if (inRO && !wcsDone) { + // Copy astrometry over + wcsDone = true; + pmHDU *inHDU = pmHDUFromCell(inRO->parent); // Template HDU + pmHDU *outHDU = pmHDUFromCell(outRO->parent); // Output HDU + pmChip *outChip = outRO->parent->parent; // Output chip + pmFPA *outFPA = outChip->parent; // Output FPA + if (!outHDU || !inHDU) { + psWarning("Unable to find HDU at FPA level to copy astrometry."); + } else { + if (!pmAstromReadWCS(outFPA, outChip, inHDU->header, 1.0)) { + psErrorClear(); + psWarning("Unable to read WCS astrometry from input FPA."); + wcsDone = false; + } else { + if (!outHDU->header) { + outHDU->header = psMetadataAlloc(); + } + if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_TOLERANCE)) { + psErrorClear(); + psWarning("Unable to write WCS astrometry to output FPA."); + wcsDone = false; + } + } + } } - float weight = options->weightings->data.F32[i]; // Weight to apply - psBinaryOp(covar->image, covar->image, "*", psScalarAlloc(weight, PS_TYPE_F32)); - sumWeights += weight; - } - if (sumWeights > 0.0) { - pmReadout *outRO = options->outRO; // Output readout - outRO->covariance = psImageCovarianceSum(options->covariances); - psBinaryOp(outRO->covariance->image, outRO->covariance->image, "/", - psScalarAlloc(sumWeights, PS_TYPE_F32)); - psFree(options->covariances); options->covariances = NULL; - psImageCovarianceTransfer(outRO->variance, outRO->covariance); } + // Set exposure time correctly + { + float exptime = 0.0; // Summed exposure time + for (int i = 0; i < options->num; i++) { + if (options->inputMask->data.U8[i]) { + continue; + } + exptime += options->exposures->data.F32[i]; + } + + { + psMetadataItem *item = psMetadataLookup(outRO->parent->concepts, "CELL.EXPOSURE"); + item->data.F32 = exptime; + } + { + psMetadataItem *item = psMetadataLookup(outRO->parent->parent->parent->concepts, "FPA.EXPOSURE"); + item->data.F32 = exptime; + } + } + + // Put version information into the header + pmHDU *hdu = pmHDUFromCell(outRO->parent); + if (!hdu) { + psError(PPSTACK_ERR_PROG, false, "Unable to find HDU for output."); + return false; + } + if (!hdu->header) { + hdu->header = psMetadataAlloc(); + } + ppStackVersionHeader(hdu->header); + + #ifdef TESTING - pmStackVisualPlotTestImage(outRO->image, "combined_initial.fits"); - ppStackWriteImage("combined_final.fits", NULL, outRO->image, config); + static int pass = 0; // Pass through + psString name = NULL; // Name of file + psStringAppend(&name, "combined_image_final_%d.fits", pass); + pass++; + ppStackWriteImage(name, NULL, outRO->image, config); + psStringSubstitute(&name, "mask", "image"); + ppStackWriteImage(name, NULL, outRO->mask, config); + psStringSubstitute(&name, "variance", "mask"); + ppStackWriteImage(name, NULL, outRO->variance, config); + psFree(name); + + pmStackVisualPlotTestImage(outRO->image, "combined_image_final.fits"); #endif Index: branches/tap_branches/ppStack/src/ppStackCombineInitial.c =================================================================== --- branches/tap_branches/ppStack/src/ppStackCombineInitial.c (revision 25900) +++ branches/tap_branches/ppStack/src/ppStackCombineInitial.c (revision 27838) @@ -10,4 +10,6 @@ #include "ppStackLoop.h" +//#define TESTING // Enable test output + bool ppStackCombineInitial(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config) { @@ -18,4 +20,5 @@ if (!options->convolve) { // No need to do initial combination when we haven't convolved + // XXX either allocate inspect and rejected here, or do not require them downstream return true; } @@ -34,5 +37,5 @@ if (!status) { // Something went wrong - psError(PS_ERR_UNKNOWN, false, "Unable to read chunk %d", numChunk); + psError(psErrorCodeLast(), false, "Unable to read chunk %d", numChunk); return false; } @@ -54,5 +57,5 @@ if (!psThreadPoolWait(false)) { - psError(PS_ERR_UNKNOWN, false, "Unable to do initial combination."); + psError(psErrorCodeLast(), false, "Unable to do initial combination."); return false; } @@ -60,4 +63,5 @@ // Harvest the jobs, gathering the inspection lists options->inspect = psArrayAlloc(options->num); + options->rejected = psArrayAlloc(options->num); for (int i = 0; i < options->num; i++) { if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) { @@ -65,4 +69,5 @@ } options->inspect->data[i] = psArrayAllocEmpty(numChunk); + options->rejected->data[i] = psArrayAllocEmpty(numChunk); } psThreadJob *job; // Completed job @@ -71,9 +76,13 @@ "Job has incorrect type: %s", job->type); psArray *results = job->results; // Results of job + psAssert(results->n == 2, "Results array has wrong size!"); + psArray *inspect = results->data[0]; // Pixels to inspect + psArray *reject = results->data[1]; // Pixels to reject for (int i = 0; i < options->num; i++) { if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) { continue; } - options->inspect->data[i] = psArrayAdd(options->inspect->data[i], 1, results->data[i]); + options->inspect->data[i] = psArrayAdd(options->inspect->data[i], 1, inspect->data[i]); + options->rejected->data[i] = psArrayAdd(options->rejected->data[i], 1, reject->data[i]); } psFree(job); @@ -83,10 +92,10 @@ #ifdef TESTING - ppStackWriteImage("combined_initial.fits", NULL, outRO->image, config); - pmStackVisualPlotTestImage(outRO->image, "combined_initial.fits"); + ppStackWriteImage("combined_image_initial.fits", NULL, options->outRO->image, config); + ppStackWriteImage("combined_mask_initial.fits", NULL, options->outRO->mask, config); + ppStackWriteImage("combined_variance_initial.fits", NULL, options->outRO->variance, config); + + pmStackVisualPlotTestImage(options->outRO->image, "combined_image_initial.fits"); #endif - - psFree(options->matchChi2); options->matchChi2 = NULL; - if (options->stats) { Index: branches/tap_branches/ppStack/src/ppStackCombinePrepare.c =================================================================== --- branches/tap_branches/ppStack/src/ppStackCombinePrepare.c (revision 25900) +++ branches/tap_branches/ppStack/src/ppStackCombinePrepare.c (revision 27838) @@ -10,5 +10,7 @@ #include "ppStackLoop.h" -bool ppStackCombinePrepare(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config) +bool ppStackCombinePrepare(const char *outName, const char *expName, + ppStackFileList files, ppStackThreadData *stack, + ppStackOptions *options, pmConfig *config) { psAssert(stack, "Require stack"); @@ -20,10 +22,10 @@ ppStackThread *thread = stack->threads->data[0]; // Representative thread if (!pmReadoutStackSetOutputSize(&col0, &row0, &numCols, &numRows, thread->readouts)) { - psError(PS_ERR_UNKNOWN, false, "problem setting output readout size."); + psError(PPSTACK_ERR_ARGUMENTS, false, "problem setting output readout size."); return false; } pmFPAfileActivate(config->files, false, NULL); - ppStackFileActivation(config, PPSTACK_FILES_COMBINE, true); + ppStackFileActivation(config, files, true); pmFPAview *view = ppStackFilesIterateDown(config); // View to readout if (!view) { @@ -31,6 +33,10 @@ } - pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT"); // Output cell - options->outRO = pmReadoutAlloc(outCell); // Output readout + pmCell *cell = pmFPAfileThisCell(config->files, view, outName); // Output cell + options->outRO = pmReadoutAlloc(cell); // Output readout + + pmCell *expCell = pmFPAfileThisCell(config->files, view, expName); // Exposure cell + options->expRO = pmReadoutAlloc(expCell); // Output readout + psFree(view); @@ -39,6 +45,12 @@ psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels + if (!pmReadoutStackDefineOutput(options->outRO, col0, row0, numCols, numRows, true, true, maskBad)) { - psError(PS_ERR_UNKNOWN, false, "Unable to prepare output."); + psError(PPSTACK_ERR_ARGUMENTS, false, "Unable to prepare output."); + return false; + } + + if (!pmReadoutStackDefineOutput(options->expRO, col0, row0, numCols, numRows, true, true, 0)) { + psError(PPSTACK_ERR_ARGUMENTS, false, "Unable to prepare output."); return false; } Index: branches/tap_branches/ppStack/src/ppStackConvolve.c =================================================================== --- branches/tap_branches/ppStack/src/ppStackConvolve.c (revision 25900) +++ branches/tap_branches/ppStack/src/ppStackConvolve.c (revision 27838) @@ -9,4 +9,7 @@ #include "ppStack.h" #include "ppStackLoop.h" + +//#define TESTING + // Update the value of a concept @@ -39,5 +42,9 @@ options->weightings = psVectorAlloc(num, PS_TYPE_F32); // Combination weightings for images (1/noise^2) psVectorInit(options->weightings, 0.0); - options->covariances = psArrayAlloc(num); // Covariance matrices + options->origCovars = psArrayAlloc(num); + options->convCovars = psArrayAlloc(num); // Covariance matrices + + psVector *renorms = psVectorAlloc(num, PS_TYPE_F32); // Renormalisation values for variances + psVectorInit(renorms, NAN); psList *fpaList = psListAlloc(NULL); // List of input FPAs, for concept averaging @@ -52,4 +59,10 @@ pmFPAfileActivate(config->files, false, NULL); ppStackFileActivationSingle(config, PPSTACK_FILES_CONVOLVE, true, i); + if (options->convolve) { + // XXX PPSTACK.CONV.KERNEL not defined unless convolve + // pmFPAfileActivate(config->files, true, "PPSTACK.CONV.KERNEL"); + pmFPAfileActivateSingle(config->files, true, "PPSTACK.CONV.KERNEL", i); // Activated file + } + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File of interest pmFPAview *view = ppStackFilesIterateDown(config); @@ -69,5 +82,5 @@ } else if (options->numCols != readout->image->numCols || options->numRows != readout->image->numRows) { - psError(PS_ERR_UNKNOWN, true, "Sizes of input images don't match: %dx%d vs %dx%d", + psError(PPSTACK_ERR_ARGUMENTS, true, "Sizes of input images don't match: %dx%d vs %dx%d", readout->image->numCols, readout->image->numRows, options->numCols, options->numRows); psFree(rng); @@ -79,11 +92,33 @@ // Background subtraction, scaling and normalisation is performed automatically by the image matching psTimerStart("PPSTACK_MATCH"); + options->origCovars->data[i] = psMemIncrRefCounter(readout->covariance); if (!ppStackMatch(readout, options, i, config)) { - psErrorStackPrint(stderr, "Unable to match image %d --- ignoring.", i); - options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PPSTACK_MASK_MATCH; - psErrorClear(); - continue; - } - options->covariances->data[i] = psMemIncrRefCounter(readout->covariance); + // XXX many things can cause a failure of ppStackMatch -- should some be handled differently? + psErrorCode error = psErrorCodeLast(); // Error code + switch (error) { + // Fatal errors + case PM_ERR_CONFIG: + case PPSTACK_ERR_CONFIG: + case PPSTACK_ERR_IO: + psError(error, false, "Unable to match image %d due to fatal error.", i); + return false; + // Non-fatal errors + case PM_ERR_STAMPS: + case PM_ERR_SMALL_AREA: + case PPSTACK_ERR_DATA: + default: + psErrorStackPrint(stderr, "Unable to match image %d --- ignoring.", i); + options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PPSTACK_MASK_MATCH; + psErrorClear(); + continue; + } + } + options->convCovars->data[i] = psMemIncrRefCounter(readout->covariance); + + float renorm = psMetadataLookupF32(NULL, readout->analysis, PM_READOUT_ANALYSIS_RENORM); + if (!isfinite(renorm)) { + renorm = 1.0; + } + renorms->data.F32[i] = renorm; if (options->stats) { @@ -114,11 +149,29 @@ pmHDU *hdu = readout->parent->parent->parent->hdu; // HDU for convolved image assert(hdu); - ppStackWriteImage(options->imageNames->data[i], hdu->header, readout->image, config); + if (!ppStackWriteImage(options->convImages->data[i], hdu->header, readout->image, config)) { + psError(PPSTACK_ERR_IO, false, "Unable to write convolved image %d", i); + psFree(fpaList); + psFree(cellList); + psFree(rng); + return false; + } psMetadata *maskHeader = psMetadataCopy(NULL, hdu->header); // Copy of header, for mask pmConfigMaskWriteHeader(config, maskHeader); - ppStackWriteImage(options->maskNames->data[i], maskHeader, readout->mask, config); + if (!ppStackWriteImage(options->convMasks->data[i], maskHeader, readout->mask, config)) { + psError(PPSTACK_ERR_IO, false, "Unable to write convolved mask %d", i); + psFree(fpaList); + psFree(cellList); + psFree(rng); + psFree(maskHeader); + return false; + } psFree(maskHeader); - psImageCovarianceTransfer(readout->variance, readout->covariance); - ppStackWriteImage(options->varianceNames->data[i], hdu->header, readout->variance, config); + if (!ppStackWriteImage(options->convVariances->data[i], hdu->header, readout->variance, config)) { + psError(PPSTACK_ERR_IO, false, "Unable to write convolved variance %d", i); + psFree(fpaList); + psFree(cellList); + psFree(rng); + return false; + } #ifdef TESTING { @@ -129,4 +182,19 @@ psFree(name); } + { + int numCols = readout->image->numCols, numRows = readout->image->numRows; + psImage *sn = psImageAlloc(numCols, numRows, PS_TYPE_F32); + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + sn->data.F32[y][x] = readout->image->data.F32[y][x] / + sqrtf(readout->variance->data.F32[y][x]); + } + } + psString name = NULL; + psStringAppend(&name, "signoise_%d.fits", i); + ppStackWriteImage(name, hdu->header, sn, config); + psFree(name); + psFree(sn); + } #endif @@ -136,4 +204,17 @@ psListAdd(fpaList, PS_LIST_TAIL, inCell->parent->parent); + // Correct ZP + if (options->matchZPs) { + // I think I need to take off the exposure time because we're going to set the new exposure time + psMetadataItem *zpItem = psMetadataLookup(inCell->parent->parent->concepts, "FPA.ZP"); + zpItem->data.F32 += options->norm->data.F32[i] + 2.5*log10(options->sumExposure); + + psMetadataItem *expItem = psMetadataLookup(inCell->parent->parent->concepts, "FPA.EXPOSURE"); + expItem->data.F32 = options->sumExposure; + + expItem = psMetadataLookup(inCell->concepts, "CELL.EXPOSURE"); + expItem->data.F32 = options->sumExposure; + } + options->cells->data[i] = psMemIncrRefCounter(inCell); if (!ppStackFilesIterateUp(config)) { @@ -149,5 +230,4 @@ psFree(rng); - psFree(options->norm); options->norm = NULL; psFree(options->sourceLists); options->sourceLists = NULL; psFree(options->psf); options->psf = NULL; @@ -155,5 +235,5 @@ if (numGood == 0) { - psError(PS_ERR_UNKNOWN, false, "No good images to combine."); + psError(PPSTACK_ERR_REJECTED, false, "No good images to combine."); psFree(fpaList); psFree(cellList); @@ -167,8 +247,17 @@ pmFPAview view; // View for output view.chip = view.cell = view.readout = 0; + pmCell *outCell = pmFPAfileThisCell(config->files, &view, "PPSTACK.OUTPUT"); // Output cell pmFPA *outFPA = outCell->parent->parent; // Output FPA + + pmCell *unconvCell = pmFPAfileThisCell(config->files, &view, "PPSTACK.UNCONV"); // Unconvolved cell + pmFPA *unconvFPA = unconvCell->parent->parent; // Unconvolved FPA + pmConceptsAverageFPAs(outFPA, fpaList); pmConceptsAverageCells(outCell, cellList, NULL, NULL, true); + + pmConceptsAverageFPAs(unconvFPA, fpaList); + pmConceptsAverageCells(unconvCell, cellList, NULL, NULL, true); + psFree(fpaList); psFree(cellList); @@ -177,4 +266,10 @@ UPDATE_CONCEPT(outCell, "CELL.EXPOSURE", options->sumExposure); UPDATE_CONCEPT(outCell, "CELL.DARKTIME", NAN); + UPDATE_CONCEPT(outFPA, "FPA.ZP", options->zp); + + UPDATE_CONCEPT(unconvFPA, "FPA.EXPOSURE", options->sumExposure); + UPDATE_CONCEPT(unconvCell, "CELL.EXPOSURE", options->sumExposure); + UPDATE_CONCEPT(unconvCell, "CELL.DARKTIME", NAN); + UPDATE_CONCEPT(unconvFPA, "FPA.ZP", options->zp); } @@ -191,5 +286,5 @@ assert(values->n == numGood); if (!psVectorSortInPlace(values)) { - psError(PS_ERR_UNKNOWN, false, "Unable to sort vector."); + psError(PPSTACK_ERR_PROG, false, "Unable to sort vector."); psFree(values); return false; @@ -211,5 +306,5 @@ numGood = 0; // Number of good images for (int i = 0; i < num; i++) { - if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_ALL) { + if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_ALL) { continue; } @@ -229,7 +324,15 @@ if (numGood == 0) { - psError(PS_ERR_UNKNOWN, false, "No good images to combine."); + psError(PPSTACK_ERR_REJECTED, false, "No good images to combine."); return false; } + + // Correct chi^2 for renormalisation + psBinaryOp(options->matchChi2, options->matchChi2, "/", renorms); + for (int i = 0; i < num; i++) { + psLogMsg("ppStack", PS_LOG_INFO, "Additional variance for image %d: %f\n", + i, options->matchChi2->data.F32[i]); + } + psFree(renorms); return true; Index: branches/tap_branches/ppStack/src/ppStackErrorCodes.c.in =================================================================== --- branches/tap_branches/ppStack/src/ppStackErrorCodes.c.in (revision 27838) +++ branches/tap_branches/ppStack/src/ppStackErrorCodes.c.in (revision 27838) @@ -0,0 +1,37 @@ +/** @file ppStackErrorCodes.c.in + * + * @brief + * + * @ingroup ppStack + * + * @author IfA + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-05 20:44:04 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#include "pslib.h" +#include "ppStackErrorCodes.h" + +/* + * The line + { PPSTACK_ERR_$X{ErrorCode}, "$X{ErrorDescription}"}, + * (without the Xs) + * will be replaced by values from errorCodes.dat + */ +void ppStackErrorRegister(void) +{ + static psErrorDescription errors[] = { + { PPSTACK_ERR_BASE, "First value we use; lower values belong to psLib" }, + { PPSTACK_ERR_${ErrorCode}, "${ErrorDescription}"}, + }; + static int nerror = PPSTACK_ERR_NERROR - PPSTACK_ERR_BASE; ///< number of values in enum + + for (int i = 0; i < nerror; i++) { + psErrorDescription *tmp = psAlloc(sizeof(psErrorDescription)); + *tmp = errors[i]; + psErrorRegister(tmp, 1); + psFree(tmp); /* it's on the internal list */ + } + nerror = 0; // don't register more than once +} Index: branches/tap_branches/ppStack/src/ppStackErrorCodes.dat =================================================================== --- branches/tap_branches/ppStack/src/ppStackErrorCodes.dat (revision 27838) +++ branches/tap_branches/ppStack/src/ppStackErrorCodes.dat (revision 27838) @@ -0,0 +1,13 @@ +# +# This file is used to generate ppStackErrorClasses.h +# +BASE = 13000 First value we use; lower values belong to psLib +UNKNOWN Unknown ppSub error code +NOT_IMPLEMENTED Desired feature is not yet implemented +ARGUMENTS Incorrect arguments +CONFIG Problem in configure files +IO Problem in FITS I/O +DATA Problem in data values +PSF Problem determining PSF +REJECTED All images rejected +PROG Programming error Index: branches/tap_branches/ppStack/src/ppStackErrorCodes.h.in =================================================================== --- branches/tap_branches/ppStack/src/ppStackErrorCodes.h.in (revision 27838) +++ branches/tap_branches/ppStack/src/ppStackErrorCodes.h.in (revision 27838) @@ -0,0 +1,30 @@ +/** @file ppStackErrorCodes.h.in + * + * @brief + * + * @ingroup ppStack + * + * @author IfA + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-05 20:44:04 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#if !defined(PPSTACK_ERROR_CODES_H) +#define PPSTACK_ERROR_CODES_H +/* + * The line + * PPSTACK_ERR_$X{ErrorCode}, + * (without the X) + * + * will be replaced by values from errorCodes.dat + */ +typedef enum { + PPSTACK_ERR_BASE = 14000, + PPSTACK_ERR_${ErrorCode}, + PPSTACK_ERR_NERROR +} ppStackErrorCode; + +void ppStackErrorRegister(void); + +#endif Index: branches/tap_branches/ppStack/src/ppStackFiles.c =================================================================== --- branches/tap_branches/ppStack/src/ppStackFiles.c (revision 25900) +++ branches/tap_branches/ppStack/src/ppStackFiles.c (revision 27838) @@ -17,10 +17,17 @@ /// Files required for the convolution -static char *filesConvolve[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.VARIANCE", - "PPSTACK.CONV.KERNEL", NULL }; - -/// Output files for the combination -static char *filesCombine[] = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.VARIANCE", - "PPSTACK.OUTPUT.JPEG1", "PPSTACK.OUTPUT.JPEG2", NULL }; +static char *filesConvolve[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.VARIANCE", NULL }; + +// "PPSTACK.CONV.KERNEL", NULL }; + +/// Regular (convolved) stack files +static char *filesStack[] = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.VARIANCE", + "PPSTACK.OUTPUT.EXP", "PPSTACK.OUTPUT.EXPNUM", "PPSTACK.OUTPUT.EXPWT", + "PPSTACK.OUTPUT.JPEG1", "PPSTACK.OUTPUT.JPEG2", + NULL }; +/// Unconvolved stack files +static char *filesUnconv[] = { "PPSTACK.UNCONV", "PPSTACK.UNCONV.MASK", "PPSTACK.UNCONV.VARIANCE", + "PPSTACK.UNCONV.EXP", "PPSTACK.UNCONV.EXPNUM", "PPSTACK.UNCONV.EXPWT", + NULL }; /// Files for photometry @@ -35,5 +42,6 @@ case PPSTACK_FILES_PREPARE: return filesPrepare; case PPSTACK_FILES_CONVOLVE: return filesConvolve; - case PPSTACK_FILES_COMBINE: return filesCombine; + case PPSTACK_FILES_STACK: return filesStack; + case PPSTACK_FILES_UNCONV: return filesUnconv; case PPSTACK_FILES_PHOT: return filesPhot; default: @@ -113,16 +121,20 @@ pmFPAview *view = pmFPAviewAlloc(0);// Pointer into FPA hierarchy if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(PPSTACK_ERR_IO, false, "File checks failed."); return NULL; } view->chip = 0; if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(PPSTACK_ERR_IO, false, "File checks failed."); return NULL; } view->cell = 0; if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(PPSTACK_ERR_IO, false, "File checks failed."); return NULL; } view->readout = 0; if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(PPSTACK_ERR_IO, false, "File checks failed."); return NULL; } @@ -139,16 +151,20 @@ view->chip = view->cell = view->readout = 0; if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(PPSTACK_ERR_IO, false, "File checks failed."); return false; } view->readout = -1; if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(PPSTACK_ERR_IO, false, "File checks failed."); return false; } view->cell = -1; if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(PPSTACK_ERR_IO, false, "File checks failed."); return false; } view->chip = -1; if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(PPSTACK_ERR_IO, false, "File checks failed."); return false; } @@ -170,15 +186,19 @@ psFits *fits = psFitsOpen(resolved, "w"); if (!fits) { - psError(PS_ERR_IO, false, "Unable to open FITS file %s to write image.", resolved); + psError(PPSTACK_ERR_IO, false, "Unable to open FITS file %s to write image.", resolved); psFree(resolved); return false; } if (!psFitsWriteImage(fits, header, image, 0, NULL)) { - psError(PS_ERR_IO, false, "Unable to write FITS image %s.", resolved); + psError(PPSTACK_ERR_IO, false, "Unable to write FITS image %s.", resolved); psFitsClose(fits); psFree(resolved); return false; } - psFitsClose(fits); + if (!psFitsClose(fits)) { + psError(PPSTACK_ERR_IO, false, "Unable to close FITS image %s.", resolved); + psFree(resolved); + return false; + } psFree(resolved); return true; Index: branches/tap_branches/ppStack/src/ppStackFinish.c =================================================================== --- branches/tap_branches/ppStack/src/ppStackFinish.c (revision 25900) +++ branches/tap_branches/ppStack/src/ppStackFinish.c (revision 27838) @@ -4,7 +4,8 @@ #include +#include #include #include -#include +#include #include "ppStack.h" @@ -21,55 +22,92 @@ psAssert(recipe, "We've thrown an error on this before."); - // Statistics on output - if (options->stats) { - psTrace("ppStack", 1, "Gathering statistics on stacked image....\n"); - psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits for bad - psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels + bool mdok; // Status of MD lookup + bool tempDelete = psMetadataLookupBool(&mdok, recipe, "TEMP.DELETE"); // Delete temporary files? - pmFPAview *view = pmFPAviewAlloc(0); // View to readout - view->chip = view->cell = view->readout = 0; + // Delete temporary images + if (tempDelete && options->convolve) { + for (int i = 0; i < options->num; i++) { + if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) { + continue; + } - ppStatsFPA(options->stats, options->outRO->parent->parent->parent, view, maskBad, config); - - psFree(view); - } - - ppStackMemDump("stats"); - - psFree(options->outRO); options->outRO = NULL; - - // Write out the output files - if (!ppStackFilesIterateUp(config)) { - return false; + psString imageResolved = pmConfigConvertFilename(options->convImages->data[i], + config, false, false); + psString maskResolved = pmConfigConvertFilename(options->convMasks->data[i], + config, false, false); + psString varianceResolved = pmConfigConvertFilename(options->convVariances->data[i], + config, false, false); + if (unlink(imageResolved) == -1 || unlink(maskResolved) == -1 || + unlink(varianceResolved) == -1) { + psWarning("Unable to delete temporary files for image %d", i); + } + psFree(imageResolved); + psFree(maskResolved); + psFree(varianceResolved); + } } - // Write out summary statistics - if (options->stats) { - psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_STACK", 0, - "Total time", psTimerClear("PPSTACK_TOTAL")); + return true; +} - const char *statsMDC = psMetadataConfigFormat(options->stats); - if (!statsMDC || strlen(statsMDC) == 0) { - psError(PS_ERR_IO, false, "Unable to get statistics MDC file.\n"); - } else { - fprintf(options->statsFile, "%s", statsMDC); - } - psFree((void *)statsMDC); - fclose(options->statsFile); options->statsFile = NULL; - pmConfigRunFilenameAddWrite(config, "STATS", psMetadataLookupStr(NULL, config->arguments, "STATS")); + +psExit ppStackExitCode(psExit exitValue) +{ + if (exitValue != PS_EXIT_SUCCESS) { + return exitValue; } - - // Dump configuration - bool mdok; // Status of MD lookup - psString dump = psMetadataLookupStr(&mdok, config->arguments, "DUMP_CONFIG"); // File for config - if (dump) { - if (!pmConfigDump(config, dump)) { - psError(PS_ERR_IO, false, "Unable to dump configuration."); - return false; + psErrorCode errorCode = psErrorCodeLast(); // Error code + if (errorCode != PS_ERR_NONE) { + psErrorStackPrint(stderr, "Unable to perform stack."); + switch (errorCode) { + case PPSTACK_ERR_UNKNOWN: + case PS_ERR_UNKNOWN: + psLogMsg("ppStack", 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 PPSTACK_ERR_IO: + psLogMsg("ppStack", 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 PPSTACK_ERR_ARGUMENTS: + case PPSTACK_ERR_CONFIG: + psLogMsg("ppStack", PS_LOG_WARN, "Configuration error code: %x", errorCode); + exitValue = PS_EXIT_CONFIG_ERROR; + break; + case PPSTACK_ERR_PSF: + case PSPHOT_ERR_PSF: + case PM_ERR_STAMPS: + case PM_ERR_SMALL_AREA: + case PPSTACK_ERR_REJECTED: + case PPSTACK_ERR_DATA: + psLogMsg("ppStack", PS_LOG_WARN, "Data error code: %x", errorCode); + exitValue = PS_EXIT_DATA_ERROR; + break; + case PS_ERR_UNEXPECTED_NULL: + case PS_ERR_PROGRAMMING: + case PPSTACK_ERR_NOT_IMPLEMENTED: + case PPSTACK_ERR_PROG: + psLogMsg("ppStack", 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("ppStack", PS_LOG_WARN, "Unrecognised error code: %x", errorCode); + exitValue = PS_EXIT_PROG_ERROR; + break; } } - - return true; + return exitValue; } Index: branches/tap_branches/ppStack/src/ppStackLoop.c =================================================================== --- branches/tap_branches/ppStack/src/ppStackLoop.c (revision 25900) +++ branches/tap_branches/ppStack/src/ppStackLoop.c (revision 27838) @@ -10,5 +10,45 @@ #include "ppStackLoop.h" -bool ppStackLoop(pmConfig *config) +/// Print a summary of the inputs, and return the number of good inputs +static int stackSummary(const ppStackOptions *options, // Stack options, with input mask + const char *place // Place in code + ) +{ + int numGood = 0; // Number of good inputs + psString summary = NULL; // Summary of images + for (int i = 0; i < options->num; i++) { + psString reason = NULL; // Reason for rejecting + if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] == 0) { + psStringAppend(&reason, " Good."); + numGood++; + } else { + if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_CAL) { + psStringAppend(&reason, " Calibration failed."); + } + if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_PSF) { + psStringAppend(&reason, " PSF measurement failed."); + } + if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_MATCH) { + psStringAppend(&reason, " PSF matching failed."); + } + if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_CHI2) { + psStringAppend(&reason, " PSF matching chi^2 deviant."); + } + if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_REJECT) { + psStringAppend(&reason, " Rejection exceeded threshold."); + } + } + psStringAppend(&summary, "Image %d: %s\n", i, reason); + psFree(reason); + } + psLogMsg("ppStack", PS_LOG_INFO, "Summary of images for %s:\n%s", place, summary); + psFree(summary); + + return numGood; +} + + + +bool ppStackLoop(pmConfig *config, ppStackOptions *options) { assert(config); @@ -16,12 +56,9 @@ psTimerStart("PPSTACK_TOTAL"); psTimerStart("PPSTACK_STEPS"); - - ppStackOptions *options = ppStackOptionsAlloc(); // Options for stacking // Setup psTrace("ppStack", 1, "Setup....\n"); if (!ppStackSetup(options, config)) { - psError(PS_ERR_UNKNOWN, false, "Unable to setup."); - psFree(options); + psError(psErrorCodeLast(), false, "Unable to setup."); return false; } @@ -33,6 +70,5 @@ psTrace("ppStack", 1, "Preparation for stacking: merging sources, determining target PSF....\n"); if (!ppStackPrepare(options, config)) { - psError(PS_ERR_UNKNOWN, false, "Unable to prepare for stacking."); - psFree(options); + psError(psErrorCodeLast(), false, "Unable to prepare for stacking."); return false; } @@ -40,11 +76,13 @@ psTimerClear("PPSTACK_STEPS")); ppStackMemDump("prepare"); - + if (options->quality) { + // Can't do anything else + return true; + } // Convolve inputs psTrace("ppStack", 1, "Convolving inputs to target PSF....\n"); if (!ppStackConvolve(options, config)) { - psError(PS_ERR_UNKNOWN, false, "Unable to convolve images."); - psFree(options); + psError(psErrorCodeLast(), false, "Unable to convolve images."); return false; } @@ -53,20 +91,28 @@ ppStackMemDump("convolve"); + // Ensure sufficient inputs + { + int numGood = stackSummary(options, "initial combination"); + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe + bool safe = psMetadataLookupBool(NULL, recipe, "SAFE"); // Be safe when combining + if (safe && numGood <= 1) { + psError(PPSTACK_ERR_REJECTED, true, "Insufficient inputs for combination with safety on"); + return false; + } + } // Start threading ppStackThreadInit(); - ppStackThreadData *stack = ppStackThreadDataSetup(options, config); + ppStackThreadData *stack = ppStackThreadDataSetup(options, config, true); if (!stack) { - psError(PS_ERR_IO, false, "Unable to initialise stack threads."); - psFree(options); - return false; - } - psFree(options->cells); options->cells = NULL; + psError(psErrorCodeLast(), false, "Unable to initialise stack threads."); + return false; + } // Prepare for combination - if (!ppStackCombinePrepare(stack, options, config)) { - psError(PS_ERR_UNKNOWN, false, "Unable to prepare for combination."); - psFree(stack); - psFree(options); + if (!ppStackCombinePrepare("PPSTACK.OUTPUT", "PPSTACK.OUTPUT.EXP", PPSTACK_FILES_STACK, + stack, options, config)) { + psError(psErrorCodeLast(), false, "Unable to prepare for combination."); + psFree(stack); return false; } @@ -75,31 +121,46 @@ psTrace("ppStack", 1, "Initial stack of convolved images....\n"); if (!ppStackCombineInitial(stack, options, config)) { - psError(PS_ERR_UNKNOWN, false, "Unable to perform initial combination."); - psFree(stack); - psFree(options); - return false; - } - ppStackMemDump("convolve"); + psError(psErrorCodeLast(), false, "Unable to perform initial combination."); + psFree(stack); + return false; + } + ppStackMemDump("initial"); psLogMsg("ppStack", PS_LOG_INFO, "Stage 3: Make Initial Stack: %f sec", psTimerClear("PPSTACK_STEPS")); + // Done with stack inputs for now + for (int i = 0; i < options->num; i++) { + pmCellFreeData(options->cells->data[i]); + } + psFree(stack); // Pixel rejection psTrace("ppStack", 1, "Reject pixels....\n"); if (!ppStackReject(options, config)) { - psError(PS_ERR_UNKNOWN, false, "Unable to reject pixels."); - psFree(stack); - psFree(options); + psError(psErrorCodeLast(), false, "Unable to reject pixels."); + psFree(stack); return false; } psLogMsg("ppStack", PS_LOG_INFO, "Stage 4: Pixel Rejection: %f sec", psTimerClear("PPSTACK_STEPS")); - ppStackMemDump("reject"); - + + // Check inputs + { + int numGood = stackSummary(options, "final combination"); + if (numGood <= 0) { + psError(PPSTACK_ERR_REJECTED, true, "Insufficient inputs for combination"); + return false; + } + } + + stack = ppStackThreadDataSetup(options, config, true); + if (!stack) { + psError(psErrorCodeLast(), false, "Unable to initialise stack threads."); + return false; + } // Final combination psTrace("ppStack", 2, "Final stack of convolved images....\n"); - if (!ppStackCombineFinal(stack, options, config)) { - psError(PS_ERR_UNKNOWN, false, "Unable to perform final combination."); - psFree(stack); - psFree(options); + if (!ppStackCombineFinal(stack, options->convCovars, options, config, false, false, true)) { + psError(psErrorCodeLast(), false, "Unable to perform final combination."); + psFree(stack); return false; } @@ -107,41 +168,88 @@ ppStackMemDump("final"); + // Photometry + psTrace("ppStack", 1, "Photometering stacked image....\n"); + if (!ppStackPhotometry(options, config)) { + psError(psErrorCodeLast(), false, "Unable to perform photometry."); + return false; + } + psLogMsg("ppStack", PS_LOG_INFO, "Stage 6: Photometry Analysis: %f sec", psTimerClear("PPSTACK_STEPS")); + ppStackMemDump("photometry"); // Clean up psTrace("ppStack", 2, "Cleaning up after combination....\n"); if (!ppStackCleanup(stack, options, config)) { - psError(PS_ERR_UNKNOWN, false, "Unable to clean up."); - psFree(stack); - psFree(options); - return false; - } - psLogMsg("ppStack", PS_LOG_INFO, "Stage 6: WCS & JPEGS: %f sec", psTimerClear("PPSTACK_STEPS")); + psError(psErrorCodeLast(), false, "Unable to clean up."); + psFree(stack); + return false; + } + for (int i = 0; i < options->num; i++) { + pmCellFreeData(options->cells->data[i]); + } + psFree(stack); + psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: Cleanup, WCS & JPEGS: %f sec", psTimerClear("PPSTACK_STEPS")); ppStackMemDump("cleanup"); - psFree(stack); - - - // Photometry - psTrace("ppStack", 1, "Photometering stacked image....\n"); - if (!ppStackPhotometry(options, config)) { - psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry."); - psFree(options); - return false; - } - psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: Photometry Analysis: %f sec", psTimerClear("PPSTACK_STEPS")); - ppStackMemDump("photometry"); - +#if 1 + // Unconvolved stack --- it's cheap to calculate, compared to everything else! + if (options->convolve) { + // Start threading + ppStackThreadData *stack = ppStackThreadDataSetup(options, config, false); + if (!stack) { + psError(psErrorCodeLast(), false, "Unable to initialise stack threads."); + return false; + } + + // Prepare for combination + if (!ppStackCombinePrepare("PPSTACK.UNCONV", "PPSTACK.UNCONV.EXP", PPSTACK_FILES_UNCONV, + stack, options, config)) { + psError(psErrorCodeLast(), false, "Unable to prepare for combination."); + psFree(stack); + return false; + } + + psTrace("ppStack", 2, "Stack of unconvolved images....\n"); + if (!ppStackCombineFinal(stack, options->origCovars, options, config, + false, true, false)) { + psError(psErrorCodeLast(), false, "Unable to perform unconvolved combination."); + psFree(stack); + return false; + } + psLogMsg("ppStack", PS_LOG_INFO, "Stage 8: Unconvolved Stack: %f sec", psTimerClear("PPSTACK_STEPS")); + ppStackMemDump("unconv"); + + if (!ppStackFilesIterateUp(config)) { + psError(psErrorCodeLast(), false, "Unable to close files."); + psFree(stack); + return false; + } + ppStackFileActivation(config, PPSTACK_FILES_UNCONV, false); + options->outRO->data_exists = false; + options->outRO->parent->data_exists = false; + options->outRO->parent->parent->data_exists = false; + psFree(options->outRO); + options->outRO = NULL; + options->expRO->data_exists = false; + options->expRO->parent->data_exists = false; + options->expRO->parent->parent->data_exists = false; + psFree(options->expRO); + options->expRO = NULL; + + for (int i = 0; i < options->num; i++) { + pmCellFreeData(options->cells->data[i]); + } + psFree(stack); + } + psFree(options->cells); options->cells = NULL; +#endif // Finish up psTrace("ppStack", 1, "Finishing up....\n"); if (!ppStackFinish(options, config)) { - psError(PS_ERR_UNKNOWN, false, "Unable to finish up."); - psFree(options); - return false; - } - psLogMsg("ppStack", PS_LOG_INFO, "Stage 8: Final output: %f sec", psTimerClear("PPSTACK_STEPS")); + psError(psErrorCodeLast(), false, "Unable to finish up."); + return false; + } ppStackMemDump("finish"); - psFree(options); return true; } Index: branches/tap_branches/ppStack/src/ppStackLoop.h =================================================================== --- branches/tap_branches/ppStack/src/ppStackLoop.h (revision 25900) +++ branches/tap_branches/ppStack/src/ppStackLoop.h (revision 27838) @@ -7,9 +7,11 @@ #include "ppStackOptions.h" #include "ppStackThread.h" +#include "ppStack.h" // Loop over the inputs, doing the combination bool ppStackLoop( - pmConfig *config // Configuration + pmConfig *config, // Configuration + ppStackOptions *options // Options for stacking ); @@ -36,4 +38,7 @@ // Prepare for combination bool ppStackCombinePrepare( + const char *outName, // Name of output file + const char *expName, // Name of exposure file + ppStackFileList files, // Files of interest ppStackThreadData *stack, // Stack ppStackOptions *options, // Options @@ -57,6 +62,10 @@ bool ppStackCombineFinal( ppStackThreadData *stack, // Stack + psArray *covariances, // Covariances ppStackOptions *options, // Options - pmConfig *config // Configuration + pmConfig *config, // Configuration + bool safe, // Allow safe combination? + bool norm, // Normalise images? + bool grow // Grow rejection masks? ); Index: branches/tap_branches/ppStack/src/ppStackMatch.c =================================================================== --- branches/tap_branches/ppStack/src/ppStackMatch.c (revision 25900) +++ branches/tap_branches/ppStack/src/ppStackMatch.c (revision 27838) @@ -18,5 +18,5 @@ #define COVAR_FRAC 0.01 // Truncation fraction for covariance matrix -//#define TESTING // Enable debugging output +// #define TESTING // Enable debugging output #ifdef TESTING @@ -31,10 +31,10 @@ psFree(resolved); if (!fits) { - psError(PS_ERR_IO, false, "Unable to open previously produced image: %s", name); + psError(PPSTACK_ERR_IO, false, "Unable to open previously produced image: %s", name); return false; } psImage *image = psFitsReadImage(fits, psRegionSet(0,0,0,0), 0); // Image of interest if (!image) { - psError(PS_ERR_IO, false, "Unable to read previously produced image: %s", name); + psError(PPSTACK_ERR_IO, false, "Unable to read previously produced image: %s", name); psFitsClose(fits); return false; @@ -115,4 +115,6 @@ psFree(coords); psFree(tree); + psFree(x); + psFree(y); psLogMsg("ppStack", PS_LOG_INFO, "Filtered out %d of %d sources", numFiltered, numGood); @@ -144,5 +146,5 @@ psMetadataAddU8(psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "user-defined mask", maskBad); - psImage *binned = psphotBackgroundModel(ro, config); // Binned background model + psImage *binned = psphotModelBackgroundReadoutNoFile(ro, config); // Binned background model psImageBinning *binning = psMetadataLookupPtr(NULL, ro->analysis, "PSPHOT.BACKGROUND.BINNING"); // Binning for model @@ -150,14 +152,12 @@ psImage *unbinned = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Unbinned background model if (!psImageUnbin(unbinned, binned, binning)) { - psError(PS_ERR_UNKNOWN, false, "Unable to unbin background model"); + psError(PPSTACK_ERR_DATA, false, "Unable to unbin background model"); + psFree(binned); psFree(unbinned); return NULL; } - - // XXX should these really be here?? (probably not...) - // pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL"); - // pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL.STDEV"); - - return unbinned; + psFree(binned); + + return unbinned; } @@ -167,4 +167,5 @@ ) { +#if 1 bool mdok; // Status of metadata lookups @@ -176,15 +177,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(PPSTACK_ERR_CONFIG, 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(PPSTACK_ERR_CONFIG, 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(PPSTACK_ERR_CONFIG, true, "RENORM.MAX is not set in the recipe"); return false; } @@ -192,5 +193,9 @@ psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask + psImageCovarianceTransfer(readout->variance, readout->covariance); return pmReadoutVarianceRenormalise(readout, maskBad, num, minValid, maxValid); +#else + return true; +#endif } @@ -212,5 +217,5 @@ int size = psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Kernel half-size - psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in + psString maskValStr = psMetadataLookupStr(NULL, ppsub, "MASK.VAL"); // Name of bits to mask going in psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch psString maskPoorStr = psMetadataLookupStr(NULL, recipe, "MASK.POOR"); // Name of bits to mask for poor @@ -224,5 +229,5 @@ if (!pmReadoutMaskNonfinite(readout, maskVal)) { - psError(PS_ERR_UNKNOWN, false, "Unable to mask non-finite pixels in readout."); + psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels in readout."); return false; } @@ -251,5 +256,5 @@ psFree(resolved); if (!fits || !pmReadoutReadSubtractionKernels(conv, fits)) { - psError(PS_ERR_IO, false, "Unable to read previously produced kernel"); + psError(PPSTACK_ERR_IO, false, "Unable to read previously produced kernel"); psFitsClose(fits); return false; @@ -257,8 +262,8 @@ psFitsClose(fits); - if (!readImage(&readout->image, options->imageNames->data[index], config) || - !readImage(&readout->mask, options->maskNames->data[index], config) || - !readImage(&readout->variance, options->varianceNames->data[index], config)) { - psError(PS_ERR_IO, false, "Unable to read previously produced image."); + if (!readImage(&readout->image, options->convImages->data[index], config) || + !readImage(&readout->mask, options->convMasks->data[index], config) || + !readImage(&readout->variance, options->convVariances->data[index], config)) { + psError(PPSTACK_ERR_IO, false, "Unable to read previously produced image."); return false; } @@ -269,9 +274,11 @@ PM_SUBTRACTION_ANALYSIS_KERNEL); - pmSubtractionAnalysis(readout->analysis, NULL, kernels, region, + pmSubtractionAnalysis(conv->analysis, NULL, kernels, region, readout->image->numCols, readout->image->numRows); psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel + bool oldThreads = psImageCovarianceSetThreads(true); // Old thread setting psKernel *covar = psImageCovarianceCalculate(kernel, readout->covariance); // Covariance matrix + psImageCovarianceSetThreads(oldThreads); psFree(readout->covariance); readout->covariance = covar; @@ -292,5 +299,10 @@ int iter = psMetadataLookupS32(NULL, ppsub, "ITER"); // Rejection iterations float rej = psMetadataLookupF32(NULL, ppsub, "REJ"); // Rejection threshold - float sysError = psMetadataLookupF32(NULL, ppsub, "SYS"); // Relative systematic error in kernel + float kernelError = psMetadataLookupF32(NULL, ppsub, "KERNEL.ERR"); // Relative systematic error in kernel + float normFrac = psMetadataLookupF32(NULL, ppsub, "NORM.FRAC"); // Fraction of window for normalisn windw + float sysError = psMetadataLookupF32(NULL, ppsub, "SYS.ERR"); // Relative systematic error in images + float skyErr = psMetadataLookupF32(NULL, ppsub, "SKY.ERR"); // Additional error in sky + float covarFrac = psMetadataLookupF32(NULL, ppsub, "COVAR.FRAC"); // Fraction for covariance calculation + const char *typeStr = psMetadataLookupStr(NULL, ppsub, "KERNEL.TYPE"); // Kernel type pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Kernel type @@ -309,4 +321,16 @@ float poorFrac = psMetadataLookupF32(&mdok, ppsub, "POOR.FRACTION"); // Fraction for "poor" + bool scale = psMetadataLookupBool(NULL, ppsub, "SCALE"); // Scale kernel parameters? + float scaleRef = psMetadataLookupF32(NULL, ppsub, "SCALE.REF"); // Reference for scaling + float scaleMin = psMetadataLookupF32(NULL, ppsub, "SCALE.MIN"); // Minimum for scaling + float scaleMax = psMetadataLookupF32(NULL, ppsub, "SCALE.MAX"); // Maximum for scaling + if (!isfinite(scaleRef) || !isfinite(scaleMin) || !isfinite(scaleMax)) { + psError(PPSTACK_ERR_CONFIG, false, + "Scale parameters (SCALE.REF=%f, SCALE.MIN=%f, SCALE.MAX=%f) not set in PPSUB recipe.", + scaleRef, scaleMin, scaleMax); + return false; + } + + // These values are specified specifically for stacking const char *stampsName = psMetadataLookupStr(NULL, config->arguments, "STAMPS");// Stamps filename @@ -322,5 +346,5 @@ psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) { - psError(PS_ERR_UNKNOWN, false, "Can't measure background for image."); + psError(PPSTACK_ERR_DATA, false, "Can't measure background for image."); psFree(fake); psFree(optWidths); @@ -338,8 +362,9 @@ footprint); // Filtered list of sources + bool oldThreads = pmReadoutFakeThreads(true); // Old threading state if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows, stampSources, SOURCE_MASK, NULL, NULL, options->psf, minFlux, footprint + size, false, true)) { - psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image with target PSF."); + psError(PPSTACK_ERR_DATA, false, "Unable to generate fake image with target PSF."); psFree(fake); psFree(optWidths); @@ -347,11 +372,14 @@ return false; } + pmReadoutFakeThreads(oldThreads); fake->mask = psImageCopy(NULL, readout->mask, PS_TYPE_IMAGE_MASK); +#if 1 // Add the background into the target image psImage *bgImage = stackBackgroundModel(readout, config); // Image of background psBinaryOp(fake->image, fake->image, "+", bgImage); psFree(bgImage); +#endif #ifdef TESTING @@ -379,5 +407,5 @@ if (threads > 0) { - pmSubtractionThreadsInit(readout, fake); + pmSubtractionThreadsInit(); } @@ -387,7 +415,7 @@ if (kernel) { if (!pmSubtractionMatchPrecalc(NULL, conv, fake, readout, readout->analysis, - stride, sysError, maskVal, maskBad, maskPoor, + stride, kernelError, covarFrac, maskVal, maskBad, maskPoor, poorFrac, badFrac)) { - psError(PS_ERR_UNKNOWN, false, "Unable to convolve images."); + psError(psErrorCodeLast(), false, "Unable to convolve images."); psFree(fake); psFree(optWidths); @@ -395,26 +423,46 @@ psFree(conv); if (threads > 0) { - pmSubtractionThreadsFinalize(readout, fake); + pmSubtractionThreadsFinalize(); } return false; } } else { - if (!pmSubtractionMatch(NULL, conv, fake, readout, footprint, stride, regionSize, spacing, - threshold, stampSources, stampsName, type, size, order, widths, - orders, inner, ringsOrder, binning, penalty, - optimum, optWidths, optOrder, optThresh, iter, rej, sysError, - maskVal, maskBad, maskPoor, poorFrac, badFrac, - PM_SUBTRACTION_MODE_2)) { - psError(PS_ERR_UNKNOWN, false, "Unable to match images."); + // Scale the input parameters + psVector *widthsCopy = psVectorCopy(NULL, widths, PS_TYPE_F32); // Copy of kernel widths + if (scale && !pmSubtractionParamsScale(&size, &footprint, widthsCopy, + options->inputSeeing->data.F32[index], + options->targetSeeing, scaleRef, scaleMin, scaleMax)) { + psError(psErrorCodeLast(), false, "Unable to scale kernel parameters"); psFree(fake); psFree(optWidths); psFree(stampSources); psFree(conv); + psFree(widthsCopy); if (threads > 0) { - pmSubtractionThreadsFinalize(readout, fake); + pmSubtractionThreadsFinalize(); } return false; } - } + + if (!pmSubtractionMatch(NULL, conv, fake, readout, footprint, stride, regionSize, spacing, + threshold, stampSources, stampsName, type, size, order, widthsCopy, + orders, inner, ringsOrder, binning, penalty, + optimum, optWidths, optOrder, optThresh, iter, rej, normFrac, + sysError, skyErr, kernelError, covarFrac, maskVal, maskBad, maskPoor, + poorFrac, badFrac, PM_SUBTRACTION_MODE_2)) { + psError(psErrorCodeLast(), false, "Unable to match images."); + psFree(fake); + psFree(optWidths); + psFree(stampSources); + psFree(conv); + psFree(widthsCopy); + if (threads > 0) { + pmSubtractionThreadsFinalize(); + } + return false; + } + psFree(widthsCopy); + } + #ifdef TESTING @@ -447,5 +495,5 @@ if (threads > 0) { - pmSubtractionThreadsFinalize(readout, fake); + pmSubtractionThreadsFinalize(); } @@ -487,10 +535,5 @@ while ((item = psMetadataGetAndIncrement(iter))) { assert(item->type == PS_DATA_UNKNOWN); - // Set the normalisation dimensions, since these will be otherwise unavailable when reading - // the images by scans. pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction - kernel->numCols = readout->image->numCols; - kernel->numRows = readout->image->numRows; - kernels = psArrayAdd(kernels, ARRAY_BUFFER, kernel); } @@ -518,10 +561,35 @@ } + // Kernel normalisation + { + double sum = 0.0; // Sum of chi^2 + int num = 0; // Number of measurements of chi^2 + psString regex = NULL; // Regular expression + psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_NORM); + psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex); + psFree(regex); + psMetadataItem *item = NULL;// Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + assert(item->type == PS_TYPE_F32); + float norm = item->data.F32; // Normalisation + sum += norm; + num++; + } + psFree(iter); + float conv = sum/num; // Mean normalisation from convolution + float stars = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation from stars + float renorm = stars / conv; // Renormalisation to apply + psLogMsg("ppStack", PS_LOG_INFO, "Renormalising image %d by %f (kernel: %f, stars: %f)\n", + index, renorm, conv, stars); + psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(renorm, PS_TYPE_F32)); + psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(renorm), PS_TYPE_F32)); + } + // Reject image completely if the maximum deconvolution fraction exceeds the limit float deconv = psMetadataLookupF32(NULL, conv->analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Max deconvolution fraction if (deconv > deconvLimit) { - psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting\n", - deconv, deconvLimit); + psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting image %d\n", + deconv, deconvLimit, index); psFree(conv); return NULL; @@ -542,13 +610,13 @@ psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) { - psWarning("Can't measure background for image."); - psErrorClear(); + psWarning("Can't measure background for image."); + psErrorClear(); } else { - if (!psMetadataLookupBool(NULL, config->arguments, "PPSTACK.SKIP.BG.SUB")) { - psLogMsg("ppStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)", - psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), psStatsGetValue(bg, PS_STAT_ROBUST_STDEV)); - (void)psBinaryOp(readout->image, readout->image, "-", - psScalarAlloc(psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), PS_TYPE_F32)); - } + if (!psMetadataLookupBool(NULL, config->arguments, "PPSTACK.SKIP.BG.SUB")) { + psLogMsg("ppStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)", + psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), psStatsGetValue(bg, PS_STAT_ROBUST_STDEV)); + (void)psBinaryOp(readout->image, readout->image, "-", + psScalarAlloc(psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), PS_TYPE_F32)); + } } @@ -560,14 +628,18 @@ // Measure the variance level for the weighting - if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) { - psError(PS_ERR_UNKNOWN, false, "Can't measure mean variance for image."); - psFree(rng); - psFree(bg); - return false; - } - options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) * - psImageCovarianceFactor(readout->covariance)); - psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "PPSTACK.WEIGHTING", 0, - "Weighting by 1/noise^2 for stack", options->weightings->data.F32[index]); + if (psMetadataLookupBool(NULL, recipe, "WEIGHTS")) { + if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) { + psError(PPSTACK_ERR_DATA, false, "Can't measure mean variance for image."); + psFree(rng); + psFree(bg); + return false; + } + options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) * + psImageCovarianceFactor(readout->covariance)); + } else { + options->weightings->data.F32[index] = 1.0; + } + psLogMsg("ppStack", PS_LOG_INFO, "Weighting for image %d is %f\n", + index, options->weightings->data.F32[index]); psFree(rng); Index: branches/tap_branches/ppStack/src/ppStackOptions.c =================================================================== --- branches/tap_branches/ppStack/src/ppStackOptions.c (revision 25900) +++ branches/tap_branches/ppStack/src/ppStackOptions.c (revision 27838) @@ -12,11 +12,18 @@ fclose(options->statsFile); } - psFree(options->imageNames); - psFree(options->maskNames); - psFree(options->varianceNames); + psFree(options->origImages); + psFree(options->origMasks); + psFree(options->origVariances); + psFree(options->origCovars); + psFree(options->convImages); + psFree(options->convMasks); + psFree(options->convVariances); psFree(options->psf); + psFree(options->inputSeeing); + psFree(options->exposures); psFree(options->inputMask); psFree(options->sourceLists); psFree(options->norm); + psFree(options->sources); psFree(options->cells); psFree(options->kernels); @@ -24,6 +31,7 @@ psFree(options->matchChi2); psFree(options->weightings); - psFree(options->covariances); + psFree(options->convCovars); psFree(options->outRO); + psFree(options->expRO); psFree(options->inspect); psFree(options->rejected); @@ -37,14 +45,27 @@ options->convolve = true; + options->matchZPs = true; + options->photometry = false; options->stats = NULL; options->statsFile = NULL; - options->imageNames = NULL; - options->maskNames = NULL; - options->varianceNames = NULL; + options->origImages = NULL; + options->origMasks = NULL; + options->origVariances = NULL; + options->origCovars = NULL; + options->convImages = NULL; + options->convMasks = NULL; + options->convVariances = NULL; options->num = 0; + options->quality = 0; options->psf = NULL; + options->sumExposure = NAN; + options->zp = NAN; + options->inputSeeing = NULL; + options->exposures = NULL; + options->targetSeeing = NAN; options->inputMask = NULL; options->sourceLists = NULL; options->norm = NULL; + options->sources = NULL; options->cells = NULL; options->kernels = NULL; @@ -54,6 +75,7 @@ options->matchChi2 = NULL; options->weightings = NULL; - options->covariances = NULL; + options->convCovars = NULL; options->outRO = NULL; + options->expRO = NULL; options->inspect = NULL; options->rejected = NULL; Index: branches/tap_branches/ppStack/src/ppStackOptions.h =================================================================== --- branches/tap_branches/ppStack/src/ppStackOptions.h (revision 25900) +++ branches/tap_branches/ppStack/src/ppStackOptions.h (revision 27838) @@ -9,14 +9,24 @@ // Setup bool convolve; // Convolve images? + bool matchZPs; // Adjust relative fluxes based on transparency analysis? + bool photometry; // Perform photometry? psMetadata *stats; // Statistics for output FILE *statsFile; // File to which to write statistics - psArray *imageNames, *maskNames, *varianceNames; // Filenames for the temporary convolved images + psArray *origImages, *origMasks, *origVariances; // Filenames of the original images + psArray *convImages, *convMasks, *convVariances; // Filenames for the temporary convolved images + psArray *origCovars; // Original covariances matrices int num; // Number of inputs + int quality; // Bad data quality flag // Prepare pmPSF *psf; // Target PSF + psVector *inputSeeing; // Input seeing FWHMs + float targetSeeing; // Target seeing FWHM + psVector *exposures; // Exposure times float sumExposure; // Sum of exposure times + float zp; // Zero point for output psVector *inputMask; // Mask for inputs psArray *sourceLists; // Individual lists of sources for matching psVector *norm; // Normalisation for each image + psArray *sources; // Matched sources // Convolve psArray *cells; // Cells for convolved images --- a handle for reading again @@ -26,7 +36,8 @@ psVector *matchChi2; // chi^2 for stamps from matching psVector *weightings; // Combination weightings for images (1/noise^2) - psArray *covariances; // Covariance matrices + psArray *convCovars; // Convolved covariance matrices // Combine initial pmReadout *outRO; // Output readout + pmReadout *expRO; // Exposure readout psArray *inspect; // Array of arrays of pixels to inspect // Rejection Index: branches/tap_branches/ppStack/src/ppStackPSF.c =================================================================== --- branches/tap_branches/ppStack/src/ppStackPSF.c (revision 25900) +++ branches/tap_branches/ppStack/src/ppStackPSF.c (revision 27838) @@ -14,4 +14,6 @@ const psArray *psfs, const psVector *inputMask) { + bool mdok = false; + #ifndef TESTING // Get the recipe values @@ -24,4 +26,11 @@ int psfOrder = psMetadataLookupS32(NULL, recipe, "PSF.ORDER"); // Spatial order for PSF + psString maskValStr = psMetadataLookupStr(&mdok, recipe, "MASK.VAL"); // Name of bits to mask going in + if (!mdok || !maskValStr) { + psError(PPSTACK_ERR_CONFIG, false, "Unable to find MASK.VAL in recipe"); + return false; + } + psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask + for (int i = 0; i < psfs->n; i++) { if (inputMask->data.U8[i]) { @@ -32,8 +41,7 @@ // Solve for the target PSF - pmPSF *psf = pmPSFEnvelope(numCols, numRows, psfs, psfInstances, psfRadius, psfModel, - psfOrder, psfOrder); + pmPSF *psf = pmPSFEnvelope(numCols, numRows, psfs, psfInstances, psfRadius, psfModel, psfOrder, psfOrder, maskVal); if (!psf) { - psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF."); + psError(PPSTACK_ERR_PSF, false, "Unable to determine output PSF."); return NULL; } @@ -42,5 +50,5 @@ pmPSF *psf = pmPSFBuildSimple("PS_MODEL_PS1_V1", 4.0, 4.0, 0.0, 1.0); if (!psf) { - psError(PS_ERR_UNKNOWN, false, "Unable to build dummy PSF."); + psError(PPSTACK_ERR_PSF, false, "Unable to build dummy PSF."); return NULL; } Index: branches/tap_branches/ppStack/src/ppStackPhotometry.c =================================================================== --- branches/tap_branches/ppStack/src/ppStackPhotometry.c (revision 25900) +++ branches/tap_branches/ppStack/src/ppStackPhotometry.c (revision 27838) @@ -19,6 +19,5 @@ psAssert(recipe, "We've thrown an error on this before."); - bool mdok; // Status of MD lookup - if (!psMetadataLookupBool(&mdok, recipe, "PHOTOMETRY")) { + if (!options->photometry) { // Nothing to do return true; @@ -27,8 +26,8 @@ psTimerStart("PPSTACK_PHOT"); - ppStackFileActivation(config, PPSTACK_FILES_COMBINE, false); + pmFPAfileActivate(config->files, false, NULL); ppStackFileActivation(config, PPSTACK_FILES_PHOT, true); pmFPAview *photView = ppStackFilesIterateDown(config); // View to readout - ppStackFileActivation(config, PPSTACK_FILES_COMBINE, true); + ppStackFileActivation(config, PPSTACK_FILES_STACK, true); pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT"); // File for photometry @@ -60,12 +59,12 @@ "Bits to use for mark", markValue); - pmCell *sourcesCell = pmFPAfileThisCell(config->files, photView, "PPSTACK.OUTPUT"); - psArray *inSources = psMetadataLookupPtr(NULL, sourcesCell->analysis, "PSPHOT.SOURCES"); + psArray *inSources = options->sources; if (!inSources) { - psError(PS_ERR_UNKNOWN, false, "Unable to find input sources"); + psError(PPSTACK_ERR_PROG, false, "Unable to find input sources"); psFree(photView); return false; } + pmModelClassSetLimits(PM_MODEL_LIMITS_LAX); if (!psphotReadoutKnownSources(config, photView, inSources)) { // This is likely a data quality issue @@ -73,7 +72,6 @@ psErrorStackPrint(stderr, "Unable to perform photometry on image"); psWarning("Unable to perform photometry on image --- suspect bad data quality."); - if (options->stats && psMetadataLookupS32(NULL, options->stats, "QUALITY") == 0) { - psMetadataAddS32(options->stats, PS_LIST_TAIL, "QUALITY", PS_META_REPLACE, - "Unable to perform photometry on image", psErrorCodeLast()); + if (options->quality == 0) { + options->quality = psErrorCodeLast(); } psErrorClear(); @@ -84,5 +82,5 @@ !pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV") || !pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND")) { - psError(PS_ERR_UNKNOWN, false, "Unable to drop PSPHOT internal files."); + psError(PPSTACK_ERR_PROG, false, "Unable to drop PSPHOT internal files."); return false; } @@ -92,10 +90,13 @@ if (options->stats) { pmReadout *photRO = pmFPAviewThisReadout(photView, photFile->fpa); // Readout with the sources - psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources - psMetadataAddS32(options->stats, PS_LIST_TAIL, "NUM_SOURCES", 0, - "Number of sources detected", sources ? sources->n : 0); - psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_PHOT", 0, - "Time to do photometry", psTimerMark("PPSTACK_PHOT")); + pmDetections *detections = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.DETECTIONS"); // detections + if (detections && detections->allSources) { + psMetadataAddS32(options->stats, PS_LIST_TAIL, "NUM_SOURCES", 0, "Number of sources detected", detections->allSources->n); + } else { + psMetadataAddS32(options->stats, PS_LIST_TAIL, "NUM_SOURCES", 0, "Number of sources detected", 0); + } + psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_PHOT", 0, "Time to do photometry", psTimerMark("PPSTACK_PHOT")); } + psFree(photView); Index: branches/tap_branches/ppStack/src/ppStackPrepare.c =================================================================== --- branches/tap_branches/ppStack/src/ppStackPrepare.c (revision 25900) +++ branches/tap_branches/ppStack/src/ppStackPrepare.c (revision 27838) @@ -26,15 +26,15 @@ float zpRadius = psMetadataLookupS32(&mdok, recipe, "PHOT.RADIUS"); // Radius for PHOT measurement if (!mdok) { - psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.RADIUS in recipe"); + psError(PPSTACK_ERR_CONFIG, true, "Unable to find PHOT.RADIUS in recipe"); return false; } float zpSigma = psMetadataLookupF32(&mdok, recipe, "PHOT.SIGMA"); // Gaussian sigma for photometry if (!mdok) { - psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.SIGMA in recipe"); + psError(PPSTACK_ERR_CONFIG, true, "Unable to find PHOT.SIGMA in recipe"); return false; } float zpFrac = psMetadataLookupF32(&mdok, recipe, "PHOT.FRAC"); // Fraction of good pixels for photometry if (!mdok) { - psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.FRAC in recipe"); + psError(PPSTACK_ERR_CONFIG, true, "Unable to find PHOT.FRAC in recipe"); return false; } @@ -42,5 +42,5 @@ psString maskValStr = psMetadataLookupStr(&mdok, recipe, "MASK.VAL"); // Name of bits to mask going in if (!mdok || !maskValStr) { - psError(PS_ERR_UNKNOWN, false, "Unable to find MASK.VAL in recipe"); + psError(PPSTACK_ERR_CONFIG, false, "Unable to find MASK.VAL in recipe"); return false; } @@ -53,5 +53,5 @@ psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics if (!psImageBackground(stats, NULL, image, mask, maskVal, rng)) { - psError(PS_ERR_UNKNOWN, false, "Unable to measure background for image"); + psError(PPSTACK_ERR_DATA, false, "Unable to measure background for image"); psFree(stats); psFree(rng); @@ -126,4 +126,6 @@ options->inputMask = psVectorAlloc(num, PS_TYPE_VECTOR_MASK); // Mask for inputs psVectorInit(options->inputMask, 0); + options->exposures = psVectorAlloc(options->num, PS_TYPE_F32); + psVectorInit(options->exposures, NAN); pmFPAfileActivate(config->files, false, NULL); @@ -134,28 +136,28 @@ } - psMetadataIterator *fileIter = psMetadataIteratorAlloc(config->files, PS_LIST_HEAD, "^PPSTACK.INPUT$"); - psMetadataItem *fileItem; // Item from iteration psArray *psfs = psArrayAlloc(num); // PSFs for PSF envelope int numCols = 0, numRows = 0; // Size of image - int index = 0; // Index for file - while ((fileItem = psMetadataGetAndIncrement(fileIter))) { - assert(fileItem->type == PS_DATA_UNKNOWN); - pmFPAfile *inputFile = fileItem->data.V; // An input file + options->sumExposure = 0.0; + for (int i = 0; i < num; i++) { + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File of interest + pmCell *cell = pmFPAviewThisCell(view, file->fpa); // Cell of interest + + options->exposures->data.F32[i] = psMetadataLookupF32(NULL, cell->concepts, "CELL.EXPOSURE"); + options->sumExposure += options->exposures->data.F32[i]; // Get list of PSFs, to determine target PSF if (options->convolve) { - pmChip *chip = pmFPAviewThisChip(view, inputFile->fpa); // The chip: holds the PSF + pmChip *chip = pmFPAviewThisChip(view, file->fpa); // The chip: holds the PSF pmPSF *psf = psMetadataLookupPtr(NULL, chip->analysis, "PSPHOT.PSF"); // PSF if (!psf) { - psError(PS_ERR_UNKNOWN, false, "Unable to find PSF."); - psFree(view); - psFree(fileIter); + psError(PPSTACK_ERR_PROG, false, "Unable to find PSF."); + psFree(view); psFree(psfs); return false; } - psfs->data[index] = psMemIncrRefCounter(psf); + psfs->data[i] = psMemIncrRefCounter(psf); psMetadataRemoveKey(chip->analysis, "PSPHOT.PSF"); - pmCell *cell = pmFPAviewThisCell(view, inputFile->fpa); // Cell of interest + pmCell *cell = pmFPAviewThisCell(view, file->fpa); // Cell of interest pmHDU *hdu = pmHDUFromCell(cell); assert(hdu && hdu->header); @@ -163,7 +165,6 @@ int naxis2 = psMetadataLookupS32(NULL, hdu->header, "NAXIS2"); // Number of rows if (naxis1 <= 0 || naxis2 <= 0) { - psError(PS_ERR_UNKNOWN, false, "Unable to determine size of image from PSF."); - psFree(view); - psFree(fileIter); + psError(PPSTACK_ERR_PROG, false, "Unable to determine size of image from PSF."); + psFree(view); psFree(psfs); return false; @@ -176,18 +177,29 @@ - pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Readout with sources - psArray *sources = psMetadataLookupPtr(NULL, ro->analysis, "PSPHOT.SOURCES"); // Sources - if (!sources) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find sources in readout."); - return NULL; - } - options->sourceLists->data[index] = psMemIncrRefCounter(sources); + bool redoPhot = psMetadataLookupBool(NULL, recipe, "PHOT"); + + pmDetections *detections = NULL; + if (options->convolve || options->matchZPs || options->photometry || redoPhot) { + pmReadout *ro = pmFPAviewThisReadout(view, file->fpa); // Readout with sources + detections = psMetadataLookupPtr(NULL, ro->analysis, "PSPHOT.DETECTIONS"); // Sources + if (!detections || !detections->allSources) { + psWarning("No detections found for image %d --- rejecting.", i); + options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_CAL; + continue; + } + psAssert (detections->allSources, "missing sources?"); + + options->sourceLists->data[i] = psMemIncrRefCounter(detections->allSources); + } // Re-do photometry if we don't trust the source lists - if (psMetadataLookupBool(NULL, recipe, "PHOT")) { - psTrace("ppStack", 2, "Photometering input %d of %d....\n", index, num); + if (redoPhot) { + psTrace("ppStack", 2, "Photometering input %d of %d....\n", i, num); pmFPAfileActivate(config->files, false, NULL); - ppStackFileActivationSingle(config, PPSTACK_FILES_CONVOLVE, true, index); - pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", index); // File + ppStackFileActivationSingle(config, PPSTACK_FILES_CONVOLVE, true, i); + if (options->convolve) { + pmFPAfileActivate(config->files, true, "PPSTACK.CONV.KERNEL"); + } + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File pmFPAview *photView = ppStackFilesIterateDown(config); if (!photView) { @@ -198,6 +210,6 @@ pmReadout *ro = pmFPAviewThisReadout(view, file->fpa); // Readout of interest - if (!ppStackInputPhotometer(ro, sources, config)) { - psError(PS_ERR_UNKNOWN, false, "Unable to do photometry on input sources"); + if (!ppStackInputPhotometer(ro, detections->allSources, config)) { + psError(psErrorCodeLast(), false, "Unable to do photometry on input sources"); psFree(view); psFree(photView); @@ -213,8 +225,45 @@ ppStackFileActivation(config, PPSTACK_FILES_PREPARE, true); } - - index++; - } - psFree(fileIter); + } + + psString log = psStringCopy("Input seeing FWHMs:\n"); // Log message + bool havePSFs = false; // Do we have any PSFs? + options->inputSeeing = psVectorAlloc(num, PS_TYPE_F32); + psVectorInit(options->inputSeeing, NAN); + for (int i = 0; i < num; i++) { + pmPSF *psf = psfs->data[i]; // PSF for image + if (!psf) { + continue; + } + havePSFs = true; + + int xNum = PS_MAX(psf->trendNx, 1), yNum = PS_MAX(psf->trendNy, 1); // Number of realisations + double sumFWHM = 0.0; // FWHM for image + int numFWHM = 0; // Number of FWHM measurements + for (int y = 0; y < yNum; y++) { + float yPos = ((float)y + 0.5) / (float)yNum * numRows; + for (int x = 0; x < xNum; x++) { + float xPos = ((float)x + 0.5) / (float)xNum * numCols; + float fwhm = pmPSFtoFWHM(psf, xPos, yPos); // FWHM for image + if (isfinite(fwhm)) { + sumFWHM += fwhm; + numFWHM++; + } + } + } + if (numFWHM == 0) { + options->inputSeeing->data.F32[i] = NAN; + options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_PSF; + psLogMsg("ppStack", PS_LOG_INFO, "Unable to measure PSF FWHM for image %d --- rejected.", i); + } else { + options->inputSeeing->data.F32[i] = sumFWHM / (float)numFWHM; + } + + psStringAppend(&log, "Input %d: %f\n", i, options->inputSeeing->data.F32[i]); + } + if (havePSFs) { + psLogMsg("ppStack", PS_LOG_INFO, "%s", log); + } + psFree(log); // Generate target PSF @@ -223,5 +272,5 @@ psFree(psfs); if (!options->psf) { - psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF."); + psError(psErrorCodeLast(), false, "Unable to determine output PSF."); psFree(view); return false; @@ -229,6 +278,8 @@ psMetadataAddPtr(config->arguments, PS_LIST_TAIL, "PSF.TARGET", PS_DATA_UNKNOWN, "Target PSF for stack", options->psf); - - pmChip *outChip = pmFPAfileThisChip(config->files, view, "PPSTACK.OUTPUT"); // Output chip + options->targetSeeing = pmPSFtoFWHM(options->psf, 0.5 * numCols, 0.5 * numRows); // FWHM for target + psLogMsg("ppStack", PS_LOG_INFO, "Target seeing FWHM: %f\n", options->targetSeeing); + + pmChip *outChip = pmFPAfileThisChip(config->files, view, "PPSTACK.TARGET.PSF"); // Output chip psMetadataAddPtr(outChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN, "Target PSF", options->psf); @@ -238,5 +289,5 @@ // Zero point calibration if (!ppStackSourcesTransparency(options, view, config)) { - psError(PS_ERR_UNKNOWN, false, "Unable to calculate transparency differences"); + psError(PPSTACK_ERR_DATA, false, "Unable to calculate transparency differences"); psFree(view); return false; Index: branches/tap_branches/ppStack/src/ppStackReadout.c =================================================================== --- branches/tap_branches/ppStack/src/ppStackReadout.c (revision 25900) +++ branches/tap_branches/ppStack/src/ppStackReadout.c (revision 27838) @@ -23,13 +23,12 @@ psVector *mask = options->inputMask; // Mask for inputs psVector *weightings = options->weightings; // Weightings (1/noise^2) for each image + psVector *exposures = options->exposures; // Exposure times for each image psVector *addVariance = options->matchChi2; // Additional variance when rejecting - psArray *inspect = ppStackReadoutInitial(config, outRO, thread->readouts, mask, - weightings, addVariance); - - job->results = inspect; + job->results = ppStackReadoutInitial(config, outRO, thread->readouts, mask, + weightings, exposures, addVariance); thread->busy = false; - return inspect ? true : false; + return job->results ? true : false; } @@ -40,18 +39,23 @@ psArray *args = job->args; // Arguments ppStackThread *thread = args->data[0]; // Thread - ppStackOptions *options = args->data[1]; // Options - pmConfig *config = args->data[2]; // Configuration - - pmReadout *outRO = options->outRO; // Output readout + psArray *reject = args->data[1]; // Rejected pixels for each image + ppStackOptions *options = args->data[2]; // Options + pmConfig *config = args->data[3]; // Configuration + bool safety = PS_SCALAR_VALUE(args->data[4], U8); // Safety switch on? + bool normalise = PS_SCALAR_VALUE(args->data[5], U8); // Normalise images? + psVector *mask = options->inputMask; // Mask for inputs - psArray *rejected = options->rejected; // Rejected pixels psVector *weightings = options->weightings; // Weightings (1/noise^2) for each image + psVector *exposures = options->exposures; // Exposure times for each image psVector *addVariance = options->matchChi2; // Additional variance when rejecting - - bool status = ppStackReadoutFinal(config, outRO, thread->readouts, mask, rejected, - weightings, addVariance); // Status of operation + psVector *norm = normalise ? options->norm : NULL; // Normalisations to apply to images + + bool status = ppStackReadoutFinal(config, options->outRO, options->expRO, thread->readouts, mask, reject, + weightings, exposures, addVariance, safety, norm); // Status of operation thread->busy = false; + psAssert(status, "Stacking failed."); + return status; } @@ -63,24 +67,35 @@ psArray *args = job->args; // Input arguments - psArray *inspect = args->data[0]; // Array of pixel arrays - int index = PS_SCALAR_VALUE(args->data[1], S32); // Index of interest - - psArray *inputs = inspect->data[index]; // Array of interest - psPixels *output = NULL; // Output pixel list - for (int i = 0; i < inputs->n; i++) { - psPixels *input = inputs->data[i]; // Input pixel list - if (!input || input->n == 0) { - continue; - } - output = psPixelsConcatenate(output, input); - } - - if (!output) { - // If there are no pixels to inspect, then just fake it - output = psPixelsAllocEmpty(0); - } - - psFree(inputs); - inspect->data[index] = output; + psArray *inspects = args->data[0]; // Array of pixel arrays + psArray *rejects = args->data[1]; // Array of pixel arrays + int index = PS_SCALAR_VALUE(args->data[2], S32); // Index of interest + + psArray *inInspects = inspects->data[index]; // Array of interest + psArray *inRejects = rejects->data[index]; // Array of interest + psAssert(inInspects->n == inRejects->n, "Size should be the same"); + psPixels *outInspect = NULL, *outReject = NULL; // Output pixel lists + for (int i = 0; i < inInspects->n; i++) { + psPixels *inInspect = inInspects->data[i]; // Input pixel list + if (inInspect && inInspect->n > 0) { + outInspect = psPixelsConcatenate(outInspect, inInspect); + } + psPixels *inReject = inRejects->data[i]; // Input pixel list + if (inReject && inReject->n > 0) { + outReject = psPixelsConcatenate(outReject, inReject); + } + } + + // If there are no pixels to inspect, then just fake it + if (!outInspect) { + outInspect = psPixelsAllocEmpty(0); + } + if (!outReject) { + outReject = psPixelsAllocEmpty(0); + } + + psFree(inspects->data[index]); + inspects->data[index] = outInspect; + psFree(rejects->data[index]); + rejects->data[index] = outReject; return true; @@ -89,5 +104,6 @@ psArray *ppStackReadoutInitial(const pmConfig *config, pmReadout *outRO, const psArray *readouts, - const psVector *mask, const psVector *weightings, const psVector *addVariance) + const psVector *mask, const psVector *weightings, const psVector *exposures, + const psVector *addVariance) { assert(config); @@ -104,5 +120,5 @@ bool mdok; // Status of MD lookup - int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations + float iter = psMetadataLookupF32(NULL, recipe, "COMBINE.ITER"); // Rejection iterations float combineRej = psMetadataLookupF32(NULL, recipe, "COMBINE.REJ"); // Combination threshold float combineSys = psMetadataLookupF32(NULL, recipe, "COMBINE.SYS"); // Combination systematic error @@ -116,4 +132,6 @@ psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch + psString maskSuspectStr = psMetadataLookupStr(NULL, recipe, "MASK.SUSPECT"); // Name of suspect mask bits + psImageMaskType maskSuspect = pmConfigMaskGet(maskSuspectStr, config); // Suspect bits psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels @@ -129,15 +147,4 @@ } -#if 0 - // This doesn't seem to work, so getting the weightings directly from a vector - float weighting = psMetadataLookupF32(&mdok, ro->analysis, "PPSTACK.WEIGHTING"); // Relative weight - if (!mdok || !isfinite(weighting)) { - psWarning("No weighting supplied for image %d --- set to unity.", i); - weighting = 1.0; - } else { - psLogMsg("ppStack", PS_LOG_INFO, "Weighting for image %d is %f", i, weighting); - } -#endif - // Ensure there is a mask, or pmStackCombine will complain if (!ro->mask) { @@ -146,9 +153,10 @@ } - stack->data[i] = pmStackDataAlloc(ro, weightings->data.F32[i], addVariance->data.F32[i]); - } - - if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, kernelSize, iter, - combineRej, combineSys, combineDiscard, true, useVariance, safe, false)) { + stack->data[i] = pmStackDataAlloc(ro, weightings->data.F32[i], exposures->data.F32[i], + addVariance->data.F32[i]); + } + + if (!pmStackCombine(outRO, NULL, stack, maskVal | maskBad, maskSuspect, maskBad, kernelSize, iter, + combineRej, combineSys, combineDiscard, useVariance, safe, false)) { psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection."); psFree(stack); @@ -156,6 +164,7 @@ } - // Save list of pixels to inspect + // Save lists of pixels psArray *inspect = psArrayAlloc(num); // List of pixels to inspect + psArray *reject = psArrayAlloc(num); // List of pixels rejected for (int i = 0; i < num; i++) { pmStackData *data = stack->data[i]; // Data for this image @@ -168,4 +177,5 @@ } inspect->data[i] = psMemIncrRefCounter(data->inspect); + reject->data[i] = psMemIncrRefCounter(data->reject); } psFree(stack); @@ -173,15 +183,21 @@ sectionNum++; - return inspect; -} - - - -bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, const psArray *readouts, + psArray *results = psArrayAlloc(2); // Array of results + results->data[0] = inspect; + results->data[1] = reject; + + return results; +} + + + +bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, pmReadout *expRO, const psArray *readouts, const psVector *mask, const psArray *rejected, const psVector *weightings, - const psVector *addVariance) + const psVector *exposures, const psVector *addVariance, bool safety, + const psVector *norm) { assert(config); assert(outRO); + assert(expRO); assert(readouts); assert(!rejected || readouts->n == rejected->n); @@ -196,8 +212,4 @@ bool mdok; // Status of MD lookup - int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations - float combineRej = psMetadataLookupF32(NULL, recipe, "COMBINE.REJ"); // Combination threshold - float combineSys = psMetadataLookupF32(NULL, recipe, "COMBINE.SYS"); // Combination systematic error - float combineDiscard = psMetadataLookupF32(NULL, recipe, "COMBINE.DISCARD"); // Olympic discard fraction bool useVariance = psMetadataLookupBool(&mdok, recipe, "VARIANCE"); // Use variance for rejection? bool safe = psMetadataLookupBool(&mdok, recipe, "SAFE"); // Be safe when combining small numbers of pixels @@ -205,4 +217,6 @@ psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch + psString maskSuspectStr = psMetadataLookupStr(NULL, recipe, "MASK.SUSPECT"); // Name of suspect mask bits + psImageMaskType maskSuspect = pmConfigMaskGet(maskSuspectStr, config); // Suspect bits psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels @@ -211,34 +225,17 @@ psArray *stack = psArrayAlloc(num); // Array for stacking - bool entire = true; // Combine entire image? - if (rejected) { - // We have rejection from a previous combination: combine without flagging pixels to inspect - entire = false; - safe = false; - iter = 0; - combineRej = NAN; - combineSys = NAN; - } + // We have rejection from a previous combination: combine without flagging pixels to inspect + safe &= safety; + int iter = 0; + float combineRej = NAN; + float combineSys = NAN; + float combineDiscard = NAN; for (int i = 0; i < num; i++) { pmReadout *ro = readouts->data[i]; - if (mask->data.U8[i] & (PPSTACK_MASK_REJECT | PPSTACK_MASK_BAD)) { - // Image completely rejected since previous combination - entire = true; - continue; - } else if (mask->data.U8[i]) { - // Image completely rejected before original combination - continue; - } - -#if 0 - // This doesn't seem to work, so getting the weightings directly from a vector - bool mdok; // Status of MD lookup - float weighting = psMetadataLookupF32(&mdok, ro->analysis, "PPSTACK.WEIGHTING"); // Relative weight - if (!mdok || !isfinite(weighting)) { - psWarning("No WEIGHTING supplied for image %d --- set to unity.", i); - weighting = 1.0; - } -#endif + if (mask->data.U8[i]) { + // Image completely rejected + continue; + } // Ensure there is a mask, or pmStackCombine will complain @@ -248,13 +245,18 @@ } - pmStackData *data = pmStackDataAlloc(ro, weightings->data.F32[i], + pmStackData *data = pmStackDataAlloc(ro, weightings->data.F32[i], exposures->data.F32[i], addVariance ? addVariance->data.F32[i] : NAN); data->reject = rejected ? psMemIncrRefCounter(rejected->data[i]) : NULL; stack->data[i] = data; - } - - if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, 0, - iter, combineRej, combineSys, combineDiscard, - entire, useVariance, safe, !rejected)) { + + if (norm) { + float normalise = powf(10.0, -0.4 * norm->data.F32[i]); // Normalisation + psBinaryOp(ro->image, ro->image, "*", psScalarAlloc(normalise, PS_TYPE_F32)); + psBinaryOp(ro->variance, ro->variance, "*", psScalarAlloc(PS_SQR(normalise), PS_TYPE_F32)); + } + } + + if (!pmStackCombine(outRO, expRO, stack, maskVal | maskBad, maskSuspect, maskBad, 0, iter, combineRej, + combineSys, combineDiscard, useVariance, safe, rejected)) { psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts."); psFree(stack); @@ -264,9 +266,14 @@ pmCell *outCell = outRO->parent; // Output cell pmChip *outChip = outCell->parent; // Output chip - outRO->data_exists = true; outCell->data_exists = true; outChip->data_exists = true; + pmCell *expCell = expRO->parent; // Exposure cell + pmChip *expChip = expCell->parent; // Exposure chip + expRO->data_exists = true; + expCell->data_exists = true; + expChip->data_exists = true; + psFree(stack); Index: branches/tap_branches/ppStack/src/ppStackReject.c =================================================================== --- branches/tap_branches/ppStack/src/ppStackReject.c (revision 25900) +++ branches/tap_branches/ppStack/src/ppStackReject.c (revision 27838) @@ -10,5 +10,5 @@ #include "ppStackLoop.h" -//#define TESTING +// #define TESTING bool ppStackReject(ppStackOptions *options, pmConfig *config) @@ -23,5 +23,4 @@ int num = options->num; // Number of inputs - options->rejected = psArrayAlloc(num); psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe @@ -32,5 +31,4 @@ float threshold = psMetadataLookupF32(NULL, recipe, "THRESHOLD.MASK"); // Threshold for mask deconvolution - float poorFrac = psMetadataLookupF32(NULL, recipe, "POOR.FRACTION"); // Fraction for "poor" float imageRej = psMetadataLookupF32(NULL, recipe, "IMAGE.REJ"); // Maximum fraction of image to reject // before rejecting entire image @@ -53,4 +51,5 @@ psThreadJob *job = psThreadJobAlloc("PPSTACK_INSPECT"); // Job to start psArrayAdd(job->args, 1, options->inspect); + psArrayAdd(job->args, 1, options->rejected); PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32); if (!psThreadJobAddPending(job)) { @@ -61,5 +60,5 @@ } if (!psThreadPoolWait(true)) { - psError(PS_ERR_UNKNOWN, false, "Unable to concatenate inspection lists."); + psError(psErrorCodeLast(), false, "Unable to concatenate inspection lists."); return false; } @@ -92,15 +91,11 @@ #endif - psPixels *reject = pmStackReject(options->inspect->data[i], options->numCols, options->numRows, - threshold, poorFrac, stride, options->regions->data[i], - options->kernels->data[i]); // Rejected pixels - #ifdef TESTING { - psImage *mask = psPixelsToMask(NULL, reject, + psImage *mask = psPixelsToMask(NULL, options->rejected->data[i], psRegionSet(0, options->numCols - 1, 0, options->numRows - 1), 0xff); // Mask image psString name = NULL; // Name of image - psStringAppend(&name, "reject_%03d.fits", i); + psStringAppend(&name, "pre_reject_%03d.fits", i); pmStackVisualPlotTestImage(mask, name); psFits *fits = psFitsOpen(name, "w"); @@ -111,4 +106,8 @@ } #endif + + psPixels *reject = pmStackReject(options->inspect->data[i], options->numCols, options->numRows, + threshold, stride, options->regions->data[i], + options->kernels->data[i]); // Rejected pixels psFree(options->inspect->data[i]); @@ -127,13 +126,29 @@ "exceeds limit (%.3f)", i, frac, imageRej); psFree(reject); - // reject == NULL means reject image completely reject = NULL; options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PPSTACK_MASK_BAD; numRejected++; + } else { + // Add to list of pixels already rejected + reject = psPixelsConcatenate(reject, options->rejected->data[i]); + options->rejected->data[i] = psPixelsDuplicates(options->rejected->data[i], reject); } } - // Images without a list of rejected pixels (the list may be empty) are rejected completely - options->rejected->data[i] = reject; +#ifdef TESTING + { + psImage *mask = psPixelsToMask(NULL, options->rejected->data[i], + psRegionSet(0, options->numCols - 1, 0, options->numRows - 1), + 0xff); // Mask image + psString name = NULL; // Name of image + psStringAppend(&name, "reject_%03d.fits", i); + pmStackVisualPlotTestImage(mask, name); + psFits *fits = psFitsOpen(name, "w"); + psFree(name); + psFitsWriteImage(fits, NULL, mask, 0, NULL); + psFree(mask); + psFitsClose(fits); + } +#endif if (options->stats) { @@ -143,4 +158,6 @@ "Number of pixels rejected", reject ? reject->n : 0); } + + psFree(reject); psLogMsg("ppStack", PS_LOG_INFO, "Time to perform rejection on image %d: %f sec", i, psTimerClear("PPSTACK_REJECT")); @@ -148,9 +165,7 @@ psFree(options->inspect); options->inspect = NULL; - psFree(options->kernels); options->kernels = NULL; - psFree(options->regions); options->regions = NULL; - if (numRejected >= num - 1) { - psError(PS_ERR_UNKNOWN, true, "All inputs completely rejected; unable to proceed."); + if (numRejected >= num) { + psError(PPSTACK_ERR_REJECTED, true, "All inputs completely rejected; unable to proceed."); return false; } Index: branches/tap_branches/ppStack/src/ppStackSetup.c =================================================================== --- branches/tap_branches/ppStack/src/ppStackSetup.c (revision 25900) +++ branches/tap_branches/ppStack/src/ppStackSetup.c (revision 27838) @@ -22,4 +22,9 @@ psAssert(recipe, "We've thrown an error on this before."); + // XXX : switch to this name? options->matchZPs = psMetadataLookupBool(NULL, recipe, "MATCH.ZERO.POINTS"); // Adjust zero points based on tranparency analysis? + options->matchZPs = psMetadataLookupBool(NULL, recipe, "ZP"); // Adjust zero points? + + options->photometry = psMetadataLookupBool(NULL, recipe, "PHOTOMETRY"); // Perform photometry? + options->convolve = psMetadataLookupBool(NULL, recipe, "CONVOLVE"); // Convolve images? if (!psMetadataLookupBool(NULL, config->arguments, "HAVE.PSF")) { @@ -37,5 +42,5 @@ options->statsFile = fopen(resolved, "w"); if (!options->statsFile) { - psError(PS_ERR_IO, true, "Unable to open statistics file %s for writing.\n", resolved); + psError(PPSTACK_ERR_IO, true, "Unable to open statistics file %s for writing.\n", resolved); psFree(resolved); return false; @@ -43,11 +48,13 @@ psFree(resolved); options->stats = psMetadataAlloc(); - psMetadataAddS32(options->stats, PS_LIST_TAIL, "QUALITY", 0, "No problems", 0); } // Generate temporary names for convolved images - const char *tempDir = psMetadataLookupStr(NULL, recipe, "TEMP.DIR"); // Directory for temporary images + const char *tempDir = psMetadataLookupStr(NULL, config->arguments, "-temp-dir"); // Directory for temps if (!tempDir) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find TEMP.DIR in recipe"); + tempDir = psMetadataLookupStr(NULL, config->site, "TEMP.DIR"); + } + if (!tempDir) { + psError(PPSTACK_ERR_CONFIG, false, "Unable to find TEMP.DIR in site configuration"); return false; } @@ -57,5 +64,5 @@ const char *tempName = basename(outputName); if (!tempName) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to construct basename for temporary files."); + psError(PPSTACK_ERR_ARGUMENTS, false, "Unable to construct basename for temporary files."); psFree(outputName); return false; @@ -66,5 +73,5 @@ const char *tempVariance = psMetadataLookupStr(NULL, recipe, "TEMP.VARIANCE"); // Suffix for var maps if (!tempImage || !tempMask || !tempVariance) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, + psError(PPSTACK_ERR_CONFIG, false, "Unable to find TEMP.IMAGE, TEMP.MASK and TEMP.VARIANCE in recipe"); psFree(outputName); @@ -72,7 +79,7 @@ } - options->imageNames = psArrayAlloc(num); - options->maskNames = psArrayAlloc(num); - options->varianceNames = psArrayAlloc(num); + options->convImages = psArrayAlloc(num); + options->convMasks = psArrayAlloc(num); + options->convVariances = psArrayAlloc(num); for (int i = 0; i < num; i++) { psString imageName = NULL, maskName = NULL, varianceName = NULL; // Names for convolved images @@ -81,12 +88,33 @@ psStringAppend(&varianceName, "%s/%s.%d.%s", tempDir, tempName, i, tempVariance); psTrace("ppStack", 5, "Temporary files: %s %s %s\n", imageName, maskName, varianceName); - options->imageNames->data[i] = imageName; - options->maskNames->data[i] = maskName; - options->varianceNames->data[i] = varianceName; + options->convImages->data[i] = imageName; + options->convMasks->data[i] = maskName; + options->convVariances->data[i] = varianceName; } psFree(outputName); + // Original images + options->origImages = psArrayAlloc(num); + options->origMasks = psArrayAlloc(num); + options->origVariances = psArrayAlloc(num); + pmFPAview *view = pmFPAviewAlloc(0); + for (int i = 0; i < num; i++) { + { + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); + options->origImages->data[i] = pmFPAfileName(file, view, config); + } + { + // We want the convolved mask, since that defines the area that has been tested for outliers + options->origMasks->data[i] = psMemIncrRefCounter(options->convMasks->data[i]); + } + { + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT.VARIANCE", i); + options->origVariances->data[i] = pmFPAfileName(file, view, config); + } + } + psFree(view); + if (!pmConfigMaskSetBits(NULL, NULL, config)) { - psError(PS_ERR_UNKNOWN, false, "Unable to determine mask value."); + psError(PPSTACK_ERR_CONFIG, false, "Unable to determine mask value."); return false; } Index: branches/tap_branches/ppStack/src/ppStackSources.c =================================================================== --- branches/tap_branches/ppStack/src/ppStackSources.c (revision 25900) +++ branches/tap_branches/ppStack/src/ppStackSources.c (revision 27838) @@ -8,4 +8,6 @@ //#define TESTING // Enable debugging output + +//#define ASTROMETRY // Correct astrometry? #ifdef TESTING @@ -61,4 +63,10 @@ PS_ASSERT_PTR_NON_NULL(config, false); + if (!options->matchZPs && !options->photometry) { + options->norm = psVectorAlloc(options->num, PS_TYPE_F32); + psVectorInit(options->norm, 0.0); + return true; + } + psArray *sourceLists = options->sourceLists; // Source lists for each input psVector *inputMask = options->inputMask; // Mask for inputs @@ -79,4 +87,14 @@ } source->psfMag += 1.0; +#ifdef ASTROMETRY + if (source->modelPSF) { + source->modelPSF->params->data.F32[PM_PAR_XPOS] += 1.0; + source->modelPSF->params->data.F32[PM_PAR_YPOS] += 1.0; + } + if (source->peak) { + source->peak->xf += 1.0; + source->peak->yf += 1.0; + } +#endif } } @@ -105,18 +123,34 @@ psMetadata *airmassZP = psMetadataLookupMetadata(NULL, recipe, "ZP.AIRMASS"); // Airmass terms if (!airmassZP) { - psError(PS_ERR_UNKNOWN, false, "Unable to find ZP.AIRMASS in recipe."); + psError(PPSTACK_ERR_CONFIG, false, "Unable to find ZP.AIRMASS in recipe."); return false; } - - int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs + psMetadata *zpTargetMenu = psMetadataLookupMetadata(NULL, recipe, "ZP.TARGET"); // Target zero point terms + if (!zpTargetMenu) { + psError(PPSTACK_ERR_CONFIG, false, "Unable to find ZP.TARGET in recipe."); + return false; + } + + int num = options->num; // Number of inputs psAssert(num == sourceLists->n, "Wrong number of source lists: %ld\n", sourceLists->n); - psVector *zp = psVectorAlloc(num, PS_TYPE_F32); // Zero points for each image + psVector *zp = psVectorAlloc(num, PS_TYPE_F32); // Relative zero points for each image + psVector *zpExp = psVectorAlloc(num, PS_TYPE_F32); // Measured zero points for each image (maybe) + int zpExpNum = 0; // Number of measured zero points const char *filter = NULL; // Filter name float airmassTerm = NAN; // Airmass term - float sumExpTime = 0.0; // Sum of the exposure time + float zpTarget = NAN; // Target zero point + int numGoodImages = 0; // Number of good images for (int i = 0; i < num; i++) { + psArray *sources = sourceLists->data[i]; // Source list + if (!sources || sources->n == 0) { + psLogMsg("ppStack", PS_LOG_WARN, "Image %d has no sources for transparency measurement.", i); + options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PPSTACK_MASK_CAL; + zp->data.F32[i] = NAN; + continue; + } + numGoodImages++; + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File of interest - pmCell *cell = pmFPAviewThisCell(view, file->fpa); // Cell of interest #if defined(TESTING) && 0 @@ -135,11 +169,14 @@ #endif - - float exptime = psMetadataLookupF32(NULL, cell->concepts, "CELL.EXPOSURE"); // Exposure time + float exptime = options->exposures->data.F32[i]; // Exposure time float airmass = psMetadataLookupF32(NULL, file->fpa->concepts, "FPA.AIRMASS"); // Airmass const char *expFilter = psMetadataLookupStr(NULL, file->fpa->concepts, "FPA.FILTER"); // Filter name + zpExp->data.F32[i] = psMetadataLookupF32(NULL, file->fpa->concepts, "FPA.ZP"); // Exposure zero point + psLogMsg("ppStack", PS_LOG_INFO, + "Image %d: %.2f sec exposure in %s at airmass %.2f with zero point %.2f", + i, exptime, expFilter, airmass, zpExp->data.F32[i]); if (!isfinite(exptime) || exptime == 0 || !isfinite(airmass) || airmass == 0 || !expFilter || strlen(expFilter) == 0) { - psError(PS_ERR_UNEXPECTED_NULL, false, + psError(PPSTACK_ERR_CONFIG, false, "Unable to find exposure time (%f), airmass (%f) or filter (%s)", exptime, airmass, expFilter); @@ -147,29 +184,70 @@ return false; } + if (isfinite(zpExp->data.F32[i])) { + zpExp->data.F32[i] += 2.5 * log10(exptime); + zpExpNum++; + } if (!filter) { filter = expFilter; - airmassTerm = psMetadataLookupF32(NULL, airmassZP, filter); - if (!isfinite(airmassTerm)) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, + bool mdok; + airmassTerm = psMetadataLookupF32(&mdok, airmassZP, filter); + if (!mdok || !isfinite(airmassTerm)) { + psError(PPSTACK_ERR_CONFIG, false, "Unable to find airmass term (ZP.AIRMASS) for filter %s", filter); psFree(zp); return false; } + zpTarget = psMetadataLookupF32(&mdok, zpTargetMenu, filter); + if (!mdok || !isfinite(zpTarget)) { + psError(PPSTACK_ERR_CONFIG, false, + "Unable to find target zero point (ZP.TARGET) for filter %s", filter); + psFree(zp); + return false; + } } else if (strcmp(filter, expFilter) != 0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Filters don't match: %s vs %s", filter, expFilter); + psError(PPSTACK_ERR_CONFIG, false, "Filters don't match: %s vs %s", filter, expFilter); psFree(zp); return false; } - zp->data.F32[i] = airmassTerm * airmass - 2.5 * log10(exptime); - sumExpTime += exptime; - } - - options->sumExposure = sumExpTime; + zp->data.F32[i] = airmassTerm * airmass + 2.5 * log10(exptime); + } + + if (numGoodImages == 0) { + psLogMsg("ppStack", PS_LOG_WARN, "No images with sources to measure transparency."); + options->quality = PPSTACK_ERR_REJECTED; + psFree(zp); + psFree(zpExp); + return true; + } + if (numGoodImages == 1) { + psArray *sources = NULL; // Sources + for (int i = 0; i < num && !sources; i++) { + if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) { + continue; + } + sources = sourceLists->data[i]; + } + options->quality = PPSTACK_ERR_REJECTED; + options->sources = psMemIncrRefCounter(sources); + options->norm = psVectorAlloc(num, PS_TYPE_F32); + psVectorInit(options->norm, 1.0); + options->zp = NAN; + psLogMsg("ppStack", PS_LOG_WARN, "Single image with sources --- cannot match transparency."); + psFree(zp); + psFree(zpExp); + return true; + } + + if (zpExpNum == numGoodImages) { + for (int i = 0; i < num; i++) { + zp->data.F32[i] = zpExp->data.F32[i]; + } + } psArray *matches = pmSourceMatchSources(sourceLists, radius, true); // List of matches if (!matches) { - psError(PS_ERR_UNKNOWN, false, "Unable to match sources"); + psError(PPSTACK_ERR_DATA, false, "Unable to match sources"); psFree(zp); return false; @@ -180,92 +258,171 @@ #endif - psVector *trans = pmSourceMatchRelphot(matches, zp, tol, iter1, starRej1, starSys1, - iter2, starRej2, starSys2, starLimit, - transIter, transRej, transThresh); // Transparencies for each image - if (!trans) { - psError(PS_ERR_UNKNOWN, false, "Unable to measure transparencies"); - return false; - } + if (options->matchZPs) { + psVector *trans = pmSourceMatchRelphot(matches, zp, tol, iter1, starRej1, starSys1, + iter2, starRej2, starSys2, starLimit, + transIter, transRej, transThresh); // Transparencies per image + if (!trans) { + psError(PPSTACK_ERR_DATA, false, "Unable to measure transparencies"); + return false; + } + for (int i = 0; i < trans->n; i++) { + if (!isfinite(trans->data.F32[i])) { + inputMask->data.U8[i] |= PPSTACK_MASK_CAL; + } + } + // M = m + c0 + c1 * airmass - 2.5log(t) + transparency + // Want sources to have m corresponding to airmass = 1 and t = sumExpTime and transparency = 0 + // m_0 + c1 * airmass_0 + 2.5log(t_0) - trans_0 = m_1 + c1 * airmass_1 + 2.5log(t_1) - trans_1 + // m_0 = m_1 + zp_1 + trans_1 - c1 * airmass_0 - 2.5log(t_0) + // We don't need to know the magnitude zero point for the filter, since it cancels out + if (options->matchZPs) { + options->norm = psVectorAlloc(num, PS_TYPE_F32); + for (int i = 0; i < num; i++) { + if (inputMask->data.U8[i] || !isfinite(trans->data.F32[i])) { + continue; + } + psArray *sources = sourceLists->data[i]; // Sources of interest + float magCorr = zp->data.F32[i] + trans->data.F32[i] - 2.5*log10(options->sumExposure); + if (zpExpNum == numGoodImages) { + // Using measured zero points, so attempt to set target zero point + magCorr -= zpTarget; + } + options->norm->data.F32[i] = -magCorr; + psLogMsg("ppStack", PS_LOG_INFO, + "Applying scale correction to image %d: %f mag (%f)\n", + i, magCorr, trans->data.F32[i]); + + for (int j = 0; j < sources->n; j++) { + pmSource *source = sources->data[j]; // Source of interest + if (!source) { + continue; + } + source->psfMag -= magCorr; + } + } + } + + if (zpExpNum == numGoodImages) { + // Producing image with target zero point + options->zp = zpTarget; + } else { + options->zp = NAN; + } + #ifdef TESTING - dumpMatches("source_mags.dat", num, matches, zp, trans); -#endif - - for (int i = 0; i < trans->n; i++) { - if (!isfinite(trans->data.F32[i])) { - inputMask->data.U8[i] |= PPSTACK_MASK_CAL; - } - } - - // Save best matches SOMEWHERE for future photometry - // XXX this is a really poor output location; clean up the pmFPAfiles used in ppStack - pmCell *sourcesCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT"); - psArray *sourcesBest = psArrayAllocEmpty(matches->n); - - // XXX something of a hack: require at least 2 detections or the nominated fraction of the max possible - int minMatches = PS_MAX(2, fracMatch * num);// Minimum number of matches required - for (int i = 0; i < matches->n; i++) { - pmSourceMatch *match = matches->data[i]; // Match of interest - if (match->num < minMatches) { - continue; - } - - // We need to grab a single instance of this source: just take the first available - int image = match->image->data.S32[0]; // Index of image - int index = match->index->data.S32[0]; // Index of source within image - psArray *sources = sourceLists->data[image]; // Sources for image - pmSource *source = sources->data[index]; // Source of interest - - psArrayAdd(sourcesBest, sourcesBest->n, source); - } - psMetadataAdd(sourcesCell->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY | PS_META_REPLACE, - "psphot sources", sourcesBest); - psLogMsg("ppStack", PS_LOG_INFO, "Selected %ld sources for photometry analysis", sourcesBest->n); - psFree(sourcesBest); + dumpMatches("source_mags.dat", num, matches, zp, trans); +#endif + psFree(trans); + +#ifdef TESTING + // Double check: all transparencies should be zero + { + psArray *matches = pmSourceMatchSources(sourceLists, radius, true); // List of matches + if (!matches) { + psError(PPSTACK_ERR_DATA, false, "Unable to match sources"); + psFree(zp); + return false; + } + psVector *trans = pmSourceMatchRelphot(matches, zp, tol, iter1, starRej1, starSys1, + iter2, starRej2, starSys2, starLimit, + transIter, transRej, transThresh); // Transparencies + for (int i = 0; i < num; i++) { + fprintf(stderr, "Transparency of image %d: %f\n", i, trans->data.F32[i]); + } + psFree(trans); + psFree(matches); + } +#endif + } + + psFree(zp); + psFree(zpExp); + +#ifdef ASTROMETRY + // Position offsets + { + psArray *offsets = pmSourceMatchRelastro(matches, num, tol, iter1, starRej1, + iter2, starRej2, starLimit); // Shifts for each image + if (!offsets) { + psError(PPSTACK_ERR_DATA, false, "Unable to measure offsets"); + return false; + } + for (int i = 0; i < num; i++) { + if (options->inputMask->data.U8[i]) { + continue; + } + psArray *sources = sourceLists->data[i]; // Sources of interest + psVector *offset = offsets->data[i]; // Offsets for image + float dx = offset->data.F32[0], dy = offset->data.F32[1]; // Offsets to apply + if (!isfinite(dx) || !isfinite(dy)) { + continue; + } + psLogMsg("ppStack", PS_LOG_INFO, "Applying astrometric correction to image %d: %f,%f\n", + i, dx, dy); + for (int j = 0; j < sources->n; j++) { + pmSource *source = sources->data[j]; // Source of interest + if (!source) { + continue; + } + if (source->modelPSF) { + source->modelPSF->params->data.F32[PM_PAR_XPOS] -= dx; + source->modelPSF->params->data.F32[PM_PAR_YPOS] -= dy; + } + if (source->peak) { + source->peak->xf -= dx; + source->peak->yf -= dy; + } + } + } + psFree(offsets); + } +#endif + +#if (defined TESTING && defined ASTROMETRY) + // Double check: all offsets should be zero + { + psArray *matches = pmSourceMatchSources(sourceLists, radius, true); // List of matches + if (!matches) { + psError(PPSTACK_ERR_DATA, false, "Unable to match sources"); + psFree(zp); + return false; + } + psArray *offsets = pmSourceMatchRelastro(matches, num, tol, iter1, starRej1, + iter2, starRej2, starLimit); // Shifts for each image + for (int i = 0; i < num; i++) { + psVector *offset = offsets->data[i]; // Offsets for image + fprintf(stderr, "Offset of image %d: %f,%f\n", i, offset->data.F32[0], offset->data.F32[1]); + } + psFree(offsets); + psFree(matches); + } +#endif + + + if (options->photometry) { + // Save best matches for future photometry + psArray *sourcesBest = options->sources = psArrayAllocEmpty(matches->n); + + // XXX something of a hack: require at least 2 detections or the nominated fraction of the max possible + int minMatches = PS_MAX(2, fracMatch * num);// Minimum number of matches required + for (int i = 0; i < matches->n; i++) { + pmSourceMatch *match = matches->data[i]; // Match of interest + if (match->num < minMatches) { + continue; + } + + // We need to grab a single instance of this source: just take the first available + int image = match->image->data.S32[0]; // Index of image + int index = match->index->data.S32[0]; // Index of source within image + psArray *sources = sourceLists->data[image]; // Sources for image + pmSource *source = sources->data[index]; // Source of interest + + psArrayAdd(sourcesBest, sourcesBest->n, source); + } + psLogMsg("ppStack", PS_LOG_INFO, "Selected %ld sources for photometry analysis", sourcesBest->n); + } psFree(matches); - - // M = m + c0 + c1 * airmass - 2.5log(t) + transparency - // Want sources to have m corresponding to airmass = 1 and t = sumExpTime and transparency = 0 - // m_0 + c1 * airmass_0 + 2.5log(t_0) - trans_0 = m_1 + c1 * airmass_1 + 2.5log(t_1) - trans_1 - // m_0 = m_1 + zp_1 + trans_1 - c1 * airmass_0 - 2.5log(t_0) - // We don't need to know the magnitude zero point for the filter, since it cancels out - options->norm = psVectorAlloc(num, PS_TYPE_F32); - for (int i = 0; i < num; i++) { - if (!isfinite(trans->data.F32[i])) { - continue; - } - psArray *sources = sourceLists->data[i]; // Sources of interest - float magCorr = airmassTerm - 2.5*log10(sumExpTime) - zp->data.F32[i] - trans->data.F32[i]; - options->norm->data.F32[i] = magCorr; - psLogMsg("ppStack", PS_LOG_INFO, "Applying magnitude correction to image %d: %f\n", i, magCorr); - - for (int j = 0; j < sources->n; j++) { - pmSource *source = sources->data[j]; // Source of interest - if (!source) { - continue; - } - source->psfMag += magCorr; - } - } - psFree(trans); - -#ifdef TESTING - // Double check: all transparencies should be zero - { - psArray *matches = pmSourceMatchSources(sourceLists, radius); // List of matches - if (!matches) { - psError(PS_ERR_UNKNOWN, false, "Unable to match sources"); - psFree(zp); - return false; - } - psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej, - transThresh, starRej, starSys); - for (int i = 0; i < num; i++) { - fprintf(stderr, "Transparency of image %d: %f\n", i, trans->data.F32[i]); - } - psFree(trans); - } -#endif return true; Index: branches/tap_branches/ppStack/src/ppStackThread.c =================================================================== --- branches/tap_branches/ppStack/src/ppStackThread.c (revision 25900) +++ branches/tap_branches/ppStack/src/ppStackThread.c (revision 27838) @@ -56,5 +56,5 @@ } -ppStackThreadData *ppStackThreadDataSetup(const ppStackOptions *options, const pmConfig *config) +ppStackThreadData *ppStackThreadDataSetup(const ppStackOptions *options, const pmConfig *config, bool conv) { psAssert(options, "Require options"); @@ -62,14 +62,22 @@ const psArray *cells = options->cells; // Array of input cells - const psArray *imageNames = options->imageNames; // Names of images to read - const psArray *maskNames = options->maskNames; // Names of masks to read - const psArray *varianceNames = options->varianceNames; // Names of variance maps to read - const psArray *covariances = options->covariances; // Covariance matrices (already read) + const psArray *imageNames = conv ? options->convImages : options->origImages; // Names of images to read + const psArray *maskNames = conv ? options->convMasks : options->origMasks; // Names of masks to read + const psArray *varianceNames = conv ? options->convVariances : options->origVariances; // Variance names + const psArray *covariances = conv ? options->convCovars : options->origCovars; // Covariance matrices PS_ASSERT_ARRAY_NON_NULL(cells, NULL); - PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, imageNames, NULL); - PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, maskNames, NULL); - PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, varianceNames, NULL); - PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, covariances, NULL); + if (imageNames) { + PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, imageNames, NULL); + } + if (maskNames) { + PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, maskNames, NULL); + } + if (varianceNames) { + PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, varianceNames, NULL); + } + if (covariances) { + PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, covariances, NULL); + } ppStackThreadData *stack = psAlloc(sizeof(ppStackThreadData)); // Thread data, to return @@ -87,19 +95,21 @@ continue; } - // Resolved names - psString imageResolved = pmConfigConvertFilename(imageNames->data[i], config, false, false); - psString maskResolved = pmConfigConvertFilename(maskNames->data[i], config, false, false); - psString varianceResolved = pmConfigConvertFilename(varianceNames->data[i], config, false, false); - stack->imageFits->data[i] = psFitsOpen(imageResolved, "r"); - stack->maskFits->data[i] = psFitsOpen(maskResolved, "r"); - stack->varianceFits->data[i] = psFitsOpen(varianceResolved, "r"); - psFree(imageResolved); - psFree(maskResolved); - psFree(varianceResolved); - if (!stack->imageFits->data[i] || !stack->maskFits->data[i] || !stack->varianceFits->data[i]) { - psError(PS_ERR_UNKNOWN, false, "Unable to open convolved files %s, %s, %s", - (char*)imageNames->data[i], (char*)maskNames->data[i], (char*)varianceNames->data[i]); - return NULL; - } + +// Open an image +#define IMAGE_OPEN(NAMES, FITS, INDEX) \ + if (NAMES) { \ + psString resolved = pmConfigConvertFilename((NAMES)->data[INDEX], config, false, false); \ + (FITS)->data[INDEX] = psFitsOpen(resolved, "r"); \ + if (!(FITS)->data[INDEX]) { \ + psError(PPSTACK_ERR_IO, false, "Unable to open file %s", (char*)(NAMES)->data[INDEX]); \ + psFree(resolved); \ + return NULL; \ + } \ + psFree(resolved); \ + } + + IMAGE_OPEN(imageNames, stack->imageFits, i); + IMAGE_OPEN(maskNames, stack->maskFits, i); + IMAGE_OPEN(varianceNames, stack->varianceFits, i); } @@ -116,5 +126,7 @@ } pmReadout *ro = pmReadoutAlloc(cell); // Readout for thread - ro->covariance = psMemIncrRefCounter(covariances->data[j]); + if (covariances) { + ro->covariance = psMemIncrRefCounter(covariances->data[j]); + } readouts->data[j] = ro; } @@ -140,10 +152,10 @@ psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // Recipe if (!recipe) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find recipe %s", PPSTACK_RECIPE); + psError(PPSTACK_ERR_CONFIG, false, "Unable to find recipe %s", PPSTACK_RECIPE); return NULL; } int rows = psMetadataLookupS32(NULL, recipe, "ROWS"); // Number of rows to read per chunk if (rows <= 0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, "ROWS is not set in the recipe."); + psError(PPSTACK_ERR_CONFIG, false, "ROWS is not set in the recipe."); return NULL; } @@ -186,11 +198,12 @@ psFits *varianceFits = stack->varianceFits->data[i]; // FITS file for variance - - int zMax = 0; + int zMax = 0; bool keepReading = false; - if (pmReadoutMore(ro, imageFits, 0, &zMax, rows, config)) { + + if (imageFits && pmReadoutMore(ro, imageFits, 0, &zMax, rows, config)) { keepReading = true; if (!pmReadoutReadChunk(ro, imageFits, 0, NULL, rows, overlap, config)) { - psError(PS_ERR_IO, false, "Unable to read chunk %d for file PPSTACK.INPUT %d", + psError(PPSTACK_ERR_IO, false, + "Unable to read chunk %d for file PPSTACK.INPUT %d", numChunk, i); *status = false; @@ -199,8 +212,9 @@ } - if (pmReadoutMoreMask(ro, maskFits, 0, &zMax, rows, config)) { + if (maskFits && pmReadoutMoreMask(ro, maskFits, 0, &zMax, rows, config)) { keepReading = true; if (!pmReadoutReadChunkMask(ro, maskFits, 0, NULL, rows, overlap, config)) { - psError(PS_ERR_IO, false, "Unable to read chunk %d for file PPSTACK.INPUT.MASK %d", + psError(PPSTACK_ERR_IO, false, + "Unable to read chunk %d for file PPSTACK.INPUT.MASK %d", numChunk, i); *status = false; @@ -209,8 +223,8 @@ } - if (pmReadoutMoreVariance(ro, varianceFits, 0, &zMax, rows, config)) { + if (varianceFits && pmReadoutMoreVariance(ro, varianceFits, 0, &zMax, rows, config)) { keepReading = true; if (!pmReadoutReadChunkVariance(ro, varianceFits, 0, NULL, rows, overlap, config)) { - psError(PS_ERR_IO, false, + psError(PPSTACK_ERR_IO, false, "Unable to read chunk %d for file PPSTACK.INPUT.VARIANCE %d", numChunk, i); @@ -263,5 +277,5 @@ { - psThreadTask *task = psThreadTaskAlloc("PPSTACK_INSPECT", 2); + psThreadTask *task = psThreadTaskAlloc("PPSTACK_INSPECT", 3); task->function = &ppStackInspect; psThreadTaskAdd(task); @@ -270,5 +284,5 @@ { - psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 3); + psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 6); task->function = &ppStackReadoutFinalThread; psThreadTaskAdd(task); Index: branches/tap_branches/ppStack/src/ppStackThread.h =================================================================== --- branches/tap_branches/ppStack/src/ppStackThread.h (revision 25900) +++ branches/tap_branches/ppStack/src/ppStackThread.h (revision 27838) @@ -35,5 +35,6 @@ ppStackThreadData *ppStackThreadDataSetup( const ppStackOptions *options, // Options - const pmConfig *config // Configuration + const pmConfig *config, // Configuration + bool conv // Use convolved products? );