Index: trunk/ippconfig/recipes/filerules-mef.mdc =================================================================== --- trunk/ippconfig/recipes/filerules-mef.mdc (revision 26075) +++ trunk/ippconfig/recipes/filerules-mef.mdc (revision 26076) @@ -263,4 +263,6 @@ PPSTACK.OUTPUT.MASK OUTPUT {OUTPUT}.mk.fits MASK COMP_MASK FPA TRUE NONE PPSTACK.OUTPUT.VARIANCE OUTPUT {OUTPUT}.wt.fits VARIANCE COMP_WT FPA TRUE NONE +PPSTACK.UNCONV.MASK OUTPUT {OUTPUT}.unconv.mask.fits MASK COMP_MASK FPA TRUE NONE +PPSTACK.UNCONV.VARIANCE OUTPUT {OUTPUT}.unconv.wt.fits VARIANCE COMP_WT FPA TRUE NONE PPSTACK.TARGET.PSF OUTPUT {OUTPUT}.target.psf PSF NONE CHIP TRUE NONE PPSTACK.CONV.KERNEL OUTPUT {OUTPUT}.{FILE.INDEX}.kernel SUBKERNEL NONE FPA TRUE NONE Index: trunk/ippconfig/recipes/filerules-simple.mdc =================================================================== --- trunk/ippconfig/recipes/filerules-simple.mdc (revision 26075) +++ trunk/ippconfig/recipes/filerules-simple.mdc (revision 26076) @@ -212,4 +212,6 @@ PPSTACK.OUTPUT.MASK OUTPUT {OUTPUT}.mask.fits MASK NONE FPA TRUE NONE PPSTACK.OUTPUT.VARIANCE OUTPUT {OUTPUT}.weight.fits VARIANCE NONE FPA TRUE NONE +PPSTACK.UNCONV.MASK OUTPUT {OUTPUT}.unconv.mask.fits MASK COMP_MASK FPA TRUE NONE +PPSTACK.UNCONV.VARIANCE OUTPUT {OUTPUT}.unconv.wt.fits VARIANCE COMP_WT FPA TRUE NONE PPSTACK.TARGET.PSF OUTPUT {OUTPUT}.target.psf PSF NONE CHIP TRUE NONE PPSTACK.CONV.KERNEL OUTPUT {OUTPUT}.{FILE.INDEX}.kernel SUBKERNEL NONE FPA TRUE NONE Index: trunk/ippconfig/recipes/filerules-split.mdc =================================================================== --- trunk/ippconfig/recipes/filerules-split.mdc (revision 26075) +++ trunk/ippconfig/recipes/filerules-split.mdc (revision 26076) @@ -229,4 +229,6 @@ PPSTACK.OUTPUT.MASK OUTPUT {OUTPUT}.mask.fits MASK COMP_MASK FPA TRUE NONE PPSTACK.OUTPUT.VARIANCE OUTPUT {OUTPUT}.wt.fits VARIANCE COMP_WT FPA TRUE NONE +PPSTACK.UNCONV.MASK OUTPUT {OUTPUT}.unconv.mask.fits MASK COMP_MASK FPA TRUE NONE +PPSTACK.UNCONV.VARIANCE OUTPUT {OUTPUT}.unconv.wt.fits VARIANCE COMP_WT FPA TRUE NONE PPSTACK.TARGET.PSF OUTPUT {OUTPUT}.target.psf PSF NONE CHIP TRUE NONE PPSTACK.CONV.KERNEL OUTPUT {OUTPUT}.{FILE.INDEX}.kernel SUBKERNEL NONE FPA TRUE NONE Index: trunk/ippconfig/recipes/ppStack.config =================================================================== --- trunk/ippconfig/recipes/ppStack.config (revision 26075) +++ trunk/ippconfig/recipes/ppStack.config (revision 26076) @@ -2,9 +2,10 @@ CONVOLVE BOOL TRUE # Convolve images when stacking? -ITER S32 1 # Number of rejection iterations +COMBINE.ITER F32 0.5 # Number of rejection iterations per input COMBINE.REJ F32 2.5 # Rejection threshold in combination (sigma) -COMBINE.SYS F32 0.03 # Relative systematic error in combination +COMBINE.SYS F32 0.1 # Relative systematic error in combination COMBINE.DISCARD F32 0.2 # Discard fraction for Olympic weighted mean -MASK.VAL STR MASK.VALUE,CONV.BAD,BURNTOOL # Mask value of input bad pixels +MASK.IN STR MASK.VALUE,CONV.BAD # Mask value of input bad pixels +MASK.SUSPECT STR SUSPECT,CONV.POOR # Mask value of suspect pixels MASK.BAD STR BLANK # Mask value to give bad pixels MASK.POOR STR CONV.POOR # Mask value to give poor pixels @@ -33,5 +34,5 @@ PHOT.FRAC F32 0.9 # Minimum fraction of good pixels -ZP.RADIUS F32 1.0 # Radius (pixels) for matching sources +ZP.RADIUS F32 4.0 # Radius (pixels) for matching sources ZP.ITER.1 S32 5 # Iterations for zero point; pass 1 ZP.ITER.2 S32 1000 # Iterations for zero point; pass 2 Index: trunk/ippconfig/recipes/ppSub.config =================================================================== --- trunk/ippconfig/recipes/ppSub.config (revision 26075) +++ trunk/ippconfig/recipes/ppSub.config (revision 26076) @@ -15,5 +15,5 @@ SYS.ERR F32 0.1 # Relative systematic error in images -MASK.IN STR MASK.VALUE,CONV.BAD # Mask value for input +MASK.VAL STR MASK.VALUE,CONV.BAD # Mask value for input MASK.BAD STR CONV.BAD # Mask value to give bad pixels MASK.POOR STR CONV.POOR # Mask value to give poor pixels Index: trunk/ppStack/src/ppStack.h =================================================================== --- trunk/ppStack/src/ppStack.h (revision 26075) +++ trunk/ppStack/src/ppStack.h (revision 26076) @@ -59,5 +59,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 @@ -83,5 +83,7 @@ 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 *addVariance, // Additional variance for rejection + bool safety, // Enable safety switch? + const psVector *norm // Normalisations to apply ); Index: trunk/ppStack/src/ppStackArguments.c =================================================================== --- trunk/ppStack/src/ppStackArguments.c (revision 26075) +++ trunk/ppStack/src/ppStackArguments.c (revision 26076) @@ -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); @@ -250,5 +250,5 @@ } - 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 +260,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); Index: trunk/ppStack/src/ppStackCamera.c =================================================================== --- trunk/ppStack/src/ppStackCamera.c (revision 26075) +++ trunk/ppStack/src/ppStackCamera.c (revision 26076) @@ -233,11 +233,11 @@ // Output image - pmFPA *fpa = pmFPAConstruct(config->camera, config->cameraName); // FPA to contain the output - if (!fpa) { + pmFPA *outFPA = pmFPAConstruct(config->camera, config->cameraName); // FPA to contain the output + if (!outFPA) { 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 + 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")); @@ -250,7 +250,6 @@ output->save = true; - if (!pmFPAAddSourceFromFormat(fpa, "Stack", output->format)) { + if (!pmFPAAddSourceFromFormat(outFPA, "Stack", output->format)) { psError(PS_ERR_UNKNOWN, false, "Unable to generate output FPA."); - psFree(fpa); return false; } @@ -294,4 +293,55 @@ targetPSF->save = true; } + +#if 1 + // Unconvolved stack + pmFPA *unconvFPA = pmFPAConstruct(config->camera, config->cameraName); // FPA to contain unconvolved output + if (!unconvFPA) { + psError(PS_ERR_UNEXPECTED_NULL, 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(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.UNCONV")); + return false; + } + if (unConv->type != PM_FPA_FILE_IMAGE) { + psError(PS_ERR_IO, true, "PPSTACK.UNCONV is not of type IMAGE"); + return false; + } + unConv->save = true; + + if (!pmFPAAddSourceFromFormat(unconvFPA, "Stack", unConv->format)) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate output FPA."); + return false; + } + + // Unconvolved mask + pmFPAfile *unconvMask = pmFPAfileDefineOutput(config, unconvFPA, "PPSTACK.UNCONV.MASK"); + if (!unconvMask) { + psError(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.UNCONV.MASK")); + return false; + } + if (unconvMask->type != PM_FPA_FILE_MASK) { + psError(PS_ERR_IO, 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(PS_ERR_IO, false, _("Unable to generate output file from PPSTACK.UNCONV.VARIANCE")); + return false; + } + if (unconvVariance->type != PM_FPA_FILE_VARIANCE) { + psError(PS_ERR_IO, true, "PPSTACK.UNCONV.VARIANCE is not of type VARIANCE"); + return false; + } + unconvVariance->save = true; + } +#endif // Output JPEGs Index: trunk/ppStack/src/ppStackCleanup.c =================================================================== --- trunk/ppStack/src/ppStackCleanup.c (revision 26075) +++ trunk/ppStack/src/ppStackCleanup.c (revision 26076) @@ -74,9 +74,9 @@ if (tempDelete) { - psString imageResolved = pmConfigConvertFilename(options->imageNames->data[i], + psString imageResolved = pmConfigConvertFilename(options->convImages->data[i], config, false, false); - psString maskResolved = pmConfigConvertFilename(options->maskNames->data[i], + psString maskResolved = pmConfigConvertFilename(options->convMasks->data[i], config, false, false); - psString varianceResolved = pmConfigConvertFilename(options->varianceNames->data[i], + psString varianceResolved = pmConfigConvertFilename(options->convVariances->data[i], config, false, false); if (unlink(imageResolved) == -1 || unlink(maskResolved) == -1 || Index: trunk/ppStack/src/ppStackCombineFinal.c =================================================================== --- trunk/ppStack/src/ppStackCombineFinal.c (revision 26075) +++ trunk/ppStack/src/ppStackCombineFinal.c (revision 26076) @@ -10,9 +10,16 @@ #include "ppStackLoop.h" -bool ppStackCombineFinal(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config) +//#define TESTING // Enable test output + +bool ppStackCombineFinal(pmReadout *target, ppStackThreadData *stack, psArray *covariances, + ppStackOptions *options, pmConfig *config, bool safe, bool normalise) { psAssert(stack, "Require stack"); psAssert(options, "Require options"); psAssert(config, "Require configuration"); + + if (!target->mask) { + target->mask = psImageAlloc(target->image->numCols, target->image->numRows, PS_TYPE_IMAGE_MASK); + } stack->lastScan = 0; // Reset read @@ -30,9 +37,12 @@ } - // call: ppStackReadoutFinal(config, outRO, readouts, rejected) + // call: ppStackReadoutFinal(config, target, readouts, rejected) psThreadJob *job = psThreadJobAlloc("PPSTACK_FINAL_COMBINE"); // Job to start + psArrayAdd(job->args, 1, target); psArrayAdd(job->args, 1, thread); 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); @@ -48,31 +58,43 @@ // 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; - continue; + if (covariances) { + double sumWeights = 0.0; // Sum of weights + for (int i = 0; i < options->num; i++) { + if (options->inputMask->data.U8[i]) { + psFree(covariances->data[i]); + covariances->data[i] = NULL; + continue; + } + psKernel *covar = covariances->data[i]; // Covariance matrix + if (!covar) { + continue; + } + float weight = options->weightings->data.F32[i]; // Weight to apply + psBinaryOp(covar->image, covar->image, "*", psScalarAlloc(weight, PS_TYPE_F32)); + sumWeights += weight; } - psKernel *covar = options->covariances->data[i]; // Covariance matrix - if (!covar) { - continue; + if (sumWeights > 0.0) { + target->covariance = psImageCovarianceSum(covariances); + psBinaryOp(target->covariance->image, target->covariance->image, "/", + psScalarAlloc(sumWeights, PS_TYPE_F32)); + psImageCovarianceTransfer(target->variance, target->covariance); } - 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); + } else { + target->covariance = psImageCovarianceNone(); } #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, target->image, config); + psStringSubstitute(&name, "mask", "image"); + ppStackWriteImage(name, NULL, target->mask, config); + psStringSubstitute(&name, "variance", "mask"); + ppStackWriteImage(name, NULL, target->variance, config); + psFree(name); + + pmStackVisualPlotTestImage(target->image, "combined_image_final.fits"); #endif Index: trunk/ppStack/src/ppStackCombineInitial.c =================================================================== --- trunk/ppStack/src/ppStackCombineInitial.c (revision 26075) +++ trunk/ppStack/src/ppStackCombineInitial.c (revision 26076) @@ -9,4 +9,6 @@ #include "ppStack.h" #include "ppStackLoop.h" + +//#define TESTING // Enable test output bool ppStackCombineInitial(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config) @@ -60,4 +62,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 +68,5 @@ } options->inspect->data[i] = psArrayAllocEmpty(numChunk); + options->rejected->data[i] = psArrayAllocEmpty(numChunk); } psThreadJob *job; // Completed job @@ -71,9 +75,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 +91,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: trunk/ppStack/src/ppStackCombinePrepare.c =================================================================== --- trunk/ppStack/src/ppStackCombinePrepare.c (revision 26075) +++ trunk/ppStack/src/ppStackCombinePrepare.c (revision 26076) @@ -33,4 +33,8 @@ pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT"); // Output cell options->outRO = pmReadoutAlloc(outCell); // Output readout + + pmCell *unconvCell = pmFPAfileThisCell(config->files, view, "PPSTACK.UNCONV"); // Unconvolved cell + options->unconvRO = pmReadoutAlloc(unconvCell); // Unconvolved readout + psFree(view); @@ -39,4 +43,5 @@ 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."); @@ -44,4 +49,11 @@ } + options->unconvRO->image = psImageCopy(NULL, options->outRO->image, PS_TYPE_F32); +// options->unconvRO->mask = psImageCopy(NULL, options->outRO->mask, PS_TYPE_IMAGE_MASK); + options->unconvRO->col0 = options->outRO->col0; + options->unconvRO->row0 = options->outRO->row0; + + + return true; } Index: trunk/ppStack/src/ppStackConvolve.c =================================================================== --- trunk/ppStack/src/ppStackConvolve.c (revision 26075) +++ trunk/ppStack/src/ppStackConvolve.c (revision 26076) @@ -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 @@ -79,4 +86,5 @@ // 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); @@ -85,5 +93,11 @@ continue; } - options->covariances->data[i] = psMemIncrRefCounter(readout->covariance); + 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 +128,10 @@ pmHDU *hdu = readout->parent->parent->parent->hdu; // HDU for convolved image assert(hdu); - ppStackWriteImage(options->imageNames->data[i], hdu->header, readout->image, config); + ppStackWriteImage(options->convImages->data[i], hdu->header, readout->image, config); psMetadata *maskHeader = psMetadataCopy(NULL, hdu->header); // Copy of header, for mask pmConfigMaskWriteHeader(config, maskHeader); - ppStackWriteImage(options->maskNames->data[i], maskHeader, readout->mask, config); + ppStackWriteImage(options->convMasks->data[i], maskHeader, readout->mask, config); psFree(maskHeader); - psImageCovarianceTransfer(readout->variance, readout->covariance); - ppStackWriteImage(options->varianceNames->data[i], hdu->header, readout->variance, config); + ppStackWriteImage(options->convVariances->data[i], hdu->header, readout->variance, config); #ifdef TESTING { @@ -129,4 +142,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 @@ -149,5 +177,4 @@ psFree(rng); - psFree(options->norm); options->norm = NULL; psFree(options->sourceLists); options->sourceLists = NULL; psFree(options->psf); options->psf = NULL; @@ -211,5 +238,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; } @@ -233,4 +260,12 @@ } + // 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: trunk/ppStack/src/ppStackFiles.c =================================================================== --- trunk/ppStack/src/ppStackFiles.c (revision 26075) +++ trunk/ppStack/src/ppStackFiles.c (revision 26076) @@ -22,4 +22,5 @@ /// Output files for the combination static char *filesCombine[] = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.VARIANCE", + "PPSTACK.UNCONV", "PPSTACK.UNCONV.MASK", "PPSTACK.UNCONV.VARIANCE", "PPSTACK.OUTPUT.JPEG1", "PPSTACK.OUTPUT.JPEG2", NULL }; Index: trunk/ppStack/src/ppStackLoop.c =================================================================== --- trunk/ppStack/src/ppStackLoop.c (revision 26075) +++ trunk/ppStack/src/ppStackLoop.c (revision 26076) @@ -56,5 +56,5 @@ // 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."); @@ -62,5 +62,4 @@ return false; } - psFree(options->cells); options->cells = NULL; // Prepare for combination @@ -98,5 +97,5 @@ // Final combination psTrace("ppStack", 2, "Final stack of convolved images....\n"); - if (!ppStackCombineFinal(stack, options, config)) { + if (!ppStackCombineFinal(options->outRO, stack, options->convCovars, options, config, false, false)) { psError(PS_ERR_UNKNOWN, false, "Unable to perform final combination."); psFree(stack); @@ -106,5 +105,4 @@ psLogMsg("ppStack", PS_LOG_INFO, "Stage 5: Final Stack: %f sec", psTimerClear("PPSTACK_STEPS")); ppStackMemDump("final"); - // Clean up @@ -121,4 +119,28 @@ psFree(stack); +#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(PS_ERR_IO, false, "Unable to initialise stack threads."); + psFree(options); + return false; + } + psTrace("ppStack", 2, "Stack of unconvolved images....\n"); + if (!ppStackCombineFinal(options->unconvRO, stack, options->origCovars, options, config, false, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to perform unconvolved combination."); + psFree(stack); + psFree(options); + return false; + } + psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: Unconvolved Stack: %f sec", psTimerClear("PPSTACK_STEPS")); + ppStackMemDump("unconv"); + + psFree(stack); + } + psFree(options->cells); options->cells = NULL; +#endif // Photometry @@ -129,7 +151,6 @@ return false; } - psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: Photometry Analysis: %f sec", psTimerClear("PPSTACK_STEPS")); + psLogMsg("ppStack", PS_LOG_INFO, "Stage 8: Photometry Analysis: %f sec", psTimerClear("PPSTACK_STEPS")); ppStackMemDump("photometry"); - // Finish up @@ -140,8 +161,9 @@ return false; } - psLogMsg("ppStack", PS_LOG_INFO, "Stage 8: Final output: %f sec", psTimerClear("PPSTACK_STEPS")); + psLogMsg("ppStack", PS_LOG_INFO, "Stage 9: Final output: %f sec", psTimerClear("PPSTACK_STEPS")); ppStackMemDump("finish"); psFree(options); + return true; } Index: trunk/ppStack/src/ppStackLoop.h =================================================================== --- trunk/ppStack/src/ppStackLoop.h (revision 26075) +++ trunk/ppStack/src/ppStackLoop.h (revision 26076) @@ -56,7 +56,11 @@ // Final combination bool ppStackCombineFinal( + pmReadout *target, // Target readout ppStackThreadData *stack, // Stack + psArray *covariances, // Covariances ppStackOptions *options, // Options - pmConfig *config // Configuration + pmConfig *config, // Configuration + bool safe, // Allow safe combination? + bool norm // Normalise images? ); Index: trunk/ppStack/src/ppStackMatch.c =================================================================== --- trunk/ppStack/src/ppStackMatch.c (revision 26075) +++ trunk/ppStack/src/ppStackMatch.c (revision 26076) @@ -167,4 +167,5 @@ ) { +#if 1 bool mdok; // Status of metadata lookups @@ -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 @@ -257,7 +262,7 @@ 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)) { + 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(PS_ERR_IO, false, "Unable to read previously produced image."); return false; @@ -378,4 +383,7 @@ } #endif + + fprintf(stderr, "vf = %f\n", psImageCovarianceFactor(readout->covariance)); + if (threads > 0) { @@ -517,4 +525,5 @@ psFree(iter); options->matchChi2->data.F32[index] = sum / (psImageCovarianceFactor(readout->covariance) * num); + fprintf(stderr, "chi2 = %f ; vf = %f\n", sum/num, psImageCovarianceFactor(readout->covariance)); } @@ -543,13 +552,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)); + } } @@ -569,4 +578,6 @@ options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) * psImageCovarianceFactor(readout->covariance)); + psLogMsg("ppStack", PS_LOG_INFO, "Weighting for image %d is %f\n", + index, options->weightings->data.F32[index]); psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "PPSTACK.WEIGHTING", 0, "Weighting by 1/noise^2 for stack", options->weightings->data.F32[index]); Index: trunk/ppStack/src/ppStackOptions.c =================================================================== --- trunk/ppStack/src/ppStackOptions.c (revision 26075) +++ trunk/ppStack/src/ppStackOptions.c (revision 26076) @@ -12,7 +12,11 @@ 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->inputMask); @@ -24,6 +28,7 @@ psFree(options->matchChi2); psFree(options->weightings); - psFree(options->covariances); + psFree(options->convCovars); psFree(options->outRO); + psFree(options->unconvRO); psFree(options->inspect); psFree(options->rejected); @@ -39,7 +44,11 @@ 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->psf = NULL; @@ -54,6 +63,7 @@ options->matchChi2 = NULL; options->weightings = NULL; - options->covariances = NULL; + options->convCovars = NULL; options->outRO = NULL; + options->unconvRO = NULL; options->inspect = NULL; options->rejected = NULL; Index: trunk/ppStack/src/ppStackOptions.h =================================================================== --- trunk/ppStack/src/ppStackOptions.h (revision 26075) +++ trunk/ppStack/src/ppStackOptions.h (revision 26076) @@ -11,5 +11,7 @@ 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 // Prepare @@ -26,7 +28,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 *unconvRO; // Unconvolved readout psArray *inspect; // Array of arrays of pixels to inspect // Rejection Index: trunk/ppStack/src/ppStackReadout.c =================================================================== --- trunk/ppStack/src/ppStackReadout.c (revision 26075) +++ trunk/ppStack/src/ppStackReadout.c (revision 26076) @@ -25,11 +25,9 @@ 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, addVariance); thread->busy = false; - return inspect ? true : false; + return job->results ? true : false; } @@ -39,16 +37,19 @@ 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 + pmReadout *target = args->data[0]; // Output readout + ppStackThread *thread = args->data[1]; // Thread + 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 *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, target, thread->readouts, mask, rejected, + weightings, addVariance, safety, norm); // Status of operation thread->busy = false; @@ -63,24 +64,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; @@ -104,5 +116,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 @@ -114,6 +126,8 @@ int kernelSize = 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, recipe, "MASK.IN"); // 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 @@ -149,6 +163,6 @@ } - if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, kernelSize, iter, - combineRej, combineSys, combineDiscard, true, useVariance, safe, false)) { + if (!pmStackCombine(outRO, 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 +170,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 +183,5 @@ } inspect->data[i] = psMemIncrRefCounter(data->inspect); + reject->data[i] = psMemIncrRefCounter(data->reject); } psFree(stack); @@ -173,5 +189,9 @@ sectionNum++; - return inspect; + psArray *results = psArrayAlloc(2); // Array of results + results->data[0] = inspect; + results->data[1] = reject; + + return results; } @@ -180,5 +200,5 @@ bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, const psArray *readouts, const psVector *mask, const psArray *rejected, const psVector *weightings, - const psVector *addVariance) + const psVector *addVariance, bool safety, const psVector *norm) { assert(config); @@ -196,13 +216,11 @@ 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 - psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in + psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.IN"); // 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,22 +229,15 @@ 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 + if (mask->data.U8[i]) { + // Image completely rejected continue; } @@ -252,9 +263,14 @@ 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, stack, maskVal | maskBad, maskSuspect, maskBad, 0, iter, combineRej, + combineSys, combineDiscard, useVariance, safe, true)) { psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts."); psFree(stack); Index: trunk/ppStack/src/ppStackReject.c =================================================================== --- trunk/ppStack/src/ppStackReject.c (revision 26075) +++ trunk/ppStack/src/ppStackReject.c (revision 26076) @@ -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 @@ -53,4 +52,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)) { @@ -92,15 +92,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 +107,8 @@ } #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 psFree(options->inspect->data[i]); @@ -127,5 +127,4 @@ "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; @@ -134,6 +133,25 @@ } - // Images without a list of rejected pixels (the list may be empty) are rejected completely - options->rejected->data[i] = reject; + if (reject) { + // Add to list of pixels already rejected + reject = psPixelsConcatenate(reject, options->rejected->data[i]); + options->rejected->data[i] = psPixelsDuplicates(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 +161,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")); Index: trunk/ppStack/src/ppStackSetup.c =================================================================== --- trunk/ppStack/src/ppStackSetup.c (revision 26075) +++ trunk/ppStack/src/ppStackSetup.c (revision 26076) @@ -72,7 +72,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,9 +81,30 @@ 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)) { Index: trunk/ppStack/src/ppStackSources.c =================================================================== --- trunk/ppStack/src/ppStackSources.c (revision 26075) +++ trunk/ppStack/src/ppStackSources.c (revision 26076) @@ -188,4 +188,21 @@ } +#if 0 + // Position offsets + { + psArray *offsets = pmSourceMatchRelastro(matches, num, tol, iter1, starRej1, + iter2, starRej2, starLimit); // Shifts for each image + if (!offsets) { + psError(PS_ERR_UNKNOWN, false, "Unable to measure offsets"); + return false; + } + for (int i = 0; i < num; i++) { + psVector *offset = offsets->data[i]; + fprintf(stderr, "Offset %d: %f,%f\n", i, offset->data.F32[0], offset->data.F32[1]); + } + psFree(offsets); + } +#endif + #ifdef TESTING dumpMatches("source_mags.dat", num, matches, zp, trans); Index: trunk/ppStack/src/ppStackThread.c =================================================================== --- trunk/ppStack/src/ppStackThread.c (revision 26075) +++ trunk/ppStack/src/ppStackThread.c (revision 26076) @@ -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(PS_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; } @@ -186,8 +198,8 @@ 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)) { @@ -199,5 +211,5 @@ } - 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)) { @@ -209,5 +221,5 @@ } - 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)) { @@ -263,5 +275,5 @@ { - psThreadTask *task = psThreadTaskAlloc("PPSTACK_INSPECT", 2); + psThreadTask *task = psThreadTaskAlloc("PPSTACK_INSPECT", 3); task->function = &ppStackInspect; psThreadTaskAdd(task); @@ -270,5 +282,5 @@ { - psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 3); + psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 6); task->function = &ppStackReadoutFinalThread; psThreadTaskAdd(task); Index: trunk/ppStack/src/ppStackThread.h =================================================================== --- trunk/ppStack/src/ppStackThread.h (revision 26075) +++ trunk/ppStack/src/ppStackThread.h (revision 26076) @@ -35,5 +35,6 @@ ppStackThreadData *ppStackThreadDataSetup( const ppStackOptions *options, // Options - const pmConfig *config // Configuration + const pmConfig *config, // Configuration + bool conv // Use convolved products? ); Index: trunk/ppSub/src/ppSubArguments.c =================================================================== --- trunk/ppSub/src/ppSubArguments.c (revision 26075) +++ trunk/ppSub/src/ppSubArguments.c (revision 26076) @@ -86,4 +86,5 @@ psMetadataAddS32(arguments, PS_LIST_TAIL, "-convolve", 0, "Image to convolve [1 or 2]", 0); psMetadataAddBool(arguments, PS_LIST_TAIL, "-photometry", 0, "Perform photometry?", NULL); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp", 0, "Zero point for photometry", NAN); psMetadataAddBool(arguments, PS_LIST_TAIL, "-inverse", 0, "Generate inverse subtractions?", NULL); psMetadataAddBool(arguments, PS_LIST_TAIL, "-visual", 0, "Show diagnostic plots", NULL); Index: trunk/psModules/src/camera/pmFPAMaskWeight.c =================================================================== --- trunk/psModules/src/camera/pmFPAMaskWeight.c (revision 26075) +++ trunk/psModules/src/camera/pmFPAMaskWeight.c (revision 26076) @@ -444,5 +444,6 @@ } - return true; + return psMetadataAddF32(readout->analysis, PS_LIST_TAIL, PM_READOUT_ANALYSIS_RENORM, 0, + "Renormalisation of variance", PS_SQR(correction)); } Index: trunk/psModules/src/camera/pmFPAMaskWeight.h =================================================================== --- trunk/psModules/src/camera/pmFPAMaskWeight.h (revision 26075) +++ trunk/psModules/src/camera/pmFPAMaskWeight.h (revision 26076) @@ -15,4 +15,6 @@ /// @addtogroup Camera Camera Layout /// @{ + +#define PM_READOUT_ANALYSIS_RENORM "READOUT.RENORM" // Name on analysis metadata for renormalisation /// Set a temporary readout mask using CELL.SATURATION and CELL.BAD Index: trunk/psModules/src/imcombine/pmStack.c =================================================================== --- trunk/psModules/src/imcombine/pmStack.c (revision 26075) +++ trunk/psModules/src/imcombine/pmStack.c (revision 26076) @@ -33,7 +33,8 @@ #define NUM_DIRECT_STDEV 5 // For less than this number of values, measure stdev directly + //#define TESTING // Enable test output -//#define TEST_X 1085 // x coordinate to examine -//#define TEST_Y 3371 // y coordinate to examine +//#define TEST_X 3122-1 // x coordinate to examine +//#define TEST_Y 1028-1 // y coordinate to examine @@ -42,9 +43,9 @@ typedef struct { psVector *pixels; // Pixel values - psVector *masks; // Pixel masks psVector *variances; // Pixel variances psVector *weights; // Pixel weightings psVector *sources; // Pixel sources (which image did they come from?) psVector *limits; // Rejection limits + psVector *suspects; // Pixel is suspect? psVector *sort; // Buffer for sorting (to get a robust estimator of the standard dev) } combineBuffer; @@ -53,9 +54,9 @@ { psFree(buffer->pixels); - psFree(buffer->masks); psFree(buffer->variances); psFree(buffer->weights); psFree(buffer->sources); psFree(buffer->limits); + psFree(buffer->suspects); psFree(buffer->sort); return; @@ -69,9 +70,9 @@ buffer->pixels = psVectorAlloc(numImages, PS_TYPE_F32); - buffer->masks = psVectorAlloc(numImages, PS_TYPE_VECTOR_MASK); buffer->variances = psVectorAlloc(numImages, PS_TYPE_F32); buffer->weights = psVectorAlloc(numImages, PS_TYPE_F32); buffer->sources = psVectorAlloc(numImages, PS_TYPE_U16); buffer->limits = psVectorAlloc(numImages, PS_TYPE_F32); + buffer->suspects = psVectorAlloc(numImages, PS_TYPE_U8); buffer->sort = psVectorAlloc(numImages, PS_TYPE_F32); return buffer; @@ -143,23 +144,16 @@ float *stdev, // Standard deviation value, to return const psVector *values, // Values to combine - const psVector *masks, // Mask to apply psVector *sortBuffer // Buffer for sorting ) { assert(values); - assert(!masks || values->n == masks->n); assert(values->type.type == PS_TYPE_F32); - assert(!masks || masks->type.type == PS_TYPE_VECTOR_MASK); assert(sortBuffer && sortBuffer->nalloc >= values->n && sortBuffer->type.type == PS_TYPE_F32); - // Need to filter out clipped values - int num = 0; // Number of valid values - for (int i = 0; i < values->n; i++) { - if (!masks || !masks->data.PS_TYPE_VECTOR_MASK_DATA[i]) { - sortBuffer->data.F32[num++] = values->data.F32[i]; - } - } - sortBuffer->n = num; - if (!psVectorSortInPlace(sortBuffer)) { + int num = values->n; // Number of values + sortBuffer = psVectorSortIndex(sortBuffer, values); + if (!sortBuffer) { + *median = NAN; + *stdev = NAN; return false; } @@ -167,14 +161,15 @@ if (num == 3) { // Attempt to measure standard deviation with only three values (and one of those possibly corrupted) - *median = sortBuffer->data.F32[1]; + *median = values->data.F32[sortBuffer->data.S32[1]]; if (stdev) { - float diff1 = sortBuffer->data.F32[0] - *median; - float diff2 = sortBuffer->data.F32[2] - *median; + float diff1 = values->data.F32[sortBuffer->data.S32[0]] - *median; + float diff2 = values->data.F32[sortBuffer->data.S32[2]] - *median; // This factor of sqrt(2) might not be exact, but it's about right *stdev = M_SQRT2 * PS_MIN(fabsf(diff1), fabsf(diff2)); } } else { - *median = num % 2 ? sortBuffer->data.F32[num / 2] : - (sortBuffer->data.F32[num / 2 - 1] + sortBuffer->data.F32[num / 2]) / 2.0; + *median = num % 2 ? values->data.F32[sortBuffer->data.S32[num / 2]] : + (values->data.F32[sortBuffer->data.S32[num / 2 - 1]] + + values->data.F32[sortBuffer->data.S32[num / 2]]) / 2.0; if (stdev) { if (num <= NUM_DIRECT_STDEV) { @@ -182,11 +177,11 @@ double sum = 0.0; for (int i = 0; i < num; i++) { - sum += PS_SQR(sortBuffer->data.F32[i] - *median); + sum += PS_SQR(values->data.F32[sortBuffer->data.S32[i]] - *median); } *stdev = sqrt(sum / (double)(num - 1)); } else { // Standard deviation from the interquartile range - *stdev = 0.74 * (sortBuffer->data.F32[(int)(0.75 * num)] - - sortBuffer->data.F32[(int)(0.25 * num)]); + *stdev = 0.74 * (values->data.F32[sortBuffer->data.S32[(int)(0.75 * num)]] - + values->data.F32[sortBuffer->data.S32[(int)(0.25 * num)]]); } } @@ -195,72 +190,13 @@ return true; } - -#if 0 -// Return the weighted median for the pixels -// This does not appear to produce as clean images as the weighted Olympic mean -static float combinationWeightedMedian(const psVector *values, // Values to combine - const psVector *weights, // Weights to combine - const psVector *masks, // Mask to apply - psVector *sortBuffer // Buffer for sorting - ) -{ - double sumWeight = 0.0; // Sum of weights - for (int j = 0; j < values->n; j++) { - if (masks->data.PS_TYPE_VECTOR_MASK_DATA[j]) { - continue; - } - sumWeight += weights->data.F32[j]; - } - - sortBuffer = psVectorSortIndex(sortBuffer, values); - double target = sumWeight / 2.0; // Target weight - - int dominant = -1; // Index of dominant value, if any - double cumulativeWeight = 0.0; // Sum of weights - for (int j = 0; j < values->n; j++) { - if (masks->data.PS_TYPE_VECTOR_MASK_DATA[j]) { - continue; - } - int index = sortBuffer->data.S32[j]; // Index of value of interest - float weight = weights->data.F32[index]; // Weight for value of interest - if (weight >= target) { - // Get the weighted median of the rest - dominant = index; - sumWeight -= weight; - target = sumWeight / 2.0; - continue; - } - cumulativeWeight += weight; - if (cumulativeWeight >= target) { - float median = values->data.F32[index]; // Weighted median median - if (dominant != -1) { - // In the case that a single value contains a disproportionate weight compared to the rest, - // we use a weighted mean between that dominant value and the weighted median of the rest. - return (values->data.F32[dominant] * weights->data.F32[dominant] + median * sumWeight) / - (weights->data.F32[dominant] + sumWeight); - } - return median; - } - } - - return NAN; -} -#endif // Return the weighted Olympic mean for the pixels static float combinationWeightedOlympic(const psVector *values, // Values to combine const psVector *weights, // Weights to combine - const psVector *masks, // Mask to apply float frac, // Fraction to discard psVector *sortBuffer // Buffer for sorting ) { - int numGood = 0; // Number of good values - for (int i = 0; i < values->n; i++) { - if (masks->data.PS_TYPE_VECTOR_MASK_DATA[i]) { - continue; - } - numGood++; - } + int numGood = values->n; // Number of good values int numBad = frac * numGood + 0.5; // Number of bad values @@ -272,7 +208,4 @@ for (int i = 0, j = 0; i < values->n; i++) { int index = sortBuffer->data.S32[i]; // Index of interest - if (masks->data.PS_TYPE_VECTOR_MASK_DATA[index]) { - continue; - } j++; if (j > high) { @@ -290,9 +223,15 @@ // Mark a pixel for inspection -static inline void combineInspect(const psArray *inputs, // Stack data - int x, int y, // Pixel - int source // Source image index - ) -{ +// Value in pixel doesn't seem to agree with the stack, so need to look closer +static inline void combineMarkInspect(const psArray *inputs, // Stack data + int x, int y, // Pixel + int source // Source image index + ) +{ +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Marking image %d for inspection\n", source); + } +#endif pmStackData *data = inputs->data[source]; // Stack data of interest if (!data) { @@ -304,45 +243,62 @@ } -// Given a stack of images, combine with optional rejection. -// Pixels in the stack that are rejected are marked for subsequent inspection -static void combinePixels(psImage *image, // Combined image, for output - psImage *mask, // Combined mask, for output - psImage *variance, // Combined variance map, for output - const psArray *inputs, // Stack data - const psVector *weights, // Global (single value) weights for data, or NULL - const psVector *addVariance, // Additional variance for data - const psVector *reject, // Indices of pixels to reject, or NULL - int x, int y, // Coordinates of interest; frame of output image - psImageMaskType maskVal, // Value to mask - psImageMaskType bad, // Value to give bad pixels - int numIter, // Number of rejection iterations - float rej, // Number of standard deviations at which to reject - float sys, // Relative systematic error - float discard,// Fraction of values to discard (Olympic weighted mean) - bool useVariance, // Use variance for rejection when combining? - bool safe, // Combine safely? - bool rejectInspect, // Reject values marked for inspection from combination? - combineBuffer *buffer // Buffer for combination; to avoid multiple allocations - ) +// Mark a pixel for rejection +// Cannot possibly inspect this pixel and confirm that it's good. +// e.g., Only a single input +static inline void combineMarkReject(const psArray *inputs, // Stack data + int x, int y, // Pixel + int source // Source image index + ) +{ +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Marking pixel image %d for rejection\n", source); + } +#endif + pmStackData *data = inputs->data[source]; // Stack data of interest + if (!data) { + psWarning("Can't find input data for source %d", source); + return; + } + data->reject = psPixelsAdd(data->reject, data->reject->nalloc, x, y); + return; +} + + +// Extract vectors for simple combination/rejection operations +static void combineExtract(int *num, // Number of good pixels + bool *suspect, // Any suspect pixels? + combineBuffer *buffer, // Buffer with vectors + psImage *image, // Combined image, for output + psImage *mask, // Combined mask, for output + psImage *variance, // Combined variance map, for output + const psArray *inputs, // Stack data + const psVector *weights, // Global (single value) weights for data, or NULL + const psVector *addVariance, // Additional variance for data + const psVector *reject, // Indices of pixels to reject, or NULL + int x, int y, // Coordinates of interest; frame of output image + psImageMaskType maskVal, // Value to mask + psImageMaskType maskSuspect // Value to suspect + ) { // Rudimentary error checking + assert(buffer); assert(image); assert(mask); assert(inputs); - assert(numIter >= 0); - assert(buffer); - assert(addVariance); - assert((useVariance && variance) || !useVariance); psVector *pixelData = buffer->pixels; // Values for the pixel of interest - psVector *pixelMasks = buffer->masks; // Masks for the pixel of interest psVector *pixelVariances = variance ? buffer->variances : NULL; // Variances for the pixel of interest psVector *pixelWeights = buffer->weights; // Image weights for the pixel of interest psVector *pixelSources = buffer->sources; // Sources for the pixel of interest psVector *pixelLimits = buffer->limits; // Limits for the pixel of interest - psVector *sort = buffer->sort; // Sort buffer + psVector *pixelSuspects = buffer->suspects; // Is the pixel suspect? + + if (suspect) { + *suspect = false; + } // Extract the pixel and mask data - int num = 0; // Number of good images + int numGood = 0; // Number of good pixels for (int i = 0, j = 0; i < inputs->n; i++) { // Check if this pixel has been rejected. Assumes that the rejection pixel list is sorted --- it @@ -364,45 +320,69 @@ } + pixelSuspects->data.U8[numGood] = mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & maskSuspect ? + true : false; + psImage *image = data->readout->image; // Image of interest psImage *variance = data->readout->variance; // Variance map of interest - pixelData->data.F32[num] = image->data.F32[yIn][xIn]; + pixelData->data.F32[numGood] = image->data.F32[yIn][xIn]; if (variance) { - pixelVariances->data.F32[num] = variance->data.F32[yIn][xIn] * addVariance->data.F32[i]; - } - pixelWeights->data.F32[num] = data->weight; - pixelSources->data.U16[num] = i; - num++; - } - pixelData->n = num; + pixelVariances->data.F32[numGood] = variance->data.F32[yIn][xIn] * addVariance->data.F32[i]; + } + pixelWeights->data.F32[numGood] = data->weight; + pixelSources->data.U16[numGood] = i; + numGood++; + } + pixelData->n = numGood; if (variance) { - pixelVariances->n = num; - } - pixelWeights->n = num; - pixelSources->n = num; - pixelLimits->n = num; + pixelVariances->n = numGood; + } + pixelWeights->n = numGood; + pixelSources->n = numGood; + pixelLimits->n = numGood; + pixelSuspects->n = numGood; + *num = numGood; #ifdef TESTING if (x == TEST_X && y == TEST_Y) { - for (int i = 0; i < num; i++) { - fprintf(stderr, "Input %d (%" PRIu16 "): %f %f %f\n", + for (int i = 0; i < numGood; i++) { + fprintf(stderr, "Input %d (%" PRIu16 "): %f %f (%f) %f %d\n", i, pixelSources->data.U16[i], pixelData->data.F32[i], pixelVariances->data.F32[i], - pixelWeights->data.F32[i]); - } - } -#endif - - // The sensible thing to do varies according to how many good pixels there are. + addVariance->data.F32[i], pixelWeights->data.F32[i], pixelSuspects->data.U8[i]); + } + } +#endif + + return; +} + + +// Combine pixels +static void combinePixels(psImage *image, // Combined image, for output + psImage *mask, // Combined mask, for output + psImage *variance, // Combined variance map, for output + int num, // Number of good pixels + combineBuffer *buffer, // Buffer with vectors + int x, int y, // Coordinates of interest; frame of output image + psImageMaskType bad, // Value for bad pixels + bool safe // Safe combination? + ) +{ + psVector *pixelData = buffer->pixels; // Values for the pixel of interest + psVector *pixelVariances = variance ? buffer->variances : NULL; // Variances for the pixel of interest + psVector *pixelWeights = buffer->weights; // Image weights for the pixel of interest + // Default option is that the pixel is bad float imageValue = NAN, varianceValue = NAN; // Value for combined image and variance map psImageMaskType maskValue = bad; // Value for combined mask switch (num) { - case 0: - // Nothing to combine: it's bad -#ifdef TESTING - if (x == TEST_X && y == TEST_Y) { - fprintf(stderr, "No inputs to combine, pixel is bad.\n"); - } -#endif - break; + case 0: { + // Nothing to combine: it's bad +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "No inputs to combine, pixel is bad.\n"); + } +#endif + break; + } case 1: { // Accept the single pixel unless we have to be safe @@ -427,7 +407,6 @@ } case 2: { - // Accept the mean of the pixels only if we're going to reject based on the variance, or we're not - // playing safe - if (useVariance || !safe) { + // Automatically accept the mean of the pixels only if we're not playing safe + if (!safe) { float mean, var; // Mean and variance from combination if (combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) { @@ -437,43 +416,20 @@ #ifdef TESTING if (x == TEST_X && y == TEST_Y) { - fprintf(stderr, "Two inputs to combine using variance/unsafe --> %f %f\n", - mean, var); + fprintf(stderr, "Two inputs to combine using unsafe --> %f %f\n", mean, var); } #endif } } - if (useVariance && numIter > 0) { - // Use variance to check that the two are consistent - float diff = 0.5 * (pixelData->data.F32[0] - pixelData->data.F32[1]); // Mean flux difference - float var1 = pixelVariances->data.F32[0]; // Variance of first - float var2 = pixelVariances->data.F32[1]; // Variance of second - // Systematic error contributes to the rejection level - var1 += PS_SQR(sys * pixelData->data.F32[0]); - var2 += PS_SQR(sys * pixelData->data.F32[1]); - - float sigma2 = var1 + var2; // Combined variance - if (PS_SQR(diff) > PS_SQR(rej) * sigma2) { - // Not consistent: mark both for inspection - if (rejectInspect) { - imageValue = NAN; - varianceValue = NAN; - maskValue = bad; - } else { - combineInspect(inputs, x, y, pixelSources->data.U16[0]); - combineInspect(inputs, x, y, pixelSources->data.U16[1]); - } -#ifdef TESTING - if (x == TEST_X && y == TEST_Y) { - fprintf(stderr, "Both pixels marked for inspection (%f > %f x %f\n)", - diff, rej, sqrtf(sigma2)); - } -#endif +#ifdef TESTING + else { + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Two inputs to combine, safety on, pixel is bad\n"); } } +#endif break; } default: { - // Record the value derived with no clipping, because pixels rejected using the harsh clipping - // applied in the first pass might later be accepted. + // Can combine without too much worrying float mean, var; // Mean and variance of the combination if (!combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) { @@ -488,139 +444,4 @@ } #endif - - // Prepare for clipping iteration - if (numIter > 0) { - pixelMasks->n = num; - psVectorInit(pixelMasks, 0); - if (useVariance) { - // Convert to rejection limits --- saves doing it later. - // Using squared rejection limit because it's cheaper than sqrts - float rej2 = PS_SQR(rej); // Rejection level squared - double sumWeights = 0.0; - for (int i = 0; i < num; i++) { - sumWeights += pixelWeights->data.F32[i]; - } - for (int i = 0; i < num; i++) { - // Systematic error contributes to the rejection level - float sysVar = PS_SQR(sys * pixelData->data.F32[i]); -#ifdef TESTING - // Correct variance for comparison against weighted mean including itself - float compare = 1.0 - pixelWeights->data.F32[i] / sumWeights; - if (x == TEST_X && y == TEST_Y) { - fprintf(stderr, "Variance %d (%d): %f %f %f\n", i, pixelSources->data.U16[i], - pixelVariances->data.F32[i], sysVar, compare); - } -#endif - - pixelLimits->data.F32[i] = rej2 * (pixelVariances->data.F32[i] + sysVar); - } - } - } - - // The clipping that follows is solely to identify suspect pixels. - // These suspect pixels will be inspected in more detail by other functions. - int numClipped = INT_MAX; // Number of pixels clipped per iteration - int totalClipped = 0; // Total number of pixels clipped - for (int i = 0; i < numIter && numClipped > 0 && num - totalClipped > 2; i++) { - numClipped = 0; - float median = NAN; // Middle of distribution - float limit = NAN; // Rejection limit - if (!useVariance) { - float stdev; // Median and stdev of the combination, for rejection - if (!combinationMedianStdev(&median, useVariance ? NULL : &stdev, - pixelData, pixelMasks, sort)) { - psWarning("Bad median/stdev at %d,%d", x, y); - break; - } - limit = rej * stdev; -#ifdef TESTING - if (x == TEST_X && y == TEST_Y) { - fprintf(stderr, "Rejecting without variance; rejection limit: %f\n", limit); - } -#endif - } else { -#ifdef TESTING - if (x == TEST_X && y == TEST_Y) { - fprintf(stderr, "Rejecting with variance...\n"); - } -#endif - median = combinationWeightedOlympic(pixelData, pixelWeights, pixelMasks, discard, sort); - } - -#ifdef TESTING - if (x == TEST_X && y == TEST_Y) { - fprintf(stderr, "Median: %f\n", median); - } -#endif - - -// Mask a pixel for inspection -#define MASK_PIXEL_FOR_INSPECTION() \ - pixelMasks->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xff; \ - if (!rejectInspect) { \ - combineInspect(inputs, x, y, pixelSources->data.U16[j]); \ - } \ - numClipped++; \ - totalClipped++; - - for (int j = 0; j < num; j++) { - if (pixelMasks->data.PS_TYPE_VECTOR_MASK_DATA[j]) { - continue; - } - float diff = pixelData->data.F32[j] - median; // Difference from expected - if (useVariance) { - // Comparing squares --- cheaper than lots of sqrts - // pixelVariances includes the rejection limit, from above - if (PS_SQR(diff) > pixelLimits->data.F32[j]) { - MASK_PIXEL_FOR_INSPECTION(); -#ifdef TESTING - if (x == TEST_X && y == TEST_Y) { - fprintf(stderr, "Rejecting input %d based on variance: %f > %f\n", - j, diff, sqrtf(pixelLimits->data.F32[j])); - } -#endif - } - } else if (fabsf(diff) > limit) { - MASK_PIXEL_FOR_INSPECTION(); -#ifdef TESTING - if (x == TEST_X && y == TEST_Y) { - fprintf(stderr, "Rejecting input %d based on distribution: %f > %f\n", - j, diff, limit); - } -#endif - } - } - } - - if (rejectInspect && totalClipped > 0) { - // Get rid of the masked values - // The alternative to this is to make combinationMeanVariance() accept a mask - int good = 0; // Index of good value - for (int i = 0; i < num; i++) { - if (pixelMasks->data.PS_TYPE_VECTOR_MASK_DATA[i]) { - continue; - } - if (i != good) { - pixelData->data.F32[good] = pixelData->data.F32[i]; - pixelVariances->data.F32[good] = pixelVariances->data.F32[i]; - pixelWeights->data.F32[good] = pixelWeights->data.F32[i]; - pixelData->data.F32[good] = pixelData->data.F32[i]; - } - good++; - } - pixelData->n = good; - pixelVariances->n = good; - pixelWeights->n = good; - if (combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) { - imageValue = mean; - varianceValue = var; - maskValue = 0; - } else { - imageValue = NAN; - varianceValue = NAN; - maskValue = bad; - } - } - break; } @@ -633,5 +454,341 @@ } +#ifdef TESTING + mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = num; +#endif + + return; +} + + +// Test pixels to be combined +// Returns false to repeat without suspect pixels +static bool combineTest(int num, // Number of good pixels + bool suspect, // Does the stack contain suspect pixels? + psArray *inputs, // Original inputs (for flagging) + combineBuffer *buffer, // Buffer with vectors + int x, int y, // Coordinates of interest; frame of output image + float iter, // Number of rejection iterations per input + float rej, // Number of standard deviations at which to reject + float sys, // Relative systematic error + float olympic,// Fraction of values to discard (Olympic weighted mean) + bool useVariance, // Use variance for rejection when combining? + bool safe // Combine safely? + ) +{ + if (iter <= 0) { + return true; + } + + int numIter = PS_MAX(iter * num, 1); // Number of iterations + +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Testing pixel %d,%d: %d %f %f %f %d %d\n", + x, y, numIter, rej, sys, olympic, useVariance, safe); + } +#endif + + psVector *pixelData = buffer->pixels; // Values for the pixel of interest + psVector *pixelWeights = buffer->weights; // Is the pixel suspect? + psVector *pixelVariances = buffer->variances; // Variances for the pixel of interest + psVector *pixelSources = buffer->sources; // Sources for the pixel of interest + psVector *pixelSuspects = buffer->suspects; // Is the pixel suspect? + psVector *pixelLimits = buffer->limits; // Is the pixel suspect? + + // Set up rejection limits + float rej2 = PS_SQR(rej); // Rejection level squared + if (num > 2 && useVariance) { + // Convert rejection limits --- saves doing it later multiple times + // Using squared rejection limit because it's cheaper than sqrts + double sumWeights = 0.0; + for (int i = 0; i < num; i++) { + sumWeights += pixelWeights->data.F32[i]; + } + for (int i = 0; i < num; i++) { + // Systematic error contributes to the rejection level + float sysVar = PS_SQR(sys * pixelData->data.F32[i]); +#ifdef TESTING + // Correct variance for comparison against weighted mean including itself + float compare = 1.0 - pixelWeights->data.F32[i] / sumWeights; + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Variance %d (%d): %f %f %f\n", i, pixelSources->data.U16[i], + pixelVariances->data.F32[i], sysVar, compare); + } +#endif + pixelLimits->data.F32[i] = rej2 * (pixelVariances->data.F32[i] + sysVar); + } + } + + int maskIndex = 0; // Index of pixel to mask + int totalClipped = 0; // Total number of pixels clipped + for (int i = 0; i < numIter && maskIndex >= 0; i++) { + maskIndex = -1; // Nothing to reject + + switch (num) { + case 0: + break; + case 1: + if (i == 0 && safe) { + combineMarkReject(inputs, x, y, pixelSources->data.U16[0]); + } + break; + case 2: { + if (useVariance) { + // Use variance to check that the two are consistent + float diff = 0.5 * (pixelData->data.F32[0] - pixelData->data.F32[1]); // Mean flux difference + float var1 = pixelVariances->data.F32[0]; // Variance of first + float var2 = pixelVariances->data.F32[1]; // Variance of second + // Systematic error contributes to the rejection level + var1 += PS_SQR(sys * pixelData->data.F32[0]); + var2 += PS_SQR(sys * pixelData->data.F32[1]); + + float sigma2 = var1 + var2; // Combined variance + if (PS_SQR(diff) > rej2 * sigma2) { + // Not consistent: don't believe either! + if (i == 0 && suspect) { + combineMarkReject(inputs, x, y, pixelSources->data.U16[0]); + combineMarkReject(inputs, x, y, pixelSources->data.U16[1]); + } else { + combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]); + combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]); + } +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Flagged both pixels (%f > %f x %f\n)", + diff, rej, sqrtf(sigma2)); + } +#endif + } + } else if (i == 0 && safe) { + // Can't test them, and we want to be safe, so reject + combineMarkReject(inputs, x, y, pixelSources->data.U16[0]); + combineMarkReject(inputs, x, y, pixelSources->data.U16[1]); + } + break; + } +#if 0 + case 3: { + // Want to be a bit careful on the rejection than for a larger number of inputs + if (!useVariance) { + return combineTestGeneral(num, suspect, inputs, buffer, x, y, numIter, rej, sys, + olympic, useVariance, safe, allowSuspect); + } + + // Differences between pixel values + float diff01 = pixelData->data.F32[0] - pixelData->data.F32[1]; + float diff12 = pixelData->data.F32[1] - pixelData->data.F32[2]; + float diff20 = pixelData->data.F32[2] - pixelData->data.F32[0]; + // Variance for each pixel + float var0 = pixelVariances->data.F32[0] + PS_SQR(sys * pixelData->data.F32[0]); + float var1 = pixelVariances->data.F32[1] + PS_SQR(sys * pixelData->data.F32[1]); + float var2 = pixelVariances->data.F32[2] + PS_SQR(sys * pixelData->data.F32[2]); + // Errors in pixel differences + float err01 = var0 + var1; + float err12 = var1 + var2; + float err20 = var2 + var0; + +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Diff 0-1: %f %f\n", diff01, err01); + fprintf(stderr, "Diff 1-2: %f %f\n", diff12, err12); + fprintf(stderr, "Diff 2-0: %f %f\n", diff20, err20); + } +#endif + + int badPairs = 0; // Number of bad pairs + bool bad01 = false, bad12 = false, bad20 = false; // Pair is bad? + if (PS_SQR(diff01) > rej2 * err01) { + bad01 = true; + badPairs++; + } + if (PS_SQR(diff12) > rej2 * err12) { + bad12 = true; + badPairs++; + } + if (PS_SQR(diff20) > rej2 * err20) { + bad20 = true; + badPairs++; + } + + if (badPairs > 0 && allowSuspect && suspect) { + return false; + } + + switch (badPairs) { + case 0: + // Nothing to worry about! + break; + case 1: + // Can't tell which image is bad, so be sure to get it + if (bad01) { + combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]); + combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]); + break; + } + if (bad12) { + combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]); + combineMarkInspect(inputs, x, y, pixelSources->data.U16[2]); + break; + } + if (bad20) { + combineMarkInspect(inputs, x, y, pixelSources->data.U16[2]); + combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]); + break; + } + psAbort("Should never get here"); + case 2: + if (bad01 && bad12) { + // 2 and 0 are good + combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]); + break; + } + if (bad12 && bad20) { + // 0 and 1 are good + combineMarkInspect(inputs, x, y, pixelSources->data.U16[2]); + break; + } + if (bad20 && bad01) { + // 1 and 2 are good + combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]); + break; + } + psAbort("Should never get here"); + case 3: + // Everything's bad + combineMarkInspect(inputs, x, y, pixelSources->data.U16[0]); + combineMarkInspect(inputs, x, y, pixelSources->data.U16[1]); + combineMarkInspect(inputs, x, y, pixelSources->data.U16[2]); + break; + } + break; + } +#endif + default: { + if (useVariance) { + float median = combinationWeightedOlympic(pixelData, pixelWeights, + olympic, buffer->sort); // Median for stack +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Flag with variance, median = %f\n", median); + } +#endif + float worst = -INFINITY; // Largest deviation + for (int j = 0; j < num; j++) { + float diff = pixelData->data.F32[j] - median; // Difference from expected +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Testing input %d: %f\n", j, diff); + } +#endif + + // Comparing squares --- cheaper than lots of sqrts + // pixelVariances includes the rejection limit, from above + float diff2 = PS_SQR(diff); // Square difference + if (diff2 > pixelLimits->data.F32[j]) { + float dev = diff2 / pixelLimits->data.F32[j]; // Deviation + if (dev > worst) { + worst = dev; + maskIndex = j; + } + } + } + } else { + float median = NAN, stdev = NAN; // Median and stdev of the combination, for rejection + combinationMedianStdev(&median, &stdev, pixelData, buffer->sort); + float limit = rej * stdev; // Rejection limit +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Flag without variance; median = %f, stdev = %f, limit = %f\n", + median, stdev, limit); + } +#endif + float worst = -INFINITY; // Largest deviation + for (int j = 0; j < num; j++) { + float diff = fabsf(pixelData->data.F32[j] - median); // Difference from expected + + if (diff > limit) { + float dev = diff / limit; // Deviation + if (dev > worst) { + worst = dev; + maskIndex = j; + } + } + } + } + } + } + + // Do the actual rejection of the pixel + if (maskIndex >= 0) { + if (suspect) { +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Throwing out all suspect pixels\n"); + } +#endif + // Throw out all suspect pixels + int numGood = 0; // Number of good pixels + for (int j = 0; j < num; j++) { + if (pixelSuspects->data.U8[j]) { + combineMarkReject(inputs, x, y, pixelSources->data.U16[j]); + continue; + } + if (numGood == j) { + numGood++; + continue; + } + pixelData->data.F32[numGood] = pixelData->data.F32[j]; + pixelWeights->data.F32[numGood] = pixelWeights->data.F32[j]; + pixelSources->data.U16[numGood] = pixelSources->data.U16[j]; + pixelLimits->data.F32[numGood] = pixelLimits->data.F32[j]; + pixelVariances->data.F32[numGood] = pixelVariances->data.F32[j]; + numGood++; + } + pixelData->n = numGood; + pixelWeights->n = numGood; + pixelSources->n = numGood; + pixelLimits->n = numGood; + pixelVariances->n = numGood; + totalClipped += num - numGood; + num = numGood; + suspect = false; + } else { + // Throw out masked pixel +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Throwing out input %d\n", maskIndex); + } +#endif + combineMarkInspect(inputs, x, y, pixelSources->data.U16[maskIndex]); + int numGood = 0; // Number of good pixels + for (int j = 0; j < num; j++) { + if (j == maskIndex) { + continue; + } + if (numGood == j) { + numGood++; + continue; + } + pixelData->data.F32[numGood] = pixelData->data.F32[j]; + pixelWeights->data.F32[numGood] = pixelWeights->data.F32[j]; + pixelSources->data.U16[numGood] = pixelSources->data.U16[j]; + pixelLimits->data.F32[numGood] = pixelLimits->data.F32[j]; + pixelVariances->data.F32[numGood] = pixelVariances->data.F32[j]; + numGood++; + } + pixelData->n = numGood; + pixelWeights->n = numGood; + pixelSources->n = numGood; + pixelLimits->n = numGood; + pixelVariances->n = numGood; + totalClipped++; + num--; + } + } + } + + return true; } @@ -639,8 +796,8 @@ // Ensure the input array of pmStackData is valid, and get some details out of it static bool validateInputData(bool *haveVariances, // Do we have variance maps in the sky images? - bool *haveRejects, // Do we have lists of rejected pixels? int *num, // Number of inputs int *numCols, int *numRows, // Size of (sky) images - psArray *input // Input array of pmStackData to validate + const psArray *input, // Input array of pmStackData to validate + const pmReadout *output // Output readout ) { @@ -668,5 +825,5 @@ PS_ASSERT_IMAGE_TYPE(data->readout->variance, PS_TYPE_F32, false); } - *haveRejects = (data->reject != NULL); + bool haveRejects = (data->reject != NULL); // Do we have rejected pixels? // Make sure the rest correspond with the first @@ -681,5 +838,5 @@ return false; } - if ((*haveRejects && !data->reject) || (data->reject && !*haveRejects)) { + if ((haveRejects && !data->reject) || (data->reject && !haveRejects)) { psError(PS_ERR_UNEXPECTED_NULL, true, "The rejected pixels are specified in some but not all inputs."); @@ -697,4 +854,13 @@ PS_ASSERT_IMAGE_TYPE(data->readout->variance, PS_TYPE_F32, false); } + } + + PM_ASSERT_READOUT_NON_NULL(output, false); + if (output->image) { + PS_ASSERT_IMAGE_NON_NULL(output->image, false); + PS_ASSERT_IMAGE_TYPE(output->image, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_NON_NULL(output->mask, false); + PS_ASSERT_IMAGE_TYPE(output->mask, PS_TYPE_IMAGE_MASK, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(output->image, output->mask, false); } @@ -782,31 +948,21 @@ /// Stack input images -bool pmStackCombine(pmReadout *combined, psArray *input, psImageMaskType maskVal, psImageMaskType bad, - int kernelSize, int numIter, float rej, float sys, float discard, - bool entire, bool useVariance, bool safe, bool rejectInspect) -{ - PS_ASSERT_PTR_NON_NULL(combined, false); +bool pmStackCombine(pmReadout *combined, psArray *input, psImageMaskType maskVal, psImageMaskType maskSuspect, + psImageMaskType bad, int kernelSize, float iter, float rej, float sys, float olympic, + bool useVariance, bool safe, bool rejection) +{ bool haveVariances; // Do we have the variance maps? - bool haveRejects; // Do we have lists of rejected pixels? int num; // Number of inputs int numCols, numRows; // Size of (sky) images - if (!validateInputData(&haveVariances, &haveRejects, &num, &numCols, &numRows, input)) { + if (!validateInputData(&haveVariances, &num, &numCols, &numRows, input, combined)) { return false; } PS_ASSERT_INT_NONNEGATIVE(kernelSize, false); PS_ASSERT_INT_POSITIVE(bad, false); - PS_ASSERT_INT_NONNEGATIVE(numIter, false); if (isnan(rej)) { - PS_ASSERT_INT_EQUAL(numIter, 0, false); + PS_ASSERT_FLOAT_EQUAL(iter, 0, false); } else { + PS_ASSERT_FLOAT_LARGER_THAN(iter, 0, false); PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false); - } - if (haveRejects) { - // This is a subsequent combination, so expect that the image and mask already exist - PS_ASSERT_IMAGE_NON_NULL(combined->image, false); - PS_ASSERT_IMAGE_TYPE(combined->image, PS_TYPE_F32, false); - PS_ASSERT_IMAGE_NON_NULL(combined->mask, false); - PS_ASSERT_IMAGE_TYPE(combined->mask, PS_TYPE_IMAGE_MASK, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(combined->image, combined->mask, false); } if (useVariance && !haveVariances) { @@ -833,6 +989,12 @@ } #endif - if (!haveRejects && !data->inspect) { - data->inspect = psPixelsAllocEmpty(PIXEL_LIST_BUFFER); + if (!rejection) { + // Ensure pixels can be put on the appropriate list + if (!data->inspect) { + data->inspect = psPixelsAllocEmpty(PIXEL_LIST_BUFFER); + } + if (!data->reject) { + data->reject = psPixelsAllocEmpty(PIXEL_LIST_BUFFER); + } } } @@ -863,91 +1025,85 @@ combineBuffer *buffer = combineBufferAlloc(num); - if (haveRejects) { - psImage *combinedImage = combined->image; // Combined image - psImage *combinedMask = combined->mask; // Combined mask - psImage *combinedVariance = combined->variance; // Combined variance map - - psArray *pixelMap = pixelMapGenerate(input, minInputCols, maxInputCols, - minInputRows, maxInputRows); // Map of pixels to source - psPixels *pixels = NULL; // Total list of pixels, with no duplicates + // Pull the products out, allocate if necessary + psImage *combinedImage = combined->image; // Combined image + if (!combinedImage) { + combined->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); + combinedImage = combined->image; + } + psImage *combinedMask = combined->mask; // Combined mask + if (!combinedMask) { + combined->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); + combinedMask = combined->mask; + } + + psImage *combinedVariance = combined->variance; // Combined variance map + if (haveVariances && !combinedVariance) { + combined->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); + combinedVariance = combined->variance; + } + + // Set up rejection list + psArray *pixelMap = NULL; // Map of pixels to source + if (rejection) { + pixelMap = pixelMapGenerate(input, minInputCols, maxInputCols, minInputRows, maxInputRows); + } + + // Combine each pixel + for (int y = minInputRows; y < maxInputRows; y++) { + for (int x = minInputCols; x < maxInputCols; x++) { +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Combining pixel %d,%d: %x %x %f %f %f %f %d %d %d\n", + x, y, maskVal, bad, iter, rej, sys, olympic, useVariance, safe, rejection); + } +#endif + psVector *reject = NULL; // Images to reject for this pixel + if (rejection) { + reject = pixelMapQuery(pixelMap, minInputCols, minInputRows, x, y); +#ifdef TESTING + if (x == TEST_X && y == TEST_Y) { + fprintf(stderr, "Rejected inputs: "); + if (!reject) { + fprintf(stderr, "\n"); + } else { + for (int i = 0; i < reject->n; i++) { + fprintf(stderr, "%d ", reject->data.U16[i]); + } + fprintf(stderr, "\n"); + } + } +#endif + } + + int num; // Number of good pixels + bool suspect; // Suspect pixels in stack? + combineExtract(&num, &suspect, buffer, combinedImage, combinedMask, combinedVariance, + input, weights, addVariance, reject, x, y, maskVal, maskSuspect); + combinePixels(combinedImage, combinedMask, combinedVariance, num, buffer, x, y, bad, safe); + + if (iter > 0) { + combineTest(num, suspect, input, buffer, x, y, iter, rej, sys, olympic, + useVariance, safe); + } + } + } + + psFree(pixelMap); + psFree(weights); + psFree(buffer); + + + +#ifndef PS_NO_TRACE + if (!rejection && psTraceGetLevel("psModules.imcombine") >= 5) { for (int i = 0; i < num; i++) { - pmStackData *data = input->data[i]; // Stacking data; contains the list of pixels - if (!data) { + pmStackData *data = input->data[i]; // Stack data for this input + if (!data || !data->inspect) { continue; } - pixels = psPixelsConcatenate(pixels, data->reject); - } - pixels = psPixelsDuplicates(pixels, pixels); - - if (entire) { - // Combine entire image - for (int y = minInputRows; y < maxInputRows; y++) { - for (int x = minInputCols; x < maxInputCols; x++) { - psVector *reject = pixelMapQuery(pixelMap, minInputCols, minInputRows, - x, y); // Reject these images - combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, - addVariance, reject, x, y, maskVal, bad, numIter, rej, sys, discard, - useVariance, safe, rejectInspect, buffer); - } - } - } else { - // Only combine previously rejected pixels - for (int i = 0; i < pixels->n; i++) { - // Pixel coordinates are in the frame of the original image - int x = pixels->data[i].x, y = pixels->data[i].y; // Coordinates of interest - if (x < minInputCols || x >= maxInputCols || y < minInputRows || y >= maxInputRows) { - continue; - } - psVector *reject = pixelMapQuery(pixelMap, minInputCols, minInputRows, - x, y); // Reject these images - combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, - addVariance, reject, x, y, maskVal, bad, numIter, rej, sys, discard, - useVariance, safe, rejectInspect, buffer); - } - } - psFree(pixels); - psFree(pixelMap); - } else { - // Pull the products out, allocate if necessary - psImage *combinedImage = combined->image; // Combined image - if (!combinedImage) { - combined->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); - combinedImage = combined->image; - } - psImage *combinedMask = combined->mask; // Combined mask - if (!combinedMask) { - combined->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); - combinedMask = combined->mask; - } - - psImage *combinedVariance = combined->variance; // Combined variance map - if (haveVariances && !combinedVariance) { - combined->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); - combinedVariance = combined->variance; - } - - for (int y = minInputRows; y < maxInputRows; y++) { - for (int x = minInputCols; x < maxInputCols; x++) { - combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, - addVariance, NULL, x, y, maskVal, bad, numIter, rej, sys, discard, - useVariance, safe, rejectInspect, buffer); - } - } - -#ifndef PS_NO_TRACE - if (psTraceGetLevel("psModules.imcombine") >= 5) { - for (int i = 0; i < num; i++) { - pmStackData *data = input->data[i]; // Stack data for this input - if (!data || !data->inspect) { - continue; - } - psTrace("psModules.imcombine", 5, "Image %d: %ld pixels to inspect.\n", i, data->inspect->n); - } - } -#endif - } - - psFree(weights); - psFree(buffer); + psTrace("psModules.imcombine", 5, "Image %d: %ld pixels to inspect.\n", i, data->inspect->n); + } + } +#endif return true; Index: trunk/psModules/src/imcombine/pmStack.h =================================================================== --- trunk/psModules/src/imcombine/pmStack.h (revision 26075) +++ trunk/psModules/src/imcombine/pmStack.h (revision 26076) @@ -44,11 +44,11 @@ psArray *input, ///< Input array of pmStackData psImageMaskType maskVal, ///< Mask value of bad pixels + psImageMaskType suspect, ///< Mask value of suspect pixels psImageMaskType bad, ///< Mask value to give rejected pixels int kernelSize, ///< Half-size of the convolution kernel - int numIter, ///< Number of iterations + float iter, ///< Number of iterations per input float rej, ///< Rejection limit (standard deviations) float sys, ///< Relative systematic error float discard, ///< Fraction of values to discard for Olympic weighted mean - bool entire, ///< Combine entire image even if rejection lists provided? bool useVariance, ///< Use variance values for rejection? bool safe, ///< Play safe with small numbers of input pixels (mask if N <= 2)? Index: trunk/psModules/src/objects/pmSourceMatch.c =================================================================== --- trunk/psModules/src/objects/pmSourceMatch.c (revision 26075) +++ trunk/psModules/src/objects/pmSourceMatch.c (revision 26076) @@ -17,4 +17,5 @@ #define SOURCES_MAX_LEAF 2 // Maximum number of points on a tree leaf #define ARRAY_BUFFER 16 // Buffer for array + ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -124,4 +125,6 @@ match->mag = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_F32); match->magErr = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_F32); + match->x = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_F32); + match->y = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_F32); match->image = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_U32); match->index = psVectorAllocEmpty(ARRAY_BUFFER, PS_TYPE_U32); @@ -133,4 +136,5 @@ void pmSourceMatchAdd(pmSourceMatch *match, // Match data float mag, float magErr, // Magnitude and error + float x, float y, // Position int image, // Image index int index // Source index @@ -141,4 +145,6 @@ match->mag = psVectorExtend(match->mag, match->mag->nalloc, 1); match->magErr = psVectorExtend(match->magErr, match->magErr->nalloc, 1); + match->x = psVectorExtend(match->x, match->x->nalloc, 1); + match->y = psVectorExtend(match->y, match->y->nalloc, 1); match->image = psVectorExtend(match->image, match->image->nalloc, 1); match->index = psVectorExtend(match->index, match->index->nalloc, 1); @@ -147,4 +153,6 @@ match->mag->data.F32[num] = mag; match->magErr->data.F32[num] = magErr; + match->x->data.F32[num] = x; + match->y->data.F32[num] = y; match->image->data.S32[num] = image; match->index->data.S32[num] = index; @@ -192,5 +200,5 @@ pmSourceMatch *match = pmSourceMatchAlloc(); // Match data pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], - i, indices->data.S32[j]); + xImage->data.F32[j], yImage->data.F32[j], i, indices->data.S32[j]); matches->data[j] = match; } @@ -212,5 +220,5 @@ pmSourceMatch *match = pmSourceMatchAlloc(); // Match data pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], - i, indices->data.S32[j]); + xImage->data.F32[j], yImage->data.F32[j], i, indices->data.S32[j]); matches->data[k] = match; } @@ -238,5 +246,5 @@ pmSourceMatch *match = matches->data[index]; // Match data pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], - i, indices->data.S32[j]); + xImage->data.F32[j], yImage->data.F32[j], i, indices->data.S32[j]); numMatch++; } else { @@ -244,5 +252,5 @@ pmSourceMatch *match = pmSourceMatchAlloc(); // Match data pmSourceMatchAdd(match, magImage->data.F32[j], magErrImage->data.F32[j], - i, indices->data.S32[j]); + xImage->data.F32[j], yImage->data.F32[j], i, indices->data.S32[j]); xMaster->data.F32[numMaster] = xImage->data.F32[j]; yMaster->data.F32[numMaster] = yImage->data.F32[j]; @@ -304,5 +312,6 @@ pmSourceMatch *match = matches->data[i]; // Match of interest for (int j = 0; j < match->num && !source; j++) { - if (!isfinite(match->mag->data.F32[j]) || !isfinite(match->magErr->data.F32[j])) { + if (!isfinite(match->mag->data.F32[j]) || !isfinite(match->magErr->data.F32[j]) || + !isfinite(match->x->data.F32[j]) || !isfinite(match->y->data.F32[j])) { continue; } @@ -365,5 +374,5 @@ double star = 0.0, starErr = 0.0; // Accumulators for star for (int j = 0; j < match->num; j++) { - if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) { + if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_PHOT) { continue; } @@ -396,5 +405,5 @@ pmSourceMatch *match = matches->data[i]; // Matched stars for (int j = 0; j < match->num; j++) { - if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) { + if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_PHOT) { continue; } @@ -424,5 +433,5 @@ pmSourceMatch *match = matches->data[i]; // Matched stars for (int j = 0; j < match->num; j++) { - if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) { + if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_PHOT) { continue; } @@ -543,5 +552,5 @@ pmSourceMatch *match = matches->data[i]; // Matched stars for (int j = 0; j < match->num; j++) { - if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) { + if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_PHOT) { continue; } @@ -564,5 +573,5 @@ if (PS_SQR(dev) > starClip * (PS_SQR(magErr) + sysErr2)) { numRejected++; - match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xFF; + match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] |= PM_SOURCE_MATCH_MASK_PHOT; } } @@ -627,8 +636,8 @@ return NULL; } - psTrace("psModules.objects", 3, "Pass 1: Determined %d/%d are photometric", numPhoto, numImages); + psTrace("psModules.objects", 3, "Pass 1: Determined %d/%d are photometric\n", numPhoto, numImages); fracRej = sourceMatchRelphotReject(trans, stars, matches, zp, photo, badImage, rej1, sys1); - psTrace("psModules.objects", 3, "Pass 1: %f%% of measurements rejected", fracRej * 100); + psTrace("psModules.objects", 3, "Pass 1: %f%% of measurements rejected\n", fracRej * 100); chi2 = sourceMatchRelphotIterate(trans, stars, badImage, matches, zp, photo, sys1); @@ -649,8 +658,8 @@ return NULL; } - psTrace("psModules.objects", 3, "Pass 2: Determined %d/%d are photometric", numPhoto, numImages); + psTrace("psModules.objects", 3, "Pass 2: Determined %d/%d are photometric\n", numPhoto, numImages); fracRej = sourceMatchRelphotReject(trans, stars, matches, zp, photo, badImage, rej2, sys2); - psTrace("psModules.objects", 3, "Pass 2: %f%% of measurements rejected", fracRej * 100); + psTrace("psModules.objects", 3, "Pass 2: %f%% of measurements rejected\n", fracRej * 100); chi2 = sourceMatchRelphotIterate(trans, stars, badImage, matches, zp, photo, sys2); @@ -668,2 +677,220 @@ return trans; } + + +// Iterate on the star positions and image shifts +// Returns the solution chi^2 +static float sourceMatchRelastroIterate(psVector *xShift, psVector *yShift, // Shift for image + psVector *xStar, psVector *yStar, // Position for star + const psArray *matches // Array of matches + ) +{ + psAssert(matches, "Need list of matches"); + + int numImages = xShift->n; // Number of images + int numStars = matches->n; // Number of stars + + psAssert(xShift && xShift->type.type == PS_TYPE_F32 && yShift && yShift->type.type == PS_TYPE_F32, + "Need shifts"); + psAssert(yShift->n == numImages, "Not enough shifts: %ld\n", yShift->n); + psAssert(xStar && xStar->type.type == PS_TYPE_F32 && yStar && yStar->type.type == PS_TYPE_F32, + "Need star positions"); + psAssert(xStar->n == numStars && yStar->n == numStars, "Not enough stars\n"); + + // Solve the star positions + psVectorInit(xStar, NAN); + psVectorInit(yStar, NAN); + psVector *starMask = psVectorAlloc(numStars, PS_TYPE_U8); // Mask for stars + psVectorInit(starMask, 0xFF); + int numGoodStars = 0; // Number of stars with good measurements + for (int i = 0; i < numStars; i++) { + pmSourceMatch *match = matches->data[i]; // Matched stars + int numMeasurements = 0; // Number of unmasked measurements for star + double xSum = 0.0, ySum = 0.0; // Accumulators for star + for (int j = 0; j < match->num; j++) { + if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_ASTRO) { + continue; + } + numMeasurements++; + int index = match->image->data.U32[j]; // Image index + + xSum += match->x->data.F32[j] - xShift->data.F32[index]; + ySum += match->y->data.F32[j] - yShift->data.F32[index]; + } + if (numMeasurements > 1) { + // It's only a good star (contributing to the chi^2) if there's more than 1 measurement + numGoodStars++; + xStar->data.F32[i] = xSum / numMeasurements; + yStar->data.F32[i] = ySum / numMeasurements; + starMask->data.U8[i] = 0; + } + } + + // Solve for the shifts + psVectorInit(xShift, 0.0); + psVectorInit(yShift, 0.0); + psVector *num = psVectorAlloc(numImages, PS_TYPE_S32); // Number of stars + psVectorInit(num, 0); + for (int i = 0; i < numStars; i++) { + if (starMask->data.U8[i]) { + continue; + } + pmSourceMatch *match = matches->data[i]; // Matched stars + for (int j = 0; j < match->num; j++) { + if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_ASTRO) { + continue; + } + int index = match->image->data.U32[j]; // Image index + + xShift->data.F32[index] += match->x->data.F32[j] - xStar->data.F32[i]; + yShift->data.F32[index] += match->y->data.F32[j] - yStar->data.F32[i]; + num->data.S32[index]++; + } + } + for (int i = 0; i < numImages; i++) { + xShift->data.F32[i] /= num->data.S32[i]; + yShift->data.F32[i] /= num->data.S32[i]; + psTrace("psModules.objects", 3, "Shift for image %d: %f,%f\n", + i, xShift->data.F32[i], yShift->data.F32[i]); + } + psFree(num); + + // Once more through to evaluate chi^2 + float chi2 = 0.0; // chi^2 for iteration + int dof = 0; // Degrees of freedom + for (int i = 0; i < numStars; i++) { + pmSourceMatch *match = matches->data[i]; // Matched stars + if (starMask->data.U8[i]) { + continue; + } + for (int j = 0; j < match->num; j++) { + if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j]) { + continue; + } + + int index = match->image->data.U32[j]; // Image index + float dx = match->x->data.F32[j] - xShift->data.F32[index] - xStar->data.F32[i]; + float dy = match->y->data.F32[j] - yShift->data.F32[index] - yStar->data.F32[i]; + + chi2 += PS_SQR(dx) + PS_SQR(dy); + dof++; + } + } + dof -= numGoodStars + numImages; + chi2 /= dof; + + return chi2; +} + +// Reject star measurements +// Returns the fraction of measurements that were rejected +static float sourceMatchRelastroReject(const psVector *xShift, const psVector *yShift, // Shifts for each image + const psVector *xStar, const psVector *yStar, // Positions for each star + const psArray *matches, // Array of matches + float chi2, // chi^2 from fit + float rej // Rejection threshold + ) +{ + psAssert(matches, "Need list of matches"); + + int numImages = xShift->n; // Number of images + int numStars = matches->n; // Number of stars + + psAssert(xShift && xShift->type.type == PS_TYPE_F32 && yShift && yShift->type.type == PS_TYPE_F32, + "Need shifts"); + psAssert(yShift->n == numImages, "Not enough shifts: %ld\n", yShift->n); + psAssert(xStar && xStar->type.type == PS_TYPE_F32 && yStar && yStar->type.type == PS_TYPE_F32, + "Need star positions"); + psAssert(xStar->n == numStars && yStar->n == numStars, "Not enough stars\n"); + + int numRejected = 0; // Number rejected + int numMeasurements = 0; // Number of measurements + + float thresh = PS_SQR(rej) * chi2; // Threshold for rejection + + for (int i = 0; i < numStars; i++) { + pmSourceMatch *match = matches->data[i]; // Matched stars + for (int j = 0; j < match->num; j++) { + if (match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_SOURCE_MATCH_MASK_ASTRO) { + continue; + } + numMeasurements++; + int index = match->image->data.U32[j]; // Image index + + float dx = match->x->data.F32[j] - xShift->data.F32[index] - xStar->data.F32[i]; + float dy = match->y->data.F32[j] - yShift->data.F32[index] - yStar->data.F32[i]; + + if (PS_SQR(dx) + PS_SQR(dy) > thresh) { + numRejected++; + match->mask->data.PS_TYPE_VECTOR_MASK_DATA[j] |= PM_SOURCE_MATCH_MASK_ASTRO; + } + } + } + + return (float)numRejected / (float)numMeasurements; +} + +psArray *pmSourceMatchRelastro(const psArray *matches, // Array of matches + int numImages, // Number of images + float tol, // Relative tolerance for convergence + int iter1, // Number of iterations for pass 1 + float rej1, // Limit on rejection between iterations for pass 1 + int iter2, // Number of iterations for pass 2 + float rej2, // Limit on rejection between iterations for pass 2 + float rejLimit // Limit on rejection between iterations + ) +{ + PS_ASSERT_ARRAY_NON_NULL(matches, NULL); + + int numStars = matches->n; // Number of stars + psVector *xShift = psVectorAlloc(numImages, PS_TYPE_F32); // x shift for each image + psVector *yShift = psVectorAlloc(numImages, PS_TYPE_F32); // y shift for each image + psVectorInit(xShift, 0.0); + psVectorInit(yShift, 0.0); + psVector *xStar = psVectorAlloc(numStars, PS_TYPE_F32); // x position for each star + psVector *yStar = psVectorAlloc(numStars, PS_TYPE_F32); // y position for each star + + float chi2 = sourceMatchRelastroIterate(xShift, yShift, xStar, yStar, matches); // chi^2 for solution + psTrace("psModules.objects", 1, "Initial: chi^2 = %f\n", chi2); + float lastChi2 = INFINITY; // chi^2 on last iteration + float fracRej = INFINITY; // Fraction of measurements rejected + + // In the first passes, the shifts are not well deteremined: use high systematic error and + // rejection thresholds + for (int i = 0; i < iter1; i++) { + fracRej = sourceMatchRelastroReject(xShift, yShift, xStar, yStar, matches, chi2, rej1); + psTrace("psModules.objects", 3, "Pass 1: %f%% of measurements rejected\n", fracRej * 100); + + chi2 = sourceMatchRelastroIterate(xShift, yShift, xStar, yStar, matches); + psTrace("psModules.objects", 1, "Pass 1: iter = %d: chi^2 = %f rejected = %f\n", i, chi2, fracRej); + } + + for (int i = 0; i < iter2 && (fabsf(lastChi2 - chi2) > tol * chi2 || fracRej > rejLimit); i++) { + lastChi2 = chi2; + + fracRej = sourceMatchRelastroReject(xShift, yShift, xStar, yStar, matches, chi2, rej2); + psTrace("psModules.objects", 3, "Pass 2: %f%% of measurements rejected\n", fracRej * 100); + + chi2 = sourceMatchRelastroIterate(xShift, yShift, xStar, yStar, matches); + psTrace("psModules.objects", 1, "Pass 2: iter = %d: chi^2 = %f rejected = %f\n", i, chi2, fracRej); + } + + psFree(xStar); + psFree(yStar); + + if (fabsf(lastChi2 - chi2) > tol * chi2 || fracRej > rejLimit) { + psWarning("Unable to converge to relphot solution (%f,%f)", (lastChi2 - chi2) / chi2, fracRej); + } + + psArray *results = psArrayAlloc(numImages); // Array of results + for (int i = 0; i < numImages; i++) { + psVector *offset = results->data[i] = psVectorAlloc(2, PS_TYPE_F32); // Offset for image + offset->data.F32[0] = xShift->data.F32[i]; + offset->data.F32[1] = yShift->data.F32[i]; + } + psFree(xShift); + psFree(yShift); + + return results; +} + Index: trunk/psModules/src/objects/pmSourceMatch.h =================================================================== --- trunk/psModules/src/objects/pmSourceMatch.h (revision 26075) +++ trunk/psModules/src/objects/pmSourceMatch.h (revision 26076) @@ -1,4 +1,12 @@ #ifndef PM_SOURCE_MATCH_H #define PM_SOURCE_MATCH_H + +#include + +/// Mask values for matched sources +typedef enum { + PM_SOURCE_MATCH_MASK_PHOT = 0x01, // Source was rejected from photometry fit + PM_SOURCE_MATCH_MASK_ASTRO = 0x02, // Source was rejected from astrometry fit +} pmSourceMatchMask; /// Matched sources @@ -12,4 +20,5 @@ psVector *mag; // Magnitudes psVector *magErr; // Magnitude errors + psVector *x, *y; // Positions psVector *mask; // Mask for measurements } pmSourceMatch; @@ -26,8 +35,12 @@ PS_ASSERT_VECTOR_NON_NULL((MATCH)->mag, RVAL); \ PS_ASSERT_VECTOR_NON_NULL((MATCH)->magErr, RVAL); \ + PS_ASSERT_VECTOR_NON_NULL((MATCH)->x, RVAL); \ + PS_ASSERT_VECTOR_NON_NULL((MATCH)->y, RVAL); \ PS_ASSERT_VECTOR_SIZE((MATCH)->image, (MATCH)->num, RVAL); \ PS_ASSERT_VECTOR_SIZE((MATCH)->index, (MATCH)->num, RVAL); \ PS_ASSERT_VECTOR_SIZE((MATCH)->mag, (MATCH)->num, RVAL); \ PS_ASSERT_VECTOR_SIZE((MATCH)->magErr, (MATCH)->num, RVAL); \ + PS_ASSERT_VECTOR_SIZE((MATCH)->x, (MATCH)->num, RVAL); \ + PS_ASSERT_VECTOR_SIZE((MATCH)->y, (MATCH)->num, RVAL); \ } @@ -38,4 +51,5 @@ void pmSourceMatchAdd(pmSourceMatch *match, // Match data float mag, float magErr, // Magnitude and error + float x, float y, // Position int image, // Image index int index // Source index @@ -75,3 +89,14 @@ ); +/// Perform relative astrometry to calibrate images +psArray *pmSourceMatchRelastro(const psArray *matches, // Array of matches + int numImages, // Number of images + float tol, // Relative tolerance for convergence + int iter1, // Number of iterations for pass 1 + float rej1, // Limit on rejection between iterations for pass 1 + int iter2, // Number of iterations for pass 2 + float rej2, // Limit on rejection between iterations for pass 2 + float rejLimit // Limit on rejection between iterations + ); + #endif