Index: /trunk/ppStack/src/ppStack.h =================================================================== --- /trunk/ppStack/src/ppStack.h (revision 21476) +++ /trunk/ppStack/src/ppStack.h (revision 21477) @@ -46,4 +46,5 @@ const psArray *maskNames, // Names of masks to read const psArray *varianceNames, // Names of variance maps to read + const psArray *covariances, // Covariance matrices (already read) const pmConfig *config // Configuration ); @@ -97,6 +98,5 @@ pmReadout *outRO, // Output readout const psArray *readouts, // Input readouts - const psArray *regions, // Array with array of regions used in each PSF match - const psArray *kernels, // Array with array of kernels used in each PSF match + const psVector *mask, // Mask for input readouts const psVector *weightings, // Weighting factors for each image const psVector *addVariance // Additional variance for rejection @@ -115,6 +115,8 @@ pmReadout *outRO, // Output readout const psArray *readouts, // Input readouts + const psVector *mask, // Mask for input readouts const psArray *rejected, // Array with pixels rejected in each image - const psVector *weightings // Weighting factors for each image + const psVector *weightings, // Weighting factors for each image + const psVector *addVariance // Additional variance for rejection ); Index: /trunk/ppStack/src/ppStackArguments.c =================================================================== --- /trunk/ppStack/src/ppStackArguments.c (revision 21476) +++ /trunk/ppStack/src/ppStackArguments.c (revision 21477) @@ -147,4 +147,6 @@ psMetadataAddF32(arguments, PS_LIST_TAIL, "-combine-sys", 0, "Relative systematic error in combination", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-combine-discard", 0, + "Discard fraction for Olympic weighted mean", NAN); psMetadataAddStr(arguments, PS_LIST_TAIL, "-mask-val", 0, "Mask value of input bad pixels", NULL); psMetadataAddStr(arguments, PS_LIST_TAIL, "-mask-bad", 0, "Mask value to give bad pixels", NULL); @@ -232,12 +234,13 @@ } - VALUE_ARG_RECIPE_INT("-iter", "ITER", S32, 0); - VALUE_ARG_RECIPE_FLOAT("-combine-rej", "COMBINE.REJ", F32); - VALUE_ARG_RECIPE_FLOAT("-combine-sys", "COMBINE.SYS", F32); - VALUE_ARG_RECIPE_FLOAT("-threshold-mask", "THRESHOLD.MASK", F32); - VALUE_ARG_RECIPE_FLOAT("-image-rej", "IMAGE.REJ", F32); - VALUE_ARG_RECIPE_FLOAT("-deconv-limit", "DECONV.LIMIT", F32); - VALUE_ARG_RECIPE_INT("-rows", "ROWS", S32, 0); - VALUE_ARG_RECIPE_FLOAT("-poor-frac", "POOR.FRACTION", F32); + VALUE_ARG_RECIPE_INT("-iter", "ITER", S32, 0); + VALUE_ARG_RECIPE_FLOAT("-combine-rej", "COMBINE.REJ", F32); + VALUE_ARG_RECIPE_FLOAT("-combine-sys", "COMBINE.SYS", F32); + VALUE_ARG_RECIPE_FLOAT("-combine-discard", "COMBINE.DISCARD", F32); + VALUE_ARG_RECIPE_FLOAT("-threshold-mask", "THRESHOLD.MASK", F32); + VALUE_ARG_RECIPE_FLOAT("-image-rej", "IMAGE.REJ", F32); + VALUE_ARG_RECIPE_FLOAT("-deconv-limit", "DECONV.LIMIT", F32); + VALUE_ARG_RECIPE_INT("-rows", "ROWS", S32, 0); + VALUE_ARG_RECIPE_FLOAT("-poor-frac", "POOR.FRACTION", F32); valueArgRecipeStr(arguments, recipe, "-mask-val", "MASK.VAL", recipe); Index: /trunk/ppStack/src/ppStackLoop.c =================================================================== --- /trunk/ppStack/src/ppStackLoop.c (revision 21476) +++ /trunk/ppStack/src/ppStackLoop.c (revision 21477) @@ -480,5 +480,6 @@ psMetadataAddF32(stats, PS_LIST_TAIL, "STAMP.NUM", PS_META_DUPLICATE_OK, "Number of stamps", kernels->numStamps); - + psMetadataAddF32(stats, PS_LIST_TAIL, "PPSTACK.WEIGHTING", PS_META_DUPLICATE_OK, + "Weighting for image", weightings->data.F32[i]); } psLogMsg("ppStack", PS_LOG_INFO, "Time to match image %d: %f sec", i, psTimerClear("PPSTACK_MATCH")); @@ -495,5 +496,14 @@ writeImage(maskNames->data[i], maskHeader, readout->mask, config); psFree(maskHeader); + psImageCovarianceTransfer(readout->variance, readout->covariance); writeImage(varianceNames->data[i], hdu->header, readout->variance, config); +#ifdef TESTING + { + psString name = NULL; + psStringAppend(&name, "covariance_%d.fits", i); + writeImage(name, hdu->header, readout->covariance->image, config); + psFree(name); + } +#endif pmCell *inCell = readout->parent; // Input cell @@ -605,4 +615,6 @@ i, matchChi2->data.F32[i]); } else { + psLogMsg("ppStack", PS_LOG_INFO, "Image %d has matching chi^2: %f", + i, matchChi2->data.F32[i]); numGood++; } @@ -630,5 +642,6 @@ // Start threading ppStackThreadInit(); - ppStackThreadData *stack = ppStackThreadDataSetup(cells, imageNames, maskNames, varianceNames, config); + ppStackThreadData *stack = ppStackThreadDataSetup(cells, imageNames, maskNames, varianceNames, + covariances, config); if (!stack) { psError(PS_ERR_IO, false, "Unable to initialise stack threads."); @@ -729,6 +742,5 @@ psArrayAdd(job->args, 1, config); psArrayAdd(job->args, 1, outRO); - psArrayAdd(job->args, 1, subRegions); - psArrayAdd(job->args, 1, subKernels); + psArrayAdd(job->args, 1, inputMask); psArrayAdd(job->args, 1, weightings); psArrayAdd(job->args, 1, matchChi2); @@ -761,5 +773,4 @@ return false; } - psFree(matchChi2); // Harvest the jobs, gathering the inspection lists @@ -787,4 +798,8 @@ memDump("initial"); } + +#ifdef TESTING + writeImage("combined_initial.fits", NULL, outRO->image, config); +#endif if (stats) { @@ -828,4 +843,5 @@ psFree(rejected); psFree(covariances); + psFree(matchChi2); return false; } @@ -844,4 +860,5 @@ psFree(rejected); psFree(covariances); + psFree(matchChi2); return false; } @@ -939,4 +956,5 @@ psFree(outRO); psFree(covariances); + psFree(matchChi2); return false; } @@ -969,4 +987,5 @@ psFree(outRO); psFree(covariances); + psFree(matchChi2); return false; } @@ -981,6 +1000,8 @@ psArrayAdd(job->args, 1, config); psArrayAdd(job->args, 1, outRO); + psArrayAdd(job->args, 1, inputMask); psArrayAdd(job->args, 1, rejected); psArrayAdd(job->args, 1, weightings); + psArrayAdd(job->args, 1, matchChi2); if (!psThreadJobAddPending(job)) { psFree(job); @@ -991,4 +1012,5 @@ psFree(outRO); psFree(covariances); + psFree(matchChi2); return false; } @@ -1004,4 +1026,5 @@ psFree(outRO); psFree(covariances); + psFree(matchChi2); return false; } @@ -1009,4 +1032,8 @@ memDump("final"); + +#ifdef TESTING + writeImage("combined_final.fits", NULL, outRO->image, config); +#endif // Sum covariance matrices @@ -1019,4 +1046,5 @@ outRO->covariance = psImageCovarianceSum(covariances); psFree(covariances); + psFree(matchChi2); if (stats) { Index: /trunk/ppStack/src/ppStackMatch.c =================================================================== --- /trunk/ppStack/src/ppStackMatch.c (revision 21476) +++ /trunk/ppStack/src/ppStackMatch.c (revision 21477) @@ -170,4 +170,5 @@ assert(kernels && !*kernels); assert(config); + *weighting = 0.0; psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe @@ -260,4 +261,11 @@ pmSubtractionAnalysis(readout->analysis, kernels, region, readout->image->numCols, readout->image->numRows); + + psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel + psKernel *covar = psImageCovarianceCalculate(kernel, readout->covariance); // New covariance matrix + psFree(readout->covariance); + readout->covariance = covar; + psFree(kernel); + } else { #endif @@ -541,5 +549,5 @@ } psFree(iter); - *chi2 /= num; + *chi2 /= psImageCovarianceFactor(readout->covariance) * num; } @@ -579,9 +587,17 @@ (void)psBinaryOp(readout->image, readout->image, "-", psScalarAlloc(psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), PS_TYPE_F32)); - *weighting = 1.0 / PS_SQR(psStatsGetValue(bg, PS_STAT_ROBUST_STDEV)); - psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "PPSTACK.WEIGHTING", 0, - "Weighting by 1/noise^2 for stack", - 1.0 / PS_SQR(psStatsGetValue(bg, PS_STAT_ROBUST_STDEV))); - } + } + + // Measure the variance level for the weighting + if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) { + psError(PS_ERR_UNKNOWN, false, "Can't measure mean variance for image."); + psFree(output); + return false; + } + *weighting = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) * + psImageCovarianceFactor(readout->covariance)); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "PPSTACK.WEIGHTING", 0, + "Weighting by 1/noise^2 for stack", *weighting); + psFree(bg); Index: /trunk/ppStack/src/ppStackReadout.c =================================================================== --- /trunk/ppStack/src/ppStackReadout.c (revision 21476) +++ /trunk/ppStack/src/ppStackReadout.c (revision 21477) @@ -10,6 +10,4 @@ #include "ppStack.h" -//#define TESTING // Write debugging output? - bool ppStackReadoutInitialThread(psThreadJob *job) { @@ -20,10 +18,9 @@ pmConfig *config = args->data[1]; // Configuration pmReadout *outRO = args->data[2]; // Output readout - psArray *subRegions = args->data[3]; // Regions for PSF-matching - psArray *subKernels = args->data[4]; // Kernels for PSF-matching - psVector *weightings = args->data[5]; // Weightings (1/noise^2) for each image - psVector *addVariance = args->data[6]; // Additional variance when rejecting - - psArray *inspect = ppStackReadoutInitial(config, outRO, thread->readouts, subRegions, subKernels, + psVector *mask = args->data[3]; // Mask for inputs + psVector *weightings = args->data[4]; // Weightings (1/noise^2) for each image + psVector *addVariance = args->data[5]; // Additional variance when rejecting + + psArray *inspect = ppStackReadoutInitial(config, outRO, thread->readouts, mask, weightings, addVariance); @@ -42,9 +39,11 @@ pmConfig *config = args->data[1]; // Configuration pmReadout *outRO = args->data[2]; // Output readout - psArray *rejected = args->data[3]; // Rejected pixels - psVector *weightings = args->data[4]; // Weighting (1/noise^2) for each image - - bool status = ppStackReadoutFinal(config, outRO, thread->readouts, rejected, - weightings); // Status of operation + psVector *mask = args->data[3]; // Mask for inputs + psArray *rejected = args->data[4]; // Rejected pixels + psVector *weightings = args->data[5]; // Weighting (1/noise^2) for each image + psVector *addVariance = args->data[6]; // Weighting (1/noise^2) for each image + + bool status = ppStackReadoutFinal(config, outRO, thread->readouts, mask, rejected, + weightings, addVariance); // Status of operation thread->busy = false; @@ -85,14 +84,10 @@ psArray *ppStackReadoutInitial(const pmConfig *config, pmReadout *outRO, const psArray *readouts, - const psArray *regions, const psArray *kernels, const psVector *weightings, - const psVector *addVariance) + const psVector *mask, const psVector *weightings, const psVector *addVariance) { assert(config); assert(outRO); assert(readouts); - assert(regions); - assert(kernels); - assert(readouts->n == regions->n); - assert(regions->n == kernels->n); + assert(mask && mask->n == readouts->n && mask->type.type == PS_TYPE_VECTOR_MASK); assert(weightings && weightings->n == readouts->n && weightings->type.type == PS_TYPE_F32); assert(addVariance && addVariance->n == readouts->n && addVariance->type.type == PS_TYPE_F32); @@ -107,4 +102,5 @@ 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 @@ -123,5 +119,5 @@ for (int i = 0; i < num; i++) { pmReadout *ro = readouts->data[i]; - if (!ro) { + if (!ro || mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) { // Bad image continue; @@ -148,21 +144,10 @@ } - if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, kernelSize, iter, combineRej, combineSys, - true, useVariance, safe)) { + if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, kernelSize, iter, + combineRej, combineSys, combineDiscard, true, useVariance, safe)) { psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection."); psFree(stack); return false; } - -#ifdef TESTING - { - psString name = NULL; // Name of image - psStringAppend(&name, "combined_initial_%03d.fits", sectionNum); - psFits *fits = psFitsOpen(name, "w"); - psFree(name); - psFitsWriteImage(fits, NULL, outRO->image, 0, NULL); - psFitsClose(fits); - } -#endif // Save list of pixels to inspect @@ -189,5 +174,6 @@ bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, const psArray *readouts, - const psArray *rejected, const psVector *weightings) + const psVector *mask, const psArray *rejected, const psVector *weightings, + const psVector *addVariance) { assert(config); @@ -196,4 +182,5 @@ assert(rejected); assert(readouts->n == rejected->n); + assert(mask && mask->n == readouts->n && mask->type.type == PS_TYPE_VECTOR_MASK); assert(weightings && weightings->n == readouts->n && weightings->type.type == PS_TYPE_F32); @@ -217,12 +204,10 @@ int numGood = num; // Number of good inputs: images that haven't been completely rejected for (int i = 0; i < num; i++) { - if (!rejected->data[i]) { + pmReadout *ro = readouts->data[i]; + if (!ro || !rejected->data[i] || mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) { // Image completely rejected numGood--; continue; } - - pmReadout *ro = readouts->data[i]; - assert(ro); #if 0 @@ -242,10 +227,10 @@ } - pmStackData *data = pmStackDataAlloc(ro, weightings->data.F32[i], NAN); + pmStackData *data = pmStackDataAlloc(ro, weightings->data.F32[i], addVariance->data.F32[i]); data->reject = psMemIncrRefCounter(rejected->data[i]); stack->data[i] = data; } - if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, 0, 0, NAN, NAN, + if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, 0, 0, NAN, NAN, NAN, numGood != num, useVariance, false)) { psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts."); @@ -254,15 +239,4 @@ } -#ifdef TESTING - { - psString name = NULL; // Name of image - psStringAppend(&name, "combined_final_%03d.fits", sectionNum); - psFits *fits = psFitsOpen(name, "w"); - psFree(name); - psFitsWriteImage(fits, NULL, outRO->image, 0, NULL); - psFitsClose(fits); - } -#endif - pmCell *outCell = outRO->parent; // Output cell pmChip *outChip = outCell->parent; // Output chip Index: /trunk/ppStack/src/ppStackThread.c =================================================================== --- /trunk/ppStack/src/ppStackThread.c (revision 21476) +++ /trunk/ppStack/src/ppStackThread.c (revision 21477) @@ -49,5 +49,5 @@ ppStackThreadData *ppStackThreadDataSetup(const psArray *cells, const psArray *imageNames, const psArray *maskNames, const psArray *varianceNames, - const pmConfig *config) + const psArray *covariances, const pmConfig *config) { PS_ASSERT_ARRAY_NON_NULL(cells, NULL); @@ -55,4 +55,5 @@ PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, maskNames, NULL); PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, varianceNames, NULL); + PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, covariances, NULL); ppStackThreadData *stack = psAlloc(sizeof(ppStackThreadData)); // Thread data, to return @@ -98,5 +99,7 @@ continue; } - readouts->data[j] = pmReadoutAlloc(cell); + pmReadout *ro = pmReadoutAlloc(cell); // Readout for thread + ro->covariance = psMemIncrRefCounter(covariances->data[j]); + readouts->data[j] = ro; } threads->data[i] = ppStackThreadAlloc(readouts); @@ -236,5 +239,5 @@ { - psThreadTask *task = psThreadTaskAlloc("PPSTACK_INITIAL_COMBINE", 7); + psThreadTask *task = psThreadTaskAlloc("PPSTACK_INITIAL_COMBINE", 6); task->function = &ppStackReadoutInitialThread; psThreadTaskAdd(task); @@ -250,5 +253,5 @@ { - psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 5); + psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 7); task->function = &ppStackReadoutFinalThread; psThreadTaskAdd(task);