Index: branches/pap/ppStack/src/ppStackArguments.c =================================================================== --- branches/pap/ppStack/src/ppStackArguments.c (revision 26006) +++ branches/pap/ppStack/src/ppStackArguments.c (revision 26007) @@ -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: branches/pap/ppStack/src/ppStackCombineFinal.c =================================================================== --- branches/pap/ppStack/src/ppStackCombineFinal.c (revision 26006) +++ branches/pap/ppStack/src/ppStackCombineFinal.c (revision 26007) @@ -10,5 +10,5 @@ #include "ppStackLoop.h" -#define TESTING // Enable test output +//#define TESTING // Enable test output bool ppStackCombineFinal(pmReadout *target, ppStackThreadData *stack, psArray *covariances, Index: branches/pap/ppStack/src/ppStackCombineInitial.c =================================================================== --- branches/pap/ppStack/src/ppStackCombineInitial.c (revision 26006) +++ branches/pap/ppStack/src/ppStackCombineInitial.c (revision 26007) @@ -10,5 +10,5 @@ #include "ppStackLoop.h" -#define TESTING // Enable test output +//#define TESTING // Enable test output bool ppStackCombineInitial(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config) Index: branches/pap/ppStack/src/ppStackConvolve.c =================================================================== --- branches/pap/ppStack/src/ppStackConvolve.c (revision 26006) +++ branches/pap/ppStack/src/ppStackConvolve.c (revision 26007) @@ -9,4 +9,7 @@ #include "ppStack.h" #include "ppStackLoop.h" + +//#define TESTING + // Update the value of a concept @@ -130,5 +133,4 @@ ppStackWriteImage(options->convMasks->data[i], maskHeader, readout->mask, config); psFree(maskHeader); - psImageCovarianceTransfer(readout->variance, readout->covariance); ppStackWriteImage(options->convVariances->data[i], hdu->header, readout->variance, config); #ifdef TESTING @@ -139,4 +141,19 @@ pmStackVisualPlotTestImage(readout->covariance->image, name); 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 @@ -221,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; } @@ -245,4 +262,8 @@ // 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); Index: branches/pap/ppStack/src/ppStackLoop.c =================================================================== --- branches/pap/ppStack/src/ppStackLoop.c (revision 26006) +++ branches/pap/ppStack/src/ppStackLoop.c (revision 26007) @@ -130,5 +130,5 @@ } psTrace("ppStack", 2, "Stack of unconvolved images....\n"); - if (!ppStackCombineFinal(options->unconvRO, stack, options->origCovars, options, config, true, true)) { + if (!ppStackCombineFinal(options->unconvRO, stack, options->origCovars, options, config, false, true)) { psError(PS_ERR_UNKNOWN, false, "Unable to perform unconvolved combination."); psFree(stack); Index: branches/pap/ppStack/src/ppStackMatch.c =================================================================== --- branches/pap/ppStack/src/ppStackMatch.c (revision 26006) +++ branches/pap/ppStack/src/ppStackMatch.c (revision 26007) @@ -18,5 +18,5 @@ #define COVAR_FRAC 0.01 // Truncation fraction for covariance matrix -#define TESTING // Enable debugging output +//#define TESTING // Enable debugging output #ifdef TESTING @@ -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 @@ -377,4 +382,7 @@ } #endif + + fprintf(stderr, "vf = %f\n", psImageCovarianceFactor(readout->covariance)); + if (threads > 0) { @@ -516,4 +524,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)); } @@ -542,13 +551,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)); + } } Index: branches/pap/ppStack/src/ppStackReadout.c =================================================================== --- branches/pap/ppStack/src/ppStackReadout.c (revision 26006) +++ branches/pap/ppStack/src/ppStackReadout.c (revision 26007) @@ -116,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 @@ -126,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 @@ -161,5 +163,5 @@ } - if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, kernelSize, iter, + 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."); @@ -217,6 +219,8 @@ 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 @@ -267,5 +271,5 @@ } - if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, 0, iter, combineRej, + 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."); Index: branches/pap/ppStack/src/ppStackReject.c =================================================================== --- branches/pap/ppStack/src/ppStackReject.c (revision 26006) +++ branches/pap/ppStack/src/ppStackReject.c (revision 26007) @@ -10,5 +10,5 @@ #include "ppStackLoop.h" -#define TESTING +//#define TESTING bool ppStackReject(ppStackOptions *options, pmConfig *config)