IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 21477


Ignore:
Timestamp:
Feb 13, 2009, 1:53:56 PM (17 years ago)
Author:
Paul Price
Message:

Using covariance factor to set the correct variance level for the combination. Adding additional recipe parameter to provide the fraction of values to discard when performing a 'olympic weighted mean'.

Location:
trunk/ppStack/src
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppStack/src/ppStack.h

    r21366 r21477  
    4646                                          const psArray *maskNames, // Names of masks to read
    4747                                          const psArray *varianceNames, // Names of variance maps to read
     48                                          const psArray *covariances, // Covariance matrices (already read)
    4849                                          const pmConfig *config // Configuration
    4950    );
    … …  
    9798                               pmReadout *outRO,   // Output readout
    9899                               const psArray *readouts, // Input readouts
    99                                const psArray *regions, // Array with array of regions used in each PSF match
    100                                const psArray *kernels, // Array with array of kernels used in each PSF match
     100                               const psVector *mask, // Mask for input readouts
    101101                               const psVector *weightings, // Weighting factors for each image
    102102                               const psVector *addVariance // Additional variance for rejection
    … …  
    115115                         pmReadout *outRO,   // Output readout
    116116                         const psArray *readouts, // Input readouts
     117                         const psVector *mask, // Mask for input readouts
    117118                         const psArray *rejected, // Array with pixels rejected in each image
    118                          const psVector *weightings // Weighting factors for each image
     119                         const psVector *weightings, // Weighting factors for each image
     120                         const psVector *addVariance // Additional variance for rejection
    119121    );
    120122

  • trunk/ppStack/src/ppStackArguments.c

    r21366 r21477  
    147147    psMetadataAddF32(arguments, PS_LIST_TAIL, "-combine-sys", 0,
    148148                     "Relative systematic error in combination", NAN);
     149    psMetadataAddF32(arguments, PS_LIST_TAIL, "-combine-discard", 0,
     150                     "Discard fraction for Olympic weighted mean", NAN);
    149151    psMetadataAddStr(arguments, PS_LIST_TAIL, "-mask-val", 0, "Mask value of input bad pixels", NULL);
    150152    psMetadataAddStr(arguments, PS_LIST_TAIL, "-mask-bad", 0, "Mask value to give bad pixels", NULL);
    … …  
    232234    }
    233235
    234     VALUE_ARG_RECIPE_INT("-iter",             "ITER",           S32, 0);
    235     VALUE_ARG_RECIPE_FLOAT("-combine-rej",    "COMBINE.REJ",    F32);
    236     VALUE_ARG_RECIPE_FLOAT("-combine-sys",    "COMBINE.SYS",    F32);
    237     VALUE_ARG_RECIPE_FLOAT("-threshold-mask", "THRESHOLD.MASK", F32);
    238     VALUE_ARG_RECIPE_FLOAT("-image-rej",      "IMAGE.REJ",      F32);
    239     VALUE_ARG_RECIPE_FLOAT("-deconv-limit",   "DECONV.LIMIT",   F32);
    240     VALUE_ARG_RECIPE_INT("-rows",             "ROWS",           S32, 0);
    241     VALUE_ARG_RECIPE_FLOAT("-poor-frac",      "POOR.FRACTION",  F32);
     236    VALUE_ARG_RECIPE_INT("-iter",              "ITER",            S32, 0);
     237    VALUE_ARG_RECIPE_FLOAT("-combine-rej",     "COMBINE.REJ",     F32);
     238    VALUE_ARG_RECIPE_FLOAT("-combine-sys",     "COMBINE.SYS",     F32);
     239    VALUE_ARG_RECIPE_FLOAT("-combine-discard", "COMBINE.DISCARD", F32);
     240    VALUE_ARG_RECIPE_FLOAT("-threshold-mask",  "THRESHOLD.MASK",  F32);
     241    VALUE_ARG_RECIPE_FLOAT("-image-rej",       "IMAGE.REJ",       F32);
     242    VALUE_ARG_RECIPE_FLOAT("-deconv-limit",    "DECONV.LIMIT",    F32);
     243    VALUE_ARG_RECIPE_INT("-rows",              "ROWS",            S32, 0);
     244    VALUE_ARG_RECIPE_FLOAT("-poor-frac",       "POOR.FRACTION",   F32);
    242245
    243246    valueArgRecipeStr(arguments, recipe, "-mask-val",  "MASK.VAL",  recipe);

  • trunk/ppStack/src/ppStackLoop.c

    r21376 r21477  
    480480            psMetadataAddF32(stats, PS_LIST_TAIL, "STAMP.NUM", PS_META_DUPLICATE_OK,
    481481                             "Number of stamps", kernels->numStamps);
    482 
     482            psMetadataAddF32(stats, PS_LIST_TAIL, "PPSTACK.WEIGHTING", PS_META_DUPLICATE_OK,
     483                             "Weighting for image", weightings->data.F32[i]);
    483484        }
    484485        psLogMsg("ppStack", PS_LOG_INFO, "Time to match image %d: %f sec", i, psTimerClear("PPSTACK_MATCH"));
    … …  
    495496        writeImage(maskNames->data[i], maskHeader, readout->mask, config);
    496497        psFree(maskHeader);
     498        psImageCovarianceTransfer(readout->variance, readout->covariance);
    497499        writeImage(varianceNames->data[i], hdu->header, readout->variance, config);
     500#ifdef TESTING
     501        {
     502            psString name = NULL;
     503            psStringAppend(&name, "covariance_%d.fits", i);
     504            writeImage(name, hdu->header, readout->covariance->image, config);
     505            psFree(name);
     506        }
     507#endif
    498508
    499509        pmCell *inCell = readout->parent; // Input cell
    … …  
    605615                             i, matchChi2->data.F32[i]);
    606616                } else {
     617                    psLogMsg("ppStack", PS_LOG_INFO, "Image %d has matching chi^2: %f",
     618                             i, matchChi2->data.F32[i]);
    607619                    numGood++;
    608620                }
    … …  
    630642    // Start threading
    631643    ppStackThreadInit();
    632     ppStackThreadData *stack = ppStackThreadDataSetup(cells, imageNames, maskNames, varianceNames, config);
     644    ppStackThreadData *stack = ppStackThreadDataSetup(cells, imageNames, maskNames, varianceNames,
     645                                                      covariances, config);
    633646    if (!stack) {
    634647        psError(PS_ERR_IO, false, "Unable to initialise stack threads.");
    … …  
    729742            psArrayAdd(job->args, 1, config);
    730743            psArrayAdd(job->args, 1, outRO);
    731             psArrayAdd(job->args, 1, subRegions);
    732             psArrayAdd(job->args, 1, subKernels);
     744            psArrayAdd(job->args, 1, inputMask);
    733745            psArrayAdd(job->args, 1, weightings);
    734746            psArrayAdd(job->args, 1, matchChi2);
    … …  
    761773            return false;
    762774        }
    763         psFree(matchChi2);
    764775
    765776        // Harvest the jobs, gathering the inspection lists
    … …  
    787798        memDump("initial");
    788799    }
     800
     801#ifdef TESTING
     802    writeImage("combined_initial.fits", NULL, outRO->image, config);
     803#endif
    789804
    790805    if (stats) {
    … …  
    828843                psFree(rejected);
    829844                psFree(covariances);
     845                psFree(matchChi2);
    830846                return false;
    831847            }
    … …  
    844860            psFree(rejected);
    845861            psFree(covariances);
     862            psFree(matchChi2);
    846863            return false;
    847864        }
    … …  
    939956            psFree(outRO);
    940957            psFree(covariances);
     958            psFree(matchChi2);
    941959            return false;
    942960        }
    … …  
    969987                psFree(outRO);
    970988                psFree(covariances);
     989                psFree(matchChi2);
    971990                return false;
    972991            }
    … …  
    9811000            psArrayAdd(job->args, 1, config);
    9821001            psArrayAdd(job->args, 1, outRO);
     1002            psArrayAdd(job->args, 1, inputMask);
    9831003            psArrayAdd(job->args, 1, rejected);
    9841004            psArrayAdd(job->args, 1, weightings);
     1005            psArrayAdd(job->args, 1, matchChi2);
    9851006            if (!psThreadJobAddPending(job)) {
    9861007                psFree(job);
    … …  
    9911012                psFree(outRO);
    9921013                psFree(covariances);
     1014                psFree(matchChi2);
    9931015                return false;
    9941016            }
    … …  
    10041026            psFree(outRO);
    10051027            psFree(covariances);
     1028            psFree(matchChi2);
    10061029            return false;
    10071030        }
    … …  
    10091032
    10101033    memDump("final");
     1034
     1035#ifdef TESTING
     1036    writeImage("combined_final.fits", NULL, outRO->image, config);
     1037#endif
    10111038
    10121039    // Sum covariance matrices
    … …  
    10191046    outRO->covariance = psImageCovarianceSum(covariances);
    10201047    psFree(covariances);
     1048    psFree(matchChi2);
    10211049
    10221050    if (stats) {

  • trunk/ppStack/src/ppStackMatch.c

    r21366 r21477  
    170170    assert(kernels && !*kernels);
    171171    assert(config);
     172    *weighting = 0.0;
    172173
    173174    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
    … …  
    260261        pmSubtractionAnalysis(readout->analysis, kernels, region,
    261262                              readout->image->numCols, readout->image->numRows);
     263
     264        psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel
     265        psKernel *covar = psImageCovarianceCalculate(kernel, readout->covariance); // New covariance matrix
     266        psFree(readout->covariance);
     267        readout->covariance = covar;
     268        psFree(kernel);
     269
    262270    } else {
    263271#endif
    … …  
    541549        }
    542550        psFree(iter);
    543         *chi2 /= num;
     551        *chi2 /= psImageCovarianceFactor(readout->covariance) * num;
    544552    }
    545553
    … …  
    579587        (void)psBinaryOp(readout->image, readout->image, "-",
    580588                         psScalarAlloc(psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), PS_TYPE_F32));
    581         *weighting = 1.0 / PS_SQR(psStatsGetValue(bg, PS_STAT_ROBUST_STDEV));
    582         psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "PPSTACK.WEIGHTING", 0,
    583                          "Weighting by 1/noise^2 for stack",
    584                          1.0 / PS_SQR(psStatsGetValue(bg, PS_STAT_ROBUST_STDEV)));
    585     }
     589    }
     590
     591    // Measure the variance level for the weighting
     592    if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) {
     593        psError(PS_ERR_UNKNOWN, false, "Can't measure mean variance for image.");
     594        psFree(output);
     595        return false;
     596    }
     597    *weighting = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) *
     598                        psImageCovarianceFactor(readout->covariance));
     599    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "PPSTACK.WEIGHTING", 0,
     600                     "Weighting by 1/noise^2 for stack", *weighting);
     601
    586602    psFree(bg);
    587603

  • trunk/ppStack/src/ppStackReadout.c

    r21183 r21477  
    1010#include "ppStack.h"
    1111
    12 //#define TESTING                  // Write debugging output?
    13 
    1412bool ppStackReadoutInitialThread(psThreadJob *job)
    1513{
    … …  
    2018    pmConfig *config = args->data[1];   // Configuration
    2119    pmReadout *outRO = args->data[2];   // Output readout
    22     psArray *subRegions = args->data[3]; // Regions for PSF-matching
    23     psArray *subKernels = args->data[4]; // Kernels for PSF-matching
    24     psVector *weightings = args->data[5]; // Weightings (1/noise^2) for each image
    25     psVector *addVariance = args->data[6]; // Additional variance when rejecting
    26 
    27     psArray *inspect = ppStackReadoutInitial(config, outRO, thread->readouts, subRegions, subKernels,
     20    psVector *mask = args->data[3];     // Mask for inputs
     21    psVector *weightings = args->data[4]; // Weightings (1/noise^2) for each image
     22    psVector *addVariance = args->data[5]; // Additional variance when rejecting
     23
     24    psArray *inspect = ppStackReadoutInitial(config, outRO, thread->readouts, mask,
    2825                                             weightings, addVariance);
    2926
    … …  
    4239    pmConfig *config = args->data[1];   // Configuration
    4340    pmReadout *outRO = args->data[2];   // Output readout
    44     psArray *rejected = args->data[3];  // Rejected pixels
    45     psVector *weightings = args->data[4];  // Weighting (1/noise^2) for each image
    46 
    47     bool status = ppStackReadoutFinal(config, outRO, thread->readouts, rejected,
    48                                       weightings); // Status of operation
     41    psVector *mask = args->data[3];     // Mask for inputs
     42    psArray *rejected = args->data[4];  // Rejected pixels
     43    psVector *weightings = args->data[5];  // Weighting (1/noise^2) for each image
     44    psVector *addVariance = args->data[6];  // Weighting (1/noise^2) for each image
     45
     46    bool status = ppStackReadoutFinal(config, outRO, thread->readouts, mask, rejected,
     47                                      weightings, addVariance); // Status of operation
    4948
    5049    thread->busy = false;
    … …  
    8584
    8685psArray *ppStackReadoutInitial(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
    87                                const psArray *regions, const psArray *kernels, const psVector *weightings,
    88                                const psVector *addVariance)
     86                               const psVector *mask, const psVector *weightings, const psVector *addVariance)
    8987{
    9088    assert(config);
    9189    assert(outRO);
    9290    assert(readouts);
    93     assert(regions);
    94     assert(kernels);
    95     assert(readouts->n == regions->n);
    96     assert(regions->n == kernels->n);
     91    assert(mask && mask->n == readouts->n && mask->type.type == PS_TYPE_VECTOR_MASK);
    9792    assert(weightings && weightings->n == readouts->n && weightings->type.type == PS_TYPE_F32);
    9893    assert(addVariance && addVariance->n == readouts->n && addVariance->type.type == PS_TYPE_F32);
    … …  
    107102    float combineRej = psMetadataLookupF32(NULL, recipe, "COMBINE.REJ"); // Combination threshold
    108103    float combineSys = psMetadataLookupF32(NULL, recipe, "COMBINE.SYS"); // Combination systematic error
     104    float combineDiscard = psMetadataLookupF32(NULL, recipe, "COMBINE.DISCARD"); // Olympic discard fraction
    109105    bool useVariance = psMetadataLookupBool(&mdok, recipe, "VARIANCE"); // Use variance for rejection?
    110106    bool safe = psMetadataLookupBool(&mdok, recipe, "SAFE"); // Be safe when combining small numbers of pixels
    … …  
    123119    for (int i = 0; i < num; i++) {
    124120        pmReadout *ro = readouts->data[i];
    125         if (!ro) {
     121        if (!ro || mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
    126122            // Bad image
    127123            continue;
    … …  
    148144    }
    149145
    150     if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, kernelSize, iter, combineRej, combineSys,
    151                         true, useVariance, safe)) {
     146    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, kernelSize, iter,
     147                        combineRej, combineSys, combineDiscard, true, useVariance, safe)) {
    152148        psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection.");
    153149        psFree(stack);
    154150        return false;
    155151    }
    156 
    157 #ifdef TESTING
    158     {
    159         psString name = NULL;           // Name of image
    160         psStringAppend(&name, "combined_initial_%03d.fits", sectionNum);
    161         psFits *fits = psFitsOpen(name, "w");
    162         psFree(name);
    163         psFitsWriteImage(fits, NULL, outRO->image, 0, NULL);
    164         psFitsClose(fits);
    165     }
    166 #endif
    167152
    168153    // Save list of pixels to inspect
    … …  
    189174
    190175bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, const psArray *readouts,
    191                          const psArray *rejected, const psVector *weightings)
     176                         const psVector *mask, const psArray *rejected, const psVector *weightings,
     177                         const psVector *addVariance)
    192178{
    193179    assert(config);
    … …  
    196182    assert(rejected);
    197183    assert(readouts->n == rejected->n);
     184    assert(mask && mask->n == readouts->n && mask->type.type == PS_TYPE_VECTOR_MASK);
    198185    assert(weightings && weightings->n == readouts->n && weightings->type.type == PS_TYPE_F32);
    199186
    … …  
    217204    int numGood = num;                  // Number of good inputs: images that haven't been completely rejected
    218205    for (int i = 0; i < num; i++) {
    219         if (!rejected->data[i]) {
     206        pmReadout *ro = readouts->data[i];
     207        if (!ro || !rejected->data[i] || mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
    220208            // Image completely rejected
    221209            numGood--;
    222210            continue;
    223211        }
    224 
    225         pmReadout *ro = readouts->data[i];
    226         assert(ro);
    227212
    228213#if 0
    … …  
    242227        }
    243228
    244         pmStackData *data = pmStackDataAlloc(ro, weightings->data.F32[i], NAN);
     229        pmStackData *data = pmStackDataAlloc(ro, weightings->data.F32[i], addVariance->data.F32[i]);
    245230        data->reject = psMemIncrRefCounter(rejected->data[i]);
    246231        stack->data[i] = data;
    247232    }
    248233
    249     if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, 0, 0, NAN, NAN,
     234    if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, 0, 0, NAN, NAN, NAN,
    250235                        numGood != num, useVariance, false)) {
    251236        psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts.");
    … …  
    254239    }
    255240
    256 #ifdef TESTING
    257     {
    258         psString name = NULL;           // Name of image
    259         psStringAppend(&name, "combined_final_%03d.fits", sectionNum);
    260         psFits *fits = psFitsOpen(name, "w");
    261         psFree(name);
    262         psFitsWriteImage(fits, NULL, outRO->image, 0, NULL);
    263         psFitsClose(fits);
    264     }
    265 #endif
    266 
    267241    pmCell *outCell = outRO->parent;    // Output cell
    268242    pmChip *outChip = outCell->parent;  // Output chip

  • trunk/ppStack/src/ppStackThread.c

    r21377 r21477  
    4949ppStackThreadData *ppStackThreadDataSetup(const psArray *cells, const psArray *imageNames,
    5050                                          const psArray *maskNames, const psArray *varianceNames,
    51                                           const pmConfig *config)
     51                                          const psArray *covariances, const pmConfig *config)
    5252{
    5353    PS_ASSERT_ARRAY_NON_NULL(cells, NULL);
    … …  
    5555    PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, maskNames, NULL);
    5656    PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, varianceNames, NULL);
     57    PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, covariances, NULL);
    5758
    5859    ppStackThreadData *stack = psAlloc(sizeof(ppStackThreadData)); // Thread data, to return
    … …  
    9899                continue;
    99100            }
    100             readouts->data[j] = pmReadoutAlloc(cell);
     101            pmReadout *ro = pmReadoutAlloc(cell); // Readout for thread
     102            ro->covariance = psMemIncrRefCounter(covariances->data[j]);
     103            readouts->data[j] = ro;
    101104        }
    102105        threads->data[i] = ppStackThreadAlloc(readouts);
    … …  
    236239
    237240    {
    238         psThreadTask *task = psThreadTaskAlloc("PPSTACK_INITIAL_COMBINE", 7);
     241        psThreadTask *task = psThreadTaskAlloc("PPSTACK_INITIAL_COMBINE", 6);
    239242        task->function = &ppStackReadoutInitialThread;
    240243        psThreadTaskAdd(task);
    … …  
    250253
    251254    {
    252         psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 5);
     255        psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 7);
    253256        task->function = &ppStackReadoutFinalThread;
    254257        psThreadTaskAdd(task);
Note: See TracChangeset for help on using the changeset viewer.